Propagation Methods¶

QDYN has a build-in propagation routine, which allows for propagating a state or a density matrix. The following page will explain

How to propagate (using QDYN)
Which method to use?
Methods implemented
Bibliography

How to propagate (using QDYN)¶

Basics ¶

The main module for propagation is implemented in the prop module. For standard propagation, particularly the routines prop() and init_prop() are important, which are described in the following.

Initializing the propagation¶

Before being able to propagate, one needs to initialize the propagation by:

call init_prop(...).

The function requires some parameters, which are listed and described in the corresponding subroutine description. The most important parameters are explained here:

gen:: This variable contains the so called dynamical generator. It contains the Hamiltonian or the Liouvillian of the system and should be initialized before, by e.g. calling the init() routine.
grid:: The spatial grid can also be initialized using the init() routine.
prop_work:: This is initialized by init_prop() and used as a temporary storage variable.
para:: These are the config file parameters, which can be read from the config file using read_para(). Most important is the variable of the para%prop, which contains the information on the propagation. For further details see How to set the method.
pulses(:):: This is the array of pulses needed for the propagation, which should be linked to their corresponding generators in gen. Make sure this is done by, for example, using a proper config file and using the init(). Even for a time-independent generator an empty array must be created for this purpose.

The method, which should be used for the propagation, can be set by specifying the method attribute of the function. See How to set the method for details.

Propagation¶

After having initialized the propagation, one can propagate the state by using

call prop(...).

This function also requires some parameters, which are listed and described in the corresponding subroutine description. In addition to the paramters already explained in Initializing the propagation, the state, which needs to be propagated needs to be given in the prop() routine:

state:: State to be propagated. This can for example also be set using the config file.

If you call the function with just the necessary parameters, the variable state will be changed to the propagated one. However, if you want to obtain the state at any time step you need to specify the argument storage. This is of type prop_seg_storage_t and is the variable into which the states at any times are stored. The variable of type prop_seg_storage_t is initialized by the init_prop_seg_storage() routine. For further details on how to retrieve the data from this type, see the prop_seg_storage module.

The method, which should be used for the propagation, can be set by specifying the method attribute of the function. See How to set the method for details.

The update and prop info hook¶

The prop() routine has two optional external interfaces, which can be used to interfere with the propagation process.

prop_info_hook:: The first interface is used for analysis purposes. It should mainly be used to print and analyze the data in the propagation process. However, it can also be used to manipulate the propagated state.
upd_hook:: The second interface is used for updating the dynamical generator and the grid. This can be used, for example, to realize moving grids.

For the structure of the required interfaces see the documentation of the prop() routine.

How to set the method (and further propagation properties)¶

In general there are three ways to set the propagation method and further properties for each method. The individual properties are explained in the corresponding methods section of the Methods implemented chapter.

The first and by far the most comfortable option is to set the method (and all the other properties) of the propagation in the config file (especially prop). As an example:
prop: method='cheby'
which sets the method to be the Chebyshev propagation. For a detailed list on all the available methods, see the section Which method to use? or the documentation of the type prop_pt. Also further example config files can be found in the examples folder, which is contained in QDYN.
The second way of setting the method (and further properties) is to add them, when calling
call init_prop(...,method='<method_to_use>',...)
or
call prop(...,method='<method_to_use>',...).
(substitute <method_to_use> with the method you want to use). Although this works well for setting only the method, nearly none of the method specific properties can be set this way!
The last and probably most cumbersome way of setting the propagation variables, is to manually define a variable of type para_t and also manually set the element para_t%prop with type prop_pt. This is not recommended!

How to propagate inhomogeneous equations?¶

The propagation of inhomogeneous equations is similar to the propagation of homogeneous equations of motion. The main difference here is the usage of the inhom_prop module instead of the prop module. Further, some additional parameters can be set in the config file (or directly in the para_t%prop variable).

Initializing the inhomogeneous propagation¶

Similar to the standard propagation, the inhomogeneous propagation needs to be initialized by

call init_inhom_prop(...).

The main arguments of the routine are similar to the init_prop() routine described in detail in the section Initializing the propagation. For a detailed list of all the arguments see the doumentation of the init_inhom_prop() routine.

Inhomogeneous propagation¶

The inhomogeneous propagation is also very similar to the homogeneous one:

call inhom_prop(...).

However, in addition to the parameters from the prop routine, the inhomogeneity needs to be specified. Therefore, the function has two required arguments and a required external interface. The arguments are

inhom_states: This argument of type prop_seg_storage_t contains the inhomogeneities at each time step. This variable is in general “filled” by a propagation, but can also be filled manually (demanding). See the prop_seg_storage module for more information.
lambda_b: This parameter is intended to be used to determine the strength of the inhomogeneity. However, it is just used as parameter for the external interface. Although it can therefore hypothetically be used as additional parameter in the external interface, this might cause problems in future versions of QDYN, due to the implementation of lambda_b as an inhomogeneity strength parameter only.

An additional external interface is required for the inhomogeneous propagation. This interface is called get_inhom_phi and allows to modify the inhom_states variable given to the routine. For the expected structure of the interface see the documentation of inhom_prop().

For further arguments and interfaces see the documentation of inhom_prop() and/or the explanation on the arguments of the homogeneous (init_)prop routine.

Additional config file parameters for the inhomogeneous propagation¶

For the inhomogeneous propagation, there exist some more parameters, which can be set in the config file.

inhom_expan_err:: Absolute error for coefficients in expansion series within inhomogeneous and ITO propagation; coefficient series will be truncated onces falling below this value (default 1.0d-12).
n_taylor:: Order of Taylor expansion in the inhomogeneous propagation for small arguments (default is 20), C.f. [17, 18].

Which method to use?¶

Methods implemented ¶

QDYN provides various methods to propagate states or density matrices. In the following, a detailed description of all the implemented propagation methods is given. If you don’t know which method is the right one for your purposes, see the “Which method to use?” decision tree.

Exact ¶

The method exact can be chosen by setting the method attribute to exact.

The calculations are carried out by the routine prop_exact().

Why to choose (and why not)¶

As the name suggests, this method will give the most accurate results. However, due to the diagonalization used during the process, this method is also by far the most expansive one. While this does not matter for problems with only a few Hilbert space dimensions, it leads to a very slow propagation with larger systems.

Idea¶

The idea of the exact method is to diagonalize the Hamiltonian \(\mathcal{H}(t)\) at every time \(t_n\) and propagate it via

\[\ket{\psi(t_{n+1})} = \sum_\lambda e^{i\lambda\cdot\mathrm{d}t} \ket{v_\lambda}\bra{v_\lambda} \ket{\psi(t_{n})},\]

where \(\ket{v_\lambda}\) are the eigenvectors of \(\mathcal{H}(t_n)\) to the according eigenvalue \(\lambda\).

Homogeneous Cheby ¶

The method cheby can be chosen by setting the method attribute to ‘cheby’.

The calculations are carried out by the routines in the module qdyn_cheby_mod.

Why to choose (and why not)¶

The Chebyshev propagator is in general used, due to the fact, that the expansion in Chebyshev polynomials lead to the fastest converging polynomial series expansion of a function [19]. Because of this fast convergence, it is the best choice to express the function

\[\hat{\mathcal{H}} = \frac{i}{\hbar} \mathcal{H} \mathrm{d}t\]

in a series expansion. This representation of the function represents the propagation step of the linear Schrödinger or non-dissipative Liouville-von Neumann equation (or any other linear/homogeneous equation of motion).

However, as prerequisite for using this expansion, the operator for of the homogeneous equation of motion (i.e. the Hamiltonian or the Liouvillian) needs to be hermitian. This is not the case for a general dissipative Liouvillian.

Basic idea¶

The idea behind the Chebyshev propagation is to expand the propagation step of, for example the Schrödinger equation

\[\ket{\Psi(t_{n+1})} = e^{-\frac{i}{\hbar}\mathcal{H}\mathrm{d}t} \ket{\Psi(t_{n})},\]

in a series of Chebyshev polynomials

\[e^{-\frac{i}{\hbar}\mathcal{H}\mathrm{d}t} \ket{\Psi(t_{n})} \approx \sum^{N-1}\ *{n} a_n P_n(-i\hat{\mathcal{H}}*\ {\mathrm{norm}}) \ket{\Psi(t_{n})},\]

which is truncated at certain order \(N\) [20]. Since the Chebyshev polynomials \(P(x)\) are only defined in the interval \(x\in[-1,1]\) the Hamiltonian needs to be normalized to \(\hat{\mathcal{H}}_{\mathrm{norm}}\).

For a detailed description see chapter 3.2.1 of reference [20].

Inhomogeneous Cheby ¶

If your equation of motion is inhomogeneous, you need to set the inhom_method attribute to cheby.

To set the maximum order of the expansion of the inhomogeneous term, use the parameter prop:inhom_max_order in the config file (default is 5). For details on how to do this, see the How to set the method section.

Basic idea¶

For the idea of the inhomogeneous Chebyshev propagator, see the documentation of the inhom_cheby module, which contains the subroutines for the calculations of the inhomogeneous Chebyshev propagator.

Homogeneous Newton ¶

The method newton can be chosen by setting the method attribute to newton.

The calculations are carried out by the routines in the module qdyn_newton_mod.

Using newton in QDYN¶

For using the newton propagator, there are four parameters in the prop_pt variable (for a description on how these can be set, see the How to set the method section), which can be used to modify the behavior of the algorithm

newton_relerr: Minimum Norm of the solution of the Newton series in relation to the correct one. Determines when the series is truncated (default is 1.0d-13).
newton_norm_min: Minimum norm of a state before the state is assumed to be identical to zero. This primarily affects the detection of eigenstates. You should only change this (and probably only to a lower value) if the Newton propagator gives wrong results due to assuming a state is an eigenstate when it is not (default is 1.0d-15).
newton_arnoldi_order: Dimension of the Hessenberg matrix to be constructed in each restart of the Newton propagator (default is 4).
newton_max_restarts: Maximum number of restarts in the Newton propagator (default is 1000).

The init_prop() routine has a `state0`-argument, which gives the average state for the newton propagation and is necessary for the initialization of the propagator. If this is not specified, a random state is used. For more detail see the argument description in the init_prop() routine.

Why to choose (and why not)¶

The Newton propagator is in general used to propagate systems with dissipative Liouvillians or in general non-hermitian Hamiltonians. This makes it a universal propagator for systems with linear equations of motion. However, if the Hamiltonian is hermitian or the Liouvillian is non-dissipative, the Chebyshev-propagator performs better.

Basic idea¶

The basic idea of the Newton propagation is to use Newton polynomials, to expand the propagation step of, for example the Liouville-von Neumann equation

\[\ket{\Psi(t_{n+1})} = e^{-\frac{i}{\hbar}\mathcal{L}\mathrm{d}t} \ket{\Psi(t_{n})},\]

in a series of Newton polynomials

\[e^{-\frac{i}{\hbar}\mathcal{H}\mathrm{d}t} \ket{\Psi(t_{n})} \approx \sum^{N-1}_{n} a_n R_n(\mathcal{L}) \ket{\Psi(t_{n})}.\]

For a detailed description see chapter 3.2.2 of reference [20].

Inhomogeneous Newton ¶

If your equation of motion is inhomogeneous, you need to set the inhom_method attribute to newton.

Using inhom_newton in QDYN¶

The usage of the inhomogeneous inhom_newton() propagation is similar to the usage of the homogeneous newton propagation. However, the inhom_prop() routine, needs to be used for the inhomogeneous propagation. On how to use this function, see the How to propagate inhomogeneous equations? section for more details.

Basic idea¶

For the idea of the inhomogeneous Newton propagator, see the documentation of the inhom_newton module, which contains the subroutines for the calculations of the inhomogeneous Newton propagator.

ITO ¶

Iterative-Time-Ordering (ITO) is a method which accounts for the time ordering effects and is able to propagate under rapidly oscillating fields.

The current implementation of ITO is called nito (Newton/new - ITO) and can be chosen by setting the method attribute nITO.

The calculations are carried out by the routines in the module qdyn_nito_mod.

Using ITO in QDYN¶

For using the ITO propagator, it can be chosen, whether the solution of the time-dependent Schrödinger equation is solved iteratively, this can be controlled by the parameter in the prop_pt-variable (for a description on how these can be set, see the How to set the method section)

ITO_do_iter If .false., iterative algorithm in ITO propagation is omitted and convergence after one iteration step is assumed (default is .true.).

In addition to that, the inhomogeneous expansion can be controlled by the parameters

inhom_max_order: Maximum order M in the expansion of the inhomogeneous term in the inhomogeneous Chebyshev or the ITO propagation (default is 5).
inhom_expan_err: Absolute error for coefficients in expansion series within inhomogeneous and ITO propagation; coefficient series will be truncated onces falling below this value (default is 1.0d-12).

and since Newton polynomials are used for the expansion, the parameters of the prop_pt-variable used to set the Newton-method also affect the algorithm

newton_relerr: Minimum Norm of the solution of the Newton series in relation to the correct one. Determines when the series is truncated (default is 1.0d-13).
newton_norm_min: Minimum norm of a state before the state is assumed to be identical to zero. This primarily affects the detection of eigenstates. You should only change this (and probably only to a lower value) if the Newton propagator gives wrong results due to assuming a state is an eigenstate when it is not (default is 1.0d-15).
newton_arnoldi_order: Dimension of the Hessenberg matrix to be constructed in each restart of the Newton propagator (default is 4).
newton_max_restarts: Maximum number of restarts in the Newton propagator (default is 1000).

Why to choose (and why not)¶

As mentioned above, the main reason to use ITO, is due to its ability to propagate with very high accuracy for rapidly changing pulses. However, this comes along with rather high numerical costs, which makes it unsuitable if you have sufficiently smooth pulses.

Idea¶

The main idea of the ITO method is to subdivide the time interval and solve the time dependent Schrödinger equation, which is turned into a inhomogeneous equation of motion, self-consistently and therefore also accounts for time ordering effects.

For further explanations of the idea behind ITO, see the documentation of the nito module and the paper referred to there.

Homogeneous Runge-Kutta 45 ¶

The method rk45 can be chosen by setting the method attribute to rk45.

The calculations are carried out by the routines in the module qdyn_rk45_mod.

Using rk45 in QDYN¶

For using rk45 you need to specify a function, which returns the derivative of the state (c.f. \(f(t,y(t))\) in the equation in the ‘Why to choose (and why not)’ section below). This function can be passed by an external interface rk_deriv, which is given in the description of rk45() under external interfaces. The interface rk_deriv is then given to the prop() routine. Several examples of these are provided by the rk45 module.

In addition, there are two parameters in the prop_pt-variable (For a description on how these can be set, see the How to set the method section)

rk45_relerr
rk45_abserr

These two variables give the highest relative and absolute error allowed, respectively. By default these values are \(10^{-11}\) and \(10^{-15}\), respectively.

Why to choose (and why not)¶

The Runge-Kutta-Fehlberg (rk45) method provides a general algorithm for solving ordinary differential equations. More specifically it gives a numerical solution to the first order differential equation

\[y’(x) = f(t,y(t)), \qquad\qquad \qquad\qquad y(t_0)=x_0.\]

The advantage of this method is, that nearly any equation of motion in the form of the equation above can be solved, such as for example the Gross-Pitaevskii equation, which has a non-linear term in the equation of motion.

However, with the possibility of solving general problems, there is a loss of efficiency and accuracy. Generally speaking, the Runge-Kutta-Fehlberg method accumulates the inevitable errors of each time step. Therefore, this method of propagation should only be used if you have a non-linear equation of motion. Otherwise, the expansion in polynomials, such as Chebyshev- or Newton-polynomials should be used, which distribute the error over the whole interval. See the decision tree to choose another method.

Inhomogeneous Runge-Kutta 45 ¶

If your equation of motion is inhomogeneous, you need to set the inhom_method attribute to rk45.

The calculations are carried out by the routines in the module qdyn_inhom_rk45_mod.

Using inhom_rk45 in QDYN¶

The usage of the inhomogeneous rk45 propagation is similar to the usage of the homogeneous rk45 propagation. However, the rk_deriv external interface is given to the inhom_prop() routine, which is used for inhomogeneous propagation. On how to use this function, see the How to propagate inhomogeneous equations? section for more details.

Idea¶

Runge-Kutta methods are a generalization of the Euler-method, which recursively solves the first order equation of motion, via

\[y(t_{n+1}) = y(t_{n})+ \Delta t\cdot f(t_n,y(t_n)),\]

which means, that one follows the gradient at of \(y(t)\) at \(t_{n}\).

In contrast to this, the idea of the Runge-Kutta methods is to use auxiliary steps to evaluate the function and use the gradient information to refine the approximation for the overall step from \(y(t_n)\) to \(y(t_{n+1})\).

The Runge-Kutta-Fehlberg algorithm (rk45) is an embedded method, that is, an algorithm that combines two Runge-Kutta methods to estimate the local truncate error, allowing to control it by using an adaptive stepsize. The methods used have errors of order \(O (h^4)\) and \(O (h^5)\), thus the acronym rk45.

Bibliography ¶

[17] Ndong et al, J. Chem. Phys. 130, 124108 (2009)
[18] Ndong et al, J. Chem. Phys. 132, 064105 (2010)
[19] A. Gil, J. Segura, and N. Temme, Numerical Methods for Special Functions, Society for Industrial and Applied Mathematics, 2007
[20] M. Goerz, Optimizing Robust Quantum Gates in Open Quantum Systems, Ph.D. Thesis. Universität Kassel, 2015