Multireference Methods are Realistic and Useful Tools for Modeling Catalysis

Highly correlated systems, in particular those that include transition metals, are ubiquitous in catalysis. The significant static correlation found in such systems is often poorly accounted for using Kohn Sham density functional theory methods, as they are single determinantal in nature. Applications to catalysis of more rigorous and appropriate multiconfigurational methods have been reported in select instances, but their use remains rare. We discuss obstacles that hinder the routine application of multireference (MR) wave function theoretical calculations to catalytic systems and the current state of the art with respect to removing those


Introduction
Computational chemistry methods are an indispensable tool for the rational design of catalysts, offering a far more efficient alternative to "trial and error"-based heuristic approaches. [1] By identifying the reaction mechanism(s) as well as structureactivity and property-activity relationships involving computed descriptors, they can be employed to screen, in silico, large sets of chemical systems for new catalysts candidates. [2] The results of the screening can then be used to direct synthetic efforts, saving both time and physical and chemical resources. [3] We assert that we have reached an era in which computational catalysis does not run behind experiment, but rather it works in close partnership to drive design.
Properties of reaction intermediates that are quantitatively correlated with a property of the transition state (e. g., the free energy of activation) are called descriptors. Descriptors themselves can vary widely in character: they may be associated with structural features, distributions of charge, or energetics along the reaction pathway(s). In the final category, while it is a given that computed free energies of activation associated with a set of related rate-or stereo-determining reaction steps will directly correlate with activity, it may also be the case that the relative energies of related reaction intermediates will correlate with those free energies of activation, and hence with activity. Finding and optimizing reaction intermediates is typically much less timeconsuming than the analogous task for transition-state structures, so recourse to property-activity relationships involving the former can make the discovery of new catalysts faster and more costeffective.

DFT Methods
Among computational methods today, Kohn-Sham density functional theory (KS-DFT) remains the most commonly employed for catalytic studies. This is due to three key advantages it often enjoys compared to other methods: (1) it is user-friendly to apply (albeit perhaps misleadingly); [4] (2) it is relatively undemanding of computational resources; (3) modern functionals have evolved to achieve near chemical accuracy [5] in systems amenable to characterization by a single Kohn-Sham determinant. [6] Although the potential of DFT (in the following we will use DFT as a shorthand for KS-DFT, noting that the vast majority of DFT applications adopt the Kohn-Sham formalism) remains far from having been fully exploited, the approximations on which DFT methods rely (because the exact exchange-correlation functional is unknown) make DFT for systems containing unpaired electrons (open-shell systems) at least problematic. [6c,7] These systems often have multi-configurational character, i. e., more than a single electronic configuration is important for the accurate description of their electronic structure. Such multiconfigurational character is common when there is a narrow separation in energies of nominally frontier orbitals, which is typically the case for transition-metal-containing species lacking a complete set of frontier d electrons. Indeed, it is this very feature of flexibility in electronic structure that makes transition metals so readily tuned for catalysis, but it also complicates the accurate representation of their electronic structure, [8] and this becomes particularly acute in the case of transition-state structures, which by their nature include partial bonds and substantial charge delocalization. [9] As noted above, KS-DFT represents a system's electron density using only a single Slater determinant (sometimes call the KS "wave function", although this is somewhat sloppy and misleading since the energy is computed from the density and a particular functional, not from application of the Hamiltonian operator to the KS determinant). For systems characterized by significant multiconfigurational character, the suitability of DFT modeling needs to be evaluated case-by-case through comparisons among different functionals (which often show enormous variation for the most challenging cases) and/or through benchmarks with experiment or higher-level computations. [1,10] Much has been learned about this situation from many years of studies. Many compounds incorporating one or more transition metals are characterized by two or more unpaired electrons. Single-determinantal solutions, as from KS-DFT, are usually not eigenfunctions of the S 2 operator for spin multiplicities (2S + 1) lower than the maximum multiplicity. [1] As a result, DFT methods can behave poorly when used to compute spin ladder energetics for open shell systems, [7] and indeed such analysis can be one means to assess the suitability of a given DFT functional for a particular chemical problem where such a spin ladder exists. [4] Reasonable energetics for open-shell systems can sometimes be calculated from KS-DFT using so called "broken symmetry" solutions, although these are neither spin eigenfunctions nor do they generally produce correct spin densities. [1,10e,11] In such situations, even though electronic energies may be computed with surprisingly good accuracy, other features of the system may not be physically meaningful. Other approximate strategies include adoption of an ad hoc percentage of Hartree-Fock exchange [12] or of a Hubbard term [13] tuned to reproduce reference values. [11a] For many open-shell systems, then -as noted above -DFT is only seemingly a user-friendly method: while it will straightforwardly optimize a KS determinant for a given structure, considerable experience is required to assess whether the Christopher J. Cramer As senior vice president and chief research officer for Underwriters Laboratories Inc., Dr. Christopher J. Cramer oversees a broad portfolio of projects all sharing a focus on fundamental and applied safety science. His guidance enhances the breadth and impact of the organization's capabilities as a world-class safety science research and standards development institution, and follows an academic career as a theoretical chemist at the University of Minnesota that culminated in his service as the institution's vice president for research. Cramer continues to hold the title of Distinguished McKnight and University Teaching Professor Emeritus.

Laura Gagliardi is the Richard and Kathy
Leventhal Professor at the University of Chicago. She received her Ph.D. degree in theoretical chemistry from the University of Bologna in 1997, and then was a postdoc for two years at Cambridge University. She began her independent career as an assistant professor in 2002 at the University of Palermo, moving to the University of Geneva as associate professor in 2005. In 2009, she moved to the University of Minnesota as a professor, where she remained prior to her move to Chicago in 2020. Laura is a quantum chemist who works on the development of electronic structure methods and their use for understanding complex energy-relevant chemical systems and materials. She is an expert in porous materials, homogeneous and heterogeneous catalysis, heavy-element chemistry and magnetic systems.
results from that functional are adequate for the chemistry under examination.
That said, DFT methods can be quite successfully applied across a broad range of systems, and DFT is used routinely for the production of large training data sets for machine learning in catalysis. In particular, meta-GGA functionals [14] have been shown generally to offer a good compromise for many systems, offering good accuracy not only for reaction barriers [11b,15] but also for more subtle features like the vibrational shift of molecular probes. [16] And, even if the "right answer for the wrong reason" rule applies to the particular application, [17] it may be offset by the robust character of catalyst descriptors when looking for correlations between computed properties and experimental activities. [18] Thus, DFT methods will certainly continue to play an important role in the design of new catalysts (directly or in combination with machine learning) and in the interpretation of experiments.
However, in this article we are especially interested in those multiconfigurational systems that prove pathological to DFT, and for which recourse to a multireference (MR) treatment is essentially required. With that as a focus, we will briefly review the most important MR methods adopted in catalytic studies (i. e., active-space-based and excitationbased truncation), [1] highlighting those factors that hinder their routine application. We will also summarize current strategies to overcome these obstacles, and we will explore the feasibility of the routinely use of MR methods for catalysis.

MR Methods
The most accurate wave function for a general chemical system, in a given basis set, and neglecting relativistic effects, is the full configuration interaction (FCI) wave function. In an FCI wave function, one includes all the configuration state functions that can be formed within a given spin symmetry considering all of the electrons distributed across all of the orbitals. As is routinely noted in textbooks, FCI calculations are not computationally feasible except for systems having a very small number of electrons, so approximations are employed. A useful definition of electron correlation energy is the difference between the FCI energy (even if it is not known) and the single-reference Hartree-Fock energy. [19] For pedagogical reasons, electron correlation is typically classified as deriving from two contributions, "dynamic" and "nondynamic" (also commonly referred to as "static"), although a practical division may not be straightforward. [20] Static correlation is associated with the multireference character of an electronic system. It is called "static" because it is not directly associated with the dynamical character of electron-electron interaction dynamics (which are ignored in mean-field approximations), but instead relates to the inadequacy of using only a single Slater determinant to describe the wave function. In the expansion of the MR wave function as a sum of configurational states, the relative weights of alternative configuration state functions as they contribute to the overall wave function are a measure of the importance of static correlation. Jiang et al. [21] have studies other diagnostics as well that quantify the importance of static correlation relative to a single-reference description for 3d transition-metal compounds, including spin contamination, as alluded to in the DFT discussion above.
Dynamic correlation, by contrast, accounts for the difference between approximate (mean-field, in the case of Hartree-Fock) and accurate electron-electron interactions, encompassing such important effects as short-range exchange repulsion and long-range dispersion interactions. This pedagogical distinction will prove useful in further discussion.
Among multiconfiguration self-consistent-field (MCSCF) methods, the complete active space self-consistent field (CASSCF) method [22] involves a FCI of n electrons within a limited number N of orbitals that form the active space. By convention, CASSCF active spaces are labeled as (n, N). Besides the active space in which the FCI takes place, there are additional orbitals outside the active space that are optimized under the constraint of fixed occupancy (0 or 2) for all configurations included in the CASSCF wave function. In general, active space orbitals will be those spanning the frontier of the occupied/virtual orbital space, while inactive occupied orbitals will be lower in the energy manifold and inactive virtual orbitals will be higher in the energy manifold.
CASSCF is designed mainly to capture static correlation, although increasing amounts of dynamic correlation are accounted for with growing size of the active space (as must be the case, since having every orbital in the space is the FCI limit, which captures all correlation). However, recovery of the dynamic correlation (also called external correlation) [23] converges slowly with respect to active space size, [24] so it is generally more efficient to pursue post-CASCCF treatments to address this quantity. One such example is perturbation theory to second order, i. e., complete active space second-order perturbation theory (CASPT2), [25] which is often used as a "gold-standard" benchmark (although it has some flaws as well). [26] Secondorder n-electron valence state perturbation theory (NEVPT2) [27] is the most common alternative to CASPT2 and it has a comparable computational cost. [28] Computational cost currently limits the utilization of CASSCF and CASPT2 methods, as it increases exponentially with active space size. This is a special concern for transition-metal-based catalysis because active spaces generally need to include at least the d orbital manifold of the metal(s) as well as additional orbitals associated with making and breaking bonds of substrates. To date, a (20,20) active space appears to be the largest that has been successfully employed, using massive parallelism. [29] On conventional high-performance computers, active spaces tend not to exceed (14,15) and within the context of that limit the size of the chemical model [30] and/or the basis set [31] is sometimes reduced.
Some approximations to CASSCF and CASPT2 scale more favorably with active space size (Section 3.1). Nevertheless, geometry optimization in particular is out of reach for catalytic systems when using MR methods. Some exceptions include Wang et al. [28b] reporting optimization at the CASSCF level of a transition-state structure for CO dissociation catalyzed by a dicopper cluster, but with a very small (3,3) active space. More significantly, Ma et al. [32] accomplished CASSCF geometry optimizations for all 12 excited states of a nickel bipyridine complex using (10,10) and (10,11) active spaces.
More commonly, MR methods may be used to compute the energetics of various points along reaction coordinates that are defined by DFT calculations. [33] That is, MR energies are computed as single-point calculations on structures optimized at the DFT level (and MR spin ladders may simultaneously be used to inform the choice of DFT functional used to generate the reaction coordinate). This approach has been applied to study catalysis by vanadium oxide [34] and bimetallic cobalt oxide clusters, [35] iron fluorides, [36] polyoxometalates, [37] transition metal organometallic complexes (e. g. MnN4py, [28c] [(Me 2 Py 2 TACN)FeN] 2 + , [38] porphyrins, [26a,33c,e,39] iron and ruthenium imido [40] and aquo complexes, [41] ruthenium diolefin diazadienes, [42] iridium-based systems, [43] and complexes featuring short metal-metal bonds such as Mn 2 (C 6 H 6 ) 2 [44] and N,N,N-tri(2-(2pyridylamino)ethyl)amine complexes [45] ), metal-organic frameworks (MIL-100(Fe), [10b] Fe 0.1 Mg 1.9 (dobdc) 2 , [10a] ZIF-8(Fe) [31] ), and iron-based zeolites (Fe-BEA*, [46] Fe-CHA, [30] Fe-FER, [46b] Fe-ZSM5). [46b] For any given study, a significant complication is the choice of orbitals to include in the active space. This is often far from straightforward as experience and chemical intuition is required to appreciate the full range of orbitals whose partial occupation may be expected to lead to highly weighted configuration state functions in the final MR wave function; extensive testing of alternative choices is typically required to be certain of accuracy. However, algorithms for automated active space selections have begun to emerge and offer considerable promise with respecting to ameliorating this complication (Section 3.2).

Reducing the Cost of MR Calculations
To reduce the cost of MR calculations applied to reaction coordinates for catalysis one can pursue any or all of: (1) using an approximation to the CASSCF wave function, (2) simplifying the post-CASSCF step to compute dynamic correlation energy, or (3) simplifying structure optimization.

Lowering the Computational Cost of Accounting for Static Correlation
The cost of CASSCF calculations scales exponentially with the active space size because that is how the number of configurations in the CI expansion grows. Many of these configurations contribute only minimally to the wave function, or, as Ruedenberg said, "they are deadwood". [47] In order to reduce the number of configurations, one can use the so-called restricted active space self-consistent field (RASSCF) method. [48] In a RASSCF calculation, the active space is divided into three subspaces, RAS1, RAS2, and RAS3. The RAS1 subspace contains doubly occupied orbitals from which excitations are allowed only up to a user-defined level (e. g., singles, doubles, etc., but less than would be the case for a full CAS). In RAS2, a FCI is performed. Finally, the RAS3 subspace contains virtual orbitals to which user-defined excitations are allowed (e. g., singles, doubles, etc., but again, less than would be allowed in a CAS calculation). [1,48] The concept behind RASSCF is taken still further with the generalized active space SCF (GASSCF), [49] the occupation-restricted multiple active space (ORMAS), [50] and the SplitGAS methods. [51] These are all different flavors of truncated CI methods, as are the separated pair (SP) approximation [52] and the generalized valence bond (GVB) [53] method. Utilization of these approaches is more limited than CASSCF because they tend to be less general, and the quality of their results tends to depend on expert choices of active subspaces and levels of excitations. In addition, as they are truncated-CI methods, they are not size extensive.
A separate large cost in wave-function calculations compared to DFT in general is associated with computation of the four-index two-electron integrals. This can be reduced by using density fitting (DF) or resolution-of-theidentity (RI) approximations, [54] in which product densities are expanded in auxiliary basis sets. [55] For example, Cholesky decomposition decreases the cost of CASSCF calculations by one-to-two orders of magnitude, making it possible to study systems of medium size with basis sets including diffuse functions, for instance. [55] Nevertheless, the active space size can remain a bottleneck.
One promising approach is the density matrix renormalization group (DMRG) algorithm, [56] which is an approximation to an FCI (or CASSCF) solver. [56e] A review of the DMRG algorithm is reported in Refs. [57]. In the DMRG approach, each orbital is treated as a ''site'' in a onedimensional lattice. In the guess, the orbitals are ordered following various criteria, e. g. by their Hartree-Fock energy.
[56e] Once the orbitals are ordered, they are divided in collections or "blocks", based on one or more properties. The possible states for the many-body blocks are sorted in an optimal way by DMRG as the eigenstates of a manyparticle density matrix. During each step, two blocks are considered. At each iteration, the orbitals from one block are moved (or swept) to the other block one at a time, always from the right to the left block, until the right block has only one orbital, and then reversing the direction. [56e] The states defining each block are chosen after each sweep. [56e] The cost of DMRG depends on the number of states and on the number of electrons considered. DMRG applied to CASSCF has a polynomial scaling with respect to the size of the active space, as opposite to the exponential scaling of standard CASSCF. As a result, DMRG-CASSCF allows practical computation of numerically well-converged solutions for active spaces up to four times larger than standard CASSCF. [58] DMRG has been used for MR calculations of the electronic structures of various enzyme cofactors implicated in key biological processes. [59] Similarly, DMRG has been used in conjunction with active spaces of (30,32) and (54,36) by Sharma et al. [59a] to study [2FeÀ 2S] and [4FeÀ 4S] clusters, respectively. The larger dimensions of active spaces that can be handled by DMRG reduce the responsibility of the user with respect to selecting a balanced collection of orbitals. However, DMRG itself is not yet a black-box method as selection of userdefined parameters involving shapes and orders of orbitals remain important, these being associated mainly with the one-dimensional nature of the DMRG wavefunction ansatz. [24,56f,58] An alternative strategy to DMRG to avoid the exponential scaling of the CASSCF cost on the active space size is to combine [60] MC-SCF methods with Full Configuration Interaction Quantum Monte Carlo (FCIQMC). [61] FCIQMC uses a probabilistic approach to solve the CI problem while accumulating one-and two-body reduced density matrices. [60b] The method combining CASSCF with FCIQMC is indicated in the literature as FCIQMC-CASSCF or stochastic-CASSCF. [60a] Application of stochastic-CASSCF to iron-sulfur dimer models has included active spaces as large as (22,26), [62] which exceeds the current limit of (24,24) achieved by CASSCF, [29] without the need of massive parallel implementations.
As another alternative, density matrix embedding theory (DMET), [63] is an efficient wave-function-based embedding approach to treat strongly correlated systems. [64] In DMET, a system is divided into fragments where each fragment is treated using a different quantum mechanical level. This means that a higher accuracy can be used for some specific fragment (e. g., in a catalytic system, the metal center and its first two coordination shells), while a significant decrease in computational resources can be achieved by treating the rest of the system with a much lower level theory. In the simplest case, the structure is divided into two parts, the "impurity" and the "environment", where the impurity is the most important fragment, i. e., the one to be treated using the highest level of theory. Phan et al. [64] have described CASSCF as a correlated impurity solver for DMET (CAS-DMET), and CAS-DMET has been successfully used to describe single bond breaking in the H 10 ring model system as well as N=N double bond breaking in azomethane and pentyldiazene. For the successful applica-tion of CAS-DMET, a proper choice of the impurity and of the environment must be accomplished, and that can be less straightforward in arbitrary systems compared to those in Ref. [64]. Applications to realistic chemical systems call for improvements in the convergence rate of both the chemical potential and correlation potential, which are key components linking the DMET subsystems. [64]

Lowering the Computational Cost of Computing External Correlation
A perennial challenge in quantum chemistry is the design of methods that compute the external correlation energy efficiently. [24,65] Second-order restricted active space perturbation theory (RASPT2) [66] is more affordable than CASPT2, but also less general. DMRG-CASPT2 and DMRG-NEVPT2 have also been reported, although they are not yet practical for large systems as the computational cost remains mostly determined by the active space size as for CASPT2 and RASPT2. [56f,57,58,67] This cost can be significantly reduced by avoiding or reducing the cost of the computation of the four-particle reduced density matrix (4-RDM). [57,67f] This can be accomplished through different strategies, [33d,56f,57,58,67b-h] including the employment of pseudo-canonical orbitals, [67f] or deriving the 4-RDM from 3-RDM calculations. [57,68] Phung and Pierloot showed that the latter approach is particularly promising in terms of both accuracy of the results and computational cost in a study on iron-oxo porphyrin. [69] Overall, the cost of recovering external correlation remains prohibitive in many cases. Multiconfiguration pair-density functional theory (MC-PDFT) [65,70] has considerable promise, since it is less expensive than CASPT2. [70b,71] Importantly, the accuracy in the calculation of reaction energetics involving small and medium-size organic molecules is similar to CASPT2. [72] Considering spin-state energetics, MC-PDFT has been used in combination with a DMRG wave function, showing very promising results for singlet-triplet gaps in polyacenes and polyacetylenes. [24]

Lowering the Computational Cost of a Structure Optimization
MR methods are often employed to perform single-point energy calculations on KS-DFT optimized geometries. In many cases, KS-DFT predicts reliable ground-state equilibrium structures. Non-equilibrium geometries and excitedstate structures can be more challenging, as they may need an accurate accounting of spin localization(s), which are not necessarily well predicted by single-reference methods. [73] By contrast, while optimization with MR methods can avoid such deficiencies, a new challenge is posed in so far as the important frontier orbitals may change along a reaction coordinate, rendering the active space choice decidedly non-trivial. In fact, an important challenge in the use of MR methods to model chemical reactions in general (and catalytic reactions in particular) is to define a common active space for the reactant(s), transition state(s), and product(s). In addition, the computational cost associated with gradient calculations is usually an order of magnitude greater than for a single point calculation, creating a resource burden. Unsurprisingly, most structures optimized at MR levels tend to contain only light elements, with a few exceptions. [74] In recent years, however, a significant improvement in the efficiency of minimization algorithms has been attained, in part because of the implementation of analytical gradients for many MR methods (e. g. CASSCF, [74b] CASPT2, [75] DMRG-CASSCF, [67a] MC-PDFT [76] ). The most efficient algorithms combine analytical gradients with density fitting techniques for twoelectron integrals. [74b,76c] However, the calculation of second derivatives (to compute thermochemical quantities or vibrational frequencies of molecular probes in catalysis) remains a challenge.
As already alluded to above, one simplifying approach that can be used irrespective of quantum methodology involves the application of less complete levels of theory to regions well separated from the active site. For this purpose, QM/MM methods can be quite powerful. [77] QM/MM methods in which the QM level is a MR model have been used to study the nature of FeÀ O 2 bonding in oxymyoglobin [78] and in non-heme iron enzymes. [79] An embedded MR/DFT approach has been recently adopted by Carter and co-workers to study adsorbates on surfaces [80] and heterogeneous catalysis, [33a] in particular plasmoninduced catalysis in metal nanoparticles [81] including methane dry reforming [82] and plasmon-induced ammonia decomposition on copper nanoparticles. [83]

Towards the Automatic Selection of the Active Space
What constitutes a good active space? In a nutshell, a good active space generates a qualitatively correct wave function that can be used in subsequent post-SCF calculations to predict quantitatively accurate energetics. Having at hand infinite computational resources, and time, one would test more and more active spaces of increasing size and composed of different sets of orbitals. One could also compare the different wave functions and the energetics obtained by PT2 or PDFT treatments on top of the MR wave functions. Of course, such a procedure is not practical with finite resources, except for very small systems. A more efficient approach typically relies on user experience, with criteria for confirming the quality of the optimized active orbitals often strongly system dependent. While this makes it difficult to generalize rules for active space selection across all chemical systems, one can still define some transferable rules for systems having certain chemical similarities. For example, for systems based on transition metals, the 3d orbitals should be active, especially those that are only partially occupied because of the large static correlation among the electrons in these orbitals. It is generally advisable to include another shell of correlating d orbitals, 3d' (unless the metal has very few d electrons to begin). Some relevant methods designed to standardize the (non-automatic) active orbitals selection for multiconfigurational calculations can be found in Refs. [84], and a discussion and comparison among them can be found in Refs. [85]. Obviously, it would be helpful to minimize human intervention and develop automated schemes that make these decisions.
One algorithm for a transition-metal system could launch the initial guess of the SCF cycle with an active space containing the 3d metal orbitals and force them to remain in the active space for some number of iterations, considering also the energetic ranking expected based on crystal field theory or natural bond occupations. [86] This would involve an alternation of CASSCF and configuration interaction (CI) steps. The constraint on the active space would be removed after a certain number of iterations, permitting energy-based updates. Several schemes for automated active space selection have been reported in the literature, [84a,c,87] based on imposed [86,88] and ranked orbital [85b] approaches. While these schemes show promise, they lack generality, i. e., they work only for specific sets of systems, with limited transferable prescriptions. [89] Stein and Reiher [90] have taken another approach to this challenge. Instead of imposing a set of rules, their AutoCAS exploits the use of DMRG for the approximate CI solver in order to use very large active spaces. Smaller active spaces are then selected based on orbital entropy information and used to perform calculations at an increased level of accuracy. In this way, the selection of the active space is automated. Nevertheless, this approach is only seemingly black box because, as noted in Section 3.1.1., the DMRG calculations require substantial practitioner oversight. Moreover, they come with a significant computational cost. [90d] It is intriguing to consider the application of artificial intelligence to this problem. Machine-learning algorithms should be able to learn the same chemical principles that permit a human to choose a proper active space depending on the system. Acquisition of "chemical knowledge" by machine learning algorithms has been demonstrated by Jablonka et al. in a study designed to assign automatically the correct oxidation states of atoms. [91] More recently, Jeong et al. [92] have reported a proof-of-concept study using machine learning for automated active space selection. Their results, although limited to small chemical spaces, demonstrate the feasibility of an automated procedure for active space selection that can in principle be enlarged by increasing the dimensions and variability of the training set. Golub et al. [93] in another recent study developed an alternative machine learning algorithm for the automatic selection of active spaces, verifying its transferability to systems outside the training set. In particular they successfully tested it for Fe(II)-porphyrin, a model system for the active site of several enzymes.

Outlook
Many strategies are now actively pursued to reduce the cost of MR computational methods and to make them more black-box-like in character. While selection of the active space has historically been considered a step requiring human intervention, it is encouraging that promising automated procedures have started to appear, particularly those exploiting powerful machine learning techniques. [92][93] This should lead to MR methods becoming more routine, particularly for single point calculations. Future studies should target enlarging the set of systems to which these approaches can be systematically applied. A further development of these methods that would be particularly rewarding would be the selection of active spaces that can be used all along a reaction pathway.
With respect to routine application of MR methods to catalysis, some problems will likely remain intractable for a longer time. For example, MR geometry optimizations are likely to require human supervision more than single-point calculations because of how the active space can change along an optimization (or reaction) coordinate. Of course, as machine learning becomes more robust, improved recognition of such variation will be possible. To date, optimizations of copper, [74a] NiFe, [94] iron, [74c] and ruthenium complexes [74b] have been reported, but the resources required are too high to make them routine. As noted above, multiscalar methods that combine KS-DFT geometry optimizations with single-point DMRG-PDFT steps remain an interim option, and these are certainly amenable to automation.
With respect to avoiding the computation of wave functions altogether, machine learning methods have already been used to avoid DFT SCF calculations for very large systems [95] either by directly training on energies [17,96] or by developing interatomic potentials to replace the DFT calculation. [1,97] By analogy, machine learning could replace MR calculations as well, but appropriate, large MR data sets that would be required for training are not presently available.
Finally, advances in quantum computing offer the potential to speed dramatically MR calculations. Recently Tilly et al. [98] reported the optimization of CO using CASSCF on a 4-qubit system. Although current "noisy intermediate-scale quantum" technologies are plagued by random error, [99] Tilly et al. achieved convergence of the active space orbitals through a light-touch mitigation strategy that compensated for gate errors. While a single CO molecule is hardly a catalyst surface, continued progress on the development of algorithms and theory suitable to the quantum computing technologies of today and tomorrow has the potential for transformational impact.