63 THE TDHF AND TDKS PROGRAMS

Real-time electronic dynamics using time-dependent Hartree-Fock and time-dependent Kohn-Sham theories can be performed using the commands TDHF and TDKS respectively, which have to be preceded by a HF and KS command. Unrestriced versions are available through TDUHF and TDUKS and should be preceded by UHF and UKS runs respectively. All methods require symmetry to be switched off. For details on the theory and methods see H. Eshuis, G. G. Balint-Kurti and F. R. Manby, J. Chem. Phys. $\mathbf{128}$, 114113 (2008), and references therein. The commands take several options:

command,t=,dt=,ns=,ng=,grsize=,print=;
PULSE,options

The total propagation time (in au) is set by t; dt sets the timestep and ns the number of steps, where two of the three have to be provided. ng sets the number of grid points in one dimension (default = 0) and grsize the grid size in bohr (default = 10 bohr). Setting ng $>2$ switches on the calculation of quantum currents (see below). The option print determines the level of output (0=normal output, 1=object linear in matrices, 2=matrices as well, $>2$ debug). The subcommand PULSE determines the envelope used, and takes several options depending on the envelope selected. Possibilities are:

NONE: no pulse, no options

STEP, $e_x, e_y, e_z$, length
a DC-field of strength $e_q$ is applied in the $q$th direction for a total time length (in au)

DC, $e_x, e_y, e_z, \alpha$
a DC-field of strength $e_q$ is applied in the $q$th direction. The field is switched on exponentially with a rate determined by $\alpha$.

TRAP, $e_x, e_y, e_z, \omega, \alpha$
an oscillating field of strength $e_q$ is applied in the $q$th direction. The field oscillates with angular frequency $\omega$. The envelope reaches $e_q$ in one period of the field, stays constant for $\alpha$ periods and then decays to zero in one period.

CW, $e_x, e_y, e_z, \omega$
an oscillating field of strength $e_q$ is applied in the $q$th direction. The field oscillates with angular frequency $\omega$. No envelope present.

CWSIN, $e_x, e_y, e_z, \omega, \alpha$
an oscillating field of strength $e_q$ is applied in the $q$th direction. The field oscillates with angular frequency $\omega$ and is switched on using a $\sin^2$ envelope, where $\alpha$ determines how fast the field is switched on.

CWGAUSS, $e_x, e_y, e_z, \omega, \alpha$
like CWSIN, but with a Gaussian envelope

The finite pulses TRAP and STEP produce an absorption spectrum obtained from the time-dependent dipole moment sampled after the field is switched off. It is located in ($input).spec, or in molpro.spec when running interactively, and contains of 4 columns, which contain energy (eV) and the absorption in the $x, y, z$ direction respectively.

All runs produce a file ($input).dat (or molpro.dat) which contains time-dependent properties, like the components of the field, components of the dipole moment, total energy, orbital occupation numbers, orbital energies and total number of electrons. The Molpro output file specifies the order of the data in the file. For very long runs the size of the file is restricted by printing only every twentieth set of data.

In case of an unrestricted run three files are produced: ($input).dat, ($input)a.dat, ($input)b.dat. The first file contains the generic data about the field, total energy, dipole moments and the expectation value of the total spin operator. The other two files contain the orbital energies and occupation numbers for the $\alpha$ and $\beta$ electrons respectively.

TDKS and TDUKS use the options provided by the KS and UKS comands. At the moment TDKS/TDUKS suffers from numerical instabilities when using strong fields. Divergence of the energy is observed, possibly due to the use of quadrature for the evaluation of the potential.

Quantum currents will be calculated when choosing $ng>2$. A cubic grid will be computed of size grsize with a total of $ng^3$ gridpoints. The imaginary part of the density will be summed at every timestep and after the dynamics the total current will be evaluated at every gridpoint. The user can extract the required data from this array by printing out parts of it, or by integrating over point, but this requires actual coding, as this has not been implemented sufficiently neat. It is also straightforward to evaluate the currents at every timestep or at selected timesteps. This is not done automatically, because it slows down the dynamics considerably.