Improving the Efficiency of the Multireference Driven Similarity Renormalization Group via Sequential Transformation, Density Fitting, and the Noninteracting Virtual Orbital Approximation

Forte: A suite of advanced multireference quantum chemistry methods

Forte is an open-source library specialized in multireference electronic structure theories for molecular systems and the rapid prototyping of new methods. This paper gives an overview of the capabilities of Forte, its software architecture, and examples of applications enabled by the methods it implements.

Frozen Natural Orbitals for the State-Averaged Driven Similarity Renormalization Group

We present a reduced-cost implementation of the state-averaged driven similarity renormalization group (SA-DSRG) based on the frozen natural orbital (FNO) approach. The natural orbitals (NOs) are obtained by diagonalizing the one-body reduced density matrix from SA-DSRG second-order perturbation theory (SA-DSRG-PT2). We consider three criteria to truncate the virtual NOs for the subsequent electron correlation treatment beyond SA-DSRG-PT2. An additive second-order correction is applied to the SA-DSRG Hamiltonian to reintroduce correlation effects from the discarded orbitals. The FNO SA-DSRG method is benchmarked on 35 small organic molecules in the QUEST database. When keeping 98-99% of the cumulative occupation numbers, the mean absolute error in the vertical transition energies due to FNO is less than 0.01 eV. Using the same FNO threshold, we observe a speedup of 9 times compared to the conventional SA-DSRG implementation for nickel carbonyl with a quadruple-ζ basis set. The FNO approach enables nonperturbative SA-DSRG computations on chloroiron corrole [FeCl(C19H11N4)] with more than 1000 basis functions, surpassing the current limit of a conventional implementation.

Multicomponent Cholesky Decomposition: Application to Nuclear-Electronic Orbital Theory

The Cholesky decomposition technique is commonly used to reduce the memory requirement for storing two-particle repulsion integrals in quantum chemistry calculations that use atomic orbital bases. However, when quantum methods use multicomponent bases, such as nuclear-electronic orbitals, additional challenges are introduced due to asymmetric two-particle integrals. This work proposes several multicomponent Cholesky decomposition methods for calculations using nuclear-electronic orbital density functional theory. To analyze the errors in different Cholesky decomposition components, benchmark calculations using water clusters are carried out. The largest benchmark calculation is a water cluster (H2O)27 where all 54 protons are treated quantum mechanically. This study provides energetic and complexity analyses to demonstrate the accuracy and performance of the proposed multicomponent Cholesky decomposition method.

Relativistic resolution-of-the-identity with Cholesky integral decomposition

In this study, we present an efficient integral decomposition approach called the restricted-kinetic-balance resolution-of-the-identity (RKB-RI) algorithm, which utilizes a tunable RI method based on the Cholesky integral decomposition for in-core relativistic quantum chemistry calculations. The RKB-RI algorithm incorporates the restricted-kinetic-balance condition and offers a versatile framework for accurate computations. Notably, the Cholesky integral decomposition is employed not only to approximate symmetric large-component electron repulsion integrals but also those involving small-component basis functions. In addition to comprehensive error analysis, we investigate crucial conditions, such as the kinetic balance condition and variational stability, which underlie the applicability of Dirac relativistic electronic structure theory. We compare the computational cost of the RKB-RI approach with the full in-core method to assess its efficiency. To evaluate the accuracy and reliability of the RKB-RI method proposed in this work, we employ actinyl oxides as benchmark systems, leveraging their properties for validation purposes. This investigation provides valuable insights into the capabilities and performance of the RKB-RI algorithm and establishes its potential as a powerful tool in the field of relativistic quantum chemistry.

Unraveling the Spin-State Energetics of FeN4 Complexes with Ab Initio Methods

A systematic analysis was conducted to explore the spin-state energetics of a series of 19 FeN4 complexes. The performance of a large number of multireference methods was assessed, highlighting the significant challenges associated with accurately describing the spin-state energetics of FeN4 complexes. Most multireference methods were found to be susceptible to errors originating from the reference CASSCF wavefunction, leading to an overstabilization of high-spin states. Nonetheless, a few multireference methods, namely, CASPT2/CC, DSRG-MRPT3, and LDSRG(2), demonstrated promising performance compared to the benchmark CCSD(T) method. Furthermore, our study revealed that FeN4 complexes having a quintet ground state are exceedingly rare. Accordingly, only one specific model (Fe(L2)) and one synthesized complex (Fe(OTBP)) have the quintet ground state among the studied complexes. This scarcity of quintet FeN4 complexes highlights the unique nature of these systems and raises intriguing questions regarding the factors influencing spin states, such as the size of the macrocycle cavity, the introduction of substituents, or the induction of out-of-plane deformation.

A Similarity Renormalization Group Approach to Green's Function Methods

The family of Green's function methods based on the GW approximation has gained popularity in the electronic structure theory thanks to its accuracy in weakly correlated systems combined with its cost-effectiveness. Despite this, self-consistent versions still pose challenges in terms of convergence. A recent study [Monino and Loos J. Chem. Phys. 2022, 156, 231101.] has linked these convergence issues to the intruder-state problem. In this work, a perturbative analysis of the similarity renormalization group (SRG) approach is performed on Green's function methods. The SRG formalism enables us to derive, from first-principles, the expression of a naturally static and Hermitian form of the self-energy that can be employed in quasiparticle self-consistent GW (qsGW) calculations. The resulting SRG-based regularized self-energy significantly accelerates the convergence of qsGW calculations, slightly improves the overall accuracy, and is straightforward to implement in existing code.

Assessment of State-Averaged Driven Similarity Renormalization Group on Vertical Excitation Energies: Optimal Flow Parameters and Applications to Nucleobases

We present a comprehensive excited-state benchmark for the state-averaged (SA) driven similarity renormalization group (DSRG) [Li, C.; Evangelista, F. A. J. Chem. Phys. 2018, 148, 124106]. Following the QUEST database [Véril, M.; Scemama, A.; Caffarel, M.; Lipparini, F.; Boggio-Pasqua, M.; Jacquemin, D.; Loos, P.-F. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2021, 11, e1517], 280 vertical transition energies of 35 medium-sized molecules are computed using the SA-DSRG derived second- and third-order perturbation theories (PT2/PT3) along with a nonperturbative approach [sq-LDSRG(2)]. Comparing to the theoretical best estimates, the optimal flow parameter is found to be 0.35 and 2.0 Eh-2 for SA-DSRG-PT2 and SA-DSRG-PT3, respectively. For SA-sq-LDSRG(2), a flow parameter of 1.5 Eh-2 provides converged equations without compromising the accuracy. We then assess the accuracy of the SA-DSRG hierarchy using these parameters. The SA-DSRG-PT2 scheme outperforms the level-shifted CASPT2 by 0.10 eV in mean absolute error (MAE), yet this accuracy is slightly inferior than that of CASPT2 with the ionization-potential-electron-affinity shift. Both SA-DSRG-PT3 and SA-sq-LDSRG(2) yield a MAE of 0.10 eV, which is comparable to that of CASPT3 (0.09 eV). Finally, we compute vertical excitation energies of several low-lying singlet states of nucleobases. The SA-sq-LDSRG(2) approach provides highly accurate results for π → π* excitations, while n → π* transitions are better described by SA-DSRG-PT3.

Reference Energies for Cyclobutadiene: Automerization and Excited States

Cyclobutadiene is a well-known playground for theoretical chemists and is particularly suitable to test ground- and excited-state methods. Indeed, due to its high spatial symmetry, especially at the D_4h square geometry but also in the D_2h rectangular arrangement, the ground and excited states of cyclobutadiene exhibit multiconfigurational characters and single-reference methods, such as standard adiabatic time-dependent density-functional theory (TD-DFT) or standard equation-of-motion coupled cluster (EOM-CC), are notoriously known to struggle in such situations. In this work, using a large panel of methods and basis sets, we provide an extensive computational study of the automerization barrier (defined as the difference between the square and rectangular ground-state energies) and the vertical excitation energies at D_2h and D_4h equilibrium structures. In particular, selected configuration interaction (SCI), multireference perturbation theory (CASSCF, CASPT2, and NEVPT2), and coupled-cluster (CCSD, CC3, CCSDT, CC4, and CCSDTQ) calculations are performed. The spin-flip formalism, which is known to provide a qualitatively correct description of these diradical states, is also tested within TD-DFT (combined with numerous exchange-correlation functionals) and the algebraic diagrammatic construction [ADC(2)-s, ADC(2)-x, and ADC(3)]. A theoretical best estimate is defined for the automerization barrier and for each vertical transition energy.

Second-Order Active-Space Embedding Theory

Quantum embedding schemes are a promising way to extend multireference computations to large molecules with strong correlation effects localized on a small number of atoms. This work introduces a second-order active-space embedding theory [ASET(2)] which improves upon mean-field frozen embedding by treating fragment-environment interactions via an approximate canonical transformation. The canonical transformation employed in ASET(2) is formulated using the driven similarity renormalization group. The ASET(2) scheme is benchmarked on the N═N bond dissociation in pentyldiazene, the S₀ to S₁ excitation in 1-octene, and the interaction energy of the O₂-benzene complex. The ASET(2) explicit treatment of fragment-environment interactions beyond the mean-field level generally improves the accuracy of embedded computations, and it becomes necessary to achieve an accurate description of excitation energies of 1-octene and the singlet-triplet gap of the O₂-benzene complex.

Assessing the Accuracy of Tailored Coupled Cluster Methods Corrected by Electronic Wave Functions of Polynomial Cost

Tailored coupled cluster theory represents a computationally inexpensive way to describe static and dynamical electron correlation effects. In this work, we scrutinize the performance of various coupled cluster methods tailored by electronic wave functions of polynomial cost. Specifically, we focus on frozen-pair coupled cluster (fpCC) methods, which are tailored by pair-coupled cluster doubles (pCCD), and coupled cluster theory tailored by matrix product state wave functions optimized by the density matrix renormalization group (DMRG) algorithm. As test system, we selected a set of various small- and medium-sized molecules containing diatomics (N₂, F₂, C₂, CN⁺, CO, BN, BO⁺, and Cr₂) and molecules (ammonia, ethylene, cyclobutadiene, benzene, hydrogen chains, rings, and cuboids) for which the conventional single-reference coupled cluster singles and doubles (CCSD) method is not able to produce accurate results for spectroscopic constants, potential energy surfaces, and barrier heights. Most importantly, DMRG-tailored and pCCD-tailored approaches yield similar errors in spectroscopic constants and potential energy surfaces compared to accurate theoretical and/or experimental reference data. Although fpCC methods provide a reliable description for the dissociation pathway of molecules featuring single and quadruple bonds, they fail in the description of triple or hextuple bond-breaking processes or avoided crossing regions.

Theoretical Calculation of Core-Excited States along Dissociative Pathways beyond Second-Order Perturbation Theory

We extend the multireference driven similarity renormalization (MR-DSRG) method to compute core-excited states by combining it with a GASSCF treatment of orbital relaxation and static electron correlation effects. We consider MR-DSRG treatments of dynamical correlation truncated at the level of perturbation theory (DSRG-MRPT2/3) and iterative linearized approximations with one- and two-body operators [MR-LDSRG(2)] in combination with a spin-free exact-two-component (X2C) one-electron treatment of scalar relativistic effects. This approach is calibrated and tested on a series of 16 core-excited states of five closed- and open-shell diatomic molecules containing first-row elements (C, N, and O). All GASSCF-MR-DSRG theories show excellent agreement with experimental adiabatic transitions energies, with mean absolute errors ranging between 0.17 and 0.35 eV, even for the challenging partially doubly excited states of the N₂⁺ molecule. The vibrational structure of all these transitions, obtained from using a full potential energy scan, shows a mean absolute error as low as 25 meV for DSRG-MRPT2 and 12/13 meV for DSRG-MRPT3 and MR-LDSRG(2). We generally find that a treatment of dynamical correlation that goes beyond the second-order level in perturbation theory improves the accuracy of the potential energy surface, especially in the bond-dissociation region.

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.

High-level coupled-cluster energetics by merging moment expansions with selected configuration interaction

Inspired by our earlier semi-stochastic work aimed at converging high-level coupled-cluster (CC) energetics [J. E. Deustua, J. Shen, and P. Piecuch, Phys. Rev. Lett. 119, 223003 (2017) and J. E. Deustua, J. Shen, and P. Piecuch, J. Chem. Phys. 154, 124103 (2021)], we propose a novel form of the CC(P; Q) theory in which the stochastic Quantum Monte Carlo propagations, used to identify dominant higher-than-doubly excited determinants, are replaced by the selected configuration interaction (CI) approach using the perturbative selection made iteratively (CIPSI) algorithm. The advantages of the resulting CIPSI-driven CC(P; Q) methodology are illustrated by a few molecular examples, including the dissociation of F₂ and the automerization of cyclobutadiene, where we recover the electronic energies corresponding to the CC calculations with a full treatment of singles, doubles, and triples based on the information extracted from compact CI wave functions originating from relatively inexpensive Hamiltonian diagonalizations.

Spin-free formulation of the multireference driven similarity renormalization group: A benchmark study of first-row diatomic molecules and spin-crossover energetics

We report a spin-free formulation of the multireference (MR) driven similarity renormalization group (DSRG) based on the ensemble normal ordering of Mukherjee and Kutzelnigg [J. Chem. Phys. 107, 432 (1997)]. This ensemble averages over all microstates of a given total spin quantum number, and therefore, it is invariant with respect to SU(2) transformations. As such, all equations may be reformulated in terms of spin-free quantities and they closely resemble those of spin-adapted closed-shell coupled cluster (CC) theory. The current implementation is used to assess the accuracy of various truncated MR-DSRG methods (perturbation theory up to third order and iterative methods with single and double excitations) in computing the constants of 33 first-row diatomic molecules. The accuracy trends for these first-row diatomics are consistent with our previous benchmark on a small subset of closed-shell diatomic molecules. We then present the first MR-DSRG application on transition-metal complexes by computing the spin splittings of the [Fe(H₂O)₆]²⁺ and [Fe(NH₃)₆]²⁺ molecules. A focal point analysis (FPA) shows that third-order perturbative corrections are essential to achieve reasonably converged energetics. The FPA based on the linearized MR-DSRG theory with one- and two-body operators and up to a quintuple-ζ basis set predicts the spin splittings of [Fe(H₂O)₆]²⁺ and [Fe(NH₃)₆]²⁺ to be -35.7 and -17.1 kcal mol^-1, respectively, showing good agreement with the results of local CC theory with singles, doubles, and perturbative triples.

Toward the Minimal Floating Operation Count Cholesky Decomposition of Electron Repulsion Integrals

As quantum chemistry calculations deal with molecular systems of increasing size, the memory requirement to store electron-repulsion integrals (ERIs) greatly outpaces the physical memory available in computing hardware. The Cholesky decomposition of ERIs provides a convenient yet accurate technique to reduce the storage requirement of integrals. Recent developments of a two-step algorithm have drastically reduced the memory operation (MOP) count, leaving the floating operation (FLOP) count as the last frontier of cost reduction in the Cholesky ERI algorithm. In this report, we introduce a dynamic integral tracking, reusing, and compression/elimination protocol embedded in the two-step Cholesky ERI method. Benchmark studies suggest that this technique becomes particularly advantageous when the basis set consists of many computationally expensive high-angular-momentum basis functions. With this dynamic-ERI improvement, the Cholesky ERI approach proves to be a highly efficient algorithm with minimal FLOP and MOP count.

High-level coupled-cluster energetics by Monte Carlo sampling and moment expansions: Further details and comparisons

We recently proposed a novel approach to converging electronic energies equivalent to high-level coupled-cluster (CC) computations by combining the deterministic CC(P;Q) formalism with the stochastic configuration interaction (CI) and CC Quantum Monte Carlo (QMC) propagations. This article extends our initial study [J. E. Deustua, J. Shen, and P. Piecuch, Phys. Rev. Lett. 119, 223003 (2017)], which focused on recovering the energies obtained with the CC method with singles, doubles, and triples (CCSDT) using the information extracted from full CI QMC and CCSDT-MC, to the CIQMC approaches truncated at triples and quadruples. It also reports our first semi-stochastic CC(P;Q) calculations aimed at converging the energies that correspond to the CC method with singles, doubles, triples, and quadruples (CCSDTQ). The ability of the semi-stochastic CC(P;Q) formalism to recover the CCSDT and CCSDTQ energies, even when electronic quasi-degeneracies and triply and quadruply excited clusters become substantial, is illustrated by a few numerical examples, including the F-F bond breaking in F₂, the automerization of cyclobutadiene, and the double dissociation of the water molecule.

Coupled Cluster Externally Corrected by Adaptive Configuration Interaction

An externally corrected coupled cluster (CC) method, where an adaptive configuration interaction (ACI) wave function provides the external cluster amplitudes, named ACI-CC, is presented. By exploiting the connection between configuration interaction and CC through cluster analysis, the higher-order T₃ and T₄ terms obtained from ACI are used to augment the T₁ and T₂ amplitude equations from traditional CC. These higher-order contributions are kept frozen during the CC iterations and do not contribute to an increased cost with respect to coupled cluster including the single and double excitations (CCSD). We have benchmarked this method on three closed-shell systems: beryllium dimer, carbonyl oxide, and cyclobutadiene, with good results compared to other corrected CC methods. In all cases, the inclusion of these external corrections improved upon the "gold standard" CCSD(T) results, indicating that ACI-CCSD(T) can be used to assess strong correlation effects in a system and as an inexpensive starting point for more complex external corrections.

Connected three-body terms in single-reference unitary many-body theories: Iterative and perturbative approximations

This work introduces various approaches to include connected three-body terms in unitary many-body theories, focusing on the driven similarity renormalization group (DSRG). Starting from the least approximate method-the linearized DSRG truncated to one-, two-, and three-body operators [LDSRG(3)]-we develop several approximate LDSRG(3) models with reduced computational cost. Through a perturbative analysis, we motivate a family of iterative LDSRG(3)-n and -n' (n = 1, 2, 3, 4) methods that contain a subset of the LDSRG(3) diagrams. Among these variants, the LDSRG(3)-2 scheme has the same computational complexity of coupled cluster theory with singles, doubles, and triples (CCSDT), but it outperforms CCSDT in the accuracy of the predicted correlation energies. We also propose and implement two perturbative triples corrections based on the linearized DSRG truncated to one- and two-body operators augmented with recursive semi-quadratic commutators [qDSRG(2)]. The resulting qDSRG(2)+(T) approach matches the accuracy of the "gold-standard" coupled cluster theory with singles, doubles, and perturbative triples model on the energetics of twenty-eight closed-shell atoms and small molecules.

Psi4 1.4: Open-source software for high-throughput quantum chemistry

PSI4 is a free and open-source ab initio electronic structure program providing implementations of Hartree-Fock, density functional theory, many-body perturbation theory, configuration interaction, density cumulant theory, symmetry-adapted perturbation theory, and coupled-cluster theory. Most of the methods are quite efficient, thanks to density fitting and multi-core parallelism. The program is a hybrid of C++ and Python, and calculations may be run with very simple text files or using the Python API, facilitating post-processing and complex workflows; method developers also have access to most of PSI4's core functionalities via Python. Job specification may be passed using The Molecular Sciences Software Institute (MolSSI) QCSCHEMA data format, facilitating interoperability. A rewrite of our top-level computation driver, and concomitant adoption of the MolSSI QCARCHIVE INFRASTRUCTURE project, makes the latest version of PSI4 well suited to distributed computation of large numbers of independent tasks. The project has fostered the development of independent software components that may be reused in other quantum chemistry programs.