I am trying to create MPOs (long-range Hamiltonians) via their finite state automaton representations. The hamiltonian I have is given by:
$$H = - \sum_{ij, s} c^{\dagger}{i,s} c{j,s} + V_0 \sum_{i,j>i, s, s’} \frac{e^{-r/D}}{\sqrt{1 + r^2/d^2}} n_{i,s} n_{j,s}$$
I understand that there are two ways to create this MPO in ITensor (1) using the OpSum() method or (2) define it as a linear combination of exponential decays and then define each tensor for each site. The later is supposed to speed things up. Since I am studying the above model for different (D,d) I will end to create many such MPOs. The problem I am facing is that the above MPO created in the two methods yield different groundstates via DMRG. I am not sure what’s causing the issues. I am comparing the two MPOs and they seem to be the same.
I implement the first method as:
function construct_HMPO_screened(sites, D, d, V_int_strength)
N = length(sites); os = OpSum()
for j=1:N-1
os += -1,"Cdagup",j,"Cup",j+1
os += -1, "Cdagdn",j,"Cdn",j+1
os += -1, "Cdagup",j+1,"Cup",j
os += -1, "Cdagdn",j+1,"Cdn",j
end
for i in 1:N-1
for j in i+1:N
os += V_int_strength*V_int(abs(i-j), D, d),"Ntot",i,"Ntot",j
end
end
return MPO(os,sites)
end
N=240
sites = siteinds("Electron",conserve_nf=true,conserve_sz=true, N);
V_int_strength = 0.29
D = 30; d = 4
H_MPO = construct_HMPO_screened(sites, D, d, V_int_strength)
nelectrons = 24
### create initial state
init_state = ["Emp" for i in 1:N]
for i in 1:nelectrons
if i % 2 == 0
init_state[2i] = "Up"
elseif i%2 == 1
init_state[end - 2i+1] = "Dn"
end
end
psi_init = MPS(sites, init_state)
### perform DMRG
energy, psi_WC = dmrg(H_MPO,psi_init;nsweeps=12, maxdim=[650], cutoff=1e-10,
noise=[1e-8, 1e-10], eigsolve_krylovdim=16)
When I create the MPO via its Finite state Automaton framework, I first estimate the parameters \{ \phi_k, \lambda_k\} such that (\frac{1}{N} \sum_{k=1}^N \left | 1 - \frac{\frac{e^{-r/D}}{\sqrt{1 + r^2/d^2}}}{\sum_k \phi_k e^{-\lambda_k r_k}}\right|) < 0.01 where N=4. This gives me the parameters:
best_params = [0.03968434, 0.11529795, 0.28910056, 0.73098612], [0.03656105, 0.0539406, 0.09868307, 0.22186] for \{ \phi_k, \lambda_k\}.
Then, I create the FSA as follows:
function get_HMPO_manybody(sites, V_int_strength, best_params)
other_params = [V_int_strength,length(sites)]
L = length(sites);
# define QNs for each column of FSA
# First column is initial state
# then is followed by the 4 (length of best_params[1]) columns of interaction strengths and
# the sixth column is the final state
# the last 4 columns are the Kinetic energy terms
αα = [Index(
QN(("Nf", 0, -1), ("Sz", 0)) => 2+n_phi,
QN(("Nf", 1, -1), ("Sz", 1)) => 1,
QN(("Nf", -1, -1), ("Sz", -1)) => 1,
QN(("Nf", 1, -1), ("Sz", -1)) => 1,
QN(("Nf", -1, -1), ("Sz", 1)) => 1) for j = 1:L+1]
A = [single_mpo_tensor_manybody(ss[j], αα[j], dag(αα[j+1]), best_params, j, other_params) for j = 1:L]
A[1] *= bcL_manybody(dag(αα[1]), initial_state_index, total_bonddim )
A[end] *= bcR_manybody(αα[end], final_state_index, total_bonddim)
H_MPO = exp_mpo_manybody(sites, best_params, other_params)
return MPO(A)
end
function single_mpo_tensor_manybody(s :: Index, αL :: Index, αR :: Index, best_params, index_num, other_params)
U, nsites = other_params
initial_state_index = 1
n_phi = length(best_params[1])
final_state_index = 2+n_phi
KE_term1_index = 3+n_phi
KE_term2_index = 4+n_phi
KE_term3_index = 5+n_phi
KE_term4_index = 6+n_phi
# Initial and final sites
A = onehot(αL => initial_state_index, αR => initial_state_index) * op("I",s)
A += onehot(αL => final_state_index, αR => final_state_index) * op("I",s)
# KE ; spin up
# notation --> K(state_start, state_end) = operator
A += onehot(αL => initial_state_index, αR => KE_term1_index) * -1 * op("Cdagup",s)
A += onehot(αL => KE_term1_index, αR => final_state_index) * op("Cup",s)
# hermitian conjugate
A += onehot(αL => initial_state_index, αR => KE_term2_index) * -1 * op("Cup",s)
A += onehot(αL => KE_term2_index, αR => final_state_index) * op("Cdagup",s)
# KE ; spin down
A += onehot(αL => initial_state_index, αR => KE_term3_index) * -1 * op("Cdagdn",s)
A += onehot(αL => KE_term3_index, αR => final_state_index) * op("Cdn",s)
# hermitian conjugate
A += onehot(αL => initial_state_index, αR => KE_term4_index) * -1 * op("Cdn",s)
A += onehot(αL => KE_term4_index, αR => final_state_index) * op("Cdagdn",s)
# Coloumb interaction
for i in 1:length(n_phi)
A += onehot(αL => initial_state_index, αR => i+1) * U * best_params[1][i] * op("Ntot", s)
A += onehot(αL => i+1, αR => final_state_index) * exp(-best_params[2][i]) * op("Ntot", s)
A += onehot(αL => i+1, αR => i+1) * exp(-best_params[2][i]) * op("I",s)
end
return ITensor(A |> array, inds(A))
end
bcL_manybody(α, wanted_dim, total_dim) = (@assert dim(α) == total_dim; onehot(α => wanted_dim))
bcR_manybody(α, wanted_dim, total_dim) = (@assert dim(α) == total_dim; onehot(α => wanted_dim))
## this is called before DMRG as follows:
H_MPO = get_HMPO_manybody(sites, V_int_strength, best_params)
However, both these MPOs yield different ground states. The right answer is given by the first method, where I expect spins to dimerize, whereas in the later case, the electrons never dimerize. Where are am I going wrong? Is their something wrong with the MPO creation method in approach 2 ?
For sanity checks, I checked the overlap of the two MPOs with different states. The following seems to be a little different:
state3 = [(i % 2 == 0) ? "Up" : "Dn" for i in 1:240]
psi3 = MPS(siteids_to_use, state3)
inner(psi3', H_MPO1, psi3) |> println # gives 591.9147815238581
inner(psi3', H_MPO2, psi3) |> println # gives 592.8440939055168
Thanks in advance!