Manual:Time Dependent

From OctopusWiki

Jump to: navigation, search

Contents

[edit]

Time evolution

When CalculationMode = td, the code performs the time propagation of the electronic orbitals and – if required – of the ionic positions. The 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 \psi_j(t)\,\!:

i{\partial \over \partial t}\psi_j(t) = H(t)\psi_j(t)\,\!

\psi_j(t=0) = \psi_j^{(0)}\,\!

Because this equation is linear (the time derivative and the Hamiltonian are both linear operators), one may formally define a linear “evolution” operator, U(T, t)\,\!, which transforms the initial vector into the solution at time T\,\!:

\psi_j(T) = U(T, 0)\psi_j^{(0)}

Moreover, there is the formal exact expression for the evolution operator

\psi_j(T) = \mathcal{T}\!\!\exp\left\{ -i\!\!\int_0^{T}d\tau H(\tau)\right\} \psi_j^{(0)}\,\!,

where \mathcal{T}\!\!\exp\,\! is the time-ordered exponential, which is a short-hand for:

\psi_j(T) = \left\{\sum_{n=0}^{\infty} \frac{\left(-i\right)^n}{n!} \int_0^t d\tau_1 \cdots \int_0^t d\tau_n H(\tau_1) \cdots H(\tau_n)\right\}\psi_j^{(0)}\,\!

If the Hamiltonian commutes 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 simply written as:

\psi_j(T) = \exp\left\{ -iTH\right\} \psi_j^{(0)}\,\!.

Unfortunately, this is not the case for TDDFT when the system is exposed to external time-dependent perturbations like electric and magnetic fields or pulsed lasers. But even without an external time-dependency, there remains the intrinsic time-dependency of the Kohn-Sham Hamiltonian, which is built “self-consistenly” from the varying electronic density.

The first step to tackle this problem is to split the propagation of the long intervall [0, T]\,\! into N\,\! smaller steps by utilizing the group property

U(T, t) = U(T, t')U(t', t)\,\!

of the time evolution operator. This yields the following time discretization:

U(T, 0) = \prod_{i=0}^{N-1}U(t_i+\Delta t, t_i)\,\!,

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:

\psi_j(t+\Delta t) = U(t+\Delta t, t)\psi_j(t) = \mathcal{T}\!\!\exp\left\{ -i\int_{t}^{t+\Delta t}d\tau H(\tau)\right\} \psi_j(t)\,\!.

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_{\mathrm{max}}\,\! is the maximum frequency that we want to discern, \Delta t\,\! should be no larger than \approx 1/\omega_{\mathrm{max}}\,\!. Below this \Delta t_{\mathrm{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 need 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 TDEvolutionMethod – and, if this algorithm requires it, we will also need an algorithm to approximate the exponential of a matrix operator \exp\left\{ A\right\}\,\! – which will be specified by variable TDExponentialMethod.


[edit]

Propagators for the time dependent Kohn-Sham equations

In the following, we describe the propagators for the time-dependent Kohn-Sham equations, which are available in Octopus. These propagators solve the problem of approximating the orbitals \psi_j(t+\Delta t)\,\! by the knowledge of \psi_j(\tau)\,\! and H(\tau)\,\!for 0\le\tau\le t\,\!. Some methods require the knowledge of the Hamiltonian at some points \tau\,\! in time between t\,\! and t+\Delta t\,\!.

This quantity can be approximated by the following procedure:

  1. Approximate H(\tau)\,\! through extrapolation by a polynomial fit to a certain number of previous time steps.
  2. Propagate \psi_j(t)\,\! to get \psi_j(t+\Delta t)\,\!.
  3. Calculate H(t+\Delta t)\,\! from the orbitals \psi_j(t+\Delta t)\,\!.
  4. Interpolate the required H(\tau)\,\! from H(t)\,\! and H(t+\Delta t)\,\!.
  5. Repeat the steps 2 to 4 until self-consistency is reached.

In Octopus, however, the above scheme is dropped for performance reasons and only step 1 is implemented by a second order extrapolation, except for the first two steps, where the extrapolation obviously cannot be trusted. Instead, we rely on a sufficiently small \Delta t\,\!.

[edit]
Splitting techniques

Splitting techniques exploit the fact that the Hamiltonian is composed of two terms, H=T+V_{\mathrm{KS}}(t)\,\!, the first kinetic term being diagonal in Fourier space, and the second potential term diagonal (or almost diagonal) in real space.

The Split Operator (TDEvolutionMethod=split) is defined as

U_{\mathrm{SO}}(t+\Delta t, t) = \exp\left\{-i\frac{\Delta t}{2}T\right\}\exp\left\{-i\Delta t V'_{\mathrm{KS}}(t)\right\}\exp\left\{-i\frac{\Delta t}{2} T\right\}\,\!

with

V'_{\mathrm{KS}}(t) = v_{\mathrm{ext}}(t+\Delta t/2) + v_H[n'](t) + v_{\mathrm{xc}}[n'](t)\,\!.

In the above expression, the density n'\,\! is calculated after applying the first kinetic term of U_{\mathrm{SO}}(t+\Delta t, t)\,\! to the orbitals \psi_j(t)\,\!. Only after calculation of the new potential V'_{\mathrm{KS}}(t)\,\! the potential term and second kinetic term are applied to yield \psi_j(t+\Delta t)\,\!.

The two kinetic terms are evaluated in Fourier space by the aid of the Fast Fourier Transform algorithm.

The Split Operator method in this realization is of second order in \Delta t\,\!.

The Suzuki-Trotter method (TDEvolutionMethod=suzuki_trotter) is a fourth order propagator defined by

U_{\mathrm{ST}}(t+\Delta t, t) = \prod_{k=1}^5\exp\left\{-i p_k\frac{\Delta t}{2}T(t_k)\right\}\exp\left\{-i p_k\Delta t V(t_k)\right\}\exp\left\{-i p_k\frac{\Delta t}{2}T(t_k)\right\}\,\!.

The p_k\,\! are a suitably chosen set of real numbers related to the t_k\,\! by t_k = t + \left(p_1 + \cdots = p_k/2\right)\Delta t\,\!. The potential at times t_k\,\! is extrapolated in order to get a proper fourth order method. This method allows for comparably large timesteps but each time step takes roughly five times longer than for the second order Split Operator.

[edit]
Midpoint rules

The implicit midpoint rule or Crank-Nicholson method (TDEvolutionMethod=crank_nicholson) calculates the exponential of the Hamiltonian by a first order Padé approximation. The Hamiltonian is evaluated at t + \Delta t/2\,\! and the integral is dropped:

U_{\mathrm{CN}}(t + \Delta t, t) = \frac{1- i \frac{\Delta t}{2}H(t+\Delta t/2)}{1+ \frac{\Delta t}{2}H(t+\Delta t/2)}\,\!

The calculation of the matrix fraction is transformed into the solution of the linear system

L \psi_j(t+\Delta t) = b\,\!

with the known quantities

L = 1+i\frac{\Delta t}{2}H(t+\Delta t/2)\,\!

and

b = \left\{1-i\frac{\Delta t}{2}H(t+\Delta t/2)\right\}\psi_j(t).\,\!

The Crank-Nicholson scheme is unitary and preserves time-reversal symmetry.

Likewise, exponential midpoint rule (TDEvolutionMethod=exp_mid) is unitary and time-reversal symmetry preserving but simpler with the drawback that it requires fairly small time steps:

U_{\mathrm{EM}}(t + \Delta t, t) = \exp\left\{-i\Delta t H(t + \Delta t/2)\right\}\,\!

[edit]
Magnus expansions

Magnus expansion(TDEvolutionMethod=magnus) is the most sophisticated propagation implemented in Octopus. According to Magnus, there exists an operator \Omega(t+\Delta t, t)\,\! for which holds

U(t+\Delta t, t) = \exp\left\{\Omega(t+\Delta t, t)\right\}\,\!

given by the series

\Omega(t + \Delta t, t) = \sum_{k=1}^\infty \Omega_k(t + \Delta t, t)\,\!

that converges at least for some local environment of t\,\!.

The operators \Omega_k(t+\Delta t, t)\,\! are generated with the recursive procedure

\Omega_k(t+\Delta t, t) = \sum_{l=0}^{k-1} \frac{B_l}{l!}\int_t^{t+\Delta t} S^l_k(\tau)d\tau,\,\!

S_1^0(\tau) = -iH(\tau), \quad S_k^0=0\ \mathrm{for}\ k>1,\,\!

S_k^j(\tau) = \sum_{m=1}^{k-l}\left[\Omega_m(t+\Delta t, t), S_{k-m}^{l-1}(\tau)\right]\ \mathrm{for}\ 1\le l\le k-1,\,\!

B_l\,\! Bernoulli numbers.

In Octopus we have implemented a fourth order Magnus expansion, which means that the operator series is truncated to second order and the integrals are calculated by a second order quadrature formula. The resulting operator is

\Omega_{M(4)}(t + \Delta t, t) = -i\frac{\Delta t}{2}\left[H(t_1)+ H(t_2)\right] - \frac{\sqrt{3}\Delta t^2}{12}\left[H(t_2), H(t_1)\right]\,\!

with the quadrature sampling points t_{1,2} = t + \left[\frac{1}{2}\mp \frac{\sqrt{3}}{6}\right]\Delta t\,\!.

With a modified Hamiltonian

H_{M(4)}(t, \Delta t) = \overline{H}(t, \Delta t) + i \left[T + v_{\mathrm{ext}}^{\mathrm{nonlocal}}(t), \overline{\Delta V}(t, \Delta t)\right],\,\!

\overline{H}(t, \Delta t) = T + \frac{1}{2}\left[V_{\mathrm{KS}}(t_1)+V_{\mathrm{KS}}(t_2)\right],\,\!

\overline{\Delta V}(t, \Delta t) = \frac{\sqrt{3}}{12}\Delta t\left[V_{\mathrm{KS}}(t_1)-V_{\mathrm{KS}}(t_1)\right]\,\!

the propagator takes the form

U_{M(4)}(t+\Delta t, t) = \exp\left\{-i \Delta t H_{M(4)}(t, \Delta t)\right\}.\,\!

Note, that only the nonlocal components of the Kohn-Sham Hamilton contribute to the commutator and, furthermore, that we make the assumption that the nonlocal potential v_{\mathrm{ext}}^{\mathrm{nonlocal}}(t)\,\! does not vary significantly in the interval [t, t+\Delta t]\,\!. This nonlocal component stems from the ionic pseudopotentials. Consequently, its variation is caused by the ionic movement, which is negligible on the electronic time scale determining \Delta t\,\!.

[edit]
Time reversal symmetry based propagation

Due to time reversal symmetry, propagating backwards by \Delta t/2\,\! starting from \psi_j(t+\Delta t)\,\! should lead to the same result as propagating forwards by \Delta t/2\,\! starting from \psi_j(t)\,\!. Using the simplest approximation to the evolution operator, we end up with the condition

\exp\left\{+i\frac{\Delta t}{2}H(t +\Delta t)\right\} \psi_j(t+\Delta t) = \exp\left\{-i\frac{\Delta t}{2}H(t)\right\} \psi_j(t)\,\!

Rearrangeing the terms gives an approximation for the propagator:

U_{\mathrm{ETRS}}(t + \Delta t, t) = \exp\left\{-i \frac{\Delta t}{2} H(t + \Delta t)\right\} \exp\left\{-i \frac{\Delta t}{2}H(t)\right\}\,\!

This enforced time reversal symmetry method is available in Octopus by setting (TDEvolutionMethod=etrs).

It is worthwhile to give a sidenote on the approximation of H(t+\Delta t)\,\!: for the etrs method the Hamiltonian at time t+\Delta t\,\! is not extrapolated by build from the density n' = \sum_j|\psi'_j(t+\Delta t)|^2\,\!, which, in turn, is calculated from the orbital estimate

\psi'_j(t+\Delta t) = \exp\left\{-i\Delta t H(t)\right\}\psi_j(t).\,\!

There exists, however, also a variant TDEvolutionMethod=aetrs (approximated enforced time reversal symmetry) that performs a second order polynomial extrapolation of H(t-k\Delta t)\,\!, k=0,1,2\,\!, to get H(t+\Delta t)\,\!, which is abuot 40 % faster than etrs.

[edit]

Approximations to the exponential of an operator

Most of the evolution methods described in the previous section require the calculation of the exponential of a matrix. Before going into the details of the TDExponentialMethod that Octopus provides, one general difficulty deserves mentioning.

In principle, we would like to calculate \exp\left\{A\right\}v\,\! by first calculating \exp\left\{A\right\}\,\! and then applying this result to any vector v\,\!. Unfortunately, the size of our Hamiltonian matrix is of the order \approx 10^5\,\! which forbids its full storage in matrix form. Apart from that, methods to calculate the exponential of a matrix explicitly are limited to a few thousand matrix elements. For this reason, we have to turn to iterative solutions that calculate \exp\left\{A\right\}v\,\!, for a particular vector v\,\!.

[edit]
Polynomial expansions

The simplest possibility is to use the definition of the exponential of a matrix

\exp\left\{A\right\} = \sum_{k=0}^\infty \frac{1}{k!}A^k\,\!

to get the approximation of order k\,\!(TDExponentialMethod=standard):

\mathrm{taylor}_N\left\{A, v\right\} = \sum_{k=0}^N\frac{1}{k!}A^kv\,\!

\mathrm{taylor}_N\left\{A, v\right\}\,\! amounts to the expansion of the exponential in the standard polynomial base \left\{1, x, x^2, \ldots\right\}\,\!.

Experience shows that the choice k=4\,\! gives particularly good results for our TDDFT implementation.

The standard polynomial base is not the only possibility, Octopus also implements the expansion of the exponential in the Chebychev base:

\mathrm{cheb}_N\left\{A, v\right\} = \sum_{k=0}^Nc_kT_k(A)kv\,\!

with T_k\,\! being the Chebychev polynomial of order k\,\!.

The truncation N\,\! for both these expansions can be set by the input variable ExpOrder.

[edit]
Krylov subspace projection

The Krylov subspace of order N\,\! for a given Operator A\,\! and vector v\,\! is defined as

\mathcal{K}_N\left\{A, v\right\} = \mathrm{span}\left\{v, Av, A^2v, \ldots, A^{N-1}v\right\}.\,\!

The idea is to find the element of \mathcal{K}_N\left\{A, v\right\}\,\! that optimally approximates \exp\left\{A\right\}\,\! in a least square sense. It can be proved that this is given by

\beta V_N(V_N^T\exp\left\{A\right\}V_N)e_1\,\!

with V_N = [v_1, \ldots, v_N]\,\! being an orthonormal base of \mathcal{K}_N\left\{A, v\right\}\,\! calculated with the Lanczos procedure. To get rid of the exponential, the approximation

V_N^T\exp\left\{A\right\}V_N \approx \exp\left\{V_N^TAV_N\right\}\,\!

is introduced. The object H_N=V_N^TAV_N\,\! is also a result of the Lanczos procedure and of the small dimension k\times k\,\! which allows for direct calculation of the remaining exponential. So, we have

\mathrm{lanczos}_N\left\{A, v\right\} = V_N\exp\left\{H_N\right\}e_1,\,\!

which can be selected by setting (TDExponentialMethod=lanczos). When actually calculating the exponential, Octopus recursively generates the quantities H_N\,\! and V_N\,\! until some convergence criterion – controlled by TDLanczosTol – is met or the order exceed the maximum allowed order controlled by . In the latter case, a warning is emitted. Care should be taken in chosing the convergence parameter small enough to avoid faulty evolutions.

[edit]
Splitting techniques

Splitting techniques can also be used to calculate the exponential of the Hamiltonian (TDExponentialMethod=split) which work in the same way as described in the previous section for the propgator. Higher order splitting schemes (TDExponentialMethod=suzuki_trotter) are also available.

[edit]

Performing a time propagation

To actually perform a time propagation, a necessary prerequisite is to have the ground state density (see Ground State of the system under consideration.

Then, the CalculationMode has to be set to td.

Besides the already mentioned input variables concerning the propagation method, two more important parameters have to be set: TDTimeStep and TDMaximumIter. The first one sets \Delta t\,\!, the second sets the number of iterations. It is convenient to define these two variables like this: 

T  = 0.1     # Length of propagation.
dt = 0.002   # Length od time step.

TDMaximumIter = T/dt
TDTimeStep = dt

Without perturbing the system, its total energy should be a constant. This check to determine a reasonable time step has to be done before any other calculation.

The relevant parts of Octopus' output might look like 

  Iter           Time        Energy     Elapsed Time
      1       0.002000   -541.542073         4.800
      2       0.004000   -541.542073         4.560
      3       0.006000   -541.542073         3.880
      .
      .
      .
     48       0.096000   -541.542073         7.420
     49       0.098000   -541.542073         7.620
     50       0.100000   -541.542073         7.490

The energy is printed in the third column and should remain constant. In general, larger time steps are desireable to shorten computational propagation time but keep in mind that the size of the time step is related to the maximum frequency you will be able to observe, and, consequently, also to any external perturbation you might wish to expose the system to.

[edit]

External fields

[edit]
Delta kick - Calculating an Absorption Spectrum

To obtain the linear optical absoption spectrum of the system, we follow the scheme proposed by Yabana and Bertsch, and excite all frequencies of the system by giving some small momemtum ({\mathcal K}) to the electrons. This is achieved by transforming the ground-state wave-functions according to:

\varphi_i(r, \delta t) = e^{i {\mathcal K} z} \varphi_i(r, 0) \,,

and then propagating these wave-functions for some (finite) time. The spectrum can then be obtained from the expression of the dipole strength function S(ω):

S(\omega) = \frac{2 \omega}{\pi} {\mathcal Im} \alpha(\omega) \,,

where the dynamical polarizability, α(ω), is essentially the Fourier transform of the dipole moment of the system d(t):

\alpha(\omega) = \frac{1}{\mathcal{K}}\int\!\! dt \; e^{i\omega t} \left[d(t)-d(0)\right]\,.

With this definition, the Thomas-Reiche-Kuhn f-sum rule for the number of electrons, N, is given by the integral:

N = \int\!\! d\omega\;S(\omega)\,.

This sum rule can be used to check the quality of the calculations. Another check is energy conservation, which the TDDFT respects when no external field is applied.

To obtain a spectrum with such recipe in octopus one has to follow the steps:

  1. Choose a system of independent (can be non-orthogonal) axes (defined by the variable TDPolarization). The default are the 3 Cartesian axes.
  2. Add an electric kick using the variable TDDeltaStrength. Its value should be small so that we remain in the linear response, but large enough to avoid numerical errors.
  3. Run a time-dependent simulation for the 3 independent directions. This involves running Octopus 3 times changing the value of TDPolarizationDirection. After each run, move the file td.general/multipoles to td.general/multipoles.X, where X is 1, 2, or 3.
  4. Run the utility oct-cross-section
[edit]
Lasers

For the purpose of obtaining nonlinear optical properties, we follow the evolution of the system under the influence of a laser field that is treated in the dipole approximation (although this constraint can be removed). The emitted harmonic spectra can then be calculated from the acceleration of the dipole moment\cite{nlo}:

H(\omega) \propto \left| \int\!\! dt \; e^{i\omega t} \frac{d^2}{dt^2} d(t) \right|^2\,.

During the propagation, charge density is absorbed at the boundaries of the simulation region, either by an imaginary absorbing potential or a mask function. In the first case, we add to the Kohn-Sham potential:

V_{\rm eff}(r, t) = V_{\rm KS}(r,t) - i V_{\rm abs}(r) \,,

where Vabs is zero in the inner region of the simulation box, and rises smoothly till the edges. By adjusting both the height and the shape of the potential, we can select which momenta are absorbed and prevent the unwanted reflections at the boundary. When using a mask, the wave-function is multiplied in each time-step by a function which is 1 in the inner simulation region and gradually goes to 0 at the borders:

\varphi(r, t) \rightarrow M(r) \varphi(r, t) \,.

The absorbed charge can be interpreted as an ionization probability and can be used to estimate the photo-electron spectra. The box size has to be big enough so that the physical system is not perturbed by the absorbing boundaries. Note that the wavefunctions are no longer normalized as the system slowly gets charged.

Previous Manual:Ground_State - Next Manual:Casida

Back to Manual

Retrieved from "http://www.tddft.org/programs/octopus/wiki/index.php/Manual:Time_Dependent"
Personal tools