QuantumPropagators Package

Initialize a Propagator.

propagator = initprop(
    state, generator, tlist;
    method=:auto,
    backward=false,
    inplace=true,
    piecewise=nothing,
    pwc=nothing,
    kwargs...
)

initializes a propagator for the time propagation of the given state over a time grid tlist under the time-dependent generator (Hamiltonian/Liouvillian) generator.

Arguments

• state: The "initial" state for the propagation. For backward=false, this state is taken to be at initial time (tlist[begin]); and for backward=true, at the final time (tlist[end])
• generator: The time-dependent generator of the dynamics
• tlist: The time grid over which which the propagation is defined. This may or may not be equidistant.

Keyword arguments

• method: The propagation method to use. The default value of :auto attempts to choose the best method available, based on the properties of the given state, tlist, and generator, cf. choose_propmethod
• backward: If true, initialize the propagator for a backward propagation. The resulting propagator.t will be tlist[end], and subsequent calls to propstep! will move backward on tlist.
• inplace: If true, the state property of the resulting propagator will be changed in-place by any call to propstep!. If false, each call to propstep! changes the reference for propgator.state, and the progation will not use any in-place operations. Not all propagation methods may support both in-place and not-in-place propagation. In-place propagation is generally more efficient but may not be compatible, e.g., with automatic differentiation.
• piecewise: If given a a boolean, true enforces that the resulting propagator is a PiecewisePropagator, and false enforces is not to be a PiecewisePropagator
• pwc: Like piecewise, for for the stronger PWCPropagator

All other kwargs are method-dependent and are ignored for methods that do not support them.

The type of the returned propagator is a sub-type of AbstractPropagator, respectively a sub-type of PiecewisePropagator if piecewise=true or a sub-type of PWCPropagator if pwc=true.

See also

• reinitprop! — Re-initialize a propagator
• propagate — Higher-level propagation interface

cheby_propagator = initprop(
    state,
    generator,
    tlist;
    method=:cheby,
    inplace=true,
    backward=false,
    verbose=false,
    parameters=nothing,
    control_ranges=nothing,
    specrange_method=:auto,
    specrange_buffer=0.01,
    cheby_coeffs_limit=1e-12,
    check_normalization=false,
    specrange_kwargs...
)

initializes a ChebyPropagator.

Method-specific keyword arguments

• control_ranges: a dict the maps the controls in generator (see getcontrols) to a tuple of min/max values. The Chebychev coefficients will be calculated based on a spectral envelope that assumes that each control can take arbitrary values within the min/max range. If not given, the ranges are determined automatically. Specifying manual control ranges can be useful when the the control amplitudes (parameters) may change during the propagation, e.g. in a sequential-update control scheme.
• specrange_method: Method to pass to the specrange function
• specrange_buffer: An additional factor by which to enlarge the estimated spectral range returned by specrange, in order to ensure that Chebychev coefficients are based on an overestimation of the spectral range.
• cheby_coeffs_limit: The maximum magnitude of Chebychev coefficients that should be treated as non-zero
• check_normalization: Check whether the Hamiltonian has been properly normalized, i.e., that the spectral range of generator has not been underestimated. This slowes down the propagation, but is advisable for novel generators.
• specrange_kwargs: All further keyword arguments are passed to the specrange function

newton_propagator = initprop(
    state,
    generator,
    tlist,
    method::Val{:newton};
    inplace=true,
    backward=false,
    verbose=false,
    parameters=nothing,
    m_max=10,
    func=(z -> exp(-1im * z)),
    norm_min=1e-14,
    relerr=1e-12,
    max_restarts=50,
    _...
)

initializes a NewtonPropagator.

Method-specific keyword arguments

• m_max: maximum Krylov dimension, cf. NewtonWrk
• func, norm_min, relerr, max_restarts: parameter to pass to newton!

exp_propagator = initprop(
    state,
    generator,
    tlist,
    method::Val{:expprop};
    inplace=true,
    backward=false,
    verbose=false,
    parameters=nothing,
    func=(H_dt -> exp(-1im * H_dt))
    convert_state=_exp_prop_convert_state(state),
    convert_operator=_exp_prop_convert_operator(generator),
    _...
)

initializes an ExpPropagator.

Method-specific keyword arguments

• func: The function to evaluate. The argument H_dt is obtained by constructing an operator H from generator via the evalcontrols function and the multiplied with the time step dt for the current time interval. The propagation then simply multiplies the return value of func with the current state
• convert_state: Type to which to temporarily convert states before multiplying the return value of func.
• convert_operator: Type to which to convert the operator H before multiplying it with dt and plugging the result into func

The convert_state and convert_operator parameters are useful for when the generator and or state are unusual data structures for which the relevant methods to calculate func are not defined. Often, it is easier to temporarily convert them to standard complex matrices and vectors than to implement the missing methods.

Propagate a state over an entire time grid.

state = propagate(
    state,
    generator,
    tlist;
    method=:auto,
    backward=false,
    inplace=true,
    verbose=false,
    piecewise=nothing,
    pwc=nothing,
    storage=nothing,
    observables=(<store state>, ),
    callback=nothing,
    showprogress=false,
    initprop_kwargs...)

propagates state of the entire time grid and returns the propagates states, or a storage array of data collected during the propagation.

Arguments

• state: The "initial" state for the propagation. For backward=false, this state is taken to be at initial time (tlist[begin]); and for backward=true, at the final time (tlist[end])
• generator: The time-dependent generator of the dynamics
• tlist: The time grid over which which the propagation is defined. This may or may not be equidistant.

Keyword arguments

• method: The propagation method to use. The default value of :auto attempts to choose the best method available, based on the properties of the given state, tlist, and generator.
• backward: If true, propagate backward in time
• inplace: If true, propagate using in-place operations. If false, avoid in-place operations. Not all propagation methods support both in-place and not-in-place propagation.
• piecewise: If given a a boolean, limit the propagation to "piecewise" methods, respectively disallow piecewise methods
• pwc: If given a a boolean, limit the propagation to piecewise-constant methods, respectively disallow piecewise-constant methods
• storage: Flag whether to store and return the propagated states / observables, or pre-allocated storage array. See Notes below.
• observables: Converters for data to be stored in storage. See Notes below.
• callback: Function to call after each propagation step. See Notes below.
• showprogess: Whether to show a progress bar. See Notes below.

All remaining keyword arguments are passed to initprop to initialize the Propagator that is used internally to drive the optimization. Unknown keyword arguments will be ignored.

Notes

In general, there is no requirement that tlist has a constant time step, although some propagation methods (most notably cheby!) only support a uniform time grid.

If storage is given as an Array, it will be filled with data determined by the observables. The default "observable" results in the propagated states at every point in time being stored. The storage array should be created with init_storage. See its documentation for details.

The storage parameter may also be given as true, and a new storage array will be created internally with init_storage and returned instead of the propagated state:

data = propagate(
    state, generator, tlist; method=:auto
    backward=false; storage=true, observables=observables,
    callback=nothing, showprogress=false, kwargs...)

If backward is true, the input state is assumed to be at time tlist[end], and the propagation progresses backward in time (with a negative time step dt). If storage is given, it will be filled back-to-front during the backward propagation.

If callback is given as a callable, it will be called after each propagation step, as callback(propagator, observables) where propagator is Propagator object driving the propagation. The callback is called before calculating any observables. Example usage includes writing data to file, or modifying state via set_state!, e.g., removing amplitude from the lowest and highest level to mitigate "truncation error".

If showprogress is given as true, a progress bar will be shown for long-running propagationn. In order to customize the progress bar, showprogress may also be a function that receives length(tlist) and returns a ProgressMeter.Progress instance.

If in_place=false is given, the propagation avoids in-place operations. This is slower than inplace=true, but is often required in the context of automatic differentiation (AD), e.g., with Zygote. That is, use in_place=false if propagate is called inside a function to be passed to Zygote.gradient, Zygote.pullback, or a similar function. In an AD context, storage and showprogress should not be used.

The propagate routine returns the propagated state at tlist[end], respectively tlist[1] if backward=true, or a storage array with the stored states / observable data if storage=true.

See also

• initprop — Propagate via a Propagator object

Advance the propagator by a single time step.

state = propstep!(propagator)

returns the state obtained from propagating to the next point on the time grid from propagator.t, respectively the previous point if propagator.backward is true.

When the propagation would lead out of the time grid, propstep! leaves propagator unchanged and returns nothing. Thus, a return value of nothing may be used to signal that a propagation has completed.

Re-initialize a propagator.

reinitprop!(propagator, state; kwargs...)

resets the propagator to state at the beginning of the time grid, respectively the end of the time grid if propagator.backward is true.

At a minimum, this is equivalent to a call to set_state! follow by a call to set_t!, but some propagators may have additional requirements on re-initialization, such as refreshing expansion coefficients for ChebyPropagator. In this case, the kwargs may be additional keyword arguments specific to the concrete type of propagator.

reinitprop!(
    propagator::ChebyPropagator,
    state;
    transform_control_ranges=((c, ϵ_min, ϵ_max, check) => (ϵ_min, ϵ_max)),
    kwargs...
)

re-initializes an existing ChebyPropagator. This may or may not involve recalculating the Chebychev coefficients based on the current control amplitudes in propagator.parameters.

Method-specific keyword arguments

• transform_control_ranges: a function (c, ϵ_min, ϵ_max, check) => (ϵ_min′, ϵ_max′). For each control c, the function is called with check=true and ϵ_min (ϵ_max) the current minimum (maximum) values for the control from propagator.parameters). The Chebychev coefficients will be recalculated if the existing coefficients were obtained assuming a range for c outside the returned ϵ_min′, ϵ_max′.

  If the coefficients do need to be recalculated, transform_control_ranges is called a second time with check=false, and the returned (ϵ_min′, ϵ_max′) are used for estimating the new spectral range.

  For example,

  function transform_control_ranges(c, ϵ_min, ϵ_max, check)
      if check
          return (min(ϵ_min, 2 * ϵ_min), max(ϵ_max, 2 * ϵ_max))
      else
          return (min(ϵ_min, 5 * ϵ_min), max(ϵ_max, 5 * ϵ_max))
      end
  end

  will re-calculate the Chebychev coefficients only if the current amplitudes differ by more than a factor of two from the ranges that were used when initializing the propagator (control_ranges parameter in initprop, which would have had to overestimate the actual amplitudes by at least a factor of two). When re-calculating, the control_ranges will overestimate the amplitudes by a factor of five. With this transform_control_ranges, the propagation will be stable as long as the amplitudes do not change dynamically by more than a factor of 2.5 from their original range, while also not re-calculating coefficients unnecessarily in each pass because of modest changes in the amplitudes.

  The transform_control_ranges argument is only relevant in the context of optimal control, where the same propagator will be used for many iterations with changing control field amplitudes.

All other keyword arguments are ignored.

Set the current state of the propagator.

set_state!(propagator, state)

sets the propagator.state property. In order to mutate the current state after a call to propstep!, the following pattern is recommended:

Ψ = propagator.state
mutate!(Ψ)
set_state!(propagator, Ψ)

This is guaranteed to work efficiently both for in-place and not-in-place propagators, without incurring unnecessary copies.

mutate!(propagator.state)

set_state!(propagator, propagator.state)

See also

• set_t! — set propagator.t.

AbstractPropagator sub-type for piecewise propagators.

A piecewise propagator is determined by a single parameter per control and time grid interval. Consequently, the propagator.parameters are a dictionary mapping the controls found in the generator via getcontrols to a vector of values defined on the intervals of the time grid, see discretize_on_midpoints. This does not necessarily imply that these values are the piecewise-constant amplitudes for the intervals. A general piecwise propagatore might use interpolation to obtain actual amplitudes within any given time interval.

When the amplitudes are piecewise-constant, the propagator should be a concrete intantiation of a PWCPropagator.

PiecewisePropagator sub-type for piecewise-constant propagators.

Like the more general PiecewisePropagator, this is characterized by propagator.parameters mapping the controls in the generator to a vector of amplitude value on the midpoints of the time grid intervals.

The propagation will use these values as constant within each interval.

Abstract base type for all Propagator objects.

All Propagator objects must be instantiated via initprop and implement the following interface.

Properties

• state (read-only): The current quantum state in the propagation
• tlist (read-only): The time grid for the propatation
• t (read-only): The time at which state is defined. An element of tlist.
• parameters: parameters that determine the dynamics. The structure of the parameters depends on the concrete Propagator type (i.e., the propagation method). Mutating the parameters affects subsequent propagation steps.
• backward: Boolean flag to indicate whether the propagation moves forward or backward in time
• inplace: Boolean flag to indicate whether propagator.state is modified in-place or is recreated by every call of propstep! or set_state!. For inplace=true, we find Ψ = propagator.state; propstep!(propagator); propagator.state === Ψ to be true, while for inplace=false it is false.

Concrete Propagator types may have additional properties or fields, but these should be considered private.

Methods

• reinitprop! — reset the propagator to a new initial state at the beginning of the time grid (or the end, for backward propagation)
• propstep! – advance the propagator by one step forward (or backward) on the time grid.
• set_state! — safely mutate the current quantum state of the propagation. Note that directly mutating the state property is not safe. However, Ψ = propagator.state; mutate!(Ψ), set_state!(propagator, Ψ) is guaranteed to be safe and efficient for both in-place and not-in-place propagators.
• set_t! — safely mutate the current time (propagator.t), snapping to the values of tlist.

See also

• PiecewisePropagator — specialization of AbstractPropagator for piecewise propagation methods.
• PWCPropagator — specialization of PiecewisePropagator for piecewise-constant propagation methods.

Choose a suitable propagation method.

method = choose_propmethod(generator, state, tlist;
                           pwc=nothing, piecewise=nothing, inplace=true)

identifies a suitable propagation method for the given generator, state and tlist. If piecewise or pwc are given as true, only consider methods that result in in a PiecewisePropagator or PWCPropagator, respectively. If piecewise or pwc are given as false, disregard any methods that result in these propagators. Only propagators that support the given inplace are taken into account.

Determine the spectral envelope of a generator.

E_min, E_max = cheby_get_spectral_envelope(
    generator, tlist, control_ranges, method; kwargs...
)

estimates a lower bound E_min the lowest eigenvalue of the generator for any values of the controls specified by control_ranges, and an upper bound E_max for the highest eigenvalue.

This is done by constructing operators from the extremal values for the controls as specified in control_ranges and taking the smallest/largest return values from specrange for those operators.

Arguments

• generator: dynamical generator, e.g. a time-dependent
• tlist: The time grid for the propagation
• control_ranges: a dict that maps controls that occur in generator (cf. getcontrols to a tuple of mimimum and maximum amplitude for that control
• method: method name to pass to specrange
• kwargs: Any remaining keyword arguments are passed to specrange

Set the current time for the propagation.

set_t!(propagator, t)

Sets propagator.t to the given value of t, where t must be an element of propagator.tlist.

See also

• set_state! — set propagator.state

Propagator for Chebychev propagation (method=:cheby).

This is a PWCPropagator.

Propagator for propagation via direct exponentiation (method=:expprop)

This is a PWCPropagator.

Propagator for Newton propagation (method=:newton).

This is a PWCPropagator.

m = arnoldi!(Hess, q, m, Ψ, H, dt=1.0; extended=true, norm_min=1e-15)

Calculate the Hessenberg matrix and Arnoldi vectors of H dt, from Ψ.

For a given order m, the m×m Hessemberg matrix is calculated and stored in in the pre-allocated Hess. Further an array of m normalized Arnoldi vectors is stored in in the pre-allocated q, plus one additional unnormalized Arnoldi vector. The unnormalized m+1st vector could be used to easily extend a given m×m Hessenberg matrix to a (m+1)×(m+1) matrix.

If the extended Hessenberg matrix is requested (extended=true, default), the m+1st Arnoldi vector is also normalized, and it's norm will be stored in m+1, m entry of the (extended) Hessenberg matrix, which is an (m+1)×(m+1) matrix.

Return the size m of the calculated Hessenberg matrix. This will usually be the input m, except when the Krylov dimension of H starting from Ψ is less then m. E.g., if Ψ is an eigenstate of H, the returned m will be 1.

See http://en.wikipedia.org/wiki/Arnoldi_iteration for a description of the algorithm.

Arguments

• Hess::Matrix{ComplexF64}: Pre-allocated storage for the Hessemberg matrix. Can be uninitialized on input. The matrix must be at least of size m×m, or (m+1)×(m+1) if extended=true. On output, the m×m sub-matrix of Hess (with the returned output m) will contain the Hessenberg matrix, and all other elements of Hess be be set to zero.
• q: Pre-allocated array of states similar to Ψ, as storage for the calculated Arnoldi vectors. These may be un-initialized on input. Must be at least of length m+1
• m: The requested dimensions of the output Hessenberg matrix.
• Ψ: The starting vector for the Arnoldi procedure. This can be of any type, as long as Φ = H * Ψ results in a vector similar to Ψ, there is an inner products of Φ and Ψ (Ψ⋅Φ is defined), and norm(Ψ) is defined.
• H: The operator (up to dt) for which to calculate the Arnoldi procedure. Can be of any type, as long as H * Ψ is defined.
• dt: The implicit time step; the total operator for which to calculate the Arnoldi procedure is H * dt
• extended: If true (default), calculate the extended Hessenberg matrix, and normalized the final Arnoldi vector
• norm_min: the minimum value of the norm of Ψ at which Ψ should be considered the zero vector

Extend dimension of Hessenberg matrix by one.

extend_arnoldi!(Hess, q, m, H, dt; norm_min=1e-15)

extends the entries in Hess from size (m-1)×(m-1) to size m×m, and the list q of Arnoldi vectors from m to (m+1). It is assumed that the input Hess was created by a call to arnoldi! with extended=false or a previous call to extend_arnoldi!. Note that Hess itself is not resized, so it must be allocated to size m×m or greater on input.

diagonalize_hessenberg_matrix(Hess, m; accumulate=false)

Diagonalize the m × m top left submatrix of the given Hessenberg matrix.

If accumulate is true, return the concatenated eigenvalues for Hess[1:1,1:1] to Hess[1:m,1:m], that is, all sumatrices of size 1 through m.

Workspace for the Chebychev propagation routine.

ChebyWrk(Ψ, Δ, E_min, dt; limit=1e-12)

initializes the workspace for the propagation of a state similar to Ψ under a Hamiltonian with eigenvalues between E_min and E_min + Δ, and a time step dt. Chebychev coefficients smaller than the given limit are discarded.

Evaluate Ψ = exp(i- H dt) Ψ.

Ψ_out = cheby(Ψ, H, dt, wrk; E_min=nothing, check_normalization=false)

acts like cheby! but does not modify Ψ in-place.

Evaluate Ψ = exp(-i H dt) Ψ in-place.

cheby!(Ψ, H, dt, wrk; E_min=nothing, check_normalization=false)

Arguments

• Ψ: on input, initial vector. Will be overwritten with result.
• H: Hermitian operator
• dt: time step
• wrk: internal workspace
• E_min: minimum eigenvalue of H, to be used instead of the E_min from the initialization of wrk. The same wrk may be used for different values E_min, as long as the spectra radius Δ and the time step dt are the same as those used for the initialization of wrk.
• check_normalizataion: perform checks that the H does not exceed the spectral radius for which the the workspace was initialized.

The routine will not allocate any internal storage. This implementation requires copyto! lmul!, and axpy! to be implemented for Ψ, and the three-argument mul! for Ψ and H.

Calculate Chebychev coefficients.

a::Vector{Float64} = cheby_coeffs(Δ, dt; limit=1e-12)

return an array of coefficiencts larger than limit.

Arguments

• Δ: the spectral radius of the underlying operator
• dt: the time step

See also cheby_coeffs! for an in-place version.

Calculate Chebychev coefficients in-place.

n::Int = cheby_coeffs!(coeffs, Δ, dt, limit=1e-12)

overwrites the first n values in coeffs with new coefficients larger than limit for the given new spectral radius Δ and time step dt. The coeffs array will be resized if necessary, and may have a length > n on exit.

See also cheby_coeffs for an non-in-place version.

Evaluate control at every point of tlist.

values = discretize(control, tlist; via_midpoints=true)

discretizes the given control to a Vector of values defined on the points of tlist.

If control is a function, it will will first be evaluated at the midpoint of tlist, see discretize_on_midpoints, and then the values on the midpoints are converted to values on tlist. This discretization is more stable than directly evaluationg the control function at the values of tlist, and ensures that repeated round-trips between discretize and discretize_on_midpoints can be done safely, see the note in the documentation of discretize_on_midpoints.

The latter can still be achieved by passing via_midpoints=false. While such a direct discretization is suitable e.g. for plotting, but it is unsuitable for round-trips between discretize and discretize_on_midpoints (constant controls on tlist may result in a zig-zag on the intervals of tlist).

If control is a vector, it will be returned un-modified if it is of the same length as tlist. Otherwise, control must have one less value than tlist, and is assumed to be defined on the midpoins of tlist. In that case, discretize acts as the inverse of discretize_on_midpoints. See discretize_on_midpoints for how control values on tlist and control values on the intervals of tlist are related.

Evaluate control at the midpoints of tlist.

values = discretize_on_midpoints(control, tlist)

discretizes the given control to a Vector of values on the midpoints of tlist. Hence, the resulting values will contain one less value than tlist.

If control is a vector of values defined on tlist (i.e., of the same length as tlist), it will be converted to a vector of values on the intervals of tlist. The value for the first and last "midpoint" will remain the original values at the beginning and end of tlist, in order to ensure exact bounary conditions. For all other midpoints, the value for that midpoint will be calculated by "un-averaging".

For example, for a control and tlist of length 5, consider the following diagram:

tlist index:       1   2   3   4   5
tlist:             ⋅   ⋅   ⋅   ⋅   ⋅   input values cᵢ (i ∈ 1..5)
                   |̂/ ̄ ̄ ̂\ / ̂\ / ̂ ̄ ̄\|̂
midpoints:         x     x   x     x   output values pᵢ (i ∈ 1..4)
midpoints index:   1     2   3     4

We will have $p₁=c₁$ for the first value, $p₄=c₅$ for the last value. For all other points, the control values $cᵢ = \frac{p_{i-1} + p_{i}}{2}$ are the average of the values on the midpoints. This implies the "un-averaging" for the midpoint values $pᵢ = 2 c_{i} - p_{i-1}$.

If control is a vector of one less length than tlist, it will be returned unchanged, under the assumption that the input is already properly discretized.

If control is a function, the function will be directly evaluated at the midpoints marked as x in the above diagram..

Replace the controls in generator with static values.

G = evalcontrols(generator, vals_dict)

replaces the time-dependent controls in generator with the values in vals_dict and returns the static operator G.

The vals_dict is a dictionary (IdDict) mapping controls as returned by getcontrols(generator) to values.

In-place version of evalcontrols.

evalcontrols!(G, generator, vals_dict)

acts as evalcontrols, but modifies G in-place.

Shift time grid values the interval midpoints

tlist_midpoints = get_tlist_midpoints(tlist)

takes a vector tlist of length $n$ and returns a vector of length $n-1$ containing the midpoint values of each interval. The intervals in tlist are not required to be uniform.

Get the derivative of the generator $G$ w.r.t. the control $ϵ(t)$.

μ  = getcontrolderiv(generator, control)

returns nothing if the generator (Hamiltonian or Liouvillian) does not depend on control, or a function μ(v) that evaluates

\[μ(v) = \left.\frac{∂G}{∂ϵ(t)}\right\vert_{ϵ(t)=v}\]

otherwise. That is, a call μ(v) will return the static operator resulting from evaluating the derivative of the dynamical generator $G$ with respect to the control filed $ϵ(t)$ at a particular point in time where the control field takes the value $v$.

Note that for the common case of linear control terms, e.g., $Ĥ = Ĥ_0 + \sum_l ϵ_l(t) Ĥ_l$, the derivative $∂Ĥ/∂ϵ_l(t)$ is simply the control Hamiltonian $Ĥ_l$. Thus, the resulting function μ will simply return $Ĥ_l$, ignoring the argument v.

Get a vector of the derivatives of generator w.r.t. each control.

getcontrolderivs(generator, controls)

return as vector containing the derivative of generator with respect to each control in controls. The elements of the vector are either nothing if generator does not depend on that particular control, or a function μ(α) that evaluates the derivative for a particular value of the control, see getcontrolderiv.

Extract a Tuple of controls.

controls = getcontrols(generator)

extracts the controls from a single dynamical generator.

By default, assumes that any generator is a nested Tuple, e.g. (H0, (H1, ϵ1), (H2, ϵ2), ...) and extracts (ϵ1, ϵ2)

Each control must be a valid argument for discretize.

getcontrols(operator)

for a static operator (matrix) returns an empty tuple.

controls = getcontrols(objectives)

extracts the controls from a list of objectives (i.e., from each objective's generator). Controls that occur multiple times in the different objectives will occur only once in the result.

Substitute the controls inside a generator with different controls.

new_generator = substitute_controls(generator, controls_map)

Creates a new generator from generator by replacing any control that is in the dict controls_map with controls_map[control]. Controls that are not in controls_map are kept unchanged.

The substituted controls must be time-dependent; to substitute static values for the controls, converting the time-depdentned generator into a static operator, use evalcontrols.

ExpPropWrk(v0)

Workspace for propagation via direct matrix exponentiation.

Initializes the workspace for the propagation of a vector v0

expprop!(Ψ, H, dt, wrk; func=(H_dt -> exp(-1im * H_dt)), _...)

Evaluate Ψ = func(H*dt) Ψ by directly evaluating U = func(H*dt), i.e. by matrix exponentiation for the default func, and then multiplying U and Ψ in-place with mul!.

The workspace wrk must be initialized with ExpPropWrk to provide storage for a temporary state.

Keyword arguments besides func are ignored.

Ψ_out = expprop(Ψ, H, dt, wrk; func=(H_dt -> exp(-1im * H_dt)), _...)

evaluates Ψ_out = func(H*dt) Ψ as in expprop!, but not acting in-place.

NewtonWrk(v0, m_max=10)

Workspace for the Newton-with-restarted-Arnoldi propagation routine.

Initializes the workspace for the propagation of a vector v0, using a maximum Krylov dimension of m_max in each restart iteration. Note that m_max should be smaller than the length of v0.

newton!(Ψ, H, dt, wrk; func=(z -> exp(-1im*z)), norm_min=1e-14, relerr=1e-12,
        max_restarts=50, _...)

Evaluate Ψ = func(H*dt) Ψ using a Newton-with-restarted-Arnoldi scheme.

Arguments

• Ψ: The state to propagate, will be overwritten in-place with the propagated state
• H: Operator acting on Ψ. Together with dt, this is the argument to func
• dt: Implicit time step. Together with H, this is the argument to func
• wkr: Work array, initialized with NewtonWrk
• func: The function to apply to H dt, taking a single (scalar) complex-valued argument z in place of H dt. The default func is to evaluate the time evoluation operator for the Schrödinger equation
• norm_min: the minium norm at which to consider a state similar to Ψ as zero
• relerr: The relative error defining the convergence condition for the restart iteration. Propagation stops when the norm of the accumulated Ψ is stable up to the given relative error
• max_restarts: The maximum number of restart iterations. Exceeding max_restarts will throw an AssertionError.

All other keyword arguments are ignored.

extend_leja!(leja, n, newpoints, n_use)

Given an array of n (ordered) Leja points, extract n_use points from newpoints, and append them to the existing Leja points. The array leja should be sufficiently large to hold the new Leja points, which are appended after index n_old. It will be re-allocated if necessary and may have a size of up to 2*(n+n_use).

Arguments

• leja: Array of leja values. Must contain the "old" leja values to be kept in leja(0:n-1). On output, n_use new leja points will be in leja(n+:n+n_use-1), for the original value of n. The leja array must use zero-based indexing.
• n: On input, number of "old" leja points in leja. On output, total number of leja points (i.e. n=n+n_use)
• newpoints: On input, candidate points for new leja points. The n_use best values will be chosen and added to leja. On output, the values of new_points are undefined.
• n_use: Number of points that should be added to leja

extend_newton_coeffs!(a, n_a, leja, func, n_leja, radius)

Extend the array a of existing Newton coefficients for the expansion of the func from n_a coefficients to n_leja coefficients. Return a new value n_a=n_a+n_leja with the total number of Newton coefficients in the updated a.

Arguments

• a: On input, a zero-based array of length n_a or greater, containing Newton coefficients. On output, array containing a total n_leja coefficients. The array a will be resized if necessary, and may have a length greater than n_leja on output
• n_a: The number of Newton coefficients in a, on input. Elements of a beyond the first n_a elements will be overwritten.
• leja: Array of normalized Leja points, containing at least n_leja elements.
• func: Function for which to calcluate Newton coeffiecients
• n_leja: The number of elements in leja to use for calculating new coefficients, and the total number of Newton coefficients on output
• radius: Normalization radius for divided differences

Calculate the spectral range of a Hamiltonian H on the real axis.

E_min, E_max = specrange(H; method=:auto, kwargs...)

calculates the approximate lowest and highest eigenvalues of H. Any imaginary part in the eigenvalues is ignored: the routine is intended for (although not strictly limited to) a Hermitian H.

This delegates to

specrange(H, method; kwargs...)

for the different methods.

The default method=:auto chooses the best method for the given H. This is :diag for small matrices, and :arnoldi otherwise. Keyword arguments not relevant to the underlying implementation will be ignored.

E_min, E_max = specrange(H, :arnoldi; state=random_state(H), m_min=20,
                         m_max=60, prec=1e-3, norm_min=1e-15, enlarge=true)

uses Arnoldi iteration with state as the starting vector. It approximates the eigenvalues of H with between m_min and m_max Ritz values, until the lowest and highest eigenvalue are stable to a releative precision of prec. The norm_min parameter is passed to the underlying arnoldi!.

If enlarge=true (default) the returned E_min and E_max will be enlarged via a heuristic to slightly over-estimate the spectral radius instead of under-estimating it.

E_min, E_max = specrange(H, :diag)

uses exact diagonization via the standard eigvals function to obtain the smallest and largest eigenvalue. This should only be used for relatively small matrices.

Calculate a vector for Ritz values converged to a given precision.

R = ritzvals(G, state, m_min, m_max=2*m_min; prec=1e-5, norm_min=1e-15)

calculates a complex vector R of at least m_min and at most m_max Ritz values.

Random normalized quantum state.

    Ψ = random_state(H)

returns a random normalized state compatible with the Hamiltonian H. This is intended to provide a starting vector for estimating the spectral radius of H via an Arnoldi method.

Obtain data from storage.

get_from_storage!(state, storage, i)

extracts data from the storage for the i'th time slot. Invese of write_to_storage!

Create a storage array for propagation.

storage = init_storage(state, tlist)

creates a storage array suitable for storing a state for each point in tlist.

storage = init_storage(state, tlist, observables))

creates a storage array suitable for the data generated by the observables applied to state, see map_observables, for each point in tlist.

storage = init_storage(data, nt))

creates a storage arrays suitable for storing data nt times, where nt=length(tlist). By default, this will be a vector of typeof(data) and length nt, or a n × nt Matrix with the same eltype as data if data is a Vector of length n.

Apply a single observable to state.

data = map_observable(observable, state)

By default, observable is assumed to be callable, and the above is equivalent to data = observable(state).

If observable is a matrix and state is a vector evaluate the expectation value of the observable as dot(state, observable, state).

Obtain "observable" data from state.

data = map_observables(observables, state)

calculates the data for a tuple of observables applied to state. For a single observable (tuple of length 1), simply return the result of map_observable.

For multiple observables, return the tuple resulting from applying map_observable for each observable. If the tuple is "uniform" (all elements are of the same type, e.g. if each observable calculates the expectation value of a Hermitian operator), it is converted to a Vector. This allows for compact storage in a storage array, see init_storage.

Place data into storage for time slot i.

write_to_storage!(storage, i, state, observables)

For a storage array created by init_storage, store the data obtains from map_observables into the storage for time slot i. This delegates to the more general

write_to_storage!(storage, i, data)

Conceptually, this corresponds roughly to storage[i] = data, but storage may have its own idea on how to store data for a specific time slot. For example, with the default init_storage Vector data will be stored in a matrix, and write_to_storage! will in this case write data to the i'th column of the matrix.

For a given type of storage and data, it is the developer's responsibility that init_storage and write_to_storage! are compatible.