8.2 Theoretical Background

We briefly state the basic working equations in the following, as far as required to understand the program output. For a more detailed treatment of the theory see refs. [24,101104] and refs. therein. The following discussion is restriced to the one-component (nonrelativistic) treatment, for the sake of convenience.

The first-order frequency dependent response of the density matrix can be expanded as

     ′   ∑           *  ′           *  ′
γ(x,x) =   {Xaiφi(x)φa(x) +Yaiφa(x)φi(x)}.

The (real) expansion coefficients Xai and Y ai are conveniently gathered in a “super-vector”

        (    )
|X,Y ⟩ =   Y

on L, the linear space of products of occupied and virtual ground state MOs φi(x)φa*(x) plus their complex conjugates. X and Y describe the first-order change of the ground state MOs due to an external perturbation which is represented by |P,Qon L. For example, if an oscillating electric dipole perturbation along the z axis is applied, |P,Q= |μz, where μ is the electric dipole operator.

Next we define the 2 × 2 “super-matrices”

   (  A  B  )        ( 1   0 )
Λ =           , Δ =            ,
      B  A             0  - 1

where the four-index quantities A and B are the so-called “orbital rotation Hessians”. Explicit expressions for the standard A and B can be found, e.g., in ref.  [24]. For MGGA functionals, the linear response of the paramagnetic current density leads to additional XC kernel matrix elements, and subsequently to modified defintions of A and B [105]. The vector |X,Y is determined as the solution of the TDHF/TDDFT response problem,

(Λ - ωΔ )|X,Y ⟩ = - |P,Q ⟩.

If |Xα,Y αarises from an electric dipole perturbation |μα, the electronic dipole polarizability at frequency ω is

α αβ(ω ) = - ⟨X α,Yα|μβ⟩,

α,β ∈{x,y,z}. Similarly, if |mαis a component of the magnetic dipole moment operator, the optical rotation is [106]

δαβ(ω) = - c3ωIm⟨Xα,Y α|m β⟩,

where c is the light velocity.

Excitation energies Ωn are the poles of the frequency-dependent density matrix response. They are thus the zeros of the operator on the left-hand side of Eq. (8.4),

(Λ - Ωn Δ)|Xn, Yn⟩ = 0.

The corresponding eigenvectors |Xn,Y nare the transition density matrices for a given excitation (also called “excitation vectors” in the following). They are required to be normalized according to

⟨Xn,Yn|Δ|Xn,Yn ⟩ = 1.

Transition moments are evaluated by taking the trace with one-particle operators, e.g.,

μ0n = ⟨Xn,Yn|μ⟩

for the electric and

m   = ⟨Xn,Yn|m ⟩

for the magnetic transition dipole moments.

The full TDHF/TDDFT formalism is gauge-invariant, i.e., the dipole-length and dipole-velocity gauges lead to the same transition dipole moments in the basis set limit. This can be used as a check for basis set quality in excited state calculations. The TDA can formally be derived as an approximation to full TDHF/TDDFT by constraining the Y vectors to zero. For TDHF, the TDA is equivalent to configuration interaction including all single excitations from the HF reference (CIS). The TDA is not gauge invariant and does not satisfy the usual sum rules [102], but it is somewhat less affected by stability problems (see below). For MGGA functionals, the response of the paramagnetic current density is required to ensure gauge invariance and is included by default.

Stability analysis of closed-shell electronic wavefunctions amounts to computing the lowest eigenvalues of the electric orbital rotation Hessian A + B, which decomposes into a singlet and a triplet part, and of the magnetic orbital rotation Hessian A-B. Note that A-B is diagonal for non-hybrid and non-MGGA DFT, while A + B generally is not. See refs. [23,105,107] for further details.

The two-component relativistic TDDFT eigenvalue problem for excitations of (Kramers-restricted) closed-shell systems taking (approximately) into account the effect of spin-orbit coupling is

M  (X  + Y)n = Ω2n(X + Y)n .

M is a Hermitian (complex) matrix containing spinor energy differences, Coulomb matrix elements as well as matrix elements of the two-component noncollinear exchange-correlation kernel. An explicit expression for M can be found in ref. [108]. (X + Y )n are complex two-component excitation vectors.

Properties of excited states are defined as derivatives of the excited state energy with respect to an external perturbation. It is advantageous to consider a fully variational Lagrangian of the excited state energy [24],

L [X, Y,Ω,C,Z,W ] = EGS + ⟨X, Y|Λ|X, Y⟩ - Ω(⟨X,Y |Δ |X, Y⟩- 1)
                    ∑          ∑
                  +    ZiaFia -   Wpq (Spq - δpq).
                     ia          pq

Here EGS denotes the ground state energy, F and S are the Fock and overlap matrices, respectively, and indices p,q run over all, occupied and virtual MOs.

First, L is made stationary with respect to all its parameters. The additional Lagrange multipliers Z and W enforce that the MOs satisfy the ground state HF/KS equations and are orthonormal. Z is the so-called Z-vector, while W turns out to be the excited state energy-weighted density matrix. Computation of Z and W requires the solution of a single static TDHF/TDKS response equation (8.4), also called coupled and perturbed HF/KS equation. Once the relaxed densities have been computed, excited state properties are obtained by simple contraction with derivative integrals in the atomic orbital (AO) basis. Thus, computation of excited state gradients is more expensive than that of ground state gradients only by a constant factor which is usually in the range of 14.

TDHF/TDDFT expressions for components of the frequency-dependent polarizability ααβ(ω) can also be reformulated as variational polarizability Lagrangians [109]

Lαβ[Xα,Yα,X β,Yβ,C,Zαβ,W αβ](ω )

    = ⟨X α,Yα|(Λ - ωΔ )|X β,Y β⟩+ ⟨X α,Yα|μβ⟩+ ⟨μα|X β,Yβ⟩
   + ∑  Z αβFiaσ - ∑    W αβ(Spqσ - δpq).
      iaσ  iaσ      pqσ,p≤q  pqσ

The stationary point of Lαβ(ω) equals to -ααβ(ω). The requirement that Lαβ(ω) be stationary with respect to all variational parameters determines the Lagrange multipliers Zαβ and Wαβ. All polarizability components αβ are processed simultaneously which allows for computation of polarizability derivatives at the computational cost which is only 2–3 higher than for the electronic polarizability itself.

Within TDDFT and TDHF, the X and Y coefficients are normalized as follows:

  (|Xia|2 - |Yia|2) = 1,

where i and a label occupied and virtual MOs, respectively. Thus, the squared "coefficient" of a single electron excitation from orbital i to orbital a can be defined as

|cia|2 = |Xia |2 - |Yia|2.

escf prints out |cia|2 * 100 starting with the largest coefficient, until the sum of the coefficients is 0.9 or greater. TDA is contained as special case with Y ia=0.