From 2883dc19badaa78d5e0b66fb3cddcbc8ae18f1ff Mon Sep 17 00:00:00 2001 From: DaanMaertens Date: Tue, 10 Oct 2023 15:55:54 +0200 Subject: [PATCH] new standards --- src/MPSKit.jl | 2 + src/algorithms/expval.jl | 43 ++++-------- src/algorithms/timestep/integrators.jl | 31 +++++--- src/algorithms/timestep/tdvp.jl | 17 ++--- src/algorithms/timestep/time_evolve.jl | 97 ++++++++++++++++++++++++++ src/environments/multipleenv.jl | 5 ++ src/operators/multipliedoperator.jl | 14 ---- src/operators/sumofoperators.jl | 10 +-- src/utility/linearcombination.jl | 7 ++ 9 files changed, 152 insertions(+), 74 deletions(-) create mode 100644 src/algorithms/timestep/time_evolve.jl create mode 100644 src/utility/linearcombination.jl diff --git a/src/MPSKit.jl b/src/MPSKit.jl index 0a9d2d0f..607f80b6 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -63,6 +63,7 @@ module Defaults const tol = 1e-12 const verbose = true _finalize(iter, state, opp, envs) = (state, envs) + _time_finalize(iter, state, opp, envs) = (state, envs, tobesaved) import KrylovKit: GMRES, Arnoldi const linearsolver = GMRES(; tol, maxiter) @@ -72,6 +73,7 @@ end include("utility/periodicarray.jl") include("utility/utility.jl") #random utility functions include("utility/plotting.jl") +include("utility/linearcombination.jl") #maybe we should introduce an abstract state type include("states/window.jl") diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl index 742b0ad3..16b8d1e6 100644 --- a/src/algorithms/expval.jl +++ b/src/algorithms/expval.jl @@ -192,58 +192,39 @@ end expectation_value(Ψ, op, t, args...) = expectation_value(Ψ, op, args...) # define expectation_value for MultipliedOperator as scalar multiplication of the non-multiplied result, instead of multiplying the operator itself -function expectation_value(Ψ, op::TimedOperator, t::Number, args...) - return op.f(t) * expectation_value(Ψ, op.op, args...) -end -function expectation_value(Ψ, op::UntimedOperator, t::Number, args...) - return expectation_value(Ψ, op::UntimedOperator, args...) -end + function expectation_value(Ψ, op::UntimedOperator, args...) return op.f * expectation_value(Ψ, op.op, args...) end # define expectation_value for SumOfOperators -function expectation_value(Ψ, ops::SumOfOperators, t::Number, at::Int64) - return sum(op -> expectation_value(Ψ, op, t, at), ops) +function expectation_value(Ψ, ops::SumOfOperators, at::Int64) + return sum(op -> expectation_value(Ψ, op, at), ops) end function expectation_value( - Ψ, ops::SumOfOperators, t::Number, envs::MultipleEnvironments=environments(Ψ, ops) + Ψ, ops::SumOfOperators, envs::MultipleEnvironments=environments(Ψ, ops) ) - return sum(map((op, env) -> expectation_value(Ψ, op, t, env), ops.ops, envs)) + return sum(map((op, env) -> expectation_value(Ψ, op, env), ops.ops, envs)) end # define expectation_value for Window - -# when no time is given -function expectation_value(Ψ::WindowMPS, windowH::Window, at::Int64) - return expectation_value(Ψ, windowH, 0.0, at) -end - -function expectation_value( - Ψ::WindowMPS, windowH::Window, windowEnvs=environments(Ψ, windowH) -) - return expectation_value(Ψ, windowH, 0.0, windowEnvs) -end - -# with time argument -function expectation_value(Ψ::WindowMPS, windowOp::Window, t::Number, at::Int64) +function expectation_value(Ψ::WindowMPS, windowOp::Window, at::Int64) if at < 1 - return expectation_value(Ψ.left_gs, windowOp.left, t, at) + return expectation_value(Ψ.left_gs, windowOp.left, at) elseif 1 <= at <= length(Ψ.window) - return expectation_value(Ψ, windowOp.middle, t, at) + return expectation_value(Ψ, windowOp.middle, at) else - return expectation_value(Ψ.right_gs, windowOp.right, t, at) + return expectation_value(Ψ.right_gs, windowOp.right, at) end end function expectation_value( Ψ::WindowMPS, windowH::Window, - t::Number, windowEnvs::Window{C,D,C}=environments(Ψ, windowH), ) where {C<:Union{MultipleEnvironments,Cache},D<:Union{MultipleEnvironments,Cache}} - left = expectation_value(Ψ.left_gs, windowH.left, t, windowEnvs.left) - middle = expectation_value(Ψ.window, windowH.middle, t, windowEnvs.middle) - right = expectation_value(Ψ.right_gs, windowH.right, t, windowEnvs.right) + 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/integrators.jl b/src/algorithms/timestep/integrators.jl index 32c0edcc..af3f8df0 100644 --- a/src/algorithms/timestep/integrators.jl +++ b/src/algorithms/timestep/integrators.jl @@ -1,8 +1,7 @@ """ integrate(f, y₀, t₀, a, dt, algorithm) -Integrate the differential equation ``dy/dt = a f(y,t)`` over a time step 'dt' starting from ``y(t₀)=y₀``, using the provided algorithm. -For time-independent operators (i.e. not a TimedOperator) t₀ is ingored. +Integrate the differential equation ``dy/dt = a*f(y,t)`` over a time step 'dt' starting from ``y(t₀)=y₀``, using the provided algorithm. # Arguments - `f::Function`: driving function @@ -10,20 +9,32 @@ For time-independent operators (i.e. not a TimedOperator) t₀ is ingored. - `t₀`: time f is evaluated at - `a` : scalar prefactor - `dt`: timestep -- `method`: method to integrate TDVP equations +- `algorithm`: integration scheme """ function integrate end -# for backwards compatibility -integrate(f, y₀, a, dt, method) = integrate(f, y₀, 0.0, a, dt, method) +# make time evaluation dispatchable +# user provides f(y,t) that can be called with two arguments +Eval(f,x,t::Number) = f(x) +Eval(f::F,x,t::Number) where {O<:TimedOperator,F<:Union{O,SumOfOperators{O}}} = f(x,t) -# wrap function into UntimedOperator by default -integrate(f, y₀, t₀, a, dt, method) = integrate(UntimedOperator(f), y₀, t₀, a, dt, method) + +""" + ExpIntegrator + +Method that solves ``dy/dt = a*f(y,t)`` by exponentiating the action of f(y,t). + +# Fields +- `krylovmethod::Union{Arnoldi,Lanczos}`` KrylovKit method for exponentiation. For options such as tolerance, see Lanczos/Arnoldi in KrylovKit. +""" +@kwdef struct ExpIntegrator + krylovmethod::Union{Arnoldi,Lanczos} = Lanczos(); +end #original integrator in iTDVP, namely exponentiation function integrate( - f::F, y₀, t₀, a, dt, method::Union{Arnoldi,Lanczos} -) where {F<:Union{MultipliedOperator,SumOfOperators}} - sol, convhist = exponentiate(x -> f(x, t₀ + dt / 2), a * dt, y₀, method) + f, y₀, t₀::Number, a::Number, dt::Number, method::ExpIntegrator +) + sol, convhist = exponentiate( x->Eval_t(f,x,t₀), a*dt, y₀, method.krylovmethod) return sol, convhist.converged, convhist end diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl index 69124c27..f6231293 100644 --- a/src/algorithms/timestep/tdvp.jl +++ b/src/algorithms/timestep/tdvp.jl @@ -22,12 +22,12 @@ algorithm for time evolution. # Fields - `integrator::A`: integration algorithm (defaults to Lanczos exponentiation) - `tolgauge::Float64`: tolerance for gauging algorithm -- `maxiter::Int`: maximum amount of gauging iterations +- `gaugemaxiter::Int`: maximum amount of gauging iterations """ @kwdef struct TDVP{A} <: Algorithm integrator::A = Lanczos(; tol=Defaults.tol) tolgauge::Float64 = Defaults.tolgauge - maxiter::Int = Defaults.maxiter + gaugemaxiter::Int = Defaults.maxiter end function timestep( @@ -67,7 +67,7 @@ function timestep( Qc, _ = leftorth!(temp_CRs[loc]; alg=TensorKit.QRpos()) @plansor temp_ACs[loc][-1 -2; -3] = QAc[-1 -2; 1] * conj(Qc[-3; 1]) end - newΨ = InfiniteMPS(temp_ACs, Ψ.CR[end]; tol=alg.tolgauge, maxiter=alg.maxiter) + newΨ = InfiniteMPS(temp_ACs, Ψ.CR[end]; tol=alg.tolgauge, maxiter=alg.gaugemaxiter) else for loc in 1:length(Ψ) @@ -76,7 +76,7 @@ function timestep( _, Qc = rightorth!(temp_CRs[mod1(loc - 1, end)]; alg=TensorKit.LQpos()) temp_ACs[loc] = _transpose_front(Qc' * QAc) end - newΨ = InfiniteMPS(Ψ.CR[0], temp_ACs; tol=alg.tolgauge, maxiter=alg.maxiter) + newΨ = InfiniteMPS(Ψ.CR[0], temp_ACs; tol=alg.tolgauge, maxiter=alg.gaugemaxiter) end recalculate!(envs, newΨ) @@ -148,13 +148,13 @@ algorithm for time evolution. # Fields - `integrator::A`: integrator algorithm (defaults to Lanczos exponentiation) - `tolgauge::Float64`: tolerance for gauging algorithm -- `maxiter::Int`: maximum amount of gauging iterations +- `gaugemaxiter::Int`: maximum amount of gauging iterations - `trscheme`: truncation algorithm for [tsvd][TensorKit.tsvd](@ref) """ @kwdef struct TDVP2{A} <: Algorithm integrator::A = Lanczos(; tol=Defaults.tol) tolgauge::Float64 = Defaults.tolgauge - maxiter::Int = Defaults.maxiter + gaugemaxiter::Int = Defaults.maxiter trscheme = truncerr(1e-3) end @@ -206,11 +206,6 @@ function timestep!( return Ψ, envs end -# time-independent version -function timestep(Ψ, H, dt, alg, env=environments(Ψ, H); kwargs...) - return timestep(Ψ, H, 0.0, dt, alg, env; kwargs...) -end - #copying version function timestep( Ψ::AbstractFiniteMPS, diff --git a/src/algorithms/timestep/time_evolve.jl b/src/algorithms/timestep/time_evolve.jl new file mode 100644 index 00000000..3fe65a1b --- /dev/null +++ b/src/algorithms/timestep/time_evolve.jl @@ -0,0 +1,97 @@ +#put in defaults +_time_finalize(iter, state, opp, envs) = (state, envs, nothing) + +""" + SimpleScheme{F} <: Algorithm + +The simplest numerical scheme to do time evolution by doing consecutive individual timesteps. + +# Fields +- `stepalg::S`: algorithm used to do the individual timesteps +- `ts`::AbstractVector: time step points +- `dts`::AbstractVector: dt of each time step +- `time_finalize::F`: user-supplied function which is applied after each iteration which can be used to save objects during the time evolution, with + signature `finalize(iter, Ψ, H, envs) -> Ψ, envs, saved` + +""" +struct SimpleScheme{S,F,N} <: MPSKit.Algorithm + stepalg::S + ts::AbstractVector{N} + dts::AbstractVector{N} + time_finalize::F + + function SimpleScheme(stepalg::S,ts::AbstractVector{N},dts::AbstractVector{N},time_finalize::F) where {N<:Number,S,F} + length(ts) == length(dts)+1 || throw(ArgumentError("times and timesteps length need to be compatible")) + all(isapprox.(test.ts[1:end-1] .+ test.dts,test.ts[2:end])) || throw(ArgumentError("times and timesteps need to be compatible")) + return new{S,F,N}(stepalg,ts,dts,time_finalize) + end +end +function SimpleScheme(stepalg,ts::AbstractVector,time_finalize) + return SimpleScheme(stepalg,ts,ts[2:end] .- ts[1:end-1],time_finalize) +end + +function SimpleScheme(stepalg,ts::AbstractRange,time_finalize) + return SimpleScheme(stepalg,ts,repeat([step(ts)],length(ts)-1),time_finalize) +end + +#implement iteration +function Base.iterate(x::SimpleScheme,state=1) + if length(x.ts) == 0 ||state == length(x.ts) + return nothing + else + return ((x.ts[state],x.dts[state]),state+1) + end +end + +""" + time_evolve(Ψ, H, tspan, [environments]; kwargs...) + time_evolve(Ψ, H, scheme, [environments]; kwargs...) + +Time evolve the initial state `Ψ` with Hamiltonian `H` over a time span `tspan`. If not specified the +time step `scheme` defaults to `SimpleScheme` which is just a consecutive iteration over `tspan`. +If not specified, an algorithm for the individual time steps will be attempted based on the supplied keywords. + +## Arguments +- `Ψ::AbstractMPS`: initial state +- `H::AbstractMPO`: operator that generates the time evolution (can be time-dependent). +- `[environments]`: MPS environment manager +- `tspan::AbstractVector`: time points over which the time evolution is stepped +- `scheme`: time step scheme + +## Keywords +- `tol::Float64`: tolerance for convergence criterium +- `maxiter::Int`: maximum amount of iterations +- `verbose::Bool`: display progress information +""" +function time_evolve!(Ψ,H,tspan::AbstractVector{<:Number},envs::Cache=environments(Ψ,H); + integrator= Lanczos(; tol=tol), + tol=Defaults.tol, + tolgauge=Defaults.tolgauge, + gaugemaxiter=Defaults.maxiter, + verbose=Defaults.verbose, + trscheme=nothing, + stepalg = TDVP(;integrator = integrator, tol = tol, tolgauge=tolgauge, gaugemaxiter = gaugemaxiter, verbose=verbose), + time_finalize=Defaults._time_finalize, + saved = [] + ) + if !isnothing(trscheme) + stepalg = TDVP2(;integrator = integrator, tol = tol, tolgauge=tolgauge, gaugemaxiter = gaugemaxiter, verbose=verbose, trscheme=trscheme) + end + scheme = SimpleScheme(stepalg,tspan,time_finalize) + return time_evolve!(Ψ,H,scheme,envs;saved=saved) +end + +time_evolve(Ψ,H,tspan,envs::Cache=environments(Ψ,H); kwargs...) = time_evolve(copy(Ψ),H,tspan,envs; kwargs...) + +function time_evolve!(Ψ,H,scheme::SimpleTimeScheme,envs=environments(Ψ,H);saved=[]) + _,_, tobesaved = alg.finalize(1, Ψ, H, envs) + isnothing(tobesaved) || push!(saved,tobesaved) + for (iter,(t,dt)) in enumerate(scheme) + Ψ,envs = timestep!(Ψ,H,t,dt,scheme.stepalg,envs) + + Ψ, envs, tobesaved = scheme.finalize(iter, Ψ, H, envs)::Tuple{typeof(Ψ),typeof(envs),eltype(saved)} + isnothing(tobesaved) || push!(saved,tobesaved) + end + return Ψ, envs, saved +end +time_evolve(Ψ,H,scheme,envs::Cache=environments(Ψ,H); kwargs...) = time_evolve!(copy(Ψ),H,scheme,envs; kwargs...) diff --git a/src/environments/multipleenv.jl b/src/environments/multipleenv.jl index 59c5c17e..c7832c6c 100644 --- a/src/environments/multipleenv.jl +++ b/src/environments/multipleenv.jl @@ -15,6 +15,11 @@ function environments(st, ham::SumOfOperators) return MultipleEnvironments(ham, map(op -> environments(st, op), ham.ops)) end +#broadcast vs map? +function environments(state, ham::LinearCombination) + return MultipleEnvironments(ham, broadcast(o -> environments(state, o), ham.opps)) +end; + function environments( st::WindowMPS, ham::SumOfOperators; diff --git a/src/operators/multipliedoperator.jl b/src/operators/multipliedoperator.jl index 42ad81ac..18f006a7 100644 --- a/src/operators/multipliedoperator.jl +++ b/src/operators/multipliedoperator.jl @@ -29,15 +29,6 @@ UntimedOperator(x::O, c::C) where {C<:Real,O} = MultipliedOperator(x, c) TimedOperator(x) = TimedOperator(x, t -> 1) UntimedOperator(x) = UntimedOperator(x, 1) -TimedOperator(x::UntimedOperator) = TimedOperator(x.op, t -> x.f) -function UntimedOperator(x::TimedOperator) - throw( - ArgumentError( - "Cannot make a UntimedOperator from a TimedOperator without a time argument" - ), - ) -end - # what to do when we multiply by a scalar function Base.:*(op::UntimedOperator, b::Number) return UntimedOperator(op.op, b * op.f) @@ -47,11 +38,6 @@ function Base.:*(op::TimedOperator, b::Number) end Base.:*(b::Number, op::MultipliedOperator) = op * b -(x::UntimedOperator)() = x.f * x.op - -#ignore time-dependence by default -(x::MultipliedOperator)(t::Number) = x - (x::TimedOperator)(t::Number) = UntimedOperator(x.op, x.f(t)) # logic for derivatives diff --git a/src/operators/sumofoperators.jl b/src/operators/sumofoperators.jl index 89fd3a48..28802d8e 100644 --- a/src/operators/sumofoperators.jl +++ b/src/operators/sumofoperators.jl @@ -19,14 +19,13 @@ SumOfOperators(x) = SumOfOperators([x]) function SumOfOperators(ops::AbstractVector, fs::AbstractVector) return SumOfOperators(map((op, f) -> MultipliedOperator(op, f), ops, fs)) end -LinearCombination(ops::Tuple, fs::Tuple) = SumOfOperators(collect(ops), collect(fs)) # we can add operators to SumOfOperators by using + function Base.:+(op1::MultipliedOperator{O,F}, op2::MultipliedOperator{O,G}) where {O,F,G} return SumOfOperators([op1, op2]) end -# this we can also do with promote +# this we could also do with promote function Base.:+(op1::TimedOperator{O}, op2::Union{O,UntimedOperator{O}}) where {O} return SumOfOperators(TimedOperator{O}[op1, TimedOperator(op2)]) end @@ -58,12 +57,7 @@ function Base.:+(op::O, SumOfOps::SumOfOperators{T}) where {O,T<:MultipliedOpera return SumOfOperators(UntimedOperator(op)) + SumOfOps end -(x::SumOfOperators{<:UntimedOperator})() = sum(op -> op(), x) - -#ignore time-dependence by default -(x::SumOfOperators)(t::Number) = x - -(x::SumOfOperators{<:MultipliedOperator})(t::Number) = SumOfOperators(map(op -> op(t), x)) #will convert to SumOfOperators{UnTimedOperator} +(x::SumOfOperators{<:TimedOperator})(t::Number) = SumOfOperators(map(op -> op(t), x)) #will convert to SumOfOperators{UnTimedOperator} # logic for derivatives Base.:*(x::SumOfOperators, v) = x(v); diff --git a/src/utility/linearcombination.jl b/src/utility/linearcombination.jl new file mode 100644 index 00000000..9fb3d375 --- /dev/null +++ b/src/utility/linearcombination.jl @@ -0,0 +1,7 @@ +struct LinearCombination{O<:Tuple,C<:Tuple} + opps::O + coeffs::C +end + +Base.:*(h::LinearCombination, v) = sum((c * (o * v) for (o, c) in zip(h.opps, h.coeffs))) +(h::LinearCombination)(v) = sum((c * o(v) for (o, c) in zip(h.opps, h.coeffs)))