7.14 Real-Time SCF Methods (RT-TDDFT, RT-HF, OSCF2)

Linear response calculations are the most efficient way to predict resonant electronic response frequencies and intensities. Q-Chem can also propagate mean-field theories (HF, DFT) in real-time with and without dissipation and finite-temperature effects. These methods time-evolve the Liouville-von Neumann equation of motion for a one-particle density operator, $\widehat{\rho }(t)$:

$$\begin{array}{c}\hfill (i/\mathrm{\hslash })\widehat{\rho }(t)=[\widehat{F}\{\rho \},\widehat{\rho }(t)]\end{array}$$

(7.98)

These real-time methods are useful for simulating attosecond to picosecond timescales, density response to multiple or strong applied fields, and timescales of energy relaxation. Symmetry, non-singlet densities, Meta-GGAs and 5th-rung functionals are not supported for real-time propagation, and nuclei are fixed. The CPU cost of a single time step is comparable to that of a single SCF cycle, if the modified-midpoint unitary transformation (MMUT) propagation algorithm⁵⁸⁰ is employed, and no more than a few times the cost of a single SCF cycle if predictor/corrector algorithms (which facilitate stable time propagation using much larger time steps) are used.¹¹¹⁹ Memory costs are similar to those of a ground-state SCF calculation, which for large systems, or those with a large density of states, represents a dramatic reduction in the memory requirement relative to linear-reponse TDDFT. The number of required Fock-build steps can be estimated from the default electronic time-step, which is 0.02 atomic units ($4.8\times {10}^{-4}$ fs). The real-time code exploits real-time shared-memory parallelism, and the use of 8 or more cores (–nt 8) is suggested. The input file of a basic real-time propagation is relatively simple, but sophisticated jobs require additional input files, and generate additional output.

$molecule
   0 1
   B    0.000000     0.000000     0.000000
   H    1.115609    -0.322048     0.260309
   H   -0.332368     1.137184     0.111382
   H   -0.782205    -0.814474    -0.375294
$end

$rem
   SCF_GUESS         core
   GEN_SCFMAN        true
   RTTDDFT           1
   BASIS             3-21g
   THRESH            10
   SCF_CONVERGENCE   9
   EXCHANGE          b3lyp
   SYMMETRY          false
   SYM_IGNORE        TRUE
   SCF_ALGORITHM     DIIS
   UNRESTRICTED      FALSE
   MAX_SCF_CYCLES    200
$end

Before executing this job, a directory (/logs) must be made in the output directory of the job to collect most of the results of the propagation. By default, a real-time job begins from the ground state SCF density, applies a weak, brief (0.07 atomic time units) oscillating electric field to the $y$-axis of a molecule which excites a superposition of all $y$-polarized excited states, and outputs the resulting time-dependent dipole moment, energy and other quantities to the logs directory. The parameters of the time-dependent run are printed when the propagation begins. During the propagation the state of the molecule and propagation is summarized in the Q-Chem output file, including estimates of the elapsed and remaining simulation wall-time. Propagation stops when MaxIter time steps are exceeded, or will stop prematurely if the density matrix becomes nonphysical. The logs directory will be populated by several text, white space delimited tables by default. The ./logs/Pol.csv file is a table consisting of time (in a.u.), ${\mu }_x$, ${\mu }_y$, ${\mu }_z$, total energy, $\text{tr}\widehat{\rho }$, the gap, the electronic energy, and two unused columns. The ./logs/Fourier_x.csv file contains Fourier transform $ℱ[{\mu }_x(t)]$, in the format of its total value in energy units (eV) followed by its negative imaginary part and then its real part. Analogous files ./logs/Fourier_y.csv and ./logs/Fourier_z.csv are also created.

By adjusting the options of the propagation with a file TDSCF.prm, a much larger amount of output can be generated, including the electron and hole densities at all times as sequential MolDen files (*.mol) viewable with the free package Gabedit (gabedit.sourceforge.net). The *.mol files generated by the real-time code have fractional orbital occupation numbers, and do not render properly in viewers other than Gabedit to our knowledge. So long as the applied field is weak and has short duration the positions at which peaks appear in the Fourier are the same as a linear-response TDDFT-RPA job (note: not the same as a LR-TDDFT-TDA job), as shown in the sample. In an energy-conserving job, the width of the peak is inversely proportional to the duration of the signal sample used to construct the Fourier transform. If the applied field pulse is long ($>1$ a.u.), or has strong intensity ($>0.01$ a.u.) non-linear response can be studied. Non-linear effects in real-time SCF are an area of active investigation and development.

The absorption cross-section ${\sigma }_{ii}(\omega )$ in the $i$ direction, $i\in \{x,y,z\}$, can be obtained from the imaginary part of the frequency-dependent polarizability ${\alpha }_{ii}(\omega )$:

$${\sigma }_{ii}(\omega )=\left(\frac{4\pi \omega }{c}\right)\mathrm{\Im }\left[{\alpha }_{ii}(\omega )\right],$$

(7.99)

and the rotationally-averaged absorption spectrum (within the dipole approximation) is

$$A(\omega )=\frac{1}{3}\left[{\sigma }_{xx}(\omega )+{\sigma }_{yy}(\omega )+{\sigma }_{zz}(\omega )\right].$$

(7.100)

The frequency-dependent polarizability ${\alpha }_{ij}(\omega )$ is obtained from the Fourier transform of ${\mu }_i(t)$, for a perturbing field in the $j$ direction:¹¹¹⁹

$${\alpha }_{ij}(\omega )=\frac{ℱ\left[{\mu }_i(t)\right]}{ℱ\left[E_j(t)\right]}.$$

(7.101)

Thus to compute the spectrum in Eq. (7.100), three RT-TDDFT simulations are required, perturbing separately in the $x$, $y$, and $z$ directions. A script to obtain the spectrum after the RT-TDDFT calculation is provided and can be find at bin/sp3d.m. A script to obtain the spectrum after the RT-TDDFT calculation is provided and can be find at bin/sp3d.m.

Because of the large number of floating point arguments used to control a real-time job, a separate input file TDSCF.prm in the same directory as the Q-Chem input file is used for parameters. The file is two columns of plain text. The first column is a string naming the parameter (which must match the case and spelling printed in the output exactly to be interpreted correctly), and the second column is a floating point number. The number must only be followed by a new line. Several inputs are interpreted as true = 1 or false = 0. The most useful parameters are discussed in this section. The code will signal if this file is not found and default values will be supplied instead.

Q-Chem Typical TDSCF.prm for a real-time B3LYP calculation.

dt               0.02
Stabilize        0
TCLOn            0
MaxIter          80000
ApplyImpulse     1
ApplyCw          0
FieldFreq        0.7
Tau              0.07
FieldAmplitude   0.001
ExDir            1.0
EyDir            1.0
EzDir            1.0
Print            0
StatusEvery      10
SaveDipoles      1
DipolesEvery     2
SavePopulations  0
SaveFockEnergies 0
WriteDensities   0
SaveEvery        500
FourierEvery     5000
MMUT             1
LFLPPC           0