In Frozen-Density Embedding Theory (FDET), an $N_{\mathit{AB}}$ electron density is represented by means of two independent variables: i) an $N_A$-electron wavefunction (${\mathrm{\Psi }}_A$) and ii) an non-negative real function (${\rho }_B(r)$) such that $\mathrm{\int \; <!--\; nolimits\; -->}<!--\; FUNCTION\; APPLICATION\; -->{\rho }_B(r)dr=N_B$, where $N_B$ is an integer such that $N_{\mathit{AB}}=N_A+N_B$.

The FDET total energy functional ($E_{v_{\mathit{AB}}}^{\mathit{FDET}}[{\mathrm{\Psi }}_A,{\rho }_B]$) satisfies the following relation with the Hohenberg-Kohn energy functional ($E_{v_{\mathit{AB}}}^{\mathit{HK}}[\rho ]$):

$$\begin{array}{rcll}\underset{{\mathrm{\Psi }}_A\to N_A}{\min <!--\; FUNCTION\; APPLICATION\; -->}{E}_{v_{\mathit{AB}}}^{\mathit{FDET}}[{\mathrm{\Psi }}_A,{\rho }_B]=E_{v_{\mathit{AB}}}^{\mathit{FDET}}[{\mathrm{\Psi }}_A^o,{\rho }_B]=E_{v_{\mathit{AB}}}^{\mathit{HK}}[{\rho }_A^o+{\rho }_B],& & & \end{array}$$

where $v_{\mathit{AB}}(r)$ is a given external potential and ${\rho }_A^o(r)=<{\mathrm{\Psi }}_A^o\left|{\sum <!--\; FUNCTION\; APPLICATION\; -->}_{i=1}^{N_A}\delta \left(r_i-r\right)\right|{\mathrm{\Psi }}_A^o>$. The optimal embedded wavefunction is obtained from the FDET eigenvalue equation:

$$\begin{array}{rcll}\left({Ĥ}_{v_{\mathit{AB}}}+{\widehat{v}}_{\mathit{emb}}^{\mathit{FDET}}[{\rho }_A^o,{\rho }_B]\right){\mathrm{\Psi }}_A^o={\lambda }^o{\mathrm{\Psi }}_A^o& & & \end{array}$$

The expressions for $E_{v_{\mathit{AB}}}^{\mathit{FDET}}[{\mathrm{\Psi }}_A,{\rho }_B]$ and the corresponding embedding potential are given in [Wesolowski, Phys. Rev. A 77, 012504 (2008)] for various choices for the form of the embedded wavefunction.

For the embedded non-interacting reference system of electrons, see
[Wesolowski & Warshel, J. Phys. Chem. 97, 8050 (1993)].

For the embedded one-particle density matrix, see
[Pernal & Wesolowski, Int. J. Quant. Chem. 109, 2520 (2009)]

FDET beyond ground-state:

1) For non-interacting reference embedded wavefunctions (LR-TDDFT way):
[Wesolowski, J. Am. Acad. Sci. 126, 11444 (2004)]

2) For interaction embedded wavefunctions (orthogonal stationary-state way):
[Wesolowski, J. Chem. Phys. 140, 18A530 (2014)]

FDET equalities for non-variational embedded wavefunctions:
[ Wesolowski, J. Chem. Theor. & Comput. 16, 6880-6885 (2020) ]

In FDET, the constrained optimisation of the Hohenberg-Kohn energy functional E_vAB^HK[ ρ ] is performed using Ψ_A and ρ_B as independent variables. FDET provides, therefore, the formal basis of various possible multi-scale/multi-level simulation methods combining the quantum mechanical level of description for the embedded species (Ψ_A) with methods applying various physical laws to generate ρ_B.

• FDET with ρ_B and Ψ_A from different quantum mechanical methods
  

  • – ρ_B as the superposition of electron densities of isolated molecules in the environment
    [Wesolowski & Warshel, J. Phys. Chem., 98, 5183 (1994)]; [Zhou et al, J. Am. Chem. Soc., 136, 2723 (2014)] for instance
    

  • – ρ_B as the superposition of electron densities of isolated atoms in the environment
    [Zbiri et al., Chem. Phys. Lett. , 397, 441 (2004)] or [Humbert-Droz et al., Theor. Chem. Acc., 132, 1405 (2013)] for instance

  • – ρ_B as the density of the whole environment
    [Zech et al., J. Chem. Theor. & Comput., 14 4028 (2018)] for instance

  • – ρ_B as the density of the whole environment perturbed by the interaction with the embedded species
    [Zbiri et al., Chem. Phys. Lett. , 397, 441 (2004)] or [Humbert-Droz et al., Theor. Chem. Acc., 132, 1405 (2013)] for instance

• FDET with ρ_B from methods beyond quantum mechanics

  The FDET relations hold for any admissible ρ_B obtained from physical laws appropriate for the characteristic length- and time-scale for the environment.

  • A continuum function (< ρ_B > (r) corresponding to the electron density averaged over configurations of nuclei in a statistical ensemble can be also used as ρ_B in all FDET equations. We have used various techniques (3D-RISM, Molecular DFT, and explicit molecular dynamics simulations) to generate < ρ_B > and use it to estimate the solvent effect on absorption- (see [Kaminski et al., J. Phys. Chem., 114, 6082 (2010); Zhou et al., Phys. Chem. Chem. Phys., 13 10565 (2011) ; Laktionov et al, Phys. Chem. Chem. Phys., 18 21069 (2016) ; Gonzalez-Espinoza et al., J. Chem. Theor. & Comput., 18 1072 (2022)] ) or emission [Shedge & Wesolowski, ChemPhysChem, 15, 3291 (2014)] spectra.

  • Recently, we have shown that the experimental electron density obtained from X-ray diffraction data on molecular crystals used as ρ_B in FDET equations yields excitation energies of embedded chromophores comparable to the ones obtained using ρ_B from ab initio calculations [Ricardi et al. Acta Crystallographica A - Foundations and Advances, 76, 571 (2020)].

FDET vs Orbital-Free DFT:

The density functional for the kinetic energy (T_s[ρ]) is a key quantity in both, the embedding methods based on FDET and in methods based on the orbital-free formulation of DFT (OF-DFT). (For an overview of OF-DFT methods, see Recent Progress in Orbital-Free Density-Functional Theory, T.A. Wesolowski & Y.A. Wang, Eds., World Scientific, 2013). Any approximation used for T_s[ρ] and its functional derivative applicable in OF-DFT methods can be also used to approximate the non-additive kinetic energy:

T_s^nad[ρ_A,ρ_B]= T_s[ρ_A+ρ_B]-T_s[ρ_A]-T_s[ρ_B]

and its functional derivative with respect to ρ_A (non-additive kinetic otential), which is one of the components of the FDET embedding potential. Such a "top-down" strategy to approximate this derivative proved to be prone to surprising failures in FDET. For example, although addition of the second-order expansion correction to Thomas-Fermi functional is known to improve the kinetic energy, the addition of the corresponsing correction to the non-additive kinetic energy and potential worsens the FDET results as shown already in the our original work reporting the first FDET based embedding method [Wesolowski & Warshel, J. Phys. Chem. 97, 8050 (1993)]. Subsequent dedicated studies, have shown that this failure is due to the incorrect behaviour of the derivative of T_s^nad[ρ_A,ρ_B] rather than the T_s^nad[ρ_A,ρ_B] contribution to the energy [Wesolowski & Weber, Intl. J. Quant. Chem. 61, 303 (1997)]. We pursued two strategies in developeding approximations for the non-additive kinetic potential and energy. The decomposable one, in which the second order contributions are smoothly dupmed (GGA97 approximation - see below) or the non-decomposable strategy ("bottom-up") in which the exact mathematical properties are imposed on the non-additive kinetic potential without constructing the "parent" approximation for T_s[ρ] (NDCS abd NDSD functionals - see below).

Approximations to the bi-functional of the non-additive kinetic energy and/or potential:

- Generalized Gradient Approximation for the non-additive kinetic potential (GGA97): [Wesolowski et al., J. Chem. Phys. 105, 9182 (1996)] and [Wesolowski, J. Chem. Phys. 106, 8516 (1997) ],
- non-decomposable approximation for the non-additive kinetic potential (NDSD): [Garcia Lastra et al., J. Chem. Phys. 129, 074107 (2008)],
- improved non-decomposable approximation for the non-additive kinetic potential (NDCS): [Polak et al., J. Chem. Phys. 156, 044103 (2022)]

I. The origins:

Subsystem formulation of DFT traces its origin to the methods used in the solid-state physics community to model ionic and rare gas crystals [Senatore & Subbaswamy, Phys. Rev. B 34, 5754 (1986)]. In 1991, Cortona gave such methods a formal foundations by relating the underlying relations for energy and density to the Hohenberg-Kohn theorems [Cortona, Phys. Rev. B 44, 8454 (1991)] and applied a method based on subsystem DFT to model structure and elastic properties of ionic solids.

II. Our contributions:

In 1996, we started exploring the applicability of methods based on subsystem DFT in the field of intermolecular interactions [Wesolowski & Weber, Chem. Phys. Lett. 248, 71 (1996)]. This lead us to the development of:
Efficient algorithms:
- "freeze-and-thaw" optimisation of subsystem densities: [Wesolowski & Weber, Chem. Phys. Lett. 248, 71 (1996)]
- simultaneous optimisation of density and nuclear coordinates:[Dulak et al., J. Chem. Theor. & Comput. 3, 735 (2007)];
- linearisation of the non-electrostatic components of the FDET embedding potential [Dulak & Wesolowski, J. Chem. Theory. & Comput. 2 1538 (2006)],
and
Approximations to the bi-functional of the non-additive kinetic energy and/or potential:
- Generalized Gradient Approximation for the non-additive kinetic potential (GGA97): [Wesolowski et al., J. Chem. Phys. 105, 9182 (1996)] and [Wesolowski, J. Chem. Phys. 106, 8516 (1997) ],
- non-decomposable approximation for the non-additive kinetic potential (NDSD): [Garcia Lastra et al., J. Chem. Phys. 129, 074107 (2008)],
- improved non-decomposable approximation for the non-additive kinetic potential (NDCS): [Polak et al., J. Chem. Phys. 156, 044103 (2022)]

In 1996-2007 we applied subsytem DFT based mathod to various types of intermolecular complexes which allowed us to determine the domain of applicability of the semi-local approximations for the density functionals approximated in subsystem DFT . For a benchmarking studies, see [Dulak & Wesolowski, J. Molecular Modeling 13, 631 (2007)] for interaction energies, and [Dulak et al., J. Chem. Theor. & Comput. 3, 735 (2007)] for equilibrium complexes.

III. Extension of subsystem DFT to excited states:

In 2004, in collaboration with Mark E. Casida, we generalized subsystem formulation of DFT as a ground state formalism to excited states using the general framework of linear-response time-dependent DFT. [M.E. Casida & T.A. Wesolowski, Intl. J. Quant. Chem. 96 577 (2004)]

IV. Subsystem DFT: current interests

The overview of such developments till 2006 is given in the review [Tomasz A. Wesolowski, One-electron Equations for Embedded Electron Density: Challenge for Theory and Practical Payoffs in Multi-Level Modeling of Complex Polyatomic Systems, in: Computational Chemistry: Reviews of Current Trends - Vol. 10, World Scientific, 1-82 (2006)]. Since the publication of the benchmarking results in 2006 -2007, we are using subsystem DFT based methods only occasionally: either as a tool to determine the accuracy of the developed approximations for the non-additive kinetic energy bi-functional (and its derivative) [Polak et al., J. Chem. Phys. 156, 044103 (2022)] or as one of many possible techniques to generate the frozen density in Frozen-Density Embedding Theory based multi-level/multi-scale simulations.