MPS after imaginary time evolution.

Hi!

I am trying to get a finite temperature state using imaginary time evolution. One thing to note that my Hamiltonian is next nearest neighbour. I have followed the example in ITensor.jl/purification closely for the test. But after imaginary time evolution,the MPS I get have all the link indices with tag “Link,n=1” irrespective of their actual site. The id though is correct. Here is a small snippet to reproduce the issue.

using ITensors

function ITensors.op(::OpName"expτXXZ", ::SiteType"S=1/2", s1::Index, s2::Index; τ,J,Δ)
  
    h =
        J * op("Sx", s1) * op("Sx", s2) +
        J * op("Sy", s1) * op("Sy", s2) +
        (J*Δ)*op("Sz", s1) * op("Sz", s2)
      return exp(τ * h)
  end

function inf_temp_mps(sites)

    num_sites=length(sites)
    if(num_sites %2 !=0)
        throw(DomainError(num_sites,"Expects even number of sites for ancilla-physical singlets."))
    else
        
        ψ=MPS(sites)
        for j=1:2:num_sites-1

            
            s1=sites[j]
            s2=sites[j+1]
            
            if(j==1)
                

                rightlink=commonind(ψ[j+1],ψ[j+2])
                A=ITensor(ComplexF64 ,s1,s2,rightlink )

                A[s1=>1,s2=>2,rightlink=>1]= 1/sqrt(2)
                A[s1=>2,s2=>1,rightlink=>1]= -1/sqrt(2)

                U,S,V=svd(A,(s1),cutoff=1e-15)
                ψ[j]= replacetags(U,"u", "l=$(j)")
                ψ[j+1]= replacetags(S*V,"u", "l=$(j)")

            elseif(j==num_sites-1)
                
                leftlink =commonind(ψ[j-1],ψ[j])
                A=ITensor(ComplexF64 ,s1,s2,leftlink )

                A[s1=>1,s2=>2,leftlink =>1]= 1/sqrt(2)
                A[s1=>2,s2=>1,leftlink =>1]= -1/sqrt(2)

                U,S,V=svd(A,(s1,leftlink ),cutoff=1e-15)
                ψ[j]= replacetags(U,"u", "l=$(j)")
                ψ[j+1]= replacetags(S*V,"u", "l=$(j)")
            else
                

                rightlink=commonind(ψ[j+1],ψ[j+2])
                leftlink =commonind(ψ[j-1],ψ[j])
            
                A=ITensor(ComplexF64 ,s1,s2,rightlink,leftlink )

                A[s1=>1,s2=>2,rightlink=>1,leftlink =>1]= 1/sqrt(2)
                A[s1=>2,s2=>1,rightlink=>1,leftlink =>1]= -1/sqrt(2)

                U,S,V=svd(A,(s1,leftlink ),cutoff=1e-15)
                ψ[j]= replacetags(U,"u", "l=$(j)")
                ψ[j+1]= replacetags(S*V,"u", "l=$(j)")
            end
        end

        return ψ

    end
end

# Configure system
δτ=0.1# trotter step
# Model parameters
J=0.5
Δ=0.1 
#Inverse temp
beta_max=2.0
N=6
s=siteinds("S=1/2",N)
ψ=inf_temp_mps(s)

os = OpSum()
for j in 1:2:(N - 2)
    
    os += J,"Sx",j,"Sx",j+2
    os += J,"Sy",j,"Sy",j+2
    os += J*Δ,"Sz",j,"Sz",j+2
    
end
#Physical hamiltonian
H_physical = MPO(os, s)

# Create gates for imaginary time evolution (Apply only on physical sites)

  # Make gates (1,3),(3,5),(5,7),...
  gates_imag_timeevo = ops([("expτXXZ", (n, n + 2), (τ=-δτ / 2,J,Δ)) for n in 1:2:(N - 2)], s)
  # Include gates in reverse order to complete Trotter formula
  append!(gates_imag_timeevo, reverse(gates_imag_timeevo))

for β in 0:δτ:beta_max
  
  energy=inner(ψ',H_physical,ψ)
  push!(energy_psi,energy)

  @printf("β = %.2f energy = %.8f\n", β, energy)
  ψ= apply(gates_imag_timeevo, ψ; cutoff=1e-15)
  normalize!(ψ)
end
# Check the link index tags in the output ψ
println(ψ)

The output is given below.You may notice the MPS printed with multiple linkindices having tag for “Link,n=1”

β = 0.00 energy = 0.00000000
β = 0.10 energy = -0.01260422
β = 0.20 energy = -0.02526980
β = 0.30 energy = -0.03796322
β = 0.40 energy = -0.05065010
β = 0.50 energy = -0.06329558
β = 0.60 energy = -0.07586472
β = 0.70 energy = -0.08832296
β = 0.80 energy = -0.10063654
β = 0.90 energy = -0.11277288
β = 1.00 energy = -0.12470100
β = 1.10 energy = -0.13639189
β = 1.20 energy = -0.14781876
β = 1.30 energy = -0.15895735
β = 1.40 energy = -0.16978611
β = 1.50 energy = -0.18028634
β = 1.60 energy = -0.19044226
β = 1.70 energy = -0.20024107
β = 1.80 energy = -0.20967287
β = 1.90 energy = -0.21873063
β = 2.00 energy = -0.22741003
MPS
[1] ((dim=2|id=198|"S=1/2,Site,n=1"), (dim=2|id=105|"Link,n=1"))
[2] ((dim=2|id=105|"Link,n=1"), (dim=2|id=450|"S=1/2,Site,n=2"), (dim=4|id=667|"Link,n=1"))
[3] ((dim=4|id=667|"Link,n=1"), (dim=2|id=834|"S=1/2,Site,n=3"), (dim=8|id=923|"Link,n=1"))
[4] ((dim=2|id=656|"S=1/2,Site,n=4"), (dim=4|id=522|"Link,n=1"), (dim=8|id=923|"Link,n=1"))
[5] ((dim=2|id=790|"S=1/2,Site,n=5"), (dim=2|id=302|"Link,l=5"), (dim=4|id=522|"Link,n=1"))
[6] ((dim=2|id=395|"S=1/2,Site,n=6"), (dim=2|id=302|"Link,l=5"))

They are associated with proper sites as it seems from the link ids.

To be precise, my even sites are ancilla & odd sites are physical qubits.
I guess using apply for non-consecutive sites have some restrictions in C++. I checked the apply function definition in Julia version & the docstring suggests this has been taken care of using movesite. I am not sure what is going wrong here. Though as the link ids are correct,I guess the results should not get affected in subsequent calculations.

I would really appreciate any suggestion.

Hi, so I believe the issue here is your initial state being all down (“Dn”) spins. If you think about the action of the S^x_i S^x_j + S^y_i S^y_j term as being equal to S^+_i S^-_j + S^-_i S^+_j (here I am leaving out a factor of 1/2), then this term tries to flip a pair of spins. But if both spins in the pair are down, it gives zero. So your initial state will exactly be an eigenstate of your Hamiltonian and it can’t act in a non-trivial way.

Please try a different initial state. In fact if you are trying to do the purification method, the initial state should be a product state of entangled “Bell pairs” between the even and odd spins.

Hope that helps!

Thanks a lot for the response.
Sorry @miles for this . I have now updated the code with infinite temperature state & the issue persists.

Ok glad you tried that. Could you be more specific though about what you changed?

Also note that we have an updated purification code in the ITensor source repo which uses an MPO now instead of an MPS. This makes the gate pattern simpler since the gates can be nearest-neighbor and just act on one side of the MPO. (If the MPO starts as an identity operator it is equivalent to an MPS of Bell pairs.)

Thank you for the suggestion. In the code snippet above now ,I have an MPS made of Bell pairs between physical & auxiliary sites now (which represents an infinite temperature state). I do an imaginary time evolution with apply function (the Hamiltonian acts only on physical sites). The state that comes out have wrong tags for link indices in bulk sites. All of the has “Link,n=1” .(The output is now added in the post).

Thanks for the reference to the new purification code. It works great but I was interested in MPS approach first.

Really appreciate your kind support.

Oh I see– the issue was with the tags? In this case the tags of the virtual indices of an MPS are not really so important… the important thing is that the dimensions of the indices are growing as they should automatically througout the time evolution process.

Are the dimensions growing to values greater than 1?

Yeah, the bond dimension grows as they should & link ids are alright. The issue was with the tags only.Anyway, that doesn’t seem to affect any calculation.

Thanks a lot for helping with this.Have a nice day.

Thanks for the quick & helpful feedback and discussion. It sounds like this is a minor bug with the apply function and its tagging behavior – we will look into fixing that, thanks!