Local Energy Decomposition of Open-Shell Molecular Systems in the Domain-Based Local Pair Natural Orbital Coupled Cluster Framework

Local energy decomposition (LED) analysis decomposes the interaction energy between two fragments calculated at the domain-based local pair natural orbital CCSD(T) (DLPNO-CCSD(T)) level of theory into a series of chemically meaningful contributions and has found widespread applications in the study of noncovalent interactions. Herein, an extension of this scheme that allows for the analysis of interaction energies of open-shell molecular systems calculated at the UHF-DLPNO-CCSD(T) level is presented. The new scheme is illustrated through applications to the CH2···X (X = He, Ne, Ar, Kr, and water) and heme···CO interactions in the low-lying singlet and triplet spin states. The results are used to discuss the mechanism that governs the change in the singlet–triplet energy gap of methylene and heme upon adduct formation.


INTRODUCTION
Weak intermolecular interactions play a major role in virtually all areas of chemical research. 1−6 Perturbative and supermolecular approaches can be used to evaluate weak interaction energies between two or more fragments. Within a perturbative approach the Hamiltonian is partitioned into contributions of noninteracting fragment terms plus a series of perturbing potentials representing the interaction between the fragments. The most popular approach of this type is the symmetryadapted perturbation theory (SAPT), which also provides a decomposition of the interaction energy into a series of physically meaningful contributions, including electrostatics, induction, London dispersion, and exchange-repulsion terms. 7 Although these interaction energy components have no unique definition, their quantification has been instrumental for the rationalization of the underlying mechanism that gives rise to the interaction. 8 Within a supermolecular approach, the interaction energy is computed as the difference between the total energy of the adduct and the energy of the separated monomers. The decomposition of the resulting interaction energy into physical terms is then obtained via energy decomposition analysis (EDA) schemes, which are mainly based on a seminal work of Morokuma. 9 In these schemes the interaction energy is typically partitioned into electrostatics, charge transfer and/ or polarization, and Pauli repulsion (also called exchangerepulsion) terms. 10−14 The large majority of EDA schemes applied in mainstream computational chemistry relies on density functional theory (DFT) for the calculation of interaction energies. Hence, the range of applicability of such schemes is limited by the accuracy of the chosen functional. This is especially limiting in the context of noncovalent interactions, as DFT does not properly describe London dispersion. Although several strategies have been suggested 15−19 in order to practically deal with this shortcoming, a quantitative understanding of weak interactions typically requires the use of highly correlated wave function-based ab initio methods in conjunction with large basis sets. 8,20−23 In particular, the coupled-cluster method with single, double, and perturbative treatment of triple excitations, i.e., CCSD(T), 22,24 has proven its reliability in a wide range of contexts and typically allows one to compute relative energies with chemical accuracy (defined here as 1 kcal/mol). Unfortunately, the computational cost of CCSD(T) increases as the seventh power of the molecular size. Hence, the calculation of CCSD(T) energies is only possible for benchmark studies involving relatively small systems. To overcome this limitation, the domain-based local pair natural orbital CCSD(T) method, i.e., DLPNO-CCSD(T), was developed. 25−34 This technique typically retains the accuracy and reliability of CCSD(T), as shown on many benchmark data sets, 34−37 while allowing at the same time for the calculation of single-point energies for systems with thousands of basis functions.
In order to aid in the interpretation of DLPNO-CCSD(T) results, we recently introduced an EDA approach called local

Journal of Chemical Theory and Computation
Article and heme upon the formation of the adducts are discussed in terms of the corresponding LED interaction energy components for the low-lying singlet and triplet states. The dependence of the results on the various technical aspects of the calculations is discussed in the Supporting Information (Tables S2−S13). The distance dependence of the open-shell LED terms is discussed in section 4 for the PFe-( 1 A 1g / 3 A 2g / 3 E g )···CO interaction and given in Figure S1 for the CH 2 ···H 2 O interaction. The last section of the paper is devoted to discussion of the results and concluding remarks of the study.
where i and j are the occupied orbitals of the reference determinant, n i,α and n i,β denote the occupation numbers (1 or 0) for α and β electrons, respectively, n i = n i,α + n i,β ; (ii|jj) and (ij|ij) are two electrons integrals in Mulliken notation, and Z A (Z B ) is the nuclear charge of atom A (B) at position R A (R B ). The DLPNO-CCSD(T) correlation energy can be written essentially as a sum of electron-pair correlation energy (ε ij , where i and j denote the localized orbitals) contributions plus the perturbative triples correction (E C−(T) ). Local second-order many-body perturbation theory is used to divide the ε ij terms into "weak pairs", with an expected negligible contribution to the correlation energy, and "strong pairs". The contribution coming from the weak pairs is kept at the second-order level, whereas the strong pairs are treated at the coupled cluster level. Hence, the overall correlation energy reads where E C−SP denotes the energy contribution from the strong pairs, E C−WP is the correlation contribution from the weak pairs, and the last term is associated with the perturbative triples correction E C−(T) .
The dominant strong pair contribution reads where the indices with an overbar denote β spin−orbitals and those without overbar are used for α spin−orbitals. The first two terms represent the contribution from the single excitations (a i , the singles PNOs; t a i i , the singles amplitudes).
These terms vanish if the Brillouin's theorem is satisfied. 73 This is not generally the case when the quasi-restricted orbitals (QROs) 75 or restricted open-shell HF (ROHF) determinant is used, the former of which is the standard choice in open-shell DLPNO-CCSD(T) calculations. 73 The ε ij , ε ̅ ̅ i j and ε ̅ ij terms denote α−α, β−β, and α−β pair correlation energies, defined as where a ij and b ij are the PNOs belonging to the ij pair, (ia ij | jb ij ) are the two-electrons integrals, are the cluster amplitudes in the PNO basis, and t a ij b ij ij , t a ij i , t b ij i are the doubles and singles amplitudes of the coupled cluster equations. 2.2. Open-Shell Variant of the Local Energy Decomposition. Within a supermolecular approach, the energy of a molecular adduct XY relative to the total energies of noninteracting fragments X and Y, i.e, the binding energy of the fragments (ΔE), can be written as where ΔE geo-prep is the geometric preparation energy needed to distort the fragments X and Y from their structures at infinite separation to their in-adduct geometry. ΔE int is the interaction Scheme 1. Simple Schematic Representation of the Singlet and Low-Lying Triplet Electronic Configurations of PFe and PFe··· CO a a One of the degenerate pair of 3 E g configurations is shown. The other configuration is clearly (d

Journal of Chemical Theory and Computation
Article energy between the fragments X and Y frozen in the geometry they have in the adduct XY, which is defined as In the DLPNO-CCSD(T) framework, the orbitals of the reference determinant are initially localized. Hence, each orbital can be usually assigned to the fragment where it is dominantly localized. By exploiting the localization of the occupied orbitals, it is possible to regroup the terms of the reference energy of the XY adduct E ref XY (eq 1) into intra-and intermolecular contributions in exactly the same way as it was described for the closed-shell case 38,40 The intermolecular reference energy can be further partitioned into electrostatic (E elstat ) and exchange (E exch ) interactions. Accordingly, the ΔE int ref term can be written as The electronic preparation energy ΔE el-prep ref is positive and thus repulsive. It corresponds to the energy needed to bring the electronic structures of the isolated fragments into the one that is optimal for the interaction. E elstat and E exch are the electrostatic and exchange interactions between the interacting fragments, respectively (note that the "ref" superscript in these terms is omitted for the sake of simplicity from now on). It is worth noting here that the intermolecular exchange describes a stabilizing component of the interaction, lowering the repulsion between electrons of the same spin. Note that E elstat incorporates the Coulomb interaction between the distorted electronic clouds of the fragments. Hence, it accounts for both induced and permanent electrostatics. We recently proposed an approach for disentangling these two contributions in E elstat in the closed-shell case, 76 which is in principle also applicable to the open-shell case.
The correlation contribution to the interaction energy E int C can thus be expressed as a sum of three contributions where ΔE int C-SP , ΔE int C-WP , and ΔE int C-(T) are the strong pairs, weak pairs, and triples correction components of the correlation contribution to the interaction energy, respectively. The ΔE int C-SP , ΔE int C-WP , and ΔE int C-(T) terms can be further divided into electronic preparation and interfragment interaction energies based on the localization of the occupied orbitals, as already described previously. 31 For the (typically) dominant ΔE int C-SP contribution, a more sophisticated approach is used in order to provide a clear-cut definition of London dispersion. The E C-SP XY term can be written as a sum of the contributions of single and double excitations according to eqs 3−6. By exploiting the localization of both the occupied and the virtual orbitals in the DLPNO-CCSD(T) framework, E C-SP XY can be divided into five different contributions where each energy term contains contributions from the αα, ββ, and αβ excitations constituting the correlation energy, as given in eqs 4−6. The corresponding energy expressions of these terms are reported in the following. Note that the subscript ij in the PNOs is omitted for clarity and that the indices X and Y represent the fragment on which the orbital is dominantly localized. For the sake of simplicity, only contributions from the αα pairs (and singles excitations of α spin−orbitals) are shown as an example The E C-SP (Y) and E C-SP CT(X←Y) terms can be obtained from eqs 14 and 15 by exchanging the X and Y labels in all terms. Note that the last term in E C-SP (X) denotes the contribution from the singles excitations, whose physical meaning will be discussed later in this section. The relevant pair excitation contributions constituting these terms are shown pictorially in Figure 2.
E C-SP (X) and E C-SP (Y) describe the correlation contribution from excitations occurring within the same fragment (also called intrafragment correlation). The dynamic charge transfer (CT) terms E C-SP CT(X→Y) and E C-SP CT(X←Y) represent the correlation energy contributions that arise from instantaneous cation−anion pair formations. In contrast to the Coulombic interaction energy decaying with r −1 , they occur with a small probability decaying exponentially with distance. These terms are essential for correcting overpolarized electron densities at the reference level. E C-SP DISP(X,Y) describes the energy contribution associated with genuine and exchange dispersion (see Figure 2).
For the sake of simplicity, it may be useful to combine several terms. For example, ΔE el-prep C-SP (i.e., the difference between the intrafragment contributions E C-SP (X/Y) and the strong

Journal of Chemical Theory and Computation
Article pair correlation energy of the isolated fragments) is always positive, while E C-SP CT(X→Y) and E C-SP CT(X←Y) are negative. These two terms typically compensate for each other to a large extent. 38,40 Hence, they can be combined to give the SP correlation contribution to the interaction energy excluding dispersion contribution (ΔE no-disp C-SP ) where ΔE no-disp C-SP describes the correlation correction for the interaction terms approximately accounted for at the reference level (e.g., permanent and induced electrostatics).
Note that the intermolecular part of the weak-pair contribution has mainly dispersive character as it describes the correlation energy of very distant pairs of electrons. Hence, it can be added to the strong pair dispersion term E DISP C-SP in order to obtain the total dispersion contribution at the DLPNO-CCSD level. We label this summation as E disp C-CCSD . The remaining part of the correlation interaction energy at the DLPNO-CCSD level is labeled as ΔE no-disp C-CCSD . Therefore Collecting all the terms we obtain for the binding energy Equation 19 is the base of our LED scheme and applies identically to the closed-shell and to the open-shell cases. As a final remark, it is worth noting that the magnitude of the single excitations term of the dynamic electronic preparation energy is typically very close to zero for closedshell fragments due to the Brillouin's theorem. Hence, this term provides a useful tool for localizing the fragments in which unpaired electrons are present (see Tables S5, S6, and S12).
2.3. LED for the Analysis of Singlet−Triplet Energy Gap. As an illustrative example of the usefulness and applicability of LED scheme, we discuss in this work the influence of the different LED components on the singlet− triplet energy gap (E S−T ) of a molecule Y upon formation of a weak intermolecular interaction with X. The E S−T gap of the corresponding 1,3 Y···X adduct can be written as The same quantity for the free Y reads as The variation in the singlet−triplet gap of Y upon interaction with a molecule X (Δ) can be obtained by subtracting eqs 20 and 21 Thus, Δ equals the difference between the ΔE( 1 Y···X) and ΔE( 3 Y···X) binding energies. In the following, the LED scheme is used to decompose ΔE( 1 Y···X) and ΔE( 1 Y···X). This approach provides insights into the physical mechanism responsible for Δ in different systems.

COMPUTATIONAL DETAILS
All calculations were carried out with a development version of ORCA software package. 43,44 The presently described developments were made accessible in the ORCA 4.1 release, which is free of charge to the scientific community. 3.1. Geometries. The CH 2 molecule and its adducts (see Figure 1) were fully optimized numerically at the DLPNO-CCSD(T) level 25−34 by employing the aug-cc-pVTZ basis set and its matching auxiliary counterparts. 77−80 We have shown on the CH 2 ···He adduct that DLPNO-CCSD(T) gives almost identical geometries as its parent method, i.e., the canonical CCSD(T), with nonbonded interatomic distances having a maximum deviation of 0.007 Å (see Table S1 and the associated coordinates in section S1.1). The fully optimized structures of the adducts of both 1 CH 2 and 3 CH 2 at the DLPNO-CCSD(T)/aug-cc-pVTZ level are shown in Figure 1 (see section S1.1 for their coordinates).
The ferrous heme species and their CO adducts (see Figure  1) were fully optimized using the RIJCOSX approximation 81,82 at the B3LYP 83−85 level incorporating the atom-pairwise dispersion correction (D3) with Becke−Johnson damping (BJ), i.e., B3LYP-D3. 86 The def2-TZVP basis set was used in conjunction with its matching auxiliary JK counterpart. 87,88 Relaxed potential energy surface (PES) scans were performed at the B3LYP-D3/def2-TZVP level for the PFe( 1 A 1g )···CO, PFe( 3 A 2g )···CO, and PFe( 3 E g )···CO interactions. The optimized coordinates of the 1 A 1g , 3 E g , and 3 A 2g states of PFe and PFe···CO (see Figure 1) at the B3LYP-D3/def2-TZVP level are given in section S2.1. It should be noted here that in the optimized structures the CO moiety is almost perpendicular to the heme plane on both singlet and triplet surfaces, consistent with previous experimental and computational studies. 89−92 Figure 2. Schematic representation of strong pair excitations from occupied orbitals to virtual orbitals (PNOs) in the framework of the DLPNO-CCSD(T)/LED method. Intramolecular excitations occur within the same fragment. For the sake of simplicity, they are shown only for the excitations on X. Dynamic electronic preparation energy is obtained by subtracting the strong pairs energy of the isolated fragments X and Y (frozen in their in-adduct geometries) from the corresponding intramolecular terms. Only the charge transfer excitations from X to Y are shown. Analogous charge transfer excitations also exist from Y to X.  Tables  S3 and S4). The distance dependence of the open-shell LED terms is given in Figure S1 for CH 2 ···H 2 O and in section 4 for PFe···CO.

Journal of Chemical Theory and Computation
Unless otherwise specified, DLPNO-CCSD(T) and LED calculations were performed using "TightPNO" settings. 33,34 For the heme species "NormalPNO" settings 33,34 with conservative TCutPairs thresholds (10 −5 ) were also used for comparison. Hence, the so-called "NormalPNO" settings used in this study are slightly different than the standard "NormalPNO" settings. For CH 2 ···H 2 O and heme species, the resolution-of-the-identity (RI) approach was utilized in the SCF part for both the Coulomb and the exchange terms (RIJK). For all other adducts, the RI approach was only used for the Coulomb term (RIJONX, called also as RIJDX). 81,82 As default in ORCA, 93 both geometry optimizations and singlepoint energy calculations at the DLPNO-CCSD(T) level were performed using the frozen-core approximation by excluding only the 1s orbital of C, N, and O, and the 1s, 2s, and 2p orbitals of Fe from the correlation treatment.
For CH 2 adducts, DLPNO-CCSD(T) and LED energies were obtained with the aug-cc-pV5Z basis set. The corresponding LED terms are not affected much by the basis set size. As shown in Tables S2−S4, aug-cc-pVnZ basis sets, where n = D, T, Q, and 5, provide similar results. Recently, it has been shown that the DLPNO-CCSD(T) method can be used to compute accurate singlet−triplet gaps for aryl carbenes. 94 On a set of 12 aryl carbenes, the mean absolute error (MAE) is only 0.2 kcal/mol. 94 Analogously, the MAE of the DLPNO-CCSD(T) method compared with its canonical counterpart is only 0.07 kcal/mol for binding energies and 0.04 kcal/mol for singlet−triplet energy gaps for the CH 2 adducts studied in this work (see Table S5).
DLPNO-CCSD(T) and LED energies for heme species were obtained using TightPNO settings at the extrapolated complete basis set (CBS) limit by using def2-TZVP and def2-QZVP basis sets, as described previously. 39,95 However, for the sake of simplicity, the distance dependence of the DLPNO-CCSD(T) and LED energy terms was studied at the DLPNO-CCSD(T)/NormalPNO/def2-TZVP level. This methodology provides results that are in qualitative agreement with those obtained at the DLPNO-CCSD(T)/TightPNO/CBS level at the equilibrium geometry (see Tables S7−S12). Note that a large number of benchmark and application studies has already shown that the DLPNO-CCSD(T)/def2-TZVP methodology typically provides a good balance between accuracy and computational cost for relative energies. 35,40,96,97 For triplet heme species, DLPNO-CCSD(T)/def2-TZVP binding energies are in fact very close to the estimated CBS limit (see Table  S9). However, it was found that the PFe( 1 A 1g )···CO binding energy converges slowly with the basis set size and is also quite sensitive with respect to DLPNO thresholds used (see Table  S9).
For heme species, scalar relativistic effects calculated at the DLPNO-CCSD(T)/NormalPNO/def2-TZVP level by utilizing the DKH2 Hamiltonian 98,99 are quite small and nearly cancel out with the zero-point vibrational energy (ZPVE) calculated at the B3LYP-D3/def2-TZVP level (see Table S13). Therefore, the energies reported in the main paper are not corrected for these effects.
In DLPNO-CCSD(T) single-point energy calculations and geometry optimizations, the occupied orbitals were localized through the Foster-Boys and augmented Hessian Foster-Boys schemes, respectively. 100 Consistent with the closed-shell formalism, 38 the LED terms discussed demonstrate only a slight dependence to the localization scheme used for the occupied orbitals (see Table S6). In the LED calculations, pair natural orbitals (PNOs) were localized with the Pipek-Mezey 101 scheme.
Finally, it is worth mentioning that for open-shell species, QRO and UHF absolute energies differ significantly. However, the difference typically cancels out in relative energies. In the present case, QRO and UHF binding energies are essentially identical for both spin states (see Tables S5 and S7).

RESULTS AND DISCUSSION
According to eq 22, the variation (Δ) in the singlet−triplet gap (E S−T ) of methylene upon intermolecular interaction with a molecule X is equal to the difference in the 1 CH 2 ···X and 3 CH 2 ···X binding energies. Analogously, Δ in E S−T of ferrous heme (PFe) upon interacting with CO is equal to the difference in the PFe( 1 A 1g )···CO and PFe( 3 A 2g / 3 E g )···CO binding energies. In the following, these binding energies are decomposed by means of the open-shell DLPNO-CCSD(T)/ LED scheme. The results are then used to rationalize the different values of Δ obtained for the different CH 2 ···X (X = H 2 O, He, Ne, Ar, and Kr) and PFe···CO adducts studied in this work. A positive (negative) E S−T implies that the triplet state is more (less) stable, while the positive (negative) sign of Δ implies the increase (decrease) in E S−T of CH 2 upon interacting with water and rare gases. b The corresponding canonical CCSD(T) value is 9.34 kcal/mol. Experiment including zero-point vibrational energy: 9.023 kcal/mol. 102,103 Journal of Chemical Theory and Computation

. Singlet−Triplet Energy
Gap of CH 2 Adducts. The singlet−triplet energy gap of the isolated CH 2 and of CH 2 ···X adducts calculated at the DLPNO-CCSD(T)/TightPNO/aug-cc-pV5Z level is given in Table 1 with the computed Δ values obtained at the QRO, DLPNO-CCSD, and DLPNO-CCSD(T) levels of theory.
The triplet state of the bare CH 2 is calculated to be 9.28 and 9.34 kcal/mol lower in energy than its singlet state at the DLPNO-CCSD(T) and CCSD(T) levels, respectively. The inclusion of the harmonic ZPVE correction (−0.64 kcal/mol) reduces the singlet triplet gap to 8.7 kcal/mol, which is reasonably close to the experimental value of 9.023 kcal/ mol. 102,103 When interacting with water, the singlet−triplet gap reduces to 5.6 kcal/mol (Δ = −3.7 kcal/mol). Hence, the interaction stabilizes the singlet state more than the triplet. The largest contribution to this differential spin state stabilization comes from the reference energy, with almost no net contribution from the electron correlation term.
The interaction of CH 2 with rare gases is weaker than with water. At the reference level, it stabilizes the triplet more than the singlet state (Δ > 0). However, when electron correlation is incorporated, the opposite trend is observed for all adducts (Δ < 0). Importantly, Δ increases in absolute value with the polarizability of the rare gases.
The different role that electron correlation plays in the two situations above indicates that the physical mechanism responsible for the change in the singlet−triplet gap is different for methylene adducts with water and rare gases. A deeper insight into the origin of this difference comes by analyzing the 1 CH 2 ···X and 3 CH 2 ···X binding energies, which determine the overall Δ through eq 22. An in-depth analysis of this aspect is reported in the following by means of the LED scheme.
4.1.2. Binding Energies of CH 2 Adducts. The 1 CH 2 ···X and 3 CH 2 ···X binding energies are reported in Table 2. In the same table their decomposition into geometric preparation and interaction energy is also given. The latter is further decomposed into its reference and correlation energy contributions. Consistent with the trend of Δ previously discussed, water and rare gases interact more strongly with 1 CH 2 than with 3 CH 2 , thus resulting in a lowering of the singlet−triplet gap in CH 2 ···X compounds (Δ < 0).  Table 3). Thus, ΔE int ref is practically negligible, demonstrating that dynamic electron correlation is responsible for the stability of the 3 CH 2 ···H 2 O adduct. In particular, London dispersion is one of the most important energy components of the interaction, as demonstrated by the large E disp C-CCSD /ΔE ratio of 0.57.
A different picture emerges in 1 CH 2 ···H 2 O. In this case, the interacting species are closer (see Figure 1)  Even though the correlation contribution to the 1,3 CH 2 ··· H 2 O interaction is largely dominated by the London dispersion term E disp C-CCSD (see Table 3) and significantly contributes to the intermolecular interaction, its magnitude is almost identical for both spin states. Hence, the decrease in the singlet−triplet gap of CH 2 while interacting with water is driven by the fact that 1 CH 2 ···H 2 O shows a much larger electrostatic interaction than 3  It is worth mentioning here that the singles term included in ΔE no-disp C-CCSD amounts to −2.4 kcal/mol for the 3 CH 2 fragment of all of the adducts studied in this work (see Table S5). In contrast, the adducts containing 1 CH 2 features a negligible singles term (the Brillouin's theorem is satisfied in this case). In general, the magnitude of the singles term can be used to identify the fragment in which the unpaired electron is localized.
4.1.4. LED Analysis of the CH 2 ···Rg Interaction. The interaction of CH 2 with rare gases is relatively weak. In both singlet and triplet states (see Tables 2 and 3

Journal of Chemical Theory and Computation
Article energies due to the distortion of electron clouds of both CH 2 and rare gases dominate over the sum of the attractive electrostatic and exchange interactions. Hence, the overall ΔE is dominated by electron correlation and, in particular, by the London dispersion contribution (E disp C-CCSD /ΔE > 1.5). It is worth mentioning here that the first explanation of the attraction between two nonpolar molecules was given by F. London. 106,107 An approximated expression for the London dispersion energy between two atoms (X and Y), i.e., E disp,L can be written where C 6 is the atom pairwise induced dipole−induced dipole interaction coefficient, I X and I Y are the first ionization potential of the interacting X and Y molecules, α X and α Y are the polarizabilities of X and Y, and r is the distance between X and Y. It can be assessed how well this simple London equation quantifies the dispersion interaction energy between CH 2 and rare gases by comparing the London dispersion energies obtained from eq 23 with those derived from the LED method. To do that we computed ionization potential and numerical polarizabilities of the isolated CH 2 and rare gas atoms at the DLPNO-CCSD(T)/aug-cc-pV5Z level. The computed values are reasonably close to the available experimental values (see Table 4). 108−114 As an approximation we take r as the distance between the rare gas and the carbon atom (see Figure 1) at the equilibrium geometry of CH 2 ···Rg adducts.
The correlation between the dispersion energies obtained with London formula and the corresponding LED values for the CH 2 ···Rg systems studied in this work is given in Figure 4. Despite the simplicity of the London equation and the various approximations adopted, the dispersion energy estimates obtained with the two methods show a reasonably good linear correlation, with an R 2 larger than 0.98 for both the 1 CH 2 ···Rg and the 3 CH 2 ···Rg series. In all cases, London dispersion increases with the increase of the polarizability of the rare gas atoms.   Figure 3. Valence electron configurations and molecular electrostatic potential maps of 1 CH 2 (left) and 3 CH 2 (right) projected onto the corresponding molecular electron densities calculated at the UHF/ aug-cc-pV5Z level in the unit of E h /e. Red region identifies lowest electrostatic potential and thus highest electron density, while blue region identifies the opposite. A 2g ) energy gaps of the free and CO-bound heme are given in Table 5 at various levels of theory. The corresponding Δ values are also reported in the same table.
Experimental studies on bare iron−porphyrin (see ref 115 for a review) agree on its ground-state multiplicity (triplet) but differ in the interpretation of the ground-state electronic configuration. In particular, Mossbauer data of iron(II)− tetraphenylporphyrln (FeTPP) were interpreted as an indication of either a 3 E g 116 or a 3 A 2g 70,117 ground state. Ligand field calculations that are consistent with the magnetic susceptibility measurements predict an 3 A 2g ground state. 118 Consistently, previous hybrid density functional, CASPT2, and several CI studies find the 3 A 2g state more stable than the 3 E g state by 2 kcal/mol or less. 89,115 Consistently, the present B3LYP-D3/def2-TZVP and DLPNO-CCSD(T)/TightPNO/ CBS calculations find the 3 A 2g state of heme to be only 0.96 and 1.97 kcal/mol more stable than the 3 E g state, respectively (see Table 5).
The d z 2 orbital of Fe(II), which is doubly occupied in the 3 A 2g state, is destabilized upon the binding of an axial ligand to the ferrous heme. 89,119 The amount of this destabilization and thus spin-state energetics depends strongly on the nature of the axial ligand. For example, upon imidazole binding, the d z 2 orbital raises in energy and becomes singly occupied. Thus, the ground state changes from triplet to quintet, and the lowest triplet changes from 3 A 2g to 3 E g . 89 As CO is a strong field ligand, the destabilization of the d z 2 orbital upon its binding to heme is even larger. In fact, PFe−CO systems having side chains feature a singlet ground state with a formally empty d z 2 orbital. 67 Consequently, the identity of the most stable triplet state of CO-bound heme complexes is rarely studied. The present calculations find the 3 E g state of the CO-bound heme to be 2.16 and 0.57 kcal/mol more stable than 3 A 2g at the B3LYP-D3/def2-TZVP and DLPNO-CCSD(T)/TightPNO/ CBS levels, respectively. Hence, CO binding probably reverses the energetic order of the 3 E g and 3 A 2g states, consistent with the imidazole-bound heme case. These changes can be rationalized in terms of the destabilization of the d z 2 orbital of Fe(II).
Unfortunately, no direct experimental measure of E S−T gaps exists for these complexes. At the B3LYP-D3/def2-TZVP and DLPNO-CCSD(T)/TightPNO/CBS levels, the E S−T gap of the free heme is 33.32 and 32.13 kcal/mol relative to 3 E g , while  A positive (negative) E rel implies that the 3 A 2g state is more (less) stable than the 3 E g , while the negative (positive) sign of Δ implies the increase (decrease) in the stability of the 3 E g state relative to the 3 A 2g state of PFe upon interacting with CO. b In this case, E rel corresponds to E S−T . A positive (negative) E rel implies that the triplet state is more (less) stable than the singlet state, while the positive (negative) sign of Δ implies the increase (decrease) in E S−T of PFe upon interacting with CO.

Journal of Chemical Theory and Computation
Article it is 34.88 and 34.10 kcal/mol relative to 3 A 2g (see Table 5). As detailed in Table S7, these figures are only weakly affected by the technical parameters of the calculations and are consistent with those obtained previously at the CASPT2 level (∼35 kcal/mol). 89 As seen in Table 5, upon interacting with CO, at the DLPNO-CCSD(T)/TightPNO/CBS level, the E S−T gap relative to the 3 E g and 3 A 2g states reduces to −6.50 and −7.07 kcal/mol (Δ = −38.63 and −41.17 kcal/mol), respectively. Hence, CO binding significantly stabilizes the 1 A 1g state (which becomes the ground state), consistent with the above-mentioned experimental findings on related systems. An in-depth discussion of the physical mechanism behind this differential spin state stabilization is reported in the following sections.
4.2.2. Binding Energy of the Heme Adducts. The calculated binding energies of PFe( 1 A 1g )···CO, PFe( 3 E g )··· CO, and PFe( 3 A 2g )···CO adducts at the DLPNO-CCSD(T)/ TightPNO/CBS level are reported in Table 6. In the same table, their decomposition into geometric preparation and interaction energy is also given. The latter is further decomposed into reference and correlation energy contributions.
The PFe( 1 A 1g )···CO interaction is much stronger than the PFe( 3 A 2g / 3 E g )···CO ones, as apparent from their large negative Δ values of −41.17/−38.63 kcal/mol (see Table 6). This leads to the observed reduction of the singlet−triplet gap of heme derivatives shown in Table 5. However, at the reference level, the PFe( 1 A 1g )···CO interaction is strongly repulsive (ΔE int ref > 0). Thus, electron correlation counteracts this repulsion in the singlet state and makes the overall interaction significantly attractive.
On the technical side, PFe( 3 E g / 3 A 2g )···CO binding energies do not depend significantly on the computational settings used in the DLPNO-CCSD(T) calculations. In contrast, the binding energy of the PFe( 1 A 1g )···CO adduct is more sensitive to the DLPNO thresholds and basis sets (see Table S8) adopted. For example, the PFe( 1 A 1g )···CO binding energy is −42.84 kcal/ mol with TightPNO/CBS and −36.90 kcal/mol with Normal-PNO/def2-TZVP settings. Hence, NormalPNO/def2-TZVP settings should only be used if one is interested in mechanistic tendencies rather than in accurate quantitative estimates of binding energies of heme species.
The PFe( 1 A 1g )···CO binding energy calculated at the DLPNO-CCSD(T)/TightPNO/CBS (−42.84 kcal/mol), B3LYP-D3/def2-TZVP (−46.95 kcal/mol), and CASPT2 89 (−51.3 kcal/mol) levels varies significantly. Relative to the ground PFe( 3 A 2g ) state, these estimates are −8.74, −12.67, and −16.4 kcal/mol, respectively. As already mentioned, experimental data are not available for this system. The only experimentally determined binding energy on a related system was obtained for a heme derivative in which four tetrakis(4sulfonatophenyl) anions (tpps) are bound to the four meso carbons of porphyrin connecting pyrrole moieties. In this case, the experimental gas-phase binding energy of the PFe( 1 A 1g )− CO adduct was measured to be −15.85 kcal/mol relative to the ground state of PFe. 120 It is worth stressing here that the four tpps anions (i.e., a total charge of −4) bend the porphyrin moiety significantly and thus are expected to alter the electronic structure of the system as well. In fact, the binding energy of this system is predicted to be 2.05 kcal/mol larger in absolute value than that of side chain free PFe( 1 A 1g )···CO at the B3LYP-D3/def2-TZVP level of theory.
4.2.3. LED Analysis of the Interaction between Heme and CO. LED terms of the PFe( 1 A 1g / 3 E g / 3 A 2g )···CO interactions are given in Table 7 at the DLPNO-CCSD(T)/TightPNO/ CBS level. For an in-depth analysis of these interactions, we also performed DLPNO-CCSD(T)/NormalPNO/def2-TZVP single-point energy calculations on a series of structures obtained from constrained B3LYP-D3 geometry optimizations in which the Fe−C distance (r Fe ··· C ) was varied from ∼1.7 to ∼4.5 Å with an increment of 0.2 Å. The resulting energy profiles are reported in Figure 5 together with their decomposition into the various LED terms.
As mentioned above, PFe( 3 A 2g ) features a doubly occupied d z 2 orbital that points toward the CO lone pair located on the carbon atom, while in PFe( 3 E g ) the d z 2 orbital is singly occupied. For this reason, PFe( 3 A 2g )···CO features a steeper repulsive wall than PFe( 3 E g )···CO. This leads to a longer Fe− C bond distance in PFe( 3 A 2g )···CO (3.27 Å) than in PFe( 3 E g )···CO (2.29 Å). Accordingly, the PFe( 3 A 2g )···CO interaction (−1.67 kcal/mol) is weaker than the PFe( 3 E g )··· CO one (−4.21 kcal/mol), and all of the LED terms of PFe( 3 A 2g )···CO are smaller than those of PFe( 3 E g )···CO at the

Journal of Chemical Theory and Computation
Article B3LYP-D3 equilibrium geometry (see Tables 6 and 7, and  Tables S9−S11). For both triplet states, the interaction between PFe-( 3 E g / 3 A 2g ) and CO is repulsive at the QRO level (top panels of Figure 5 and Table 7) in the short range. Consistent with this picture, the LED decomposition of the reference PFe( 3 E g / 3 A 2g )···CO interaction energies (central panels of Figure 5 and Table 7) shows that the repulsive electronic preparation energies due to the distortion of the electron clouds of PFe( 3 E g / 3 A 2g ) and CO are always larger or approximately equal in magnitude than the sum of the attractive electrostatic and exchange interaction energies.
Hence, dynamic electron correlation is essential to the stability of PFe( 3 E g / 3 A 2g )···CO adducts. Note that the B3LYP-D3 equilibrium geometries for the PFe( 3 E g / 3 A 2g )···CO adducts feature shorter intermolecular distances (of about 0.6/0.2 Å) than the DLPNO-CCSD(T) minima. However, this has only a small impact on the DLPNO-CCSD(T) binding energy due to the flatness of the PES.
The decomposition of the correlation binding energy (see bottom panels of Figure 5 and Table 7) shows that the dispersion energy dominates the PFe( 3 E g / 3 A 2g )···CO interaction. The corresponding E disp C-CCSD value amounts to −7.53/ −2.61 kcal/mol at the B3LYP-D3 equilibrium geometry. In fact, both PFe( 3 E g )···CO and PFe( 3 A 2g )···CO adducts have a E disp C-CCSD /ΔE ratio larger than 1.5 (see Tables 7 and S11) and can be both described as a van der Waals adduct mainly stabilized by London dispersion forces.
By contrast, the nature of the PFe( 1 A 1g )···CO interaction is completely different. The energy profile computed at the QRO level (top right panel of Figure 5) shows a shallow minimum at about 2.5 Å. The corresponding DLPNO-CCSD(T) energy profile shows a very deep minimum (NormalPNO/def2-TZVP, − 36.90 kcal/mol; TightPNO/CBS, − 42.84 kcal/mol) at ∼1.7 Å, where the QRO interaction energy is repulsive. Thus, electron correlation affects significantly both the stability and the geometry of the PFe( 1 A 1g )···CO adduct.
The LED decomposition of the reference PFe( 1 A 1g )···CO interaction energy (central right panel of Figure 5) shows that the electrostatic interaction is larger than the corresponding electronic preparation at the QRO minimum. Importantly, the PFe( 1 A 1g )···CO interaction at the QRO level is associated with larger LED components for both repulsive and attractive terms than the PFe( 3 E g / 3 A 2g )···CO interaction at all Fe−C distances. Hence, the electronic clouds of the interacting fragments are distorted more significantly and interact stronger in PFe-( 1 A 1g )···CO than in PFe( 3 E g / 3 A 2g )···CO, indicating that a significant polarization of the fragments occurs upon CO binding already at the QRO level. , respectively. Hence, these three correlation terms are responsible for the significant strength of the PFe( 1 A 1g )··· CO interaction, which in turn determines the reduction observed in the singlet−triplet gap of PFe upon CO binding. In particular, the large magnitude of the nondispersive component suggests that electron correlation significantly affects the polarization of the interacting fragments. 76 To analyze this aspect in more detail, we obtained 3D contour plots of the one electron density difference function Δρ(x,y,z) describing the electron density rearrangement taking place upon CO binding for each spin state. The Δρ(x,y,z) function is computed as the difference between the electron density of the PFe−CO adduct (for a given spin state) and the sum of the electron densities of PFe (at the same spin state) and CO frozen at their in-adduct geometries. In order to investigate electron correlation effects to the PFe···CO binding, Δρ was divided into a reference (Δρ QRO ) and an unrelaxed DLPNO-CCSD correlation (Δρ C ) contribution calculated via the solution of Λ equations. 121,122 The corresponding contour plots are shown in Figure 6.
Consistent with the fact that the LED describes the PFe( 3 E g / 3 A 2g ) ··· CO adducts as van der Waals complexes held together by London dispersion, the corresponding Δρ QRO and Δρ C contour plots show that only a small charge rearrangement takes place upon bond formation in the triplet states. In contrast, the PFe( 1 A 1g )···CO interaction is accompanied by a significant charge accumulation in the region of the bond, consistent with the fact that the PFe( 1 A 1g )···CO interaction is associated with larger electrostatic, exchange, and electronic preparation energy components. This is evident from the contour plots of the corresponding Δρ QRO and Δρ C functions and consistent with the fact that the PFe( 1 A 1g )···CO interaction is associated with a large nondispersive correlation term.
The contour plots just discussed are consistent with the variations in d-orbital populations occurring upon CO binding (Δq) reported in Table 8. The Δq values are negligible for all d orbitals in PFe( 3 E g / 3 A 2g )−CO adducts, while they are significant in the PFe( 1 A 1g )−CO adduct. In particular, the population of the formally empty d z 2 orbital of the singlet state increases significantly upon CO binding. At the same time, the population of the formally doubly occupied d xz and d yz orbitals decreases. On the basis of these findings, it is clear that the Fe(II)···CO interaction in the singlet state can be understood in terms of the well-known Deward−Chatt−Duncanson

Journal of Chemical Theory and Computation
Article (DCD) bonding model, which was originally introduced to discuss the binding of olefins to transition metals. 123 A significant σ-donation of charge takes place from the lone pair located on the CO ligand to the empty d z 2 orbital located on the Fe(II) of the singlet adduct, while at the same time πbackdonation occurs from the d xz and d yz orbitals to the empty π* orbitals of CO, consistent with natural bond orbital (NBO) 124 results (see Scheme S1) and the elongation of the C−O bond length 125 in the PFe( 1 A 1g )···CO adduct (1.141 Å) with respect to that of free CO (1.125 Å, which is the same as in the PFe( 3 E g / 3 A 2g )···CO complexes). The importance of πbackdonation has been already pointed out on related compounds based on vibrational spectroscopy measurements and DFT computations. 126−128 A schematic representation of the binding modes just discussed, along with the associated amount of the electron transfer (Δq) based on the Mulliken population analysis, is reported in Scheme 2. Note that electron correlation significantly affects the magnitude of the DCD components, consistent with our previous findings on agostic complexes. 41

CONCLUSIONS
The LED analysis in the DLPNO-CCSD(T) framework decomposes the interaction energy of an arbitrary number of fragments of a closed-shell molecular adduct into repulsive electronic preparation and interfragment electrostatic, exchange, and London dispersion contributions. In this study, we presented an extension of the LED scheme to open-shell molecular systems that was implemented in the ORCA program package. As a first illustrative case study this scheme was applied to investigate the mechanism that governs the change in the singlet−triplet energy gap of CH 2 upon interaction with water and rare gases (Rg) and of heme (PFe) upon interaction with CO.
The CH 2 ···X interaction (X = H 2 O, He, Ne, Ar, and Kr) was found to be attractive for both the triplet 3 CH 2 and the singlet 1 CH 2 methylene. The interaction is stronger with 1 CH 2 than with 3 CH 2 , resulting in a lowering of the singlet−triplet energy gap (E S−T ). The LED analysis of the CH 2 ···H 2 O interaction showed that electrostatics dominates the interaction for both spin states of methylene. The lowering of E S−T can thus be understood in terms of the larger electrostatic interaction in 1 CH 2 ···H 2 O than that in 3 CH 2 ···H 2 O, consistent with chemical intuition. In contrast, the interaction of methylene with rare gases (Rg) is dominated by London dispersion forces for both spin states. In this case, the lowering of E S−T arises from the stronger London dispersion in 1 CH 2 ···Rg than that in 3 CH 2 ··· Rg. This is consistent with the larger polarizability of 1 CH 2 than that of 3 CH 2 .
As regards heme systems, it was found that the PFe···CO interaction is dominated by the correlation contribution for all the spin states of PFe investigated in this work ( 1 A 1g , 3 E g , and 3 A 2g ). Although PFe is a triplet in its ground state, the PFe( 1 A 1g )···CO interaction is significantly stronger than the interaction of PFe in the two lowest-lying triplets ( 3 E g and 3 A 2g ) with CO, leading to a lowering of E S−T in the adduct. In fact, the LED analysis demonstrated that PFe( 3 E g )···CO and the PFe( 3 A 2g )···CO adducts can be described as van der Waals complexes stabilized by weak dispersion forces. In contrast, the PFe( 1 A 1g )···CO bond can be described as a strong donor/ acceptor interaction in which electron density is transferred from the carbon lone pair into the formally empty d z 2 orbital on the central iron atom (Fe(II) ← CO σ-donation) and from the filled d xz and d yz orbitals into the π* antibonding orbitals on the CO molecule (Fe(II) → CO π-backdonation).
It should be stressed here again that the results just discussed refer to the minimal porphyrin model in the gas phase. The presence of heme side chains, solvent, or protein matrices may drastically change the nature of the Fe···CO interaction, especially for triplet heme adducts. However, the tools presented in this work appear to be particularly powerful for the in-depth understanding of intermolecular interactions in adducts with complicated electronic structures, such as heme species.

Journal of Chemical Theory and Computation
Article