From bc6a4662f96fc1090e14adee87fa1e4d3938eb0c Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Thu, 23 Nov 2023 16:30:41 +0100 Subject: [PATCH 01/62] window stuff --- src/algorithms/expval.jl | 185 +++++++++++-------- src/algorithms/timestep/windowtdvp.jl | 246 ++++++++++++++++++++++++++ src/states/window.jl | 6 +- src/states/windowmps.jl | 5 + 4 files changed, 361 insertions(+), 81 deletions(-) create mode 100644 src/algorithms/timestep/windowtdvp.jl diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl index 2cea5c151..330e98798 100644 --- a/src/algorithms/expval.jl +++ b/src/algorithms/expval.jl @@ -1,13 +1,30 @@ -#works for general tensors +""" + expectation_value(ψ, O, [location], [environments]) + +Compute the expectation value of an operator `O` on a state `ψ`. If `location` is given, the +operator is applied at that location. If `environments` is given, the expectation value +might be computed more efficiently by re-using previously calculated environments. + +# Arguments +* `ψ::AbstractMPS` : the state on which to compute the expectation value +* `O` : the operator to compute the expectation value of. This can either be an `AbstractMPO`, a single `AbstractTensorMap` or an array of `AbstractTensorMap`s. +* `location::Union{Int,AbstractRange{Int}}` : the location at which to apply the operator. Only applicable for operators that act on a subset of all sites. +* `environments::Cache` : the environments to use for the calculation. If not given, they will be calculated. +""" +function expectation_value end + +# Local operators +# --------------- function expectation_value( - state::Union{InfiniteMPS,WindowMPS,FiniteMPS}, opp::AbstractTensorMap -) - return expectation_value(state, fill(opp, length(state))) + ψ::Union{InfiniteMPS,WindowMPS,FiniteMPS}, O::AbstractTensorMap{S,1,1} +) where {S} + return expectation_value(ψ, fill(O, length(ψ))) end function expectation_value( - state::Union{InfiniteMPS,WindowMPS,FiniteMPS}, opps::AbstractArray{<:AbstractTensorMap} -) - map(zip(state.AC, opps)) do (ac, opp) + ψ::Union{InfiniteMPS,WindowMPS,FiniteMPS}, + opps::AbstractArray{<:AbstractTensorMap{S,1,1}}, +) where {S} + return map(zip(ψ.AC, opps)) do (ac, opp) tr( ac' * transpose( opp * transpose( @@ -22,86 +39,78 @@ function expectation_value( end end -""" -calculates the expectation value of op, where op is a plain tensormap where the first index works on site at -""" +# Multi-site operators +# -------------------- +# TODO: replace Vector{MPOTensor} with FiniteMPO function expectation_value( - state::Union{FiniteMPS{T},WindowMPS{T},InfiniteMPS{T}}, op::AbstractTensorMap, at::Int -) where {T<:MPSTensor} - return expectation_value(state, decompose_localmpo(add_util_leg(op)), at) + ψ::Union{FiniteMPS{T},WindowMPS{T},InfiniteMPS{T}}, O::AbstractTensorMap{S,N,N}, at::Int +) where {S,N,T<:MPSTensor{S}} + return expectation_value(ψ, decompose_localmpo(add_util_leg(O)), at) end - -""" -calculates the expectation value of op = op1*op2*op3*... (ie an N site operator) starting at site at -""" function expectation_value( - state::Union{FiniteMPS{T},WindowMPS{T},InfiniteMPS{T}}, - op::AbstractArray{<:AbstractTensorMap}, + ψ::Union{FiniteMPS{T},WindowMPS{T},InfiniteMPS{T}}, + O::AbstractArray{<:MPOTensor{S}}, at::Int, -) where {T<:MPSTensor} - firstspace = _firstspace(first(op)) - (firstspace == oneunit(firstspace) && _lastspace(last(op)) == firstspace') || throw( +) where {S,T<:MPSTensor{S}} + firstspace = _firstspace(first(O)) + (firstspace == oneunit(firstspace) && _lastspace(last(O)) == firstspace') || throw( ArgumentError( "localmpo should start and end in a trivial leg, not with $(firstspace)" ), ) - ut = fill_data!(similar(op[1], firstspace), one) + ut = fill_data!(similar(O[1], firstspace), one) @plansor v[-1 -2; -3] := isomorphism( - storagetype(T), - left_virtualspace(state, at - 1), - left_virtualspace(state, at - 1), + storagetype(T), left_virtualspace(ψ, at - 1), left_virtualspace(ψ, at - 1) )[ -1 -3 ] * conj(ut[-2]) tmp = - v * TransferMatrix( - state.AL[at:(at + length(op) - 1)], op, state.AL[at:(at + length(op) - 1)] - ) + v * TransferMatrix(ψ.AL[at:(at + length(O) - 1)], O, ψ.AL[at:(at + length(O) - 1)]) return @plansor tmp[1 2; 3] * ut[2] * - state.CR[at + length(op) - 1][3; 4] * - conj(state.CR[at + length(op) - 1][1; 4]) + ψ.CR[at + length(O) - 1][3; 4] * + conj(ψ.CR[at + length(O) - 1][1; 4]) end -""" - calculates the expectation value for the given operator/hamiltonian -""" -expectation_value(state, envs::Cache) = expectation_value(state, envs.opp, envs); -function expectation_value(state, ham::MPOHamiltonian) - return expectation_value(state, ham, environments(state, ham)) +# MPOHamiltonian +# -------------- +# TODO: add section in documentation to explain convention +expectation_value(ψ, envs::Cache) = expectation_value(ψ, envs.opp, envs) +function expectation_value(ψ, ham::MPOHamiltonian) + return expectation_value(ψ, ham, environments(ψ, ham)) end -function expectation_value(state::WindowMPS, ham::MPOHamiltonian, envs::FinEnv) - vals = expectation_value_fimpl(state, ham, envs) + +#under review +function expectation_value(ψ::WindowMPS, ham::MPOHamiltonian, envs::FinEnv) + vals = expectation_value_fimpl(ψ, ham, envs) tot = 0.0 + 0im for i in 1:(ham.odim), j in 1:(ham.odim) - tot += @plansor leftenv(envs, length(state), state)[i][1 2; 3] * - state.AC[end][3 4; 5] * - rightenv(envs, length(state), state)[j][5 6; 7] * - ham[length(state)][i, j][2 8; 4 6] * - conj(state.AC[end][1 8; 7]) + tot += @plansor leftenv(envs, length(ψ), ψ)[i][1 2; 3] * + ψ.AC[end][3 4; 5] * + rightenv(envs, length(ψ), ψ)[j][5 6; 7] * + ham[length(ψ)][i, j][2 8; 4 6] * + conj(ψ.AC[end][1 8; 7]) end - return vals, tot / (norm(state.AC[end])^2) + return vals, tot / (norm(ψ.AC[end])^2) end -function expectation_value(state::FiniteMPS, ham::MPOHamiltonian, envs::FinEnv) - return expectation_value_fimpl(state, ham, envs) +function expectation_value(ψ::FiniteMPS, ham::MPOHamiltonian, envs::FinEnv) + return expectation_value_fimpl(ψ, ham, envs) end -function expectation_value_fimpl( - state::AbstractFiniteMPS, ham::MPOHamiltonian, envs::FinEnv -) - ens = zeros(scalartype(state), length(state)) - for i in 1:length(state), (j, k) in keys(ham[i]) +function expectation_value_fimpl(ψ::AbstractFiniteMPS, ham::MPOHamiltonian, envs::FinEnv) + ens = zeros(scalartype(ψ), length(ψ)) + for i in 1:length(ψ), (j, k) in keys(ham[i]) !((j == 1 && k != 1) || (k == ham.odim && j != ham.odim)) && continue - cur = @plansor leftenv(envs, i, state)[j][1 2; 3] * - state.AC[i][3 7; 5] * - rightenv(envs, i, state)[k][5 8; 6] * - conj(state.AC[i][1 4; 6]) * + cur = @plansor leftenv(envs, i, ψ)[j][1 2; 3] * + ψ.AC[i][3 7; 5] * + rightenv(envs, i, ψ)[k][5 8; 6] * + conj(ψ.AC[i][1 4; 6]) * ham[i][j, k][2 4; 7 8] if !(j == 1 && k == ham.odim) cur /= 2 @@ -110,7 +119,7 @@ function expectation_value_fimpl( ens[i] += cur end - n = norm(state.AC[end])^2 + n = norm(ψ.AC[end])^2 return ens ./ n end @@ -164,28 +173,50 @@ function expectation_value( return tot end -function expectation_value(st::InfiniteMPS, mpo::DenseMPO) - return expectation_value(convert(MPSMultiline, st), convert(MPOMultiline, mpo)) -end; -function expectation_value(st::MPSMultiline, mpo::MPOMultiline) - return expectation_value(st, environments(st, mpo)) -end; -function expectation_value(st::InfiniteMPS, ca::PerMPOInfEnv) - return expectation_value(convert(MPSMultiline, st), ca) -end; -function expectation_value(st::MPSMultiline, opp::MPOMultiline, ca::PerMPOInfEnv) - retval = PeriodicArray{scalartype(st.AC[1, 1]),2}(undef, size(st, 1), size(st, 2)) - for (i, j) in product(1:size(st, 1), 1:size(st, 2)) - retval[i, j] = @plansor leftenv(ca, i, j, st)[1 2; 3] * - opp[i, j][2 4; 6 5] * - st.AC[i, j][3 6; 7] * - rightenv(ca, i, j, st)[7 5; 8] * - conj(st.AC[i + 1, j][1 4; 8]) +# Transfer matrices +# ----------------- +function expectation_value(ψ::InfiniteMPS, mpo::DenseMPO) + return expectation_value(convert(MPSMultiline, ψ), convert(MPOMultiline, mpo)) +end +function expectation_value(ψ::MPSMultiline, mpo::MPOMultiline) + return expectation_value(ψ, environments(ψ, mpo)) +end +function expectation_value(ψ::InfiniteMPS, ca::PerMPOInfEnv) + return expectation_value(convert(MPSMultiline, ψ), ca) +end +function expectation_value(ψ::MPSMultiline, O::MPOMultiline, ca::PerMPOInfEnv) + retval = PeriodicMatrix{scalartype(ψ)}(undef, size(ψ, 1), size(ψ, 2)) + for (i, j) in product(1:size(ψ, 1), 1:size(ψ, 2)) + retval[i, j] = @plansor leftenv(ca, i, j, ψ)[1 2; 3] * + O[i, j][2 4; 6 5] * + ψ.AC[i, j][3 6; 7] * + rightenv(ca, i, j, ψ)[7 5; 8] * + conj(ψ.AC[i + 1, j][1 4; 8]) end return retval end -expectation_value(state::FiniteQP, opp) = expectation_value(convert(FiniteMPS, state), opp) -function expectation_value(state::FiniteQP, opp::MPOHamiltonian) - return expectation_value(convert(FiniteMPS, state), opp) +expectation_value(ψ::FiniteQP, O) = expectation_value(convert(FiniteMPS, ψ), O) +function expectation_value(ψ::FiniteQP, O::MPOHamiltonian) + return expectation_value(convert(FiniteMPS, ψ), O) +end + +# define expectation_value for Window +function expectation_value(Ψ::WindowMPS, windowOp::Window, at::Int) + if at < 1 + return expectation_value(Ψ.left_gs, windowOp.left, at) + elseif 1 <= at <= length(Ψ.window) + return expectation_value(Ψ, windowOp.middle, at) + else + return expectation_value(Ψ.right_gs, windowOp.right, at) + end +end + +function expectation_value( + Ψ::WindowMPS, windowH::Window, windowEnvs::Window{C,D,C}=environments(Ψ, windowH) +) where {C<:Union{MultipleEnvironments,Cache},D<:Union{MultipleEnvironments,Cache}} + left = expectation_value(Ψ.left_gs, windowH.left, windowEnvs.left) + middle = expectation_value(Ψ.window, windowH.middle, windowEnvs.middle) + right = expectation_value(Ψ.right_gs, windowH.right, windowEnvs.right) + return [left.data..., middle..., right.data...] end diff --git a/src/algorithms/timestep/windowtdvp.jl b/src/algorithms/timestep/windowtdvp.jl new file mode 100644 index 000000000..9b1557bbe --- /dev/null +++ b/src/algorithms/timestep/windowtdvp.jl @@ -0,0 +1,246 @@ + +function _update_leftEnv!( + nleft::InfiniteMPS, WindowEnv::Window{E,F,E} +) where {E<:Cache,F<:Cache} + l = leftenv(WindowEnv.left, 1, nleft) + WindowEnv.middle.ldependencies[:] = similar.(WindowEnv.middle.ldependencies) # forget the old left dependencies - this forces recalculation whenever leftenv is called + WindowEnv.middle.leftenvs[1] = l + + return Window(WindowEnv.left, WindowEnv.middle, WindowEnv.right) +end + +function _update_leftEnv!( + nleft::InfiniteMPS, WindowEnv::Window{E,F,E} +) where {E<:MultipleEnvironments,F<:MultipleEnvironments} + @assert length(WindowEnv.middle) == length(WindowEnv.left) + + for (subEnvLeft, subEnvMiddle) in zip(WindowEnv.left, WindowEnv.middle) + l = leftenv(subEnvLeft, 1, nleft) + subEnvMiddle.ldependencies[:] = similar.(subEnvMiddle.ldependencies) # forget the old left dependencies - this forces recalculation whenever leftenv is called + subEnvMiddle.leftenvs[1] = l + end + return Window(WindowEnv.left, WindowEnv.middle, WindowEnv.right) +end + +function _update_rightEnv!( + nright::InfiniteMPS, WindowEnv::Window{E,F,E} +) where {E<:Cache,F<:Cache} + r = rightenv(WindowEnv.right, length(nright), nright) + WindowEnv.middle.rdependencies[:] = similar.(WindowEnv.middle.rdependencies) # forget the old right dependencies - this forces recalculation + WindowEnv.middle.rightenvs[end] = r + + return Window(WindowEnv.left, WindowEnv.middle, WindowEnv.right) +end + +function _update_rightEnv!( + nright::InfiniteMPS, WindowEnv::Window{E,F,E} +) where {E<:MultipleEnvironments,F<:MultipleEnvironments} + @assert length(WindowEnv.middle) == length(WindowEnv.right) + + for (subEnvMiddle, subEnvRight) in zip(WindowEnv.middle, WindowEnv.right) + r = rightenv(subEnvRight, length(nright), nright) + subEnvMiddle.rdependencies[:] = similar.(subEnvMiddle.rdependencies) # forget the old right dependencies - this forces recalculation + subEnvMiddle.rightenvs[end] = r + end + return Window(WindowEnv.left, WindowEnv.middle, WindowEnv.right) +end + +function leftexpand( + st::WindowMPS, + H::Union{<:MultipliedOperator,<:SumOfOperators}, + Envs; + singval=1e-2, + growspeed=10, +) + (U, S, V) = tsvd(st.left_gs.CR[1]; alg=TensorKit.SVD()) + + if minimum([minimum(diag(v)) for (k, v) in blocks(S)]) > singval + (nst, _) = changebonds( + st.left_gs, H, OptimalExpand(; trscheme=truncbelow(singval, growspeed)), Envs + ) + + # the AL-bond dimension changed, and therefore our window also needs updating + v = TensorMap( + rand, ComplexF64, left_virtualspace(nst, 0), left_virtualspace(st.left_gs, 0) + ) + (vals, vecs) = eigsolve( + flip(TransferMatrix(st.left_gs.AL, nst.AL)), v, 1, :LM, Arnoldi() + ) + rho = pinv(nst.CR[0]) * vecs[1] * st.left_gs.CR[0] #CR[0] == CL[1] + st.AC[1] = _transpose_front(normalize!(rho * st.CR[0]) * _transpose_tail(st.AR[1])) + + recalculate!(Envs, nst) #updates left infinite env based on expanded state + return nst, Envs + end + return st.left_gs, Envs +end +function leftexpand( + st::WindowMPS, H::Union{MPOHamiltonian,DenseMPO,SparseMPO}, t::Number, Envs; kwargs... +) + return leftexpand(st, UntimedOperator(H), t, Envs; kwargs...) +end + +function rightexpand( + st::WindowMPS, + H::Union{<:MultipliedOperator,<:SumOfOperators}, + Envs; + singval=1e-2, + growspeed=10, +) + (U, S, V) = tsvd(st.left_gs.CR[1]; alg=TensorKit.SVD()) + + if minimum([minimum(diag(v)) for (k, v) in blocks(S)]) > singval + (nst, _) = changebonds( + st.right_gs, + H, + OptimalExpand(; trscheme=truncbelow(singval, growspeed)), + Envs, + ) + + v = TensorMap( + rand, ComplexF64, right_virtualspace(st.right_gs, 0), right_virtualspace(nst, 0) + ) + (vals, vecs) = eigsolve( + TransferMatrix(st.right_gs.AR, nst.AR), v, 1, :LM, Arnoldi() + ) + rho = st.right_gs.CR[0] * vecs[1] * pinv(nst.CR[0]) + st.AC[end] = st.AL[end] * normalize!(st.CR[end] * rho) + + recalculate!(Envs, nst) #updates right infinite env based on expanded state + return nst, Envs + end + return st.right_gs, Envs +end +function rightexpand( + st::WindowMPS, H::Union{MPOHamiltonian,DenseMPO,SparseMPO}, t::Number, Envs; kwargs... +) + return rightexpand(st, UntimedOperator(H), t, Envs; kwargs...) +end + +function timestep!( + Ψ::WindowMPS, + H::Window, + t::Number, + dt::Number, + alg::TDVP, + env::Window=environments(Ψ, H); + leftevolve=true, + rightevolve=true, +) + #first evolve left state + if leftevolve + nleft, _ = timestep(Ψ.left_gs, H.left, t, dt, alg, env.left; leftorthflag=true) #env gets updated in place + _update_leftEnv!(nleft, env) + else + nleft = Ψ.left_gs + end + + # left to right sweep on window + for i in 1:(length(Ψ) - 1) + h_ac = ∂∂AC(i, Ψ, H.middle, env.middle) + Ψ.AC[i] = integrate(h_ac, Ψ.AC[i], t, dt / 2, alg.integrator) + + h_c = ∂∂C(i, Ψ, H.middle, env.middle) + Ψ.CR[i] = integrate(h_c, Ψ.CR[i], t, -dt / 2, alg.integrator) + end + + h_ac = ∂∂AC(length(Ψ), Ψ, H.middle, env.middle) + Ψ.AC[end] = integrate(h_ac, Ψ.AC[end], t, dt / 2, alg.integrator) + + if rightevolve + nright, _ = timestep(Ψ.right_gs, H.right, t + dt, dt, alg, env.right; leftorthflag=false) # env gets updated in place + _update_rightEnv!(nright, env) + else + nright = Ψ.right_gs + end + + # right to left sweep on window + for i in length(Ψ):-1:2 + h_ac = ∂∂AC(i, Ψ, H.middle, env.middle) + Ψ.AC[i] = integrate(h_ac, Ψ.AC[i], t + dt / 2, dt / 2, alg.integrator) + + h_c = ∂∂C(i - 1, Ψ, H.middle, env.middle) + Ψ.CR[i - 1] = integrate(h_c, Ψ.CR[i - 1], t + dt / 2, -dt / 2, alg.integrator) + end + + h_ac = ∂∂AC(1, Ψ, H.middle, env.middle) + Ψ.AC[1] = integrate(h_ac, Ψ.AC[1], t, dt / 2, alg.integrator) + + return WindowMPS(nleft, Ψ.window, nright), env +end + +function timestep!( + Ψ::WindowMPS, + H::Window, + t::Number, + dt::Number, + alg::TDVP2, + env::Window=environments(Ψ, H); + leftevolve=false, + rightevolve=false, + kwargs..., +) + singleTDVPalg = TDVP(; + integrator=alg.integrator, tolgauge=alg.tolgauge, maxiter=alg.maxiter + ) + + # first evolve left state + if leftevolve + # expand the bond dimension using changebonds + nleft, _ = leftexpand(Ψ, H.left(t), env.left; kwargs...) + # fill it by doing regular TDVP + nleft, _ = timestep( + nleft, H.left, t, dt, singleTDVPalg, env.left; leftorthflag=true + ) + _update_leftEnv!(nleft, env) + else + nleft = Ψ.left_gs + end + + # left to right sweep on window + for i in 1:(length(Ψ) - 1) + h_ac2 = ∂∂AC2(i, Ψ, H.middle, env.middle) + ac2 = Ψ.AC[i] * _transpose_tail(Ψ.AR[i + 1]) + ac2 = integrate(h_ac2, ac2, t, dt / 2, alg.integrator) + + U, S, V, = tsvd!(ac2; alg=TensorKit.SVD(), trunc=alg.trscheme) + + Ψ.AC[i] = (U, S) + Ψ.AC[i + 1] = (S, _transpose_front(V)) + + if i < length(Ψ) - 1 + h_ac = ∂∂AC(i + 1, Ψ, H.middle, env.middle) + Ψ.AC[i + 1] = integrate(h_ac, Ψ.AC[i + 1], t, -dt / 2, alg.integrator) + end + end + + if rightevolve + # expand the bond dimension using changebonds + nright, _ = rightexpand(Ψ, H.right(t), env.right; kwargs...) + # fill it by doing regular TDVP + nright, _ = timestep( + nright, H.right, t + dt, dt, singleTDVPalg, env.right; leftorthflag=false + ) #env gets updated in place + _update_rightEnv!(nright, env) + else + nright = Ψ.right_gs + end + + # right to left sweep on window + for i in length(Ψ):-1:2 + h_ac2 = ∂∂AC2(i - 1, Ψ, H.middle, env.middle) + ac2 = Ψ.AL[i - 1] * _transpose_tail(Ψ.AC[i]) + ac2 = integrate(h_ac2, ac2, t + dt / 2, dt / 2, alg.integrator) + U, S, V, = tsvd!(ac2; alg=TensorKit.SVD(), trunc=alg.trscheme) + + Ψ.AC[i - 1] = (U, S) + Ψ.AC[i] = (S, _transpose_front(V)) + + if i > 2 + h_ac = ∂∂AC(i - 1, Ψ, H.middle, env.middle) + Ψ.AC[i - 1] = integrate(h_ac, Ψ.AC[i - 1], t + dt / 2, -dt / 2, alg.integrator) + end + end + + return WindowMPS(nleft, Ψ.window, nright), env +end diff --git a/src/states/window.jl b/src/states/window.jl index 8e0b89546..415bdb453 100644 --- a/src/states/window.jl +++ b/src/states/window.jl @@ -3,7 +3,7 @@ " Window(leftstate,window,rightstate) - general struct an object with a left, middle and right part. + general struct of an object with a left, middle and right part. " struct Window{L,M,R} left::L @@ -11,7 +11,5 @@ struct Window{L,M,R} right::R end -Base.length(win::Window) = length(win.middle) - # do we need copy? -# Base.copy(win::Window) = Window(copy(win.left),copy(win.middle),copy(win.right)) +# Base.copy(win::Window) = Window(copy(win.left),copy(win.middle),copy(win.right)) \ No newline at end of file diff --git a/src/states/windowmps.jl b/src/states/windowmps.jl index a0403ed94..5b46d8fdb 100644 --- a/src/states/windowmps.jl +++ b/src/states/windowmps.jl @@ -53,6 +53,11 @@ end #=========================================================================================== Constructors ===========================================================================================# +function WindowMPS_copied( + Ψₗ::InfiniteMPS{A,B}, Ψₘ::FiniteMPS{A,B}, Ψᵣ::InfiniteMPS{A,B}=Ψₗ +) where {A<:GenericMPSTensor,B<:MPSBondTensor} + return WindowMPS{A,B}(copy(Ψₗ), Ψₘ, copy(Ψᵣ)) +end function WindowMPS( Ψₗ::InfiniteMPS, site_tensors::AbstractVector{<:GenericMPSTensor}, Ψᵣ::InfiniteMPS=Ψₗ From 6b9359f7ed52ce05e14726db49f8fc1c29ecac65 Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Sat, 25 Nov 2023 13:27:01 +0100 Subject: [PATCH 02/62] first test added --- src/algorithms/expval.jl | 2 +- test/algorithms.jl | 10 ++++++++++ 2 files changed, 11 insertions(+), 1 deletion(-) diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl index 330e98798..94f1b53b5 100644 --- a/src/algorithms/expval.jl +++ b/src/algorithms/expval.jl @@ -214,7 +214,7 @@ end function expectation_value( Ψ::WindowMPS, windowH::Window, windowEnvs::Window{C,D,C}=environments(Ψ, windowH) -) where {C<:Union{MultipleEnvironments,Cache},D<:Union{MultipleEnvironments,Cache}} +) where {C<:Cache,D<:Cache} left = expectation_value(Ψ.left_gs, windowH.left, windowEnvs.left) middle = expectation_value(Ψ.window, windowH.middle, windowEnvs.middle) right = expectation_value(Ψ.right_gs, windowH.right, windowEnvs.right) diff --git a/test/algorithms.jl b/test/algorithms.jl index 2c79badcd..9fa1ab496 100644 --- a/test/algorithms.jl +++ b/test/algorithms.jl @@ -75,6 +75,16 @@ end E = expectation_value(ψ, H, envs) @test sum(E₀) ≈ sum(E) atol = 1e-2 end + + Ψwindow₀ = WindowMPS(ψ₀, 10) + Hwindow = Window(H, repeat(H, 10), H) + E₀ = expectation_value(Ψwindow₀, Hwindow) + + @testset "WindowMPS TDVP" begin + ψwindow, envs = timestep(Ψwindow₀, H, dt, TDVP()) + E = expectation_value(ψ, H, envs) + @test sum(E₀) ≈ sum(E) atol = 1e-2 + end end @testset "leading_boundary" verbose = true begin From aed4a2c60639c6c75ff1356cd7f8c96d8b50e35d Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Sat, 25 Nov 2023 13:39:20 +0100 Subject: [PATCH 03/62] untested update to windowtdvp --- src/algorithms/timestep/windowtdvp.jl | 41 ++------------------------- 1 file changed, 3 insertions(+), 38 deletions(-) diff --git a/src/algorithms/timestep/windowtdvp.jl b/src/algorithms/timestep/windowtdvp.jl index 9b1557bbe..91d289665 100644 --- a/src/algorithms/timestep/windowtdvp.jl +++ b/src/algorithms/timestep/windowtdvp.jl @@ -9,19 +9,6 @@ function _update_leftEnv!( return Window(WindowEnv.left, WindowEnv.middle, WindowEnv.right) end -function _update_leftEnv!( - nleft::InfiniteMPS, WindowEnv::Window{E,F,E} -) where {E<:MultipleEnvironments,F<:MultipleEnvironments} - @assert length(WindowEnv.middle) == length(WindowEnv.left) - - for (subEnvLeft, subEnvMiddle) in zip(WindowEnv.left, WindowEnv.middle) - l = leftenv(subEnvLeft, 1, nleft) - subEnvMiddle.ldependencies[:] = similar.(subEnvMiddle.ldependencies) # forget the old left dependencies - this forces recalculation whenever leftenv is called - subEnvMiddle.leftenvs[1] = l - end - return Window(WindowEnv.left, WindowEnv.middle, WindowEnv.right) -end - function _update_rightEnv!( nright::InfiniteMPS, WindowEnv::Window{E,F,E} ) where {E<:Cache,F<:Cache} @@ -32,22 +19,10 @@ function _update_rightEnv!( return Window(WindowEnv.left, WindowEnv.middle, WindowEnv.right) end -function _update_rightEnv!( - nright::InfiniteMPS, WindowEnv::Window{E,F,E} -) where {E<:MultipleEnvironments,F<:MultipleEnvironments} - @assert length(WindowEnv.middle) == length(WindowEnv.right) - - for (subEnvMiddle, subEnvRight) in zip(WindowEnv.middle, WindowEnv.right) - r = rightenv(subEnvRight, length(nright), nright) - subEnvMiddle.rdependencies[:] = similar.(subEnvMiddle.rdependencies) # forget the old right dependencies - this forces recalculation - subEnvMiddle.rightenvs[end] = r - end - return Window(WindowEnv.left, WindowEnv.middle, WindowEnv.right) -end - +#Note: this needs to be tested function leftexpand( st::WindowMPS, - H::Union{<:MultipliedOperator,<:SumOfOperators}, + H, Envs; singval=1e-2, growspeed=10, @@ -74,15 +49,10 @@ function leftexpand( end return st.left_gs, Envs end -function leftexpand( - st::WindowMPS, H::Union{MPOHamiltonian,DenseMPO,SparseMPO}, t::Number, Envs; kwargs... -) - return leftexpand(st, UntimedOperator(H), t, Envs; kwargs...) -end function rightexpand( st::WindowMPS, - H::Union{<:MultipliedOperator,<:SumOfOperators}, + H, Envs; singval=1e-2, growspeed=10, @@ -111,11 +81,6 @@ function rightexpand( end return st.right_gs, Envs end -function rightexpand( - st::WindowMPS, H::Union{MPOHamiltonian,DenseMPO,SparseMPO}, t::Number, Envs; kwargs... -) - return rightexpand(st, UntimedOperator(H), t, Envs; kwargs...) -end function timestep!( Ψ::WindowMPS, From 78d439277d2a769a306634f2f74f5ba17b5d294b Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Tue, 13 Feb 2024 17:05:42 +0100 Subject: [PATCH 04/62] first version of new windowenv --- src/environments/FinEnv.jl | 29 +++++--------- src/environments/windowenv.jl | 65 +++++++++++++++++++++++++++++++ src/states/windowmps.jl | 73 ++++++++++++++++++++++++----------- 3 files changed, 125 insertions(+), 42 deletions(-) create mode 100644 src/environments/windowenv.jl diff --git a/src/environments/FinEnv.jl b/src/environments/FinEnv.jl index b45eaa597..e60a282b3 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,17 +95,18 @@ function environments(state::Union{FiniteMPS,WindowMPS}, opp::ProjectionOperator end #notify the cache that we updated in-place, so it should invalidate the dependencies +# this forces the transfers to be recalculated lazily function poison!(ca::FinEnv, ind) ca.ldependencies[ind] = similar(ca.ldependencies[ind]) - return ca.rdependencies[ind] = similar(ca.rdependencies[ind]) + ca.rdependencies[ind] = similar(ca.rdependencies[ind]) 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/windowenv.jl b/src/environments/windowenv.jl new file mode 100644 index 000000000..82829ac3f --- /dev/null +++ b/src/environments/windowenv.jl @@ -0,0 +1,65 @@ +""" + 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} <: Cache + left_env::A + fin_env::B + right_env::C +end + +#automatically construct the correct leftstart/rightstart for a WindowMPS +function WindowEnv(state::WindowMPS, O::Union{SparseMPO,MPOHamiltonian,DenseMPO}, above=nothing; lenvs=environments(state.left_gs, O),renvs=environments(state.right_gs, O)) + fin_env = environments(state, O, above, copy(leftenv(lenvs, 1, state.left_gs)), + copy(rightenv(renvs, length(state), state.right_gs))) + return WindowEnv(lenvs,fin_env,renvs) +end + +function environments(below::WindowMPS, above::WindowMPS) + above.left_gs == below.left_gs || throw(ArgumentError("left gs differs")) + above.right_gs == below.right_gs || throw(ArgumentError("right gs differs")) + + opp = fill(nothing, length(below)) + return environments(below, opp, above, l_LL(above), r_RR(above)) +end + +#notify the cache that we updated in-place, so it should invalidate the dependencies +poison!(ca::WindowEnv, ind) = poison!(ca.fin_env,ind) + +# check and recalculate the left and right start +function check_rightinfenv!(ca::WindowEnv) + if ca.right_env.rw[:,0] != ca.fin_env.rightstart[end]# !== doesn't work as intended + + poison!(ca, lastindex(ca.fin_env.rightstart)) #forces transfers to be recalculated lazily + + ca.fin_env.rightstart = ca.right_env.rw[:,0] #do some other checks and recalcs for bonddimensions? + + end +end + +function check_leftinfenv!(ca::WindowEnv) + if ca.left_env.lw[:,end+1] != ca.fin_env.leftstart[1]# !== doesn't work as intended + + poison!(ca, firstindex(ca.fin_env.leftstart)) #forces transfers to be recalculated lazily + + ca.fin_env.leftstart = ca.left_env.rw[:,end+1] #do some other checks and recalcs for bonddimensions? + + end +end + +#rightenv[ind] will be contracteable with the tensor on site [ind] +function rightenv(ca::WindowEnv, ind, state) + check_rightinfenv!(ca) + return rightenv(ca.fin_env, ind, state) +end + +function leftenv(ca::FinEnv, ind, state) + check_leftinfenv!(ca) + return leftenv(ca.fin_env, ind, state) +end + +function fix_infinite(ψ::WindowMPS,env::WindowEnv) + newenv = FinEnv(env.above,env.opp,env.ldependencies,env.rdependencies,copy(env.leftenvs),copy(env.rightenvs)) + return fix_infinite(ψ),newenv +end diff --git a/src/states/windowmps.jl b/src/states/windowmps.jl index 831da9a12..e77411c68 100644 --- a/src/states/windowmps.jl +++ b/src/states/windowmps.jl @@ -1,7 +1,7 @@ """ 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 @@ -26,64 +26,63 @@ a window state or a vector of tensors to construct the window. Alternatively, it 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. - WindowMPS(state::InfiniteMPS, L::Int) Construct a WindowMPS from an InfiniteMPS, by promoting a region of length `L` to a `FiniteMPS`. """ -struct WindowMPS{A<:GenericMPSTensor,B<:MPSBondTensor} <: AbstractFiniteMPS +struct WindowMPS{A<:GenericMPSTensor,B<:MPSBondTensor,Vₗ,Vᵣ} <: AbstractMPS left_gs::InfiniteMPS{A,B} window::FiniteMPS{A,B} right_gs::InfiniteMPS{A,B} - function WindowMPS(ψₗ::InfiniteMPS{A,B}, ψₘ::FiniteMPS{A,B}, - ψᵣ::InfiniteMPS{A,B}=copy(ψₗ)) where {A<:GenericMPSTensor, - B<:MPSBondTensor} + function WindowMPS{A,B,Vₗ,Vᵣ}(ψₗ::InfiniteMPS{A,B}, ψₘ::FiniteMPS{A,B}, + ψᵣ::InfiniteMPS{A,B}) where {A<:GenericMPSTensor, + B<:MPSBondTensor,Vₗ,Vᵣ} 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}(ψₗ, ψₘ, ψᵣ) + return new{A,B,Vₗ,Vᵣ}(ψₗ, ψₘ, ψᵣ) end + +end + +function WindowMPS(ψₗ::InfiniteMPS{A,B}, ψₘ::FiniteMPS{A,B}, + ψᵣ::InfiniteMPS{A,B}) where {A<:GenericMPSTensor, + B<:MPSBondTensor} + return WindowMPS{A,B,:V,:V}(ψₗ, ψₘ, ψᵣ) end +# alias to help dispatching? + #=========================================================================================== Constructors ===========================================================================================# -function WindowMPS_copied( - Ψₗ::InfiniteMPS{A,B}, Ψₘ::FiniteMPS{A,B}, Ψᵣ::InfiniteMPS{A,B}=Ψₗ -) where {A<:GenericMPSTensor,B<:MPSBondTensor} - return WindowMPS{A,B}(copy(Ψₗ), Ψₘ, copy(Ψᵣ)) -end - function WindowMPS(ψₗ::InfiniteMPS, site_tensors::AbstractVector{<:GenericMPSTensor}, - ψᵣ::InfiniteMPS=ψₗ) + ψᵣ::InfiniteMPS) return WindowMPS(ψₗ, FiniteMPS(site_tensors), ψᵣ) end function WindowMPS(f, elt, physspaces::Vector{<:Union{S,CompositeSpace{S}}}, maxvirtspace::S, ψₗ::InfiniteMPS, - ψᵣ::InfiniteMPS=ψₗ) where {S<:ElementarySpace} + ψᵣ::InfiniteMPS) where {S<:ElementarySpace} ψₘ = FiniteMPS(f, elt, physspaces, maxvirtspace; left=left_virtualspace(ψₗ, 0), right=right_virtualspace(ψᵣ, length(physspaces))) return WindowMPS(ψₗ, ψₘ, ψᵣ) end function WindowMPS(physspaces::Vector{<:Union{S,CompositeSpace{S}}}, maxvirtspace::S, - ψₗ::InfiniteMPS, ψᵣ::InfiniteMPS=ψₗ) where {S<:ElementarySpace} + ψₗ::InfiniteMPS, ψᵣ::InfiniteMPS) where {S<:ElementarySpace} return WindowMPS(rand, Defaults.eltype, physspaces, maxvirtspace, ψₗ, ψᵣ) end function WindowMPS(f, elt, physspaces::Vector{<:Union{S,CompositeSpace{S}}}, virtspaces::Vector{S}, ψₗ::InfiniteMPS, - ψᵣ::InfiniteMPS=ψₗ) where {S<:ElementarySpace} + ψᵣ::InfiniteMPS) where {S<:ElementarySpace} ψₘ = FiniteMPS(f, elt, physspaces, virtspaces) return WindowMPS(ψₗ, ψₘ, ψᵣ) end function WindowMPS(physspaces::Vector{<:Union{S,CompositeSpace{S}}}, virtspaces::Vector{S}, - ψₗ::InfiniteMPS, ψᵣ::InfiniteMPS=ψₗ) where {S<:ElementarySpace} + ψₗ::InfiniteMPS, ψᵣ::InfiniteMPS) where {S<:ElementarySpace} return WindowMPS(rand, Defaults.eltype, physspaces, virtspaces, ψₗ, ψᵣ) end @@ -116,13 +115,41 @@ function WindowMPS(ψ::InfiniteMPS{A,B}, L::Int) where {A,B} end #=========================================================================================== -Utility +Left and right variable status changes (usefull for time evolution) ===========================================================================================# +function fix_left(ψ::WindowMPS{A,B,Vₗ,Vᵣ}) where {A,B,Vₗ,Vᵣ} + return WindowMPS{A,B,:F,Vᵣ}(ψ.left_gs,ψ.window,ψ.right_gs) +end +function fix_left(ψ::WindowMPS{A,B,:V,:F}) where {A,B} + return copy(ψ.window) +end + +function fix_right(ψ::WindowMPS{A,B,Vₗ,Vᵣ}) where {A,B,Vₗ,Vᵣ} + return WindowMPS{A,B,Vₗ,:F}(ψ.left_gs,ψ.window,ψ.right_gs) +end +function fix_right(ψ::WindowMPS{A,B,:F,:V}) where {A,B} + return copy(ψ.window) +end + +function fix_infinite(ψ::WindowMPS) + return copy(ψ.window) +end -function Base.copy(ψ::WindowMPS) +#=========================================================================================== +Utility +===========================================================================================# +function Base.copy(ψ::WindowMPS{A,B,:V,:V}) return WindowMPS(copy(ψ.left_gs), copy(ψ.window), copy(ψ.right_gs)) end +function Base.copy(ψ::WindowMPS{A,B,:F,:V}) + return WindowMPS(ψ.left_gs, copy(ψ.window), copy(ψ.right_gs)) +end + +function Base.copy(ψ::WindowMPS{A,B,:V,:F}) + return WindowMPS(copy(ψ.left_gs), copy(ψ.window), ψ.right_gs) +end + # not sure about the underlying methods... Base.length(ψ::WindowMPS) = length(ψ.window) Base.size(ψ::WindowMPS, i...) = size(ψ.window, i...) From eacefe34218c8ac0abf1a3dd071f917636e78bcd Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Tue, 13 Feb 2024 17:07:11 +0100 Subject: [PATCH 05/62] some name changes --- src/MPSKit.jl | 3 ++- src/environments/FinEnv.jl | 2 +- src/environments/windowenv.jl | 6 +++--- 3 files changed, 6 insertions(+), 5 deletions(-) diff --git a/src/MPSKit.jl b/src/MPSKit.jl index f5ec8c508..795bb384c 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -91,10 +91,11 @@ include("transfermatrix/transfer.jl") abstract type Cache end # cache "manages" environments -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") diff --git a/src/environments/FinEnv.jl b/src/environments/FinEnv.jl index e60a282b3..b08737cb2 100644 --- a/src/environments/FinEnv.jl +++ b/src/environments/FinEnv.jl @@ -96,7 +96,7 @@ end #notify the cache that we updated in-place, so it should invalidate the dependencies # this forces the transfers to be recalculated lazily -function poison!(ca::FinEnv, ind) +function invalidate!(ca::FinEnv, ind) ca.ldependencies[ind] = similar(ca.ldependencies[ind]) ca.rdependencies[ind] = similar(ca.rdependencies[ind]) end diff --git a/src/environments/windowenv.jl b/src/environments/windowenv.jl index 82829ac3f..7d82fe671 100644 --- a/src/environments/windowenv.jl +++ b/src/environments/windowenv.jl @@ -25,13 +25,13 @@ function environments(below::WindowMPS, above::WindowMPS) end #notify the cache that we updated in-place, so it should invalidate the dependencies -poison!(ca::WindowEnv, ind) = poison!(ca.fin_env,ind) +invalidate!(ca::WindowEnv, ind) = invalidate!(ca.fin_env,ind) # check and recalculate the left and right start function check_rightinfenv!(ca::WindowEnv) if ca.right_env.rw[:,0] != ca.fin_env.rightstart[end]# !== doesn't work as intended - poison!(ca, lastindex(ca.fin_env.rightstart)) #forces transfers to be recalculated lazily + invalidate!(ca, lastindex(ca.fin_env.rightstart)) #forces transfers to be recalculated lazily ca.fin_env.rightstart = ca.right_env.rw[:,0] #do some other checks and recalcs for bonddimensions? @@ -41,7 +41,7 @@ end function check_leftinfenv!(ca::WindowEnv) if ca.left_env.lw[:,end+1] != ca.fin_env.leftstart[1]# !== doesn't work as intended - poison!(ca, firstindex(ca.fin_env.leftstart)) #forces transfers to be recalculated lazily + invalidate!(ca, firstindex(ca.fin_env.leftstart)) #forces transfers to be recalculated lazily ca.fin_env.leftstart = ca.left_env.rw[:,end+1] #do some other checks and recalcs for bonddimensions? From 667fb8a596faf24503239e82eedcade31022220f Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Tue, 13 Feb 2024 17:51:56 +0100 Subject: [PATCH 06/62] first fixes --- src/MPSKit.jl | 1 + src/states/windowmps.jl | 6 +++--- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/MPSKit.jl b/src/MPSKit.jl index 795bb384c..43aa49e71 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -14,6 +14,7 @@ using Base: @kwdef # bells and whistles for mpses export InfiniteMPS, FiniteMPS, WindowMPS, MPSMultiline export PeriodicArray, Window +export fix_left,fix_right,fix_infinite export MPSTensor export QP, LeftGaugedQP, RightGaugedQP export leftorth, diff --git a/src/states/windowmps.jl b/src/states/windowmps.jl index e77411c68..3c19271bb 100644 --- a/src/states/windowmps.jl +++ b/src/states/windowmps.jl @@ -138,15 +138,15 @@ end #=========================================================================================== Utility ===========================================================================================# -function Base.copy(ψ::WindowMPS{A,B,:V,:V}) +function Base.copy(ψ::WindowMPS{A,B,:V,:V}) where {A,B} return WindowMPS(copy(ψ.left_gs), copy(ψ.window), copy(ψ.right_gs)) end -function Base.copy(ψ::WindowMPS{A,B,:F,:V}) +function Base.copy(ψ::WindowMPS{A,B,:F,:V}) where {A,B} return WindowMPS(ψ.left_gs, copy(ψ.window), copy(ψ.right_gs)) end -function Base.copy(ψ::WindowMPS{A,B,:V,:F}) +function Base.copy(ψ::WindowMPS{A,B,:V,:F}) where {A,B} return WindowMPS(copy(ψ.left_gs), copy(ψ.window), ψ.right_gs) end From ef56ad917ae39432b7dfcaef2d7332e1525673c3 Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Wed, 14 Feb 2024 21:11:25 +0100 Subject: [PATCH 07/62] Update windowenv.jl --- src/environments/windowenv.jl | 77 +++++++++++++++++++++++------------ 1 file changed, 51 insertions(+), 26 deletions(-) diff --git a/src/environments/windowenv.jl b/src/environments/windowenv.jl index 7d82fe671..e7a4da2b6 100644 --- a/src/environments/windowenv.jl +++ b/src/environments/windowenv.jl @@ -3,19 +3,24 @@ 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} <: Cache +struct WindowEnv{A,B,C,D} <: Cache left_env::A fin_env::B right_env::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 -function WindowEnv(state::WindowMPS, O::Union{SparseMPO,MPOHamiltonian,DenseMPO}, above=nothing; lenvs=environments(state.left_gs, O),renvs=environments(state.right_gs, O)) - fin_env = environments(state, O, above, copy(leftenv(lenvs, 1, state.left_gs)), - copy(rightenv(renvs, length(state), state.right_gs))) - return WindowEnv(lenvs,fin_env,renvs) +# 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::Union{SparseMPO,MPOHamiltonian,DenseMPO}, above=nothing; lenvs=environments(ψ.left_gs, O),renvs=environments(ψ.right_gs, O)) + fin_env = environments(ψ, O, above, leftenv(lenvs, 1, ψ.left_gs), + rightenv(renvs, length(ψ), ψ.right_gs)) + return WindowEnv(lenvs,fin_env,renvs,copy(ψ.left_gs.AL),copy(ψ.right_gs.AR)) end +# is this intended for overlaps? we already have dot for this. function environments(below::WindowMPS, above::WindowMPS) above.left_gs == below.left_gs || throw(ArgumentError("left gs differs")) above.right_gs == below.right_gs || throw(ArgumentError("right gs differs")) @@ -27,39 +32,59 @@ end #notify the cache that we updated in-place, so it should invalidate the dependencies invalidate!(ca::WindowEnv, ind) = invalidate!(ca.fin_env,ind) -# check and recalculate the left and right start -function check_rightinfenv!(ca::WindowEnv) - if ca.right_env.rw[:,0] != ca.fin_env.rightstart[end]# !== doesn't work as intended - - invalidate!(ca, lastindex(ca.fin_env.rightstart)) #forces transfers to be recalculated lazily +# Check the infinite envs and recalculate the left and right start +function check_rightinfenv!(ca::WindowEnv, ψ::InfiniteMPS) + println("Doing right check") + if !all(ca.rinfdeps .=== ψ.AR) + println("changing right env") + invalidate!(ca, length(ψ)) #forces transfers to be recalculated lazily - ca.fin_env.rightstart = ca.right_env.rw[:,0] #do some other checks and recalcs for bonddimensions? - + ca.fin_env.rightenvs[end] = rightenv(ca.right_env, 0, ψ) #automatic recalculate of right_env + ca.rinfdeps .= ψ.AR + #do some other checks and recalcs for bonddimensions? end end -function check_leftinfenv!(ca::WindowEnv) - if ca.left_env.lw[:,end+1] != ca.fin_env.leftstart[1]# !== doesn't work as intended - - invalidate!(ca, firstindex(ca.fin_env.leftstart)) #forces transfers to be recalculated lazily +function check_leftinfenv!(ca::WindowEnv, ψ::InfiniteMPS) + println("Doing left check") + if !all(ca.linfdeps .=== ψ.AL) + println("changing left env") + invalidate!(ca, 1) #forces transfers to be recalculated lazily - ca.fin_env.leftstart = ca.left_env.rw[:,end+1] #do some other checks and recalcs for bonddimensions? - + # replace this line with a function to do this for lazy environments + ca.fin_env.leftenvs[1] = leftenv(ca.right_env, length(ψ)+1, ψ) + ca.linfdeps .= ψ.AL + #do some other checks and recalcs for bonddimensions? end end -#rightenv[ind] will be contracteable with the tensor on site [ind] -function rightenv(ca::WindowEnv, ind, state) - check_rightinfenv!(ca) - return rightenv(ca.fin_env, ind, state) +# only does the check when the env is variable +function check_infenv!(ca::WindowEnv, ψ::WindowMPS) + check_leftinfenv!(ca,ψ.left_gs) + check_rightinfenv!(ca,ψ.right_gs) +end + +function check_infenv!(ca::WindowEnv, ψ::WindowMPS{A,B,:F,Vᵣ}) where {A,B,Vᵣ} + check_rightinfenv!(ca,ψ.right_gs) end -function leftenv(ca::FinEnv, ind, state) - check_leftinfenv!(ca) - return leftenv(ca.fin_env, ind, state) +function check_infenv!(ca::WindowEnv, ψ::WindowMPS{A,B,Vₗ,:F}) where {A,B,Vₗ} + check_leftinfenv!(ca,ψ.left_gs) +end + +# for LazySum and the like, we do not want to wrap every subenv in a WindowEnv, +# so instead we will just put in a check before the derivatives are called +# we could consider something similar for expectation_value +for der = (:∂∂AC,:∂∂C,:∂∂AC2) + @eval begin + function $der(pos::Int,mps::WindowMPS,opp,ca::WindowEnv) + check_infenv!(ca, mps) + return $der(pos,mps.window,opp,ca.fin_env) + end + end end function fix_infinite(ψ::WindowMPS,env::WindowEnv) - newenv = FinEnv(env.above,env.opp,env.ldependencies,env.rdependencies,copy(env.leftenvs),copy(env.rightenvs)) + newenv = environments(ψ.window,env.fin_env.opp,env.fin_env.above,copy(leftenv(env.left_env, 1, ψ.left_gs)),copy(rightenv(env.right_env, length(ψ), ψ.right_gs))) return fix_infinite(ψ),newenv end From 86ff2c1f2b5f9de455a82617f1134a329ad8ce10 Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Mon, 19 Feb 2024 16:57:15 +0100 Subject: [PATCH 08/62] some renaming --- src/algorithms/toolbox.jl | 4 +- src/environments/multipleenv.jl | 4 +- src/states/orthoview.jl | 16 ++-- src/states/windowmps.jl | 147 +++++++++++++------------------- 4 files changed, 72 insertions(+), 99 deletions(-) diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl index 38b942939..f7da741ad 100644 --- a/src/algorithms/toolbox.jl +++ b/src/algorithms/toolbox.jl @@ -160,11 +160,11 @@ function variance(state::InfiniteQP, H::MPOHamiltonian, envs=environments(state, state.trivial || throw(ArgumentError("variance of domain wall excitations is not implemented")) - rescaled_H = H - expectation_value(state.left_gs, H) + rescaled_H = H - expectation_value(state.left, H) #I don't remember where the formula came from E_ex = dot(state, effective_excitation_hamiltonian(H, state, envs)) - E_f = expectation_value(state.left_gs, rescaled_H, 1:0) + E_f = expectation_value(state.left, rescaled_H, 1:0) H2 = rescaled_H * rescaled_H diff --git a/src/environments/multipleenv.jl b/src/environments/multipleenv.jl index 108ea99ef..e0888f0fe 100644 --- a/src/environments/multipleenv.jl +++ b/src/environments/multipleenv.jl @@ -32,8 +32,8 @@ end function environments(st::WindowMPS, H::LazySum; - lenvs=environments(st.left_gs, H), - renvs=environments(st.right_gs, H)) + lenvs=environments(st.left, H), + renvs=environments(st.right, H)) return MultipleEnvironments(H, map((op, sublenv, subrenv) -> environments(st, op; lenvs=sublenv, diff --git a/src/states/orthoview.jl b/src/states/orthoview.jl index 981e64fca..4273a9555 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/windowmps.jl b/src/states/windowmps.jl index 3c19271bb..f98a1ba91 100644 --- a/src/states/windowmps.jl +++ b/src/states/windowmps.jl @@ -1,3 +1,6 @@ +const WINDOW_FIXED = :F +const WINDOW_VARIABLE = :V + """ WindowMPS{A<:GenericMPSTensor,B<:MPSBondTensor} <: AbstractFiniteMPS @@ -5,56 +8,51 @@ Type that represents a finite Matrix Product State embedded in an infinite Matri ## 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, + 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. +left and right state to construct the WindowMPS in one step. WindowMPS(state::InfiniteMPS, L::Int) Construct a WindowMPS from an InfiniteMPS, by promoting a region of length `L` to a `FiniteMPS`. + +Options for fixing the left and right infinite state (i.e. so they don't get time evolved) +can be done via the keyword arguments `fixleft` and `fixright`. """ -struct WindowMPS{A<:GenericMPSTensor,B<:MPSBondTensor,Vₗ,Vᵣ} <: AbstractMPS - left_gs::InfiniteMPS{A,B} - window::FiniteMPS{A,B} - right_gs::InfiniteMPS{A,B} - - function WindowMPS{A,B,Vₗ,Vᵣ}(ψₗ::InfiniteMPS{A,B}, ψₘ::FiniteMPS{A,B}, - ψᵣ::InfiniteMPS{A,B}) where {A<:GenericMPSTensor, - B<:MPSBondTensor,Vₗ,Vᵣ} +struct WindowMPS{A<:GenericMPSTensor,B<:MPSBondTensor,VL,VR} <: AbstractMPS + window::Window{InfiniteMPS{A,B},FiniteMPS{A,B},InfiniteMPS{A,B}} + + function WindowMPS(ψₗ::InfiniteMPS{A,B}, ψₘ::FiniteMPS{A,B}, + ψᵣ::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,Vₗ,Vᵣ}(ψₗ, ψₘ, ψᵣ) + + VL = fixleft ? WINDOW_FIXED : WINDOW_VARIABLE + VR = fixright ? WINDOW_FIXED : WINDOW_VARIABLE + return new{A,B,VL,VR}(Window(ψₗ, ψₘ, ψᵣ)) end end -function WindowMPS(ψₗ::InfiniteMPS{A,B}, ψₘ::FiniteMPS{A,B}, - ψᵣ::InfiniteMPS{A,B}) where {A<:GenericMPSTensor, - B<:MPSBondTensor} - return WindowMPS{A,B,:V,:V}(ψₗ, ψₘ, ψᵣ) -end - -# alias to help dispatching? - #=========================================================================================== Constructors ===========================================================================================# @@ -100,7 +98,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) @@ -111,103 +109,78 @@ 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), ψ) -end - -#=========================================================================================== -Left and right variable status changes (usefull for time evolution) -===========================================================================================# -function fix_left(ψ::WindowMPS{A,B,Vₗ,Vᵣ}) where {A,B,Vₗ,Vᵣ} - return WindowMPS{A,B,:F,Vᵣ}(ψ.left_gs,ψ.window,ψ.right_gs) -end -function fix_left(ψ::WindowMPS{A,B,:V,:F}) where {A,B} - return copy(ψ.window) -end - -function fix_right(ψ::WindowMPS{A,B,Vₗ,Vᵣ}) where {A,B,Vₗ,Vᵣ} - return WindowMPS{A,B,Vₗ,:F}(ψ.left_gs,ψ.window,ψ.right_gs) -end -function fix_right(ψ::WindowMPS{A,B,:F,:V}) where {A,B} - return copy(ψ.window) -end - -function fix_infinite(ψ::WindowMPS) - return copy(ψ.window) + return WindowMPS(ψ, FiniteMPS(ALs, ARs, ACs, CLs), copyright ? copy(ψ) : ψ; kwargs...) end #=========================================================================================== Utility ===========================================================================================# -function Base.copy(ψ::WindowMPS{A,B,:V,:V}) where {A,B} - return WindowMPS(copy(ψ.left_gs), copy(ψ.window), copy(ψ.right_gs)) -end - -function Base.copy(ψ::WindowMPS{A,B,:F,:V}) where {A,B} - return WindowMPS(ψ.left_gs, copy(ψ.window), copy(ψ.right_gs)) +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{A,B,:V,:F}) where {A,B} - return WindowMPS(copy(ψ.left_gs), copy(ψ.window), ψ.right_gs) +function Base.copy(ψ::WindowMPS{A,B,VL,VR}) where {A,B,VL,VR} + left = VL === WINDOW_VARIABLE ? copy(ψ.left) : ψ.left + right = VR === WINDOW_VARIABLE ? copy(ψ.right) : ψ.right + return WindowMPS(left, copy(ψ.middle), right) 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 site_type(::Type{<:WindowMPS{A}}) where {A} = A bond_type(::Type{<:WindowMPS{<:Any,B}}) where {B} = B TensorKit.space(ψ::WindowMPS, n::Integer) = space(ψ.AC[n], 2) -left_virtualspace(ψ::WindowMPS, n::Integer) = left_virtualspace(ψ.window, n); -right_virtualspace(ψ::WindowMPS, n::Integer) = right_virtualspace(ψ.window, n); +left_virtualspace(ψ::WindowMPS, n::Integer) = left_virtualspace(ψ.middle, n); +right_virtualspace(ψ::WindowMPS, n::Integer) = right_virtualspace(ψ.middle, n); -r_RR(ψ::WindowMPS) = r_RR(ψ.right_gs, length(ψ)) -l_LL(ψ::WindowMPS) = l_LL(ψ.left_gs, 1) - -function Base.getproperty(ψ::WindowMPS, prop::Symbol) - if prop == :AL - return ALView(ψ) - elseif prop == :AR - return ARView(ψ) - elseif prop == :AC - return ACView(ψ) - elseif prop == :CR - return CRView(ψ) - else - return getfield(ψ, prop) - end -end +r_RR(ψ::WindowMPS) = r_RR(ψ.right, length(ψ)) +l_LL(ψ::WindowMPS) = l_LL(ψ.left, 1) -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(ψ)) From 5e0ab0246b80d08eba66cda62f6f38f6aafda75c Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Mon, 19 Feb 2024 16:59:46 +0100 Subject: [PATCH 09/62] first version of windowtdvp --- src/MPSKit.jl | 1 + src/algorithms/timestep/windowtdvp.jl | 137 ++++++++------------- src/environments/windowenv.jl | 169 +++++++++++++++++++------- 3 files changed, 176 insertions(+), 131 deletions(-) diff --git a/src/MPSKit.jl b/src/MPSKit.jl index 43aa49e71..8ec314ffb 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -120,6 +120,7 @@ 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/groundstate/vumps.jl") include("algorithms/groundstate/idmrg.jl") diff --git a/src/algorithms/timestep/windowtdvp.jl b/src/algorithms/timestep/windowtdvp.jl index 91d289665..b0b626392 100644 --- a/src/algorithms/timestep/windowtdvp.jl +++ b/src/algorithms/timestep/windowtdvp.jl @@ -1,103 +1,61 @@ +# make a struct for WindowTDVP, finalize for bondexpansion and glue after that -function _update_leftEnv!( - nleft::InfiniteMPS, WindowEnv::Window{E,F,E} -) where {E<:Cache,F<:Cache} - l = leftenv(WindowEnv.left, 1, nleft) - WindowEnv.middle.ldependencies[:] = similar.(WindowEnv.middle.ldependencies) # forget the old left dependencies - this forces recalculation whenever leftenv is called - WindowEnv.middle.leftenvs[1] = l - - return Window(WindowEnv.left, WindowEnv.middle, WindowEnv.right) -end +# we should use finenv here explicitly probably +function timestep!(Ψ::WindowMPS{A,B,VL,VR},H,t::Number,dt::Number,alg::TDVP,env=environments(Ψ, H)) where {A,B,VL,VR} + #first evolve left state + if VL === WINDOW_VARIABLE + println("Doing left") + nleft, _ = timestep(Ψ.left, H, t, dt, alg, env.left; leftorthflag=true) #env gets updated in place + Ψ = WindowMPS(nleft,Ψ.middle,Ψ.right) + end -function _update_rightEnv!( - nright::InfiniteMPS, WindowEnv::Window{E,F,E} -) where {E<:Cache,F<:Cache} - r = rightenv(WindowEnv.right, length(nright), nright) - WindowEnv.middle.rdependencies[:] = similar.(WindowEnv.middle.rdependencies) # forget the old right dependencies - this forces recalculation - WindowEnv.middle.rightenvs[end] = r + # left to right sweep on window + for i in 1:(length(Ψ) - 1) + h_ac = ∂∂AC(i, Ψ, H, env) + Ψ.AC[i] = integrate(h_ac, Ψ.AC[i], t, dt / 2, alg.integrator) - return Window(WindowEnv.left, WindowEnv.middle, WindowEnv.right) -end + h_c = ∂∂C(i, Ψ, H, env) + Ψ.CR[i] = integrate(h_c, Ψ.CR[i], t, -dt / 2, alg.integrator) + end -#Note: this needs to be tested -function leftexpand( - st::WindowMPS, - H, - Envs; - singval=1e-2, - growspeed=10, -) - (U, S, V) = tsvd(st.left_gs.CR[1]; alg=TensorKit.SVD()) + h_ac = ∂∂AC(length(Ψ), Ψ, H, env) + Ψ.AC[end] = integrate(h_ac, Ψ.AC[end], t, dt / 2, alg.integrator) - if minimum([minimum(diag(v)) for (k, v) in blocks(S)]) > singval - (nst, _) = changebonds( - st.left_gs, H, OptimalExpand(; trscheme=truncbelow(singval, growspeed)), Envs - ) + if VR === WINDOW_VARIABLE + println("Doing right") + nright, _ = timestep(Ψ.right, H, t + dt, dt, alg, env.right; leftorthflag=false) # env gets updated in place + Ψ = WindowMPS(Ψ.left,Ψ.middle,nright) + end - # the AL-bond dimension changed, and therefore our window also needs updating - v = TensorMap( - rand, ComplexF64, left_virtualspace(nst, 0), left_virtualspace(st.left_gs, 0) - ) - (vals, vecs) = eigsolve( - flip(TransferMatrix(st.left_gs.AL, nst.AL)), v, 1, :LM, Arnoldi() - ) - rho = pinv(nst.CR[0]) * vecs[1] * st.left_gs.CR[0] #CR[0] == CL[1] - st.AC[1] = _transpose_front(normalize!(rho * st.CR[0]) * _transpose_tail(st.AR[1])) + # right to left sweep on window + for i in length(Ψ):-1:2 + h_ac = ∂∂AC(i, Ψ, H, env) + Ψ.AC[i] = integrate(h_ac, Ψ.AC[i], t + dt / 2, dt / 2, alg.integrator) - recalculate!(Envs, nst) #updates left infinite env based on expanded state - return nst, Envs + h_c = ∂∂C(i - 1, Ψ, H, env) + Ψ.CR[i - 1] = integrate(h_c, Ψ.CR[i - 1], t + dt / 2, -dt / 2, alg.integrator) end - return st.left_gs, Envs -end -function rightexpand( - st::WindowMPS, - H, - Envs; - singval=1e-2, - growspeed=10, -) - (U, S, V) = tsvd(st.left_gs.CR[1]; alg=TensorKit.SVD()) - - if minimum([minimum(diag(v)) for (k, v) in blocks(S)]) > singval - (nst, _) = changebonds( - st.right_gs, - H, - OptimalExpand(; trscheme=truncbelow(singval, growspeed)), - Envs, - ) - - v = TensorMap( - rand, ComplexF64, right_virtualspace(st.right_gs, 0), right_virtualspace(nst, 0) - ) - (vals, vecs) = eigsolve( - TransferMatrix(st.right_gs.AR, nst.AR), v, 1, :LM, Arnoldi() - ) - rho = st.right_gs.CR[0] * vecs[1] * pinv(nst.CR[0]) - st.AC[end] = st.AL[end] * normalize!(st.CR[end] * rho) + h_ac = ∂∂AC(1, Ψ, H, env) + Ψ.AC[1] = integrate(h_ac, Ψ.AC[1], t, dt / 2, alg.integrator) - recalculate!(Envs, nst) #updates right infinite env based on expanded state - return nst, Envs - end - return st.right_gs, Envs + return Ψ, env end +#= + function timestep!( - Ψ::WindowMPS, - H::Window, + Ψ::WindowMPS{A,B,VL,VR}, + H, t::Number, dt::Number, alg::TDVP, - env::Window=environments(Ψ, H); - leftevolve=true, - rightevolve=true, -) + env=environments(Ψ, H); +) where {A,B,VL,VR} #first evolve left state - if leftevolve - nleft, _ = timestep(Ψ.left_gs, H.left, t, dt, alg, env.left; leftorthflag=true) #env gets updated in place - _update_leftEnv!(nleft, env) - else - nleft = Ψ.left_gs + if VL === WINDOW_VARIABLE + nleft, _ = timestep(Ψ.left, H.left, t, dt, alg, env.left; leftorthflag=true) #env gets updated in place + Ψ = WindowMPS(nleft,Ψ.middle,Ψ.right) end # left to right sweep on window @@ -112,11 +70,9 @@ function timestep!( h_ac = ∂∂AC(length(Ψ), Ψ, H.middle, env.middle) Ψ.AC[end] = integrate(h_ac, Ψ.AC[end], t, dt / 2, alg.integrator) - if rightevolve + if VR === WINDOW_VARIABLE nright, _ = timestep(Ψ.right_gs, H.right, t + dt, dt, alg, env.right; leftorthflag=false) # env gets updated in place - _update_rightEnv!(nright, env) - else - nright = Ψ.right_gs + Ψ = WindowMPS(Ψ.left,Ψ.middle,nright) end # right to left sweep on window @@ -131,9 +87,11 @@ function timestep!( h_ac = ∂∂AC(1, Ψ, H.middle, env.middle) Ψ.AC[1] = integrate(h_ac, Ψ.AC[1], t, dt / 2, alg.integrator) - return WindowMPS(nleft, Ψ.window, nright), env + return Ψ, env end +=# +#= function timestep!( Ψ::WindowMPS, H::Window, @@ -159,7 +117,7 @@ function timestep!( ) _update_leftEnv!(nleft, env) else - nleft = Ψ.left_gs + nleft = Ψ.left end # left to right sweep on window @@ -209,3 +167,4 @@ function timestep!( return WindowMPS(nleft, Ψ.window, nright), env end +=# diff --git a/src/environments/windowenv.jl b/src/environments/windowenv.jl index e7a4da2b6..f7b06fef7 100644 --- a/src/environments/windowenv.jl +++ b/src/environments/windowenv.jl @@ -4,87 +4,172 @@ """ struct WindowEnv{A,B,C,D} <: Cache - left_env::A - fin_env::B - right_env::C + window::Window{A,B,C} linfdeps::PeriodicArray{D} #the data we used to calculate leftenvs/rightenvs rinfdeps::PeriodicArray{D} end +const WindowEnvUnion = Union{T,Window{A,T,B}} where {T<:WindowEnv,A,B} + #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::Union{SparseMPO,MPOHamiltonian,DenseMPO}, above=nothing; lenvs=environments(ψ.left_gs, O),renvs=environments(ψ.right_gs, O)) - fin_env = environments(ψ, O, above, leftenv(lenvs, 1, ψ.left_gs), - rightenv(renvs, length(ψ), ψ.right_gs)) - return WindowEnv(lenvs,fin_env,renvs,copy(ψ.left_gs.AL),copy(ψ.right_gs.AR)) +function environments(ψ::WindowMPS, O::Union{SparseMPO,MPOHamiltonian,DenseMPO}, above=nothing; lenvs=environments(ψ.left, O),renvs=environments(ψ.right, O)) + fin_env = environments(ψ, O, above, leftenv(lenvs, 1, ψ.left), + rightenv(renvs, length(ψ), ψ.right)) + return WindowEnv(Window(lenvs,fin_env,renvs),copy(ψ.left.AL),copy(ψ.right.AR)) +end + +# when supplied with a Window Hamiltonian, we cannot assume H.left/H.right is the same as H.middle +function environments(ψ::WindowMPS, O::Window, above=nothing; lenvs=environments(ψ.left, O.left),renvs=environments(ψ.right, O.right)) + window_env = Window(environments(ψ.left, O.middle),environments(ψ.middle, O.middle),environments(ψ.left, O.middle)) + return Window(lenvs,WindowEnv(window_env,copy(ψ.left.AL),copy(ψ.right.AR)),renvs) end -# is this intended for overlaps? we already have dot for this. function environments(below::WindowMPS, above::WindowMPS) - above.left_gs == below.left_gs || throw(ArgumentError("left gs differs")) - above.right_gs == below.right_gs || throw(ArgumentError("right gs differs")) + 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) + 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 +finenv(ca::Window{A,<:WindowEnv,B},ψ::WindowMPS) where {A,B} = finenv(ca.middle,ψ) + +left_of_finenv(ca::WindowEnv) = ca.left +right_of_finenv(ca::WindowEnv) = ca.right +left_of_finenv(ca::Window{A,<:WindowEnv,B}) where {A,B} = ca.middle.left +right_of_finenv(ca::Window{A,<:WindowEnv,B}) where {A,B} = ca.middle.right + #notify the cache that we updated in-place, so it should invalidate the dependencies -invalidate!(ca::WindowEnv, ind) = invalidate!(ca.fin_env,ind) +invalidate!(ca::WindowEnv, ind) = invalidate!(ca.middle,ind) -# Check the infinite envs and recalculate the left and right start -function check_rightinfenv!(ca::WindowEnv, ψ::InfiniteMPS) +# Check the infinite envs and recalculate the right start +function check_rightinfenv!(ca::WindowEnv, ψ::WindowMPS) println("Doing right check") - if !all(ca.rinfdeps .=== ψ.AR) + if !all(ca.rinfdeps .=== ψ.right.AR) println("changing right env") - invalidate!(ca, length(ψ)) #forces transfers to be recalculated lazily + invalidate!(ca, length(ψ.middle)) #forces transfers to be recalculated lazily - ca.fin_env.rightenvs[end] = rightenv(ca.right_env, 0, ψ) #automatic recalculate of right_env - ca.rinfdeps .= ψ.AR + update_rightstart!(ca.middle,ca.right,ψ.right) + glue_right!(ψ,ca.rinfdeps) # + ca.rinfdeps .= ψ.right.AR #do some other checks and recalcs for bonddimensions? end end -function check_leftinfenv!(ca::WindowEnv, ψ::InfiniteMPS) +function check_leftinfenv!(ca::WindowEnv, ψ::WindowMPS) println("Doing left check") - if !all(ca.linfdeps .=== ψ.AL) + if !all(ca.linfdeps .=== ψ.left.AL) println("changing left env") invalidate!(ca, 1) #forces transfers to be recalculated lazily # replace this line with a function to do this for lazy environments - ca.fin_env.leftenvs[1] = leftenv(ca.right_env, length(ψ)+1, ψ) - ca.linfdeps .= ψ.AL - #do some other checks and recalcs for bonddimensions? + update_leftstart!(ca.middle,ca.left,ψ.left) + glue_left!(ψ,ca.linfdeps) # + ca.linfdeps .= ψ.left.AL + end +end + +#replace ψ.left by appropriate site index +function glue_left!(ψ::WindowMPS,oldALs) + #do we need C of finite part too? + newD = left_virtualspace(ψ.left, 0) + oldD = left_virtualspace(ψ, 0) + if newD == oldD + return nothing + end + v = TensorMap(rand, ComplexF64, newD, oldD) + (vals, vecs) = eigsolve( + flip(TransferMatrix(oldAls, ψ.left.AL)), v, 1, :LM, Arnoldi() + ) + rho = pinv(ψ.left.CR[0]) * vecs[1] * ψ.CR[0] #CR[0] == CL[1] + ψ.AC[1] = _transpose_front(normalize!(rho * ψ.CR[0]) * _transpose_tail(ψ.AR[1])) +end + +function glue_right!(ψ::WindowMPS,oldARs) + newD = right_virtualspace(ψ.right, 0) + oldD = right_virtualspace(ψ, 0) + if newD == oldD + return nothing end + v = TensorMap(rand, ComplexF64, oldD, newD) + (vals, vecs) = eigsolve( + TransferMatrix(oldARs, ψ.right.AR), v, 1, :LM, Arnoldi() + ) + rho = ψ.CR[0] * vecs[1] * pinv(ψ.right.CR[0]) + ψ.AC[end] = ψ.AL[end] * normalize!(ψ.CR[end] * rho) end -# only does the check when the env is variable -function check_infenv!(ca::WindowEnv, ψ::WindowMPS) - check_leftinfenv!(ca,ψ.left_gs) - check_rightinfenv!(ca,ψ.right_gs) +# these are to be extended for LazySum and the like +function update_leftstart!(ca_fin::FinEnv,ca_infin::MPOHamInfEnv, ψ::InfiniteMPS) + ca_fin.leftenvs[1] = leftenv(ca_infin, length(ψ)+1, ψ) end -function check_infenv!(ca::WindowEnv, ψ::WindowMPS{A,B,:F,Vᵣ}) where {A,B,Vᵣ} - check_rightinfenv!(ca,ψ.right_gs) +function update_rightstart!(ca_fin::FinEnv,ca_infin::MPOHamInfEnv, ψ::InfiniteMPS) + ca_fin.rightenvs[end] = rightenv(ca_infin, 0, ψ) end -function check_infenv!(ca::WindowEnv, ψ::WindowMPS{A,B,Vₗ,:F}) where {A,B,Vₗ} - check_leftinfenv!(ca,ψ.left_gs) +# under review +function leftenv(ca::WindowEnvUnion, ind, ψ::WindowMPS) + if ind < 1 + return leftenv(left_of_finenv(ca),ind,ψ.left) + elseif ind > length(ψ) + return leftenv(right_of_finenv(ca),ind,ψ.right) + else + return leftenv(finenv(ca,ψ),ind,ψ) + end end -# for LazySum and the like, we do not want to wrap every subenv in a WindowEnv, -# so instead we will just put in a check before the derivatives are called -# we could consider something similar for expectation_value -for der = (:∂∂AC,:∂∂C,:∂∂AC2) - @eval begin - function $der(pos::Int,mps::WindowMPS,opp,ca::WindowEnv) - check_infenv!(ca, mps) - return $der(pos,mps.window,opp,ca.fin_env) - end +function rightenv(ca::WindowEnvUnion, ind, ψ::WindowMPS) + if ind < 1 + return rightenv(left_of_finenv(ca),ind,ψ.left) + elseif ind > length(ψ) + return rightenv(right_of_finenv(ca),ind,ψ.right) + else + return rightenv(finenv(ca,ψ),ind,ψ) end end -function fix_infinite(ψ::WindowMPS,env::WindowEnv) - newenv = environments(ψ.window,env.fin_env.opp,env.fin_env.above,copy(leftenv(env.left_env, 1, ψ.left_gs)),copy(rightenv(env.right_env, length(ψ), ψ.right_gs))) - return fix_infinite(ψ),newenv +# to be moved +function expectation_value(Ψ::WindowMPS, windowH::Window, windowEnvs::Window{A,<:WindowEnv,B}=environments(Ψ, windowH)) where {A,B} + left = expectation_value(Ψ.left, windowH.left, windowEnvs.left) + middle = expectation_value(Ψ.middle, windowH.middle, finenv(windowEnvs,Ψ)) + right = expectation_value(Ψ.right, windowH.right, windowEnvs.right) + return left,middle,right +end + +function expectation_value(Ψ::WindowMPS, H, windowEnvs::WindowEnv=environments(Ψ, H)) + left = expectation_value(Ψ.left, H, windowEnvs.left) + middle = expectation_value(Ψ.middle, H, finenv(windowEnvs,Ψ)) + right = expectation_value(Ψ.right, H, windowEnvs.right) + return left,middle,right +end + + +# do we need this? I find it convenvient +function forget_infinite(ψ::WindowMPS,env::WindowEnvUnion) + fin_env = finenv(env,ψ) + newenv = environments(ψ.middle,fin_env.opp,fin_env.above,copy(leftenv(env, 1, ψ)),copy(rightenv(env, length(ψ), ψ))) + return copy(ψ.middle),newenv end From a855e3f7d346123d12d0eebba0b0304d5d0737cb Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Mon, 19 Feb 2024 20:54:12 +0100 Subject: [PATCH 10/62] fix WindowMPS subtyping --- src/states/windowmps.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/states/windowmps.jl b/src/states/windowmps.jl index f98a1ba91..1ae63c650 100644 --- a/src/states/windowmps.jl +++ b/src/states/windowmps.jl @@ -37,7 +37,7 @@ Construct a WindowMPS from an InfiniteMPS, by promoting a region of length `L` t Options for fixing the left and right infinite state (i.e. so they don't get time evolved) can be done via the keyword arguments `fixleft` and `fixright`. """ -struct WindowMPS{A<:GenericMPSTensor,B<:MPSBondTensor,VL,VR} <: AbstractMPS +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}, From e4242e28764143c1bf2b93697bf763f07e3a399a Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Mon, 19 Feb 2024 20:54:59 +0100 Subject: [PATCH 11/62] note change --- src/states/window.jl | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/src/states/window.jl b/src/states/window.jl index 415bdb453..2c7028683 100644 --- a/src/states/window.jl +++ b/src/states/window.jl @@ -1,4 +1,4 @@ -# Note : this is intended to be a template for windowmps and windows of operators/environments but this clashes with abstractfinitemps. +# Note : this is intended to be a template for windowmps and windows of operators/environments " Window(leftstate,window,rightstate) @@ -9,7 +9,4 @@ struct Window{L,M,R} left::L middle::M right::R -end - -# do we need copy? -# Base.copy(win::Window) = Window(copy(win.left),copy(win.middle),copy(win.right)) \ No newline at end of file +end \ No newline at end of file From 393786a7a1778684b7aef560d1ac5a9b821e71b0 Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Mon, 19 Feb 2024 20:55:15 +0100 Subject: [PATCH 12/62] Update windowtdvp.jl --- src/algorithms/timestep/windowtdvp.jl | 25 +++++++++++++++---------- 1 file changed, 15 insertions(+), 10 deletions(-) diff --git a/src/algorithms/timestep/windowtdvp.jl b/src/algorithms/timestep/windowtdvp.jl index b0b626392..9734b8efc 100644 --- a/src/algorithms/timestep/windowtdvp.jl +++ b/src/algorithms/timestep/windowtdvp.jl @@ -1,42 +1,47 @@ # make a struct for WindowTDVP, finalize for bondexpansion and glue after that -# we should use finenv here explicitly probably -function timestep!(Ψ::WindowMPS{A,B,VL,VR},H,t::Number,dt::Number,alg::TDVP,env=environments(Ψ, H)) where {A,B,VL,VR} +function timestep!(Ψ::WindowMPS{A,B,VL,VR},H,t::Number,dt::Number,alg::TDVP,env::Cache=environments(Ψ, H)) where {A,B,VL,VR} + if H isa Window + Hl,Hm,Hr = (H.left,H.middle,H.right) + else + Hl = Hm = Hr = H + end + #first evolve left state if VL === WINDOW_VARIABLE println("Doing left") - nleft, _ = timestep(Ψ.left, H, t, dt, alg, env.left; leftorthflag=true) #env gets updated in place + nleft, _ = timestep(Ψ.left, Hl, t, dt, alg, env.left; leftorthflag=true) #env gets updated in place Ψ = WindowMPS(nleft,Ψ.middle,Ψ.right) end # left to right sweep on window for i in 1:(length(Ψ) - 1) - h_ac = ∂∂AC(i, Ψ, H, env) + h_ac = ∂∂AC(i, Ψ, Hm, finenv(env,Ψ)) Ψ.AC[i] = integrate(h_ac, Ψ.AC[i], t, dt / 2, alg.integrator) - h_c = ∂∂C(i, Ψ, H, env) + h_c = ∂∂C(i, Ψ, Hm, finenv(env,Ψ)) Ψ.CR[i] = integrate(h_c, Ψ.CR[i], t, -dt / 2, alg.integrator) end - h_ac = ∂∂AC(length(Ψ), Ψ, H, env) + h_ac = ∂∂AC(length(Ψ), Ψ, Hm, finenv(env,Ψ)) Ψ.AC[end] = integrate(h_ac, Ψ.AC[end], t, dt / 2, alg.integrator) if VR === WINDOW_VARIABLE println("Doing right") - nright, _ = timestep(Ψ.right, H, t + dt, dt, alg, env.right; leftorthflag=false) # env gets updated in place + nright, _ = timestep(Ψ.right, Hr, t + dt, dt, alg, env.right; leftorthflag=false) # env gets updated in place Ψ = WindowMPS(Ψ.left,Ψ.middle,nright) end # right to left sweep on window for i in length(Ψ):-1:2 - h_ac = ∂∂AC(i, Ψ, H, env) + h_ac = ∂∂AC(i, Ψ, Hm, finenv(env,Ψ)) Ψ.AC[i] = integrate(h_ac, Ψ.AC[i], t + dt / 2, dt / 2, alg.integrator) - h_c = ∂∂C(i - 1, Ψ, H, env) + h_c = ∂∂C(i - 1, Ψ, Hm, finenv(env,Ψ)) Ψ.CR[i - 1] = integrate(h_c, Ψ.CR[i - 1], t + dt / 2, -dt / 2, alg.integrator) end - h_ac = ∂∂AC(1, Ψ, H, env) + h_ac = ∂∂AC(1, Ψ, Hm, finenv(env,Ψ)) Ψ.AC[1] = integrate(h_ac, Ψ.AC[1], t, dt / 2, alg.integrator) return Ψ, env From fe96dd7af10fa3164e5a7a152bd1c8112ccefc87 Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Mon, 19 Feb 2024 20:55:27 +0100 Subject: [PATCH 13/62] tmp version for windowenv.jl --- src/environments/windowenv.jl | 35 ++++++++++++++++++----------------- 1 file changed, 18 insertions(+), 17 deletions(-) diff --git a/src/environments/windowenv.jl b/src/environments/windowenv.jl index f7b06fef7..66afda5bd 100644 --- a/src/environments/windowenv.jl +++ b/src/environments/windowenv.jl @@ -56,11 +56,6 @@ function finenv(ca::WindowEnv,ψ::WindowMPS{A,B,VL,VR}) where {A,B,VL,VR} end finenv(ca::Window{A,<:WindowEnv,B},ψ::WindowMPS) where {A,B} = finenv(ca.middle,ψ) -left_of_finenv(ca::WindowEnv) = ca.left -right_of_finenv(ca::WindowEnv) = ca.right -left_of_finenv(ca::Window{A,<:WindowEnv,B}) where {A,B} = ca.middle.left -right_of_finenv(ca::Window{A,<:WindowEnv,B}) where {A,B} = ca.middle.right - #notify the cache that we updated in-place, so it should invalidate the dependencies invalidate!(ca::WindowEnv, ind) = invalidate!(ca.middle,ind) @@ -72,9 +67,8 @@ function check_rightinfenv!(ca::WindowEnv, ψ::WindowMPS) invalidate!(ca, length(ψ.middle)) #forces transfers to be recalculated lazily update_rightstart!(ca.middle,ca.right,ψ.right) - glue_right!(ψ,ca.rinfdeps) # + glue_right!(ψ,ca.rinfdeps) ca.rinfdeps .= ψ.right.AR - #do some other checks and recalcs for bonddimensions? end end @@ -84,9 +78,8 @@ function check_leftinfenv!(ca::WindowEnv, ψ::WindowMPS) println("changing left env") invalidate!(ca, 1) #forces transfers to be recalculated lazily - # replace this line with a function to do this for lazy environments update_leftstart!(ca.middle,ca.left,ψ.left) - glue_left!(ψ,ca.linfdeps) # + glue_left!(ψ,ca.linfdeps) ca.linfdeps .= ψ.left.AL end end @@ -96,28 +89,29 @@ function glue_left!(ψ::WindowMPS,oldALs) #do we need C of finite part too? newD = left_virtualspace(ψ.left, 0) oldD = left_virtualspace(ψ, 0) - if newD == oldD - return nothing - end + #if newD == oldD + # return nothing + #end v = TensorMap(rand, ComplexF64, newD, oldD) (vals, vecs) = eigsolve( - flip(TransferMatrix(oldAls, ψ.left.AL)), v, 1, :LM, Arnoldi() + flip(TransferMatrix(oldALs, ψ.left.AL)), v, 1, :LM, Arnoldi() ) rho = pinv(ψ.left.CR[0]) * vecs[1] * ψ.CR[0] #CR[0] == CL[1] ψ.AC[1] = _transpose_front(normalize!(rho * ψ.CR[0]) * _transpose_tail(ψ.AR[1])) end +#space mismatch here still function glue_right!(ψ::WindowMPS,oldARs) newD = right_virtualspace(ψ.right, 0) oldD = right_virtualspace(ψ, 0) - if newD == oldD - return nothing - end + #if newD == oldD + # return nothing + #end v = TensorMap(rand, ComplexF64, oldD, newD) (vals, vecs) = eigsolve( TransferMatrix(oldARs, ψ.right.AR), v, 1, :LM, Arnoldi() ) - rho = ψ.CR[0] * vecs[1] * pinv(ψ.right.CR[0]) + rho = ψ.CR[end] * vecs[1] * pinv(ψ.right.CR[0]) ψ.AC[end] = ψ.AL[end] * normalize!(ψ.CR[end] * rho) end @@ -131,6 +125,12 @@ function update_rightstart!(ca_fin::FinEnv,ca_infin::MPOHamInfEnv, ψ::InfiniteM end # under review +#= +left_of_finenv(ca::WindowEnv) = ca.left +right_of_finenv(ca::WindowEnv) = ca.right +left_of_finenv(ca::Window{A,<:WindowEnv,B}) where {A,B} = ca.middle.left +right_of_finenv(ca::Window{A,<:WindowEnv,B}) where {A,B} = ca.middle.right + function leftenv(ca::WindowEnvUnion, ind, ψ::WindowMPS) if ind < 1 return leftenv(left_of_finenv(ca),ind,ψ.left) @@ -150,6 +150,7 @@ function rightenv(ca::WindowEnvUnion, ind, ψ::WindowMPS) return rightenv(finenv(ca,ψ),ind,ψ) end end +=# # to be moved function expectation_value(Ψ::WindowMPS, windowH::Window, windowEnvs::Window{A,<:WindowEnv,B}=environments(Ψ, windowH)) where {A,B} From 07d32692b3ebad05f02b45fd9025bf22095c05b2 Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Wed, 21 Feb 2024 17:17:03 +0100 Subject: [PATCH 14/62] some refactoring --- src/algorithms/timestep/windowtdvp.jl | 25 ++++----- src/environments/multipleenv.jl | 11 ---- src/environments/windowenv.jl | 74 ++++++++++++++++----------- 3 files changed, 54 insertions(+), 56 deletions(-) diff --git a/src/algorithms/timestep/windowtdvp.jl b/src/algorithms/timestep/windowtdvp.jl index 9734b8efc..4365eb0fa 100644 --- a/src/algorithms/timestep/windowtdvp.jl +++ b/src/algorithms/timestep/windowtdvp.jl @@ -1,47 +1,42 @@ # make a struct for WindowTDVP, finalize for bondexpansion and glue after that -function timestep!(Ψ::WindowMPS{A,B,VL,VR},H,t::Number,dt::Number,alg::TDVP,env::Cache=environments(Ψ, H)) where {A,B,VL,VR} - if H isa Window - Hl,Hm,Hr = (H.left,H.middle,H.right) - else - Hl = Hm = Hr = H - end - +function timestep(Ψ::WindowMPS{A,B,VL,VR},H::Union{Window,LazySum},t::Number,dt::Number,alg::TDVP,env::Cache=environments(Ψ, H)) where {A,B,VL,VR} #first evolve left state if VL === WINDOW_VARIABLE println("Doing left") - nleft, _ = timestep(Ψ.left, Hl, t, dt, alg, env.left; leftorthflag=true) #env gets updated in place + nleft, _ = timestep(Ψ.left, H.left, t, dt, alg, env.left; leftorthflag=true) #env gets updated in place Ψ = WindowMPS(nleft,Ψ.middle,Ψ.right) end # left to right sweep on window for i in 1:(length(Ψ) - 1) - h_ac = ∂∂AC(i, Ψ, Hm, finenv(env,Ψ)) + h_ac = ∂∂AC(i, Ψ, H.middle, env) Ψ.AC[i] = integrate(h_ac, Ψ.AC[i], t, dt / 2, alg.integrator) - h_c = ∂∂C(i, Ψ, Hm, finenv(env,Ψ)) + h_c = ∂∂C(i, Ψ, H.middle, env) Ψ.CR[i] = integrate(h_c, Ψ.CR[i], t, -dt / 2, alg.integrator) end - h_ac = ∂∂AC(length(Ψ), Ψ, Hm, finenv(env,Ψ)) + h_ac = ∂∂AC(length(Ψ), Ψ, H.middle, env) Ψ.AC[end] = integrate(h_ac, Ψ.AC[end], t, dt / 2, alg.integrator) + #then evolve right state if VR === WINDOW_VARIABLE println("Doing right") - nright, _ = timestep(Ψ.right, Hr, t + dt, dt, alg, env.right; leftorthflag=false) # env gets updated in place + nright, _ = timestep(Ψ.right, H.right, t + dt, dt, alg, env.right; leftorthflag=false) # env gets updated in place Ψ = WindowMPS(Ψ.left,Ψ.middle,nright) end # right to left sweep on window for i in length(Ψ):-1:2 - h_ac = ∂∂AC(i, Ψ, Hm, finenv(env,Ψ)) + h_ac = ∂∂AC(i, Ψ, H.middle, env) Ψ.AC[i] = integrate(h_ac, Ψ.AC[i], t + dt / 2, dt / 2, alg.integrator) - h_c = ∂∂C(i - 1, Ψ, Hm, finenv(env,Ψ)) + h_c = ∂∂C(i - 1, Ψ, H.middle, env) Ψ.CR[i - 1] = integrate(h_c, Ψ.CR[i - 1], t + dt / 2, -dt / 2, alg.integrator) end - h_ac = ∂∂AC(1, Ψ, Hm, finenv(env,Ψ)) + h_ac = ∂∂AC(1, Ψ, H.middle, env) Ψ.AC[1] = integrate(h_ac, Ψ.AC[1], t, dt / 2, alg.integrator) return Ψ, env diff --git a/src/environments/multipleenv.jl b/src/environments/multipleenv.jl index e0888f0fe..bb125e2b6 100644 --- a/src/environments/multipleenv.jl +++ b/src/environments/multipleenv.jl @@ -30,17 +30,6 @@ end # return MultipleEnvironments(H, broadcast(o -> environments(state, o), H.opps)) # end; -function environments(st::WindowMPS, - H::LazySum; - lenvs=environments(st.left, H), - renvs=environments(st.right, H)) - return MultipleEnvironments(H, - map((op, sublenv, subrenv) -> environments(st, op; - lenvs=sublenv, - renvs=subrenv), - H.ops, lenvs, renvs)) -end - # we need to define how to recalculate """ Recalculate in-place each sub-env in MultipleEnvironments diff --git a/src/environments/windowenv.jl b/src/environments/windowenv.jl index 66afda5bd..a68f70df4 100644 --- a/src/environments/windowenv.jl +++ b/src/environments/windowenv.jl @@ -10,22 +10,14 @@ struct WindowEnv{A,B,C,D} <: Cache rinfdeps::PeriodicArray{D} end -const WindowEnvUnion = Union{T,Window{A,T,B}} where {T<:WindowEnv,A,B} - #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::Union{SparseMPO,MPOHamiltonian,DenseMPO}, above=nothing; lenvs=environments(ψ.left, O),renvs=environments(ψ.right, O)) - fin_env = environments(ψ, O, above, leftenv(lenvs, 1, ψ.left), +function environments(ψ::WindowMPS, O::Window, above=nothing; lenvs=environments(ψ.left, O.left),renvs=environments(ψ.right, O.right)) + fin_env = environments(ψ, O.middle, above, leftenv(lenvs, 1, ψ.left), rightenv(renvs, length(ψ), ψ.right)) return WindowEnv(Window(lenvs,fin_env,renvs),copy(ψ.left.AL),copy(ψ.right.AR)) end -# when supplied with a Window Hamiltonian, we cannot assume H.left/H.right is the same as H.middle -function environments(ψ::WindowMPS, O::Window, above=nothing; lenvs=environments(ψ.left, O.left),renvs=environments(ψ.right, O.right)) - window_env = Window(environments(ψ.left, O.middle),environments(ψ.middle, O.middle),environments(ψ.left, O.middle)) - return Window(lenvs,WindowEnv(window_env,copy(ψ.left.AL),copy(ψ.right.AR)),renvs) -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")) @@ -38,7 +30,7 @@ end #=========================================================================================== Utility ===========================================================================================# -function Base.getproperty(ca::WindowEnv,sym::Symbol) +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? @@ -48,13 +40,35 @@ function Base.getproperty(ca::WindowEnv,sym::Symbol) end end +#to be tested +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 + +function Base.getproperty(ca::MultipleEnvironments{<:LazySum,<:WindowEnv},sym::Symbol) + if sym === :left || 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 + # 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 -finenv(ca::Window{A,<:WindowEnv,B},ψ::WindowMPS) where {A,B} = finenv(ca.middle,ψ) +#to be tested +function finenv(ca::MultipleEnvironments{<:WindowEnv},ψ::WindowMPS) + return MultipleEnvironments(ca.opp.middle,map(x->finenv(x.middle,ψ),ca)) +end #notify the cache that we updated in-place, so it should invalidate the dependencies invalidate!(ca::WindowEnv, ind) = invalidate!(ca.middle,ind) @@ -84,14 +98,13 @@ function check_leftinfenv!(ca::WindowEnv, ψ::WindowMPS) end end -#replace ψ.left by appropriate site index function glue_left!(ψ::WindowMPS,oldALs) #do we need C of finite part too? newD = left_virtualspace(ψ.left, 0) oldD = left_virtualspace(ψ, 0) - #if newD == oldD - # return nothing - #end + if newD == oldD + return nothing + end v = TensorMap(rand, ComplexF64, newD, oldD) (vals, vecs) = eigsolve( flip(TransferMatrix(oldALs, ψ.left.AL)), v, 1, :LM, Arnoldi() @@ -100,13 +113,12 @@ function glue_left!(ψ::WindowMPS,oldALs) ψ.AC[1] = _transpose_front(normalize!(rho * ψ.CR[0]) * _transpose_tail(ψ.AR[1])) end -#space mismatch here still function glue_right!(ψ::WindowMPS,oldARs) newD = right_virtualspace(ψ.right, 0) - oldD = right_virtualspace(ψ, 0) - #if newD == oldD - # return nothing - #end + oldD = right_virtualspace(ψ, length(ψ)) + if newD == oldD + return nothing + end v = TensorMap(rand, ComplexF64, oldD, newD) (vals, vecs) = eigsolve( TransferMatrix(oldARs, ψ.right.AR), v, 1, :LM, Arnoldi() @@ -130,29 +142,30 @@ left_of_finenv(ca::WindowEnv) = ca.left right_of_finenv(ca::WindowEnv) = ca.right left_of_finenv(ca::Window{A,<:WindowEnv,B}) where {A,B} = ca.middle.left right_of_finenv(ca::Window{A,<:WindowEnv,B}) where {A,B} = ca.middle.right +=# -function leftenv(ca::WindowEnvUnion, ind, ψ::WindowMPS) +# automatic recalculation +function leftenv(ca::WindowEnv, ind, ψ::WindowMPS) if ind < 1 - return leftenv(left_of_finenv(ca),ind,ψ.left) + return leftenv(ca.left,ind,ψ.left) elseif ind > length(ψ) - return leftenv(right_of_finenv(ca),ind,ψ.right) + return leftenv(ca.right,ind,ψ.right) else return leftenv(finenv(ca,ψ),ind,ψ) end end -function rightenv(ca::WindowEnvUnion, ind, ψ::WindowMPS) +function rightenv(ca::WindowEnv, ind, ψ::WindowMPS) if ind < 1 - return rightenv(left_of_finenv(ca),ind,ψ.left) + return rightenv(ca.left,ind,ψ.left) elseif ind > length(ψ) - return rightenv(right_of_finenv(ca),ind,ψ.right) + return rightenv(ca.right,ind,ψ.right) else return rightenv(finenv(ca,ψ),ind,ψ) end end -=# -# to be moved +# remove optional env like the rest and use doorgever function expectation_value(Ψ::WindowMPS, windowH::Window, windowEnvs::Window{A,<:WindowEnv,B}=environments(Ψ, windowH)) where {A,B} left = expectation_value(Ψ.left, windowH.left, windowEnvs.left) middle = expectation_value(Ψ.middle, windowH.middle, finenv(windowEnvs,Ψ)) @@ -160,7 +173,7 @@ function expectation_value(Ψ::WindowMPS, windowH::Window, windowEnvs::Window{A, return left,middle,right end -function expectation_value(Ψ::WindowMPS, H, windowEnvs::WindowEnv=environments(Ψ, H)) +function expectation_value(Ψ::WindowMPS, H, windowEnvs::WindowEnv) left = expectation_value(Ψ.left, H, windowEnvs.left) middle = expectation_value(Ψ.middle, H, finenv(windowEnvs,Ψ)) right = expectation_value(Ψ.right, H, windowEnvs.right) @@ -169,6 +182,7 @@ end # do we need this? I find it convenvient +# currently broken though function forget_infinite(ψ::WindowMPS,env::WindowEnvUnion) fin_env = finenv(env,ψ) newenv = environments(ψ.middle,fin_env.opp,fin_env.above,copy(leftenv(env, 1, ψ)),copy(rightenv(env, length(ψ), ψ))) From f8c307166d40f36da5203e02e32627d9cc28ca1a Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Wed, 21 Feb 2024 21:34:24 +0100 Subject: [PATCH 15/62] constructor for MultipliedOperator{Window} --- src/operators/multipliedoperator.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/operators/multipliedoperator.jl b/src/operators/multipliedoperator.jl index 6f9100e69..e51b474c3 100644 --- a/src/operators/multipliedoperator.jl +++ b/src/operators/multipliedoperator.jl @@ -50,3 +50,7 @@ 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 From 8938f1bf150c860f71784a36f36810e626c169ae Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Wed, 21 Feb 2024 22:56:15 +0100 Subject: [PATCH 16/62] refactoring of tdvp sweep --- src/MPSKit.jl | 1 + src/algorithms/timestep/tdvp.jl | 38 +++--- src/algorithms/timestep/timestep.jl | 24 ++++ src/algorithms/timestep/windowtdvp.jl | 181 ++++---------------------- 4 files changed, 71 insertions(+), 173 deletions(-) create mode 100644 src/algorithms/timestep/timestep.jl diff --git a/src/MPSKit.jl b/src/MPSKit.jl index 8ec314ffb..3046decdb 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -119,6 +119,7 @@ include("algorithms/changebonds/randexpand.jl") include("algorithms/timestep/tdvp.jl") include("algorithms/timestep/timeevmpo.jl") include("algorithms/timestep/integrators.jl") +include("algorithms/timestep/timestep.jl") include("algorithms/timestep/time_evolve.jl") include("algorithms/timestep/windowtdvp.jl") diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl index 57c2bd489..6a8e2e110 100644 --- a/src/algorithms/timestep/tdvp.jl +++ b/src/algorithms/timestep/tdvp.jl @@ -58,34 +58,40 @@ function timestep(Ψ::InfiniteMPS, H, t::Number, dt::Number, alg::TDVP, 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) + Ψ.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) + Ψ.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) + Ψ.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) + Ψ.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) + Ψ.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) + Ψ.AC[1] = integrate(h_ac, Ψ.AC[1], t + dt, dt, alg.integrator) return Ψ, envs end @@ -112,8 +118,8 @@ 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) @@ -131,6 +137,12 @@ function timestep!(Ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP2, 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]) @@ -149,9 +161,3 @@ function timestep!(Ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP2, 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 diff --git a/src/algorithms/timestep/timestep.jl b/src/algorithms/timestep/timestep.jl new file mode 100644 index 000000000..22a5c4449 --- /dev/null +++ b/src/algorithms/timestep/timestep.jl @@ -0,0 +1,24 @@ +# 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 + +# 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 index 4365eb0fa..147b1393f 100644 --- a/src/algorithms/timestep/windowtdvp.jl +++ b/src/algorithms/timestep/windowtdvp.jl @@ -1,170 +1,37 @@ -# make a struct for WindowTDVP, finalize for bondexpansion and glue after that - -function timestep(Ψ::WindowMPS{A,B,VL,VR},H::Union{Window,LazySum},t::Number,dt::Number,alg::TDVP,env::Cache=environments(Ψ, H)) where {A,B,VL,VR} - #first evolve left state - if VL === WINDOW_VARIABLE - println("Doing left") - nleft, _ = timestep(Ψ.left, H.left, t, dt, alg, env.left; leftorthflag=true) #env gets updated in place - Ψ = WindowMPS(nleft,Ψ.middle,Ψ.right) - end - - # left to right sweep on window - for i in 1:(length(Ψ) - 1) - h_ac = ∂∂AC(i, Ψ, H.middle, env) - Ψ.AC[i] = integrate(h_ac, Ψ.AC[i], t, dt / 2, alg.integrator) - - h_c = ∂∂C(i, Ψ, H.middle, env) - Ψ.CR[i] = integrate(h_c, Ψ.CR[i], t, -dt / 2, alg.integrator) - end - - h_ac = ∂∂AC(length(Ψ), Ψ, H.middle, env) - Ψ.AC[end] = integrate(h_ac, Ψ.AC[end], t, dt / 2, alg.integrator) - - #then evolve right state - if VR === WINDOW_VARIABLE - println("Doing right") - nright, _ = timestep(Ψ.right, H.right, t + dt, dt, alg, env.right; leftorthflag=false) # env gets updated in place - Ψ = WindowMPS(Ψ.left,Ψ.middle,nright) - end - - # right to left sweep on window - for i in length(Ψ):-1:2 - h_ac = ∂∂AC(i, Ψ, H.middle, env) - Ψ.AC[i] = integrate(h_ac, Ψ.AC[i], t + dt / 2, dt / 2, alg.integrator) - - h_c = ∂∂C(i - 1, Ψ, H.middle, env) - Ψ.CR[i - 1] = integrate(h_c, Ψ.CR[i - 1], t + dt / 2, -dt / 2, alg.integrator) - end - - h_ac = ∂∂AC(1, Ψ, H.middle, env) - Ψ.AC[1] = integrate(h_ac, Ψ.AC[1], t, dt / 2, alg.integrator) - - return Ψ, env +""" + WindowTDVP{A} <: Algorithm + +[Mixed TDVP](https://arxiv.org/abs/2007.15035) algorithm for time evolution. + +# Fields +- `finite_alg::A`: algorithm to do the timestep of the finite part of the WindowMPS. This can be `TDVP2` to expand the bonddimension. +- `infinite_alg::A`: algorithm to do the timestep of the infinite part of the WindowMPS +- `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,F} <: Algorithm + finite_alg::A = TDVP() + infinite_alg::B = TDVP() + finalize::F = Defaults._finalize end -#= - -function timestep!( - Ψ::WindowMPS{A,B,VL,VR}, - H, - t::Number, - dt::Number, - alg::TDVP, - env=environments(Ψ, H); -) where {A,B,VL,VR} +function timestep(Ψ::WindowMPS{A,B,VL,VR},H::Union{Window,LazySum{<:Window}},t::Number,dt::Number,alg::WindowTDVP,env::Cache=environments(Ψ, H)) where {A,B,VL,VR} + #first evolve left state if VL === WINDOW_VARIABLE - nleft, _ = timestep(Ψ.left, H.left, t, dt, alg, env.left; leftorthflag=true) #env gets updated in place + nleft, _ = timestep(Ψ.left, H.left, t, dt, alg.infinite_alg, env.left; leftorthflag=true) #env gets updated in place Ψ = WindowMPS(nleft,Ψ.middle,Ψ.right) end - # left to right sweep on window - for i in 1:(length(Ψ) - 1) - h_ac = ∂∂AC(i, Ψ, H.middle, env.middle) - Ψ.AC[i] = integrate(h_ac, Ψ.AC[i], t, dt / 2, alg.integrator) - - h_c = ∂∂C(i, Ψ, H.middle, env.middle) - Ψ.CR[i] = integrate(h_c, Ψ.CR[i], t, -dt / 2, alg.integrator) - end - - h_ac = ∂∂AC(length(Ψ), Ψ, H.middle, env.middle) - Ψ.AC[end] = integrate(h_ac, Ψ.AC[end], t, dt / 2, alg.integrator) + Ψ, env = ltr_sweep!(Ψ, H.middle, t, dt / 2, alg.finite_alg, env) + #then evolve right state if VR === WINDOW_VARIABLE - nright, _ = timestep(Ψ.right_gs, H.right, t + dt, dt, alg, env.right; leftorthflag=false) # env gets updated in place + nright, _ = timestep(Ψ.right, H.right, t + dt, dt, alg.infinite_alg, env.right; leftorthflag=false) # env gets updated in place Ψ = WindowMPS(Ψ.left,Ψ.middle,nright) end - # right to left sweep on window - for i in length(Ψ):-1:2 - h_ac = ∂∂AC(i, Ψ, H.middle, env.middle) - Ψ.AC[i] = integrate(h_ac, Ψ.AC[i], t + dt / 2, dt / 2, alg.integrator) - - h_c = ∂∂C(i - 1, Ψ, H.middle, env.middle) - Ψ.CR[i - 1] = integrate(h_c, Ψ.CR[i - 1], t + dt / 2, -dt / 2, alg.integrator) - end - - h_ac = ∂∂AC(1, Ψ, H.middle, env.middle) - Ψ.AC[1] = integrate(h_ac, Ψ.AC[1], t, dt / 2, alg.integrator) - + Ψ, env = rtl_sweep!(Ψ, H.middle, t + dt / 2, dt / 2, alg.finite_alg, env) + return Ψ, env -end -=# - -#= -function timestep!( - Ψ::WindowMPS, - H::Window, - t::Number, - dt::Number, - alg::TDVP2, - env::Window=environments(Ψ, H); - leftevolve=false, - rightevolve=false, - kwargs..., -) - singleTDVPalg = TDVP(; - integrator=alg.integrator, tolgauge=alg.tolgauge, maxiter=alg.maxiter - ) - - # first evolve left state - if leftevolve - # expand the bond dimension using changebonds - nleft, _ = leftexpand(Ψ, H.left(t), env.left; kwargs...) - # fill it by doing regular TDVP - nleft, _ = timestep( - nleft, H.left, t, dt, singleTDVPalg, env.left; leftorthflag=true - ) - _update_leftEnv!(nleft, env) - else - nleft = Ψ.left - end - - # left to right sweep on window - for i in 1:(length(Ψ) - 1) - h_ac2 = ∂∂AC2(i, Ψ, H.middle, env.middle) - ac2 = Ψ.AC[i] * _transpose_tail(Ψ.AR[i + 1]) - ac2 = integrate(h_ac2, ac2, t, dt / 2, alg.integrator) - - U, S, V, = tsvd!(ac2; alg=TensorKit.SVD(), trunc=alg.trscheme) - - Ψ.AC[i] = (U, S) - Ψ.AC[i + 1] = (S, _transpose_front(V)) - - if i < length(Ψ) - 1 - h_ac = ∂∂AC(i + 1, Ψ, H.middle, env.middle) - Ψ.AC[i + 1] = integrate(h_ac, Ψ.AC[i + 1], t, -dt / 2, alg.integrator) - end - end - - if rightevolve - # expand the bond dimension using changebonds - nright, _ = rightexpand(Ψ, H.right(t), env.right; kwargs...) - # fill it by doing regular TDVP - nright, _ = timestep( - nright, H.right, t + dt, dt, singleTDVPalg, env.right; leftorthflag=false - ) #env gets updated in place - _update_rightEnv!(nright, env) - else - nright = Ψ.right_gs - end - - # right to left sweep on window - for i in length(Ψ):-1:2 - h_ac2 = ∂∂AC2(i - 1, Ψ, H.middle, env.middle) - ac2 = Ψ.AL[i - 1] * _transpose_tail(Ψ.AC[i]) - ac2 = integrate(h_ac2, ac2, t + dt / 2, dt / 2, alg.integrator) - U, S, V, = tsvd!(ac2; alg=TensorKit.SVD(), trunc=alg.trscheme) - - Ψ.AC[i - 1] = (U, S) - Ψ.AC[i] = (S, _transpose_front(V)) - - if i > 2 - h_ac = ∂∂AC(i - 1, Ψ, H.middle, env.middle) - Ψ.AC[i - 1] = integrate(h_ac, Ψ.AC[i - 1], t + dt / 2, -dt / 2, alg.integrator) - end - end - - return WindowMPS(nleft, Ψ.window, nright), env -end -=# +end \ No newline at end of file From 06c856c2ac2c59d5b4afaabd0acf0c3c6f3ebb6e Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Wed, 21 Feb 2024 23:10:37 +0100 Subject: [PATCH 17/62] fixed wrong typing in windowtdvp.jl --- src/algorithms/timestep/windowtdvp.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algorithms/timestep/windowtdvp.jl b/src/algorithms/timestep/windowtdvp.jl index 147b1393f..ee2a71573 100644 --- a/src/algorithms/timestep/windowtdvp.jl +++ b/src/algorithms/timestep/windowtdvp.jl @@ -15,7 +15,7 @@ finalize::F = Defaults._finalize end -function timestep(Ψ::WindowMPS{A,B,VL,VR},H::Union{Window,LazySum{<:Window}},t::Number,dt::Number,alg::WindowTDVP,env::Cache=environments(Ψ, H)) where {A,B,VL,VR} +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 From 556b1b0adba7892b83809790a2b2437282aa8bff Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Wed, 21 Feb 2024 23:13:22 +0100 Subject: [PATCH 18/62] cleanup of windowenv.jl --- src/environments/windowenv.jl | 85 ++++------------------------------- 1 file changed, 9 insertions(+), 76 deletions(-) diff --git a/src/environments/windowenv.jl b/src/environments/windowenv.jl index a68f70df4..c0e9227ed 100644 --- a/src/environments/windowenv.jl +++ b/src/environments/windowenv.jl @@ -40,7 +40,6 @@ function Base.getproperty(ca::WindowEnv,sym::Symbol) #under review end end -#to be tested function Base.getproperty(sumops::LazySum{<:Window},sym::Symbol) if sym === :left || sym === :middle || sym === :right #extract the left/right parts @@ -65,7 +64,7 @@ function finenv(ca::WindowEnv,ψ::WindowMPS{A,B,VL,VR}) where {A,B,VL,VR} VR === WINDOW_FIXED || check_rightinfenv!(ca,ψ) return ca.middle end -#to be tested + function finenv(ca::MultipleEnvironments{<:WindowEnv},ψ::WindowMPS) return MultipleEnvironments(ca.opp.middle,map(x->finenv(x.middle,ψ),ca)) end @@ -80,8 +79,10 @@ function check_rightinfenv!(ca::WindowEnv, ψ::WindowMPS) println("changing right env") invalidate!(ca, length(ψ.middle)) #forces transfers to be recalculated lazily - update_rightstart!(ca.middle,ca.right,ψ.right) - glue_right!(ψ,ca.rinfdeps) + #update_rightstart!(ca.middle,ca.right,ψ.right) + #glue_right!(ψ,ca.rinfdeps) + + ca_fin.rightenvs[end] = rightenv(ca_infin, 0, ψ) ca.rinfdeps .= ψ.right.AR end end @@ -92,58 +93,14 @@ function check_leftinfenv!(ca::WindowEnv, ψ::WindowMPS) println("changing left env") invalidate!(ca, 1) #forces transfers to be recalculated lazily - update_leftstart!(ca.middle,ca.left,ψ.left) - glue_left!(ψ,ca.linfdeps) - ca.linfdeps .= ψ.left.AL - end -end + #update_leftstart!(ca.middle,ca.left,ψ.left) + #glue_left!(ψ,ca.linfdeps) -function glue_left!(ψ::WindowMPS,oldALs) - #do we need C of finite part too? - newD = left_virtualspace(ψ.left, 0) - oldD = left_virtualspace(ψ, 0) - if newD == oldD - return nothing - end - v = TensorMap(rand, ComplexF64, newD, oldD) - (vals, vecs) = eigsolve( - flip(TransferMatrix(oldALs, ψ.left.AL)), v, 1, :LM, Arnoldi() - ) - rho = pinv(ψ.left.CR[0]) * vecs[1] * ψ.CR[0] #CR[0] == CL[1] - ψ.AC[1] = _transpose_front(normalize!(rho * ψ.CR[0]) * _transpose_tail(ψ.AR[1])) -end - -function glue_right!(ψ::WindowMPS,oldARs) - newD = right_virtualspace(ψ.right, 0) - oldD = right_virtualspace(ψ, length(ψ)) - if newD == oldD - return nothing + ca_fin.leftenvs[1] = leftenv(ca_infin, length(ψ)+1, ψ) + ca.linfdeps .= ψ.left.AL end - v = TensorMap(rand, ComplexF64, oldD, newD) - (vals, vecs) = eigsolve( - TransferMatrix(oldARs, ψ.right.AR), v, 1, :LM, Arnoldi() - ) - rho = ψ.CR[end] * vecs[1] * pinv(ψ.right.CR[0]) - ψ.AC[end] = ψ.AL[end] * normalize!(ψ.CR[end] * rho) -end - -# these are to be extended for LazySum and the like -function update_leftstart!(ca_fin::FinEnv,ca_infin::MPOHamInfEnv, ψ::InfiniteMPS) - ca_fin.leftenvs[1] = leftenv(ca_infin, length(ψ)+1, ψ) end -function update_rightstart!(ca_fin::FinEnv,ca_infin::MPOHamInfEnv, ψ::InfiniteMPS) - ca_fin.rightenvs[end] = rightenv(ca_infin, 0, ψ) -end - -# under review -#= -left_of_finenv(ca::WindowEnv) = ca.left -right_of_finenv(ca::WindowEnv) = ca.right -left_of_finenv(ca::Window{A,<:WindowEnv,B}) where {A,B} = ca.middle.left -right_of_finenv(ca::Window{A,<:WindowEnv,B}) where {A,B} = ca.middle.right -=# - # automatic recalculation function leftenv(ca::WindowEnv, ind, ψ::WindowMPS) if ind < 1 @@ -164,27 +121,3 @@ function rightenv(ca::WindowEnv, ind, ψ::WindowMPS) return rightenv(finenv(ca,ψ),ind,ψ) end end - -# remove optional env like the rest and use doorgever -function expectation_value(Ψ::WindowMPS, windowH::Window, windowEnvs::Window{A,<:WindowEnv,B}=environments(Ψ, windowH)) where {A,B} - left = expectation_value(Ψ.left, windowH.left, windowEnvs.left) - middle = expectation_value(Ψ.middle, windowH.middle, finenv(windowEnvs,Ψ)) - right = expectation_value(Ψ.right, windowH.right, windowEnvs.right) - return left,middle,right -end - -function expectation_value(Ψ::WindowMPS, H, windowEnvs::WindowEnv) - left = expectation_value(Ψ.left, H, windowEnvs.left) - middle = expectation_value(Ψ.middle, H, finenv(windowEnvs,Ψ)) - right = expectation_value(Ψ.right, H, windowEnvs.right) - return left,middle,right -end - - -# do we need this? I find it convenvient -# currently broken though -function forget_infinite(ψ::WindowMPS,env::WindowEnvUnion) - fin_env = finenv(env,ψ) - newenv = environments(ψ.middle,fin_env.opp,fin_env.above,copy(leftenv(env, 1, ψ)),copy(rightenv(env, length(ψ), ψ))) - return copy(ψ.middle),newenv -end From 7f305d548e322331c48bf6367f6309c2e1c37a24 Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Wed, 21 Feb 2024 23:15:29 +0100 Subject: [PATCH 19/62] changed windowmps related expvals --- src/algorithms/expval.jl | 33 +++++++++------------------------ 1 file changed, 9 insertions(+), 24 deletions(-) diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl index c3e8f1bad..c61b10262 100644 --- a/src/algorithms/expval.jl +++ b/src/algorithms/expval.jl @@ -64,32 +64,11 @@ end # -------------- # TODO: add section in documentation to explain convention expectation_value(ψ, envs::Cache) = expectation_value(ψ, envs.opp, envs) -function expectation_value(ψ, H::MPOHamiltonian) +function expectation_value(ψ, H::Union{MPOHamiltonian,Window,LazySum}) return expectation_value(ψ, H, environments(ψ, H)) end -""" - expectation_value(ψ::WindowMPS, H::MPOHAmiltonian, envs) -> vals, tot - -TODO -""" -function expectation_value(ψ::WindowMPS, H::MPOHamiltonian, envs::FinEnv) - vals = expectation_value_fimpl(ψ, H, envs) - - tot = 0.0 + 0im - for i in 1:(H.odim), j in 1:(H.odim) - tot += @plansor leftenv(envs, length(ψ), ψ)[i][1 2; 3] * ψ.AC[end][3 4; 5] * - rightenv(envs, length(ψ), ψ)[j][5 6; 7] * - H[length(ψ)][i, j][2 8; 4 6] * conj(ψ.AC[end][1 8; 7]) - end - - return vals, tot / (norm(ψ.AC[end])^2) -end - -function expectation_value(ψ::FiniteMPS, H::MPOHamiltonian, envs::FinEnv) - return expectation_value_fimpl(ψ, H, envs) -end -function expectation_value_fimpl(ψ::AbstractFiniteMPS, H::MPOHamiltonian, envs::FinEnv) +function expectation_value(ψ::AbstractFiniteMPS, H::MPOHamiltonian, envs::FinEnv) ens = zeros(scalartype(ψ), length(ψ)) for i in 1:length(ψ), (j, k) in keys(H[i]) !((j == 1 && k != 1) || (k == H.odim && j != H.odim)) && continue @@ -187,7 +166,7 @@ end function expectation_value(ψ, ops::LazySum, at::Int) return sum(op -> expectation_value(ψ, op, at), ops) end -function expectation_value(ψ, ops::LazySum, envs::MultipleEnvironments=environments(ψ, ops)) +function expectation_value(ψ, ops::LazySum, envs::MultipleEnvironments) return sum(((op, env),) -> expectation_value(ψ, op, env), zip(ops.ops, envs)) end @@ -210,3 +189,9 @@ function expectation_value(ψ::FiniteMPS, O::ProjectionOperator, n = norm(ψ.AC[end])^2 return ens ./ n end + +# Window +# ------ +function expectation_value(Ψ::WindowMPS, windowH::Union{Window,LazySum{<:Window}}, windowenvs) = + return expectation_value(Ψ, windowH.middle, windowenvs.middle) +end \ No newline at end of file From e98416a35399a42b873e3f172a12ac52d7e0a321 Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Thu, 22 Feb 2024 12:58:51 +0100 Subject: [PATCH 20/62] now it will precompile --- src/algorithms/expval.jl | 2 +- src/environments/multipleenv.jl | 36 +++++++++++++++++++++++++-------- src/environments/windowenv.jl | 22 -------------------- src/operators/lazysum.jl | 9 +++++++++ src/states/windowmps.jl | 20 +++++++++--------- 5 files changed, 48 insertions(+), 41 deletions(-) diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl index c61b10262..3164694ea 100644 --- a/src/algorithms/expval.jl +++ b/src/algorithms/expval.jl @@ -192,6 +192,6 @@ end # Window # ------ -function expectation_value(Ψ::WindowMPS, windowH::Union{Window,LazySum{<:Window}}, windowenvs) = +function expectation_value(Ψ::WindowMPS, windowH::Union{Window,LazySum{<:Window}}, windowenvs) return expectation_value(Ψ, windowH.middle, windowenvs.middle) end \ No newline at end of file diff --git a/src/environments/multipleenv.jl b/src/environments/multipleenv.jl index bb125e2b6..b615371f5 100644 --- a/src/environments/multipleenv.jl +++ b/src/environments/multipleenv.jl @@ -30,6 +30,31 @@ end # return MultipleEnvironments(H, broadcast(o -> environments(state, o), H.opps)) # end; + +#=========================================================================================== +Utility +===========================================================================================# +function Base.getproperty(ca::MultipleEnvironments{<:LazySum,<:WindowEnv},sym::Symbol) + if sym === :left || 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{<:WindowEnv},ψ::WindowMPS) + return MultipleEnvironments(ca.opp.middle,map(x->finenv(x.middle,ψ),ca)) +end + # we need to define how to recalculate """ Recalculate in-place each sub-env in MultipleEnvironments @@ -41,11 +66,6 @@ function recalculate!(env::MultipleEnvironments, args...; kwargs...) 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 index c0e9227ed..ec428ffb6 100644 --- a/src/environments/windowenv.jl +++ b/src/environments/windowenv.jl @@ -40,24 +40,6 @@ function Base.getproperty(ca::WindowEnv,sym::Symbol) #under review end end -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 - -function Base.getproperty(ca::MultipleEnvironments{<:LazySum,<:WindowEnv},sym::Symbol) - if sym === :left || 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 - # 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,ψ) @@ -65,10 +47,6 @@ function finenv(ca::WindowEnv,ψ::WindowMPS{A,B,VL,VR}) where {A,B,VL,VR} return ca.middle end -function finenv(ca::MultipleEnvironments{<:WindowEnv},ψ::WindowMPS) - return MultipleEnvironments(ca.opp.middle,map(x->finenv(x.middle,ψ),ca)) -end - #notify the cache that we updated in-place, so it should invalidate the dependencies invalidate!(ca::WindowEnv, ind) = invalidate!(ca.middle,ind) diff --git a/src/operators/lazysum.jl b/src/operators/lazysum.jl index 5d911cbf4..ad5602e39 100644 --- a/src/operators/lazysum.jl +++ b/src/operators/lazysum.jl @@ -57,3 +57,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/states/windowmps.jl b/src/states/windowmps.jl index 1ae63c650..7fda4f892 100644 --- a/src/states/windowmps.jl +++ b/src/states/windowmps.jl @@ -57,31 +57,31 @@ end Constructors ===========================================================================================# function WindowMPS(ψₗ::InfiniteMPS, site_tensors::AbstractVector{<:GenericMPSTensor}, - ψᵣ::InfiniteMPS) - return WindowMPS(ψₗ, FiniteMPS(site_tensors), ψᵣ) + ψᵣ::InfiniteMPS; kwargs...) + return WindowMPS(ψₗ, FiniteMPS(site_tensors), ψᵣ; kwargs...) end 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...) From 97e54386c1fbce064b8a89f96566894a73dcd760 Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Thu, 22 Feb 2024 14:50:20 +0100 Subject: [PATCH 21/62] now timestep actually works --- src/MPSKit.jl | 4 ++-- src/algorithms/timestep/windowtdvp.jl | 2 +- src/environments/windowenv.jl | 4 ++-- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/MPSKit.jl b/src/MPSKit.jl index 3046decdb..9d61fee62 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -40,7 +40,7 @@ export VUMPS, DMRG, DMRG2, IDMRG1, IDMRG2, GradientGrassmann export excitations, FiniteExcited, QuasiparticleAnsatz export marek_gap, correlation_length, correlator export time_evolve, timestep!, timestep -export TDVP, TDVP2, make_time_mpo, WI, WII, TaylorCluster +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 @@ -119,9 +119,9 @@ include("algorithms/changebonds/randexpand.jl") include("algorithms/timestep/tdvp.jl") include("algorithms/timestep/timeevmpo.jl") include("algorithms/timestep/integrators.jl") -include("algorithms/timestep/timestep.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/timestep/windowtdvp.jl b/src/algorithms/timestep/windowtdvp.jl index ee2a71573..f0c66ec2d 100644 --- a/src/algorithms/timestep/windowtdvp.jl +++ b/src/algorithms/timestep/windowtdvp.jl @@ -15,7 +15,7 @@ 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} +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 diff --git a/src/environments/windowenv.jl b/src/environments/windowenv.jl index ec428ffb6..58134a7a9 100644 --- a/src/environments/windowenv.jl +++ b/src/environments/windowenv.jl @@ -60,7 +60,7 @@ function check_rightinfenv!(ca::WindowEnv, ψ::WindowMPS) #update_rightstart!(ca.middle,ca.right,ψ.right) #glue_right!(ψ,ca.rinfdeps) - ca_fin.rightenvs[end] = rightenv(ca_infin, 0, ψ) + ca.middle.rightenvs[end] = rightenv(ca.right, 0, ψ.right) ca.rinfdeps .= ψ.right.AR end end @@ -74,7 +74,7 @@ function check_leftinfenv!(ca::WindowEnv, ψ::WindowMPS) #update_leftstart!(ca.middle,ca.left,ψ.left) #glue_left!(ψ,ca.linfdeps) - ca_fin.leftenvs[1] = leftenv(ca_infin, length(ψ)+1, ψ) + ca.middle.leftenvs[1] = leftenv(ca.left, length(ψ.left)+1, ψ.left) ca.linfdeps .= ψ.left.AL end end From 7d0adc591801052e78e3e7a67f43127af25c1225 Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Fri, 23 Feb 2024 16:25:34 +0100 Subject: [PATCH 22/62] time_evolve! was not exported --- src/MPSKit.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/MPSKit.jl b/src/MPSKit.jl index 9d61fee62..a7d1f477c 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -39,7 +39,7 @@ export find_groundstate!, find_groundstate, leading_boundary export VUMPS, DMRG, DMRG2, IDMRG1, IDMRG2, GradientGrassmann export excitations, FiniteExcited, QuasiparticleAnsatz export marek_gap, correlation_length, correlator -export time_evolve, timestep!, timestep +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 From 02ec8c1ba6c324f6e7047a609f0248e44eb90427 Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Mon, 26 Feb 2024 16:21:45 +0100 Subject: [PATCH 23/62] cleanup --- src/environments/windowenv.jl | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/src/environments/windowenv.jl b/src/environments/windowenv.jl index 58134a7a9..8c5346840 100644 --- a/src/environments/windowenv.jl +++ b/src/environments/windowenv.jl @@ -52,28 +52,18 @@ invalidate!(ca::WindowEnv, ind) = invalidate!(ca.middle,ind) # Check the infinite envs and recalculate the right start function check_rightinfenv!(ca::WindowEnv, ψ::WindowMPS) - println("Doing right check") if !all(ca.rinfdeps .=== ψ.right.AR) - println("changing right env") invalidate!(ca, length(ψ.middle)) #forces transfers to be recalculated lazily - #update_rightstart!(ca.middle,ca.right,ψ.right) - #glue_right!(ψ,ca.rinfdeps) - ca.middle.rightenvs[end] = rightenv(ca.right, 0, ψ.right) ca.rinfdeps .= ψ.right.AR end end function check_leftinfenv!(ca::WindowEnv, ψ::WindowMPS) - println("Doing left check") if !all(ca.linfdeps .=== ψ.left.AL) - println("changing left env") invalidate!(ca, 1) #forces transfers to be recalculated lazily - #update_leftstart!(ca.middle,ca.left,ψ.left) - #glue_left!(ψ,ca.linfdeps) - ca.middle.leftenvs[1] = leftenv(ca.left, length(ψ.left)+1, ψ.left) ca.linfdeps .= ψ.left.AL end From 8b4737a57fad232a97c834ea1292af09f06c5710 Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Mon, 26 Feb 2024 16:22:04 +0100 Subject: [PATCH 24/62] cleanup --- src/environments/multipleenv.jl | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/environments/multipleenv.jl b/src/environments/multipleenv.jl index b615371f5..48370b817 100644 --- a/src/environments/multipleenv.jl +++ b/src/environments/multipleenv.jl @@ -35,7 +35,7 @@ end Utility ===========================================================================================# function Base.getproperty(ca::MultipleEnvironments{<:LazySum,<:WindowEnv},sym::Symbol) - if sym === :left || sym === :right + 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 @@ -51,10 +51,6 @@ function Base.getproperty(envs::MultipleEnvironments, prop::Symbol) end end -function finenv(ca::MultipleEnvironments{<:WindowEnv},ψ::WindowMPS) - return MultipleEnvironments(ca.opp.middle,map(x->finenv(x.middle,ψ),ca)) -end - # we need to define how to recalculate """ Recalculate in-place each sub-env in MultipleEnvironments From 58e51927c901a91207d3b38801e7079bad6e6584 Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Mon, 26 Feb 2024 16:22:42 +0100 Subject: [PATCH 25/62] Update windowmps.jl --- src/states/windowmps.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/states/windowmps.jl b/src/states/windowmps.jl index 7fda4f892..4c582663e 100644 --- a/src/states/windowmps.jl +++ b/src/states/windowmps.jl @@ -35,7 +35,7 @@ Construct a WindowMPS from an InfiniteMPS, by promoting a region of length `L` t `FiniteMPS`. Options for fixing the left and right infinite state (i.e. so they don't get time evolved) -can be done via the keyword arguments `fixleft` and `fixright`. +can be done via the Boolean keyword arguments `fixleft` and `fixright`. """ struct WindowMPS{A<:GenericMPSTensor,B<:MPSBondTensor,VL,VR} <: AbstractFiniteMPS window::Window{InfiniteMPS{A,B},FiniteMPS{A,B},InfiniteMPS{A,B}} @@ -53,6 +53,8 @@ struct WindowMPS{A<:GenericMPSTensor,B<:MPSBondTensor,VL,VR} <: AbstractFiniteMP end + + #=========================================================================================== Constructors ===========================================================================================# From b6b0953dc5b8564f159e8825511fba41db57a8e1 Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Mon, 26 Feb 2024 16:23:09 +0100 Subject: [PATCH 26/62] new timestep for WindowMPS --- src/algorithms/timestep/timestep.jl | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/algorithms/timestep/timestep.jl b/src/algorithms/timestep/timestep.jl index 22a5c4449..da134fe76 100644 --- a/src/algorithms/timestep/timestep.jl +++ b/src/algorithms/timestep/timestep.jl @@ -1,6 +1,6 @@ # TDVP -function timestep!(Ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::Union{TDVP,TDVP2}, +function timestep!(Ψ::FiniteMPS, H, t::Number, dt::Number, alg::Union{TDVP,TDVP2}, envs::Union{Cache,MultipleEnvironments}=environments(Ψ, H)) Ψ,envs = ltr_sweep!(Ψ, H, t, dt / 2, alg, envs) @@ -10,12 +10,18 @@ function timestep!(Ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::Union{T end #copying version -function timestep(Ψ::AbstractFiniteMPS, H, time::Number, timestep::Number, +function timestep(Ψ::FiniteMPS, 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, From f2f4774592bf6571b1ff94ae3ed7221020cb4437 Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Mon, 26 Feb 2024 17:00:30 +0100 Subject: [PATCH 27/62] better dispatching for window in expval --- src/algorithms/expval.jl | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl index 3164694ea..5055de5ab 100644 --- a/src/algorithms/expval.jl +++ b/src/algorithms/expval.jl @@ -162,7 +162,8 @@ function expectation_value(ψ, op::UntimedOperator, args...) return op.f * expectation_value(ψ, op.op, args...) end -# define expectation_value for LazySum +# LazySum +# ------- function expectation_value(ψ, ops::LazySum, at::Int) return sum(op -> expectation_value(ψ, op, at), ops) end @@ -192,6 +193,11 @@ end # Window # ------ -function expectation_value(Ψ::WindowMPS, windowH::Union{Window,LazySum{<:Window}}, windowenvs) - return expectation_value(Ψ, windowH.middle, windowenvs.middle) -end \ No newline at end of file +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 From fd4883f0c80a738d5e095daa9448759007c1e9ab Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Mon, 26 Feb 2024 17:44:45 +0100 Subject: [PATCH 28/62] forgot about time dependent stuff --- src/operators/lazysum.jl | 6 ------ src/operators/timedependence.jl | 8 ++++++++ src/states/window.jl | 16 ++++++++++++---- 3 files changed, 20 insertions(+), 10 deletions(-) diff --git a/src/operators/lazysum.jl b/src/operators/lazysum.jl index ad5602e39..13b8952a7 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 diff --git a/src/operators/timedependence.jl b/src/operators/timedependence.jl index 363b9eb2f..208d39c55 100644 --- a/src/operators/timedependence.jl +++ b/src/operators/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/src/states/window.jl b/src/states/window.jl index 2c7028683..ea5ff49fb 100644 --- a/src/states/window.jl +++ b/src/states/window.jl @@ -1,7 +1,5 @@ -# Note : this is intended to be a template for windowmps and windows of operators/environments - " - Window(leftstate,window,rightstate) + Window(left,middle,right) general struct of an object with a left, middle and right part. " @@ -9,4 +7,14 @@ struct Window{L,M,R} left::L middle::M right::R -end \ No newline at end of file +end + +# Holy traits +TimeDependence(x::Window) = istimed(x) ? TimeDependent() : NotTimeDependent() +istimed(x::Window) = istimed(x.left) || istimed(x.middle) || istimed(x.right) + +_eval_at(x::Window, args...) = Window(_eval_at(x.left,args...),_eval_at(x.middle,args...),_eval_at(x.right,args...)) +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 From 21be93c96e3199b8d5299a85e7e6745aeb0dabba Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Mon, 26 Feb 2024 22:43:09 +0100 Subject: [PATCH 29/62] moved timedependence.jl --- src/{operators => utility}/timedependence.jl | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/{operators => utility}/timedependence.jl (100%) diff --git a/src/operators/timedependence.jl b/src/utility/timedependence.jl similarity index 100% rename from src/operators/timedependence.jl rename to src/utility/timedependence.jl From 512f33626ebe6b98d7b1bb43b3e8509ae9c820db Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Mon, 26 Feb 2024 23:19:54 +0100 Subject: [PATCH 30/62] struggled with variance --- src/algorithms/expval.jl | 12 ++++++++ src/algorithms/propagator/corvector.jl | 40 ++++++++++++++------------ src/algorithms/toolbox.jl | 13 ++++----- 3 files changed, 40 insertions(+), 25 deletions(-) diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl index 5055de5ab..040d60ce8 100644 --- a/src/algorithms/expval.jl +++ b/src/algorithms/expval.jl @@ -130,6 +130,18 @@ function expectation_value(st::InfiniteMPS, H::MPOHamiltonian, range::UnitRange{ return tot end +# maybe rename this so there is no confusion on what it does? +function expectation_value(ψ::WindowMPS, H::MPOHamiltonian, envs::WindowEnv) + tot = 0.0 + 0im + for i in 1:(H.odim), j in 1:(H.odim) + tot += @plansor leftenv(envs, length(ψ), ψ)[i][1 2; 3] * ψ.AC[end][3 4; 5] * + rightenv(envs, length(ψ), ψ)[j][5 6; 7] * + H[length(ψ)][i, j][2 8; 4 6] * conj(ψ.AC[end][1 8; 7]) + end + + return tot / (norm(ψ.AC[end])^2) +end + # Transfer matrices # ----------------- function expectation_value(ψ::InfiniteMPS, mpo::DenseMPO) diff --git a/src/algorithms/propagator/corvector.jl b/src/algorithms/propagator/corvector.jl index 8b5df2676..8f0523591 100644 --- a/src/algorithms/propagator/corvector.jl +++ b/src/algorithms/propagator/corvector.jl @@ -133,19 +133,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 @@ -153,22 +145,34 @@ 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(Ψ::FiniteMPS, H::MPOHamiltonian,envs=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=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/toolbox.jl b/src/algorithms/toolbox.jl index f7da741ad..f0ea188ad 100644 --- a/src/algorithms/toolbox.jl +++ b/src/algorithms/toolbox.jl @@ -144,13 +144,6 @@ function variance(state::FiniteMPS, H::MPOHamiltonian, envs=environments(state, sum(expectation_value(state, H, envs))^2) end -function variance(state::WindowMPS, H::MPOHamiltonian, envs=environments(state, H)) - #tricky to define - H2, nenvs = squaredenvs(state, H, envs) - return real(expectation_value(state, H2, nenvs)[2] - - expectation_value(state, H, envs)[2]^2) -end - function variance(state::FiniteQP, H::MPOHamiltonian, args...) return variance(convert(FiniteMPS, state), H) end; @@ -178,6 +171,12 @@ function variance(ψ, H::LazySum, envs=environments(ψ, sum(H))) return variance(ψ, sum(H), envs) end +function variance(ψ::WindowMPS, H::Window, envs::WindowEnv=environments(ψ,H)) + #tricky to define + H2, nenvs = squaredenvs(ψ, H, envs) + return real(expectation_value(ψ, H2.middle, nenvs) - expectation_value(ψ, H.middle, 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 From b33aa568d7592d83d09c44a4cfc7904aa1386efa Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Mon, 26 Feb 2024 23:20:18 +0100 Subject: [PATCH 31/62] handy constructor for WindowMPS --- src/environments/windowenv.jl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/environments/windowenv.jl b/src/environments/windowenv.jl index 8c5346840..ca058b76d 100644 --- a/src/environments/windowenv.jl +++ b/src/environments/windowenv.jl @@ -13,11 +13,14 @@ 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, leftenv(lenvs, 1, ψ.left), - rightenv(renvs, length(ψ), ψ.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::MPOHamiltonian, above, lenvs::MPOHamInfEnv,renvs::MPOHamInfEnv) + return fin_env = environments(ψ, O, above, leftenv(lenvs, 1, ψ.left),rightenv(renvs, length(ψ), ψ.right)) +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")) From 11321ed27c80f7a8173abf6bcef0859e433bdb03 Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Mon, 26 Feb 2024 23:20:40 +0100 Subject: [PATCH 32/62] timedependence move to utility --- src/MPSKit.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/MPSKit.jl b/src/MPSKit.jl index a7d1f477c..06c3492c8 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -18,7 +18,7 @@ export fix_left,fix_right,fix_infinite 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? @@ -65,6 +65,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/window.jl") @@ -83,7 +84,6 @@ 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") From c348a68b095775f4e00010dc8dd7c3ad9a8b409d Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Mon, 26 Feb 2024 23:20:55 +0100 Subject: [PATCH 33/62] + defined for window --- src/states/window.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/states/window.jl b/src/states/window.jl index ea5ff49fb..853af199e 100644 --- a/src/states/window.jl +++ b/src/states/window.jl @@ -9,6 +9,10 @@ struct Window{L,M,R} 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) From fd157cfb75a817c0b8935ba8db52d9349a388972 Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Tue, 27 Feb 2024 12:56:06 +0100 Subject: [PATCH 34/62] find_groundstate for WindowMPS --- .../groundstate/find_groundstate.jl | 20 ++++++++++++++++++- src/environments/windowenv.jl | 2 +- 2 files changed, 20 insertions(+), 2 deletions(-) diff --git a/src/algorithms/groundstate/find_groundstate.jl b/src/algorithms/groundstate/find_groundstate.jl index 0ce90984b..90cd17cc1 100644 --- a/src/algorithms/groundstate/find_groundstate.jl +++ b/src/algorithms/groundstate/find_groundstate.jl @@ -16,7 +16,7 @@ optimization algorithm will be attempted based on the supplied keywords. - `maxiter::Int`: maximum amount of iterations - `verbose::Bool`: 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, verbose=Defaults.verbose, trscheme=nothing) if isa(ψ, InfiniteMPS) @@ -37,5 +37,23 @@ function find_groundstate(ψ::AbstractMPS, H, envs::Cache=environments(ψ, H); else throw(ArgumentError("Unknown input state type")) end + if isa(ψ, WindowMPS) + alg_infin = VUMPS(; tol_galerkin=tol, verbose=verbose, 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::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) + end + if VR === WINDOW_VARIABLE + (gs_right, _) = find_groundstate(state.right, H.right, alg.right, envs.right) + end + #set infinite parts + state = WindowMPS(gs_left,state.middle,gs_right) + return find_groundstate(state, H.middle, alg.middle, envs) +end diff --git a/src/environments/windowenv.jl b/src/environments/windowenv.jl index ca058b76d..feddbf83e 100644 --- a/src/environments/windowenv.jl +++ b/src/environments/windowenv.jl @@ -18,7 +18,7 @@ function environments(ψ::WindowMPS, O::Window, above=nothing; lenvs=environment end function environments(ψ::WindowMPS, O::MPOHamiltonian, above, lenvs::MPOHamInfEnv,renvs::MPOHamInfEnv) - return fin_env = environments(ψ, O, above, leftenv(lenvs, 1, ψ.left),rightenv(renvs, length(ψ), ψ.right)) + return environments(ψ, O, above, leftenv(lenvs, 1, ψ.left),rightenv(renvs, length(ψ), ψ.right)) end function environments(below::WindowMPS, above::WindowMPS) From 6760faee77e6057c06abd9d696278e54dc019133 Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Tue, 27 Feb 2024 13:01:58 +0100 Subject: [PATCH 35/62] Update find_groundstate.jl --- src/algorithms/groundstate/find_groundstate.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/algorithms/groundstate/find_groundstate.jl b/src/algorithms/groundstate/find_groundstate.jl index 90cd17cc1..cdcc367d6 100644 --- a/src/algorithms/groundstate/find_groundstate.jl +++ b/src/algorithms/groundstate/find_groundstate.jl @@ -45,7 +45,7 @@ function find_groundstate(ψ::AbstractMPS, H, envs::Union{Cache,MultipleEnvironm return find_groundstate(ψ, H, alg, envs) end -function find_groundstate(state::WindowMPS{A,B,VL,VR}, H::Window, alg::Window, envs=environments(state, H)) where {A,B,VL,VR} +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) @@ -57,3 +57,7 @@ function find_groundstate(state::WindowMPS{A,B,VL,VR}, H::Window, alg::Window, e state = WindowMPS(gs_left,state.middle,gs_right) return find_groundstate(state, H.middle, alg.middle, envs) end + +function find_groundstate(ψ::WindowMPS, H::Union{Window,LazySum{<:Window}}, alg::Window, envs...) + return find_groundstate!(copy(ψ), H, alg, envs...) +end From 52e66054d9c370fcf16a853f53c9c94d92432700 Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Tue, 27 Feb 2024 13:12:44 +0100 Subject: [PATCH 36/62] changed comment in windowenv.jl --- src/environments/windowenv.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/environments/windowenv.jl b/src/environments/windowenv.jl index feddbf83e..5cb0f424c 100644 --- a/src/environments/windowenv.jl +++ b/src/environments/windowenv.jl @@ -72,7 +72,7 @@ function check_leftinfenv!(ca::WindowEnv, ψ::WindowMPS) end end -# automatic recalculation +#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) From 5e34180250728638a916c199a3d9a3197464b532 Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Tue, 27 Feb 2024 13:28:19 +0100 Subject: [PATCH 37/62] Update windowtdvp.jl --- src/algorithms/timestep/windowtdvp.jl | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/src/algorithms/timestep/windowtdvp.jl b/src/algorithms/timestep/windowtdvp.jl index f0c66ec2d..64c49e959 100644 --- a/src/algorithms/timestep/windowtdvp.jl +++ b/src/algorithms/timestep/windowtdvp.jl @@ -9,9 +9,10 @@ - `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,F} <: Algorithm - finite_alg::A = TDVP() - infinite_alg::B = TDVP() +@kwdef struct WindowTDVP{A,B,C,F} <: Algorithm + left::A = TDVP() + middle::B = TDVP() + right::C = left finalize::F = Defaults._finalize end @@ -19,19 +20,19 @@ function timestep!(Ψ::WindowMPS{A,B,VL,VR},H::Union{Window,LazySum{<:Window}},t #first evolve left state if VL === WINDOW_VARIABLE - nleft, _ = timestep(Ψ.left, H.left, t, dt, alg.infinite_alg, env.left; leftorthflag=true) #env gets updated in place + 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.finite_alg, env) + Ψ, 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.infinite_alg, env.right; leftorthflag=false) # env gets updated in place + 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.finite_alg, env) + Ψ, env = rtl_sweep!(Ψ, H.middle, t + dt / 2, dt / 2, alg.middle, env) return Ψ, env end \ No newline at end of file From 7c832061c0aff1fff0863d32666b1a076c665a55 Mon Sep 17 00:00:00 2001 From: Daan Maertens <44669533+DaanMaertens@users.noreply.github.com> Date: Tue, 27 Feb 2024 13:43:38 +0100 Subject: [PATCH 38/62] Formatter Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/algorithms/expval.jl | 2 +- src/algorithms/groundstate/find_groundstate.jl | 14 ++++++++------ src/algorithms/propagator/corvector.jl | 8 ++++---- src/algorithms/timestep/tdvp.jl | 6 +++--- 4 files changed, 16 insertions(+), 14 deletions(-) diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl index 040d60ce8..db7ad0db6 100644 --- a/src/algorithms/expval.jl +++ b/src/algorithms/expval.jl @@ -209,7 +209,7 @@ 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) + 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 cdcc367d6..6a012dd17 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 - `verbose::Bool`: display progress information """ -function find_groundstate(ψ::AbstractMPS, H, envs::Union{Cache,MultipleEnvironments}=environments(ψ, H); +function find_groundstate(ψ::AbstractMPS, H, + envs::Union{Cache,MultipleEnvironments}=environments(ψ, H); tol=Defaults.tol, maxiter=Defaults.maxiter, verbose=Defaults.verbose, trscheme=nothing) if isa(ψ, InfiniteMPS) @@ -39,13 +40,13 @@ function find_groundstate(ψ::AbstractMPS, H, envs::Union{Cache,MultipleEnvironm end if isa(ψ, WindowMPS) alg_infin = VUMPS(; tol_galerkin=tol, verbose=verbose, maxiter=maxiter) - alg = Window(alg_infin,alg,alg_infin) + 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} +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) @@ -54,10 +55,11 @@ function find_groundstate!(state::WindowMPS{A,B,VL,VR}, H::Union{Window,LazySum{ (gs_right, _) = find_groundstate(state.right, H.right, alg.right, envs.right) end #set infinite parts - state = WindowMPS(gs_left,state.middle,gs_right) + state = WindowMPS(gs_left, state.middle, gs_right) return find_groundstate(state, H.middle, alg.middle, envs) end -function find_groundstate(ψ::WindowMPS, H::Union{Window,LazySum{<:Window}}, alg::Window, envs...) +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 8f0523591..2a71ea4b5 100644 --- a/src/algorithms/propagator/corvector.jl +++ b/src/algorithms/propagator/corvector.jl @@ -162,7 +162,7 @@ function squaredenvs(Ψ, H, envs, H2, envs2) return H2, envs2 end -function squaredenvs(Ψ::FiniteMPS, H::MPOHamiltonian,envs=environments(Ψ, H)) +function squaredenvs(Ψ::FiniteMPS, H::MPOHamiltonian, envs=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 @@ -170,9 +170,9 @@ function squaredenvs(Ψ::FiniteMPS, H::MPOHamiltonian,envs=environments(Ψ, H)) return squaredenvs(Ψ, H, envs, nH, environments(Ψ, nH)) end -function squaredenvs(Ψ::WindowMPS, H::Window,envs=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)) +function squaredenvs(Ψ::WindowMPS, H::Window, envs=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 6a8e2e110..3e2f27a32 100644 --- a/src/algorithms/timestep/tdvp.jl +++ b/src/algorithms/timestep/tdvp.jl @@ -59,7 +59,7 @@ function timestep(Ψ::InfiniteMPS, H, t::Number, dt::Number, alg::TDVP, end function ltr_sweep!(Ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP, - envs::Union{Cache,MultipleEnvironments}=environments(Ψ, H)) + envs::Union{Cache,MultipleEnvironments}=environments(Ψ, H)) # sweep left to right for i in 1:(length(Ψ) - 1) @@ -78,7 +78,7 @@ function ltr_sweep!(Ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP, end function rtl_sweep!(Ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP, - envs::Union{Cache,MultipleEnvironments}=environments(Ψ, H)) + envs::Union{Cache,MultipleEnvironments}=environments(Ψ, H)) # sweep right to left for i in length(Ψ):-1:2 @@ -119,7 +119,7 @@ algorithm for time evolution. end function ltr_sweep!(Ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP2, - envs=environments(Ψ, H)) + envs=environments(Ψ, H)) # sweep left to right for i in 1:(length(Ψ) - 1) From 479018acd958a4601cb5d8867497bca9a6cad13f Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Tue, 27 Feb 2024 13:48:13 +0100 Subject: [PATCH 39/62] removed obselete exports --- src/MPSKit.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/MPSKit.jl b/src/MPSKit.jl index 17e6a46a5..3592b2736 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -14,7 +14,6 @@ using Base: @kwdef # bells and whistles for mpses export InfiniteMPS, FiniteMPS, WindowMPS, MPSMultiline export PeriodicArray, Window -export fix_left,fix_right,fix_infinite export MPSTensor export QP, LeftGaugedQP, RightGaugedQP export leftorth, From 6063406ea3ce039f9c743bdf875ecd811433ccce Mon Sep 17 00:00:00 2001 From: Daan Maertens <44669533+DaanMaertens@users.noreply.github.com> Date: Tue, 27 Feb 2024 15:20:42 +0100 Subject: [PATCH 40/62] Formatter Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/operators/lazysum.jl | 4 ++-- src/operators/multipliedoperator.jl | 5 +++-- src/states/window.jl | 10 +++++++--- src/states/windowmps.jl | 11 ++++------- 4 files changed, 16 insertions(+), 14 deletions(-) diff --git a/src/operators/lazysum.jl b/src/operators/lazysum.jl index 13b8952a7..4439988f0 100644 --- a/src/operators/lazysum.jl +++ b/src/operators/lazysum.jl @@ -52,10 +52,10 @@ 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) +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) + return map(x -> getproperty(x, sym), sumops) else return getfield(sumops, sym) end diff --git a/src/operators/multipliedoperator.jl b/src/operators/multipliedoperator.jl index e51b474c3..a51c241eb 100644 --- a/src/operators/multipliedoperator.jl +++ b/src/operators/multipliedoperator.jl @@ -51,6 +51,7 @@ 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)) +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/window.jl b/src/states/window.jl index 853af199e..92954d418 100644 --- a/src/states/window.jl +++ b/src/states/window.jl @@ -10,15 +10,19 @@ struct Window{L,M,R} end function Base.:+(window1::Window, window2::Window) - return Window(window1.left+window2.left,window1.middle+window2.middle,window1.right+window2.right) + 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) -_eval_at(x::Window, args...) = Window(_eval_at(x.left,args...),_eval_at(x.middle,args...),_eval_at(x.right,args...)) -safe_eval(::TimeDependent, x::Window, t::Number) = _eval_at(x,t) +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 4c582663e..3abc3187f 100644 --- a/src/states/windowmps.jl +++ b/src/states/windowmps.jl @@ -41,20 +41,17 @@ struct WindowMPS{A<:GenericMPSTensor,B<:MPSBondTensor,VL,VR} <: AbstractFiniteMP window::Window{InfiniteMPS{A,B},FiniteMPS{A,B},InfiniteMPS{A,B}} function WindowMPS(ψₗ::InfiniteMPS{A,B}, ψₘ::FiniteMPS{A,B}, - ψᵣ::InfiniteMPS{A,B}; fixleft::Bool=false, fixright::Bool=false) where {A,B} + ψᵣ::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")) - VL = fixleft ? WINDOW_FIXED : WINDOW_VARIABLE VR = fixright ? WINDOW_FIXED : WINDOW_VARIABLE return new{A,B,VL,VR}(Window(ψₗ, ψₘ, ψᵣ)) end - end - - #=========================================================================================== Constructors ===========================================================================================# @@ -117,8 +114,8 @@ end #=========================================================================================== Utility ===========================================================================================# -function Base.getproperty(ψ::WindowMPS,sym::Symbol) - if sym === :left || sym === :middle || sym === :right +function Base.getproperty(ψ::WindowMPS, sym::Symbol) + if sym === :left || sym === :middle || sym === :right return getfield(ψ.window, sym) elseif sym === :AL return ALView(ψ) From 3dfb7ce3c1d626f0e612b47f2cd04505396d2247 Mon Sep 17 00:00:00 2001 From: Daan Maertens <44669533+DaanMaertens@users.noreply.github.com> Date: Tue, 27 Feb 2024 15:21:37 +0100 Subject: [PATCH 41/62] Formatter Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/algorithms/timestep/timestep.jl | 16 +++++++--------- src/algorithms/toolbox.jl | 5 +++-- src/environments/FinEnv.jl | 2 +- src/environments/multipleenv.jl | 6 +++--- 4 files changed, 14 insertions(+), 15 deletions(-) diff --git a/src/algorithms/timestep/timestep.jl b/src/algorithms/timestep/timestep.jl index da134fe76..773eeee0a 100644 --- a/src/algorithms/timestep/timestep.jl +++ b/src/algorithms/timestep/timestep.jl @@ -1,24 +1,22 @@ # TDVP function timestep!(Ψ::FiniteMPS, 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) + 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(Ψ::FiniteMPS, H, time::Number, timestep::Number, - alg::Union{TDVP,TDVP2}, envs=environments(Ψ, H); kwargs...) - + 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...) - +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 diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl index f0ea188ad..77c4033f7 100644 --- a/src/algorithms/toolbox.jl +++ b/src/algorithms/toolbox.jl @@ -171,10 +171,11 @@ function variance(ψ, H::LazySum, envs=environments(ψ, sum(H))) return variance(ψ, sum(H), envs) end -function variance(ψ::WindowMPS, H::Window, envs::WindowEnv=environments(ψ,H)) +function variance(ψ::WindowMPS, H::Window, envs::WindowEnv=environments(ψ, H)) #tricky to define H2, nenvs = squaredenvs(ψ, H, envs) - return real(expectation_value(ψ, H2.middle, nenvs) - expectation_value(ψ, H.middle, envs)^2) + return real(expectation_value(ψ, H2.middle, nenvs) - + expectation_value(ψ, H.middle, envs)^2) end """ diff --git a/src/environments/FinEnv.jl b/src/environments/FinEnv.jl index b08737cb2..cefedee9f 100644 --- a/src/environments/FinEnv.jl +++ b/src/environments/FinEnv.jl @@ -98,7 +98,7 @@ end # this forces the transfers to be recalculated lazily function invalidate!(ca::FinEnv, ind) ca.ldependencies[ind] = similar(ca.ldependencies[ind]) - ca.rdependencies[ind] = similar(ca.rdependencies[ind]) + return ca.rdependencies[ind] = similar(ca.rdependencies[ind]) end #rightenv[ind] will be contracteable with the tensor on site [ind] diff --git a/src/environments/multipleenv.jl b/src/environments/multipleenv.jl index 48370b817..434790441 100644 --- a/src/environments/multipleenv.jl +++ b/src/environments/multipleenv.jl @@ -30,14 +30,14 @@ end # return MultipleEnvironments(H, broadcast(o -> environments(state, o), H.opps)) # end; - #=========================================================================================== Utility ===========================================================================================# -function Base.getproperty(ca::MultipleEnvironments{<:LazySum,<:WindowEnv},sym::Symbol) +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)) + return MultipleEnvironments(getproperty(ca.opp, sym), + map(x -> getproperty(x, sym), ca)) else return getfield(ca, sym) end From a131db4bf212b670fe492f289599a9a84324b1ca Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Tue, 27 Feb 2024 15:55:14 +0100 Subject: [PATCH 42/62] Formatter (I don't know what I'm doing) --- src/algorithms/timestep/tdvp.jl | 2 +- src/algorithms/timestep/windowtdvp.jl | 14 +++++----- src/environments/windowenv.jl | 39 ++++++++++++++------------- 3 files changed, 30 insertions(+), 25 deletions(-) diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl index 3e2f27a32..dd93f7818 100644 --- a/src/algorithms/timestep/tdvp.jl +++ b/src/algorithms/timestep/tdvp.jl @@ -141,7 +141,7 @@ function ltr_sweep!(Ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP2, end function rtl_sweep!(Ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP2, - envs=environments(Ψ, H)) + envs=environments(Ψ, H)) # sweep right to left for i in length(Ψ):-1:2 diff --git a/src/algorithms/timestep/windowtdvp.jl b/src/algorithms/timestep/windowtdvp.jl index 64c49e959..8652e8f2c 100644 --- a/src/algorithms/timestep/windowtdvp.jl +++ b/src/algorithms/timestep/windowtdvp.jl @@ -16,23 +16,25 @@ 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} - +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) + Ψ = 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) + 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/environments/windowenv.jl b/src/environments/windowenv.jl index 5cb0f424c..f65962417 100644 --- a/src/environments/windowenv.jl +++ b/src/environments/windowenv.jl @@ -12,13 +12,17 @@ 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)) +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)) + return WindowEnv(Window(lenvs, fin_env, renvs), copy(ψ.left.AL), copy(ψ.right.AR)) end -function environments(ψ::WindowMPS, O::MPOHamiltonian, above, lenvs::MPOHamInfEnv,renvs::MPOHamInfEnv) - return environments(ψ, O, above, leftenv(lenvs, 1, ψ.left),rightenv(renvs, length(ψ), ψ.right)) +function environments(ψ::WindowMPS, O::MPOHamiltonian, above, lenvs::MPOHamInfEnv, + renvs::MPOHamInfEnv) + return environments(ψ, O, above, leftenv(lenvs, 1, ψ.left), + rightenv(renvs, length(ψ), ψ.right)) end function environments(below::WindowMPS, above::WindowMPS) @@ -29,12 +33,11 @@ function environments(below::WindowMPS, above::WindowMPS) 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 +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) @@ -44,14 +47,14 @@ function Base.getproperty(ca::WindowEnv,sym::Symbol) #under review 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,ψ) +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) +invalidate!(ca::WindowEnv, ind) = invalidate!(ca.middle, ind) # Check the infinite envs and recalculate the right start function check_rightinfenv!(ca::WindowEnv, ψ::WindowMPS) @@ -67,7 +70,7 @@ 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.middle.leftenvs[1] = leftenv(ca.left, length(ψ.left) + 1, ψ.left) ca.linfdeps .= ψ.left.AL end end @@ -75,20 +78,20 @@ 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) + return leftenv(ca.left, ind, ψ.left) elseif ind > length(ψ) - return leftenv(ca.right,ind,ψ.right) + return leftenv(ca.right, ind, ψ.right) else - return leftenv(finenv(ca,ψ),ind,ψ) + return leftenv(finenv(ca, ψ), ind, ψ) end end function rightenv(ca::WindowEnv, ind, ψ::WindowMPS) if ind < 1 - return rightenv(ca.left,ind,ψ.left) + return rightenv(ca.left, ind, ψ.left) elseif ind > length(ψ) - return rightenv(ca.right,ind,ψ.right) + return rightenv(ca.right, ind, ψ.right) else - return rightenv(finenv(ca,ψ),ind,ψ) + return rightenv(finenv(ca, ψ), ind, ψ) end end From 6e1979ba70b6a9fcd26391f6fd28b1c7689c560d Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Wed, 28 Feb 2024 16:03:31 +0100 Subject: [PATCH 43/62] fix copy for WindowMPS --- src/states/windowmps.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/states/windowmps.jl b/src/states/windowmps.jl index 3abc3187f..043c9db78 100644 --- a/src/states/windowmps.jl +++ b/src/states/windowmps.jl @@ -132,8 +132,10 @@ end 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 - return WindowMPS(left, copy(ψ.middle), right) + fixright = VR === WINDOW_VARIABLE ? false : true + return WindowMPS(left, copy(ψ.middle), right; fixleft=fixleft, fixright=fixright) end # not sure about the underlying methods... From 493772860f35cd39bd62a1dd65c601daf35a8762 Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Wed, 28 Feb 2024 16:04:04 +0100 Subject: [PATCH 44/62] backwards compability --- src/environments/windowenv.jl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/environments/windowenv.jl b/src/environments/windowenv.jl index f65962417..6d55f8686 100644 --- a/src/environments/windowenv.jl +++ b/src/environments/windowenv.jl @@ -25,6 +25,12 @@ function environments(ψ::WindowMPS, O::MPOHamiltonian, above, lenvs::MPOHamInfE rightenv(renvs, length(ψ), ψ.right)) end +#Backwards compability +function environments(ψ::WindowMPS{A,B,WINDOW_FIXED,WINDOW_FIXED}, H::MPOHamiltonian, + above=nothing) where {A,B} + 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")) From da044f158621267097b03cd94b9f307a6edead30 Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Wed, 28 Feb 2024 16:04:24 +0100 Subject: [PATCH 45/62] fix multiplenv of lazysum --- src/environments/multipleenv.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/environments/multipleenv.jl b/src/environments/multipleenv.jl index 434790441..3c75ee0e4 100644 --- a/src/environments/multipleenv.jl +++ b/src/environments/multipleenv.jl @@ -51,6 +51,10 @@ function Base.getproperty(envs::MultipleEnvironments, prop::Symbol) 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 """ Recalculate in-place each sub-env in MultipleEnvironments @@ -61,7 +65,3 @@ function recalculate!(env::MultipleEnvironments, args...; kwargs...) end return env end - - - - From ae5e07070c7b054f33c434f08b646f1d883ab671 Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Wed, 28 Feb 2024 16:05:41 +0100 Subject: [PATCH 46/62] old part of expval for windowmps split off --- src/algorithms/expval.jl | 27 +++++++++++++++++++-------- 1 file changed, 19 insertions(+), 8 deletions(-) diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl index db7ad0db6..45add5a76 100644 --- a/src/algorithms/expval.jl +++ b/src/algorithms/expval.jl @@ -130,16 +130,27 @@ function expectation_value(st::InfiniteMPS, H::MPOHamiltonian, range::UnitRange{ return tot end -# maybe rename this so there is no confusion on what it does? -function expectation_value(ψ::WindowMPS, H::MPOHamiltonian, envs::WindowEnv) - tot = 0.0 + 0im - for i in 1:(H.odim), j in 1:(H.odim) - tot += @plansor leftenv(envs, length(ψ), ψ)[i][1 2; 3] * ψ.AC[end][3 4; 5] * - rightenv(envs, length(ψ), ψ)[j][5 6; 7] * - H[length(ψ)][i, j][2 8; 4 6] * conj(ψ.AC[end][1 8; 7]) +function expectation_value(ψ::WindowMPS, H::MPOHamiltonian, range::UnitRange{Int}, + envs::FinEnv=environments(ψ, H)) + if range == 1:length(ψ) + tot = 0.0 + 0im + for i in 1:(H.odim), j in 1:(H.odim) + tot += @plansor leftenv(envs, length(ψ), ψ)[i][1 2; 3] * ψ.AC[end][3 4; 5] * + rightenv(envs, length(ψ), ψ)[j][5 6; 7] * + H[length(ψ)][i, j][2 8; 4 6] * conj(ψ.AC[end][1 8; 7]) + end + + return tot / (norm(ψ.AC[end])^2) + else + # probably doesn't make sense for any other range + # does it even make sense for unequal left and right states? + throw(ArgumentError("range $range not supported")) end +end - return tot / (norm(ψ.AC[end])^2) +function expectation_value(ψ::WindowMPS, H::Window, range::UnitRange{Int}, + envs::WindowEnv=environments(ψ, H)) + return expectation_value(ψ, H.middle, range, finenv(envs, ψ)) end # Transfer matrices From 8b5e5972c567f74847174689e3fb4859eff73238 Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Wed, 28 Feb 2024 16:06:06 +0100 Subject: [PATCH 47/62] typing in corvector.jl --- src/algorithms/propagator/corvector.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/algorithms/propagator/corvector.jl b/src/algorithms/propagator/corvector.jl index 2a71ea4b5..da7843e7c 100644 --- a/src/algorithms/propagator/corvector.jl +++ b/src/algorithms/propagator/corvector.jl @@ -162,7 +162,8 @@ function squaredenvs(Ψ, H, envs, H2, envs2) return H2, envs2 end -function squaredenvs(Ψ::FiniteMPS, H::MPOHamiltonian, envs=environments(Ψ, H)) +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 @@ -170,7 +171,7 @@ function squaredenvs(Ψ::FiniteMPS, H::MPOHamiltonian, envs=environments(Ψ, H)) return squaredenvs(Ψ, H, envs, nH, environments(Ψ, nH)) end -function squaredenvs(Ψ::WindowMPS, H::Window, envs=environments(Ψ, H)) +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 From 76a77558ca4e9eca06537ec2476caa153a57b6de Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Wed, 28 Feb 2024 16:06:22 +0100 Subject: [PATCH 48/62] Update timestep.jl --- src/algorithms/timestep/timestep.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/algorithms/timestep/timestep.jl b/src/algorithms/timestep/timestep.jl index 773eeee0a..2dac6929b 100644 --- a/src/algorithms/timestep/timestep.jl +++ b/src/algorithms/timestep/timestep.jl @@ -1,6 +1,6 @@ # TDVP -function timestep!(Ψ::FiniteMPS, H, t::Number, dt::Number, alg::Union{TDVP,TDVP2}, +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) @@ -9,7 +9,7 @@ function timestep!(Ψ::FiniteMPS, H, t::Number, dt::Number, alg::Union{TDVP,TDVP end #copying version -function timestep(Ψ::FiniteMPS, H, time::Number, timestep::Number, +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 From ba48826ca3d82109649e90e53706a62e39a2d70a Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Wed, 28 Feb 2024 16:06:35 +0100 Subject: [PATCH 49/62] too much replace --- src/algorithms/toolbox.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl index 77c4033f7..3c561cf1b 100644 --- a/src/algorithms/toolbox.jl +++ b/src/algorithms/toolbox.jl @@ -153,11 +153,11 @@ function variance(state::InfiniteQP, H::MPOHamiltonian, envs=environments(state, state.trivial || throw(ArgumentError("variance of domain wall excitations is not implemented")) - rescaled_H = H - expectation_value(state.left, H) + rescaled_H = H - expectation_value(state.left_gs, H) #I don't remember where the formula came from E_ex = dot(state, effective_excitation_hamiltonian(H, state, envs)) - E_f = expectation_value(state.left, rescaled_H, 1:0) + E_f = expectation_value(state.left_gs, rescaled_H, 1:0) H2 = rescaled_H * rescaled_H @@ -171,11 +171,11 @@ function variance(ψ, H::LazySum, envs=environments(ψ, sum(H))) return variance(ψ, sum(H), envs) end -function variance(ψ::WindowMPS, H::Window, envs::WindowEnv=environments(ψ, H)) +function variance(ψ::WindowMPS, H::Union{MPOHamiltonian,Window}, envs=environments(ψ, H)) #tricky to define H2, nenvs = squaredenvs(ψ, H, envs) - return real(expectation_value(ψ, H2.middle, nenvs) - - expectation_value(ψ, H.middle, envs)^2) + return real(expectation_value(ψ, H2, 1:length(ψ), nenvs) - + expectation_value(ψ, H, 1:length(ψ), envs)^2) end """ From b0ff4e9ca361342ac93232208f5f1f2a9829406d Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Wed, 28 Feb 2024 16:06:50 +0100 Subject: [PATCH 50/62] Update find_groundstate.jl --- src/algorithms/groundstate/find_groundstate.jl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/algorithms/groundstate/find_groundstate.jl b/src/algorithms/groundstate/find_groundstate.jl index 6a012dd17..b5a5664aa 100644 --- a/src/algorithms/groundstate/find_groundstate.jl +++ b/src/algorithms/groundstate/find_groundstate.jl @@ -54,9 +54,12 @@ function find_groundstate!(state::WindowMPS{A,B,VL,VR}, H::Union{Window,LazySum{ if VR === WINDOW_VARIABLE (gs_right, _) = find_groundstate(state.right, H.right, alg.right, envs.right) end - #set infinite parts + #what do we if bond dimension has changed? + #set infinite parts and put through to optimize finite part state = WindowMPS(gs_left, state.middle, gs_right) - return find_groundstate(state, H.middle, alg.middle, envs) + 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, From 7514ae30f1a58014f33dab4ddec19c56403a41df Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Wed, 28 Feb 2024 16:06:55 +0100 Subject: [PATCH 51/62] tests --- test/algorithms.jl | 115 ++++++++++++++++++++++++++++++++++++++++++++- test/states.jl | 84 ++++++++++++++++++++++++++++----- 2 files changed, 185 insertions(+), 14 deletions(-) diff --git a/test/algorithms.jl b/test/algorithms.jl index 0515dd78d..1cc236a3e 100644 --- a/test/algorithms.jl +++ b/test/algorithms.jl @@ -92,6 +92,63 @@ using TensorKit: ℙ @test expectation_value(ψ, Hlazy, envs) ≈ expectation_value(ψ_nolazy, sum(Hlazy), envs_nolazy) atol = 1 - 06 end + + infinite_algs = [VUMPS(; tol_galerkin=tol, verbose=verbosity > 0), + IDMRG1(; tol_galerkin=tol, verbose=verbosity > 0), + GradientGrassmann(; tol=tol, verbosity=verbosity)] + finite_algs = [DMRG(; verbose=verbosity > 0), + DMRG2(; verbose=verbosity > 0, trscheme=truncdim(D)), + DMRG(; verbose=verbosity > 0)] + 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 @@ -116,6 +173,25 @@ end @test (3 + 1.55 - 0.1) * sum(E₀) ≈ sum(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 + Ht = MultipliedOperator(H, t -> 4) + MultipliedOperator(H, 1.45) @testset "Finite TimeDependent LazySum $(alg isa TDVP ? "TDVP" : "TDVP2")" for alg in @@ -156,6 +232,38 @@ end Et = expectation_value(ψt, Ht(1.0), envst) @test sum(E) ≈ sum(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 @@ -186,11 +294,13 @@ end Hwindow = Window(H, repeat(H, 10), H) E₀ = expectation_value(Ψwindow₀, Hwindow) + #= @testset "WindowMPS TDVP" begin ψwindow, envs = timestep(Ψwindow₀, H, dt, TDVP()) E = expectation_value(ψ, H, envs) @test sum(E₀) ≈ sum(E) atol = 1e-2 end + =# end @testset "leading_boundary" verbose = true begin @@ -348,13 +458,14 @@ 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(; verbose=false)) - 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(TensorMap(ComplexF64[0.5 0; 0 -0.5], ℂ^2 ← ℂ^2)) @test expectation_value(gs, szd)[1] ≈ expectation_value(window, szd)[1] atol = 1e-10 polepos = expectation_value(gs, ham, 10) - @test polepos ≈ expectation_value(window, ham)[2] + @test polepos ≈ expectation_value(window, ham, 1:length(window)) vals = (-0.5:0.2:0.5) .+ polepos eta = 0.3im diff --git a/test/states.jl b/test/states.jl index c551c4f97..92e04279a 100644 --- a/test/states.jl +++ b/test/states.jl @@ -96,7 +96,62 @@ end end end -@testset "WindowMPS" begin +@testset "Fixed WindowMPS" begin + ham = force_planar(transverse_field_ising(; g=8.0)) + (gs, _, _) = find_groundstate(InfiniteMPS([ℙ^2], [ℙ^10]), ham, VUMPS(; verbose=false)) + + #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) + + #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(; verbose=false)) + 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(; verbose=false)) @@ -116,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 @@ -127,24 +182,29 @@ end @test 9 * 9 ≈ norm(window)^2 normalize!(window) - e1 = expectation_value(window, ham) + ham = Window(ham, ham, ham) + edens1 = expectation_value(window, ham) + e1 = expectation_value(window, ham, 1:length(window)) v1 = variance(window, ham) - (window, envs, _) = find_groundstate(window, ham, DMRG(; verbose=false)) + gs_alg = Window(VUMPS(; verbose=false), DMRG(; verbose=false), VUMPS(; verbose=false)) + (window, envs, _) = find_groundstate(window, ham, gs_alg) v2 = variance(window, ham) - e2 = expectation_value(window, ham) + edens2 = expectation_value(window, ham) + e2 = expectation_value(window, ham, 1:length(window)) @test v2 < v1 - @test real(e2[2]) ≤ real(e1[2]) + @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, ham) + edens3 = expectation_value(window, ham) + e3 = expectation_value(window, ham, 1:length(window), envs) - @test e2[1] ≈ e3[1] atol = 1e-4 - @test e2[2] ≈ e3[2] atol = 1e-4 + @test edens2 ≈ edens3 atol = 1e-4 + @test e2 ≈ e3 atol = 1e-4 end @testset "Quasiparticle state" verbose = true begin From c07a308ead1be37e8b62d969778ae03b8fd77d1a Mon Sep 17 00:00:00 2001 From: Daan Maertens <44669533+DaanMaertens@users.noreply.github.com> Date: Wed, 28 Feb 2024 16:17:28 +0100 Subject: [PATCH 52/62] Rename FinEnv.jl to finenv.jl --- src/environments/{FinEnv.jl => finenv.jl} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/environments/{FinEnv.jl => finenv.jl} (100%) diff --git a/src/environments/FinEnv.jl b/src/environments/finenv.jl similarity index 100% rename from src/environments/FinEnv.jl rename to src/environments/finenv.jl From 4f356b8ba7333dc038eb1e97b176fb008520f4d6 Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Wed, 28 Feb 2024 22:28:40 +0100 Subject: [PATCH 53/62] WindowMPS docstring update --- src/states/windowmps.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/states/windowmps.jl b/src/states/windowmps.jl index 043c9db78..84921d8ff 100644 --- a/src/states/windowmps.jl +++ b/src/states/windowmps.jl @@ -20,7 +20,7 @@ Type that represents a finite Matrix Product State embedded in an infinite Matri WindowMPS(left::InfiniteMPS, middle_tensors::AbstractVector, right::InfiniteMPS) WindowMPS([f, eltype], physicalspaces::Vector{<:Union{S, CompositeSpace{S}}, virtualspaces::Vector{<:Union{S, CompositeSpace{S}}, left::InfiniteMPS, - [right_gs::InfiniteMPS]) + right_gs::InfiniteMPS) WindowMPS([f, eltype], physicalspaces::Vector{<:Union{S,CompositeSpace{S}}}, maxvirtualspace::S, left::InfiniteMPS, right_gs::InfiniteMPS) @@ -29,10 +29,10 @@ a middle state or a vector of tensors to construct the middle. Alternatively, it to supply the same arguments as for the constructor of [`FiniteMPS`](@ref), followed by a 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`. @@ -60,6 +60,8 @@ function WindowMPS(ψₗ::InfiniteMPS, site_tensors::AbstractVector{<:GenericMPS 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; kwargs...) where {S<:ElementarySpace} From 3c66740fcf6895552d757dc2fe01244eacd2c262 Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Wed, 28 Feb 2024 22:29:19 +0100 Subject: [PATCH 54/62] WindowMPS env extra check --- src/environments/windowenv.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/environments/windowenv.jl b/src/environments/windowenv.jl index 6d55f8686..77c7a5ac7 100644 --- a/src/environments/windowenv.jl +++ b/src/environments/windowenv.jl @@ -28,6 +28,8 @@ 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 From bab69251d9d334b62d27fbe388b16399e4c949b0 Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Thu, 29 Feb 2024 11:39:22 +0100 Subject: [PATCH 55/62] Formatter --- src/environments/windowenv.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/environments/windowenv.jl b/src/environments/windowenv.jl index 77c7a5ac7..38bdeb930 100644 --- a/src/environments/windowenv.jl +++ b/src/environments/windowenv.jl @@ -28,8 +28,10 @@ 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")) + 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 From 2e23fc54b8601f974c2a37145a8fe568fb25ea73 Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Tue, 12 Mar 2024 12:04:03 +0100 Subject: [PATCH 56/62] updated docstring for windowtdvp --- src/algorithms/timestep/windowtdvp.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/algorithms/timestep/windowtdvp.jl b/src/algorithms/timestep/windowtdvp.jl index 8652e8f2c..5e1b38d33 100644 --- a/src/algorithms/timestep/windowtdvp.jl +++ b/src/algorithms/timestep/windowtdvp.jl @@ -4,8 +4,9 @@ [Mixed TDVP](https://arxiv.org/abs/2007.15035) algorithm for time evolution. # Fields -- `finite_alg::A`: algorithm to do the timestep of the finite part of the WindowMPS. This can be `TDVP2` to expand the bonddimension. -- `infinite_alg::A`: algorithm to do the timestep of the infinite part of the WindowMPS +- `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. """ From ca4bf09d818e8875edf4c93710976860ddaa4ac9 Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Tue, 12 Mar 2024 12:04:11 +0100 Subject: [PATCH 57/62] updated windowmps examples --- examples/windowmps.jl | 180 ++++++++++++++++++++++++++++++++++++------ 1 file changed, 154 insertions(+), 26 deletions(-) diff --git a/examples/windowmps.jl b/examples/windowmps.jl index 807c8c807..ad3ff6135 100644 --- a/examples/windowmps.jl +++ b/examples/windowmps.jl @@ -1,37 +1,165 @@ +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")) From 04b10e993e43f5419360c12d81fafefeb9dbe5a0 Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Tue, 12 Mar 2024 15:18:49 +0100 Subject: [PATCH 58/62] small fix in find_groundstate.jl --- src/algorithms/groundstate/find_groundstate.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/algorithms/groundstate/find_groundstate.jl b/src/algorithms/groundstate/find_groundstate.jl index b5a5664aa..748c6901f 100644 --- a/src/algorithms/groundstate/find_groundstate.jl +++ b/src/algorithms/groundstate/find_groundstate.jl @@ -50,13 +50,13 @@ function find_groundstate!(state::WindowMPS{A,B,VL,VR}, H::Union{Window,LazySum{ # 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 - #what do we if bond dimension has changed? - #set infinite parts and put through to optimize finite part - state = WindowMPS(gs_left, state.middle, gs_right) + # then find finite groundstate state, _, delta = find_groundstate(state, H.middle, alg.middle, finenv(envs, state)) return state, envs, delta From 4cca009a07d7731668df37485c7ebb6f5aea6f05 Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Thu, 4 Apr 2024 12:59:05 +0200 Subject: [PATCH 59/62] fixed wrong merge for tests --- .../groundstate/find_groundstate.jl | 2 +- test/algorithms.jl | 71 ++++++++++--------- test/states.jl | 8 +-- 3 files changed, 42 insertions(+), 39 deletions(-) diff --git a/src/algorithms/groundstate/find_groundstate.jl b/src/algorithms/groundstate/find_groundstate.jl index e5eb59e34..aafb6ee1d 100644 --- a/src/algorithms/groundstate/find_groundstate.jl +++ b/src/algorithms/groundstate/find_groundstate.jl @@ -39,7 +39,7 @@ function find_groundstate(ψ::AbstractMPS, H, throw(ArgumentError("Unknown input state type")) end if isa(ψ, WindowMPS) - alg_infin = VUMPS(; tol_galerkin=tol, verbose=verbose, maxiter=maxiter) + alg_infin = VUMPS(; tol=tol, verbosity=verbosity, maxiter=maxiter) alg = Window(alg_infin, alg, alg_infin) end return find_groundstate(ψ, H, alg, envs) diff --git a/test/algorithms.jl b/test/algorithms.jl index 1cc236a3e..d734e5791 100644 --- a/test/algorithms.jl +++ b/test/algorithms.jl @@ -14,13 +14,12 @@ using TensorKit: ℙ @testset "find_groundstate" verbose = true begin tol = 1e-8 verbosity = 0 - infinite_algs = [VUMPS(; tol_galerkin=tol, verbose=verbosity > 0), - IDMRG1(; tol_galerkin=tol, verbose=verbosity > 0), - IDMRG2(; trscheme=truncdim(12), tol_galerkin=tol, - verbose=verbosity > 0), - GradientGrassmann(; tol=tol, verbosity=verbosity), - VUMPS(; tol_galerkin=100 * tol, verbose=verbosity > 0) & - GradientGrassmann(; tol=tol, verbosity=verbosity)] + infinite_algs = [VUMPS(; tol, verbosity), + IDMRG1(; tol, verbosity), + IDMRG2(; trscheme=truncdim(12), tol, verbosity), + GradientGrassmann(; tol, verbosity), + VUMPS(; tol=100 * tol, verbosity) & + GradientGrassmann(; tol, verbosity)] g = 4.0 D = 6 @@ -59,8 +58,8 @@ using TensorKit: ℙ expectation_value(ψ_nolazy, sum(Hlazy), envs_nolazy) atol = 1 - 06 end - finite_algs = [DMRG(; verbose=verbosity > 0), - DMRG2(; verbose=verbosity > 0, trscheme=truncdim(D)), + finite_algs = [DMRG(; verbosity), + DMRG2(; verbosity, trscheme=truncdim(D)), GradientGrassmann(; tol, verbosity, maxiter=300)] H = force_planar(transverse_field_ising(; g)) @@ -93,12 +92,12 @@ using TensorKit: ℙ expectation_value(ψ_nolazy, sum(Hlazy), envs_nolazy) atol = 1 - 06 end - infinite_algs = [VUMPS(; tol_galerkin=tol, verbose=verbosity > 0), - IDMRG1(; tol_galerkin=tol, verbose=verbosity > 0), + infinite_algs = [VUMPS(; tol, verbosity), + IDMRG1(; tol, verbosity), GradientGrassmann(; tol=tol, verbosity=verbosity)] - finite_algs = [DMRG(; verbose=verbosity > 0), - DMRG2(; verbose=verbosity > 0, trscheme=truncdim(D)), - DMRG(; verbose=verbosity > 0)] + 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 @@ -304,7 +303,10 @@ end end @testset "leading_boundary" verbose = true begin - algs = [VUMPS(; tol_galerkin=1e-5, verbose=false), GradientGrassmann(; verbosity=0)] + tol = 1e-5 + verbosity = 0 + algs = [VUMPS(; tol, verbosity), VOMPS(; tol, verbosity), + GradientGrassmann(; verbosity)] mpo = force_planar(classical_ising()) ψ₀ = InfiniteMPS([ℙ^2], [ℙ^10]) @@ -322,7 +324,7 @@ end @testset "infinite (ham)" begin H = repeat(force_planar(heisenberg_XXX()), 2) ψ = InfiniteMPS([ℙ^3, ℙ^3], [ℙ^48, ℙ^48]) - ψ, envs, _ = find_groundstate(ψ, H; maxiter=400, verbose=false, tol=1e-11) + ψ, envs, _ = find_groundstate(ψ, H; maxiter=400, verbosity=0, tol=1e-11) energies, ϕs = excitations(H, QuasiparticleAnsatz(), Float64(pi), ψ, envs) @test energies[1] ≈ 0.41047925 atol = 1e-4 @test variance(ϕs[1], H) < 1e-8 @@ -330,23 +332,23 @@ end @testset "infinite (mpo)" begin H = repeat(sixvertex(), 2) ψ = InfiniteMPS([ℂ^2, ℂ^2], [ℂ^10, ℂ^10]) - ψ, envs, _ = leading_boundary(ψ, H, VUMPS(; maxiter=400, verbose=false)) + ψ, envs, _ = leading_boundary(ψ, H, VUMPS(; maxiter=400, verbosity=0)) energies, ϕs = excitations(H, QuasiparticleAnsatz(), [0.0, Float64(pi / 2)], ψ, - envs; verbose=false) + envs; verbosity=0) @test abs(energies[1]) > abs(energies[2]) # has a minima at pi/2 end @testset "finite" begin - verbose = false + verbosity = 0 H = force_planar(transverse_field_ising()) ψ = InfiniteMPS([ℙ^2], [ℙ^10]) - ψ, envs, _ = find_groundstate(ψ, H; maxiter=400, verbose, tol=1e-9) + ψ, envs, _ = find_groundstate(ψ, H; maxiter=400, verbosity, tol=1e-9) energies, ϕs = excitations(H, QuasiparticleAnsatz(), 0.0, ψ, envs) inf_en = energies[1] fin_en = map([20, 10]) do len ψ = FiniteMPS(rand, ComplexF64, len, ℙ^2, ℙ^15) - (ψ, envs, _) = find_groundstate(ψ, H; verbose) + (ψ, envs, _) = find_groundstate(ψ, H; verbosity) #find energy with quasiparticle ansatz energies_QP, ϕs = excitations(H, QuasiparticleAnsatz(), ψ, envs) @@ -355,7 +357,7 @@ end #find energy with normal dmrg energies_dm, _ = excitations(H, FiniteExcited(; - gsalg=DMRG(; verbose, + gsalg=DMRG(; verbosity, tol=1e-6)), ψ) @test energies_dm[1] ≈ energies_QP[1] + sum(expectation_value(ψ, H, envs)) atol = 1e-4 @@ -457,7 +459,7 @@ 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(; verbose=false)) + 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; fixleft=true, fixright=true) @@ -473,7 +475,7 @@ end predicted = [1 / (v + eta - polepos) for v in vals] @testset "Flavour $f" for f in (Jeckelmann(), NaiveInvert()) - alg = DynamicalDMRG(; flavour=f, verbose=false, tol=1e-8) + alg = DynamicalDMRG(; flavour=f, verbosity=0, tol=1e-8) data = map(vals) do v result, = propagator(window, v + eta, ham, alg) return result @@ -495,7 +497,7 @@ end for λ in [1.05, 2.0, 4.0] H = hamiltonian(λ) ψ = InfiniteMPS([ℂ^2], [ℂ^16]) - ψ, envs, = find_groundstate(ψ, H, VUMPS(; maxiter=100, verbose=false)) + ψ, envs, = find_groundstate(ψ, H, VUMPS(; maxiter=100, verbosity=0)) numerical_scusceptibility = fidelity_susceptibility(ψ, H, [H_X], envs; maxiter=10) @test numerical_scusceptibility[1, 1] ≈ analytical_susceptibility(λ) atol = 1e-2 @@ -503,7 +505,7 @@ end # test if the finite fid sus approximates the analytical one with increasing system size fin_en = map([20, 15, 10]) do L ψ = FiniteMPS(rand, ComplexF64, L, ℂ^2, ℂ^16) - ψ, envs, = find_groundstate(ψ, H, DMRG(; verbose=false)) + ψ, envs, = find_groundstate(ψ, H, DMRG(; verbosity=0)) numerical_scusceptibility = fidelity_susceptibility(ψ, H, [H_X], envs; maxiter=10) return numerical_scusceptibility[1, 1] / L @@ -516,12 +518,12 @@ end @testset "correlation length / entropy" begin st = InfiniteMPS([ℙ^2], [ℙ^10]) th = force_planar(transverse_field_ising()) - (st, _) = find_groundstate(st, th, VUMPS(; verbose=false)) + (st, _) = find_groundstate(st, th, VUMPS(; verbosity=0)) len_crit = correlation_length(st)[1] entrop_crit = entropy(st) th = force_planar(transverse_field_ising(; g=4)) - (st, _) = find_groundstate(st, th, VUMPS(; verbose=false)) + (st, _) = find_groundstate(st, th, VUMPS(; verbosity=0)) len_gapped = correlation_length(st)[1] entrop_gapped = entropy(st) @@ -532,7 +534,7 @@ end @testset "expectation value / correlator" begin st = InfiniteMPS([ℂ^2], [ℂ^10]) th = transverse_field_ising(; g=4) - (st, _) = find_groundstate(st, th, VUMPS(; verbose=false)) + (st, _) = find_groundstate(st, th, VUMPS(; verbosity=0)) sz_mpo = TensorMap([1.0 0; 0 -1], ℂ^1 * ℂ^2, ℂ^2 * ℂ^1) sz = TensorMap([1.0 0; 0 -1], ℂ^2, ℂ^2) @@ -560,6 +562,7 @@ end end @testset "approximate" verbose = true begin + verbosity = 0 @testset "mpo * infinite ≈ infinite" begin st = InfiniteMPS([ℙ^2, ℙ^2], [ℙ^10, ℙ^10]) th = force_planar(repeat(transverse_field_ising(; g=4), 2)) @@ -570,9 +573,9 @@ end W1 = convert(DenseMPO, sW1) W2 = convert(DenseMPO, sW2) - st1, _ = approximate(st, (sW1, st), VUMPS(; verbose=false)) - st2, _ = approximate(st, (W2, st), VUMPS(; verbose=false)) - st3, _ = approximate(st, (W1, st), IDMRG1(; verbose=false)) + st1, _ = approximate(st, (sW1, st), VUMPS(; verbosity)) + st2, _ = approximate(st, (W2, st), VUMPS(; verbosity)) + st3, _ = approximate(st, (W1, st), IDMRG1(; verbosity)) st4, _ = approximate(st, (sW2, st), IDMRG2(; trscheme=truncdim(20), verbose=false)) st5, _ = timestep(st, th, 0.0, dt, TDVP()) st6 = changebonds(W1 * st, SvdCut(; trscheme=truncdim(10))) @@ -586,7 +589,7 @@ end @test dim(space(nW1.opp[1, 1], 1)) == 1 end - finite_algs = [DMRG(; verbose=false), DMRG2(; verbose=false, trscheme=truncdim(10))] + finite_algs = [DMRG(; verbosity), DMRG2(; verbosity, trscheme=truncdim(10))] @testset "finitemps1 ≈ finitemps2" for alg in finite_algs a = FiniteMPS(10, ℂ^2, ℂ^10) b = FiniteMPS(10, ℂ^2, ℂ^20) @@ -634,7 +637,7 @@ end ψ = FiniteMPS(len, ℂ^2, ℂ^10) - gs, envs = find_groundstate(ψ, th, DMRG(; verbose=false)) + gs, envs = find_groundstate(ψ, th, DMRG(; verbosity=0)) #translation mpo: @tensor bulk[-1 -2; -3 -4] := isomorphism(ℂ^2, ℂ^2)[-2, -4] * diff --git a/test/states.jl b/test/states.jl index 92e04279a..5955f24a1 100644 --- a/test/states.jl +++ b/test/states.jl @@ -98,7 +98,7 @@ end @testset "Fixed WindowMPS" begin ham = force_planar(transverse_field_ising(; g=8.0)) - (gs, _, _) = find_groundstate(InfiniteMPS([ℙ^2], [ℙ^10]), ham, VUMPS(; verbose=false)) + (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; fixleft=true, @@ -132,7 +132,7 @@ end e1 = expectation_value(window, ham, 1:length(window)) v1 = variance(window, ham) - (window, envs, _) = find_groundstate(window, ham, DMRG(; verbose=false)) + (window, envs, _) = find_groundstate(window, ham, DMRG(; verbosity=0)) v2 = variance(window, ham) edens2 = expectation_value(window, ham) @@ -153,7 +153,7 @@ end @testset "Variable WindowMPS" begin ham = force_planar(transverse_field_ising(; g=8.0)) - (gs, _, _) = find_groundstate(InfiniteMPS([ℙ^2], [ℙ^10]), ham, VUMPS(; verbose=false)) + (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) @@ -187,7 +187,7 @@ end e1 = expectation_value(window, ham, 1:length(window)) v1 = variance(window, ham) - gs_alg = Window(VUMPS(; verbose=false), DMRG(; verbose=false), VUMPS(; verbose=false)) + gs_alg = Window(VUMPS(; verbosity=0), DMRG(; verbosity=0), VUMPS(; verbosity=0)) (window, envs, _) = find_groundstate(window, ham, gs_alg) v2 = variance(window, ham) From 142e4092985d72aa98e3bbe542841cfc24bc90b5 Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Thu, 4 Apr 2024 13:59:25 +0200 Subject: [PATCH 60/62] forgot to fix everything last commit --- test/algorithms.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/algorithms.jl b/test/algorithms.jl index d734e5791..4fc5d459e 100644 --- a/test/algorithms.jl +++ b/test/algorithms.jl @@ -573,10 +573,10 @@ end W1 = convert(DenseMPO, sW1) W2 = convert(DenseMPO, sW2) - st1, _ = approximate(st, (sW1, st), VUMPS(; verbosity)) - st2, _ = approximate(st, (W2, st), VUMPS(; verbosity)) + st1, _ = approximate(st, (sW1, st), VOMPS(; verbosity)) + st2, _ = approximate(st, (W2, st), VOMPS(; verbosity)) st3, _ = approximate(st, (W1, st), IDMRG1(; verbosity)) - st4, _ = approximate(st, (sW2, st), IDMRG2(; trscheme=truncdim(20), verbose=false)) + st4, _ = approximate(st, (sW2, st), IDMRG2(; trscheme=truncdim(20), verbosity)) st5, _ = timestep(st, th, 0.0, dt, TDVP()) st6 = changebonds(W1 * st, SvdCut(; trscheme=truncdim(10))) From 2ea9ee14039957d666c5df5f844ad7cdddfe9adc Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Thu, 28 Nov 2024 14:37:56 +0100 Subject: [PATCH 61/62] made environments(::WindowMPS) less strict --- src/environments/windowenv.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/environments/windowenv.jl b/src/environments/windowenv.jl index 38bdeb930..ff34e3e0a 100644 --- a/src/environments/windowenv.jl +++ b/src/environments/windowenv.jl @@ -19,8 +19,9 @@ function environments(ψ::WindowMPS, O::Window, above=nothing; return WindowEnv(Window(lenvs, fin_env, renvs), copy(ψ.left.AL), copy(ψ.right.AR)) end -function environments(ψ::WindowMPS, O::MPOHamiltonian, above, lenvs::MPOHamInfEnv, - renvs::MPOHamInfEnv) +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 From 7fcf8a1553120d0c1e7a1b01fdf72a03fc89e39b Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Thu, 28 Nov 2024 14:42:36 +0100 Subject: [PATCH 62/62] Format of windowmps.jl --- examples/windowmps.jl | 40 ++++++++++++++++++++++------------------ 1 file changed, 22 insertions(+), 18 deletions(-) diff --git a/examples/windowmps.jl b/examples/windowmps.jl index ad3ff6135..e93b855b8 100644 --- a/examples/windowmps.jl +++ b/examples/windowmps.jl @@ -12,9 +12,9 @@ function my_transverse_field_ising(gs) return a + b end -function my_timedependent_ising(gl,gs,gr,f) +function my_timedependent_ising(gl, gs, gr, f) L = length(gs) - lattice = InfiniteChain(1) + lattice = InfiniteChain(1) latticeL = InfiniteChain(L) ZZ = rmul!(σᶻᶻ(), -1) X = rmul!(σˣ(), -1) @@ -27,16 +27,18 @@ function my_timedependent_ising(gl,gs,gr,f) Xm = @mpoham sum(gs[i] * X{i} for i in vertices(latticeL)) Xr = @mpoham sum(gr * X{i} for i in vertices(lattice)) - H1 = Window(ZZl,ZZm,ZZr) - H2 = Window(Xl,Xm,Xr) - return LazySum([H1,MultipliedOperator(H2,f)]) + H1 = Window(ZZl, ZZm, ZZr) + H2 = Window(Xl, Xm, Xr) + return LazySum([H1, MultipliedOperator(H2, f)]) end -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) +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) + right = expectation_value(Ψwindow.right, O.right) + return vcat(left, middle, right) end function my_finalize(t, Ψ, H, envs, si, tosave) @@ -75,7 +77,7 @@ 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)); + 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); @@ -143,23 +145,25 @@ Ts = fill(TensorMap(rand, ComplexF64, ℂ^D * ℂ^2, ℂ^D), L); (Ψwindow, envs_gs, _) = find_groundstate(Ψwindow, Hgs); #define the quench Hamiltonian -f(t) = 0.1*t #we take a linear ramp +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 = 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); +envs = environments(Ψwindow, Hquench); -sdat = [my_expectation_value(Ψwindow, Window(sx,sx,sx))] +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)); +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")) +display(heatmap(t_span, 0:(L + 1), real.(reduce((a, b) -> [a b], sdat)); xlabel="t", + ylabel="i"))