Time dependent TDVP

.gitignore

@@ -1,3 +1,4 @@
 docs/build
 Manifest.toml
 .vscode
+.DS_Store

examples/timeevolution.jl

@@ -0,0 +1,73 @@
+using MPSKit, MPSKitModels, TensorKit, Test
+
+#local redefinition
+function nonsym_ising(; J=1, g=0.5)
+    (sx, _, sz) = nonsym_spintensors(1//2)
+
+    return MPOHamiltonian(
+        LocalOperator(-4 * J * sz ⊗ sz, (1, 2)) + LocalOperator(-2 * g * sx, (1,))
+    )
+end
+
+# for now this uses the old MPSKitModels code
+
+(sx, _, sz) = nonsym_spintensors(1//2)
+σₓ = 2 * sx;
+σz = 2 * sz;  #factor 2 to get pauli matrices
+J = 1;
+g = 0.5;
+
+HJ = MPOHamiltonian(LocalOperator(-J * σz ⊗ σz, (1, 2)))
+Hg = MPOHamiltonian(LocalOperator(-g * σₓ, (1,)))
+H₀ = HJ + Hg
+
+gs = InfiniteMPS([2], [50]); #create MPS with physical bond dimension d=2 and virtual D=50
+(gs, envs, _) = find_groundstate(gs, H₀, VUMPS(; maxiter=400));
+
+sx_gs = sum(expectation_value(gs, σₓ)) / length(gs)
+@show sx_gs
+E_gs = sum(expectation_value(gs, H₀, envs)) / length(gs)
+@show E_gs
+
+# some time function that slowly ramps
+f = t -> t == 0.0 ? 1 : 1 + 2 * min(0.1 * t, 2.0)
+
+# time dependent Hamiltonian = HJ+f(t)*Hg
+Hₜ = TimeDepProblem((HJ, Hg), (t -> 1, f));
+
+# the corresponding environment
+envs = environments(gs, Hₜ);
+
+# the algorithm we will use for time evolution
+# expalg performs the integration of the tdvp equations
+# here we choose the implicit midpoint method IM() (others are availible/ easily implemented see integrators.jl)
+alg = TDVP(; expalg=IM())
+
+# timestep and how many steps
+dt = 0.001
+N = 100
+
+# containers for observables
+sxs = zeros(ComplexF64, N + 1)
+Es = zeros(ComplexF64, N + 1)
+
+# the actual time evolution
+let nstate = copy(gs), nenvs = environments(gs, Hₜ), t = 0.0
+    for i in 1:(N + 1)
+        sxs[i] = expectation_value(nstate, σₓ, 1)
+        Es[i] = expectation_value(nstate, Hₜ, t, nenvs)[1]
+        nstate, nenvs = timestep(nstate, Hₜ, t, dt, alg, nenvs)
+        t += dt
+    end
+end
+
+#check that for t=0 before timestep everything is same as in gs
+@test real(sxs[1]) ≈ real(sx_gs)
+@test real(Es[1]) ≈ real(E_gs)
+
+# numerical free fermion solution @ t=0.1
+sxt = 0.25869106309840484
+Et = -1.0661310695057724
+
+@show abs(real(sxs[end]) / sxt - 1)
+@show abs(real(Es[end]) / Et - 1)

src/MPSKit.jl

@@ -27,7 +27,8 @@ export entanglementplot, transferplot
 
 # hamiltonian things
 export Cache
-export SparseMPO, MPOHamiltonian, DenseMPO, MPOMultiline
+export SparseMPO,
+    MPOHamiltonian, DenseMPO, MPOMultiline, UntimedOperator, TimedOperator, LazySum, ConvertOperator
 export ∂C, ∂AC, ∂AC2, environments, expectation_value, effective_excitation_hamiltonian
 export leftenv, rightenv
 
@@ -37,6 +38,7 @@ export VUMPS, DMRG, DMRG2, IDMRG1, IDMRG2, GradientGrassmann
 export excitations, FiniteExcited, QuasiparticleAnsatz
 export marek_gap, correlation_length, correlator
 export timestep!, timestep, TDVP, TDVP2, make_time_mpo, WI, WII, TaylorCluster
+export time_evolve, time_evolve!
 export splitham, infinite_temperature, entanglement_spectrum, transfer_spectrum, variance
 export changebonds!, changebonds, VUMPSSvdCut, OptimalExpand, SvdCut, UnionTrunc, RandExpand
 export entropy
@@ -79,6 +81,8 @@ include("operators/sparsempo/sparsempo.jl")
 include("operators/mpohamiltonian.jl") # the mpohamiltonian objects
 include("operators/mpomultiline.jl")
 include("operators/projection.jl")
+include("operators/multipliedoperator.jl")
+include("operators/lazysum.jl")
 
 include("transfermatrix/transfermatrix.jl")
 include("transfermatrix/transfer.jl")
@@ -89,6 +93,7 @@ include("environments/FinEnv.jl")
 include("environments/abstractinfenv.jl")
 include("environments/permpoinfenv.jl")
 include("environments/mpohaminfenv.jl")
+include("environments/multipleenv.jl")
 include("environments/qpenv.jl")
 include("environments/idmrgenv.jl")
 include("environments/lazylincocache.jl")
@@ -108,7 +113,10 @@ include("algorithms/changebonds/svdcut.jl")
 include("algorithms/changebonds/randexpand.jl")
 
 include("algorithms/timestep/tdvp.jl")
+include("algorithms/timestep/windowtdvp.jl")
 include("algorithms/timestep/timeevmpo.jl")
+include("algorithms/timestep/integrators.jl")
+include("algorithms/timestep/time_evolve.jl")
 
 include("algorithms/groundstate/vumps.jl")
 include("algorithms/groundstate/idmrg.jl")

src/algorithms/derivatives.jl

@@ -1,5 +1,6 @@
 # Given a state and it's environments, we can act on it
 
+#make this into one struct?
 """
     Draft operators
 """
@@ -21,12 +22,16 @@ struct MPO_∂∂AC2{O,L,R}
     rightenv::R
 end
 
-Base.:*(h::Union{MPO_∂∂C,MPO_∂∂AC,MPO_∂∂AC2}, v) = h(v);
+const DerivativeOperator = Union{MPO_∂∂C,MPO_∂∂AC,MPO_∂∂AC2}
+
+Base.:*(h::DerivativeOperator, v) = h(v);
 
 (h::MPO_∂∂C)(x) = ∂C(x, h.leftenv, h.rightenv);
 (h::MPO_∂∂AC)(x) = ∂AC(x, h.o, h.leftenv, h.rightenv);
 (h::MPO_∂∂AC2)(x) = ∂AC2(x, h.o1, h.o2, h.leftenv, h.rightenv);
+(h::DerivativeOperator)(v, ::Number) = h(v)
 
+# can we reduce code duplication for ∂∂C,∂∂AC,∂∂AC2 by defining constructors for DerivativeOperator?
 # draft operator constructors
 function ∂∂C(pos::Int, mps, opp::Union{MPOHamiltonian,SparseMPO,DenseMPO}, cache)
     return MPO_∂∂C(leftenv(cache, pos + 1, mps), rightenv(cache, pos, mps))
@@ -326,3 +331,36 @@ function ∂∂AC2(pos::Int, state, opp::ProjectionOperator, env)
         rightenv(env, pos + 1, state),
     )
 end
+
+# MultipliedOperator and SumOfOperators
+(x::TimedOperator{<:DerivativeOperator})(y, t::Number) = x.f(t) * x.op(y)
+(x::UntimedOperator{<:DerivativeOperator})(y, ::Number) = x.f * x.op(y)
+(x::LazySum{<:Union{MPSKit.MultipliedOperator{D}, D} where {D<:MPSKit.DerivativeOperator}})(y, t::Number) = sum(O -> O(y, t), x) #not strict enough?
+
+function ∂∂C(pos::Int, mps, opp::MultipliedOperator, cache)
+    return MultipliedOperator(∂∂C(pos::Int, mps, opp.op, cache), opp.f)
+end
+
+function ∂∂AC(pos::Int, mps, opp::MultipliedOperator, cache)
+    return MultipliedOperator(∂∂AC(pos::Int, mps, opp.op, cache), opp.f)
+end
+
+function ∂∂AC2(pos::Int, mps, opp::MultipliedOperator, cache)
+    return MultipliedOperator(∂∂AC2(pos::Int, mps, opp.op, cache), opp.f)
+end
+
+function ∂∂C(pos::Int, mps, opp::LazySum, cache::MultipleEnvironments)
+    return LazySum{Union{MPO_∂∂C,MultipliedOperator{<:MPO_∂∂C}}}(map((op, openv) -> ∂∂C(pos, mps, op, openv), opp.ops, cache.envs))
+end
+
+function ∂∂AC(pos::Int, mps, opp::LazySum, cache::MultipleEnvironments)
+    return LazySum{Union{MPO_∂∂AC,MultipliedOperator{<:MPO_∂∂AC}}}(
+        map((op, openv) -> ∂∂AC(pos, mps, op, openv), opp.ops, cache.envs)
+    )
+end
+
+function ∂∂AC2(pos::Int, mps, opp::LazySum, cache::MultipleEnvironments)
+    return LazySum{Union{MPO_∂∂AC2,MultipliedOperator{<:MPO_∂∂AC2}}}(
+        map((op, openv) -> ∂∂AC2(pos, mps, op, openv), opp.ops, cache.envs)
+    )
+end