Linear-Scaling Open-Shell MP2 Approach: Algorithm, Benchmarks, and Large-Scale Applications

Overview of Developments in the MRCC Program System

mrcc is a versatile suite of quantum chemistry programs designed for accurate ab initio and density functional theory (DFT) calculations. This contribution outlines the general features and recent developments of the package. The most popular features include the open-ended coupled-cluster (CC) code, state-of-the-art CC singles and doubles with perturbative triples [CCSD(T)], second-order algebraic-diagrammatic construction, and combined wave function theory-DFT approaches. Cost-reduction techniques are implemented, such as natural orbital (NO), local NO (LNO), and natural auxiliary function approximations, which significantly decrease the computational demands of these methods. This paper also details the method developments made over the past five years, including efficient schemes to approach the complete basis set limit for CCSD(T) and the extension of our LNO-CCSD(T) method to open-shell systems. Additionally, we discuss the new approximations introduced to accelerate the self-consistent field procedure and the cost-reduction techniques elaborated for analytic gradient calculations at various levels. Furthermore, embedding techniques and novel range-separated double-hybrid functionals are presented for excited-state calculations, while the extension of the theories established to describe core excitations and ionized states is also discussed. For academic purposes, the program and its source code are available free of charge, and its commercial use is also facilitated.

Advancing Non-Atom-Centered Basis Methods for More Accurate Interaction Energies: Benchmarks and Large-Scale Applications

Recent advances in local electron correlation approaches have enabled the relatively routine access to CCSD(T) [that is, coupled cluster (CC) with single, double, and perturbative triple excitations] computations for molecules of a hundred or more atoms. Here, approaching their complete basis set (CBS) limit becomes more challenging due to extensive basis set superposition errors, often necessitating the use of large atomic orbital (AO) basis sets with diffuse functions. Here, we study a potential remedy in the form of non-atom-centered or floating orbitals (FOs). FOs are still rarely employed even for small molecules due to the practical complication of defining their position, number, exponents, etc. The most frequently used FO method thus simply places a single FO center with a large number of FOs toward the middle of noncovalent dimers; however, a single FO center for larger complexes can soon become insufficient. A recent alternative uses a grid of FO centers around the monomers with a single s function per center, which is currently applicable only for H, C, N, and O atoms. Here, we build on the above advantages and mitigate some drawbacks of previous FO approaches by using a layer of FO centers and 4-9 FOs/center for each monomer. Thus, a double layer of FOs is placed between the interacting subsystems. When extending the double-ζ AO basis with this double layer of FOs, the quality of conventional augmented double-ζ or conventional triple-ζ AO bases can be reached or surpassed with less orbitals, leading to few tenths of a kcal/mol basis set errors for medium-sized dimers. This good performance extends to larger molecules (shown here up to 72 atoms), as efficient local natural orbital (LNO) CCSD(T) computations with only double-ζ AO and 4 FOs/center FO bases match our LNO-CCSD(T)/CBS reference within ca. 0.1 kcal/mol. These developments introduce FO methods to the accurate modeling of large molecular complexes without limitations to atom types by further accelerating efficient correlation calculations, like LNO-CCSD(T).

Basis-Set Limit CCSD(T) Energies for Large Molecules with Local Natural Orbitals and Reduced-Scaling Basis-Set Corrections

The calculation of density-based basis-set correction (DBBSC), which remedies the basis-set incompleteness (BSI) error of the correlation energy, is combined with local approximations. Aiming at large-scale applications, the procedure is implemented in our efficient local natural orbital-based coupled-cluster singles and doubles with perturbative triples [LNO-CCSD(T)] scheme. To this end, the range-separation function, which characterizes the one-electron BSI in space, is decomposed into the sum of contributions from individual localized molecular orbitals (LMOs). A compact domain is constructed around each LMO, and the corresponding contributions are evaluated only within these restricted domains. Furthermore, for the calculation of the complementary auxiliary basis set (CABS) correction, which significantly improves the Hartree-Fock (HF) energy, the local density fitting approximation is utilized. The errors arising from the local approximations are examined in detail, efficient prescreening techniques are introduced to compress the numerical quadrature used for DBBSC, and conservative default thresholds are selected for the truncation parameters. The efficiency of the DBBSC-LNO-CCSD(T) method is demonstrated through representative examples of up to 1000 atoms. Based on the numerical results, we conclude that the corrections drastically reduce the BSI error using double-ζ basis sets, often to below 1 kcal/mol compared to the reliable LNO-CCSD(T) complete basis set references, while significant improvements are also achieved with triple-ζ basis sets. Considering that the calculation of the DBBSC and CABS corrections only moderately increases the wall-clock time required for the post-HF steps in practical applications, the proposed DBBSC-LNO-CCSD(T) method offers a highly efficient and robust tool for large-scale calculations.

State-of-the-art local correlation methods enable affordable gold standard quantum chemistry for up to hundreds of atoms

In this feature, we review the current capabilities of local electron correlation methods up to the coupled cluster model with single, double, and perturbative triple excitations [CCSD(T)], which is a gold standard in quantum chemistry. The main computational aspects of the local method types are assessed from the perspective of applications, but the focus is kept on how to achieve chemical accuracy (i.e., <1 kcal mol-1 uncertainty), as well as on the broad scope of chemical problems made accessible. The performance of state-of-the-art methods is also compared, including the most employed DLPNO and, in particular, our local natural orbital (LNO) CCSD(T) approach. The high accuracy and efficiency of the LNO method makes chemically accurate CCSD(T) computations accessible for molecules of hundreds of atoms with resources affordable to a broad computational community (days on a single CPU and 10-100 GB of memory). Recent developments in LNO-CCSD(T) enable systematic convergence and robust error estimates even for systems of complicated electronic structure or larger size (up to 1000 atoms). The predictive power of current local CCSD(T) methods, usually at about 1-2 order of magnitude higher cost than hybrid density functional theory (DFT), has become outstanding on the palette of computational chemistry applicable for molecules of practical interest. We also review more than 50 LNO-based and other advanced local-CCSD(T) applications for realistic, large systems across molecular interactions as well as main group, transition metal, bio-, and surface chemistry. The examples show that properly executed local-CCSD(T) can contribute to binding, reaction equilibrium, rate constants, etc. which are able to match measurements within the error estimates. These applications demonstrate that modern, open-access, and broadly affordable local methods, such as LNO-CCSD(T), already enable predictive computations and atomistic insight for complicated, real-life molecular processes in realistic environments.

Chemical Potentials and the One-Electron Hamiltonian of the Second-Order Perturbation Theory from the Functional Derivative Approach

We develop a functional derivative approach to calculate the chemical potentials of second-order perturbation theory (MP2). In the functional derivative approach, the correlation part of the MP2 chemical potential, which is the derivative of the MP2 correlation energy with respect to the occupation number of frontier orbitals, is obtained from the chain rule via the noninteracting Green's function. First, the MP2 correlation energy is expressed in terms of the noninteracting Green's function, and its functional derivative to the noninteracting Green's function is the second-order self-energy. Then, the derivative of the noninteracting Green's function to the occupation number is obtained by including the orbital relaxation effect. We show that the MP2 chemical potentials obtained from the functional derivative approach agree with that obtained from the finite difference approach. The one-electron Hamiltonian, defined as the derivative of the MP2 energy with respect to the one particle density matrix, is also derived using the functional derivative approach, which can be used in the self-consistent calculations of MP2 and double-hybrid density functionals. The developed functional derivative approach is promising for calculating the chemical potentials and the one-electron Hamiltonian of approximate functionals and many-body perturbation approaches dependent explicitly on the noninteracting Green's function.

Toward Real Chemical Accuracy on Current Quantum Hardware Through the Transcorrelated Method

Quantum computing is emerging as a new computational paradigm with the potential to transform several research fields including quantum chemistry. However, current hardware limitations (including limited coherence times, gate infidelities, and connectivity) hamper the implementation of most quantum algorithms and call for more noise-resilient solutions. We propose an explicitly correlated Ansatz based on the transcorrelated (TC) approach to target these major roadblocks directly. This method transfers, without any approximation, correlations from the wave function directly into the Hamiltonian, thus reducing the resources needed to achieve accurate results with noisy quantum devices. We show that the TC approach allows for shallower circuits and improves the convergence toward the complete basis set limit, providing energies within chemical accuracy to experiment with smaller basis sets and, thus, fewer qubits. We demonstrate our method by computing bond lengths, dissociation energies, and vibrational frequencies close to experimental results for the hydrogen dimer and lithium hydride using two and four qubits, respectively. To demonstrate our approach's current and near-term potential, we perform hardware experiments, where our results confirm that the TC method paves the way toward accurate quantum chemistry calculations already on today's quantum hardware.

Linear-Scaling Local Natural Orbital CCSD(T) Approach for Open-Shell Systems: Algorithms, Benchmarks, and Large-Scale Applications

The extension of the highly optimized local natural orbital (LNO) coupled cluster (CC) with single-, double-, and perturbative triple excitations [LNO-CCSD(T)] method is presented for high-spin open-shell molecules based on restricted open-shell references. The techniques enabling the outstanding efficiency of the closed-shell LNO-CCSD(T) variant are adopted, including the iteration- and redundancy-free second-order Møller-Plesset and (T) formulations as well as the integral-direct, memory- and disk use-economic, and OpenMP-parallel algorithms. For large molecules, the efficiency of our open-shell LNO-CCSD(T) method approaches that of its closed-shell parent method due to the application of restricted orbital sets for demanding integral transformations and a novel approximation for higher-order long-range spin-polarization effects. The accuracy of open-shell LNO-CCSD(T) is extensively tested for radicals and reactions thereof, ionization processes, as well as spin-state splittings, and transition-metal compounds. At the size range where the canonical CCSD(T) reference is accessible (up to 20-30 atoms), the average open-shell LNO-CCSD(T) correlation energies are found to be 99.9 to 99.95% accurate, which translates into average absolute deviations of a few tenths of kcal/mol in the investigated energy differences already with the default settings. For more extensive molecules, the local errors may grow, but they can be estimated and decreased via affordable systematic convergence studies. This enables the accurate modeling of large systems with complex electronic structures, as illustrated on open-shell organic radicals and transition-metal complexes of up to 179 atoms as well as on challenging biochemical systems, including up to 601 atoms and 11,000 basis functions. While the protein models involve difficulties for local approximations, such as the spin states of a bounded iron ion or an extremely delocalized singly occupied orbital, the corresponding single-node LNO-CCSD(T) computations were feasible in a matter of days with 10s to 100 GB of memory use. Therefore, the new LNO-CCSD(T) implementation enables highly accurate computations for open-shell systems of unprecedented size and complexity with widely accessible hardware.

Local Second-Order Møller-Plesset Theory with a Single Threshold Using Orthogonal Virtual Orbitals: Theory, Implementation, and Assessment

It has long been clear that electron correlation methods exhibit unphysical compute scalings with molecular size, which has motivated the development of local correlation methods to discard effectively zero contributions in a controlled way to yield an approximate correlation energy. The ideal local correlation method should have a single numerical threshold that controls the dropping of terms with the ability to have that threshold set small enough so that the correlation energy is reproduced to enough significant figures such that the result is chemically identical. This work reports such a method for the second-order Møller-Plesset (MP2) theory. The theory, implementation, and testing of this local MP2 theory are reported. Thresholds ranging from 10-5 to 10-8 and basis sets ranging from split valence plus polarization through to quadruple-ζ are assessed for local MP2 calculations on a range of molecules, including linear chains and molecules with two- and three-dimensional character. The implementation is shared memory parallel via OpenMP and yields roughly 50% parallel efficiency with 16 cores for a large job. Considerable efforts were made to minimize memory demands, which increased as thresholds were tightened. A variety of relative energy calculations are presented as a function of threshold to provide some guidance to users on how to obtain adequate precision at a low compute cost. It is particularly clear that derivative properties require tighter thresholds in order to achieve an adequate precision.

Corrigendum: Coupled cluster theory on modern heterogeneous supercomputers

[This corrects the article DOI: 10.3389/fchem.2023.1154526.].

Coupled cluster theory on modern heterogeneous supercomputers

This study examines the computational challenges in elucidating intricate chemical systems, particularly through ab-initio methodologies. This work highlights the Divide-Expand-Consolidate (DEC) approach for coupled cluster (CC) theory-a linear-scaling, massively parallel framework-as a viable solution. Detailed scrutiny of the DEC framework reveals its extensive applicability for large chemical systems, yet it also acknowledges inherent limitations. To mitigate these constraints, the cluster perturbation theory is presented as an effective remedy. Attention is then directed towards the CPS (D-3) model, explicitly derived from a CC singles parent and a doubles auxiliary excitation space, for computing excitation energies. The reviewed new algorithms for the CPS (D-3) method efficiently capitalize on multiple nodes and graphical processing units, expediting heavy tensor contractions. As a result, CPS (D-3) emerges as a scalable, rapid, and precise solution for computing molecular properties in large molecular systems, marking it an efficient contender to conventional CC models.

Stability of Carbocyclic Phosphinyl Radicals: Effect of Ring Size, Delocalization, and Sterics

In this computational study, we report on the stability of cyclic phosphinyl radicals with an aim for a systematical assessment of stabilization effects. The radical stabilization energies (RSEs) were calculated using isodesmic reactions for a large number of carbocyclic radicals possessing different ring sizes and grades of unsaturation. In general, the RSE values range from −1.2 to −14.0 kcal·mol^–1, and they show practically no correlation with the spin populations at the P-centers. The RSE values correlate with the reaction Gibbs free energies calculated for the dimerization of the studied simple radicals. Therefore, the more easily accessible RSE values offer a cost-effective estimation of global stability in a straightforward manner. To explore the effect of unsaturation on the RSE values, delocalization energies were determined using appropriate isodesmic reactions. Introducing unsaturations beside the P-center into the backbone of the rings leads to an additive increase in the magnitude of the delocalization energy (∼10, 20, and 30 kcal·mol^–1, respectively, for radicals with one, two, and three C=C bonds in the conjugation). Parallelly, the spin populations at the P-centers also dwindle gradually by ∼0.1 e in the same order, indicating that the lone electron delocalizes over the π-system. Radicals containing exocyclic C=C π-bonds were also investigated, and all of these radicals have rather similar stabilities independently of the ring size, outlining the primary importance of the two exocyclic π-bonds in the conjugation. Among the radicals involved in our study, those with the best electronic stabilization are the unsaturated three-, five-, six-, and seven-membered rings containing the maximum number of conjugated vinyl fragments. The largest delocalization energy of 31.5 kcal·mol^–1 and the lowest obtained spin population of 0.665 e were found for the fully unsaturated seven-membered radical (phosphepin derivative). Importantly, the electronic stabilization effects alone are insufficient for stabilizing the radicals in monomeric forms epitomized by the exothermic dimerization energies (−40 to −58 kcal·mol^–1). Therefore, it is essential to apply sterically demanding bulky substituents on the α-C-atoms. Tweaking the steric congestion enabled us to propose radicals that are expected to be stable against dimerization and, consequently, may be realistic target species for synthetic investigations. The effects contributing to the stability of radicals having sterically encumbered substituents have also been explored.

To systematically evaluate the stabilization effects, the radical stabilization energies of various carbocyclic phosphinyl radicals having saturated backbones or unsaturation(s) in either endocyclic or exocyclic manner have been determined and analyzed. As the electronic stabilization is alone insufficient to hamper the possible dimerization of these species, the effect of several sterically demanding substituents has been explored for the congeners with best electronic stabilizations, thus enabling us to propose synthetically accessible candidates in the future.

Quasi-Restricted Orbital Description of the Copper(I) Photoredox Catalytic Cycle

In this computational study, the electronic structure changes along the oxidative and reductive quenching cycles of a homoleptic and a heteroleptic prototype Cu(I) photoredox catalyst, namely [Cu(dmp)2]+ (dmp = 2,9-dimethyl-1,10-phenanthroline) and [Cu(phen)(POP)]+ (POP = bis[2-(diphenylphosphino)phenyl]ether) are scrutinized and characterized using quasi-restricted orbitals (QRO), electron density differences and spin densities. After validating our density functional theory-based computational protocol, the equilibrium geometries and wavefunctions (using QROs and atom/fragment compositions) of the four states involved in photoredox cycle (S0, T1, Dox and Dred) are systematically and thoroughly described. The formal ground and excited state ligand- and metal-centered redox events are substantiated by the QRO description of the open-shell triplet 3MLCT (d9L-1), Dox (d9L0) and Dred (d10L-1) species and the corresponding structural changes, e.g., flattening distortion, shortening/elongation of Cu-N/Cu-P bonds, are rationalized in terms of the underlying electronic structure transformations. Amongst others, we reveal the molecular-scale delocalization of the ligand-centered radical in the a 3MLCT (d9L-1) and Dred (d9L-1) states of homoleptic [Cu(dmp)2]+ and its localization to the redox-active phenanthroline ligand in the case of heteroleptic [Cu(phen)(POP)]+.

Orbital Pair Selection for Relative Energies in the Domain-Based Local Pair Natural Orbital Coupled-Cluster Method

For the accurate computation of relative energies, domain-based local pair natural orbital coupled-cluster [DLPNO-CCSD(T0)] has become increasingly popular. Even though DLPNO-CCSD(T0) shows a formally linear scaling of the computational effort with the system size, accurate predictions of relative energies remain costly. Therefore, multi-level approaches are attractive that focus the available computational resources on a minor part of the molecular system, e.g., a reaction center, where changes in the correlation energy are expected to be the largest. We present a pair-selected multi-level DLPNO-CCSD(T0) ansatz that automatically partitions the orbital pairs according to their contribution to the overall correlation energy change in a chemical reaction. To this end, the localized orbitals are mapped between structures in the reaction; all pair energies are approximated through computationally efficient semi-canonical second-order Møller--Plesser perturbation theory, and the orbital pairs for which the pair energies change significantly are identified. This multi-level approach is significantly more robust than our previously suggested, orbital selection-based multi-level DLPNO-CCSD(T0) ansatz [ J. Chem. Phys. 2021, 155, 224102] for reactions showing only small changes in the occupied orbitals. At the same time, it is even more efficient without added input complexity or accuracy loss compared to the full DLPNO-CCSD(T0) calculation. We demonstrate the accuracy of the multi-level approach for a total of 128 chemical reactions and potential energy curves of weakly interacting complexes from the S66x8 benchmark set.

A combined first- and second-order optimization method for improving convergence of Hartree–Fock and Kohn–Sham calculations

In this work, we investigate the optimization of Hartree–Fock (HF) orbitals with our recently proposed combined first- and second-order (SO-SCI) method, which was originally developed for multi-configuration self-consistent field (MCSCF) and complete active space SCF (CASSCF) calculations. In MCSCF/CASSCF, it unites a second-order optimization of the active orbitals with a Fock-based first-order treatment of the remaining closed-virtual orbital rotations. In the case of the single-determinant wavefunctions, the active space is replaced by a preselected “second-order domain,” and all rotations involving orbitals in this subspace are treated at second-order. The method has been implemented for spin-restricted and spin-unrestricted Hartree–Fock (RHF, UHF), configuration-averaged Hartree–Fock (CAHF), as well as Kohn–Sham (KS) density functional theory (RKS, UKS). For each of these cases, various choices of the second-order domain have been tested, and appropriate defaults are proposed. The performance of the method is demonstrated for several transition metal complexes. It is shown that the SO-SCI optimization provides faster and more robust convergence than the standard SCF procedure but requires, in many cases, even less computation time. In difficult cases, the SO-SCI method not only speeds up convergence but also avoids convergence to saddle-points. Furthermore, it helps to find spin-symmetry broken solutions in the cases of UHF or UKS. In the case of CAHF, convergence can also be significantly improved as compared to a previous SCF implementation. This is particularly important for multi-center cases with two or more equal heavy atoms. The performance is demonstrated for various two-center complexes with different lanthanide atoms.

Recent Advances in First-Principles Based Molecular Dynamics

First-principles molecular dynamics (FPMD) and its quantum mechanical-molecular mechanical (QM/MM) extensions are powerful tools to follow the real-time dynamics of a broad variety of systems in their ground as well as electronically excited states. The continued advances in computational power have enabled simulations of QM regions of larger sizes for more extended time scales. In addition, development of the parallel algorithms has boosted the performance of QM/MM methods even on existing computer architectures. In the case of density functional-based FPMD, systems of several hundreds to thousands of atoms can now be customarily simulated for tens to hundreds of picoseconds. In spite of this progress, the time scale limitations remain severe, especially when high-rung exchange-correlation functionals or high-level wave function based quantum mechanical methods are used. To ameliorate this, a large number of enhanced sampling methods have been introduced but most of the approaches that have been developed to increase the efficiency of FPMD based simulations sacrifice the real-time dynamics in favor of enhancing sampling. Here, we present some recent advances in boosting the efficiency of FPMD based simulations while keeping the full dynamic information. These include a highly efficient recent implementation of FPMD-based QM/MM simulations that not only enables fully flexible combinations of different electronic structure methods and force fields via a highly efficient communication library, it also fully exploits parallelism for both quantum and classical descriptions. The second type of acceleration methods we discuss is a large family of specially devised multiple-time-step algorithms that make use of suitable breakups of the total nuclear forces into fast components that can be calculated via lower level methods and slowly varying correction forces evaluated with a high-level method at long time intervals. The computational gain of this scheme mostly depends on the cost difference between the two methods and advantageous combinations can yield large speedups without compromising the accuracy of the high-level method. And finally, the third class of FPMD acceleration methods presented here are machine learning models to accelerated FPMD and their powerful combinations with multiple-time-step techniques. The combination of all the approaches enables substantial speedups of FPMD simulations of several orders of magnitude while fully preserving the real-time dynamics and accuracy.

Direct orbital selection within the domain-based local pair natural orbital coupled-cluster method

Domain-based local pair natural orbital coupled cluster (DLPNO-CC) has become increasingly popular to calculate relative energies (e.g., reaction energies and reaction barriers). It can be applied within a multi-level DLPNO-CC-in-DLPNO-CC ansatz to reduce the computational cost and focus the available computational resources on a specific subset of the occupied orbitals. We demonstrate how this multi-level DLPNO-CC ansatz can be combined with our direct orbital selection (DOS) approach [M. Bensberg and J. Neugebauer, J. Chem. Phys. 150, 214106 (2019)] to automatically select orbital sets for any multi-level calculation. We find that the parameters for the DOS procedure can be chosen conservatively such that they are transferable between reactions. The resulting automatic multi-level DLPNO-CC method requires no user input and is extremely robust and accurate. The computational cost is easily reduced by a factor of 3 without sacrificing accuracy. We demonstrate the accuracy of the method for a total of 61 reactions containing up to 174 atoms and use it to predict the relative stability of conformers of a Ru-based catalyst.

Coupled Cluster Benchmark of New DFT and Local Correlation Methods: Mechanisms of Hydroarylation and Oxidative Coupling Catalyzed by Ru(II, III) Chloride Carbonyls

We have evaluated a set of accurate canonical CCSD(T) energies for stationary points on the potential energy surface for Ru(II, III) chloride carbonyl catalysis of two competing reactions between benzene and methyl acrylate (MA), namely, hydroarylation and oxidative coupling. We have then applied this set to evaluate the performance of localized orbital coupled-cluster methods and several new and common density functionals. We find that (a) DLPNO-CCSD(T) with TightPNO cutoffs is an acceptable substitute for full canonical CCSD(T) calculations on this system; (b) for the closed-shell systems where it could be applied, LNO-CCSD(T) with tight convergence criteria is very close to the canonical results; (c) the recent ωB97X-V and ωB97M-V functionals exhibit superior performance to commonly used DFT functionals in both closed- and open-shell calculations; (d) the revDSD-PBEP86 revision of the DSD-PBEP86 double hybrid represents an improvement over the original, even though transition metals were not involved in its parametrization; and (e) DSD-SCAN and DOD-SCAN show comparable efficiency. Most tested (meta)-GGA and hybrid density functionals perform better for open-shell than for closed-shell complexes; this is not the case for the double hybrids considered.