Search Results


11.6 Projection-Based Density Embedding

11.6.1 Introduction & Theory

(May 7, 2024)

The formally exact density embedding method developed by Manby, Miller and coworkers 826 Manby F. R. et al.
J. Chem. Theory Comput.
(2012), 8, pp. 2564.
, 748 Lee S. J. R. et al.
Acc. Chem. Res.
(2019), 52, pp. 1359.
affords a step further than electrostatic embedding in the description of chemical environment. This embedding scheme allows for the partition of a system into two interacting subsystems that are treated at two different levels of QM theories, e.g., coupled cluster embedded in DFT. This type of embedding fully accounts for polarization as well as quantum mechanical exchange, as calculated from the super-molecular embedding density and the exchange correlation functional used. The goal of this embedding theory is to perform a higher-level QM (DFT or WFT) calculation in an environment described by a lower-level QM theory. Working Equations

For a system whose total density can be represented as 𝜸=𝜸A+𝜸B, i.e., the molecular orbitals on fragments A and B are orthogonal to each other, its total energy can be expressed as

Etot[𝜸A+𝜸B]=tr(𝜸A𝐡)+J[𝜸A]+Exc[𝜸A]Energy of A+tr(𝜸B𝐡)+J[𝜸B]+Exc[𝜸B]Energy of B+J[𝜸A,𝜸B]+Exc[𝜸A,𝜸B]Non-additive terms (11.99)

where 𝐡 is the core-Hamiltonian matrix, J and Exc are the Coulomb and XC energies (exact exchange K will also be present if a hybrid functional is employed). In an embedding calculation where the electron density of A (denoted as γ~A) is optimized at a high-level theory in the presence of γB (which is obtained at a low-level theory and then fixed), the total energy can be expressed as

Etot[γ~A,γB]=Ehigh[γ~A]+Elow[γ~A+γB]-Elow[γ~A]. (11.100)

Differentiating Eq. (11.100) with respect to γ~A gives the Fock matrix for the embedding calculation:

𝐟A-in-B=δδγ~AEtot[γ~A,γB]=𝐅high[γ~A]+𝐯emb[γ~A,γB], (11.101)

where the embedding potential

𝐯emb[γ~A,γB]=γ~A(Elow[γ~A+γB]-Elow[γ~A])=𝐅low[γ~A+γB]-𝐅low[γ~A]=(𝐉[𝜸A+𝜸B]-𝐉[𝜸A])+(𝐯xclow[𝜸A+𝜸B]-𝐯xclow[𝜸A])+μ𝐏B, (11.102)

where 𝐏B=𝐒𝜸B𝐒 is the projector formed by B’s orbitals and μ is a large enough constant (e.g. 106 a.u.) that effectively elevates the energies of B’s orbitals and thereby enforces the orthogonality between MOs on A and B. In the Q-Chem implementation (with GEN_SCFMAN_EMBED = TRUE), an alternative approach is employed, where one diagonalizes a modified Fock matrix from which the variational degrees of freedom spanned by B’s orbitals are projected out:

𝐟~A-in-B=(𝐈-𝐒𝜸B)𝐟A-in-B(𝐈-𝜸B𝐒), (11.103)

where 𝐒 is the AO overlap matrix.

One disadvantage of using the energy expression given by Eq. (11.100) is that the embedding potential [Eq. (11.102)] depends on γ~A, which needs to be updated in every iteration of the embedding calculation. Thus this scheme affords no savings in computational cost for DFT-in-DFT calculations. To address this limitation, a linearized approximation was suggested, 748 Lee S. J. R. et al.
Acc. Chem. Res.
(2019), 52, pp. 1359.
which is based on the following expansion (to the first order of γ~A-γA):

Elow[γ~A+γB]-Elow[γ~A]Elow[γA+γB]-Elow[γA]+Tr((γ~A-γA)𝐯emb[γA,γB]), (11.104)


𝐯emb[γA,γB]=δδγA(Elow[γA+γB]-Elow[γA])=𝐅low[γA+γB]-𝐅[γA] (11.105)

Based on this linearized approximation, the total energy of the entire system is approximately given by

Etot[γ~A;γA,γB]=Ehigh[γ~A]+Elow[γA+γB]-Elow[γA]+Tr((γ~A-γA)𝐯emb[γA,γB]), (11.106)

and the Fock matrix for embedding calculation:

𝐟A-in-B=δδγ~AEtot[γ~A;γA,γB]=𝐅high[γ~A]+𝐯emb[γA,γB] (11.107)

where the embedding potential stays unchanged during the embedding calculation. In Q-Chem 5.4.1 and versions after, this linearized approximation is used by default in projection-based embedding calculations.

For WFT-in-DFT calculations, one can absorb the embedding potential 𝐯emb[γA,γB] into the Hamiltonian that is used for the correlated WFT calculation, converting Eq. (11.106) into

E[ΨA;γA,γB]=ΨA|H^A-in-B|ΨA+Elow[γA+γB]-Elow[γA]-Tr(γA𝐯emb[γA,γB]). (11.108) Partition of the Occupied Space

An embedding calculation usually starts from an SCF calculation of the full system at the lower level of theory, which yields canonical MOs. Therefore, it is necessary to partition the occupied space and assign orbitals to fragments A and B (without losing generality, assuming A is the embedded fragment). In the original work by Manby et al., 826 Manby F. R. et al.
J. Chem. Theory Comput.
(2012), 8, pp. 2564.
this was achieved by

  • Performing a Pipek-Mezey localization 1019 Pipek J., Mezey P. G.
    J. Chem. Phys.
    (1989), 90, pp. 4916.
    of the canonical occupied orbitals;

  • Assigning a PM-localized orbital to the “active” fragment A if its Mulliken population on A is greater than 0.4.

This approach has been adopted as the default occupied space partition method in Q-Chem.

Recently a parameter-free and more robust partition scheme was proposed by Claudino and Mayhall, which is known as the Subsystem Projected AO Decomposition (SPADE) procedure. 236 Claudino D., Mayhall N. J.
J. Chem. Theory Comput.
(2019), 15, pp. 1053.
In this approach, one first transforms the occupied orbitals into the symmetrically orthogonalized AO basis:

𝐂¯occ=𝐒1/2𝐂occ (11.109)

and then denotes the rows in 𝐂¯occ that correspond to fragment A as 𝐂¯occA. A singular value decomposition (SVD) is then applied to 𝐂¯occA: 𝐂¯occA=𝐔AΣA(𝐕A)T, and the SPADE orbitals are then obtained by rotating the original 𝐂occ:

𝐂occSPADE=𝐂occ𝐕A. (11.110)

The largest gap in the singular value spectrum determines the most appropriate occupied orbital partition under the given fragmentation. Truncation of the Virtual Space Using Concentric Localization

A WFT-in-DFT calculation requires not only the occupied orbitals on the “active” fragment but also the virtual orbitals. Unlike the occupied orbitals, the virtual orbitals obtained from a projection-based embedding calculation are not assigned to fragments but stay delocalized. If the full virtual space is used in the post-SCF calculation, the savings on computational cost will be rather limited since only the number of occupied orbitals is reduced. Therefore, it is desirable to further truncate the virtual space so that one can significantly reduce the computational cost of WFT-in-DFT calculations.

Claudino and Mayhall recently proposed a simple and efficient approach to truncate the virtual space based on concentric localization (CL), 237 Claudino D., Mayhall N. J.
J. Chem. Theory Comput.
(2019), 15, pp. 6085.
which shares the same spirit as the SPADE partition scheme for occupied space. As the first step, the original set of delocalized virtual orbitals (𝐂vir) represented in the working basis (WB) are projected onto the embedded fragment A in a user-specified projection basis (PB):

𝐂virA=(𝐒PB,A-1)(A𝐒PB,WB)𝐂vir (11.111)

where the superscript “A” indicates that only the rows corresponding to fragment A’s basis functions are included and 𝐒PB,A denotes the overlap matrix for PB functions on fragment A only. One can choose PB to be the same as WB or even a smaller basis set. A particular set of virtual orbitals denoted as 𝐂0 can then be selected by performing an SVD on the overlap between 𝐂virA and 𝐂vir:

(A𝐂vir)T(A𝐒PB,WB)𝐂vir=𝐔𝚺𝐕T (11.112)


𝐕 =[𝐕span|𝐕null] (11.113a)
𝐂0 =𝐂vir𝐕span (11.113b)
𝐂null,0 =𝐂vir𝐕null (11.113c)

By construction, 𝐂0 should consist of the virtual valence shell of the WB. In order to achieve higher accuracy for the embedded correlated method, one can select more virtual orbitals from 𝐂null,0 in a stepwise fashion. A recommended way 237 Claudino D., Mayhall N. J.
J. Chem. Theory Comput.
(2019), 15, pp. 6085.
is to singular value decompose the matrix 𝐂T𝐅𝐂null,0, i.e., the coupling between 𝐂0 and 𝐂null,0 through the Fock operator

𝐂0T𝐅𝐂null,0=𝐔0𝚺0𝐕0T (11.114)


𝐕0 =[𝐕span,0|𝐕null,0] (11.115a)
𝐂1 =𝐂null,0𝐕span,0 (11.115b)
𝐂null,1 =𝐂null,0𝐕null,0 (11.115c)
𝐂vir, act =[𝐂0|𝐂1] (11.115d)

As the size of 𝐂1 is the same as 𝐂0, going through the procedure given by Eq. (11.114) doubles the number of active virtual orbitals. This procedure can be carried on iteratively, rendering the accuracy of this method tunable:

𝐂nT𝐅𝐂null,n=𝐔n𝚺n𝐕nT (11.116)


𝐕n =[𝐕span,n|𝐕null,n] (11.117a)
𝐂n+1 =𝐂null,n𝐕span,n (11.117b)
𝐂null,n+1 =𝐂null,n𝐕null,n (11.117c)
𝐂vir, act =[𝐂0|𝐂1||𝐂n+1] (11.117d)

The virtual orbitals that span the null space, 𝐂null,n+1, will remain inactive in the post-SCF calculations. In practice, one is often able to obtain sub-kcal/mol accuracy by only including 𝐂0 and 𝐂1, which is known as the “double-ζ” CL shell model.