When= td, the code performs the time propagation of the electronic orbitals and – if required – of the ionic positions. This latter task does not pose major algorithmical problems (the usual Verlet algorithms deal with that task); however the best way to propagate a Schrödinger-like equation is still unclear. Due to this fact, we provide with a rather excessive selection of possibilities for that purpose. Before describing the set of variables necessary to specify the way in which the time evolution is to be performed, it is worth making a brief introduction to the problem.
We are concerned with a set of Schrödinger-like equations for the electronic orbitals:
Being the equation linear, one may formally define a linear “evolution” operator, which trasforms the initial vector into the solution at time T:
Moreover, there is a formal exact expression for the evolution operator:
where T\exp is the time-ordered exponential. If the Hamiltonian conmutes with itself at different times, we can drop the time-ordering product, and leave a simple exponential. If the Hamiltonian is time-independent – which makes it trivially self commuting, the solution is then simply written as:
Unfortunately, this is not the case in general. We have to find an algorithm able to cope with time-dependent Hamiltonians, such as the self-consistent time-dependent Kohn-Sham operator, which is built “self consistently” from the varying electronic density.
The first step is to perform a time-discretization: the full propagation between 0 and T is decomposed as:
where t_0=0, t_N=T, \delta t = T/N. So at each time step we are dealing with the problem of performing the short-time propagation:
In this way, one can monitor the evolution in the interior of [0,t]. In fact, the possibility of monitoring the evolution is generally a requirement; this requirement imposes a natural restriction on the maximum size of \delta t: if \omega_\rm max is the maximum frequency that we want to discern, \delta t should be no larger than \approx 1/\omega_\rm max. Below this \delta t_\rm max, we are free to choose \delta t considering performance reasons: Technically, the reason for the discretization is twofold: the time-dependence of H is alleviated, and the norm of the exponential argument is reduced (the norm increases linearly with \delta t).
Since we cannot drop the time-ordering product, the desired algorithm cannot be reduced, in principle, to the calculation of the action of the exponential of an operator over the initial vector. Some algorithms tailored to approximate the evolution operator, in fact, do not even require to peform such operator exponentials. Most of them, however, do rely on the calculation of one or more exponentials, such as the ones used by octopus. This is why in principle we need to specify two different issues: the “evolution method”, and the “exponential method”. In other words: we need an algorithm to approximate the evolution operator U(t+\delta t, t) – which will be specified by variable– and, if this algorithm requires it, we will also need an algorithm to approximate the exponential of a matrix operator \exp\lbrace A\rbrace – which will be specified by variable .