Tensor-Hypercontracted MP2 First Derivatives: Runtime and Memory Efficient Computation of Hyperfine Coupling Constants

Authors: Felix H. Bangerter, Michael Glasbrenner, Christian Ochsenfeld ⚬

Affiliations: †Chair of Theoretical Chemistry, Department of Chemistry, University of Munich (LMU), D-81377 Munich, Germany; ‡Max Planck Institute for Solid State Research, D-70569 Stuttgart, Germany

License: © 2022 The Authors. Published by American Chemical Society CC BY 4.0 Permits the broadest form of re-use including for commercial purposes, provided that author attribution and integrity are maintained (https://creativecommons.org/licenses/by/4.0/).

Article links: DOI: 10.1021/acs.jctc.2c00118  |  PubMed: 35943450  |  PMC: PMC9476664

Relevance: Relevant: mentioned in keywords or abstract

Full text: PDF (548 KB)

Introduction

Ever since the advent of modern computers in the 1990s, Møller–Plesset perturbation theory (MPn)¹ has been a good compromise in the family of quantum chemical methods, being sufficiently accurate for many applications while still being computationally affordable.² As opposed to coupled cluster theory (CC), MPn lacks infinite-order corrections present in the cluster operator expansion of the CC models, which generally makes MPn less accurate.² However, when going from energies to gradients and molecular properties, MPn, especially second-order MPn (MP2), was shown to yield accurate hyperfine coupling constants (HFCCs)³⁻⁷ and relative nuclear magnetic resonance (NMR) shifts.⁸⁻¹² However, MP2 is sensitive to spin-contamination in the Hartree–Fock wave function, which can be improved upon when used in its orbital-optimized variant⁴ or as part of double-hybrid density functionals.¹³

Furthermore, when comparing MPn and CC at the same expansion orders, for example, MP2 and singles and doubles CC (CCSD), MPn comes with a scaling advantage, both in the prefactor and the scaling exponent. Nonetheless, canonical MP2, as well as the associated first and second derivative,^6,11 still scale with the fifth power of the molecule size, thereby severely restricting the accessible chemical space. To alleviate this limitation several formulations of the MP2 derivatives have been proposed.

Early work from Pulay and Sæbø^14,15 on local correlation was applied to the computation of MP2 gradients.^16,17 Following up on this, in recent years the domain-based local pair natural orbital (DLPNO) formulation of MP2 by Neese and co-workers¹⁸⁻²⁰ was extended to first and second derivatives. Conceptually related is the divide-expand-consolidate (DEC) formulation of MP2 by Jørgensen and co-workers,^21,22 which likewise was extended to the computation of molecular gradients in a linear-scaling and massively parallel manner.^23,24 Instead of exploiting locality in the correlation space, the equations of the MP2 first^6,25 and second derivative¹¹ can be reformulated entirely in terms of atomic orbitals (AOs). However, for an efficient implementation and the reduction of the scaling prefactor, an orbital localization by pivoted Cholesky decomposition (PCD) of the associated pseudodensities is essential.⁶

Besides ensuring low-scaling and efficiency, for derivative calculations of electron correlation methods in general, it is pivotal to efficiently manage the available memory and disk space. Compared to the energy equations, the associated gradients and higher derivatives often not only include electron repulsion integrals (ERIs), in either atomic or molecular orbital basis, but can also include partially transformed ERIs and derivatives thereof, for example, with respect to the magnetic field in NMR calculations. Since canonical ERIs are fourth-order tensors, their memory requirements prohibitively scale with the forth power of the number of basis functions and saving multiple such tensors quickly becomes unfeasible. To reduce the memory footprint of the ERIs, tensor decomposition methods, particularly the resolution-of-the-identity (RI) ansatz,²⁶⁻²⁹ have been broadly applied in the context of MP2 derivatives.^5,6,12,30 However, with increasing molecule size, even the third-order RI tensors eventually exceed the available disk space of conventional high performance computing nodes. To overcome this storage limitation, further reduction of the dimensionality of the ERIs is desirable. This can be achieved by the recently introduced tensor hypercontraction (THC) factorization of Martínez and co-workers,³¹⁻³⁴ which in general terms approximates a fourth-order integral tensor (μν|ô|λσ), where ô is a singular two electron interaction kernel, by five second-order tensors.³⁵ In the least-squares formulation of THC (LS-THC),³¹ four of these tensors are simply obtained by evaluation of the basis functions at real-space grid nodes and the singular ô operator is replaced by the LS-fitted Z matrix. If ô is the Coulombic 1/r operator, a factorization of the regular ERIs is achieved, which has been employed in reduced scaling formulations of exact exchange,³⁶ different orders of MPn,^31,33−41 the random phase approximation (RPA),⁴² complete active space perturbation theory (CASPT2),⁴³ and various flavors of CC theory,⁴⁴⁻⁴⁶ as well as equation-of-motion CC (EOM-CC) theory.⁴⁷ Recently Matthews³⁵ thoroughly investigated amplitude factorizations within MP3, as a stepping stone toward CCSD, and noted that the LS-THC factorization of nonlocal integrals, such as the exchange integrals, incurs an additional error.

The applicability of THC for the computation of molecular gradients is largely unexplored. Song et al.³⁹ derived equations for the analytical gradient of THC-AO-MP2 with application to geometry optimizations and ab initio molecular dynamics (AIMD) simulations. As commonly the derivation of gradient equations of electron correlation methods is rather involved, Song et al.⁴⁸ also proposed an automatic differentiation scheme for the automated generation of working equations for gradients of THC-based correlation methods.

When it comes to applying tensor factorizations, be it RI or THC, two possible routes to molecular gradients can be taken: The first approach is to differentiate the RI- or THC-approximated energy equation, then the associated gradient describes the slope of an approximated potential energy surface (PES). The second approach is to take the gradient of the canonical energy and insert the approximation into the exact gradient; this way an approximate gradient is used to move along the exact PES.

The latter can lead to an unwanted buildup of errors during the course of a simulation, when used in molecular dynamics simulations. However, when used in conjunction with thermostats constant energy is traded for constant temperature, and depending on the extent of the gradient error, this approach can still be applicable but must be tested for the chosen gradient. The first approach was taken by Song et al.³⁹ for their THC-based MP2 molecular gradient, whereas in the present work, the gradient equations with respect to a perturbation, that the basis functions are independent of, are derived by inserting the THC factorization into the equations of the exact gradient. As will be discussed in section 2.4, the resulting equations are identical to the ones obtained by differentiating the THC-AO-MP2 energy equation. As a representative case of these kinds of perturbations, we apply our recently developed low-scaling LS-THC algorithm⁴¹ to the AO-MP2 energy derivative with respect to the nuclear magnetic moment, for the computation of HFCCs on the MP2 level of theory.

We present ways to efficiently treat the Coulomb- and the exchange-like part of the most expensive intermediate of the derivative within the LS-THC approximation. We demonstrate the low-scaling behavior of our THC-ω-RI-CDD-MP2 derivative method for various chemically relevant systems. We also show that THC-ω-RI-CDD-MP2 significantly outperforms our previous implementation of ω-RI-CDD-MP2^6,30 for the computation of isotropic HFCCs.

Theory

Notation

Throughout this publication we make use of the following indices:

• μ, ν, λ, σ: atomic orbital indices belonging to the AO basis {χ_μ} of size N_bf.
• α, β, γ, δ: auxiliary basis indices belonging to the density fitting basis {χ_α} of size N_aux (usually N_aux ≈ 3·N_bf).
• P, Q, R, S: grid point indices belonging to the LS-THC grid of size N_grid (usually N_grid ≈ 3·N_aux).
• i, j: occupied molecular orbital indices belonging to the MO basis {ϕ_i} of size N_occ.
• a, b: virtual molecular orbital indices belonging to the MO basis {ϕ_a} of size N_virt (N_virt ≫ N_occ).
• η, η′: spin indices, for either α- or β-electrons, with η′ ≠ η.
• κ: index of the Laplace quadrature points for the MP2 energy denominator with weights ω_κ (usually integration with 5–8 points is sufficiently accurate).
• k: index of the nucleus under consideration.

Review of the AO-MP2 Gradient

The unrestricted AO-MP2 energy with a Laplace transformation⁴⁹⁻⁵¹ for the energy denominator is given as

with

for which the transformed ERIs are given by

and P̲ and P̅ are the usual pseudodensities, given by

To obtain the AO-MP2 gradient, eq 1 has to be differentiated with respect to a perturbation ξ. Since the focus of this work is on HFCCs, as an example for a property for which the basis functions are independent of the perturbation—here ξ′—the following derivation is restricted to this special case. An in-depth derivation of the AO-MP2 gradient equations, as well as a comparison to the MO-MP2 gradient, is available in refs (ref. 6) and (ref. 25).

In order to obtain the gradient of the AO-MP2 energy, eq 1 is differentiated with respect to ξ′

with

and intermediates R̲ and R̅ given by

The above R intermediates can be thought of as the contraction of all perturbation-independent parts of the gradient. Note that eq 5 only involves the derivative of the pseudodensities and no integral derivatives, as the basis functions are taken to be independent of ξ′. In order to avoid the evaluation of the derivatives of the perturbed occupied and virtual pseudodensities, these intermediates are expanded in terms of the regular occupied and virtual densities P_occ and P_virt as

Thus, the derivatives of eq 8 are given by

By further making use of the identity

the perturbed virtual density can be related to the perturbed occupied density as

Note again that eq 11 does not contain the derivative of the overlap matrix S, due to S being independent of ξ′. By making use of the above relations for the densities, eq 6 becomes

Cyclic permutation under the trace was applied above to obtain terms of the general form

with

which can be solved for Y by recursion, as detailed in refs (ref. 25) and (ref. 12). Let Y̅^η be the solution of eq 13 with A ≡ τ_κP_occ^ηF^η and B ≡ P_occ^ηR̅^η, and let Y̲^η be the solution with A ≡ – τ_κP_virt^ηF^η and B ≡ P_virt^ηR̲^η. Then, by making use of the other transformations outlined above, the derivative from eq 6 can be rewritten as

with

where intermediate resembles the Fock matrix with substituting for the density matrix; that is,

Note that the first term in eq 15 only includes the derivative of the core Hamiltonian matrix as the integrals in the Fock matrix, from which this term originates from, are independent of ξ′.

Equation 15 permits an elegant solution for the perturbed density, avoiding the need to solve coupled-perturbed self-consistent field (CPSCF) equations for all perturbations ξ′, by means of applying a AO-based Z-vector-like method²⁵ originally proposed by Handy and Schaefer.⁵² The implicit first derivative of the occupied density can therefore be efficiently obtained by applying the density matrix-based Laplace-transform unrestricted CPSCF (DL-UCPSCF) method by Beer and Ochsenfeld.⁵³ The intricacies of this method are detailed in refs (ref. 53, ref. 25), and (ref. 6).

To conclude, eqs 5 and 15 yield the first derivative of the AO-MP2 energy with respect to a perturbation ξ′. More specifically, if ξ′ was an external electric field then eq 5 would yield permanent dipole moments, and if ξ′ was equal to the nuclear magnetic moment M_k of a given nucleus k then the isotropic contribution to the HFCC of nucleus k in the absence of spin–orbit coupling would be obtained. The latter property will be used as a sample property for the newly developed THC-ω-RI-CDD-MP2 derivative method presented in section 2.4.

RI-CDD-MP2 HFCCs

The computational bottleneck of obtaining the AO-MP2 gradient in eqs 5 and 15 are the integral contractions in the formation of the R matrix intermediates given by eq 7. Because forming the R-matrices involves the same contraction with the pseudodensities as the AO-MP2 energy, the computation of the gradient will a priori also have quintic scaling. As is common practice when dealing with these kinds of integral contractions, the RI approximation can be inserted into eq 7 to lower the computational cost as well as the memory requirements by avoiding the fourth-order ERI tensors. To further lower the prefactor of the integral transformations a PCD of the pseudodensities can be used, which is known as the Cholesky-decomposed pseudodensity (CDD) approach. The CDD method produces a set of so-called Cholesky pseudo-MO coefficient matrices L̲ and L̅ according to

which are used analogously to the regular MO coefficients. The combined RI-CDD approach leads to a reformulation of the integrals incorporated in the R intermediates as

with

In combination with QQR-type integral screening the RI-CDD-MP2 gradient method was shown to yield cubic scaling for the computation of molecular gradients and HFCCs.⁶ To further reduce the scaling, Vogler et al.³⁰ employed the attenuated Coulomb metric^54,55 in the RI approximation as well as the scaled-opposite spin (SOS)⁵⁶ approximation, which removes the same spin contribution entirely. Still, the expensive formation of the R intermediates has to be done for every Laplace point and thus constitutes the predominant part of the wall time for the evaluation of the MP2 gradient with respect to a perturbation ξ′, even with the ω-RI-CDD-SOS-MP2 method.^6,30

THC-CDD-MP2 HFCCs

ERIs are ubiquitous in electron correlation methods and their transformation and contraction usually represents the bottleneck of the calculation. This is especially the case when many different ERIs, that is, fully and partially transformed into the MO space or contracted with a perturbed density matrix, are needed, and their formation has to be carried out repeatedly, such as inside a Laplace expansion or during the iterative solution of amplitude equations. Since the scaling behavior of these operations is dependent on the dimensionality of the representation of the ERI tensor, a most compact representation is desirable. Of particular interest is thus the THC factorization, which in its AO formulation approximates an ERI as

and—in LS-THC—the X-matrices are simply obtained by evaluation of basis functions at the THC grid.³¹ The analytical expression of the Z-matrix can be shown to be the solution to the normal equations associated with the least-squares equation of finding the THC factorization (see the Supporting Information). As has been shown previously,^{33,34,36,39,41−43} the THC factorization can achieve major savings in computation time for intermediates involving ERI contractions, by reducing the representation of the ERIs to only second-order tensors. In this work, we make use of our recently reported low-scaling THC method⁴¹ based on the ω-RI approximation for the ERIs contained in Z and natural blocking (NB).^57,58 However, in contrast to our work on THC-MP2 energies,⁴¹ in the present work the AO ERIs are fitted, since the gradient for calculating the MP2 HFCCs is based on our RI-CDD approach to the computation of AO-MP2 energy gradients.⁶ It is important to note here, that while the equations for the THC-based gradient method are derived by inserting the THC factorization into the equations of the RI-CDD-MP2 gradient with respect to ξ′, there is no difference to directly differentiating the THC-AO-MP2 energy equations. This is the case, because neither X nor Z depend on the perturbation ξ′ in the AO-THC formulation. This independence would not be given if either the MO-THC approach was used or for the more general derivative with respect to ξ, which would necessitate the derivatives of the THC tensors. In other words, by simply inserting the THC factorization into eq 7 the derivative of the THC-ω-RI-CDD-MP2 energy with respect to ξ′ is obtained.

To reduce the computation time needed for forming the expensive R-matrices of the UMP2 gradient, the AO-THC factorization is inserted into eq 7 in its RI-CDD formulation to yield

where X̲ and X̅ are the collocation matrices X transformed into the occupied and virtual Cholesky pseudo-MO basis, respectively. In MP2 it is often advisable to treat Coulomb- and exchange-like contributions separately,⁵⁸ thus the R-matrices are partitioned into R^C, the Coulomb-like contribution, and R^X, the exchange-like contribution.

THC R-Matrices: Coulomb-like Contribution

The Coulomb-like parts R̲^C and R̅^C are given by

and can—in analogy to the THC-MP2 energy—be efficiently computed using sparse linear algebra.⁴¹ In doing so, R̲^C and R̅^C are especially efficient to compute, as their formation only involves BLAS level 3 operations. For an efficient implementation it is important to realize, that for large enough molecules and appropriate ordering of the THC grid points, the collocation matrices X become sparse, while the Z-matrix, being the representation of the long-ranged 1/r operator, will always remain dense. By closer inspection of eq 23, it can be noticed, that the X-matrices corresponding to the MOs of the ket of the decomposed ERI and the Z-matrix are identical for all terms. Therefore, these matrices can be collected in an intermediate D^η given by Algorithm 1, where intermediates A and B are given by

and represent the occupied and virtual pseudodensity in the grid basis, respectively. Additionally, the Λ-factorization of the Z-matrix, i.e., Z ≡ ΛΛ^T, is used to lower the prefactor of this step.^34,41

As forming intermediate D requires a series of dense matrix–matrix-multiplications, this step will require the majority of the computation time for R̲^C and R̅^C. However, Algorithm 1 has to be performed only once per Laplace point and electron spin. The final contribution to the R-matrices can then simply be obtained by a Schur product and two matrix–matrix-multiplications given by Algorithm 2.

The effects of Algorithms 1 and 2 can be best understood by visualizing the underlying tensor contractions in a tensor network diagram, as given in Figure . For an introduction of tensor network diagrams, also in the context of THC, refer to the work by Schutski et al.⁴⁶ Algorithms 1 and 2 are then pieced together with the algorithm for the exchange-like part, detailed in the next section, for the final Algorithm 4 in section 2.4.3.

THC R-Matrices: Exchange-like Contribution

Usually when higher than second-order tensors arise, the associated tensor contractions are either carried out with the tensors reshaped into matrices or batched over the dimensions exceeding matrix dimensionality in so-called tensor slices. For an efficient contraction when iterating over tensor slices it is advisable to make use of an underlying structure to reduce the dimensions of the slices, either by matrix decomposition of the slice or by neglecting noncontributing elements. If whole rows and columns are excluded based on some significance criterion, one arrives at the natural blocking (NB) formalism^57,58 for tensor contractions. NB relies on significance lists, which in general terms describe which pairs, of a general index pair i and j, contribute to a tensor contraction involving these indices. Mathematically speaking these lists are sets and thus can, in set-builder notation, be represented as

where A is a screening matrix involving indices i and j, and ε_NB is the NB screening threshold. To avoid confusion, we use j_i as a shorthand notation for the set of all significant j for a particular index i and {j_i} for the set of all j_i. If two elements in a set {j_i} are identical, it technically becomes a multiset, as elements in sets are only allowed to have a multiplicity of 1. In set terminology, the set {i}_j is the transpose of {j}_i and can analogously determined from A^T as

or directly from {j_i}. Another important quantity is the number of significant pairs N_ij, which is defined as

NB and THC work particularly well together for exchange-like contractions of ERIs, as within the THC formalism the necessary screening matrices can easily be constructed as outlined in the following. The following prototypical exchange-like ERI contraction from the THC-CDD-MP2 energy expression will serve as an example:

Two types of indices are present in the above equation, the orbital indices i and j (occupied space) as well as a and b (virtual space) and the THC auxiliary indices P, Q, R, and S. In the following discussion we use the LS formulation of THC, but the statements also hold for other THC variants. Thus, there are three general types of index pairs: orbital–orbital, orbital–grid, and grid–grid index pairs. The set of significance lists for orbital-grid type pairs is especially easy to construct, since it can be directly derived from the collocation matrices X. For example, the set of all significant grid points P for a given occupied orbital i can be constructed as

Significant orbital–orbital pairs are also easily obtained from the collocation matrices; for example, the set of significant virtual orbitals a for a given occupied orbital i can be built as

The screening criterion from eq 30 can also be interpreted to yield only ia pairs, for which the orbitals have significant overlap and which produce non-negligible charge densities. For the development of low-scaling exchange-type contractions it is important to make use of the exponential coupling between all orbital indices. Orbital i couples to orbital a in an exponentially decaying fashion in the bra of the first ERI in eq 28. Likewise, orbital j also couples to orbital a in an exponentially decaying fashion in the ket of the second ERI. Thus, there is indirect exponential coupling between orbitals i and j via orbital a and the set of significant orbitals i for a given orbital j can be derived from the sets {a_j}, which is identical to {a_i}, and {a_i}:

With the significance lists given by eqs 29, 30, and 31, an asymptotically linear scaling algorithm for the exchange-like energy contribution to MP2 can be devised. The algorithm is detailed in the Supporting Information and close to linear scaling is demonstrated. Here, however, the focus lies on the exchange-like parts of the R-matrices, which are conceptually similar, but different in that two AO indices remain uncontracted. The tensor contractions necessary for forming the R^X parts can again best be understood from the tensor network diagram in Figure .

While Figure only shows the R̅^X part, the R̲^X contribution can be constructed analogously. First, intermediates A or B are formed from the collocation matrices X̲ or X̅, respectively. Second, the remaining collocation matrices, the Z-matrices and intermediates A or B are contracted to the third-order tensor intermediate D, where the symmetry of the tensors is used to reduce the operation count. The idea for an efficient implementation is then as follows: to avoid storing third-order tensors, the contraction is batched over the occupied orbital index common to both R̲ and R̅, and to reduce the cost of the dgemm operations within the loop NB is applied. In this way separate algorithms for R̲^X and R̅^X can be formulated, which are presented in the Supporting Information. By closer inspection, however, it can be seen, that both R^X-matrices share the most expensive to compute intermediate, which incorporates the Z-matrix. Therefore, a joint computation as given by Algorithm 3 is preferred.

First, according to eqs 29, 30, and 31 all necessary significance lists are computed and then the precursor intermediates to R̲^X and R̅^X, that is, the matrices E and F, are accumulated in a loop over the common occupied orbital index j. Lines 12 and 21 represent the bottleneck of the algorithm as their evaluation formally scales as . However, since all involved grid indices P, Q, and R are connected to j through the significance lists {S_j} (identical to {Q_j}) and {R_j} (identical to {P_j}), the size of the involved matrices in NB format is asymptotically constant. Therefore, as N_occ grows linearly with the molecule size, the formation of the R^X-parts can be done in linear scaling time. For this, appropriate thresholds have to be chosen for the formation of the significance lists. Instead of using the same threshold for all pairs, we identify three different types of lists: (1) grid-occupied orbital index pairs, e.g., {S_j}, (2) grid-virtual orbital index pairs, e.g., {b_R}, where, in both cases, the orbital index is coupled to the grid index directly by a collocation matrix, and (3) virtual-occupied orbital index pairs, e.g., {b_j}, which are coupled in real-space over grid points. While the final algorithm only requires the {S_j} and {R_j}/{P_j} lists, the selection of index pairs included in {b_R} and {b_j} is still important, as the {R_j}/{P_j} lists are built from these lists analogously to eq 31.

THC R-Matrices: Final Algorithm

Piecing together Algorithms 1 and 2 for the Coulomb-like part and Algorithm 3 for the exchange-like part, the final algorithm for the formation of the R-matrices is given by Algorithm 4.

Overall, Algorithm 4 has to be executed once per Laplace point and is separated into a precontraction phase and the phase for the actual formation of the contributions to R̲ and R̅. In the precontraction phase intermediates D^η are formed, which represent the most expensive part for the Coulomb-like terms, as the subsequent calls to Algorithm 2 only contribute a Schur product and two dgemm operations. These dgemm calls contribute a negligible overhead compared with line 3 of Algorithm 1, as the dimensions of the matrices involved are reduced. In total, however, Algorithm 3 will dominate the runtime for forming the R-matrices due to its formal scaling.

Computational Details

The above-described THC-ω-RI-CDD-MP2 HFCC code is implemented within our quantum chemistry package FermiONs++.⁵⁹⁻⁶¹ For the THC-based HFCC calculations the hand-optimized grids by Martínez and co-workers³⁸ were used together with the Dunning cc-pVXZ (X ∈ {D, T}) basis sets⁶² and the corresponding auxiliary basis sets. For the phosphorus atoms in the DNA backbone, the fluorine grids were used without loss of accuracy as reported in our work on THC-MP2 energies.⁴¹ All calculations were carried out without the frozen-core approximation. All preceding SCF calculations were converged to an energy difference of 10^–8 H and a FPS−SPF commutator difference of 10^–7 using DIIS acceleration.⁶³ For the gradient calculations seven Laplace points were used in the expansion and the DL-UCPSCF algorithm was converged to an error of 10^–4 for all molecules of the benchmark set in section 4. For all subsequent calculations on larger molecules, a threshold of 10^–3 was used. These settings were shown to yield errors below 1 MHz.^6,30 For the assessment of the accuracy of the THC-ω-RI-CDD-MP2 HFCCs against other methods, the standard orientation was used. For the THC factorization of the ERIs an attenuation strength of 0.1 in the attenuated Coulomb metric was used^54,55 and the same general protocol for screening based on integral partition bounds (IPBs)⁶⁴ and NB, as in our work on THC-MP2 energies, was followed, although adjusted to fit ERIs in the AO basis.⁴¹ All timings are done on an AMD EPYC 7302 (3.30 GHz) CPU node with 256 GB RAM and 1.7 TB of SSD disk space.

Results and Discussion

First, the accuracy of the newly developed THC-ω-RI-CDD-MP2 method for the calculation of isotropic HFCCs is assessed against our reference ω-RI-CDD-MP2 implementation.^6,30 Next, the thresholds necessary for the screening in the expensive exchange-like contribution R^X are optimized on a set of medium-sized organic radicals. Finally, the scaling of the THC-ω-RI-CDD-MP2 method is analyzed and timings are compared to the ω-RI-CDD-MP2 reference for a set of representative radicals.

Accuracy of THC-ω-RI-CDD-MP2 HFCCs

Throughout this publication, our ω-RI-CDD-MP2 implementation for the computation of HFCCs by Vogler et al.,^6,30 which was verified against the RI-MP2 implementation in the ORCA program package,⁶⁵ will serve as reference. The original implementation, however, made use of the SOS approximation and excluded the exchange-like terms. To enable a fair comparison, the exchange-like terms were added analogously to the earlier implemented RI-CDD-based AO-MP2 gradient,⁶ albeit with the attenuated Coulomb metric for the RI integrals. For the comparison, first the accuracy of the presented THC-ω-RI-CDD-MP2 method for HFCCs is assessed. The THC-ω-RI-CDD-MP2 method is benchmarked using a set of 12 organic radicals from a recent study⁷ on the effects of electron correlation, molecular dynamic contributions, and solvation effects on HFCCs. Mean absolute deviations (MAD), root-mean-square deviations (RMSD), and absolute maximum deviations (MAX) are given in Table 1. We note that we used all 12 radicals for the comparison, even though, as Vogler et al.⁷ pointed out, some molecules are spin contaminated. While the latter certainly has an effect on the reliability of the results, when comparing to experiment, it should not influence the comparison of different methods.

Table 1 shows that the mean errors for the THC-ω-RI-CDD-MP2 method are below 1 MHz for both basis sets, while the MAX error corresponds to atoms with high spin density, that is, ¹⁹F in the CF₃ radical for the double-ζ basis set and ¹¹B in the BH₃ radical for the triple-ζ basis set. While for these nuclei the absolute error is larger, the relative error is still below 1%, as due to their high spin density the HFCCs are large in magnitude. The origin of these errors is mainly based on the following two shortcomings: (1) Using the hand-optimized THC grids by Martínez and co-workers³⁸ for AO-THC incurs additional errors over MO-THC as these grids were optimized for fitting MO-based ERIs, for which the orbital space is much more compact compared to ERIs in the AO basis. This phenomenon was observed similarly in our recent work on THC-ω-RI-CDD-MP2 energies.⁴¹ (2) Since grid-based THC uses DFT-like integration grids, THC-based properties are likewise prone to not being rotationally invariant. In DFT the problem is alleviated through larger grids, which THC cannot make use of without forfeiting the reduction in computational cost compared to the respective canonical method. However, we found that these rotational errors are on the order of 0.01−0.05 MHz, depending on the magnitude of the spin density on the respective nucleus. Overall, we consider a mean deviation of less than 1 MHz to be less than the method error of RI-MP2 and certainly accurate enough, when compared against experimental results, where, as shown by Vogler et al.,⁷ other effects like dynamic contributions or solvation effects contribute significantly.

Threshold Optimization

After having established that the presented THC-ω-RI-CDD-MP2 method provides reliable HFCCs, the focus is now on optimizing the time complexity of the underlying algorithm while preserving the accuracy. An obvious point for optimization is the formation of the exchange-like parts R^X, for which Algorithm 3 has quartic scaling if no screening is applied. However, as discussed in section 2.4.2, by applying NB with carefully chosen thresholds, the formation of R^X should be possible with linear time complexity. As outlined in section 2.4.2, different thresholds will be used in the screening process for the different types of significance lists. We associate the thresholds ε_Sj, ε_bR, and ε_bj with the significance lists {S_j}, {b_R}, and {b_j}, respectively. Since ε_Sj directly determines the pairs included in {S_j}, and ε_bR/ε_bj only indirectly determine the pairs in {R_j}, the optimization of these thresholds is simplified by separately optimizing ε_Sj and ε_bR/ε_bj. The thresholds are first optimized for the smaller double-ζ basis set and later transferred to the triple-ζ basis. For this, we chose a set of six medium-sized radical molecules and supramolecular assemblies, which are large enough for the screening to have effect. For further information on this benchmark set see the Supporting Information. Table 2 summarizes the errors and the resulting average number of significant pairs N̅_Sj in {S_j} for the optimization of ε_Sj.

While for threshold values of 10^–6 through 10^–4 the error remains negligible, the screening shows an effect in that N̅_Sj indicates that only roughly a third of the pairs in {S_j} are necessary for this accuracy. The mean errors grow roughly linearly with loosening thresholds and remain sufficiently small for a range of threshold values. For ε_Sj ≥ 10^–3 especially the MAX error deteriorates above 1 MHz, while the MAD remains below 0.1 MHz. For the optimization of ε_bR and ε_bj, combinations of thresholds have to be considered, since they determine, through eq 31, the significant pairs in {R_j}. The results of this optimization are shown in Figure as a heatmap.

From Figure it can be seen, that the MAD is stable through a wide range of threshold values, only significantly worsening when choosing ε_bR > 10^–2, irrespective of the value chosen for ε_bj. The observation that even for looser thresholds, for example, ε_bR = ε_bj = 10^–2, still roughly 80% of the pairs in {R_j} are significant stems from the fact that indices R and j are coupled indirectly over a virtual orbital b. According to eq 31 a pair of indices R and j is only considered insignificant, if they do not share significant overlap with any virtual orbital b. Therefore, N̅_Rj will always be greater than N̅_Sj for any sensibly chosen combination of thresholds.

With the separate optimization of ε_Sj and ε_bR/ε_bj as a starting point, different combinations of these three thresholds were tested. The best trade-off between accuracy and the number of significant pairs appeared to be for the thresholds ε_Sj = 10^–3, ε_bR = 10^–2, and ε_bj = 10^–2, for which the errors are summarized in Table 3.

The chosen combination of thresholds provides good accuracy for both basis sets, with the errors for the triple-ζ basis being somewhat larger. The latter could be improved through a separate optimization of the thresholds for the triple-ζ basis. However, in view of the fact that the mean errors are still below 0.1 MHz, the thresholds optimized for the double-ζ basis set seem to be suitable for the larger basis set as well.

Timings and Scaling

With the optimized screening thresholds at hand, the scaling behavior of the THC-ω-RI-CDD-MP2 method for the computation of HFCCs is analyzed and compared to the previous ω-RI-CDD-MP2 implementation.^6,30 For the assessment of the asymptotic scaling behavior, HFCCs for a series of linear alkyl radicals C_nH_2n+1 are computed. In Figure the timings are shown together with the underlying contributions from the most significant steps for both basis sets.

As is evident from Figure , the overall scaling behavior is governed by the contribution from the exchange-like parts of R (Algorithm 3), while the Coulomb-like terms (Algorithms 1 and 2) and the overhead from obtaining the THC factorization only contribute marginally. For both basis sets, the THC-ω-RI-CDD-MP2 method reaches subquadratic scaling, while the scaling is also partly influenced by the Fock matrix builds in the Z-vector step. The increased cost of the Z-vector step for larger fragment sizes in the case of the double-ζ basis set is also the reason for the overall scaling exponent slightly deteriorating beyond C₁₄₀H₂₈₁ to 1.62. To prevent this unfavorable scaling for the larger triple-ζ basis set, the recommendations by Laqua et al.,⁶⁶ which are default settings in FermiONs++, are followed, and the recently presented seminumerical exchange method (sn-LinK) is used for the exchange part of the Fock matrices. The latter makes the overall scaling for the cc-pVTZ basis set almost entirely be governed by Algorithm 3 (blue bars). Therefore, the scaling reduces to close to linear for the largest fragment size considered. The same is true for the cc-pVDZ basis set, for which Algorithm 3 reaches an apparent asymptotic scaling of 1.3.

To go toward more chemically relevant systems and beyond what was possible with our previous ω-RI-CDD-MP2 implementation, the scaling behavior for spin-labeled adenine–thymine base pair stacks (AT)_n is assessed.

Figure (left) shows the scaling behavior for spin-labeled DNA fragments up to seven repetition units or 5101 basis functions. As expected, the onset for subquadratic scaling is for greater fragment sizes compared to the alkyl radicals. Nonetheless, the scaling exponent decreases to 1.87 for (AT)₆ → (AT)₇. Further reduction with increasing fragment size can be expected based on the growth rate of the number of significant pairs in {S_j} and {R_j}. In Figure (right) the logarithm of the total number of significant pairs is shown for all index pairs relevant for Algorithm 3. The runtime of Algorithm 3 is mainly governed by N_Sj and N_Rj and relies on only a constant number of grid points being significant for a given occupied orbital, see section 2.4.2. If only a constant number of grid points is significant for a given occupied orbital, then N_Sj and N_Rj will grow linearly with increasing molecule size. The latter is demonstrated for N_Sj, that is, the grid point-occupied orbital pair directly coupled by a collocation matrix. N_Rj is inherently greater than N_Sj, since {R_j} is formed from eq 31 with coupling of R and j over a virtual orbital b. The latter is also the reason for the scaling exponent not quite reducing to 1.0 and there being an onset for the close to linear scaling. This also explains why Algorithm 3 reaches subquadratic scaling for the DNA system, but not quite linear scaling. Nonetheless, the THC-ω-RI-CDD-MP2 method also reaches subquadratic scaling for the spin-labeled DNA fragments under consideration and allows for the computation of HFCCs for almost 500 atoms and more than 5000 basis functions.

While the asymptotic scaling exponent of a quantum chemical method is certainly important for the treatment of large (bio)chemical systems, the prefactor oftentimes determines the applicability of a method for a certain problem. In other words, a method can be linear scaling, but have a prefactor so large, that calculations still remain unfeasible. Another aspect which determines feasibility are memory/storage requirements, oftentimes governed by the necessity to store ERIs or amplitude tensors present in electron correlation methods. These aspects are considered in the following for a comparison of the ω-RI-CDD-MP2 method in its all-nuclei variant^6,30 and the presented THC-ω-RI-CDD-MP2 method for a collection of organic radicals. An additional comparison against the selected-nuclei variant of ω-RI-CDD-MP2,³⁰ which was proposed to alleviate some of the shortcomings of the ω-RI-CDD-MP2 method, is given in the Supporting Information. Here, the focus is on typical applications on medium-sized organic radicals. Table 4 summarizes the estimated storage requirements M_est and the wall times t for the radicals, as well as the speedups S relative to the ω-RI-CDD-MP2 implementation. More information on how the storage requirements are estimated is given in the Supporting Information.

THC has a natural advantage over RI-based methods when it comes to storage requirements. In THC-based methods the largest tensors necessary to keep in memory or store on disk are second-order tensors of dimension N_grid × N_grid. For RI-based methods the dimensionality of the factorized representation of the fourth-order ERI tensor increases to three with dimensions N_occ × N_virt × N_aux in the MO basis. Therefore, THC-ω-RI-CDD-MP2 is superior with an order of magnitude less storage requirements, as can be seen from Table 4. Furthermore, the THC-based method considerably outperforms ω-RI-CDD-MP2 in terms of computation time with speedups up to roughly 600 for the double-ζ basis set and 800–1200 for the triple-ζ basis set. The relative speedups are somewhat reduced when comparing the methods for predominantly linear and very sparse systems like linear alkyl radicals, signifying that the ω-RI-CDD-MP2 uses the sparsity well in these model cases. For larger and more globular structures, the computational savings with the THC-ω-RI-CDD-MP2 method are significantly greater, indicating that Algorithm 4 utilizes the sparsity well, even in nonlinear systems. The latter, and the reduced memory requirements make the THC-ω-RI-CDD-MP2 method attractive for the computation of HFCCs of large, globular systems, as they are commonly encountered in proteins and enzymes. Furthermore, due to the greatly reduced computational cost, the method is applicable in double-hybrid functionals, where usually the MP2 part is the computational bottleneck, while also allowing for sampling of multiple points on the PES.

Conclusions and Outlook

In this work, we presented the THC-ω-RI-CDD-MP2 method for the efficient and accurate computation of isotropic HFCCs for large organic radicals. As usual for MP2 methods, the exchange-like terms, here R^X, govern the scaling and runtime of the method. This issue was addressed through screening based on the THC collocation matrices X in combination with natural blocking for the tensor contractions. An asymptotically linear scaling recipe for the contraction of exchange-like terms in THC format is provided. This recipe is applied to the R^X terms, reducing the formal quartic scaling to effectively subquadratic, as shown for linear alkyl radicals and spin-labeled DNA strands. The THC-ω-RI-CDD-MP2 method furthermore highlights the attractiveness of THC-based methods for derivatives of electron correlation methods. Derivative calculations usually involve more types of ERIs, for example, half-transformed integrals, integrals with mixed spin in openshell calculations, or integrals contracted with perturbed densities. This generally increases storage requirements compared to single point calculations. Furthermore, the computational cost of the method is strongly increased due to additional ERI contractions. Both challenges are overcome by using THC-factorized ERIs, for which only second-order tensors have to be stored and integral transformations and contractions can be easily reduced to simple dgemm operations. The advantages of THC-based gradient methods are demonstrated for HFCC calculations on a range of medium-sized organic radicals and spin-labeled DNA strands with more than 5000 basis functions.

For future applications, the availability of THC grids has to be improved, as ideally EPR-specific basis sets,⁶⁷ and corresponding THC grids, should be used for the calculations presented. The grids used throughout this publications were hand optimized³⁸ and are only viable for the cc-pVXZ (X ∈{D,T}) basis sets. Furthermore, these grids were optimized based on MO-THC-MP2 and not for AO-THC, as used throughout this work. This incurs additional errors due to the larger fitting space. Different techniques for the on-the-fly generation of THC grids have been proposed, either based on PCD of the THC metric in a larger parent grid basis,⁴⁰ or based on centroidal Voronoi tesselation (CVT),³⁶ and could be applied for the calculation of THC-MP2 HFCCs in future work.

Finally, the developed THC-ω-RI-CDD-MP2 method can easily be used for the MP2 part of double-hybrid functionals and, once appropriate grids are available, enable accurate HFCC predictions for large molecules. Furthermore, due to the reduced computational complexity and memory requirements, THC-ω-RI-CDD-MP2, in conjunction with an appropriate double-hybrid functional, is an attractive candidate for QM/MM HFCC calculations.

References

1. C. Møller, M. S. Plesset. Note on an approximation treatment for many-electron systems.. Phys. Rev., 1934. [DOI]
2. D. Cremer. Møller-Plesset perturbation theory: From small molecule methods to methods for thousands of atoms.. Wiley Interdiscip. Rev. Comput. Mol. Sci., 2011. [DOI]
3. G. D. Bent. Many-body perturbation theory electronic structure calculations for the methoxy radical. II. Hyperfine coupling coefficients.. J. Chem. Phys., 1994. [DOI]
4. S. Kossmann, F. Neese. Correlated ab initio spin densities for larger molecules: Orbital-optimized spin-component-scaled MP2 method.. J. Phys. Chem. A, 2010. [DOI | PubMed]
5. B. Sandhoefer, S. Kossmann, F. Neese. Derivation and assessment of relativistic hyperfine-coupling tensors on the basis of orbital-optimized second-order Møller-Plesset perturbation theory and the second-order Douglas-Kroll-Hess transformation.. J. Chem. Phys., 2013. [DOI | PubMed]
6. S. Vogler, M. Ludwig, M. Maurer, C. Ochsenfeld. Low-scaling first-order properties within second-order Møller-Plesset perturbation theory using Cholesky decomposed density matrices.. J. Chem. Phys., 2017. [DOI | PubMed]
7. S. Vogler, J. C. B. Dietschreit, L. D. M. Peters, C. Ochsenfeld. Important components for accurate hyperfine coupling constants: electron correlation, dynamic contributions, and solvation effects.. Mol. Phys., 2020. [DOI]
8. J. Gauss. Calculation of NMR chemical shifts at second-order many-body perturbation theory using gauge-including atomic orbitals.. Chem. Phys. Lett., 1992. [DOI]
9. J. Gauss. Effects of electron correlation in the calculation of nuclear magnetic resonance chemical shifts.. J. Chem. Phys., 1993. [DOI]
10. D. Flaig, M. Maurer, M. Hanni, K. Braunger, L. Kick, M. Thubauville, C. Ochsenfeld. Benchmarking Hydrogen and Carbon NMR Chemical Shifts at HF, DFT, and MP2 Levels.. J. Chem. Theory Comput., 2014. [DOI | PubMed]
11. M. Maurer, C. Ochsenfeld. A linear- and sublinear-scaling method for calculating NMR shieldings in atomic orbital-based second-order Møller-Plesset perturbation theory.. J. Chem. Phys., 2013. [DOI | PubMed]
12. M. Glasbrenner, S. Vogler, C. Ochsenfeld. Efficient low-scaling computation of NMR shieldings at the second-order Møller-Plesset perturbation theory level with Cholesky-decomposed densities and an attenuated Coulomb metric.. J. Chem. Phys., 2021. [DOI | PubMed]
13. S. Kossmann, B. Kirchner, F. Neese. Performance of modern density functional theory for the prediction of hyperfine structure: Meta-GGA and double hybrid functionals.. Mol. Phys., 2007. [DOI]
14. S. Sæbø, P. Pulay. Local configuration interaction: An efficient approach for larger molecules.. Chem. Phys. Lett., 1985. [DOI]
15. S. Sæbø, P. Pulay. Local treatment of electron correlation.. Annu. Rev. Phys. Chem., 1993. [DOI]
16. A. El Azhary, G. Rauhut, P. Pulay, H.-J. Werner. Analytical energy gradients for local second-order Møller-Plesset perturbation theory.. J. Chem. Phys., 1998. [DOI]
17. M. Schütz, H. J. Werner, R. Lindh, F. R. Manby. Analytical energy gradients for local second-order Møller-Plesset perturbation theory using density fitting approximations.. J. Chem. Phys., 2004. [DOI | PubMed]
18. P. Pinski, F. Neese. Communication: Exact analytical derivatives for the domain-based local pair natural orbital MP2 method (DLPNO-MP2).. J. Chem. Phys., 2018. [DOI | PubMed]
19. P. Pinski, F. Neese. Analytical gradient for the domain-based local pair natural orbital second order Møller-Plesset perturbation theory method (DLPNO-MP2).. J. Chem. Phys., 2019. [DOI | PubMed]
20. G. L. Stoychev, A. A. Auer, J. Gauss, F. Neese. DLPNO-MP2 second derivatives for the computation of polarizabilities and NMR shieldings.. J. Chem. Phys., 2021. [DOI | PubMed]
21. K. Kristensen, I. M. Høyvik, B. Jansik, P. Jørgensen, T. Kjærgaard, S. Reine, J. Jakowski. MP2 energy and density for large molecular systems with internal error control using the Divide-Expand-Consolidate scheme.. Phys. Chem. Chem. Phys., 2012. [DOI | PubMed]
22. P. Baudin, P. Ettenhuber, S. Reine, K. Kristensen, T. Kjærgaard. Efficient linear-scaling second-order Møller-Plesset perturbation theory: The divide-expand-consolidate RI-MP2 model.. J. Chem. Phys., 2016. [DOI | PubMed]
23. K. Kristensen, P. Jorgensen, B. Jansík, T. Kjærgaard, S. Reine. Molecular gradient for second-order Møller-Plesset perturbation theory using the divide-expand-consolidate (DEC) scheme.. J. Chem. Phys., 2012. [DOI | PubMed]
24. D. Bykov, K. Kristensen, T. Kjærgaard. The molecular gradient using the divide-expand-consolidate resolution of the identity second-order Møller-Plesset perturbation theory: The DEC-RI-MP2 gradient.. J. Chem. Phys., 2016. [DOI | PubMed]
25. S. Schweizer, B. Doser, C. Ochsenfeld. An atomic orbital-based reformulation of energy gradients in second-order Møller-Plesset perturbation theory.. J. Chem. Phys., 2008. [DOI | PubMed]
26. J. L. Whitten. Coulombic potential energy integrals and approximations.. J. Chem. Phys., 1973. [DOI]
27. B. I. Dunlap, J. W. Connolly, J. R. Sabin. On some approximations in applications of Xα theory.. J. Chem. Phys., 1979. [DOI]
28. O. Vahtras, J. Almlöf, M. W. Feyereisen. Integral approximations for LCAO-SCF calculations.. Chem. Phys. Lett., 1993. [DOI]
29. F. Weigend, M. Häser, H. Patzelt, R. Ahlrichs. RI-MP2: Optimized auxiliary basis sets and demonstration of efficiency.. Chem. Phys. Lett., 1998. [DOI]
30. S. Vogler, G. Savasci, M. Ludwig, C. Ochsenfeld. Selected-Nuclei Method for the Computation of Hyperfine Coupling Constants within Second-Order Møller-Plesset Perturbation Theory.. J. Chem. Theory Comput., 2018. [DOI | PubMed]
31. R. M. Parrish, E. G. Hohenstein, T. J. Martínez, C. D. Sherrill. Tensor hypercontraction. II. Least-squares renormalization.. J. Chem. Phys., 2012. [DOI | PubMed]
32. R. M. Parrish, E. G. Hohenstein, N. F. Schunck, C. D. Sherrill, T. J. Martínez. Exact tensor hypercontraction: A universal technique for the resolution of matrix elements of local finite-range N-body potentials in many-body quantum problems.. Phys. Rev. Lett., 2013. [DOI]
33. C. Song, T. J. Martínez. Atomic orbital-based SOS-MP2 with tensor hypercontraction. I. GPU-based tensor construction and exploiting sparsity.. J. Chem. Phys., 2016. [DOI | PubMed]
34. C. Song, T. J. Martínez. Atomic orbital-based SOS-MP2 with tensor hypercontraction. II. Local tensor hypercontraction.. J. Chem. Phys., 2017. [DOI | PubMed]
35. D. A. Matthews. A critical analysis of least-squares tensor hypercontraction applied to MP3.. J. Chem. Phys., 2021. [DOI | PubMed]
36. J. Lee, L. Lin, M. Head-Gordon. Systematically Improvable Tensor Hypercontraction: Interpolative Separable Density-Fitting for Molecules Applied to Exact Exchange, Second- and Third-Order Møller-Plesset Perturbation Theory.. J. Chem. Theory Comput., 2020. [DOI | PubMed]
37. E. G. Hohenstein, R. M. Parrish, T. J. Martínez. Tensor hypercontraction density fitting. I. Quartic scaling second- and third-order Møller-Plesset perturbation theory.. J. Chem. Phys., 2012. [DOI | PubMed]
38. S. I. Kokkila Schumacher, E. G. Hohenstein, R. M. Parrish, L. P. Wang, T. J. Martínez. Tensor Hypercontraction Second-Order Møller-Plesset Perturbation Theory: Grid Optimization and Reaction Energies.. J. Chem. Theory Comput., 2015. [DOI | PubMed]
39. C. Song, T. J. Martínez. Analytical gradients for tensor hyper-contracted MP2 and SOS-MP2 on graphical processing units.. J. Chem. Phys., 2017. [DOI | PubMed]
40. D. A. Matthews. Improved Grid Optimization and Fitting in Least Squares Tensor Hypercontraction.. J. Chem. Theory Comput., 2020. [DOI | PubMed]
41. F. H. Bangerter, M. Glasbrenner, C. Ochsenfeld. Low-Scaling Tensor Hypercontraction in the Cholesky Molecular Orbital Basis Applied to Second-Order Møller-Plesset Perturbation Theory.. J. Chem. Theory Comput., 2021. [DOI | PubMed]
42. N. Shenvi, H. Van Aggelen, Y. Yang, W. Yang. Tensor hypercontracted ppRPA: Reducing the cost of the particle-particle random phase approximation from O(r6) to O(r4).. J. Chem. Phys., 2014. [DOI | PubMed]
43. C. Song, T. J. Martínez. Reduced scaling CASPT2 using supporting subspaces and tensor hyper-contraction.. J. Chem. Phys., 2018. [DOI | PubMed]
44. E. G. Hohenstein, S. I. Kokkila, R. M. Parrish, T. J. Martínez. Quartic scaling second-order approximate coupled cluster singles and doubles via tensor hypercontraction: THC-CC2.. J. Chem. Phys., 2013. [DOI | PubMed]
45. R. M. Parrish, C. D. Sherrill, E. G. Hohenstein, S. I. Kokkila, T. J. Martínez. Communication: Acceleration of coupled cluster singles and doubles via orbital-weighted least-squares tensor hypercontraction.. J. Chem. Phys., 2014. [DOI | PubMed]
46. R. Schutski, J. Zhao, T. M. Henderson, G. E. Scuseria. Tensor-structured coupled cluster theory.. J. Chem. Phys., 2017. [DOI | PubMed]
47. E. G. Hohenstein, S. I. Kokkila, R. M. Parrish, T. J. Martínez. Tensor hypercontraction equation-of-motion second-order approximate coupled cluster: Electronic excitation energies in O(N4) time.. J. Phys. Chem. B, 2013. [DOI | PubMed]
48. C. Song, T. J. Martínez, J. B. Neaton. A diagrammatic approach for automatically deriving analytical gradients of tensor hyper-contracted electronic structure methods.. J. Chem. Phys., 2021. [DOI | PubMed]
49. J. Almlöf. Elimination of energy denominators in Møller-Plesset perturbation theory by a Laplace transform approach.. Chem. Phys. Lett., 1991. [DOI]
50. M. Häser, J. Almlöf. Laplace transform techniques in Møller-Plesset perturbation theory.. J. Chem. Phys., 1992. [DOI]
51. M. Häser. Møller-Plesset (MP2) perturbation theory for large molecules.. Theor. Chim. Acta, 1993. [DOI]
52. N. C. Handy, H. F. Schaefer. On the evaluation of analytic energy derivatives for correlated wave functions.. J. Chem. Phys., 1984. [DOI]
53. M. Beer, C. Ochsenfeld. Efficient linear-scaling calculation of response properties: Density matrix-based Laplace-transformed coupled-perturbed self-consistent field theory.. J. Chem. Phys., 2008. [DOI | PubMed]
54. S. Reine, E. Tellgren, A. Krapp, T. Kjærgaard, T. Helgaker, B. Jansik, S. Høst, P. Salek. Variational and robust density fitting of four-center two-electron integrals in local metrics.. J. Chem. Phys., 2008. [DOI | PubMed]
55. A. Luenser, H. F. Schurkus, C. Ochsenfeld. Vanishing-Overhead Linear-Scaling Random Phase Approximation by Cholesky Decomposition and an Attenuated Coulomb-Metric.. J. Chem. Theory Comput., 2017. [DOI | PubMed]
56. Y. Jung, R. C. Lochan, A. D. Dutoi, M. Head-Gordon. Scaled opposite-spin second order Møller-Plesset correlation energy: An economical electronic structure method.. J. Chem. Phys., 2004. [DOI | PubMed]
57. Y. Jung, Y. Shao, M. Head-Gordon. Fast Evaluation of Scaled Opposite Spin Second-Order Møller-Plesset Correlation Energies Using Auxiliary Basis Expansions and Exploiting Sparsity.. J. Comput. Chem., 2007. [DOI | PubMed]
58. M. Glasbrenner, D. Graf, C. Ochsenfeld. Efficient Reduced-Scaling Second-Order Møller-Plesset Perturbation Theory with Cholesky-Decomposed Densities and an Attenuated Coulomb Metric.. J. Chem. Theory Comput., 2020. [DOI | PubMed]
59. J. Kussmann, C. Ochsenfeld. Pre-selective screening for matrix elements in linear-scaling exact exchange calculations.. J. Chem. Phys., 2013. [DOI | PubMed]
60. J. Kussmann, C. Ochsenfeld. Preselective screening for linear-scaling exact exchange-gradient calculations for graphics processing units and general strong-scaling massively parallel calculations.. J. Chem. Theory Comput., 2015. [DOI | PubMed]
61. J. Kussmann, C. Ochsenfeld. Hybrid CPU/GPU Integral Engine for Strong-Scaling Ab Initio Methods.. J. Chem. Theory Comput., 2017. [DOI | PubMed]
62. T. H. Dunning. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen.. J. Chem. Phys., 1989. [DOI]
63. P. Pulay. Convergence acceleration of iterative sequences. the case of scf iteration.. Chem. Phys. Lett., 1980. [DOI]
64. T. H. Thompson, C. Ochsenfeld. Integral partition bounds for fast and effective screening of general one-, two-, and many-electron integrals.. J. Chem. Phys., 2019. [DOI | PubMed]
65. F. Neese. The ORCA program system.. Wiley Interdiscip. Rev. Comput. Mol. Sci., 2012. [DOI]
66. H. Laqua, T. H. Thompson, J. Kussmann, C. Ochsenfeld. Highly Efficient, Linear-Scaling Seminumerical Exact-Exchange Method for Graphic Processing Units.. J. Chem. Theory Comput., 2020. [DOI | PubMed]
67. P. Jakobsen, F. Jensen. Probing basis set requirements for calculating hyperfine coupling constants.. J. Chem. Phys., 2019. [DOI | PubMed]