Real time evolution of Hubbard chain

ITensor Julia Questions
,
1

Dear professors,

I test a very simple Hubbard chain with ITensor.
Firstly, I calculate the ground state and energy.
And then I calculate the real time evolution of it.
However, I meet a mistake. During the time evolution, the total energy is not conserved.
In fact, it should not be changes.

using ITensors, ITensorMPS, DelimitedFiles
let
  N = 10
  Npart = 10
  t = -0.5
  U = 4.0
  sites = siteinds("Electron", N, conserve_qns=true)

#################
  os = OpSum()
  for j=1:N-1
     os += -t, "Cdagup", j, "Cup", j + 1;
     os += -t, "Cdagup", j + 1, "Cup", j;
     os += -t, "Cdagdn", j, "Cdn", j + 1;
     os += -t, "Cdagdn", j + 1, "Cdn", j;
  end
  for j=1:N
     os += U, "Nupdn", j;
  end
  H = MPO(os , sites)







#######################################################################
#Set the initial wavefunction matrix product state to be a Neel state.
p = Npart;
state = Vector{String}(undef, N);
for i in range(start=N, step=-1, stop=1)
if p > i
println("Doubly occupying site ",i);
state[i] = "UpDn";
p -= 2;
elseif p > 0
println("Singly occupying site ",i);
state[i] = (isodd(i) ? "Up" : "Dn");
p -= 1;
else
state[i] = "Emp";
end
end








#######################################################################
  nsweeps = 50
  maxdim = [10,20,100,100,200]
  cutoff = [1E-10]
  noise = [1E-7,1E-8,0.0]
  psi0_init = random_mps(sites , state ; linkdims=2)
  energy0 , psi0 = dmrg(H,psi0_init;nsweeps,maxdim,cutoff,noise)
  println("ground energy = ", energy0, ", total magnetization = ", sum(expect(psi0, "Sz")))




Cuu = correlation_matrix(psi0, "Nup", "Ndn")
println("Cuu = ", round.(Cuu,digits=3))
sum_diag = sum(diag(Cuu))
println("double occupation = ", sum_diag) 


#######################################################################
  cutoff = 1E-8
  tau = 0.01
  ttotal = 10.0
  U2 = 0


#######################################################################
  gates = ITensor[]
  for j in 1:(N - 1)
    s1 = sites[j]
    s2 = sites[j + 1]
    hj = - t * op("Cdagup", s1) * op("Cup", s2)  - t * op("Cdagup", s2) * op("Cup", s1)  - t * op("Cdagdn", s1) * op("Cdn", s2)  - t * op("Cdagdn", s2) * op("Cdn", s1)
    Gj = exp(-im * tau / 2 * hj)
    push!(gates, Gj)
  end
#######################################################################
    for j in 1:N
        sj = sites[j]
        hj = (U+U2) * op("Nupdn",sj)
        Gj = exp(-im * tau / 2 * hj)
        push!(gates,Gj)
    end

append!(gates, reverse(gates))
#################################################
nsteps = Int(ttotal/tau) + 1
data = Matrix{Float64}(undef, nsteps, 2)

for (i, ti) in enumerate(0.0:tau:ttotal)
    data[i, 1] = ti
    
    ti≈ttotal && break
    psi0 = apply(gates, psi0; cutoff)

normalize!(psi0)
print("\nenergy = ", real(inner(psi0', H, psi0)))
Cuu = correlation_matrix(psi0, "Nup", "Ndn")
sum_diag = sum(diag(Cuu))
data[i, 2] = real(sum_diag)


end
writedlm("method1.txt", data)
end

Could you help me find my mistakes?
Thank you very much for your attention
Look forward to your reply.
Regards.
YD Shen.

2

Hi YD,
Generally the “TEBD” method of applying gates is known not to conserve energy in an precise sense. However, in the limit of the method being done to higher and higher accuracy it should approximately conserve energy, at least over short ranges of time.

I would suggest adjusting the various parameters that control the errors of this method. These are mainly:

  1. the time-step size (smaller is more accurate, though after too small of a time step then you are making more truncations per unit of time so it actually makes total accuracy worse to make the time step too small)
  2. the MPS truncation parameters like cutoff and `maxdim’

I would suggest you re-run your code with different choices of the above three parameters and see if the energy remains conserved better as you adjust them.

Also I would strongly suggest in your code that you include the interaction “Nupdn” terms into the same gates as the ones where you include the hopping term, instead of as separate gates. This may take some care e.g. to not over count the interactions and make sure they are included the correct number of times.