Hi, Ryan. Thanks so much for looking over the code. It was silly of me to compute the purity instead of the trace. I believe the trace error is now fixed. Thanks also for directing me to the entropy tools GitHub. It was very useful. I lifted a piece of code from one of your modules.
I am still unsure if the reduced density matrix I obtain from using your code instead of outer(.,.) is the correct one. Maybe I am making a mistake that I can’t seem to spot? I have computed the average energy of S in two ways. (1) from the purified \psi_{\rm{pure}} as \rm{Tr} [\bra{\psi_{\rm{pure}}} H \ket{\psi_{\rm{pure}}}] where H acts only S's degrees of freedom and (2) from the reduced density matrix as \rm{Tr} [\rho_S H]. Both these energies are not equal, shouldn’t they be?
Here is a sample code, where H is the Hamiltonian for long-range TFIM.
N = 4 # Number of spins
J = 1.0 # Coupling strength
h = 1.0 # Transverse field
alpha = 2.0 #exponential decay strength
cutoff = 1e-10
#-----------------------------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------------------------
function get_local_bell_state(N::Int,sites)
@assert iseven(length(sites))
N = length(sites)
ψ0 = MPS(sites)
for j in 1:2:N
s1 = sites[j]
s2 = sites[j+1]
t = ITensor(s1,s2)
t[s1 => "Up", s2 => "Up",] = 1/sqrt(2)
t[s1 => "Dn", s2 => "Dn",] = 1/sqrt(2)
U,S,V = svd(t,s1)
ψ0[j] = U*S
ψ0[j+1] = V
end
return ψ0
end
function getH(N, J, h, alpha, sites)
os = OpSum()
#Interaction energy
for j in 1:2:2N-3
for k in j+2:2:2N
dist = abs(j-k)
coupling = -J / (dist^alpha)
os += coupling, "X", j, "X", k
end
end
#Transverse field
for j in 1:2:2N
os += -h, "Z", j
end
H = MPO(os, sites);
return H
end
function getH_sys_sites(N, J, h, alpha, sys_sites)
os = OpSum()
#Interaction energy
for j in 1:N-1
for k in j+1:N
dist = abs(j-k)
coupling = -J / (dist^alpha)
os += coupling, "X", j, "X", k
end
end
#Transverse field
for j in 1:N
os += -h, "Z", j
end
H = MPO(os, sys_sites);
return H
end
#----------------------------------------------------------------------------------
#----------------------------------------------------------------------------------
sites = siteinds("S=1/2", 2N; conserve_qns=false)
sys_sites = sites[1:2:2N]
ψ = get_local_bell_state(2N,sites)
H = getH(N, J, h, alpha, sites)
# TDVP parameters
nsteps = 5
δβ = 0.1
β = 2.0
num_β_steps = Int(round(β/δβ))
β_0 = 0.0
sweeps = Sweeps(5)
maxdim!(sweeps, 64)
cutoff!(sweeps, 1e-10)
println("---------- TDVP Evolution ----------")
ψ1 = nothing
for bstep in 0:num_β_steps
β_current = β_0 + bstep * δβ
ψ1 = ψ # reset or reuse initial state
τ = -(β_current) / (2nsteps)
for n in 1:nsteps
ψ1 = tdvp(
H,
τ,
ψ1;
updater=exponentiate_updater,
nsweeps=5,
maxdim=64,
cutoff=1e-10,
normalize=true
)
end
end
#Ryan's code to get the reduced density matrix of S
function density_matrix_sites(ψ_::AbstractMPS, region;)
"""
Obtain a density matrix, leaving some indices uncontracted
Assumes very little about tag structure, only site/linkinds
Currently assume ordered and orthogonality center is within region
"""
region = sort(region)
start, stop = region[1], region[end]
s = siteinds(ψ_)
ψ = orthogonalize(ψ_, start)
if length(region) == 1
return ψ[start] * prime(dag(ψ[start]), s[start])
end
ψH = dag(prime(ψ, linkinds(ψ)..., s[region]...))
lᵢ₁ = linkinds(ψ, start - 1) #this is n,n+1 link
ρ = prime(ψ[start], lᵢ₁)
ρ *= ψH[start]
for k in (start + 1):(stop - 1) # replace this with ITensorNetworks iterator
ρ *= ψ[k]
ρ *= ψH[k]
end
lⱼ = linkinds(ψ, stop)
ρ *= prime(ψ[stop], lⱼ)
ρ *= ψH[stop]
return ρ
end
ψ_pure= ψ1
E_pure_state = inner(ψ_pure',H,ψ_pure)
sites_to_keep = collect(1:2:2N)
ρ_S_p = density_matrix_sites(ψ_pure, sites_to_keep)
ρ_S = MPO(ρ_S_p,sys_sites)
H_sys = getH_sys_sites(N,J,h,alpha,sys_sites)
E_rdm = inner(H_sys,ρ_S)
println("energy from rdm:", E_rdm)
println("energy from pure state:", E_pure_state)
println("diff in energy:", E_rdm-E_pure_state)
Thanks again for your help!
