diff --git a/examples/windowmps.jl b/examples/windowmps.jl index 807c8c80..e93b855b 100644 --- a/examples/windowmps.jl +++ b/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 -(∑_{} Z_i Z_j + f(t) * ∑_{} 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")) diff --git a/src/MPSKit.jl b/src/MPSKit.jl index 882e4974..b5c251f7 100644 --- a/src/MPSKit.jl +++ b/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") diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl index 120f9880..583f4ca1 100644 --- a/src/algorithms/expval.jl +++ b/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 \ No newline at end of file diff --git a/src/algorithms/groundstate/find_groundstate.jl b/src/algorithms/groundstate/find_groundstate.jl index 75a9afba..aafb6ee1 100644 --- a/src/algorithms/groundstate/find_groundstate.jl +++ b/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 diff --git a/src/algorithms/propagator/corvector.jl b/src/algorithms/propagator/corvector.jl index 6e6e0f27..fbe2a0e5 100644 --- a/src/algorithms/propagator/corvector.jl +++ b/src/algorithms/propagator/corvector.jl @@ -153,19 +153,11 @@ 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 @@ -173,22 +165,35 @@ function squaredenvs(state::AbstractFiniteMPS, H::MPOHamiltonian, 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) diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl index 9f46898d..9a56381c 100644 --- a/src/algorithms/timestep/tdvp.jl +++ b/src/algorithms/timestep/tdvp.jl @@ -18,76 +18,82 @@ algorithm for time evolution. finalize::F = Defaults._finalize end -function timestep(ψ::InfiniteMPS, H, t::Number, dt::Number, alg::TDVP, - envs::Union{Cache,MultipleEnvironments}=environments(ψ, H); +function timestep(Ψ::InfiniteMPS, H, t::Number, dt::Number, alg::TDVP, + envs::Union{Cache,MultipleEnvironments}=environments(Ψ, H); leftorthflag=true) - temp_ACs = similar(ψ.AC) - temp_CRs = similar(ψ.CR) - - @static if Defaults.parallelize_sites - @sync for (loc, (ac, c)) in enumerate(zip(ψ.AC, ψ.CR)) - Threads.@spawn begin - h_ac = ∂∂AC(loc, ψ, H, envs) - temp_ACs[loc] = integrate(h_ac, ac, t, dt, alg.integrator) - end - - Threads.@spawn begin - h_c = ∂∂C(loc, ψ, H, envs) - temp_CRs[loc] = integrate(h_c, c, t, dt, alg.integrator) - end - end - else - for (loc, (ac, c)) in enumerate(zip(ψ.AC, ψ.CR)) - h_ac = ∂∂AC(loc, ψ, H, envs) + temp_ACs = similar(Ψ.AC) + temp_CRs = similar(Ψ.CR) + @sync for (loc, (ac, c)) in enumerate(zip(Ψ.AC, Ψ.CR)) + Threads.@spawn begin + h_ac = ∂∂AC(loc, Ψ, H, envs) temp_ACs[loc] = integrate(h_ac, ac, t, dt, alg.integrator) - h_c = ∂∂C(loc, ψ, H, envs) + end + + Threads.@spawn begin + h_c = ∂∂C(loc, Ψ, H, envs) temp_CRs[loc] = integrate(h_c, c, t, dt, alg.integrator) end end if leftorthflag - regauge!.(temp_ACs, temp_CRs; alg=TensorKit.QRpos()) - ψ′ = InfiniteMPS(temp_ACs, ψ.CR[end]; tol=alg.tolgauge, maxiter=alg.gaugemaxiter) + for loc in 1:length(Ψ) + # find AL that best fits these new Acenter and centers + QAc, _ = leftorth!(temp_ACs[loc]; alg=TensorKit.QRpos()) + Qc, _ = leftorth!(temp_CRs[loc]; alg=TensorKit.QRpos()) + @plansor temp_ACs[loc][-1 -2; -3] = QAc[-1 -2; 1] * conj(Qc[-3; 1]) + end + newΨ = InfiniteMPS(temp_ACs, Ψ.CR[end]; tol=alg.tolgauge, maxiter=alg.gaugemaxiter) + else - circshift!(temp_CRs, 1) - regauge!.(temp_CRs, temp_ACs; alg=TensorKit.LQpos()) - ψ′ = InfiniteMPS(ψ.CR[0], temp_ACs; tol=alg.tolgauge, maxiter=alg.gaugemaxiter) + for loc in 1:length(Ψ) + # find AR that best fits these new Acenter and centers + _, QAc = rightorth!(_transpose_tail(temp_ACs[loc]); alg=TensorKit.LQpos()) + _, Qc = rightorth!(temp_CRs[mod1(loc - 1, end)]; alg=TensorKit.LQpos()) + temp_ACs[loc] = _transpose_front(Qc' * QAc) + end + newΨ = InfiniteMPS(Ψ.CR[0], temp_ACs; tol=alg.tolgauge, maxiter=alg.gaugemaxiter) end - recalculate!(envs, ψ′) - return ψ′, envs + recalculate!(envs, newΨ) + return newΨ, envs end -function timestep!(ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP, - envs::Union{Cache,MultipleEnvironments}=environments(ψ, H)) +function ltr_sweep!(Ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP, + envs::Union{Cache,MultipleEnvironments}=environments(Ψ, H)) # sweep left to right - for i in 1:(length(ψ) - 1) - h_ac = ∂∂AC(i, ψ, H, envs) - ψ.AC[i] = integrate(h_ac, ψ.AC[i], t, dt / 2, alg.integrator) + for i in 1:(length(Ψ) - 1) + h_ac = ∂∂AC(i, Ψ, H, envs) + Ψ.AC[i] = integrate(h_ac, Ψ.AC[i], t, dt, alg.integrator) - h_c = ∂∂C(i, ψ, H, envs) - ψ.CR[i] = integrate(h_c, ψ.CR[i], t, -dt / 2, alg.integrator) + h_c = ∂∂C(i, Ψ, H, envs) + Ψ.CR[i] = integrate(h_c, Ψ.CR[i], t, -dt, alg.integrator) end # edge case - h_ac = ∂∂AC(length(ψ), ψ, H, envs) - ψ.AC[end] = integrate(h_ac, ψ.AC[end], t, dt / 2, alg.integrator) + h_ac = ∂∂AC(length(Ψ), Ψ, H, envs) + Ψ.AC[end] = integrate(h_ac, Ψ.AC[end], t, dt, alg.integrator) + + return Ψ, envs +end + +function rtl_sweep!(Ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP, + envs::Union{Cache,MultipleEnvironments}=environments(Ψ, H)) # sweep right to left - for i in length(ψ):-1:2 - h_ac = ∂∂AC(i, ψ, H, envs) - ψ.AC[i] = integrate(h_ac, ψ.AC[i], t + dt / 2, dt / 2, alg.integrator) + for i in length(Ψ):-1:2 + h_ac = ∂∂AC(i, Ψ, H, envs) + Ψ.AC[i] = integrate(h_ac, Ψ.AC[i], t, dt, alg.integrator) - h_c = ∂∂C(i - 1, ψ, H, envs) - ψ.CR[i - 1] = integrate(h_c, ψ.CR[i - 1], t + dt / 2, -dt / 2, alg.integrator) + h_c = ∂∂C(i - 1, Ψ, H, envs) + Ψ.CR[i - 1] = integrate(h_c, Ψ.CR[i - 1], t, -dt, alg.integrator) end # edge case - h_ac = ∂∂AC(1, ψ, H, envs) - ψ.AC[1] = integrate(h_ac, ψ.AC[1], t + dt / 2, dt / 2, alg.integrator) + h_ac = ∂∂AC(1, Ψ, H, envs) + Ψ.AC[1] = integrate(h_ac, Ψ.AC[1], t + dt, dt, alg.integrator) - return ψ, envs + return Ψ, envs end """ @@ -112,46 +118,46 @@ algorithm for time evolution. finalize::F = Defaults._finalize end -function timestep!(ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP2, - envs=environments(ψ, H)) +function ltr_sweep!(Ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP2, + envs=environments(Ψ, H)) # sweep left to right - for i in 1:(length(ψ) - 1) - ac2 = _transpose_front(ψ.AC[i]) * _transpose_tail(ψ.AR[i + 1]) - h_ac2 = ∂∂AC2(i, ψ, H, envs) + for i in 1:(length(Ψ) - 1) + ac2 = _transpose_front(Ψ.AC[i]) * _transpose_tail(Ψ.AR[i + 1]) + h_ac2 = ∂∂AC2(i, Ψ, H, envs) nac2 = integrate(h_ac2, ac2, t, dt / 2, alg.integrator) nal, nc, nar = tsvd!(nac2; trunc=alg.trscheme, alg=TensorKit.SVD()) - ψ.AC[i] = (nal, complex(nc)) - ψ.AC[i + 1] = (complex(nc), _transpose_front(nar)) + Ψ.AC[i] = (nal, complex(nc)) + Ψ.AC[i + 1] = (complex(nc), _transpose_front(nar)) - if i != (length(ψ) - 1) - ψ.AC[i + 1] = integrate(∂∂AC(i + 1, ψ, H, envs), ψ.AC[i + 1], t, -dt / 2, + if i != (length(Ψ) - 1) + Ψ.AC[i + 1] = integrate(∂∂AC(i + 1, Ψ, H, envs), Ψ.AC[i + 1], t, -dt / 2, alg.integrator) end end + return Ψ, envs +end + +function rtl_sweep!(Ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP2, + envs=environments(Ψ, H)) + # sweep right to left - for i in length(ψ):-1:2 - ac2 = _transpose_front(ψ.AL[i - 1]) * _transpose_tail(ψ.AC[i]) - h_ac2 = ∂∂AC2(i - 1, ψ, H, envs) + for i in length(Ψ):-1:2 + ac2 = _transpose_front(Ψ.AL[i - 1]) * _transpose_tail(Ψ.AC[i]) + h_ac2 = ∂∂AC2(i - 1, Ψ, H, envs) nac2 = integrate(h_ac2, ac2, t + dt / 2, dt / 2, alg.integrator) nal, nc, nar = tsvd!(nac2; trunc=alg.trscheme, alg=TensorKit.SVD()) - ψ.AC[i - 1] = (nal, complex(nc)) - ψ.AC[i] = (complex(nc), _transpose_front(nar)) + Ψ.AC[i - 1] = (nal, complex(nc)) + Ψ.AC[i] = (complex(nc), _transpose_front(nar)) if i != 2 - ψ.AC[i - 1] = integrate(∂∂AC(i - 1, ψ, H, envs), ψ.AC[i - 1], t + dt / 2, + Ψ.AC[i - 1] = integrate(∂∂AC(i - 1, Ψ, H, envs), Ψ.AC[i - 1], t + dt / 2, -dt / 2, alg.integrator) end end - return ψ, envs -end - -#copying version -function timestep(ψ::AbstractFiniteMPS, H, time::Number, timestep::Number, - alg::Union{TDVP,TDVP2}, envs=environments(ψ, H); kwargs...) - return timestep!(copy(ψ), H, time, timestep, alg, envs; kwargs...) + return Ψ, envs end diff --git a/src/algorithms/timestep/timestep.jl b/src/algorithms/timestep/timestep.jl new file mode 100644 index 00000000..2dac6929 --- /dev/null +++ b/src/algorithms/timestep/timestep.jl @@ -0,0 +1,28 @@ +# TDVP + +function timestep!(Ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::Union{TDVP,TDVP2}, + envs::Union{Cache,MultipleEnvironments}=environments(Ψ, H)) + Ψ, envs = ltr_sweep!(Ψ, H, t, dt / 2, alg, envs) + Ψ, envs = rtl_sweep!(Ψ, H, t + dt / 2, dt / 2, alg, envs) + + return Ψ, envs +end + +#copying version +function timestep(Ψ::AbstractFiniteMPS, H, time::Number, timestep::Number, + alg::Union{TDVP,TDVP2}, envs=environments(Ψ, H); kwargs...) + return timestep!(copy(Ψ), H, time, timestep, alg, envs; kwargs...) +end + +function timestep(Ψ::WindowMPS, H::Union{Window,LazySum{<:Window}}, time::Number, + timestep::Number, + alg::WindowTDVP, envs=environments(Ψ, H); kwargs...) + return timestep!(copy(Ψ), H, time, timestep, alg, envs; kwargs...) +end + +# Time MPO +#= +function timestep(Ψ::FiniteMPS, H, t::Number, dt::Number, alg, + envs::Union{Cache,MultipleEnvironments}=environments(Ψ, H)) +end +=# \ No newline at end of file diff --git a/src/algorithms/timestep/windowtdvp.jl b/src/algorithms/timestep/windowtdvp.jl new file mode 100644 index 00000000..5e1b38d3 --- /dev/null +++ b/src/algorithms/timestep/windowtdvp.jl @@ -0,0 +1,41 @@ +""" + WindowTDVP{A} <: Algorithm + +[Mixed TDVP](https://arxiv.org/abs/2007.15035) algorithm for time evolution. + +# Fields +- `left::A`: algorithm to do the timestep of the infinite part of the WindowMPS. +- `middle::B`: algorithm to do the timestep of the finite part of the WindowMPS. This can be `TDVP2` to expand the bonddimension. +- `right::C`: algorithm to do the timestep of the right part of the WindowMPS. By default the same as left +- `finalize::F`: user-supplied function which is applied after each timestep, with + signature `finalize(t, Ψ, H, envs) -> Ψ, envs`. Can be used to enlarge the bond dimension of the infinite part. +""" +@kwdef struct WindowTDVP{A,B,C,F} <: Algorithm + left::A = TDVP() + middle::B = TDVP() + right::C = left + finalize::F = Defaults._finalize +end + +function timestep!(Ψ::WindowMPS{A,B,VL,VR}, H::Union{Window,LazySum{<:Window}}, t::Number, + dt::Number, alg::WindowTDVP, env=environments(Ψ, H)) where {A,B,VL,VR} + + #first evolve left state + if VL === WINDOW_VARIABLE + nleft, _ = timestep(Ψ.left, H.left, t, dt, alg.left, env.left; leftorthflag=true) #env gets updated in place + Ψ = WindowMPS(nleft, Ψ.middle, Ψ.right) + end + + Ψ, env = ltr_sweep!(Ψ, H.middle, t, dt / 2, alg.middle, env) + + #then evolve right state + if VR === WINDOW_VARIABLE + nright, _ = timestep(Ψ.right, H.right, t + dt, dt, alg.right, env.right; + leftorthflag=false) # env gets updated in place + Ψ = WindowMPS(Ψ.left, Ψ.middle, nright) + end + + Ψ, env = rtl_sweep!(Ψ, H.middle, t + dt / 2, dt / 2, alg.middle, env) + + return Ψ, env +end \ No newline at end of file diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl index adc486a8..4e967ebd 100644 --- a/src/algorithms/toolbox.jl +++ b/src/algorithms/toolbox.jl @@ -210,6 +210,13 @@ function variance(ψ, H::LazySum, envs=environments(ψ, sum(H))) return variance(ψ, sum(H), envs) end +function variance(ψ::WindowMPS, H::Union{MPOHamiltonian,Window}, envs=environments(ψ, H)) + #tricky to define + H2, nenvs = squaredenvs(ψ, H, envs) + return real(expectation_value(ψ, H2, 1:length(ψ), nenvs) - + expectation_value(ψ, H, 1:length(ψ), envs)^2) +end + """ You can impose periodic boundary conditions on an mpo-hamiltonian (for a given size) That creates a new mpo-hamiltonian with larger bond dimension diff --git a/src/environments/FinEnv.jl b/src/environments/finenv.jl similarity index 82% rename from src/environments/FinEnv.jl rename to src/environments/finenv.jl index b45eaa59..cefedee9 100644 --- a/src/environments/FinEnv.jl +++ b/src/environments/finenv.jl @@ -1,5 +1,5 @@ """ - FinEnv keeps track of the environments for FiniteMPS / WindowMPS + FinEnv keeps track of the environments for FiniteMPS It automatically checks if the queried environment is still correctly cached and if not - recalculates if above is set to nothing, above === below. @@ -18,6 +18,9 @@ struct FinEnv{A,B,C,D} <: Cache rightenvs::Vector{D} end +#do we want to call these function environments? or do we want to restrict it to +# only the constructor environments(state,ham)? + function environments(below, t::Tuple, args...; kwargs...) return environments(below, t[1], t[2], args...; kwargs...) end; @@ -80,20 +83,7 @@ function MPSKit.environments(below::FiniteMPS{S}, O::DenseMPO, above=nothing) wh return environments(below, O, above, leftstart, rightstart) end -#extract the correct leftstart/rightstart for WindowMPS -function environments(state::WindowMPS, O::Union{SparseMPO,MPOHamiltonian,DenseMPO}, - above=nothing; lenvs=environments(state.left_gs, O), - renvs=environments(state.right_gs, O)) - return environments(state, O, above, copy(leftenv(lenvs, 1, state.left_gs)), - copy(rightenv(renvs, length(state), state.right_gs))) -end - -function environments(below::S, above::S) where {S<:Union{FiniteMPS,WindowMPS}} - S isa WindowMPS && - (above.left_gs == below.left_gs || throw(ArgumentError("left gs differs"))) - S isa WindowMPS && - (above.right_gs == below.right_gs || throw(ArgumentError("right gs differs"))) - +function environments(below::FiniteMPS, above::FiniteMPS) opp = fill(nothing, length(below)) return environments(below, opp, above, l_LL(above), r_RR(above)) end @@ -105,7 +95,8 @@ function environments(state::Union{FiniteMPS,WindowMPS}, opp::ProjectionOperator end #notify the cache that we updated in-place, so it should invalidate the dependencies -function poison!(ca::FinEnv, ind) +# this forces the transfers to be recalculated lazily +function invalidate!(ca::FinEnv, ind) ca.ldependencies[ind] = similar(ca.ldependencies[ind]) return ca.rdependencies[ind] = similar(ca.rdependencies[ind]) end @@ -113,9 +104,9 @@ end #rightenv[ind] will be contracteable with the tensor on site [ind] function rightenv(ca::FinEnv, ind, state) a = findfirst(i -> !(state.AR[i] === ca.rdependencies[i]), length(state):-1:(ind + 1)) - a = a == nothing ? nothing : length(state) - a + 1 + a = isnothing(a) ? nothing : length(state) - a + 1 - if a != nothing + if !isnothing(a) #we need to recalculate for j in a:-1:(ind + 1) above = isnothing(ca.above) ? state.AR[j] : ca.above.AR[j] @@ -131,7 +122,7 @@ end function leftenv(ca::FinEnv, ind, state) a = findfirst(i -> !(state.AL[i] === ca.ldependencies[i]), 1:(ind - 1)) - if a != nothing + if !isnothing(a) #we need to recalculate for j in a:(ind - 1) above = isnothing(ca.above) ? state.AL[j] : ca.above.AL[j] diff --git a/src/environments/multipleenv.jl b/src/environments/multipleenv.jl index 108ea99e..3c75ee0e 100644 --- a/src/environments/multipleenv.jl +++ b/src/environments/multipleenv.jl @@ -30,15 +30,29 @@ end # return MultipleEnvironments(H, broadcast(o -> environments(state, o), H.opps)) # end; -function environments(st::WindowMPS, - H::LazySum; - lenvs=environments(st.left_gs, H), - renvs=environments(st.right_gs, H)) - return MultipleEnvironments(H, - map((op, sublenv, subrenv) -> environments(st, op; - lenvs=sublenv, - renvs=subrenv), - H.ops, lenvs, renvs)) +#=========================================================================================== +Utility +===========================================================================================# +function Base.getproperty(ca::MultipleEnvironments{<:LazySum,<:WindowEnv}, sym::Symbol) + if sym === :left || sym === :middle || sym === :right + #extract the left/right parts + return MultipleEnvironments(getproperty(ca.opp, sym), + map(x -> getproperty(x, sym), ca)) + else + return getfield(ca, sym) + end +end + +function Base.getproperty(envs::MultipleEnvironments, prop::Symbol) + if prop === :solver + return map(env -> env.solver, envs) + else + return getfield(envs, prop) + end +end + +function finenv(ca::MultipleEnvironments{<:LazySum,<:WindowEnv}, ψ::WindowMPS) + return MultipleEnvironments(ca.opp.middle, map(x -> finenv(x, ψ), ca)) end # we need to define how to recalculate @@ -51,12 +65,3 @@ function recalculate!(env::MultipleEnvironments, args...; kwargs...) end return env end - -#maybe this can be used to provide compatibility with existing code? -function Base.getproperty(envs::MultipleEnvironments, prop::Symbol) - if prop === :solver - return map(env -> env.solver, envs) - else - return getfield(envs, prop) - end -end diff --git a/src/environments/windowenv.jl b/src/environments/windowenv.jl new file mode 100644 index 00000000..ff34e3e0 --- /dev/null +++ b/src/environments/windowenv.jl @@ -0,0 +1,108 @@ +""" + WindowEnv keeps track of the environments for WindowMPS + Changes in the infinite environments are checked and taken into account whenever the environment of the finite part is queried. + +""" +struct WindowEnv{A,B,C,D} <: Cache + window::Window{A,B,C} + + linfdeps::PeriodicArray{D} #the data we used to calculate leftenvs/rightenvs + rinfdeps::PeriodicArray{D} +end + +#automatically construct the correct leftstart/rightstart for a WindowMPS +# copying the vector where the tensors reside in makes it have another memory adress, while keeping the references of the elements +function environments(ψ::WindowMPS, O::Window, above=nothing; + lenvs=environments(ψ.left, O.left), + renvs=environments(ψ.right, O.right)) + fin_env = environments(ψ, O.middle, above, lenvs, renvs) + return WindowEnv(Window(lenvs, fin_env, renvs), copy(ψ.left.AL), copy(ψ.right.AR)) +end + +function environments(ψ::WindowMPS, O::Union{DenseMPO,SparseMPO,MPOHamiltonian}, above, + lenvs::Union{PerMPOInfEnv,MPOHamInfEnv}, + renvs::Union{PerMPOInfEnv,MPOHamInfEnv}) + return environments(ψ, O, above, leftenv(lenvs, 1, ψ.left), + rightenv(renvs, length(ψ), ψ.right)) +end + +#Backwards compability +function environments(ψ::WindowMPS{A,B,WINDOW_FIXED,WINDOW_FIXED}, H::MPOHamiltonian, + above=nothing) where {A,B} + length(H) == length(ψ.left) || + throw(ArgumentError("length(ψ.left) != length(H), use H=Window(Hleft,Hmiddle,Hright) instead")) + length(H) == length(ψ.right) || + throw(ArgumentError("length(ψ.right) != length(H), use H=Window(Hleft,Hmiddle,Hright) instead")) + return environments(ψ, H, above, environments(ψ.left, H), environments(ψ.right, H)) +end + +function environments(below::WindowMPS, above::WindowMPS) + above.left == below.left || throw(ArgumentError("left gs differs")) + above.right == below.right || throw(ArgumentError("right gs differs")) + + opp = fill(nothing, length(below)) + return environments(below, opp, above, l_LL(above), r_RR(above)) +end + +#=========================================================================================== +Utility +===========================================================================================# +function Base.getproperty(ca::WindowEnv, sym::Symbol) #under review + if sym === :left || sym === :middle || sym === :right + return getfield(ca.window, sym) + elseif sym === :opp #this is for derivatives. Could we remove opp field from derivatives? + return getfield(ca.window.middle, sym) + else + return getfield(ca, sym) + end +end + +# when accesing the finite part of the env, use this function +function finenv(ca::WindowEnv, ψ::WindowMPS{A,B,VL,VR}) where {A,B,VL,VR} + VL === WINDOW_FIXED || check_leftinfenv!(ca, ψ) + VR === WINDOW_FIXED || check_rightinfenv!(ca, ψ) + return ca.middle +end + +#notify the cache that we updated in-place, so it should invalidate the dependencies +invalidate!(ca::WindowEnv, ind) = invalidate!(ca.middle, ind) + +# Check the infinite envs and recalculate the right start +function check_rightinfenv!(ca::WindowEnv, ψ::WindowMPS) + if !all(ca.rinfdeps .=== ψ.right.AR) + invalidate!(ca, length(ψ.middle)) #forces transfers to be recalculated lazily + + ca.middle.rightenvs[end] = rightenv(ca.right, 0, ψ.right) + ca.rinfdeps .= ψ.right.AR + end +end + +function check_leftinfenv!(ca::WindowEnv, ψ::WindowMPS) + if !all(ca.linfdeps .=== ψ.left.AL) + invalidate!(ca, 1) #forces transfers to be recalculated lazily + + ca.middle.leftenvs[1] = leftenv(ca.left, length(ψ.left) + 1, ψ.left) + ca.linfdeps .= ψ.left.AL + end +end + +#should be used to calculate expvals for operators over the boundary, also need views to work for this +function leftenv(ca::WindowEnv, ind, ψ::WindowMPS) + if ind < 1 + return leftenv(ca.left, ind, ψ.left) + elseif ind > length(ψ) + return leftenv(ca.right, ind, ψ.right) + else + return leftenv(finenv(ca, ψ), ind, ψ) + end +end + +function rightenv(ca::WindowEnv, ind, ψ::WindowMPS) + if ind < 1 + return rightenv(ca.left, ind, ψ.left) + elseif ind > length(ψ) + return rightenv(ca.right, ind, ψ.right) + else + return rightenv(finenv(ca, ψ), ind, ψ) + end +end diff --git a/src/operators/lazysum.jl b/src/operators/lazysum.jl index 5d911cbf..4439988f 100644 --- a/src/operators/lazysum.jl +++ b/src/operators/lazysum.jl @@ -34,13 +34,7 @@ LazySum(ops::AbstractVector, fs::AbstractVector) = LazySum(map(MultipliedOperato # wrapper around _eval_at safe_eval(::TimeDependent, x::LazySum, t::Number) = map(O -> _eval_at(O, t), x) -function safe_eval(::TimeDependent, x::LazySum) - throw(ArgumentError("attempting to evaluate time-dependent LazySum without specifiying a time")) -end safe_eval(::NotTimeDependent, x::LazySum) = sum(_eval_at, x) -function safe_eval(::NotTimeDependent, x::LazySum, t::Number) - throw(ArgumentError("attempting to evaluate time-independent LazySum at time")) -end # For users # using (t) should return NotTimeDependent LazySum @@ -57,3 +51,12 @@ Base.:+(op1, op2::LazySum) = LazySum(op1) + op2 Base.:+(op1::MultipliedOperator, op2::MultipliedOperator) = LazySum([op1, op2]) Base.repeat(x::LazySum, args...) = LazySum(repeat.(x, args...)) + +function Base.getproperty(sumops::LazySum{<:Window}, sym::Symbol) + if sym === :left || sym === :middle || sym === :right + #extract the left/right parts + return map(x -> getproperty(x, sym), sumops) + else + return getfield(sumops, sym) + end +end diff --git a/src/operators/multipliedoperator.jl b/src/operators/multipliedoperator.jl index 6f9100e6..a51c241e 100644 --- a/src/operators/multipliedoperator.jl +++ b/src/operators/multipliedoperator.jl @@ -50,3 +50,8 @@ Base.:*(op::UntimedOperator, g::Function) = MultipliedOperator(op.op, t -> g(t) function environments(st, x::MultipliedOperator, args...; kwargs...) return environments(st, x.op, args...; kwargs...) end + +function MultipliedOperator(x::Window, f) + return Window(MultipliedOperator(x.left, f), MultipliedOperator(x.middle, f), + MultipliedOperator(x.right, f)) +end \ No newline at end of file diff --git a/src/states/orthoview.jl b/src/states/orthoview.jl index 981e64fc..4273a955 100644 --- a/src/states/orthoview.jl +++ b/src/states/orthoview.jl @@ -10,8 +10,8 @@ end function Base.getindex(v::ALView{<:WindowMPS,E}, i::Int)::E where {E} i <= length(v.parent) || throw(ArgumentError("out of bounds")) - i < 1 && return v.parent.left_gs.AL[i] - return ALView(v.parent.window)[i] + i < 1 && return v.parent.left.AL[i] + return ALView(v.parent.middle)[i] end Base.getindex(v::ALView{<:Multiline}, i::Int, j::Int) = v.parent[i].AL[j] @@ -32,8 +32,8 @@ end function Base.getindex(v::ARView{<:WindowMPS,E}, i::Int)::E where {E} i >= 1 || throw(ArgumentError("out of bounds")) - i > length(v.parent) && return v.parent.right_gs.AR[i] - return ARView(v.parent.window)[i] + i > length(v.parent) && return v.parent.right.AR[i] + return ARView(v.parent.middle)[i] end Base.getindex(v::ARView{<:Multiline}, i::Int, j::Int) = v.parent[i].AR[j] @@ -78,9 +78,9 @@ function Base.setindex!(v::CRView{<:FiniteMPS}, vec, i::Int) return setindex!(v.parent.CLs, vec, i + 1) end -Base.getindex(v::CRView{<:WindowMPS}, i::Int) = CRView(v.parent.window)[i] +Base.getindex(v::CRView{<:WindowMPS}, i::Int) = CRView(v.parent.middle)[i] function Base.setindex!(v::CRView{<:WindowMPS}, vec, i::Int) - return setindex!(CRView(v.parent.window), vec, i) + return setindex!(CRView(v.parent.middle), vec, i) end Base.getindex(v::CRView{<:Multiline}, i::Int, j::Int) = v.parent[i].CR[j] function Base.setindex!(v::CRView{<:Multiline}, vec, i::Int, j::Int) @@ -144,10 +144,10 @@ end function Base.getindex(v::ACView{<:WindowMPS,E}, i::Int)::E where {E} (i >= 1 && i <= length(v.parent)) || throw(ArgumentError("out of bounds")) - return ACView(v.parent.window)[i] + return ACView(v.parent.middle)[i] end function Base.setindex!(v::ACView{<:WindowMPS}, vec, i::Int) - return setindex!(ACView(v.parent.window), vec, i) + return setindex!(ACView(v.parent.middle), vec, i) end Base.getindex(v::ACView{<:Multiline}, i::Int, j::Int) = v.parent[i].AC[j] diff --git a/src/states/window.jl b/src/states/window.jl new file mode 100644 index 00000000..92954d41 --- /dev/null +++ b/src/states/window.jl @@ -0,0 +1,28 @@ +" + Window(left,middle,right) + + general struct of an object with a left, middle and right part. +" +struct Window{L,M,R} + left::L + middle::M + right::R +end + +function Base.:+(window1::Window, window2::Window) + return Window(window1.left + window2.left, window1.middle + window2.middle, + window1.right + window2.right) +end + +# Holy traits +TimeDependence(x::Window) = istimed(x) ? TimeDependent() : NotTimeDependent() +istimed(x::Window) = istimed(x.left) || istimed(x.middle) || istimed(x.right) + +function _eval_at(x::Window, args...) + return Window(_eval_at(x.left, args...), _eval_at(x.middle, args...), + _eval_at(x.right, args...)) +end +safe_eval(::TimeDependent, x::Window, t::Number) = _eval_at(x, t) + +# For users +(x::Window)(t::Number) = safe_eval(x, t) \ No newline at end of file diff --git a/src/states/windowmps.jl b/src/states/windowmps.jl index 7225ddb2..88e2a066 100644 --- a/src/states/windowmps.jl +++ b/src/states/windowmps.jl @@ -1,85 +1,88 @@ +const WINDOW_FIXED = :F +const WINDOW_VARIABLE = :V + """ WindowMPS{A<:GenericMPSTensor,B<:MPSBondTensor} <: AbstractFiniteMPS -Type that represents a finite Matrix Product State embedded in an infinte Matrix Product State. +Type that represents a finite Matrix Product State embedded in an infinite Matrix Product State. ## Fields -- `left_gs::InfiniteMPS` -- left infinite environment -- `window::FiniteMPS` -- finite window Matrix Product State -- `right_gs::InfiniteMPS` -- right infinite environment +- `left::InfiniteMPS` -- left infinite state +- `middle::FiniteMPS` -- finite state in the middle +- `right::InfiniteMPS` -- right infinite state --- ## Constructors - WindowMPS(left_gs::InfiniteMPS, window_state::FiniteMPS, [right_gs::InfiniteMPS]) - WindowMPS(left_gs::InfiniteMPS, window_tensors::AbstractVector, [right_gs::InfiniteMPS]) + WindowMPS(left::InfiniteMPS, middle::FiniteMPS, right::InfiniteMPS) + WindowMPS(left::InfiniteMPS, middle_tensors::AbstractVector, right::InfiniteMPS) WindowMPS([f, eltype], physicalspaces::Vector{<:Union{S, CompositeSpace{S}}, - virtualspaces::Vector{<:Union{S, CompositeSpace{S}}, left_gs::InfiniteMPS, - [right_gs::InfiniteMPS]) + virtualspaces::Vector{<:Union{S, CompositeSpace{S}}, left::InfiniteMPS, + right_gs::InfiniteMPS) WindowMPS([f, eltype], physicalspaces::Vector{<:Union{S,CompositeSpace{S}}}, - maxvirtualspace::S, left_gs::InfiniteMPS, [right_gs::InfiniteMPS]) + maxvirtualspace::S, left::InfiniteMPS, right_gs::InfiniteMPS) -Construct a WindowMPS via a specification of left and right infinite environment, and either -a window state or a vector of tensors to construct the window. Alternatively, it is possible +Construct a WindowMPS via a specification of left and right infinite state, and either +a middle state or a vector of tensors to construct the middle. Alternatively, it is possible to supply the same arguments as for the constructor of [`FiniteMPS`](@ref), followed by a -left (and right) environment to construct the WindowMPS in one step. - -!!! note - By default, the right environment is chosen to be equal to the left, however no copy is - made. In this case, changing the left state will also affect the right state. +left and right state to construct the WindowMPS in one step. - WindowMPS(state::InfiniteMPS, L::Int) + WindowMPS(state::InfiniteMPS, L::Int; copyright=false) Construct a WindowMPS from an InfiniteMPS, by promoting a region of length `L` to a -`FiniteMPS`. +`FiniteMPS`. Note that by default the right state is not copied (and thus .left === .right). + +Options for fixing the left and right infinite state (i.e. so they don't get time evolved) +can be done via the Boolean keyword arguments `fixleft` and `fixright`. """ -struct WindowMPS{A<:GenericMPSTensor,B<:MPSBondTensor} <: AbstractFiniteMPS - left_gs::InfiniteMPS{A,B} - window::FiniteMPS{A,B} - right_gs::InfiniteMPS{A,B} +struct WindowMPS{A<:GenericMPSTensor,B<:MPSBondTensor,VL,VR} <: AbstractFiniteMPS + window::Window{InfiniteMPS{A,B},FiniteMPS{A,B},InfiniteMPS{A,B}} function WindowMPS(ψₗ::InfiniteMPS{A,B}, ψₘ::FiniteMPS{A,B}, - ψᵣ::InfiniteMPS{A,B}=copy(ψₗ)) where {A<:GenericMPSTensor, - B<:MPSBondTensor} + ψᵣ::InfiniteMPS{A,B}; fixleft::Bool=false, + fixright::Bool=false) where {A,B} left_virtualspace(ψₗ, 0) == left_virtualspace(ψₘ, 0) && right_virtualspace(ψₘ, length(ψₘ)) == right_virtualspace(ψᵣ, length(ψₘ)) || throw(SpaceMismatch("Mismatch between window and environment virtual spaces")) - return new{A,B}(ψₗ, ψₘ, ψᵣ) + VL = fixleft ? WINDOW_FIXED : WINDOW_VARIABLE + VR = fixright ? WINDOW_FIXED : WINDOW_VARIABLE + return new{A,B,VL,VR}(Window(ψₗ, ψₘ, ψᵣ)) end end #=========================================================================================== Constructors ===========================================================================================# - function WindowMPS(ψₗ::InfiniteMPS, site_tensors::AbstractVector{<:GenericMPSTensor}, - ψᵣ::InfiniteMPS=ψₗ) - return WindowMPS(ψₗ, FiniteMPS(site_tensors), ψᵣ) + ψᵣ::InfiniteMPS; kwargs...) + return WindowMPS(ψₗ, FiniteMPS(site_tensors), ψᵣ; kwargs...) end +#perhaps we want to consider not using the FiniteMPS constructor since I believe these constructs +#the spaces so that the edge virtual sapces are one dimensional. function WindowMPS(f, elt, physspaces::Vector{<:Union{S,CompositeSpace{S}}}, maxvirtspace::S, ψₗ::InfiniteMPS, - ψᵣ::InfiniteMPS=ψₗ) where {S<:ElementarySpace} + ψᵣ::InfiniteMPS; kwargs...) where {S<:ElementarySpace} ψₘ = FiniteMPS(f, elt, physspaces, maxvirtspace; left=left_virtualspace(ψₗ, 0), right=right_virtualspace(ψᵣ, length(physspaces))) - return WindowMPS(ψₗ, ψₘ, ψᵣ) + return WindowMPS(ψₗ, ψₘ, ψᵣ; kwargs...) end function WindowMPS(physspaces::Vector{<:Union{S,CompositeSpace{S}}}, maxvirtspace::S, - ψₗ::InfiniteMPS, ψᵣ::InfiniteMPS=ψₗ) where {S<:ElementarySpace} - return WindowMPS(rand, Defaults.eltype, physspaces, maxvirtspace, ψₗ, ψᵣ) + ψₗ::InfiniteMPS, ψᵣ::InfiniteMPS; kwargs...) where {S<:ElementarySpace} + return WindowMPS(rand, Defaults.eltype, physspaces, maxvirtspace, ψₗ, ψᵣ; kwargs...) end function WindowMPS(f, elt, physspaces::Vector{<:Union{S,CompositeSpace{S}}}, virtspaces::Vector{S}, ψₗ::InfiniteMPS, - ψᵣ::InfiniteMPS=ψₗ) where {S<:ElementarySpace} + ψᵣ::InfiniteMPS; kwargs...) where {S<:ElementarySpace} ψₘ = FiniteMPS(f, elt, physspaces, virtspaces) - return WindowMPS(ψₗ, ψₘ, ψᵣ) + return WindowMPS(ψₗ, ψₘ, ψᵣ; kwargs...) end function WindowMPS(physspaces::Vector{<:Union{S,CompositeSpace{S}}}, virtspaces::Vector{S}, - ψₗ::InfiniteMPS, ψᵣ::InfiniteMPS=ψₗ) where {S<:ElementarySpace} - return WindowMPS(rand, Defaults.eltype, physspaces, virtspaces, ψₗ, ψᵣ) + ψₗ::InfiniteMPS, ψᵣ::InfiniteMPS; kwargs...) where {S<:ElementarySpace} + return WindowMPS(rand, Defaults.eltype, physspaces, virtspaces, ψₗ, ψᵣ; kwargs...) end function WindowMPS(f, elt, P::ProductSpace, args...; kwargs...) @@ -96,7 +99,7 @@ function WindowMPS(N::Int, V::VectorSpace, args...; kwargs...) return WindowMPS(fill(V, N), args...; kwargs...) end -function WindowMPS(ψ::InfiniteMPS{A,B}, L::Int) where {A,B} +function WindowMPS(ψ::InfiniteMPS{A,B}, L::Int; copyright=false, kwargs...) where {A,B} CLs = Vector{Union{Missing,B}}(missing, L + 1) ALs = Vector{Union{Missing,A}}(missing, L) ARs = Vector{Union{Missing,A}}(missing, L) @@ -107,20 +110,39 @@ function WindowMPS(ψ::InfiniteMPS{A,B}, L::Int) where {A,B} ACs .= ψ.AC[1:L] CLs .= ψ.CR[0:L] - return WindowMPS(ψ, FiniteMPS(ALs, ARs, ACs, CLs), ψ) + return WindowMPS(ψ, FiniteMPS(ALs, ARs, ACs, CLs), copyright ? copy(ψ) : ψ; kwargs...) end #=========================================================================================== Utility ===========================================================================================# +function Base.getproperty(ψ::WindowMPS, sym::Symbol) + if sym === :left || sym === :middle || sym === :right + return getfield(ψ.window, sym) + elseif sym === :AL + return ALView(ψ) + elseif sym === :AR + return ARView(ψ) + elseif sym === :AC + return ACView(ψ) + elseif sym === :CR + return CRView(ψ) + else + return getfield(ψ, sym) + end +end -function Base.copy(ψ::WindowMPS) - return WindowMPS(copy(ψ.left_gs), copy(ψ.window), copy(ψ.right_gs)) +function Base.copy(ψ::WindowMPS{A,B,VL,VR}) where {A,B,VL,VR} + left = VL === WINDOW_VARIABLE ? copy(ψ.left) : ψ.left + fixleft = VL === WINDOW_VARIABLE ? false : true + right = VR === WINDOW_VARIABLE ? copy(ψ.right) : ψ.right + fixright = VR === WINDOW_VARIABLE ? false : true + return WindowMPS(left, copy(ψ.middle), right; fixleft=fixleft, fixright=fixright) end # not sure about the underlying methods... -Base.length(ψ::WindowMPS) = length(ψ.window) -Base.size(ψ::WindowMPS, i...) = size(ψ.window, i...) +Base.length(ψ::WindowMPS) = length(ψ.middle) +Base.size(ψ::WindowMPS, i...) = size(ψ.middle, i...) Base.eltype(::Type{<:WindowMPS{A}}) where {A} = A Base.checkbounds(::Type{Bool}, ψ::WindowMPS, i::Integer) = true @@ -130,13 +152,13 @@ bond_type(::Type{<:WindowMPS{<:Any,B}}) where {B} = B TensorKit.space(ψ::WindowMPS, n::Integer) = space(ψ.AC[n], 2) for f in (:physicalspace, :left_virtualspace, :right_virtualspace) - @eval $f(ψ::WindowMPS, n::Integer) = n < 1 ? $f(ψ.left_gs, n) : - n > length(ψ) ? $f(ψ.right_gs, n - length(ψ)) : - $f(ψ.window, n) + @eval $f(ψ::WindowMPS, n::Integer) = n < 1 ? $f(ψ.left, n) : + n > length(ψ) ? $f(ψ.right, n - length(ψ)) : + $f(ψ.middle, n) end function physicalspace(ψ::WindowMPS) - return WindowArray(physicalspace(ψ.left_gs), physicalspace(ψ.window), - physicalspace(ψ.right_gs)) + return WindowArray(physicalspace(ψ.left), physicalspace(ψ.middle), + physicalspace(ψ.right)) end r_RR(ψ::WindowMPS) = r_RR(ψ.right_gs, length(ψ)) l_LL(ψ::WindowMPS) = l_LL(ψ.left_gs, 1) @@ -155,35 +177,35 @@ function Base.getproperty(ψ::WindowMPS, prop::Symbol) end end -max_Ds(ψ::WindowMPS) = max_Ds(ψ.window) +max_Ds(ψ::WindowMPS) = max_Ds(ψ.middle) Base.:*(ψ::WindowMPS, a::Number) = rmul!(copy(ψ), a) Base.:*(a::Number, ψ::WindowMPS) = lmul!(a, copy(ψ)) function TensorKit.lmul!(a::Number, ψ::WindowMPS) - lmul!(a, ψ.window) + lmul!(a, ψ.middle) return ψ end function TensorKit.rmul!(ψ::WindowMPS, a::Number) - rmul!(ψ.window, a) + rmul!(ψ.middle, a) return ψ end function TensorKit.dot(ψ₁::WindowMPS, ψ₂::WindowMPS) length(ψ₁) == length(ψ₂) || throw(ArgumentError("MPS with different length")) - ψ₁.left_gs == ψ₂.left_gs || - dot(ψ₁.left_gs, ψ₂.left_gs) ≈ 1 || + ψ₁.left == ψ₂.left || + dot(ψ₁.left, ψ₂.left) ≈ 1 || throw(ArgumentError("left InfiniteMPS are different")) - ψ₁.right_gs == ψ₂.right_gs || - dot(ψ₁.right_gs, ψ₂.right_gs) ≈ 1 || + ψ₁.right == ψ₂.right || + dot(ψ₁.right, ψ₂.right) ≈ 1 || throw(ArgumentError("right InfiniteMPS are different")) ρr = TransferMatrix(ψ₂.AR[2:end], ψ₁.AR[2:end]) * r_RR(ψ₂) return tr(_transpose_front(ψ₁.AC[1])' * _transpose_front(ψ₂.AC[1]) * ρr) end -TensorKit.norm(ψ::WindowMPS) = norm(ψ.window) +TensorKit.norm(ψ::WindowMPS) = norm(ψ.middle) -TensorKit.normalize!(ψ::WindowMPS) = normalize!(ψ.window) +TensorKit.normalize!(ψ::WindowMPS) = normalize!(ψ.middle) TensorKit.normalize(ψ::WindowMPS) = normalize!(copy(ψ)) diff --git a/src/operators/timedependence.jl b/src/utility/timedependence.jl similarity index 63% rename from src/operators/timedependence.jl rename to src/utility/timedependence.jl index 363b9eb2..208d39c5 100644 --- a/src/operators/timedependence.jl +++ b/src/utility/timedependence.jl @@ -11,5 +11,13 @@ istimed(x) = istimed(TimeDependence(x)) safe_eval(x, args...) = safe_eval(TimeDependence(x), x, args...) +# wrapper around _eval_at +function safe_eval(::TimeDependent, x) + throw(ArgumentError("attempting to evaluate time-dependent object without specifiying a time")) +end +function safe_eval(::NotTimeDependent, x, t::Number) + throw(ArgumentError("attempting to evaluate time-independent object at a time")) +end + # Internal use only, works always _eval_at(x, args...) = x # -> this is what you should define for custom structs diff --git a/test/algorithms.jl b/test/algorithms.jl index 4d0ca284..6e0490c5 100644 --- a/test/algorithms.jl +++ b/test/algorithms.jl @@ -287,6 +287,63 @@ end @test abs(dot(ψ₀, ψ_lazy)) ≈ 1 atol = atol end + + infinite_algs = [VUMPS(; tol, verbosity), + IDMRG1(; tol, verbosity), + GradientGrassmann(; tol=tol, verbosity=verbosity)] + finite_algs = [DMRG(; verbosity), + DMRG2(; verbosity, trscheme=truncdim(D)), + DMRG(; verbosity)] + window_algs = map((x, y) -> Window(x, y, x), infinite_algs, finite_algs) + + L = 10 + H = force_planar(transverse_field_ising(; g)) + H = Window(H, repeat(H, L), H) + + @testset "Window $i" for (i, alg) in enumerate(window_algs) + ψl = InfiniteMPS(rand, ComplexF64, [ℙ^2], [ℙ^D]) + ψr = InfiniteMPS(rand, ComplexF64, [ℙ^2], [ℙ^D]) + ψ₀ = WindowMPS(rand, ComplexF64, L, ℙ^2, ℙ^D, ψl, ψr) + + v₀ = variance(ψ₀, H) + ψ, envs, δ = find_groundstate(ψ₀, H, alg) + v = variance(ψ, H, envs) + + @test sum(δ) < 1e-3 + @test v₀ > v && v < 1e-2 # energy variance should be low + + #make sure that it is uniform + edens = expectation_value(ψ, H, envs) + @test edens[1] ≈ edens[2] atol = 1e-06 + @test edens[end - 1] ≈ edens[end] atol = 1e-06 + end + + g += 1 + H2 = force_planar(transverse_field_ising(; g)) + H2 = Window(H2, repeat(H2, L), H2) + Hlazy = LazySum([H, H2]) + + @testset "LazySum Window $i" for (i, alg) in enumerate(window_algs) + ψl = InfiniteMPS(rand, ComplexF64, [ℙ^2], [ℙ^D]) + ψr = InfiniteMPS(rand, ComplexF64, [ℙ^2], [ℙ^D]) + ψ₀ = WindowMPS(rand, ComplexF64, 10, ℙ^2, ℙ^D, ψl, ψr) + + v₀ = variance(ψ₀, Hlazy) + ψ, envs, δ = find_groundstate(ψ₀, Hlazy, alg) + v = variance(ψ, Hlazy) + + @test sum(δ) < 1e-3 + @test v₀ > v && v < 1e-2 # energy variance should be low + + ψ_nolazy, envs_nolazy, _ = find_groundstate(ψ₀, sum(Hlazy), alg) + @test expectation_value(ψ, Hlazy, envs) ≈ + expectation_value(ψ_nolazy, sum(Hlazy), envs_nolazy) atol = 1 - 06 + + #make sure that it is uniform + edens = expectation_value(ψ, Hlazy, envs) + @test edens[1] ≈ edens[2] atol = 1e-06 + @test edens[end - 1] ≈ edens[end] atol = 1e-06 + end end @testset "timestep" verbose = true begin @@ -311,6 +368,44 @@ end @test (3 + 1.55 - 0.1) * E₀ ≈ E atol = 1e-2 end + ψl₀ = InfiniteMPS(rand, ComplexF64, [ℙ^2], [ℙ^1]) + ψr₀ = InfiniteMPS(rand, ComplexF64, [ℙ^2], [ℙ^1]) + ψ₀ = WindowMPS(ψl₀, ψ₀, ψr₀; fixleft=true, fixright=true) + E₀ = expectation_value(ψ₀, H) + + @testset "Fixed Window $(alg isa TDVP ? "TDVP" : "TDVP2")" for alg in algs + ψ, envs = timestep(ψ₀, H, 0.0, dt, alg) + E = expectation_value(ψ, H, envs) + @test sum(E₀) ≈ sum(E) atol = 1e-2 + end + + Hlazy = LazySum([3 * H, 1.55 * H, -0.1 * H]) + + @testset "Fixed Window LazySum $(alg isa TDVP ? "TDVP" : "TDVP2")" for alg in algs + ψ, envs = timestep(ψ₀, Hlazy, 0.0, dt, alg) + E = expectation_value(ψ, Hlazy, envs) + @test (3 + 1.55 - 0.1) * sum(E₀) ≈ sum(E) atol = 5e-2 + end + + ψl₀ = InfiniteMPS(rand, ComplexF64, [ℙ^2], [ℙ^1]) + ψr₀ = InfiniteMPS(rand, ComplexF64, [ℙ^2], [ℙ^1]) + ψ₀ = WindowMPS(ψl₀, ψ₀, ψr₀; fixleft=true, fixright=true) + E₀ = expectation_value(ψ₀, H) + + @testset "Fixed Window $(alg isa TDVP ? "TDVP" : "TDVP2")" for alg in algs + ψ, envs = timestep(ψ₀, H, 0.0, dt, alg) + E = expectation_value(ψ, H, envs) + @test sum(E₀) ≈ sum(E) atol = 1e-2 + end + + Hlazy = LazySum([3 * H, 1.55 * H, -0.1 * H]) + + @testset "Fixed Window LazySum $(alg isa TDVP ? "TDVP" : "TDVP2")" for alg in algs + ψ, envs = timestep(ψ₀, Hlazy, 0.0, dt, alg) + E = expectation_value(ψ, Hlazy, envs) + @test (3 + 1.55 - 0.1) * sum(E₀) ≈ sum(E) atol = 5e-2 + end + Ht = MultipliedOperator(H, t -> 4) + MultipliedOperator(H, 1.45) @testset "Finite TimeDependent LazySum $(alg isa TDVP ? "TDVP" : "TDVP2")" for alg in @@ -351,6 +446,38 @@ end Et = expectation_value(ψt, Ht(1.0), envst) @test E ≈ Et atol = 1e-8 end + + ψl₀ = InfiniteMPS(rand, ComplexF64, [ℙ^2], [ℙ^10]) + ψr₀ = InfiniteMPS(rand, ComplexF64, [ℙ^2], [ℙ^10]) + ψ₀ = FiniteMPS(fill(TensorMap(rand, ComplexF64, ℙ^10 * ℙ^2, ℙ^10), 5)) + ψ₀ = WindowMPS(ψl₀, ψ₀, ψr₀; fixleft=false, fixright=false) + H = force_planar(transverse_field_ising(; g=3)) + Hwindow = Window(H, H, H) + ψ₀, envs, _ = find_groundstate(ψ₀, Hwindow) + E₀ = expectation_value(ψ₀, Hwindow, envs) + Hquench = force_planar(transverse_field_ising(; g=4)) + Hquench = Window(Hquench, Hquench, Hquench) + + @testset "Variable Window $(alg isa TDVP ? "TDVP" : "TDVP2")" for alg in algs + ψ, envs = timestep(ψ₀, Hquench, 0.0, dt, WindowTDVP(; middle=alg)) + E = expectation_value(ψ, Hquench, envs) + @test sum(E₀) ≈ sum(expectation_value(ψ, Hwindow)) atol = 0.1 + #test if still uniform + @test E[1] ≈ E[2] atol = 1e-04 + @test E[end - 1] ≈ E[end] atol = 1e-04 + end + + Hquench = LazySum([MultipliedOperator(Hwindow, 3.5), MultipliedOperator(Hquench, 1.7)]) + + @testset "Variable Window LazySum/MultipliedOperator $(alg isa TDVP ? "TDVP" : "TDVP2")" for alg in + algs + ψ, envs = timestep(ψ₀, Hquench, 0.0, dt, WindowTDVP(; middle=alg)) + E = expectation_value(ψ, Hquench, envs) + @test sum(E₀) ≈ sum(expectation_value(ψ, Hwindow)) atol = 0.1 + #test if still uniform + @test E[1] ≈ E[2] atol = 1e-03 + @test E[end - 1] ≈ E[end] atol = 1e-03 + end end @testset "time_evolve" verbose = true begin @@ -376,6 +503,14 @@ end E = expectation_value(ψ, H, envs) @test E₀ ≈ E atol = 1e-2 end + + #= + @testset "WindowMPS TDVP" begin + ψwindow, envs = timestep(Ψwindow₀, H, dt, TDVP()) + E = expectation_value(ψ, H, envs) + @test E₀ ≈ E atol = 1e-2 + end + =# end @testset "leading_boundary" verbose = true begin @@ -458,15 +593,13 @@ end # random nn interaction nn = TensorMap(rand, ComplexF64, pspace * pspace, pspace * pspace) nn += nn' - H = MPOHamiltonian(nn) - Δt = 0.1 - expH = make_time_mpo(H, Δt, WII()) - O = convert(DenseMPO, expH) - Op = periodic_boundary_conditions(O, 10) - Op′ = changebonds(Op, SvdCut(; trscheme=truncdim(5))) + mpo1 = periodic_boundary_conditions(convert(DenseMPO, + make_time_mpo(MPOHamiltonian(nn), 0.1, + WII())), 10) + mpo2 = changebonds(mpo1, SvdCut(; trscheme=truncdim(5))) - @test dim(space(Op′[5], 1)) < dim(space(Op[5], 1)) + @test dim(space(mpo2[5], 1)) < dim(space(mpo1[5], 1)) end @testset "infinite mps" begin @@ -559,7 +692,8 @@ end @testset "Dynamical DMRG" verbose = true begin ham = force_planar(-1.0 * transverse_field_ising(; g=-4.0)) gs, = find_groundstate(InfiniteMPS([ℙ^2], [ℙ^10]), ham, VUMPS(; verbosity=0)) - window = WindowMPS(gs, copy.([gs.AC[1]; [gs.AR[i] for i in 2:10]]), gs) + window = WindowMPS(gs, copy.([gs.AC[1]; [gs.AR[i] for i in 2:10]]), gs; fixleft=true, + fixright=true) szd = force_planar(S_z()) @test [expectation_value(gs, i => szd) for i in 1:length(window)] ≈ diff --git a/test/states.jl b/test/states.jl index 5bfac853..b4587b5b 100644 --- a/test/states.jl +++ b/test/states.jl @@ -96,16 +96,64 @@ end end end -@testset "WindowMPS" begin - g = 8.0 - ham = force_planar(transverse_field_ising(; g)) +@testset "Fixed WindowMPS" begin + ham = force_planar(transverse_field_ising(; g=8.0)) + (gs, _, _) = find_groundstate(InfiniteMPS([ℙ^2], [ℙ^10]), ham, VUMPS(; verbosity=0)) - # operator for testing expectation_value - X = S_x(; spin=1 // 2) - E = TensorMap(ComplexF64[1 0; 0 1], ℂ^2 ← ℂ^2) - O = force_planar(-(S_zz(; spin=1 // 2) + (g / 2) * (X ⊗ E + E ⊗ X))) + #constructor 1 - give it a plain array of tensors + window_1 = WindowMPS(gs, copy.([gs.AC[1]; [gs.AR[i] for i in 2:10]]), gs; fixleft=true, + fixright=true) - gs, = find_groundstate(InfiniteMPS([ℙ^2], [ℙ^10]), ham, VUMPS(; verbosity=0)) + #constructor 2 - used to take a "slice" from an infinite mps + window_2 = WindowMPS(gs, 10; fixleft=true, fixright=true) + + # we should logically have that window_1 approximates window_2 + ovl = dot(window_1, window_2) + @test ovl ≈ 1 atol = 1e-8 + + #constructor 3 - random initial tensors + window = WindowMPS(rand, ComplexF64, 10, ℙ^2, ℙ^10, gs, gs; fixleft=true, fixright=true) + normalize!(window) + + for i in 1:length(window) + @test window.AC[i] ≈ window.AL[i] * window.CR[i] + @test window.AC[i] ≈ MPSKit._transpose_front(window.CR[i - 1] * + MPSKit._transpose_tail(window.AR[i])) + end + + @test norm(window) ≈ 1 + window = window * 3 + @test 9 ≈ norm(window)^2 + window = 3 * window + @test 9 * 9 ≈ norm(window)^2 + normalize!(window) + + edens1 = expectation_value(window, ham) + e1 = expectation_value(window, ham, 1:length(window)) + + v1 = variance(window, ham) + (window, envs, _) = find_groundstate(window, ham, DMRG(; verbosity=0)) + v2 = variance(window, ham) + + edens2 = expectation_value(window, ham) + e2 = expectation_value(window, ham, 1:length(window)) + + @test v2 < v1 + @test real(e2) ≤ real(e1) + + (window, envs) = timestep(window, ham, 0.1, 0.0, TDVP2(), envs) + (window, envs) = timestep(window, ham, 0.1, 0.0, TDVP(), envs) + + edens3 = expectation_value(window, ham) + e3 = expectation_value(window, ham, 1:length(window), envs) + + @test edens2 ≈ edens3 atol = 1e-4 + @test e2 ≈ e3 atol = 1e-4 +end + +@testset "Variable WindowMPS" begin + ham = force_planar(transverse_field_ising(; g=8.0)) + (gs, _, _) = find_groundstate(InfiniteMPS([ℙ^2], [ℙ^10]), ham, VUMPS(; verbosity=0)) # constructor 1 - give it a plain array of tensors window_1 = WindowMPS(gs, copy.([gs.AC[1]; [gs.AR[i] for i in 2:10]]), gs) @@ -123,8 +171,8 @@ end for i in 1:length(window) @test window.AC[i] ≈ window.AL[i] * window.CR[i] - @test window.AC[i] ≈ - _transpose_front(window.CR[i - 1] * _transpose_tail(window.AR[i])) + @test window.AC[i] ≈ MPSKit._transpose_front(window.CR[i - 1] * + MPSKit._transpose_tail(window.AR[i])) end @test norm(window) ≈ 1 @@ -136,17 +184,21 @@ end e1 = expectation_value(window, (1, 2) => O) - window, envs, _ = find_groundstate(window, ham, DMRG(; verbosity=0)) + v1 = variance(window, ham) + gs_alg = Window(VUMPS(; verbosity=0), DMRG(; verbosity=0), VUMPS(; verbosity=0)) + (window, envs, _) = find_groundstate(window, ham, gs_alg) + v2 = variance(window, ham) e2 = expectation_value(window, (1, 2) => O) @test real(e2) ≤ real(e1) - window, envs = timestep(window, ham, 0.1, 0.0, TDVP2(), envs) - window, envs = timestep(window, ham, 0.1, 0.0, TDVP(), envs) + (window, envs) = timestep(window, ham, 0.1, 0.0, WindowTDVP(; middle=TDVP2()), envs) + (window, envs) = timestep(window, ham, 0.1, 0.0, WindowTDVP(), envs) e3 = expectation_value(window, (1, 2) => O) + @test edens2 ≈ edens3 atol = 1e-4 @test e2 ≈ e3 atol = 1e-4 end @@ -181,7 +233,7 @@ end @testset "Infinite" for (th, D, d) in [(force_planar(transverse_field_ising()), ℙ^10, ℙ^2), - (heisenberg_XXX(SU2Irrep; spin=1), Rep[SU₂](1 => 3, 0 => 2), + (heisenberg_XXX(SU2Irrep; spin=1), Rep[SU₂](1 => 1, 0 => 3), Rep[SU₂](1 => 1))] period = rand(1:4) ψ = InfiniteMPS(fill(d, period), fill(D, period))