0 votes
asked by (360 points)
edited by


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
ampo += -t1exp(ii(+0.5)chi),"Cdagdn",b,"Cdn",b+1;
ampo += -t1


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,


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);
            p -= 2;
        if(p > 0)
            println("Singly occupying site ",i);
            state.set(i,(i%2==1 ? "Up" : "Dn"));
            p -= 1;

    auto psi0 = MPS(state);


    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;

    // 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,


1 Answer

+1 vote
answered by (70.1k points)

(See answer to question in comments above)

Welcome to ITensor Support Q&A, where you can ask questions and receive answers from other members of the community.

Formatting Tips:
  • To format code, indent by four spaces
  • To format inline LaTeX, surround it by @@ on both sides
  • To format LaTeX on its own line, surround it by $$ above and below
  • For LaTeX, it may be necessary to backslash-escape underscore characters to obtain proper formatting. So for example writing \sum\_i to represent a sum over i.
If you cannot register due to firewall issues (e.g. you cannot see the capcha box) please email Miles Stoudenmire to ask for an account.

To report ITensor bugs, please use the issue tracker.