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 symmetry-adapted 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 difChemical">ference 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 (<Chemical">span class="Chemical">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 exchange-repulsion) terms.[10−14]

The large majority of Chemical">EDA sc<Chemical">span class="Chemical">hemes 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 Chemical">EDA approach called local energy decomposition (LED).[38] In this sc<Chemical">span class="Chemical">heme, the interaction energy between two or more fragments is calculated at the DLPNO-CCSD(T) level and is then decomposed into a repulsive intramolecular energy term called electronic preparation, plus a series of intermolecular energy terms such as electrostatic, quantum mechanical exchange, and London dispersion interactions.[38] This scheme has been already applied in the context of H-bond interactions,[38,39] frustrated Lewis pairs,[40] agostic interactions,[41] and the interactions in dipnictenes.[42] Quantitative comparisons of the terms of LED and SAPT have been also reported.[38,39]

Herein, we present an extension of the LED scChemical">heme to open-shell systems in the framework of the DLPNO-CCSD(T) method. The present developments were made accessible through the recent ORCA 4.1 program release.[43,44] This opens an unprecedented opportunity for understanding and thus controlling a wide range of chemical and biological processes, eChemical">specially in view of the limited number of previously developed <Chemical">span class="Chemical">EDA and SAPT approaches for open-shell molecular systems.[45−52] The new scheme is used to investigate the mechanism that governs the change in the singlet–triplet energy gap (ES–T) of carbenes upon their noncovalent interactions with various chemical species and of heme upon CO binding. The systems studied are as shown in Figure .

Figure 1

Investigated CH2···X and PFe···CO adducts and the associated equilibrium intermolecular distances obtained at the DLPNO-CCSD(T)/TightPNO/aug-cc-pVTZ and B3LYP-D3/def2-TZVP levels, respectively.

Investigated CH2···X and PChemical">Fe···CO adducts and the associated equilibrium intermolecular distances obtained at the DLPNO-CCSD(T)/TightPNO/aug-cc-pVTZ and B3LYP-D3/def2-TZVP levels, reChemical">spectively.

Chemical">Carbenes are highly reactive neutral divalent molecules that exist in triplet and <Chemical">span class="Chemical">singlet spin states.[53−62] The singlet and triplet states of carbenes yield different reaction products. For example, the simplest carbene CH2 (methylene), being one of the most reactive molecules, undergoes stereospecific reactions in the singlet state and stereoselective reactions in the triplet state.[55−62] Although many carbenes have a triplet ground state in the gas phase, the singlet state is typically stabilized significantly in several solvents, causing an equilibrium between the two spin states or even leading to a singlet ground state.[55−62]

Herein, the presently introduced open-shell variant of the DLPNO-CCSD(T)/LED scChemical">heme is applied to a series of prototypical molecular systems (Figure ) of the type CH2···X, in which methylene (in both its <Chemical">span class="Chemical">singlet 1CH2 and its triplet 3CH2 spin states) interacts through noncovalent interactions either with water or with rare gases (Rg = He, Ne, Ar, and Kr). The variation in the ES–T of methylene upon adduct formation is discussed in terms of the 1CH2···X and 3CH2···X interaction energy components extracted from the LED scheme.

The other illustrative case study discussed here involves Chemical">heme, which is an <Chemical">span class="Chemical">iron-coordinated porphyrin (P) derivative that is essential for the function of all aerobic cells.[63] It serves as a chromophore for many proteins including hemoglobin, myoglobin, and cytochromes, and thus regulates diverse biological functions.[64,65]

Chemical">Heme owes its functional diversity in the protein matrix to its dif<Chemical">span class="Chemical">fering side chain environments, bound axial ligands and their environments, coordination number, and oxidation and spin states of the iron center.[65] In particular, the binding and dissociation reactions of CO to the ferrous heme in the protein matrix have significant impacts on its role in respiration and regulation processes.[66] Experimentally, it is known that ferrous iron-coordinated porphyrin (PFe, also denoted hereafter as free heme) derivatives and their CO-bound PFe···CO analogs have triplet and singlet ground states, respectively.[67−70] Thus, the CO binding reverses the spin-state ordering of heme derivatives. However, the chemical mechanism behind this experimental observation is not clear.

Herein, the LED scChemical">heme is used to provide a first insight into the factors contributing to the stabilization of the lowest-lying <Chemical">span class="Chemical">singlet state (1A1g) and two close-lying lowest-energy triplets (3Eg and 3A2g) of PFe (see Scheme ) upon CO binding, i.e., the PFe(1A1g/3Eg/3A2g)···CO interactions. As seen in Scheme , the 1A1g singlet has an electronic configuration of (dx)2(d)2(d)2(dz)0(d)0 while the 3A2g and 3Eg triplets have the (d)2(d)2(d)1(d)1(d)0 and (d)2(d)2(d)1(d)1(d)0 configurations, respectively.

Scheme 1

Simple Schematic Representation of the Singlet and Low-Lying Triplet Electronic Configurations of PFe and PFe···CO

One of the degenerate pair of 3Eg configurations is shown. The other configuration is clearly (d)2(d)1(d)2(d)1(d)0.

One of the degenerate pair of Chemical">3Eg configurations is shown. The other configuration is clearly (d)2(d)1(d)2(d)1(d)0.

As a prototype case study, we focus on the bare porhyrin that has no peripheral or other axial substituents (see Figure ). It should be emphasized here that spin-state ordering of <Chemical">span class="Chemical">heme species is dependent on such substituents as well as on their surrounding media.[71] For example, if a single imidazole ligand is bound to the free heme rather than CO then the quintet state becomes the ground state both in the solid phase and in protein matrices.[72] A thorough analysis of the role of different substituents and surrounding media on the singlet–triplet energy gap of heme species would require an extensive systematic study, which is beyond the scope of this study.

The paper is organized as follows. In section the theory of the open-shell LED scChemical">heme is given, while computational details are reported in section . In section , the ES–T <Chemical">span class="Chemical">singlet–triplet energy gap of methylene, heme, and their investigated adducts is reported at the DLPNO-CCSD(T) level. Then in the framework of the presently introduced open-shell DLPNO-CCSD(T)/LED method, the variations in ES–T of methylene 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 for the PFe(1A1g/3A2g/3Eg)···CO interaction and given in Figure S1 for the CH2···H2O interaction. The last section of the paper is devoted to discussion of the results and concluding remarks of the study.

Theory

Theoretical Background

The theory and implementation of the DLPNO-CCSD(T) method in both its closed-shell and its open-shell variants have been described in detail in a series of recent publications.[25−34,38,73,74] We thus only recall here the energy expression for the open-shell DLPNO-CCSD(T) method, which is decomposed in the present variant of the LED scChemical">heme.

The total DLPNO-CCSD(T) energy can be expressed as a sum of reChemical">ference and correlation energies E = Eref + EC. In the closed-shell variant of DLPNO-CCSD(T), the re<Chemical">span class="Chemical">ference determinant is typically the Hartree–Fock (HF) determinant, and hence, Eref corresponds to the HF energy. In the open-shell variant, Eref is the energy of a high-spin open-shell single determinantal function, consisting of a single set of molecular orbitals. The energy of such a determinant can be expressed aswhere i and j are the occupied orbitals of the reference determinant, n and n denote the occupation numbers (1 or 0) for α and β electrons, respectively, n = n + n; (ii|jj) and (ij|ij) are two electrons integrals in Mulliken notation, and Z (Z) is the nuclear charge of atom A (B) at position ().

The DLPNO-CCSD(T) correlation energy can be written essentially as a sum of electron-pair correlation energy (ε, where i and j denote the localized orbitals) contributions plus the perturbative triples correction (EC–(T)). Local second-order many-body perturbation theory is used to divide the ε 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 readswhere EC–SP denotes the energy contribution from the strong pairs, EC–WP is the correlation contribution from the weak pairs, and the last term is associated with the perturbative triples correction EC–(T).

The dominant strong pair contribution readswhere the indices with an overbar denote β spin–orbitals and those without overbar are used for α <Chemical">span class="Gene">spin–orbitals. The first two terms represent the contribution from the single excitations (a, the singles PNOs; t, 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 ε, and terms denote α–α, β–β, and α–β pair correlation energies, defined aswhere a and b are the PNOs belonging to the ij pair, (ia | jb) are the two-electrons integrals, τ = t + tt – tt and are the cluster amplitudes in the PNO basis, and t, t, t are the doubles and singles amplitudes of the coupled cluster equations.

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 aswhere ΔEgeo-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. ΔEint is the interaction energy between the fragments X and Y frozen in the geometry they have in the adduct XY, which is defined aswhere EXY is energy of the XY adduct while EX (EY) is the energy of the isolated X (Y) fragment frozen at the same geometry as that in the adduct. ΔEint can be decomposed into a reChemical">ference contribution ΔEintref (i.e., the contribution to the interaction energy from the QRO determinant in the open shell case) and a correlation contribution ΔEintC

In the DLPNO-CCSD(T) framework, the orbitals of the reChemical">ference 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 re<Chemical">span class="Chemical">ference energy of the XY adduct ErefXY (eq ) into intra- and intermolecular contributions in exactly the same way as it was described for the closed-shell case[38,40]

The intermolecular reChemical">ference energy can be further partitioned into electrostatic (Eelstat) and exchange (Eexch) interactions. Accordingly, the ΔEintref term can be written as

The electronic preparation energy ΔEel-prepref 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. Eelstat and Eexch 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 Eelstat 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 Eelstat in the closed-shell case,[76] which is in principle also applicable to the open-shell case.

The correlation contribution to the interaction energy EintC can thus be expressed as a sum of three contributionswhere ΔEintC-SP, ΔEintC-WP, and ΔEC-(T) are the strong pairs, weak pairs, and triples correction components of the correlation contribution to the interaction energy, respectively. The ΔEintC-SP, ΔEintC-WP, and ΔEintC-(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 ΔEintC-SP contribution, a more sophisticated approach is used in order to provide a clear-cut definition of London dispersion. The Chemical">EC-SPXY term can be written as a sum of the contributions of single and double excitations according to eqs –6. By exploiting the localization of both the occupied and the virtual orbitals in the DLPNO-CCSD(T) framework, <Chemical">span class="Chemical">EC-SPXY can be divided into five different contributionswhere each energy term contains contributions from the αα, ββ, and αβ excitations constituting the correlation energy, as given in eqs –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 Chemical">EC-SP(Y) and <Chemical">span class="Chemical">EC-SPCT(X←Y) terms can be obtained from eqs and 15 by exchanging the X and Y labels in all terms. Note that the last term in EC-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 .

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.

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 transChemical">fer excitations from X to Y are shown. Analogous charge trans<Chemical">span class="Chemical">fer excitations also exist from Y to X.

EC-SP(X) and EC-SP(Y) describe the correlation contribution from excitations occurring within the same fragment (also called intrafragment correlation). The dynamic charge trans<Chemical">span class="Chemical">fer (CT) terms EC-SPCT(X→Y) and EC-SPCT(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. EC-SPDISP(X,Y) describes the energy contribution associated with genuine and exchange dispersion (see Figure ).

For the sake of simplicity, it may be useful to combine several terms. For example, ΔEel-prepC-SP (i.e., the difChemical">ference between the intrafragment contributions <Chemical">span class="Chemical">EC-SP(X/Y) and the strong pair correlation energy of the isolated fragments) is always positive, while EC-SPCT(X→Y) and EC-SPCT(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 (ΔEno-dispC-SP)where ΔEno-dispC-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 Chemical">EDISPC-SP in order to obtain the total diChemical">spersion contribution at the DLPNO-CCSD level. We label this summation as EdiChemical">spC-CCSD. The remaining part of the correlation interaction energy at the DLPNO-CCSD level is labeled as ΔEno-diChemical">spC-CCSD. Therefore

Collecting all the terms we obtain for the binding energyEquation is the base of our LED scChemical">heme 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 closed-shell 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).

LED for the Analysis of Singlet–Triplet Energy Gap

As an illustrative example of the usefulness and applicability of LED scChemical">heme, we discuss in this work the influence of the dif<Chemical">span class="Chemical">ferent LED components on the singlet–triplet energy gap (ES–T) of a molecule Y upon formation of a weak intermolecular interaction with X. The ES–T gap of the corresponding 1,3Y···X adduct can be written asThe same quantity for the free Y reads asThe variation in the singlet–triplet gap of Y upon interaction with a molecule X (Δ) can be obtained by subtracting eqs and 21Thus, Δ equals the difference between the ΔE(1Y···X) and ΔE(3Y···X) binding energies. In the following, the LED scheme is used to decompose ΔE(1Y···X) and ΔE(1Y···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.

Geometries

The CH2 molecule and its adducts (see Figure ) 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 CH2···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 1CH2 and 3CH2 at the DLPNO-CCSD(T)/aug-cc-pVTZ level are shown in Figure (see section S1.1 for their coordinates).

The Chemical">ferrous heme Chemical">species and their CO adducts (see Figure ) were fully optimized using the RIJCOSX approximation[81,82] at the B3LYP[83−85] level incorporating the atom-pairwise diChemical">spersion 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 P<Chemical">span class="Chemical">Fe(1A1g)···CO, PFe(3A2g)···CO, and PFe(3Eg)···CO interactions. The optimized coordinates of the 1A1g, 3Eg, and 3A2g states of PFe and PFe···CO (see Figure ) 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]

Single-Point DLPNO-CCSD(T) and LED Calculations

DLPNO-CCSD(T) and LED calculations of both Chemical">singlet and triplet states were performed with the open-shell DLPNO-CCSD(T) algorithm. Note that open-shell and closed-shell DLPNO-CCSD(T) implementations give almost identical LED energy terms for closed-shell Chemical">species (see Tables S3 and S4). The distance dependence of the open-shell LED terms is given in Figure S1 for CH2···<Chemical">span class="Chemical">H2O and in section for PFe···CO.

Unless otherwise specified, DLPNO-CCSD(T) and LED calculations were performed using “TightPNO” settings.[33,34] For the Chemical">heme Chemical">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 dif<Chemical">span class="Chemical">ferent than the standard “NormalPNO” settings. For CH2···H2O 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 single-point 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 CH2 adducts, DLPNO-CCSD(T) and LED energies were obtained with the aug-cc-pV5Z basis set. The corresponding LED terms are not afChemical">fected much by the basis set size. As shown in Tables S2–S4, aug-cc-<Chemical">span class="Species">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 CH2 adducts studied in this work (see Table S5).

DLPNO-CCSD(T) and LED energies for Chemical">heme Chemical">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 <Chemical">span class="Chemical">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(1A1g)···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 Chemical">heme Chemical">species, scalar relativistic ef<Chemical">span class="Chemical">fects 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-<Chemical">span class="Species">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 difChemical">fer significantly. However, the dif<Chemical">span class="Chemical">ference 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 , the variation (Δ) in the Chemical">singlet–triplet gap (ES–T) of methylene upon intermolecular interaction with a molecule X is equal to the dif<Chemical">span class="Chemical">ference in the 1CH2···X and 3CH2···X binding energies. Analogously, Δ in ES–T of ferrous heme (PFe) upon interacting with CO is equal to the difference in the PFe(1A1g)···CO and PFe(3A2g/3Eg)···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 CH2···X (X = H2O, He, Ne, Ar, and Kr) and PFe···CO adducts studied in this work.

CH2···X Interaction

Singlet–Triplet Energy Gap of CH2 Adducts

The Chemical">singlet–triplet energy gap of the isolated CH2 and of CH2···X adducts calculated at the DLPNO-CCSD(T)/TightPNO/aug-cc-pV5Z level is given in Table with the computed Δ values obtained at the QRO, DLPNO-CCSD, and DLPNO-CCSD(T) levels of theory.

Table 1

Calculated Singlet-Triplet Energy Gap (ES–T, in kcal/mol) of the Free CH2 and Its Water- and Rare Gas-Interacted Structures Together with the Variation (Δ) of ES–T Relative to the Free CH2 at the Reference QRO/aug-cc-pV5Z, DLPNO-CCSD/TightPNO/aug-cc-pV5Z, and DLPNO-CCSD(T)/TightPNO/aug-cc-pV5Z Levelsa

	ES–T	Δ
molecule	DLPNO-CCSD(T)	DLPNO-CCSD(T)	DLPNO-CCSD	reference (QRO)
CH2	9.28b	0.00	0.00	0.00
CH2···H2O	5.57	–3.71	–3.62	–3.67
CH2···He	8.97	–0.31	–0.26	0.18
CH2···Ne	8.99	–0.29	–0.21	0.25
CH2···Ar	8.23	–1.05	–0.76	–0.06
CH2···Kr	7.40	–1.88	–1.31	2.05

A positive (negative) ES–T implies that the triplet state is more (less) stable, while the positive (negative) sign of Δ implies the increase (decrease) in ES–T of CH2 upon interacting with water and rare gases.

The corresponding canonical CCSD(T) value is 9.34 kcal/mol. Experiment including zero-point vibrational energy: 9.023 kcal/mol.[102,103]

A positive (negative) ES–T implies that the triplet state is more (less) stable, while the positive (negative) sign of Δ implies the increase (decrease) in ES–T of CH2 upon interacting with Chemical">water and rare gases.

The triplet state of the bare CH2 is calculated to be 9.28 and 9.34 kcal/mol lower in energy than its Chemical">singlet state at the DLPNO-CCSD(T) and CCSD(T) levels, reChemical">spectively. The inclusion of the harmonic ZPVE correction (−0.64 kcal/mol) reduces the <Chemical">span class="Chemical">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 Chemical">water, the <Chemical">span class="Chemical">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 CH2 with rare gases is weaker than with Chemical">water. At the re<Chemical">span class="Chemical">ference 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 difChemical">ferent role that electron correlation plays in the two situations above indicates that the physical mechanism reChemical">sponsible for the change in the <Chemical">span class="Chemical">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 1CH2···X and 3CH2···X binding energies, which determine the overall Δ through eq . An in-depth analysis of this aspect is reported in the following by means of the LED scheme.

Binding Energies of CH2 Adducts

The 1CH2···X and 3CH2···X binding energies are reported in Table . In the same table their decomposition into geometric preparation and interaction energy is also given. The latter is further decomposed into its reChemical">ference and correlation energy contributions.

Table 2

Calculated Equilibrium ΔE Binding Energies (kcal/mol) of the Studied CH2 Adducts and Their Decomposition into the Reference (QRO/aug-cc-pV5Z) and DLPNO-CCSD(T)/TightPNO/aug-cc-pV5Z Correlation Energies Together with the Contribution Δ of Each Term to the Singlet–Triplet Gap

			ΔE terms		ΔEint terms
molecule	state	ΔE	ΔEint	ΔEgeo-prep	ΔEintref	ΔEintC
CH2···H2O	T0	–1.56	–1.58	0.02	–0.30	–1.28
	S1	–5.27	–5.36	0.09	–3.96	–1.41
	Δ	–3.71	–3.78	0.07	–3.66	–0.13
CH2···He	T0	–0.02	–0.02	0.00	0.02	–0.04
	S1	–0.33	–0.33	0.00	0.19	–0.53
	Δ	–0.31	–0.31	0.00	0.17	–0.49
CH2···Ne	T0	–0.06	–0.06	0.00	0.06	–0.12
	S1	–0.35	–0.35	0.00	0.30	–0.65
	Δ	–0.29	–0.29	0.00	0.24	–0.53
CH2···Ar	T0	–0.20	–0.20	0.00	0.21	–0.41
	S1	–1.25	–1.25	0.00	1.04	–2.30
	Δ	–1.05	–1.05	0.00	0.83	–1.89
CH2···Kr	T0	–0.31	–0.31	0.00	0.31	–0.62
	S1	–2.19	–2.19	0.00	2.39	–4.58
	Δ	–1.88	–1.88	0.00	2.08	–3.96

Consistent with the trend of Δ previously discussed, Chemical">water and rare gases interact more strongly with 1CH2 than with 3CH2, thus resulting in a lowering of the <Chemical">span class="Chemical">singlet–triplet gap in CH2···X compounds (Δ < 0). The extent of the ΔEintref and ΔEintC contributions to the 1CH2···X and 3CH2···X interactions varies depending on the system. In particular, ΔEintref dominates the CH2···H2O interaction, while the CH2···Rg interaction is dominated by ΔEint. The ΔEintref and ΔEintC contributions of the 1,3CH2···X interactions are decomposed with the LED scheme into physically meaningful terms. The corresponding LED energy terms at the equilibrium geometries are reported in Table and discussed in the next section.

Table 3

LED of the Reference (QRO/aug-cc-pV5Z) and DLPNO-CCSD(T)/TightPNO/aug-cc-pV5Z Correlation Contributions to Interaction Energies (kcal/mol) for the Studied CH2 Adducts and the Contribution Δ of Each Term to the Singlet–Triplet Gap

		reference energy decomposition				correlation energy decomposition
molecule	state	ΔEintref	ΔEel-prepref	Eelstat	Eexch	ΔEintC	EdispC-CCSD	ΔEno-dispC-CCSD	ΔEintC-(T)
CH2···H2O	T0	–0.30	20.71	–15.88	–5.13	–1.28	–0.89	–0.23	–0.16	0.57
	S1	–3.96	29.71	–28.70	–4.96	–1.41	–1.22	0.04	–0.23	0.23
	Δ	–3.66	9.00	–12.82	0.17	–0.13	–0.33	0.27	–0.07	0.09
CH2···He	T0	0.02	0.09	–0.06	–0.03	–0.04	–0.05	0.01	0.00	2.50
	S1	0.19	1.39	–0.78	–0.43	–0.53	–0.44	–0.04	–0.05	1.33
	Δ	0.17	1.30	–0.72	–0.40	–0.49	–0.39	–0.05	–0.05	1.26
CH2···Ne	T0	0.06	0.69	–0.54	–0.13	–0.12	–0.12	0.01	–0.01	2.00
	S1	0.30	2.63	–1.85	–0.52	–0.65	–0.58	0.01	–0.08	1.66
	Δ	0.24	1.94	–1.31	–0.39	–0.53	–0.46	0.00	–0.07	1.59
CH2···Ar	T0	0.21	3.18	–2.36	–0.55	–0.41	–0.39	0.02	–0.04	1.95
	S1	1.04	11.69	–8.00	–2.60	–2.30	–1.65	–0.31	–0.34	1.32
	Δ	0.83	8.51	–5.64	–2.05	–1.89	–1.26	–0.33	–0.30	1.20
CH2···Kr	T0	0.31	1.68	–4.27	–0.88	–0.62	–0.48	–0.09	–0.05	1.55
	S1	2.39	27.78	–22.53	–6.63	–4.58	–2.67	–1.30	–0.61	1.22
	Δ	2.08	26.1	–18.26	–5.75	–3.96	–2.19	–1.21	–0.56	1.16

LED Analysis of the CH2···H2O Interaction

For the 3CH2···Chemical">H2O case, the sum of attractive electrostatic Eelstat and exchange Eexch terms of ΔEintref is almost entirely compensated by its repulsive static electronic preparation term ΔEel-prepref (see Table ). Thus, ΔEintref is practically negligible, demonstrating that dynamic electron correlation is reChemical">sponsible for the stability of the 3CH2···<Chemical">span class="Chemical">H2O adduct. In particular, London dispersion is one of the most important energy components of the interaction, as demonstrated by the large EdispC-CCSD/ΔE ratio of 0.57.

A difChemical">ferent picture emerges in 1CH2···<Chemical">span class="Chemical">H2O. In this case, the interacting species are closer (see Figure ) and ΔEintref consists of larger ΔEel-prepref and Eelstat values compared with those of 3CH2···H2O. These two contributions largely cancel each other. Hence, the overall ΔEintref is on the order of the remaining attractive exchange interaction.

Even though the correlation contribution to the Chemical">1,3CH2···<Chemical">span class="Chemical">H2O interaction is largely dominated by the London dispersion term EdispC-CCSD (see Table ) 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 CH2 while interacting with water is driven by the fact that 1CH2···H2O shows a much larger electrostatic interaction than 3CH2···H2O.

This efChemical">fect can be simply rationalized by looking at the schematic representation of the electronic configuration of <Chemical">span class="Chemical">singlet (Figure , top left) and triplet (Figure , top right) carbenes,[104] reported in Figure . At the bottom of the same figure, their molecular electrostatic potential (MEP)[105] maps projected onto the corresponding one electron densities at the reference level are also reported. The MEP map of the singlet methylene (Figure , bottom left) shows a negative electrostatic potential in the region of the lone pair. This favors H-bond interactions. Conversely, the MEP of the triplet carbene is much more isotropic, which leads to weaker H-bonding interactions.

Figure 3

Valence electron configurations and molecular electrostatic potential maps of 1CH2 (left) and 3CH2 (right) projected onto the corresponding molecular electron densities calculated at the UHF/aug-cc-pV5Z level in the unit of Eh/e. Red region identifies lowest electrostatic potential and thus highest electron density, while blue region identifies the opposite.

It is worth mentioning here that the singles term included in ΔEno-dispC-CCSD amounts to −2.4 kcal/mol for the 3CH2 fragment of all of the adducts studied in this work (see Table S5). In contrast, the adducts containing 1CH2 Chemical">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.

LED Analysis of the CH2···Rg Interaction

The interaction of CH2 with rare gases is relatively weak. In both Chemical">singlet and triplet states (see Tables and 3), the ΔEintref values are positive, meaning that the repulsive electronic preparation energies due to the distortion of electron clouds of both CH2 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 diChemical">spersion contribution (EdiChemical">spC-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., Edisp,L can be writtenwhere C6 is the atom pairwise induced Chemical">dipole–induced <Chemical">span class="Chemical">dipole interaction coefficient, IX and IY 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 CH2 and rare gases by comparing the London dispersion energies obtained from eq with those derived from the LED method. To do that we computed ionization potential and numerical polarizabilities of the isolated CH2 and rare <Chemical">span class="Disease">gas atoms at the DLPNO-CCSD(T)/aug-cc-pV5Z level. The computed values are reasonably close to the available experimental values (see Table ).[108−114] As an approximation we take r as the distance between the rare gas and the carbon atom (see Figure ) at the equilibrium geometry of CH2···Rg adducts.

Table 4

Numerical DLPNO-CCSD(T)/TightPNO/aug-cc-pV5Z Polarizabilities (α) and the Ionization Potentials (I) Compared with Available Experimental Data

	state	αcalc (Å3)	αexp. (Å3)[108−110]	Icalc (eV)	Iexp (eV)[111−114]
CH2	T0	2.186		10.393	10.386
CH2	S1	2.391		10.599
H2O	S0	1.449	1.456	12.753	12.621
He	S0	0.205	0.205	24.442	24.587
Ne	S0	0.393	0.395	21.572	21.565
Ar	S0	1.642	1.643	15.784	15.760
Kr	S0	2.475	2.486	14.200	14.000

The correlation between the dispersion energies obtained with London formula and the corresponding LED values for the CH2···Rg systems studied in this work is given in Figure . 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 R2 larger than 0.98 for both the 1CH2···Rg and the 3CH2···Rg series. In all cases, London dispersion increases with the increase of the polarizability of the rare gas atoms.

Figure 4

Plot of dispersion energies of CH2 interacting with rare gases: EdispC-CCSD of DLPNO-CCSD(T)/TightPNO/LED/aug-cc-pV5Z vs Edisp,L of London equation.

The 1CH2···Rg dispersion energies obtained by the London formula are very similar to the LED values, with variations of 13% in average. On the other hand, the London formula underestimates (−40%) the 3CH2···Rg dispersion interaction. The incorporation of the higher order terms seems thus necessary for open-shell systems in order to estimate dispersion energies more accurately.

Ferrous Heme···CO Interaction

Relative Spin State Energies of Heme Species

The calculated Chemical">singlet–triplet (1A1g – <Chemical">span class="Chemical">3Eg and 1A1g – 3A2g) and triplet–triplet (3Eg – 3A2g) energy gaps of the free and CO-bound heme are given in Table at various levels of theory. The corresponding Δ values are also reported in the same table.

Table 5

Calculated Relative Spin State Energies (Erel) of the Free and CO-Bound Heme Together with the Variation (Δ) of Erel Relative to the Free Heme at the Reference QRO/CBS, DLPNO-CCSD/TightPNO/CBS, and DLPNO-CCSD(T)/TightPNO/CBS Levels (in kcal/mol)a

	Erel	Δ
	DLPNO-CCSD(T)	DLPNO-CCSD(T)	DLPNO-CCSD	reference (QRO)
3Eg – 3A2ga
PFe	1.97	0.00	0.00	0.00
PFe···CO	–0.57	–2.54	–1.24	6.67
1A1g – 3Egb
PFe	32.13	0.00	0.00	0.00
PFe···CO	–6.50	–38.63	–31.78	9.10
1A1g – 3A2gb
PFe	34.10	0.00	0.00	0.00
PFe···CO	–7.07	–41.17	–33.02	15.77

A positive (negative) Erel implies that the 3A2g state is more (less) stable than the 3Eg, while the negative (positive) sign of Δ implies the increase (decrease) in the stability of the 3Eg state relative to the 3A2g state of PFe upon interacting with CO.

In this case, Erel corresponds to ES–T. A positive (negative) Erel implies that the triplet state is more (less) stable than the singlet state, while the positive (negative) sign of Δ implies the increase (decrease) in ES–T of PFe upon interacting with CO.

A positive (negative) Erel implies that the 3A2g state is more (less) stable than the Chemical">3Eg, while the negative (positive) sign of Δ implies the increase (decrease) in the stability of the <Chemical">span class="Chemical">3Eg state relative to the 3A2g state of PFe upon interacting with CO.

In this case, Erel corresponds to ES–T. A positive (negative) Erel implies that the triplet state is more (less) stable than the Chemical">singlet state, while the positive (negative) sign of Δ implies the increase (decrease) in ES–T of P<Chemical">span class="Chemical">Fe upon interacting with CO.

Experimental studies on bare Chemical">iron–<Chemical">span class="Chemical">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, Mössbauer data of iron(II)–tetraphenylporphyrln (FeTPP) were interpreted as an indication of either a 3Eg[116] or a 3A2g[70,117] ground state. Ligand field calculations that are consistent with the magnetic susceptibility measurements predict an 3A2g ground state.[118] Consistently, previous hybrid density functional, CASPT2, and several CI studies find the 3A2g state more stable than the 3Eg 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 3A2g state of heme to be only 0.96 and 1.97 kcal/mol more stable than the 3Eg state, respectively (see Table ).

The d orbital of Chemical">Fe(II), which is doubly occupied in the 3A2g state, is destabilized upon the binding of an axial ligand to the <Chemical">span class="Chemical">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 orbital raises in energy and becomes singly occupied. Thus, the ground state changes from triplet to quintet, and the lowest triplet changes from 3A2g to 3Eg.[89] As CO is a strong field ligand, the destabilization of the d 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 orbital.[67] Consequently, the identity of the most stable triplet state of CO-bound heme complexes is rarely studied. The present calculations find the 3Eg state of the CO-bound heme to be 2.16 and 0.57 kcal/mol more stable than 3A2g at the B3LYP-D3/def2-TZVP and DLPNO-CCSD(T)/TightPNO/CBS levels, respectively. Hence, CO binding probably reverses the energetic order of the 3Eg and 3A2g states, consistent with the imidazole-bound heme case. These changes can be rationalized in terms of the destabilization of the d orbital of Fe(II).

Unfortunately, no direct experimental measure of ES–T gaps exists for these complexes. At the B3LYP-D3/def2-TZVP and DLPNO-CCSD(T)/TightPNO/CBS levels, the ES–T gap of the free Chemical">heme is 33.32 and 32.13 kcal/mol relative to <Chemical">span class="Chemical">3Eg, while it is 34.88 and 34.10 kcal/mol relative to 3A2g (see Table ). 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 , upon interacting with CO, at the DLPNO-CCSD(T)/TightPNO/CBS level, the ES–T gap relative to the Chemical">3Eg and 3A2g states reduces to −6.50 and −7.07 kcal/mol (Δ = −38.63 and −41.17 kcal/mol), reChemical">spectively. Hence, CO binding significantly stabilizes the 1A1g 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 dif<Chemical">span class="Chemical">ferential spin state stabilization is reported in the following sections.

Binding Energy of the Heme Adducts

The calculated binding energies of PChemical">Fe(1A1g)···CO, P<Chemical">span class="Chemical">Fe(3Eg)···CO, and PFe(3A2g)···CO adducts at the DLPNO-CCSD(T)/TightPNO/CBS level are reported in Table . 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.

Table 6

Calculated Equilibrium ΔE Binding Energies (kcal/mol) of the CO-Bound Heme Adducts and Their Decomposition into the Reference (QRO/CBS) and DLPNO-CCSD(T)/TightPNO/CBS Correlation Energies Together with the Contribution Δ of Each Term to the Singlet–Triplet Energy Gap

		ΔE terms		ΔEint terms
	ΔE	ΔEint	ΔEgeo-prep	ΔEintref	ΔEintC
3Eg	–4.21	–5.94	1.74	7.31	–13.25
3A2g	–1.67	–1.77	0.10	1.06	–2.84
1A1g	–42.84	–46.28	3.44	12.28	–58.55
Δ1A1g–3Eg	–38.63	–40.34	1.70	4.97	–45.30
Δ1A1g–3A2g	–41.17	–44.51	3.34	11.22	–55.71

The PChemical">Fe(1A1g)···CO interaction is much stronger than the P<Chemical">span class="Chemical">Fe(3A2g/3Eg)···CO ones, as apparent from their large negative Δ values of −41.17/–38.63 kcal/mol (see Table ). This leads to the observed reduction of the singlet–triplet gap of heme derivatives shown in Table . However, at the reference level, the PFe(1A1g)···CO interaction is strongly repulsive (Δintref > 0). Thus, electron correlation counteracts this repulsion in the singlet state and makes the overall interaction significantly attractive.

On the technical side, PChemical">Fe(<Chemical">span class="Chemical">3Eg/3A2g)···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(1A1g)···CO adduct is more sensitive to the DLPNO thresholds and basis sets (see Table S8) adopted. For example, the PFe(1A1g)···CO binding energy is −42.84 kcal/mol with TightPNO/CBS and −36.90 kcal/mol with NormalPNO/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 PChemical">Fe(1A1g)···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 P<Chemical">span class="Chemical">Fe(3A2g) 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(4-sulfonatophenyl) 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(1A1g)–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(1A1g)···CO at the B3LYP-D3/def2-TZVP level of theory.

LED Analysis of the Interaction between Heme and CO

LED terms of the PChemical">Fe(1A1g/<Chemical">span class="Chemical">3Eg/3A2g)···CO interactions are given in Table 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 (rFe···C) was varied from ∼1.7 to ∼4.5 Å with an increment of 0.2 Å. The resulting energy profiles are reported in Figure together with their decomposition into the various LED terms.

Table 7

LED of the Reference (QRO/CBS) and DLPNO-CCSD(T)/TightPNO/CBS Correlation Interaction Energies (kcal/mol) for the CO-Bound Heme Adducts and the Contribution Δ of Each Term to the Singlet–Triplet Gap

	reference energy decomposition				correlation energy decomposition
	ΔEintref	ΔEel-prepref	Eelstat	Eexch	ΔEintC	EdispC-CCSD	ΔEno-dispC- CCSD	ΔEintC-(T)
3Eg	7.31	153.39	–121.54	–24.54	–13.25	–7.53	–4.11	–1.61	1.79
3A2g	1.06	12.94	–9.67	–2.21	–2.84	–2.61	0.07	–0.29	1.56
1A1g	12.28	714.55	–601.76	–100.51	–58.55	–22.56	–27.88	–8.11	0.53
Δ1A1g–3Eg	4.97	561.16	–480.22	–75.97	–46.30	–15.03	–23.77	–6.50	0.39
Δ1A1g–3A2g	11.22	701.61	–592.09	–98.3	–55.71	–19.95	–27.95	–7.82	0.48

Figure 5

Distance dependence of DLPNO-CCSD(T)/NormalPNO/LED/def2-TZVP terms of the interaction of CO with the 3Eg triplet (left), 3A2g triplet (middle), and 1A1g singlet (right) states of PFe at the B3LYP-D3/def2-TZVP geometries. Relaxed PES scans were performed at the B3LYP-D3/def2-TZVP level. Vertical dotted lines correspond to the B3LYP-D3 minima. Energy of the dissociated fragments is used as reference in all cases.

Distance dependence of DLPNO-CCSD(T)/NormalPNO/LED/def2-TZVP terms of the interaction of CO with the Chemical">3Eg triplet (left), 3A2g triplet (middle), and 1A1g <Chemical">span class="Chemical">singlet (right) states of PFe at the B3LYP-D3/def2-TZVP geometries. Relaxed PES scans were performed at the B3LYP-D3/def2-TZVP level. Vertical dotted lines correspond to the B3LYP-D3 minima. Energy of the dissociated fragments is used as reference in all cases.

As mentioned above, PChemical">Fe(3A2g) <Chemical">span class="Chemical">features a doubly occupied d orbital that points toward the CO lone pair located on the carbon atom, while in PFe(3Eg) the d orbital is singly occupied. For this reason, PFe(3A2g)···CO features a steeper repulsive wall than PFe(3Eg)···CO. This leads to a longer Fe–C bond distance in PFe(3A2g)···CO (3.27 Å) than in PFe(3Eg)···CO (2.29 Å). Accordingly, the PFe(3A2g)···CO interaction (−1.67 kcal/mol) is weaker than the PFe(3Eg)···CO one (−4.21 kcal/mol), and all of the LED terms of PFe(3A2g)···CO are smaller than those of PFe(3Eg)···CO at the B3LYP-D3 equilibrium geometry (see Tables and 7, and Tables S9–S11).

For both triplet states, the interaction between PChemical">Fe(<Chemical">span class="Chemical">3Eg/3A2g) and CO is repulsive at the QRO level (top panels of Figure and Table ) in the short range. Consistent with this picture, the LED decomposition of the reference PFe(3Eg/3A2g)···CO interaction energies (central panels of Figure and Table ) shows that the repulsive electronic preparation energies due to the distortion of the electron clouds of PFe(3Eg/3A2g) 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 PChemical">Fe(<Chemical">span class="Chemical">3Eg/3A2g)···CO adducts. Note that the B3LYP-D3 equilibrium geometries for the PFe(3Eg/3A2g)···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 and Table ) shows that the dispersion energy dominates the PChemical">Fe(<Chemical">span class="Chemical">3Eg/3A2g)···CO interaction. The corresponding EdispC-CCSD value amounts to −7.53/–2.61 kcal/mol at the B3LYP-D3 equilibrium geometry. In fact, both PFe(3Eg)···CO and PFe(3A2g)···CO adducts have a EdispC-CCSD/ΔE ratio larger than 1.5 (see Tables 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 PChemical">Fe(1A1g)···CO interaction is completely dif<Chemical">span class="Chemical">ferent. The energy profile computed at the QRO level (top right panel of Figure ) 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(1A1g)···CO adduct.

The LED decomposition of the reChemical">ference P<Chemical">span class="Chemical">Fe(1A1g)···CO interaction energy (central right panel of Figure ) shows that the electrostatic interaction is larger than the corresponding electronic preparation at the QRO minimum. Importantly, the PFe(1A1g)···CO interaction at the QRO level is associated with larger LED components for both repulsive and attractive terms than the PFe(3Eg/3A2g)···CO interaction at all Fe–C distances. Hence, the electronic clouds of the interacting fragments are distorted more significantly and interact stronger in PFe(1A1g)···CO than in PFe(3Eg/3A2g)···CO, indicating that a significant polarization of the fragments occurs upon CO binding already at the QRO level.

The analysis of the LED terms of the correlation energy reported in the bottom right panel of Figure is illuminating. Dispersive, nondispersive, and triples contribution are all significant in the short range. The corresponding TightPNO/CBS values at the B3LYP-D3 equilibrium geometry are −22.56, −27.88, and −8.11 kcal/mol for EdispC-CCSD, ΔEno-dispC-CCSD, and ΔEintC-(T), respectively. Hence, these three correlation terms are responsible for the significant strength of the PChemical">Fe(1A1g)···CO interaction, which in turn determines the reduction observed in the <Chemical">span class="Chemical">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 difChemical">ference function Δρ(x,y,z) describing the electron density rearrangement taking place upon CO binding for each <Chemical">span class="Gene">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 .

Figure 6

For (a) 3Eg, (b) 3A2g, and (c) 1A1g states, the contour plots of the reference QRO (up) and the DLPNO-CCSD/NormalPNO/def2-TZVP correlation (bottom) contributions to the electron density rearrangements taking place upon CO binding. Plots are all given with the density isosurface contour value of ±0.002 e/Bohr3. Blue and red surfaces identify regions of electron density accumulation and depletion, respectively.

For (a) Chemical">3Eg, (b) 3A2g, and (c) 1A1g states, the contour plots of the re<Chemical">span class="Chemical">ference QRO (up) and the DLPNO-CCSD/NormalPNO/def2-TZVP correlation (bottom) contributions to the electron density rearrangements taking place upon CO binding. Plots are all given with the density isosurface contour value of ±0.002 e/Bohr3. Blue and red surfaces identify regions of electron density accumulation and depletion, respectively.

Consistent with the fact that the LED describes the PChemical">Fe(<Chemical">span class="Chemical">3Eg/3A2g)···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(1A1g)···CO interaction is accompanied by a significant charge accumulation in the region of the bond, consistent with the fact that the PFe(1A1g)···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(1A1g)···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 . The Δq values are negligible for all d orbitals in PChemical">Fe(<Chemical">span class="Chemical">3Eg/3A2g)–CO adducts, while they are significant in the PFe(1A1g)–CO adduct. In particular, the population of the formally empty d orbital of the singlet state increases significantly upon CO binding. At the same time, the population of the formally doubly occupied d and d 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 (DCD) bonding model, which was originally introduced to discuss the binding of olefins to transition metals.[123]

Table 8

UHF/def2-TZVP and DLPNO-CCSD/NormalPNO/def2-TZVP Mulliken d-Orbital Populations (q in Units of e) for the 3Eg, 3A2g, and 1A1g States of PFe and PFe···CO Together with Their Variations Upon CO Binding (Δq)

	UHF			DLPNO-CCSD
	qPFe	qPFe···CO	Δq	qPFe	qPFe···CO	Δq
3Eg
dz2	1.15	1.15	0.00	1.13	1.14	0.01
dxz	1.02	1.02	0.00	1.03	1.03	0.00
dyz	2.00	2.00	0.00	1.98	1.98	0.00
dx2–iy2	1.89	1.87	–0.02	1.91	1.87	–0.04
dxy	0.36	0.36	0.00	0.48	0.48	0.00
3A2g
dz2	1.96	1.97	0.01	1.93	1.94	0.01
dxz	1.02	1.02	0.00	1.03	1.02	–0.01
dyz	1.02	1.02	0.00	1.03	1.02	–0.01
dx2–y2	2.00	2.00	0.00	1.99	1.99	0.00
dxy	0.36	0.36	0.00	0.49	0.49	0.00
1A1g
dz2	0.07	0.22	0.15	0.11	0.40	0.29
dxz	2.01	1.93	–0.08	1.99	1.84	–0.15
dyz	2.01	1.93	–0.08	1.99	1.84	–0.15
dx2–y2	2.00	2.02	0.02	1.99	2.01	0.02
dxy	0.32	0.35	0.03	0.45	0.50	0.05

A significant σ-donation of charge takes place from the lone pair located on the CO ligand to the empty d orbital located on the Chemical">Fe(II) of the <Chemical">span class="Chemical">singlet adduct, while at the same time π-backdonation occurs from the d and d 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(1A1g)···CO adduct (1.141 Å) with respect to that of free CO (1.125 Å, which is the same as in the PFe(3Eg/3A2g)···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 . Note that electron correlation significantly affects the magnitude of the DCD components, consistent with our previous findings on agostic complexes.[41]

Scheme 2

Schematic Diagram of the σ-Donation and π-Backdonation Taking Place upon CO Binding to the (a) 3Eg, (b) 3A2g, and (c) 1A1g States of PFe

The corresponding amount of electron transfer (Δq) based on Mulliken population analysis at the DLPNO-CCSD/NormalPNO/def2-TZVP and HF (in parentheses) levels is also reported.

The corresponding amount of electron transChemical">fer (Δq) based on Mulliken population analysis at the DLPNO-CCSD/NormalPNO/def2-TZVP and HF (in parentheses) levels is also reported.

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 scChemical">heme to open-shell molecular systems that was implemented in the ORCA program package. As a first illustrative case study this sc<Chemical">span class="Chemical">heme was applied to investigate the mechanism that governs the change in the singlet–triplet energy gap of CH2 upon interaction with water and rare gases (Rg) and of heme (PFe) upon interaction with CO.

The CH2···X interaction (X = Chemical">H2O, He, Ne, Ar, and Kr) was found to be attractive for both the triplet 3CH2 and the <Chemical">span class="Chemical">singlet 1CH2 methylene. The interaction is stronger with 1CH2 than with 3CH2, resulting in a lowering of the singlet–triplet energy gap (ES–T). The LED analysis of the CH2···H2O interaction showed that electrostatics dominates the interaction for both spin states of methylene. The lowering of ES–T can thus be understood in terms of the larger electrostatic interaction in 1CH2···H2O than that in 3CH2···H2O, 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 ES–T arises from the stronger London dispersion in 1CH2···Rg than that in 3CH2···Rg. This is consistent with the larger polarizability of 1CH2 than that of 3CH2.

As regards Chemical">heme systems, it was found that the P<Chemical">span class="Chemical">Fe···CO interaction is dominated by the correlation contribution for all the spin states of PFe investigated in this work (1A1g, 3Eg, and 3A2g). Although PFe is a triplet in its ground state, the PFe(1A1g)···CO interaction is significantly stronger than the interaction of PFe in the two lowest-lying triplets (3Eg and 3A2g) with CO, leading to a lowering of ES–T in the adduct. In fact, the LED analysis demonstrated that PFe(3Eg)···CO and the PFe(3A2g)···CO adducts can be described as van der Waals complexes stabilized by weak dispersion forces. In contrast, the PFe(1A1g)···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 orbital on the central iron atom (Fe(II) ← CO σ-donation) and from the filled d and d orbitals into the π* antibonding orbitals on the CO molecule (Fe(II) → CO π-backdonation).

It should be stressed here again that the results just discussed reChemical">fer to the minimal <Chemical">span class="Chemical">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.