We present a derivation and efficient implementation of the formally complete analytic second derivatives for the domain-based local pair natural orbital second order Møller–Plesset perturbation theory (MP2) method, applicable to electric or magnetic field-response properties but not yet to harmonic frequencies. We also discuss the occurrence and avoidance of numerical instability issues related to singular linear equation systems and near linear dependences in the projected atomic orbital domains. A series of benchmark calculations on medium-sized systems is performed to assess the effect of the local approximation on calculated nuclear magnetic resonance shieldings and the static dipole polarizabilities. Relative deviations from the resolution of the identity-based MP2 (RI-MP2) reference for both properties are below 0.5% with the default truncation thresholds. For large systems, our implementation achieves quadratic effective scaling, is more efficient than RI-MP2 starting at 280 correlated electrons, and is never more than 5–20 times slower than the equivalent Hartree–Fock property calculation. The largest calculation performed here was on the vancomycin molecule with 176 atoms, 542 correlated electrons, and 4700 basis functions and took 3.3 days on 12 central processing unit cores.

## I. INTRODUCTION

The early development of local electron correlation methods is due to Pulay and Sæbø,^{1–3} and thanks to recent advancements in several research groups,^{4–41} the popularity and applicability of these methods for the calculation of relative energies have grown tremendously. With modern approximations, the “gold standard” method of computational chemistry—coupled cluster theory with single, double, and perturbative triple [CCSD(T)] excitations^{42}—has become available at a computational cost only a few times higher than that of density functional theory (DFT).^{43,44} Considering the massive progress that has been made in this field for the calculation of electronic energies, there are comparatively few works that make use of local correlation approximations to compute molecular properties. This is, in part, because many observable properties—such as the dipole moment and polarizability, nuclear magnetic resonance (NMR) shielding and spin–spin coupling, and electron paramagnetic resonance (EPR) g-tensor and hyperfine coupling (HFC)—are related to derivatives of the energy, and analytic derivatives of local correlation methods are challenging due to the complexity of the theories.

This is not to say that nothing has been accomplished on this subject. Werner, Schütz, and co-workers derived and implemented analytic nuclear gradients for several local correlation methods^{45–50} based on projected atomic orbitals (PAOs),^{1,2} including second order Møller–Plesset perturbation theory (MP2), and also used them to semi-numerically calculate vibrational frequencies.^{51–53} Gauss and Werner also presented a pilot implementation for NMR shielding calculations with PAO-based local MP2 (LMP2).^{54} The first efficient such implementation, also employing the resolution of the identity (RI) approximation, was later reported by Loibl and Schütz^{55} and also adapted to magnetizabilities.^{56} Another study by Werner and co-workers examines the accuracy of PAO-based local correlation methods for polarizability calculations via finite differences,^{57} an approach that is difficult to apply to magnetic properties as it requires an implementation based on complex algebra. Maurer and Ochsenfeld developed an AO-based Laplace-transformed MP2 method for the computation of NMR shieldings, which can asymptotically achieve linear scaling or even sub-linear if only a few nuclei are of interest.^{58} Static and dynamic polarizabilities have also been implemented for MP2 and the approximate coupled cluster method CC2 using both RI and a Laplace-transformation.^{59} Frank *et al.* reported analytic gradients for a local MP2 variant, based on pair natural orbitals (PNOs), although they neglected pair natural orbitals (PNOs) relaxation.^{60} The complete analytic gradient for an orbital-specific virtual (OSV) local MP2 method was also published recently by Zhou *et al.*^{61} We should also mention here the work by Crawford and co-workers in the field of local coupled cluster linear response theory for the calculation of molecular properties, including frequency-dependent polarizabilities and specific optical rotations.^{62–65} Work by our group in Mülheim has focused on the domain-based local pair natural orbital (DLPNO) approximation and has resulted in orbital-unrelaxed first derivatives of DLPNO-CCSD, which can be used to calculate, e.g., dipole moments and HFCs.^{66,67} Fully orbital-relaxed first derivatives and nuclear gradients were also implemented for the DLPNO-MP2 method,^{68,69} and those results highlight the importance of the PNO relaxation contributions, e.g., for electric field gradients.^{68}

In the current work, we present a derivation of analytic second derivatives of DLPNO-MP2, applicable to field-response properties, such as dipole polarizabilities and NMR shieldings. Section II focuses on the theory and working equations, while Sec. V provides a detailed benchmark study of the accuracy and efficiency of the method. To our knowledge, this is the first implementation of formally complete analytic second derivatives of a PNO-based local correlation method. On the one hand, MP2 is the simplest post-Hartree–Fock (HF) correlation method, and the present work is a stepping stone toward DLPNO-CCSD second derivatives, which we intend to pursue in the future. The pitfalls we encountered and their proposed solutions, discussed in Sec. III, are also likely relevant for other local correlation methods. On the other hand, MP2 itself and even more so the related double-hybrid DFT (DHDFT) have been shown to be very accurate for both polarizabilities and NMR shieldings.^{70–73} The present work extends their applicability to much larger systems and should therefore be viewed as part of an ongoing effort to reduce the computational cost of MP2 response property calculations. Works in this context are the integral-direct gauge-including atomic orbital (GIAO) MP2 implementation for shieldings of Kollwitz, Häser, and Gauss,^{74,75} the derivation of RI-MP2 second derivatives in combination with chain-of-spheres exchange (COSX),^{76} the Laplace-based approaches of Ochsenfeld, Hättig, and their co-workers,^{25,59} the proof-of-concept GIAO-LMP2 implementation,^{54} and, in particular, the efficient RI-based version of the latter.^{55,56} Note that the present derivation and implementation are not applicable to harmonic vibrational frequencies, as the efficient evaluation of the nuclear Hessian requires a substantially different algorithm.

## II. THEORY

Both the polarizability tensor **α** and the shielding tensor for nucleus *A*, **σ**^{A}, can be obtained as second derivatives of the energy,^{77,78}

where *E*, **D**, and **h** are the method-specific energy expression, density matrix, and one-electron part of the Fock matrix, respectively; *x* and *y* denote arbitrary Cartesian directions; and **F**^{ext}, **B**^{ext}, and **M**^{A} are the external electric and magnetic fields and the nuclear magnetic moment of *A*, respectively. Note that the expressions in Eqs. (1) and (2) are not symmetric with respect to the perturbation and thus not the only possible expressions for these properties. Note also that **h** only contains terms that are linear in **F**^{ext}. In the following discussion, **λ** will denote a generic multi-dimensional perturbation (such as **F**^{ext} or **B**^{ext}), and derivatives with respect to it will be signified as upper indices (e.g., **D**^{λ}). We note here that in the case of magnetic field perturbations, GIAOs^{77,79–82} are used to tackle the gauge-origin problem.^{83,84}

Analytic derivatives of MP2 and variants thereof have been thoroughly discussed in the literature,^{85–94} so we will provide rather limited background here. We will focus on DLPNO-MP2 second derivatives, in particular, and remark on similarities and differences to the DLPNO-MP2 gradient.^{68,69} One point of note is that using a Lagrangian formulation together with the second derivative expressions above, which are “asymmetric” with respect to the two perturbations, is particularly advantageous for magnetic shielding calculations (though not necessarily for the polarizability tensor and other properties) as it only requires derivatives of the wave function parameters and the two-electron integrals with respect to the three magnetic field components.^{78} In contrast, the “symmetric” formulation, commonly used for vibrational frequency calculations,^{88–90} would also require derivatives with respect to the 3*N*_{atoms} components of the nuclear magnetic moments. The downside of the “asymmetric” approach is that field derivatives of all other parameters of the Lagrangian are necessary, and thus, perturbed (i.e., first-order) Z-vector equations must be solved.

### A. Notation

Due to the many basis transformations and domain truncations involved in DLPNO methods, mathematical notations can quickly become cumbersome, imprecise, or both. Here, we try to give a brief explanation of the symbols used in the present derivations. In the Nomenclature section, the symbols used to denote different sets of orbitals are listed. Note that virtual orbitals belonging to a specific orbital or pair domain are usually denoted with corresponding lower indices, e.g., $\mu \u0303i$ or $\xe3ij$. However, the occupied orbital indices are often obvious from context or are added as upper indices of a matrix, e.g., $F\mu \u0303i$ or $T\xe3b\u0303ij$, and thus, the repetition is avoided to simplify the notation.

Be aware that despite the similar notation, the symbols $K\xe3b\u0303ij$, $T\xe3b\u0303ij$, $F\xe3b\u0303ij$, and $S\xe3b\u0303ij$ have rather different meanings: the first denotes the two-electron repulsion integral $(i\xe3ij\u2009|\u2009jb\u0303ij)$, the second parameterizes excitations from orbitals *i* and *j* into PNOs $\xe3$ and $b\u0303$, while the last two are simply the Fock and overlap matrix elements in the PNO basis of pair *ij*.

When expressing derivatives in matrix form, we will use the following:

In the case of magnetic properties, special attention must be paid to complex conjugation: we use the superscripts “*,” “T,” and “†” to denote the complex conjugate, matrix transpose, and complex transpose of a quantity. A matrix trace is denoted as $trA=\u2211iAii$.

### B. Overview of DLPNO-MP2

Here, we provide a brief review of the closed-shell (spin-restricted) DLPNO-MP2 method, which has been described in significant detail elsewhere.^{19} The following steps are taken to obtain the DLPNO-MP2 energy:

The self-consistent field (SCF) equations are solved, giving the HF energy,

*E*_{HF}, and canonical molecular orbitals (CMOs).The occupied orbitals are localized according to the Foster–Boys (FB) criterion, as modified by Boys (referred to below as FB localization).

^{95,96}Core and valence orbitals are localized separately.PAOs are obtained as

and then normalized using the factor $N\mu \u0303\u2032=\mu \u0303u\u2032\u2009|\u2009\mu \u0303u\u2032\u221212$,

For each occupied orbital, a correlation domain of PAOs is selected based on the differential overlap criterion,

with the default value of the parameter *T*_{CutDOPre} = 0.03. These domains are then extended to all PAOs coming from the included atoms.

Non-redundant orthonormal PAOs (NPAOs) are constructed by diagonalizing the overlap matrix of each PAO domain, discarding eigenvectors with eigenvalues below

*T*_{S}(10^{−8}by default) and normalizing the rest. These are then transformed into pseudo-canonical NPAOs ${\mu \u0303}$ by diagonalizing the NPAO Fock matrix.Strongly correlated electron pairs are determined using a three-step procedure. A pair

*ij*is kept if DOI_{ij}is greater than*T*_{CutDOij}= 10^{−5}and the semi-canonical pair energy estimate [Eq. (8)] using the collinear dipole approximation to the exchange integrals^{97}[Eq. (8) computed with $M\mu \u0303\nu \u0303ij$ from Eq. (10) instead of Eq. (9)] is greater than*T*_{CutPre}= 10^{−6}. The energy contribution from the screened-out pairs, Δ*E*_{Pre}, is estimated using the more accurate dipole approximation to the exchange integrals^{11}[Eq. (9)],

where $F\mu \u0303i$ and $F\nu \u0303j$ are NPAO orbital energies in the domain of orbitals *i* and *j*, respectively, while *F*_{ii} and *F*_{jj} are elements of the Fock matrix in the localized MO (LMO) basis. The operator **r** in Eq. (12) is the electronic coordinate with respect to the global origin.

The correlation domains are reconstructed as in step 4 but with a tighter threshold

*T*_{CutDO}= 0.01. Pair domains are then made as the union of orbital domains. The pair domain PAOs are later orthonormalized as in step 5.PNOs are constructed by diagonalizing the semi-canonical pair density $D\u2323ij$ in the basis of NPAOs for the given orbital pair,

where $F\mu \u0303ij$ and $F\nu \u0303ij$ are NPAO orbital energies in the *ij* pair domain. PNOs with an occupation number greater than *T*_{CutPNO} = 10^{−8} are kept, while the rest form the “complementary” PNO (CPNO) basis referred to below. A tighter threshold, *T*_{CutPNO} × *T*_{ScalePNOCore} = 10^{−10}, is used for pairs that include core orbitals.^{98} The PAO to PNO transformation matrices are stored on a disk. A correction to the energy, Δ*E*_{PNO}, is computed as the difference between the semi-canonical correlation energy, calculated in the NPAO and PNO basis sets.

The full amplitude equations $\u2202E2DLPNO/\u2202T\u0303\xe3b\u0303ij*=0$ (see Sec. II C for definitions) are solved iteratively in the PNO basis, and the final DLPNO-MP2 energy is calculated,

### C. DLPNO-MP2 Lagrangian

Because the DLPNO-MP2 energy expression is non-variational, a Lagrangian formulation is used for analytic derivatives,^{99} which is slightly different from the one used previously,^{68,69}

The first four terms are the energy contributions, the next five are the Brillouin, localization, core–valence (CV), semicanonical, and PNO Lagrangian constraints, and the last three are orthonormality conditions for the MOs, PNOs, and NPAOs. The various terms will be explained as follows. *E*_{HF} is the HF energy expression,

with two-electron integrals in the Mulliken (1^{*}1|2^{*}2) notation. $E2DLPNO$ is the PNO-basis Hylleraas functional in the spin-adapted form^{100}

Δ*E*_{PNO} is an estimate for the error due to the PNO truncation, calculated as the difference between the correlation energy computed from the semicanonical amplitudes in the full NPAO basis and in the truncated PNO basis,

The energy contribution of the screened-out pairs takes the form of an approximate Hylleraas functional,

where 1_{S}(*ij*) = 1 if pair *ij* was screened-out, 0 otherwise. The dipole approximation to the exchange integrals is given in Eq. (9). The Brillouin condition is formulated as

where *z*_{ai} is the required Lagrange parameter. Introducing a magnetic field means that all quantities throughout the formalism are necessarily complex. Therefore, for any Lagrangian condition, it is necessary to constrain both its real and imaginary part or, equivalently, the condition and its complex conjugate. We favor the latter approach here, which also leads to complex Lagrange multipliers and their conjugates as separate optimization parameters, which can be arranged in Hermitian matrices, e.g., $zia\u2261zai*$. For infinitesimal perturbations, the final working equations in first order of the perturbation are always fully real (e.g., for electric fields) or fully imaginary (e.g., for magnetic fields) and are therefore over-parameterized. Thus, linear combinations are formed to obtain separate equations for the real and imaginary components of the Lagrange multipliers, of which one provides the required solution and the other trivially gives zero. The factors of one half in Eq. (27) as well as Eqs. (30), (32), (37), and (41) (see below) are introduced for consistency with the definitions of the Lagrange multipliers in Ref. 69. We also note in passing that Eq. (17) is formulated such that all terms are real.

In the present derivation, the molecular orbital (MO) response to perturbations is parameterized as^{101}

where **c** denotes the LMO coefficients. In contrast to the exponential parameterization used in Refs. 68 and 69, **U** is not constrained to be unitary but is instead constrained by the MO orthonormality condition,

which is fulfilled by the CMO coefficients and the unperturbed LMO coefficients **c**^{(0)}. Here, *δ*_{pq} is the Kronecker delta, and *x*_{pq} are the associated Lagrange multipliers. The term *C*_{Loc} enforces the FB localization condition separately on the core and valence orbitals and introduces the required Lagrange multipliers *l*_{ij},

The final MO constraint is the core–valence separation condition with an associated Lagrange parameter $z\xafim$,

Recall that the core and valence orbitals are localized separately because different PNO truncation thresholds are used for pairs that include core orbitals.^{98} Therefore, the core–valence separation condition is necessary even if all electrons are included in the correlation treatment. On the other hand, the localization condition for core orbitals is not needed if the latter are not localized, which is the case in the frozen-core approximation.

We denote the PNOs by $\xe3,b\u0303,c\u0303,d\u0303$, the CPNOs by $\xe3\u2032,b\u0303\u2032,c\u0303\u2032,d\u0303\u2032$, and the union of the two spaces by $\xe3\u2033,b\u0303\u2033,c\u0303\u2033,d\u0303\u2033$. The transformation matrices from redundant PAOs to each of these sets (for orbital pair *ij*) are denoted by $d\mu \u0303\u2032\xe3ij$, $d\mu \u0303\u2032\xe3\u2032\u2032ij$, and $d\mu \u0303\u2032\xe3\u2033\u2032\u2032ij$, respectively. The PNOs and CPNOs combined span the same space as the NPAOs. In the presence of a perturbation, we allow the PNOs to relax within the entire PNO + CPNO space. Rotations within each domain are thus parameterized as

Analogous to **U**, the matrices **θ**^{ij} are not necessarily unitary but rather are constrained by the PNO (and CPNO) orthonormality condition,

Note that **S**^{ij} is expanded as

A further PNO constraint ensures that the semi-canonical pair densities $D\u2323ij$ remain block-diagonal,

where $v\xe3b\u0303\u2032ij$ are unknown Lagrange multipliers. An orthonormality condition also applies to the orbital domain NPAOs used for the prescreening correction, parameterized as

where $\pi \mu \u0303\u2032\nu \u0303i$ are the transformation matrices from redundant to non-redundant orbital domain PAOs, **S**^{i} is the overlap matrix in the latter basis, and $x\mu \u0303\nu \u0303i$ are the necessary Lagrange multipliers. Finally, the semi-canonical amplitudes are parameterized as

and the semi-canonical residual conditions *C*_{SC} are included in the Lagrangian,

where $w\xe3\u2033b\u0303\u2033ij$ are Lagrange multipliers and the semi-canonical residuals $R\u2323\xe3\u2033b\u0303\u2033ij$ are defined as

The optimization parameters in the present DLPNO-MP2 Lagrangian are thus *U*_{pq}, $T\xe3b\u0303ij$, $T\u2322\xe3\u2033b\u0303\u2033ij$, $\theta \xe3\u2033b\u0303\u2033ij$, $\theta \mu \u0303\nu \u0303i$, $t\mu \u0303\nu \u0303ij$, *z*_{ai}, $z\xafim$, *l*_{ij}, *x*_{pq}, $v\xe3b\u0303\u2032ij$, $x\xe3\u2033b\u0303\u2033ij$, $x\mu \u0303\nu \u0303i$, $w\xe3\u2033b\u0303\u2033ij$, and their complex conjugates.

### D. Equations for first derivatives

The derivative of the DLPNO-MP2 energy with respect to the first external perturbation **κ** (in this work either an electric field or nuclear magnetic moment, on which only the one-electron part of the Hamiltonian depends) is equal to the respective partial derivative of the Lagrangian,

if and only if the Lagrangian is made stationary with respect to all parameters. The equations for the necessary stationarity conditions will be presented in this section. Despite the slightly different formulation of the Lagrangian, the final working equations are completely equivalent to those derived in Ref. 69. Note that Eq. (43) is only valid if the basis functions do not depend on the perturbation **κ**. Therein, **P**^{SCF} is the SCF density matrix,

while **D** is the orbital-relaxed DLPNO-MP2 difference density matrix,

where, as stated above, only the occupied–virtual and virtual–occupied blocks of **z** are non-zero with $zia\u2261zai*$. Analogously, $z\xaf$ only has non-zero core–valence blocks with $z\xafmi\u2261z\xafim*$. **D**′ is the “PNO-relaxed” difference density matrix,

The prescreening contributions are defined as

The values of the required variables are obtained by solving the stationarity conditions equations (and their complex conjugates),

The first eight are already fulfilled at the end of the DLPNO-MP2 energy calculation. The PNO amplitudes are obtained by iteratively solving the amplitude equations,

Note that the equations obtained from the derivative with respect to the contravariant amplitudes, $\u2202L/\u2202T\u0303ij$, are equivalent to a linear combination of $\u2202L/\u2202Tij$ and $\u2202L/\u2202Tji$ but produce a more convenient expression. As explained in Refs. 5 and 19, terms in the sum over *k* in Eq. (53) are skipped if the corresponding element of the Fock matrix has an absolute value below a threshold *F*_{Cut} = 10^{−5}.

The semi-canonical amplitude stationarity conditions take the form

The matrix $N\xe3\u2033b\u0303\u2033ij=\delta c\u0303\xe3\u2033\delta c\u0303b\u0303\u2033$ is introduced to “select” the truncated PNO basis. Equation (54) can be transformed to the pseudo-canonical NPAO basis using the unitary matrix $S\xe3\u2033\mu \u0303ij(0)$ and solved for **w**^{ij} at **λ** = 0 (i.e., **θ**^{ij} = **I**),

with

Equations for **v**^{ij} are obtained from the Lagrangian stationarity conditions with respect to $\theta ab\u2032ij$,

These can be solved in the PNO + CPNO basis, which diagonalizes $D\u2323ij$. A detailed derivation, as well as expressions for **x**^{ij}, is given in the Appendix. The remaining Lagrange multipliers **z**, $z\xaf$, and **x** are determined from the MO rotation stationarity conditions. As expected, the Lagrangian is invariant to rotations among virtual orbitals, but occupied–virtual and occupied–occupied (core–core or valence–valence) rotations lead to the so-called Z-vector coupled-perturbed self-consistent field (Z-CPSCF)^{102} and Z-vector coupled-perturbed localization (Z-CPL)^{45} equations for **z** and **l**, respectively, while rotations between core and valence orbitals result in the Z-vector core–valence (Z-CV) equations for $z\xaf$. For rotations within a block of localized orbitals (either core or valence), we obtain the Z-CPL equations

where *i*, *j*, *k* are all within the same (core or valence) block. The localization response operator is defined as follows (note the change in sign and indexing compared to Ref. 69):

with

$\Gamma i\mu \u0303\u2032K$ is the two-body density matrix

$LijPre$ contains contributions from the approximate exchange integrals in the pre-screening correction,

where

The Z-CPL equation (59) must be solved iteratively. The stationarity condition with respect to rotations between core and valence MOs gives the Z-CV equations,

These are solved in the CMO basis. Note that in frozen-core calculations, the terms involving $\Gamma m\mu \u0303\u2032K$, $ym\mu \u0303\u2032\alpha ,x$, and $ym\beta ,x$ drop out.

Finally, the Z-CPSCF equations are obtained from the occupied–virtual block of Eq. (A22),

with the two-electron operator defined as

and likewise for matrices in the AO basis. The right-hand side of Eq. (69) contains several terms,

Note that the sum in Eq. (72) assumes a one-to-one correspondence between PAOs and their parent AOs and that $Sa\mu \u0303\u2032=Sa\mu $. The intermediates **Y**^{α}, **Y**^{β}, and **Y**^{γ} are defined as follows:

As discussed in the Appendix, **x**^{ij} is already Hermitian if Eqs. (56) and (58) are fulfilled, so the first two terms in Eq. (76) are actually equivalent. The same applies to the multipliers **x**^{i}, which are obtained analogously by solving $\u2202L/\u2202\theta i=\u2202L/\u2202\theta i\u2020=0$. We note that the first term in Eq. (79) is missing from Eq. (56) of Ref. 69, which appears to be an error in the publication as the term was, in fact, implemented in the DLPNO-MP2 gradient code. This is the final major step of the DLPNO-MP2 density calculation, after which **D** can be completed.

### E. Equations for second derivatives

The second derivative of the DLPNO-MP2 energy with respect to perturbations **κ** and **λ**, the latter being an electric or magnetic field in this work, is equal to the derivative of Eq. (43),

Apart from the second derivative integrals **h**^{κ,λ} (which are zero for the polarizability, i.e., when **κ** = **λ** = **F**^{ext}), evaluating this expression requires derivatives of the Lagrange parameters and unknown multipliers with respect to **λ**, which are obtained by taking derivatives of the constraints and stationarity conditions discussed in Sec. II D. Note that all **λ**-derivatives in this section are real for electric perturbations and imaginary for magnetic ones. The SCF response density **P**^{SCF,λ} requires the MO response coefficients **U**^{λ},

where in the second line, we have used the perturbed MO orthonormality constraint

We note in passing that the AOs depend on the external magnetic field when GIAOs are used, but not on the external electric field; thus, $S\mu \nu Fext=0$. $Uai\lambda $ are obtained from the perturbed Brillouin condition, which in conjunction with Eq. (82) gives the CPSCF equations,^{101}

where **F**^{(λ)} excludes derivatives of the MO coefficients. In the electric (magnetic) field case, the right-hand side is real (imaginary) and therefore so are $Uai\lambda $. The equations are solved in the CMO basis as usual, and the solution is transformed to the LMO basis. The remaining blocks of **U**^{λ}, which are needed below, are obtained from the perturbed localization, core–valence separation, and MO orthonormality constraints as follows. The core–valence block is obtained from the condition $Fim\lambda =0$ [in conjunction with Eq. (82)],

which can be solved in the CMO basis. The valence–valence block of **U**^{λ} is obtained from the perturbed localization condition $sij\lambda =0$, which gives the CPL equations,^{45}

where the plus/minus signs are used in the magnetic/electric field case, respectively, and

**s**^{(λ)} excludes derivatives of the MO coefficients and is only non-zero in the GIAO case, where it contains the perturbed dipole integrals,

which can be calculated as linear combinations of quadrupole integrals. If all electrons are included in the MP2 treatment, an analogous set of CPL equations must be solved for the core–core block, which can be obtained from Eq. (85) by substituting core for valence indices and vice versa. In the frozen-core approximation, the core–core block of **U**^{λ} is only constrained by Eq. (82), which can be solved as

The same solution is used for the virtual–virtual block of **U**^{λ} . Note that the left-hand side of CPL equations in the electric field case is equivalent to that of the Z-CPL equations (59), which can be written as

and $Aij,klloc,\xb1=Akl,ijloc,\xb1$ if *s*_{ij} = 0. Therefore, the same solver can be used for both equations (with some modifications for magnetic perturbations).

The expression for the DLPNO-MP2 response difference density **D**^{λ} is obtained by straightforward differentiation of Eqs. (45)–(51) and requires the perturbed Lagrange parameters **U**^{λ}, **z**^{λ}, $z\xaf\lambda $, **T**^{ij,λ}, $T\u2322ij,\lambda $, **w**^{ij,λ}, **θ**^{ij,λ}, **t**^{ij,λ}, and **θ**^{i,λ}. The PAO coefficients $P\u0303$ also depend on the perturbation through **U**,

For magnetic perturbations, $N\mu \u0303\u2032\lambda =0$. The perturbed (C)PNO orthonormality condition gives a relationship analogous to Eq. (82),

with

Note that the perturbed PAO overlap $S\u0303\mu \u0303\u2032\nu \u0303\u2032\lambda $ includes the PAO response [Eq. (92)] and is non-zero even for electric perturbations, unlike the perturbed AO overlap. We use Eq. (95) to remove the dependence on **θ**^{ij,λ} in the perturbed SC residual equations,

where

Equation (96) is once again solved in the pseudo-canonical NPAO basis. Next, the perturbed PNO constraint gives the PNO–CPNO block of **θ**^{ij,λ},

These equations are solved in the basis that diagonalizes $D\u2323ij$. The PNO–PNO and CPNO–CPNO blocks of **θ**^{ij,λ} are only constrained by Eq. (95) and can be assigned “symmetrically” analogously to Eq. (90), although an alternative choice will be discussed in Sec. III B. The perturbed PNO amplitudes **T**^{ij,λ} are then obtained from the perturbed PNO residual equations,

which can be solved iteratively with the same solver as Eq. (53) The prescreening contributions to the response density, **D**^{o(Pre),λ} and **D**^{v(Pre),λ}, require the NPAO response coefficients **θ**^{i,λ}, which are only constrained by the perturbed NPAO orthonormality condition with the possible solution,

Also required are the perturbed prescreening amplitudes, obtained from the equations $dd\lambda \u2202L\u2202tij*=0$ giving

For magnetic perturbations, $Rij\lambda =0$. The perturbed Lagrange multipliers **w**^{ij,λ} are obtained from the equations

which are solved in the pseudo-canonical NPAO basis with

The equations for the perturbed multipliers **v**^{ij,λ} are derived in the Appendix and have the form

with

The perturbed localization condition multiplier **l**^{λ} is obtained from the perturbed Z-CPL equations,

where $Lij\lambda $ and $Rijloc,\lambda [l(0)]$ are straightforwardly derived from Eqs. (60) and (61). The derivatives of the prescreening contributions $yi\mu \u0303\u2032\alpha ,x$ and $yi\beta ,x$ are somewhat more involved, so we provide them here,

The perturbed Z-CV equations have the form

where $LmjCV,\lambda $ is the derivative of Eq. (68). Finally, the perturbed Z-CPSCF equations must be solved for **z**^{λ},

Thus, the full response density can be assembled and contracted with the property integrals to obtain the second derivative.

At the end of this section, we briefly summarize the order of computational steps needed to obtain the DLPNO-MP2 response density in Table I. The actual procedure used in our implementation is given in more detail in Sec. IV.

In . | . | Equations . | . | Out . |
---|---|---|---|---|

SCF, FB-loc. | → | c, E_{HF} | ||

c | → | LMP2 guess (13)–(15) | → | $T\u2323ij,d\u2032\u2032ij,\Delta EPNO$ |

d^{ij} | → | LMP2 iter. (53) | → | $Tij,E2DLPNO$ |

$T\u2323ij,d\u2032\u2032ij,Tij$ | → | PNO cond. (58) | → | v^{ij} |

$T\u2323ij,d\u2032\u2032ij,vij$ | → | SC cond. (56) | → | w^{ij} |

$T\u2323ij,d\u2032\u2032ij,Tij,wij$ | → | Equation (63) | → | Γ^{K} |

$c,T\u2323ij,d\u2032\u2032ij,Tij,wij$ | → | Equation (46) | → | D′ |

D′, Γ^{K} | → | Z-CPL (59) | → | l |

l, Γ^{K} | → | Z-CV (67) | → | $z\xaf$ |

$c,l,\Gamma K,D\u2032,z\xaf$ | → | Z-CPSCF (69) | → | z |

$c,D\u2032,z\xaf,z$ | → | Equation (45) | → | D |

c | → | CPSCF (83), CV (84), CPL (85), $CMOO\lambda $ (82) | → | U^{λ}, c^{λ} |

$c\lambda ,T\u2323ij,d\u2032\u2032ij$ | → | PNO, respectively, (95), (96), and (99) | → | $T\u2323ij,\lambda ,\theta ij,\lambda ,d\u2032\u2032ij,\lambda $ |

d^{ij,λ}, T^{ij} | → | Equation (100) | → | T^{ij,λ} |

$T\u2323ij,\lambda ,d\u2032\u2032ij,\lambda ,Tij,\lambda ,vij$ | → | Pert. PNO cond. (107) | → | v^{ij,λ} |

$T\u2323ij,\lambda ,d\u2032\u2032ij,\lambda ,vij,\lambda ,wij$ | → | Pert. SC cond. (105) | → | w^{ij,λ} |

$T\u2323ij,\lambda ,d\u2032\u2032ij,\lambda ,Tij,\lambda ,wij,\lambda $ | → | Deriv. of Eq. (63) | → | Γ^{K,λ} |

$c\lambda ,T\u2323ij,\lambda ,d\u2032\u2032ij,\lambda ,Tij,\lambda ,wij,\lambda $ | → | Deriv. of Eq. (46) | → | D′^{λ} |

D′^{λ}, Γ^{K,λ}, l | → | Pert. Z-CPL (109) | → | l^{λ} |

$l\lambda ,\Gamma K,\lambda ,z\xaf$ | → | Pert. Z-CV (112) | → | $z\xaf\lambda $ |

$c\lambda ,l\lambda ,\Gamma K,\lambda ,D\u2032\lambda ,z\xaf\lambda ,z$ | → | Pert. Z-CPSCF (113) | → | z^{λ} |

$c\lambda ,D\u2032\lambda ,z\xaf\lambda ,z\lambda $ | → | Deriv. of Eq. (45) | → | D^{λ} |

c^{λ}, D^{λ} | → | Property equations (1), (2), and (80) | → | α, σ^{A} |

In . | . | Equations . | . | Out . |
---|---|---|---|---|

SCF, FB-loc. | → | c, E_{HF} | ||

c | → | LMP2 guess (13)–(15) | → | $T\u2323ij,d\u2032\u2032ij,\Delta EPNO$ |

d^{ij} | → | LMP2 iter. (53) | → | $Tij,E2DLPNO$ |

$T\u2323ij,d\u2032\u2032ij,Tij$ | → | PNO cond. (58) | → | v^{ij} |

$T\u2323ij,d\u2032\u2032ij,vij$ | → | SC cond. (56) | → | w^{ij} |

$T\u2323ij,d\u2032\u2032ij,Tij,wij$ | → | Equation (63) | → | Γ^{K} |

$c,T\u2323ij,d\u2032\u2032ij,Tij,wij$ | → | Equation (46) | → | D′ |

D′, Γ^{K} | → | Z-CPL (59) | → | l |

l, Γ^{K} | → | Z-CV (67) | → | $z\xaf$ |

$c,l,\Gamma K,D\u2032,z\xaf$ | → | Z-CPSCF (69) | → | z |

$c,D\u2032,z\xaf,z$ | → | Equation (45) | → | D |

c | → | CPSCF (83), CV (84), CPL (85), $CMOO\lambda $ (82) | → | U^{λ}, c^{λ} |

$c\lambda ,T\u2323ij,d\u2032\u2032ij$ | → | PNO, respectively, (95), (96), and (99) | → | $T\u2323ij,\lambda ,\theta ij,\lambda ,d\u2032\u2032ij,\lambda $ |

d^{ij,λ}, T^{ij} | → | Equation (100) | → | T^{ij,λ} |

$T\u2323ij,\lambda ,d\u2032\u2032ij,\lambda ,Tij,\lambda ,vij$ | → | Pert. PNO cond. (107) | → | v^{ij,λ} |

$T\u2323ij,\lambda ,d\u2032\u2032ij,\lambda ,vij,\lambda ,wij$ | → | Pert. SC cond. (105) | → | w^{ij,λ} |

$T\u2323ij,\lambda ,d\u2032\u2032ij,\lambda ,Tij,\lambda ,wij,\lambda $ | → | Deriv. of Eq. (63) | → | Γ^{K,λ} |

$c\lambda ,T\u2323ij,\lambda ,d\u2032\u2032ij,\lambda ,Tij,\lambda ,wij,\lambda $ | → | Deriv. of Eq. (46) | → | D′^{λ} |

D′^{λ}, Γ^{K,λ}, l | → | Pert. Z-CPL (109) | → | l^{λ} |

$l\lambda ,\Gamma K,\lambda ,z\xaf$ | → | Pert. Z-CV (112) | → | $z\xaf\lambda $ |

$c\lambda ,l\lambda ,\Gamma K,\lambda ,D\u2032\lambda ,z\xaf\lambda ,z$ | → | Pert. Z-CPSCF (113) | → | z^{λ} |

$c\lambda ,D\u2032\lambda ,z\xaf\lambda ,z\lambda $ | → | Deriv. of Eq. (45) | → | D^{λ} |

c^{λ}, D^{λ} | → | Property equations (1), (2), and (80) | → | α, σ^{A} |

## III. TREATMENT OF NUMERICAL INSTABILITIES

### A. Localization response singularities

It was discussed at some length in Ref. 69 that for some systems, the Z-CPL (as well as the CPL and perturbed Z-CPL) equations can be (near-)singular. In other words, the left-hand side matrix **A**^{loc,±} can have (near-)zero eigenvalues. This happens when a continuous degeneracy exists among the LMOs, i.e., when some subset of the localized orbitals can be arbitrarily rotated without changing the value of the localization sum. Formally, continuous degeneracies are only possible for certain point groups;^{103–105} however, even local pseudosymmetry can result in small eigenvalues of **A**^{loc,±} and thus prevent convergence of the CPL/Z-CPL equations or negatively impact the results, as shown for the (CH_{3})_{3}P(CO)_{4}Os⋯Cr(CO)_{5} molecule (reference code KAMDOR)^{106} in Ref. 69. To systematically treat this problem, here we take the same approach as for the DLPNO-MP2 gradient, namely, to modify the Lagrangian so as to obtain a nonsingular problem with the nullspace projected out. Note that the matrix **A**^{loc,−}, which occurs in the Z-CPL equation (59) and in the CPL (85) and perturbed Z-CPL (109) equations for electric perturbations, is different from the matrix **A**^{loc,+}, which occurs in the CPL and perturbed Z-CPL equations for magnetic perturbations. Thus, the two matrices have different eigenvalues and eigenvectors. In order to ensure that both types of equations are properly treated, it is most convenient to reformulate the localization constraint in the Lagrangian using separate real and imaginary terms,

where the equivalence of the two formulations is demonstrated for the valence orbitals. The core orbital terms are analogous, so for simplicity, we will omit them in this section. Note that both of the newly introduced Lagrange multiplier matrices $lR$ and $lI$ are real by definition. Using this definition of *C*_{Loc}, Eq. (91) becomes

where

The real and imaginary parts of Eq. (116) can be separated, and the latter yields $lI=0$ because *L*_{ij} is real. Note that as the unperturbed quantities are purely real, $Aloc,R$ coincides with *A*^{loc,−} and $Aloc,I$ coincides with *A*^{loc,+}. We then define the eigendecompositions of these matrices as

where the index *n*″ denotes the eigenvectors of either matrix (determined from context), *n*′ denotes those with (near-)zero eigenvalues, and *n* denotes the rest. We change the Lagrangian constraint to

where *δ*_{ij,kl} = *δ*_{ik}*δ*_{jl} and the projectors onto the nullspace of $Aloc,R/I$ are

Thus, the Z-CPL equations become

where

Once again, $lI$ is zero, while the contributions of $lR$ to Eqs. (68), (73), (112), and (113) must be calculated with the projected matrix,

The electric/magnetic CPL equations become

i.e., the nullspace must first be projected out of the right-hand side in Eq. (127), the equations are solved, and afterward, the nullspace is also projected out of the solution before the perturbed orthonormality constraint is used to obtain the lower triangle of the matrix,

The perturbed Z-CPL equations become

with

Equation (130) is separated into real and imaginary parts, of which only the former/latter yields a nonzero solution for electric/magnetic perturbations, respectively. Note that for magnetic perturbations, $rii\u2212rjj\lambda =0$ and $rij+rij*\lambda =0$. The contributions of $lR/I,\lambda $ to Eqs. (112) and (113) must be calculated with the projected matrices

Finally, we demonstrate the numerical instabilities that can occur if the localization response matrix, $Aloc,R$ or $Aloc,I$, is singular, as well as the effectiveness of the procedure to remove the near-zero eigenvalues. We chose the octahedral molecules SF_{6} and SeF_{6} with bond lengths 1.560 and 1.688 Å, respectively. SF_{6} has six eigenvalues of $Aloc,R$ smaller than 10^{−4} due to a continuous degeneracy in the FB LMOs among the 18 fluorine sp^{3}-hybridized orbitals (3 on each F atom) directed away from the S atom. Three more singular eigenvalues are associated with the sulfur 2sp shell, which is treated as core orbitals. In addition to these, SeF_{6} has ten more singular eigenvalues in the valence and three in the core region, associated with the 3d and 3sp shells, respectively. The magnetic response matrix $Aloc,I$ only has ten singular eigenvalues for SeF_{6}, due to the 3d orbitals, and none for SF_{6}. Table II shows that if no measures are taken to remove the singularities, the iterative solution of the perturbed CPL, Z-CPL, or PNO amplitude equations may not converge or large errors may occur in the calculated properties, which is arguably worse. However, if the singular eigenvectors are projected out, the iterative solutions are stable and the results are in line with the expected accuracy for DLPNO-MP2.

System . | SF_{6}
. | SeF_{6}
. | ||
---|---|---|---|---|

. | FC . | AE . | FC . | AE . |

$Aloc,R$ nullspace | 6 | 9 | 16 | 24 |

$Aloc,I$ nullspace | 0 | 0 | 10 | 10 |

Errors without projection: | ||||

|Δα| (Å^{3}) | 8 × 10^{5} | …^{a} | 17 | …^{a} |

max (|Δσ_{F}|) (ppm) | 0.26 | 0.28 | 2.1 | 29 |

Errors with projection: | ||||

|Δα| (Å^{3}) | <0.001 | <0.001 | <0.001 | <0.001 |

max (|Δσ_{F}|) (ppm) | 0.26 | 0.28 | 0.47 | 0.55 |

System . | SF_{6}
. | SeF_{6}
. | ||
---|---|---|---|---|

. | FC . | AE . | FC . | AE . |

$Aloc,R$ nullspace | 6 | 9 | 16 | 24 |

$Aloc,I$ nullspace | 0 | 0 | 10 | 10 |

Errors without projection: | ||||

|Δα| (Å^{3}) | 8 × 10^{5} | …^{a} | 17 | …^{a} |

max (|Δσ_{F}|) (ppm) | 0.26 | 0.28 | 2.1 | 29 |

Errors with projection: | ||||

|Δα| (Å^{3}) | <0.001 | <0.001 | <0.001 | <0.001 |

max (|Δσ_{F}|) (ppm) | 0.26 | 0.28 | 0.47 | 0.55 |

^{a}

The calculation did not converge.

### B. PAO domain redundancy

While testing our initial implementation, we discovered large errors due to numerical instabilities related to the domain truncation. After thorough investigation, we found that the source was the appearance of near linear dependencies in the “non-redundant” PAO domains. Recall that the full PAO space is linearly dependent because it is of the same size as the full AO basis, *N*_{AO}, but it only represents the *N*_{virt} virtual orbitals. Therefore, there are exactly *N*_{occ} redundant PAOs. However, any subset (domain) of PAOs is not *exactly* linearly dependent, which becomes obvious when inspecting the eigenvalue spectrum of the overlap matrix of a domain of PAOs—see Fig. 1. While the spectrum of the full overlap matrix clearly reveals the *N*_{occ} redundant eigenvectors, all smaller domains have virtually continuous eigenvalue spectra, i.e., the matrix is ill-conditioned. By default, eigenvectors with eigenvalues greater than *T*_{S} = 10^{−8} form the NPAO domain. It is apparent that the lowest remaining eigenvalue *s*_{n} will be very close to the threshold, and thus, the normalization coefficient of that eigenvector will be $sn\u221212\u2248104$. This is not a cause for concern, and indeed, the energy and gradient calculations do not suffer. It is, however, a problem for second derivatives due to the perturbed PNO/CPNO coefficients **d**″^{ij,λ}. To illustrate this, we express the latter as

where **U**^{S} is the unitary matrix that diagonalizes the PAO domain overlap matrix, **s** is the diagonal matrix of eigenvalues, and **U**^{SD} is a unitary transformation from the overlap eigenbasis to the SC pair density eigenbasis. We will drop the *ij* indices in the rest of this section for simplicity of notation. Let us also assume for now that no PNO truncation is used; thus, using Eq. (95), we can express

and therefore,

The largest elements of the matrices are estimated based on $s\u221212\u2248TS\u221212=104$. These large numbers lead to a loss of numerical precision, which gets amplified during the multiple iterative procedures, especially the solution of the CPSCF equations. An example of this problem in its purest form is illustrated in Table III for a toy model of two water molecules in a triangular prism configuration with C_{2v} symmetry at a distance of 10.5 bohrs between the two molecular planes. No truncations were used other than the domain truncation with *T*_{CutDO} = 10^{−2}. Ignoring the line “Analytic v2” for now, the large error in the analytic polarizability at *T*_{S} = 10^{−8} is apparent. This disappears at higher values of the cutoff, and at *T*_{S} = 10^{−5}, there is very good agreement with both the semi-numeric polarizability (calculated from analytic dipole moment calculations with an explicit electric field of 10^{−4} a.u.) and the RI-MP2 value. Note that even the errors in the relaxed MP2 density get larger with lower values of *T*_{S}, hinting at the underlying numerical instability, although, of course, these errors are negligible.

. | T_{S} = 10^{−8}
. | T_{S} = 10^{−7}
. | T_{S} = 10^{−6}
. | T_{S} = 10^{−5}
. | RI-MP2 . |
---|---|---|---|---|---|

Numeric | 19.1393 | 19.0770 | 19.1034 | 19.1146 | 19.1190 |

Analytic v1 | 17.5287 | 19.0675 | 19.0973 | 19.1149 | 19.1190 |

Analytic v2 | 19.1033 | 19.0688 | 19.0974 | 19.1149 | |

max(|ΔD_{μν}|) | 0.0117 | 0.0110 | 0.0055 | 0.0010 |

. | T_{S} = 10^{−8}
. | T_{S} = 10^{−7}
. | T_{S} = 10^{−6}
. | T_{S} = 10^{−5}
. | RI-MP2 . |
---|---|---|---|---|---|

Numeric | 19.1393 | 19.0770 | 19.1034 | 19.1146 | 19.1190 |

Analytic v1 | 17.5287 | 19.0675 | 19.0973 | 19.1149 | 19.1190 |

Analytic v2 | 19.1033 | 19.0688 | 19.0974 | 19.1149 | |

max(|ΔD_{μν}|) | 0.0117 | 0.0110 | 0.0055 | 0.0010 |

After some experimentation, we found an alternative solution for Eq. (95), which partially alleviates the problem and results in the line “Analytic v2” in Table III. First, we express **θ**^{λ} as

where both **θ**^{1,λ} and **θ**^{2,λ} are defined as skew-Hermitian matrices and the latter only has elements in the PNO–CPNO block. Thus, the first and last terms in Eq. (136) are used to ensure that the perturbed orthonormality and SC conditions [Eqs. (95) and (99)] are fulfilled, respectively, while the second term is arbitrary. We got the best results by choosing it so as to minimize the norm of the perturbed PNO/CPNO coefficients, excluding the **θ**^{2,λ} part,

which is solved in the overlap eigenbasis as

All of the above also applies to the orbital-specific PAO domains used for the prescreening contributions and their response parameterized by **θ**^{i,λ}. Of course, in that case, **θ**^{i,2,λ} = 0.

This approach basically serves to better distribute the large floating point numbers throughout **d**″^{λ} and seems to reduce the error in the calculated properties, as evidenced by the toy example in Table III. However, agreement with the numeric polarizability is still not perfect, and in fact, large, unpredictable errors still occur in calculations for realistic systems when the default value of *T*_{S} is used, as will be shown in Sec. V. Therefore, in addition to the “Analytic v2” ansatz presented above, we have settled on a cutoff of *T*_{S} = 10^{−5}, with which the issue no longer occurs. The reduced NPAO domains do lead to a slightly increased error in the correlation energy, but the difference is negligible compared to the deviations due to the domain and PNO approximations, which are themselves sufficiently small by design, so as to obtain results with chemical accuracy.

We note in passing that the entire problem disappears if a linearly independent set of virtual orbitals is used, such as FB-localized MOs. In fact, we implemented this option and confirmed that the instability does not occur. However, on the one hand, localizing virtual orbitals to a true maximum, rather than a saddle point, of the FB condition is very difficult and time-consuming, even with the augmented Hessian solver used in our implementation. On the other hand, much larger domains were necessary to obtain accurate results, corroborating recent observations.^{107,108} Therefore, we will not discuss this approach in the rest of this paper.

### C. PNO response singularities

To solve Eqs. (58), (99), and (107), which follow from the PNO constraint $D\u2323\xe3b\u0303\u2032ij$, it is necessary to divide by the difference $(n\xe3ij\u2212nb\u0303\u2032ij)$ between occupation numbers of PNOs and CPNOs. Evidently, if this difference approaches zero, the equations become near-singular, which may lead to numeric instability. Usually, the eigenvalue spectrum of $D\u2323ij$ is nearly continuous, and thus, both the smallest PNO and the largest CPNO occupation numbers are of the same order of magnitude as *T*_{CutPNO}, and it is not uncommon for the threshold to fall between two nearly degenerate PNOs. In our tests, we found rare cases, where this instability caused spurious, sometimes very large, errors in the calculated properties. In all of these cases, the smallest relative difference $(n\xe3ij\u2212nb\u0303\u2032ij)/n\xe3ij$ was smaller than 10^{−2}, which can be taken as a necessary, but not sufficient, condition for problems to occur. Clearly, it is necessary to treat this problem systematically, since it is not possible to know when the results will be, or indeed have been, affected, without reference calculations, e.g., at other values of *T*_{CutPNO}.

The approach we propose is to introduce a level shift *ε* in the denominator of Eqs. (58), (99), and (107), thus ensuring that $(n\xe3ij\u2212nb\u0303\u2032ij+\epsilon )>\epsilon $. To that end, we modify the PNO condition as follows:

In effect, we shift all kept PNO occupation numbers up by *ε*. This modification has no effect on the DLPNO-MP2 energy, as the occupation numbers are only used for PNO truncation. It does, however, lead to extra terms in the right-hand sides of Eqs. (99) and (107) (given below in the pair density-diagonalizing basis),

An appropriate value of the level shift parameter *ε* must be chosen: if it is too low, it will have no effect. However, if it is too large, it can introduce additional round-off errors. It makes sense to choose a value proportional to *T*_{CutPNO}, i.e., *ε* = *ε*_{scale} × *T*_{CutPNO}. After some experimentation (see below), we found that *ε*_{scale} = 0.1 provides stable results without introducing additional errors. Note that because the value of *T*_{CutPNO} is different for occupied pairs, which include core orbitals, the value of *ε* reflects that. In principle, it is also possible for the level shift to be pair-specific, e.g., only applied to “problematic” domains, but we did not find that to be necessary.

We demonstrate a rather extreme example of the instability discussed above in Fig. 2. In shielding calculations for caffeine with LoosePNO default thresholds, *T*_{S} = 10^{−5} and *T*_{CutPNO} = 10^{−7}, as well as with *T*_{CutPNO} = 5.96 × 10^{−7}, the maximal absolute errors with respect to RI-MP2 are huge (40.5 and 106.3 ppm, respectively). The level shift scheme with *ε*_{scale} in the range 0.05–0.1 is successful in resolving the problem and producing a smooth curve. A value of 0.01 is not sufficient to remove the singularity at *T*_{CutPNO} = 10^{−7}, while *ε*_{scale} = 1.0 results in a smooth curve but with an additional small deviation for all values of *T*_{CutPNO}. Thus, we consider *ε*_{scale} = 0.1 to be an appropriate value.

Initially, no adjustment or level shift was applied to the calculations in Secs. V B and V C. Consequently, we encountered several calculations with conspicuously large errors, which are discussed whenever they occur in Secs. V B and V C. In all of these cases, the level shift approach (with *ε*_{scale} = 0.1) was successful in reducing the errors to within the expected bounds.

## IV. IMPLEMENTATION

We have implemented the DLPNO-MP2 response density calculation in the ORCA package,^{109} closely following the algorithm for the DLPNO-MP2 gradient.^{69} The code is structured as follows:

First, the SCF equations are solved.

Then, the DLPNO-MP2 energy and relaxed density calculation is performed. The PNO coefficients

**d**^{ij}and amplitudes**T**^{ij}are kept on disk, as well as the list of active orbital pairs and the domain information.Next, the SCF-level property calculation is performed, and the CPSCF solution $Uai\lambda $ is stored.

The remaining blocks of

**U**^{λ}are calculated, according to Eqs. (82), (84), (85), and (90).**F**^{λ}and**c**^{λ}are completed, as well as $P\u0303\lambda $.The three-index integrals $(i\mu \u0303\u2032\u2009|\u2009K)$ and $(i\mu \u0303\u2032\u2009|\u2009K)\lambda $ are calculated and stored.

In a loop over orbital pairs, the CPNO coefficients

**d**′^{ij}are reconstructed, and $T\u2323ij$ and $T\u2322ij,\lambda $ are calculated according to Eqs. (14) and (96).**θ**^{ij,λ}are computed, and the perturbed PNO coefficients**d**^{ij,λ}are stored on disk.The perturbed PNO amplitude equations, [Eq. (100)] are solved iteratively, separately for each perturbation.

The contributions to the response density matrix from

*D*_{ij},**D**^{ij}, and their derivatives are calculated.The PNO-dependent terms are evaluated in a loop over orbital pairs:

- 9.1.
The full

**d**″^{ij}, $T\u2323ij$,**d**″^{ij,λ}, and $T\u2323ij,\lambda $ are reconstructed. - 9.2.
The intermediates

**τ**^{ij}and**τ**^{ij,λ}are calculated. - 9.3.
Equations (58) and (107) are solved to obtain

**v**^{ij}and**v**^{ij,λ}. - 9.4.
Equations (56) and (105) are solved to obtain

**w**^{ij}and**w**^{ij,λ}. Their contributions to**D**′^{o},**D**′^{ij}, and their derivatives are evaluated. - 9.5.
**G**^{i(ij)},**G**^{j(ij)},**G**^{i(ij),λ}, and**G**^{j(ij),λ}and their contributions to**Γ**^{K},**ξ**^{ij}, and**Y**^{γ}, and their derivatives are computed. - 9.6.
Terms related to

**ξ**^{ij}and**τ**^{ij}and their derivatives are added to**Y**^{δ}and**Y**^{δ,λ}.

- 9.1.
**Γ**^{K}and**Γ**^{K,λ}are read from disk, and their contributions to $Lij\lambda $, $Lmj\lambda $,**Y**^{α,λ},**Y**^{β}, and**Y**^{β,λ}are evaluated in the AO basis by generating the required three-index integrals on the fly.The prescreening contributions to $Lij\lambda $, $Lmj\lambda $, $Lia\lambda $, and

**D**′ are calculated.$Lmj\lambda $ is completed, and the Z-CV equations (112) are solved to obtain $zmj\lambda $.

$Lij\lambda $ is completed, and the Z-CPL equations (109) are solved to obtain $zij\lambda $. In our implementation, a conjugate gradient solver is used.

$Lia\lambda $ is completed, and the Z-CPSCF equations (69) are solved to obtain $zia\lambda $. This is most conveniently done in the CMO basis using a standard CPSCF solver.

The full response density

**D**^{λ}is assembled, and the DLPNO-MP2-level properties are calculated.

As the goal here is the calculation of field-response properties, the algorithm has been optimized for a small number of perturbations. The loop over the three Cartesian components of the external field is the innermost one in most cases such that the unperturbed and a single perturbed quantity are kept in memory at any given time. In some situations, e.g., step 9.2, it is faster to calculate all perturbed quantities at the same time, at the cost of higher memory usage. Likewise, all intermediates that are stored on a disk—such as $(i\mu \u0303\u2032\u2009|\u2009K)\lambda $, **d**^{ij,λ}, and **T**^{ij,λ}—are kept for all perturbations simultaneously. For these reasons, the implemented algorithm is not suitable for perturbations whose number grows with the system size, e.g., nuclear Hessians. It is also worth pointing out that the NPAO domains and the SC amplitudes are calculated four times in total—in steps 2 (twice), 6, and 9.1—while their response is calculated twice—in steps 6, and 9.1. This is not an insignificant computational overhead, as shown in Sec. V D, but is done intentionally to avoid the storage and I/O of quantities in the full NPAO domains, which are much larger than the PNO domains.

The asymptotic scaling behavior of the DLPNO-MP2 gradient algorithm was discussed in Ref. 69, and that discussion applies here as well. Briefly, the number of orbital pairs that survive the dipole-based prescreening scales linearly with the system size in the asymptotic limit, while the PAO correlation domains reach a finite size for large systems. This allows for asymptotic linear scaling of most major steps in the algorithm, provided sparsity is properly exploited using the sparse map infrastructure introduced in Ref. 19. In particular, the RI integral transformation, PNO generation, and iterative solution of the amplitude equations are asymptotically linear scaling. Two-electron contributions coming from the Fock matrix, e.g., in the CPSCF and Z-CPSCF equations, can be calculated with effective quadratic scaling using the RIJCOSX (RI for Coulomb and chain-of-spheres exchange) approximation.^{110} The calculation of terms involving (*K*|*ip*) integrals in the Z-CPSCF equations’ [Eq. (74)] right-hand sides asymptotically scales as $O(N2)$. Together with the $(i\mu \u0303\u2032\u2009|\u2009K)$ contributions [Eq. (75)], these terms take about 10%–15% of the computation time for the largest systems discussed here. Contributions from the screened-out pairs scale quadratically, while localization, CPL, and Z-CPL equations scale as $O(N3)$, but these steps constitute a very small part of the overall computation time.

## V. RESULTS AND DISCUSSION

### A. Computational details

To examine how quickly DLPNO-MP2 response properties converge toward the RI-MP2 results as the various thresholds are tightened, we performed a series of polarizability and NMR shielding calculations on the molecules depicted in Fig. 3.

For the shielding calculations, the pcSseg-2^{113} basis set was used, along with the cc-pwCVTZ/C^{114–116} auxiliary basis set, which was shown to be an appropriate combination in our previous study.^{72} RI was used only for Coulomb integrals in all Fock matrix builds with no approximation for the exchange integrals, while RIJK was used for the perturbed two-electron integrals over GIAOs. In both cases, the def2/JK^{117} auxiliary basis set was employed. All orbitals were included in the MP2 treatment.

For the polarizability calculations, we used the def2-TZVPD^{118,119} basis set together with the aug-cc-pVTZ/C^{114} auxiliary basis set. Once again, def2/JK^{117} was used for Coulomb fitting, and all electrons were correlated.

In order to isolate the influences of the different DLPNO parameters, we performed several sets of calculations as follows:

Varying

*T*_{CutPNO}for the valence orbitals with very conservative values of*T*_{CutPNO}= 10^{−12}for the core orbitals. We set*F*_{Cut}= 10^{−8},*T*_{S}= 10^{−8}, and all other DLPNO thresholds to 0.Varying

*T*_{CutPNO}for the core orbitals with*T*_{CutPNO}= 10^{−10}for the valence orbitals. We set*F*_{Cut}= 10^{−8},*T*_{S}= 10^{−8}, and all other DLPNO thresholds equal to 0.Varying

*F*_{Cut}with*T*_{CutPNO}= 10^{−12}for the core orbitals and*T*_{CutPNO}= 10^{−10}for the valence orbitals. We set*T*_{S}= 10^{−8}and all other DLPNO thresholds equal to 0.Same as above but also fixing

*F*_{Cut}= 10^{−5}and varying*T*_{CutDO}. Two values were used for*T*_{S}to examine the numerical stability: 10^{−8}and 10^{−5}.Same as above but also fixing

*T*_{CutDO}= 10^{−3},*T*_{CutDOPre}= 0.03, and*T*_{S}= 10^{−5}and varying*T*_{CutPre}. We set*T*_{CutDOij}= 0.05, a large value, so as to allow as many pairs as possible to be screened out without the dipole approximation breaking down.Setting all thresholds to the default (“NormalPNO”) values used for energy and gradient calculations:

*T*_{CutPNO}= 10^{−8}(10^{−10}for core orbitals),*F*_{Cut}= 10^{−5},*T*_{CutDO}= 10^{−2},*T*_{CutDOPre}= 0.03,*T*_{CutDOij}= 10^{−5},*T*_{CutPre}= 10^{−6},*T*_{CutMKN}= 10^{−3},*T*_{CutC}= 10^{−3}, but with*T*_{S}= 10^{−5}and a PNO level shift*ε*_{scale}= 0.1.Same as above with the “LoosePNO” default thresholds:

*T*_{CutPNO}= 10^{−7}(10^{−9}for core orbitals) and*T*_{CutDO}= 2 × 10^{−2}.We also performed shielding calculations with RI-MP2, NormalPNO, LoosePNO, and TightPNO (

*T*_{CutPNO}= 10^{−9}for valence and 10^{−11}for core orbitals,*T*_{CutDO}= 5 × 10^{−3}), together with the RIJCOSX approximation to see how the COSX and DLPNO errors stack up. The keyword DefGrid3 was used to choose finer COSX integration grids than the default. Note that the grids in ORCA were updated for version 5 of the program, and the DefGrid3 settings are similar to the “RIJCOSX-XL” settings used in Ref. 72.

To investigate the efficiency of the implementation for larger system sizes, we used idealized glycine chains in either a linear (*ϕ* = 180°, *ψ* = 180°) or *α*-helix (*ϕ* = −57°, *ψ* = −47°) conformation. The neutral form of the molecules was used, rather than the zwitterion. The structures were generated using Gabedit.^{120} As above, these calculations were performed with NormalPNO (and RI-MP2) and RIJCOSX. These calculations ran in a node-exclusive fashion on 8 Intel Xeon E5-2690v2@3.00GHz central processing unit (CPU) cores with 6 GB of random-access memory (RAM) per core.

To showcase all the capabilities of the implementation, another large shielding calculation was performed on the vancomycin molecule using the DSD-PBEP86 DHDFT functional,^{121} NormalPNO thresholds, RIJCOSX, frozen-core orbitals, and CPCM for implicit solvation in dimethylsulfoxide (DMSO). This calculation ran on 12 Intel Xeon E5-2687Wv4@3.00GHz CPU cores with 15 GB RAM per core.

Finally, a series of DLPNO-MP2 polarizability calculations were performed on [*n*]helicene (*n* = 7, 9, 11, 13, 15) molecules using the aug-cc-pVTZ,^{122,123} aug-cc-pVTZ/C, and def2/JK basis sets, NormalPNO thresholds, frozen-core orbitals, and RIJCOSX. The geometries were optimized with RI-MP2 (for *n* = 7–11) and DLPNO-MP2 (for *n* = 13–15) using the cc-pVTZ^{123,124} and cc-pVTZ/C^{114} basis sets. These calculations ran on 12 Intel Xeon E5-2687Wv3@3.10GHz CPU cores with 8 GB of RAM per core.

### B. Accuracy: NMR shieldings

In this section, we use some of the systems in Fig. 3 to examine how much the DLPNO-MP2 isotropic shielding constants deviate from the RI-MP2 values, as we tighten the DLPNO thresholds. The mean absolute errors (MAEs) vs RI-MP2 for different values of *T*_{CutPNO} are shown in Fig. 4(a) for the different nuclides in the systems ATP^{4−}, caffeine, ebselen, and penicillin. We note that the errors diminish rapidly with decreasing *T*_{CutPNO}, and at the default value of 10^{−8} for valence orbitals, they are much lower than the inherent error in MP2. For example, the MAE for C is 0.1 ppm, which is an order of magnitude smaller than the average deviation from CCSD(T).^{71,72} It is also noteworthy that the results are insensitive to core PNO truncation, except for the heavier elements, particularly P and Se. This may be, in part, because the pcSseg-2 basis set does not fully capture core correlation effects.^{72}

For simplicity in the discussion, we want a single measure of the DLPNO error; however, the shielding constants have very different magnitudes for different nuclides. Thus, we define an absolute weighted error (AWE), whereby we scale the error in the calculated isotropic shielding (e.g., vs RI-MP2) by a factor w^{A}, derived from the experimental range of chemical shifts for the nuclide *A*,

We use the values for $\delta minA$ and $\delta maxA$ shown in Table IV. A reasonable target accuracy is a mean AWE (MAWE) below 0.5%, which corresponds to 0.5 ppm for ^{13}C, 0.03 ppm for ^{1}H, and about 1 ppm for ^{15}N and ^{17}O. As can be seen from Fig. 4(b), the MAWE faithfully represents the trends in Fig. 4(a) while also allowing us to combine the data for different nuclides, as they now have the same order of magnitude. The same data are displayed in Fig. 5, separating the different systems and averaging over all nuclei. Here, we can see that ATP^{4−} is slightly more sensitive to core PNO truncation probably because it contains all the phosphorus atoms in the dataset. The same system also suffers from a PNO response instability at valence *T*_{CutPNO} = 10^{−9} with MaxAE = 1.51 ppm for ^{15}N vs 0.59 and 0.04 ppm for the same nucleus at *T*_{CutPNO} = 10^{−8} and *T*_{CutPNO} = 10^{−10}, respectively. With the level shift-based adjustment (see Sec. III C), this error drops to 0.17 ppm at valence *T*_{CutPNO} = 10^{−9}.

A
. | $\delta minA$ . | $\delta maxA$ . | w_{A}
. |
---|---|---|---|

^{1}H | −1 | 12 | 6.5 |

^{13}C | 0 | 200 | 100 |

^{15}N | 0 | 900 | 450 |

^{17}O | −40 | 1120 | 580 |

^{19}F | −300 | 400 | 350 |

^{33}S | −290 | 670 | 480 |

^{31}P | −180 | 250 | 215 |

^{77}Se | −1000 | 2000 | 1500 |

A
. | $\delta minA$ . | $\delta maxA$ . | w_{A}
. |
---|---|---|---|

^{1}H | −1 | 12 | 6.5 |

^{13}C | 0 | 200 | 100 |

^{15}N | 0 | 900 | 450 |

^{17}O | −40 | 1120 | 580 |

^{19}F | −300 | 400 | 350 |

^{33}S | −290 | 670 | 480 |

^{31}P | −180 | 250 | 215 |

^{77}Se | −1000 | 2000 | 1500 |

Setting *T*_{CutPNO} to a conservative value, we examine the effect of *F*_{Cut} in Fig. 6. Apparently, it imparts negligible errors on the shieldings, especially at the default value of 10^{−5}.

The influence of *T*_{CutDO} (at a conservative value of *T*_{CutPNO}) is shown in Fig. 7 for the same set of molecules. In the left subplot, the redundant PAO domains were truncated with the default eigenvalue threshold *T*_{S} = 10^{−8}, which leads to some numerical issues, as discussed in Sec. III B, seen in the figure as large kinks in the curves for some molecules. In the right subplot, a value of *T*_{S} = 10^{−5} was used, which leads to smooth convergence or the results toward the RI-MP2 reference. At the default value of *T*_{CutDO} = 10^{−2}, the MAWE over the whole dataset is under 0.05%.

Turning to the pair prescreening threshold *T*_{CutPre}, its effect on the shielding MAWE is shown in Fig. 8. Despite the relatively small systems, using a large *T*_{DOij} = 0.05 leads to a significant number of pairs to be screened-out: 50%–80% at *T*_{CutPre} = 10^{−4} and 20–50% at *T*_{CutPre} = 10^{−6}. Even so, the effect on the shieldings is minor. Neglecting the contributions from the prescreening correction Δ*E*_{Pre} also does not worsen the results significantly for *T*_{CutPre} ≤ 10^{−4}; however, calculating these contributions takes a very small proportion of the computation time, so there are no savings to be made by skipping them. On the other hand, the default threshold values of *T*_{DOij} = 10^{−5} and *T*_{CutPre} = 10^{−6} might be too conservative for NMR shielding calculations.

Having separately examined the effects of the most important DLPNO parameters on the shieldings, we now see what happens when all parameters are set to their default (NormalPNO) values. Figure 9 shows the error distributions for different nuclides over all eight test systems in Fig. 3. The NormalPNO/RIJONX distributions are fairly narrow, and both the medians and the ranges (which are more relevant for chemical shift calculations) are much smaller than the inherent error in MP2. There is an overall bias toward overestimation of the shielding constants, so some error compensation can be expected in chemical shift calculations. The same figure also contains results with the LoosePNO threshold settings, which produce both larger systematic deviations and wider error distributions, at the borderline of what we may consider tolerable.

All calculations discussed so far were performed using the RIJONX approximation. In practice, it is preferable to also approximate the exchange part of the Fock matrix, e.g., using the COSX approximation, as this significantly reduces the computation time. However, the semi-numerical integration used in COSX introduces additional errors and numerical noise to the calculation, which may get amplified in DLPNO calculations. Therefore, we have included in Fig. 9 the error distributions from calculations employing the RIJCOSX approximation with a rather large set of integration grids (corresponding to the DefGrid3 keyword). Although the DLPNO and COSX errors are basically cumulative, the latter are fairly small, so the conclusions from the previous paragraph still apply for the most part. One exception is the ^{1}H shieldings, where the combined NormalPNO and COSX errors lead to a maximum deviation of 0.2 ppm and a spread of almost 0.4 ppm, which is larger than the method error and could, in principle, lead to a wrong assignment, compared to RI-MP2. This can be remedied by using larger COSX grids, in particular, for the CPSCF and Z-CPSCF equations, which appear to be the main source of the deviations. As expected, tightening *T*_{CutPNO} and *T*_{CutDO} from LoosePNO to NormalPNO to TightPNO systematically reduces the errors toward the RI-MP2 reference.

As the effect of the local approximations becomes more significant with the increasing system size, the relative error (RE) in the properties can also be expected to increase, up to the point where all approximations are fully active: orbital domains have reached a constant size, etc. For an intrinsic (i.e., local) property such as NMR shielding, this means that absolute errors should also reach an upper limit for large enough systems. To examine this relationship, we use the calculations on the glycine chain systems, discussed in Sec. V D. Figure 10 shows the maximum absolute errors in the shieldings of different nuclides, as a function of the chain length. Apparently, the errors do increase slightly with the system size, and the effect is more pronounced for the less sparse *α*-helix conformation. However, the errors do not grow indefinitely, leveling off around (gly)_{10}, and are within the ranges shown in Fig. 9. A notable outlier is an oxygen nucleus in linear (gly)_{3} with a rather large error of 0.97 ppm for ^{17}O vs 0.25 and 0.41 ppm for (gly)_{2} and (gly)_{4}, respectively. Using the level shift-based adjustment, discussed in Sec. III C, the MaxAE for ^{17}O in linear (gly)_{3} is reduced to 0.42 ppm, consistent with the slightly oscillating trend.

### C. Accuracy: Dipole polarizabilities

As evident from Sec. II, the same general procedure can be used to calculate both magnetic shielding and dipole polarizability tensors. It is worth noting that the polarizability is an extensive global property of the system and as such is less useful to obtain for large molecules, which are the main targets of local correlation methods. Nonetheless, we have included the capability in our implementation and will evaluate the accuracy of the DLPNO-MP2 polarizabilities in this section. It is most convenient to use the relative error (RE) in the isotropic polarizability,

and its absolute value (ARE). This should be compared to the inherent error of MP2 vs more accurate methods such as CCSD(T), which is around 2% for closed-shell species.^{73} DHDFT performs slightly better but still over 1%. Therefore, we consider a DLPNO error of up to 0.5% to be reasonable.

Figure 11 shows the dependence of the error on *T*_{CutPNO} for several molecules. The DLPNO-MP2 polarizabilities correctly converge toward the RI-MP2 reference values and are accurate enough for *T*_{CutPNO} ≤ 10^{−7}. The influence of PNO truncation for core electrons is almost negligible, which is to be expected, since they do not contribute significantly to the polarizability, especially since we did not include core-polarization basis functions in these calculations. We also note a number of kinks in the curves on the left—one for coronene at *T*_{CutPNO} = 10^{−12} and two for ATP^{4−} at 10^{−7} and 10^{−10} (the latter also leads to the slightly larger error in the right subplot). We associate these discrepancies with numerical instabilities due to two factors. The first one is near linear dependencies in the AO basis: even though *T*_{S} was set to 10^{−5} for these calculations, the analogous threshold for the AO basis was 10^{−8} and the lowest overlap eigenvalues for caffeine, penicillin, coronene, and ATP^{4−} were 2.7 × 10^{−7}, 4.2 × 10^{−7}, 2.2 × 10^{−8}, and 4.6 × 10^{−7}, respectively. The second factor is small differences in the occupation numbers of PNOs and CPNOs, as discussed in Sec. V B. Indeed, either using a cutoff of 10^{−6} for the MO orthonormalization or applying the level shift, introduced in Sec. III C, both eliminate these discrepancies (see the supplementary material).

The effect of domain truncation on the polarizability is shown in Fig. 12. As for the shieldings, a low value of *T*_{S} results in erratic behavior and huge errors, especially for coronene, which already has linear dependency issues in the AO basis. At *T*_{S} = 10^{−5}, however, convergence toward the RI-MP2 reference is much smoother, and the errors are already sufficiently small at *T*_{CutDO} = 3 × 10^{−2}.

The influence of *F*_{Cut} on the polarizability, shown in Fig. 13, is rather small. Apparently, the default value of 10^{−5} is more than conservative enough.

Convergence with the pair prescreening threshold *T*_{CutPre}, as well as the significance of the terms derived from Δ*E*_{Pre}, is shown in Fig. 14. With the default parameters, only 18% of pairs are skipped for the largest system, ATP^{4−}; however with the high value of *T*_{CutDOij} used here, between 35% and 67% of pairs are screened out already at *T*_{CutPre} = 10^{−5}, and the errors in the polarizability are below 0.2%. As for the shieldings, the dipole-based correction terms only make a difference for very high values of the threshold. The somewhat larger errors for coronene at *T*_{CutPre} ≥ 10^{−4} are probably due to the strong delocalization in the system, together with the large value of *T*_{CutDOij}.

The combined effect of all DLPNO approximations with the default NormalPNO parameters (and *T*_{S} = 10^{−5}, *ε*_{scale} = 0) is shown in Fig. 15. The signed error is plotted to illustrate that DLPNO-MP2 consistently underestimates the RI-MP2 polarizabilities. In turn, MP2 tends toward overestimation with respect to CCSD(T),^{73} so some fortuitous error compensation can be expected with DLPNO-MP2. In any case, the DLPNO-MP2 values only deviate from the RI-MP2 reference by 0.1%–0.2%, which is well within our 0.5% tolerance.

### D. Computational efficiency

To compare the efficiency of our RI-MP2 and DLPNO-MP2 implementations, we performed calculations on a series of polyglycine chains in both linear and *α*-helical conformations. The latter are much denser and expected to be a bigger challenge for the local approximations. The total wall-clock times for the shielding calculations, including SCF and CPSCF solutions, are presented in Fig. 16. The steep scaling of the RI-MP2 calculations is apparent, and there is no difference between the two conformations, which is expected, as sparsity is not exploited apart from basic and rather conservative pre-screening. In contrast, both series of DLPNO calculations achieve quadratic effective scaling, although the prefactor for the helical conformation is about 4.5 times higher. The crossover point between RI- and DLPNO-MP2 is thus surpassed at (gly)_{5} (38 atoms, 160 electrons, and 994 basis functions) and (gly)_{9} (66 atoms, 280 electrons, and 1738 basis functions) in the linear and helical series, respectively. In other words, calculations that take more than a day or two (on eight cores) with RI-MP2 can be performed more efficiently with DLPNO-MP2. For the largest system reported here—(gly)_{15} (108 atoms, 460 electrons, and 3040 basis functions)—the difference is substantial: 1–3 days with DLPNO-MP2 (depending on the conformation) vs 36 days with RI-MP2. Before the crossover point, the DLPNO-MP2 calculation is at most two times slower [for the helical (gly)_{5}] than the RI-MP2 one. Another worthwhile comparison is given in Fig. 17, namely, how much more expensive a shielding calculation is at the MP2-level vs the SCF-level (HF or hybrid DFT). Apparently, with NormalPNO and RIJCOSX, the DLPNO-MP2 shieldings are obtained at only 5–20 times the cost of the preceding HF shielding calculation (the latter is included in the DLPNO-MP2 timings), and the ratio is still going down at (gly)_{15}. This is in stark contrast to RI-MP2, where the calculations with 400 correlated electrons are already more than 100 times more expensive than HF or hybrid DFT and hardly feasible beyond that. Based on the calculations performed for Fig. 9, the TightPNO (LoosePNO) settings increase (decrease) the total computation time by a factor of 1.3–2.2, compared to NormalPNO.

Detailed information on the most time-consuming steps for an even larger system—the vancomycin molecule—is given in Table V. In this case, a frozen-core calculation was performed with a DHDFT functional (DSD-PBEP86), also including implicit solvent DMSO effects using CPCM. Some key takeaways from the data are as follows: (1) The total time for the DLPNO-DHDFT calculation is only eight times longer than an equivalent hybrid DFT calculation (SCF plus SCF-level NMR timings). (2) The cost to calculate the DLPNO-MP2 response densities is roughly six times higher than the DLPNO-MP2 energy and density calculation, as for RI-MP2.^{72} (3) Reconstructing the SC amplitudes, NPAO space, and the perturbed (C)PNO coefficients twice takes about a third of the DLPNO-MP2 response time. (4) The equations related to orbital localization make up less than 1% of the computational time. (5) Timings for the evaluation of DFT functional, CPCM, and dipole-based prescreening terms are almost negligible and therefore not shown.

Atoms . | 176 . |
---|---|

Basis size (pcSseg-2) | 4 700 |

AuxJ size (def2/JK) | 9 097 |

AuxC size (cc-pwCVTZ/C) | 13 591 |

Core electrons | 218 |

Valence electrons | 542 |

Kept orbital pairs | 22 265 |

Mean PAO domain size | 1 013 |

Mean PNO domain size | 22 |

Mean AuxC domain size | 1 852 |

Total time | 4 791 |

SCF (15 cycles) | 355 |

SCF-level NMR | 242 |

→CPSCF right-hand size | 56 |

→CPSCF (8 cycles) | 175 |

DLPNO-MP2 energy + density | 579 |

$\u2192(i\mu \u0303\u2032\u2009|\u2009K)$ transformation | 23 |

$\u2192T\u2323ij$ and d^{ij} construction | 88 |

→PNO amplitude iterations | 35 |

→NPAO space reconstruction | 163 |

→PNO-specific terms | 70 |

$\u2192\Gamma i\mu \u0303\u2032K$ contributions | 84 |

→Fock response terms | 10 |

→Z-CPSCF (10 cycles) | 56 |

→ Localization + Z-CPL | 5 |

DLPNO-MP2 response density | 3 603 |

$\u2192(i\mu \u0303\u2032\u2009|\u2009K)\lambda $ transformation | 276 |

$\u2192T\u2322ij,\lambda $ and d″^{ij,λ} construction | 809 |

→T^{ij,λ} iterations | 160 |

→NPAO space reconstruction | 270 |

→PNO-specific terms | 671 |

$\u2192\Gamma i\mu \u0303\u2032K$ and $\Gamma i\mu \u0303\u2032K,\lambda $ contributions | 748 |

→Fock response terms | 127 |

→Perturbed Z-CPSCF (8 cycles) | 176 |

→CPL + perturbed Z-CPL | 27 |

Atoms . | 176 . |
---|---|

Basis size (pcSseg-2) | 4 700 |

AuxJ size (def2/JK) | 9 097 |

AuxC size (cc-pwCVTZ/C) | 13 591 |

Core electrons | 218 |

Valence electrons | 542 |

Kept orbital pairs | 22 265 |

Mean PAO domain size | 1 013 |

Mean PNO domain size | 22 |

Mean AuxC domain size | 1 852 |

Total time | 4 791 |

SCF (15 cycles) | 355 |

SCF-level NMR | 242 |

→CPSCF right-hand size | 56 |

→CPSCF (8 cycles) | 175 |

DLPNO-MP2 energy + density | 579 |

$\u2192(i\mu \u0303\u2032\u2009|\u2009K)$ transformation | 23 |

$\u2192T\u2323ij$ and d^{ij} construction | 88 |

→PNO amplitude iterations | 35 |

→NPAO space reconstruction | 163 |

→PNO-specific terms | 70 |

$\u2192\Gamma i\mu \u0303\u2032K$ contributions | 84 |

→Fock response terms | 10 |

→Z-CPSCF (10 cycles) | 56 |

→ Localization + Z-CPL | 5 |

DLPNO-MP2 response density | 3 603 |

$\u2192(i\mu \u0303\u2032\u2009|\u2009K)\lambda $ transformation | 276 |

$\u2192T\u2322ij,\lambda $ and d″^{ij,λ} construction | 809 |

→T^{ij,λ} iterations | 160 |

→NPAO space reconstruction | 270 |

→PNO-specific terms | 671 |

$\u2192\Gamma i\mu \u0303\u2032K$ and $\Gamma i\mu \u0303\u2032K,\lambda $ contributions | 748 |

→Fock response terms | 127 |

→Perturbed Z-CPSCF (8 cycles) | 176 |

→CPL + perturbed Z-CPL | 27 |

A final test is shown in Table VI. Here, DLPNO-MP2 polarizability calculations were carried out for a series of [*n*]helicene molecules. The results are compared to those obtained by Friese *et al.* with their Laplace-based RI-MP2 implementation.^{59} The calculated isotropic polarizabilities agree to within 0.3%, as expected. It is not fair to compare the timings directly due to the different computer architectures, but we can look at the trends and effective scaling of both implementations with the system size. Based on the calculations for *n* = 7–11, the RI-MP2 algorithm scales as $O(Nbas4.5)$, roughly equivalent to the $O(Nbas4.3)$ scaling of our DLPNO-MP2 code for the same systems. The poor performance of the DLPNO approximations here can be explained by the high degree of delocalization in the helicenes, combined with the diffuse basis sets, which leads to large correlation domains and to practically all occupied orbital pairs surviving the dipole-based prescreening. Regardless, the performance of DLPNO-MP2 improves for the larger systems, and the effective scaling for *n* = 11–15 drops to $O(Nbas2.8)$.

. | DLPNO-MP2 . | RI-MP2^{59}
. | |||
---|---|---|---|---|---|

N
. | N_{bas}
. | α_{iso}
. | t^{a}
. | α_{iso}
. | t^{b}
. |

7 | 1794 | 357.7 | 11 | 357 | 65 |

9 | 2254 | 436.4 | 30 | 437 | 134 |

11 | 2714 | 516.5 | 66 | 518 | 425 |

13 | 3174 | 596.5 | 104 | ||

15 | 3634 | 675.7 | 150 |

. | DLPNO-MP2 . | RI-MP2^{59}
. | |||
---|---|---|---|---|---|

N
. | N_{bas}
. | α_{iso}
. | t^{a}
. | α_{iso}
. | t^{b}
. |

7 | 1794 | 357.7 | 11 | 357 | 65 |

9 | 2254 | 436.4 | 30 | 437 | 134 |

11 | 2714 | 516.5 | 66 | 518 | 425 |

13 | 3174 | 596.5 | 104 | ||

15 | 3634 | 675.7 | 150 |

^{a}

On 12 Intel Xeon E5-2687Wv3@3.10GHz CPU cores with 8 GB of RAM per core.

^{b}

On 12 Intel Xeon X5670@2.93GHz CPU cores.

## VI. CONCLUSIONS AND OUTLOOK

In this work, we have presented a derivation of analytic second derivatives for the DLPNO-MP2 method, applicable to cases where the AOs are independent of at least one of the two perturbations. Our main focus was NMR shieldings, but we have also discussed dipole polarizabilities. An efficient implementation of the method in the ORCA package was reported and used to benchmark the accuracy of the DLPNO approximations on a test set of medium-sized molecules.

Errors in the calculated properties were found to be sufficiently small. Setting the local approximation thresholds to the default (“NormalPNO”) values used in energy and gradient calculations results in errors in the calculated shieldings an order of magnitude smaller than the inherent error of MP2. For example, the maximum absolute deviations from the RI-MP2 reference for ^{13}C and ^{1}H shieldings in our test set are 0.3 and 0.04 ppm, respectively. The DLPNO-MP2 isotropic polarizabilities with default thresholds were also within 0.2% of the RI-MP2 reference values.

Comparing the computational cost of DLPNO-MP2 and RI-MP2 property calculations for linear and helical polyglycine, we found that the DLPNO-MP2 implementation becomes more efficient for systems larger than 38 atoms/160 electrons and 66 atoms/280 electrons, respectively. Due to the asymptotic linear scaling of the major post-SCF steps, DLPNO-MP2 property calculations are never more than about 20 times more expensive than equivalent HF or hybrid DFT calculations. The largest DLPNO-MP2 NMR shielding calculation reported here—on vancomycin (176 atoms, 542 correlated electrons, and 4700 basis functions)—took 80 h on 12 CPU cores, compared to 10 h for hybrid DFT.

The reported DLPNO-MP2 second derivative implementation is valuable on its own as it allows accurate MP2 and DHDFT property calculations for virtually all systems that could be treated at the hybrid DFT level. However, the derivation presented here also serves as a stepping stone toward analytic second derivatives of higher-level local correlation methods, such as DLPNO-CCSD, which we intend to pursue in the future. It is important to note that MP2, DHDFT, and any approximate variants thereof can fail fundamentally for certain systems, such as those with a small HOMO–LUMO gap (e.g., transition metal complexes) or high degree of static correlation, whereas many of these cases can be treated accurately with coupled cluster methods.

The present work also opens the door to multi-layer schemes such as those reported in Ref. 125, in which some orbital pair interactions are treated with lower accuracy in the DLPNO approximations or completely neglected in the correlation treatment. Such an approach may be particularly well suited to NMR shielding calculations of solvated systems or cluster models of molecular crystals.

## SUPPLEMENTARY MATERIAL

See the supplementary material for Cartesian coordinates of all test systems (.xyz); calculated isotropic polarizabilities, isotropic shieldings, and shielding anisotropies (.csv); and example ORCA input files for DLPNO-MP2 NMR shielding and dipole polarizability calculations (.txt).

## DATA AVAILABILITY

The data that support the findings of this study are available within the article and its supplementary material. Additional data, such as complete ORCA output files, are available from the corresponding author upon reasonable request.

## ACKNOWLEDGMENTS

The authors would like to thank Frank Wennmohs and Ute Becker for their help in implementation and parallelization of the code. G.L.S. thanks Ida-Marie Høyvik and Fred Manby for helpful discussions during the 2019 Molecular Quantum Mechanics conference in Heidelberg and especially Peter Pinski for an in-depth introduction to DLPNO-MP2 and its derivatives. G.L.S., A.A.A., and F.N. gratefully acknowledge the Max Planck Society for financial support.

## NOMENCLATURE

*a*,*b*,*c*,*d*canonical virtual orbitals

*i*,*j*,*k*,*l*(FB-localized) occupied orbitals or, in particular, valence orbitals

*m*,*n*(FB-localized) core orbitals

*p*,*q*,*r*,*s*occupied or virtual molecular orbitals

*μ*,*ν*,*η*atomic orbitals (AOs)

- $\xe3,b\u0303,c\u0303,d\u0303$
pair natural orbitals (PNOs)

- $\xe3\u2032,b\u0303\u2032,c\u0303\u2032,d\u0303\u2032$
discarded/“complementary” pair natural orbitals (CPNOs)

- $\xe3\u2033,b\u0303\u2033,c\u0303\u2033,d\u0303\u2033$
union of PNOs and CPNOs

- $\mu \u0303,\nu \u0303,\eta \u0303$
orbital/pair domain of orthonormal pseudo-canonical non-redundant projected atomic orbitals (NPAOs)

- $\mu \u0303\u2032,\nu \u0303\u2032,\eta \u0303\u2032$
normalized redundant projected atomic orbitals (PAOs)

### APPENDIX: DERIVATION OF THE PNO ROTATION STATIONARITY CONDITIONS

In the present work, the MO and PNO response is parameterized using full, complex-valued matrices, constrained with an explicit orthonormality condition, in contrast to the exponential parameterization used previously for DLPNO-MP2.^{69} This has the consequence that the final working equations for electric/magnetic perturbations, which are purely real/imaginary, must be obtained as linear combinations of several Lagrangian stationarity conditions. In this section, we show a detailed derivation of the equations related to PNO rotations and briefly discuss the analogous treatment of MO rotations.

The PNO rotation stationarity conditions have the following form:

The separate terms are as follows:

where quantities in the full PNO + CPNO basis are expanded as

The index $\mu \u0303\u2032$ of the intermediate $\tau \xe3\mu \u0303\u2032ij$ spans the PAO domains of all pairs *ik* and *jk* included in the sum, i.e., those for which $\u2009Fkj\u2009$ or $\u2009Fki\u2009$, respectively, is greater than *F*_{Cut}. Note also that the amplitudes **T**^{ij} are implicitly non-zero only in the PNO block. To obtain equations for **v**^{ij}, we take the following linear combination, where the contributions from *C*_{SC} and *C*_{PNOO} drop out:

Only rotations between PNOs and CPNOs result in the non-zero equations shown in Eq. (58). In addition, we use Eq. (55) (at **λ** = 0) to express Eq. (A5) as

and after expanding the semi-canonical residuals in Eq. (A7), we obtain a different expression for Eq. (A1),

where, padding **D**^{ij} with zeros up to the full PNO + CPNO basis, we have substituted

If we define $L\u2212PNOO=L\u2212CPNOO$, we obtain a solution for **x**^{ij},

Note that we obtain a different expression from the Hermitian conjugate,

but both equations are equivalent, provided that Eq. (A12) is fulfilled.

To obtain the perturbed multipliers **v**^{ij,λ}, we first express the perturbed PNO/CPNO stationarity conditions as

in which **x**^{ij,λ} cancels out to obtain Eq. (107). A solution for **x**^{ij,λ} can also be obtained from Eq. (A18),

This is equal to the Hermitian conjugate, obtained from Eq. (A19), provided that Eq. (A20) is fulfilled.

The derivation of the MO rotation stationarity conditions in analogous—the contributions from *C*_{MOO} cancel in the linear combination,

which gives the Z-CPL, Z-CV, and Z-CPSCF equations (59), (67), and (69), respectively, while solutions for **x** can be obtained from either equation separately,

which are equivalent, provided that Eq. (A22) is fulfilled. The second derivatives are also treated analogously—we first express the stationarity conditions as

and then obain the perturbed Z-CPL, Z-CV, and Z-CPSCF equations (109), (112), and (113), respectively, from the linear combination in which **x**^{λ} cancels,