new standards

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")

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

src/algorithms/timestep/integrators.jl

@@ -1,29 +1,40 @@
 """
     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
 - `y₀`: object to integrate
 - `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

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,

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...)

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;

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

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);

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)))