DLPNO-MP2 second derivatives for the computation of polarizabilities and NMR shieldings

Stoychev, Georgi L.; Auer, Alexander A.; Gauss, Jürgen; Neese, Frank

doi:10.1063/5.0047125

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⁴²—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).⁵⁴ The first efficient such implementation, also employing the resolution of the identity (RI) approximation, was later reported by Loibl and Schütz⁵⁵ and also adapted to magnetizabilities.⁵⁶ Another study by Werner and co-workers examines the accuracy of PAO-based local correlation methods for polarizability calculations via finite differences,⁵⁷ 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.⁵⁸ Static and dynamic polarizabilities have also been implemented for MP2 and the approximate coupled cluster method CC2 using both RI and a Laplace-transformation.⁵⁹ 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.⁶⁰ The complete analytic gradient for an orbital-specific virtual (OSV) local MP2 method was also published recently by Zhou et al.⁶¹ 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.⁶⁸

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),⁷⁶ the Laplace-based approaches of Ochsenfeld, Hättig, and their co-workers,^25,59 the proof-of-concept GIAO-LMP2 implementation,⁵⁴ 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

α_{x y} = \frac{d^{2} E}{d F_{x}^{ext} d F_{y}^{ext}} = \sum_{μ ν} \frac{d D_{μ ν}}{d F_{x}^{ext}} \frac{d h_{μ ν}}{d F_{y}^{ext}},

(1)

σ_{y x}^{A} = \frac{d^{2} E}{d B_{x}^{ext} d M_{y}^{A}} = \sum_{μ ν} [D_{μ ν} \frac{d^{2} h_{μ ν}}{d B_{x}^{ext} d M_{y}^{A}} + \frac{d D_{μ ν}}{d B_{x}^{ext}} \frac{d h_{μ ν}}{d M_{y}^{A}}],

(2)

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.⁷⁸ In contrast, the “symmetric” formulation, commonly used for vibrational frequency calculations,^88–90 would also require derivatives with respect to the 3N_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., ${\tilde{μ}}_{i}$ or $ã_{i j}$ ⁠. However, the occupied orbital indices are often obvious from context or are added as upper indices of a matrix, e.g., $F_{\tilde{μ}}^{i}$ or $T_{ã \tilde{b}}^{i j}$ ⁠, and thus, the repetition is avoided to simplify the notation.

Be aware that despite the similar notation, the symbols $K_{ã \tilde{b}}^{i j}$ ⁠, $T_{ã \tilde{b}}^{i j}$ ⁠, $F_{ã \tilde{b}}^{i j}$ ⁠, and $S_{ã \tilde{b}}^{i j}$ have rather different meanings: the first denotes the two-electron repulsion integral $(i ã_{i j} | j {\tilde{b}}_{i j})$ ⁠, the second parameterizes excitations from orbitals i and j into PNOs $ã$ and $\tilde{b}$ ⁠, 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:

\frac{\partial A}{\partial b_{i j}} = C_{i j} \Leftrightarrow \frac{\partial A}{\partial b} = C,

(3)

A_{i j}^{b} \equiv \frac{d A_{i j}}{d b} = C_{i j} \Leftrightarrow A^{b} \equiv \frac{d A}{d b} = C .

(4)

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 $tr (A) = \sum_{i} A_{i i}$ ⁠.

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.¹⁹ 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

|{\tilde{μ}}_{u}^{'}⟩ = (1 - \sum_{i} |i⟩ ⟨i|) |μ⟩ \equiv \sum_{ν} {\tilde{P}}_{ν {\tilde{μ}}_{u}^{'}}^{u} |ν⟩

(5)

and then normalized using the factor $N_{{\tilde{μ}}^{'}} = {⟨{\tilde{μ}}_{u}^{'} | {\tilde{μ}}_{u}^{'}⟩}^{- \frac{1}{2}}$ ⁠,

|{\tilde{μ}}^{'}⟩ = N_{{\tilde{μ}}^{'}} |{\tilde{μ}}_{u}^{'}⟩ \equiv \sum_{ν} {\tilde{P}}_{ν {\tilde{μ}}^{'}} |ν⟩ .

(6)

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

{DOI}_{i {\tilde{μ}}^{'}} = \sqrt{⟨{|i|}^{2} | {|{\tilde{μ}}^{'}|}^{2}⟩} > T_{CutDOPre},

(7)

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⁻⁸ by default) and normalizing the rest. These are then transformed into pseudo-canonical NPAOs ${\tilde{μ}}$ 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⁻⁵ and the semi-canonical pair energy estimate [Eq. (8)] using the collinear dipole approximation to the exchange integrals⁹⁷ [Eq. (8) computed with $M_{\tilde{μ} \tilde{ν}}^{i j}$ from Eq. (10) instead of Eq. (9)] is greater than T_CutPre = 10⁻⁶. The energy contribution from the screened-out pairs, ΔE_Pre, is estimated using the more accurate dipole approximation to the exchange integrals¹¹ [Eq. (9)],

ε_{i j}^{dip} = - 4 \sum_{\tilde{μ} \tilde{ν}} \frac{M_{\tilde{μ} \tilde{ν}}^{i j} M_{\tilde{μ} \tilde{ν}}^{i j *}}{F_{\tilde{μ}}^{i} + F_{\tilde{ν}}^{j} - F_{i i} - F_{j j}},

(8)

M_{\tilde{μ} \tilde{ν}}^{i j} = \frac{r_{i \tilde{μ}} r_{j \tilde{ν}}}{{|R_{i j}|}^{3}} - 3 \frac{(r_{i \tilde{μ}} R_{i j}) (r_{j \tilde{ν}} R_{i j})}{{|R_{i j}|}^{5}},

(9)

M_{\tilde{μ} \tilde{ν}}^{i j, col} = 4 \frac{r_{i \tilde{μ}} r_{j \tilde{ν}}}{{|R_{i j}|}^{3}},

(10)

R_{i j} = r_{i i} - r_{j j},

(11)

r_{i \tilde{μ}} = ⟨i| r |{\tilde{μ}}_{i}⟩,

(12)

where $F_{\tilde{μ}}^{i}$ and $F_{\tilde{ν}}^{j}$ 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 ${\overset{⌣}{D}}^{i j}$ in the basis of NPAOs for the given orbital pair,

{\overset{⌣}{D}}_{\tilde{μ} \tilde{ν}}^{i j} = \sum_{\tilde{η}} [{\tilde{\overset{⌣}{T}}}_{\tilde{μ} \tilde{η}}^{i j *} {\overset{⌣}{T}}_{\tilde{ν} \tilde{η}}^{i j} + {\tilde{\overset{⌣}{T}}}_{\tilde{η} \tilde{μ}}^{i j *} {\overset{⌣}{T}}_{\tilde{η} \tilde{ν}}^{i j}],

(13)

{\overset{⌣}{T}}_{\tilde{μ} \tilde{ν}}^{i j} = - \frac{({\tilde{μ}}_{i j} i | {\tilde{ν}}_{i j} j)}{F_{\tilde{μ}}^{i j} + F_{\tilde{ν}}^{i j} - F_{i i} - F_{j j}},

(14)

{\tilde{\overset{⌣}{T}}}_{\tilde{μ} \tilde{ν}}^{i j} = \frac{1}{1 + δ_{i j}} (4 {\overset{⌣}{T}}_{\tilde{μ} \tilde{ν}}^{i j} - 2 {\overset{⌣}{T}}_{\tilde{ν} \tilde{μ}}^{i j}),

(15)

where $F_{\tilde{μ}}^{i j}$ and $F_{\tilde{ν}}^{i j}$ are NPAO orbital energies in the ij pair domain. PNOs with an occupation number greater than T_CutPNO = 10⁻⁸ are kept, while the rest form the “complementary” PNO (CPNO) basis referred to below. A tighter threshold, T_CutPNO × T_ScalePNOCore = 10⁻¹⁰, is used for pairs that include core orbitals.⁹⁸ 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 $\partial E_{2}^{D L P N O} / \partial {\tilde{T}}_{ã \tilde{b}}^{i j *} = 0$ (see Sec. II C for definitions) are solved iteratively in the PNO basis, and the final DLPNO-MP2 energy is calculated,

E_{DLPNO - MP 2} = E_{HF} + E_{2}^{D L P N O} + Δ E_{P N O} + Δ E_{P r e} .

(16)

Sections II C–II E list the steps necessary to compute analytic derivatives based on the DLPNO-MP2 energy.

C. DLPNO-MP2 Lagrangian

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

\begin{align} L & = E_{HF} + E_{2}^{DLPNO} + Δ E_{PNO} + Δ E_{Pre} + C_{Bri} + C_{Loc} + C_{CV} \\ + C_{SC} + C_{PNO} + C_{MOO} + C_{PNOO} + C_{NPAOO} . \end{align}

(17)

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,

E_{HF} = 2 \sum_{i} h_{i i} + \sum_{i j} [2 (i i | j j) - (i j | j i)],

(18)

with two-electron integrals in the Mulliken (1^*1|2^*2) notation. $E_{2}^{D L P N O}$ is the PNO-basis Hylleraas functional in the spin-adapted form¹⁰⁰

E_{2}^{D L P N O} = \sum_{i \geq j} \sum_{ã \tilde{b}} [K_{ã \tilde{b}}^{i j} {\tilde{T}}_{ã \tilde{b}}^{i j} + K_{ã \tilde{b}}^{i j *} {\tilde{T}}_{ã \tilde{b}}^{i j *} + D_{ã \tilde{b}}^{i j} F_{ã \tilde{b}}^{i j}] - \sum_{i j} D_{i j} F_{i j},

(19)

D_{ã \tilde{b}}^{i j} = \sum_{\tilde{c}} [{\tilde{T}}_{ã \tilde{c}}^{i j *} T_{\tilde{b} \tilde{c}}^{i j} + {\tilde{T}}_{\tilde{c} ã}^{i j *} T_{\tilde{c} \tilde{b}}^{i j}],

(20)

D_{i j} = \sum_{k} (1 + δ_{j k}) \sum_{ã \tilde{b} \tilde{c} \tilde{d}} S_{ã \tilde{b}}^{i k, k j *} {\tilde{T}}_{\tilde{b} \tilde{c}}^{k j *} S_{\tilde{c} \tilde{d}}^{k j, i k} T_{\tilde{d} ã}^{i k},

(21)

K_{ã \tilde{b}}^{i j} = (i ã_{i j} | j {\tilde{b}}_{i j}),

(22)

S_{ã \tilde{b}}^{i j, k l} = ⟨ã_{i j} | {\tilde{b}}_{k l}⟩,

(23)

{\tilde{T}}_{ã \tilde{b}}^{i j} = \frac{1}{1 + δ_{i j}} (4 T_{ã \tilde{b}}^{i j} - 2 T_{\tilde{b} ã}^{i j}) .

(24)

Δ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,

\begin{align} Δ E_{P N O} & = \frac{1}{2} \sum_{i \geq j} [\sum_{ã^{^{″}} {\tilde{b}}^{^{″}}} (K_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j} {\tilde{\overset{⌣}{T}}}_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j} + K_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j *} {\tilde{\overset{⌣}{T}}}_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j *})) \\ - (\sum_{ã \tilde{b}} (K_{ã \tilde{b}}^{i j} {\tilde{\overset{⌣}{T}}}_{ã \tilde{b}}^{i j} + K_{ã \tilde{b}}^{i j *} {\tilde{\overset{⌣}{T}}}_{ã \tilde{b}}^{i j *})] . \end{align}

(25)

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

\begin{align} Δ E_{P r e} & = 4 \sum_{i > j} 1_{S} (i j) \sum_{\tilde{μ} \tilde{ν}} [M_{\tilde{μ} \tilde{ν}}^{i j} t_{\tilde{μ} \tilde{ν}}^{i j} + M_{\tilde{μ} \tilde{ν}}^{i j *} t_{\tilde{μ} \tilde{ν}}^{i j *} - (F_{i i} + F_{j j})) \\ \times (t_{\tilde{μ} \tilde{ν}}^{i j *} t_{\tilde{μ} \tilde{ν}}^{i j} + \sum_{\tilde{η}} (F_{\tilde{μ} \tilde{ν}}^{i} t_{\tilde{μ} \tilde{η}}^{i j *} t_{\tilde{ν} \tilde{η}}^{i j} + F_{\tilde{μ} \tilde{ν}}^{j} t_{\tilde{η} \tilde{μ}}^{i j *} t_{\tilde{η} \tilde{ν}}^{i j})], \end{align}

(26)

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

C_{Bri} = \frac{1}{2} \sum_{a i} (z_{a i} F_{a i} + z_{a i}^{*} F_{a i}^{*}),

(27)

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., $z_{i a} \equiv z_{a i}^{*}$ ⁠. 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¹⁰¹

c (λ) = c^{(0)} U (λ), U (0) = I,

(28)

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,

C_{MOO} = \sum_{p q} x_{p q} (S_{p q} - δ_{p q}),

(29)

which is fulfilled by the CMO coefficients and the unperturbed LMO coefficients c⁽⁰⁾. 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,

C_{Loc} = \frac{1}{2} \sum_{i < j}^{v a l e n c e} (l_{i j} s_{i j} + l_{i j}^{*} s_{i j}^{*}) + \frac{1}{2} \sum_{m < n}^{c o r e} (l_{i j} s_{i j} + l_{i j}^{*} s_{i j}^{*}),

(30)

s_{i j} = (r_{i i} - r_{j j}) r_{i j} = - s_{j i}^{*} .

(31)

The final MO constraint is the core–valence separation condition with an associated Lagrange parameter ${\bar{z}}_{i m}$ ⁠,

C_{CV} = \frac{1}{2} \sum_{i}^{v a l e n c e} \sum_{m}^{c o r e} ({\bar{z}}_{i m} F_{i m} + {\bar{z}}_{i m}^{*} F_{i m}^{*}) .

(32)

Recall that the core and valence orbitals are localized separately because different PNO truncation thresholds are used for pairs that include core orbitals.⁹⁸ 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 $ã, \tilde{b}, \tilde{c}, \tilde{d}$ ⁠, the CPNOs by $ã^{'}, {\tilde{b}}^{'}, {\tilde{c}}^{'}, {\tilde{d}}^{'}$ ⁠, and the union of the two spaces by $ã^{^{″}}, {\tilde{b}}^{^{″}}, {\tilde{c}}^{^{″}}, {\tilde{d}}^{^{″}}$ ⁠. The transformation matrices from redundant PAOs to each of these sets (for orbital pair ij) are denoted by $d_{{\tilde{μ}}^{'} ã}^{i j}$ ⁠, $d_{{\tilde{μ}}^{'} ã^{'}}^{' i j}$ ⁠, and $d_{{\tilde{μ}}^{'} ã^{^{″}}}^{'' i j}$ ⁠, 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

d^{'' i j} (λ) = d^{'' i j (0)} θ^{i j} (λ), θ^{i j} (0) = I .

(33)

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

C_{PNOO} = \sum_{i \geq j} \sum_{ã^{^{″}} {\tilde{b}}^{^{″}}} x_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j} (S_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j} - δ_{ã^{^{″}} {\tilde{b}}^{^{″}}}) with S^{i j} \equiv S^{i j, i j} .

(34)

Note that S^ij is expanded as

S^{i j} = θ^{i j †} S^{i j} θ^{i j} = θ^{i j †} d^{'' i j (0) †} \tilde{S} d^{'' i j (0)} θ^{i j},

(35)

{\tilde{S}}_{{\tilde{μ}}^{'} {\tilde{ν}}^{'}} = ⟨{\tilde{μ}}^{'} | {\tilde{ν}}^{'}⟩ .

(36)

A further PNO constraint ensures that the semi-canonical pair densities ${\overset{⌣}{D}}^{i j}$ remain block-diagonal,

C_{PNO} = \frac{1}{2} \sum_{i \geq j} \sum_{ã {\tilde{b}}^{'}} (v_{ã {\tilde{b}}^{'}}^{i j} {\overset{⌣}{D}}_{ã {\tilde{b}}^{'}}^{i j} + v_{ã {\tilde{b}}^{'}}^{i j *} {\overset{⌣}{D}}_{ã {\tilde{b}}^{'}}^{i j *}),

(37)

where $v_{ã {\tilde{b}}^{'}}^{i j}$ are unknown Lagrange multipliers. An orthonormality condition also applies to the orbital domain NPAOs used for the prescreening correction, parameterized as

π^{i} (λ) = π^{i (0)} θ^{i} (λ), θ^{i} (0) = I,

(38)

C_{NPAOO} = \sum_{i} \sum_{\tilde{μ} \tilde{ν}} x_{\tilde{μ} \tilde{ν}}^{i} (S_{\tilde{μ} \tilde{ν}}^{i} - δ_{\tilde{μ} \tilde{ν}}),

(39)

where $π_{{\tilde{μ}}^{'} \tilde{ν}}^{i}$ are the transformation matrices from redundant to non-redundant orbital domain PAOs, Sⁱ is the overlap matrix in the latter basis, and $x_{\tilde{μ} \tilde{ν}}^{i}$ are the necessary Lagrange multipliers. Finally, the semi-canonical amplitudes are parameterized as

{\overset{⌣}{T}}^{i j} = θ^{i j †} {\overset{⌢}{T}}^{i j} θ^{i j *},

(40)

and the semi-canonical residual conditions C_SC are included in the Lagrangian,

C_{SC} = \frac{1}{2} \sum_{i \geq j} \sum_{ã^{^{″}} {\tilde{b}}^{^{″}}} (w_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j} {\overset{⌣}{R}}_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j} + w_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j *} {\overset{⌣}{R}}_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j *}),

(41)

where $w_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j}$ are Lagrange multipliers and the semi-canonical residuals ${\overset{⌣}{R}}_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j}$ are defined as

{\overset{⌣}{R}}_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j} = K_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j *} + \sum_{{\tilde{c}}^{^{″}}} ({\overset{⌣}{T}}_{ã^{^{″}} {\tilde{c}}^{^{″}}}^{i j} F_{{\tilde{b}}^{^{″}} {\tilde{c}}^{^{″}}}^{i j} + F_{ã^{^{″}} {\tilde{c}}^{^{″}}}^{i j} {\overset{⌣}{T}}_{{\tilde{c}}^{^{″}} {\tilde{b}}^{^{″}}}^{i j}) - {\overset{⌣}{T}}_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j} (F_{i i} + F_{j j}) .

(42)

The optimization parameters in the present DLPNO-MP2 Lagrangian are thus U_pq, $T_{ã \tilde{b}}^{i j}$ ⁠, ${\overset{⌢}{T}}_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j}$ ⁠, $θ_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j}$ ⁠, $θ_{\tilde{μ} \tilde{ν}}^{i}$ ⁠, $t_{\tilde{μ} \tilde{ν}}^{i j}$ ⁠, z_ai, ${\bar{z}}_{i m}$ ⁠, l_ij, x_pq, $v_{ã {\tilde{b}}^{'}}^{i j}$ ⁠, $x_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j}$ ⁠, $x_{\tilde{μ} \tilde{ν}}^{i}$ ⁠, $w_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j}$ ⁠, 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,

\frac{d E_{DLPNO - MP 2}}{d κ} = \frac{\partial L}{\partial κ} = \sum_{μ ν} (P_{μ ν}^{SCF} + D_{μ ν}) h_{μ ν}^{κ},

(43)

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,

P_{μ ν}^{SCF} = 2 \sum_{i} c_{μ i}^{*} c_{ν i},

(44)

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

D = D^{'} + \frac{1}{2} c^{*} (z + \bar{z}) c^{T},

(45)

where, as stated above, only the occupied–virtual and virtual–occupied blocks of z are non-zero with $z_{i a} \equiv z_{a i}^{*}$ ⁠. Analogously, $\bar{z}$ only has non-zero core–valence blocks with ${\bar{z}}_{m i} \equiv {\bar{z}}_{i m}^{*}$ ⁠. D′ is the “PNO-relaxed” difference density matrix,

D^{'} = {\tilde{P}}^{*} D^{' v} {\tilde{P}}^{T} - c^{*} D^{' o} c^{T},

(46)

D_{i j}^{' o} = D_{i j} + \frac{1}{2} \sum_{k \geq l} δ_{i j} (δ_{i k} + δ_{i l}) tr (w^{k l} {\overset{⌣}{T}}^{l k} + w^{k l *} {\overset{⌣}{T}}^{l k *}) + D_{i j}^{o (P r e)},

(47)

D_{{\tilde{μ}}^{'} {\tilde{ν}}^{'}}^{' v} = \sum_{i \geq j} D_{{\tilde{μ}}^{'} {\tilde{ν}}^{'}}^{' i j} + D_{{\tilde{μ}}^{'} {\tilde{ν}}^{'}}^{v (P r e)},

(48)

\begin{align} D^{' i j} & = d^{i j *} D^{i j} d^{i j, T} \\ + \frac{1}{2} d^{'' i j *} (w^{i j} {\overset{⌣}{T}}^{j i} + w^{i j, T} {\overset{⌣}{T}}^{i j} + {\overset{⌣}{T}}^{j i *} w^{i j *} + {\overset{⌣}{T}}^{i j *} w^{i j †}) d^{'' i j, T} . \end{align}

(49)

The prescreening contributions are defined as

D_{i j}^{o (P r e)} = 4 δ_{i j} \sum_{k} 1_{S} (i k) tr (t^{i k *} t^{k i}),

(50)

D^{v (P r e)} = 4 \sum_{i j} 1_{S} (i j) π^{i *} t^{i j *} t^{j i} π^{i, T} .

(51)

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

\begin{align} 0 & \equiv \frac{\partial L}{\partial T_{ã \tilde{b}}^{i j}} = \frac{\partial L}{\partial w_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j}} = \frac{\partial L}{\partial v_{ã {\tilde{b}}^{'}}^{i j}} = \frac{\partial L}{\partial z_{a i}} = \frac{\partial L}{\partial {\bar{z}}_{i m}} = \frac{\partial L}{\partial l_{i j}} \\ = \frac{\partial L}{\partial x_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j}} = \frac{\partial L}{\partial x_{p q}} = \frac{\partial L}{\partial {\overset{⌢}{T}}_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j}} = \frac{\partial L}{\partial θ_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j}} = \frac{\partial L}{\partial U_{p q}} . \end{align}

(52)

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,

\begin{align} 0 & \equiv \frac{\partial L}{\partial {\tilde{T}}_{ã \tilde{b}}^{i j *}} = \frac{\partial E_{2}^{DLPNO}}{\partial {\tilde{T}}_{ã \tilde{b}}^{i j *}} \\ = K_{ã \tilde{b}}^{i j *} + \sum_{\tilde{c}} (T_{ã \tilde{c}}^{i j} F_{\tilde{b} \tilde{c}}^{i j} + F_{ã \tilde{c}}^{i j} T_{\tilde{c} \tilde{b}}^{i j}) \\ - \sum_{k \tilde{c} \tilde{d}} (F_{k j} S_{ã \tilde{c}}^{i j, i k} T_{\tilde{c} \tilde{d}}^{i k} S_{\tilde{d} \tilde{b}}^{i k, i j *} + F_{k i} S_{ã \tilde{c}}^{i j, k j} T_{\tilde{c} \tilde{d}}^{k j} S_{\tilde{d} \tilde{b}}^{k j, i j *}) . \end{align}

(53)

Note that the equations obtained from the derivative with respect to the contravariant amplitudes, $\partial L / \partial {\tilde{T}}^{i j}$ ⁠, are equivalent to a linear combination of $\partial L / \partial T^{i j}$ and $\partial L / \partial T^{j i}$ 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⁻⁵.

The semi-canonical amplitude stationarity conditions take the form

\begin{align} 0 & \equiv \frac{\partial L}{\partial {\overset{⌢}{T}}^{i j}} = θ^{i j *} [\frac{1}{1 + δ_{i j}} (2 K^{i j} - K^{j i} - 2 N^{i j} K^{i j} N^{i j} + N^{i j} K^{j i} N^{i j})) \\ + (\frac{1}{2} (v^{i j, T} {\tilde{\overset{⌣}{T}}}^{i j *} + {\tilde{\overset{⌣}{T}}}^{i j *} v^{i j} + F^{i j *} w^{i j} + w^{i j} F^{i j} - w^{i j} (F_{j j} + F_{i i}))] θ^{i j †}, \end{align}

(54)

0 \equiv \frac{\partial L}{\partial {\overset{⌢}{T}}^{i j *}} = {(\frac{\partial L}{\partial {\overset{⌢}{T}}^{i j}})}^{*} .

(55)

The matrix $N_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j} = δ_{\tilde{c} ã^{^{″}}} δ_{\tilde{c} {\tilde{b}}^{^{″}}}$ is introduced to “select” the truncated PNO basis. Equation (54) can be transformed to the pseudo-canonical NPAO basis using the unitary matrix $S_{ã^{^{″}} \tilde{μ}}^{i j (0)}$ and solved for w^ij at λ = 0 (i.e., θ^ij = I),

\begin{align} (F_{\tilde{μ}}^{i j} + F_{\tilde{ν}}^{i j} - F_{j j} - F_{i i}) w_{\tilde{μ} \tilde{ν}}^{i j} \\ = - \sum_{ã^{^{″}} {\tilde{b}}^{^{″}}} S_{ã^{^{″}} \tilde{μ}}^{i j (0)} [\frac{1 - δ_{ã ã^{^{″}}} δ_{\tilde{b} {\tilde{b}}^{^{″}}}}{1 + δ_{i j}} (4 K^{i j} - 2 K^{j i})) \\ + {((v^{i j, T} {\tilde{\overset{⌣}{T}}}^{i j *} + {\tilde{\overset{⌣}{T}}}^{i j *} v^{i j})]}_{ã^{^{″}} {\tilde{b}}^{^{″}}} S_{{\tilde{b}}^{^{″}} \tilde{ν}}^{i j (0)}, \end{align}

(56)

with

w_{\tilde{μ} \tilde{ν}}^{i j} = \sum_{ã^{^{″}} {\tilde{b}}^{^{″}}} S_{ã^{^{″}} \tilde{μ}}^{i j (0)} w_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j} S_{{\tilde{b}}^{^{″}} \tilde{ν}}^{i j (0)} .

(57)

Equations for v^ij are obtained from the Lagrangian stationarity conditions with respect to $θ_{a b^{'}}^{i j}$ ⁠,

\begin{align} \sum_{\tilde{c}} v_{\tilde{c} {\tilde{b}}^{'}}^{i j} {\overset{⌣}{D}}_{\tilde{c} ã}^{i j} - \sum_{{\tilde{c}}^{'}} {\overset{⌣}{D}}_{{\tilde{b}}^{'} {\tilde{c}}^{'}}^{i j} v_{ã {\tilde{c}}^{'}}^{i j} \\ = 2 {(\frac{\partial E_{2}^{D L P N O}}{\partial θ_{{\tilde{b}}^{'} ã}^{i j}}|}_{λ = 0} - \sum_{\tilde{c}} [K_{{\tilde{b}}^{'} \tilde{c}}^{j i} {\tilde{\overset{⌢}{T}}}_{\tilde{c} ã}^{i j} + K_{{\tilde{b}}^{'} \tilde{c}}^{i j} {\tilde{\overset{⌢}{T}}}_{ã \tilde{c}}^{i j} + K_{ã \tilde{c}}^{i j *} {\tilde{\overset{⌢}{T}}}_{{\tilde{b}}^{'} \tilde{c}}^{i j *} + K_{\tilde{c} ã}^{i j *} {\tilde{\overset{⌢}{T}}}_{\tilde{c} {\tilde{b}}^{'}}^{i j *}] . \end{align}

(58)

These can be solved in the PNO + CPNO basis, which diagonalizes ${\overset{⌣}{D}}^{i j}$ ⁠. A detailed derivation, as well as expressions for x^ij, is given in the Appendix. The remaining Lagrange multipliers z, $\bar{z}$ ⁠, 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)¹⁰² and Z-vector coupled-perturbed localization (Z-CPL)⁴⁵ equations for z and l, respectively, while rotations between core and valence orbitals result in the Z-vector core–valence (Z-CV) equations for $\bar{z}$ ⁠. For rotations within a block of localized orbitals (either core or valence), we obtain the Z-CPL equations

R_{i j}^{loc} [l] - {(R_{j i}^{loc} [l])}^{*} = - L_{i j},

(59)

\begin{align} L_{i j} & = 2 \sum_{k} D_{i k}^{' o} F_{j k} - 2 \sum_{k} F_{k i} D_{k j}^{' o} \\ + \sum_{K {\tilde{μ}}^{'}} Γ_{j {\tilde{μ}}^{'}}^{K *} (K | i {\tilde{μ}}^{'}) - \sum_{K {\tilde{μ}}^{'}} Γ_{i {\tilde{μ}}^{'}}^{K} (j {\tilde{μ}}^{'} | K) + L_{i j}^{P r e}, \end{align}

(60)

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):

R_{p i}^{l o c} [l] = \sum_{k} [l_{k i} r_{k p} (r_{k k} - r_{i i}) - r_{i p} (r_{k i} l_{k i} + r_{k i}^{*} l_{k i}^{*})]

with

l^{†} = - l .

(61)

$Γ_{i {\tilde{μ}}^{'}}^{K}$ is the two-body density matrix

Γ_{i {\tilde{μ}}^{'}}^{K} = \sum_{j} (1 + δ_{i j}) G_{K {\tilde{μ}}^{'}}^{i (i j)},

(62)

G_{K {\tilde{μ}}^{'}}^{i (i j)} = \sum_{L ã^{^{″}} {\tilde{b}}^{^{″}}} d_{{\tilde{μ}}^{'} ã^{^{″}}}^{'' i j} g_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j} (j {\tilde{b}}^{^{″}} | L) {(V^{- 1})}_{L K}, with V_{K L} = (K | L),

(63)

g_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j} = δ_{ã^{^{″}} ã} δ_{{\tilde{b}}^{^{″}} \tilde{b}} (2 {\tilde{T}}_{ã \tilde{b}}^{i j} - {\tilde{\overset{⌣}{T}}}_{ã \tilde{b}}^{i j}) + {\tilde{\overset{⌣}{T}}}_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j} + w_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j *} .

(64)

$L_{i j}^{P r e}$ contains contributions from the approximate exchange integrals in the pre-screening correction,

\begin{align} L_{i j}^{P r e} & = \sum_{x} [\sum_{{\tilde{μ}}^{'}} (r_{i {\tilde{μ}}^{'}}^{x *} y_{j {\tilde{μ}}^{'}}^{α, x *} - r_{j {\tilde{μ}}^{'}}^{x} y_{i {\tilde{μ}}^{'}}^{α, x})) \\ + (\frac{1}{2} r_{i j}^{x *} (y_{j}^{β, x} + y_{j}^{β, x *} - y_{i}^{β, x} - y_{i}^{β, x *})], \end{align}

(65)

where

y_{i {\tilde{μ}}^{'}}^{α, x} = 8 \sum_{j} 1_{S} (i j) \sum_{\tilde{ν} \tilde{η}} π_{{\tilde{μ}}^{'} \tilde{ν}}^{i} t_{\tilde{ν} \tilde{η}}^{i j} [\frac{r_{j \tilde{η}}^{x}}{{|R_{i j}|}^{3}} - 3 \frac{R_{i j}^{x} (r_{j \tilde{η}} R_{i j})}{{|R_{i j}|}^{5}}],

(66a)

\begin{matrix} y_{i}^{β, x} & = & \sum_{j} 1_{S} (i j) \sum_{\tilde{μ} \tilde{ν}} t_{\tilde{μ} \tilde{ν}}^{i j} [(48 \frac{r_{i \tilde{μ}} r_{j \tilde{ν}}}{{|R_{i j}|}^{5}} - 240 \frac{(r_{i \tilde{μ}} R_{i j}) (r_{j \tilde{ν}} R_{i j})}{{|R_{i j}|}^{7}})) \\ \times (R_{j i}^{x} - 48 \frac{r_{i \tilde{μ}}^{x} (r_{j \tilde{ν}} R_{i j}) + r_{j \tilde{ν}}^{x} (r_{i \tilde{μ}} R_{i j})}{{|R_{i j}|}^{5}}] . \end{matrix}

(66b)

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,

\sum_{k} {\bar{z}}_{m k} F_{j k} - \sum_{n} F_{n m} {\bar{z}}_{n j} = L_{m j}^{CV},

(67)

\begin{align} L_{m j}^{CV} & = \sum_{K {\tilde{μ}}^{'}} Γ_{j {\tilde{μ}}^{'}}^{K *} (K | m {\tilde{μ}}^{'}) - \sum_{K {\tilde{μ}}^{'}} Γ_{m {\tilde{μ}}^{'}}^{K} (j {\tilde{μ}}^{'} | K) \\ + R_{m j}^{loc} [l] - {(R_{j m}^{loc} [l])}^{*} + L_{m j}^{P r e} . \end{align}

(68)

These are solved in the CMO basis. Note that in frozen-core calculations, the terms involving $Γ_{m {\tilde{μ}}^{'}}^{K}$ ⁠, $y_{m {\tilde{μ}}^{'}}^{α, x}$ ⁠, and $y_{m}^{β, x}$ drop out.

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

\sum_{b} z_{i b} F_{a b} - \sum_{j} F_{j i} z_{j a} + 2 g_{a i} [z] = L_{i a},

(69)

with the two-electron operator defined as

g_{a i} [z] \equiv \sum_{p q} z_{p q} (a i ‖ p q) = \sum_{j b} [z_{j b} (a i ‖ j b) + z_{j b}^{*} (a i ‖ b j)],

(70)

(r s ‖ p q) = (r s | p q) - \frac{1}{2} (r q | p s)

(71)

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

L_{i a} = Y_{i a}^{α} + \sum_{{\tilde{μ}}^{'}} [S_{i μ}^{*} N_{{\tilde{μ}}^{'}} Y_{a {\tilde{μ}}^{'}}^{β} + S_{a μ} N_{{\tilde{μ}}^{'}} Y_{i {\tilde{μ}}^{'}}^{β *} - S_{i μ}^{*} N_{{\tilde{μ}}^{'}} Y_{{\tilde{μ}}^{'}}^{γ} S_{a {\tilde{μ}}^{'}}] .

(72)

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

\begin{align} Y_{i a}^{α} & = - \sum_{K {\tilde{μ}}^{'}} Γ_{i {\tilde{μ}}^{'}}^{K} (a {\tilde{μ}}^{'} | K) - 4 g_{a i} [D^{'} + \frac{1}{2} c^{*} \bar{z} c^{T}] - {(R_{a i}^{loc} [l])}^{*} \\ - \sum_{x {\tilde{μ}}^{'}} r_{a {\tilde{μ}}^{'}}^{x} y_{i {\tilde{μ}}^{'}}^{α, x} - \frac{1}{2} \sum_{x} r_{a i}^{x} (y_{i}^{β, x} + y_{i}^{β, x *}), \end{align}

(73)

\begin{align} Y_{p {\tilde{μ}}^{'}}^{β} = \sum_{K i} Γ_{i {\tilde{μ}}^{'}}^{K *} (K | i p) + 2 \sum_{{\tilde{ν}}^{'}} D_{{\tilde{μ}}^{'} {\tilde{ν}}^{'}}^{' v} F_{p {\tilde{ν}}^{'}} + \sum_{{\tilde{ν}}^{'}} Y_{{\tilde{μ}}^{'} {\tilde{ν}}^{'}}^{δ} S_{p {\tilde{ν}}^{'}} + \sum_{x i} r_{p i}^{x} y_{i {\tilde{μ}}^{'}}^{α, x *}, \end{align}

(74)

\begin{align} Y_{{\tilde{μ}}^{'}}^{γ} & = \frac{1}{2} \sum_{K i} Γ_{i {\tilde{μ}}^{'}}^{K} (i {\tilde{μ}}^{'} | K) + \sum_{{\tilde{ν}}^{'}} D_{{\tilde{μ}}^{'} {\tilde{ν}}^{'}}^{' v} F_{{\tilde{μ}}^{'} {\tilde{ν}}^{'}} + \frac{1}{2} \sum_{{\tilde{ν}}^{'}} Y_{{\tilde{μ}}^{'} {\tilde{ν}}^{'}}^{δ} S_{{\tilde{μ}}^{'} {\tilde{ν}}^{'}} \\ + \frac{1}{2} \sum_{x i} y_{i {\tilde{μ}}^{'}}^{α, x} r_{i {\tilde{μ}}^{'}}^{x} + complex conjugate, \end{align}

(75)

Y_{{\tilde{μ}}^{'} {\tilde{ν}}^{'}}^{δ} = \sum_{i \geq j} (ξ_{{\tilde{μ}}^{'} {\tilde{ν}}^{'}}^{i j} + ξ_{{\tilde{ν}}^{'} {\tilde{μ}}^{'}}^{i j *} - 2 \sum_{ã} d_{{\tilde{μ}}^{'} ã}^{i j *} τ_{ã {\tilde{ν}}^{'}}^{i j}) - y_{{\tilde{μ}}^{'} {\tilde{ν}}^{'}}^{δ} - y_{{\tilde{ν}}^{'} {\tilde{μ}}^{'}}^{δ *},

(76)

\begin{align} ξ_{{\tilde{μ}}^{'} {\tilde{ν}}^{'}}^{i j} & = {(d^{'' i j *} x^{i j} d^{'' i j, T})}_{{\tilde{μ}}^{'} {\tilde{ν}}^{'}} \\ = \sum_{{\tilde{b}}^{^{″}} {\tilde{λ}}^{'}} d_{{\tilde{μ}}^{'} {\tilde{b}}^{^{″}}}^{'' i j *} d_{{\tilde{λ}}^{'} {\tilde{b}}^{^{″}}}^{'' i j} [- \frac{1}{2} \sum_{K} (G_{K {\tilde{ν}}^{'}}^{i (i j)} (i {\tilde{λ}}^{'} | K) + G_{K {\tilde{ν}}^{'}}^{j (i j)} (j {\tilde{λ}}^{'} | K))) \\ - (\sum_{{\tilde{κ}}^{'}} F_{{\tilde{κ}}^{'} {\tilde{λ}}^{'}} D_{{\tilde{κ}}^{'} {\tilde{ν}}^{'}}^{' v} + \sum_{{\tilde{κ}}^{'} ã} S_{{\tilde{κ}}^{'} {\tilde{λ}}^{'}} τ_{ã {\tilde{κ}}^{'}}^{i j *} d_{{\tilde{ν}}^{'} ã}^{i j}], \end{align}

(77)

\begin{align} τ_{ã {\tilde{κ}}^{'}}^{i j} & = \sum_{k} [F_{k j} ({\tilde{T}}^{i j *} S^{i j, k i} T^{k i} + {\tilde{T}}^{j i *} S^{i j, k i} T^{i k}) d^{k i, T}) \\ {(+ F_{k i} ({\tilde{T}}^{i j *} S^{i j, j k} T^{j k} + {\tilde{T}}^{j i *} S^{i j, j k} T^{k j}) d^{j k, T}]}_{ã {\tilde{κ}}^{'}}, \end{align}

(78)

\begin{align} y^{δ} & = - \sum_{i} π^{i *} x^{i} π^{i, T} \\ = 4 \sum_{i j} 1_{S} (i j) π^{i *} (M^{i j} + F^{i *} t^{i j *}) t^{j i} π^{i, T} . \end{align}

(79)

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ⁱ, which are obtained analogously by solving $\partial L / \partial θ^{i} = \partial L / \partial θ^{i †} = 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),

\begin{align} \frac{d^{2} E_{DLPNO - MP 2}}{d λ d κ} & = \frac{d}{d λ} \frac{\partial L}{\partial κ} = \sum_{μ ν} (P_{μ ν}^{SCF} + D_{μ ν}) h_{μ ν}^{κ, λ} \\ + \sum_{μ ν} (P_{μ ν}^{SCF, λ} + D_{μ ν}^{λ}) h_{μ ν}^{κ} . \end{align}

(80)

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^λ,

\begin{align} P_{μ ν}^{SCF, λ} & = 2 \sum_{p i} (c_{μ p}^{*} U_{p i}^{λ *} c_{ν i} + c_{μ i}^{*} U_{p i}^{λ} c_{ν p}) \\ = 2 \sum_{a i} (c_{μ a}^{*} U_{a i}^{λ *} c_{ν i} + c_{μ i}^{*} U_{a i}^{λ} c_{ν a}) - 2 \sum_{i j} c_{μ j}^{*} S_{i j}^{(λ)} c_{ν i}, \end{align}

(81)

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

U_{p q}^{λ} + U_{q p}^{λ *} + S_{p q}^{(λ)} = 0 with S_{p q}^{(λ)} = \sum_{μ ν} c_{μ p}^{(0) *} S_{μ ν}^{λ} c_{ν q}^{(0)} .

(82)

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_{μ ν}^{F^{ext}} = 0$ ⁠. $U_{a i}^{λ}$ are obtained from the perturbed Brillouin condition, which in conjunction with Eq. (82) gives the CPSCF equations,¹⁰¹

\begin{align} \sum_{b} F_{a b} U_{b i}^{λ} - \sum_{j} U_{a j}^{λ} F_{j i} + 2 \sum_{b j} [U_{b j}^{λ *} (a i ‖ b j) + U_{b j}^{λ} (a i ‖ j b)] \\ = \sum_{j} S_{a j}^{(λ)} F_{j i} - F_{a i}^{(λ)} + 2 \sum_{k j} S_{j k}^{(λ)} (a i ‖ k j), \end{align}

(83)

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 $U_{a i}^{λ}$ ⁠. 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 $F_{i m}^{λ} = 0$ [in conjunction with Eq. (82)],

\sum_{j} F_{i j} U_{j m}^{λ} - \sum_{n} U_{i n}^{λ} F_{n m} = \sum_{n} S_{i n}^{(λ)} F_{n m} - F_{i m}^{(λ)} - g_{i m} [P^{λ}],

(84)

which can be solved in the CMO basis. The valence–valence block of U^λ is obtained from the perturbed localization condition $s_{i j}^{λ} = 0$ ⁠, which gives the CPL equations,⁴⁵

\sum_{k < l} U_{k l}^{λ} A_{k l, i j}^{l o c, \pm} = - B_{i j}^{λ, \pm},

(85)

\begin{align} B_{i j}^{λ, \pm} & = s_{i j}^{(λ)} - \sum_{k < l} S_{l k}^{(λ)} (\frac{d s_{i j}}{d U_{l k}} \pm \frac{d s_{j i}}{d U_{l k}}) \\ + \sum_{p k} (δ_{p m} + δ_{p a}) U_{p k}^{λ} (\frac{d s_{i j}}{d U_{p k}} \pm \frac{d s_{j i}}{d U_{p k}}), \end{align}

(86)

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

A_{p q, r s}^{l o c, \pm} = \frac{d s_{r s}}{d U_{p q}} + \frac{d s_{s r}}{d U_{q p}} \pm \frac{d s_{s r}}{d U_{p q}} \pm \frac{d s_{r s}}{d U_{q p}},

(87)

\frac{d s_{r s}}{d U_{p q}} = r_{r s} (r_{r p} δ_{q r} - r_{s p} δ_{q s}) + δ_{q s} r_{r p} (r_{r r} - r_{s s}) .

(88)

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

r_{μ ν, y}^{B_{x}^{ext}} = \frac{i}{2} ⟨μ| {(R_{M N} \times r)}_{x} r_{y} |ν⟩,

(89)

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

U_{m n}^{λ} = U_{n m}^{λ *} = - \frac{1}{2} S_{m n}^{(λ)} .

(90)

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

\sum_{k < l} l_{k l} A_{i j, k l}^{l o c, -} = - L_{i j},

(91)

and $A_{i j, k l}^{l o c, \pm} = A_{k l, i j}^{l o c, \pm}$ 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^λ, ${\bar{z}}^{λ}$ ⁠, T^ij,λ, ${\overset{⌢}{T}}^{i j, λ}$ ⁠, w^ij,λ, θ^ij,λ, t^ij,λ, and θ^i,λ. The PAO coefficients $\tilde{P}$ also depend on the perturbation through U,

{\tilde{P}}_{ν {\tilde{μ}}^{'}}^{λ} = {\tilde{P}}_{ν {\tilde{μ}}_{u}^{'}}^{u, λ} N_{{\tilde{μ}}^{'}} + {\tilde{P}}_{ν {\tilde{μ}}_{u}^{'}}^{u} N_{{\tilde{μ}}^{'}}^{λ},

(92)

{\tilde{P}}_{ν {\tilde{μ}}_{u}^{'}}^{u, λ} = - \frac{1}{2} \sum_{η} (P_{ν η}^{λ} S_{η μ} + P_{ν η} S_{η μ}^{λ}),

(93)

N_{{\tilde{μ}}^{'}}^{λ} = - \frac{1}{2} N_{{\tilde{μ}}^{'}}^{3} S_{{\tilde{μ}}_{u}^{'} {\tilde{μ}}_{u}^{'}}^{λ} .

(94)

For magnetic perturbations, $N_{{\tilde{μ}}^{'}}^{λ} = 0$ ⁠. The perturbed (C)PNO orthonormality condition gives a relationship analogous to Eq. (82),

θ_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j, λ} + θ_{{\tilde{b}}^{^{″}} ã^{^{″}}}^{i j, λ *} + S_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j, λ} = 0,

with

S_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j, λ} = \sum_{{\tilde{μ}}^{'} {\tilde{ν}}^{'}} d_{{\tilde{μ}}^{'} ã^{^{″}}}^{'' i j (0) *} {\tilde{S}}_{{\tilde{μ}}^{'} {\tilde{ν}}^{'}}^{λ} d_{{\tilde{ν}}^{'} {\tilde{b}}^{^{″}}}^{'' i j (0)} .

(95)

Note that the perturbed PAO overlap ${\tilde{S}}_{{\tilde{μ}}^{'} {\tilde{ν}}^{'}}^{λ}$ 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,

\begin{align} (F_{i i} + F_{j j}) {\overset{⌢}{T}}^{i j, λ} - {\overset{⌢}{T}}^{i j, λ} ℱ^{i j *} - ℱ^{i j} {\overset{⌢}{T}}^{i j, λ} \\ = K^{i j, λ *} - {\overset{⌢}{T}}^{i j} (F_{i i}^{λ} + F_{j j}^{λ}) + {\overset{⌢}{T}}^{i j} (F^{i j, λ *} - S^{i j, λ *} F^{i j *}) \\ + (F^{i j, λ} - F^{i j} S^{i j, λ}) {\overset{⌢}{T}}^{i j}, \end{align}

(96)

where

K_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j, λ} = \sum_{{\tilde{μ}}^{'} {\tilde{ν}}^{'}} d_{{\tilde{μ}}^{'} ã^{^{″}}}^{'' i j (0)} K_{{\tilde{μ}}^{'} {\tilde{ν}}^{'}}^{i j, λ} d_{{\tilde{ν}}^{'} {\tilde{b}}^{^{″}}}^{'' i j (0)},

(97)

F_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j, λ} = \sum_{{\tilde{μ}}^{'} {\tilde{ν}}^{'}} d_{{\tilde{μ}}^{'} ã^{^{″}}}^{'' i j (0) *} F_{{\tilde{μ}}^{'} {\tilde{ν}}^{'}}^{λ} d_{{\tilde{ν}}^{'} {\tilde{b}}^{^{″}}}^{'' i j (0)} .

(98)

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

\begin{align} \sum_{\tilde{c}} {\overset{⌣}{D}}_{ã \tilde{c}}^{i j} θ_{{\tilde{b}}^{'} \tilde{c}}^{i j, λ} - \sum_{{\tilde{c}}^{'}} θ_{{\tilde{c}}^{'} ã}^{i j, λ} {\overset{⌣}{D}}_{{\tilde{c}}^{'} {\tilde{b}}^{'}}^{i j} \\ = ({\tilde{\overset{⌢}{T}}}^{i j, λ *} {\overset{⌢}{T}}^{j i} + {\tilde{\overset{⌢}{T}}}^{j i, λ *} {\overset{⌢}{T}}^{i j} + {\tilde{\overset{⌢}{T}}}^{i j *} {\overset{⌢}{T}}^{j i, λ} + {\tilde{\overset{⌢}{T}}}^{j i *} {\overset{⌢}{T}}^{i j, λ}) \\ - {({\overset{⌣}{D}}^{i j} S^{i j, λ *} - {\tilde{\overset{⌢}{T}}}^{i j *} S^{i j, λ} {\overset{⌢}{T}}^{j i} - {\tilde{\overset{⌢}{T}}}^{j i *} S^{i j, λ} {\overset{⌢}{T}}^{i j})}_{ã {\tilde{b}}^{'}} . \end{align}

(99)

These equations are solved in the basis that diagonalizes ${\overset{⌣}{D}}^{i j}$ ⁠. 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,

\begin{align} 0 & \equiv \frac{d}{d λ} \frac{\partial L}{\partial {\tilde{T}}_{ã \tilde{b}}^{i j *}} = X_{ã \tilde{b}}^{i j, λ} + \sum_{\tilde{c}} (T_{ã \tilde{c}}^{i j, λ} F_{\tilde{b} \tilde{c}}^{i j} + F_{ã \tilde{c}}^{i j} T_{\tilde{c} \tilde{b}}^{i j, λ}) \\ - \sum_{k \tilde{c} \tilde{d}} (F_{k j} S_{ã \tilde{c}}^{i j, i k} T_{\tilde{c} \tilde{d}}^{i k, λ} S_{\tilde{d} \tilde{b}}^{i k, i j *} + F_{k i} S_{ã \tilde{c}}^{i j, k j} T_{\tilde{c} \tilde{d}}^{k j, λ} S_{\tilde{d} \tilde{b}}^{k j, i j *}), \end{align}

(100)

\begin{matrix} X_{ã \tilde{b}}^{i j, λ} & = & K_{ã \tilde{b}}^{i j, λ *} + \sum_{\tilde{c}} (T_{ã \tilde{c}}^{i j} F_{\tilde{b} \tilde{c}}^{i j, λ} + F_{ã \tilde{c}}^{i j, λ} T_{\tilde{c} \tilde{b}}^{i j}) \\ - \sum_{k \tilde{c} \tilde{d}} [{(F_{k j} S_{ã \tilde{c}}^{i j, i k} S_{\tilde{d} \tilde{b}}^{i k, i j *})}^{λ} T_{\tilde{c} \tilde{d}}^{i k} + {(F_{k i} S_{ã \tilde{c}}^{i j, k j} S_{\tilde{d} \tilde{b}}^{k j, i j *})}^{λ} T_{\tilde{c} \tilde{d}}^{k j}], \end{matrix}

(101)

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,

θ_{\tilde{μ} \tilde{ν}}^{i, λ} = θ_{\tilde{ν} \tilde{μ}}^{i, λ *} = - \frac{1}{2} S_{\tilde{μ} \tilde{ν}}^{i, λ} with S_{\tilde{μ} \tilde{ν}}^{i, λ} = \sum_{{\tilde{μ}}^{'} {\tilde{ν}}^{'}} π_{{\tilde{μ}}^{'} \tilde{μ}}^{i (0) *} {\tilde{S}}_{{\tilde{μ}}^{'} {\tilde{ν}}^{'}}^{λ} π_{{\tilde{ν}}^{'} \tilde{ν}}^{i (0)} .

(102)

Also required are the perturbed prescreening amplitudes, obtained from the equations $\frac{d}{d λ} \frac{\partial L}{\partial t^{i j *}} = 0$ giving

\begin{align} t_{\tilde{μ} \tilde{ν}}^{i j, λ} & = {(F_{i i} + F_{j j} - F_{\tilde{μ}}^{i} - F_{\tilde{ν}}^{i})}^{- 1} \\ \times [M_{\tilde{μ} \tilde{ν}}^{i j, λ *} - (F_{i i}^{λ} + F_{j j}^{λ}) t_{\tilde{μ} \tilde{ν}}^{i j} + \sum_{\tilde{η}} (F_{\tilde{μ} \tilde{η}}^{i, λ} t_{\tilde{η} \tilde{ν}}^{i j} + t_{\tilde{μ} \tilde{η}}^{i j} F_{\tilde{ν} \tilde{η}}^{j, λ})], \end{align}

(103)

\begin{align} M_{\tilde{μ} \tilde{ν}}^{i j, λ_{x}} & = \frac{r_{i \tilde{μ}}^{λ_{x}} r_{j \tilde{ν}} + r_{i \tilde{μ}} r_{j \tilde{ν}}^{λ_{x}}}{{|R_{i j}|}^{3}} - 3 \frac{(r_{i \tilde{μ}}^{λ_{x}} R_{i j}) (r_{j \tilde{ν}} R_{i j}) + (r_{i \tilde{μ}} R_{i j}) (r_{j \tilde{ν}}^{λ_{x}} R_{i j})}{{|R_{i j}|}^{5}} \\ - 3 \frac{(r_{i \tilde{μ}} R_{i j}^{λ_{x}}) (r_{j \tilde{ν}} R_{i j}) + (r_{i \tilde{μ}} R_{i j}) (r_{j \tilde{ν}} R_{i j}^{λ_{x}})}{{|R_{i j}|}^{5}} \\ - 3 (R_{i j} R_{i j}^{λ_{x}}) (\frac{r_{i \tilde{μ}} r_{j \tilde{ν}}}{{|R_{i j}|}^{5}} - 5 \frac{(r_{i \tilde{μ}} R_{i j}) (r_{j \tilde{ν}} R_{i j})}{{|R_{i j}|}^{7}}) . \end{align}

(104)

For magnetic perturbations, $R_{i j}^{λ} = 0$ ⁠. The perturbed Lagrange multipliers w^ij,λ are obtained from the equations

\sum_{{\tilde{c}}^{^{″}}} (w_{ã^{^{″}} {\tilde{c}}^{^{″}}}^{i j, λ} F_{{\tilde{c}}^{^{″}} {\tilde{b}}^{^{″}}}^{i j} + F_{{\tilde{c}}^{^{″}} ã^{^{″}}}^{i j} w_{{\tilde{c}}^{^{″}} {\tilde{b}}^{^{″}}}^{i j, λ}) - (F_{i i} + F_{j j}) w_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j, λ} = - W_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j, λ},

(105)

which are solved in the pseudo-canonical NPAO basis with

\begin{align} W^{i j, λ} & = \frac{d}{d λ} \{θ^{i j *} [F^{i j *} w^{i j (0)} + w^{i j (0)} F^{i j} - (F_{i i} + F_{j j}) w^{i j (0)})) \\ + \frac{2}{1 + δ_{i j}} (2 K^{i j} - K^{j i} - 2 N^{i j} K^{i j} N^{i j} + N^{i j} K^{j i} N^{i j}) \\ (+ ({\tilde{\overset{⌣}{T}}}^{i j *} v^{i j, T} + v^{i j} {\tilde{\overset{⌣}{T}}}^{i j *}] θ^{i j †}\} . \end{align}

(106)

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

\sum_{\tilde{c}} v_{\tilde{c} {\tilde{b}}^{'}}^{i j, λ} {\overset{⌣}{D}}_{\tilde{c} ã}^{i j} - \sum_{{\tilde{c}}^{'}} {\overset{⌣}{D}}_{{\tilde{b}}^{'} {\tilde{c}}^{'}}^{i j} v_{ã {\tilde{c}}^{'}}^{i j, λ} = - V_{{\tilde{b}}^{'} ã}^{i j, λ},

(107)

with

\begin{align} V_{{\tilde{b}}^{'} ã}^{i j, λ} & = \frac{d}{d λ} [v^{i j (0), T} {\overset{⌣}{D}}^{i j} - {\overset{⌣}{D}}^{i j} v^{i j (0), T} + K^{j i} N^{i j} {\tilde{\overset{⌣}{T}}}^{i j}) \\ + K^{i j} N^{i j} {\tilde{\overset{⌣}{T}}}^{j i} + {\tilde{\overset{⌣}{T}}}^{i j *} N^{i j} K^{j i *} + {\tilde{\overset{⌣}{T}}}^{j i *} N^{i j} K^{i j *} \\ - {(2 (K^{j i} {\tilde{T}}^{i j} + K^{i j} {\tilde{T}}^{j i} + F^{i j *} D^{i j} - d^{'' i j, T} {\tilde{S}}^{*} τ^{i j †})]}_{{\tilde{b}}^{'} ã} . \end{align}

(108)

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

R_{i j}^{loc} [l^{λ}] - {(R_{j i}^{loc} [l^{λ}])}^{*} = {(R_{j i}^{loc, λ} [l^{(0)}])}^{*} - R_{i j}^{loc, λ} [l^{(0)}] - L_{i j}^{λ},

(109)

where $L_{i j}^{λ}$ and $R_{i j}^{loc, λ} [l^{(0)}]$ are straightforwardly derived from Eqs. (60) and (61). The derivatives of the prescreening contributions $y_{i {\tilde{μ}}^{'}}^{α, x}$ and $y_{i}^{β, x}$ are somewhat more involved, so we provide them here,

\begin{align} y_{i {\tilde{μ}}^{'}}^{α, x, λ} & = 8 \sum_{j} 1_{S} (i j) \sum_{\tilde{ν} \tilde{η}} \{{(π_{{\tilde{μ}}^{'} \tilde{ν}}^{i} t_{\tilde{ν} \tilde{η}}^{i j})}^{λ} [\frac{r_{j \tilde{η}}^{x}}{{|R_{i j}|}^{3}} - 3 \frac{R_{i j}^{x} (r_{j \tilde{η}} R_{i j})}{{|R_{i j}|}^{5}}]) \\ (+ π_{{\tilde{μ}}^{'} \tilde{ν}}^{i} t_{\tilde{ν} \tilde{η}}^{i j} [\frac{r_{j \tilde{η}}^{x, λ}}{{|R_{i j}|}^{3}} - 3 \frac{R_{i j}^{x, λ} (r_{j \tilde{η}} R_{i j}) + R_{i j}^{x} (r_{j \tilde{η}}^{λ} R_{i j} + r_{j \tilde{η}} R_{i j}^{λ})}{{|R_{i j}|}^{5}}]) \\ (- π_{{\tilde{μ}}^{'} \tilde{ν}}^{i} t_{\tilde{ν} \tilde{η}}^{i j} [3 \frac{r_{j \tilde{η}}^{x}}{{|R_{i j}|}^{5}} - 15 \frac{R_{i j}^{x} (r_{j \tilde{η}} R_{i j})}{{|R_{i j}|}^{7}}] (R_{i j}^{λ} R_{i j})\}, \end{align}

(110)

\begin{align} y_{i}^{β, x, λ} & = 48 \sum_{j} 1_{S} (i j) \sum_{\tilde{μ} \tilde{ν}} \{t_{\tilde{μ} \tilde{ν}}^{i j, λ} [(\frac{r_{i \tilde{μ}} r_{j \tilde{ν}}}{{|R_{i j}|}^{5}} - 5 \frac{(r_{i \tilde{μ}} R_{i j}) (r_{j \tilde{ν}} R_{i j})}{{|R_{i j}|}^{7}}) R_{j i}^{x} - \frac{r_{i \tilde{μ}}^{x} (r_{j \tilde{ν}} R_{i j}) + r_{j \tilde{ν}}^{x} (r_{i \tilde{μ}} R_{i j})}{{|R_{i j}|}^{5}}] + t_{\tilde{μ} \tilde{ν}}^{i j} [(\frac{r_{i \tilde{μ}} r_{j \tilde{ν}}}{{|R_{i j}|}^{5}} - 5 \frac{(r_{i \tilde{μ}} R_{i j}) (r_{j \tilde{ν}} R_{i j})}{{|R_{i j}|}^{7}}) R_{j i}^{x, λ})) \\ - \frac{r_{i \tilde{μ}}^{x, λ} (r_{j \tilde{ν}} R_{i j}) + r_{j \tilde{ν}}^{x, λ} (r_{i \tilde{μ}} R_{i j})}{{|R_{i j}|}^{5}} + (\frac{{(r_{i \tilde{μ}} r_{j \tilde{ν}})}^{λ}}{{|R_{i j}|}^{5}} - 5 \frac{{(r_{i \tilde{μ}} R_{i j})}^{λ} (r_{j \tilde{ν}} R_{i j}) + (r_{i \tilde{μ}} R_{i j}) {(r_{j \tilde{ν}} R_{i j})}^{λ}}{{|R_{i j}|}^{7}}) R_{j i}^{x} \\ (- (\frac{r_{i \tilde{μ}}^{x} {(r_{j \tilde{ν}} R_{i j})}^{λ} + r_{j \tilde{ν}}^{x} {(r_{i \tilde{μ}} R_{i j})}^{λ}}{{|R_{i j}|}^{5}}] - t_{\tilde{μ} \tilde{ν}}^{i j} [(5 \frac{r_{i \tilde{μ}} r_{j \tilde{ν}}}{{|R_{i j}|}^{7}} - 35 \frac{(r_{i \tilde{μ}} R_{i j}) (r_{j \tilde{ν}} R_{i j})}{{|R_{i j}|}^{9}}) R_{j i}^{x} - 5 \frac{r_{i \tilde{μ}}^{x} (r_{j \tilde{ν}} R_{i j}) + r_{j \tilde{ν}}^{x} (r_{i \tilde{μ}} R_{i j})}{{|R_{i j}|}^{7}}] (R_{j i}^{λ} R_{j i})\} . \end{align}

(111)

The perturbed Z-CV equations have the form

\sum_{k} {\bar{z}}_{m k}^{λ} F_{j k} - \sum_{n} F_{n m} {\bar{z}}_{n j}^{λ} = \sum_{n} F_{n m}^{λ} {\bar{z}}_{n j} - \sum_{k} {\bar{z}}_{m k} F_{j k}^{λ} + L_{m j}^{CV, λ},

(112)

where $L_{m j}^{CV, λ}$ is the derivative of Eq. (68). Finally, the perturbed Z-CPSCF equations must be solved for z^λ,

\begin{align} \sum_{b} z_{i b}^{λ} F_{a b} - \sum_{j} F_{j i} z_{j a}^{λ} + 2 g_{a i} [z^{λ}] \\ = L_{i a}^{λ} - \sum_{b} z_{i b} F_{a b}^{λ} + \sum_{j} F_{j i}^{λ} z_{j a} - 2 {(g_{a i} [z^{(0)}])}^{λ} . \end{align}

(113)

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.

TABLE I.

Interdependencies between the main Lagrange multipliers and target quantities with references to relevant equations. Some terms are excluded for brevity, e.g., those related to the prescreening correction. Wherever a perturbed quantity is required as input (e.g., c^λ), so is the unperturbed one (i.e., c).

In		Equations		Out
		SCF, FB-loc.	→	c, E_HF
c	→	LMP2 guess (13)–(15)	→	${\overset{⌣}{T}}^{i j}, d^{'' i j}, Δ E_{PNO}$
d^ij	→	LMP2 iter. (53)	→	$T^{i j}, E_{2}^{DLPNO}$
${\overset{⌣}{T}}^{i j}, d^{'' i j}, T^{i j}$	→	PNO cond. (58)	→	v^ij
${\overset{⌣}{T}}^{i j}, d^{'' i j}, v^{i j}$	→	SC cond. (56)	→	w^ij
${\overset{⌣}{T}}^{i j}, d^{'' i j}, T^{i j}, w^{i j}$	→	Equation (63)	→	Γ^K
$c, {\overset{⌣}{T}}^{i j}, d^{'' i j}, T^{i j}, w^{i j}$	→	Equation (46)	→	D′
D′, Γ^K	→	Z-CPL (59)	→	l
l, Γ^K	→	Z-CV (67)	→	$\bar{z}$
$c, l, Γ^{K}, D^{'}, \bar{z}$	→	Z-CPSCF (69)	→	z
$c, D^{'}, \bar{z}, z$	→	Equation (45)	→	D
c	→	CPSCF (83), CV (84), CPL (85), $C_{MOO}^{λ}$ (82)	→	U^λ, c^λ
$c^{λ}, {\overset{⌣}{T}}^{i j}, d^{'' i j}$	→	PNO, respectively, (95), (96), and (99)	→	${\overset{⌣}{T}}^{i j, λ}, θ^{i j, λ}, d^{'' i j, λ}$
d^ij,λ, T^ij	→	Equation (100)	→	T^ij,λ
${\overset{⌣}{T}}^{i j, λ}, d^{'' i j, λ}, T^{i j, λ}, v^{i j}$	→	Pert. PNO cond. (107)	→	v^ij,λ
${\overset{⌣}{T}}^{i j, λ}, d^{'' i j, λ}, v^{i j, λ}, w^{i j}$	→	Pert. SC cond. (105)	→	w^ij,λ
${\overset{⌣}{T}}^{i j, λ}, d^{'' i j, λ}, T^{i j, λ}, w^{i j, λ}$	→	Deriv. of Eq. (63)	→	Γ^K,λ
$c^{λ}, {\overset{⌣}{T}}^{i j, λ}, d^{'' i j, λ}, T^{i j, λ}, w^{i j, λ}$	→	Deriv. of Eq. (46)	→	D′^λ
D′^λ, Γ^K,λ, l	→	Pert. Z-CPL (109)	→	l^λ
$l^{λ}, Γ^{K, λ}, \bar{z}$	→	Pert. Z-CV (112)	→	${\bar{z}}^{λ}$
$c^{λ}, l^{λ}, Γ^{K, λ}, D^{' λ}, {\bar{z}}^{λ}, z$	→	Pert. Z-CPSCF (113)	→	z^λ
$c^{λ}, D^{' λ}, {\bar{z}}^{λ}, z^{λ}$	→	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)	→	${\overset{⌣}{T}}^{i j}, d^{'' i j}, Δ E_{PNO}$
d^ij	→	LMP2 iter. (53)	→	$T^{i j}, E_{2}^{DLPNO}$
${\overset{⌣}{T}}^{i j}, d^{'' i j}, T^{i j}$	→	PNO cond. (58)	→	v^ij
${\overset{⌣}{T}}^{i j}, d^{'' i j}, v^{i j}$	→	SC cond. (56)	→	w^ij
${\overset{⌣}{T}}^{i j}, d^{'' i j}, T^{i j}, w^{i j}$	→	Equation (63)	→	Γ^K
$c, {\overset{⌣}{T}}^{i j}, d^{'' i j}, T^{i j}, w^{i j}$	→	Equation (46)	→	D′
D′, Γ^K	→	Z-CPL (59)	→	l
l, Γ^K	→	Z-CV (67)	→	$\bar{z}$
$c, l, Γ^{K}, D^{'}, \bar{z}$	→	Z-CPSCF (69)	→	z
$c, D^{'}, \bar{z}, z$	→	Equation (45)	→	D
c	→	CPSCF (83), CV (84), CPL (85), $C_{MOO}^{λ}$ (82)	→	U^λ, c^λ
$c^{λ}, {\overset{⌣}{T}}^{i j}, d^{'' i j}$	→	PNO, respectively, (95), (96), and (99)	→	${\overset{⌣}{T}}^{i j, λ}, θ^{i j, λ}, d^{'' i j, λ}$
d^ij,λ, T^ij	→	Equation (100)	→	T^ij,λ
${\overset{⌣}{T}}^{i j, λ}, d^{'' i j, λ}, T^{i j, λ}, v^{i j}$	→	Pert. PNO cond. (107)	→	v^ij,λ
${\overset{⌣}{T}}^{i j, λ}, d^{'' i j, λ}, v^{i j, λ}, w^{i j}$	→	Pert. SC cond. (105)	→	w^ij,λ
${\overset{⌣}{T}}^{i j, λ}, d^{'' i j, λ}, T^{i j, λ}, w^{i j, λ}$	→	Deriv. of Eq. (63)	→	Γ^K,λ
$c^{λ}, {\overset{⌣}{T}}^{i j, λ}, d^{'' i j, λ}, T^{i j, λ}, w^{i j, λ}$	→	Deriv. of Eq. (46)	→	D′^λ
D′^λ, Γ^K,λ, l	→	Pert. Z-CPL (109)	→	l^λ
$l^{λ}, Γ^{K, λ}, \bar{z}$	→	Pert. Z-CV (112)	→	${\bar{z}}^{λ}$
$c^{λ}, l^{λ}, Γ^{K, λ}, D^{' λ}, {\bar{z}}^{λ}, z$	→	Pert. Z-CPSCF (113)	→	z^λ
$c^{λ}, D^{' λ}, {\bar{z}}^{λ}, z^{λ}$	→	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₃)₃P(CO)₄Os⋯Cr(CO)₅ molecule (reference code KAMDOR)¹⁰⁶ 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,

\begin{align} C_{Loc}^{valence} & = \sum_{i < j}^{v a l e n c e} (l_{i j}^{R} R [s_{i j}] + l_{i j}^{I} I [s_{i j}]) \\ = \frac{1}{2} \sum_{i < j}^{v a l e n c e} [l_{i j}^{R} (s_{i j} + s_{i j}^{*}) + i l_{i j}^{I} (s_{i j} - s_{i j}^{*})] \end{align}

(114)

= \frac{1}{2} \sum_{i < j}^{v a l e n c e} [\underset{l_{i j}}{\underset{︸}{(l_{i j}^{R} + i l_{i j}^{I})}} s_{i j} + \underset{l_{i j}^{*}}{\underset{︸}{(l_{i j}^{R} - i l_{i j}^{I})}} s_{i j}^{*}],

(115)

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 $l^{R}$ and $l^{I}$ are real by definition. Using this definition of C_Loc, Eq. (91) becomes

\sum_{k < l} l_{i j}^{R} A_{i j, k l}^{l o c, R} + i \sum_{k < l} l_{k l}^{I} A_{i j, k l}^{l o c, I} = - L_{i j},

(116)

where

A_{i j, k l}^{l o c, R} = [\frac{d s_{k l}}{d U_{i j}} + {(\frac{d s_{l k}}{d U_{j i}})}^{*} - \frac{d s_{l k}}{d U_{i j}} - {(\frac{d s_{k l}}{d U_{j i}})}^{*}],

(117)

A_{i j, k l}^{l o c, I} = [\frac{d s_{k l}}{d U_{i j}} + {(\frac{d s_{l k}}{d U_{j i}})}^{*} + \frac{d s_{l k}}{d U_{i j}} + {(\frac{d s_{k l}}{d U_{j i}})}^{*}] .

(118)

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

\sum_{i < j} A_{k l, i j}^{l o c, R} u_{i j, n^{^{″}}}^{R} = u_{k l, n^{^{″}}}^{R} ω_{n^{^{″}}}^{R},

(119)

\sum_{i < j} A_{k l, i j}^{l o c, I} u_{i j, n^{^{″}}}^{I} = u_{k l, n^{^{″}}}^{I} ω_{n^{^{″}}}^{I},

(120)

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

\begin{align} C_{Loc}^{valence} & = \sum_{n}^{R} l_{n}^{R} R [s_{n}] + \sum_{n}^{I} l_{n}^{I} I [s_{n}] + \sum_{n^{'}}^{R} l_{n^{'}}^{R} R [U_{n^{'}}] + \sum_{n^{'}}^{I} l_{n^{'}}^{I} I [U_{n^{'}}] \\ = \frac{1}{2} \sum_{i < j}^{valence} \sum_{k < l}^{valence} [l_{i j}^{R} (δ_{i j, k l} - O_{i j, k l}^{R}) (s_{k l} + s_{k l}^{*})) \\ + (i l_{i j}^{I} (δ_{i j, k l} - O_{i j, k l}^{I}) (s_{k l} - s_{k l}^{*})) \\ (+ l_{i j}^{R} O_{i j, k l}^{R} (U_{k l} + U_{k l}^{*}) + i l_{i j}^{I} O_{i j, k l}^{I} (U_{k l} - U_{k l}^{*})], \end{align}

(121)

where δ_ij,kl = δ_ikδ_jl and the projectors onto the nullspace of $A^{l o c, R / I}$ are

O_{i j, k l}^{R / I} = \sum_{n^{'}} u_{i j, n^{'}}^{R / I} u_{k l, n^{'}}^{R / I} .

(122)

Thus, the Z-CPL equations become

\sum_{k < l} l_{k l}^{R} Ã_{i j, k l}^{l o c, R} + i \sum_{k < l} l_{k l}^{I} Ã_{i j, k l}^{l o c, I} = - L_{i j},

(123)

where

{\tilde{A}}^{l o c, R} = (I - 𝒪^{R}) A^{l o c, R} (I - 𝒪^{R}) + 𝒪^{R},

(124)

{\tilde{A}}^{l o c, I} = (I - 𝒪^{I}) A^{l o c, I} (I - 𝒪^{I}) + 𝒪^{I} .

(125)

Once again, $l^{I}$ is zero, while the contributions of $l^{R}$ to Eqs. (68), (73), (112), and (113) must be calculated with the projected matrix,

{\tilde{l}}_{i j}^{R} = \sum_{k < l} (δ_{i j, k l} - O_{i j, k l}^{R}) l_{k l}^{R} .

(126)

The electric/magnetic CPL equations become

\sum_{k < l} U_{k l}^{λ} Ã_{k l, i j}^{l o c, R / I} = - \sum_{k < l} {(I - O^{R / I})}_{i j, k l} B_{k l}^{l o c, λ},

(127)

\begin{align} B_{i j}^{l o c, λ} & = s_{i j}^{(λ)} - \sum_{k < l} [\frac{d s_{i j}}{d U_{l k}} S_{l k}^{(λ)} - {(\frac{d s_{j i}}{d U_{l k}})}^{*} S_{l k}^{(λ) *}] \\ + \sum_{p k} (δ_{p m} + δ_{p a}) [\frac{d s_{i j}}{d U_{p k}} U_{p k}^{λ} - {(\frac{d s_{j i}}{d U_{p k}})}^{*} U_{p k}^{λ *}], \end{align}

(128)

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,

Ũ_{i j}^{λ} = \sum_{k < l} (δ_{i j, k l} - O_{i j, k l}^{R / I}) U_{k l}^{λ} .

(129)

The perturbed Z-CPL equations become

\sum_{k < l} l_{k l}^{R, λ} Ã_{i j, k l}^{l o c, R} + i \sum_{k < l} l_{k l}^{I, λ} Ã_{i j, k l}^{l o c, I} = - \sum_{k < l} {\tilde{l}}_{k l}^{R} A_{i j, k l}^{l o c, R, λ} - L_{i j}^{λ},

(130)

with

\begin{align} A_{p q, i j}^{l o c, R, λ} & = [(r_{i j} + r_{i j}^{*}) (δ_{p j} r_{q j} + δ_{q i} r_{i p} - δ_{p i} r_{q i} - δ_{q j} r_{j p})) \\ {(+ (r_{i i} - r_{j j}) (δ_{q j} r_{i p} + δ_{q i} r_{j p} - δ_{p i} r_{q j} - δ_{p j} r_{q i})]}^{λ} . \end{align}

(131)

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, ${(r_{i i} - r_{j j})}^{λ} = 0$ and ${(r_{i j} + r_{i j}^{*})}^{λ} = 0$ ⁠. The contributions of $l^{R / I, λ}$ to Eqs. (112) and (113) must be calculated with the projected matrices

{\tilde{l}}_{i j}^{R / I, λ} = \sum_{k < l} (δ_{i j, k l} - O_{i j, k l}^{R / I}) l_{k l}^{R / I, λ} .

(132)

Finally, we demonstrate the numerical instabilities that can occur if the localization response matrix, $A^{l o c, R}$ or $A^{l o c, I}$ ⁠, is singular, as well as the effectiveness of the procedure to remove the near-zero eigenvalues. We chose the octahedral molecules SF₆ and SeF₆ with bond lengths 1.560 and 1.688 Å, respectively. SF₆ has six eigenvalues of $A^{l o c, R}$ smaller than 10⁻⁴ due to a continuous degeneracy in the FB LMOs among the 18 fluorine sp³-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₆ 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 $A^{l o c, I}$ only has ten singular eigenvalues for SeF₆, due to the 3d orbitals, and none for SF₆. 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.

TABLE II.

(Maximum) absolute errors vs RI-MP2 in the isotropic polarizability (α) and ¹⁹F shielding (σ_F) for the systems SF₆ and SeF₆ from frozen-core (FC) and all-electron (AE) NormalPNO calculations with and without projection of the localization response nullspace. The first two rows show the number of near-zero eigenvalues of the localization response matrices, see text for details.

System	SF₆		SeF₆
	FC	AE	FC	AE
$A^{l o c, R}$ nullspace	6	9	16	24
$A^{l o c, I}$ nullspace	0	0	10	10
	Errors without projection:
\|Δα\| (Å³)	8 × 10⁵	…^a	17	…^a
max (\|Δσ_F\|) (ppm)	0.26	0.28	2.1	29
	Errors with projection:
\|Δα\| (Å³)	<0.001	<0.001	<0.001	<0.001
max (\|Δσ_F\|) (ppm)	0.26	0.28	0.47	0.55

System	SF₆		SeF₆
	FC	AE	FC	AE
$A^{l o c, R}$ nullspace	6	9	16	24
$A^{l o c, I}$ nullspace	0	0	10	10
	Errors without projection:
\|Δα\| (Å³)	8 × 10⁵	…^a	17	…^a
max (\|Δσ_F\|) (ppm)	0.26	0.28	2.1	29
	Errors with projection:
\|Δα\| (Å³)	<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⁻⁸ 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 $s_{n}^{- \frac{1}{2}} \approx 1 0^{4}$ ⁠. 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

d^{'' i j, λ} = U^{S} s^{- \frac{1}{2}} U^{S D} θ^{i j, λ},

(133)

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

\begin{align} θ^{λ} & = - \frac{1}{2} d^{'' (0) †} {\tilde{S}}^{λ} d^{'' (0)} \\ = - \frac{1}{2} U^{S D †} s^{- \frac{1}{2}} U^{S †} {\tilde{S}}^{λ} U^{S} s^{- \frac{1}{2}} U^{S D} \propto s^{- 1} \approx 1 0^{8}, \end{align}

(134)

and therefore,

d^{'' λ} \propto s^{- \frac{3}{2}} \approx 1 0^{12},

(135)

The largest elements of the matrices are estimated based on $s^{- \frac{1}{2}} \approx T_{S}^{- \frac{1}{2}} = 1 0^{4}$ ⁠. 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⁻². Ignoring the line “Analytic v2” for now, the large error in the analytic polarizability at T_S = 10⁻⁸ is apparent. This disappears at higher values of the cutoff, and at T_S = 10⁻⁵, there is very good agreement with both the semi-numeric polarizability (calculated from analytic dipole moment calculations with an explicit electric field of 10⁻⁴ 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.

FIG. 1.

View large Download slide

Eigenvalues of the PAO domain overlap matrix for different domains in a pentane molecule with the def2-TZVP basis set. The atoms included in each domain are depicted graphically with matching colors.

TABLE III.

Isotropic polarizability (bohr³) of the water dimer model (see text) calculated semi-numerically and analytically using RI-MP2 and DLPNO-MP2 with the “symmetric” solution for $θ_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j}$ (v1) or the alternative solution discussed in the text (v2). Only the domain truncation was used. The maximal absolute deviation of the MP2 density matrix elements [ max(|ΔD_μν|)] is also given.

	T_S = 10⁻⁸	T_S = 10⁻⁷	T_S = 10⁻⁶	T_S = 10⁻⁵	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⁻⁸	T_S = 10⁻⁷	T_S = 10⁻⁶	T_S = 10⁻⁵	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

θ^{λ} = - \frac{1}{2} d^{'' (0) †} {\tilde{S}}^{λ} d^{'' (0)} + θ^{1, λ} + θ^{2, λ},

(136)

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,

{|U^{S} s^{- \frac{1}{2}} U^{S D} (- \frac{1}{2} d^{'' (0) †} {\tilde{S}}^{λ} d^{'' (0)} + θ^{1, λ})|}^{2} \to min,

(137)

which is solved in the overlap eigenbasis as

{(U^{S D} θ^{1, λ} U^{S D †})}_{\tilde{μ} \tilde{ν}} = \frac{(s_{\tilde{ν}} - s_{\tilde{μ}})}{2 \sqrt{s_{\tilde{μ}} s_{\tilde{ν}}} (s_{\tilde{ν}} + s_{\tilde{μ}})} {(U^{S †} {\tilde{S}}^{λ} U^{S})}_{\tilde{μ} \tilde{ν}} .

(138)

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⁻⁵, 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 ${\overset{⌣}{D}}_{ã {\tilde{b}}^{'}}^{i j}$ ⁠, it is necessary to divide by the difference $(n_{ã}^{i j} - n_{{\tilde{b}}^{'}}^{i j})$ 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 ${\overset{⌣}{D}}^{i j}$ 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_{ã}^{i j} - n_{{\tilde{b}}^{'}}^{i j}) / n_{ã}^{i j}$ was smaller than 10⁻², 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_{ã}^{i j} - n_{{\tilde{b}}^{'}}^{i j} + ε) > ε$ ⁠. To that end, we modify the PNO condition as follows:

C_{PNO} = \frac{1}{2} \sum_{i \geq j} \sum_{ã {\tilde{b}}^{'}} [v_{ã {\tilde{b}}^{'}}^{i j} ({\overset{⌣}{D}}_{ã {\tilde{b}}^{'}}^{i j} + Ξ_{ã {\tilde{b}}^{'}}^{i j}) + v_{ã {\tilde{b}}^{'}}^{i j *} ({\overset{⌣}{D}}_{ã {\tilde{b}}^{'}}^{i j *} + Ξ_{ã {\tilde{b}}^{'}}^{i j *})],

(139)

Ξ^{i j} = θ^{i j *} Ξ^{i j} θ^{i j, T},

(140)

Ξ_{ã^{^{″}} {\tilde{b}}^{^{″}}}^{i j} = δ_{\tilde{c} ã^{^{″}}} δ_{\tilde{c} {\tilde{b}}^{^{″}}} ε .

(141)

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),

\begin{align} (n_{ã}^{i j} - n_{{\tilde{b}}^{'}}^{i j} + ε) θ_{{\tilde{b}}^{'} ã}^{i j, λ} \\ = - (n_{ã}^{i j} + ε) S_{{\tilde{b}}^{'} ã}^{i j, λ} - ({\tilde{\overset{⌢}{T}}}^{i j *} S^{i j, λ} {\overset{⌢}{T}}^{j i} + {\tilde{\overset{⌢}{T}}}^{j i *} S^{i j, λ} {\overset{⌢}{T}}^{i j}) \\ - {({\tilde{\overset{⌢}{T}}}^{i j, λ *} {\overset{⌢}{T}}^{j i} - {\tilde{\overset{⌢}{T}}}^{j i, λ *} {\overset{⌢}{T}}^{i j} - {\tilde{\overset{⌢}{T}}}^{i j *} {\overset{⌢}{T}}^{j i, λ} - {\tilde{\overset{⌢}{T}}}^{j i *} {\overset{⌢}{T}}^{i j, λ})}_{ã {\tilde{b}}^{'}}, \end{align}

(142)

(n_{ã}^{i j} - n_{{\tilde{b}}^{'}}^{i j} + ε) v_{ã {\tilde{b}}^{'}}^{i j, λ} = - V_{{\tilde{b}}^{'} ã}^{i j, λ} + ε \sum_{\tilde{c}} S_{ã \tilde{c}}^{i j, λ} v_{\tilde{c} {\tilde{b}}^{'}}^{i j} .

(143)

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⁻⁵ and T_CutPNO = 10⁻⁷, as well as with T_CutPNO = 5.96 × 10⁻⁷, 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⁻⁷, 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.

FIG. 2.

View large Download slide

Maximum absolute errors vs RI-MP2 in the shieldings of caffeine from LoosePNO calculations with varying T_CutPNO and different level shift scaling parameters ε_scale in the PNO response equations (ε_scale = 0 is equivalent to no level shift). The legend applies to both subplots.

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,¹⁰⁹ closely following the algorithm for the DLPNO-MP2 gradient.⁶⁹ 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 $U_{a i}^{λ}$ 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 ${\tilde{P}}^{λ}$ ⁠.
The three-index integrals $(i {\tilde{μ}}^{'} | K)$ and ${(i {\tilde{μ}}^{'} | K)}^{λ}$ are calculated and stored.
In a loop over orbital pairs, the CPNO coefficients d′^ij are reconstructed, and ${\overset{⌣}{T}}^{i j}$ and ${\overset{⌢}{T}}^{i j, λ}$ 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, ${\overset{⌣}{T}}^{i j}$ ⁠, d″^ij,λ, and ${\overset{⌣}{T}}^{i j, λ}$ 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^δ,λ.
Γ^K and Γ^K,λ are read from disk, and their contributions to $L_{i j}^{λ}$ ⁠, $L_{m j}^{λ}$ ⁠, Y^α,λ, Y^β, and Y^β,λ are evaluated in the AO basis by generating the required three-index integrals on the fly.
The prescreening contributions to $L_{i j}^{λ}$ ⁠, $L_{m j}^{λ}$ ⁠, $L_{i a}^{λ}$ ⁠, and D′ are calculated.
$L_{m j}^{λ}$ is completed, and the Z-CV equations (112) are solved to obtain $z_{m j}^{λ}$ ⁠.
$L_{i j}^{λ}$ is completed, and the Z-CPL equations (109) are solved to obtain $z_{i j}^{λ}$ ⁠. In our implementation, a conjugate gradient solver is used.
$L_{i a}^{λ}$ is completed, and the Z-CPSCF equations (69) are solved to obtain $z_{i a}^{λ}$ ⁠. 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 {\tilde{μ}}^{'} | K)}^{λ}$ ⁠, 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.¹¹⁰ The calculation of terms involving (K|ip) integrals in the Z-CPSCF equations’ [Eq. (74)] right-hand sides asymptotically scales as $O (N^{2})$ ⁠. Together with the $(i {\tilde{μ}}^{'} | K)$ 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 (N^{3})$ ⁠, but these steps constitute a very small part of the overall computation time.

Alongside MP2, the implementation can be used for DHDFT—the relevant modifications are identical to the RI-MP2 case, discussed in Ref. 72. Spin-component scaling is also implemented.¹¹¹ Implicit solvent effects can be included using the conductor-like polarizable continuum model (CPCM).¹¹²

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.

FIG. 3.

View large Download slide

Structures of the benchmark systems with the number of electrons and pcSseg-2 basis functions.

For the shielding calculations, the pcSseg-2¹¹³ 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.⁷² 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¹¹⁷ 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¹¹⁴ auxiliary basis set. Once again, def2/JK¹¹⁷ 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⁻¹² for the core orbitals. We set F_Cut = 10⁻⁸, T_S = 10⁻⁸, and all other DLPNO thresholds to 0.
Varying T_CutPNO for the core orbitals with T_CutPNO = 10⁻¹⁰ for the valence orbitals. We set F_Cut = 10⁻⁸, T_S = 10⁻⁸, and all other DLPNO thresholds equal to 0.
Varying F_Cut with T_CutPNO = 10⁻¹² for the core orbitals and T_CutPNO = 10⁻¹⁰ for the valence orbitals. We set T_S = 10⁻⁸ and all other DLPNO thresholds equal to 0.
Same as above but also fixing F_Cut = 10⁻⁵ and varying T_CutDO. Two values were used for T_S to examine the numerical stability: 10⁻⁸ and 10⁻⁵.
Same as above but also fixing T_CutDO = 10⁻³, T_CutDOPre = 0.03, and T_S = 10⁻⁵ 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⁻⁸ (10⁻¹⁰ for core orbitals), F_Cut = 10⁻⁵, T_CutDO = 10⁻², T_CutDOPre = 0.03, T_CutDOij = 10⁻⁵, T_CutPre = 10⁻⁶, T_CutMKN = 10⁻³, T_CutC = 10⁻³, but with T_S = 10⁻⁵ and a PNO level shift ε_scale = 0.1.
Same as above with the “LoosePNO” default thresholds: T_CutPNO = 10⁻⁷ (10⁻⁹ for core orbitals) and T_CutDO = 2 × 10⁻².
We also performed shielding calculations with RI-MP2, NormalPNO, LoosePNO, and TightPNO (T_CutPNO = 10⁻⁹ for valence and 10⁻¹¹ for core orbitals, T_CutDO = 5 × 10⁻³), 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.¹²⁰ 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,¹²¹ 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¹¹⁴ 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⁴⁻, caffeine, ebselen, and penicillin. We note that the errors diminish rapidly with decreasing T_CutPNO, and at the default value of 10⁻⁸ 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.⁷²

FIG. 4.

View large Download slide

Mean absolute error (a) and mean absolute weighted error (b) vs RI-MP2 in the shieldings of different nuclei in the systems ATP⁴⁻, caffeine, ebselen, and penicillin. T_CutPNO was varied for either the valence (left) or core (right) orbitals, with the other parameter set to 10⁻¹² and 10⁻¹⁰, respectively. The line and marker styles in the legend apply to all subplots.

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,

AWE = \frac{|σ_{i}^{DLPNO - MP 2} - σ_{i}^{RI - MP 2}|}{w_{A}} \times 100 %, i \in A,

(144)

w_{A} = \frac{1}{2} (δ_{max}^{A} - δ_{min}^{A}) .

(145)

We use the values for $δ_{min}^{A}$ and $δ_{max}^{A}$ shown in Table IV. A reasonable target accuracy is a mean AWE (MAWE) below 0.5%, which corresponds to 0.5 ppm for ¹³C, 0.03 ppm for ¹H, and about 1 ppm for ¹⁵N and ¹⁷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⁴⁻ 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⁻⁹ with MaxAE = 1.51 ppm for ¹⁵N vs 0.59 and 0.04 ppm for the same nucleus at T_CutPNO = 10⁻⁸ and T_CutPNO = 10⁻¹⁰, respectively. With the level shift-based adjustment (see Sec. III C), this error drops to 0.17 ppm at valence T_CutPNO = 10⁻⁹.

TABLE IV.

Experimental chemical shift ranges (ppm) for different nuclides used for scaling in MAWE calculations.

A	$δ_{min}^{A}$	$δ_{max}^{A}$	w_A
¹H	−1	12	6.5
¹³C	0	200	100
¹⁵N	0	900	450
¹⁷O	−40	1120	580
¹⁹F	−300	400	350
³³S	−290	670	480
³¹P	−180	250	215
⁷⁷Se	−1000	2000	1500

A	$δ_{min}^{A}$	$δ_{max}^{A}$	w_A
¹H	−1	12	6.5
¹³C	0	200	100
¹⁵N	0	900	450
¹⁷O	−40	1120	580
¹⁹F	−300	400	350
³³S	−290	670	480
³¹P	−180	250	215
⁷⁷Se	−1000	2000	1500

FIG. 5.

View large Download slide

Mean absolute weighted error vs RI-MP2 in the shieldings of the systems ATP⁴⁻, caffeine, ebselen, and penicillin. T_CutPNO was varied for either the valence (top) or core (bottom) orbitals, with the other parameter set to 10⁻¹² and 10⁻¹⁰, respectively. The line and marker styles in the legend apply to both subplots.

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⁻⁵.

FIG. 6.

View large Download slide

Mean absolute weighted error vs RI-MP2 in the shieldings of the systems ATP⁴⁻, caffeine, ebselen, and penicillin with varying F_Cut.

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⁻⁸, 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⁻⁵ was used, which leads to smooth convergence or the results toward the RI-MP2 reference. At the default value of T_CutDO = 10⁻², the MAWE over the whole dataset is under 0.05%.

FIG. 7.

View large Download slide

Mean absolute weighted error vs RI-MP2 in the shieldings of the systems ATP⁴⁻, caffeine, ebselen, and penicillin. T_CutDO was varied, and T_S was set to either 10⁻⁸ (left) or 10⁻⁵ (right). The line and marker styles in the legend apply to both subplots.

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⁻⁴ and 20–50% at T_CutPre = 10⁻⁶. 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⁻⁴; 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⁻⁵ and T_CutPre = 10⁻⁶ might be too conservative for NMR shielding calculations.

FIG. 8.

View large Download slide

Mean absolute weighted error vs RI-MP2 in the shieldings of the systems ATP⁴⁻, caffeine, ebselen, and penicillin. T_CutPre was varied, and dipole estimate-based contributions from the screened-out pairs were either included (solid lines) or excluded (dotted lines).

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.

FIG. 9.

View large Download slide

Distribution of errors in the shieldings for the systems: (anth)₂, ATP⁴⁻, caffeine, coronene, ebselen, penicillin, and tweezer complex. Calculations were performed (from left to right, respectively) using LoosePNO thresholds without and with the COSX approximation, using NormalPNO thresholds without and with COSX, and using TightPNO and RI-MP2 with COSX, and the errors were calculated against RI-MP2 without COSX for all datasets. The lines, boxes, and whiskers denote the median, interquartile range, and minimum/maximum values, respectively.

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 ¹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)₁₀, and are within the ranges shown in Fig. 9. A notable outlier is an oxygen nucleus in linear (gly)₃ with a rather large error of 0.97 ppm for ¹⁷O vs 0.25 and 0.41 ppm for (gly)₂ and (gly)₄, respectively. Using the level shift-based adjustment, discussed in Sec. III C, the MaxAE for ¹⁷O in linear (gly)₃ is reduced to 0.42 ppm, consistent with the slightly oscillating trend.

FIG. 10.

View large Download slide

Maximum absolute errors in the shieldings of different nuclides from NormalPNO calculations vs RI-MP2 for polyglycine chains in either linear (dotted lines) or α-helix conformation (solid lines). RIJCOSX was used throughout.

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,

RE = \frac{α^{DLPNO - MP 2} - α^{RI - MP 2}}{α^{RI - MP 2}} \times 100 %,

(146)

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.⁷³ 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⁻⁷. 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⁻¹² and two for ATP⁴⁻ at 10⁻⁷ and 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⁻⁵ for these calculations, the analogous threshold for the AO basis was 10⁻⁸ and the lowest overlap eigenvalues for caffeine, penicillin, coronene, and ATP⁴⁻ were 2.7 × 10⁻⁷, 4.2 × 10⁻⁷, 2.2 × 10⁻⁸, and 4.6 × 10⁻⁷, 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⁻⁶ for the MO orthonormalization or applying the level shift, introduced in Sec. III C, both eliminate these discrepancies (see the supplementary material).

FIG. 11.

View large Download slide

Absolute relative errors vs RI-MP2 in the isotropic polarizability of the systems ATP⁴⁻, caffeine, coronene, and penicillin. T_CutPNO was varied for either the valence (left) or core (right) orbitals, with the other parameter set to 10⁻¹² and 10⁻¹⁰, respectively. The line and marker styles in the legend apply to both subplots.

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⁻⁵, however, convergence toward the RI-MP2 reference is much smoother, and the errors are already sufficiently small at T_CutDO = 3 × 10⁻².

FIG. 12.

View large Download slide

Absolute relative errors vs RI-MP2 in the isotropic polarizability of the systems ATP⁴⁻, caffeine, coronene, and penicillin. T_CutDO was varied, and T_S was set to either 10⁻⁸ (top) or 10⁻⁵ (bottom). The line and marker styles in the legend apply to both subplots.

The influence of F_Cut on the polarizability, shown in Fig. 13, is rather small. Apparently, the default value of 10⁻⁵ is more than conservative enough.

FIG. 13.

View large Download slide

Absolute relative errors vs RI-MP2 in the isotropic polarizability of the systems ATP⁴⁻, caffeine, coronene, and penicillin with varying F_Cut.

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⁴⁻; however with the high value of T_CutDOij used here, between 35% and 67% of pairs are screened out already at T_CutPre = 10⁻⁵, 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⁻⁴ are probably due to the strong delocalization in the system, together with the large value of T_CutDOij.

FIG. 14.

View large Download slide

Absolute relative errors vs RI-MP2 in the isotropic polarizability of the systems ATP⁴⁻, caffeine, coronene, and penicillin. T_CutPre was varied, and dipole estimate-based contributions from the screened-out pairs were either included (solid lines) or excluded (dotted lines).

The combined effect of all DLPNO approximations with the default NormalPNO parameters (and T_S = 10⁻⁵, ε_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),⁷³ 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.

FIG. 15.

View large Download slide

Relative errors in the isotropic polarizability from NormalPNO calculations vs RI-MP2 for the systems caffeine, ebselen, penicillin, coronene, ATP⁴⁻, and (anth)₂, ordered by increasing polarizability.

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)₅ (38 atoms, 160 electrons, and 994 basis functions) and (gly)₉ (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)₁₅ (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)₅] 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)₁₅. 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.

FIG. 16.

View large Download slide

Wall-clock times (on 8 CPU cores) for shielding calculations on polyglycine chains in either linear (dotted lines) or α-helix (solid lines) conformation using either DLPNO-MP2 with NormalPNO or RI-MP2. RIJCOSX was used throughout.

FIG. 17.

View large Download slide

Ratio of the wall-clock times (on 8 CPU cores) for MP2 vs HF shielding calculations on polyglycine chains in either linear (dotted lines) or α-helix (solid lines) conformation using either DLPNO-MP2 with NormalPNO or RI-MP2. RIJCOSX was used throughout.

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.⁷² (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.

TABLE V.

System size indicators and wall-clock times (min) for a frozen-core DLPNO-DSD-PBEP86 NMR shielding calculation on vancomycin using RIJCOSX. The calculation was performed on 12 Intel Xeon E5-2687Wv4@3.00GHz CPU cores with 15 GB RAM per core. Subparts of calculation steps are marked with “→.”

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
$\to (i {\tilde{μ}}^{'} \| K)$ transformation	23
$\to {\overset{⌣}{T}}^{i j}$ and d^ij construction	88
→PNO amplitude iterations	35
→NPAO space reconstruction	163
→PNO-specific terms	70
$\to Γ_{i {\tilde{μ}}^{'}}^{K}$ contributions	84
→Fock response terms	10
→Z-CPSCF (10 cycles)	56
→ Localization + Z-CPL	5
DLPNO-MP2 response density	3 603
$\to {(i {\tilde{μ}}^{'} \| K)}^{λ}$ transformation	276
$\to {\overset{⌢}{T}}^{i j, λ}$ and d″^ij,λ construction	809
→T^ij,λ iterations	160
→NPAO space reconstruction	270
→PNO-specific terms	671
$\to Γ_{i {\tilde{μ}}^{'}}^{K}$ and $Γ_{i {\tilde{μ}}^{'}}^{K, λ}$ 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
$\to (i {\tilde{μ}}^{'} \| K)$ transformation	23
$\to {\overset{⌣}{T}}^{i j}$ and d^ij construction	88
→PNO amplitude iterations	35
→NPAO space reconstruction	163
→PNO-specific terms	70
$\to Γ_{i {\tilde{μ}}^{'}}^{K}$ contributions	84
→Fock response terms	10
→Z-CPSCF (10 cycles)	56
→ Localization + Z-CPL	5
DLPNO-MP2 response density	3 603
$\to {(i {\tilde{μ}}^{'} \| K)}^{λ}$ transformation	276
$\to {\overset{⌢}{T}}^{i j, λ}$ and d″^ij,λ construction	809
→T^ij,λ iterations	160
→NPAO space reconstruction	270
→PNO-specific terms	671
$\to Γ_{i {\tilde{μ}}^{'}}^{K}$ and $Γ_{i {\tilde{μ}}^{'}}^{K, λ}$ 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.⁵⁹ 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 (N_{bas}^{4.5})$ ⁠, roughly equivalent to the $O (N_{bas}^{4.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 (N_{bas}^{2.8})$ ⁠.

TABLE VI.

Isotropic polarizability (α_iso/bohr³) and wall-clock time (t/h) for calculations on [n]helicene using DLPNO-MP2, compared to RI-MP2 data from Ref. 59. The number of (aug-cc-pVTZ) basis functions (N_bas) is also given.

	DLPNO-MP2		RI-MP2⁵⁹
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⁵⁹
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 ¹³C and ¹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)
$ã, \tilde{b}, \tilde{c}, \tilde{d}$: pair natural orbitals (PNOs)
$ã^{'}, {\tilde{b}}^{'}, {\tilde{c}}^{'}, {\tilde{d}}^{'}$: discarded/“complementary” pair natural orbitals (CPNOs)
$ã^{^{″}}, {\tilde{b}}^{^{″}}, {\tilde{c}}^{^{″}}, {\tilde{d}}^{^{″}}$: union of PNOs and CPNOs
$\tilde{μ}, \tilde{ν}, \tilde{η}$: orbital/pair domain of orthonormal pseudo-canonical non-redundant projected atomic orbitals (NPAOs)
${\tilde{μ}}^{'}, {\tilde{ν}}^{'}, {\tilde{η}}^{'}$: 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.⁶⁹ 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:

0 \equiv \frac{\partial L}{\partial θ^{i j}} = \frac{\partial E_{2}^{D L P N O}}{\partial θ^{i j}} + \frac{\partial Δ E_{P N O}}{\partial θ^{i j}} + \frac{\partial C_{PNO}}{\partial θ^{i j}} + \frac{\partial C_{PNOO}}{\partial θ^{i j}} + \frac{\partial C_{SC}}{\partial θ^{i j}},

(A1)

0 \equiv \frac{\partial L}{\partial θ^{i j †}} = {(\frac{\partial L}{\partial θ^{i j}})}^{†} .

(A2)

The separate terms are as follows:

\frac{\partial E_{2}^{D L P N O}}{\partial θ^{i j}} = 𝒦^{j i} θ^{i j} {\tilde{T}}^{i j} + 𝒦^{i j} θ^{i j} {\tilde{T}}^{j i} + F^{i j *} θ^{i j *} D^{i j} - d^{'' i j (0), T} {\tilde{S}}^{*} τ^{i j †},

(A3)

\begin{align} \frac{\partial Δ E_{P N O}}{\partial θ^{i j}} & = \frac{1}{2} (K^{j i} θ^{i j} {\tilde{\overset{⌣}{T}}}^{i j} + K^{i j} θ^{i j} {\tilde{\overset{⌣}{T}}}^{j i} + {\tilde{\overset{⌢}{T}}}^{i j *} θ^{i j} K^{j i *} + {\tilde{\overset{⌢}{T}}}^{j i *} θ^{j i} K^{i j *}) \\ - \frac{1}{2} (K^{j i} θ^{i j} N^{i j} {\tilde{\overset{⌣}{T}}}^{i j} + K^{i j} θ^{i j} N^{i j} {\tilde{\overset{⌣}{T}}}^{j i} + {\tilde{\overset{⌢}{T}}}^{i j *} θ^{i j} N^{i j} K^{j i *}) \\ + ({\tilde{\overset{⌢}{T}}}^{j i *} θ^{j i} N^{i j} K^{i j *}) N^{i j}, \end{align}

(A4)

\frac{\partial C_{PNO}}{\partial θ^{i j}} = \frac{1}{2} {\tilde{\overset{⌢}{T}}}^{i j *} θ^{i j} ({\overset{⌣}{T}}^{j i} v^{i j, T} + v^{i j} {\overset{⌣}{T}}^{j i}) + \frac{1}{2} {\tilde{\overset{⌢}{T}}}^{j i *} θ^{i j} ({\overset{⌣}{T}}^{i j} v^{i j, T} + v^{i j} {\overset{⌣}{T}}^{i j}),

(A5)

\frac{\partial C_{PNOO}}{\partial θ^{i j}} = 𝒮^{i j *} θ^{i j *} x^{i j},

(A6)

\begin{align} \frac{\partial C_{SC}}{\partial θ^{i j}} & = \frac{1}{2} [F^{i j *} θ^{i j *} w^{i j, T} {\overset{⌣}{T}}^{i j} + F^{i j *} θ^{i j *} w^{i j} {\overset{⌣}{T}}^{j i} + {\overset{⌢}{T}}^{j i *} θ^{i j} w^{i j *} F^{i j *}) \\ + ({\overset{⌢}{T}}^{i j *} θ^{i j} w^{i j †} F^{i j *} + {\overset{⌢}{R}}^{i j *} θ^{i j} w^{i j †} + {\overset{⌢}{R}}^{i j †} θ^{i j} w^{i j *}], \end{align}

(A7)

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

K^{i j} = θ^{i j, T} K^{i j} θ^{i j},

(A8)

{\overset{⌣}{R}}^{i j} = θ^{i j †} {\overset{⌢}{R}}^{i j} θ^{i j *},

(A9)

S^{i j} = θ^{i j †} S^{i j} θ^{i j} = θ^{i j †} d^{'' i j (0) †} \tilde{S} d^{'' i j (0)} θ^{i j},

(A10)

F^{i j} = θ^{i j †} F^{i j} θ^{i j} .

(A11)

The index ${\tilde{μ}}^{'}$ of the intermediate $τ_{ã {\tilde{μ}}^{'}}^{i j}$ spans the PAO domains of all pairs ik and jk included in the sum, i.e., those for which $|F_{k j}|$ or $|F_{k i}|$ ⁠, 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:

\begin{align} 0 & \equiv {([\frac{\partial L}{\partial θ^{i j}} - \frac{\partial L}{\partial θ^{i j †}}]|}_{λ = 0} = [\frac{\partial E_{2}^{D L P N O}}{\partial θ^{i j}} - \frac{\partial E_{2}^{D L P N O}}{\partial θ^{i j †}}) \\ {(+ (\frac{\partial Δ E_{P N O}}{\partial θ^{i j}} - \frac{\partial Δ E_{P N O}}{\partial θ^{i j †}} + \frac{\partial C_{PNO}}{\partial θ^{i j}} - \frac{\partial C_{PNO}}{\partial θ^{i j †}}]|}_{λ = 0} . \end{align}

(A12)

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

\begin{align} \frac{\partial C_{PNO}}{\partial θ^{i j}} & = \frac{1}{2} [{\tilde{\overset{⌢}{T}}}^{i j *} θ^{i j} (N^{i j} K^{j i *} N^{i j} - K^{j i *})) \\ + ({\tilde{\overset{⌢}{T}}}^{j i *} θ^{i j} (N^{i j} K^{i j *} N^{i j} - K^{i j *})) \\ + ({\overset{⌢}{T}}^{i j *} θ^{i j} (w^{i j †} (F_{j j} + F_{i i}) - F^{i j} w^{i j †} - w^{i j †} F^{i j *})) \\ + ({\overset{⌢}{T}}^{j i *} θ^{i j} (w^{i j *} (F_{j j} + F_{i i}) - F^{i j} w^{i j *} - w^{i j *} F^{i j *})], \end{align}

(A13)

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

\begin{align} 0 & \equiv \frac{\partial L}{\partial θ^{i j}} = \frac{1}{2} K^{j i} θ^{i j} g^{i j} + \frac{1}{2} K^{i j} θ^{i j} g^{i j, T} + F^{i j *} θ^{i j *} D^{'' i j} \\ - d^{'' i j (0), T} {\tilde{S}}^{*} τ^{i j †} + S^{i j *} θ^{i j *} x^{i j}, \end{align}

(A14)

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

D^{'' i j} = D^{i j} + \frac{1}{2} (w^{i j} {\overset{⌣}{T}}^{j i} + w^{i j, T} {\overset{⌣}{T}}^{i j} + {\overset{⌣}{T}}^{j i *} w^{i j *} + {\overset{⌣}{T}}^{i j *} w^{i j †}) .

(A15)

If we define $L_{- PNOO} = L - C_{PNOO}$ ⁠, we obtain a solution for x^ij,

\begin{align} x^{i j} & = - {(\frac{\partial L_{- PNOO}}{\partial θ^{i j}}|}_{λ = 0} \\ = - \frac{1}{2} K^{j i} g^{i j} - \frac{1}{2} K^{i j} g^{i j, T} - F^{i j *} D^{'' i j} + d^{'' i j (0), T} {\tilde{S}}^{*} τ^{i j †} . \end{align}

(A16)

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

\begin{align} x^{i j} & = - {(\frac{\partial L_{- PNOO}}{\partial θ^{i j †}}|}_{λ = 0} \\ = - \frac{1}{2} g^{i j †} K^{i j *} - \frac{1}{2} g^{i j *} K^{j i *} - D^{'' i j} F^{i j *} + τ^{i j} {\tilde{S}}^{*} d^{'' i j (0) *}, \end{align}

(A17)

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

\begin{align} 0 & \equiv \frac{d}{d λ} \frac{\partial L}{\partial θ^{i j}} = \frac{d}{d λ} (\frac{\partial L_{- PNOO}}{\partial θ^{i j}} + S^{i j *} θ^{i j *} x^{i j}) \\ = \frac{d}{d λ} \frac{\partial L_{- PNOO}}{\partial θ^{i j}} + θ^{i j, λ, T} \frac{\partial L_{- PNOO}}{\partial θ^{i j}} + x^{i j, λ}, \end{align}

(A18)

0 \equiv \frac{d}{d λ} \frac{\partial L}{\partial θ^{i j †}} = \frac{d}{d λ} \frac{\partial L_{- PNOO}}{\partial θ^{i j †}} + \frac{\partial L_{- PNOO}}{\partial θ^{i j †}} θ^{i j, λ *} + x^{i j, λ},

(A19)

where we have used Eqs. (A16), (A17), and (95). We then take the linear combination

0 \equiv \frac{d}{d λ} (\frac{\partial L}{\partial θ^{i j}} - \frac{\partial L}{\partial θ^{i j †}}),

(A20)

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

x^{i j, λ} = \frac{d}{d λ} (- \frac{1}{2} K^{j i} g^{i j} - \frac{1}{2} K^{i j} g^{i j, T} - F^{i j *} D^{'' i j} + d^{'' i j, T} {\tilde{S}}^{*} τ^{i j †}) .

(A21)

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,

0 \equiv {([\frac{\partial L}{\partial U} - \frac{\partial L}{\partial U^{†}}]|}_{λ = 0},

(A22)

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,

x = {(\frac{\partial L_{- MOO}}{\partial U}|}_{λ = 0},

(A23)

x = - {(\frac{\partial L_{- MOO}}{\partial U^{†}}|}_{λ = 0},

(A24)

L_{- MOO} \equiv L - C_{MOO},

(A25)

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

0 \equiv \frac{d}{d λ} \frac{\partial L}{\partial U} = \frac{d}{d λ} \frac{\partial L_{- MOO}}{\partial U} + U^{λ, T} \frac{\partial L_{- MOO}}{\partial U} + x^{λ},

(A26)

0 \equiv \frac{d}{d λ} \frac{\partial L}{\partial U^{†}} = \frac{d}{d λ} \frac{\partial L_{- MOO}}{\partial U^{†}} + \frac{\partial L_{- MOO}}{\partial U^{†}} U^{λ *} + x^{λ}

(A27)

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,

0 \equiv \frac{d}{d λ} (\frac{\partial L}{\partial U} - \frac{\partial L}{\partial U^{†}}),

(A28)

Equivalent solutions for x^λ can be obtained from either Eq. (A26) or (A27); however, this is not necessary for the properties discussed in the present work.

REFERENCES

1.

P.

Pulay

,

Chem. Phys. Lett.

100

,

151

(

1983

).

https://doi.org/10.1016/0009-2614(83)80703-9

Google Scholar

Crossref

2.

S.

Sæbø

and

P.

Pulay

,

Chem. Phys. Lett.

113

,

13

(

1985

).

https://doi.org/10.1016/0009-2614(85)85003-X

Google Scholar

Crossref

3.

S.

Sæbø

and

P.

Pulay

,

Annu. Rev. Phys. Chem.

44

,

213

(

1993

).

https://doi.org/10.1146/annurev.pc.44.100193.001241

Google Scholar

Crossref

4.

C.

Hampel

and

H. J.

Werner

,

J. Chem. Phys.

104

,

6286

(

1996

).

https://doi.org/10.1063/1.471289

Google Scholar

Crossref

5.

M.

Schütz

,

G.

Hetzer

, and

H.-J.

Werner

,

J. Chem. Phys.

111

,

5691

(

1999

).

https://doi.org/10.1063/1.479957

Google Scholar

Crossref

6.

M.

Schütz

and

H.-J.

Werner

,

Chem. Phys. Lett.

318

,

370

(

2000

).

https://doi.org/10.1016/s0009-2614(00)00066-x

Google Scholar

Crossref

7.

M.

Schütz

and

H.-J.

Werner

,

J. Chem. Phys.

114

,

661

(

2001

).

https://doi.org/10.1063/1.1330207

Google Scholar

Crossref

8.

H.-J.

Werner

,

F. R.

Manby

, and

P. J.

Knowles

,

J. Chem. Phys.

118

,

8149

(

2003

).

https://doi.org/10.1063/1.1564816

Google Scholar

Crossref

9.

H.-J.

Werner

and

M.

Schütz

,

J. Chem. Phys.

135

,

144116

(

2011

).

https://doi.org/10.1063/1.3641642

Google Scholar

Crossref

PubMed

10.

M.

Schwilk

,

D.

Usvyat

, and

H.-J.

Werner

,

J. Chem. Phys.

142

,

121102

(

2015

).

https://doi.org/10.1063/1.4916316

Google Scholar

Crossref

PubMed

11.

H.-J.

Werner

,

J. Chem. Phys.

145

,

201101

(

2016

).

https://doi.org/10.1063/1.4968595

Google Scholar

Crossref

PubMed

12.

H.-J.

Werner

,

G.

Knizia

,

C.

Krause

,

M.

Schwilk

, and

M.

Dornbach

,

J. Chem. Theory Comput.

11

,

484

(

2015

).

https://doi.org/10.1021/ct500725e

Google Scholar

Crossref

PubMed

13.

M.

Schwilk

,

Q.

Ma

,

C.

Köppl

, and

H.-J.

Werner

,

J. Chem. Theory Comput.

13

,

3650

(

2017

).

https://doi.org/10.1021/acs.jctc.7b00554

Google Scholar

Crossref

PubMed

14.

Q.

Ma

and

H.-J.

Werner

,

J. Chem. Theory Comput.

14

,

198

(

2018

).

https://doi.org/10.1021/acs.jctc.7b01141

Google Scholar

Crossref

PubMed

15.

H.-J.

Werner

,

C.

Köppl

,

Q.

Ma

, and

M.

Schwilk

, in

Fragmentation

, edited by

M. S.

Gordon

(

John Wiley & Sons, Ltd.

,

Chichester, UK

,

2017

), pp.

1

–

79

.

Google Scholar

16.

F.

Neese

,

A.

Hansen

, and

D. G.

Liakos

,

J. Chem. Phys.

131

,

064103

(

2009

).

https://doi.org/10.1063/1.3173827

Google Scholar

Crossref

PubMed

17.

F.

Neese

,

F.

Wennmohs

, and

A.

Hansen

,

J. Chem. Phys.

130

,

114108

(

2009

).

https://doi.org/10.1063/1.3086717

Google Scholar

Crossref

PubMed

18.

C.

Riplinger

and

F.

Neese

,

J. Chem. Phys.

138

,

034106

(

2013

).

https://doi.org/10.1063/1.4773581

Google Scholar

Crossref

PubMed

19.

P.

Pinski

,

C.

Riplinger

,

E. F.

Valeev

, and

F.

Neese

,

J. Chem. Phys.

143

,

034108

(

2015

).

https://doi.org/10.1063/1.4926879

Google Scholar

Crossref

PubMed

20.

C.

Riplinger

,

P.

Pinski

,

U.

Becker

,

E. F.

Valeev

, and

F.

Neese

,

J. Chem. Phys.

144

,

024109

(

2016

).

https://doi.org/10.1063/1.4939030

Google Scholar

Crossref

PubMed

21.

M.

Saitow

,

U.

Becker

,

C.

Riplinger

,

E. F.

Valeev

, and

F.

Neese

,

J. Chem. Phys.

146

,

164105

(

2017

).

https://doi.org/10.1063/1.4981521

Google Scholar

Crossref

PubMed

22.

Y.

Guo

,

C.

Riplinger

,

U.

Becker

,

D. G.

Liakos

,

Y.

Minenkov

,

L.

Cavallo

, and

F.

Neese

,

J. Chem. Phys.

148

,

011101

(

2018

).

https://doi.org/10.1063/1.5011798

Google Scholar

Crossref

PubMed

23.

P. Y.

Ayala

and

G. E.

Scuseria

,

J. Chem. Phys.

110

,

3660

(

1999

).

https://doi.org/10.1063/1.478256

Google Scholar

Crossref

24.

G. E.

Scuseria

and

P. Y.

Ayala

,

J. Chem. Phys.

111

,

8330

(

1999

).

https://doi.org/10.1063/1.480174

Google Scholar

Crossref

25.

S. A.

Maurer

,

D. S.

Lambrecht

,

J.

Kussmann

, and

C.

Ochsenfeld

,

J. Chem. Phys.

138

,

014101

(

2013

).

https://doi.org/10.1063/1.4770502

Google Scholar

Crossref

PubMed

26.

S. A.

Maurer

,

L.

Clin

, and

C.

Ochsenfeld

,

J. Chem. Phys.

140

,

224112

(

2014

).

https://doi.org/10.1063/1.4881144

Google Scholar

Crossref

PubMed

27.

J.

Yang

,

G. K.-L.

Chan

,

F. R.

Manby

,

M.

Schütz

, and

H.-J.

Werner

,

J. Chem. Phys.

136

,

144105

(

2012

).

https://doi.org/10.1063/1.3696963

Google Scholar

Crossref

PubMed

28.

Y.

Kurashige

,

J.

Yang

,

G. K.-L.

Chan

, and

F. R.

Manby

,

J. Chem. Phys.

136

,

124106

(

2012

).

https://doi.org/10.1063/1.3696962

Google Scholar

Crossref

PubMed

29.

M.

Schütz

,

J.

Yang

,

G. K.-L.

Chan

,

F. R.

Manby

, and

H.-J.

Werner

,

J. Chem. Phys.

138

,

054109

(

2013

).

https://doi.org/10.1063/1.4789415

Google Scholar

Crossref

PubMed

30.

G.

Schmitz

,

B.

Helmich

, and

C.

Hättig

,

Mol. Phys.

111

,

2463

(

2013

).

https://doi.org/10.1080/00268976.2013.794314

Google Scholar

Crossref

31.

G.

Schmitz

,

C.

Hättig

, and

D. P.

Tew

,

Phys. Chem. Chem. Phys.

16

,

22167

(

2014

).

https://doi.org/10.1039/c4cp03502j

Google Scholar

Crossref

PubMed

32.

G.

Schmitz

and

C.

Hättig

,

J. Chem. Phys.

145

,

234107

(

2016

).

https://doi.org/10.1063/1.4972001

Google Scholar

Crossref

PubMed

33.

G.

Schmitz

and

C.

Hättig

,

J. Chem. Theory Comput.

13

,

2623

(

2017

).

https://doi.org/10.1021/acs.jctc.7b00180

Google Scholar

Crossref

PubMed

34.

Y.

Guo

,

W.

Li

, and

S.

Li

,

J. Phys. Chem. A

118

,

8996

(

2014

).

https://doi.org/10.1021/jp501976x

Google Scholar

Crossref

PubMed

35.

Y.

Guo

,

U.

Becker

, and

F.

Neese

,

J. Chem. Phys.

148

,

124117

(

2018

).

https://doi.org/10.1063/1.5021898

Google Scholar

Crossref

PubMed

36.

J. J.

Eriksen

,

P.

Baudin

,

P.

Ettenhuber

,

K.

Kristensen

,

T.

Kjærgaard

, and

P.

Jørgensen

,

J. Chem. Theory Comput.

11

,

2984

(

2015

).

https://doi.org/10.1021/acs.jctc.5b00086

Google Scholar

Crossref

PubMed

37.

P.

Baudin

,

P.

Ettenhuber

,

S.

Reine

,

K.

Kristensen

, and

T.

Kjærgaard

,

J. Chem. Phys.

144

,

054102

(

2016

).

https://doi.org/10.1063/1.4940732

Google Scholar

Crossref

PubMed

38.

T.

Kjærgaard

,

P.

Baudin

,

D.

Bykov

,

K.

Kristensen

, and

P.

Jørgensen

,

Wiley Interdiscip. Rev.: Comput. Mol. Sci.

7

,

e1319

(

2017

).

https://doi.org/10.1002/wcms.1319

Google Scholar

Crossref

39.

Z.

Rolik

,

L.

Szegedy

,

I.

Ladjánszki

,

B.

Ladóczki

, and

M.

Kállay

,

J. Chem. Phys.

139

,

094105

(

2013

).

https://doi.org/10.1063/1.4819401

Google Scholar

Crossref

PubMed

40.

P. R.

Nagy

,

G.

Samu

, and

M.

Kállay

,

J. Chem. Theory Comput.

12

,

4897

(

2016

).

https://doi.org/10.1021/acs.jctc.6b00732

Google Scholar

Crossref

PubMed

41.

P. R.

Nagy

,

G.

Samu

, and

M.

Kállay

,

J. Chem. Theory Comput.

14

,

4193

(

2018

).

https://doi.org/10.1021/acs.jctc.8b00442

Google Scholar

Crossref

PubMed

42.

K.

Raghavachari

,

G. W.

Trucks

,

J. A.

Pople

, and

M.

Head-Gordon

,

Chem. Phys. Lett.

157

,

479

(

1989

).

https://doi.org/10.1016/s0009-2614(89)87395-6

Google Scholar

Crossref

43.

D. G.

Liakos

and

F.

Neese

,

J. Chem. Theory Comput.

11

,

4054

(

2015

).

https://doi.org/10.1021/acs.jctc.5b00359

Google Scholar

Crossref

PubMed

44.

D. G.

Liakos

,

Y.

Guo

, and

F.

Neese

,

J. Phys. Chem. A

124

,

90

(

2020

).

https://doi.org/10.1021/acs.jpca.9b05734

Google Scholar

Crossref

PubMed

45.

A.

El Azhary

,

G.

Rauhut

,

P.

Pulay

, and

H.-J.

Werner

,

J. Chem. Phys.

108

,

5185

(

1998

).

https://doi.org/10.1063/1.475955

Google Scholar

Crossref

46.

G.

Rauhut

and

H.-J.

Werner

,

Phys. Chem. Chem. Phys.

3

,

4853

(

2001

).

https://doi.org/10.1039/b105126c

Google Scholar

Crossref

47.

M.

Schütz

,

H.-J.

Werner

,

R.

Lindh

, and

F. R.

Manby

,

J. Chem. Phys.

121

,

737

(

2004

).

https://doi.org/10.1063/1.1760747

Google Scholar

Crossref

PubMed

48.

K.

Ledermüller

,

D.

Kats

, and

M.

Schütz

,

J. Chem. Phys.

139

,

084111

(

2013

).

https://doi.org/10.1063/1.4818586

Google Scholar

Crossref

PubMed

49.

K.

Ledermüller

and

M.

Schütz

,

J. Chem. Phys.

140

,

164113

(

2014

).

https://doi.org/10.1063/1.4872169

Google Scholar

Crossref

PubMed

50.

M.

Dornbach

and

H.-J.

Werner

,

Mol. Phys.

117

,

1252

(

2019

).

https://doi.org/10.1080/00268976.2018.1537529

Google Scholar

Crossref

51.

G.

Rauhut

,

A.

El Azhary

,

F.

Eckert

,

U.

Schumann

, and

H.-J.

Werner

,

Spectrochim. Acta, Part A

55

,

647

(

1999

).

https://doi.org/10.1016/s1386-1425(98)00268-6

Google Scholar

Crossref

52.

G.

Rauhut

and

H.-J.

Werner

,

Phys. Chem. Chem. Phys.

5

,

2001

(

2003

).

https://doi.org/10.1039/b212590k

Google Scholar

Crossref

53.

T.

Hrenar

,

H.-J.

Werner

, and

G.

Rauhut

,

J. Chem. Phys.

126

,

134108

(

2007

).

https://doi.org/10.1063/1.2718951

Google Scholar

Crossref

PubMed

54.

J.

Gauss

and

H.-J.

Werner

,

Phys. Chem. Chem. Phys.

2

,

2083

(

2000

).

https://doi.org/10.1039/b000024h

Google Scholar

Crossref

55.

S.

Loibl

and

M.

Schütz

,

J. Chem. Phys.

137

,

084107

(

2012

).

https://doi.org/10.1063/1.4744102

Google Scholar

Crossref

PubMed

56.

S.

Loibl

and

M.

Schütz

,

J. Chem. Phys.

141

,

024108

(

2014

).

https://doi.org/10.1063/1.4884959

Google Scholar

Crossref

PubMed

57.

T.

Korona

,

K.

Pflüger

, and

H.-J.

Werner

,

Phys. Chem. Chem. Phys.

6

,

2059

(

2004

).

https://doi.org/10.1039/b316709g

Google Scholar

Crossref

58.

M.

Maurer

and

C.

Ochsenfeld

,

J. Chem. Phys.

138

,

174104

(

2013

).

https://doi.org/10.1063/1.4801084

Google Scholar

Crossref

PubMed

59.

D. H.

Friese

,

N. O. C.

Winter

,

P.

Balzerowski

,

R.

Schwan

, and

C.

Hättig

,

J. Chem. Phys.

136

,

174106

(

2012

).

https://doi.org/10.1063/1.4704788

Google Scholar

Crossref

PubMed

60.

M. S.

Frank

,

G.

Schmitz

, and

C.

Hättig

,

Mol. Phys.

115

,

343

(

2017

).

https://doi.org/10.1080/00268976.2016.1263762

Google Scholar

Crossref

61.

R.

Zhou

,

Q.

Liang

, and

J.

Yang

,

J. Chem. Theory Comput.

16

,

196

(

2020

).

https://doi.org/10.1021/acs.jctc.9b00806

Google Scholar

Crossref

PubMed

62.

N. J.

Russ

and

T. D.

Crawford

,

Phys. Chem. Chem. Phys.

10

,

3345

(

2008

).

https://doi.org/10.1039/b804119a

Google Scholar

Crossref

PubMed

63.

H. R.

McAlexander

,

T. J.

Mach

, and

T. D.

Crawford

,

Phys. Chem. Chem. Phys.

14

,

7830

(

2012

).

https://doi.org/10.1039/c2cp23797k

Google Scholar

Crossref

PubMed

64.

H. R.

McAlexander

and

T. D.

Crawford

,

J. Chem. Theory Comput.

12

,

209

(

2016

).

https://doi.org/10.1021/acs.jctc.5b00898

Google Scholar

Crossref

PubMed

65.

R.

D’Cunha

and

T. D.

Crawford

,

J. Chem. Theory Comput.

17

,

290

(

2021

).

https://doi.org/10.1021/acs.jctc.0c01086

Google Scholar

Crossref