Time-evolution of electron system and the difference Julia and C++

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:

\hat{H} = \sum_{i, \sigma} c^{\dagger}_{i,\sigma} c_{i+1,\sigma} + c^{\dagger}_{i,\sigma} c_{i-1,\sigma}

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.

Hi Riku,
It’s a good question. There could be many reasons, though they could be subtle. Hopefully there isn’t an issue with the Julia version but if there is we would be happy to learn about it.

How well is your ground state converged at the beginning in the Julia code? Is the energy as low as in the C++ code?

Also do you use the same cutoff during the time evolution part of both codes?

Miles

Dear Miles,

Thank you for your prompt and kind reply.

Both cases have the same Sweep parameters in DMRG :

#Julia
sweeps = Sweeps(16)
maxdim!(sweeps, 50, 80, 100, 120, 160, 200, 240, 320, 400, 480, 560, 640)
cutoff!(sweeps, 1E-8, 1E-9, 1E-9, 1E-9, 1E-10, 1E-10, 1E-10, 1E-10, 1E-11,
1E-11, 1E-12, 1E-12, 1E-13)
noise!(sweeps, 1E-6, 1E-7, 1E-8, 1E-9, 1E-10, 1E-10, 1E-11, 1E-12, 1E-13)
//C++
auto sweeps = Sweeps(16);
sweeps.maxdim() = 50, 80, 100, 120, 160, 200, 240, 320, 400, 480, 560, 640;
sweeps.noise() = 1E-6, 1E-7, 1E-8, 1E-9, 1E-10, 1E-10, 1E-11, 1E-12, 1E-13;
sweeps.cutoff() = 1E-8, 1E-9, 1E-9, 1E-9, 1E-10, 1E-10, 1E-10, 1E-10, 1E-11,
    1E-11, 1E-12, 1E-12, 1E-13;

Their convergences of the ground state energy in DMRG are almost same in both codes:

# Julia
After sweep 14 energy=-126.6023540652746  maxlinkdim=640 maxerr=1.04E-07 time=31.505
After sweep 15 energy=-126.60235406475219  maxlinkdim=640 maxerr=1.04E-07 time=31.412
After sweep 16 energy=-126.60235406472395  maxlinkdim=640 maxerr=1.04E-07 time=31.282
# C++
    vN Entropy at center bond b=50 = 2.096864691904
    Eigs at center bond b=50: 0.4266 0.0974 0.0974 0.0974 0.0974 0.0222 0.0222 0.0222 0.0222 0.0222 
    Largest link dim during sweep 14/16 was 640
    Largest truncation error: 4.20922e-08
    Energy after sweep 14/16 is -126.602354065272
    Sweep 14/16 CPU time = 33.28s (Wall time = 35.98s)

    vN Entropy at center bond b=50 = 2.096867084767
    Eigs at center bond b=50: 0.4266 0.0974 0.0974 0.0974 0.0974 0.0222 0.0222 0.0222 0.0222 0.0222 
    Largest link dim during sweep 15/16 was 640
    Largest truncation error: 4.21051e-08
    Energy after sweep 15/16 is -126.602354064750
    Sweep 15/16 CPU time = 33.11s (Wall time = 35.76s)

    vN Entropy at center bond b=50 = 2.096867295358
    Eigs at center bond b=50: 0.4266 0.0974 0.0974 0.0974 0.0974 0.0222 0.0222 0.0222 0.0222 0.0222 
    Largest link dim during sweep 16/16 was 640
    Largest truncation error: 4.21055e-08
    Energy after sweep 16/16 is -126.602354064722
    Sweep 16/16 CPU time = 33.06s (Wall time = 35.55s)

In addition, the cutoff during the time-evolution and the time step are also same in the both codes. (cutoff=1E-8 and \tau = 0.02)

I got the results above in the case where the system size N was 100. When I tried the same for the small system (N=12), the truncation error in Julia was \sim 10^{-14}, but even in this case, the decay happened. The decay of the overlap |\braket{\psi(t)|\psi(0)}|^2 seems to have the system-size dependence.

# N = 12

# the convergence in DMRG
After sweep 14 energy=-14.59245962111034  maxlinkdim=500 maxerr=9.97E-14 time=0.885
After sweep 15 energy=-14.592459621112537  maxlinkdim=501 maxerr=9.96E-14 time=0.848
After sweep 16 energy=-14.592459621112821  maxlinkdim=504 maxerr=9.97E-14 time=0.773

# the overlap|<\psi(t)|\psi(0)>| in time evolution 
  t		overlap
0.44 0.9018382448143883
0.46 0.8933082138497178
0.48 0.884501345720261
0.5 0.875429446645427 
1 Like