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

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.

Algorithm for analytic nuclear energy gradients of state averaged DMRG-CASSCF theory with newly derived coupled-perturbed equations.

We present an algorithm for evaluating analytic nuclear energy gradients of the state-averaged density matrix renormalization group complete-active-space self-consistent field (SA-DMRG-CASSCF) theory based on the newly derived coupled-perturbed (CP) DMRG-CASSCF equations. The Lagrangian for the conventional SA-CASSCF analytic gradient theory is extended to the SA-DMRG-CASSCF variant that can fully consider a whole set of constraints on the parameters of multi-root canonical matrix product states formed at all the DMRG block configurations. An efficient algorithm to solve the CP-DMRG-CASSCF equations for determining the multipliers was developed. The complexity of the resultant analytic gradient algorithm is overall the same as that of the unperturbed SA-DMRG-CASSCF algorithm. In addition, a reduced-scaling approach was developed to directly compute the SA reduced density matrices (SA-RDMs) and their perturbed ones without calculating separate state-specific RDMs. As part of our implementation scheme, we neglect the term associated with the constraint on the active orbitals in terms of the active-active rotation in the Lagrangian. Thus, errors from the true analytic gradients may be caused in this scheme. The proposed gradient algorithm was tested with the spin-adapted implementation by checking how accurately the computed analytic energy gradients reproduce numerical gradients of the SA-DMRG-CASSCF energies using a common number of renormalized bases. The illustrative applications show that the errors are sufficiently small when using a typical number of the renormalized bases, which is required to attain adequate accuracy in DMRG's total energies.

Multireference Electron Correlation Methods: Journeys along Potential Energy Surfaces.

Multireference electron correlation methods describe static and dynamical electron correlation in a balanced way and, therefore, can yield accurate and predictive results even when single-reference methods or multiconfigurational self-consistent field theory fails. One of their most prominent applications in quantum chemistry is the exploration of potential energy surfaces. This includes the optimization of molecular geometries, such as equilibrium geometries and conical intersections and on-the-fly photodynamics simulations, both of which depend heavily on the ability of the method to properly explore the potential energy surface. Because such applications require nuclear gradients and derivative couplings, the availability of analytical nuclear gradients greatly enhances the scope of quantum chemical methods. This review focuses on the developments and advances made in the past two decades. A detailed account of the analytical nuclear gradient and derivative coupling theories is presented. Emphasis is given to the software infrastructure that allows one to make use of these methods. Notable applications of multireference electron correlation methods to chemistry, including geometry optimizations and on-the-fly dynamics, are summarized at the end followed by a discussion of future prospects.

Computational Spectroscopy of the Cr-Cr Bond in Coordination Complexes.

We report the accurate computational vibrational analysis of the Cr-Cr bond in dichromium complexes using second-order multireference complete active space methods (CASPT2), allowing direct comparison with experimental spectroscopic data both to facilitate interpreting the low-energy region of the spectra and to provide insights into the nature of the bonds themselves. Recent technological development by the authors has realized such computation for the first time. Accurate simulation of the vibrational structure of these compounds has been hampered by their notorious multiconfigurational electronic structure that yields bond distances that do not correlate with bond order. Some measured Cr-Cr vibrational stretching modes, ν(Cr2), have suggested weaker bonding, even for so-called ultrashort Cr-Cr bonds, while others are in line with the bond distance. Here, we optimize geometries and compute ν(Cr2) with CASPT2 for three well-characterized complexes, Cr2(O2CCH3)4(H2O)2, Cr2(mhp)4, and Cr2(dmp)4. We obtain CASPT2 harmonic ν(Cr2) modes in good agreement with experiment at 282 cm-1 for Cr2(mhp)4 and 353 cm-1 for Cr2(dmp)4, compute 50Cr and 54Cr isotope shifts, and demonstrate that the use of the so-called IPEA shift leads to improved Cr-Cr distances. Additionally, normal mode sampling was used to estimate anharmonicity along ν(Cr2), leading to an anharmonic mode of 272 cm-1 for Cr2(mhp)4 and 333 cm-1 for Cr2(dmp)4.

Efficient and stochastic multireference perturbation theory for large active spaces within a full configuration interaction quantum Monte Carlo framework.

Full Configuration Interaction Quantum Monte Carlo (FCIQMC) has been effectively applied to very large configuration interaction (CI) problems and was recently adapted for use as an active space solver and combined with orbital optimization. In this work, we detail an approach within FCIQMC to allow for efficient sampling of fully internally contracted multireference perturbation theories within the same stochastic framework. Schemes are described to allow for the close control over the resolution of stochastic sampling of the effective higher-body intermediates within the active space. It is found that while complete active space second-order perturbation theory seems less amenable to a stochastic reformulation, strongly contracted N-Electron Valence second-order Perturbation Theory (NEVPT2) is far more stable, requiring a similar number of walkers to converge the sc-NEVPT2 expectation values as to converge the underlying CI problem. We demonstrate the application of the stochastic approach to the computation of sc-NEVPT2 within a (24, 24) active space in a biologically relevant system and show that small numbers of walkers are sufficient for a faithful sampling of the sc-NEVPT2 energy to chemical accuracy, despite the active space already exceeding the limits of practicality for traditional approaches. This raises prospects of an efficient stochastic solver for multireference chemical problems requiring large active spaces, with an accurate treatment of external orbitals.

Visualizing Complex-Valued Molecular Orbitals.

We report an implementation of a program for visualizing complex-valued molecular orbitals. The orbital phase information is encoded on each of the vertices of triangle meshes using the standard color wheel. Using this program, we visualized the molecular orbitals for systems with spin-orbit couplings, external magnetic fields, and complex absorbing potentials. Our work has not only created visually attractive pictures but also clearly demonstrated that the phases of the complex-valued molecular orbitals carry rich chemical and physical information on the system, which has often been unnoticed or overlooked.

Large-scale relativistic complete active space self-consistent field with robust convergence.

We report an efficient algorithm using density fitting for the relativistic complete active space self-consistent field (CASSCF) method, which is significantly more stable than the algorithm previously reported by one of the authors [J. E. Bates and T. Shiozaki, J. Chem. Phys. 142, 044112 (2015)]. Our algorithm is based on the second-order orbital update scheme with an iterative augmented Hessian procedure, in which the density-fitted orbital Hessian is directly contracted to the trial vectors. Using this scheme, each microiteration is made less time consuming than one Dirac-Hartree-Fock iteration, and macroiterations converge quadratically. In addition, we show that the CASSCF calculations with the Gaunt and full Breit interactions can be efficiently performed by means of approximate orbital Hessians computed with the Dirac-Coulomb Hamiltonian. It is demonstrated that our algorithm can also be applied to systems under an external magnetic field, for which all of the molecular integrals are computed using gauge-including atomic orbitals.

On the difference between variational and unitary coupled cluster theories.

There have been assertions in the literature that the variational and unitary forms of coupled cluster theory lead to the same energy functional. Numerical evidence from previous authors was inconsistent with this claim, yet the small energy differences found between the two methods and the relatively large number of variational parameters precluded an unequivocal conclusion. Using the Lipkin Hamiltonian, we here present conclusive numerical evidence that the two theories yield different energies. The ambiguities arising from the size of the cluster parameter space are absent in the Lipkin model, particularly when truncating to double excitations. We show that in the symmetry adapted basis under strong correlation, the differences between the variational and unitary models are large, whereas they yield quite similar energies in the weakly correlated regime previously explored. We also provide a qualitative argument rationalizing why these two models cannot be the same. Additionally, we study a generalized non-unitary and non-hermitian variant that contains excitation, de-excitation, and mixed operators with different amplitudes and show that it works best when compared to the traditional, variational, unitary, and extended forms of coupled cluster doubles theories.

Electronic structure of SmO and SmO- via slow photoelectron velocity-map imaging spectroscopy and spin-orbit CASPT2 calculations.

The chemi-ionization reaction of atomic samarium, Sm + O → SmO+ + e-, has been investigated by the Air Force Research Laboratory as a means to modify local electron density in the ionosphere for reduction of scintillation of high-frequency radio waves. Neutral SmO is a likely unwanted byproduct. The spectroscopy of SmO is of great interest to aid in interpretation of optical emission spectra recorded following atmospheric releases of Sm as part of the Metal Oxide Space Cloud (MOSC) observations. Here, we report a joint experimental and theoretical study of SmO using slow photoelectron velocity-map imaging spectroscopy of cryogenically cooled SmO- anions (cryo-SEVI) and high-level spin-orbit complete active space calculations with corrections from second order perturbation theory (CASPT2). With cryo-SEVI, we measure the electron affinity of SmO to be 1.0581(11) eV and report electronic and vibrational structure of low-lying electronic states of SmO in good agreement with theory and prior experimental work. We also obtain spectra of higher-lying excited states of SmO for direct comparison to the MOSC results.

Fully relativistic self-consistent field under a magnetic field.

We present a gauge-invariant implementation of the four-component Dirac-Hartree-Fock method for simulating the electronic structure of heavy element complexes in magnetic fields. The additional cost associated with the magnetic field is shown to be only 10-13% of that at zero field. The Dirac-Hartree-Fock wave function is constructed from gauge-including atomic orbitals. The so-called restricted magnetic balance is used to generate 2-spinor basis functions for the small component. The molecular integrals for the Coulomb and Gaunt interactions are computed using density fitting. Our efficient, parallel implementation allows for simulating the electronic structure of molecules containing more than 100 atoms with a few heavy elements under magnetic fields.

Communication: automatic code generation enables nuclear gradient computations for fully internally contracted multireference theory.

Analytical nuclear gradients for fully internally contracted complete active space second-order perturbation theory (CASPT2) are reported. This implementation has been realized by an automated code generator that can handle spin-free formulas for the CASPT2 energy and its derivatives with respect to variations of molecular orbitals and reference coefficients. The underlying complete active space self-consistent field and the so-called Z-vector equations are solved using density fitting. The implementation has been applied to the vertical and adiabatic ionization potentials of the porphin molecule to illustrate its capability.

Fully relativistic complete active space self-consistent field for large molecules: quasi-second-order minimax optimization.

We develop an efficient algorithm for four-component complete active space self-consistent field (CASSCF) methods on the basis of the Dirac equation that takes into account spin-orbit and other relativistic effects self-consistently. Orbitals are optimized using a trust-region quasi-Newton method with Hessian updates so that energies are minimized with respect to rotations among electronic orbitals and maximized with respect to rotations between electronic and positronic orbitals. Utilizing density fitting and parallel computation, we demonstrate that Dirac-Coulomb CASSCF calculations can be routinely performed on systems with 100 atoms and a few heavy-elements. The convergence behavior and wall times for octachloridodirhenate(III) and a tungsten methylidene complex are presented. In addition, the excitation energies of octachloridodirhenate(III) are reported using a state-averaged variant.

Communication: Active space decomposition with multiple sites: density matrix renormalization group algorithm.

We extend the active space decomposition method, recently developed by us, to more than two active sites using the density matrix renormalization group algorithm. The fragment wave functions are described by complete or restricted active-space wave functions. Numerical results are shown on a benzene pentamer and a perylene diimide trimer. It is found that the truncation errors in our method decrease almost exponentially with respect to the number of renormalization states M, allowing for numerically exact calculations (to a few µE(h) or less) with M = 128 in both cases. This rapid convergence is because the renormalization steps are used only for the interfragment electron correlation.

Fast Emulation of Fermionic Circuits with Matrix Product States.

We describe a matrix product state (MPS) extension for the Fermionic Quantum Emulator (FQE) software library. We discuss the theory behind symmetry-adapted MPSs for approximating many-body wave functions of spin-1/2 Fermions, and we present an open-source, MPS-enabled implementation of the FQE interface (MPS-FQE). The software uses the open-source pyblock3 and block2 libraries for most elementary tensor operations, and it can largely be used as a drop-in replacement for FQE that allows for more efficient but approximate emulation of larger Fermionic circuits. Finally, we show several applications relevant to both near-term and fault-tolerant quantum algorithms where approximate emulation of larger systems is expected to be useful: characterization of state preparation strategies for quantum phase estimation, the testing of different variational quantum eigensolver ansätze, the numerical evaluation of Trotter errors, and the simulation of general quantum dynamics problems. In all these examples, approximate emulation with MPS-FQE allows us to treat systems that are significantly larger than those accessible with a full statevector emulator.

Pyrazine excited states revisited using the extended multi-state complete active space second-order perturbation method.

We demonstrate that the recently developed extended multi-state complete active space second-order perturbation theory (XMS-CASPT2) [Shiozaki et al., J. Chem. Phys., 2011, 135, 081106] provides qualitatively correct potential energy surfaces for low-lying excited singlet states of pyrazine, while the potential energy surfaces of the standard MS-CASPT2 methods are ill-behaved near the crossing point of two reference potential energy surfaces. The XMS-CASPT2 method is based on the extended multi-configuration quasi-degenerate perturbation theory proposed earlier by Granovsky [J. Chem. Phys., 2011, 134, 214113]. We show that the conical intersection at the XMS-CASPT2 level can be described without artifacts if the entire method is invariant with respect to any unitary rotations of the reference functions. The photoabsorption spectra of the 1(1)B(3u) and 1(1)B(2u) states of pyrazine are simulated, based on a vibronic-coupling model Hamiltonian. The XMS-CASPT2 spectrum of the 1(1)B(3u) band is found to be comparable to the one computed by a more expensive multireference configuration interaction (MRCI) method, while the XMS-CASPT2 simulation of the 1(1)B(2u) band is slightly inferior to the MRCI one.

Communication: An efficient algorithm for evaluating the Breit and spin-spin coupling integrals.

We present an efficient algorithm for evaluating a class of two-electron integrals of the form r12⊗r12/r12(n) over one-electron Gaussian basis functions. The full Breit interaction in four-component relativistic theories beyond the Gaunt term is such an operator with n = 3. Another example is the direct spin-spin coupling term in the quasi-relativistic Breit-Pauli Hamiltonian (n = 5). These integrals have been conventionally evaluated by expensive derivative techniques. Our algorithm is based on tailored Gaussian quadrature, similar to the Rys quadrature for electron repulsion integrals (ERIs), and can utilize the so-called horizontal recurrence relation to reduce the computational cost. The CPU time for computing all six Cartesian components of the Breit or spin-spin coupling integrals is found to be only 3 to 4 times that of the ERI evaluation.

Large-scale Dirac-Fock-Breit method using density fitting and 2-spinor basis functions.

We present an efficient theory and algorithm for computing four-component relativistic Dirac-Fock wave functions using the Coulomb, Gaunt, and full Breit interactions. Our implementation is based on density fitting, and is routinely applicable to systems with 100 atoms and a few heavy elements. The small components are expanded using 2-spinor basis functions. We show that the factorization of 3-index half-transformed integrals before building Coulomb and exchange matrices is essential for efficient evaluation of the Fock matrix. With the Coulomb interaction, the computational cost for evaluating the Fock operator has been found to be only 70-90 times that in the non-relativistic density-fitted Hartree-Fock method. The prefactors have been 170 and 350-450 for the Gaunt and Breit interactions, respectively. The largest molecule to which we have applied the Dirac-Fock-Coulomb method is an Ac(III) motexafin complex (130 atoms, 556 electrons, 1289 basis functions), for which one self-consistent iteration takes around 1100 s using 1024 CPU cores. In addition, we have found that, while the standard fitting basis sets are accurate for Dirac-Fock-Coulomb calculations, their accuracy is very poor for Dirac-Fock-Gaunt and Breit calculations. We report a prototype of accurate fitting basis sets for these cases.

Analytical energy gradients for second-order multireference perturbation theory using density fitting.

We present algorithms for computing analytical energy gradients for multi-configuration self-consistent field methods and partially internally contracted complete active space second-order perturbation theory (CASPT2) using density fitting (DF). Our implementation is applicable to both single-state and multi-state CASPT2 analytical gradients. The accuracy of the new methods is demonstrated for structures and excitation energies of valence and Rydberg states of pyrrole, as well as for structures and adiabatic singlet-triplet energy splittings for the hydro-, the O,O(')-formato-, and the N,N(')-diiminato-copper-dioxygen complexes. It is shown that the effects of density fitting on optimized structures and relative energies are negligible. For cases in which the total cost is dominated by the integral evaluations and transformations, the DF-CASPT2 gradient calculations are found to be faster than the corresponding conventional calculations by typically a factor of three to five using triple-ζ basis sets, and by about a factor of ten using quadruple-ζ basis sets.

Communication: Active-space decomposition for molecular dimers.

We have developed an active-space decomposition strategy for molecular dimers that allows for the efficient computation of the dimer's complete-active-space wavefunction while only constructing the monomers' active-space wavefunctions. Dimer states are formed from linear combinations of direct products of localized orthogonal monomer states and Hamiltonian matrix elements are computed directly without explicitly constructing the product space. This decomposition is potentially exact in the limit where a full set of monomer states is included. The adiabatic states are then found by diagonalizing the dimer Hamiltonian matrix. We demonstrate the convergence of our method to a complete-active-space calculation of the full dimer with two test cases: the benzene and naphthalene dimers.