In MPSKit v0.13.2, ``timestep!`` yields unphysical results, whereas in v0.12.6, it does not

[ZongYongyue]

For a L=8 1d Hubbard model with 1/4 hole dopping, the correct result is Fig.1 but now with MPSKit v0.13.2 it's like Fig.2. The only variable between these two different results is the MPSKit version.

[lkdvos]

Is it possible that this is the result of the change in convention for sector keywords etc?

[ZongYongyue]

I don't know if this is related to what you mentioned, but I used both the method of adjusting filling with shift MPS and MPO quantum numbers (as mentioned in Mortier, Q. et al. Fermionic tensor network methods. arXiv (2024)), and the method of injecting quantum numbers into the MPS from one end, and compared DMRG with ED, obtaining consistent results. The ground state, ground state energy, and static correlations are all correct under both versions of MPSKit. Only during time evolution does v0.12 remain normal while v0.13 contains errors. Therefore I believe the issue must lie solely in the error within timestep!.

[ZongYongyue]

If the error is in the sector level, the DMRG results might also be wrong I thought?

[lkdvos]

Yeah, that seems like a fair assessment. Do you have some script so I can investigate a bit further? I don't really have a great idea of what might have changed so I'd have to check in more detail. As always, thanks for reporting in any case!

[ZongYongyue]

MPSKit@v0.12.6:

using TensorKit
using MPSKit
using MPSKitModels: FiniteChain
using DynamicalCorrelators: hubbard, e_min, chargedMPS

elt = Float64
I = FermionParity ⊠ SU2Irrep ⊠ U1Irrep 
L = 4
H = hubbard(elt, SU2Irrep, U1Irrep, FiniteChain(L); t=0.5, U=3, μ=0)
dop = -2
spin = 0
Vs = repeat([Vect[I]((Int(!isinteger(i)), i, j) => 1  for i in 0:1//2:1 for j in -(abs(dop) + 1):(abs(dop)+1)),], L - 1)
Ps = Vect[I]((0,0,-1)=>1, (1,1//2,0)=>1, (0,0,1)=>1)
st = FiniteMPS(rand, elt, [Ps for _ in 1:L], Vs; right=Vect[I]((Int(!isinteger(spin)), spin, dop) => 1))
gs, envs, delta = find_groundstate(st, H, DMRG2(trscheme = truncerr(1e-6)));
E0 = expectation_value(gs, H)
# -1.261122521051645

cm = e_min(elt, SU2Irrep, U1Irrep; side=:L)
gf = zeros(Float64, 4, 4) 
for i in 1:4, j in 1:4
    gf[i, j] = dot(chargedMPS(cm, gs, i), chargedMPS(cm, gs, j))
end
julia> gf
4×4 Matrix{Float64}:
 0.421264  0.466667  0.357925  0.30294
 0.466667  0.578736  0.44343   0.357925
 0.357925  0.44343   0.578736  0.466667
 0.30294   0.357925  0.466667  0.421264

function dpropagator(H::MPOHamiltonian, bras::Vector{<:FiniteMPS}, ket::FiniteMPS, dt::Number, ft::Number, n::Integer)
    times = collect(0:dt:ft)
    propagators = zeros(ComplexF64, length(bras), length(times))
    propagators[:,1] = [dot(bras[i], ket) for i in 1:length(bras)]
    envs = environments(ket, H)
    for (i, t) in enumerate(times[2:end])
        alg = t > n * dt ? TDVP() : TDVP2(;trscheme=truncerr(1e-4))
        ket, envs = timestep(ket, H, 0, dt, alg, envs)
        for j in eachindex(bras)
            propagators[j,i+1] = dot(bras[j], ket)
        end
    end
    return propagators
end
gs = complex(gs)
gfrt = dpropagator(H, [chargedMPS(cm, gs, i) for i in 1:4], chargedMPS(cm, gs, 1), 0.1, 5, 3)
julia> gfrt
4×51 Matrix{ComplexF64}:
 0.421264+0.0im  0.420052+0.0310577im  …  -0.242805-0.206401im  -0.222403-0.225351im
 0.466667+0.0im  0.465286+0.0354041im      -0.27928-0.246929im  -0.255387-0.268725im
 0.357925+0.0im  0.356593+0.0319219im     -0.286244-0.355789im  -0.258893-0.377668im
  0.30294+0.0im  0.301785+0.0274059im     -0.249478-0.324526im  -0.225426-0.343628im

MPSKit@v0.13.2:

using TensorKit
using MPSKit
using MPSKitModels: FiniteChain
using DynamicalCorrelators: hubbard, e_min, chargedMPS

elt = Float64
I = FermionParity ⊠ SU2Irrep ⊠ U1Irrep 
L = 4
H = hubbard(Float64, SU2Irrep, U1Irrep, FiniteChain(L); t=0.5, U=3, μ=0)
dop = -2
spin = 0
Vs = repeat([Vect[I]((Int(!isinteger(i)), i, j) => 1  for i in 0:1//2:1 for j in -(abs(dop) + 1):(abs(dop)+1)),], L - 1)
Ps = Vect[I]((0,0,-1)=>1, (1,1//2,0)=>1, (0,0,1)=>1)
st = FiniteMPS(rand, elt, [Ps for _ in 1:L], Vs; right=Vect[I]((Int(!isinteger(spin)), spin, dop) => 1))
gs, envs, delta = find_groundstate(st, H, DMRG2(trscheme = truncerr(1e-6)));
E0 = expectation_value(gs, H)
# -1.261122521051645
cm = e_min(elt, SU2Irrep, U1Irrep; side=:L)
gf = zeros(Float64, 4, 4) 
for i in 1:4, j in 1:4
    gf[i, j] = dot(chargedMPS(cm, gs, i), chargedMPS(cm, gs, j))
end
julia> gf
4×4 Matrix{Float64}:
 0.421264  0.466667  0.357925  0.30294
 0.466667  0.578736  0.44343   0.357925
 0.357925  0.44343   0.578736  0.466667
 0.30294   0.357925  0.466667  0.421264

function dpropagator(H::MPOHamiltonian, bras::Vector{<:FiniteMPS}, ket::FiniteMPS, dt::Number, ft::Number, n::Integer)
    times = collect(0:dt:ft)
    propagators = zeros(ComplexF64, length(bras), length(times))
    propagators[:,1] = [dot(bras[i], ket) for i in 1:length(bras)]
    envs = environments(ket, H)
    for (i, t) in enumerate(times[2:end])
        alg = t > n * dt ? TDVP() : TDVP2(;trscheme=truncerr(1e-4))
        ket, envs = timestep(ket, H, 0, dt, alg, envs)
        for j in eachindex(bras)
            propagators[j,i+1] = dot(bras[j], ket)
        end
    end
    return propagators
end
gs = complex(gs)
gfrt = dpropagator(H, [chargedMPS(cm, gs, i) for i in 1:4], chargedMPS(cm, gs, 1), 0.1, 5, 3)

julia> gfrt
4×51 Matrix{ComplexF64}:
 0.421264+0.0im  0.417849-0.0530129im  0.407666-0.105129im   …  0.301554+0.103058im  0.309208+0.0680936im
 0.466667+0.0im  0.463045-0.0577396im  0.452235-0.114566im      0.350343+0.127399im   0.36048+0.0865438im
 0.357925+0.0im  0.355826-0.0395584im  0.349546-0.0787068im     0.396718+0.226133im  0.422689+0.176047im
  0.30294+0.0im  0.301214-0.0330958im  0.296048-0.0658634im     0.351053+0.210519im  0.376089+0.165692im

benchmark with ED:

using QuantumLattices
using ExactDiagonalization

unitcell  = Lattice([0,0]; vectors = [[1, 0], ])
lattice  = Lattice(unitcell, (4,), ('o',))

hilbert = Hilbert(site=>Fock{:f}(1, 2) for site=1:length(lattice))
quantumnumber = ℕ(2)
t = Hopping(:t, -0.5, 1)
U = Hubbard(:U, 3.0)
ed = ED(lattice, hilbert, (t, tzu, tzd, U), quantumnumber)
eigensystem = eigen(ed; nev=1)
gs = eigensystem.vectors[1];
E0 = eigensystem.values
print(E0)
#[-1.2611225210516321]

gf = zeros(Float64, 4, 4)
bs1 = BinaryBases(8, ℕ(2))
bs2 = BinaryBases(8, ℕ(1))
for i in 1:4, j in 1:4
    opt = Operator(1, 𝕔(i, 1, -1//2, 1)) + Operator(1, 𝕔(i, 1, 1//2, 1)) 
    opm = matrix(opt,(bs2, bs1),Table(hilbert, Metric(EDKind(hilbert), hilbert)))
    psi1 = opm*gs

    opt = Operator(1, 𝕔(j, 1, -1//2, 1)) + Operator(1, 𝕔(j, 1, 1//2, 1)) 
    opm = matrix(opt,(bs2, bs1),Table(hilbert, Metric(EDKind(hilbert), hilbert)))
    psi2 = opm*gs
    gf[i, j] = dot(psi1, psi2)
end
julia> gf
4×4 Matrix{Float64}:
 0.421264  0.466667  0.357925  0.30294
 0.466667  0.578736  0.44343   0.357925
 0.357925  0.44343   0.578736  0.466667
 0.30294   0.357925  0.466667  0.421264

[ZongYongyue]

You can see that for both the ground state energy and static correlations, the results from both versions match the exact diagonalization results. However, as time evolves, the dynamical correlation function gfrt calculated by MPSKit v0.13 shows increasingly significant deviations compared to the results from v0.12.

[ZongYongyue]

gfrt.zip

[ZongYongyue]

This is even clearer when examining the time evolution from the ground state. Theoretically, the inner product should consistently remain 1, as seen in the results from v0.12:

julia> gfrt = dpropagator(H, [gs,], gs, 0.1, 5, 3)
1×51 Matrix{ComplexF64}:
 1.0+0.0im  0.992058+0.125778im  0.96836+0.249559im  …  0.99463-0.103499im  0.999749+0.0224254im

However, in v0.13, the inner product results clearly violate physical principles:

julia> gfrt = dpropagator(H, [gs,], gs, 0.1, 5, 3)
1×51 Matrix{ComplexF64}:
 1.0+0.0im  0.997272-0.0738205im  0.989101-0.147238im  0.975533-0.219852im  …  -0.887499+0.460809im  -0.85106+0.525068im

[ZongYongyue]

This is even clearer when examining the time evolution from the ground state. Theoretically, the inner product should consistently remain 1, as seen in the results from v0.12:

julia> gfrt = dpropagator(H, [gs,], gs, 0.1, 5, 3)
1×51 Matrix{ComplexF64}:
1.0+0.0im 0.992058+0.125778im 0.96836+0.249559im … 0.99463-0.103499im 0.999749+0.0224254im
However, in v0.13, the inner product results clearly violate physical principles:

julia> gfrt = dpropagator(H, [gs,], gs, 0.1, 5, 3)
1×51 Matrix{ComplexF64}:
1.0+0.0im 0.997272-0.0738205im 0.989101-0.147238im 0.975533-0.219852im … -0.887499+0.460809im -0.85106+0.525068im

To clarify, the modulus of the inner product between the ground state and the time-evolved ground state should be 1, but there should exist a phase factor of exp(iE_0 t). Both versions of MPSKit can produce inner products with modulus 1, but v0.13 gives completely wrong phases, while v0.12 provides precisely correct phases.

[exp(E0*im*t) for t in 0:0.1:5]
51-element Vector{ComplexF64}:
                1.0 - 0.0im
 0.9920583837970895 - 0.12577822998399443im
 0.9683596737241869 - 0.2495586951095603im
                    ⋮
 0.9946295290436121 + 0.10349927513988884im
 0.9997485186915878 - 0.02242541803347092im

[lkdvos]

Thanks for the clarifications! I am currently out of office but will be back next week to hopefully fix this asap

[lkdvos]

Apologies for taking a while to get to this. I looked it over and I think I discovered what was going wrong, I linked a PR with a fix, but won't be able to merge and properly test this until next week. It would be great if you could try it out and check if the new code actually resolves the issues you are running into.

[ZongYongyue]

I apologize for disturbing you multiple times outside of work hours, and thank you very much for fixing the bug. The results are now correct, but there’s a minor issue: timestep requires both the MPS and MPO to have complex eltype, otherwise it throws an error.

[sumansuhag]

Hi!!!!!
Thanks for bringing this up, ZongYongyue! 🚨 It looks like there's a problem with timestep! in MPSKit v0.13.2 that's causing issues with the 1D Hubbard model, especially since it worked fine in v0.12.6. This seems like a regression bug.

It shows just how important it is to have solid testing and version control for scientific computing tools to keep everything accurate and trustworthy. Unphysical results can really mess with research, so bug reports like yours are crucial for keeping these tools reliable.

By the way, have you noticed any strange behaviors in other models or settings with v0.13.2, aside from the L=8 1D Hubbard model? Your detailed findings are super helpful for the MPSKit team. #MPSKit #ComputationalPhysics #QuantumManyBodyPhysics #HubbardModel #RegressionBug #ScientificComputing #Reproducibility #SoftwareQuality