In this work, we present a relativistic quantum embedding formalism capable of variationally treating relativistic effects, including scalar-relativity and spin–orbit coupling. We extend density functional theory (DFT)-in-DFT projection-based quantum embedding to a relativistic two-component formalism, where the full spin magnetization vector form is retained throughout the embedding treatment. To benchmark various relativistic embedding schemes, spin–orbit splitting of the nominally *t*_{2g} valence manifold of W(CO)_{6}, exchange coupling of [(H_{3}N)_{4}Cr(OH)_{2}Cr(NH_{3})_{4}]^{4+}, and the dissociation potential curve of WF_{6} are investigated. The relativistic embedding formalism introduced in this work is well suited for efficient modeling of open-shell systems containing late transition metal, lanthanide, and actinide molecular complexes.

## I. INTRODUCTION

As the number of electrons in a system increases, one must make a compromise between cost and accuracy for the electronic structure method used, opening the door for quantum embedding methods. Quantum embedding requires the separation of the total system into a number of subsystems where the subsystems of interest are treated with a more accurate, oftentimes more expensive, method. Developments of quantum embedding in the past 20 years have been numerous and excellent reviews are available.^{1–6}

Embedding is often used as a tool to improve the electron–electron correlation of a subsystem by exploiting the locality of the interaction. Since relativistic effects, such as the scalar relativity and spin–orbit coupling, are more local than electron–electron correlation,^{7–9} embedding should be an effective tool for introducing new physics within the Dirac framework into the Hamiltonian of a quantum subsystem. Relativistic embedding has been previously explored;^{10–15} however, the handling of open-shell fragments in a relativistic framework is yet to be studied. Open-shell environments give rise to a scalar embedding potential and an additional vector potential that corresponds to spin magnetization. A relativistic embedding framework that retains the full magnetization vector form of the spin degrees of freedom can be particularly advantageous for studies of late transition metal, lanthanide, and actinide molecular complexes where relativistic effects play an essential role in their chemical properties and reactivities.

In this work, we seek to introduce a relativistic projection-based embedding scheme within the Dirac formalism, where the full vector form of magnetization density and non-collinear spin are retained throughout the embedding treatment. We considered three tests cases: the spin–orbit splitting of the valence *t*_{2g} orbitals of W(CO)_{6}, the exchange coupling constant of [(H_{3}N)_{4}Cr(OH)_{2}Cr(NH_{3})_{4}]^{4+}, and the potential curve of WF_{6} dissociation. We first consider the splitting of the *t*_{2g} manifold of W(CO)_{6} as a case where only scalar-potential embedding is needed due to the closed-shell nature, but inclusion of spin–orbit coupling is required. The exchange coupling constant of [(H_{3}N)_{4}Cr(OH)_{2}Cr(NH_{3})_{4}]^{4+} tests how well a local treatment of a complex two-spin-center interaction is modeled. We keep the spin centers in the embedded subsystem and use scalar-potential embedding. Finally, we consider the dissociation of WF_{6}, where the spin centers are placed in two different subsystems. Proper handling of relativistic effects and embedding magnetic field are needed for this test. The method introduced here is capable of treating open-shell systems with strong relativistic effects.

## II. METHODOLOGY

### A. Two-component relativistic Dirac-Fock approach

In the spinor basis, the matrix representation of the kinetically balanced modified Dirac equation (in a.u.) is^{7,8}

where *c* is the speed of light. **V**, **T**, and **S** are the two-component nonrelativistic potential energy, kinetic energy, and overlap matrices, respectively. **W** is the relativistic potential matrix, defined as $(\sigma \u2192\u22c5p\u2192)\u2009V(\sigma \u2192\u22c5p\u2192)$, where $p\u2192$ is the linear momentum operator and $\sigma \u2192$ contains Pauli spin matrices. {*ϵ*^{+}} and {*ϵ*^{−}} are the sets of positive/negative eigenvalues with corresponding molecular orbital coefficients $(CL+\u2003CS+)T$ for the positive and $(CL\u2212\u2003CS\u2212)T$ for the negative energy solutions.

In this work, we use the exact-two-component (X2C)^{16–31} technique to “fold” the pseudo-large component coefficients $CS\u2212$ into the large component by a one-step unitary transformation, resulting in an electron-only Hamiltonian (see Ref. 24 for numerical implementations). To lower the computational cost, here we only use the bare Coulomb operator instead of transforming from the four-component two-electron operator, leading to the one-electron X2C approach. In the one-electron X2C framework, the transformation (or “picture change”) is independent of the two-electron operator. This simplification leads to a one-step procedure to construct the transformation matrix through the diagonalization of the one-electron four-component core Hamiltonian. This approximation neglects the “picture-change” of the two-electron and embedding operators. To empirically correct the neglect of two-electron spin–orbit coupling, we make use of the Boettger scaling factor^{32}

where *h*^{SNSO} is the scaling-factor-corrected spin–orbit Hamiltonian (screened-nuclear-spin–orbit Hamiltonian), *i* and *j* are basis functions, *h*^{X2C-SO} is the spin–orbit part of the one-electron X2C Hamiltonian, *Q* is the number of electrons in the filled shell, *l* is the angular momentum, and *Z* is the nuclear charge where the basis function is centered. We work in a general Kramers-unrestricted formalism for representation of the one-electron wave functions.

### B. Embedding relativistic density functional theory

The X2C transformation yields a set of spinor molecular orbitals, which are complex linear combinations of AO spin-orbitals. Spinor molecular orbitals in two-component relativistic electronic structure method can be written as

where the spatial functions ${\varphi k\alpha (r)}$, ${\varphi k\beta (r)}$ are expanded in terms of a common set of real atomic orbital basis functions with complex-valued coefficients, ${c\mu i\alpha}$ and ${c\mu i\beta}$. One-body density matrix ** γ** can be defined as

In one-component nonrelativistic electronic structure methods, due to spin collinearity, only one of the magnetization densities can be non-zero,^{30,33} which is often conveniently chosen to be the *z*-component of the magnetization density vector. In the relativistic regime, the complete vector form of the magnetization density or the current density must be retained in any embedding framework in order to properly simulate open-shell systems and spin-driven reactions, such as intersystem crossing and homogeneous bond dissociation. We will work in the Pauli matrix basis where the spin-blocked Kohn–Sham matrix **F** and density matrix ** γ** are cast as

^{30,34–36}

where

The resulting scalar (**F**_{0}) and spin/magnetic parts (**F**_{1}, **F**_{2}, **F**_{3}) of the Kohn–Sham matrix are

where **h** is the core Hamiltonian, **J** and **K** are the Coulomb and exchange matrices, respectively, *ζ* is the coefficient of the exact exchange in hybrid density functional theory (DFT), and **V**_{xc} is the torque-free non-collinear density functional kernel.^{29,30,37}

In the Pauli matrix basis, ** γ** now have a well-defined physical meaning:

*γ*_{0}is the scalar charge density and

*γ*_{1–3}are magnetization densities that give rise to the magnetic moments. A major advantage of developing embedding schemes in this framework is that each magnetization component can be independently aligned between the high-level and low-level subsystems, and components of vector embedding potentials are also decoupled. In addition, the Pauli matrix basis makes it easier to include external electromagnetic fields related to spectroscopies.

^{34–36}

We partition the quantum system into two quantum subsystems, the cluster and the environment, where it is assumed that the cluster is the region of interest. The charge density and magnetization density vector are partitioned as

leading to a partition of the total energy as

where *E*_{int} is the interaction energy. The interaction energy accounts for the non-additive contribution to the energy upon explicit cluster–environment interactions. The embedding perturbation matrices arising from the interaction energy term in Eq. (12) are

where $Vextenv$ is the external potential due to nuclei in the environment.

In Eq. (13), **V**′ is the matrix representation of the scalar embedding potential

and $(VT\gamma 0tot\u2212VT\gamma 0clu)$ is the matrix representation of the scalar potential due to the non-interacting electronic kinetic energy (**T**_{S}),

In Eq. (14), $B\u2192\u2032$ is the matrix representation of the embedding magnetic field vector,

and $VTm\u2192tot\u2212VTm\u2192clu$ is the matrix representation of the magnetic field due to the non-interacting electronic kinetic energy

The non-additive kinetic energy potential **V**_{T} is difficult to compute analytically since it is challenging to obtain a good analytic model for the kinetic energy as a function of density. See Ref. 4 for more details regarding approximations of **V**_{T}.

The partitioning above is, in principle, exact for a given level of theory. One can use the embedding perturbations in Eqs. (13) and (14), computed at the low level of theory, in a calculation on the cluster at the high level of theory

where $Bn\u2032$ is the *n*th component of the magnetic field vector, *μ*_{B} is the Bohr magneton, and the embedding magnetic field is coupled to the high-level spin degrees of freedom via a spin-Zeeman interaction.^{7,36,38}

### C. Projection-based embedding

We will focus on density functional theory (DFT) based embedding methods in this work.^{39–51} An embedding approach that avoids the non-additive kinetic energy contribution to Eq. (13) is the projection-based scheme.^{6,52} In projection-based embedding, a projection operator is introduced to enforce orthogonality between the orbitals of the cluster and environment. Thus, the non-additive kinetic energy contribution to Eqs. (13) and (14) is zero by construction,^{52} which gives rise to the following embedding perturbations:

The resulting embedding two-component Kohn–Sham matrix in Pauli matrix representation for the high-level calculation is

where

is the renormalized density for the environment subsystem. *μ*_{P} (usually about 10^{6} a.u.) is a parameter for the projection term that elevates the energy of orbitals in the environment to maintain the wave function orthogonality between different subsystems.^{6,52}

The separation of the total system into subsystems in projection-based quantum embedding is done by assigning the molecular orbitals of the total system to different subsystems. Note that schemes involving the choice of nuclei for projection-based embedding have been used to facilitate an automated choice of molecular orbitals.^{52} Projection-based embedding has been shown to be able to handle strongly bonded systems since the relevant bonding orbitals can be embedded as opposed to cutting across a bond.^{53–55}

As pointed out in our previous work,^{38} there is ambiguity in Eq. (19) in the absence of a preferred spin axis (e.g., in an external magnetic field) since the high- and low-level calculations may differ in spin orientation. In addition to the spin-dependent embedding magnetic field, the projection operator also depends upon the spin frame of the low-level calculation [see Eqs. (22) and (23)]. We propose using the low-level cluster density matrix as an initial guess for the embedded high-level cluster calculation. This strategy biases the high-level embedded calculation toward a solution with the same spin frame.

## III. RESULTS AND DISCUSSION

For projection-based embedding, the low- and high-level calculations were performed in a locally modified version of the Chronus Quantum software package.^{56} The triple-*ζ* Sapporo basis set with diffuse functions^{57} was used for W, the 6-311G(d) basis set^{58} was used for C, O, and F in the W-containing complexes, and 6-31G^{59–61} was used for [(H_{3}N)_{4}Cr(OH)_{2}Cr(NH_{3})_{4}]^{4+} in both high- and low-level calculations. A value of *μ*_{P} = 10^{6} a.u. was used for the projection operator.

### A. Relativistic-in-nonrelativistic embedding for closed-shell systems

In the experimental photoelectron spectrum of W(CO)_{6}, two peaks were observed with a 0.29 eV separation and a 2:1 ratio in intensity.^{62,63} These features can be understood at the orbital level,^{64} where the degenerate *b*_{2g} and *e*_{g} orbitals (of the *t*_{2g} manifold) split into *e*′ and $e\u2033$ upon inclusion of spin–orbit coupling. To correctly model the W(CO)_{6} spectrum, electronic structure methods that include spin–orbit effects should be used even though ground-state W(CO)_{6} is a closed-shell system.

When the fine-structure splitting of the high-level system is the subject of interest, the low-level method can be chosen to be nonrelativistic because the splitting is mainly due to the spin–orbit coupling localized at the high-level site. The photoelectron spectrum of W(CO)_{6} is used to benchmark the relativistic-in-nonrelativistic embedding scheme.

The experimental structure of W(CO)_{6} is used in this work.^{65} In projection-based embedding, the partitioning formally occurs with respect to molecular orbital. Since the molecular orbitals for W(CO)_{6} have well-defined symmetries, we did not perform any localization for the selection of projected orbital space. In Table I, we report several different ways that we partition the orbital space of W(CO)_{6}. The embedding schemes were motivated by choosing the environment orbitals with little spatial overlap with the W atom. In embedding scheme **A**, we chose the core orbitals of C and O, scheme **B** adds CO *σ* orbitals, and **C** includes the CO *π* orbitals. See supplementary material for images of orbitals referenced in Table I.

. | . | Environment orbitals . | . |
---|---|---|---|

A | $4a1g,2eg,3t1u,3eg,5t1u,6a1g$ | ||

B | $A+8t1u,5eg,8a1g,9t1u,11t1u$ | ||

C | $B+10t1u,2t2u,1t1g$ |

. | . | Environment orbitals . | . |
---|---|---|---|

A | $4a1g,2eg,3t1u,3eg,5t1u,6a1g$ | ||

B | $A+8t1u,5eg,8a1g,9t1u,11t1u$ | ||

C | $B+10t1u,2t2u,1t1g$ |

Table II presents the splitting of the highest-occupied *t*_{2g}-derived orbitals of W(CO)_{6} as a function of electronic structure method and embedding scheme. In this benchmark study, the nonrelativistic one-component PBE0 exchange-correlation functional^{66} is used in the low level of W(CO)_{6}. Because the splitting is due to spin–orbit coupling, when nonrelativistic PBE0 without spin–orbit is used for the high-level subsystem, it failed to predict the correct fine-structure splitting. In contrast, when the high-level subsystem is modeled with the relativistic X2C-PBE0 approach,^{30} the predicted splitting (0.22 eV) agrees well with the experimental value of 0.29 eV.^{63} We find that the splitting is relatively independent of the partition choice; removing 42% of the electrons from the relativistic X2C-PBE0 calculation yields an error less than 0.01 eV with respect to a calculation on the total system. Because the schemes primarily remove orbitals originating from the CO molecules, this result is reasonable since the spin–orbit is localized to the W atom.

Method . | . | . | X2C-PBE0 . | X2C-PBE0 . | X2C-PBE0 . | . |
---|---|---|---|---|---|---|

scheme . | PBE0 . | X2C-PBE0 . | A
. | B
. | C
. | Expt.^{63}
. |

Δt_{2g} | 0.00 | 0.22 | 0.22 | 0.22 | 0.22 | 0.29 |

$Neclu$ | 158 | 158 | 134 | 110 | 92 |

Method . | . | . | X2C-PBE0 . | X2C-PBE0 . | X2C-PBE0 . | . |
---|---|---|---|---|---|---|

scheme . | PBE0 . | X2C-PBE0 . | A
. | B
. | C
. | Expt.^{63}
. |

Δt_{2g} | 0.00 | 0.22 | 0.22 | 0.22 | 0.22 | 0.29 |

$Neclu$ | 158 | 158 | 134 | 110 | 92 |

### B. Relativistic-in-nonrelativistic embedding for exchange coupling between spin centers

The [(H_{3}N)_{4}Cr(OH)_{2}Cr(NH_{3})_{4}]^{4+} complex, denoted di-Cr(III), is an interesting test case due to the antiferromagnetic superexchange-coupled Cr(III) centers.^{36,67–69} In this second benchmark study, we compute the exchange-coupling constant of di-Cr(III) as

where *E*_{AFM} is the energy of the antiferromagnetic state and *E*_{FM} is the energy of the ferromagnetic state. Note that Eq. (24) assumes an isotropic, nonrelativistic spin Hamiltonian. In this benchmark system, we focus on testing if we can remove a significant amount of the electrons while maintaining the correct spin physics.

We used the computed geometry from Ref. 36 for di-Cr(III) and chose nonrelativistic BH&HLYP^{70} for the low-level calculation and X2C–BH&HLYP as the high-level calculation. For the partitioning, we removed 64 electrons primarily from NH_{3} orbitals (see the supplementary material for more details), leaving the unpaired electrons on each Cr(III) in the high-level calculation. Since the electrons in the environment are paired, the embedding magnetic field is zero, and scalar-potential embedding is used.

Table III lists the exchange-coupling constant of di-Cr(III) from an embedded X2C–BH&HLYP calculation, compared to those from the full X2C calculation. In comparing the X2C–BH&HLYP calculation on the total system with the embedded calculation, we find that we can remove 45% of the electrons from the system and introduce an error of only 0.2 cm^{−1} in *J*. This result suggests that the methodology presented can locally model couplings between multiple spin centers.

### C. Relativistic-in-scalar-relativistic embedding for dissociation curve

In the third benchmark study, we investigate the dissociation energy (Δ*E*) and ground-state potential curve for the WF_{6} → WF_{5} + F reaction. Specifically, we considered singlet WF_{6} in an octahedral geometry dissociating to triplet WF_{5} + F, where the WF_{5} has a trigonal bipyrimidal geometry (see supplementary material for structures). In order to correctly predict the energetics, the low-level system is treated with the spin-free X2C-PBE0 approach in which only scalar-relativistic effects are included, resulting in a relativistic-in-scalar-relativistic embedding scheme.

In this study, equilibrium geometries of WF_{6} and WF_{5} were optimized using two-component PBE0 in the development version of Gaussian.^{71} CRENBL effective core potential^{72} was used for the core electrons and spin–orbit of W; double-*ζ* Sapporo basis set^{57} was used for the 5*s*, 5*p*, 5*d*, and 6*p* electrons of W, and def2-SVP basis set^{73} for F. In order to obtain the dissociation curve, we carried out relaxed potential energy scan in which all molecular degrees of freedom except the W–F bond are optimized.

Canonical orbitals will not be as well localized for the distorted WF_{6} geometries along the potential curve. Thus, we partitioned the system based on localized orbitals. We used Pipek–Mezey localization^{74} due to previous good performance for projection-based embedding.^{6} We consider two embedding schemes for the WF_{6} → WF_{5} + F reaction. For WF_{6}, we projected out six *π* orbitals, and for WF_{5} + F, we projected out the dissociated F *p* orbitals and three *π* orbitals from the high-level subsystem (see the supplementary material images of orbitals). Note that the two schemes do not map to each other since the dissociated F *p*_{z} orbital forms a *σ* bond that has a significant contribution from W basis functions. At the equilibrium bond length, the partitioning scheme is designed to keep the bonding electrons in the same subspace; at the long distance, one of the bonding electrons is projected out to create two product radicals. Since the unpaired electrons are in separate subsystems, both the embedding scalar potential and embedding magnetic field are included.

Table IV compares dissociation energies computed with spin-free PBE0 (SF-PBE0), X2C-PBE0, and embedded X2C-PBE0 in SF-PBE0. The dissociation energies were computed as the difference between the energy of WF_{5} 20 Å away from a fluorine atom and the energy of WF_{6} at equilibrium. We include an experimental value for reference in Table IV, but a direct comparison between computed and experimental values is not recommended due to our neglect of thermal corrections and zero-point energy in the calculations. The addition of spin–orbit coupling lowers the dissociation energy from 5.23 eV to 5.11 eV in going from SF-PBE0 to X2C-PBE0. The embedded X2C-in-SF-PBE0 dissociation energy is in quantitative agreement (0.03 eV difference) with X2C-PBE0, indicating that projection-based embedding is an effective tool for including localized spin–orbit effects.

Figure 1 shows the potential energy curve corresponding to the dissociation of WF_{6} computed using a relaxed potential energy scan. During the singlet region of the potential curve, the bonding electrons are kept in the same subspace. Starting from 3.00 Å, one of the bonding electrons is projected out in order to simulate the homogeneous dissociation path. Inevitably, this creates an inconsistent treatment throughout the potential energy curve scan, but it offers an opportunity to analyze the importance of embedding magnetic field in open-shell systems. We find a good agreement between X2C-in-SF-PBE0 and X2C-PBE0 except in the 2.75 Å–3.25 Å region. This discrepancy is due to the presence of a singlet/triplet crossing and a choice of partitioned MOs (see the supplementary material) that does not continuously map between the singlet and triplet. The former problem can be addressed by moving to an electronic structure method capable of a balanced treatment of ground and excited states, and the latter problem is an active area of recent research.^{76–79}

In order to analyze the importance of the magnetic embedding component, we also report in Fig. 1 the WF_{6} potential curve computed using X2C-in-SF-PBE0 without the embedding magnetic field. Owing to the closed-shell nature in the bonding region, the embedding magnetic field is zero in the range of 1.75 Å–2.75 Å. In the intermediate interaction distance (2.75 Å–3.25 Å), the embedding magnetic field plays an important role of introducing as much as 0.2 eV interaction energy to the system compared to those computed with scalar-potential embedding only. As the interaction distance increases, the magnetic embedding effects decrease quickly and the potential energy curves computed with and without embedding magnetic fields are indistinguishable.

## IV. CONCLUSION

We introduced a relativistic quantum embedding formalism capable of variational treatments of relativistic effects, including scalar relativity and spin–orbit coupling. Projection-based embedding was extended to a relativistic two-component formalism for the first time, where the full spin magnetization vector is explicitly modeled. Two projection-based relativistic embedding protocols, relativistic-in-nonrelativistic and relativistic-in-scalar-relativistic schemes, were benchmarked and applied to studies of molecular complexes. The spin–orbit splitting for the valence *t*_{2g}-derived orbitals of W(CO)_{6} was computed using ground-state relativistic Kohn–Sham density functional theory and the result agrees quantitatively with experimental photoelectron spectroscopy. Using embedding, we can remove almost half of the electrons of the W(CO)_{6} complex from the high-level treatment and introduce an error less than 0.01 eV with respect to a calculation on the total system. We also applied the relativistic-in-nonrelativistic embedding scheme to the exchange coupling constant of [(H_{3}N)_{4}Cr(OH)_{2}Cr(NH_{3})_{4}]^{4+} and found that almost half of the electrons can be removed and only introduce an error of 0.2 cm^{−1}, indicating that local treatments of complex interactions between spin centers can be achieved. The dissociation energy and dissociation potential curve of WF_{6} were studied with the relativistic-in-scalar-relativistic scheme. We find a quantitative agreement of 0.03 eV in the dissociation energy with a relativistic calculation on the total system and a potential curve that is well behaved in all regions except the singlet-triplet crossing. The embedding presented in this work shows promise for treatments of heavy-element derived spin phenomena since it can accurately and efficiently include spin–orbit coupling in systems with multiple spin centers.

## SUPPLEMENTARY MATERIAL

Images of the projected orbitals for W(CO)_{6}, [(H_{3}N)_{4}Cr(OH)_{2}Cr(NH_{3})_{4}]^{4+}, and WF_{6} can be found in the supplementary material along with the geometries of WF_{6} and WF_{5}.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## ACKNOWLEDGMENTS

This work was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Award No. LAB 17-1775, as part of the Computational Chemical Sciences Program. The development of the two-component method is funded by the U.S. Department of Energy (Grant No. DE-SC0006863). The development of the Chronus Quantum open source software package was supported by the National Science Foundation (Grant No. OAC-1663636). Computations were facilitated through the use of advanced computational, storage, and networking infrastructure provided by the Hyak supercomputer system at the University of Washington, funded by the Student Technology Fee.

## REFERENCES

_{2}UO

_{2}Cl

_{4}as a test case

_{2}and Tc

_{2}

_{3})-nZP (n = D, T, Q) sets for the sixth period s-, d-, and p-block atoms

^{*}basis set for atoms K through Zn

^{2}E

_{g}

^{4}A

_{2g}pair states

_{3})

_{4}Cr(OH)

_{2}Cr(NH

_{3})

_{4}]

^{4+}and [(en)

_{2}Cr(OH)

_{2}Cr(en)

_{2}]

^{4+}