TDSCF with Psi4¶

Roberto Di Remigio (@robertodr)

European National Competence Centre Sweden

Find this notebook at https://tinyurl.com/tddft-psicon2020¶

Psi4 1.4 enables you to run TDSCF calculations:

1. using the random-phase appoximation (RPA). Credits: Andrew James (@amjames) and Daniel Smith (@dgasmith).
2. using the Tamm-Dancoff approximation (TDA). Credits: Ruhee D'Cunha (@rdcunha).
3. from Psithon input files. Credits: Jeff Schriber (@jeffschriber) and Lori Burns (@loriab).
4. writing Python scripts and PsiAPI. Credits: Andrew James (@amjames) and yours truly.
5. including solvent effects with either the polarizable continuum model (PCM) or polarizable embedding (PE). Credits: Holger Kruse (@hokru) and Maximilian Scheurer (@maxscheurer).
6. and easily plot the corresponding spectra. Credits: yours truly.

Shift + Enter

Using Psithon¶

The bridge to Psithon was built by Jeff Schriber (@jeffschriber) and Lori Burns (@loriab). A TDSCF calculation can be invoked with two styles:

• From the energy driver

energy("td-pbe0/cc-pvdz")

• With a wavefunction object

e, wfn = energy("pbe0/cc-pvdz", return_wfn=True)
tdscf(wfn)

We will ask for 10 roots (states) and to use the Tamm-Dancoff approximation (TDA):

set {
  tdscf_states 10
  tdscf_tda true
}

other options are detailed in the manual. The number of roots can be specified either as an integer or as a list. The latter form is mostly useful when dealing with symmetry: you can ask to solve for roots belonging to a specific irrep. This is the complete input file:

molecule moxy {
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
}

set_num_threads(4)

set {
  tdscf_states 10
  tdscf_tda true
}

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

Now we can run with:

$ psi4 -i moxy.inp -o moxy.dat

The result sections reports excitation energies and transition properties:

******************************************************************************************
**********                               WARNING                                **********
**********  Length-gauge rotatory strengths are **NOT** gauge-origin invariant  **********
******************************************************************************************

                                    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.30939         8.41882        -192.58499       0.0266          0.0185          0.0824          0.0552        
     2        A->A (1 A)       0.31237         8.50013        -192.58200       0.0052          0.0005         -0.0058         -0.0025        
     3        A->A (1 A)       0.33179         9.02838        -192.56259       0.0171          0.0077         -0.0356         -0.0294        
     4        A->A (1 A)       0.33881         9.21956        -192.55556       0.0474          0.0582         -0.0230         -0.0325        
     5        A->A (1 A)       0.34136         9.28892        -192.55301       0.0169          0.0179          0.0108          0.0181        
     6        A->A (1 A)       0.34444         9.37258        -192.54994       0.0130          0.0059         -0.0520         -0.0433        
     7        A->A (1 A)       0.36218         9.85546        -192.53219       0.0231          0.0170         -0.0006          0.0053        
     8        A->A (1 A)       0.37258         10.13855       -192.52179       0.0307          0.0252          0.0725          0.0708        
     9        A->A (1 A)       0.37828         10.29347       -192.51610       0.0041          0.0039          0.0001          0.0040        
     10       A->A (1 A)       0.37904         10.31409       -192.51534       0.1130          0.0521          0.0460          0.0310

Further information on the set up of the calculation is reported a bit earlier in the output, together with the progress of the iterative subspace solver:

---------------------------------------------------------
                         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         : TDA (Davidson)


  ==> Requested Excitations <==

      10 singlet states with A symmetry


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

                         Generalized Davidson Solver                         
                               By Ruhee Dcunha                               

  ==> Options <==

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

  => Iterations <=
                           Max[D[value]]     Max[|R|]   # vectors
  DavidsonSolver iter   1:   3.84154e-01  5.04370e-02     40      
  DavidsonSolver iter   2:   2.40990e-03  3.49018e-02     50      
  DavidsonSolver iter   3:   2.59452e-03  9.67573e-03     60      
  DavidsonSolver iter   4:   1.11697e-04  1.30981e-03     70      
  DavidsonSolver iter   5:   2.39422e-06  1.65681e-04     77      
  DavidsonSolver iter   6:   3.67548e-08  6.20551e-05     78      Converged

Using PsiAPI¶

Using PsiAPI is as simple as using Psithon. The relevant function is tdscf_excitations, which is the same driver called under the hood when running TDSCF from a Psithon input file. In this example, we ran the same example as before, but using the random-phase approximation (RPA). Note that the save_jk option has to be explicitly set to True when working in PsiAPI.

The res variable contains all results from the TDSCF calculation. It is a list, with one element per computed root.

Known limitations¶

The implementation cannot currently handle the following cases:

• The use of symmetry with density-fitted two-electron integrals.
• Excited states of triplet symmetry from a restricted DFT reference.
• Excited states from an unrestricted reference other than HF or LDA.

Theory and implementation¶

You can find additional material on response theory in Norman2011-ad and Dreuw2005-wp.

The time-dependent Schrödinger equation is the appropriate equation of motion: \begin{equation}\label{eq:td-schrodinger} H|0(t)\rangle = \mathrm{i}\frac{\partial |0(t)\rangle}{\partial t}, \end{equation} where the time-dependent, semiclassical matter-field Hamiltonian is used: \begin{equation} H = H_0 + V(t). \end{equation} The perturbation is a one-electron operator that is periodic in time $V(t) = V(t+T)$. Rather than solve for the time-dependent wavefunction, we can target the order-by-order corrections of time-dependent expectation values: \begin{equation} \begin{aligned} \langle 0(t) | A | 0(t) \rangle &= \langle 0 | A | 0 \rangle

\sum_{(b, B)} \langle \langle A; B \rangle \rangle_{\omega_{b}}\epsilon_{B}(\omega_{b})\exp[-\mathrm{i}\omega_{b}t] \ &+ \frac{1}{2} \sum_{(b, B);(c, C)} \langle \langle A; B, C \rangle \rangle_{\omega_{b}, \omega_{c}}\epsilon_{B}(\omega_{b})\epsilon_{C}(\omega_{c})\exp[-\mathrm{i}(\omega_{b}+\omega_{c})t]

• \ldots

\end{aligned} \end{equation} The discrete Fourier coefficients are the response functions which we can map to observable properties of the system. For example, the linear response function: \begin{equation} \langle \langle A; B \rangle \rangle_{\omega} = \sum_{i \neq j} \frac{\langle i | A | j \rangle\langle j | B | i \rangle}{\omega - \omega_{ij}}

\frac{\langle i | B | j \rangle\langle j | A | i \rangle}{\omega + \omega_{ij}} \end{equation} in the exact eigenbasis of the unperturbed Hamiltonian and it encodes all excitation energies of the system (the poles) and transition properties (the residues). The sum-over-states (SOS) expression requires the whole eigenbasis of the Hamiltonian. However, response properties can be extracted without prior knowledge of a suitable basis, as long as one has access to the molecular electronic Hessian.

In self-consistent field theories, excitation energies of molecular systems can be obtained as the eigenvalues, $\omega_{n}$, of the response eigenvalue problem: \begin{equation} \begin{pmatrix} \mathbf{A} & \mathbf{B} \\ \mathbf{B}^{*} & \mathbf{A}^{*} \end{pmatrix} \begin{pmatrix} \mathbf{X}_{n} \\ \mathbf{Y}_{n} \end{pmatrix} = \omega_{n} \begin{pmatrix} \mathbf{1} & \mathbf{0} \\ \mathbf{0} & -\mathbf{1} \end{pmatrix} \begin{pmatrix} \mathbf{X}_{n} \\ \mathbf{Y}_{n} \end{pmatrix}. \end{equation}

The $\mathbf{A}$ and $\mathbf{B}$ matrices appearing on the left-hand side are the blocks of the molecular electronic Hessian: \begin{equation} \begin{aligned} A_{aibj} &= \delta_{ij}f_{ab} - \delta_{ab}f_{ij} + (ai|jb) - \gamma (ab|ji) + w_{\mathrm{xc}; aijb} \\ B_{aibj} &= (ai|bj) - \gamma (aj|ib) + w_{\mathrm{xc}; aibj} \\ \end{aligned} \end{equation} 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.

TThe response eigenvalue equation is not an Hermitian eigenproblem, due to the nonunit metric on the right-hand side. The eigenproblem however has Hamiltonian symmetry: the roots appear in pairs $(\omega_{n}, -\omega_{n})$, as do the eigenvectors. When using real orbitals, this symmetry can be preserved by solving an equivalent problem instead. First note that: \begin{equation} (\mathbf{A} - \mathbf{B})(\mathbf{A} + \mathbf{B})| \mathbf{X}_{n} + \mathbf{Y}_{n}\rangle = \omega^{2}_{n} | \mathbf{X}_{n} + \mathbf{Y}_{n}\rangle, \end{equation} which can be brought to the Hermitian form: \begin{equation} (\mathbf{A} - \mathbf{B})^{\frac{1}{2}}(\mathbf{A} + \mathbf{B})(\mathbf{A} - \mathbf{B})^{\frac{1}{2}} \mathbf{T}_{n} = \omega^{2}_{n} \mathbf{T}_{n}, \end{equation} assuming the SCF reference is stable, i.e. $(\mathbf{A}-\mathbf{B})$ is positive-definite. The paired vectors $| \mathbf{X}_{n} - \mathbf{Y}_{n}\rangle$ are left eigenvectors and form a biorthonormal set together with the right eigenvectors $| \mathbf{X}_{n} + \mathbf{Y}_{n}\rangle$. In the Tamm-Dancoff approximation (TDA), we solve the Hermitian eigenvalue problem: \begin{equation} \mathbf{A}\mathbf{X}_{n} = \omega_{n}\mathbf{X}_{n}, \end{equation} using a regular Davidson solver.

Transition properties¶

The excitation energies tell only one part of the story. We need to have transition moments to understand whether an excitation is allowed (bright) or not (dark).

The excitation energies and eigenvectors can be used to compute transition moments, such as electric and magnetic transition dipole moments, and spectroscopic intensities, such as oscillator strengths and rotatory strengths. For example, the transition electric dipole moment in the length gauge is:

Our implementation calculates the following transition properties:

• electric dipole moment in the length gauge
• oscillator strength in the length gauge
• electric dipole moment in the velocity gauge
• oscillator strength in the velocity gauge
• magnetic dipole moment
• rotatory strength in the length gauge

upon convergence of the solver. It is rather easy to compute new transition properties. The following is a re-implementation of the oscillator strength:

Plotting spectra¶

We can use these observables to plot spectra. We use empirical lineshapes to simulate the broadening observed in experiments. For one-photon absorption (OPA), the function to plot is: $$ \varepsilon(\omega) = \frac{4\pi^{2}N_{\mathrm{A}}\omega}{3\times 1000\times \ln(10) (4 \pi \epsilon_{0}) n \hbar c} \sum_{i \rightarrow j}g_{ij}(\omega)|\mathbf{\mu}_{ij}|^{2} $$ while for electronic circular dichroism (ECD): $$ \Delta\varepsilon(\omega) = \frac{16\pi^{2}N_{\mathrm{A}}\omega}{3\times 1000\times \ln(10) (4 \pi \epsilon_{0}) n \hbar c^{2}} \sum_{i \rightarrow j}g_{ij}(\omega)\Im(\mathbf{\mu}_{ij}\cdot\mathbf{m}_{ij}). $$ See Rizzo2011-to for details.

I decided to decouple the lineshape fitting from the plotting:

• We do not entangle Psi4 with any specific plotting library.
• The user has the freedom to run whatever visualization pipeline.

Psi4 provides the spectrum function:

def spectrum(*,
             poles: Union[List[float], np.ndarray],
             residues: Union[List[float], np.ndarray],
             kind: str = "opa",
             lineshape: str = "gaussian",
             gamma: float = 0.2,
             npoints: int = 5000,
             out_units: str = "nm") -> Dict[str, np.ndarray]

that will perform the lineshape fitting of the TDDFT results and return numpy arrays for the $x$ and $y$ values of the plot.

For both one-photon absorption (OPA) and electronic circular dichroism (ECD) the poles are the excitation energies. The residues are:

• The square of the transition electric dipole moment for OPA: $|\mu_{ij}|^{2}$.
• The rotatory strength $R_{ij}$ for ECD. This is the imaginary part of the product of transition electric and magnetic dipole moments: $\Im(\mathbf{\mu}_{ij}\cdot\mathbf{m}_{ij})$.

The implementation assumes these quantities are given in the length gauge.

We can obtain these quantities in atomic units from the list of results returned by the tdscf_excitations function.

With these at hand, we invoke the spectrum function to perform the convolution with Gaussian lineshapes (inhomogeneous broadening) with empirical linewidth $\gamma = 0.01\,\textrm{a.u.}$ (angular frequency).

All data in input is given in atomic units, but the final spectrum will use nanometers on the $x$ axis and macroscopic units $\mathrm{L}\cdot\mathrm{mol}^{-1}\cdot\mathrm{cm}^{-1}$ on the $y$ axis.

The return value is a dictionary-of-dictionaries with:

1. The lineshape convolution:
   • opa_spectrum["convolution"]["x"] a NumPy array with the $x$ axis points of the lineshape convolution.
   • opa_spectrum["convolution"]["y"] a NumPy array with the $y$ axis points of the lineshape convolution.
2. The infinitely narrow (stick) representation:
   • opa_spectrum["sticks"]["poles"] a NumPy array with the poles in the chosen unit of measure.
   • opa_spectrum["sticks"]["residues"] a NumPy array with the residues in macroscopic units.

We can now plot these spectra with any library we want. In this example, I will use Altair and pandas. The raw data is stored in NumPy arrays and any plotting library will work!

Solvent effects¶

Finally, it is possible to run these calculations including solvent effects, either with the polarizable continuum model (PCM) or with polarizable embedding (PE).

And we can plot the results, of course!