Bug in MPOs with products of operators? (Julia version)

Question

Hi,

I have codes in which I perform DMRG and then compute < S^2 >.
Before I updated ITensor to the latest version, these codes were giving correct results and now they are not. I know that the previous results were correct because I benchmarked with simple models for which I can do exact diagonalization.

One simple example that may help to identify this apparent bug is the following:

Model: DMRG for S=1 antiferromagnetic Heisenberg dimer.
The groundstate should be a singlet.
I compute < Sz > and I get 0. That's good.
I compute < S^2 > and I get different than 0. That's not right.

I defined the S^2 operator as follows:

function S2_op(Nsites,sites)
ampo = AutoMPO()

for i in 1:Nsites
    for j in 1:Nsites
        ampo += 1.,"S+",i,"S-",j
        ampo += 1.,"Sz",i,"Sz",j
    end
    ampo += -1.,"Sz",i
end

S2 = MPO(ampo,sites)

return S2
end

I also tried to define it in a different manner: (EDIT: this is wrong)

function S2_op_other(Nsites,sites)
ampo = AutoMPO()

for i in 1:Nsites
    ampo += 1.,"S+*S-",i
    ampo += 1.,"Sz*Sz",i
    ampo += -1.,"Sz",i
end

S2 = MPO(ampo,sites)

return S2
end

Both cases give values that are not 0. Also, they give different results, which I think is also telling. In the first case, I get 2/3. In the second case, I get 4. None of these numbers makes any physical sense.

The exact same thing happened when I tried to create an MPO of the Sz^2 operator and compute its average.

Best,
Gonçalo Catarina

Comments

I am quite sure that you have entered the Hamiltonian incorrectly in the second case (and possibly in the first case too, although I'm not as sure about that one). In addition, the two are different from each other.

The double for loop makes your first Hamiltonian put an interaction term between EVERY pair of sites, not just adjacent ones. Is that really what you want? Because I'm pretty sure you only want intensely nonlocal interactions in very specific cases, and I thought the Heisenberg dimer wasn't one of them. Also, generally ITensor doesn't work too great on models with a lot of nonlocal interactions (I even asked about this once: see http://itensor.org/support/2051/what-governs-runtime-and-memory-usage-of-h-mpo-ampo-sites), although that doesn't mean it's not worth trying. I will also note that, when I ran DMRG on this model on my computer using Julia, I indeed got a ground state energy extremely close to zero, which you are claiming is what you should have gotten.

In contrast, your second Hamiltonian is actually the opposite--every term in the Hamiltonian acts on just one site, so you can actually just minimize the energy of each site independently. In fact, it is just a diagonal matrix, because every term is diagonal. In the spin-1 Hilbert space, S+S- is just diag(2, 2, 0), SzSz is just diag(1, 0, 1), and Sz is just diag(1, 0, -1). So your total Hamiltonian is just diag(2, 2, 2) on each site, which is literally just 2 times the identity operator! So really every state is the ground state, and you should get a ground state energy of 2*N, where N is your number of sites. I can confirm that I did indeed get this result when I ran DMRG on this Hamiltonian in Julia.

Ultimately, I assume you are referencing some book/paper/website. If so, a reference would be extremely helpful, that we we can see exactly what equation your're using.

---

Hi,

Thanks a lot for the reply.

Let me try to clarify my point.
1) My codes are untouched. Before the latest update, they were giving different outputs. Of course, this doesn't necessarily mean there is a bug now.
2) What I do is to run DMRG with an Hamiltonian H, obtaining the groundstate wavefuntion psi0. Then, I want to compute <psi0 | S^2 |psi0>. In order to do so, I create an MPO of the S^2 operator and compute inner(psi0,S2,psi0).
3) The S^2 operator can be written as S^2 = Sp*Sm + Sz^2 - Sz. Of course, I need to sum over all sites, so that Sp = sum_i Sp(i), Sm = sum_i Sm(i) and Sz = sum_i Sz(i). Back in the day, it was not possible to multiply autoMPOs, so I was using the first option (check this reply from Miles, http://itensor.org/support/546/possible-multiply-autompo-objects-before-converting-iqmpos, where I got this trick from). The second option that I sent is *wrong*. My bad, I am really sorry! You are correct, these are not the same operators and I messed it up.
4) Below, I send a code where I compute <psi0 | S^2 |psi0> for the S=1 antiferromagnetic Heisenberg dimer and I do not obtain 0 (EDIT: it should be 0 because the groundstate is a singlet). I also compute <psi0 | Sz | psi0> and I obtain 0, as expected. However, I also compute <psi0 | Sz^2 |psi0> (using the same trick) and I do not obtain 0, which is weird.

Code:

using ITensors

# Hamiltonian: Heisenberg, first neighbor
function H_Heis_1n(J,Nsites,sites)
    ampo = AutoMPO()
    
    #J
    for i in 1:Nsites-1
        ampo += -J,"Sz",i,"Sz",i+1
        ampo += -J/2,"S+",i,"S-",i+1
        ampo += -J/2,"S+",i+1,"S-",i
    end
    
    H = MPO(ampo,sites)
    
    return H
end

# S^2 operator
function S2_op(Nsites,sites)
    ampo = AutoMPO()
    
    for i in 1:Nsites
        for j in 1:Nsites
            ampo += 1.,"S+",i,"S-",j
            ampo += 1.,"Sz",i,"Sz",j
        end
        ampo += -1.,"Sz",i
    end
    
    S2 = MPO(ampo,sites)
    
    return S2
end

# Sz operator
function Sz_op(Nsites,sites)
    ampo = AutoMPO()
    
    for i in 1:Nsites
        ampo += 1.,"Sz",i
    end
    
    Sz = MPO(ampo,sites)
    
    return Sz
end

# Sz^2 operator
function Sz2_op(Nsites,sites)
    ampo = AutoMPO()
    
    for i in 1:Nsites
        for j in 1:Nsites
            ampo += 1.,"Sz",i,"Sz",j
        end
    end
    
    Sz2 = MPO(ampo,sites)
    
    return Sz2
end

# physical parameters
N = 2
J = -1. #energy in units of |J|

# DMRG parameters
sweeps0 = Sweeps(10)
maxdim!(sweeps0, 20,50,100,100,200)
cutoff!(sweeps0, 1E-5,1E-5,1E-5,1E-8,1E-8)

# system
Nsites = N
sites = siteinds("S=1",Nsites)

# initial state
ψi = randomMPS(sites)

# Hamiltonian
H = H_Heis_1n(J,Nsites,sites)

# run DMRG
## groundstate
E0,ψ0 = dmrg(H,ψi,sweeps0)

# S^2 operator
S2 = S2_op(Nsites,sites)

# Sz operator
Sz = Sz_op(Nsites,sites)

# Sz^2 operator (TEST)
Sz2 = Sz2_op(Nsites,sites)

# outputs
println("S^2 for groundstate = ", inner(ψ0,S2,ψ0))
println("Sz for groundstate = ", inner(ψ0,Sz,ψ0))
println("Sz^2 for groundstate = ", inner(ψ0,Sz2,ψ0))

---

OK, your comment helps a lot. I was clearly very confused about exactly what you were asking, so thank you for the clarification!

I ran the above code on my computer for N = 2, 10, and 20, and in every case I indeed got <S^2>, <Sz>, and <Sz^2> very close to zero. This is reassuring, since you said this is what you should get, but this also means it's not helping me see what may have gone wrong.

I did indeed copy exactly the code you put above, with absolutely no modifications, so I imagine any discrepancy must be due to differences in Julia. My local copy of ITensors.jl was last updated on June 9, and it was the most recent version at the time. Which version were you using before the update, and which one did you switch to? If you let me know, we might be able to localize what changed that made things different.

---

Hi,

Again, thanks a lot for the reply!

I agree with you.
Indeed, for chains with even N you should get a singlet groundstate and therefore all of these average values equal to 0 (the Sz^2 is silly to compute, it was just a benchmark test that I did on my own).

At the moment, I have the latest version of ITensor, where this issue is still present.
I noticed this discrepancy in the previous version of ITensor: I updated it on Friday (June 19) and then these problems started.
However, I had not updated it for a couple of weeks, so that I cannot tell you exactly in which version this discrepancy appeared.
This puts the occurrence of this discrepancy between June 9 and June 19.

---

Hi Goncalo,
What's the reason for the ampo += -1.,"Sz",i term in the first definition? I'm not sure that's correct but I only just did a pretty quick derivation so I could have missed something too.

Best,
Miles

---

Hi Miles,

The S^2 = Sx^2 + Sy^2 + Sz^2 operator can be written as S^2 = Sp*Sm + Sz^2 - Sz.
You can check it, for instance, in this website, https://inst.eecs.berkeley.edu/~cs191/sp05/lectures/lecture10.pdf (page 3, Eq.1).
The derivation is straightforward. The -Sz term comes from [Sp,Sm] = 2Sz. The key step is to verify that Sx^2 + Sy^2 = 1/2*(Sp*Sm + Sm*Sp) = Sp*Sm - Sz.

Best,
Gonçalo Catarina

---

Thanks - yes I see then & it's a good way to write it.

When you say that you get "different from zero" for <S^2>, how different? I wouldn't be surprised at all if it's not numerically exactly zero. But how large is it?

---

I get <S^2> = 2/3.
I also get <Sz^2> = 2/3.
You can check this out by copying the code I sent above and run it in the latest version of ITensor Julia.

In previous versions, I was obtaining values very close to 0.

If it helps, in the comments above we seem to have understood that the occurrence of this discrepancy happened in some update between June 9 and June 19.

---

That does help - thanks for pointing that out again about the results changing for (apparently) the same ground state MPS. You know, I did make a pretty significant change to the AutoMPO code around that very time (it was on June 12) so it could very much be related to what you're observing. I'll look into it - I'm filing an issue right now so it doesn't get overlooked.

---

Just to let you know, there's definitely a bug in AutoMPO that's responsible for part of what you're seeing. I'm working on it and will let you know when that's done.

---

Thanks for the update!

Answer 1

This bug should be fixed now in version 0.1.14 of Julia ITensor (ITensors.jl). Please let me know if it isn't! Thanks for pointing out this issue.

Best,
Miles