Time-dependent TDVP1

0.4.0](https://github.com/ITensor/ITensorMPS.jl/releases/tag/v0.4.0]), the TimeDependentSum functionality was removed from the ITensorMPS package, which means that the following tutorial will not run with version ITensorMPS version 0.4.0 or higher.

In this tutorial we will see how to solve the Schrödinger equation with a time-dependent Hamiltonian using the TDVP1 method.

Our Hamiltonian is composed of two parts \(H\sb0\) and \(H\sb1\) such that

\[\begin{equation} H(t) = H_0 + f(t) H_1, \qquad f(t) = \min\{t/T, 1\}, \label{eq:time-dependent-hamiltonian} \end{equation}\]

where \(H\sb0\) is the Hamiltonian of a collection of free electrons

\[H_0 = \sum_{k=1}^{N} \sum_{s\in\{\spinup,\spindown\}} \varepsilon \adj{c_{k,s}} c_{k,s}\phantomadj + \sum_{k=1}^{N-1} \sum_{s\in\{\spinup,\spindown\}} \lambda (\adj{c_{k,s}} c_{k+1,s}\phantomadj + \adj{c_{k+1,s}} c_{k,s}\phantomadj)\]

and \(H\sb1\) represents on-site Coulomb repulsion,

\[H_1 = \sum_{k=1}^{N} U n_{k,\spinup} n_{k,\spindown}.\]

Let's define these operators with ITensor. We also choose as the initial condition a state where all sites contain an electron in the “down” state.

julia> using ITensorMPS, MPSTimeEvolution

julia> N = 6; s = siteinds("Electron", N);

julia> ψ₀ = MPS(s, "↓");

julia> ε = 0.1; λ = 0.2; u = -2ε;

julia> h₀ = OpSum(); h₁ = OpSum();

julia> for k in 1:N
           h₀ += ε, "Ntot", k
       end

julia> for k in 1:N-1
           h₀ += λ, "Cdagup", k, "Cup", k+1
           h₀ += λ, "Cdagdn", k, "Cdn", k+1
           h₀ += λ, "Cup", k, "Cdagup", k+1
           h₀ += λ, "Cdn", k, "Cdagdn", k+1
       end

julia> H₀ = MPO(h₀, s);

julia> for k in 1:N
           h₁ += u, "Nup * Ndn", k
       end

julia> H₁ = MPO(h₁, s);

In order to implement the time-dependent ODE problem, we will need to use a custom solver argument as the first argument of tdvp1!. This argument is normally implicitly set to the default solver, the exponentiate method from KrylovKit. First, we define a time-dependent version of exponentiate (with the same structure as exponentiate):

julia> using KrylovKit: exponentiate

julia> using ITensorMPS: TimeDependentSum

julia> function time_dependent_exp(
           H::TimeDependentSum, time_step, ψ₀; current_time=0.0, outputlevel=0, kwargs...
       )
           ψₜ, info = exponentiate(H(current_time), time_step, ψ₀; kwargs...)
           return ψₜ, info
       end;

The first argument TimeDependentSum is a struct, defined in ITensorMPS (but not exported), that contains two members called coefficients and terms. When a TimeDependentSum object S is called on a number t, it yields something like sum(c(t) * s for (c, s) in zip(coefficients, terms)): this means that it can be used to define a function like \(H(t)\). In our case, we can write \(H(t)\) from \eqref{eq:time-dependent-hamiltonian} as follows, once we choose a value for \(T\):

julia> T = 2;

julia> ramp(t) = min(t / T, one(t));

julia> fs = [one, ramp]; Hs = [H₀, H₁];

Here one is the function that returns the multiplicative identity for its argument, i.e. one(t) returns 1.0 if t is a Float64. We define the time-dependent solver as follows:

julia> function time_dependent_solver(PHs::ProjMPOSum, time_step, ψ₀; kwargs...)
           return time_dependent_exp(
               TimeDependentSum(fs, PHs), time_step, ψ₀; ishermitian=true, kwargs...
           )
       end;

This function will be used inside tdvp1!, precisely in the tdvp_site_update! method that computes the update of a single site during each sweep. We also add the ishermitian=true flag, that will be forwarded to the inner exponentiate method.

Now we set the variables related to time, and the callback operator:

julia> dt = 0.01; tmax = 2T;

julia> cb = ExpValueCallback("Nup(1,5),Ndn(1,5)", s, dt)
ExpValueCallback
Operators: Nup(1), Nup(5), Ndn(1) and Ndn(5)
No measurements performed

One last thing: during the execution of time_dependent_exp, inside the stack of function calls, the iterate method is called on the PHs object. Unfortunately, at the time this tutorial is written, no such function is defined by ITensorMPS, so we need to do a bit of “plumbing” ourselves. Specifically, we need to extend the two Base.iterate methods to the ProjMPOSum type, so that Julia can use them:

julia> Base.iterate(PH::ProjMPOSum, state) = iterate(PH.terms, state)

julia> Base.iterate(PH::ProjMPOSum) = iterate(PH.terms)

Now that iterate is properly extended, we can finally call tdvp1! with everything we defined:

julia> tdvp1!(
           time_dependent_solver,
           ψ₀,
           Hs,
           dt,
           tmax,
           callback=cb,
           ishermitian=true,
           progress=false
       );