Efficient all-electron time-dependent density functional theory calculations using an enriched finite element basis

TDDFT reduces the many-electron time dependent Schrödinger equation to an auxiliary system of non-interacting electrons which yields the same time dependent electron density, $\rho(\boldsymbol{\textbf{r}},t)$, as the interacting system and whose evolution is prescribed by a set of effective single electron equations, called the time dependent Kohn-Sham (TDKS) equations. The TDKS equations, in atomic units, are given as

$$\begin{split}i\frac{\partial{\psi_{\alpha}(\boldsymbol{\textbf{r}},t)}}{\partial{t}}&=H_{\rm{KS}}[\rho](\boldsymbol{\textbf{r}},t;\boldsymbol{\textbf{R}})\psi_{\alpha}(\boldsymbol{\textbf{r}},t)\coloneqq\left[-\frac{1}{2}\nabla^{2}+v_{\rm{KS}}[\rho](\boldsymbol{\textbf{r}},t;\boldsymbol{\textbf{R}})\right]\psi_{\alpha}(\boldsymbol{\textbf{r}},t)\\ &\coloneqq\left[-\frac{1}{2}\nabla^{2}+v_{\rm{ext}}(\boldsymbol{\textbf{r}},t;\boldsymbol{\textbf{R}})+v_{\rm{H}}[\rho](\boldsymbol{\textbf{r}},t)+v_{\rm{xc}}[\rho](\boldsymbol{\textbf{r}},t)\right]\psi_{\alpha}(\boldsymbol{\textbf{r}},t)\,.\end{split}$$

(1)

In the above, $H_{\rm{KS}}[\rho](\boldsymbol{\textbf{r}},t;\boldsymbol{\textbf{R}})$ and $\psi_{\alpha}(\boldsymbol{\textbf{r}},t)$ represent the time-dependent Kohn-Sham Hamiltonian and orbitals, respectively; $\alpha$ is the index for the $N_{e}$ electrons in the system; $\boldsymbol{\textbf{R}}=\{\boldsymbol{\textbf{R}}_{1},\boldsymbol{\textbf{R}}_{2},...,\boldsymbol{\textbf{R}}_{N_{a}}\}$ is the collective representation for the positions of the $N_{a}$ atoms in the system. The electron density, $\rho(\boldsymbol{\textbf{r}},t)$, is defined by the the Kohn-Sham orbitals as

$$\rho(\boldsymbol{\textbf{r}},t)=\sum_{\alpha=1}^{N_{e}}|\psi_{\alpha}(\boldsymbol{\textbf{r}},t)|^{2}\,.$$

(2)

For simplicity, we consider only spin-unpolarized systems. However, the ideas presented can be easily extended to spin-polarized systems. This work does not account for relativistic effects, which will be the subject of a future work. The external potential, $v_{\rm{ext}}(\boldsymbol{\textbf{r}},t;\boldsymbol{\textbf{R}})$, comprises of the nuclear potential ($v_{\rm{N}}(\boldsymbol{\textbf{r}};\boldsymbol{\textbf{R}}$)) and the potential corresponding to an external field ($v_{\rm{f}}(\boldsymbol{\textbf{r}},t)$). The nuclear potential is given as

$$v_{\rm{N}}(\boldsymbol{\textbf{r}};\boldsymbol{\textbf{R}})=-\sum\limits_{I=1}^{N_{a}}\frac{Z_{I}}{|\boldsymbol{\textbf{r}}-\boldsymbol{\textbf{R}}_{I}|}\,,$$

(3)

where $Z_{I}$ is the atomic number of the $I-$th atom. $v_{\rm{f}}(\boldsymbol{\textbf{r}},t)$ is usually defined as a monochromatic laser pulse of the form

$$v_{\rm{f}}(\boldsymbol{\textbf{r}},t)=-\mathbf{E_{0}}(t)\cdot\boldsymbol{\textbf{r}}\,,$$

(4)

where $\mathbf{E_{0}}(t)$ represents the time-dependent electric field. The Hartree potential, $v_{\rm{H}}[\rho](\boldsymbol{\textbf{r}},t)$, denotes the classical electrostatic potential corresponding to the electron density, given as

$$v_{\rm{H}}[\rho](\boldsymbol{\textbf{r}},t)=\int{\frac{\rho(\boldsymbol{\textbf{r}}^{\prime},t)}{|\boldsymbol{\textbf{r}}-\boldsymbol{\textbf{r}}^{\prime}|}\,d\boldsymbol{\textbf{r}}^{\prime}}\,.$$

(5)

The exchange-correlation potential, $v_{\rm{xc}}[\rho](\boldsymbol{\textbf{r}},t)$, is the mean-field potential that accounts for the quantum many-electron interactions. Although, in general, the exchange-correlation potential is known to be nonlocal in both space and time ^{100, 101, 102} along with initial state dependence, most widely used approximations use locality in time (adiabatic approximation) and non-dependence on the initial many-electron wavefunction, for lack of knowledge of the time nonlocality. As a result, it is common to employ the existing exchange-correlation approximations used in ground-state DFT. In this work, we use the adiabatic local-density approximation (ALDA) ¹⁰³, which is local in both space and time. Specifically, we use the Ceperley-Alder form ¹⁰⁴.

We note that both $v_{\rm{N}}$ and $v_{\rm{H}}$ are extended in real space, and for an evaluation in real-space formulations, they can be recast as solutions to the following Poisson equations ^{77, 79, 81, 84, 105}:

	$$-\frac{1}{4\pi}\nabla^{2}v_{\rm{H}}(\boldsymbol{\textbf{r}},t)=\rho(\boldsymbol{\textbf{r}},t)\,,\quad v_{\rm{H}}(\boldsymbol{\textbf{r}},t)\rvert_{\partial\Omega}=f(\boldsymbol{\textbf{r}},\boldsymbol{\textbf{R}})\,,$$		(6a)
	$$-\frac{1}{4\pi}\nabla^{2}v_{\rm{N}}(\boldsymbol{\textbf{r}};\boldsymbol{\textbf{R}})=b(\boldsymbol{\textbf{r}},\boldsymbol{\textbf{R}})\quad v_{\rm{N}}(\boldsymbol{\textbf{r}})\rvert_{\partial\Omega}=-f(\boldsymbol{\textbf{r}},\boldsymbol{\textbf{R}})\,,$$		(6b)

where $b(\boldsymbol{\textbf{r}},\boldsymbol{\textbf{R}})=-\sum\limits_{I=1}^{N_{a}}{Z_{I}\delta(\left|\boldsymbol{\textbf{r}}-\boldsymbol{\textbf{R}}_{I}\right|)}$ with $\delta(x)$ denoting a Dirac delta distribution; $f(\boldsymbol{\textbf{r}},R)=\sum_{I=1}^{N_{a}}\frac{Z_{I}}{\left|\boldsymbol{\textbf{r}}-\boldsymbol{\textbf{R}}_{I}\right|}$; $\partial\Omega$ denotes the boundary of a sufficiently large bounded domain $\Omega\in\mathcal{R}^{3}$.

In this section, we present the EFE discretization of the TDKS equations. At the heart of the EFE discretization lies the augmentation of the CFE basis with compactly supported numerical atomic functions, termed as enrichment functions. As a result, we account for the sharp variations in the electronic fields close to nuclei, largely, through the enrichment functions and mitigate the need for a refined classical finite element mesh. The initial attempt at an EFE basis ⁹⁷ for groundstate DFT used smoothly truncated single atom orbitals and potentials as the enrichment functions, which resulted in a staggering $50-100\times$ speedup over the CFE basis. However, such enrichment functions are susceptible to being linearly dependent on the CFE basis, particularly while using refined FE meshes, affecting the accuracy and robustness of the resulting EFE basis. To that end, we developed an orthogonalized EFE basis ⁹⁸, wherein the enrichment functions are orthogonalized with respect to the underlying CFE basis, without comprising on the locality of the resultant basis. This work employs the more robust orthogonalized EFE basis. Hereafter, for succinctness, we omit the qualifier orthogonalized and refer to the orthogonalized EFE basis simply as EFE basis. The EFE discretization of the Kohn-Sham orbital, $\psi_{\alpha}(\boldsymbol{\textbf{r}},t)$, is given as

$$\psi_{\alpha}^{h}(\boldsymbol{\textbf{r}},t)=\underbrace{\sum_{i=1}^{n_{h}}N^{C}_{i}(\boldsymbol{\textbf{r}})\psi_{\alpha,i}^{C}(t)}_{\text{Classical}}+\underbrace{\sum_{I=1}^{N_{a}}\sum_{j=1}^{n_{I}}N^{E,\psi}_{j,I}(\boldsymbol{\textbf{r}})\psi_{\alpha,j,I}^{E}(t)}_{\text{Enriched}}\,,$$

(7)

where the superscript ‘$h$’ denotes a discrete field; {$N^{C}_{i}(\boldsymbol{\textbf{r}})$} and $\{\psi_{\alpha,i}^{C}(t)\}$ denote the classical FE basis functions and their corresponding time-dependent coefficients; similarly, {$N^{E,\psi}_{j,I}(\boldsymbol{\textbf{r}})$} and $\{\psi_{\alpha,j,I}^{E}(t)\}$ denote the enrichment basis functions and their respective time-dependent coefficients. In the above form, the $I$^th atom at $\boldsymbol{\textbf{R}}_{I}$, contributes $n_{I}$ enrichment functions, each centered around $\boldsymbol{\textbf{R}}_{I}$. The enrichment function {$N^{E,\psi}_{j,I}(\boldsymbol{\textbf{r}})$} is expressed as

$$N^{E,\psi}_{j,I}(\boldsymbol{\textbf{r}})=N^{A,\psi}_{j,I}(\boldsymbol{\textbf{r}})-N^{B,\psi}_{j,I}(\boldsymbol{\textbf{r}})\,,$$

(8)

where $N^{A,\psi}_{j,I}(\boldsymbol{\textbf{r}})$ and $N^{B,\psi}_{j,I}(\boldsymbol{\textbf{r}})$ are the atomic and orthogonalizing parts, respectively. The atomic part, $N^{A,\psi}_{j,I}(\boldsymbol{\textbf{r}})$, is given as

$$N^{A,\psi}_{j,I}(\boldsymbol{\textbf{r}})=\psi_{nlm,I}(|\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{R}}_{I}|,\beta_{\boldsymbol{\textbf{R}}_{I}},\gamma_{\boldsymbol{\textbf{R}}_{I}})u(|\boldsymbol{\textbf{x}}-\boldsymbol{\textbf{R}}_{I}|,r_{0},s)\,,$$

(9)

where $\psi_{nlm,I}$ is an atomic groundstate Kohn-Sham orbital indexed by the principal quantum number $n$, azimuthal quantum number $l$, and magnetic quantum number $m$, for an isolated atom of the atom type of the $I^{\text{th}}$ atom, defined in spherical coordinates. The function $u(r,r_{0},s)$ is a smooth cutoff function, parameterized by a cutoff radius $r_{0}$ and smoothness factor $s$. Typically, we include all the occupied and a few unoccupied groundstate single atom orbitals as the enrichment functions. The orthogonalizing part, $N^{B,\psi}_{j,I}(\boldsymbol{\textbf{r}})$, is given as a linear combination of the underlying CFE basis, i.e.,

$$N^{B,\psi}_{j,I}(\boldsymbol{\textbf{r}})=\sum_{l=1}^{n_{h}}c_{j,I,l}N^{C}_{l}(\boldsymbol{\textbf{r}})\,,$$

(10)

where $c_{j,I,l}$ are evaluated such that they satisfy the orthogonality condition $\int_{\Omega}N^{E,\psi}_{j,I}(\boldsymbol{\textbf{r}})N^{C}_{i}(\boldsymbol{\textbf{r}})\,d\boldsymbol{\textbf{r}}=0$. We refer to ⁹⁸ for a detailed description on $N^{A,\psi}_{j,I}(\boldsymbol{\textbf{r}})$, $N^{B,\psi}_{j,I}(\boldsymbol{\textbf{r}})$, and $h(r,r_{0},s)$. Hereafter, for compact notation, we combine the $\{j,I\}$ indices in $N^{E,\psi}_{j,I}(\boldsymbol{\textbf{r}})$ and $\psi_{\alpha,j,I}$ into a single index $\nu$.

Using the EFE discretization of Eq. 7 in the TDKS equations (cf. Eq. 1) yields

$$i\boldsymbol{\textbf{M}}\dot{\boldsymbol{\psi}}_{\alpha}(t)=\boldsymbol{\textbf{H}}\boldsymbol{\psi}_{\alpha}(t)\,.$$

(11)

where H and M are the discrete Kohn-Sham Hamiltonian matrix and the overlap matrix, respectively, and $\boldsymbol{\psi}_{\alpha}$ denotes the vector containing the $\psi_{\alpha,i}^{C}(t)$ and $\psi_{\alpha,\nu}^{E}(t)$ coefficients of Eq. 7. The Kohn-Sham Hamiltonian H is given by

$$H_{mn}=\frac{1}{2}\int_{\Omega}\nabla N_{m}(\boldsymbol{\textbf{r}})\cdot\nabla N_{n}(\boldsymbol{\textbf{r}})\,d\boldsymbol{\textbf{r}}+\int_{\Omega}v_{\rm{KS}}[\rho](\boldsymbol{\textbf{r}},t;\boldsymbol{\textbf{R}})N_{m}(\boldsymbol{\textbf{r}})N_{n}(\boldsymbol{\textbf{r}})\,d\boldsymbol{\textbf{r}}\,,$$

(12)

where $N_{m}(\boldsymbol{\textbf{r}})$ and $N_{n}(\boldsymbol{\textbf{r}})$ are generic representations for $N^{C}_{i}(\boldsymbol{\textbf{r}})$ and $N^{E,\psi}_{\nu}(\boldsymbol{\textbf{r}})$. The overlap matrix M, owing to the orthogonality between $N^{C}_{i}(\boldsymbol{\textbf{r}})$ and $N^{E,\psi}_{\nu}(\boldsymbol{\textbf{r}})$, has the following block-diagonal structure

$$\boldsymbol{\textbf{M}}=\left[\begin{array}[]{c|c}\boldsymbol{\textbf{M}}^{\textbf{cc}}&0\\ \hline\cr 0&\boldsymbol{\textbf{M}}^{\textbf{ee}}\end{array}\right]\,,$$

(13)

where $\boldsymbol{\textbf{M}}^{\textbf{cc}}_{jk}=\int_{\Omega}N^{C}_{j}(\boldsymbol{\textbf{r}})N^{C}_{k}(\boldsymbol{\textbf{r}})\,d\boldsymbol{\textbf{r}}$ denotes the overlap between two CFE basis functions and $\boldsymbol{\textbf{M}}^{\textbf{ee}}_{\mu,\nu}=\int_{\Omega}N^{E,\psi}_{\mu}(\boldsymbol{\textbf{r}})N^{E,\psi}_{\nu}(\boldsymbol{\textbf{r}})$ denotes the overlap between two enrichment functions. Noting the fact that M is a positive definite matrix with a unique positive definite square root, $\boldsymbol{\textbf{M}}^{1/2}$, we reformulate Eq. 11 as

$$i\dot{\bar{\boldsymbol{\psi}}}_{\alpha}(t)=\bar{\boldsymbol{\textbf{H}}}\bar{\boldsymbol{\psi}}_{\alpha}(t)\,,$$

(14)

where $\bar{\boldsymbol{\psi}}_{\alpha}(t)=\boldsymbol{\textbf{M}}^{1/2}\boldsymbol{\psi}_{\alpha}(t)$ and $\bar{\boldsymbol{\textbf{H}}}=\boldsymbol{\textbf{M}}^{-1/2}\boldsymbol{\textbf{H}}\boldsymbol{\textbf{M}}^{-1/2}$ are the representation of $\psi_{\alpha}(\boldsymbol{\textbf{r}},t)$ and $H_{\rm{KS}}$ in a Löwdin orthonormalized basis  ¹⁰⁶. We emphasize that the above transformation of the discrete TDKS equation is crucial to our use of the Magnus propagator (see Sec. 2.5).

Similar to $\psi_{\alpha}(\boldsymbol{\textbf{r}},t)$, the EFE discretization of the electrostatic potential (both nuclear and Hartree potentials), $\phi(\boldsymbol{\textbf{r}},t)$, can be written as

$$\phi^{h}(\boldsymbol{\textbf{r}},t)=\underbrace{\sum_{j=1}^{n_{h}}N^{C}_{j}(\boldsymbol{\textbf{r}})\phi_{j}^{C}(t)}_{\text{Classical}}+\underbrace{\sum_{I=1}^{N_{a}}N^{E,\phi}_{I}(\boldsymbol{\textbf{r}})\phi_{I}^{E}(t)}_{\text{Enriched}}\,.$$

(15)

We note that for the nuclear potential, $\phi(\boldsymbol{\textbf{r}},t)$ does not have any time dependence. In the above, the enrichment function $N^{E,\phi}_{I}(\boldsymbol{\textbf{r}})$, analogous to $N^{E,\psi}_{\nu}(\boldsymbol{\textbf{r}})$, comprises of an atomic part ($N^{A,\phi}_{I}(\boldsymbol{\textbf{r}})$) and an orthogonalizing part ($N^{B,\phi}_{I}(\boldsymbol{\textbf{r}})$). The atomic part, $N^{A,\phi}_{I}(\boldsymbol{\textbf{r}})$, is given as a smoothly truncated spherically symmetric electrostatic potential (nuclear or Hartree) of an isolated atom of the same type as located at $\boldsymbol{\textbf{R}}_{I}$. On the other hand, the orthogonalizing part, $N^{B,\phi}_{I}(\boldsymbol{\textbf{r}})$, is given as a linear combination of the CFE basis {$N^{C}_{i}(\boldsymbol{\textbf{r}})$}, such that it guarantees the orthogonality condition: $\int_{\Omega}N^{E,\phi}_{I}(\boldsymbol{\textbf{r}})N^{C}_{i}(\boldsymbol{\textbf{r}})\,d\boldsymbol{\textbf{r}}=0$. Using Eq. 15 in Eq. 6 yields

$$\boldsymbol{\textbf{A}}\boldsymbol{\phi}(t)=\boldsymbol{\textbf{d}}(t)\,,$$

(16)

where $\boldsymbol{\textbf{A}}_{mn}=\frac{1}{4\pi}\int_{\Omega}\nabla N_{m}(\boldsymbol{\textbf{r}})\cdot N_{n}(\boldsymbol{\textbf{r}})$, with $N_{m}(\boldsymbol{\textbf{r}})$ and $N_{n}(\boldsymbol{\textbf{r}})$ being generic representation of $N^{C}_{i}(\boldsymbol{\textbf{r}})$ and $N^{E,\phi}_{I}(\boldsymbol{\textbf{r}})$, is the discrete Laplace operator; and $\boldsymbol{\textbf{d}}_{m}(t)=\int_{\Omega}\rho(\boldsymbol{\textbf{r}},t)N_{m}(\boldsymbol{\textbf{r}})$ for the Hartree potential or $\boldsymbol{\textbf{d}}_{m}(t)=\int_{\Omega}b(\boldsymbol{\textbf{r}};\boldsymbol{\textbf{R}})N_{m}(\boldsymbol{\textbf{r}})$ for the nuclear potential.

We remark that the Löwdin transformation in Eq. 14 warrants an efficient means of evaluating $\boldsymbol{\textbf{M}}^{1/2}$ and $\boldsymbol{\textbf{M}}^{-1/2}$. To that end, we use spectral finite elements in conjunction with Gauss-Lobatto-Legendre (GLL) quadrature rule, the combination of which renders classical-classical block ($\boldsymbol{\textbf{M}}^{\textbf{cc}}$) block of M diagonal, and, hence, allows for trivial evaluation of $\left(\boldsymbol{\textbf{M}}^{\textbf{cc}}\right)^{1/2}$ and $\left(\boldsymbol{\textbf{M}}^{\textbf{cc}}\right)^{-1/2}$. We refer to ⁸⁴ for an elaborate discussion on spectral finite elements. The enriched-enriched block ($\boldsymbol{\textbf{M}}^{\textbf{ee}}$), on the other hand, is a small diagonally dominant matrix of size $\sim N_{e}\times N_{e}$, and, hence, $\left(\boldsymbol{\textbf{M}}^{\textbf{ee}}\right)^{1/2}$ and $\left(\boldsymbol{\textbf{M}}^{\textbf{ee}}\right)^{-1/2}$ can be easily evaluated through either eigenvalue decomposition or Newton-Schultz methods  ^{107, 108, 109}. Subsequently, $\boldsymbol{\textbf{M}}^{1/2}$ and $\boldsymbol{\textbf{M}}^{-1/2}$ can be written as

$$\boldsymbol{\textbf{M}}^{1/2}=\left[\begin{array}[]{c|c}(\boldsymbol{\textbf{M}}^{\textbf{cc}})^{1/2}&0\\ \hline\cr 0&\left(\boldsymbol{\textbf{M}}^{\textbf{ee}}\right)^{1/2}\end{array}\right]\,,\quad\boldsymbol{\textbf{M}}^{-1/2}=\left[\begin{array}[]{c|c}(\boldsymbol{\textbf{M}}^{\textbf{cc}})^{-1/2}&0\\ \hline\cr 0&\left(\boldsymbol{\textbf{M}}^{\textbf{ee}}\right)^{-1/2}\end{array}\right]\ \,.$$

(17)

The enrichment functions, $N^{E,\psi}_{\nu}(\boldsymbol{\textbf{r}})$ and $N^{E,\phi}_{I}(\boldsymbol{\textbf{r}})$, being characterized by sharp gradients or oscillations near the atoms, warrant the use of a high density of quadrature points near the atoms for an accurate evaluation of the integrals involving them. On the other hand, given that the enrichment functions have a small compact support, a lower quadrature density may suffice in regions farther away from atoms. We strike a good balance of accuracy and efficiency by using an adaptive quadrature, wherein we recursively refine each finite element until a set of trial integrals, involving the enrichment functions, attain convergence. To this end, we adopt the adaptive quadrature scheme proposed in  ^{110, 97} and refer to these works for the details.

We now discuss the time discretization of Eq. 14 in the context of second-order Magnus propagator. In the Magnus ansatz, the solution to Eq.  14 is given as

$$\bar{\boldsymbol{\psi}}_{\alpha}(t)=\text{exp}(\boldsymbol{\textbf{A}}(t))\bar{\boldsymbol{\psi}}_{\alpha}(0)\,,\qquad\forall t\geq 0\,,$$

(18)

where $\text{exp}(\boldsymbol{\textbf{A}}(t))$ is called the Magnus propagator with $\boldsymbol{\textbf{A}}(t)$ given explicitly as  ^{111, 112}

$$\boldsymbol{\textbf{A}}(t)=\int_{0}^{t}-i\bar{\boldsymbol{\textbf{H}}}(\tau)d\tau-\frac{1}{2}\int_{0}^{t}\left[\int_{0}^{\tau}-i\bar{\boldsymbol{\textbf{H}}}(\sigma)d\sigma,-i\bar{\boldsymbol{\textbf{H}}}(\tau)\right]d\tau+\ldots\,,$$

(19)

with $[\boldsymbol{\textbf{X}},\boldsymbol{\textbf{Y}}]=\boldsymbol{\textbf{X}}\boldsymbol{\textbf{Y}}-\boldsymbol{\textbf{Y}}\boldsymbol{\textbf{X}}$ being the commutator. Given the difficulty in accurately resolving the implicit dependence of $\bar{\boldsymbol{\textbf{H}}}(t)$ on $\rho(\boldsymbol{\textbf{r}},t)$, for practical use, we rewrite the Magnus propagator as

$$\text{exp}(\boldsymbol{\textbf{A}}(t))=\prod_{n=1}^{N}\exp(\boldsymbol{\textbf{A}}_{n})\,,$$

(20)

where $\boldsymbol{\textbf{A}}_{n}$ is given by Eq. 19 with the limits of integration modified to $[t_{n-1},t_{n}]$. In practice, we use an approximation to the exact $\boldsymbol{\textbf{A}}_{n}$, denoted as $\tilde{\textbf{A}}_{n}$, which involves a truncation of the Magnus expansion (Eq. 19) and an approximation for the time integrals in the truncated Magnus expansion. Using the first $p$ terms in the Magnus expansion results in a time-integration scheme of order $2p$. In this work, we use the second-order Magnus propagator, i.e., obtained by truncating the Magnus expansion after the first term. Additionally, we use a mid-point integration rule to evaluate $\int_{t_{n-1}}^{t_{n}}-i\bar{\boldsymbol{\textbf{H}}}(\tau)d\tau$. As a result, given a set of Kohn-Sham orbitals $\{\bar{\boldsymbol{\psi}}_{1}(t),\bar{\boldsymbol{\psi}}_{2}(t),\ldots,\bar{\boldsymbol{\psi}}_{N_{e}}(t)\}$ defining the density $\rho(\boldsymbol{\textbf{r}},t)$, we can write

$$\bar{\boldsymbol{\psi}}_{\alpha}(\boldsymbol{\textbf{r}},t_{n-1}+\Delta t)\approx e^{\tilde{\textbf{A}}_{n}}\bar{\boldsymbol{\psi}}_{\alpha}(t_{n-1})=e^{-i\bar{\boldsymbol{\textbf{H}}}\left[\rho\left(t_{n-1}+\frac{\Delta t}{2}\right)\right]\Delta t}\bar{\boldsymbol{\psi}}_{\alpha}(t_{n-1})\,,$$

(21)

where $\Delta t=t_{n}-t_{n-1}$ and $\bar{\boldsymbol{\textbf{H}}}\left[\rho\left(t_{n-1}+\frac{\Delta t}{2}\right)\right]$ is the Kohn-Sham Hamiltonian described at the future time instance $t_{n-1}+\Delta t/2$. $\bar{\boldsymbol{\textbf{H}}}\left[\rho\left(t_{n-1}+\frac{\Delta t}{2}\right)\right]$, being dependent on a future instance of the density, is evaluated either by an extrapolation of $\bar{\boldsymbol{\textbf{H}}}$ using $m(>2)$ previous steps or by a second (or higher) order predictor-corrector scheme. In this work, we use the second-order predictor-corrector scheme presented in ¹¹³ (also see ⁹⁶ for details).

The use of the second-order Magnus propagator, as defined in Eq. 21, requires an efficient means of evaluating the action of an exponential operator ($\text{exp}({{\tilde{\textbf{A}}}_{n}})$) on a vector ($\bar{\boldsymbol{\psi}}_{\alpha}(t_{n-1})$). While direct evaluation of $\text{exp}({{\tilde{\textbf{A}}}_{n}})$ is computationally prohibitive beyond small sizes, one can employ subspace projection methods to approximate the action of the exponential operator on a vector. To that end, we use an adaptive Krylov subspace method, wherein we approximate $\text{exp}({{\tilde{\textbf{A}}}_{n}})\bar{\boldsymbol{\psi}}$ as

$$\text{exp}({{\tilde{\textbf{A}}}_{n}})\bar{\boldsymbol{\psi}}\approx\left\lVert\bar{\boldsymbol{\psi}}\right\rVert\boldsymbol{\textbf{Q}}_{k}\text{exp}({\boldsymbol{\textbf{Q}}_{k}^{\dagger}\tilde{\textbf{A}}_{n}\boldsymbol{\textbf{Q}}_{k}})\hat{e}_{1}=\left\lVert\bar{\boldsymbol{\psi}}\right\rVert\boldsymbol{\textbf{Q}}_{k}\text{exp}({\boldsymbol{\textbf{T}}_{k}})\hat{e}_{1}\,,$$

(22)

where $\boldsymbol{\textbf{Q}}_{k}=\{\boldsymbol{q}_{1},\boldsymbol{q}_{2},\ldots,\boldsymbol{q}_{k}\}$ denotes $k$ orthonormal set of vectors (known as Lanczos vectors) that span the Krylov subspace $\mathcal{K}_{k}\{\tilde{\textbf{A}}_{n},\bar{\boldsymbol{\psi}}\}=\{\bar{\boldsymbol{\psi}},\tilde{\textbf{A}}_{n}\bar{\boldsymbol{\psi}},\tilde{\textbf{A}}_{n}^{2}\bar{\boldsymbol{\psi}},\ldots,\tilde{\textbf{A}}^{k-1}\bar{\boldsymbol{\psi}}\}$; $\boldsymbol{\textbf{T}}_{k}=\boldsymbol{\textbf{Q}}_{k}^{\dagger}\tilde{\textbf{A}}_{n}\boldsymbol{\textbf{Q}}$ is a tridiagonal matrix; and $\hat{e}_{1}$ is the first unit vector in $\mathbb{C}^{k}$. Thus, the problem is now reduced to the evaluation of $\text{exp}(\boldsymbol{\textbf{T}}_{k})$, where $\boldsymbol{\textbf{T}}_{k}$ is a small matrix $k\times k$ matrix, and, hence, $\text{exp}(\boldsymbol{\textbf{T}}_{k})$ can be inexpensively computed either through Taylor expansion or exact eigenvalue decomposition of $\boldsymbol{\textbf{T}}_{k}$. A distinct advantage of the Krylov subspace approach is that the error, $\epsilon_{k}$, incurred in the above approximation can be estimated as  ¹¹⁴

$$\epsilon_{k}=\left\lVert\text{exp}({\tilde{\textbf{A}}_{n}})\bar{\boldsymbol{\psi}}-\left\lVert\bar{\boldsymbol{\psi}}\right\rVert\boldsymbol{\textbf{Q}}_{k}\text{exp}({\boldsymbol{\textbf{T}}_{k}})e_{1}\right\rVert\approx\beta_{k+1,k}\left\lVert\bar{\boldsymbol{\psi}}\right\rVert\left|\left[\text{exp}({\boldsymbol{\textbf{T}}_{k}})\right]_{k,1}\right|\,,$$

(23)

where $\beta_{k+1,k}$ is the $(k+1,k)$ entry of $\boldsymbol{\textbf{T}}_{k+1}=\boldsymbol{\textbf{Q}}_{k+1}^{\dagger}\tilde{\textbf{A}}_{n}\boldsymbol{\textbf{Q}}_{k+1}$. As a result, the above relation offers a systematically convergent, efficient and adaptive recipe to evaluate the action of the second-order Magnus propagator on the Kohn-Sham orbitals.

A salient feature of the FE basis is its adaptive spatial resolution, thereby, allowing for finer resolution near the atoms and coarser resolution away from the atoms. The benefits of adaptive resolution is even more pronounced in the case of TDDFT, which warrants the use of large simulation domains to mitigate any spurious artifacts from reflections of the boundary of the simulation domain. However, an efficient and reliable choice of an adaptive FE mesh is useful to be guided by error analysis. To that end, we employ the a priori mesh adaption techniques presented in ⁹⁶ to generate the underlying CFE basis of our EFE basis. The key idea is to minimize the semi-discrete (discrete in space, but continuous in time) error in the dipole moment of the system with respect to the mesh-size distribution ($h(\boldsymbol{\textbf{r}})$). Briefly, the mesh adaption strategy is as follows: (i) for a fixed number of finite elements ($N_{\text{elem}}$), we use the ground-state atomic orbitals and electrostatic potentials to determine the optimal mesh size distribution expression presented in  ⁹⁶ (see Eq. 34 in the reference); (ii) we increase the value of $N_{\text{elem}}$ until chemical accuracy in the quantities of interest (e.g., groundstate energy, excitation energy) is attained. Typically, we conduct the mesh adaption exercise for small representative systems (atoms or small molecules) to determine the characteristic $h(\boldsymbol{\textbf{r}})$ that attains chemical accuracy. Subsequently, for larger systems, we use the $h(\boldsymbol{\textbf{r}})$ from the nearest atom.

Similar to an efficient mesh adaption, we turn to the full-discrete (discrete in both space and time) error analysis presented in ⁹⁶ for an economic choice of the time-step ($\Delta t$). The key aspect of the error analysis is to relate, for a given level of accuracy, the optimal $\Delta t$ to the spatial discretization error. In practice, we first obtain the optimal $\Delta t$ from atomic TDDFT calculations, as per the full-discrete error estimate in ⁹⁶ (see Eq. 35 in the reference). Subsequently, for a multi-atom case we use the least $\Delta t$ obtained for each of its constituent atomic species.

We study the rates of convergence of the dipole moment afforded by the EFE basis with respect to both spatial (FE mesh size $h$) and temporal discretization (time-step $\Delta t$), using a CO molecule of bond length 2.4 a.u. as a benchmark system. In order to reliably study the rates of convergence with respect to the mesh size, we need to mimic a semi-discrete (discrete in space but continuous in time) error analysis, so as to make the spatial discretization error to be the dominant source of error. To that end, we use a small time-step of $\Delta t=10^{-3}$ and small tolerance of $10^{-12}$ for the Krylov subspace projection error (Eq. 23). A large cubical domain of side 50 a.u. is chosen to ensure that the electron density decays to zero on the domain boundary, allowing us to impose Dirichlet boundary condition on the time-dependent Kohn-Sham orbitals and the Hartree potential. We consider two different orders ($p$) of finite elements—HEX27 ($p=2$) and HEX64SPECTRAL ($p=3$). For each $p$, we consider a sequence of increasingly refined meshes using the recipe presented in  2.7 with increasing values of $N_{\text{elem}}$. For all the meshes, we first, obtain the ground-state and then excite the system using a Gaussian electric field of the form $\textbf{E}_{0}(t)=\kappa e^{(t-t_{0})^{2}/2s^{2}}\hat{\boldsymbol{\textbf{x}}}_{1}$, with $\kappa=2\times 10^{-5}$ a.u., $t_{0}=3.0$ a.u., $s=0.2$ a.u.  and $\hat{\boldsymbol{\textbf{x}}}_{1}$ denoting the unit vector along $x$-direction. We note that the semi-discrete error in the dipole moment, as a function of the mesh size ($h$), can be expressed as ⁹⁶

$$\frac{\left|\mu_{x}(t)-\mu_{x}^{h}(t)\right|}{\left|\mu_{x}(t)\right|}=Ch^{q}\,,$$

(24)

where $\mu_{x}(t)$ and $\mu_{x}^{h}(t)$ are the continuum and discrete values of the $x$-component of the dipole moment at time $t$, $C$ is a mesh-independent constant, and $q$ is the rate of convergence. We first evaluate $\mu_{x}(t)$ using the EFE basis with a highly refined HEX125SPECTRAL ($p=4$) finite-element mesh. Subsequently, $C$ and $q$ are obtained by fitting the above relation to a set of $\mu_{x}^{h}(t)$ and $h$. Fig. 1 presents the relative semi-discrete error in the dipole moment at $t=5$ a.u.  (1 a.u. = 0.024188 fs). As evident, the numerical rates of convergence are in close agreement with the theoretical rate of $\mathcal{O}(h^{p})$, where $p$ is the order of the FE basis ($p=2$ for HEX27 and $p=3$ for HEX64SPECTRAL).

We, next, turn to the study the rate of convergence of the dipole moment with respect to temporal discretization. To this end, we used a finite element mesh with a sufficiently refined HEX125SPECTRAL mesh, such that it affords $10^{-4}$ relative error with respect to spatial discretization. We use the same Gaussian electric field as above and propagate the ground-state Kohn-Sham orbitals using second-order Magnus propagator with different $\Delta t$. The full-discrete (discrete both in space and time) error in the dipole moment at time $t_{n}$ is given as ⁹⁶

$$\frac{\left|\mu^{h}_{x}(t_{n})-\mu^{h,n}_{x}\right|}{\left|\mu^{h}_{x}(t_{n})\right|}=C(\Delta t)^{q}\,,$$

(25)

where $\mu^{h}_{x}(t_{n})$ and $\mu^{h,n}_{x}$ are the semi-discrete and full-discrete values of the $x$-component of the dipole moment at time $t_{n}$, $C$ is a constant, and $q$ is the rate of convergence. Fig.  2 depicts the rate of convergence of the dipole moment with respect to $\Delta t$ at $t_{n}=5.0$ a.u. . We attain a numerical rate of convergence of $q=1.99$, which is in good agreement with the quadratic rate of convergence for second-order Magnus propagator.

We now consider methane and benzene molecules as our benchmark systems to compare the performance of the EFE and CFE basis. For both the molecules, we excite the system from their groundstate, in each of the three directions, separately, using a weak Gaussian electric field of the form $\textbf{E}_{0}(t)=\kappa e^{-(t-t_{0})^{2}/2s^{2}}\hat{\boldsymbol{\textbf{x}}}_{j}$, with $\kappa=2\times 10^{-5}$, $t_{0}=3$, $s=0.2$ (all in a.u.), and $\hat{\boldsymbol{\textbf{x}}}_{j}$ being the unit vector along $j-$th coordinate axis. Both the systems are propagated for 20 fs. We use a third-order (HEX64SPECTRAL) and fourth-order (HEX125SPECTRAL) finite element mesh for the EFE and CFE based calculations, respectively, adaptively refined as per the strategy presented in Sec. 2.7, so as to attain 0.1 mHa accuracy in the groundstate energy per atom. Furthermore, we use a large cubical domain of length 60 a.u. and 70 a.u. for methane and benzene, respectively, to ensure that the Kohn-Sham orbitals decay to zero, and thereby, avoid any reflection effects. We chose a time-step ($\Delta t)$ of 0.1 and 0.05 a.u.  for the EFE and CFE calculations, respectively, and a Krylov subspace tolerance ($\epsilon$) of $10^{-7}$ for both the basis sets, to attain $<1$ mHa accuracy in the excitation energies. To evaluate the excitation energies, we first, take the Fourier transform of the dipole moment to obtain the dynamic polarizability, $\alpha_{i,j}(\omega)$, where $i,j$ are the index of the electric field’s polarization direction and measurement direction of the dipole, respectively. Subsequently, we compute the absorption spectrum $S(\omega)=\frac{4\omega}{3\pi}\frac{\text{Tr}\left[\text{Im}\left[\boldsymbol{\alpha}(\omega)\right]\right]}{{\mathcal{F}}\left[\textbf{E}_{0}\right](\omega)}$, where ${\mathcal{F}}\left[\textbf{E}_{0}\right](\omega)$ denotes the Fourier transform of the applied electric field $\textbf{E}_{0}(t)$. Finally, the peaks in the absorption spectrum correspond to the excitation energies. We damp the dipole moment with an exponential window of the form $g(t)=e^{-t/\tau}$, with $\tau=100$ a.u. to artificially broaden the excitation peaks. Fig. 3 and Fig. 4 compare the EFE and CFE basis based absorption spectrum for methane and benzene, respectively. We attain remarkable agreement in the absorption spectrum between the two basis, underscoring the accuracy of the EFE basis. Table 1 compares the performance of the EFE and CFE basis, in terms of the number of basis functions required and the computational cost incurred. As evident, the EFE basis attains a staggering $\sim 50\times$ and $\sim 100\times$ speedup over the CFE basis for the methane and the benzene molecule, respectively. This speedup is attributed to a combination of $3-5\times$ reduction in number of basis functions, $\sim 2\times$ reduction in the stencil of the discrete matrices (HEX64SPECTRAL in EFE and HEX125SPECTRAL in CFE), $2\times$ reduction in the time-step, and $4-5\times$ reduction in Krylov subspace size required by the EFE as compared to that of the CFE basis. Given the huge computational cost associated with the CFE basis, in the remaining benchmark calculations presented in subsequent sections, we only employ the EFE basis.

In order to demonstrate the efficacy of the EFE for larger systems, we study the absorption spectrum of three sodium clusters of increasing sizes—Na₃, Na₁₁, and Na₅₀. For all the three systems, we employ a third-order (HEX64SPECTRAL) finite element mesh that is adaptively refined using the strategy presented in Sec. 2.7 and commensurate with 0.1 mHa accuracy in the groundstate energy per atom. Furthermore, we use large cubical domains—70 a.u. for Na₃ and Na₁₁, and 85 a.u. for Na₅₀—to avoid any reflection effects from the boundary. We use the same Gaussian electric field as used in the methane and benzene benchmark calculations to excite the clusters from their groundstate and propagate the TDKS equations for 20 fs. Similar to the previous examples, we use a time-step of 0.1 a.u. and Krylov subspace tolerance of $10^{-7}$ to attain $<1$ mHa accuracy in excitation energies. Fig. 5 presents the absorption spectrum of the three sodium clusters, normalized with the number of electrons, wherein we have used an exponential window of the form $g(t)=e^{-t/\tau}$, with $\tau=200$ a.u.  to artificially broaden the excitation peaks. Table 2 provides the number of basis functions and the total computation time required by the EFE basis for the three sodium clusters. Using the timings reported, we attain a scaling of $\mathcal{O}(N_{e}^{2.01})$ with respect to the number of electrons, which agrees well with the expected quadratic scaling of the TDKS problem.

The green fluorescent protein (GFP), abundantly found in jellyfish, corals, copepods, is a commonly used tool to understand a wide array of biochemical processes ranging from fluorescence microscopy ¹¹⁵ to gene therapy ¹¹⁶ to reporter gene technology ¹¹⁷. In this example, we study the linear to nonlinear response of the chromophore of the GFP—a small molecule that lends it color. We employ a second-order (HEX27) finite element mesh that is adaptively refined using the strategy presented in Sec. 2.7 and commensurate with 0.1 mHa accuracy in the groundstate energy per atom. We use a large cubucal domain of 120 a.u. to avoid any reflection effects from the boundary. First, we identify the onset of nonlinear response by subjecting the chromophore, at its groundstate, to a series of Gaussian electric field of the form provided in Sec. 3.2 with the strength ($\kappa$) taking values of $2\times 10^{-5}$, $2\times 10^{-4}$, $2\times 10^{-3}$, $2\times 10^{-2}$, and $2\times 10^{-1}$ (all in a.u.). In linear response regime, the dipole moments scaled by the strength of the electric field are expected to be identical, and hence, can be used to identify the onset of nonlinear response. For the purpose of identifying the linear to nonlinear transition, for each of the strengths, we propagate the system for a short duration of 2.5 fs and observe the scaled dipole moment scaled. As evident from  6(a), the clear deviation of the scaled dipole moment for $\kappa=2\times 10^{-1}$ from that of the rest of the strengths, indicates the onset of nonlinear regime to be around $\kappa=2\times 10^{-1}$. Having identified the linear and nonlinear regimes, we propagate the system for longer timescale of 20 fs using a weak ($\kappa=2\times 10^{-4}$) and a strong perturbation ($\kappa=2\times 10^{-1}$). Fig. 6(b) presents a comparison of the absorption spectrum for both the strengths, where the peaks are artificially broadened with an exponential window $e^{-t/\tau}$ with $\tau=300$ a.u. . We observe the first excitation peak at $3.37$ eV, which agrees well with the experimental value of 3.51 eV ¹¹⁸. As evident, the strong perturbation has a decreased height of the first absorption peak, largely owing to saturation effect (i.e., the inability of the highly excited molecule under strong perturbation to absorb subsequent radiations). Furthermore, for the strong perturbation, we observe a blue-shifting of the absorption peaks in the lower energy ($<6$ eV) regime. In the higher energy regime ($>6$ eV), for the strong perturbation, we observe fewer excitation peaks, indicating the absorption to be spread over a range of frequencies. We remark that, despite the use of strong perturbation and large domain size, we require only $\sim 18,000$ basis functions per atom, highlighting the efficiency of the EFE basis even for strong perturbations in TDDFT.

We examine the competence of the EFE basis for simulating strong perturbations by studying the higher harmonic generation in a magnesium dimer with a bond-length of $4.74$ a.u. . We use a strong sinusoidal electric field of the form $\textbf{E}_{0}(t)=\kappa\text{sin}(\pi t/T)\text{sin}(\omega t)\hat{\boldsymbol{\textbf{x}}}_{1}$, with $\kappa=0.01$, $\omega=0.03$, and $T=5\times(2\pi/\omega)$ (all in a.u.) to excite the system from its groundstate and propagate for $T$ fs. For an efficient choice of finite element mesh, we employ adaptive third-order (HEX64SPECTRAL) finite elements, obtained using the approach presented in Sec. 2.7. As with previous examples, we use a cubical domain of length 100 a.u. to eliminate any reflection effects from the boundaries. We employ a time-step of 0.1 a.u. and a Krylov subspace tolerance of $10^{-7}$. Fig. 7 shows the dipole power spectrum ($P(\omega)$) for Mg₂, defined as $P(\omega)=\left|\int_{0}^{T}e^{-i\omega t}\frac{d}{dt^{2}}\mu_{x}(t)dt\right|^{2}$ (i.e., square of the Fourier transform of the acceleration of the dipole moment). We have artificially broadened the peaks by using a gaussian window of form $e^{-t^{2}/\tau}$, with $\tau=10^{5}$ a.u. . For a system with spatial inversion symmetry, only odd multiples of the frequency of the exciting laser pulse must be emitted, which is well verified in Fig. 7 with the peaks in the power spectrum coinciding with odd harmonics (odd multiples of $\omega$). We also observe that the rate of decay of the intensity of the peaks flattens beyond the $13$-$th$ harmonic, which corroborates well with the plateau phenomenon observed in experiments ¹¹⁹. Despite the large domain size used in this calculation, we require only $\sim 32,000$ basis functions per atom, further underscoring the efficiency of the EFE basis even for strong perturbations in TDDFT.

Finally, we examine the parallel scalability of the implementation of our EFE basis for TDDFT calculations, using the Na₅₀ cluster with $~{}\sim 1.7$ million basis functions as a benchmark system. We repeat the same calculation as presented in Sec. 3.3 and measure the relative speedup in the walltime with respect to 72 processors. As evident from Fig. 8, we attain close to ideal scaling until 576 processors, at which we observe a parallel efficiency of 87%. The efficiency drops to 67% at 960 processors owing to fact that, at 960 processors, the number of basis functions per processor falls below 2000, which is low to sustain good parallel scalability.