hi
I want to calculate the correlation function \left\langle\psi\right|e^{iHt}Ae^{-iHt}B\left|\psi\right\rangle by TEBD.
To derive the thermal states, I do the imaginary time evolution on enlarged Hilbert space \left | n \right \rangle \otimes \left | \tilde{n} \right \rangle, so I apply the gates to the identity operator. The sites of ancilla space are not evolved, so I add “prime” to the ancilla sites.
After imaginary time evolution, I need to calculate the correlation at moment t during real time evolution, but the function inner(MPO, MPO, MPO) does not work, so I apply the code "inner(MPO, apply(MPO,MPO)).
But the code does not work. What is the reason?
I just transferred from C++to Julia, so the code I wrote is still quite messy. Sorry.
using ITensors
using Printf
g = 0.0075
U = 0.1
t = 0.03
E = -0.01
x = -5
N =10
function ITensors.op(::OpName"expτSS", ::SiteType"Electron", s1::Index, s2::Index; tau)
if s1 == 1
factor = 1
else
factor1 = 0.5
end
s2 == N ? factor2 = 1 : factor2 = 0.5
h = factor1 * E * op("Ntot", s1) + factor1 * U * op("Nupdn", s1)
+ factor1 * E * op("Ntot", s2) + factor1 * U * op("Nupdn", s2)
+ t * op("Cdagup", s1) * op("Cup", s2) + t * op("Cdagdn", s1) * op("Cdn", s2)
return exp(tau * h)
end
function main(; N=10, cutoff=1E-8, tstep=0.1, beta_max=4.0)
println("begin")
onsite_E = E2
# Make an array of 'site' indices
s = siteinds("Electron", N; conserve_qns=false)
# Make gates (1,2),(2,3),(3,4),...
gates_imag = ops([("expτSS", (n, n + 1), (tau = -tstep / 2,)) for n in 1:(N - 1)], s)
# Include gates in reverse order to complete Trotter formula
append!(gates_imag, reverse(gates_imag))
gates_real = ops([("expτSS", (n, n + 1), (tau = -tstep / 2,)) for n in 1:(N - 1)], s)
append!(gates_real, reverse(gates_real))
# Hamiltonian of Hubbard Model
ampoHub = OpSum()
for j in 1:(N - 1)
j == 1 ? onsite_E = E1 : onsite_E = E2
ampoHub += E, "Ntot", j
ampoHub += U, "Nupdn", j
ampoHub += "Cdagup", j, "Cup", j + 1
ampoHub += "Cdagdn", j, "Cdn", j + 1
end
ampoHub += E, "Ntot", N
ampoHub += U, "Nupdn", N
Hub = MPO(ampoHub, s)
# Dx Hamiltonian of Hubbard Model
ampoDxHub = OpSum()
ampoDxHub += √2 * g, "Ntot", 1
DxHub = MPO(ampoDxHub, s)
# Initial state is infinite-temperature mixed state
#rho = MPO(s, "Id") ./ √2
X = MPO(s, "Id")
Y = MPO(prime(s), "Id")
x = convert(MPS, X)
y = convert(MPS, Y)
rho = outer(x, y) ./ √2
# Do the imaginary time evolution by applying the gates for Nsteps steps
for beta in 0 : tstep : beta_max
rho = apply(gates_imag, rho; cutoff)
rho = rho / tr(rho)
end
orho = apply(DxHub, rho; cutoff)
for i in 0 : 0.1 : 100
rho = apply(gates_real, rho; cutoff)
orho = apply(gates_real, orho; cutoff)
correlation = inner(rho, apply(DxHub, orho))
@printf("real = %.8f imag = %.8f", real(correlation), imag(correlation))
end
return nothing
end
main()