## Time-Dependent

TDExternalFields
Section: Time-Dependent
Type: block

The block TDExternalFields describes the type and shape of time-dependent external perturbations that are applied to the system, in the form $$f(x,y,z) \cos(\omega t + \phi (t)) g(t)$$, where $$f(x,y,z)$$ is defined by by a field type and polarization or a scalar potential, as below; $$\omega$$ is defined by omega; $$g(t)$$ is defined by envelope_function_name; and $$\phi(t)$$ is the (time-dependent) phase from phase.

These perturbations are only applied for time-dependent runs. If you want the value of the perturbation at time zero to be applied for time-independent runs, use TimeZero = yes.

Each line of the block describes an external field; this way you can actually have more than one laser (e.g. a "pump" and a "probe").

There are two ways to specify $$f(x,y,z)$$ but both use the same omega | envelope_function_name [| phase] for the time-dependence. The float omega will be the carrier frequency of the pulse (in energy units). The envelope of the field is a time-dependent function whose definition must be given in a TDFunctions block. envelope_function_name is a string (and therefore it must be surrounded by quotation marks) that must match one of the function names given in the first column of the TDFunctions block. phase is optional and is taken to be zero if not provided, and is also a string specifying a time-dependent function.

(A) type = electric field, magnetic field, vector_potential

For these cases, the syntax is:

%TDExternalFields
type | nx | ny | nz | omega | envelope_function_name | phase
%

The vector_potential option (constant in space) permits us to describe an electric perturbation in the velocity gauge. The three (possibly complex) numbers (nx, ny, nz) mark the polarization direction of the field.

(B) type = scalar_potential

%TDExternalFields
scalar_potential | "spatial_expression" | omega | envelope_function_name | phase
%

The scalar potential is any expression of the spatial coordinates given by the string "spatial_expression", allowing a field beyond the dipole approximation.

A NOTE ON UNITS:

It is very common to describe the strength of a laser field by its intensity, rather than using the electric-field amplitude. In atomic units (or, more precisely, in any Gaussian system of units), the relationship between instantaneous electric field and intensity is: $$I(t) = \frac{c}{8\pi} E^2(t)$$.

It is common to read intensities in W/cm$$^2$$. The dimensions of intensities are [W]/(L$$^2$$T), where [W] are the dimensions of energy. The relevant conversion factors are:

Hartree / ($$a_0^2$$ atomic_time) = $$6.4364086 \times 10^{15} \mathrm{W/cm}^2$$

eV / ( Å$$^2 (\hbar$$/eV) ) = $$2.4341348 \times 10^{12} \mathrm{W/cm}^2$$

If, in atomic units, we set the electric-field amplitude to $$E_0$$, then the intensity is:

$$I_0 = 3.51 \times 10^{16} \mathrm{W/cm}^2 (E_0^2)$$

If, working with Units = ev_angstrom, we set $$E_0$$, then the intensity is:

$$I_0 = 1.327 \times 10^{13} (E_0^2) \mathrm{W/cm}^2$$

Options:

• electric_field: The external field is an electric field, the usual case when we want to describe a laser in the length gauge.
• magnetic_field: The external field is a (homogeneous) time-dependent magnetic field.
• vector_potential: The external field is a time-dependent homogeneous vector potential, which may describe a laser field in the velocity gauge.
• scalar_potential: The external field is an arbitrary scalar potential, which may describe an inhomogeneous electrical field.

TDFreezeDFTUOccupations
Section: Time-Dependent
Type: logical
Default: no

The occupation matrices than enters in the LDA+U potential are not evolved during the time evolution.

TDFreezeHXC
Section: Time-Dependent
Type: logical
Default: no

The electrons are evolved as independent particles feeling the Hartree and exchange-correlation potentials from the ground-state electronic configuration.

TDFreezeOrbitals
Section: Time-Dependent
Type: integer
Default: 0

(Experimental) You have the possibility of "freezing" a number of orbitals during a time-propagation. The Hartree and exchange-correlation potential due to these orbitals (which will be the lowest-energy ones) will be added during the propagation, but the orbitals will not be propagated.
Options:

• sae: Single-active-electron approximation. This option is only valid for time-dependent calculations (CalculationMode = td). Also, the nuclei should not move. The idea is that all orbitals except the last one are frozen. The orbitals are to be read from a previous ground-state calculation. The active orbital is then treated as independent (whether it contains one electron or two) -- although it will feel the Hartree and exchange-correlation potentials from the ground-state electronic configuration.

It is almost equivalent to setting TDFreezeOrbitals = N-1, where N is the number of orbitals, but not completely.

TDFreezeU
Section: Time-Dependent
Type: logical
Default: no

The effective U of LDA+U is not evolved during the time evolution.

TDFunctions
Section: Time-Dependent
Type: block

This block specifies the shape of a "time-dependent function", such as the envelope needed when using the TDExternalFields block. Each line in the block specifies one function. The first element of each line will be a string that defines the name of the function. The second element specifies which type of function we are using; in the following we provide an example for each of the possible types:

Options:

• tdf_cw:

%TDFunctions
"function-name" | tdf_cw | amplitude
%

The function is just a constant of value amplitude: $$f(t)$$ = amplitude

• tdf_gaussian:

%TDFunctions
"function-name" | tdf_gaussian | amplitude | tau0 | t0
%

The function is a Gaussian, $$f(t) = F_0 \exp( - (t-t_0)^2/(2\tau_0^2) )$$, where $$F_0$$ = amplitude.

• tdf_cosinoidal:

%TDFunctions
"function-name" | tdf_cosinoidal | amplitude | tau0 | t0
%

$$f(t) = F_0 \cos( \frac{\pi}{2} \frac{t-2\tau_0-t_0}{\tau0} )$$

If $$| t - t_0 | > \tau_0$$, then $$f(t) = 0$$.

• tdf_trapezoidal:

%TDFunctions
"function-name" | tdf_trapezoidal | amplitude | tau0 | t0 | tau1
%

This function is a trapezoidal centered around t0. The shape is determined by tau0 and tau1. The function ramps linearly for tau1 time units, stays constant for tau0 time units, and then decays to zero linearly again for tau1 time units.

• tdf_from_file:

%TDFunctions
"function-name" | tdf_from_file | "filename"
%

The temporal shape of the function is contained in a file called filename. This file should contain three columns: first column is time, second and third column are the real part and the imaginary part of the temporal function f(t).

• tdf_from_expr:

%TDFunctions
"function-name" | tdf_from_expr | "expression"
%

The temporal shape of the field is given as an expression (e.g., cos(2.0*t). The letter t means time, obviously. The expression is used to construct the function f that defines the field.

TDGlobalForce
Section: Time-Dependent
Type: string

If this variable is set, a global time-dependent force will be applied to the ions in the x direction during a time-dependent run. This variable defines the base name of the force, that should be defined in the TDFunctions block. This force does not affect the electrons directly.

TDScissor
Section: Time-Dependent
Type: float
Default: 0.0

(experimental) If set, a scissor operator will be applied in the Hamiltonian, shifting the excitation energies by the amount specified. By default, it is not applied.

## Time-Dependent::Absorbing Boundaries

ABCapHeight
Section: Time-Dependent::Absorbing Boundaries
Type: float
Default: -0.2 a.u.

When AbsorbingBoundaries = cap, this is the height of the imaginary potential.

ABShape
Section: Time-Dependent::Absorbing Boundaries
Type: block

Set the shape of the absorbing boundaries. Here you can set the inner and outer bounds by setting the block as follows:

%ABShape
inner | outer | "user-defined"
%

The optional 3rd column is a user-defined expression for the absorbing boundaries. For example, $$r$$ creates a spherical absorbing zone for coordinates with $${\tt inner} < r < {\tt outer}$$, and $$z$$ creates an absorbing plane. Note, values outer larger than the box size may lead in these cases to unexpected reflection behaviours. If no expression is given, the absorbing zone follows the edges of the box (not valid for user-defined box).

ABWidth
Section: Time-Dependent::Absorbing Boundaries
Type: float

Specifies the boundary width. For a finer control over the absorbing boundary shape use ABShape.

AbsorbingBoundaries
Section: Time-Dependent::Absorbing Boundaries
Type: flag
Default: not_absorbing

To improve the quality of the spectra by avoiding the formation of standing density waves, one can make the boundaries of the simulation box absorbing and use exterior complex scaling.
Options:

• not_absorbing: Reflecting boundaries.
• mask: Absorbing boundaries with a mask function.
• cap: Absorbing boundaries with a complex absorbing potential.
• exterior: Exterior complex scaling (not yet implemented).

## Time-Dependent::PhotoElectronSpectrum

Section: Time-Dependent::PhotoElectronSpectrum
Type: float
Default: 1.0

Mask two points enlargement factor. Enlarges the mask box by adding two points at the edges of the box in each direction (x,y,z) at a distance L=Lb*PESMask2PEnlargeFactor where Lb is the box size. This allows to run simulations with an additional void space at a price of adding few points. The Fourier space associated with the new box is restricted by the same factor.

Note: needs PESMaskPlaneWaveProjection = nfft_map or pnfft_map .

Section: Time-Dependent::PhotoElectronSpectrum
Type: float
Default: 1

Mask box enlargement level. Enlarges the mask bounding box by a PESMaskEnlargeFactor. This helps to avoid wavefunction wrapping at the boundaries.

Section: Time-Dependent::PhotoElectronSpectrum
Type: float
Default: -1

In calculation with PESMaskMode = fullmask_mode and NFFT, spurious frequencies may lead to numerical instability of the algorithm. This option gives the possibility to filter out the unwanted components by setting an energy cut-off. If PESMaskFilterCutOff = -1 no filter is applied.

Section: Time-Dependent::PhotoElectronSpectrum
Type: logical
Default: false

Add the contribution of $$\Psi_A$$ in the mask region to the photo-electron spectrum. Literally adds the Fourier components of: $$\Theta(r-R1) \Psi_A(r)$$ with $$\Theta$$ being the Heaviside step function. With this option PES will contain all the contributions starting from the inner radius $$R1$$. Use this option to improve convergence with respect to the box size and total simulation time. Note: Carefully choose $$R1$$ in order to avoid contributions from returning electrons.

Section: Time-Dependent::PhotoElectronSpectrum
Type: integer

PES calculation mode.
Options:

• fullmask_mode: Full mask method. This includes a back action of the momentum-space states on the interaction region. This enables electrons to come back from the continuum.
• passive_mode: Passive analysis of the wf. Simply analyze the plane-wave components of the wavefunctions on the region r > R1. This mode employs a step masking function by default.

Section: Time-Dependent::PhotoElectronSpectrum
Type: integer
Default: fft_map

With the mask method, wavefunctions in the continuum are treated as plane waves. This variable sets how to calculate the plane-wave projection in the buffer region. We perform discrete Fourier transforms (DFT) in order to approximate a continuous Fourier transform. The major drawback of this approach is the built-in periodic boundary condition of DFT. Choosing an appropriate plane-wave projection for a given simulation in addition to PESMaskEnlargeFactor and PESMask2PEnlargeFactorwill help to converge the results.

NOTE: depending on the value of PESMaskMode PESMaskPlaneWaveProjection, may affect not only performance but also the time evolution of the density.
Options:

• fft_out: FFT filtered in order to keep only outgoing waves. 1D only.
• fft_map: FFT transform.
• nfft_map: Non-equispaced FFT map.
• pfft_map: Use PFFT library.
• pnfft_map: Use PNFFT library.

Section: Time-Dependent::PhotoElectronSpectrum
Type: integer
Default: m_sin2

The mask function shape.
Options:

• m_sin2: sin2 mask.

Section: Time-Dependent::PhotoElectronSpectrum
Type: block

Set the size of the mask function. Here you can set the inner (R1) and outer (R2) radius by setting the block as follows:

R1 | R2 | "user-defined"
%

The optional 3rd column is a user-defined expression for the mask function. For example, r creates a spherical mask (which is the default for BoxShape = sphere). Note, values R2 larger than the box size may lead in this case to unexpected reflection behaviours.

Section: Time-Dependent::PhotoElectronSpectrum
Type: float
Default: maxval(mask%Lk)$$^2/2$$

The maximum energy for the PES spectrum.

Section: Time-Dependent::PhotoElectronSpectrum
Type: float

The PES spectrum energy step.

Section: Time-Dependent::PhotoElectronSpectrum
Type: float
Default: -1.0

The time photoelectrons start to be recorded. In pump-probe simulations, this allows getting rid of an unwanted ionization signal coming from the pump. NOTE: This will enforce the mask boundary conditions for all times.

PES_Flux_ARPES_grid
Section: Time-Dependent::PhotoElectronSpectrum
Type: logical

Use a curvilinear momentum space grid that compensates the transformation used to obtain ARPES. With this choice ARPES data is laid out on a Cartesian regular grid. By default true when PES_Flux_Shape = pln.

PES_Flux_AvoidAB
Section: Time-Dependent::PhotoElectronSpectrum
Type: logical
Default: yes

For PES_Flux_Shape = cub, checks whether surface points are inside the absorbing zone and discards them if set to yes (default).

PES_Flux_BZones
Section: Time-Dependent::PhotoElectronSpectrum
Type: block

The block PES_Flux_BZones specifies the number of Brillouin zones along each periodic dimensions. This value will be used to generate the mesh in reciprocal space.

%PES_Flux_EnergyGrid
NBZx | NBZy
%

The option can also be specified with the same value for each direction as PES_Flux_BZones = NBZ. By default PES_Flux_BZones = 2; this means that we will have, on top of the 1st zone, two additional zones at positive and negative reciprocal lattice translation for each periodic dimension.

PES_Flux_DeltaK
Section: Time-Dependent::PhotoElectronSpectrum
Type: float
Default: 0.002

Spacing of the k-mesh in |k| (equidistant).

PES_Flux_EnergyGrid
Section: Time-Dependent::PhotoElectronSpectrum
Type: block

The block PES_Flux_EnergyGrid specifies the energy grid to be used. Only for PES_Flux_Shape = pln.
%PES_Flux_EnergyGrid
Emin | Emax | DeltaE
%

PES_Flux_Gpoint_Upsample
Section: Time-Dependent::PhotoElectronSpectrum
Type: integer

Increase resolution in momentum by adding Gpoints in between adjacent Kpoints. For each additional Gpoint G an entire Kpoint grid centered at G is added to the momentum grid.

PES_Flux_Kmax
Section: Time-Dependent::PhotoElectronSpectrum
Type: float
Default: 1.0

The maximum value of |k|.

PES_Flux_Lmax
Section: Time-Dependent::PhotoElectronSpectrum
Type: integer
Default: 1

Maximum order of the spherical harmonic to be integrated on an equidistant spherical grid (to be changed to Gauss-Legendre quadrature).

PES_Flux_Lsize
Section: Time-Dependent::PhotoElectronSpectrum
Type: block

For PES_Flux_Shape = cub sets the dimensions along each direction. The syntax is:

%PES_Flux_Lsize
xsize | ysize | zsize
%

where xsize, ysize, and zsize are with respect to the origin. The surface can be shifted with PES_Flux_Offset. If PES_Flux_Shape = pln, specifies the position of two planes perpendicular to the non-periodic dimension symmetrically placed at PES_Flux_Lsize distance from the origin.

PES_Flux_Offset
Section: Time-Dependent::PhotoElectronSpectrum
Type: block

Shifts the surface for PES_Flux_Shape = cub. The syntax is:

%PES_Flux_Offset
xshift | yshift | zshift
%

Section: Time-Dependent::PhotoElectronSpectrum
Type: float

The radius of the sphere, if PES_Flux_Shape == sph.

PES_Flux_Shape
Section: Time-Dependent::PhotoElectronSpectrum
Type: integer

The shape of the surface.
Options:

• cub: Uses a parallelepiped as surface. All surface points are grid points. Choose the location of the surface with variable PES_Flux_Lsize (default for 1D and 2D).
• sph: Constructs a sphere with radius PES_Flux_Radius. Points on the sphere are interpolated by trilinear interpolation (default for 3D).
• pln: This option is for periodic systems. Constructs planes perpendicular to the non-periodic dimension symmetrically placed at positive and negative values of PES_Flux_Lsize.

PES_Flux_StepsPhiK
Section: Time-Dependent::PhotoElectronSpectrum
Type: integer
Default: 90

Number of steps in $$\phi$$ ($$0 \le \phi \le 2 \pi$$) for the spherical grid in k.

PES_Flux_StepsPhiR
Section: Time-Dependent::PhotoElectronSpectrum
Type: integer

Number of steps in $$\phi$$ ($$0 \le \phi \le 2 \pi$$) for the spherical surface.

PES_Flux_StepsThetaK
Section: Time-Dependent::PhotoElectronSpectrum
Type: integer
Default: 45

Number of steps in $$\theta$$ ($$0 \le \theta \le \pi$$) for the spherical grid in k.

PES_Flux_StepsThetaR
Section: Time-Dependent::PhotoElectronSpectrum
Type: integer

Number of steps in $$\theta$$ ($$0 \le \theta \le \pi$$) for the spherical surface.

PES_Flux_UseMemory
Section: Time-Dependent::PhotoElectronSpectrum
Type: logical

Use memory to tabulate Volkov plane-wave components on the surface. This option speeds up calculations by precomputing plane-wave phases on the suface. By default true when PES_Flux_Shape = cub.

PES_spm_DeltaOmega
Section: Time-Dependent::PhotoElectronSpectrum
Type: float

The spacing in frequency domain for the photoelectron spectrum (if PES_spm_OmegaMax > 0). The default is PES_spm_OmegaMax/500.

PES_spm_OmegaMax
Section: Time-Dependent::PhotoElectronSpectrum
Type: float
Default: 0.0

If PES_spm_OmegaMax > 0, the photoelectron spectrum is directly calculated during time-propagation, evaluated by the PES_spm method. PES_spm_OmegaMax is then the maximum frequency (approximate kinetic energy) and PES_spm_DeltaOmega the spacing in frequency domain of the spectrum.

Section: Time-Dependent::PhotoElectronSpectrum
Type: float

The radius of the sphere for the spherical grid (if no PES_spm_points are given).

PES_spm_StepsPhiR
Section: Time-Dependent::PhotoElectronSpectrum
Type: integer
Default: 90

Number of steps in $$\phi$$ ($$0 \le \phi \le 2 \pi$$) for the spherical grid (if no PES_spm_points are given).

PES_spm_StepsThetaR
Section: Time-Dependent::PhotoElectronSpectrum
Type: integer
Default: 45

Number of steps in $$\theta$$ ($$0 \le \theta \le \pi$$) for the spherical grid (if no PES_spm_points are given).

PES_spm_points
Section: Time-Dependent::PhotoElectronSpectrum
Type: block

List of points at which to calculate the photoelectron spectrum by the sample point method. If no points are given, a spherical grid is generated automatically. The exact syntax is:

%PES_spm_points
x1 | y1 | z1
%

PES_spm_recipe
Section: Time-Dependent::PhotoElectronSpectrum
Type: integer
Default: phase

The type for calculating the photoelectron spectrum in the sample point method.
Options:

• raw: Calculate the photoelectron spectrum according to A. Pohl, P.-G. Reinhard, and E. Suraud, Phys. Rev. Lett. 84, 5090 (2000).
• phase: Calculate the photoelectron spectrum by including the Volkov phase (approximately), see P. M. Dinh, P. Romaniello, P.-G. Reinhard, and E. Suraud, Phys. Rev. A. 87, 032514 (2013).

PhotoElectronSpectrum
Section: Time-Dependent::PhotoElectronSpectrum
Type: integer
Default: none

This variable controls the method used for the calculation of the photoelectron spectrum. You can specify more than one value by giving them as a sum, for example: PhotoElectronSpectrum = pes_spm + pes_mask
Options:

• none: The photoelectron spectrum is not calculated. This is the default.
• pes_spm: Store the wavefunctions at specific points in order to calculate the photoelectron spectrum at a point far in the box as proposed in A. Pohl, P.-G. Reinhard, and E. Suraud, Phys. Rev. Lett. 84, 5090 (2000).
• pes_mask: Calculate the photo-electron spectrum using the mask method. U. De Giovannini, D. Varsano, M. A. L. Marques, H. Appel, E. K. U. Gross, and A. Rubio, Phys. Rev. A 85, 062515 (2012).
• pes_flux: Calculate the photo-electron spectrum using the t-surff technique, i.e., spectra are computed from the electron flux through a surface close to the absorbing boundaries of the box. (Experimental.) L. Tao and A. Scrinzi, New Journal of Physics 14, 013021 (2012).

## Time-Dependent::Propagation

ArnoldiOrthogonalization
Section: Time-Dependent::Propagation
Type: integer

The orthogonalization method used for the Arnoldi procedure. Only for TDExponentialMethod = lanczos.
Options:

• cgs: Classical Gram-Schmidt (CGS) orthogonalization. The algorithm is defined in Giraud et al., Computers and Mathematics with Applications 50, 1069 (2005).
• drcgs: Classical Gram-Schmidt orthogonalization with double-step reorthogonalization. The algorithm is taken from Giraud et al., Computers and Mathematics with Applications 50, 1069 (2005). According to this reference, this is much more precise than CGS or MGS algorithms.

IonsConstantVelocity
Section: Time-Dependent::Propagation
Type: logical
Default: no

(Experimental) If this variable is set to yes, the ions will move with a constant velocity given by the initial conditions. They will not be affected by any forces.

IonsTimeDependentDisplacements
Section: Time-Dependent::Propagation
Type: block

(Experimental) This variable allows you to specify a time-dependent function describing the displacement of the ions from their equilibrium position: $$r(t) = r_0 + \Delta r(t)$$. Specify the displacements dx(t), dy(t), dz(t) as follows, for some or all of the atoms:

%IonsTimeDependentDisplacements
atom_index | "dx(t)" | "dy(t)" | "dz(t)"
%

The displacement functions are time-dependent functions and should match one of the function names given in the first column of the TDFunctions block. If this block is set, the ions will not be affected by any forces.

MoveIons
Section: Time-Dependent::Propagation
Type: logical

This variable controls whether atoms are moved during a time propagation run. The default is yes when the ion velocity is set explicitly or implicitly, otherwise is no.

RecalculateGSDuringEvolution
Section: Time-Dependent::Propagation
Type: logical
Default: no

In order to calculate some information about the system along the evolution (e.g. projection onto the ground-state KS determinant, projection of the TDKS spin-orbitals onto the ground-state KS spin-orbitals), the ground-state KS orbitals are needed. If the ionic potential changes -- that is, the ions move -- one may want to recalculate the ground state. You may do this by setting this variable.

The recalculation is not done every time step, but only every RestartWriteInterval time steps.

TDDynamics
Section: Time-Dependent::Propagation
Type: integer
Default: ehrenfest

Type of dynamics to follow during a time propagation. For BO, you must set MoveIons = yes.
Options:

• ehrenfest: Ehrenfest dynamics.
• bo: Born-Oppenheimer (Experimental).

TDEnergyUpdateIter
Section: Time-Dependent::Propagation
Type: integer

This variable controls after how many iterations Octopus updates the total energy during a time-propagation run. For iterations where the energy is not updated, the last calculated value is reported. If you set this variable to 1, the energy will be calculated in each step. The default value is 10, unless the ions move, in which case the default is 1.

TDExpOrder
Section: Time-Dependent::Propagation
Type: integer
Default: 4

For TDExponentialMethod = standard or chebyshev, the order to which the exponential is expanded. For the Lanczos approximation, it is the Lanczos-subspace dimension.

TDExponentialMethod
Section: Time-Dependent::Propagation
Type: integer
Default: taylor

Method used to numerically calculate the exponential of the Hamiltonian, a core part of the full algorithm used to approximate the evolution operator, specified through the variable TDPropagator. In the case of using the Magnus method, described below, the action of the exponential of the Magnus operator is also calculated through the algorithm specified by this variable.
Options:

• lanczos: Allows for larger time-steps. However, the larger the time-step, the longer the computational time per time-step. In certain cases, if the time-step is too large, the code will emit a warning whenever it considers that the evolution may not be properly proceeding -- the Lanczos process did not converge. The method consists in a Krylov subspace approximation of the action of the exponential (see M. Hochbruck and C. Lubich, SIAM J. Numer. Anal. 34, 1911 (1997) for details). Two more variables control the performance of the method: the maximum dimension of this subspace (controlled by variable TDExpOrder), and the stopping criterion (controlled by variable TDLanczosTol). The smaller the stopping criterion, the more precisely the exponential is calculated, but also the larger the dimension of the Arnoldi subspace. If the maximum dimension allowed by TDExpOrder is not enough to meet the criterion, the above-mentioned warning is emitted.
• taylor: This method amounts to a straightforward application of the definition of the exponential of an operator, in terms of its Taylor expansion.

$$\exp_{\rm STD} (-i\delta t H) = \sum_{i=0}^{k} {(-i\delta t)^i\over{i!}} H^i.$$

The order k is determined by variable TDExpOrder. Some numerical considerations from Jeff Giansiracusa and George F. Bertsch suggest the 4th order as especially suitable and stable.
• chebyshev: In principle, the Chebyshev expansion of the exponential represents it more accurately than the canonical or standard expansion. As in the latter case, TDExpOrder determines the order of the expansion.

There exists a closed analytic form for the coefficients of the exponential in terms of Chebyshev polynomials:

$$\exp_{\rm CHEB} \left( -i\delta t H \right) = \sum_{k=0}^{\infty} (2-\delta_{k0})(-i)^{k}J_k(\delta t) T_k(H),$$

where $$J_k$$ are the Bessel functions of the first kind, and H has to be previously scaled to $$[-1,1]$$. See H. Tal-Ezer and R. Kosloff, J. Chem. Phys. 81, 3967 (1984); R. Kosloff, Annu. Rev. Phys. Chem. 45, 145 (1994); C. W. Clenshaw, MTAC 9, 118 (1955).

TDIonicTimeScale
Section: Time-Dependent::Propagation
Type: float
Default: 1.0

This variable defines the factor between the timescale of ionic and electronic movement. It allows reasonably fast Born-Oppenheimer molecular-dynamics simulations based on Ehrenfest dynamics. The value of this variable is equivalent to the role of $$\mu$$ in Car-Parrinello. Increasing it linearly accelerates the time step of the ion dynamics, but also increases the deviation of the system from the Born-Oppenheimer surface. The default is 1, which means that both timescales are the same. Note that a value different than 1 implies that the electrons will not follow physical behaviour.

According to our tests, values around 10 are reasonable, but it will depend on your system, mainly on the width of the gap.

Important: The electronic time step will be the value of TDTimeStep divided by this variable, so if you have determined an optimal electronic time step (that we can call dte), it is recommended that you define your time step as:

TDTimeStep = dte * TDIonicTimeScale

so you will always use the optimal electronic time step (more details).

TDLanczosTol
Section: Time-Dependent::Propagation
Type: float
Default: 1e-5

An internal tolerance variable for the Lanczos method. The smaller, the more precisely the exponential is calculated, and also the bigger the dimension of the Krylov subspace needed to perform the algorithm. One should carefully make sure that this value is not too big, or else the evolution will be wrong.

TDMaxSteps
Section: Time-Dependent::Propagation
Type: integer
Default: 1500

Number of time-propagation steps that will be performed. You cannot use this variable together with TDPropagationTime.

TDPropagationTime
Section: Time-Dependent::Propagation
Type: float

The length of the time propagation. You cannot set this variable at the same time as TDMaxSteps. By default this variable will not be used.

The units for this variable are $$\hbar$$/Hartree (or $$\hbar$$/eV if you selected ev_angstrom as input units). The approximate conversions to femtoseconds are 1 fs = 41.34 $$\hbar$$/Hartree = 1.52 $$\hbar$$/eV.

TDPropagator
Section: Time-Dependent::Propagation
Type: integer
Default: etrs

This variable determines which algorithm will be used to approximate the evolution operator $$U(t+\delta t, t)$$. That is, given $$\psi(\tau)$$ and $$H(\tau)$$ for $$\tau \le t$$, calculate $$t+\delta t$$. Note that in general the Hamiltonian is not known at times in the interior of the interval $$[t,t+\delta t]$$. This is due to the self-consistent nature of the time-dependent Kohn-Sham problem: the Hamiltonian at a given time $$\tau$$ is built from the "solution" wavefunctions at that time.

Some methods, however, do require the knowledge of the Hamiltonian at some point of the interval $$[t,t+\delta t]$$. This problem is solved by making use of extrapolation: given a number $$l$$ of time steps previous to time $$t$$, this information is used to build the Hamiltonian at arbitrary times within $$[t,t+\delta t]$$. To be fully precise, one should then proceed self-consistently: the obtained Hamiltonian at time $$t+\delta t$$ may then be used to interpolate the Hamiltonian, and repeat the evolution algorithm with this new information. Whenever iterating the procedure does not change the solution wavefunctions, the cycle is stopped. In practice, in Octopus we perform a second-order extrapolation without a self-consistency check, except for the first two iterations, where obviously the extrapolation is not reliable.

The proliferation of methods is certainly excessive. The reason for it is that the propagation algorithm is currently a topic of active development. We hope that in the future the optimal schemes are clearly identified. In the mean time, if you do not feel like testing, use the default choices and make sure the time step is small enough.
Options:

• qoct_tddft_propagator: WARNING: EXPERIMENTAL
• caetrs: (experimental) Corrected Approximated Enforced Time-Reversal Symmetry (AETRS), this is the previous propagator but including a correction step to the exponential.
• runge_kutta4: WARNING: EXPERIMENTAL. Implicit Gauss-Legendre 4th order Runge-Kutta.
• runge_kutta2: WARNING: EXPERIMENTAL. Implicit 2nd order Runge-Kutta (trapezoidal rule). Similar, but not identical, to Crank-Nicolson method.
• expl_runge_kutta4: WARNING: EXPERIMENTAL. Explicit RK4 method.
• cfmagnus4: WARNING EXPERIMENTAL
• etrs: The idea is to make use of time-reversal symmetry from the beginning:

$$\exp \left(-i\delta t H_{n} / 2 \right)\psi_n = \exp \left(i\delta t H_{n+1} / 2 \right)\psi_{n+1},$$

and then invert to obtain:

$$\psi_{n+1} = \exp \left(-i\delta t H_{n+1} / 2 \right) \exp \left(-i\delta t H_{n} / 2 \right)\psi_{n}.$$

But we need to know $$H_{n+1}$$, which can only be known exactly through the solution $$\psi_{n+1}$$. What we do is to estimate it by performing a single exponential: $$\psi^{*}_{n+1}=\exp \left( -i\delta t H_{n} \right) \psi_n$$, and then $$H_{n+1} = H[\psi^{*}_{n+1}]$$. Thus no extrapolation is performed in this case.
• aetrs: Approximated Enforced Time-Reversal Symmetry (AETRS). A modification of previous method to make it faster. It is based on extrapolation of the time-dependent potentials. It is faster by about 40%. The only difference is the procedure to estimate $$H_{n+1}$$: in this case it is extrapolated via a second-order polynomial by making use of the Hamiltonian at time $$t-2\delta t$$, $$t-\delta t$$ and $$t$$.
• exp_mid: Exponential Midpoint Rule (EM). This is maybe the simplest method, but it is very well grounded theoretically: it is unitary (if the exponential is performed correctly) and preserves time-reversal symmetry (if the self-consistency problem is dealt with correctly). It is defined as: $$U_{\rm EM}(t+\delta t, t) = \exp \left( -i\delta t H_{t+\delta t/2}\right)\,.$$
• crank_nicolson: Classical Crank-Nicolson propagator. $$(1 + i\delta t H_{n+1/2} / 2) \psi_{n+1} = (1 - i\delta t H_{n+1/2} / 2) \psi_{n}$$
• crank_nicolson_sparskit: Classical Crank-Nicolson propagator. Requires the SPARSKIT library. $$(1 + i\delta t H_{n+1/2} / 2) \psi_{n+1} = (1 - i\delta t H_{n+1/2} / 2) \psi_{n}$$
• magnus: Magnus Expansion (M4). This is the most sophisticated approach. It is a fourth-order scheme (a feature which it shares with the ST scheme; the other schemes are in principle second-order). It is tailored for making use of very large time steps, or equivalently, dealing with problem with very high-frequency time-dependence. It is still in a experimental state; we are not yet sure of when it is advantageous.

TDSCFThreshold
Section: Time-Dependent::Propagation
Type: float
Default: 1.0e-6

Since the KS propagator is non-linear, each propagation step should be performed self-consistently. In practice, for most purposes this is not necessary, except perhaps in the first iterations. This variable holds the number of propagation steps for which the propagation is done self-consistently.

The self consistency has to be measured against some accuracy threshold. This variable controls the value of that threshold.

TDStepsWithSelfConsistency
Section: Time-Dependent::Propagation
Type: integer
Default: 0

Since the KS propagator is non-linear, each propagation step should be performed self-consistently. In practice, for most purposes this is not necessary, except perhaps in the first iterations. This variable holds the number of propagation steps for which the propagation is done self-consistently.

The special value all_steps forces self-consistency to be imposed on all propagation steps. A value of 0 means that self-consistency will not be imposed. The default is 0.
Options:

• all_steps: Self-consistency is imposed for all propagation steps.

TDTimeStep
Section: Time-Dependent::Propagation
Type: float

The time-step for the time propagation. For most propagators you want to use the largest value that is possible without the evolution becoming unstable.

The default value is the maximum value that we have found empirically that is stable for the spacing $$h$$: $$dt = 0.0426 - 0.207 h + 0.808 h^2$$ (from parabolic fit to Fig. 4 of http://dx.doi.org/10.1021/ct800518j, probably valid for 3D systems only). However, you might need to adjust this value.

TemperatureFunction
Section: Time-Dependent::Propagation
Type: integer
Default: "temperature"

If a thermostat is used, this variable indicates the name of the function in the TDFunctions block that will be used to control the temperature. The values of the temperature are given in degrees Kelvin.

Thermostat
Section: Time-Dependent::Propagation
Type: integer
Default: none

This variable selects the type of thermostat applied to control the ionic temperature.
Options:

• none: No thermostat is applied. This is the default.
• velocity_scaling: Velocities are scaled to control the temperature.
• nose_hoover: Nose-Hoover thermostat.

ThermostatMass
Section: Time-Dependent::Propagation
Type: float
Default: 1.0

This variable sets the fictitious mass for the Nose-Hoover thermostat.

## Time-Dependent::Response

TDDeltaKickTime
Section: Time-Dependent::Response
Type: float
Default: 0.0

The delta-perturbation that can be applied by making use of the TDDeltaStrength variable, can be applied at the time given by this variable. Usually, this time is zero, since one wants to apply the delta pertubation or "kick" at a system at equilibrium, and no other time-dependent external potential is used. However, one may want to apply a kick on top of a laser field, for example.

TDDeltaStrength
Section: Time-Dependent::Response
Type: float
Default: 0

When no laser is applied, a delta (in time) perturbation with strength TDDeltaStrength can be applied. This is used to calculate, e.g., the linear optical spectra. If the ions are allowed to move, the kick will affect them also. The electric field is $$-(\hbar k / e) \delta(t)$$ for a dipole with zero wavevector, where k = TDDeltaStrength, which causes the wavefunctions instantaneously to acquire a phase $$e^{ikx}$$. The unit is inverse length.

TDDeltaStrengthMode
Section: Time-Dependent::Response
Type: integer
Default: kick_density

When calculating the density response via real-time propagation, one needs to perform an initial kick on the KS system, at time zero. Depending on what kind of response property one wants to obtain, this kick may be done in several modes. For use to calculate triplet excitations, see MJT Oliveira, A Castro, MAL Marques, and A Rubio, J. Nanoscience and Nanotechnology 8, 3392 (2008).
Options:

• kick_density: The total density of the system is perturbed. This mode is appropriate for electric dipole response, as for optical absorption.
• kick_spin: The individual spin densities are perturbed oppositely. Note that this mode is only possible if the run is done in spin-polarized mode, or with spinors. This mode is appropriate for the paramagnetic dipole response, which can couple to triplet excitations.
• kick_spin_and_density: A combination of the two above. Note that this mode is only possible if the run is done in spin-polarized mode, or with spinors. This mode is intended for use with symmetries to obtain both of the responses at once, at described in the reference above.

TDDeltaUserDefined
Section: Time-Dependent::Response
Type: string

By default, the kick function will be a dipole. This will change if (1) the variable TDDeltaUserDefined is present in the inp file, or (2) if the block TDKickFunction is present in the inp file. If both are present in the inp file, the TDKickFunction block will be ignored. The value of TDDeltaUserDefined should be a string describing the function that is going to be used as delta perturbation.

TDKickFunction
Section: Time-Dependent::Response
Type: block

If the block TDKickFunction is present in the input file, and the variable TDDeltaUserDefined is not present in the input file, the kick function to be applied at time zero of the time-propagation will not be a "dipole" function (i.e. $$\phi \rightarrow e^{ikx} \phi$$, but a general multipole in the form $$r^l Y_{lm}(r)$$.

Each line has three columns: integers l and m that defines the multipole, and a weight. Any number of lines may be given, and the kick will be the sum of those multipoles with the given weights.

This feature allows calculation of quadrupole, octupole, etc., response functions.

TDMomentumTransfer
Section: Time-Dependent::Response
Type: block

Momentum-transfer vector for the calculation of the dynamic structure factor. When this variable is set, a non-dipole field is applied, and an output file ftchd is created (it contains the Fourier transform of the charge density at each time). The type of the applied external field can be set by an optional last number. Possible options are qexp (default), qcos, qsin, or qcos+qsin. In the formulae below, $$\vec{q}$$ is the momentum-transfer vector.
Options:

• qexp: External field is $$e^{i \vec{q} \cdot \vec{r}}$$.
• qcos: External field is $$\cos \left( i \vec{q} \cdot \vec{r} \right)$$.
• qsin: External field is $$\sin \left( i \vec{q} \cdot \vec{r} \right)$$.
• qbessel: External field is $$j_l \left( \vec{q} \cdot \vec{r} \right) Y_{lm} \left(\vec{r} \right)$$. In this case, the block has to include two extra values (l and m).

## Time-Dependent::Response::Dipole

TDPolarization
Section: Time-Dependent::Response::Dipole
Type: block

The (real) polarization of the delta electric field. Normally one needs three perpendicular polarization directions to calculate a spectrum (unless symmetry is used). The format of the block is:

%TDPolarization
pol1x | pol1y | pol1z
pol2x | pol2y | pol2z
pol3x | pol3y | pol3z
%

Octopus uses both this block and the variable TDPolarizationDirection to determine the polarization vector for the run. For example, if TDPolarizationDirection=2 the polarization (pol2x, pol2y, pol2z) would be used. These directions may not be in periodic directions.

The default value for TDPolarization is the three Cartesian unit vectors (1,0,0), (0,1,0), and (0,0,1).

Note that the directions do not necessarily need to be perpendicular when symmetries are used.

WARNING: If you want to obtain the cross-section tensor, the TDPolarization block must be exactly the same for the run in each direction. The direction must be selected by the TDPolarizationDirection variable.

TDPolarizationDirection
Section: Time-Dependent::Response::Dipole
Type: integer

When a delta potential is included in a time-dependent run, this variable defines in which direction the field will be applied by selecting one of the lines of TDPolarization. In a typical run (without using symmetry), the TDPolarization block would contain the three Cartesian unit vectors (the default value), and one would make 3 runs varying TDPolarization from 1 to 3. If one is using symmetry, TDPolarization should run only from 1 to TDPolarizationEquivAxes.

TDPolarizationEquivAxes
Section: Time-Dependent::Response::Dipole
Type: integer
Default: 0

Defines how many of the TDPolarization axes are equivalent. This information is stored in a file and then used by oct-propagation_spectrum to rebuild the full polarizability tensor from just the first TDPolarizationEquivAxes directions. This variable is also used by CalculationMode = vdw.

TDPolarizationWprime
Section: Time-Dependent::Response::Dipole
Type: block

This block is needed only when TDPolarizationEquivAxes is set to 3. In such a case, the three directions (pol1, pol2, and pol3) defined in the TDPolarization block should be related by symmetry operations. If A is the symmetry operation that takes you from pol1 to pol2, then TDPolarizationWprime should be set to the direction defined by A$$^{-1}$$pol3. For more information see MJT Oliveira et al., J. Nanoscience and Nanotechnology 8, 3392 (2008).

## Time-Dependent::TD Output

TDExcitedStatesToProject
Section: Time-Dependent::TD Output
Type: block

[WARNING: This is a *very* experimental feature] To be used with TDOutput = populations. The population of the excited states (as defined by where |Phi(t)> is the many-body time-dependent state at time t, and |Phi_I> is the excited state of interest) can be approximated -- it is not clear how well -- by substituting for those real many-body states the time-dependent Kohn-Sham determinant and some modification of the Kohn-Sham ground-state determinant (e.g., a simple HOMO-LUMO substitution, or the Casida ansatz for excited states in linear-response theory. If you set TDOutput to contain populations, you may ask for these approximated populations for a number of excited states, which will be described in the files specified in this block: each line should be the name of a file that contains one excited state.

This file structure is the one written by the casida run mode, in the files in the directory *_excitations. The file describes the "promotions" from occupied to unoccupied levels that change the initial Slater determinant structure specified in ground_state. These promotions are a set of electron-hole pairs. Each line in the file, after an optional header, has four columns:

i a $$\sigma$$ weight

where i should be an occupied state, a an unoccupied one, and $$\sigma$$ the spin of the corresponding orbital. This pair is then associated with a creation-annihilation pair $$a^{\dagger}_{a,\sigma} a_{i,\sigma}$$, so that the excited state will be a linear combination in the form:

$$\left|{\rm ExcitedState}\right> = \sum weight(i,a,\sigma) a^{\dagger}_{a,\sigma} a_{i,\sigma} \left|{\rm GroundState}\right>$$

where weight is the number in the fourth column. These weights should be normalized to one; otherwise the routine will normalize them, and write a warning.

TDFloquetDimension
Section: Time-Dependent::TD Output
Type: integer
Default: -1

Order of Floquet Hamiltonian. If negative number is given, downfolding is performed.

TDFloquetFrequency
Section: Time-Dependent::TD Output
Type: float
Default: 0

Frequency for the Floquet analysis, this should be the carrier frequency or integer multiples of it. Other options will work, but likely be nonsense.

TDFloquetSample
Section: Time-Dependent::TD Output
Type: integer
Default: 20

Number of points on which one Floquet cycle is sampled in the time-integral of the Floquet analysis.

TDMultipoleLmax
Section: Time-Dependent::TD Output
Type: integer
Default: 1

Maximum electric multipole of the density output to the file td.general/multipoles during a time-dependent simulation. Must be non-negative.

TDOutput
Section: Time-Dependent::TD Output
Type: flag
Default: multipoles + energy (+ others depending on other options)

Defines what should be output during the time-dependent simulation. Many of the options can increase the computational cost of the simulation, so only use the ones that you need. In most cases the default value is enough, as it is adapted to the details of the TD run. If the ions are allowed to be moved, additionally the geometry and the temperature are output. If a laser is included it will output by default.

Note: the output files generated by this option are updated every RestartWriteInterval steps.
Options:

• multipoles: Outputs the (electric) multipole moments of the density to the file td.general/multipoles. This is required to, e.g., calculate optical absorption spectra of finite systems. The maximum value of $$l$$ can be set with the variable TDMultipoleLmax.
• angular: Outputs the orbital angular momentum of the system to td.general/angular, which can be used to calculate circular dichroism.
• gauge_field: If set, outputs the vector gauge field corresponding to a spatially uniform (but time-dependent) external electrical potential. This is only useful in a time-dependent periodic run. On by default if GaugeVectorField is set.
• temperature: If set, the ionic temperature at each step is printed. On by default if MoveIons = yes.
• ftchd: Write Fourier transform of the electron density to the file ftchd.X, where X depends on the kick (e.g. with sin-shaped perturbation X=sin). This is needed for calculating the dynamic structure factor. In the case that the kick mode is qbessel, the written quantity is integral over density, multiplied by spherical Bessel function times real spherical harmonic. On by default if TDMomentumTransfer is set.
• dipole_velocity: When set, outputs the dipole velocity, calculated from the Ehrenfest theorem, in the file td.general/velocity. This file can then be processed by the utility oct-harmonic-spectrum in order to obtain the harmonic spectrum.
• eigenvalues: Write the KS eigenvalues.
• ionization_channels: Write the multiple-ionization channels using the KS orbital densities as proposed in C. Ullrich, Journal of Molecular Structure: THEOCHEM 501, 315 (2000).
• total_current: Output the total current (average of the current density over the cell).
• partial_charges: Bader and Hirshfeld partial charges. The output file is called 'td.general/partial_charges'.
• td_kpoint_occup: Project propagated Kohn-Sham states to the states at t=0 given in the directory restart_proj (see %RestartOptions). This is an alternative to the option td_occup, with a formating more suitable for k-points and works only in k- and/or state parallelization
• td_floquet: Compute non-interacting Floquet bandstructure according to further options: TDFloquetFrequency, TDFloquetSample, TDFloquetDimension. This is done only once per td-run at t=0. works only in k- and/or state parallelization
• spin: (Experimental) Outputs the expectation value of the spin, which can be used to calculate magnetic circular dichroism.
• n_excited_el: Output the number of excited electrons, based on the projections of the time evolved wave-functions on the ground-state wave-functions. The output interval of this quantity is controled by the variable TDOutputComputeInterval
• coordinates_sep: Writes geometries in a separate file.
• velocities_sep: Writes velocities in a separate file.
• forces_sep: Writes forces in a separate file.
• total_heat_current: Output the total heat current (average of the heat current density over the cell).
• populations: (Experimental) Outputs the projection of the time-dependent Kohn-Sham Slater determinant onto the ground state (or approximations to the excited states) to the file td.general/populations. Note that the calculation of populations is expensive in memory and computer time, so it should only be used if it is really needed. See TDExcitedStatesToProject.
• geometry: If set (and if the atoms are allowed to move), outputs the coordinates, velocities, and forces of the atoms to the the file td.general/coordinates. On by default if MoveIons = yes.
• dipole_acceleration: When set, outputs the acceleration of the electronic dipole, calculated from the Ehrenfest theorem, in the file td.general/acceleration. This file can then be processed by the utility oct-harmonic-spectrum in order to obtain the harmonic spectrum.
• laser: If set, outputs the laser field to the file td.general/laser. On by default if TDExternalFields is set.
• energy: If set, octopus outputs the different components of the energy to the file td.general/energy. Will be zero except for every TDEnergyUpdateIter iterations.
• td_occup: (Experimental) If set, outputs the projections of the time-dependent Kohn-Sham wavefunctions onto the static (zero-time) wavefunctions to the file td.general/projections.XXX. Only use this option if you really need it, as it might be computationally expensive. See TDProjStateStart. The output interval of this quantity is controled by the variable TDOutputComputeInterval
• local_mag_moments: If set, outputs the local magnetic moments, integrated in sphere centered around each atom. The radius of the sphere can be set with LocalMagneticMomentsSphereRadius.

TDOutputComputeInterval
Section: Time-Dependent::TD Output
Type: integer
Default: 50

The TD output requested are computed when the iteration number is a multiple of the TDOutputComputeInterval variable. Must be >= 0. If it is 0, then no output is written. Implemented only for projections and number of excited electrons for the moment.

TDOutputDFTU
Section: Time-Dependent::TD Output
Type: flag
Default: none

Defines what should be output during the time-dependent simulation, related to LDA+U.

Note: the output files generated by this option are updated every RestartWriteInterval steps.
Options:

• effective_u: Writes the effective U for each orbital set as a function of time.

TDProjStateStart
Section: Time-Dependent::TD Output
Type: integer
Default: 1

To be used with TDOutput = td_occup. Not available if TDOutput = populations. Only output projections to states above TDProjStateStart. Usually one is only interested in particle-hole projections around the HOMO, so there is no need to calculate (and store) the projections of all TD states onto all static states. This sets a lower limit. The upper limit is set by the number of states in the propagation and the number of unoccupied states available.