Summation of MPS

ITensor C++ Questions
1

Hello, I am new to DMRG and ITensor package. Currently I am very interested to calculate the spectral function based on Lanzcos method as described in this article.(https://www.cond-mat.de/events/correl11/manuscript/Koch.pdf) . However, when I run the following program.I found that the summation of H\ket{\psi} and \ket\psi are incompatibal with quantum number conserved. I print out the MPS and found H acting on the state actually doesn’t change the index structure but an error was reported.

int main()
{
    int L = 4;
    double U = 1.0;
    double J = 1.0;
    double g = 3.0;
    double K = 5.0;


    auto sites = SpinHalf(L,{"ConserveQNs=",true});

    auto ampo = AutoMPO(sites);
    for(int j = 1; j <L; ++j)
        {
        ampo += 0.5,"S+",j,"S-",j+1;
        ampo += 0.5,"S-",j,"S+",j+1;
        ampo +=     "Sz",j,"Sz",j+1;
        }
    auto H = toMPO(ampo);

    auto sweeps = Sweeps(30); //number of sweeps is 5
    sweeps.maxdim() = 10,20,100,100,200;
    sweeps.cutoff() = 1E-10;

    auto state = InitState(sites);
    for(auto i : range1(L))
    {
        if(i%2 == 1) state.set(i,"Up");
        else         state.set(i,"Dn");
    }
    itensor::MPS psi0 = MPS(state);
   
    auto [energy,psi] = dmrg(H,psi0,sweeps,{"Quite",true});

    itensor::MPS tmp=itensor::applyMPO(H,psi);
    std::cout<<"++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
    println(tmp);
    
    std::cout<<"++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
    println(psi);
    
    std::cout<<"++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
    tmp=itensor::sum(tmp,psi);
    println(tmp);
    std::cout<<"++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
}

And after I turn off the QN conservation, the error was not reported. I can’t find where the mistake is. Please help me . Plus, I an really looking forward to calculate spectrum of some quantum systems, I am grateful if any suggestions can be provided on the possible algoirthm and C++ implementation.

2

Hi thanks for the question. I found the issue with the code you were trying: the H operator has two sets of site indices: the ones with no primes which match the sites of your MPS and then a second set of site indices which have a single prime on them (a “prime level” of one). (This is our standard convention for operators in ITensor.)

So when you do tmp = applyMPO(H,psi) the unprimed indices of H contract with psi, leaving the primed indices. Then before you sum tmp with psi you have to make sure they have the same indices with the same prime level. So you can do this:

tmp.noPrime();

which will set all prime levels to zero.

Then you’ll see that the summation will now work afterward.

I’m not sure why it worked in the non-QN case but not in the QN conserving case. It should have errored for both!

3

Hi, Miles. Thank you for replying my silly question. But acually this problem remains when I add noPrime to my code. The program works well when ConserveQNs=false and segment fault when I turn ConserveQNs=true. ConserveQNs is much more efficient and I want to keep it.

int main()
{
int L = 4;
double U = 1.0;
double J = 1.0;
double g = 3.0;
double K = 5.0;
int r1=L/2+1;//zero point to avoid boundry effect
auto sites = SpinHalf(L,{“ConserveQNs=”,true});

auto ampo = AutoMPO(sites);
for(int j = 1; j <L; ++j)
    {
    ampo += 0.5,"S+",j,"S-",j+1;
    ampo += 0.5,"S-",j,"S+",j+1;
    ampo +=     "Sz",j,"Sz",j+1;
    }
auto H = toMPO(ampo);

auto sweeps = Sweeps(30); //number of sweeps is 5
sweeps.maxdim() = 10,20,100,100,200;
sweeps.cutoff() = 1E-10;

auto state = InitState(sites);
for(auto i : range1(L))
{
    if(i%2 == 1) state.set(i,"Up");
    else         state.set(i,"Dn");
}
itensor::MPS psi0 = MPS(state);

auto [energy,psi] = dmrg(H,psi0,sweeps,{"Quite",true});

for(int i=1;i<=L;i++)
{
    auto psiSz1=psi;
    auto psiSz2=psi;

    psiSz1.position(i);
    auto tmp=op(sites,"Sz",i)*psiSz1(i);
    tmp.noPrime();
    psiSz1.set(i,tmp);

    psiSz2.position(r1);
    tmp=op(sites,"Sz",r1)*psiSz2(r1);
    tmp.noPrime();
    psiSz2.set(r1,tmp);

    psiSz1=itensor::sum(psiSz1,psiSz2);

}

}
It’s really strange and may be related to some deeper principles in the code which I don’t understand. Thanks a lot for your help.