This is an implementation of Time-Dependent Density Functional Theory(TDDFT) using the Tamm-Dancoff Approximation (TDA).

I could have used a cube as a space grid and then taken Riemann sums to evaluate (there's a great YouTube series by James Johns where he develops a matlab HF program and shows how to convert it to DFT. In this he uses Riemann sums to evaluate the integrals in the DFT program.) However, I decided to try for a 'proper' atom centered spherical grid approach. A useful paper was PMW Gill, BG Johnson and JA Poples 'A standard grid for density functional theory', although I didn't use this SG-1 grid the paper helped understand the techniques involved. The grid I settled on was a coarse grid of (10,15) radial points for period 1 and period 2 elements respectively. The radial grid is a Mura-Knowles radial grid ME Mura and PJ Knowles 'Improved radial grids for quadrature in density-functional calculations' JCP 104, 9848 (1996); DOI:10.1063/1.471749. The 'coarse' angular grid is of Lebedev orders (11, 15) for period 1 and period 2 respectively. This translates into 50 and 86 points respectively arranged on a spherical shell (VI Lebedev, and DN Laikov, Doklady Mathematics, 'A Quadrature formula for the sphere of the 131st algebraic order of accuracy' Vol. 59, No. 3, (1999)). There are various sources for this data given in the external links of the wikipedia article on Lebedev integration. A pruning scheme is employed to systematically reduce the number of angular points in regions where dense angular quadrature is not necessary, such as near the nuclei where the charge density is approximately spherically symmetric and at long distance from the nucleus. The pruning scheme I employed was the Treutler-Aldrich scheme O Treutler and R Ahlrich, 'Efficient molecular numerical integration schemes' JCP 102, 346 (1995); DOI:10.1063/1.469408. The partitioning of the atomic centered grids to a molecular grid follows a Becke scheme after Stratmann RE Stratmann, GE Scuseria and MJ Frisch, 'Achieving Linear scaling in exchange-correlation density functional quadratures' CPL 257, 3-4 (1996); DOI:10.1016/009-2614(96)00600-8. Finally I have implemented a final radius adjustment during the partition (Becke suggests doing this) using the Bragg radius. A second 'close' grid is also included which is a (50, 75) radial and (29, 29) angular, the latter representing 302 points on each shell. The grid routines are in ks_grid.py.

To get the DFT SCF started we need an initial density. To do this I use a HF overlap matrix S, and an initial Fock matrix composed of the sum of the 1-electron kinetic and coulomb integrals (core Hamiltonian - T+V). This Fock is then orthogonalised (F^') as (S^-0.5)^TFS^-0.5, eigensolve the resulting orthogonal Fock for orbital coefficients C orthogonal, transform back to atomic basis as S^-0.5C^', use resulting ao coefficients to compute a density matrix D_μν = c_μic_iν where i is over occupied orbitals. This initial density can be used with initial Fock and 2-electron repulsion integrals to form the coulomb integral J (we don't want the HF exchange integral K for DFT). To get these integrals I've used a modified version of Harpy's Cython integral package aello. This version is slightly different from the version in kspy_lda in that the dipole routine returns the component matrices rather than the actual dipole, additionally the angular and nabla routines have been added. These are in ks_aello.pyx.

The molecule definition is contained in a mol object which is itself comprised of objects from an atom class. Each instance of the atom class contains the atom symbol, atomic number and the coordinates of the atom center (array[3]). The molecule is hard coded as H₂O. The basis is contained in an orb object which is itself comprised of objects from a gaussian class. Each instance of the gaussian class contains the atom the Gaussian is centered on, the momentum(array[3]), the exponents (array[primatives]), the coefficients (array[primatives]), the normalisation (array[primatives]) and a copy of the atom center coordinates (array[3]). The momenta are given as s [0,0,0] px [1,0,0] py [0,1,0] and pz [0,0,1]. The basis used is a simple STO-3G so we only require s and p orbitals. The primatives exponent and coefficient values are hard-coded in the main section. (I use the psi4 format of the basis sets from BSE which have some (small) differences from the nwchem format versions as used by eg pyscf. This might lead to numerical differences in values when using high precision).

The choice of functionals here was determined solely because they have easily determined analytic derivatives. The second derivatives of the exchange-correlation energy are needed in TDDFT to determine the orbital Hessian for the coupling matrix. The exchange functional is Slater LDA and the correlation functional is RPA. For TDDFT we are working in a molecular spin basis so will use spin polarized versions of the functional (with α = β). The derivatives used are given below

Details of TDDFT can be found in Time-dependent density-functional theory for molecules and molecular solids, ME Casida, Journal of Molecular Structure: THEOCHEM 914 (2009) 3–18 and Single-Reference ab Initio Methods for the Calculation of Excited States of Large Molecules, A Dreuw and M Head-Gordon, Chem. Rev.2005,105,4009−4037. The basic equations for the excitation (de-excitation) energies are

The calculation of singlet and triplet states is analogous to the spin-adapted CIS calculation in HF theory see. The class LR_TDA is provided to compute the excitation energies and coefficients for either singlet or triplet states. The class provides a response property which is a tuple (energy, coefficients).

An example of calculating transition properties in TDDFT is given in this psicon 2020 document. We calculate the electric transition dipoles in both length (μ) and velocity (∇) gauges together with the associated oscillator strengths. Additionally the magnetic transition dipoles are calculated in the length (L) gauge and the rotary strengths in both gauges. All the previous properties are returned by the transition_properties method of the TDA_properties class.
A basic transition natural orbital treatment is available from the transition_NO method of the TDA_properties class. Martin, R. L., Journal of Chemical Physics, 118, 4775-4777.
A spectrum method of the TDA_properties class is provided to plot the spectrum of oscillator strengths in both gauges. This is all provided in the module ks_tda. See results.md for more details.

An example of doing a real-time calculation is given in the class RT_TDA. The Magnus expansion is used to propogate a complex Hamiltonian forward in time.The Magnus expansion used here is a second-order Magnus. The observable in the propogation is the electric dipole in length gauge. The external field is applied as either a 'kick' (instantaneous pulse) or a Gaussian waveform. See results.md for more details. Analysis of the dipole propogation in the frequency domain is possible with the spectrum method which uses Pade Approximants and a more general article on Pade can be found here.

Copy all files to a directory, run

python3 setup.py build_ext --inplace install --user

then

python3 ks_main.py