Segmentation Fault in Calculation of Excited States?

Question

Hi Miles,

I had a error when running the following code to calculate the exited state of a spin-1 chain,

int N = 24;
auto sites = SpinOne(N);
Real J = 0.05;
auto ampo = AutoMPO(sites);
for(int j = 1; j <= N; ++j)
    {
    ampo += 0.5,"Sz2",j;
    ampo += -0.5,"Sz",j;
    }
for(int j = 1; j < N; ++j)
    {
    ampo += -J,"S+",j,"S-",j+1;
    ampo += -J,"S+",j+1,"S-",j;
    }
auto H = IQMPO(ampo);

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

auto state = InitState(sites);

for(int i = 1; i <= N; i++)
{
    state.set(i, i <= N/2 ? "+" : "0");
}

auto psi0 = IQMPS(state);

auto en0 = dmrg(psi0,H,sweeps,{"Quiet=",true});

auto wfs = std::vector<IQMPS>(1);
wfs.at(0) = psi0;

auto psi1 = IQMPS(state);

auto en1 = dmrg(psi1,H,wfs,sweeps,{"Quiet=",true,"Weight=",20.0});

It always stops after calculating the ground state, with a message "Segmentation fault: 11". When I run it in debug mode, it says

"
From line 87, file itensor_operators.cc

div(T1)=QN(24) must equal div(T2)=QN() when adding T1+T2

div(T1)=QN(24) must equal div(T2)=QN() when adding T1+T2
Abort trap: 6
"

This error looks wired for me because the code works for N <=22. I solved this problem by modifying the initial state

for(int i = 1; i <= N; i++)
{
    state.set(i, i%2 == 1 ? "+" : "0");
}

My questions,

1. Why the former initial state would cause this error?

2. Why this error appear only for N >= 24?

Many thanks for answering all these questions!

Jin

Answer 1

Hi Jin,
For an IQMPS to be used it must be fully initialized. The first InitState you made did not specify a state for all of the sites but only the first N/2 of them. So the result is an IQMPS that is not fully defined.

I agree though that this behavior is somewhat confusing and that there should be a clearer error message that happens earlier. I'll see if I can make that happen.

Miles

Comments

Hi Miles,

Why isn't my first InitState fully initialized? I think the first one is "1...10...0" (the first N/2 sites are "1" and the other N/2 sites are "0"), and the second one is "1010...10". Right? Thanks.

Jin

---

Oh my fault I misread your code. I had better try this code myself then to see what is going on. It may be a bug. Thanks for correcting -

---

Hi Jin,
I took a closer look at the issue you were having today and it does seem to be a bug in ITensor, which I just pushed a patch for. The sample code you provided above works quite well now. Please give it a try and let me know if it doesn't work still for some reason.

What I think was happening is that for whatever reason, psi0 was becoming perfectly orthogonal to psi1 before the second call to DMRG. So when the code was computing the overlap of psi1 and psi0 during the DMRG algorithm, the result was a zero tensor which for efficiency the code was representing as an IQTensor with no blocks (no storage data) whatsoever. This should be fine in that the code could just treat this as zero, but instead it was trying to add two IQTensors with different numbers of blocks and encountering a bug. So now if you add an IQTensor with zero blocks to an existing one it just does nothing and moves on which fixes the issue at least for this case (and hopefully doesn't cause any other new issues).

Best,
Miles

---

By the way, I added the following lines at the end to check the accuracy of the results and to confirm that both psi0 and psi1 are eigenstates by checking their energy variance:

    printfln("<psi0|psi1> = %.5E",overlap(psi0,psi1));
    
    auto var0 = overlap(psi0,H,H,psi0) - en0*en0;
    auto var1 = overlap(psi1,H,H,psi1) - en1*en1;
    
    printfln("var0 = %.5E",var0);
    printfln("var1 = %.5E",var1);

---

Oh and the reason changing the initial state sometimes fixed this issue is that it might have slightly changed things so that psi1 wasn't perfectly orthogonal to psi0 when starting the second call to DMRG.

---

Hi Miles,

Thank you so much. Yes, the sample code works quite well now. And thanks for the explanation, very helpful!

Jin