Cumulant Approximated Second-Order Perturbation Theory Based on the Density Matrix Renormalization Group for Transition Metal Complexes: A Benchmark Study

Perspective: multi-configurational methods in bio-inorganic chemistry

Transition metal ions play crucial roles in the structure and function of numerous proteins, contributing to essential biological processes such as catalysis, electron transfer, and oxygen binding. However, accurately modeling the electronic structure and properties of metalloproteins poses significant challenges due to the complex nature of their electronic configurations and strong correlation effects. Multiconfigurational quantum chemistry methods are, in principle, the most appropriate tools for addressing these challenges, offering the capability to capture the inherent multi-reference character and strong electron correlation present in bio-inorganic systems. Yet their computational cost has long hindered wider adoption, making methods such as density functional theory (DFT) the method of choice. However, advancements over the past decade have substantially alleviated this limitation, rendering multiconfigurational quantum chemistry methods more accessible and applicable to a wider range of bio-inorganic systems. In this perspective, we discuss some of these developments and how they have already been used to answer some of the most important questions in bio-inorganic chemistry. We also comment on ongoing developments in the field and how the future of the field may evolve.

Ab Initio Calculation of UV-vis Absorption of Parent Mg, Fe, Co, Ni, Cu, and Zn Metalloporphyrins

Relativistic restricted active space (RAS) second-order multireference perturbation theory (MRPT2) methods, incorporating spin-orbit (SO) coupling perturbatively via state interaction (SO-MRPT2/RASSCF), were used to reproduce the absorption spectra of parent metalloporphyrins containing the Mg2+, Zn2+, Co2+, Ni2+, Cu2+, or FeCl2+ ions in the 12,500-40,000 cm-1 region. Particular attention was paid to the interaction between the porphyrin ring and the metal 3d electrons in states of different multiplicities (we used metal 3d and double d-shell or 3d' orbitals). For this class of compounds, the N-electron valence state perturbation theory (NEVPT2) method is superior to the complete active space perturbation theory (CASPT2) and successfully reproduces the energies of all four characteristic transitions (Q, B, N, and L) of closed-shell metalloporphyrins. Inclusion of SO coupling was found to have very little effect on excitation energies and oscillator strengths. For FeCl2+ porphyrin, we treated ligand-to-metal charge-transfer (LMCT; π,d), metal ligand field (d,d), and metal-to-ligand charge-transfer (MLCT; d,π*) transitions within the same framework. The broad and intense spectral features associated with its B (Soret) band are attributed to multiconfigurational LMCT (d,π*) bands involving strong metal-ligand orbital mixing.

Can Iron-Porphyrins Behave as Single-Molecule Magnets?

The study of magnetic properties, especially the magnetic anisotropy of iron-porphyrin complexes employing multiconfigurational methods, is quite challenging due to many strongly correlated electrons in nearly degenerate orbitals. However, a prerequisite for observing the magnetic anisotropy and slow magnetization relaxation, the zero-field splitting parameter, D, was experimentally observed decades ago for halide-based axially ligated penta-coordinate Fe(III)-porphyrins. In these complexes, the signs of D were reported mostly as positive; in a few cases, inconclusive signs of the D parameter were also mentioned. However, no ab initio calculations have been reported to shed light on this. Deciphering the electronic structure of these penta-coordinated complexes employing the complete active space self-consistent field method and N-electron valence second-order perturbation theory, we confirm the positive D values. However, a negative D value is highly desired to observe the single-molecule magnet properties without an external magnetic field, which we observed in the Fe(II)-porphyrin complexes with axial imidazole ligands instead of halide ligands. The detailed analysis of the multireference wave functions unravels the role of axial ligands in determining the sign and magnitude of the D parameters.

Renormalized-Residue-Based Multireference Configuration Interaction Method for Strongly Correlated Systems

The implementation of multireference configuration interaction (MRCI) methods in quantum systems with large active spaces is hindered by the expansion of configuration bases or the intricate handling of reduced density matrices (RDMs). In this work, we present a spin-adapted renormalized-residue-based MRCI (RR-MRCI) approach that leverages renormalized residues to effectively capture the entanglement between active and inactive orbitals. This approach is reinforced by a novel efficient algorithm, which also facilitates an efficient deployment of spin-adapted matrix product state MRCI (MPS-MRCI). The RR-MRCI framework possesses several advantages: (1) It considers the orbital entanglement and utilizes highly compressed MPS structure, improving computational accuracy and efficiency compared with internally contracted (ic) MRCI. (2) Utilizing small-sized buffer environments of a few external orbitals as probes based on quantum information theory, it enhances computational efficiency over MPS-MRCI and offers potential application to large molecular systems. (3) The RR framework can be implemented in conjunction with ic-MRCI, eliminating the need for high-rank RDMs, by using distinct renormalized residues. We evaluated this method across nine diverse molecular systems, including Cu2O22+ with an active space of (24e,24o) and two complexes of lanthanide and actinide with active space (38e,36o), demonstrating the method's versatility and efficacy.

Block-Correlated Coupled Cluster Theory with up to Four-Pair Correlation for Accurate Static Correlation of Strongly Correlated Systems

A block-correlated coupled cluster method with up to four-pair correlation based on the generalized valence bond wave function (GVB-BCCC4) is first implemented, which offers an alternative method for electronic structure calculations of strongly correlated systems. We developed some techniques to derive a set of compact and cost-effective equations for GVB-BCCC4, which include the definition of n-block (n = 1-4) Hamiltonian matrices, the combination of excitation operators, and the definition of independent amplitudes. We then applied the GVB-BCCC4 method to investigate several potential energy surfaces of strongly correlated systems with singlet ground states. Our calculations demonstrate that the GVB-BCCC4 method can provide nearly exact static correlation energies as the density matrix renormalization group method (on the basis of the same GVB orbitals). This work highlights the significance of four-pair correlation in quantitative descriptions of static correlation energy for strongly correlated systems.

Toward a Stochastic Complete Active Space Second-Order Perturbation Theory

In this work, an internally contracted stochastic complete active space second-order perturbation theory, stochastic-CASPT2, is reported. The method relies on stochastically sampled reduced density matrices (RDMs) up to rank four and contractions thereof with the generalized Fock matrix. A new protocol for calculating higher-order RDMs in full configuration interaction quantum Monte Carlo (FCIQMC) has been designed based on (1) restricting sampling of the corresponding excitations to a deterministic subspace, (2) averaging the RDMs from independent dynamics and (3) projecting them onto the closest positive semi-definite matrix. Our protocol avoids previously encountered numerical conditioning problems in the orthogonalization of the perturber overlap matrix stemming from numerical noise. The chromium dimer CASSCF(12,12)/CASPT2 binding curve is computed as a proof of concept.

Block2: A comprehensive open source framework to develop and apply state-of-the-art DMRG algorithms in electronic structure and beyond

block2 is an open source framework to implement and perform density matrix renormalization group and matrix product state algorithms. Out-of-the-box it supports the eigenstate, time-dependent, response, and finite-temperature algorithms. In addition, it carries special optimizations for ab initio electronic structure Hamiltonians and implements many quantum chemistry extensions to the density matrix renormalization group, such as dynamical correlation theories. The code is designed with an emphasis on flexibility, extensibility, and efficiency and to support integration with external numerical packages. Here, we explain the design principles and currently supported features and present numerical examples in a range of applications.

Local Oxidation States in {FeNO}6-8 Porphyrins: Insights from DMRG/CASSCF-CASPT2 Calculations

A first DMRG/CASSCF-CASPT2 study of a series of paradigmatic {FeNO}6, {FeNO}7, and {FeNO}8 heme-nitrosyl complexes has led to substantial new insight as well as uncovered key shortcomings of the DFT approach. By virtue of its balanced treatment of static and dynamic correlation, the calculations have provided some of the most authoritative information available to date on the energetics of low- versus high-spin states of different classes of heme-nitrosyl complexes. Thus, the calculations indicate low doublet-quartet gaps of 1-4 kcal/mol for {FeNO}7 complexes and high singlet-triplet gaps of ≳20 kcal/mol for both {FeNO}6 and {FeNO}8 complexes. In contrast, DFT calculations yield widely divergent spin state gaps as a function of the exchange-correlation functional. DMRG-CASSCF calculations also help calibrate DFT spin densities for {FeNO}7 complexes, pointing to those obtained from classic pure functionals as the most accurate. The general picture appears to be that nearly all the spin density of Fe[P](NO) is localized on the Fe, while the axial ligand imidazole (ImH) in Fe[P](NO)(ImH) pushes a part of the spin density onto the NO moiety. An analysis of the DMRG-CASSCF wave function in terms of localized orbitals and of the resulting configuration state functions in terms of resonance forms with varying NO(π*) occupancies has allowed us to address the longstanding question of local oxidation states in heme-nitrosyl complexes. The analysis indicates NO(neutral) resonance forms [i.e., Fe(II)-NO0 and Fe(III)-NO0] as the major contributors to both {FeNO}6 and {FeNO}7 complexes. This finding is at variance with the common formulation of {FeNO}6 hemes as Fe(II)-NO+ species but is consonant with an Fe L-edge XAS analysis by Solomon and co-workers. For the {FeNO}8 complex {Fe[P](NO)}-, our analysis suggests a resonance hybrid description: Fe(I)-NO0 ↔ Fe(II)-NO-, in agreement with earlier DFT studies. Vibrational analyses of the compounds studied indicate an imperfect but fair correlation between the NO stretching frequency and NO(π*) occupancy, highlighting the usefulness of vibrational data as a preliminary indicator of the NO oxidation state.

Intruder-free cumulant-truncated driven similarity renormalization group second-order multireference perturbation theory

Accurate multireference electronic structure calculations are important for constructing potential energy surfaces. Still, even in the case of low-scaling methods, their routine use is limited by the steep growth of the computational and storage costs as the active space grows. This is primarily due to the occurrence of three- and higher-body density matrices or, equivalently, their cumulants. This work examines the effect of various cumulant truncation schemes on the accuracy of the driven similarity renormalization group second-order multireference perturbation theory. We test four different levels of three-body reduced density cumulant truncations that set different classes of cumulant elements to zero. Our test cases include the singlet-triplet gap of CH2, the potential energy curves of the XΣg+1 and AΣu+3 states of N2, and the singlet-triplet splittings of oligoacenes. Our results show that both relative and absolute errors introduced by these cumulant truncations can be as small as 0.5 kcal mol-1 or less. At the same time, the amount of memory required is reduced from O(NA6) to O(NA5), where NA is the number of active orbitals. No additional regularization is needed to prevent the intruder state problem in the cumulant-truncated second-order driven similarity renormalization group multireference perturbation theory methods.

New Density Matrix Renormalization Group Approaches for Strongly Correlated Systems Coupled with Large Environments

Thanks to the high compression of the matrix product state (MPS) form of the wave function and the efficient site-by-site iterative sweeping optimization algorithm, the density matrix normalization group (DMRG) and its time-dependent variant (TD-DMRG) have been established as powerful computational tools in accurately simulating the electronic structure and quantum dynamics of strongly correlated molecules with a large number (10^1-2) of quantum degrees of freedom (active orbitals or vibrational modes). However, the quantitative characterization of the quantum many-body behaviors of realistic strongly correlated systems requires a further consideration of the interaction between the embedded active subsystem and the remaining correlated environment, e.g., a larger number (10^2-3) of external orbitals in electronic structure or infinite condensed-phase phononic modes in nucleus dynamics. To this end, we introduced three new post-DMRG and TD-DMRG approaches, namely (1) DMRG2sCI-MRCI and DMRG2sCI-ENPT by the reconstruction of selected configuration interaction (sCI) type of compact reference function from DMRG coefficients and the use of externally contracted MRCI (multireference configuration interaction) and Epstein-Nesbet perturbation theory (ENPT), without recourse to the expensive high order n-electron reduced density matrices (n-RDMs). (2) DMRG combined with RR-MRCI (renormalized residue-based MRCI), which improves the computational accuracy and efficiency of internally contracted (ic) MRCI by renormalizing the contracted bases with small-sized buffer environment(s) of a few external orbitals as probes based on quantum information theory. (3) HM (hierarchical mapping)-TD-DMRG in which a large environment is reduced to a small number of renormalized environmental modes (which accounts for the most vital system-environment interactions) through stepwise mapping transformation. These advances extend the efficacy of highly accurate DMRG/TD-DMRG computations to the quantitative characterization of the electronic structure and quantum dynamics in realistic strongly correlated systems coupled with large environments and are reviewed in this paper.

Electron Detachments of NbSin-/0 (n = 1-3) Clusters from Density Matrix Renormalization Group-CASPT2 Calculations

The electronic states of NbSi_n^-/0/+ (n = 1-3) clusters have been explored using the state-of-the-art DMRG-CASPT2 method with relatively large active spaces. The leading configurations, bond distances, vibrational frequencies, and relative energies of the low-lying states were identified. Electron detachment energies of the anionic cluster and ionization energies of the neutral clusters were reported at the DMRG-CASPT2 level. The ground states of the NbSi_n^-/0/+ (n = 1-3) clusters were predicted as the ³Δ, ⁴Π, and ⁵Π states of the linear NbSi^-/0/+, the ³A₂, ⁴B₁, and ³B₁ states of cyclic NbSi₂^-/0/+, and the ¹A', ²A', and ³A″ states of tetrahedral NbSi₃^-/0/+ isomers. The first feature in the photoelectron spectrum of NbSi^- was attributed to the transitions from the anionic ³Δ ground state to the neutral ⁴Π, ⁴Δ, and ⁴Φ states, whereas the second feature was assigned to the transitions to the neutral ²Δ, ²Σ⁺, and ²Φ states. The first band in the photoelectron spectrum of NbSi₃^- was ascribed to the transition from the anionic ¹A' ground state to the neutral 1²A' and 1²A″ states; the second band was attributed to the transitions to 2²A', 2²A″, and 3²A' states; and the third band was assigned to the transition to 3²A' states.

Kylin 1.0: An ab-initio density matrix renormalization group quantum chemistry program

The accurate evaluation of electron correlations is highly necessary for the proper descriptions of the electronic structures in strongly correlated molecules, ranging from bond-dissociating molecules, polyradicals, to large conjugated molecules and transition metal complexes. For this purpose, in this paper, a new ab-initio quantum chemistry program Kylin 1.0 for electron correlation calculations at various quantum many-body levels, including configuration interaction (CI), perturbation theory (PT), and density matrix renormalization group (DMRG), is presented. Furthermore, fundamental quantum chemistry methods such as Hartree-Fock self-consistent field (HF-SCF) and the complete active space SCF (CASSCF) are also implemented. The Kylin 1.0 program possesses the following features: (1) a matrix product operator (MPO) formulation-based efficient DMRG implementation for describing static electron correlation within a large active space composed of more than 100 orbitals, supporting both U 1 n × U 1 S z $$ \mathrm{U}{(1)}_{\mathrm{n}}\times \mathrm{U}{(1)}_{{\mathrm{S}}_{\mathrm{z}}} $$ and U 1 n × SU 2 S $$ \mathrm{U}{(1)}_{\mathrm{n}}\times \mathrm{SU}{(2)}_{\mathrm{S}} $$ symmetries; (2) an efficient second-order DMRG-self-consistent field (SCF) implementation; (3) an externally contracted multi-reference CI (MRCI) and Epstein-Nesbet PT with DMRG reference wave functions for including the remaining dynamic electron correlation outside the large active spaces. In this paper, we introduce the capabilities and numerical benchmark examples of the Kylin 1.0 program.

ipie: A Python-Based Auxiliary-Field Quantum Monte Carlo Program with Flexibility and Efficiency on CPUs and GPUs

We report the development of a python-based auxiliary-field quantum Monte Carlo (AFQMC) program, ipie, with preliminary timing benchmarks and new AFQMC results on the isomerization of [Cu₂O₂]²⁺. We demonstrate how implementations for both central and graphical processing units (CPUs and GPUs) are achieved in ipie. We show an interface of ipie with PySCF as well as a straightforward template for adding new estimators to ipie. Our timing benchmarks against other C++ codes, QMCPACK and Dice, suggest that ipie is faster or similarly performing for all chemical systems considered on both CPUs and GPUs. Our results on [Cu₂O₂]²⁺ using selected configuration interaction trials show that it is possible to converge the ph-AFQMC isomerization energy between bis(μ-oxo) and μ-η²:η² peroxo configurations to the exact known results for small basis sets with 10⁵-10⁶ determinants. We also report the isomerization energy with a quadruple-zeta basis set with an estimated error less than a kcal/mol, which involved 52 electrons and 290 orbitals with 10⁶ determinants in the trial wave function. These results highlight the utility of ph-AFQMC and ipie for systems with modest strong correlation and large-scale dynamic correlation.

Reliably assessing the electronic structure of cytochrome P450 on today's classical computers and tomorrow's quantum computers

Chemical simulation is one of the most promising applications for future quantum computers. It is thought that quantum computers may enable accurate simulation for complex molecules that are otherwise impossible to simulate classically; that is, it displays quantum advantage. To better understand quantum advantage in chemical simulation, we explore what quantum and classical resources are required to simulate a series of pharmaceutically relevant molecules. Using classical methods, we show that reliable classical simulation of these molecules requires significant resources and therefore is a promising candidate for quantum simulation. We estimate the quantum resources, both in overall simulation time and the size. The insights from this study pave the way for future quantum simulation of complex molecules.

An accurate assessment of how quantum computers can be used for chemical simulation, especially their potential computational advantages, provides important context on how to deploy these future devices. To perform this assessment reliably, quantum resource estimates must be coupled with classical computations attempting to answer relevant chemical questions and to define the classical algorithms simulation frontier. Herein, we explore the quantum computation and classical computation resources required to assess the electronic structure of cytochrome P450 enzymes (CYPs) and thus define a classical–quantum advantage boundary. This is accomplished by analyzing the convergence of density matrix renormalization group plus n-electron valence state perturbation theory (DMRG+NEVPT2) and coupled-cluster singles doubles with noniterative triples [CCSD(T)] calculations for spin gaps in models of the CYP catalytic cycle that indicate multireference character. The quantum resources required to perform phase estimation using qubitized quantum walks are calculated for the same systems. Compilation into the surface code provides runtime estimates to compare directly to DMRG runtimes and to evaluate potential quantum advantage. Both classical and quantum resource estimates suggest that simulation of CYP models at scales large enough to balance dynamic and multiconfigurational electron correlation has the potential to be a quantum advantage problem and emphasizes the important interplay between classical computations and quantum algorithms development for chemical simulation.

Revisiting the Orbital Energy-Dependent Regularization of Orbital-Optimized Second-Order Møller-Plesset Theory

Optimizing orbitals in the presence of electron correlation, as in orbital-optimized second-order Møller-Plesset perturbation theory (OOMP2), can remove artifacts associated with mean-field orbitals such as spin contamination and artificial symmetry-breaking. However, OOMP2 is known to suffer from divergent correlation energies in regimes of small orbital energy gaps. To address this issue, several approaches to amplitude regularization have been explored, with those featuring energy-gap-dependent regularizers appearing to be most transferable and physically justifiable. For instance, κ-OOMP2 was shown to address the energy divergence issue in, for example, bond-breaking processes while offering a significant improvement in accuracy for the W4-11 thermochemistry data set, and a parameter of κ = 1.45 was recommended. A more recent investigation of regularized MP2 with Hartree-Fock orbitals revealed that stronger regularization (i.e., smaller values of κ) than what had previously been recommended for κ-OOMP2 may offer huge improvements in certain cases such as noncovalent interactions while retaining a high level of accuracy for main-group thermochemistry data sets. In this study, we investigate the transferability of those findings to κ-OOMP2 and assess the implications of stronger regularization on the ability of κ-OOMP2 to diagnose strong static correlation. We found similar results using κ-OOMP2 for several main-group thermochemistry, barrier height, and noncovalent interaction data sets including both closed shell and open shell species. However, stronger regularization yielded substantially higher accuracy for open-shell transition-metal (TM) thermochemistry and is necessary to provide qualitatively correct spin symmetry breaking behavior for several large and electrochemically relevant TM systems. We therefore find a single κ value insufficient to treat all systems using κ-OOMP2.

Efficient Implementation of Block-Correlated Coupled Cluster Theory Based on the Generalized Valence Bond Reference for Strongly Correlated Systems

An optimized implementation of block-correlated coupled cluster theory based on the generalized valence bond wave function (GVB-BCCC) for the singlet ground state of strongly correlated systems is presented. The GVB-BCCC method with two-pair correlation (GVB-BCCC2b) or up to three-pair correlation (GVB-BCCC3b) will be the focus of this work. Three major techniques have been adopted to dramatically accelerate GVB-BCCC2b and GVB-BCCC3b calculations. First, the GVB-BCCC2b and GVB-BCCC3b codes are noticeably optimized by removing redundant calculations. Second, independent amplitudes are identified by constraining excited configurations to be pure singlet states and only independent amplitudes need to be solved. Third, an incremental updating scheme for the amplitudes in solving the GVB-BCCC equations is adopted. With these techniques, accurate GVB-BCCC3b calculations are now accessible for systems with relatively large active spaces (50 electrons in 50 orbitals) and GVB-BCCC2b calculations are affordable for systems with much larger active spaces. We have applied GVB-BCCC methods to investigate three typical kinds of systems: polyacenes, pentaprismane, and [Cu₂O₂]²⁺ isomers. For polyacenes, we demonstrate that GVB-BCCC3b can capture more than 94% of the total correlation energy even for 12-acene with 50 π electrons. For the potential energy curve of simultaneously stretching 15 C-C bonds in pentaprismane, our calculations show that the GVB-BCCC3b results are quite close to the results from the density matrix renormalization group (DMRG) over the whole range. For two dinuclear copper oxide isomers, their relative energy predicted by GVB-BCCC3b is also in good accord with the DMRG result. All calculations show that the inclusion of three-pair correlation in GVB-BCCC is critical for accurate descriptions of strongly correlated systems.

Multiconfiguration Pair-Density Functional Theory Calculations of Iron(II) Porphyrin: Effects of Hybrid Pair-Density Functionals and Expanded RAS and DMRG Active Spaces on Spin-State Orderings

Iron(II) porphyrins play critical roles in enzymes and synthetic catalysts. Computationally determining the spin-state ordering for even the unsubstituted iron(II) porphyrin (FeP) is challenging due to its large size. Multiconfiguration pair-density functional theory (MC-PDFT), a method capable of accurately capturing correlation with lower cost than comparably accurate methods, was previously used to predict a triplet ground state for FeP across a wide range of active spaces up to (34e, 35o). The purpose of this present MC-PDFT study is to determine the effects of including nonlocal exchange in the energy calculation and of using a larger active space size [DMRG(40e, 42o) and RAS(40, 2, 2; 16, 6, 20)] on the calculated FeP spin-state ordering. The recently developed hybrid MC-PDFT method, which uses a weighted average of the MC-PDFT energy and the energy expectation value of the reference wave function, is applied with a weight of the reference wave function energy of λ. We find that increasing λ stabilizes the quintet relative to the triplets. The hybrid tPBE0 functional (tPBE with λ set to 0.25) consistently predicts a triplet ground state with the quintet lying above by 0.10-0.16 eV, depending on the reference wave function. These values are particularly interesting in light of tPBE0's very strong performance for a diverse set of other systems.

Comprehensive Theoretical View of the [Cu2O2] Side-on-Peroxo-/Bis-μ-Oxo Equilibria

Coupled binuclear copper (CBC) sites are employed by many metalloenzymes to catalyze a broad set of biochemical transformations. Typically, the CBC catalytic sites are activated by the O2 molecule to form various [Cu2O2] reactive species. This has also inspired synthesis and development of various biomimetic inorganic complexes featuring the CBC core. From theoretical perspective, the [Cu2O2] reactivity often hinges on the side-on-peroxo-dicopper(II) (P) vs. bis-μ-oxo-dicopper(III) (O) isomerism - an equilibrium that has become almost iconic in theoretical bioinorganic chemistry. Herein, we present a comprehensive calibration and evaluation of the performance of various composite computational protocols available in contemporary computational chemistry, involving coupled-cluster and multi-reference (relativistic) wave function methods, popular density functionals and solvation models. Starting with the well-studied reference [Cu2O2(NH3)6]2+ system, we compared the performance of electronic structure methods and discussed the relativistic effects. This allowed us to select several 'calibrated' DFT functionals that can be conveniently employed to study ten experimentally well-characterized [Cu2O2] inorganic systems. We mostly correctly predicted the lowest-energy structures (P vs. O) of the studied systems. In addition, we present calibration of the used electronic structure methods for prediction of the spectroscopic features of the [Cu2O2] core.

Post-Density Matrix Renormalization Group Methods for Describing Dynamic Electron Correlation with Large Active Spaces

The ab initio density matrix renormalization group (DMRG) method has been well-established and has become one of the most accurate numerical methods for the precise electronic structure solution of large active spaces. In the past few years, to capture the missing dynamic correlation, various post-DMRG approaches have been proposed through the combination of DMRG and multireference quantum chemical methods or density functional theory. With this in mind, this work provides a brief overview of ab initio DMRG principles and the new developments within post-DMRG methods. For clarity, post-DMRG methods are classified into two main categories depending on whether high-order n-electron reduced density matrices are used, and their merits and disadvantages are properly discussed. Finally, we conclude by discussing unsolved bottlenecks and giving development perspectives of post-DMRG approaches, which are expected to yield quantitative descriptions of complex electronic structures in large strongly correlated molecules and materials.

A DMRG/CASPT2 Investigation of Metallocorroles: Quantifying Ligand Noninnocence in Archetypal 3d and 4d Element Derivatives

Hybrid density functional theory (B3LYP) and density matrix renormalization group (DMRG) theory have been used to quantitatively compare the degree of ligand noninnocence (corrole radical character) in seven archetypal metallocorroles. The seven complexes, in decreasing order of corrole noninnocent character, are Mn[Cor]Cl > Fe[Cor]Cl > Fe[Cor](NO) > Mo[Cor]Cl2 > Ru[Cor](NO) ≈ Mn[Cor]Ph ≈ Fe[Cor]Ph ≈ 0, where [Cor] refers to the unsubstituted corrolato ligand. DMRG-based second-order perturbation theory calculations have also yielded detailed excited-state energetics data on the compounds, shedding light on periodic trends involving middle transition elements. Thus, whereas the ground state of Fe[Cor](NO) (S = 0) is best described as a locally S = 1/2 {FeNO}7 unit antiferromagnetically coupled to a corrole A' radical, the calculations confirm that Ru[Cor](NO) may be described as simply {RuNO}6-Cor3-, that is, having an innocent corrole macrocycle. Furthermore, whereas the ferromagnetically coupled S = 1{FeNO}7-Cor•2- state of Fe[Cor](NO) is only ∼17.5 kcal/mol higher than the S = 0 ground state, the analogous triplet state of Ru[Cor](NO) is higher by a far larger margin (37.4 kcal/mol) relative to the ground state. In the same vein, Mo[Cor]Cl2 exhibits an adiabatic doublet-quartet gap of 36.1 kcal/mol. The large energy gaps associated with metal-ligand spin coupling in Ru[Cor](NO) and Mo[Cor]Cl2 reflect the much greater covalent character of 4d-π interactions relative to analogous interactions involving 3d orbitals. As far as excited-state energetics is concerned, DMRG-CASPT2 calculations provide moderate validation for hybrid density functional theory (B3LYP) for qualitative purposes, but underscore the possibility of large errors (>10 kcal/mol) in interstate energy differences.

An efficient implementation of the NEVPT2 and CASPT2 methods avoiding higher-order density matrices

A factorization of the matrix elements of the Dyall Hamiltonian in N-electron valence state perturbation theory allowing their evaluation with a computational effort comparable to the one needed for the construction of the third-order reduced density matrix at the most is presented. Thus, the computational bottleneck arising from explicit evaluation of the fourth-order density matrix is avoided. It is also shown that the residual terms arising in the case of an approximate complete active space configuration interaction solution and containing even the fifth-order density matrix for two excitation classes can be evaluated with little additional effort by choosing again a favorable factorization of the corresponding matrix elements. An analogous argument is also provided for avoiding the fourth-order density matrix in complete active space second-order perturbation theory. Practical calculations indicate that such an approach leads to a considerable gain in computational efficiency without any compromise in numerical accuracy or stability.

Analytic Energy Gradients for the Driven Similarity Renormalization Group Multireference Second-Order Perturbation Theory

We derive analytic energy gradients of the driven similarity renormalization group (DSRG) multireference second-order perturbation theory (MRPT2) using the method of Lagrange multipliers. In the Lagrangian, we impose constraints for a complete-active-space self-consistent-field reference wave function and the semicanonical orthonormal molecular orbitals. Solving for the associated Lagrange multipliers is found to share the same asymptotic scaling of a single DSRG-MRPT2 energy computation. A pilot implementation of the DSRG-MRPT2 analytic gradients is used to optimize the geometry of the singlet and triplet states of p-benzyne. The equilibrium bond lengths and angles are similar to those computed via other MRPT2s and Mukherjee's multireference coupled cluster theory. An approximate DSRG-MRPT2 method that neglects the contributions of the three-body density cumulant is found to introduce negligible errors in the geometry of p-benzyne, lending itself to a promising low-cost approach for molecular geometry optimizations using large active spaces.

Stochastic Generalized Active Space Self-Consistent Field: Theory and Application

An algorithm to perform stochastic generalized active space calculations, Stochastic-GAS, is presented, that uses the Slater determinant based FCIQMC algorithm as configuration interaction eigensolver. Stochastic-GAS allows the construction and stochastic optimization of preselected truncated configuration interaction wave functions, either to reduce the computational costs of large active space wave function optimizations, or to probe the role of specific electron correlation pathways. As for the conventional GAS procedure, the preselection of the truncated wave function is based on the selection of multiple active subspaces while imposing restrictions on the interspace excitations. Both local and cumulative minimum and maximum occupation number constraints are supported by Stochastic-GAS. The occupation number constraints are efficiently encoded in precomputed probability distributions, using the precomputed heat bath algorithm, which removes nearly all runtime overhead of GAS. This strategy effectively allows the FCIQMC dynamics to a priori exclude electronic configurations that are not allowed by GAS restrictions. Stochastic-GAS reduced density matrices are stochastically sampled, allowing orbital relaxations via Stochastic-GASSCF, and direct evaluation of properties that can be extracted from density matrices, such as the spin expectation value. Three test case applications have been chosen to demonstrate the flexibility of Stochastic-GAS: (a) the Stochastic-GASSCF [5·(6, 6)] optimization of a stack of five benzene molecules, that shows the applicability of Stochastic-GAS toward fragment-based chemical systems; (b) an uncontracted stochastic MRCISD calculation that correlates 96 electrons and 159 molecular orbitals, and uses a large (32, 34) active space reference wave function for an Fe(II)-porphyrin model system, showing how GAS can be applied to systematically recover dynamic electron correlation, and how in the specific case of the Fe(II)-porphyrin dynamic correlation further differentially stabilizes the ³E_g over the ⁵A_1g spin state; (c) the study of an Fe₄S₄ cluster's spin-ladder energetics via highly truncated stochastic-GAS [4·(5, 5)] wave functions, where we show how GAS can be applied to understand the competing spin-exchange and charge-transfer correlating mechanisms in stabilizing different spin-states.

Computational Study of CO2 Reduction Catalyzed by Iron(I) Complex at Different Spin States: Cooperativity of Hydrogen Bonding and Auxiliary Group Effect

To explore alternative approaches to the CO₂ reduction to formate and provide an insight into the spin state effect on the CO₂ reduction, we theoretically designed a kind of low-valence iron(I) model complex, whose doublet, quartet, and sextet states are denoted as ² Fe(I), ⁴ Fe(I), and ⁶ Fe(I), respectively. This complex is featured with an iron(I) center, which bonds to a 1,2-ethanediamine (en) and a 2-hydroxy-biphenyl group. Reaction mechanisms for the CO₂ reduction to formate catalyzed by this iron(I) model complex were explored using density functional theory (DFT) computations. Studies showed that the univalent iron(I) compound can efficiently fix and activate a CO₂ molecule, whereas its oxidized forms with trivalent iron(III) or bivalent iron(II) cannot activate CO₂. For the iron(I) compound, it was found that the lowest spin state ² Fe(I) is the most favorable for the CO₂ reduction as the reactions barriers involving ² Fe(I), ⁴ Fe(I), and ⁶ Fe(I) are 25.6, 37.2, and 35.9 kcal/mol, respectively. Yet, a photosensitizer-free visible-light-mediated high-low spin shift from ⁴ Fe(I) and ⁶ Fe(I) to ² Fe(I) is likely through the reverse intersystem crossing (RIC) because the ⁴ Fe(I) and ⁶ Fe(I) compounds have strong absorption in the visible-light range. Notably, the synergistic interaction between the hydrogen bonding from the auxiliary hydroxyl group in the 2-hydroxy-biphenyl moiety to CO₂ and an intermediate five-membered ring promotes the proton transfer, leading to the formation of the -COOH moiety from CO₂ and the Fe-O bond. With the addition of H₂, one H₂ molecule is split by the Fe-O bond and thus serves as H atom sources for both the CO₂ reduction and the recovery of the auxiliary hydroxyl group. The present theoretical study provides a novel solution for the challenging CO₂ reduction, which calls for further experimental verifications.

Density Matrix Renormalization Group with Dynamical Correlation via Adiabatic Connection

The quantum chemical version of the density matrix renormalization group (DMRG) method has established itself as one of the methods of choice for calculations of strongly correlated molecular systems. Despite its great ability to capture strong electronic correlation in large active spaces, it is not suitable for computations of dynamical electron correlation. In this work, we present a new approach to the electronic structure problem of strongly correlated molecules, in which DMRG is responsible for a proper description of the strong correlation, whereas dynamical correlation is computed via the recently developed adiabatic connection (AC) technique which requires only up to two-body active space reduced density matrices. We report the encouraging results of this approach on typical candidates for DMRG computations, namely, n-acenes (n = 2 → 7), Fe(II)-porphyrin, and the Fe₃S₄ cluster.

Machine Learning-Assisted Selection of Active Spaces for Strongly Correlated Transition Metal Systems

Active space quantum chemical methods could provide very accurate description of strongly correlated electronic systems, which is of tremendous value for natural sciences. The proper choice of the active space is crucial but a nontrivial task. In this article, we present a neural network-based approach for automatic selection of active spaces, focused on transition metal systems. The training set has been formed from artificial systems composed of one transition metal and various ligands, on which we have performed the density matrix renormalization group and calculated the single-site entropy. On the selected set of systems, ranging from small benchmark molecules up to larger challenging systems involving two metallic centers, we demonstrate that our machine learning models could predict the active space orbitals with reasonable accuracy. We also tested the transferability on out-of-the-model systems, including bimetallic complexes and complexes with ligands, which were not involved in the training set. Also, we tested the correctness of the automatically selected active spaces on a Fe(II)-porphyrin model, where we studied the lowest states at the DMRG level and compared the energy difference between spin states or the energy difference between conformations of ferrocene with recent studies.

Modern multireference methods and their application in transition metal chemistry

Transition metal chemistry is a challenging playground for quantum chemical methods owing to the simultaneous presence of static and dynamic electron correlation effects in many systems. Wavefunction based multireference (MR) methods constitute a physically sound and systematically improvable Ansatz to deal with this complexity but suffer from some conceptual difficulties and high computational costs. The latter problem partially arises from the unfavorable scaling of the Full Configuration Interaction (Full-CI) problem which in the majority of MR methods is solved for a subset of the molecular orbital space, the so-called active space. In the last years multiple methods such as modern variants of selected CI, Full-CI Quantum Monte Carlo (FCIQMC) and the density matrix renormalization group (DMRG) have been developed that solve the Full-CI problem approximately for a fraction of the computational cost required by conventional techniques thereby significantly extending the range of applicability of modern MR methods. This perspective review outlines recent advancements in the field of MR electronic structure methods together with the resulting chances and challenges for theoretical studies in the field of transition metal chemistry. In light of its emerging importance a special focus is put on the selection of adequate active spaces and the concomitant development of numerous selection aides in recent years.

Approximations of density matrices in N-electron valence state second-order perturbation theory (NEVPT2). II. The full rank NEVPT2 (FR-NEVPT2) formulation

In Paper I, the performances of pre-screening (PS), extended PS (EPS), and cumulant (CU) approximations to the fourth-order density matrix were examined in the context of second-order N-electron valence state perturbation theory (NEVPT2). It has been found that the CU, PS, and even EPS approximations with loose thresholds may introduce intruder states. In the present work, the origin of these "false intruder" states introduced by approximated density matrices is discussed. Canonical NEVPT2 implementations employ a rank reduction trick. By analyzing its residual error, we find that the omission of the rank reduction leads to a more stable multireference perturbation theory for incomplete active space reference wave functions. Such a full rank (FR)-NEVPT2 formulation is equivalent to the conventional NEVPT2 method for the complete active space self-consistent field/complete active space configuration interaction reference wave function. A major drawback of the FR-NEVPT2 formulation is the necessity of the fifth-order density matrix. To avoid the construction of the high-order density matrices, the combination of the FR-NEVPT2 with the CU approximation is studied. However, we find that the CU approximation remains problematic as it still introduces intruder states. The question of how to robustly and efficiently perform internally contracted multireference perturbation theories with approximate densities remains a challenging field of investigation.

Approximations of density matrices in N-electron valence state second-order perturbation theory (NEVPT2). I. Revisiting the NEVPT2 construction

Over the last decade, the second-order N-electron valence state perturbation theory (NEVPT2) has developed into a widely used multireference perturbation method. To apply NEVPT2 to systems with large active spaces, the computational bottleneck is the construction of the fourth-order reduced density matrix. Both its generation and storage become quickly problematic beyond the usual maximum active space of about 15 active orbitals. To reduce the computational cost of handling fourth-order density matrices, the cumulant approximation (CU) has been proposed in several studies. A more conventional strategy to address the higher-order density matrices is the pre-screening approximation (PS), which is the default one in the ORCA program package since 2010. In the present work, the performance of the CU, PS, and extended PS (EPS) approximations for the fourth-order density matrices is compared. Following a pedagogical introduction to NEVPT2, contraction schemes, as well as the approximations to density matrices, and the intruder state problem are discussed. The CU approximation, while potentially leading to large computational savings, virtually always leads to intruder states. With the PS approximation, the computational savings are more modest. However, in conjunction with conservative cutoffs, it produces stable results. The EPS approximation to the fourth-order density matrices can reproduce very accurate NEVPT2 results without any intruder states. However, its computational cost is not much lower than that of the canonical algorithm. Moreover, we found that a good indicator of intrude states problems in any approximation to high order density matrices is the eigenspectra of the Koopmans matrices.

The Shape of Full Configuration Interaction to Come

We present a Perspective on what the future holds for full configuration interaction (FCI) theory, with an emphasis on conceptual rather than technical details. Upon revisiting the early history of FCI, a number of its key contemporary approximations are compared on as equal a footing as possible, using a recent blind challenge on the benzene molecule as a testbed [Eriksen et al., J. Phys. Chem. Lett., 2020 11, 8922]. In the process, we review the scope of applications for which FCI continues to prove indispensable, and the required traits in terms of robustness, efficacy, and reliability its modern approximations must satisfy are discussed. We close by conveying a number of general observations on the merits offered by the state-of-the-art alongside some of the challenges still faced to this day. While the field has altogether seen immense progress over the years-the past decade, in particular-it remains clear that our community as a whole has a substantial way to go in enhancing the overall applicability of near-exact electronic structure theory for systems of general composition and increasing size.

Ultrafast processes: coordination chemistry and quantum theory

Coordination compounds, characterized by fascinating and tunable electronic properties, are capable of binding easily to proteins, polymers, wires and DNA. Upon irradiation, these molecular systems develop functions finding applications in solar cells, photocatalysis, luminescent and conformational probes, electron transfer triggers and diagnostic or therapeutic tools. The control of these functions is activated by the light wavelength, the metal/ligand cooperation and the environment within the first picoseconds (ps). After a brief summary of the theoretical background, this perspective reviews case studies, from 1st row to 3rd row transition metal complexes, that illustrate how spin-orbit, vibronic coupling and quantum effects drive the photophysics of this class of molecules at the early stage of the photoinduced elementary processes within the fs-ps time scale range.

Portably parallel construction of a configuration-interaction wave function from a matrix-product state using the Charm++ framework

The construction of configuration-interaction (CI) expansions from a matrix product state (MPS) involves numerous matrix operations and the skillful sampling of important configurations in a large Hilbert space. In this work, we present an efficient procedure for constructing CI expansions from MPS employing the parallel object-oriented Charm++ programming framework, upon which automatic load-balancing and object migrating facilities can be employed. This procedure was employed in the MPS-to-CI utility (Moritz et al., J. Chem. Phys. 2007, 126, 224109), the sampling-reconstructed complete active-space algorithm (SR-CAS, Boguslawski et al., J. Chem. Phys. 2011, 134, 224101), and the entanglement-driven genetic algorithm (EDGA, Luo et al., J. Chem. Theory Comput. 2017, 13, 4699). It enhances productivity and allows the sampling programs to evolve to their population-expansion versions, for example, EDGA with population expansion (PE-EDGA). Further, examples of 1,2-dioxetanone and firefly dioxetanone anion (FDO^- ) molecules demonstrated the following: (a) parallel efficiencies can be persistently improved by simply by increasing the proportions of the asynchronous executions and (b) a sampled CAS-type CI wave function of a bi-radical-state FDO^- molecule utilizing the full valence (30e,26o) active space can be constructed within a few hours with using thousands of cores.

Towards Near-Exact Solutions of Molecular Electronic Structure: Full Coupled-Cluster Reduction with a Second-Order Perturbative Correction

We introduce a new augmented adaptation of the recently developed full coupled-cluster reduction (FCCR) with a second-order perturbative correction, abbreviated as FCCR(2). FCCR is a selected coupled-cluster expansion aimed at optimally reducing the excitation manifold and commutator expansions for high-rank excitations for obtaining accurate solutions of the electronic Schödinger equation in a size-extensive manner. The present FCCR(2) enables estimating the residual correlation of FCCR by the second-order perturbative correction E⁽²⁾ from the complementary space of the FCCR projection manifold. The linear relationship between E⁽²⁾ and the energy of FCCR(2) allows accurate estimates of near-exact energies for a wide variety of molecules with strong electron correlation. The potential of the method is demonstrated using challenging cases, the ground-state electronic energy of the benzene molecule in equilibrium and stretched geometries, and the isomerization energy of the transition metal complex [Cu(NH₃)]₂O₂²⁺.

Molecular Photochemistry: Recent Developments in Theory

Photochemistry is a fascinating branch of chemistry that is concerned with molecules and light. However, the importance of simulating light-induced processes is reflected also in fields as diverse as biology, material science, and medicine. This Minireview highlights recent progress achieved in theoretical chemistry to calculate electronically excited states of molecules and simulate their photoinduced dynamics, with the aim of reaching experimental accuracy. We focus on emergent methods and give selected examples that illustrate the progress in recent years towards predicting complex electronic structures with strong correlation, calculations on large molecules, describing multichromophoric systems, and simulating non-adiabatic molecular dynamics over long time scales, for molecules in the gas phase or in complex biological environments.

Modern quantum chemistry with [Open]Molcas

MOLCAS/OpenMolcas is an ab initio electronic structure program providing a large set of computational methods from Hartree-Fock and density functional theory to various implementations of multiconfigurational theory. This article provides a comprehensive overview of the main features of the code, specifically reviewing the use of the code in previously reported chemical applications as well as more recent applications including the calculation of magnetic properties from optimized density matrix renormalization group wave functions.

CASSCF with Extremely Large Active Spaces Using the Adaptive Sampling Configuration Interaction Method

The complete active space self-consistent field (CASSCF) method is the principal approach employed for studying strongly correlated systems. However, exact CASSCF can only be performed on small active spaces of ∼20 electrons in ∼20 orbitals due to exponential growth in the computational cost. We show that employing the Adaptive Sampling Configuration Interaction (ASCI) method as an approximate Full CI solver in the active space allows CASSCF-like calculations within chemical accuracy (<1 kcal/mol for relative energies) in active spaces with more than ∼50 active electrons in ∼50 active orbitals, significantly increasing the sizes of systems amenable to accurate multiconfigurational treatment. The main challenge with using any selected CI-based approximate CASSCF is the orbital optimization problem; they tend to exhibit large numbers of local minima in orbital space due to their lack of invariance to active-active rotations (in addition to the local minima that exist in exact CASSCF). We highlight methods that can avoid spurious local extrema as a practical solution to the orbital optimization problem. We employ ASCI-SCF to demonstrate a lack of polyradical character in moderately sized periacenes with up to 52 correlated electrons and compare against heat-bath CI on an iron porphyrin system with more than 40 correlated electrons.

Transition Metal Spin-State Energetics by MC-PDFT with High Local Exchange

The energetics of the spin states of transition metal complexes have been explored with a variety of electronic structure methods, but the calculations require a compromise between accuracy and affordability. In this work, the spin splittings of several iron complexes are studied with multiconfiguration pair-density functional theory (MC-PDFT). The results are compared to previously published results obtained by complete active space second-order perturbation theory (CASPT2) and CASPT2 with coupled-cluster semicore correlation (CASPT2/CC). In contrast to CASPT2's systematic overstabilization of high-spin states with respect to the CASPT2/CC reference, MC-PDFT with the tPBE on-top functional understabilizes high-spin states. This systematic understabilization is largely corrected by revising the exchange and correlation contributions to the on-top functional using the high local-exchange approximation (tPBE-HLE). Moreover, tPBE-HLE correctly predicts the spin of the ground state in most cases, while CASPT2 incorrectly predicts high-spin ground states in all cases. This is encouraging for practical work because tPBE and tPBE-HLE are faster than CASPT2 by a factor of 50 even in a moderately sized example.

Mononuclear Fe in N-doped carbon: computational elucidation of active sites for electrochemical oxygen reduction and oxygen evolution reactions

First principles simulations show that in Fe and N co-doped carbon, Fe coordination controls the activity for oxygen reduction and oxygen evolution reactions, and that including the electrostatic potential has a major influence at high potential.

Energetics of non-heme iron reactivity: can ab initio calculations provide the right answer?

We use a variety of computational methods to characterize and compare the hydrogen atom transfer (HAT) and epoxidation reaction pathways for oxidation of cyclohexene by an iron(iv)-oxo complex.