TDSCF: Time-dependent Hartree–Fock and density-functional theory

Code author: Andrew M. James, Daniel G. A. Smith, Ruhee Dcuhna, Roberto Di Remigio and Jeff Schriber

Section author: Roberto Di Remigio

Module: Keywords, PSI Variables, LIBSCF_SOLVER

Introduction

PSI4 provides the capability to calculate excitation energies and ground to excited state transition properties for SCF reference wavefunctions in a linear response formalism [Dreuw2005-wp].

An illustrative example of using the TDSCF functionality is as follows:

molecule {
0 1
O        0.000000    0.695000   -0.092486
O       -0.000000   -0.695000   -0.092486
H       -0.388142    0.895249    0.739888
H        0.388142   -0.895249    0.739888
symmetry c1
}

set {
tdscf_states 10
}

energy('td-scf/cc-pvdz')

This will seek to converge 10 singlet roots from a restricted Hartree–Fock reference. The roots are obtained with an iterative eigensolver and the following is the printout from the calculation:

       ---------------------------------------------------------
                       TDSCF excitation energies
               by Andrew M. James and Daniel G. A. Smith
       ---------------------------------------------------------

==> Options <==

   Residual threshold  : 1.0000e-04
   Initial guess       : denominators
   Reference           : RHF
   Solver type         : RPA (Hamiltonian)


==> Requested Excitations <==

    10 singlet states with A symmetry


==> Seeking the lowest 10 singlet states with A symmetry

                      Generalized Hamiltonian Solver
                            By Andrew M. James

==> Options <==

  Max number of iterations        = 60
  Eigenvector tolerance           = 1.0000e-04
  Max number of expansion vectors = 2000

=> Iterations <=
                            Max[D[value]]     Max[|R|]   # vectors
HamiltonianSolver iter   1:   5.64572e-01  3.65441e-01     40
HamiltonianSolver iter   2:   1.70649e-02  4.40807e-02     60
HamiltonianSolver iter   3:   2.42552e-04  6.95387e-03     80
HamiltonianSolver iter   4:   2.34146e-06  7.75689e-04    100
HamiltonianSolver iter   5:   1.75483e-08  6.17293e-05    120      Converged

When convergence is reached, PSI4 will output a report of excitation energies, oscillator strengths, and rotatory strenghts in atomic units:

                                Excitation Energy         Total Energy        Oscillator Strength             Rotatory Strength
 #   Sym: GS->ES (Trans)        au              eV              au          au (length)    au (velocity)    au (length)    au (velocity)
---- -------------------- --------------- --------------- --------------- --------------- --------------- --------------- ---------------
 1        A->A (1 A)       0.26945         7.33199        -150.50964       0.0017          0.0082         -0.0019         -0.0135
 2        A->A (1 A)       0.31534         8.58073        -150.46375       0.0000          0.0002         -0.0007         -0.0096
 3        A->A (1 A)       0.35760         9.73076        -150.42148       0.0040          0.0097          0.0227          0.0352
 4        A->A (1 A)       0.37522         10.21028       -150.40386       0.0144          0.0442          0.0729          0.1223
 5        A->A (1 A)       0.43252         11.76960       -150.34656       0.0890          0.1189         -0.1942         -0.2491
 6        A->A (1 A)       0.46952         12.77624       -150.30957       0.0640          0.1157          0.0175          0.0235
 7        A->A (1 A)       0.49186         13.38426       -150.28722       0.0016          0.0012         -0.0243         -0.0212
 8        A->A (1 A)       0.50405         13.71581       -150.27504       0.4557          0.4396         -0.0197         -0.0158
 9        A->A (1 A)       0.52971         14.41407       -150.24938       0.0799          0.0948          0.0546          0.0595
 10       A->A (1 A)       0.56083         15.26092       -150.21825       0.0497          0.0567         -0.0587         -0.0650

The solvers can be used to extract the first few roots of interest for the full time-dependent DFT (TDDFT) equations, also known as the random-phase approximation (RPA), or its Tamm–Dancoff approximation. The former is a generalized eigenvalue problem and our solver leverages the Hamiltonian structure of the equations to ensure robust convergence [stratmann:1998]. The latter is a Hermitian eigenvalue problem and we employ a Davidson solver.

Known limitations

Warning

The implementation cannot currently handle the following cases: - Excited states of triplet symmetry from a restricted DFT reference. - Excited states from an unrestricted reference other than HF or LDA.

Warning

The length-gauge rotatory strengths PSI4 computes are currently not gauge-origin invariant.

Theory

The excitation energies and corresponding states are obtained from the following generalized eigenvalue problem, also known as the response eigenvalue problem:

$$\begin{array}{r}\left(\begin{array}{cc}\mathbf{A}& \mathbf{B}\\ {\mathbf{B}}^{\ast }& {\mathbf{A}}^{\ast }\end{array}\right)\left(\begin{array}{c}{\mathbf{X}}_n\\ {\mathbf{Y}}_n\end{array}\right)={\omega }_n\left(\begin{array}{cc}\mathbf{1}& \mathbf{0}\\ \mathbf{0}& -\mathbf{1}\end{array}\right)\left(\begin{array}{c}{\mathbf{X}}_n\\ {\mathbf{Y}}_n\end{array}\right).\end{array}$$

This approach has the advantage that there is no need to explicitly parametrize the wavefunctions of the molecular excited states. Furthermore, the excitation eigenvectors, $({\mathbf{X}}_n{\mathbf{Y}}_n{)}^t$, provide information on the nature of the transitions and can be used to form spectroscopic observables, such as oscillator and rotatory strengths.

The $\mathbf{A}$ and $\mathbf{B}$ matrices appearing on the left-hand side are the blocks of the molecular electronic Hessian, [Norman2018-tn] whose dimensionality is $(OV{)}^2$, with $O$ and $V$ the number of occupied and virtual molecular orbitals, respectively. This prevents explicit formation of the full Hessian, and subspace iteration methods need to be used to extract the first few roots. In such methods, the eigenvectors are expanded in a subspace of trial vectors, whose dimensionality is greatly lower than that of the full eigenproblem. The Hessian is projected down to this subspace where conventional full diagonalization algorithms can be applied. The subspace is augmented with new trial vectors, until a suitable convergence criterion is met. The efficiency of the subspace solver is determined by the first half-projection of the Hessian in the trial subspace, that is, by the efficiency of the routines performing the matrix-vector products.

It is essential to note that, despite the hermiticity of the molecular electronic Hessian, the response eigenvalue equation is not an Hermitian eigenproblem, due to the nonunit metric on the right-hand side. Indeed the Davidson solver, the standard subspace iteration method in quantum chemistry, demonstrates very poor convergence, sometimes manifesting as spurious complex eigenvalues. The eigenproblem however has Hamiltonian symmetry: the roots appear in pairs $({\omega }_n,-{\omega }_n)$, as do the eigenvectors. A robust subspace solver should preserve the Hamiltonian symmetry, by enforcing the paired structure on the trial vectors themselves. Since PSI4 employs real orbitals, the response eigenproblem can be brought to the form:

$$(\mathbf{A}-\mathbf{B})(\mathbf{A}+\mathbf{B})|{\mathbf{X}}_n+{\mathbf{Y}}_n⟩={\omega }_n^2|{\mathbf{X}}_n+{\mathbf{Y}}_n⟩,$$

and further to the Hermitian form:

$$(\mathbf{A}-\mathbf{B}{)}^{\frac{1}{2}}(\mathbf{A}+\mathbf{B})(\mathbf{A}-\mathbf{B}{)}^{\frac{1}{2}}{\mathbf{T}}_n={\omega }_n^2{\mathbf{T}}_n,$$

assuming the SCF reference is stable, i.e. $(\mathbf{A}-\mathbf{B})$ is positive-definite. The paired vectors $|{\mathbf{X}}_n-{\mathbf{Y}}_n⟩$ are left eigenvectors and form a biorthonormal set together with the right eigenvectors $|{\mathbf{X}}_n+{\mathbf{Y}}_n⟩$.

The algorithm for the subspace iteration Hamiltonian solver implemented in PSI4 was first described by Stratmann et al. [stratmann:1998]. As already mentioned, the formation and storage of the matrix-vector products $(\mathbf{A}+\mathbf{B}){\mathbf{b}}_i$ and $(\mathbf{A}-\mathbf{B}){\mathbf{b}}_i$ for all trial vectors ${\mathbf{b}}_i$ are the most compute- and memory-intensive operations in the Hamiltonian solver. These matrix-vector products are equivalent to building generalized Fock matrices and thus use the efficient $JK$ build infrastructure of PSI4.

The excitation energies and eigenvectors can then be used to compute transition moments, such as electric and magnetic transition dipole moments, and spectroscopic intensities, such as oscillator strengths and rotatory strengths [Pedersen1995-du], [Lestrange2015-xn]. For example, PSI4 will compute compute oscillator strengths from the MO basis electric dipole moment integrals, ${\mu }_u$, and the right excitation vectors, $|{\mathbf{X}}_n+{\mathbf{Y}}_n⟩$:

$$f=\frac{2}{3}{\omega }_n\sum _{u=x,y,z}\sum _{ia}|({\mathbf{X}}_n+{\mathbf{Y}}_n{)}_{ia}{\mu }_{ai,u}{|}^2.$$

Psithon keywords

TDSCF_STATES

Number of roots (excited states) we should seek to converge. This can be either an integer (total number of states to seek) or a list (number of states per irrep). The latter is only valid if the system has symmetry. Furthermore, the total number of states will be redistributed among irreps when symmetry is used.

• Type: array

• Default: No Default

TDSCF_TRIPLETS

Controls inclusion of triplet states, which is only valid for restricted references. Valid options: - none : No triplets computed (default) - also : lowest-energy triplets and singlets included, in 50-50 ratio. Note that singlets are privileged, i.e. if seeking to converge 5 states in total, 3 will be singlets and 2 will be triplets. - only : Only triplet states computed

• Type: string

• Possible Values: NONE, ALSO, ONLY

• Default: NONE

TDSCF_TDA

Run with Tamm-Dancoff approximation (TDA), uses random-phase approximation (RPA) when false

• Type: boolean

• Default: false

TDSCF_R_CONVERGENCE

Convergence threshold for the norm of the residual vector. If unset, default based on D_CONVERGENCE

• Type: conv double

• Default: 1e-4

TDSCF_MAXITER

Maximum number of TDSCF solver iterations

• Type: integer

• Default: 60

TDSCF_GUESS

Guess type, only ‘denominators’ currently supported

• Type: string

• Default: DENOMINATORS

TDSCF_PRINT

Verbosity level in TDSCF

• Type: integer

• Default: 1

PsiAPI usage

The TDSCF functionality is also accessible from PsiAPI. The example calculation shown above can be carried out as follows:

import psi4

from psi4.driver.procrouting.response.scf_response import tdscf_excitations

psi4.core.set_output_file("h2o2.out")

h2o2 = psi4.geometry("""0 1
O        0.000000    0.695000   -0.092486
O       -0.000000   -0.695000   -0.092486
H       -0.388142    0.895249    0.739888
H        0.388142   -0.895249    0.739888
symmetry c1
""", name="H2O2")

psi4.set_options({
    'save_jk': True,
})

e, wfn = psi4.energy("HF/cc-pvdz", return_wfn=True, molecule=h2o2)
res = tdscf_excitations(wfn, states=10)

Plotting one-photon absorption and electronic circular dichroism spectra

Excitation energies and corresponding spectroscopic observables can be used to produce spectra for one-photon absorption (OPA) and electronic circular dichroism (ECD) with phenomenological line broadening.

PSI4 provides the spectrum function for this purpose implementing the recommendations of Rizzo et al. [Rizzo2011-to]. This function will not plot the spectrum, but rather return a pair of NumPy arrays containing the $x$ and $y$ values resulting from the convolution with broadening of the computed spectroscopic observables.

import numpy as np

import psi4

from psi4.driver.procrouting.response.scf_response import tdscf_excitations
from psi4.driver.p4util import spectrum

psi4.core.set_output_file("moxy.out")

moxy = psi4.geometry("""0 1
C  0.152133 -0.035800  0.485797
C -1.039475  0.615938 -0.061249
C  1.507144  0.097806 -0.148460
O -0.828215 -0.788248 -0.239431
H  0.153725 -0.249258  1.552136
H -1.863178  0.881921  0.593333
H -0.949807  1.214210 -0.962771
H  2.076806 -0.826189 -0.036671
H  2.074465  0.901788  0.325106
H  1.414895  0.315852 -1.212218
""", name="(S)-methyloxirane")

psi4.set_options({
    'save_jk': True,
})

e, wfn = psi4.energy("HF/cc-pvdz", return_wfn=True, molecule=moxy)
res = tdscf_excitations(wfn, states=8, triplets="also")

# get poles and residues to plot OPA and ECD spectra
poles = [r["EXCITATION ENERGY"] for r in res]
opa_residues = [np.linalg.norm(r["LENGTH-GAUGE ELECTRIC DIPOLE TRANSITION MOMENT"])**2 for r in res]
ecd_residues = [r["LENGTH-GAUGE ROTATORY STRENGTH"] for r in res]

opa_spectrum = spectrum(poles=poles, residues=opa_residues, gamma=0.01, out_units="nm")
ecd_spectrum = spectrum(poles=poles, residues=ecd_residues, kind="ECD", gamma=0.01, out_units="nm")

The data produced by running the above PsiAPI code can, for example, be used with the Altair plotting library to produce the desired spectra.