Functionality to time evolve WindowMPS states

examples/windowmps.jl

@@ -1,37 +1,169 @@
+using Pkg #to be removed
+Pkg.activate("/Users/daan/Desktop/TimedtdvpTest/TimedTDVP") #to be removed
 using MPSKit, MPSKitModels, TensorKit, Plots
 
-let
-    #defining the hamiltonian
-    th = nonsym_ising_ham(; lambda=0.3)
-    sx, sy, sz = nonsym_spintensors(1 // 2)
+function my_transverse_field_ising(gs)
+    L = length(gs)
+    lattice = InfiniteChain(L)
+    ZZ = rmul!(σᶻᶻ(), -1)
+    X = rmul!(σˣ(), -1)
+    a = @mpoham sum(ZZ{i,j} for (i, j) in nearest_neighbours(lattice))
+    b = @mpoham sum(gs[i] * X{i} for i in vertices(lattice))
+    return a + b
+end
+
+function my_timedependent_ising(gl, gs, gr, f)
+    L = length(gs)
+    lattice = InfiniteChain(1)
+    latticeL = InfiniteChain(L)
+    ZZ = rmul!(σᶻᶻ(), -1)
+    X = rmul!(σˣ(), -1)
+
+    ZZl = @mpoham sum(ZZ{i,j} for (i, j) in nearest_neighbours(lattice))
+    ZZm = @mpoham sum(ZZ{i,j} for (i, j) in nearest_neighbours(latticeL))
+    ZZr = @mpoham sum(ZZ{i,j} for (i, j) in nearest_neighbours(lattice))
+
+    Xl = @mpoham sum(gl * X{i} for i in vertices(lattice))
+    Xm = @mpoham sum(gs[i] * X{i} for i in vertices(latticeL))
+    Xr = @mpoham sum(gr * X{i} for i in vertices(lattice))
 
-    #initilizing a random mps
-    ts = InfiniteMPS([ℂ^2], [ℂ^12])
+    H1 = Window(ZZl, ZZm, ZZr)
+    H2 = Window(Xl, Xm, Xr)
+    return LazySum([H1, MultipliedOperator(H2, f)])
+end
 
-    #Finding the groundstate
-    (ts, envs, _) = find_groundstate(ts, th, VUMPS(; maxiter=400))
+function my_expectation_value(Ψwindow::WindowMPS,
+                              O::Window{A,A,A}) where {A<:TrivialTensorMap{ComplexSpace,1,1,
+                                                                           Matrix{ComplexF64}}}
+    left = expectation_value(Ψwindow.left, O.left)
+    middle = expectation_value(Ψwindow, O.middle)
+    right = expectation_value(Ψwindow.right, O.right)
+    return vcat(left, middle, right)
+end
 
-    len = 20
-    deltat = 0.05
-    totaltime = 3.0
-    middle = Int(round(len / 2))
+function my_finalize(t, Ψ, H, envs, si, tosave)
+    push!(tosave, my_expectation_value(Ψ, si))
+    return Ψ, envs
+end
 
-    #apply a single spinflip at the middle site
-    mpco = WindowMPS(ts, len)
-    @tensor mpco.AC[middle][-1 -2; -3] := mpco.AC[middle][-1, 1, -3] * sx[-2, 1]
-    normalize!(mpco)
+# WindowMPS as bath
+#-------------------
 
-    envs = environments(mpco, th)
+#define the hamiltonian
+H = transverse_field_ising(; g=0.3)
+sx, sy, sz = σˣ(), σʸ(), σᶻ()
 
-    szdat = [expectation_value(mpco, sz)]
+#initilizing a random mps
+Ψ = InfiniteMPS([ℂ^2], [ℂ^12])
 
-    for i in 1:(totaltime / deltat)
-        (mpco, envs) = timestep(mpco, th, deltat, TDVP2(; trscheme=truncdim(20)), envs)
-        push!(szdat, expectation_value(mpco, sz))
-    end
+#Finding the groundstate
+(Ψ, envs, _) = find_groundstate(Ψ, H, VUMPS(; maxiter=400))
 
-    display(heatmap(real.(reduce((a, b) -> [a b], szdat))))
+len = 20
+middle = round(Int, len / 2)
 
-    println("Enter to continue ...")
-    readline()
-end
+# make a WindowMPS by promoting len sites to the window part
+# by setting fixleft=fixright=true we indicate that the infinite parts will not change
+Ψwindow = WindowMPS(Ψ, len; fixleft=true, fixright=true)
+#apply a single spinflip at the middle site
+@tensor Ψwindow.AC[middle][-1 -2; -3] := Ψwindow.AC[middle][-1, 1, -3] * sx[-2, 1];
+normalize!(Ψwindow);
+
+# create the environment
+# note: this method is only defined for fixleft=true=fixright=true WindowMPS and assumes the same H for left and right infinite environments
+# if it errors for these reasons use H  = Window(Hleft,Hmiddle,Hright) instead
+envs = environments(Ψwindow, H);
+szdat = [expectation_value(Ψwindow, sz)]
+
+#setup for time_evolve
+alg = TDVP2(; trscheme=truncdim(20),
+            finalize=(t, Ψ, H, envs) -> my_finalize(t, Ψ, H, envs, sz, szdat));
+t_span = 0:0.05:3.0
+Ψwindow, envs = time_evolve!(Ψwindow, H, t_span, alg, envs);
+
+display(heatmap(real.(reduce((a, b) -> [a b], szdat))))
+
+# WindowMPS as interpolation
+#----------------------------
+
+# The Hamiltonian wil be -(∑_{<i,j>} Z_i Z_j + f(t) * ∑_{<i>} g_i X_i)
+gl = 3.0
+gr = 4.0
+L = 10
+gs = range(gl, gr; length=L); #interpolating values for g
+
+Hl = my_transverse_field_ising([gl]);
+Hr = my_transverse_field_ising([gr]);
+
+Hfin = my_transverse_field_ising(gs);
+# Note: it is important that the lattice for the finite part of the WindowMPS is infinite.
+#       For a finite lattice @mpoham is smart and does not construct the terms in the MPOHamiltonian
+#       that will not be used due to the boundary. For a WindowMPS there is no boundary and we thus
+#       we need the MPO for the infinite lattice.
+
+Hwindow = Window(Hl, Hfin, Hr);
+sx, sy, sz = σˣ(), σʸ(), σᶻ()
+
+#initilizing a random mps
+D = 12
+Ψl = InfiniteMPS([ℂ^2], [ℂ^D])
+Ψr = InfiniteMPS([ℂ^2], [ℂ^D]) #we do not want Ψr === ψl
+#Ts = map(i->TensorMap(rand,ComplexF64,ℂ^D*ℂ^2,ℂ^D),eachindex(gs))
+Ts = fill(TensorMap(rand, ComplexF64, ℂ^D * ℂ^2, ℂ^D), L);
+Ψwindow = WindowMPS(Ψl, Ts, Ψr);
+
+# finding the groundstate
+(Ψwindow, envs, _) = find_groundstate(Ψwindow, Hwindow);
+# the alg for find_groundstate has to be a Window(alg_left,alg_middle,alg_right).
+# If no alg is specified the default is Window(VUMP(),DMRG(),VUMPS())
+es = real.(expectation_value(Ψwindow, Hwindow, envs));
+# note that es[1] is the expectation_value of Ψwindow.left and es[end] that of Ψwindow.right
+scatter(0:(L + 1), es; label="")
+
+# WindowMPS for non-uniform quench
+#---------------------------------
+
+#define the hamiltonian
+g_uni = 3.0
+gl = 3.0
+gr = 4.0
+L = 40
+
+Hl = my_transverse_field_ising([g_uni]);
+Hr = my_transverse_field_ising([g_uni]);
+Hfin = my_transverse_field_ising([g_uni]);
+Hgs = Window(Hl, Hfin, Hr);
+
+D = 12
+Ψl = InfiniteMPS([ℂ^2], [ℂ^D])
+Ψr = InfiniteMPS([ℂ^2], [ℂ^D]) #we do not want Ψr === ψl
+#Ts = map(i->TensorMap(rand,ComplexF64,ℂ^D*ℂ^2,ℂ^D),eachindex(gs))
+Ts = fill(TensorMap(rand, ComplexF64, ℂ^D * ℂ^2, ℂ^D), L);
+Ψwindow = WindowMPS(Ψl, Ts, Ψr);
+
+# finding the groundstate
+(Ψwindow, envs_gs, _) = find_groundstate(Ψwindow, Hgs);
+
+#define the quench Hamiltonian
+f(t) = 0.1 * t #we take a linear ramp
+gs = range(gl, gr; length=L); #interpolating values for g
+Hquench = my_timedependent_ising(gl, gs, gr, f);
+# Hquench is a time-dependent Hamiltonian i.e. we can do H(t) to get the instantanious Hamiltonian.
+# Note: To get an expectation_value of a time-dependent Hamiltonian one needs to give H(t) tot he function.
+
+envs = environments(Ψwindow, Hquench);
+
+sdat = [my_expectation_value(Ψwindow, Window(sx, sx, sx))]
+
+#setup for time_evolve
+left_alg = rightalg = TDVP();
+middle_alg = TDVP2(; trscheme=truncdim(20));
+alg = WindowTDVP(; left=left_alg, middle=middle_alg, right=rightalg,
+                 finalize=(t, Ψ, H, envs) -> my_finalize(t, Ψ, H, envs, Window(sx, sx, sx),
+                                                         sdat));
+t_span = 0:0.005:1.0
+
+Ψwindow, envs = time_evolve!(Ψwindow, Hquench, t_span, alg, envs; verbose=true);
+
+display(heatmap(t_span, 0:(L + 1), real.(reduce((a, b) -> [a b], sdat)); xlabel="t",
+                ylabel="i"))

src/MPSKit.jl

@@ -18,7 +18,7 @@ export PeriodicArray, WindowArray
 export MPSTensor
 export QP, LeftGaugedQP, RightGaugedQP
 export leftorth,
-       rightorth, leftorth!, rightorth!, poison!, uniform_leftorth, uniform_rightorth
+       rightorth, leftorth!, rightorth!, invalidate!, uniform_leftorth, uniform_rightorth
 export r_LL, l_LL, r_RR, l_RR, r_RL, r_LR, l_RL, l_LR # should be properties
 
 # useful utility functions?
@@ -39,8 +39,8 @@ export find_groundstate!, find_groundstate, leading_boundary
 export VUMPS, VOMPS, DMRG, DMRG2, IDMRG1, IDMRG2, GradientGrassmann
 export excitations, FiniteExcited, QuasiparticleAnsatz, ChepigaAnsatz, ChepigaAnsatz2
 export marek_gap, correlation_length, correlator
-export time_evolve, timestep!, timestep
-export TDVP, TDVP2, make_time_mpo, WI, WII, TaylorCluster
+export time_evolve, time_evolve!, timestep!, timestep
+export TDVP, TDVP2, WindowTDVP, make_time_mpo, WI, WII, TaylorCluster
 export splitham, infinite_temperature, entanglement_spectrum, transfer_spectrum, variance
 export changebonds!, changebonds, VUMPSSvdCut, OptimalExpand, SvdCut, UnionTrunc, RandExpand
 export entropy
@@ -78,6 +78,7 @@ include("utility/multiline.jl")
 include("utility/utility.jl") # random utility functions
 include("utility/plotting.jl")
 include("utility/linearcombination.jl")
+include("utility/timedependence.jl")
 
 # maybe we should introduce an abstract state type
 include("states/abstractmps.jl")
@@ -95,17 +96,17 @@ include("operators/sparsempo/sparsempo.jl")
 include("operators/mpohamiltonian.jl") # the mpohamiltonian objects
 include("operators/mpomultiline.jl")
 include("operators/projection.jl")
-include("operators/timedependence.jl")
 include("operators/multipliedoperator.jl")
 include("operators/lazysum.jl")
 
 include("transfermatrix/transfermatrix.jl")
 include("transfermatrix/transfer.jl")
 
-include("environments/FinEnv.jl")
+include("environments/finenv.jl")
 include("environments/abstractinfenv.jl")
 include("environments/permpoinfenv.jl")
 include("environments/mpohaminfenv.jl")
+include("environments/windowenv.jl")
 include("environments/qpenv.jl")
 include("environments/multipleenv.jl")
 include("environments/idmrgenv.jl")
@@ -128,6 +129,8 @@ include("algorithms/timestep/tdvp.jl")
 include("algorithms/timestep/timeevmpo.jl")
 include("algorithms/timestep/integrators.jl")
 include("algorithms/timestep/time_evolve.jl")
+include("algorithms/timestep/windowtdvp.jl")
+include("algorithms/timestep/timestep.jl")
 
 include("algorithms/groundstate/vumps.jl")
 include("algorithms/groundstate/idmrg.jl")

src/algorithms/expval.jl

@@ -195,3 +195,14 @@ function expectation_value(ψ::FiniteMPS, O::ProjectionOperator,
     n = norm(ψ.AC[end])^2
     return sum(ens) / n
 end
+
+# Window
+# ------
+function expectation_value(Ψ::WindowMPS, windowH::Window, windowenvs::WindowEnv)
+    left_expval = expectation_value(Ψ.left, windowH.left, windowenvs.left)
+    middle_expval = expectation_value(Ψ.middle, windowH.middle, windowenvs.middle)
+    right_expval = expectation_value(Ψ.right, windowH.right, windowenvs.right)
+    return vcat(left_expval, middle_expval, right_expval)
+end
+
+#I need to think about expval with location

src/algorithms/groundstate/find_groundstate.jl

@@ -16,7 +16,8 @@ optimization algorithm will be attempted based on the supplied keywords.
 - `maxiter::Int`: maximum amount of iterations
 - `verbosity::Int`: display progress information
 """
-function find_groundstate(ψ::AbstractMPS, H, envs::Cache=environments(ψ, H);
+function find_groundstate(ψ::AbstractMPS, H,
+                          envs::Union{Cache,MultipleEnvironments}=environments(ψ, H);
                           tol=Defaults.tol, maxiter=Defaults.maxiter,
                           verbosity=Defaults.verbosity, trscheme=nothing)
     if isa(ψ, InfiniteMPS)
@@ -37,5 +38,31 @@ function find_groundstate(ψ::AbstractMPS, H, envs::Cache=environments(ψ, H);
     else
         throw(ArgumentError("Unknown input state type"))
     end
+    if isa(ψ, WindowMPS)
+        alg_infin = VUMPS(; tol=tol, verbosity=verbosity, maxiter=maxiter)
+        alg = Window(alg_infin, alg, alg_infin)
+    end
     return find_groundstate(ψ, H, alg, envs)
 end
+
+function find_groundstate!(state::WindowMPS{A,B,VL,VR}, H::Union{Window,LazySum{<:Window}},
+                           alg::Window, envs=environments(state, H)) where {A,B,VL,VR}
+    # first find infinite groundstates
+    if VL === WINDOW_VARIABLE
+        (gs_left, _) = find_groundstate(state.left, H.left, alg.left, envs.left)
+        state = WindowMPS(gs_left, state.middle, state.right)
+    end
+    if VR === WINDOW_VARIABLE
+        (gs_right, _) = find_groundstate(state.right, H.right, alg.right, envs.right)
+        state = WindowMPS(state.left, state.middle, gs_right)
+    end
+    # then find finite groundstate
+    state, _, delta = find_groundstate(state, H.middle, alg.middle,
+                                       finenv(envs, state))
+    return state, envs, delta
+end
+
+function find_groundstate(ψ::WindowMPS, H::Union{Window,LazySum{<:Window}}, alg::Window,
+                          envs...)
+    return find_groundstate!(copy(ψ), H, alg, envs...)
+end

src/algorithms/propagator/corvector.jl

@@ -153,42 +153,47 @@ function propagator(A::AbstractFiniteMPS, z, H::MPOHamiltonian,
     return v, init
 end
 
-function squaredenvs(state::AbstractFiniteMPS, H::MPOHamiltonian,
-                     envs=environments(state, H))
-    nH = conj(H) * H
-    L = length(state)
-
-    # to construct the squared caches we will first initialize environments
-    # then make all data invalid so it will be recalculated
-    # then initialize the right caches at the edge
-    ncocache = environments(state, nH)
-
+function squaredenvs(Ψ, H, envs, H2, envs2)
+    L = length(Ψ)
     # make sure the dependencies are incorrect, so data will be recalculated
     for i in 1:L
-        poison!(ncocache, i)
+        invalidate!(envs2, i)
     end
 
     # impose the correct boundary conditions
     # (important for comoving mps, should do nothing for finite mps)
     indmap = LinearIndices((H.odim, H.odim))
     @sync begin
         Threads.@spawn begin
-            nleft = leftenv(ncocache, 1, state)
+            nleft = leftenv(envs2, 1, Ψ)
             for i in 1:(H.odim), j in 1:(H.odim)
-                nleft[indmap[i, j]] = _contract_leftenv²(leftenv(envs, 1, state)[j],
-                                                         leftenv(envs, 1, state)[i])
+                nleft[indmap[i, j]] = _contract_leftenv²(leftenv(envs, 1, Ψ)[j],
+                                                         leftenv(envs, 1, Ψ)[i])
             end
         end
         Threads.@spawn begin
-            nright = rightenv(ncocache, L, state)
+            nright = rightenv(envs2, L, Ψ)
             for i in 1:(H.odim), j in 1:(H.odim)
-                nright[indmap[i, j]] = _contract_rightenv²(rightenv(envs, L, state)[j],
-                                                           rightenv(envs, L, state)[i])
+                nright[indmap[i, j]] = _contract_rightenv²(rightenv(envs, L, Ψ)[j],
+                                                           rightenv(envs, L, Ψ)[i])
             end
         end
     end
+    return H2, envs2
+end
+
+function squaredenvs(Ψ::AbstractFiniteMPS, H::MPOHamiltonian,
+                     envs::FinEnv=environments(Ψ, H))
+    # to construct the squared caches we will first initialize environments
+    # then make all data invalid so it will be recalculated
+    # then initialize the correct caches at the edge
+    nH = conj(H) * H
+    return squaredenvs(Ψ, H, envs, nH, environments(Ψ, nH))
+end
 
-    return nH, ncocache
+function squaredenvs(Ψ::WindowMPS, H::Window, envs::WindowEnv=environments(Ψ, H))
+    nH = Window(conj(H.left) * H.left, conj(H.middle) * H.middle, conj(H.right) * H.right)
+    return squaredenvs(Ψ, H.middle, envs, nH, environments(Ψ, nH))
 end
 
 function _contract_leftenv²(GL_top, GL_bot)