11.6.1 Theory

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, the total energy can be represented as

(11.79)

where $𝐡$ is the core-Hamiltonian matrix, $J$ and $E_{\text{xc}}$ are the Coulomb and XC energies. To evaluate the energy of subsystem $A$ in the presence of $B$, one can variationally optimize $A$’s orbitals by diagonalizing the Fock matrix

$${𝐟}^{\text{A in B}}={𝐡}^{\text{A in B}}+𝐉[{𝜸}^A]+{𝐯}_{\text{xc}}[{𝜸}^A]$$

(11.80)

where the modified (embedded) core-Hamiltonian has the following form:

$${𝐡}^{\text{A in B}}=𝐡+(𝐉[{𝜸}^A+{𝜸}^B]-𝐉[{𝜸}^A])+({𝐯}_{\text{xc}}[{𝜸}^A+{𝜸}^B]-{𝐯}_{\text{xc}}[{𝜸}^A])+\mu 𝐏,$$

(11.81)

where $𝐏=𝐒{𝜸}^B𝐒$ is the projector formed by $B$’s orbitals and $\mu$ is a large enough constant (e.g. 10${}^6$ a.u.) that enforces the orthogonality between molecular orbitals on fragments $A$ and $B$. In the Q-Chem implementation, 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:

$${\widetilde{𝐟}}^{\text{A in B}}=(𝐈-𝐒{𝜸}^B){𝐟}^{\text{A in B}}(𝐈-{𝜸}^B𝐒).$$

(11.82)

For a DFT-in-DFT case, the final energy of the full system is

(11.83)

and for a WFT-in-DFT case

(11.84)

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.,⁶¹⁵ this was achieved by

• •

  Performing a Pipek-Mezey localization ⁷⁵⁸ 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.¹⁶⁸ In this approach, one first transforms the occupied orbitals into the symmetrically orthogonalized AO basis:

$${\overline{𝐂}}_{\text{occ}}={𝐒}^{1/2}{𝐂}_{\text{occ}}$$

(11.85)

and then denotes the rows in ${\overline{𝐂}}_{\text{occ}}$ that correspond to fragment $A$ as ${\overline{𝐂}}_{\text{occ}}^A$. A singular value decomposition (SVD) is then applied to ${\overline{𝐂}}_{\text{occ}}^A$: ${\overline{𝐂}}_{\text{occ}}^A={𝐔}^A{\mathrm{\Sigma }}^A{({𝐕}^A)}^T$, and the SPADE orbitals are then obtained by rotating the original ${𝐂}_{\text{occ}}$:

$${𝐂}_{\text{occ}}^{\text{SPADE}}={𝐂}_{\text{occ}}{𝐕}^A.$$

(11.86)

The largest gap in the singular value spectrum determines the most appropriate occupied orbital partition under the given fragmentation.

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 projector-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),¹⁶⁹ which shares the same spirit as the SPADE partition scheme for occupied space. As the first step, the original set of delocalized virtual orbitals (${𝐂}_{\text{vir}}$) represented in the working basis (WB) are projected onto the embedded fragment $A$ in a user-specified projection basis (PB):

$${}^A𝐂_{\text{vir}}^{\prime }=({𝐒}_{\text{PB},A}^{-1}){(}^A{𝐒}_{\text{PB,WB}}){𝐂}_{\text{vir}}$$

(11.87)

where the superscript “$A$” indicates that only the rows corresponding to fragment $A$’s basis functions are included and ${𝐒}_{\text{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 ${}^A𝐂_{\text{vir}}^{\prime }$ and ${𝐂}_{\text{vir}}$:

(11.88)

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 ${𝐂}_{\text{null},0}$ in a stepwise fashion. A recommended way ¹⁶⁹ is to singular value decompose the matrix ${𝐂}^T{\mathrm{𝐅𝐂}}_{\text{null},0}$, i.e., the coupling between ${𝐂}_0$ and ${𝐂}_{\text{null},0}$ through the Fock operator:

(11.89)

As the size of ${𝐂}_1$ is the same as ${𝐂}_0$, going through the procedure given by Eq. (11.89) doubles the number of active virtual orbitals. This procedure can be carried on iteratively, rendering the accuracy of this method tunable:

(11.90)

The virtual orbitals that span the null space, ${𝐂}_{\text{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-$\zeta$” CL shell model.