A method for the computation of nuclear magnetic resonance (NMR) shieldings with second-order Møller–Plesset perturbation theory (MP2) is presented which allows to efficiently compute the entire set of shieldings for a given molecular structure. The equations are derived using Laplace-transformed atomic orbital second-order Møller–Plesset perturbation theory as a starting point. The Z-vector approach is employed for minimizing the number of coupled-perturbed self-consistent-field equations that need to be solved. In addition, the method uses the resolution-of-the-identity approximation with an attenuated Coulomb metric and Cholesky decomposition of pseudo-density matrices. The sparsity in the three-center integrals is exploited with sparse linear algebra approaches, leading to reduced computational cost and memory demands. Test calculations show that the deviations from NMR shifts obtained with canonical MP2 are small if appropriate thresholds are used. The performance of the method is illustrated in calculations on DNA strands and on glycine chains with up to 283 atoms and 2864 basis functions.

Nuclear magnetic resonance (NMR) is a widely used spectroscopic method in the field of chemistry. NMR spectra are highly sensitive to the molecular geometry and contain a wealth of structural information. However, especially for large molecules, their interpretation can be difficult, and it can be challenging to unambiguously assign a molecular structure to a given spectrum. In such situations, the comparison with theoretically computed spectra can be very helpful. For this reason, much work has been done on the development of quantum-chemical methods for the simulation of NMR shieldings (for reviews, see, e.g., Refs. 1–3).

Methods for computing NMR shieldings with Hartree–Fock (HF) theory4–8 or density functional theory (DFT)9–11 can often provide useful shifts at moderate computational cost. For these theories, linear and sublinear scaling implementations have been developed which allow to apply them to large molecules with 1000 atoms or more.12–14 If higher accuracy is desired, wave-function based correlation methods can be used. This includes second-order Møller–Plesset perturbation theory (MP2),15,16 the multiconfigurational self-consistent field method (MCSCF),17,18 or approaches based on coupled-cluster theory.19–21 Among the mentioned correlated methods, MP2 is particularly attractive because it has lower formal scaling and computational cost than high-level coupled cluster methods but has been shown to yield typically NMR shifts of significantly higher quality than Hartree–Fock theory or DFT.16,22 Therefore, a significant amount of work has been done on MP2-based methods for computing NMR shieldings. In 1992, Gauss presented the first method for computing NMR shifts with MP215,16 and a Z-vector approach.23,24 In later years, a more efficient integral-direct implementation by Kollwitz and Gauss25 and implementations using local-correlation approximations26–28 were developed. Maurer and Ochsenfeld29 demonstrated that a linear and sublinear scaling computation of NMR shieldings is possible with atomic-orbital MP2 (AO-MP2)30–35 and integral screening approaches. Apart from pure MP2 theory, closely related methods, such as double-hybrid density functional theory36,37 or spin-component-scaled38,39 and scaled-opposite-spin MP2,39,40 showed good accuracy in benchmark studies and proved to be promising for efficient calculations on large molecules considering that the algorithmic developments for MP2 can easily be transferred to these methods.

In this work, we present a new method for efficiently computing NMR shieldings at the MP2 level of theory. Like the related approach from the study by Maurer and Ochsenfeld,29 the method presented in this work is based on Laplace-transformed AO-MP2 but differs in several aspects. First of all, our new method employs an all-nuclei formulation, which allows us to efficiently compute the entire set of NMR shieldings for all nuclei in a given molecule. In contrast, Maurer and Ochsenfeld29 used a selected-nuclei formulation, which was first introduced by Beer et al.14 at the SCF level of theory. The selected-nuclei formulation allows to exploit locality and accelerate the computation of the shieldings of a small subset of all nuclei. However, if all shieldings need to be computed, this formulation is inferior due to the large number of coupled-perturbed self-consistent field (CPSCF) equations that need to be solved. The working equations differ significantly between the all-nuclei and selected nuclei ansatz because certain terms involving derivatives of gauge-including atomic orbitals (GIAOs)4,41,42 are present in one formulation and absent in the other one. GIAOs are used in our method in order to ensure gauge-origin independence. Compared to the method reported by Maurer and Ochsenfeld,29 our new implementation avoids four-center electron repulsion integrals (ERIs) with the help of the resolution-of-the-identity (RI)43,44 approximation. Like in our previous work on the efficient computation of MP2 energies,45 an attenuated Coulomb metric46–48 is used in order to further increase the available sparsity in the three-center integrals. In this context, we would like to mention an efficient MP2 shielding method using a Cholesky decomposition (CD) of ERIs that was very recently introduced by Burger et al.49 As RI, the CD approach allows to avoid four-center integrals. Compared to RI, it has certain advantages, such as rigorous control of the numerical error of the decomposition and independence of pre-defined auxiliary basis sets.50,51 However, the Cholesky factors have lower sparsity than three-center RI integrals with an attenuated Coulomb metric for large systems.

Here we use a pivoted Cholesky decomposition52–54 of pseudo-density matrices which allows us to utilize the reduced rank of these matrices (especially for the occupied pseudo-densities) and further lowers computational cost and memory demands. The sparsity of the three-center integrals is exploited using sparse linear algebra approaches. The transformations of the three-center integrals in the atomic orbital basis with the Cholesky molecular orbitals (Cholesky-MOs) and with pseudo-density and perturbed pseudo-density matrices are accelerated using block-sparse matrix multiplications. For the subsequent steps involving three-center integrals, we use the natural blocking format,45,47,55 which allows us to reduce the computational effort for contractions of the three-center integrals and the I/O overhead. We demonstrate in test calculations that the developed method provides accurate NMR shieldings if the same thresholds are used that also proved to be appropriate for achieving chemically accurate energies in our recent work on MP2 correlation energies.45 In calculations on glycine chains and DNA strands, we analyze the efficiency of the method, the scaling, and the requirements on memory and disk space. We also compare calculations with standard MP2 to shielding calculations with scaled-opposite spin MP2 (SOS-MP2)39,40 and show that the SOS-approximation leads to further substantial gains in efficiency.

In Sec. II A, a brief review of the AO-MP2 method for energies is provided, which provides the starting point for the derivation. Section II B contains a discussion of AO-MP2 gradients with emphasis on the gradient with respect to the nuclear magnetic moment. In Sec. II C, the mixed second derivative of the energy with respect to nuclear magnetic moment and magnetic field is derived, which allows to compute NMR shieldings. In Sec. II D, a nested Z-vector approach is described, which allows to circumvent the explicit computation of the nuclear magnetic moment derivative of the density matrix and the second derivative of the density matrix and thus minimizes the number of CPSCF equations that need to be solved. Then, in Sec. II E, we discuss strategies for accelerating the time-dominating steps of the computation.

In general, the elements of the NMR shielding tensor σA of a nucleus A can be computed by taking a mixed second derivative of the electronic energy E. One derivative needs to be taken with respect to one of the components of the magnetic field Bs with sx,y,z, while the other derivative is taken with respect to a component of the nuclear magnetic moment vector mrA,

σrsA=2mrABsEmA=0,B=0.
(1)

Different quantum-chemical methods can be used for computing E; the method presented in this work is based on atomic-orbital MP2 (AO-MP2).30,32,34,35 Starting from the closed-shell MP2 energy expression in terms of canonical, doubly occupied Hartree–Fock orbitals,

EMP2=ijabia|jb2ia|jbib|jaϵa+ϵbϵiϵj,
(2)

AO-MP2 can be derived by expressing the orbital energy denominator from Eq. (2) using a Laplace transformation, which can be approximated using numerical integration30 as follows:

1ϵa+ϵbϵiϵj=0dtetϵa+ϵbϵiϵjαωαetαϵa+ϵbϵiϵj,
(3)

where tα and ωα are the exponent and weight of Laplace point α, respectively. After inserting the Laplace-transformed denominator into Eq. (2), insertion of the basis set expansion of the MOs, and subsequent summation over all MO-indices, the energy expression of AO-MP2 is obtained as follows:

EMP2AO=αμνλσμ̲ν̄|λ̲σ̄2μν|λσμσ|λν,
(4)

where

μ̲ν̄|λ̲σ̄=μνλσP̲μμP̄ννP̲λλP̄σσμν|λσ
(5)

and

P̲μμ=ωα4iCμie+ϵitαCμi,
(6)
P̄νν=ωα4aCνaeϵatαCνa.
(7)

The matrices P̲ and P̄ defined in Eqs. (6) and (7) are called the occupied and virtual pseudo-densities, respectively. Note that the pseudo-density matrices depend on the Laplace point α; for simplicity, we omit this dependence in all formulas shown below.

In Ref. 56, Schweizer et al. derived the gradient of AO-MP2. One possible expression for the gradient of AO-MP2 with respect to a general perturbation ξ reads as follows:

ξEMP2AO=2αTr[R̲P̄ξ]+Tr[R̄P̲ξ]2αμνλσμ̲ν̄|λ̲σ̄2μν|λσξμν|λσξ,
(8)

where

R̄μμ=νλσμν̄|λ̲σ̄2μν|λσμσ|λν,
(9)
R̲νν=μλσμ̲ν|λ̲σ̄2μν|λσμσ|λν.
(10)

Due to the pseudo-densities appearing in Eqs. (9) and (10), the matrices R̲ and R̄ are also dependent on the Laplace point. For the derivation of the working equations of our new MP2-NMR method, we will take the gradient with respect to the nuclear magnetic moment mA as a starting point,

mEMP2AO=2αTr[R̲P̄m]+Tr[R̄P̲m].
(11)

Here and in the following, mA shall be abbreviated as m in order to simplify the notation. Note that for this particular perturbation, the last term from Eq. (8) is zero because the two-electron integrals do not depend on m. Equation (11) requires perturbed pseudo-density matrices. Explicit expressions for the unperturbed pseudo-density matrices in terms of MO coefficients and orbital energies are shown in Eqs. (6) and (7); these matrices can also be computed only from quantities in the AO basis33,56,57 as follows:

P̲=ωα4etPFP,
(12)
P̄=ωα4etPvirtFPvirt,
(13)

where P=iCμi*Cνi is the HF density matrix and Pvirt=aCμa*Cνa is the virtual density matrix. After inserting the expressions from Eqs. (12) and (13) into Eq. (11) and differentiating using the product rule, one obtains

mEMP2AO=2αωα4Tr[R̄etPFmP]+Tr[R̄etPFPm]+Tr[R̲etPvirtFmPvirt]+Tr[R̲etPvirtFPvirtm].
(14)

By applying cyclic permutations under the trace, Eq. (14) becomes

mEMP2AO=2αωα4Tr[R̄etPFPm]+Tr[PR̄etPFm]+Tr[R̲etPvirtFPvirtm]+Tr[PvirtR̲etPvirtFm].
(15)

In general, the perturbed virtual density matrix can be computed as follows:56 

Pvirtξ=Pξ(P+Pvirt)SξS1.
(16)

For the nuclear magnetic moment perturbation, the derivative of S is zero, and Eq. (16) therefore becomes Pvirtm=Pm. The terms involving the perturbed density matrices can be efficiently treated with a Z-vector approach, which will be described in detail in Sec. II D. Here, we will focus on the terms containing perturbed matrix exponentials, which are of the general form Tr[BeAξ]. As shown in  Appendix A, such a term can be rearranged as follows:

Tr[BeAξ]=Tr[YAξ],
(17)

where

Y=k=1i=0k11k!AiBAk1i.
(18)

Y can be computed efficiently by a recursive scheme, which is also shown in  Appendix A. In order to rearrange the terms containing perturbed matrix exponentials in Eq. (15), we define two intermediates Y̲ and Ȳ. Ȳ is computed by evaluating Eq. (18) with A=tαPF and B=PR̄. Y̲ is calculated using the same formula and setting A=tαPvirtF and B=PvirtR̲. Using Y̲ and Ȳ, the intermediates Y̲1, Y̲2, Ȳ1, and Ȳ2 defined in Ref. 56 can be computed by a small number of additional matrix multiplications as follows:

Y̲1=tαFY̲,
(19)
Y̲2=tαY̲Pvirt,
(20)
Ȳ1=tαFȲ,
(21)
Ȳ2=tαȲP.
(22)

Instead of evaluating four different recursion formulas for Y̲1, Y̲2, Ȳ1, and Ȳ2 as in Ref. 56, only two recursions for Y̲ and Ȳ need to be evaluated if the new recursion formulas presented in this work are used. This reduces the number of necessary matrix multiplications by roughly a factor of two. Note that also all “Y”-matrices from Eqs. (19)(22) implicitly depend on the Laplace point.

Using the intermediates defined so far, we are able to write the formula for the gradient with respect to m in the following way:

mEMP2AO=2αωα4Tr[F(α)hm]+Tr[P(α)Pm],
(23)
F(α)=Y̲2+Ȳ2,
(24)
P(α)=Ȳ1Y̲1+GF(α)+R̄etαPFR̲etαPvirtF.
(25)

GF(α) is computed analogously to the two-electron part of the Fock matrix with F(α) taking the role of the density matrix,

GμνF(α)=λσ2μν|λσμσ|λνFλσ(α).
(26)

Noticing that only F(α) and P(α) depend on the Laplace point in Eq. (23), we define

Fαωα4F(α),
(27)
Pαωα4P(α).
(28)

Using F and P, Eq. (23) becomes

mEMP2AO=2Tr[Fhm]+Tr[PPm].
(29)

As shown in Eq. (1), NMR shieldings can be computed as the mixed second derivative of the energy with respect to the nuclear magnetic moment and the magnetic field. The expression from Eq. (29) thus needs to be differentiated with respect to the magnetic field B as follows:

2BmEMP2AO=2Tr[FhBm]+Tr[FBhm]+Tr[PPBm]+Tr[PBPm].
(30)

Equation (30) contains the derivative of F and P, which can be computed as follows:

FB=αωα4(F(α))B=αωα4{Y̲2B+Ȳ2B},
(31)
PB=αωα4Ȳ1BY̲1B+R̄BetαPF+R̄etαPFBR̲BetαPvirtFR̲etαPvirtFB+GBF+GFB.
(32)

Equations (31) and (32) require the magnetic field derivative of the “Y”-matrices from Eqs. (19)(22),

Ȳ1B=tαFBȲ+FȲB,
(33)
Ȳ2B=tαȲBP+ȲPB,
(34)
Y̲1B=tαFBY̲+FY̲B,
(35)
Y̲2B=tαY̲BPvirt+Y̲PvirtB.
(36)

Here, Y̲B and ȲB can be computed recursively. The corresponding recursion formulas are derived in  Appendix B by differentiating the recursion formulas for Y̲ and Ȳ. For all recursion formulas, asymptotic linear scaling can be achieved by employing block-sparse matrices. The matrix GBF is computed in a similar way to GF in Eq. (26) with the difference that the magnetic field derivatives of the ERIs are used instead of the standard unperturbed ERIs,

GμνBF=λσ2μν|λσ(B)μσ|λν(B)Fλσ.
(37)

For an efficient calculation of GBF and G[FB] in Eq. (32) and of GF in Eq. (26), any approach for efficient Fock matrix construction in Hartree–Fock theory can be used, such as the Continuous Fast Multipole Method (CFMM)58,59 or the RI-J method60 for the Coulomb matrix and LinK screening61,62 or semi-numerical exchange approaches63 for the exchange matrix.

Equations (31) and (32) also require the magnetic field derivative of R̲ and R̄, which can be calculated as follows:

R̲ννB=μλσμ̲Bν|λ̲σ̄+μ̲ν|λ̲Bσ̄+μ̲ν|λ̲σ̄B×2μν|λσμσ|λν+μ̲ν|λ̲σ̄(B)2μν|λσμσ|λν+μ̲ν|λ̲σ̄2μν|λσ(B)μσ|λν(B),
(38)
R̄μμB=νλσμν̄B|λ̲σ̄+μν̄|λ̲Bσ̄+μν̄|λ̲σ̄B×2μν|λσμσ|λν+μν̄|λ̲σ̄(B)2μν|λσμσ|λν+μν̄|λ̲σ̄2μν|λσ(B)μσ|λν(B).
(39)

Both Eqs. (38) and (39) contain two types of terms: terms involving magnetic field derivatives of ERIs and terms in which the unperturbed ERIs are transformed both with perturbed and unperturbed pseudo-densities. The efficient computation of these terms is discussed in Sec. II E. The notation with (B) as in μν̄|λ̲σ̄(B) means that only the integrals are differentiated, not the pseudo-density or density matrices.

While Eq. (30) would allow us to compute the NMR shieldings of all nuclei in a molecule, these equations still contain Pm and PBm. The explicit computation of these matrices should be avoided in order to reduce the number of CPSCF equations that need to be solved. This can be achieved using a nested Z-vector approach similar to the one presented by Maurer and Ochsenfeld.29 In this section, we will present the theory behind this approach in a general way, which can be applied to any formulation of CPSCF. In our implementation, we employ the density matrix based CPSCF (D-CPSCF) method introduced by Ochsenfeld and Head-Gordon;64 detailed working equations for this formulation are derived in  Appendix D.

In general, a CPSCF equation for Pm has the following structure:

APm=bm,
(40)

where bm is the right-hand side of the CPSCF equation and A is the Hartree–Fock Hessian matrix. Formally, this equation can be solved by multiplying with the inverse of the Hessian from the left as follows:

Pm=A1bm.
(41)

Usually, an explicit inversion of A is circumvented in favor of an iterative solution. In order to compute all matrices Pm for a given molecule, 3 × Nat CPSCF equations would need to be solved where Nat denotes the number of atoms. In a term of the form Tr[XPm], however, this can avoided by applying a Z-vector approach,23,24

Tr[XPm]=Tr[XA1bm]=Tr[ZXbm],
(42)

where ZX is obtained by solving a CPSCF-like equation with X as the right-hand side as follows:

AZX=X.
(43)

The advantage here is that only a single CPSCF equation needs to be solved. The outlined strategy could be directly applied to the term Tr[PBPm]; however, as we will discuss in the following, a contribution from the term Tr[PPBm] is first added to PB in our implementation.

The term Tr[PPBm] contains the mixed second derivative of the density matrix. An equation for explicitly computing PBm can be derived from the CPSCF equation APB = bB by differentiating it with respect to m as follows:

AmPB+APBm=bBm.
(44)

PBm can thus be obtained by solving a CPSCF equation with (bBmAmPB) as a right-hand side as follows:

PBm=A1bBmAmPB.
(45)

Next, we insert Eq. (45) into the term TrPPBm as follows:

TrPPBm=TrPA1bBmAmPB=TrZPbBmAmPB.
(46)

ZP is computed by solving a CPSCF equation with P as the right-hand side as follows:

AZP=P.
(47)

All terms in TrZPbBmAm[PB] that depend on Pm need to be rearranged in the form Tr[…Pm] such that another Z-vector approach can be applied. A detailed derivation is provided in  Appendix D. It is shown that the sum of the result from Eq. (46) and the term Tr[PBPm] can be rearranged in the following form:

TrZPbBmAmPB+TrPBPm=TrOFmhm+TrOYmhBm+TrOPm,
(48)

where OFm, OYm, and O are defined in  Appendix D. The explicit computation of Pm can be avoided by using a Z-vector ZO as follows:

Tr[OPm]=Tr[ZObm],
(49)

which is obtained by solving a CPSCF equation with O as the right-hand side as follows:

AZO=O.
(50)

With this, we arrive at the final form of the working equation for computing NMR shieldings as follows:

σA=2Tr[(F+OYm)hBm]2Tr[(FB+OFm)hm]2Tr[ZObm].
(51)

Figure 1 summarizes the described algorithm for computing NMR shieldings and shows in which order the steps of the calculation are carried out in our implementation. With the described nested Z-vector approach, seven CPSCF equations need to be solved in total irrespective of the size of the molecule (one equation for ZP and three equations for the magnetic field components of PB and ZO, respectively).

FIG. 1.

Algorithm for ω-RI-CDD-MP2-NMR.

FIG. 1.

Algorithm for ω-RI-CDD-MP2-NMR.

Close modal

The computation of the matrices R̲ and R̄ and their magnetic field derivatives (including the computation of the required three-center integrals) is the most time-consuming step of the shielding calculation with our new method. In this section, we describe how these steps can be accelerated using several approximations.

In order to exploit their rank deficiency, we subject the pseudo-density matrices a pivoted Cholesky decomposition,52–54 

P̲μν=i̲L̲μi̲L̲νi̲,
(52)
P̄μν=āL̄μāL̄νā.
(53)

This decomposition has been applied before in several methods for MP245,65,66 and direct RPA,67 and it was shown to preserve sparsity. It is especially useful for the occupied pseudo-density as the number of obtained Cholesky vectors is less than or equal to the number of occupied orbitals and therefore much smaller than the number of basis functions (assuming reasonably accurate basis sets). From the Cholesky factors of Eqs. (52) and (53), local pseudo-MOs ϕi̲ and ϕā can be computed as follows:

ϕi̲(r)=μχμ(r)L̲μi̲,
(54)
ϕā(r)=νχν(r)L̄νā.
(55)

The pivoted Cholesky decomposition is not applicable to the magnetic field perturbed pseudo-densities because they are not positive semi-definite. For this reason, we leave the perturbed pseudo-densities undecomposed. After inserting the Cholesky-decomposed pseudo-densities, Eqs. (9), (10), (38), and (39) become

R̲νν=i̲j̲b̄i̲ν|j̲b̄2i̲ν|j̲b̄i̲b̄|j̲ν,
(56)
R̄μμ=āj̲b̄μā|j̲b̄2μā|j̲b̄μb̄|j̲ā,
(57)
R̲ννB=μλσμ̲Bν|λ̲σ̄+μ̲ν|λ̲Bσ̄+μ̲ν|λ̲σ̄B×2μν|λσμσ|λν+i̲j̲b̄i̲ν|j̲b̄(B)2i̲ν|j̲b̄i̲b̄|j̲ν+i̲ν|j̲b̄[2i̲ν|j̲b̄(B)i̲b̄|j̲ν(B)],
(58)
R̄μμB=νλσμν̄B|λ̲σ̄+μν̄|λ̲Bσ̄+μν̄|λ̲σ̄B×2μν|λσμσ|λν+āj̲b̄μā|j̲b̄(B)2μā|j̲b̄μb̄|j̲ā+μā|j̲b̄2μā|j̲b̄(B)μb̄|j̲ā(B).
(59)

We furthermore use the RI approximation in order to avoid four-center integrals. As a metric for the RI, we employ an attenuated Coulomb metric46–48 as in our recent work on MP2 energies.45 With this metric, four-center ERIs can be approximated as follows:

formula
(60)

where

formula
(61)
formula
(62)

and

formula
(63)

The attenuated Coulomb metric depends on the parameter ω, which determines accuracy and sparsity in the three-center integrals. For ω equal to zero, the attenuated Coulomb metric reduces to the standard Coulomb metric, which is highly accurate but has long-ranged coupling between the bra and ket of the three-center integrals. On the other hand, in the limit ω → ∞, the rather inaccurate overlap metric is obtained, which leads to a high degree of sparsity in the three-center integrals. It has been shown for MP2 and direct RPA that ω = 0.1 provides a good compromise and gives significant sparsity in the integrals in combination with only small errors compared to the Coulomb metric.45,67 Therefore, we also use ω = 0.1 for all timings shown in this work.

Using the RI approximation allows us to eliminate all transformed four-center integrals from Eqs. (56)(59) by substituting them with three- and two-center RI-integrals. In consequence, several different types of transformed three-center integrals are needed. This includes integrals that are transformed with Cholesky-MO coefficients or with perturbed pseudo-densities and transformed three-center GIAO integrals, such as i̲ā|P(B). In order to compute all of these transformed integrals, the three-center integrals in the atomic orbital basis are computed once and written to disk and then read into memory for each Laplace point. The AO three-center integrals for any auxiliary basis function are copied to a block-sparse matrix; likewise, the Cholesky-MOs and pseudo-density matrices are kept in a block-sparse matrix format. Then, the necessary integral transformations can be carried out efficiently using block-sparse matrix multiplications.

Next, the matrices R̲ and R̄, as well as their B-field derivatives, are computed from the three-center integrals. Figure 2 displays the algorithm used for efficiently computing the Coulomb-type contributions to the matrices R̲ and R̄. For this, the transformed integrals are copied into a natural blocking47,55 data format, which has also been described in Ref. 45 and is used for speeding up the subsequent contractions. As an example, let us discuss the following contraction, which is one of the formally O(N4)-scaling steps in the algorithm from Fig. 2,

formula
(64)

The integrals used in this contraction are arranged such that there is one matrix for each i̲ with initial dimensions of nvirt × naux (where nvirt is the number of virtual Cholesky MOs and naux is the number of auxiliary basis functions). Before carrying out any contraction, rows and columns in these matrices get deleted if they only contain elements that are below the chosen natural blocking threshold. By this, the dimensions of these matrices are reduced and the integral contractions are accelerated.

FIG. 2.

Algorithm for the computation of the Coulomb contributions to R̲ and R̄.

FIG. 2.

Algorithm for the computation of the Coulomb contributions to R̲ and R̄.

Close modal

For a given i̲, the indices of the significant rows and columns are collected in sparse lists {ā}i̲P and {P}i̲ā. The list {ā}i̲P, e.g., contains all virtual pseudo-MOs ϕā that are significant for the orbital ϕi̲. The superscript index P shall denote the third index of the three-center integrals, from which the sparse list is derived (in this case, ). An additional superscript index B is used to denote lists derived from GIAO integrals. The list {ā}i̲P,B, e.g., is derived from the integrals. Note also that two sparse lists are identical if the contained indices refer to the same type of orbital, e.g., {ā}i̲P={b̄}j̲Q. Such sparse lists, which have also been employed in Ref. 45, are used extensively for exploiting sparsity in the algorithms from Figs. 25.

FIG. 3.

Algorithm for the computation of BB. For the definition of C̃, see Fig. 2. In several steps, intersections of two sparse lists are employed.

FIG. 3.

Algorithm for the computation of BB. For the definition of C̃, see Fig. 2. In several steps, intersections of two sparse lists are employed.

Close modal
FIG. 4.

Algorithm for the computation of the Coulomb contributions to R̲B and R̄B. For the definitions of B and C̃, see Fig. 2. BB is computed as shown in Fig. 3.

FIG. 4.

Algorithm for the computation of the Coulomb contributions to R̲B and R̄B. For the definitions of B and C̃, see Fig. 2. BB is computed as shown in Fig. 3.

Close modal
FIG. 5.

Algorithm for the computation of an exchange-type contribution to R̲, R̄, R̲B, or R̄B, which is of the general form w,y,zwx|yzwz|yx, where w, x, y, and z can stand either for Cholesky-MO-indices or for transformed or untransformed AO indices. R can stand for R̲,R̄,R̲B, or R̄B.

FIG. 5.

Algorithm for the computation of an exchange-type contribution to R̲, R̄, R̲B, or R̄B, which is of the general form w,y,zwx|yzwz|yx, where w, x, y, and z can stand either for Cholesky-MO-indices or for transformed or untransformed AO indices. R can stand for R̲,R̄,R̲B, or R̄B.

Close modal

For the contraction in Eq. (64), the scaling can be reduced from formally O(N4) to asymptotically linear because for large systems, only a constant number of virtual pseudo-MOs and auxiliary functions are expected to be significant for a given ϕi̲.

The matrix A from Eq. (64) is then multiplied from both sides with C̃; this step scales cubically but has a very small prefactor. The resulting matrix B in Fig. 2 and the intermediates D and E are all dense matrices. As these matrices are dense, the scaling of the remaining steps in the algorithm shown in Fig. 2 cannot be reduced to asymptotically linear but only to asymptotically quadratic; in these steps, natural blocking is employed again in order to exploit the sparsity in the three-center integrals.

For MP2 shielding calculations, not only the computational cost but also the memory and disk space requirements can be prohibitive. Compared to the implementation of Maurer and Ochsenfeld,29 the memory requirements are substantially reduced by the use of the RI approximation, which allows us to avoid four-center integrals. The employed attenuated Coulomb metric increases the degree of sparsity in the three-center integrals. As in our implementation, the three-center integrals are always kept either in the natural blocking format or a block-sparse data format, the sparsity can be efficiently exploited, and the memory requirements be reduced significantly. In order to further increase the size of systems that can be treated with a given amount of memory, the three-center integrals are stored on disk and only read into memory when needed during the integral transformations and contractions, which are all done in batches. The size of the batches is chosen in accordance with the available memory. Therefore, the available disk space determines the maximum system size that can be treated. In Sec. IV B, we show data for the required disk space in illustrative calculations on glycine chains and DNA strands.

Figures 3 and 4 show the algorithm used for computing the Coulomb-type contributions to the R̲B and R̄B matrices. It is more involved than the corresponding algorithm for the contributions to the unperturbed R̲ and R̄ matrices from Fig. 2 but uses similar intermediates and computation steps. Formally, the algorithm leads to an O(N4) scaling. Again, natural blocking is used for exploiting the sparsity in the three-center integrals. The asymptotic scaling of all steps except for the matrix multiplication C̃ABC̃ can be reduced to either linear or quadratic in this way.

Finally, the algorithm for computing the exchange contributions to R̲ and R̄ and their magnetic field derivatives is shown in Fig. 5. In this algorithm, four-center integrals are built and then contracted as opposed to the algorithms used for the Coulomb contributions. As a first step, three-center integrals are multiplied with the matrix C̃ from Eq. (60) as follows:

formula
(65)

where w and x can denote either Cholesky-MO-indices or transformed or untransformed AO indices. This step scales asymptotically quadratic even if the short-range coupling of the indices in the three-center integrals is exploited because C̃ is a dense matrix. Next, the four-center integrals are built by contracting the intermediate B̃ with three-center integrals as follows:

formula
(66)

For any wy-pair, the computational effort is O(1) because only O(1) orbitals with indices x and z are significant due to exponential coupling between w and x and between y and z. In addition, the summation over Q can be restricted to O(1) elements by exploiting the short-range coupling between Q and y. Therefore, the asymptotic scaling of this step is quadratic. Likewise, also the remaining steps of the algorithm in Fig. 5 scale asymptotically quadratic.

The described method was implemented in a development version of the FermiONs++ program.68–70 Unless stated differently, the def2-SVP basis set71 was used in combination with the corresponding auxiliary basis set72 in all calculations. Two different sets of parameters are used in the calculations shown in Sec. IV A, which we denote as “Loose” and “Tight.” For the “Loose” settings, the attenuation parameter ω is set to 0.1, and 5 Laplace points are used in MP2 and 10 Laplace points in DL-CPSCF73 iterations, which are used for computing PB and solving the Z-vector equations. Furthermore, 10−6 is used as the threshold for deleting rows and columns in the natural blocking format. These values have also been used in our recent work on efficient computation of MP2 energies with the ω-RI-CDD-MP2 method, where we demonstrated that chemically accurate energies can be obtained.45 For the “Tight” settings, the natural blocking threshold and ω are set to zero (meaning that effectively a Coulomb metric is used); the number of Laplace points for both the MP2 shielding part and the DL-CPSCF is increased to 13 in both cases. All timings in Sec. IV B were carried out using the “Loose” settings.

In all calculations, shell-pairs with an overlap of less than 10−12 were neglected. Integral screening with a threshold of 10−10 is employed during the SCF. The Laplace points are determined using a minimax algorithm as described in Refs. 74 and 75; the number of integration points is reduced automatically in our implementation if the fitting interval is small and no improved accuracy can be obtained with more integration points. The pseudo-density matrices are orthogonalized prior to the pivoted Cholesky decomposition as described in Ref. 67; afterward, the orthogonalization is reverted. No frozen-core approximation is employed. For comparison, canonical MP2 shielding calculations were carried out with the Turbomole program.25,76,77 All timings were performed using 64 threads on a node with AMD EPYC 7302 processors, 256 GB RAM, and a solid-state drive (SSD) with a capacity of 1.7 TB.

In order to analyze the influence of the employed approximations on the accuracy, NMR shieldings were computed for all complexes/dimers from the S22 test set78 and for all structures from the benchmark set used by Flaig et al.22 In these calculations, both the “Loose” and “Tight” settings were employed. As a reference, NMR shielding calculations with canonical MP2 performed with the Turbomole program25,76,77 are used. The statistics on the two benchmark sets are shown in Tables I and II. The “Loose” settings, which will be employed also for the timings in Sec. IV B, give satisfactory accuracy for the NMR shieldings, as the mean absolute deviations amount to only 0.024 ppm for the S22 test set and to 0.036 ppm for the benchmark set from the study by Flaig et al.22 No clear improvement compared to the “Loose” settings is observed when using the “Tight” settings. This suggests that the dominant source of error compared to canonical MP2 is the RI approximation; the errors due to the other approximations appear to be negligible for these test sets. Indeed, for both test sets, the already small deviations to Turbomole are reduced further in calculations with “Tight” settings and a cc-pVQZ-RI auxiliary basis set79,80 instead of a def2-SVP-RI basis set. We thus conclude that the RI approximation is the largest source of error, which does not preclude the computation of shieldings, in very good agreement with canonical MP2 shieldings.

TABLE I.

Statistics for MP2 shielding calculations for all dimers from the S22 test set. “Loose” and “Tight” settings are defined in the text. “Tight*” denotes shieldings computed with “Tight” settings and a cc-pVQZ-RI auxiliary basis set instead of a def2-SVP-RI basis set. MP2 shieldings computed with Turbomole25,76,77 are used as a reference. All values are given in ppm.

LooseTightTight*
Mean error 0.016 0.024 0.001 
Mean absolute error 0.024 0.026 0.011 
Root mean squared error 0.044 0.047 0.021 
Maximum error 0.229 0.271 0.151 
LooseTightTight*
Mean error 0.016 0.024 0.001 
Mean absolute error 0.024 0.026 0.011 
Root mean squared error 0.044 0.047 0.021 
Maximum error 0.229 0.271 0.151 
TABLE II.

Statistics for MP2 shielding calculations on the benchmark set from the study by Flaig et al.22 “Loose” and “Tight” settings are defined in the text. “Tight*” denotes shieldings computed with “Tight” settings and a cc-pVQZ-RI auxiliary basis set instead of a def2-SVP-RI basis set. MP2 shieldings computed with Turbomole25,76,77 are used as a reference. All values are given in ppm.

LooseTightTight*
Mean error 0.028 0.033 −0.002 
Mean absolute error 0.036 0.035 0.004 
Root mean squared error 0.098 0.104 0.017 
Maximum error 0.662 0.604 0.150 
LooseTightTight*
Mean error 0.028 0.033 −0.002 
Mean absolute error 0.036 0.035 0.004 
Root mean squared error 0.098 0.104 0.017 
Maximum error 0.662 0.604 0.150 

However, one needs to consider that the molecules in the S22 test set and the test set from the study by Flaig et al.22 are relatively small and that the errors for several of the used approximations, such as the local RI-metric or the sparse linear algebra, might be larger for more extended systems. Therefore, we also carried out calculations on all monomers in the L7 benchmark set81 and display the results in Table III. For the purpose of computational efficiency, we use the SOS-approximation in these calculations; this means that the exchange contributions are neglected, while the Coulomb contributions are scaled with 1.3 as suggested by Jung et al.40 (note that alternatively, for application studies, the basis-set specific scaling parameters from the study by Maurer and Ochsenfeld39 might also be useful). In addition, for the L7 test set, the differences between the “Loose” and “Tight” settings are small with a mean absolute deviation of 0.028 ppm. This suggests that the “Loose” settings are appropriate for computing reliable NMR shieldings also for larger systems; therefore, we use these settings also for the timings shown in Sec. IV B.

TABLE III.

Statistics for SOS-MP2 shielding calculations for all monomers from the L7 test set. SOS-MP2 shieldings computed with FermiONs++68–70 and “Tight” settings are used as a reference. All values are given in ppm.

Loose
 Mean error −0.013  
 Mean absolute error 0.028  
 Root mean square error 0.061  
 Maximum error 0.462  
Loose
 Mean error −0.013  
 Mean absolute error 0.028  
 Root mean square error 0.061  
 Maximum error 0.462  

In Table IV, timings with the new method in shielding calculations on glycine chains are shown. The scaling decreases for larger chain lengths and amounts to 2.74 between the two largest systems, glycine-6 and glycine-10. With this, the observed scaling is well below the formal O(N5) scaling of the MO-based method. This agrees well with the expectations, as the asymptotic scaling of all O(N4) and O(N5) scaling steps is reduced to linear or quadratic if sparsity is exploited (see Sec. II).

TABLE IV.

Timings for MP2 shielding calculations (including exchange) on glycine chains. “nglyc” denotes the number of glycines, and “nbas” is the number of basis functions.

nglycnbasWall time (s)
95 415 
166 2 051 
308 14 063 
450 34 861 
10 734 133 282 
nglycnbasWall time (s)
95 415 
166 2 051 
308 14 063 
450 34 861 
10 734 133 282 

Figure 6 also displays timings on glycine chains; in contrast to the calculations from Table IV, the SOS-approximation is used in all calculations. The wall times needed for the SOS-MP2 shielding calculations are much lower than those for the full MP2 shieldings (e.g., by a factor of 32.9 for glycine-10). This shows that the exchange contributions are the major bottleneck in the MP2 shielding calculations from Fig. 4. Some acceleration for the exchange contributions might be possible by using integral screening, e.g., with Schwarz estimates,82 in conjunction with the natural blocking, as done in our recent work on MP2 energies,45 but we have not explored this so far. For SOS-MP2, the scaling is close to quadratic for the larger chain lengths and thus lower than the formal scaling of SOS-MP2 (O(N4)). Due to the significantly higher computational efficiency obtained with the SOS-approximation, it is likely preferable to use SOS-MP2 for studying very large systems. For this, we recommend to use the basis-set specific scaling factors from the study by Maurer and Ochsenfeld,39 which were shown to provide similar or even better accuracy than canonical MP2.

FIG. 6.

Timings for SOS-MP2 shielding calculations on glycine chains. The colored numbers indicate the effective scaling between two data points.

FIG. 6.

Timings for SOS-MP2 shielding calculations on glycine chains. The colored numbers indicate the effective scaling between two data points.

Close modal

As discussed, memory and/or disk space can also be a bottleneck in MP2 shielding calculations. Therefore, we show data for the used disk space during calculations on glycine chains in Fig. 7. For the largest system (glycine-40), roughly 437 GB of disk space were used. The scaling continuously decreases with the increasing chain length and amounts to 1.3 between glycine-30 and glycine-40. This is significantly lower than the formal cubic scaling for the memory demands of the three-center integrals and is enabled by the sparsity provided by the attenuated Coulomb metric and by the employed sparse data formats (block sparse matrices and natural blocking), which allow us to neglect many of the insignificant integrals.

FIG. 7.

Maximum disk usage for three-center integrals in SOS-MP2 shielding calculations on glycine chains. The colored numbers indicate the effective scaling between two data points.

FIG. 7.

Maximum disk usage for three-center integrals in SOS-MP2 shielding calculations on glycine chains. The colored numbers indicate the effective scaling between two data points.

Close modal
TABLE V.

Wall times, scaling, and needed disk space in calculations on DNA strands with up to two adenine–thymine base pairs. “nbas” is the number of basis functions.

AT1AT2
nbas 625 1332 
Wall time (h) 1.0 14.0 
Scaling ⋯ 3.44 
Disk space (GB) 59.5 531.4 
AT1AT2
nbas 625 1332 
Wall time (h) 1.0 14.0 
Scaling ⋯ 3.44 
Disk space (GB) 59.5 531.4 

In addition to the glycine chains, also timings on short DNA strands with one and two adenine–thymine base pairs were carried out (see Table V). These molecules are more representative of potential applications on large biomolecular systems. The calculation on the largest DNA strand, AT2, with 128 atoms and 1332 basis functions took 11 h, which shows that the new method allows to treat systems in this size range with moderate computational effort. The scaling of the compute time between AT1 and AT2 amounts to 3.4. This is below the formal O(N4) scaling of Laplace-transformed SOS-MP2, which indicates that there is already usable sparsity for systems in this size range. As expected, due to the less extended shape of the AT base pairs, the observed scaling is larger than for glycine chains with a similar number of basis functions. In Table VI, timings for several of the calculation steps are shown. Most of the time during the SOS-MP2 shielding calculations is used for the transformations of the three-center integrals followed by the computation of the perturbed R̲B and R̄B matrices. Significantly less compute time is needed for the unperturbed R̲ and R̄ matrices, the recursion formulas, and the nested Z-vector approach. In Table V, also the maximum amount of disk space used during the shielding calculations on the DNA strands is shown. Significantly more disk space than that for glycine chains with a similar number of basis functions is needed, which is caused by the lower sparsity. For AT2, 531.4 GB disk space is needed for storing three-center integrals. Calculations on larger DNA strands were not possible with the available disk space on the employed nodes. In order to improve the applicability to larger systems, it is therefore important to circumvent the disk space bottleneck. A potential solution would be an integral-direct approach, by which we plan to address the disk space bottleneck in a future publication.

TABLE VI.

Timings of individual steps during a SOS-MP2 shielding calculation on AT2.

Wall time (s)
 Integral transformation 15 702  
 R̲ and R̄ 1 205  
 R̲B and R̄B 10 782  
 Y̲B and ȲB 2 795  
 Nested Z-vector 4 470  
 Total 39 619  
Wall time (s)
 Integral transformation 15 702  
 R̲ and R̄ 1 205  
 R̲B and R̄B 10 782  
 Y̲B and ȲB 2 795  
 Nested Z-vector 4 470  
 Total 39 619  

We presented an efficient method for computing NMR shieldings at the MP2 level of theory, which is based on Laplace-transformed AO-MP2. It employs many of the approximations that were also used in our recently published ω-RI-CDD-MP2 method45 for correlation energies, including RI with an attenuated Coulomb metric and Cholesky decomposition of pseudo-density matrices. Furthermore, block sparse linear algebra and natural blocking are used in order to speed up the computations and reduce the memory needed for the three-center integrals. A nested Z-vector allows us to efficiently compute the entire set of NMR shieldings for all nuclei (as opposed to the selected-nuclei formulation of Maurer and Ochsenfeld29) and requires the solution of only seven CPSCF equations irrespective of the molecule size. Benchmark calculations indicated that the same thresholds and settings that were shown previously to give accurate correlation energies with the ω-RI-CDD-MP2 method also allow us to compute accurate NMR shieldings. Timings on glycine chains show close to quadratic scaling for larger systems. Particularly high efficiency is obtained in combination with an SOS-approximation. In addition, the scaling of the needed disk space decreases with molecule size and becomes sub-quadratic for larger chain lengths. Calculations on DNA strands with up to two base pairs illustrate the potential of the method for computing highly accurate NMR shieldings of biomolecular systems. The new method thus enables the computation of NMR shieldings of large molecules with significantly reduced computational cost and memory requirements.

See the supplementary material for all computed NMR shieldings from Sec. IV A and the employed xyz-structures for DNA and glycine chains.

The authors acknowledge the financial support of the “Deutsche Forschungsgemeinschaft” (DFG) under Grant No. SFB 1309-32587107 and the cluster of excellence (Grant No. EXC2111-390814868) “Munich Center for Quantum Science and Technology” (MCQST). M.G. and S.V. acknowledge the “Studienstiftung des Deutschen Volkes” for graduate fellowships. C.O. acknowledges additional financial support as a Max-Planck-Fellow at MPI-FKF Stuttgart.

The authors have no conflicts to disclose.

The data that support the findings of this study are available within the article and its supplementary material.

In this section, the recursion formulas for a matrix Y are derived such that Tr[YAξ]=TrBeAξ holds. First, we insert the Taylor series expansion of the matrix exponential

TrBeAξ=TrB1+A+12!AA+13!AAA+ξ
(A1)

and subsequently differentiate each term using the product rule

TrBeAξ=TrB0+Aξ+12!AξA+AAξ+13!AξAA+AAξA+AAAξ+.
(A2)

In each term, cyclic permutations under the trace are used to move Aξ to the far right,

TrBeAξ=TrB+12!AB+BA+13!AAB+ABA+BAA+Aξ.
(A3)

Upon comparing with Eq. (17), the term in round brackets can then be identified with Y,

Y=B+12!AB+BA+13!AAB+ABA+BAA+=k=1i=0k11k!AiBAk1i.
(A4)

Y can be efficiently computed using recursion as shown in the following equations:

Y=k=0kmaxYk,
(A5)
Yk=1kAYk1+BAk1k>0,
(A6)
Y0=0,
(A7)
Ak1k!Ak.
(A8)

In order to ensure a fast convergence and avoid numerical difficulties, we employ a modified version of the “scaling and squaring” approach,83,84 which is commonly used for computing matrix exponentials. Instead of using A and B in Eqs. (A5)(A8), two matrices A′ and B′ are employed as follows:

A=1nA,
(A9)

where n is an integer number chosen such that |A′| < 0.5, and

B=i=0n1expAiBexpAn1i.
(A10)

From the result Y′, one can then compute the matrix Y as follows:

Y=1nY.
(A11)

For computing MP2 shieldings, also derivatives of the Y matrices from  Appendix A are needed. These perturbed Y matrices can also be computed using recursions, which are obtained by differentiating the recursion from Eqs. (A5)(A8) as follows:

Yξ=k=0kmaxYkξ,
(B1)

with

Ykξ=1kAξYk1+AYk1ξ+BξAk1+BAk1ξk>0,
(B2)
Y0ξ=0.
(B3)

Ak1ξ in Eq. (B2) can also be computed recursively as follows:

Akξ=1kAAk1ξ+AξAk1k>0,
(B4)
A0ξ=0.
(B5)

Faster convergence can be obtained using a modified scaling and squaring approach. For this, A′ = A/n and Aξ = Aξ/n should be used instead of A and Aξ in the recursion from Eqs. (B2) and (B4). Furthermore, B needs to be replaced with B′ from Eq. (A10). The result Yξ, which is computed by carrying out the recursion from Eq. (B2) with A′, Aξ, B′, and Bξ instead of their respective unprimed counterparts, finally needs to be multiplied with 1n in order to obtain Yξ.

The Taylor series expansion of a matrix exponential eA is given by

eA=1+A+12!AA+13!AAA+=k=01k!Ak.
(C1)

The matrix exponential can be evaluated efficiently using recursion as follows:

eA=k=0ek,
(C2)
ek=1kek1Ak>0,
(C3)
e0=1.
(C4)

The recursion from Eq. (C2) can only be expected to converge quickly if the norm of A is small. If this is not the case, “scaling and squaring”83,84 can be applied. For this, an integer number n is chosen such that A/n<0.5. Then, the matrix exponential eA/n is computed using the recursion from Eq. (C2). The exponential of A is obtained by taking eA/n to the power n as follows:

eA=eA/nn.
(C5)

For the presented MP2-NMR method, also perturbed matrix exponentials of the form (eA)ξ need to be computed. A formula for the recursive computation of the perturbed matrix exponential can be derived by differentiating Eq. (C2) as follows:

(eA)ξ=k=0ekξ,
(C6)
ekξ=1kek1ξA+ek1Aξk>0,
(C7)
e0ξ=0.
(C8)

If scaling and squaring is used, computing eAξ requires the recursive computation of eA/nξ and eA/n as follows:

(eA)ξ=(eA/n)ξ(eA/n)n1+(eA/n)1(eA/n)ξ(eA/n)n2++(eA/n)n1(eA/n)ξ=i=0n1(eA/n)i(eA/n)ξ(eA/n)n1i.
(C9)

An efficient evaluation of this formula is possible via another recursion,

(eA)ξ=En,
(C10)
En=(eA/n)En1+(eA/n)ξ(eA/n)n1n>1,
(C11)
E1=(eA/n)ξ.
(C12)

There are different formulations of CPSCF, which lead to different equations for A and bB in Eq. (40). Here, we use the D-CPSCF formulation from Ref. 64, which results in the following expressions for A and bB:

APB=3FPBS+3SPBF2FPBSPS2SPSPBF4FPSPBS4SPBSPF+GXPS+SPGX2SPGXPS,
(D1)

where

X=PBSP+PSPB2PSPBSP
(D2)

and

bB=FPSB+SBPF+2FPSBPS+2SPSBPFYPSSPY+2SPYPS,
(D3)

with

Y=F(B)GPSBP.
(D4)

The difference between FB and F(B) is that the latter does by definition not include contributions containing the perturbed density PB. For computing the second derivative of the density matrix PBm, the nuclear magnetic moment derivatives AmPB and bBm are also needed as shown in Eq. (45),

AmPB=3FmPBS+3SPBFm2FmPBSPS2FPBSPmS2SPmSPBF2SPSPBFm4FmPSPBS4FPmSPBS4SPBSPmF4SPBSPFm+GXmPS+GXPmS+SPmGX+SPGXm2SPmGXPS2SPGXmPS2SPGXPmS,
(D5)

where

Xm=PBSPm+PmSPB2PmSPBSP2PSPBSPm
(D6)

and

bBm=FmPSB+FPmSB+SBPmF+SBPFm+2FmPSBPS+2FPmSBPS+2FPSBPmS+2SPmSBPF+2SPSBPmF+2SPSBPFmYmPSYPmSSPmYSPYm+2SPmYPS+2SPYmPS+2SPYPmS,
(D7)

where

Ym=hBm+GB[Pm]GPmSBP+PSBPm.
(D8)

As it is commonly done in implementations of NMR shieldings, we do not use complex-valued matrices explicitly, but instead we treat purely imaginary matrices using skew-symmetric, real-valued matrices. In this case, one needs to be careful to use the correct sign for each term; whenever a term contains two imaginary matrices, like the term FmPSB with Fm and SB, an additional sign change needs to be applied because both imaginary matrices carry a factor i and i2 = −1.

Instead of iteratively solving for the full PBm as suggested by Eq. (45), we only solve for the occupied-virtual (ov) and virtual-occupied (vo) subspace projections. For the occupied–occupied (oo) and virtual–virtual (vv) subspace projections, explicit expressions can be derived by differentiating the idempotency condition PSP = P twice and projecting onto the corresponding subspaces as follows:29 

Tr[PPBm]=Tr[PPovBm+PvoBm]+Tr[PPooBm+PvvBm],
(D9)

with

PooBm+PvvBm=2PSPBSPmSP+PSPmSBP+PSBPmSP+PSPmSPBSP+PBSPm+PmSPB+PmSBP+PSBPm.
(D10)

The computation of the ov- and vo-projections of PBm in the first term on the right-hand side of Eq. (D9) can be avoided using a Z-vector approach,

Tr[P(PovBm+PvoBm)]=Tr[Zov+ZvobBmAmPB].
(D11)

Here, Zov and Zvo are the ov- and vo-projections of the Z-vector ZP=A1P.

After eliminating PBm from the equations, Pm remains to be eliminated. For this, we will show how terms involving Pm can be rearranged in the form Tr[…Pm], which then allows us to apply another Z-vector approach. The terms involving PooBm + PvvBm can be brought into this form by applying cyclic permutations under the trace,

TrPPooBm+PvvBm=TrO[PooBm+PvvBm]Pm,
(D12)

where

O[PooBm+PvvBm]=2SPPPSPBS+SBPPPS+SPPPSB+SPBSPPPS+PPBS+SPBP+SBPP+PPSB.
(D13)

Next, we focus on the expression Tr[Zov+ZvobBmAmPB]. First, we insert Eqs. (D5) and (D7). After applying cyclic permutations, we arrive at

Tr[Zov+ZvobBmAmPB]=TrOPmPm+OFmFm+OYm(G[Xm]+Ym),
(D14)

where

OPm=SBZ̄F+FZ̄SB+2SBPSZ̄F+2SZ̄FPSB+2SBPFZ̄S+2FZ̄SPSBSZ̄YYZ̄S+2YPSZ̄S+2SZ̄SPY+2SZ̄FPBS+2SPBFZ̄S+4SPBSZ̄F+4FZ̄SPBSSZ̄G[X]G[X]Z̄S+2G[X]PSZ̄S+2SZ̄SPG[X],
(D15)
OFm=PSBZ̄+Z̄SBP+2PSBPSZ̄+2Z̄SPSBP3PBSZ̄3Z̄SPB+2PBSPSZ̄+2Z̄SPSPB+4PSPBSZ̄+4Z̄SPBSP,
(D16)

and

OYm=PSZ̄Z̄SP+2PSZ̄SP,
(D17)

where Z̄Zov+Zvo.

Equation (D14) contains the Fock matrix derivative Fm, which also needs to be circumvented in the all-nuclei formulation. For this, we use the following identities:

TrWFm=TrWhm+TrWGPm,
(D18)
TrWGPm=TrGWPm,
(D19)

where W is some general matrix and GW is

Gμν[W]λσWλσ2μν|λσμσ|λν.
(D20)

Using Eqs. (D18) and (D19), the contribution from Fm in Eq. (D14) can be rearranged as follows:

TrOFmFm=TrOFmhm+TrGOFmPm.
(D21)

Similarly, the terms involving GXm and Ym in Eq. (D14) can be rearranged by using Eq. (D19) and then applying cyclic permutations as follows:

TrOYmG[Xm]=TrG[OYm]Xm=TrG[OYm]PBS+SPBG[OYm]2SPBSPG[OYm]2G[OYm]PSPBSPm,
(D22)
TrOYmYm=TrOYmhBm+TrGB[OYm]PmTrG[OYm]PmSBP+PSBPm=TrOYmhBm+TrGB[OYm]PmTrSBPG[OYm]+G[OYm]PSBPm.
(D23)

Using the results from Eqs. (D18)(D23) allows us to rewrite the expression Tr[Zov+Zvo(bBmAm[PB])] in the following form:

Tr[(Zov+Zvo)bBmAmPB]=TrOFmhm+TrOYmhBm+TrÕPm,
(D24)

where

Õ=OPm+GOFm+GB[OYm]+G[OYm]PBS+SPBG[OYm]2SPBSPG[OYm]2G[OYm]PSPBSSBPG[OYm]G[OYm]PSB.
(D25)

Combining Õ, O[PooBm+PvvBm], and PB from the term Tr[PBPm] allows us to compactly express all remaining terms involving Pm in the form TrOPm,

O=Õ+O[PooBm+PvvBm]+PB.
(D26)

This term now has the proper form for applying a Z-vector approach, which allows us to circumvent Pm,

TrOPm=TrZObm,
(D27)

where

ZO=A1O.
(D28)

In our implementation, the Z-vector equation for ZO is solved iteratively using the DL-CPSCF approach from the study by Beer and Ochsenfeld.73 

1.
J.
Gauss
,
J. Chem. Phys.
99
,
3629
(
1993
).
2.
T.
Helgaker
,
M.
Jaszuński
, and
K.
Ruud
,
Chem. Rev.
99
,
293
(
1999
).
3.
L. B.
Casabianca
and
A. C.
De Dios
,
J. Chem. Phys.
128
,
052201
(
2008
).
4.
5.
K.
Wolinski
,
J. F.
Hinton
, and
P.
Pulay
,
J. Am. Chem. Soc.
112
,
8251
(
1990
).
6.
M.
Häser
,
R.
Ahlrichs
,
H.
Baron
,
P.
Weis
, and
H.
Horn
,
Theor. Chim. Acta
83
,
455
(
1992
).
7.
R. M.
Stevens
,
R. M.
Pitzer
, and
W. N.
Lipscomb
,
J. Chem. Phys.
38
,
550
(
1963
).
8.
H.
Larsen
,
T.
Helgaker
,
J.
Olsen
, and
P.
Jørgensen
,
J. Chem. Phys.
115
,
10344
(
2001
).
9.
G.
Schreckenbach
and
T.
Ziegler
,
J. Phys. Chem.
99
,
606
(
1995
).
10.
J. R.
Cheeseman
,
G. W.
Trucks
,
T. A.
Keith
, and
M. J.
Frisch
,
J. Chem. Phys.
104
,
5497
(
1996
).
11.
G.
Rauhut
,
S.
Puyear
,
K.
Wolinski
, and
P.
Pulay
,
J. Phys. Chem.
100
,
6310
(
1996
).
12.
C.
Ochsenfeld
,
J.
Kussmann
, and
F.
Koziol
,
Angew. Chem.
116
,
4585
(
2004
).
13.
J.
Kussmann
and
C.
Ochsenfeld
,
J. Chem. Phys.
127
,
054103
(
2007
).
14.
M.
Beer
,
J.
Kussmann
, and
C.
Ochsenfeld
,
J. Chem. Phys.
134
,
074102
(
2011
).
15.
16.
J.
Gauss
and
J. F.
Stanton
,
Adv. Chem. Phys.
123
,
355
(
2002
).
17.
M. R.
Hoffmann
,
D. J.
Fox
,
J. F.
Gaw
,
Y.
Osamura
,
Y.
Yamaguchi
,
R. S.
Grev
,
G.
Fitzgerald
,
H. F.
Schaefer
 III
,
P. J.
Knowles
, and
N. C.
Handy
,
J. Chem. Phys.
80
,
2660
(
1984
).
18.
K.
Ruud
,
T.
Helgaker
,
R.
Kobayashi
,
P.
Jørgensen
,
K. L.
Bak
, and
H. J. A.
Jensen
,
J. Chem. Phys.
100
,
8178
(
1994
).
19.
J.
Gauss
and
J. F.
Stanton
,
J. Chem. Phys.
102
,
251
(
1995
).
20.
J.
Gauss
and
J. F.
Stanton
,
J. Chem. Phys.
104
,
2574
(
1996
).
21.
M.
Kállay
and
J.
Gauss
,
J. Chem. Phys.
120
,
6841
(
2004
).
22.
D.
Flaig
,
M.
Maurer
,
M.
Hanni
,
K.
Braunger
,
L.
Kick
,
M.
Thubauville
, and
C.
Ochsenfeld
,
J. Chem. Theory Comput.
10
,
572
(
2014
).
23.
N. C.
Handy
and
H. F.
Schaefer
 III
,
J. Chem. Phys.
81
,
5031
(
1984
).
24.
Y.
Yamaguchi
,
H. F.
Schaefer
 III
,
Y.
Osamura
,
J. D.
Goddard
 et al,
A New Dimension to Quantum Chemistry: Analytic Derivative Methods in Ab Initio Molecular Electronic Structure Theory
(
Oxford University Press
,
1994
).
25.
M.
Kollwitz
and
J.
Gauss
,
Chem. Phys. Lett.
260
,
639
(
1996
).
26.
J.
Gauss
and
H.-J.
Werner
,
Phys. Chem. Chem. Phys.
2
,
2083
(
2000
).
27.
S.
Loibl
and
M.
Schütz
,
J. Chem. Phys.
137
,
084107
(
2012
).
28.
G. L.
Stoychev
,
A. A.
Auer
,
J.
Gauss
, and
F.
Neese
,
J. Chem. Phys.
154
,
164110
(
2021
).
29.
M.
Maurer
and
C.
Ochsenfeld
,
J. Chem. Phys.
138
,
174104
(
2013
).
30.
J.
Almlöf
,
Chem. Phys. Lett.
181
,
319
(
1991
).
31.
M.
Häser
and
J.
Almlöf
,
J. Chem. Phys.
96
,
489
(
1992
).
32.
M.
Häser
,
Theor. Chim. Acta
87
,
147
(
1993
).
33.
P. Y.
Ayala
and
G. E.
Scuseria
,
J. Chem. Phys.
110
,
3660
(
1999
).
34.
S. A.
Maurer
,
D. S.
Lambrecht
,
D.
Flaig
, and
C.
Ochsenfeld
,
J. Chem. Phys.
136
,
144107
(
2012
).
35.
S. A.
Maurer
,
D. S.
Lambrecht
,
J.
Kussmann
, and
C.
Ochsenfeld
,
J. Chem. Phys.
138
,
014101
(
2013
).
36.
S.
Grimme
,
J. Chem. Phys.
124
,
034108
(
2006
).
37.
G. L.
Stoychev
,
A. A.
Auer
, and
F.
Neese
,
J. Chem. Theory Comput.
14
,
4756
(
2018
).
38.
S.
Grimme
,
J. Chem. Phys.
118
,
9095
(
2003
).
39.
M.
Maurer
and
C.
Ochsenfeld
,
J. Chem. Theory Comput.
11
,
37
(
2015
).
40.
Y.
Jung
,
R. C.
Lochan
,
A. D.
Dutoi
, and
M.
Head-Gordon
,
J. Chem. Phys.
121
,
9793
(
2004
).
42.
W.
Kutzelnigg
,
Isr. J. Chem.
19
,
193
(
1980
).
43.
M.
Feyereisen
,
G.
Fitzgerald
, and
A.
Komornicki
,
Chem. Phys. Lett.
208
,
359
(
1993
).
44.
J. L.
Whitten
,
J. Chem. Phys.
58
,
4496
(
1973
).
45.
M.
Glasbrenner
,
D.
Graf
, and
C.
Ochsenfeld
,
J. Chem. Theory Comput.
16
,
6856
(
2020
).
46.
Y.
Jung
,
A.
Sodt
,
P. M. W.
Gill
, and
M.
Head-Gordon
,
Proc. Natl. Acad. Sci. U. S. A.
102
,
6692
(
2005
).
47.
Y.
Jung
,
Y.
Shao
, and
M.
Head-Gordon
,
J. Comput. Chem.
28
,
1953
(
2007
).
48.
S.
Reine
,
E.
Tellgren
,
A.
Krapp
,
T.
Kjærgaard
,
T.
Helgaker
,
B.
Jansik
,
S.
Høst
, and
P.
Salek
,
J. Chem. Phys.
129
,
104101
(
2008
).
49.
S.
Burger
,
F.
Lipparini
,
J.
Gauss
, and
S.
Stopkowicz
,
J. Chem. Phys.
155
,
074105
(
2021
).
50.
N. H. F.
Beebe
and
J.
Linderberg
,
Int. J. Quantum Chem.
12
,
683
(
1977
).
51.
H.
Koch
,
A.
Sánchez de Merás
, and
T. B.
Pedersen
,
J. Chem. Phys.
118
,
9481
(
2003
).
52.
F.
Aquilante
,
T. B.
Pedersen
,
A.
Sánchez de Merás
, and
H.
Koch
,
J. Chem. Phys.
125
,
174101
(
2006
).
53.
N. J.
Higham
,
Wiley Interdiscip. Rev.: Comput. Stat.
1
,
251
(
2009
).
54.
H.
Harbrecht
,
M.
Peters
, and
R.
Schneider
,
Appl. Numer. Math.
62
,
428
(
2012
).
55.
Y.
Jung
and
M.
Head-Gordon
,
Phys. Chem. Chem. Phys.
8
,
2831
(
2006
).
56.
S.
Schweizer
,
B.
Doser
, and
C.
Ochsenfeld
,
J. Chem. Phys.
128
,
154101
(
2008
).
57.
P. R.
Surján
,
Chem. Phys. Lett.
406
,
318
(
2005
).
58.
C. A.
White
,
B. G.
Johnson
,
P. M. W.
Gill
, and
M.
Head-Gordon
,
Chem. Phys. Lett.
230
,
8
(
1994
).
59.
C. A.
White
,
B. G.
Johnson
,
P. M. W.
Gill
, and
M.
Head-Gordon
,
Chem. Phys. Lett.
253
,
268
(
1996
).
60.
F.
Weigend
,
Phys. Chem. Chem. Phys.
4
,
4285
(
2002
).
61.
C.
Ochsenfeld
,
C. A.
White
, and
M.
Head-Gordon
,
J. Chem. Phys.
109
,
1663
(
1998
).
62.
C.
Ochsenfeld
,
Chem. Phys. Lett.
327
,
216
(
2000
).
63.
H.
Laqua
,
T. H.
Thompson
,
J.
Kussmann
, and
C.
Ochsenfeld
,
J. Chem. Theory Comput.
16
,
1456
(
2020
).
64.
C.
Ochsenfeld
and
M.
Head-Gordon
,
Chem. Phys. Lett.
270
,
399
(
1997
).
65.
J.
Zienau
,
L.
Clin
,
B.
Doser
, and
C.
Ochsenfeld
,
J. Chem. Phys.
130
,
204112
(
2009
).
66.
S. A.
Maurer
,
L.
Clin
, and
C.
Ochsenfeld
,
J. Chem. Phys.
140
,
224112
(
2014
).
67.
A.
Luenser
,
H. F.
Schurkus
, and
C.
Ochsenfeld
,
J. Chem. Theory Comput.
13
,
1647
(
2017
).
68.
J.
Kussmann
and
C.
Ochsenfeld
,
J. Chem. Phys.
138
,
134114
(
2013
).
69.
J.
Kussmann
and
C.
Ochsenfeld
,
J. Chem. Theory Comput.
11
,
918
(
2015
).
70.
J.
Kussmann
and
C.
Ochsenfeld
,
J. Chem. Theory Comput.
13
,
3153
(
2017
).
71.
F.
Weigend
and
R.
Ahlrichs
,
Phys. Chem. Chem. Phys.
7
,
3297
(
2005
).
72.
F.
Weigend
,
M.
Häser
,
H.
Patzelt
, and
R.
Ahlrichs
,
Chem. Phys. Lett.
294
,
143
(
1998
).
73.
M.
Beer
and
C.
Ochsenfeld
,
J. Chem. Phys.
128
,
221102
(
2008
).
74.
M.
Kaltak
,
J.
Klimeš
, and
G.
Kresse
,
J. Chem. Theory Comput.
10
,
2498
(
2014
).
75.
A.
Takatsuka
,
S.
Ten-No
, and
W.
Hackbusch
,
J. Chem. Phys.
129
,
044112
(
2008
).
76.
TURBOMOLE V7.4 2019, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989-2007, TURBOMOLE GmbH, since 2007; available from http://www.turbomole.com.
77.
M.
Kollwitz
,
M.
Häser
, and
J.
Gauss
,
J. Chem. Phys.
108
,
8295
(
1998
).
78.
P.
Jurečka
,
J.
Šponer
,
J.
Černỳ
, and
P.
Hobza
,
Phys. Chem. Chem. Phys.
8
,
1985
(
2006
).
79.
C.
Hättig
,
Phys. Chem. Chem. Phys.
7
,
59
(
2005
).
80.
F.
Weigend
,
A.
Köhn
, and
C.
Hättig
,
J. Chem. Phys.
116
,
3175
(
2002
).
81.
R.
Sedlak
,
T.
Janowski
,
M.
Pitonak
,
J.
Rezac
,
P.
Pulay
, and
P.
Hobza
,
J. Chem. Theory Comput.
9
,
3364
(
2013
).
82.
M.
Häser
and
R.
Ahlrichs
,
J. Comput. Chem.
10
,
104
(
1989
).
83.
C.
Moler
and
C.
Van Loan
,
SIAM Rev.
45
,
3
(
2003
).
84.
T.
Fung
,
Int. J. Numer. Methods Eng.
59
,
1273
(
2004
).

Supplementary Material