7.4 Theoretical Background

Detailed description of methods implemented in the riper module is provided in Refs. [204, 205, 207, 208, 214, 215, 179, 177, 206] and references therein. Here, only a short summary of the underlying theory is provided.

7.4.1 Kohn-Sham DFT for Molecular and Periodic Systems

In periodic systems translational symmetry of solids leads to Bloch orbitals \(\psi_{p\sigma}^{\mathbf{k}}\) and one-particle energies \(\varepsilon_{p\sigma}^{\mathbf{k}}\) depending on the band index \(p\), spin \(\sigma\), and the wave vector \(\mathbf{k}\) within the Brillouin zone (BZ), which is the unit cell of reciprocal space. The orbitals

\[\begin{equation} \psi_{p\sigma}^{\mathbf{k}}(\mathbf{r}) = \frac{1}{\sqrt{N_{\text{UC}}}}\sum_{\mathbf{L}} e^{i\mathbf{k}^{\text{T}} \mathbf{L}} \sum_{\mu} C_{\mu p\sigma}^{\mathbf{k}} \mu_{\mathbf{L}}(\mathbf{r}) \end{equation}\](7.5)

are expanded in GTO basis functions \(\mu(\mathbf{r}-\mathbf{R}_{\mu}-\mathbf{L})\equiv \mu_{\mathbf{L}}(\mathbf{r})\) centered at atomic positions \(\mathbf{R}_{\mu}\) in direct lattice cells \(\mathbf{L}\) over all \(N_{\text{UC}}\) unit cells. This results in unrestricted Kohn-Sham equations

\[\begin{equation} \mathbf{F}_{\sigma}^{\mathbf{k}} \mathbf{C}_{\sigma}^{\mathbf{k}} = \mathbf{S}^{\mathbf{k}} \mathbf{C}_{\sigma}^{\mathbf{k}} \mathbf{\varepsilon}_{\sigma}^{\mathbf{k}} \text{,} \end{equation}\](7.6)

which may be solved separately for each \(\mathbf{k}\) in the BZ. The same equations hold for the molecular case, where only \(\mathbf{L} = \mathbf{k} = \mathbf{0}\) is a valid choice and \(N_{\text{UC}}\) is one. Equation (7.6) contains the reciprocal space Kohn-Sham and the overlap matrices \(\mathbf{F}_{\sigma}^{\mathbf{k}}\) and \(\mathbf{S}^{\mathbf{k}}\), respectively, obtained as Fourier transforms of real space matrices

\[\begin{align} F_{\mu\nu\sigma}^{\mathbf{k}} = \sum_{\mathbf{L}} e^{i\mathbf{k}^{\text{T}}\mathbf{L}} F_{\mu\nu\sigma}^{\mathbf{L}} S_{\mu\nu}^{\mathbf{k}} = \sum_{\mathbf{L}} e^{i\mathbf{k}^{\text{T}}\mathbf{L}} S_{\mu\nu}^{\mathbf{L}} \text{.} \end{align}\](7.7)

The elements \(F_{\mu\nu\sigma}^{\mathbf{L}}\) contain three contributions: elements \(T_{\mu\nu}^{\mathbf{L}}\) of the kinetic energy matrix, elements \(J_{\mu\nu}^{\mathbf{L}}\) of the Coulomb matrix, and elements \(X_{\mu\nu\sigma}^{\mathbf{L}}\) of the exchange-correlation matrix,

\[\begin{equation} F_{\mu\nu\sigma}^{\mathbf{L}} = T_{\mu\nu}^{\mathbf{L}} + J_{\mu\nu}^{\mathbf{L}} + X_{\mu\nu\sigma}^{\mathbf{L}}\text{.} \end{equation}\](7.8)

The total energy per unit cell \(E\) is calculated as the sum of the kinetic \(T\), Coulomb \(J\), and exchange-correlation \(E_{\text{XC}}\) contributions,

\[\begin{equation} E = T + J + E_{\text{XC}}\text{.} \end{equation}\](7.9)

7.4.2 RI-CFMM Approach

The key component of riper is a combination of RI approximation and CFMM applied for the electronic Coulomb term [214, 204, 205, 214]. In the RI scheme the total crystal electron density \(\rho^{\mathrm{cryst}}\) is approximated by an auxiliary crystal electron density \(\tilde{\rho}^{\mathrm{cryst}}\)

\[\begin{equation} \rho^{\mathrm{cryst}} \approx \tilde{\rho}^{\mathrm{cryst}} = \sum_{\mathbf{L}}\tilde{\rho}_{\mathbf{L}} \text{,} \end{equation}\](7.10)

composed of unit cell auxiliary densities \(\tilde{\rho}_{\mathbf{L}}\) with

\[\begin{equation} \tilde{\rho}_{\mathbf{L}} = \sum_{\alpha}\mathbf{c}^\text{T}\mathbf{\boldsymbol\alpha}_{\mathbf{L}} \text{,} \end{equation}\](7.11)

where \(\mathbf{\boldsymbol\alpha}_{\mathbf{L}}\) denotes the vector of auxiliary basis functions translated by a direct lattice vector \(\mathbf{L}\). The vector of expansion coefficients \(\mathbf{c}\) is determined by minimizing the Coulomb repulsion \(D\) of the residual density \(\delta\rho = \rho - \tilde{\rho}\)

\[\begin{equation} D = \iint \delta\rho\left(\mathbf{r}\right) \frac{1}{\mid \mathbf{r} - \mathbf{r}^{\prime}\mid} \sum_{\mathbf{L}}\delta\rho_{\mathbf{L}}\left(\mathbf{r}^{\prime}\right) \, d\mathbf{r} \, d\mathbf{r}^{\prime} = \sum_{\mathbf{L}} \left( \delta\rho \mid \delta\rho_{\mathbf{L}} \right) = \sum_{\mathbf{L}} \left( \rho - \tilde{\rho} \mid \rho_{\mathbf{L}} - \tilde{\rho}_{\mathbf{L}} \right) \text{.} \end{equation}\](7.12)

The RI approximation allows to replace four-center electron repulsion integrals (ERIs) by two- and three-center ones. In this formalism, elements \(J_{\mu \nu}^{\mathbf{L}}\) of the Coulomb matrix are defined as

\[\begin{equation} J_{\mu \nu}^{\mathbf{L}} = \sum_{\mathbf{L}^{\prime}} \left(\mu \nu_{\mathbf{L}} \mid \tilde{\rho}_{\mathbf{L}^{\prime}} - \rho_{\textrm{n}\mathbf{L}^{\prime}}\right) \text{,} \end{equation}\](7.13)

where \(\rho_{\textrm{n}}\) denotes the unit cell nuclear charge distribution. The total Coulomb energy including the nuclear contribution is

\[\begin{equation} J = \sum_{\mu\nu\mathbf{L}}D_{\mu\nu}^{\mathbf{L}} J_{\mu \nu}^{\mathbf{L}} - \frac{1}{2} \sum_{\mathbf{L}} \left(\tilde{\rho} + \rho_{\textrm{n}} \mid \tilde{\rho}_{\mathbf{L}} - \rho_{\textrm{n}\mathbf{L}}\right) \text{,} \end{equation}\](7.14)

with the real space density matrix elements obtained by integration

\[\begin{equation} D_{\mu\nu\sigma}^\mathbf{L} = \frac{1}{V_k} \int_\text{BZ} D_{\mu\nu\sigma}^{\mathbf{k}} e^{i \mathbf{k}^\text{T}\mathbf{L}} d\mathbf{k}\text{,} \end{equation}\](7.15)

of the reciprocal space density matrix

\[\begin{equation} D_{\mu\nu\sigma}^{\mathbf{k}} = \sum_{p} f_{p\sigma}^\mathbf{k} \left( C_{\mu p\sigma}^{\mathbf{k}} \right)^\ast C_{\nu p\sigma}^\mathbf{k} \end{equation}\](7.16)

over the BZ with volume \(V_k\).

Equations (7.13) and (7.14) as well as other expressions appearing in the RI scheme require calculation of infinite lattice sums of the form \[\begin{equation} \sum_{\mathbf{L}} \left(\rho_{1} \mid \rho_{2\mathbf{L}} \right), \end{equation}\](7.17) where the distribution \(\rho_{1}\) in the central cell interacts with an infinite number of distributions \(\rho_{2 \mathbf{L}}\), i.e., \(\rho_{2}\) translated by all possible direct lattice vectors \(\mathbf{L}\). In the RI-CFMM scheme[204, 205] the sum in Eq. (7.17) is partitioned into crystal far-field (CFF) and crystal near-field (CNF) parts. The CFF part contains summation over all direct space lattice vectors \(\mathbf{L}\) for which the overlap between the distributions \(\rho_{1}\) and \(\rho_{2 \mathbf{L}}\) is negligible. This part is very efficiently calculated using multipole expansions. The CNF contribution is evaluated using an octree based algorithm. In short, a cubic parent box enclosing all distribution centers of \(\rho_{1}\) and \(\rho_{2}\) is constructed that is large enough to yield a predefined number \(n_{\textrm{targ}}\) of distribution centers per lowest level box. The parent box is successively subdivided in half along all Cartesian axes yielding the octree. In the next step, all charge distributions comprising \(\rho_{1}\) and \(\rho_{2}\) are sorted into boxes based on their extents. Interactions between charges from well-separated boxes are calculated using a hierarchy of multipole expansions. Two boxes are considered well-separated if the distance between their centers is greater than sum of their lengths times \(0.5 \times\)wsicl, where wsicl is a predefined parameter \(\geq 2\). The remaining contribution to the Coulomb term is obtained from direct integration. This approach results in nearly linear scaling of the computational effort with the system size.

7.4.3 k Point Sampling Scheme

The integral in Eq. (7.15) is evaluated approximately using a set of sampling points \(\mathbf{k}\). riper uses a \(\Gamma\)-point centered mesh of \(\mathbf{k}\) points with weights \(w_{\mathbf{k}}\), so that Eq. (7.15) can be written as

\[\begin{equation} D_{\mu\nu\sigma}^\mathbf{L} \approx \sum_\mathbf{k} w_\mathbf{k} e^{i\mathbf{k}^\text{T}\mathbf{L}} D_{\mu\nu\sigma}^{\mathbf{k}} \;\text{.} \end{equation}\](7.18)

In 3D periodic systems each sampling point is defined by its components \(k_1\), \(k_2\) and \(k_3\) along the reciprocal lattice vectors \(\mathbf{b}_1\), \(\mathbf{b}_2\) and \(\mathbf{b}_3\) as

\[\begin{equation} \mathbf{k} = k_1\mathbf{b}_1 + k_2\mathbf{b}_2 + k_3\mathbf{b}_3 \;\text{.} \end{equation}\](7.19)

For 2D periodic systems \(k_3 = 0\). In case of 1D periodicity \(k_3 = 0\) and \(k_2 = 0\). In riper the three components \(k_j\) \((j = 1, 2, 3)\) of \(\mathbf{k}\) are given as

\[\begin{equation} k_j = \frac{i}{n_j} \text{ with } i = -\frac{n_j-1}{2},-\frac{n_j-1}{2} + 1,\dots,\frac{n_j-1}{2}-1,\frac{n_j-1}{2} \;\text{.} \end{equation}\](7.20)

with \(n_j\) \((j = 1, 2, 3)\) as integer numbers. riper reduces the number of \(\mathbf{k}\) points employed in actual calculation by a factor of two using time-inversion symmetry, i.e., the vectors \(\mathbf{k}\) and \(-\mathbf{k}\) are symmetry equivalent. The \(\mathbf{k}\) point mesh can be specified providing the integer values \(n_j\) within the data group $kpoints. Note that two-component Kramers-unrestricted do not use time-inversion symmetry. See Ref. [179] for details.

The number of \(\mathbf{k}\) points required in a calculation critically depends on required accuracy. Generally, metallic systems require considerably more \(\mathbf{k}\) points than insulators to reach the same precision. For metals, the number of \(\mathbf{k}\) point also depends on parameters of the Gaussian smearing [216] used in riper. Please refer to Ref. [216] for more details.

7.4.4 Metals and Semiconductors: Gaussian Smearing

Achieving reasonable accuracy of DFT calculations for metals requires a higher number of \(\mathbf{k}\) points than for semiconductors and insulators. The convergence with respect to the number of \(\mathbf{k}\) points can be improved applying partial occupancies [216]. To achieve this, riper uses the Gaussian smearing method in which occupation numbers \(f_{n\mathbf{k}}\) are calculated as \[\begin{equation} f_{n\mathbf{k}}\left( \frac{\epsilon_{n\mathbf{k}} - \mu}{\sigma} \right) = \frac{1}{2} \left(1 - \textrm{erf}\left[\frac{\epsilon_{n\mathbf{k}} - \mu}{\sigma}\right] \right) \;\text{,} \end{equation}\](7.21) where \(\epsilon_{n\mathbf{k}}\) are band (orbital) energies, \(\mu\) is the Fermi energy and \(\sigma\) is the width of the smearing.

When smearing is applied the total energy \(E\) has to be replaced a generalized free energy \(F\) \[\begin{equation} F = E - \sum_{n\mathbf{k}} \sigma S(f_{n\mathbf{k}}) \end{equation}\](7.22) in order to obtain a variational functional. riper output file reports the values of \(F\) as “FREE ENERGY” and the term \(-\sum_{n\mathbf{k}} \sigma S(f_{n\mathbf{k}})\) as “T*S”. In addition, the value of \(E\) for \(\sigma \rightarrow 0\) is given as “ENERGY (sigma->0)”.

Gaussian smearing can be switched on for riper calculations by simply providing the value of \(\sigma > 0\) within the $riper data group using the keyword sigma, e.g., sigma 0.01. The value of \(\sigma\) should be as large as possible, but small enough to yield negligible value of the “T*S” term. Note, that the value of sigma has to be provided in atomic units. Please refer to Ref. [216] for a more detailed discussion.

The use of Gaussian smearing often requires much higher damping and orbital shifting. Please adjust the values for $scfdamp and $scforbit­al­shift if you encounter SCF convergence problems.

The optional keyword desnue can be used within the $riper data group to constrain the number of unpaired electrons. This can be used to force a certain multiplicity in case of an unrestricted calculation, e.g., desnue 0 for singlet and desnue 1 for doublet.

7.4.5 Low-Memory Iterative Density Fitting Method

For calculations on very large molecular systems a low-memory modification of the RI approximation has been implemented within the riper module[208]. In the RI approximation minimization of the Coulomb repulsion of the residual density, Eq. (7.12), yields a system of linear equations \[\begin{equation} \mathbf{Vc} = \boldsymbol{\gamma}, \end{equation}\](7.23) where \(\mathbf{V}\) is the Coulomb metric matrix with elements \(V_{\alpha \beta} = \left(\alpha \mid \beta \right)\) representing Coulomb interaction between auxiliary basis functions and vector \(\boldsymbol{\gamma}\) is defined as \[\begin{equation} \gamma_{\alpha} = \sum_{\mu\nu} \left(\alpha \mid \mu\nu \right) D_{\mu\nu} . \end{equation}\](7.24) In the LMIDF approach a conjugate gradient (CG) iterative method is used for solution of Eq. (7.23). In order to decrease the number of CG iterations a preconditioning is employed, i.e., Eq. (7.23) is transformed using a preconditioner \(\mathbf{P}\) to an equivalent problem \[\begin{equation} \left(\mathbf{P^{-1}} \mathbf{V}\right) \mathbf{c} = \mathbf{P^{-1}}\boldsymbol{\gamma} \end{equation}\](7.25) with an improved condition number resulting in faster convergence of the CG method. The iterative CG solver in riper employs one of the following preconditioners that are formed from blocks of the \(\mathbf{V}\) matrix corresponding to the strongest and most important interactions between the auxiliary basis functions such that \(\mathbf{P^{-1}V} \approx \mathbf{I}\):

• atomic block preconditioner \[\begin{align*} &P_{\alpha \beta}^{at}=\begin{cases} \left(\alpha^I \mid \beta^I \right), & I \in A_I, A_I \text{ are all atoms in molecule}\\ 0, & \text{otherwise} \end{cases}& \end{align*}\]

• ss block preconditioner: \(P_{\alpha \beta}^{at} \cup P_{\alpha \beta}^{ss}\) \[\begin{align*} &P_{\alpha \beta}^{ss}=\begin{cases} \left(\alpha \mid \beta \right), & \alpha, \beta \in \lbrace S \rbrace, S \text{ are all s auxiliary basis functions}\\ 0, & \text{otherwise} \end{cases}& \end{align*}\]

• sp block preconditioner: \(P_{\alpha \beta}^{at} \cup P_{\alpha \beta}^{sp}\) \[\begin{align*} &P_{\alpha \beta}^{sp}=\begin{cases} \left(\alpha \mid \beta \right), & \alpha, \beta \in \lbrace S, P \rbrace, P \text{ are all p auxiliary basis functions}\\ 0, & \text{otherwise} \end{cases}& \end{align*}\]

The costly matrix-vector products of the \(\mathbf{Vc}\) type that need to be evaluated in each CG iteration are not calculated directly. Instead, the linear scaling CFMM implementation presented above is applied to carry out this multiplication since the elements of the \(\mathbf{Vc}\) vector represent Coulomb interaction between auxiliary basis functions \(\alpha\) and an auxiliary density \(\tilde{\rho}\) \[\begin{equation} \left( \mathbf{Vc} \right)_\alpha = \sum_\beta \left(\alpha \mid \beta \right) c_\beta = \left( \alpha \mid \sum_\beta c_\beta \beta \right) = \left(\alpha \mid \tilde{\rho} \right). \end{equation}\](7.26) Hence, in contrast to conventional RI neither the \(\mathbf{V}\) matrix nor its Cholesky factors need to be stored and thus significant memory savings are achieved.

7.4.6 RT-TDDFT

To investigate the electron dynamics in real time, RT-TDDFT based on Magnus propagator is implemented in riper module[215]. In RT-TDDFT, the time evolution of electron density \(\rho(\mathbf{r},t)\), represented by the single particle reduced density matrix \(\mathbf{D}(t)\) with elements

\[\begin{equation} D_{\mu \nu}(t)=\sum_{m=1}^{N_{\mathrm{MO}}} f_{m} C_{\mu m}^{\dagger}(t) C_{\nu m}(t) \end{equation}\](7.27)

is governed by the von Neumann equation

\[\begin{equation} i\frac{\partial \mathbf{D}(t)}{\partial t}=[\mathbf{F}(t), \mathbf{D}(t)] \end{equation}\](7.28)

where \(\mathrm{F}(t)\) is the time-dependent KS matrix in the orthonormal basis of MO. The von Neumann equation (7.28) is efficiently integrated numerically using the Magnus expansion which evolves the density matrix in time using a unitary operator \(\mathcal{U}(t+\Delta t, t)=\mathrm{e}^{\Omega_{1}+\Omega_{2}+\Omega_{3}+\cdots}\) that conserves the idempotency of \(\mathbf{D}(t)\). Second and fourth order Magnus expansions are implemented.
External perturbation is provided in the form of an electric field \(\mathbf{E}\) which is assumed to be uniform over the whole molecule. The electric field contribution to the KS matrix can be written as

\[\begin{equation} \mathbf{F}_{\mu \nu}^{\mathrm{E}}=-\sum_{j=x, y, z} M_{\mu \nu}^{j} E_{j}, \quad j=x, y, z \end{equation}\](7.29)

with the electric field vector \(\left(E_{x}, E_{y}, E_{z}\right)\) and the dipole moment matrices \(\mathbf{M^j}\)

\[\begin{equation} M_{\mu \nu}^{j}=-\int \mu(\mathbf{r}) j \nu(\mathbf{r}) \mathrm{d} \mathbf{r} \end{equation}\](7.30)

Two time integration methods have been implemented, the self-consistent field (SCF) procedure and the predictor-corrector (PC) scheme.
In the SCF procedure, starting from the ground state electron density \(\mathbf{D}(0)\) and KS \(\mathbf{F}(0)\) matrices, a guess for \(\mathbf{F}(t+\frac{\Delta t}{2})\) is made through linearly extrapolation. Next, \(\mathbf{D}(t)\) is propagated and used to calculate \(\mathbf{F}(t+\Delta t)\), which is followed by a linear interpolation for a better guess for \(\mathbf{F}(t+\frac{\Delta t}{2})\) until convergence is achieved.
In predictor-corrector scheme, \(\mathbf{F}(t+\Delta t/4)\) is predicted by linear extrapolation from previous values and this is used to step \(\mathbf{D}\) forward by \(\Delta t/2\). In contrast to the SCF procedure, which requires multiple KS matrix builds per time step, the PC scheme requires no new KS builds and is comparatively cheaper. However, it is prone to unstability for longer time steps.

Absorption spectra is calculated using the following expression for the dipole strength function

\[\begin{equation} S(\omega)=\frac{1}{3} \cdot \frac{4 \pi \omega}{c} \operatorname{Tr}\left(\operatorname{lm}\left[\alpha_{i j}\right]\right), \quad i, j=x, y, z \end{equation}\](7.31)

where \(\alpha_{i j}\) is the complex polarizability tensor, given as

\[\begin{equation} \alpha_{ij}(\omega)=\frac{\int_{-\infty}^{\infty} \mathrm{e}^{i \omega t} \mu_{j}^{\mathrm{ind}}(t) \mathrm{e}^{-\gamma \mathrm{t}} \mathrm{d} t}{\int_{-\infty}^{\infty} \mathrm{e}^{i \omega t} E_{i}(t) \mathrm{d} t} \end{equation}\](7.32)

where \(\gamma\) is the damping factor (typically in the range of \(0.003-0.005 \mathrm{au}=124-207 \mathrm{ps}^{-1}\)), \(\mu_{j}^{\mathrm{ind}}(t)\) is the time-dependent induced dipole moment and \(E_i\) is the electric field component along \(i\) direction.

RT-TDDFT is also available with hybrid functionals and current-dependent functionals. These methods also require to propagate the imaginary and antisymmetric part of the density matrix, which is related to the current density. For meta-GGAs, the keywords $curswitchengage and $curswitchdisengage enable or disable the current-dependent generalization, respectively. By default, the current-dependent generalization is applied for consistency with escf. The respective exchange-correlation potential is given by the closed-shell restricted Kohn–Sham or open-shell unrestricted Kohn–Sham limit of Refs. [43, 177, 180].

7.4.7 DFT-based Embedding

riper can also be used to perform DFT-based embedding via the riperembed python script. The embedding implementation allows the user to partition the system into active and environment subsystems and the influence of the environment on the active subsystem is accounted for via a DFT-based embedding potential. One can perform either frozen density embedding (FDE) or projection-based embedding (PbE) calculations. Within FDE and PbE, both active and environment subsystems are treated at the DFT level, therefore, these techniques are also referred to as DFT-in-DFT embedding. Since DFT calculations can already be performed extremely efficiently with riper, there is practically no advantage in performing DFT-in-DFT with the riperembed script. However, the key advantage of DFT-based embedding is that it allows to treat the active subsystem using a different and complex method like accurate but computationally expensive wavefunction theory (WFT) method or even RT-TDDFT. The influence of the environment is accounted for via the embedding potential generated by riperembed during the DFT-in-DFT (FDE or PbE) run.

In the following, the various embedding methodologies on which the current DFT-based embedding implementation is based on are described shortly and the readers are referred to Refs. [213, 2] for more details on the implementation.

Frozen Density Embedding

Within the FDE formalism, the total electron density \(\rho^{\text{tot}}\) of the system is partitioned into active and environment subsystem densities, \(\rho^{\text{act}}\) and \(\rho^{\text{env}}\), such that \(\rho^{\text{tot}}=\rho^{\text{act}}+\rho^{\text{env}}\). For a given frozen \(\rho^{\text{env}}\), the \(\rho^{\text{act}}\) is determined by solving the Kohn-Sham constrained electron density (KSCED) equations \[\begin{equation} \left[-\frac{\nabla^{2}}{2}+v_{\mathrm{eff}}^{\mathrm{KS}}\left[\rho^{\mathrm{act}}\right](\boldsymbol{r})+v_{\mathrm{emb}}\left[\rho^{\mathrm{act}}, \rho^{\mathrm{env}}, v_{\mathrm{nuc}}^{\mathrm{env}}\right](\boldsymbol{r})\right] \phi_{i}^{(\mathrm{act})}(\boldsymbol{r})=\epsilon_{i} \phi_{i}^{(\mathrm{act})}(\boldsymbol{r}) ; \quad i=1, \ldots, N^{\mathrm{act}} / 2 \enspace, \end{equation}\](7.33) where \(v_{\mathrm{eff}}^{\mathrm{KS}}\left[\rho^{\mathrm{act}}\right](\boldsymbol{r})\) is the usual KS effective potential of the isolated active subsystem, \(\phi_i^{\text{act}}\) are the KS orbitals of the active subsystem, and \(v_{\text{emb}}\) is the embedding potential describing the effect of the environment subsystem. It is defined as \[\begin{equation} v_{\mathrm{emb}}\left[\rho^{\mathrm{act}}, \rho^{\mathrm{env}}, v_{\mathrm{nuc}}^{\mathrm{env}}\right](\boldsymbol{r})=v_{\mathrm{nuc}}^{\mathrm{env}}(\boldsymbol{r})+\int \frac{\rho^{\mathrm{env}}\left(\boldsymbol{r}^{\prime}\right)}{\left|\boldsymbol{r}-\boldsymbol{r}^{\prime}\right|} d \boldsymbol{r}^{\prime}+\frac{\delta E_{\mathrm{xc}}^{\mathrm{nadd}}\left[\rho^{\mathrm{act}}, \rho^{\mathrm{env}}\right]}{\delta \rho^{\mathrm{act}}(\boldsymbol{r})}+v_{T}\left[\rho^{\mathrm{act}}, \rho^{\mathrm{env}}\right](\boldsymbol{r}) \enspace, \end{equation}\](7.34) where the first and second potentials on the right hand side are due to the nuclei and electrons of the environment, respectively. The third term is the non-additive exchange-correlation potential and the last term is responsible for enforcing the Pauli exclusion principle between the subsystems and is known as the non-additive kinetic potential \[\begin{equation} v_{T}\left[\rho^{\mathrm{act}}, \rho^{\mathrm{env}}\right](\boldsymbol{r})=\frac{\delta T_{s}^{\operatorname{nadd}}\left[\rho^{\mathrm{act}}, \rho^{\mathrm{env}}\right]}{\delta \rho^{\mathrm{act}}(\boldsymbol{r})}=\frac{\delta T_{s}[\rho^{\text{tot}}]}{\delta \rho^{\text{tot}}(\boldsymbol{r})}-\frac{\delta T_{s}[\rho^{\text{act}}]}{\delta \rho^{\text{act}}(\boldsymbol{r})} \enspace, \end{equation}\](7.35) which is evaluated using an approximate KEDF. The embedding potential in eq 7.34 is exact in the limit of the exact KEDF. However, there are certain conditions that the densities are required to fulfill in order to recover the exact energy of the total system: (a) the environment density should be a non-negative function \(\rho^\text{env}(\boldsymbol{r})\geq0\) and (b) the environment density should never be larger than the exact ground-state total density of the system \(\forall \boldsymbol{r} \rho^\text{tot}(\boldsymbol{r}) \geq \rho^\text{env}(\boldsymbol{r})\). Therefore, FDE can give only the upper bound of the ground state energy of the total system if \(\rho^\text{env}\) violates the above conditions. As a result, the quality of the results that one obtains is dependent on the choice of \(\rho^\text{env}\). One possible choice is to use the isolated environment density without the active subsystem. Alternatively, freeze-and-thaw (FaT) cycles may be performed to obtain a relaxed \(\rho^\text{env}\). FaT refers to the procedure where the roles of the active and environment subsystem are interchanged iteratively. Although in principle, the partitioning of the total system is not unique, the FaT procedure yields unique subsystem densities as an artifact of approximate KEDFs. Furthermore, due to the approximate KEDFs, the \(v_\text{emb}\) defined in eq 7.34 is suitable only for weakly overlapping subsystem densities.

Projection-Based Embedding For Molecular Systems

The need for approximate KEDFs within the FDE formalism stems from the fact that the KS orbitals of active and environment subsystems are not orthogonal with respect to each other. Therefore, the total kinetic energy of the system (\(T_s^\text{tot}\)) cannot be written as a simple sum of the kinetic energies of the subsystems (\(T_s^\text{act}+T_s^\text{env}\)) but also requires a non-additive component (\(T_s^\text{nadd}\)). However, this orthogonality can be taken care of by employing a level-shift projection operator (LSPO). The LSPO is defined as \(\mathbf{P_\mathrm{B}}=\mu \mathbf{S}^\mathrm{AB}\mathbf{D}^\mathrm{B}\mathbf{S}^\mathrm{BA}\), where \(\mathbf{S}^\mathrm{AB}\) is the overlap matrix of the active subsystem (A) basis function with the environment (B) basis functions, \(\mathbf{D}^\mathrm{B}\) is the density matrix of the environment subsystem and \(\mu\) is the parameter that tends to infinity ideally, and is taken to be \(10^6\) in practical implementation. Essentially, the LSPO raises the energy of the \(i\)th environment orbital (\(\epsilon^\text{env}_i\)) to a very high value (\(\epsilon^\text{env}_i+\mu\)) for calculations on the active subsystem, thereby ensuring orthogonality between the subsystems and enforcing the Pauli exclusion principle. The embedding potential, in matrix form, can therefore be written as \[\begin{equation} \mathbf{V}_\mathrm{emb} = \mathbf{V}^\mathrm{env}_\mathrm{nuc} + \mathbf{J}^\mathrm{env}_\mathrm{elec} + \mathbf{X}_\mathrm{nadd} + \mathbf{P}_\mathrm{B}\enspace, \end{equation}\](7.36) with the elements \(M_{\mu \nu}\) of the first three matrices on the right hand side (\(\mathbf{M}=\mathbf{V}^\mathrm{env}_\mathrm{nuc}, \enspace \mathbf{J}^\mathrm{env}_\mathrm{elec}, \enspace \mathbf{X}_\mathrm{nadd}\)) defined as \[\begin{equation} M_{\mu\nu}= \left<\mu_{\text{act}}\right|v\left|\nu_{\text{act}}\right>\enspace , \end{equation}\](7.37) where \(v=v_{\mathrm{nuc}}^{\mathrm{env}}(\boldsymbol{r}), \enspace \int \frac{\rho^{\mathrm{env}}\left(\boldsymbol{r}^{\prime}\right)}{\left|\boldsymbol{r}-\boldsymbol{r}^{\prime}\right|} d \boldsymbol{r}^{\prime}, \enspace \frac{\delta E_{\mathrm{xc}}^{\mathrm{nadd}}\left[\rho^{\mathrm{act}}, \rho^{\mathrm{env}}\right]}{\delta \rho^{\mathrm{act}}(\boldsymbol{r})}\) and \(\mu_{\text{act}},\nu_{\text{act}}\) are the basis functions of the active subsystem. PbE combined with FaT cycles can be used to converge to the exact subsystem densities (KS orbitals). It should be noted here that PbE requires a supermolecular basis for the subsystems in order to achieve exact supermolecular DFT results and the approximate monomolecular basis results do not offer any significant improvement over the classic FDE with KEDF-based embedding potential.

Projection-Based Embedding For Periodic Systems

For PbE to work for periodic systems, the only requirement is to make sure that the occupied Bloch orbitals of the active subsystem at a particular \(k\)-point are orthogonal to the occupied Bloch orbitals of subsystem B at the same \(k\)-point. The orthogonality of Bloch orbitals \(\psi_k\) at distinct \(k\)-points is ensured inherently. Therefore, the LSPO can be written as \[\begin{equation} \mathbf{P}^{\mathbf{k}}_{\mathrm{B}}=\mu\left(\mathbf{S}^{\mathbf{k}}_{\mathrm{AB}} \mathbf{D}_{\mathrm{B}}^{\mathbf{k}} \mathbf{S}_{\mathrm{BA}}^{\mathbf{k}}\right), \end{equation}\](7.38) where \(\mathbf{S}^{\mathbf{k}}_{\mathrm{AB}}\) is the overlap matrix of the Bloch functions of A and B, with elements \[\begin{equation} S_{\mu \nu;\mathrm{AB}}^{\mathbf{k}}=\left\langle \mu_{\mathrm{A}}^{\mathbf{k}} \mid \nu_{\mathrm{B}}^{\mathbf{k}}\right\rangle =\sum_{\mathrm{L}} e^{i\mathbf{k}^{\mathrm{T}} \mathbf{L}} S_{\mu \nu;\mathrm{AB}}^{\mathbf{L}}\enspace. \end{equation}\](7.39)

To embed an active periodic subsystem in a periodic environment, the embedding potential is constructed analogously to eq 7.36, where all the terms are replaced by their periodic counterparts. In matrix form, the embedding potential at a particular \(k\)-point can be written as \[\begin{equation} \mathbf{V}^\mathbf{k}_\mathrm{emb} = \mathbf{V}^{\mathbf{k}\mathrm{,env}}_\mathrm{nuc} + \mathbf{J}^{\mathbf{k}\mathrm{,env}}_\mathrm{elec} + \mathbf{X}^\mathbf{k}_\mathrm{nadd} + \mathbf{P}^\mathbf{k}_\mathrm{B}\enspace, \end{equation}\](7.40) where \(\mathbf{V}^{\mathbf{k}\mathrm{,env}}_\mathrm{nuc}\) and \(\mathbf{J}^{\mathbf{k}\mathrm{,env}}_\mathrm{elec}\) are the Coulomb potential matrices due to the nuclei and the electrons of the periodic environment, respectively, and \(\mathbf{X}^\mathbf{k}_\mathrm{nadd}\) is the non-additive exchange-correlation matrix. These are obtained from a Fourier transform of their real-space counterparts \[\begin{equation} M_{\mu \nu }^{\mathbf{k}}=\sum_{\mathbf{L}} \mathrm{e}^{i \mathbf{k}^{\mathrm{T}} \mathbf{L}} M_{\mu \nu }^{\mathbf{L}}\enspace, \end{equation}\](7.41) where \(M_{\mu \nu}\) are the elements of the potential matrices, and \(\mathbf{L}\) is the lattice vector. Similar to the molecular case, a supersystem basis along with FaT is required for exact results. This periodic-in-periodic framework can also be used to perform molecule-in-periodic embedding by using the gamma point embedding potential as an approximation to the real-space embedding potential for the molecule-like active subsystem.

Extension to Correlated Wave Function Theory Methods (WFT-in-DFT)

The energy-error compensation ansatz is employed to perform WFT-in-DFT embedding. Using this approach, the molecular DFT embedding potentials can be used to perform WFT-in-DFT embedding (also known as high level-in-low level embedding in a straightforward and practical manner. This is done by adding the purely DFT-based \(v_{\text{emb}}\) to the Hartree-Fock (HF) core potential of the active subsystem and obtaining the converged HF reference orbitals. Any post-HF method can then be used to obtain an improved description of ground state and excited state properties. After the WFT-in-DFT calculation, the correction to the ground state DFT energy of the total system can be calculated as \[\begin{equation} E_{\text{corr}}=E_{\mathrm{WFT}}^{\mathrm{act}}-E_{\mathrm{DFT}}^{\mathrm{act}} \enspace , \end{equation}\](7.42) where \(E_{\mathrm{WFT}}^{\mathrm{act}}\) and \(E_{\mathrm{DFT}}^{\mathrm{act}}\) are the WFT and DFT ground state energies, respectively, of the embedded active subsystem, obtained self-consistently in the presence of the frozen embedding potential. The corrected energy of the total system is then given as \[\begin{equation} E_{\text {WFT-in-DFT }}^{\mathrm{tot}}=E_{\text {DFT }}^{\mathrm{tot}}+\left(E_{\mathrm{WFT}}^{\mathrm{act}}-E_{\mathrm{DFT}}^{\mathrm{act}}\right) \enspace , \end{equation}\](7.43) where \(E_{\text {DFT }}^{\mathrm{tot}}\) is the low-level DFT energy of the total system. In principle, the doubly counted term \(\int v_{\mathrm{emb}} \rho^{\text {act }}_{\mathrm{WFT/DFT}} d \boldsymbol{r}\), corresponding to the energy of interaction with the environment, should be subtracted from the energies \(E_{\mathrm{WFT/DFT}}^{\mathrm{act}}\). This is because \(E_{\text {DFT }}^{\mathrm{tot}}\) already contains the interaction energy between the subsystems at the DFT level. In this work, the embedding potential is used only as a fixed additional one-electron potential in the Hamiltonian of the active subsystem to obtain \(E_{\text {WFT }}^{\mathrm{act}}\) and \(E_{\text {DFT }}^{\mathrm{act}}\), and the explicit contribution due to the embedding potential \(\int v_{\mathrm{emb}} \rho^{\text {act }}_{\mathrm{WFT/DFT}} d \boldsymbol{r}\) is simply excluded from this quantity. Therefore, the influence of the environment is accounted for only implicitly during the optimization of HF reference orbitals.

The WFT-in-DFT excitation energies can be calculated using either response-based approaches or as the difference of the excited state energies. The latter approach is employed in this work. Therefore, the excitation energy of the embedded active subsystem is given as \[\begin{equation} \Delta E = E_{\mathrm{WFT}}^{\mathrm{act}}\left[\Psi_{\mathrm{e}}^{\mathrm{act}}\right]-E_{\mathrm{WFT}}^{\mathrm{act}}\left[\Psi_{\mathrm{g}}^{\mathrm{act}}\right], \end{equation}\](7.44) where \(E_{\mathrm{WFT}}^{\mathrm{act}}\left[\Psi_{\mathrm{e}}^{\mathrm{act}}\right]\) and \(E_{\mathrm{WFT}}^{\mathrm{act}}\left[\Psi_{\mathrm{g}}^{\mathrm{act}}\right]\) are the energies of the excited and ground state of the active subsystem in the presence of embedding potential. The major approximation in this is that the same ground state embedding potential is used for both ground and excited state energy calculations. Additionally, one more approximation is employed in this work by not including the contribution due to the embedding potential \(\operatorname{tr}\left[\mathbf{D}_{\mathrm{WFT} }^{\mathrm{act}} \mathbf{V}_{\mathrm{emb}}\left[\mathbf{D}_{\mathrm{DFT}}^{\mathrm{act}}, \mathbf{D}_{\mathrm{DFT}}^{\mathrm{env}}\right]\right]\) in the excited and ground state energies of the active subsystem \(E_{\mathrm{WFT}}^{\mathrm{act}}\left[\Psi_{\mathrm{e/g}}^{\mathrm{act}}\right]\). Including the contribution due to \(\mathbf{V}_\mathrm{emb}\) amounts to first-order correction to the energy of the embedded active subsystem’s wavefunction.

Finally, it is worth noting that WFT-in-DFT would not be exact (i.e., reproduce the WFT energies) even in the absence of the approximations introduced above. This is because the DFT and WFT descriptions of the environment as well as the active subsystems are not expected to be the same, irrespective of the kind of embedding potential employed.

Extension To Real Time-Time Dependent Density Functional Theory (RT-TDDFET)

Within RT-TDDFET, the orthonormal basis electron density matrix \(\mathbf{D\,}^\mathrm{act}(t)\), representing the electron density \(\rho^\mathrm{act}(\mathbf{r}, t)\) of the active subsystem embedded in a molecular or periodic environment, is propagated using the LvN-type equation for the active subsystem \[\begin{equation} \mathrm{i} \frac{\partial \mathbf{D\,}^\mathrm{act}(t)}{\partial t}=\mathbf{F\,}^\mathrm{act}_\mathrm{emb}(t) \mathbf{D\,}^\mathrm{act}(t)-\mathbf{D\,}^\mathrm{act}(t) \mathbf{F\,}^\mathrm{act}_\mathrm{emb}(t), \end{equation}\](7.45) where \(\mathbf{F'\,}^\mathrm{act}_\mathrm{emb}(t)\) is the effective KSCED matrix of the embedded active subsystem in the orthonormal basis. In the presence of an external electric field \(\mathbf{E}(t)\), the active subsystem’s time-dependent KSCED matrix \(\mathbf{F\,}^\mathrm{act}_\mathrm{emb}(t)\) in the atomic orbital basis is given as \[\begin{equation} \mathbf{F}^\mathrm{act}_\mathrm{emb}(t)=\mathbf{F}_0^\mathrm{act}(t)+\mathbf{V}_\mathrm{emb}(t)+\mathbf{F}^{\mathrm{E}}(t), \end{equation}\](7.46) where in addition to the contribution \(\mathbf{F}^{\mathrm{E}}(t)\) due to the electric field, the embedding potential matrix \(\mathbf{V}_\mathrm{emb}(t)\) (accounting for the environmental effects) is also added to the KS matrix \(\mathbf{F}_0^\mathrm{act}(t)\) of the active subsystem. The embedding potential can either be KEDF-based (as in FDE) or LSPO-based (as in PbE).