Operators with complex coefficients

edited

Hello,

I am trying to obtain the ground state of a Hubbard model (I am adapting the "exthubbard" example) where the hopping terms are given by

for(int b = 1; b < N; ++b)
{
ampo += -t1exp(ii(-0.5)chi),"Cdagup",b,"Cup",b+1;
ampo += -t1
exp(ii(+0.5)chi),"Cdagup",b+1,"Cup",b;
ampo += -t1exp(ii(+0.5)chi),"Cdagdn",b,"Cdn",b+1;
ampo += -t1
exp(ii(-0.5)chi),"Cdagdn",b+1,"Cdn",b;
}

where

ii is the imaginary constant, and chi and t1 are real constants.

The initial wavefunction MPS is given like in the example (Neel state). Upon running I get

libc++abi.dylib: terminating with uncaught exception of type std::runtime_error: Error condition in diagHermitian

Is there any change I should make to the MPOs to allow for these complex numbers?

Best regards,

Rafael

commented by (70.1k points)
Hi Rafael, I tried to reproduce your bug but I was unable. Could you check that you get the error consistently, and if so, post your entire code in a comment below?

Here's the code I created to attempt to reproduce, and it ran without any errors

int N = 10;
int Npart = N; //half filling
auto sites = Electron(N,{"ConserveQNs=",true});

auto ii = Complex_i;
auto t1 = 1.0;
auto chi = 0.4;
auto U = 0.1;

auto ampo = AutoMPO(sites);
for(int b = 1; b < N; ++b)
{
ampo += -(cos(0.5*chi)+ii*sin(-0.5*chi)),"Cdagup",b,"Cup",b+1;
ampo += -(cos(0.5*chi)+ii*sin(+0.5*chi)),"Cdagup",b+1,"Cup",b;
ampo += -(cos(0.5*chi)+ii*sin(+0.5*chi)),"Cdagdn",b,"Cdn",b+1;
ampo += -(cos(0.5*chi)+ii*sin(-0.5*chi)),"Cdagdn",b+1,"Cdn",b;
}
for(int i = 1; i <= N; ++i)
{
ampo += U,"Nupdn",i;
}
auto H = toMPO(ampo);

auto state = InitState(sites);
int p = Npart;
for(int i = N; i >= 1; --i)
{
if(p > i)
{
println("Doubly occupying site ",i);
state.set(i,"UpDn");
p -= 2;
}
else
if(p > 0)
{
println("Singly occupying site ",i);
state.set(i,(i%2==1 ? "Up" : "Dn"));
p -= 1;
}
else
{
state.set(i,"Emp");
}
}

auto psi0 = MPS(state);

Print(totalQN(psi0));

auto sweeps = Sweeps(5);
sweeps.maxdim() = 10,20,100,100,200;
sweeps.cutoff() = 1E-10;
sweeps.niter() = 2;
sweeps.noise() = 1E-7,1E-8,0.0;
println(sweeps);

//
// Begin the DMRG calculation
//
auto [energy,psi] = dmrg(H,psi0,sweeps,{"Quiet",true});
commented by (360 points)
Hello Miles,

Thank you for the reply. I think the issue was a silly definition mistake in the imaginary factor. I will try to run your version and I will write back if there are any further issues.

Best wishes,

Rafael