Dear ITensor developers and users
I am new to ITensors.jl, so my question might be quite naive (and might have been already solved in this discourse). I am studying how to time-evolve electron systems. For this, I think about the following simple electron systems:
First, I did DMRG for this model in order to get the ground state, and then, time-evolved this ground state by the same Hamiltonian.
##make Hamiltonian and do dmrg to get ground state.
sites = siteinds("Electron", N; conserve_qns = true)
ampo = OpSum()
for j in 1:N-1
ampo += -t, "Cdagup", j, "Cup", j + 1
ampo += -t, "Cdagdn", j, "Cdn", j + 1
ampo += -t, "Cdagup", j + 1, "Cup", j
ampo += -t, "Cdagdn", j + 1, "Cdn", j
end
H = MPO(ampo, sites)
state = [isodd(n) ? "Up" : "Dn" for n=1:N]
psi0 = productMPS(sites, state)
energy, psi = dmrg(H, psi0, sweeps)
## time-evolution part
gates = ITensor[]
for j = 1:N-1
s1 = sites[j]
s2 = sites[j+1]
hj_hop =
-t * op("Adagup * F", s1) * op("Aup", s2)
-t * op("Adagdn", s1) * op("F * Adn", s2)
+t * op("Aup * F", s1) * op("Adagup", s2)
+t * op("Adn", s1) * op("F * Adagdn", s2)
Gj_hop = exp(-im * tau / 2 * hj_hop)
push!(gates, Gj_hop)
end
# Include gates in reverse order too
# (N,N-1),(N-1,N-2),...
append!(gates, reverse(gates))
for t in 0.0:tau:ttotal
Sz = inner(psi_t0,psi) * inner(psi,psi_t0)
println("$t $Sz")
t≈ttotal && break
psi = apply(gates, psi; cutoff)
normalize!(psi)
end
So, I was supposed to get the |\braket{\psi(t)| \psi(0)}|^2 \simeq 1 for the any t\geq 0. However, the results are like following:
#result for ITensor.jl
0.00 1.000000
0.02 0.999616
0.04 0.998469
0.06 0.996563
0.08 0.993902
0.10 0.990492
0.12 0.986341
0.14 0.981460
0.16 0.975859
0.18 0.969553
0.20 0.962555
0.22 0.954883
0.24 0.946554
0.26 0.937587
0.28 0.928003
0.30 0.917825
0.32 0.907075
0.34 0.895776
0.36 0.883955
0.38 0.871637
0.40 0.858848
0.42 0.845616
0.44 0.831969
0.46 0.817936
0.48 0.803544
0.50 0.788823
My question is “Why did the absolute value of overlap |\braket{\psi(t)|\psi(0)}|^2 decays ?“ I also tried this by ITensor C++ version , but in this case, decay did not happen:
auto gates = vector<BondGate>();
// Create the gates exp(-i*tstep/2*hterm)
// and add them to gates
for (int b = 1; b <= N - 1; ++b) {
auto hterm = -1. * sites.op("Adagup*F", b) * sites.op("Aup", b + 1);
hterm += -1. * sites.op("Adagdn", b) * sites.op("F*Adn", b + 1);
hterm += +1. * sites.op("Aup*F", b) * sites.op("Adagup", b + 1);
hterm += +1. * sites.op("Adn", b) * sites.op("F*Adagdn", b + 1);
auto g = BondGate(sites, b, b + 1, BondGate::tReal, tstep / 2., hterm);
gates.push_back(g);
}
// Create the gates exp(-i*tstep/2*hterm) in reverse
// order (to get a second order Trotter breakup which
// does a time step of "tstep") and add them to gates
for (int b = N - 1; b >= 1; --b) {
auto hterm = -1. * sites.op("Adagup*F", b) * sites.op("Aup", b + 1);
hterm += -1. * sites.op("Adagdn", b) * sites.op("F*Adn", b + 1);
hterm += +1. * sites.op("Aup*F", b) * sites.op("Adagup", b + 1);
hterm += +1. * sites.op("Adn", b) * sites.op("F*Adagdn", b + 1);
auto g = BondGate(sites, b, b + 1, BondGate::tReal, tstep / 2., hterm);
gates.push_back(g);
}
#result for C++ ver.
0.02 0.999998
0.04 0.999997
0.06 0.999996
0.08 0.999995
0.10 0.999994
0.12 0.999993
0.14 0.999993
0.16 0.999992
0.18 0.999991
0.20 0.999990
0.22 0.999990
0.24 0.999989
0.26 0.999988
0.28 0.999987
0.30 0.999986
0.32 0.999986
0.34 0.999985
0.36 0.999984
0.38 0.999983
0.40 0.999982
0.42 0.999981
0.44 0.999981
0.46 0.999980
0.48 0.999979
0.50 0.999978
Thank you very much in advance.