18.1 Background Theory

In the subsystem formulation of the density-functional theory a large system is decomposed into several constituting fragments that are treated individually. This approach offers the advantage of focusing the attention and computational cost on a limited portion of the whole system while including all the remaining enviromental effects through an effective embedding potential. Here we refer in particular to the (fully-variational) Frozen Density Embedding (FDE) [168] with the Kohn-Sham Constrained Electron Density (KSCED) equations [169,170].

In the FDE/KSCED method the embedding potential required by an embedded subsystem with density ρA to account for the presence of another (frozen) subsystem with density ρB is:

         B                δTsnadd[ρA;ρB]  δEnaxcdd[ρA;ρB]
vemb(r) = vext(r)+ vJ[ρB](r)+----δρA(r)----+ ----δρA(r)----

where vextB(r) and vJ[ρB](r) are the electrostatic potentials generated by the nuclei and electron density of the subsystem B, respectively, and

Tnadd[ρ ;ρ  ] = T [ρ + ρ ]- T [ρ  ]- T [ρ ],
 s    A  B     s A    B    s A    s  B

Enxacdd[ρA;ρB] = Exc [ρA + ρB]- Exc[ρA ]- Exc[ρB]

are the non-additive non-interacting kinetic energy and exchange-correlation energy functionals, respectively. In the expressions above Ts[ρ] is the (unknown) non-interacting kinetic energy density functional and Exc[ρ] is the exchange-correlation energy functional. Note that, while the first two terms in Eq. (18.1) refer to classical electrostatics (and could be described by e.g. external point-charges), the last two terms are related to quantum-mechanical effects.

Using freeze-and-thaw [171] cycles, the role of the frozen and the embedded subsystem is iteratively exchanged, till convergence. If expressions (18.2) and (18.3) are computed exactly, then the density ρA + ρB will coincide with the exact density of the total system.

Because the FDE/KSCED was originally developed in the Kohn-Sham framework, using standard GGA approximations for Exc[ρ], the non-additive exchange-correlation potential (δExcnadd∕δρA(r)) can be computed exactly as a functional of the density, leaving the expression of the non-additive kinetic energy term as the only approximation (with respect the corresponding GGA calculation of the total system), because the exact explicit density dependence of Ts from the density is not known. Using GGA approximations for the kinetic energy functional (Ts ˜
TsGGA) we have:

 nadd         ˜GGA           ˜GGA       ˜GGA
Ts   [ρA;ρB] ≈ T s [ρA + ρB]- Ts   [ρA]- Ts   [ρB]


δT-nadd[ρA;ρB]≈ ˜vGGA [ρ  + ρ ](r) -v˜GGA [ρ ](r).
    δρA(r)       T    A    B      T     A

where TGGA(r) = δT˜sGGA∕δρ(r).

The FDE total energy of total system is:

EFDE [˜ρA, ˜ρB] =   Ts[˜ρA]+ Ts[˜ρB]+ Tnsadd[˜ρA;˜ρB ]

             +   Vext[A+ B ]+ J[˜ρA + ˜ρB]+ Exc[˜ρA + ρ˜B ].     (18.6)
Note that this energy differs from the KS total energy of the total system due to the approximation in Eq. (18.4) as well as the approximated kinetic potential (see Eq. 18.5) which lead to approximated embedded densities (˜ρ A ρA and ρ˜ B ρB). With the current state-of-the-art GGA kinetic approximations, the error in the binding energy for weakly interacting systems is close to chemical accuracy.

Using the Generalized Kohn-Sham (GKS) theory, also hybrid exchange-correlation functionals can be used in embedding calculations. To obtain a practical computational method, the obtained embedding potential must be approximated by a local expression as shown in Ref. [172]. This corresponds to performing for each subsystem hybrid calculations including the interaction with other subsystems through an embedding potential derived at a semilocal level of theory. When orbital dependent exchange-correlation functionals (e.g. hybrid functional and LHF) are considered within the FDE method, the embedding potential includes a non-additive exchange-correlation term of the form

Enadd[ρ  ;ρ  ] = E [ΦA+B [ρ + ρ  ]]- E  [ΦA[ρ ]]- E [ΦB [ρ ]]
 xc   A  B     xc       A   B     xc     A     xc    B

where ΦA+B[ρA + ρB] denotes the Slater determinant which yields the total density ρA + ρB. Since such a determinant is not easily available, the non-additive exchange-correlation contribution cannot be determined directly and the non-additive exchange-correlation term can be approximated as [173]

  nadd           GGA            GGA        GGA
E xc [ρA;ρB] ≈ E xc [ρA + ρB ]- E xc [ρA]- E xc [ρB].