Time-dependent Hartree-Fock and density functional theory

Modules: pyscf.tdscf, pyscf.pbc.tdscf

Introduction

PySCF implements the time-dependent Hartree-Fock (TDHF) and time-dependent density functional theory (TDDFT) (frequency domain) linear response theories to compute excited-state energies and transition properties in the tdscf module. A minimal example that runs a TDDFT calculation is as follows

from pyscf import gto, scf, dft, tddft
mol = gto.Mole()
mol.build(
    atom = 'H 0 0 0; F 0 0 1.1',  # in Angstrom
    basis = '631g',
    symmetry = True,
    verbose = 4,
)

mf = dft.RKS(mol)
mf.xc = 'b3lyp'
mf.kernel()

mytd = tddft.TDDFT(mf)
mytd.nstates = 10
mytd.kernel()
mytd.analyze()

The example above computes the excitation energies, oscillator strengths and transition dipole moments of the ten lowest singlet exicted states.

Theory

Using first-order time-dependent perturbation theory within HF or KS theory, one obtains the non-Hermitian TDHF or TDDFT equations for the excitation energies [51]:

\[\begin{split}\left(\!\!\begin{array}{ll} \mathbf{A} & \mathbf{B} \\ \mathbf{B}^\ast & \mathbf{A}^\ast \end{array}\!\!\right) \left(\!\!\begin{array}{c} \mathbf{X} \\ \mathbf{Y} \end{array}\!\!\right) = \omega \left(\!\!\begin{array}{lr} \mathbf{1} & \mathbf{0} \\ \mathbf{0} & -\mathbf{1} \end{array}\!\!\right) \left(\!\!\begin{array}{c} \mathbf{X} \\ \mathbf{Y} \end{array}\!\!\right) \;,\end{split}\]

where \(\mathbf{A}\) and \(\mathbf{B}\) are the orbital hessians which also appear in the stability analysis for reference states (see Stability analysis), \(\omega\) is the excitation energy, and \(\mathbf{X}\) and \(\mathbf{Y}\) represent the response of the density matrix. In cases where the system possesses a degenerate ground state or has triplet instabilities, the algorithms used to solve the above equations may be unstable. This can be solved by applying the Tamm-Dancoff approximation (TDA) [52], which simply neglects the \(\mathbf{B}\) and \(\mathbf{Y}\) matrices and leads to a Hermitian eigenvalue problem

\[\mathbf{AX} = \omega\mathbf{X} \;.\]

Methods

For TDHF or TDDFT calculations, the reference state can be either restricted or unrestricted:

mytd = mol.RKS().run().TDDFT().run()
mytd = mol.UKS().run().TDDFT().run()

By default, only singlet excited states are computed. In order to compute triplet excited states, one needs to set the attribute singlet to False:

mytd.singlet = False
mytd.kernel()

One can also perform symmetry analysis by calling the analyze() method, which also computes the oscillator strengths and dipole moments:

mytd.analyze(verbose=4)

Property calculation

Oscillator strengths

Oscillator strengths for each excited state can be computed in both length and velocity gauges:

mytd.oscillator_strength(gauge='length')
mytd.oscillator_strength(gauge='velocity')

Higher order corrections [53] to the oscillator strength can also be included:

#include corrections due to magnetic dipole and electric quadruple
mytd.oscillator_strength(gauge='velocity', order=1)
#also include corrections due to magnetic quadruple and electric octupole
mytd.oscillator_strength(gauge='velocity', order=2)

Transition moments

PySCF implements various types of transition moments between the reference SCF state and the TDHF or TDDFT excited states. These include:

  • electric dipole, quadrupole and octupole transition moments in both length and velocity gauges:

    mytd.transition_dipole()
    mytd.transition_velocity_dipole()
    mytd.transition_quadrupole()
    mytd.transition_velocity_quadrupole()
    mytd.transition_octupole()
    mytd.transition_velocity_octupole()
    
  • magnetic dipole and quadrupole transition moments:

    mytd.transition_magnetic_dipole()
    mytd.transition_magnetic_quadrupole()
    

Nuclear gradients

Analytic nuclear gradients are available for TDHF and TDDFT, and they can be computed as follows:

tdg = mytd.Gradients()
g1 = tdg.kernel() #default will compute the gradients of first excited state
g1 = tdg.kernel(state=1) #first excited state
g2 = tdg.kernel(state=2) #second excited state

Natural transition orbital analysis

Natural transition orbitals (NTOs) can be computed by singular value decomposition of the transition density matrix. In PySCF, these orbitals can be obtained as follows:

weights, nto_coeff = mytd.get_nto(state=1)

where nto_coeff are the coefficients for NTOs represented in AO basis, and they are ordered as occupied orbitals followed by virtual orbitals.