Hi there! I am trying to implement local basis optimization after this paper. I’ve done my best to follow their algorithm (figure 2 is very helpful in this respect), but I am definitely doing something wrong – my GS energy for an 8-site Hubbard Holstein system is on the order of +25, while I know that it should be closer to -2. Here is the gist of what I am doing:
function dmrg(...)
# Randomly initialize a rotation matrix (R) on every site
# All other pre-DMRG steps, like calling orthogonalize! and position! on site 1
for sw in 1:nsweep(sweeps)
for (b, ha) in sweepnext(N)
if do_LBO
# Perform LBO (see below)
end
# Perform the remainder of 2-site DMRG using ITensors code
end
end
Now, I will get into more detail and include snippets of actual code. Here is how I am randomly initializing my rotation matrices:
if LBO
Rs = []
for i in 1:N
# Start off with having the same dimension
S_ind = getfirst(x -> hastags(x, "Site"), inds(psi[i]))
Trunc_ind = settags(copy(S_ind),"Trunc")
R = randomITensor(prime(S_ind), dag(Trunc_ind))
push!(Rs,R)
end
end
And here is my LBO implementation:
## HELPER: Get the left indices of the tensor we want to SVD##
function get_Linds(O::ITensor, tag::String)
Rind = getfirst(x -> hastags(x, tag), inds(O))
return noncommoninds([Rind], inds(O))
end
# Put into canonical form
orthogonalize!(psi,b)
M = psi[b]
R = Rs[b]
H = PH.H[b]
# Perform the first SVD
X, Λ, Y = svd(M, get_Linds(M, "Site"), lefttags="svd1,u", righttags="svd1,v")
u = getfirst(x -> hastags(x, "svd1,u"), inds(Λ))
s = getfirst(x -> hastags(x, "Site"), inds(Y))
t = getfirst(x -> hastags(x, "Trunc"), inds(R))
# Make and optimize R̃
R̃ = R * Y * Λ * dag(delta(t,s))
X̃, Λ̃, R = svd(R̃, u, lefttags="svd2,u", righttags="Trunc")
# Perform updates
t = getfirst(x -> hastags(x, "Trunc"), inds(R))
if dim(t) >= 4 # Keep at least the minimum Hubbard dimension
# Make the new M̃ and H̃
M̃ = Λ̃ * X̃ * X
H̃ = dag(prime(R,t)) * H * noprime(R)
# Rearrange/relabel some indices to get back to starting arrangement
replacetags!(M̃, "Trunc", tags(s))
replacetags!(H̃, "Trunc", tags(s))
s_old = getfirst(x -> hastags(x, "Site"), inds(R))
s_new = getfirst(x -> hastags(x, "Site"), inds(H̃))
R = R * delta(dag(s_old), s_new)
# Update
psi[b] = M̃
Rs[b] = R
PH.H[b] = H̃
end
I know that my relabelling/ rearranging of indices might be a bit inelegant, but I am not sure if this is behind my problem. Would love to hear any thoughts you might have on this issue! Thanks so much, and have a great week.