Imaginary time evolution with MPS using apply function.

Hi!

I am extremely sorry if this sounds like a repetition of the question here. I am trying to find out the finite-T state of an MPS using ancilla approach. I am using the XXZ model .

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

I am aware of the MPO approach in ITensor.jl repo & it works fine (I have done the calculations using ED & they match).But,I am insterested in the MPS approach.

My trotter gates for imaginary time evolution is defined as,

# 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), (τ=-δτ_im / 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))

The imaginary time evolution looks like

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)




# Do imaginary time evolution to obtain thermal state @ inverse temp

ψ_β=ψ_0
println("--------------------------ψ_beta initial--------------------")
println(ψ_β)


for β in 0:δτ_im:beta_max
    
    
    energy=inner(ψ_β',H_physical,ψ_β)
    
    println("β = ",β, " energy = ",energy," Maxdim ",maxlinkdim(ψ_β)," overlap ",inner(ψ_β,ψ_0)," corr ",correlation_matrix(ψ_β,"Sz","Sz")[1,3])
    flush(stdout)
    ψ_β= apply(gates_imag_timeevo,ψ_β; maxdim=maxdim,cutoff=cutoff)
    noprime!(ψ_β)
   
    ψ_β=(1.0/norm(ψ_β))*ψ_β

end

I have used the same approach as real time evolution with just tau turned to imaginary. But, the energy values are way off the ED result(& result when done with MPO). Does apply function has any limitation for non -nearest neighbour interactions like this (the matrix is also not Hermitian)? The initial state is an infinite temperature singlet state (it is separately verified).

I can’t figure out at all what is going wrong here. The same gate set gives exact result for MPO approach. Is there any way to do it using MPS?

I would appreciate any help.

Note: \psi _0 is the infinite temperature state. The function used for construction is given below if required.

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

Your code and approach looks very reasonable, and looking over it I can’t see anything obviously wrong. Most likely there is just a subtle bug somewhere.

I would start with your initial state construction: after making this initial state does it have the correct properties? Can you verify that the entanglement is log(2) on odd bonds and zero on even bonds? Is the expected value of the physical energy what it should be?

If it helps, you could just start from a totally “empty” MPS and not worry about the link indices. You can just set the MPS tensors such as to make singlets on each pair of sites, then call orthogonalize!(psi) at the end and it will put in the link indices for you at the end. It could simplify your code a good bit.

Then if that part seems to work, are you sure your mapping from imaginary time step to temperature is correct? (Not off by a factor of 2?) I think it is correct but good to check these things.

Finally, I would verify certain other limits such as, does the code eventually give the correct ground state (T=0) energy for large enough imaginary time evolution? If not that could help you find the bug.

Miles

P.S. note that only the state at T=0 is the “ground state”. All other states you are computing for finite T are just finite temperature density matrices and not ground states.

1 Like

Really sorry for late reply. Thanks a lot @miles for the helpful clues. The time step to temperature was off by a factor of 2 exactly as you indicated.The results now match with ED quite well.

Just another maybe related question regarding apply function. It seems too slow to apply non-local gates (next nearest neighbour). I guess it internally swaps the qubits’ position before applying the gate & this swapping contributes to increasing linkdim. Is there any trick to make it faster? (like using MPO. I guess I can take e^{-iHt}=1-iHt+O(t^2 ) to construct the MPO as I am not sure whether it is possible to automatically construct the MPO for propagator from the Hamiltonian.This wouldn’t rewuire swapping)

Thank you for your kind help.

Hi!
I am trying to find out the finite T state of an MPS using ancilla approach for electronic system. As a beginner of TD-DMRG, I don not understand how to write this program, although I have already read the relevant papers. Could you tell me how to derive the thermal state of electronic system or why the matrix A is A=\begin{pmatrix}0&\frac{1}{\sqrt{2} } \\-\frac{1}{\sqrt{2} }&0\end{pmatrix}

Hi @evan,
Glad to hear you are learning about the ancilla method and TD-DMRG. However, a few things:

  • please post new questions as new topics, rather than at the end of an existing topic
  • regarding your question, please as much as possible try to frame your questions in terms of specific ITensor code you are trying to write, after you have given it a try yourself and attempted to debug it
  • also regarding your question about the ancilla method, I would encourage you to take some time to teach it to yourself by working some small examples (this is what I would do, and you’d find it very valuable I think). So you can write the entire state of two spins (one physical, one ‘ancila’) starting from a singlet (infinite temperature), then act with exp(-beta H/2) which is a matrix you can fully work out, then act with this on your state, trace the ancilla with two copies of the state, and verify that you get the correct thermal density matrix for the physical spin

Thanks for your suggestions.