Understanding the mechanism stabilizing intermediate spin states in Fe(II)-Porphyrin
UUnderstanding the mechanism stabilizingintermediate spin states in Fe(II)-Porphyrin
Giovanni Li Manni ∗ and Ali Alavi ∗ Max-Planck-Institut für Festkörperforschung, Heisenbergstraße 1, 70569 Stuttgart,Germany
E-mail: [email protected]; [email protected]
Abstract
Spin fluctuations in Fe(II)-porphyrins are at the heart of heme-proteins functional-ity. Despite significant progress in porphyrin chemistry, the mechanisms that rule spinstate stabilization remain elusive. Here, it is demonstrated by using multiconfigura-tional quantum chemical approaches, including the novel Stochastic-CASSCF method,that electron delocalization between the metal centre and the π system of the macrocy-cle differentially stabilizes the triplet spin states over the quintet. This delocalizationtakes place via charge-transfer excitations, involving the out-of-plane iron d orbitals,key linking orbitals between metal and macrocycle. Through a correlated breathingmechanism the 3d electrons can make transitions towards the π orbitals of the macro-cycle. This guarantees a strong coupling between the on-site radial correlation on themetal and electron delocalization. Opposite-spin 3d electrons of the triplet can effec-tively reduce electron repulsion in this manner. Constraining the out-of-plane orbitalsfrom breathing hinders delocalization and reverses the spin ordering. Our results finda qualitative analogue in Kekulé resonance structures involving also the metal centre. a r X i v : . [ phy s i c s . c h e m - ph ] D ec Introduction
Metal-porphyrins are versatile chemical species which biological systems make abundant useof, with Mg(II)-porphyrins and Fe(II)-porphyrins being the most common in nature. Thelatter are used in a number of vital functions in aerobic life, including dioxygen transport andreduction. From an electronic point of view Mg(II) porphyrins are closed shell diamagneticcompounds, while Fe(II)-porphyrins with a d configuration at the metal centre may showa multitude of low-lying electronic states. The high-spin (quintet), the intermediate-spin(triplet) and the low-spin (singlet) states are near degenerate and, depending on the coordi-nating ligands, geometry and thermodynamical conditions, their relative energy may easilychange. Spin changes are the key feature that enables enzymatic and biomimetic reactionsinvolving Fe-porphyrins. Molecular and electron transport as well as metabolic reactionstake place thanks to the facile spin and oxidation state changes of these compounds. Theoxidative oxygenation (insertion of one oxygen in a C–H bond) by the cytochrome P450s isan example. In this reaction, the Fe(II)-porphyrin represents the active species that bindsmolecular oxygen and weakens its bond, forming the actual oxo-species that proceeds to theoxygenation. It is then crucial from a mechanistic point of view to understand the electronicstructure of Fe(II)-porphyrins and the main elements that stabilize one spin state over theothers.The first ab initio calculations on the free-base porphyrin, that can be regarded as theparent compound for such systems, were done by Almlöf. Later numerous density functionaltheory calculations were carried out in order to explore the electronic properties of metalporphyrins.
In spite of a large amount of experimental and theoretical data, many ques-tions are still unanswered regarding their electronic properties and reactivity. For instance,a definitive assignment of the ground state of four-coordinated ferrous porphyrins is stillmissing and the main ingredients governing ground state electron configuration unknownto date. A A g ground state with configuration ( d x − y ) ( d xz , d yz ) ( d z ) ( d xy refers to theanti-bonding orbital pointing at the N atoms) was suggested by Mössbauer, magnetic measurements of the Fe(II)-tetraphenilporphyrin (FeTPP). A E g statewas suggested by Raman spectroscopy, with configuration ( d x − y ) ( d xz , d yz ) ( d z ) forthe Fe(II)-octaethylporphyrin (FeOEP). A high-spin state was reported for the octamethyl-tetrabenzporphyrin-iron(II). Many factors, such as functionalization of the macrocycle andsolvation, may affect the relative ordering of the low-lying states.A large number of theoretical studies are available for model systems of the ferrousporphyrin.
Density Functional Theory (DFT) predicts a triplet ground state for theFe(II)-porphyrin, although still consensus has yet to be reached for the specific symmetry ofthe state. Swart et al. found a high sensitivity of spin gaps on the type of functional usedin density functional approximations. The OPBE functional predicted a E g groundstate for the ferrous porphyrin, with the A g and the A g at 4.0 and 7.2 kcal/mol aboverespectively. Interestingly, other functionals (BP86 and B3LYP) predicted a A g groundstate with the E g at 1.8 and 6.2 kcal/mol above respectively. A completely different scenario is depicted by wave-function theory based approaches,including high-level methodologies generally regarded as being more reliable than DFT.The Restricted Open-shell Hartree-Fock method, ROHF, predicts the high-spin A g groundstate. The gap between the quintet and the triplet spin states shrinks when post-SCF meth-ods are used, however, the spin-ordering remains in favor of the high-spin state, suggesting asystematic error in the theoretical framework. Pierloot has extensively studied these systemsby multiconfigurational methods and how spin gap predictions depend on the choice of theactive space. A definitive argument on the mechanism stabilizing the triplet spin statewas not suggested. Transition metal spin chemistry has always been challenging for quan-tum chemical methods and a simple and reliable theoretical approach for spin-dependentproperties is still not available. In addition to the lack of consensus among the theoreticalmethodologies, a more fundamental question is still unanswered, what factors dictate therelative stabilization of the competing spin states?
In this report we will show in great detail the mechanisms that rule spin gaps in the bare3errous porphyrin, by analyzing the six low-lying states, B g , A g , E g , A g , B g and E g ,of the Fe(II)-porphyrin (Figure 1). The Stochastic Complete Active Space Self-ConsistentFigure 1: Dominant electron configuration of the six competing spin states of Fe(II)-porphyrin.Field method, Stochastic-CASSCF, is the method of choice for this investigation. TheCASSCF represents a simple and natural way to systematically probe correlation channelsin correlated molecular systems. The core concept of CASSCF is the active space, a list of“critical” orbitals with their electrons for which a complete many-body expansion is generatedand orbitals are variationally optimized under the field generated by the multiconfigurationalwave-function, removing any bias related to the choice of the trial orbitals. Equipped withthe Stochastic-CASSCF method we have been able to identify the main correlation effectsthat stabilize the intermediate spin-state over the quintet spin-state and have been able toestablish that important communication pathways between the aromatic macrocycle andthe metal centre exist only for the low-lying triplet states. Qualitatively these correlationchannels can be described as Kekulé resonance structures involving also the metal centre.4 Results
Energy splittings at the various levels of theory and VTZP basis set are summarized inFigure 2 (additional results available in the Supplementary Material). Energy gaps betweenthe quintet spin states are rather insensitive to the choice of the active space, basis set anddynamic correlation correction via second order perturbation theory, CASPT2.
Whenturning our attention to the triplet spin states the scenario is rather different. A strong ’background.dat’30.0-5.00.05.010.015.020.025.0 B
2g 5 E g 3 B
1g 3 A
2g 3 E g Spin States59.063.068.0 -3.1 R e l a ti v e E n e r gy [ k ca l/ m o l ] CAS(6,5)CAS(8,6)CAS(8,12)CAS(14,16) PT2(6,5)PT2(8,6)PT2(8,12)PT2(14,16) Stoch-CAS(32,34)/500MStoch-CAS(32,34)/1BRAS(32,2,2;12,6,16)
Figure 2: Energy splittings relative to the A g state in VTZP basis set. The upper reddashed line marks the lowest relative energy value that one could obtain with conventionalCASPT2(14,16). The lower red dashed lines marks the lowest relative energy value for theStochastic-CASSCF(32,34) method.dependence of the relative energy with respect to the active space and method of choice isobserved. CASPT2 spin gaps are affected both by the under-lying active space and basisset of choice. Enlarging the active space reduces the gap between the triplet states andthe A g state. Perturbative correction to the second order for any choice of active space5pproximately halves the triplet-quintet gap. The smallest CASSCF(6,5) and CASSCF(8,6)place the lowest triplet, E g , at ∼
18 kcal/mol above the A g state. When perturbativecorrection is added the A g becomes the lowest triplet state, lying at ∼ E g at 0 . A g ground state (3.2 kcal/mol with VDZP basis set). These results clearly showthat conventional CASPT2 is still not converged and the triplet might ( and will ) be furtherstabilized by more accurate methods, namely a larger under-lying active space.The most notable result summarized in Figure 2 is the energy splitting obtained by thelarge CASSCF(32,34). At this level of theory both the E g ( − . A g ( − . A g state, setting the E g as the ground state for thissystem. The A g is only 0 . E g . The energy lowering obtained by thislarge active space is substantial.Surprisingly, the restricted active space approach, RASSCF(32,34) sets the triplet statesagain above the A g state (2.9 and 4.1 kcal/mol for the E g and the A g state respectively).Truncation of the excitation level to only single and double excitations from RAS1 and toRAS3 obviously has a strong unfavorable impact on the triplet-quintet splittings. Within theStochastic-CASSCF approach the impact of the target number of walkers is nearly negligible.Variations of less than 0.1 kcal/mol were observed by enlarging the walker population from500 million to 1 billion (before convergence is reached, higher accuracy is expected for alarger walker population). The large CASSCF calculations seem to be able to circumvent thelimitations of smaller active spaces, truncated CI expansions (RAS case) and perturbatively-corrected results. Active space size limitations are substantial and RAS type of truncationsdo not represent a solution. CASPT2 energy estimates with small active space referencewave functions are not to be considered reliable for this class of compounds.From the inspection of the CASSCF(32,34) natural orbitals for the E g state a strongmixing of the 3 d yz orbital (and 3 d xz for the degenerate state) with π orbitals of the macrocycleare observed (see orbitals π y and π Hy of Figure 3 as an example). This mixing is missing in6igure 3: Stochastic-CAS(32,34) active natural orbitals for the E g state of the Fe(II)-porphyrin model system and their occupation numbers.7he A g state (natural orbitals for the quintet spin state are depicted in the SupplementaryMaterial). The mixing of the doubly occupied 3 d yz and the symmetry allowed π -orbitals forthe triplet state is not fortuitous. It means that off-diagonal elements in the one-body densitymatrix are quite large and the eigen-vectors leading to the natural orbitals will have largecontributions from the metal centre and the macrocycle. Large off-diagonal matrix elementscan arise only when the wave function is multiconfigurational, with a strong entanglementbetween orbitals related to those large off-diagonal matrix elements and, with cumulativelylarge amplitudes for charge-transfer (CT) determinants (CT determinants are those wherestarting from a given determinant electrons are excited from the occupied 3 d orbitals tothe empty π ∗ orbitals and from the occupied π orbitals to empty 3 d orbitals). Therefore,in the E g state, π orbitals are strongly correlated to the 3 d yz orbital (and 3 d xz ) via CTdeterminants. This property, that is not present in the A g state, is the driving force thatstabilizes the triplet over the quintet spin state.Adding the double-shell d orbitals into the active space allows the triplet state to dif-ferentially reduce on-site electron repulsion by exciting electrons out of the doubly occupied3 d orbitals, into the d shell. This effect is referred to as “ radial correlation ” or “ breath-ing ”. However, the d orbitals also open a “correlation pathway”, a communication channelbetween the π orbitals and the 3d orbitals (see Figure 4). Via the correlating d orbitalsvalence out-of-plane electrons are able to breathe out, expand towards the π orbitals of themacrocycle. This guarantees a strong coupling between the on-site radial correlation on themetal centre and electron delocalisation into the macrocycle. Constraining the out-of-planeorbitals from breathing hinders delocalization and leads to the reversal of the spin order-ing. This confirms the crucial role of the breathing effect in the spin chemistry of thesesystems. In the CAS(14,16) four π orbitals of the macrocycle are added into the activespace, π Hx , π Hy , π ∗ Lx and π ∗ Ly ( x and y subscripts refer to the symmetry of these orbitals)and are explicitly correlated to the 3d and d orbitals at the metal centre. They are the twodegenerate HOMO ( H subscript) and two degenerate LUMO ( L subscript) orbitals of the8igure 4: The correlating d orbitals provide a correlation pathway for connecting 3 d and π orbitals in the iron complex.free-base porphyrin. Including these π orbitals in the active space in absence of the double-shell orbitals, leads to less populated anti-bonding π orbitals. This difference is entirely dueto the d double-shell. This is a very interesting finding. Correlating d and π orbitals syn-ergistically favour electron delocalization and reduce on-site electron repulsion at the metalcentre. This synergic effect arises from the coupling through the Hamiltonian operator ofdeterminants of the type | d → d i , | π H → d i , | d → π ∗ L i , | π H → π ∗ L i and | π ∗ L → d i . Thistype of excitations are somehow included, although only up to the second order both in theRASSCF(32,34) and the CASPT2(8,12). As neither CASPT2(8,12) or RAS(32,34) providedconverged energetics, we conclude that it is not sufficient to correlate these orbitals solelyvia singly- and doubly-substituted determinants. Higher order excitations are responsiblefor the stabilization of the triplet over the quintet spin state. Also, these excitations havebeen included in the CAS(14,16), but there full delocalization is not explicitely accountedfor in the method (most of the π orbitals are not in the active space) and as a consequencethe triplet lies still above the quintet spin state. The important excitations, however, havebeen included in the large CASSCF(32,34) calculations.Inclusion of the entire π -system in the active space and the complete CI expansion of9he CASSCF(32,34) reveal another important result. The π non-bonding orbitals, π NB and π NB , have rather low occupation numbers, 1.90 and 1.72 respectively, both for the E g andthe A g states. At the same time π ∗ Lx and π ∗ Ly orbitals reach relatively high occupationnumbers (see Figure 5 and Figure 6). σ Ν σ Ν σ Ν σ M L π π π π π π H y π H x π N B π N B x - y y z x z z xy π * L x π * L y π * π * π * π * π * s x y z ’ x - y ’ y z ’ x z ’ z ’ xy Natural Orbitals O cc up a ti on N u m b e r CAS(6,5)CAS(8,6)CAS(8,12)CAS(12,13)CAS(14,16)RAS(32,34)(32,34)/500M
Figure 5: Natural orbitals occupation numbers for the E g state within the VTZP basis setfor different choices of active spaces. Bars are related to the minimum CAS(6,5) active spacethat can be considered the reference for larger active spaces.These results sensibly differ from the RAS(32,34) results and demonstrates that theCAS(32,34) wave function is substantially different from the other wave functions here an-alyzed, with some of the π orbitals heavily involved in the correlation for this system. Thelarge CASSCF(32,34) provides a more accurate description of the electron delocalization(conjugation and aromaticity) of the macrocycle. It also better accounts for CT determi-nants. As a consequence a more realistic field around the metal centre is created that,in response, differentially stabilizes the triplet over the quintet spin state. This effect is10 .000.050.100.150.200.25 σ Ν σ Ν σ Ν σ M L π π π π π π H y π H x π N B π N B z x - y y z x z xy π * L x π * L y π * π * π * π * π * s x y z ’ x - y ’ y z ’ x z ’ z ’ xy Natural Orbitals O cc up a ti on N u m b e r CAS(6,5)CAS(8,6)CAS(8,12)CAS(12,13)CAS(14,16)RAS(32,34)(32,34)/500M
Figure 6: Natural orbitals occupation numbers for the A g state within the VTZP basis setfor different choices of active space. Bars are related to the minimum CAS(6,5).11chieved only when the wave function is relaxed with respect to CI and orbital parametersand, contains higher order excitations coupled via the true Hamiltonian operator.A deeper analysis of the large Stochastic-CAS(32,34) wave function corroborates theprevious findings. Ligand-to-metal charge-transfer (LMCT) excitations, π → (3 d xz , d yz ),contribute for ∼
1% for the E g state. Numerous LMCT excitations of the type π → ( d xz , d yz ), and metal-to-ligand charge-transfer (MLCT) excitations of the type ( d xz , d yz ) → π ∗ and (3 d xz , d yz ) → π ∗ , contribute to the CAS(32,34) wave function of the E g statewith weights around 0.1-0.5%. CT determinants were already observed and reported in ourprevious work. However, in the present work π → ( d xz , d yz ) and ( d xz , d yz ) → π ∗ are alsoobserved. This finding reinforces the above statement, double-shell d orbitals contribute tothe “radial correlation” and also serve as a communication channel bridging the gap betweenmetal centre and macrocycle orbitals. They are actively involved in the delocalization of thevalence electrons. Valence out-of-plane 3d electrons expand with a breathing mechanism viathe correlating d orbitals, which have a larger overlap with the π orbitals of the macrocycle,and thus delocalize into the π system.Charge-transfer excitations are rare for the quintet spin state. The LMCT π → d xz,yz excitations contribute for only 0.3% to the wave function. MLCT excitations are even lessrepresentative of the wave function, with the 3 d xz,yz → π ∗ contributing for only 0.05%. The E g state of the ferrous porphyrin is characterized by an important interaction between π electrons and valence electrons at the metal centre. This feature is unique to porphyrinoidshosting transition metal centres. In fact, the Mg(II) porphyrin reported in our previous workdid not show any interaction (via CT excitations) between the magnesium orbitals and the π − π ∗ system of the aromatic macrocycle. In the Mg(II) compound doubly occupied orbitalsof the magnesium lie low in the energy spectrum and the virtual manifold too high in energyto mix with the orbitals of the π − π ∗ system. In the Fe(II)-compound 3d orbitals mix wellboth in terms of symmetry/overlap (considering the role of the d shell) and energy.12 Discussion
It is demonstrated by means of multiconfigurational approaches that for the bare Fe(II)-porphyrin the triplet spin state is stabilized over the quintet spin via metal-to-ligand andligand-to-metal charge-transfer excitations. These are numerous in the intermediate spinstate while being extremely rare in the quintet spin state.Previous quantum chemical simulation of the same CASSCF type, but based on smalleractive spaces, do not describe adequately the above mentioned charge-transfer configurations(they have too small an amplitude or do not exist at all). None of the smaller correlatedmethods here discussed shows 3 d xz,yz / π orbital mixing. Simply stated previous methods donot show electron delocalization between metal centre and macrocycle and, as a consequencethe high spin states are overstabilized. Previous multiconfigurational methods fail in cap-turing simultaneously ring correlation (aromaticity), correlation at the metal centre (radialcorrelation) and the interaction between metal centre and macrocycle via charge-transferexcitations which differentially stabilise the intermediate spin states over the high spin ones.The interaction between metal and macrocycle is rather complicated. In the largeCAS(32,34) wave function, we observe non-negligible charge-transfer electron configurationscoupling directly π and 3 d orbitals. Higher order interactions between macrocycle and metalcentre are also present via the double-shell d orbitals. The double-shell orbitals play a dualrole, they account for radial correlation at the metal centre and build a bridging pathwaybetween π orbitals at the macrocycle and valence orbitals at the metal centre. A breathingmechanism can be invoked for rationalize the stabilization of the triplet spin state over thequintet. Via the correlating d orbitals, which provide the necessary breathing mechanism,the valence electrons can more easily delocalize into the π system, and to a far greater extentthan the regular (non-breathing) 3d orbitals would allow.The charge-transfer configurations can be qualitatively described as Kekulé aromaticresonance structures, involving movement of the π orbitals of the macrocycle as well as thevalence electrons on the 3 d xz and 3 d yz orbitals of the metal centre as described in Figure 7.13o the best of our knowledge this is the first quantum chemical investigation that explicitlyFigure 7: Kekulé resonance structures involving movement of iron valence electrons of the3 d xz and 3 d yz orbitals.shows how the electron delocalization between metal centre and macrocycle is the key featurefor the stabilization of a spin state over others.This finding uncovers a new facet of metal-porphyrin chemistry and, to some extent, waysto control it. Our insights on the ferrous porphyrin may open new possibilities to manipulateelectron delocalization, for instance via chemical functionalization of the macrocycle and havecontrol over spin. The role of peripheral functional groups in functionalized metal-porphyrincan be related to the proposed mechanism and spin control can be achieved by ad hoc groups that enhance or hinder electron delocalization. Although the investigation focusedexclusively on the ferrous porphyrin, our manuscript creates a paradigm that could extend toother transition metal porphyrins and therefore have a broader impact. We predict that thesame delocalization mechanism is responsible for the ability of metal-porphyrins to undergoreduction and oxidation with ease in redox processes of living systems. To date the nature ofthe oxidized or reduced species is poorly understood, and so is the character of the acceptororbital in reduction reactions. According to the proposed mechanism, delocalization of theadditional electrons (in the reduced form) is expected. This delocalization of the extra chargereduces on-site repulsion making the event favorable.14 Acknowledgments
Roland Lindh and Ignacio Férnandez Galvan are acknowledged for their assistence in theevaluation of AO integrals within the Molcas chemistry package, via the Cholesky decom-position procedure. Professors Marcel Swart and Kristine Pierloot are acknowledged for theinspiring discussions on the relevance of accurate computational tools for spin chemistryinvestigations. The Max-Planck Society and the COST Action CM1305 (ECOSTBio) areacknowledged for their financial support.
The authors have contributed equally in the investigation, conceiving ideas, developing therequired theory and algorithm and planning and performing calculations. Both authors havediscussed the results and contributed to the final manuscript.
The authors declare no competing financial interest.
Correspondence and requests for materials should be addressed to [email protected] and [email protected] . Supplementary materials
The Supplementary Material contains details on: Model System. Basis Set and electronrepulsion integrals. Complete Active Space choices. uantum Monte Carlo setup. Past andPresent proposed active spaces. Tables S1: Active Orbitals in D2h point group. Energy15plittings with VDZP basis set. Figure S1: VDZP energy gaps. Listing 1: CartesianCoordinates. Figure S2: Natura orbitals of the quintet spin state. Mixing of d and π orbitals. Correlating d’ orbitals. Wave function analysis. Listing 2: Molcas input example.Listing 3: QMC input example. References cited within the Supporting Material. References (1) Sono, M.; Roach, M. P.; Coulter, E. D.; Dawson, J. H.
Chem. Rev. , , 2841–2888.(2) Meunier, B.; de Visser, S. P.; Shaik, S. Chem. Rev. , , 3947–3980.(3) Schöneboom, J. C.; Lin, H.; Reuter, N.; Thiel, W.; Cohen, S.; Ogliaro, F.; Shaik, S. J. Am. Chem. Soc. , , 8142–8151.(4) Almlöf, J. Int. J. Quantum Chem. , , 915–924.(5) Jones, D. H.; Hinman, A. S.; Ziegler, T. Inorg. Chem. , , 2092–2095.(6) Rosa, A.; Baerends, E. J. Inorg. Chem. , , 5637–5639.(7) Collman, J. P.; Hoard, J. L.; Kim, N.; Lang, G.; Reed, C. A. J. Am. Chem. Soc. , , 2676–2681.(8) Lang, G.; Spartalian, K.; Reed, C. A.; Collman, J. P. J. Chem. Phys. , ,5424–5427.(9) Boyd, P. D. W.; Buckingham, A. D.; McMeeking, R. F.; Mitra, S. Inorg. Chem. , , 3585–3591.(10) Goff, H.; La Mar, G. N.; Reed, C. A. J. Am. Chem. Soc. , , 3641–3647.(11) Mispelter, J.; Momenteau, M.; Lhoste, J. M. J. Chem. Phys. , , 1003–1012.1612) Kitagawa, T.; Teraoka, J. Chem. Phys. Lett. , , 443–445.(13) Sams, J. R.; Tsin, T. B. Chem. Phys. Lett. , , 599–601.(14) Matsuzawa, N.; Ata, M.; Dixon, D. A. J. Phys. Chem. , , 7698–7706.(15) Kozlowski, P. M.; Spiro, T. G.; Bérces, A.; Zgierski, M. Z. J. Phys. Chem. B , , 2603–2608.(16) Liao, M.-S.; Scheiner, S. J. Chem. Phys. , , 3635–3645.(17) Smith, D. M. A.; Dupuis, M.; Straatsma, T. P. Mol. Phys. , , 273–278.(18) Chen, H.; Lai, W.; Shaik, S. J. Phys. Chem. B , , 1727–1742.(19) Choe, Y.-K.; Nakajima, T.; Hirao, K.; Lindh, R. J. Chem. Phys. , , 3837–3844.(20) Choe, Y.-K.; Hashimoto, T.; Nakano, H.; Hirao, K. Chem. Phys. Lett. , ,380–388.(21) Vancoillie, S.; Zhao, H.; Tran, V. T.; Hendrickx, M. F. A.; Pierloot, K. J. Chem.Theory Comput. , , 3961–3977.(22) Pierloot, K.; Phung, Q. M.; Domingo, A. J. Chem. Theory Comput. , , 537–553.(23) Phung, Q. M.; Wouters, S.; Pierloot, K. J. Chem. Theory Comput. , , 4352–4361.(24) Radoń, M.; Pierloot, K. J. Phys. Chem. A , , 11824–11832.(25) Vancoillie, S.; Zhao, H.; Radoń, M.; Pierloot, K. J. Chem. Theory Comput. , ,576–582.(26) Radoń, M.; Broclawik, E.; Pierloot, K. J. Phys. Chem. B , , 1518–1528.1727) Radoń, M. J. Chem. Theory Comput. , , 2306–2321.(28) Swart, M.; Groenhof, A. R.; Ehlers, A. W.; Lammertsma, K. J. Phys. Chem. A , , 5479–5483.(29) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. , , 3865–3868.(30) Handy, N. C.; Cohen, A. J. Mol. Phys. , , 403–412.(31) Swart, M. J. Chem. Theory Comput. , , 2057–2066.(32) Swart, M. Int. J. Quantum Chem. , , 2–7.(33) Costas, M.; Harvey, J. N. Nat. Chem. , , 7–9.(34) Song, S.; Kim, M.-C.; Sim, E.; Benali, A.; Heinonen, O.; Burke, K. arXiv:1708.08425 (35) Li Manni, G.; Smart, S. D.; Alavi, A. J. Chem. Theory Comput. , , 1245–1258.(36) Roos, B. O.; Taylor, P. R.; Siegbahn, P. E. M. Chem. Phys. , , 157–173.(37) Roos, B. O. Int. J. Quantum Chem. , , 175–189.(38) Siegbahn, P. E. M.; Heiberg, A.; Roos, B. O.; Levy, B. Phys. Scr. , , 323–327.(39) Siegbahn, P. E. M.; Almlöf, J.; Heiberg, A.; Roos, B. O. J. Chem. Phys. , ,2384–2396.(40) Roos, B. O. Advances in Chemical Physics ; 2007; pp 399–445.(41) Andersson, K.; Malmqvist, P.-Å.; Roos, B. O.; Sadlej, A. J.; Wolinski, K.
J. Phys.Chem. , , 5483–5488.(42) Andersson, K.; Malmqvist, P.-Å.; Roos, B. O. J. Chem. Phys. , , 1218–1226.(43) Andersson, K.; Roos, B. O. Int. J. Quantum Chem. , , 591–607.1844) Finley, J.; Malmqvist, P.-Å.; Roos, B. O.; Serrano-Andrés, L. Chem. Phys. Lett. , , 299–306.(45) Ghigo, G.; Roos, B. O.; Malmqvist, P.-Å. Chem. Phys. Lett. , , 142–149.(46) Pierloot, K. Mol. Phys. , , 2083–2094.(47) Pierloot, K.; Vancoillie, S. J. Chem. Phys. , , 034104.(48) Pierloot, K.; Vancoillie, S. J. Chem. Phys. , , 124303.19 Supplementary Material.
This PDF file includes: • Model System.• Basis Set and electron repulsion integrals.• Details on the Complete Active Space choice.• Details on the Quantum Monte Carlo setup.• Past and Present proposed active spaces.• Tables S1: Active Orbitals in D2h point group.• Energy Splittings with VDZP basis set.• Figure S1: VDZP energy gaps.• Listing 1: Cartesian Coordinates.• Figure S2: Natura orbitals of the quintet spin state.• Mixing of d and π orbitals.• Correlating d’ orbitals.• Wave function analysis.• Listing 2: Molcas input example.• Listing 3: QMC input example.• References cited within the Supporting Material20 odel System. Our model system for the Fe(II)-porphyrin was derived from Pierloot’sstudy. The β -carbon atoms were removed and bonds saturated with hydrogen atoms (co-ordinates reported in Listing 1). This simplification helped us to keep the calculations simplewithout removing the most important features of the system. The Fe–N bonds were kept ata bond length of 1.989 Å. The molecule was kept planar with D h point group symmetry.Aromaticity was preserved, with ring current and complete electron delocalization in the“inner-cross”, an 18 π electrons system and 16 carbon atoms. This simplification does notintroduce bias towards the understanding of correlation effects in metallo-porphyrins.The molecule was placed on the xy plane, with the N atoms in between the x and y axes.The D h point group was used for all calculations, such that the π − π ∗ orbital system belongsto the b u , b g , b g and a u irreducible representations. Orbitals d z and d x − y belong to the a g representation and, the d xy , d xz and d yz to the b g , b g and b g respectively. Orbitals d xz and d yz and some of the π − π ∗ orbitals belong to the same irreducible representations ( b g and b g ) and, their overlap plays a major role in stabilizing the triplet spin state. Basis set and electron repulsion integrals.
Generally-contracted atomic natural or-bitals (ANO-RCC) basis sets have been employed, obtained from the Fe(21s15p10d6f4g2h),C,N(14s9p4d3f2g), H(8s4p3d1f) primitive functions. Two contraction schemes have beenemployed. In one case the primitive functions have been contracted to Fe(5S4P2D1F),C,N(3S2P1D), H(2S1P), giving a basis set of split-valence double- ζ plus polarization quality(VDZP). In another case the primitive functions have been contracted to Fe(6s5p3d2f1g),C,N(4s3p2d1f), H(3s2p1d), giving a basis set of split-valence triple- ζ plus polarization qual-ity (VTZP). This second basis set choice led to a total of 707 basis functions. Scalar rel-ativistic effects were introduced via second order Douglas-Kroll-Hess integral correction.The evaluation of the electron repulsion integrals has been greatly simplified by means ofthe resolution-of-identity Cholesky decomposition techniques as implemented in the Molcassoftware package, with a decomposition threshold of 10 − a.u. etails on the Complete Active Space choice. CASSCF is probably the simplest andmore natural method to tackle multi-configurational systems in chemistry.
The core con-cept of CASSCF is the active space, a list of “critical” orbitals with their electrons for whicha complete many-body expansion is generated and orbitals are variationally optimized underthe field generated by the multiconfigurational wave function, removing any bias related tothe choice of the trial orbitals. Three main weaknesses of the CASSCF method need tobe highlighted. (a) The larger the active space the exponentially larger the ConfigurationInteraction (CI) expansion, limiting the active space to at most 18 active electrons and 18active orbitals. (b) The active space represents a non-numerical parameter that introducesa certain level of arbitrariness and, active spaces that are smaller than the necessary mightreturn wrong energetics even upon perturbative correction. (c) Correlation outside the activespace is completely neglected at CASSCF level and post-CASSCF methods must be utilizedfor quantitative accuracy.Many methods have been proposed to reduce the exponential scaling of CAS wave func-tions. It is important to mention special forms of truncated CI expansions that can be ob-tained via the Restricted Active Space (RAS) and the Generalized Active Space (GAS)approaches.
It is also important to mention the Density Matrix Renormalization Group(DMRG) approach, the variational two-electron reduced density matrices approach and, the most recent Stochastic-CASSCF method.
The latter is the method of choicefor this report. These methods partially circumvent the exponential scaling problem andenable the investigation of larger active spaces. Using massively parallelized architecturesconventional CAS(20,20) calculations have recently been made possible. When a suffi-ciently large active space is employed the bias due to the choice of the active space is toa great extent lifted. In order to recover dynamic correlation outside the active space,perturbation theory to the second order (such as CASPT2,
NEVPT2 ) andmulti-reference configuration interaction (MRCI) using CASSCF wave functions as reference,have been employed with great success in a wide range of chemical systems. RASPT2 variants are also available. To date these methods represent the practicalstandard for transition metal chemistry. However, they become prohibitively expensive whencoupled to reference wave functions built from large active spaces, requiring in many casesfurther approximations, as discussed in great details in the literature. Additionally,a second order approximation will not account for higher-order correlation processes andorbital relaxation ( vide intra ). The Multi-Configuration Pair-Density Functional Theory(MCPDFT) method has been proposed as a cheap alternative to CASPT2.
Its com-putational cost is nearly independent of the size of the underlying reference wave functionand therefore nearly insensitive to the size of the active space.For the CASSCF calculations several active spaces have been chosen: (a) The CAS(6,5)is the smallest active space that includes solely the six valence electrons of the metal centreand its five 3d orbitals. (b) In the CAS(8,6) a doubly occupied σ orbital is added, mostlylocalized on the N atoms and pointing to the direction of the 3 d xy orbital of the iron atom.(c) The CAS(8,12) adds five empty correlating d and the Fe 4s orbitals into the activespace. The Fe 4s has been added as it could compete with the d orbitals in accounting forcorrelation effects. (d) In the CAS(14,16), the role of the frontier π orbitals was probed. Thetwo highest bonding π orbitals and the two lowest anti-bonding π ∗ orbitals have been addedto Pierloot CAS(8,11) active space together with their four electrons. We also added thedoubly occupied 3s orbital on the Fe centre to investigate its role in the electron correlationlandscape. (e) The much larger active space, CAS(32,34), consists of the 10 Fe (3 d , d ) orbitalsand their 6 electrons, the entire π system (18 electrons and 16 orbitals), the four orbitals ofthe Fe (4s4p) shell and four doubly occupied N (2 p x ,2 p y ) symmetrically combined “radial”orbitals pointing at the metal centre. The four remaining N (2 p x ,2 p y ) orbitals, symmetricallycombined to form “tangential” orbitals, were not included in the active space. The orbitalscorrelated in the CAS(32,34) are different than the ones used in our previous work. Forthe CAS(32,29), only valence orbitals on the macrocycle and the metal centre were included.In the present work double-shell orbitals, “radial” N 2p orbitals and, the Fe (4 s ,4 p ) orbitals23ave been explicitely correlated. The enlarged active space shows the breathing mechanism and provides correct energy ordering of the spin states.For the small active spaces (cases (a) to (d) above) the CASPT2 method has been usedto recover dynamic correlation outside the active space. At CASPT2 level, core orbitals (1 s on C and N atoms and 1 s , 2 s , 2 p orbitals on Fe atom) were kept frozen. The standardIPEA zeroth order Hamiltonian has been utilized with the default IPEA denominator shiftof 0.25 a.u. No method for dynamic correlation has been coupled to the Stochastic-CASSCFwave functions. Calculations at CASPT2 level show that predictions are largely affected bythe choice of the under-lying active space with the high spin state still over-stabilized overthe triplet spin state. For comparison with the Stochastic-CAS(32,34) results, also RASSCFcalculations have been performed, in which 32 electrons and 34 active orbitals have beenpartially correlated. Following Pierloot’s approach, in our RAS(32,34) twelve orbitals wereput in the RAS1 space, including three σ (N 2p) orbitals and, nine π orbitals, six orbitalsin RAS2, including one bonding σ orbital and the five 3d orbitals and, sixteen orbitals inRAS3, including five double-shell correlating orbitals and the π ∗ orbitals of the macrocycle.Only single and double excitations out of RAS1 and into RAS3 have been allowed. Details on the Quantum Monte Carlo setup.
For the FCIQMC dynamics the initiator formulation of the method has been used with a threshold value of n a = 3 . using a deterministic subspace consisting of |D| =10000 most populated determinants. In the initiator approximation, walkers populatingdeterminants with largest weight are promoted to be “initiators”. Only initiators are ableto spawn walkers on empty determinants. Non-initiators are allowed to spawn only ondeterminants that are already occupied. The calculations were run in replica mode inorder to sample the one- and two-body reduced density matrices necessary to the orbitalrotation step. CASSCF natural orbitals from smaller active spaces were used as the startingorbitals for CASSCF optimizations with larger active spaces. 5 × walkers were employed24or the initial dynamics and the first five CASSCF iterations. The number of walkers wasgradually increased to 5 × . CASSCF convergence was reached at this walker populationfor all states here investigated. The approach of increasing the walker population in stepsfollows from our initial findings, already discussed in our previous paper. A small walkerdistribution is able to generate a convenient averaged field to allow for an effective orbitaloptimization step at the early stages of the CASSCF procedure. This procedure guaranteesfast orbital rotations and a limited number of CASSCF iterations at the high-populationregime when sub-milliHartree accuracy is required. After CASSCF convergence was reached,more refined solutions were obtained by increasing the target number of walkers to 1 × .This procedure is standard to reduce the initiator error on the stochastic sampling of thewave function and, was used to confirm that no bias on the spin splitting was introduced dueto undersampling of the determinantal space. The time-step ∆ τ was found via an automaticsearch procedure for each simulation, and took typical values in the range 5 − × − a.u.A typical FCIQMC simulation, took ∼
24 hours for each CASSCF iteration on 640 cores.Orbital rotations were performed using the Super-CI method with a quasi-Newton update.The entire CASSCF procedure converged in 10-15 iterations for the states here investigated.All calculations have been performed using the OpenMolcas chemistry software package. Past and present proposed active spaces.
The smallest active space for this systemwould be a CAS(6,5) including only the six valence electrons in the five valence 3d orbitals.Pierloot pointed out that one additional σ Fe–N bonding orbital, and its two electrons, mustbe included in the active space, leading to a CAS(8,6). She also found that it is crucial toinclude five double-shell correlating orbitals in the active space when a PT2 treatment is usedfor dynamic correlation. These orbitals account for a quite strong radial electron correlationand, when included into the active space, lead to a CAS(8,11). Recently Pierloot and co-workers have investigated the role of the (3s3p) and (4s4p) shells to give a CAS(16,15) anda CAS(16,19), respectively. For the latter being prohibitively large, the RASSCF/RASPT225r the DMRG/PT2 methodologies have been used. Pierloot has also considered the roleof the π -system of the macrocycle, by means of the RASSCF/RASPT2 approach. ARAS(34,2,2;13,6,16) has been chosen containing a total of 34 electrons and 35 orbitals. RAS1contains 13 active orbitals, doubly occupied in the reference determinant, RAS2 contains only6 active orbitals and RAS3 the remaining 16 that are empty in the reference determinant.Only single and double excitations were allowed from RAS1 and to RAS3.Among the successful wave function methods able to predict a triplet ground state, it isimportant to mention Radoń’s CCSD(T) calculations, that have shown a triplet ground stateonly upon extrapolation to the complete basis set limit, a DMRG-CI study and, therecent Heat-bath Configuration Interaction Self-Consistent Field (HCISCF) by Sharma. In the last two cases a CAS(44,44) was chosen.The CAS(32,34) used in the present report differs from the one used in the past.
The CAS(44,44) reported in the literature includes 4 p x and 4 p y orbitals of the metal centrebut, does not include the 4 s or the 4 p z orbitals. Therefore, the (4s,4p) shell is somehowincomplete. Olivares-Amaya and Sharma included the entire list of eight MOs resultingfrom a symmetry adapted combination of the 2 p x and 2 p y orbitals on the N atoms. In thepresent work only four have been included, leaving the tangential ones in the inactive space.No orbitals have been frozen or deleted at the CASSCF level of theory.In the paper introducing the Stochastic-CASSCF method, we performed calculationson a model system of Fe(II)-porphyrin. Our model system had a Fe–N bond length of2.05 Å, which is closer to the one predicted for the quintet A g state. An active space of 32electrons in 29 orbitals was chosen, including the 24 orbitals of the aromatic π -system on themacrocycle with their 26 electrons and the five 3d valence orbitals of the metal centre withtheir six electrons. The Stochastic-CASSCF(32,29) led to a quintet ground state ( A g ) withthe triplet at 14 kcal/mol. Neither the double-shell orbitals, nor the bonding Fe–N σ orbitalwere added to the active space. Inspection of the CAS(32,29) wave function revealed thatthe dominant configuration (3 d z ) (3 d x − y ) (3 d yz ) (3 d xz ) (3 d xy ) , contributed for 37.3%.26he other relevant configurations (with weight > π orbitals. Older investigationsconcluded that the iron 3d orbitals and the porphyrin π system were well separated, whereasour findings showed the opposite, with the Stochastic-CASSCF wave function showing a closeinteraction between the aromatic ring and the metal centre via non-negligible charge-transferconfigurations.The distribution of orbitals for the CAS(32,34) proposed here in the eight irreduciblerepresentations of the D h point group is given in Table S1.Table S1: Distribution of molecular orbitals among inactive, active and secondary spaces foreach irreducible representation.Ag B3u B2u B1g B1u B2g B3g AuInactive 16 13 13 9 2 0 0 0Active 6 2 2 3 6 6 6 3Secondary VTZP 113 104 104 93 59 51 51 45 Energy Splittings with VDZP basis set.
Energy splitting estimates at the various levelof theory and VDZP basis set are summarized in Figure S1. A comparison with Figure 2of the main text (where a larger basis set of TZVP quality is used) shows the effect of thebasis set on the spin ordering. 27 background.dat’30.00.05.010.015.020.025.0 B
2g 5 E g 3 B
1g 3 A
2g 3 E g Spin States59.063.068.0 R e l a ti v e E n e r gy [ k ca l/ m o l ] CAS(6,5)CAS(8,6)CAS(8,12)CAS(14,16) PT2(6,5)PT2(8,6)PT2(8,12)PT2(14,16)
Figure S1: Energy splittings relative to the A g state in VDZP basis set. The red dashedline marks the lowest triplet-quintet energy gap obtained.28isting 1: Cartesian coordinates for the Fe(II)-porphyrin model system (Angstrom) A g state of the Fe(II)-porphyrin model system and their occupation numbers.30 n the mixing of metal-d and π orbitals. In order to understand the importance ofour results on the mixing of the metal and ligand natural orbitals for the triplet spin stateand not for the quintet spin state, we would like to comment on some properties of theCASSCF approach. CASSCF is invariant to orbital rotations within the active space andthus rotations within the active orbitals do not alter the energy gap between the states.Therefore any sort of unitary transformation acting on the active orbitals would not changethe physics of the states. Here we use natural orbitals. Natural orbitals follow a strict recipe– they diagonalize the reduced one-body density matrix – that leads to a uniquely definedset of rotated CASSCF active orbitals. When occupation number degeneracies arise anyrotations within orbitals with same occupancy is possible and in this case natural orbitalscan be presented in various ways. This special case does not arise in the current investigation.In computing the natural orbitals one may think of a two step procedure. In a first step,active orbitals are localized in a way that each of them is either sitting on the metal centre oron the macrocycle. This step is conceptually important to be able to clearly state in whichpart of the molecule one electron (or pair of electrons) is located; however, in practice, it isnot needed as the shape of the natural orbitals will not be affected by the starting orbitalsused to build the one-body density matrix. In a second step, starting from these localizedorbitals, the one-body density matrix is built and diagonalized. Diagonalization of the one-body density matrix is the responsible step that leads (or not) to the mixing of the localizedorbitals. If non negligible off-diagonal elements of the one-body density matrix exist, theywill have an effect on the mixing. Non-negligible off-diagonal elements can exist only forcorrelated systems. This is exactly what we observe in the system under investigation. π orbitals and out-of-plane 3d orbitals are correlated via charge-transfer excitations, they willlead to non-negligible off-diagonal elements in the one-body density matrix and thus to themixing of the metal centre and the ligand orbitals.31 otes on the correlating d’ orbitals. In our active space we added a set of five corre-lating d’ orbitals. These are known as double-shell d orbitals. They have a nodal structurethat resemble the 4d physical orbitals. Nonetheless, their radial distribution is closer to the3d orbitals and their energy is higher than the 4d orbitals. They could be defined as con-tracted 4d orbitals. This feature is very general and not specific to the present compound.CASSCF forces them to be closer to the 3d orbitals to maximize their overlap and the radialcorrelation of the 3d electrons. They are responsible for the breathing effect , according towhich valence electrons of the correlated method are more expanded than the equivalentelectron in a non correlated (or less correlated) approach. This breathing mechanism leadsto the reduction of the on-site electron repulsion of the doubly occupied 3d orbitals of themetal centre. Wave function analysis.
The E g state is dominated by the ( σ N ) ( π ) ( σ ML ) (3 d x − y ) (3 d xz , d yz ) (3 d z ) configuration with a weight of 59%. The second two most relevant con-figurations (contributing for 3% and 2% respectively) consist of double excitations fromthe π NB to the π ∗ Lx and the π ∗ Ly respectively. These two configurations are the main re-sponsible for populating the π ∗ Lx and the π ∗ Ly natural orbitals. The same configurationscontribute for less than 0.5% in the equivalent RAS wave function. π → π ∗ excitationscontribute for ∼
1% to the wave function. σ ML → d and 3 d → d excitations also appearin the leading determinant of the triplet spin state. The A g state is dominated by the( σ N ) ( π ) ( σ ML ) (3 d z ) (3 d x − y ) (3 d xz ) (3 d yz ) (3 d xy ) configuration with a weight of 61%.Other relevant configurations are the double excitations from the π NB and π NB to the π ∗ Lx and the π ∗ Ly respectively contributing for 3% each to the wave function. The 3 d z → d x − y and 3 d z → s excitations contribute for ∼ π → π ∗ (ring delocalization) and 3 d → d excitations (radialcorrelation). 32isting 2: A possible FCIQMC-CASSCF input setup for the porphyrin molecule withinOpenMolcas & GATEWAYRICDCoord3 A2g . xyzbasisANO - RCC - VTZPgroupfull& SEWARD& RASSCFNECIEXNENWAL50000000DEFD1 2 3 4 5 13 14 17 18 21 22 27 28 29 30 3132 39 40 41 42 43 51 52 53 54 55 56 63 64 65 66CYCLe200000RSPCutoff0.3TIMENeci2000RDMStart50000RDMPick1000THRS1.0 e -4 1.0 e -1 5.0 e -4spin3Symmetry6nactel32 0 0inactive16 13 13 9 2 0 0 0ras26 2 2 3 6 6 6 3deleted0 0 0 0 0 0 0 0 Listing 3: Relevant input keywords used by the FCIQMC program NECI
TitleSystem readelectrons 32nonuniformrandexcits 4 ind - weighted -2nobrillouintheoremspin - restrict 2freeformatendsyscalcdefinedet 1 2 3 4 5 13 14 17 18 21 22 27 2829 30 31 32 39 40 41 42 43 51 52 53 54 55 5663 64 65 66methodsmethod vertex fcimcendmethodstotalwalkers 1000000000diagshift 0.00 hiftdamp 0.02nmcyc 200000stepsshift 10proje - changeref 1.2truncinitiatoraddtoinitiator 3allrealcoeffrealspawncutoff 0.30jump - shifttau 0.001 searchmax - tau 0.02maxwalkerbloom 1memoryfacspawn 10.0memoryfacpart 5.0time 1400startsinglepart 100readpopswalkcontgrowsemi - stochasticpops - core 10000trial - wavefunctionpops - trial 500rdmsamplingiters 10000endcalcloggingPRINT - SPIN - RESOLVED - RDMS( READRDMSHDF5 - POPSHighlypopwrite 200( binarypopsprintonerdmdiagflyonerdmcalcrdmonfly 3 5000 500endlogend References (1) Sono, M.; Roach, M. P.; Coulter, E. D.; Dawson, J. H.
Chem. Rev. , , 2841–2888.(2) Meunier, B.; de Visser, S. P.; Shaik, S. Chem. Rev. , , 3947–3980.(3) Schöneboom, J. C.; Lin, H.; Reuter, N.; Thiel, W.; Cohen, S.; Ogliaro, F.; Shaik, S. J. Am. Chem. Soc. , , 8142–8151.(4) Almlöf, J. Int. J. Quantum Chem. , , 915–924.(5) Jones, D. H.; Hinman, A. S.; Ziegler, T. Inorg. Chem. , , 2092–2095.(6) Rosa, A.; Baerends, E. J. Inorg. Chem. , , 5637–5639.347) Collman, J. P.; Hoard, J. L.; Kim, N.; Lang, G.; Reed, C. A. J. Am. Chem. Soc. , , 2676–2681.(8) Lang, G.; Spartalian, K.; Reed, C. A.; Collman, J. P. J. Chem. Phys. , ,5424–5427.(9) Boyd, P. D. W.; Buckingham, A. D.; McMeeking, R. F.; Mitra, S. Inorg. Chem. , , 3585–3591.(10) Goff, H.; La Mar, G. N.; Reed, C. A. J. Am. Chem. Soc. , , 3641–3647.(11) Mispelter, J.; Momenteau, M.; Lhoste, J. M. J. Chem. Phys. , , 1003–1012.(12) Kitagawa, T.; Teraoka, J. Chem. Phys. Lett. , , 443–445.(13) Sams, J. R.; Tsin, T. B. Chem. Phys. Lett. , , 599–601.(14) Matsuzawa, N.; Ata, M.; Dixon, D. A. J. Phys. Chem. , , 7698–7706.(15) Kozlowski, P. M.; Spiro, T. G.; Bérces, A.; Zgierski, M. Z. J. Phys. Chem. B , , 2603–2608.(16) Liao, M.-S.; Scheiner, S. J. Chem. Phys. , , 3635–3645.(17) Smith, D. M. A.; Dupuis, M.; Straatsma, T. P. Mol. Phys. , , 273–278.(18) Chen, H.; Lai, W.; Shaik, S. J. Phys. Chem. B , , 1727–1742.(19) Choe, Y.-K.; Nakajima, T.; Hirao, K.; Lindh, R. J. Chem. Phys. , , 3837–3844.(20) Choe, Y.-K.; Hashimoto, T.; Nakano, H.; Hirao, K. Chem. Phys. Lett. , ,380–388.(21) Vancoillie, S.; Zhao, H.; Tran, V. T.; Hendrickx, M. F. A.; Pierloot, K. J. Chem.Theory Comput. , , 3961–3977. 3522) Pierloot, K.; Phung, Q. M.; Domingo, A. J. Chem. Theory Comput. , , 537–553.(23) Phung, Q. M.; Wouters, S.; Pierloot, K. J. Chem. Theory Comput. , , 4352–4361.(24) Radoń, M.; Pierloot, K. J. Phys. Chem. A , , 11824–11832.(25) Vancoillie, S.; Zhao, H.; Radoń, M.; Pierloot, K. J. Chem. Theory Comput. , ,576–582.(26) Radoń, M.; Broclawik, E.; Pierloot, K. J. Phys. Chem. B , , 1518–1528.(27) Radoń, M. J. Chem. Theory Comput. , , 2306–2321.(28) Swart, M.; Groenhof, A. R.; Ehlers, A. W.; Lammertsma, K. J. Phys. Chem. A , , 5479–5483.(29) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. , , 3865–3868.(30) Handy, N. C.; Cohen, A. J. Mol. Phys. , , 403–412.(31) Swart, M. J. Chem. Theory Comput. , , 2057–2066.(32) Swart, M. Int. J. Quantum Chem. , , 2–7.(33) Costas, M.; Harvey, J. N. Nat. Chem. , , 7–9.(34) Song, S.; Kim, M.-C.; Sim, E.; Benali, A.; Heinonen, O.; Burke, K. arXiv:1708.08425 (35) Li Manni, G.; Smart, S. D.; Alavi, A. J. Chem. Theory Comput. , , 1245–1258.(36) Roos, B. O.; Taylor, P. R.; Siegbahn, P. E. M. Chem. Phys. , , 157–173.(37) Roos, B. O. Int. J. Quantum Chem. , , 175–189.(38) Siegbahn, P. E. M.; Heiberg, A.; Roos, B. O.; Levy, B. Phys. Scr. , , 323–327.3639) Siegbahn, P. E. M.; Almlöf, J.; Heiberg, A.; Roos, B. O. J. Chem. Phys. , ,2384–2396.(40) Roos, B. O. Advances in Chemical Physics ; 2007; pp 399–445.(41) Andersson, K.; Malmqvist, P.-Å.; Roos, B. O.; Sadlej, A. J.; Wolinski, K.
J. Phys.Chem. , , 5483–5488.(42) Andersson, K.; Malmqvist, P.-Å.; Roos, B. O. J. Chem. Phys. , , 1218–1226.(43) Andersson, K.; Roos, B. O. Int. J. Quantum Chem. , , 591–607.(44) Finley, J.; Malmqvist, P.-Å.; Roos, B. O.; Serrano-Andrés, L. Chem. Phys. Lett. , , 299–306.(45) Ghigo, G.; Roos, B. O.; Malmqvist, P.-Å. Chem. Phys. Lett. , , 142–149.(46) Pierloot, K. Mol. Phys. , , 2083–2094.(47) Pierloot, K.; Vancoillie, S. J. Chem. Phys. , , 034104.(48) Pierloot, K.; Vancoillie, S. J. Chem. Phys. , , 124303.(49) Widmark, P.-O.; Malmqvist, P.-Å.; Roos, B. O. Theor. Chem. Acc. , , 291–306.(50) Roos, B. O.; Lindh, R.; Malmqvist, P.-Å.; Veryazov, V.; Widmark, P.-O. J. Phys.Chem. A , , 2851–2858.(51) Aquilante, F. et al. J. Comput. Chem. , , 506–541.(52) Aquilante, F.; Lindh, R.; Pedersen, T. B. J. Chem. Phys. , , 114107.(53) Aquilante, F.; Pedersen, T. B.; Lindh, R. J. Chem. Phys. , , 194106.(54) Aquilante, F.; Pedersen, T. B.; Lindh, R.; Roos, B. O.; de Merás, A. S.; Koch, H. J.Chem. Phys. , , 024113. 3755) Gagliardi, F. A. L.; Pedersen, T. B.; Lindh, R. J. Chem. Phys. , , 154107.(56) Pedersen, T. B.; Aquilante, F.; Lindh, R. Theor. Chem. Acc. , , 1–10.(57) Olsen, J.; Roos, B. O.; Jørgensen, P.; Jensen, H. J. A. J. Chem. Phys. , ,2185–2192.(58) Malmqvist, P.-Å.; Rendell, A.; Roos, B. O. J. Phys. Chem. , , 5477–5482.(59) Celani, P.; Werner, H. J. J. Chem. Phys. , , 5546–5557.(60) Ma, D.; Li Manni, G.; Gagliardi, L. J. Chem. Phys. , , 044128.(61) Vogiatzis, K. D.; Li Manni, G.; Stoneburner, S. J.; Ma, D.; Gagliardi, L. J. Chem.Theory Comput. , , 3010–3021.(62) Odoh, S. O.; Li Manni, G.; Carlson, R. K.; Truhlar, D. G.; Gagliardi, L. Chem. Sci. , , 2399–2413.(63) White, S. R. Phys. Rev. Lett. , , 2863–2866.(64) Chan, G. K.-L.; Head-Gordon, M. J. Chem. Phys. , , 4462–4476.(65) White, S. R. Phys. Rev. B , , 180403(R).(66) Schollwöck, U. Rev. Mod. Phys. , , 259–315.(67) Marti, K. H.; Ondík, I. M.; Moritz, G.; Reiher, M. J. Chem. Phys. , , 014104.(68) Chan, G. K.-L.; Kállay, M.; Gauss, J. J. Chem. Phys. , , 6110–6116.(69) Chan, G. K.-L.; Sharma, S. Annu. Rev. Phys. Chem. , , 465–481.(70) Yanai, T.; Kurashige, Y.; Ghosh, D.; Chan, G. K.-L. Int. J. Quantum Chem. , , 2178–2190. 3871) Fosso-Tande, J.; Nascimento, D. R.; DePrince III, A. E. Mol. Phys. , , 423–430.(72) Fosso-Tande, J.; Nguyen, T.-S.; Gidofalvi, G.; DePrince, A. E. I. J. Chem. TheoryComput. , , 2260–2271.(73) Mazziotti, D. A. Chem. Rev. , , 244–262.(74) Pelzer, K.; Greenman, L.; Gidofalvi, G.; Mazziotti, D. A. J. Phys. Chem. A , , 5632–5640.(75) Gidofalvi, G.; Mazziotti, D. A. J. Chem. Phys. , , 134108.(76) Mazziotti, D. A. Reduced-Density-Matrix Mechanics: With Application to Many-Electron Atoms and Molecules ; 2007; pp 19–59.(77) Mazziotti, D. A.
Acc. Chem. Res. , , 207–215.(78) Thomas, R. E.; Sun, Q.; Alavi, A.; Booth, G. H. J. Chem. Theory Comput. , ,5316–5325.(79) Vogiatzis, K. D.; Ma, D.; Olsen, J.; Gagliardi, L.; de Jong, W. arXiv:1707.04346 (80) Pulay, P. Int. J. Quantum Chem. , , 3273–3279.(81) Forsberg, N.; Malmqvist, P.-Å. Chem. Phys. Lett. , , 196–204.(82) Angeli, C.; Bories, B.; Cavallini, A.; Cimiraglia, R. J. Chem. Phys. , , 054108.(83) Angeli, C.; Borini, S.; Cavallini, A.; Cestari, M.; Cimiraglia, R.; Ferrighi, L.; Sparta, M. Int. J. Quantum Chem. , , 686–691.(84) Angeli, C.; Borini, S.; Cestari, M.; Cimiraglia, R. J. Chem. Phys. , , 4043–4049. 3985) Angeli, C.; Cimiraglia, R.; Evangelisti, S.; Leininger, T.; Malrieu, J.-P. J. Chem. Phys. , , 10252–10264.(86) Angeli, C.; Cimiraglia, R.; Malrieu, J.-P. J. Chem. Phys. , , 9138–9153.(87) Angeli, C.; Cimiraglia, R.; Malrieu, J.-P. Chem. Phys. Lett. , , 297–305.(88) Angeli, C.; Pastore, M.; Cimiraglia, R. Theor. Chem. Acc. , , 743–754.(89) Sauri, V.; Serrano-Andrés, L.; Shahi, A. R. M.; Gagliardi, L.; Vancoillie, S.; Pier-loot, K. J. Chem. Theory Comput. , , 153–168.(90) Malmqvist, P.-Å.; Pierloot, K.; Shahi, A. R. M.; Cramer, C. J.; Gagliardi, L. J. Chem.Phys. , , 204109.(91) Ma, D.; Li Manni, G.; Olsen, J.; Gagliardi, L. J. Chem. Theory Comput. , ,3208–3213.(92) Kurashige, Y.; Yanai, T. J. Chem. Phys. , , 094104.(93) Kurashige, Y. Mol. Phys. , , 1485–1494.(94) Yanai, T.; Kurashige, Y.; Mizukami, W.; Chalupská, J.; Lan, T. N.; Saitow, M. Int.J. Quantum Chem. , , 283–299.(95) Li Manni, G.; Carlson, R. K.; Luo, S.; Ma, D.; Olsen, J.; Truhlar, D. G.; Gagliardi, L. J. Chem. Theory Comput. , , 3669–3680.(96) Carlson, R. K.; Li Manni, G.; Sonnenberger, A. L.; Truhlar, D. G.; Gagliardi, L. J.Chem. Theory Comput. , , 82–90.(97) Gagliardi, L.; Truhlar, D. G.; Li Manni, G.; Carlson, R. K.; Hoyer, C. E.; Bao, J. L. Acc. Chem. Res. , , 66–73. 4098) Hoyer, C. E.; Ghosh, S.; Truhlar, D. G.; Gagliardi, L. J. Phys. Chem. Lett. , ,586–591.(99) Wilbraham, L.; Verma, P.; Truhlar, D. G.; Gagliardi, L.; Ciofini, I. J. Phys. Chem.Lett. , , 2026–2030.(100) Ghosh, S.; Sonnenberger, A. L.; Hoyer, C. E.; Truhlar, D. G.; Gagliardi, L. J. Chem.Theory Comput. , , 3643–3649.(101) Hoyer, C. E.; Gagliardi, L.; Truhlar, D. G. J. Phys. Chem. Lett. , , 4184–4188.(102) Ghosh, S.; Cramer, C. J.; Truhlar, D. G.; Gagliardi, L. Chem. Sci. , , 2741–2750.(103) Cleland, D.; Booth, G. H.; Alavi, A. J. Chem. Phys. , , 041103.(104) Cleland, D.; Booth, G. H.; Alavi, A. J. Chem. Phys. , , 024112.(105) Petruzielo, F. R.; Holmes, A. A.; Changlani, H. J.; Nightingale, M. P.; Umrigar, C. J. Phys. Rev. Lett. , , 230201.(106) Blunt, N. S.; Smart, S. D.; Kersten, J. A.-F.; Spencer, J. S.; Booth, G. H.; Alavi, A. J. Chem. Phys. , , 184107.(107) Overy, C.; Booth, G. H.; Blunt, N. S.; Shepherd, J. J.; Cleland, D.; Alavi, A. J. Chem.Phys. , , 244117.(108) Booth, G. H.; Smart, S. D.; Alavi, A. Mol. Phys. , , 1855–1869.(109) Olivares-Amaya, R.; Hu, W.; Nakatani, N.; Sharma, S.; Yang, J.; Chan, G. K.-L. J.Chem. Phys. , , 034102.(110) Smith, J. E.; Mussard, B.; Holmes, A. A.; Sharma, S. J. Chem. Theory Comput.2017