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)] excitations42—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 methods45–50 based on projected atomic orbitals (PAOs),1,2 including second order Møller–Plesset perturbation theory (MP2), and also used them to semi-numerically calculate vibrational frequencies.51–53 Gauss and Werner also presented a pilot implementation for NMR shielding calculations with PAO-based local MP2 (LMP2).54 The first efficient such implementation, also employing the resolution of the identity (RI) approximation, was later reported by Loibl and Schütz55 and also adapted to magnetizabilities.56 Another study by Werner and co-workers examines the accuracy of PAO-based local correlation methods for polarizability calculations via finite differences,57 an approach that is difficult to apply to magnetic properties as it requires an implementation based on complex algebra. Maurer and Ochsenfeld developed an AO-based Laplace-transformed MP2 method for the computation of NMR shieldings, which can asymptotically achieve linear scaling or even sub-linear if only a few nuclei are of interest.58 Static and dynamic polarizabilities have also been implemented for MP2 and the approximate coupled cluster method CC2 using both RI and a Laplace-transformation.59 Frank et al. reported analytic gradients for a local MP2 variant, based on pair natural orbitals (PNOs), although they neglected pair natural orbitals (PNOs) relaxation.60 The complete analytic gradient for an orbital-specific virtual (OSV) local MP2 method was also published recently by Zhou et al.61 We should also mention here the work by Crawford and co-workers in the field of local coupled cluster linear response theory for the calculation of molecular properties, including frequency-dependent polarizabilities and specific optical rotations.62–65 Work by our group in Mülheim has focused on the domain-based local pair natural orbital (DLPNO) approximation and has resulted in orbital-unrelaxed first derivatives of DLPNO-CCSD, which can be used to calculate, e.g., dipole moments and HFCs.66,67 Fully orbital-relaxed first derivatives and nuclear gradients were also implemented for the DLPNO-MP2 method,68,69 and those results highlight the importance of the PNO relaxation contributions, e.g., for electric field gradients.68
In the current work, we present a derivation of analytic second derivatives of DLPNO-MP2, applicable to field-response properties, such as dipole polarizabilities and NMR shieldings. Section II focuses on the theory and working equations, while Sec. V provides a detailed benchmark study of the accuracy and efficiency of the method. To our knowledge, this is the first implementation of formally complete analytic second derivatives of a PNO-based local correlation method. On the one hand, MP2 is the simplest post-Hartree–Fock (HF) correlation method, and the present work is a stepping stone toward DLPNO-CCSD second derivatives, which we intend to pursue in the future. The pitfalls we encountered and their proposed solutions, discussed in Sec. III, are also likely relevant for other local correlation methods. On the other hand, MP2 itself and even more so the related double-hybrid DFT (DHDFT) have been shown to be very accurate for both polarizabilities and NMR shieldings.70–73 The present work extends their applicability to much larger systems and should therefore be viewed as part of an ongoing effort to reduce the computational cost of MP2 response property calculations. Works in this context are the integral-direct gauge-including atomic orbital (GIAO) MP2 implementation for shieldings of Kollwitz, Häser, and Gauss,74,75 the derivation of RI-MP2 second derivatives in combination with chain-of-spheres exchange (COSX),76 the Laplace-based approaches of Ochsenfeld, Hättig, and their co-workers,25,59 the proof-of-concept GIAO-LMP2 implementation,54 and, in particular, the efficient RI-based version of the latter.55,56 Note that the present derivation and implementation are not applicable to harmonic vibrational frequencies, as the efficient evaluation of the nuclear Hessian requires a substantially different algorithm.
II. THEORY
Both the polarizability tensor α and the shielding tensor for nucleus A, σA, can be obtained as second derivatives of the energy,77,78
where E, D, and h are the method-specific energy expression, density matrix, and one-electron part of the Fock matrix, respectively; x and y denote arbitrary Cartesian directions; and Fext, Bext, and MA 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 Fext. In the following discussion, λ will denote a generic multi-dimensional perturbation (such as Fext or Bext), 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, GIAOs77,79–82 are used to tackle the gauge-origin problem.83,84
Analytic derivatives of MP2 and variants thereof have been thoroughly discussed in the literature,85–94 so we will provide rather limited background here. We will focus on DLPNO-MP2 second derivatives, in particular, and remark on similarities and differences to the DLPNO-MP2 gradient.68,69 One point of note is that using a Lagrangian formulation together with the second derivative expressions above, which are “asymmetric” with respect to the two perturbations, is particularly advantageous for magnetic shielding calculations (though not necessarily for the polarizability tensor and other properties) as it only requires derivatives of the wave function parameters and the two-electron integrals with respect to the three magnetic field components.78 In contrast, the “symmetric” formulation, commonly used for vibrational frequency calculations,88–90 would also require derivatives with respect to the 3Natoms 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., or . However, the occupied orbital indices are often obvious from context or are added as upper indices of a matrix, e.g., or , and thus, the repetition is avoided to simplify the notation.
Be aware that despite the similar notation, the symbols , , , and have rather different meanings: the first denotes the two-electron repulsion integral , the second parameterizes excitations from orbitals i and j into PNOs and , while the last two are simply the Fock and overlap matrix elements in the PNO basis of pair ij.
When expressing derivatives in matrix form, we will use the following:
In the case of magnetic properties, special attention must be paid to complex conjugation: we use the superscripts “*,” “T,” and “†” to denote the complex conjugate, matrix transpose, and complex transpose of a quantity. A matrix trace is denoted as .
B. Overview of DLPNO-MP2
Here, we provide a brief review of the closed-shell (spin-restricted) DLPNO-MP2 method, which has been described in significant detail elsewhere.19 The following steps are taken to obtain the DLPNO-MP2 energy:
The self-consistent field (SCF) equations are solved, giving the HF energy, EHF, and canonical molecular orbitals (CMOs).
The occupied orbitals are localized according to the Foster–Boys (FB) criterion, as modified by Boys (referred to below as FB localization).95,96 Core and valence orbitals are localized separately.
PAOs are obtained as
and then normalized using the factor ,
For each occupied orbital, a correlation domain of PAOs is selected based on the differential overlap criterion,
with the default value of the parameter TCutDOPre = 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 TS (10−8 by default) and normalizing the rest. These are then transformed into pseudo-canonical NPAOs by diagonalizing the NPAO Fock matrix.
Strongly correlated electron pairs are determined using a three-step procedure. A pair ij is kept if DOIij is greater than TCutDOij = 10−5 and the semi-canonical pair energy estimate [Eq. (8)] using the collinear dipole approximation to the exchange integrals97 [Eq. (8) computed with from Eq. (10) instead of Eq. (9)] is greater than TCutPre = 10−6. The energy contribution from the screened-out pairs, ΔEPre, is estimated using the more accurate dipole approximation to the exchange integrals11 [Eq. (9)],
where and are NPAO orbital energies in the domain of orbitals i and j, respectively, while Fii and Fjj 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 TCutDO = 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 in the basis of NPAOs for the given orbital pair,
where and are NPAO orbital energies in the ij pair domain. PNOs with an occupation number greater than TCutPNO = 10−8 are kept, while the rest form the “complementary” PNO (CPNO) basis referred to below. A tighter threshold, TCutPNO × TScalePNOCore = 10−10, is used for pairs that include core orbitals.98 The PAO to PNO transformation matrices are stored on a disk. A correction to the energy, ΔEPNO, is computed as the difference between the semi-canonical correlation energy, calculated in the NPAO and PNO basis sets.
The full amplitude equations (see Sec. II C for definitions) are solved iteratively in the PNO basis, and the final DLPNO-MP2 energy is calculated,
C. DLPNO-MP2 Lagrangian
Because the DLPNO-MP2 energy expression is non-variational, a Lagrangian formulation is used for analytic derivatives,99 which is slightly different from the one used previously,68,69
The first four terms are the energy contributions, the next five are the Brillouin, localization, core–valence (CV), semicanonical, and PNO Lagrangian constraints, and the last three are orthonormality conditions for the MOs, PNOs, and NPAOs. The various terms will be explained as follows. EHF is the HF energy expression,
with two-electron integrals in the Mulliken (1*1|2*2) notation. is the PNO-basis Hylleraas functional in the spin-adapted form100
ΔEPNO is an estimate for the error due to the PNO truncation, calculated as the difference between the correlation energy computed from the semicanonical amplitudes in the full NPAO basis and in the truncated PNO basis,
The energy contribution of the screened-out pairs takes the form of an approximate Hylleraas functional,
where 1S(ij) = 1 if pair ij was screened-out, 0 otherwise. The dipole approximation to the exchange integrals is given in Eq. (9). The Brillouin condition is formulated as
where zai 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., . 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 as101
where c denotes the LMO coefficients. In contrast to the exponential parameterization used in Refs. 68 and 69, U is not constrained to be unitary but is instead constrained by the MO orthonormality condition,
which is fulfilled by the CMO coefficients and the unperturbed LMO coefficients c(0). Here, δpq is the Kronecker delta, and xpq are the associated Lagrange multipliers. The term CLoc enforces the FB localization condition separately on the core and valence orbitals and introduces the required Lagrange multipliers lij,
The final MO constraint is the core–valence separation condition with an associated Lagrange parameter ,
Recall that the core and valence orbitals are localized separately because different PNO truncation thresholds are used for pairs that include core orbitals.98 Therefore, the core–valence separation condition is necessary even if all electrons are included in the correlation treatment. On the other hand, the localization condition for core orbitals is not needed if the latter are not localized, which is the case in the frozen-core approximation.
We denote the PNOs by , the CPNOs by , and the union of the two spaces by . The transformation matrices from redundant PAOs to each of these sets (for orbital pair ij) are denoted by , , and , respectively. The PNOs and CPNOs combined span the same space as the NPAOs. In the presence of a perturbation, we allow the PNOs to relax within the entire PNO + CPNO space. Rotations within each domain are thus parameterized as
Analogous to U, the matrices θij are not necessarily unitary but rather are constrained by the PNO (and CPNO) orthonormality condition,
Note that Sij is expanded as
A further PNO constraint ensures that the semi-canonical pair densities remain block-diagonal,
where are unknown Lagrange multipliers. An orthonormality condition also applies to the orbital domain NPAOs used for the prescreening correction, parameterized as
where are the transformation matrices from redundant to non-redundant orbital domain PAOs, Si is the overlap matrix in the latter basis, and are the necessary Lagrange multipliers. Finally, the semi-canonical amplitudes are parameterized as
and the semi-canonical residual conditions CSC are included in the Lagrangian,
where are Lagrange multipliers and the semi-canonical residuals are defined as
The optimization parameters in the present DLPNO-MP2 Lagrangian are thus Upq, , , , , , zai, , lij, xpq, , , , , and their complex conjugates.
D. Equations for first derivatives
The derivative of the DLPNO-MP2 energy with respect to the first external perturbation κ (in this work either an electric field or nuclear magnetic moment, on which only the one-electron part of the Hamiltonian depends) is equal to the respective partial derivative of the Lagrangian,
if and only if the Lagrangian is made stationary with respect to all parameters. The equations for the necessary stationarity conditions will be presented in this section. Despite the slightly different formulation of the Lagrangian, the final working equations are completely equivalent to those derived in Ref. 69. Note that Eq. (43) is only valid if the basis functions do not depend on the perturbation κ. Therein, PSCF is the SCF density matrix,
while D is the orbital-relaxed DLPNO-MP2 difference density matrix,
where, as stated above, only the occupied–virtual and virtual–occupied blocks of z are non-zero with . Analogously, only has non-zero core–valence blocks with . D′ is the “PNO-relaxed” difference density matrix,
The prescreening contributions are defined as
The values of the required variables are obtained by solving the stationarity conditions equations (and their complex conjugates),
The first eight are already fulfilled at the end of the DLPNO-MP2 energy calculation. The PNO amplitudes are obtained by iteratively solving the amplitude equations,
Note that the equations obtained from the derivative with respect to the contravariant amplitudes, , are equivalent to a linear combination of and 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 FCut = 10−5.
The semi-canonical amplitude stationarity conditions take the form
The matrix is introduced to “select” the truncated PNO basis. Equation (54) can be transformed to the pseudo-canonical NPAO basis using the unitary matrix and solved for wij at λ = 0 (i.e., θij = I),
with
Equations for vij are obtained from the Lagrangian stationarity conditions with respect to ,
These can be solved in the PNO + CPNO basis, which diagonalizes . A detailed derivation, as well as expressions for xij, is given in the Appendix. The remaining Lagrange multipliers 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)102 and Z-vector coupled-perturbed localization (Z-CPL)45 equations for z and l, respectively, while rotations between core and valence orbitals result in the Z-vector core–valence (Z-CV) equations for . For rotations within a block of localized orbitals (either core or valence), we obtain the Z-CPL equations
where i, j, k are all within the same (core or valence) block. The localization response operator is defined as follows (note the change in sign and indexing compared to Ref. 69):
with
is the two-body density matrix
contains contributions from the approximate exchange integrals in the pre-screening correction,
where
The Z-CPL equation (59) must be solved iteratively. The stationarity condition with respect to rotations between core and valence MOs gives the Z-CV equations,
These are solved in the CMO basis. Note that in frozen-core calculations, the terms involving , , and drop out.
Finally, the Z-CPSCF equations are obtained from the occupied–virtual block of Eq. (A22),
with the two-electron operator defined as
and likewise for matrices in the AO basis. The right-hand side of Eq. (69) contains several terms,
Note that the sum in Eq. (72) assumes a one-to-one correspondence between PAOs and their parent AOs and that . The intermediates Yα, Yβ, and Yγ are defined as follows:
As discussed in the Appendix, xij 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 xi, which are obtained analogously by solving . We note that the first term in Eq. (79) is missing from Eq. (56) of Ref. 69, which appears to be an error in the publication as the term was, in fact, implemented in the DLPNO-MP2 gradient code. This is the final major step of the DLPNO-MP2 density calculation, after which D can be completed.
E. Equations for second derivatives
The second derivative of the DLPNO-MP2 energy with respect to perturbations κ and λ, the latter being an electric or magnetic field in this work, is equal to the derivative of Eq. (43),
Apart from the second derivative integrals hκ,λ (which are zero for the polarizability, i.e., when κ = λ = Fext), 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 PSCF,λ requires the MO response coefficients Uλ,
where in the second line, we have used the perturbed MO orthonormality constraint
We note in passing that the AOs depend on the external magnetic field when GIAOs are used, but not on the external electric field; thus, . are obtained from the perturbed Brillouin condition, which in conjunction with Eq. (82) gives the CPSCF equations,101
where F(λ) excludes derivatives of the MO coefficients. In the electric (magnetic) field case, the right-hand side is real (imaginary) and therefore so are . 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 [in conjunction with Eq. (82)],
which can be solved in the CMO basis. The valence–valence block of Uλ is obtained from the perturbed localization condition , which gives the CPL equations,45
where the plus/minus signs are used in the magnetic/electric field case, respectively, and
s(λ) excludes derivatives of the MO coefficients and is only non-zero in the GIAO case, where it contains the perturbed dipole integrals,
which can be calculated as linear combinations of quadrupole integrals. If all electrons are included in the MP2 treatment, an analogous set of CPL equations must be solved for the core–core block, which can be obtained from Eq. (85) by substituting core for valence indices and vice versa. In the frozen-core approximation, the core–core block of Uλ is only constrained by Eq. (82), which can be solved as
The same solution is used for the virtual–virtual block of Uλ . Note that the left-hand side of CPL equations in the electric field case is equivalent to that of the Z-CPL equations (59), which can be written as
and if sij = 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λ, , Tij,λ, , wij,λ, θij,λ, tij,λ, and θi,λ. The PAO coefficients also depend on the perturbation through U,
For magnetic perturbations, . The perturbed (C)PNO orthonormality condition gives a relationship analogous to Eq. (82),
with
Note that the perturbed PAO overlap includes the PAO response [Eq. (92)] and is non-zero even for electric perturbations, unlike the perturbed AO overlap. We use Eq. (95) to remove the dependence on θij,λ in the perturbed SC residual equations,
where
Equation (96) is once again solved in the pseudo-canonical NPAO basis. Next, the perturbed PNO constraint gives the PNO–CPNO block of θij,λ,
These equations are solved in the basis that diagonalizes . 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 Tij,λ are then obtained from the perturbed PNO residual equations,
which can be solved iteratively with the same solver as Eq. (53) The prescreening contributions to the response density, Do(Pre),λ and Dv(Pre),λ, require the NPAO response coefficients θi,λ, which are only constrained by the perturbed NPAO orthonormality condition with the possible solution,
Also required are the perturbed prescreening amplitudes, obtained from the equations giving
For magnetic perturbations, . The perturbed Lagrange multipliers wij,λ are obtained from the equations
which are solved in the pseudo-canonical NPAO basis with
The equations for the perturbed multipliers vij,λ are derived in the Appendix and have the form
with
The perturbed localization condition multiplier lλ is obtained from the perturbed Z-CPL equations,
where and are straightforwardly derived from Eqs. (60) and (61). The derivatives of the prescreening contributions and are somewhat more involved, so we provide them here,
The perturbed Z-CV equations have the form
where is the derivative of Eq. (68). Finally, the perturbed Z-CPSCF equations must be solved for zλ,
Thus, the full response density can be assembled and contracted with the property integrals to obtain the second derivative.
At the end of this section, we briefly summarize the order of computational steps needed to obtain the DLPNO-MP2 response density in Table I. The actual procedure used in our implementation is given in more detail in Sec. IV.
In . | . | Equations . | . | Out . |
---|---|---|---|---|
SCF, FB-loc. | → | c, EHF | ||
c | → | LMP2 guess (13)–(15) | → | |
dij | → | LMP2 iter. (53) | → | |
→ | PNO cond. (58) | → | vij | |
→ | SC cond. (56) | → | wij | |
→ | Equation (63) | → | ΓK | |
→ | Equation (46) | → | D′ | |
D′, ΓK | → | Z-CPL (59) | → | l |
l, ΓK | → | Z-CV (67) | → | |
→ | Z-CPSCF (69) | → | z | |
→ | Equation (45) | → | D | |
c | → | CPSCF (83), CV (84), CPL (85), (82) | → | Uλ, cλ |
→ | PNO, respectively, (95), (96), and (99) | → | ||
dij,λ, Tij | → | Equation (100) | → | Tij,λ |
→ | Pert. PNO cond. (107) | → | vij,λ | |
→ | Pert. SC cond. (105) | → | wij,λ | |
→ | Deriv. of Eq. (63) | → | ΓK,λ | |
→ | Deriv. of Eq. (46) | → | D′λ | |
D′λ, ΓK,λ, l | → | Pert. Z-CPL (109) | → | lλ |
→ | Pert. Z-CV (112) | → | ||
→ | Pert. Z-CPSCF (113) | → | zλ | |
→ | Deriv. of Eq. (45) | → | Dλ | |
cλ, Dλ | → | Property equations (1), (2), and (80) | → | α, σA |
In . | . | Equations . | . | Out . |
---|---|---|---|---|
SCF, FB-loc. | → | c, EHF | ||
c | → | LMP2 guess (13)–(15) | → | |
dij | → | LMP2 iter. (53) | → | |
→ | PNO cond. (58) | → | vij | |
→ | SC cond. (56) | → | wij | |
→ | Equation (63) | → | ΓK | |
→ | Equation (46) | → | D′ | |
D′, ΓK | → | Z-CPL (59) | → | l |
l, ΓK | → | Z-CV (67) | → | |
→ | Z-CPSCF (69) | → | z | |
→ | Equation (45) | → | D | |
c | → | CPSCF (83), CV (84), CPL (85), (82) | → | Uλ, cλ |
→ | PNO, respectively, (95), (96), and (99) | → | ||
dij,λ, Tij | → | Equation (100) | → | Tij,λ |
→ | Pert. PNO cond. (107) | → | vij,λ | |
→ | Pert. SC cond. (105) | → | wij,λ | |
→ | Deriv. of Eq. (63) | → | ΓK,λ | |
→ | Deriv. of Eq. (46) | → | D′λ | |
D′λ, ΓK,λ, l | → | Pert. Z-CPL (109) | → | lλ |
→ | Pert. Z-CV (112) | → | ||
→ | Pert. Z-CPSCF (113) | → | 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 Aloc,± 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 Aloc,± and thus prevent convergence of the CPL/Z-CPL equations or negatively impact the results, as shown for the (CH3)3P(CO)4Os⋯Cr(CO)5 molecule (reference code KAMDOR)106 in Ref. 69. To systematically treat this problem, here we take the same approach as for the DLPNO-MP2 gradient, namely, to modify the Lagrangian so as to obtain a nonsingular problem with the nullspace projected out. Note that the matrix Aloc,−, 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 Aloc,+, which occurs in the CPL and perturbed Z-CPL equations for magnetic perturbations. Thus, the two matrices have different eigenvalues and eigenvectors. In order to ensure that both types of equations are properly treated, it is most convenient to reformulate the localization constraint in the Lagrangian using separate real and imaginary terms,
where the equivalence of the two formulations is demonstrated for the valence orbitals. The core orbital terms are analogous, so for simplicity, we will omit them in this section. Note that both of the newly introduced Lagrange multiplier matrices and are real by definition. Using this definition of CLoc, Eq. (91) becomes
where
The real and imaginary parts of Eq. (116) can be separated, and the latter yields because Lij is real. Note that as the unperturbed quantities are purely real, coincides with Aloc,− and coincides with Aloc,+. We then define the eigendecompositions of these matrices as
where the index n″ denotes the eigenvectors of either matrix (determined from context), n′ denotes those with (near-)zero eigenvalues, and n denotes the rest. We change the Lagrangian constraint to
where δij,kl = δikδjl and the projectors onto the nullspace of are
Thus, the Z-CPL equations become
where
Once again, is zero, while the contributions of to Eqs. (68), (73), (112), and (113) must be calculated with the projected matrix,
The electric/magnetic CPL equations become
i.e., the nullspace must first be projected out of the right-hand side in Eq. (127), the equations are solved, and afterward, the nullspace is also projected out of the solution before the perturbed orthonormality constraint is used to obtain the lower triangle of the matrix,
The perturbed Z-CPL equations become
with
Equation (130) is separated into real and imaginary parts, of which only the former/latter yields a nonzero solution for electric/magnetic perturbations, respectively. Note that for magnetic perturbations, and . The contributions of to Eqs. (112) and (113) must be calculated with the projected matrices
Finally, we demonstrate the numerical instabilities that can occur if the localization response matrix, or , is singular, as well as the effectiveness of the procedure to remove the near-zero eigenvalues. We chose the octahedral molecules SF6 and SeF6 with bond lengths 1.560 and 1.688 Å, respectively. SF6 has six eigenvalues of smaller than 10−4 due to a continuous degeneracy in the FB LMOs among the 18 fluorine sp3-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, SeF6 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 only has ten singular eigenvalues for SeF6, due to the 3d orbitals, and none for SF6. Table II shows that if no measures are taken to remove the singularities, the iterative solution of the perturbed CPL, Z-CPL, or PNO amplitude equations may not converge or large errors may occur in the calculated properties, which is arguably worse. However, if the singular eigenvectors are projected out, the iterative solutions are stable and the results are in line with the expected accuracy for DLPNO-MP2.
System . | SF6 . | SeF6 . | ||
---|---|---|---|---|
. | FC . | AE . | FC . | AE . |
nullspace | 6 | 9 | 16 | 24 |
nullspace | 0 | 0 | 10 | 10 |
Errors without projection: | ||||
|Δα| (Å3) | 8 × 105 | …a | 17 | …a |
max (|ΔσF|) (ppm) | 0.26 | 0.28 | 2.1 | 29 |
Errors with projection: | ||||
|Δα| (Å3) | <0.001 | <0.001 | <0.001 | <0.001 |
max (|ΔσF|) (ppm) | 0.26 | 0.28 | 0.47 | 0.55 |
System . | SF6 . | SeF6 . | ||
---|---|---|---|---|
. | FC . | AE . | FC . | AE . |
nullspace | 6 | 9 | 16 | 24 |
nullspace | 0 | 0 | 10 | 10 |
Errors without projection: | ||||
|Δα| (Å3) | 8 × 105 | …a | 17 | …a |
max (|ΔσF|) (ppm) | 0.26 | 0.28 | 2.1 | 29 |
Errors with projection: | ||||
|Δα| (Å3) | <0.001 | <0.001 | <0.001 | <0.001 |
max (|ΔσF|) (ppm) | 0.26 | 0.28 | 0.47 | 0.55 |
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, NAO, but it only represents the Nvirt virtual orbitals. Therefore, there are exactly Nocc 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 Nocc redundant eigenvectors, all smaller domains have virtually continuous eigenvalue spectra, i.e., the matrix is ill-conditioned. By default, eigenvectors with eigenvalues greater than TS = 10−8 form the NPAO domain. It is apparent that the lowest remaining eigenvalue sn will be very close to the threshold, and thus, the normalization coefficient of that eigenvector will be . This is not a cause for concern, and indeed, the energy and gradient calculations do not suffer. It is, however, a problem for second derivatives due to the perturbed PNO/CPNO coefficients d″ij,λ. To illustrate this, we express the latter as
where US is the unitary matrix that diagonalizes the PAO domain overlap matrix, s is the diagonal matrix of eigenvalues, and USD is a unitary transformation from the overlap eigenbasis to the SC pair density eigenbasis. We will drop the ij indices in the rest of this section for simplicity of notation. Let us also assume for now that no PNO truncation is used; thus, using Eq. (95), we can express
and therefore,
The largest elements of the matrices are estimated based on . 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 C2v symmetry at a distance of 10.5 bohrs between the two molecular planes. No truncations were used other than the domain truncation with TCutDO = 10−2. Ignoring the line “Analytic v2” for now, the large error in the analytic polarizability at TS = 10−8 is apparent. This disappears at higher values of the cutoff, and at TS = 10−5, there is very good agreement with both the semi-numeric polarizability (calculated from analytic dipole moment calculations with an explicit electric field of 10−4 a.u.) and the RI-MP2 value. Note that even the errors in the relaxed MP2 density get larger with lower values of TS, hinting at the underlying numerical instability, although, of course, these errors are negligible.
. | TS = 10−8 . | TS = 10−7 . | TS = 10−6 . | TS = 10−5 . | RI-MP2 . |
---|---|---|---|---|---|
Numeric | 19.1393 | 19.0770 | 19.1034 | 19.1146 | 19.1190 |
Analytic v1 | 17.5287 | 19.0675 | 19.0973 | 19.1149 | 19.1190 |
Analytic v2 | 19.1033 | 19.0688 | 19.0974 | 19.1149 | |
max(|ΔDμν|) | 0.0117 | 0.0110 | 0.0055 | 0.0010 |
. | TS = 10−8 . | TS = 10−7 . | TS = 10−6 . | TS = 10−5 . | RI-MP2 . |
---|---|---|---|---|---|
Numeric | 19.1393 | 19.0770 | 19.1034 | 19.1146 | 19.1190 |
Analytic v1 | 17.5287 | 19.0675 | 19.0973 | 19.1149 | 19.1190 |
Analytic v2 | 19.1033 | 19.0688 | 19.0974 | 19.1149 | |
max(|ΔDμν|) | 0.0117 | 0.0110 | 0.0055 | 0.0010 |
After some experimentation, we found an alternative solution for Eq. (95), which partially alleviates the problem and results in the line “Analytic v2” in Table III. First, we express θλ as
where both θ1,λ and θ2,λ are defined as skew-Hermitian matrices and the latter only has elements in the PNO–CPNO block. Thus, the first and last terms in Eq. (136) are used to ensure that the perturbed orthonormality and SC conditions [Eqs. (95) and (99)] are fulfilled, respectively, while the second term is arbitrary. We got the best results by choosing it so as to minimize the norm of the perturbed PNO/CPNO coefficients, excluding the θ2,λ part,
which is solved in the overlap eigenbasis as
All of the above also applies to the orbital-specific PAO domains used for the prescreening contributions and their response parameterized by θi,λ. Of course, in that case, θi,2,λ = 0.
This approach basically serves to better distribute the large floating point numbers throughout d″λ and seems to reduce the error in the calculated properties, as evidenced by the toy example in Table III. However, agreement with the numeric polarizability is still not perfect, and in fact, large, unpredictable errors still occur in calculations for realistic systems when the default value of TS 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 TS = 10−5, with which the issue no longer occurs. The reduced NPAO domains do lead to a slightly increased error in the correlation energy, but the difference is negligible compared to the deviations due to the domain and PNO approximations, which are themselves sufficiently small by design, so as to obtain results with chemical accuracy.
We note in passing that the entire problem disappears if a linearly independent set of virtual orbitals is used, such as FB-localized MOs. In fact, we implemented this option and confirmed that the instability does not occur. However, on the one hand, localizing virtual orbitals to a true maximum, rather than a saddle point, of the FB condition is very difficult and time-consuming, even with the augmented Hessian solver used in our implementation. On the other hand, much larger domains were necessary to obtain accurate results, corroborating recent observations.107,108 Therefore, we will not discuss this approach in the rest of this paper.
C. PNO response singularities
To solve Eqs. (58), (99), and (107), which follow from the PNO constraint , it is necessary to divide by the difference 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 is nearly continuous, and thus, both the smallest PNO and the largest CPNO occupation numbers are of the same order of magnitude as TCutPNO, 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 was smaller than 10−2, which can be taken as a necessary, but not sufficient, condition for problems to occur. Clearly, it is necessary to treat this problem systematically, since it is not possible to know when the results will be, or indeed have been, affected, without reference calculations, e.g., at other values of TCutPNO.
The approach we propose is to introduce a level shift ε in the denominator of Eqs. (58), (99), and (107), thus ensuring that . To that end, we modify the PNO condition as follows:
In effect, we shift all kept PNO occupation numbers up by ε. This modification has no effect on the DLPNO-MP2 energy, as the occupation numbers are only used for PNO truncation. It does, however, lead to extra terms in the right-hand sides of Eqs. (99) and (107) (given below in the pair density-diagonalizing basis),
An appropriate value of the level shift parameter ε must be chosen: if it is too low, it will have no effect. However, if it is too large, it can introduce additional round-off errors. It makes sense to choose a value proportional to TCutPNO, i.e., ε = εscale × TCutPNO. After some experimentation (see below), we found that εscale = 0.1 provides stable results without introducing additional errors. Note that because the value of TCutPNO 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, TS = 10−5 and TCutPNO = 10−7, as well as with TCutPNO = 5.96 × 10−7, the maximal absolute errors with respect to RI-MP2 are huge (40.5 and 106.3 ppm, respectively). The level shift scheme with εscale in the range 0.05–0.1 is successful in resolving the problem and producing a smooth curve. A value of 0.01 is not sufficient to remove the singularity at TCutPNO = 10−7, while εscale = 1.0 results in a smooth curve but with an additional small deviation for all values of TCutPNO. Thus, we consider εscale = 0.1 to be an appropriate value.
Initially, no adjustment or level shift was applied to the calculations in Secs. V B and V C. Consequently, we encountered several calculations with conspicuously large errors, which are discussed whenever they occur in Secs. V B and V C. In all of these cases, the level shift approach (with εscale = 0.1) was successful in reducing the errors to within the expected bounds.
IV. IMPLEMENTATION
We have implemented the DLPNO-MP2 response density calculation in the ORCA package,109 closely following the algorithm for the DLPNO-MP2 gradient.69 The code is structured as follows:
First, the SCF equations are solved.
Then, the DLPNO-MP2 energy and relaxed density calculation is performed. The PNO coefficients dij and amplitudes Tij 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 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 .
The three-index integrals and are calculated and stored.
In a loop over orbital pairs, the CPNO coefficients d′ij are reconstructed, and and are calculated according to Eqs. (14) and (96). θij,λ are computed, and the perturbed PNO coefficients dij,λ 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 Dij, Dij, and their derivatives are calculated.
The PNO-dependent terms are evaluated in a loop over orbital pairs:
- 9.1.
The full d″ij, , d″ij,λ, and are reconstructed.
- 9.2.
The intermediates τij and τij,λ are calculated.
- 9.3.
Equations (58) and (107) are solved to obtain vij and vij,λ.
- 9.4.
Equations (56) and (105) are solved to obtain wij and wij,λ. Their contributions to D′o, D′ij, and their derivatives are evaluated.
- 9.5.
Gi(ij), Gj(ij), Gi(ij),λ, and Gj(ij),λ and their contributions to ΓK, ξij, and Yγ, and their derivatives are computed.
- 9.6.
Terms related to ξij and τij and their derivatives are added to Yδ and Yδ,λ.
- 9.1.
ΓK and ΓK,λ are read from disk, and their contributions to , , Yα,λ, Yβ, and Yβ,λ are evaluated in the AO basis by generating the required three-index integrals on the fly.
The prescreening contributions to , , , and D′ are calculated.
is completed, and the Z-CV equations (112) are solved to obtain .
is completed, and the Z-CPL equations (109) are solved to obtain . In our implementation, a conjugate gradient solver is used.
is completed, and the Z-CPSCF equations (69) are solved to obtain . 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 , dij,λ, and Tij,λ—are kept for all perturbations simultaneously. For these reasons, the implemented algorithm is not suitable for perturbations whose number grows with the system size, e.g., nuclear Hessians. It is also worth pointing out that the NPAO domains and the SC amplitudes are calculated four times in total—in steps 2 (twice), 6, and 9.1—while their response is calculated twice—in steps 6, and 9.1. This is not an insignificant computational overhead, as shown in Sec. V D, but is done intentionally to avoid the storage and I/O of quantities in the full NPAO domains, which are much larger than the PNO domains.
The asymptotic scaling behavior of the DLPNO-MP2 gradient algorithm was discussed in Ref. 69, and that discussion applies here as well. Briefly, the number of orbital pairs that survive the dipole-based prescreening scales linearly with the system size in the asymptotic limit, while the PAO correlation domains reach a finite size for large systems. This allows for asymptotic linear scaling of most major steps in the algorithm, provided sparsity is properly exploited using the sparse map infrastructure introduced in Ref. 19. In particular, the RI integral transformation, PNO generation, and iterative solution of the amplitude equations are asymptotically linear scaling. Two-electron contributions coming from the Fock matrix, e.g., in the CPSCF and Z-CPSCF equations, can be calculated with effective quadratic scaling using the RIJCOSX (RI for Coulomb and chain-of-spheres exchange) approximation.110 The calculation of terms involving (K|ip) integrals in the Z-CPSCF equations’ [Eq. (74)] right-hand sides asymptotically scales as . Together with the 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 , but these steps constitute a very small part of the overall computation time.
V. RESULTS AND DISCUSSION
A. Computational details
To examine how quickly DLPNO-MP2 response properties converge toward the RI-MP2 results as the various thresholds are tightened, we performed a series of polarizability and NMR shielding calculations on the molecules depicted in Fig. 3.
For the shielding calculations, the pcSseg-2113 basis set was used, along with the cc-pwCVTZ/C114–116 auxiliary basis set, which was shown to be an appropriate combination in our previous study.72 RI was used only for Coulomb integrals in all Fock matrix builds with no approximation for the exchange integrals, while RIJK was used for the perturbed two-electron integrals over GIAOs. In both cases, the def2/JK117 auxiliary basis set was employed. All orbitals were included in the MP2 treatment.
For the polarizability calculations, we used the def2-TZVPD118,119 basis set together with the aug-cc-pVTZ/C114 auxiliary basis set. Once again, def2/JK117 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 TCutPNO for the valence orbitals with very conservative values of TCutPNO = 10−12 for the core orbitals. We set FCut = 10−8, TS = 10−8, and all other DLPNO thresholds to 0.
Varying TCutPNO for the core orbitals with TCutPNO = 10−10 for the valence orbitals. We set FCut = 10−8, TS = 10−8, and all other DLPNO thresholds equal to 0.
Varying FCut with TCutPNO = 10−12 for the core orbitals and TCutPNO = 10−10 for the valence orbitals. We set TS = 10−8 and all other DLPNO thresholds equal to 0.
Same as above but also fixing FCut = 10−5 and varying TCutDO. Two values were used for TS to examine the numerical stability: 10−8 and 10−5.
Same as above but also fixing TCutDO = 10−3, TCutDOPre = 0.03, and TS = 10−5 and varying TCutPre. We set TCutDOij = 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: TCutPNO = 10−8 (10−10 for core orbitals), FCut = 10−5, TCutDO = 10−2, TCutDOPre = 0.03, TCutDOij = 10−5, TCutPre = 10−6, TCutMKN = 10−3, TCutC = 10−3, but with TS = 10−5 and a PNO level shift εscale = 0.1.
Same as above with the “LoosePNO” default thresholds: TCutPNO = 10−7 (10−9 for core orbitals) and TCutDO = 2 × 10−2.
We also performed shielding calculations with RI-MP2, NormalPNO, LoosePNO, and TightPNO (TCutPNO = 10−9 for valence and 10−11 for core orbitals, TCutDO = 5 × 10−3), together with the RIJCOSX approximation to see how the COSX and DLPNO errors stack up. The keyword DefGrid3 was used to choose finer COSX integration grids than the default. Note that the grids in ORCA were updated for version 5 of the program, and the DefGrid3 settings are similar to the “RIJCOSX-XL” settings used in Ref. 72.
To investigate the efficiency of the implementation for larger system sizes, we used idealized glycine chains in either a linear (ϕ = 180°, ψ = 180°) or α-helix (ϕ = −57°, ψ = −47°) conformation. The neutral form of the molecules was used, rather than the zwitterion. The structures were generated using Gabedit.120 As above, these calculations were performed with NormalPNO (and RI-MP2) and RIJCOSX. These calculations ran in a node-exclusive fashion on 8 Intel Xeon E5-2690v2@3.00GHz central processing unit (CPU) cores with 6 GB of random-access memory (RAM) per core.
To showcase all the capabilities of the implementation, another large shielding calculation was performed on the vancomycin molecule using the DSD-PBEP86 DHDFT functional,121 NormalPNO thresholds, RIJCOSX, frozen-core orbitals, and CPCM for implicit solvation in dimethylsulfoxide (DMSO). This calculation ran on 12 Intel Xeon E5-2687Wv4@3.00GHz CPU cores with 15 GB RAM per core.
Finally, a series of DLPNO-MP2 polarizability calculations were performed on [n]helicene (n = 7, 9, 11, 13, 15) molecules using the aug-cc-pVTZ,122,123 aug-cc-pVTZ/C, and def2/JK basis sets, NormalPNO thresholds, frozen-core orbitals, and RIJCOSX. The geometries were optimized with RI-MP2 (for n = 7–11) and DLPNO-MP2 (for n = 13–15) using the cc-pVTZ123,124 and cc-pVTZ/C114 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 TCutPNO are shown in Fig. 4(a) for the different nuclides in the systems ATP4−, caffeine, ebselen, and penicillin. We note that the errors diminish rapidly with decreasing TCutPNO, and at the default value of 10−8 for valence orbitals, they are much lower than the inherent error in MP2. For example, the MAE for C is 0.1 ppm, which is an order of magnitude smaller than the average deviation from CCSD(T).71,72 It is also noteworthy that the results are insensitive to core PNO truncation, except for the heavier elements, particularly P and Se. This may be, in part, because the pcSseg-2 basis set does not fully capture core correlation effects.72
For simplicity in the discussion, we want a single measure of the DLPNO error; however, the shielding constants have very different magnitudes for different nuclides. Thus, we define an absolute weighted error (AWE), whereby we scale the error in the calculated isotropic shielding (e.g., vs RI-MP2) by a factor wA, derived from the experimental range of chemical shifts for the nuclide A,
We use the values for and shown in Table IV. A reasonable target accuracy is a mean AWE (MAWE) below 0.5%, which corresponds to 0.5 ppm for 13C, 0.03 ppm for 1H, and about 1 ppm for 15N and 17O. 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 ATP4− 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 TCutPNO = 10−9 with MaxAE = 1.51 ppm for 15N vs 0.59 and 0.04 ppm for the same nucleus at TCutPNO = 10−8 and TCutPNO = 10−10, respectively. With the level shift-based adjustment (see Sec. III C), this error drops to 0.17 ppm at valence TCutPNO = 10−9.
A . | . | . | wA . |
---|---|---|---|
1H | −1 | 12 | 6.5 |
13C | 0 | 200 | 100 |
15N | 0 | 900 | 450 |
17O | −40 | 1120 | 580 |
19F | −300 | 400 | 350 |
33S | −290 | 670 | 480 |
31P | −180 | 250 | 215 |
77Se | −1000 | 2000 | 1500 |
A . | . | . | wA . |
---|---|---|---|
1H | −1 | 12 | 6.5 |
13C | 0 | 200 | 100 |
15N | 0 | 900 | 450 |
17O | −40 | 1120 | 580 |
19F | −300 | 400 | 350 |
33S | −290 | 670 | 480 |
31P | −180 | 250 | 215 |
77Se | −1000 | 2000 | 1500 |
Setting TCutPNO to a conservative value, we examine the effect of FCut in Fig. 6. Apparently, it imparts negligible errors on the shieldings, especially at the default value of 10−5.
The influence of TCutDO (at a conservative value of TCutPNO) 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 TS = 10−8, which leads to some numerical issues, as discussed in Sec. III B, seen in the figure as large kinks in the curves for some molecules. In the right subplot, a value of TS = 10−5 was used, which leads to smooth convergence or the results toward the RI-MP2 reference. At the default value of TCutDO = 10−2, the MAWE over the whole dataset is under 0.05%.
Turning to the pair prescreening threshold TCutPre, its effect on the shielding MAWE is shown in Fig. 8. Despite the relatively small systems, using a large TDOij = 0.05 leads to a significant number of pairs to be screened-out: 50%–80% at TCutPre = 10−4 and 20–50% at TCutPre = 10−6. Even so, the effect on the shieldings is minor. Neglecting the contributions from the prescreening correction ΔEPre also does not worsen the results significantly for TCutPre ≤ 10−4; however, calculating these contributions takes a very small proportion of the computation time, so there are no savings to be made by skipping them. On the other hand, the default threshold values of TDOij = 10−5 and TCutPre = 10−6 might be too conservative for NMR shielding calculations.
Having separately examined the effects of the most important DLPNO parameters on the shieldings, we now see what happens when all parameters are set to their default (NormalPNO) values. Figure 9 shows the error distributions for different nuclides over all eight test systems in Fig. 3. The NormalPNO/RIJONX distributions are fairly narrow, and both the medians and the ranges (which are more relevant for chemical shift calculations) are much smaller than the inherent error in MP2. There is an overall bias toward overestimation of the shielding constants, so some error compensation can be expected in chemical shift calculations. The same figure also contains results with the LoosePNO threshold settings, which produce both larger systematic deviations and wider error distributions, at the borderline of what we may consider tolerable.
All calculations discussed so far were performed using the RIJONX approximation. In practice, it is preferable to also approximate the exchange part of the Fock matrix, e.g., using the COSX approximation, as this significantly reduces the computation time. However, the semi-numerical integration used in COSX introduces additional errors and numerical noise to the calculation, which may get amplified in DLPNO calculations. Therefore, we have included in Fig. 9 the error distributions from calculations employing the RIJCOSX approximation with a rather large set of integration grids (corresponding to the DefGrid3 keyword). Although the DLPNO and COSX errors are basically cumulative, the latter are fairly small, so the conclusions from the previous paragraph still apply for the most part. One exception is the 1H 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 TCutPNO and TCutDO from LoosePNO to NormalPNO to TightPNO systematically reduces the errors toward the RI-MP2 reference.
As the effect of the local approximations becomes more significant with the increasing system size, the relative error (RE) in the properties can also be expected to increase, up to the point where all approximations are fully active: orbital domains have reached a constant size, etc. For an intrinsic (i.e., local) property such as NMR shielding, this means that absolute errors should also reach an upper limit for large enough systems. To examine this relationship, we use the calculations on the glycine chain systems, discussed in Sec. V D. Figure 10 shows the maximum absolute errors in the shieldings of different nuclides, as a function of the chain length. Apparently, the errors do increase slightly with the system size, and the effect is more pronounced for the less sparse α-helix conformation. However, the errors do not grow indefinitely, leveling off around (gly)10, and are within the ranges shown in Fig. 9. A notable outlier is an oxygen nucleus in linear (gly)3 with a rather large error of 0.97 ppm for 17O vs 0.25 and 0.41 ppm for (gly)2 and (gly)4, respectively. Using the level shift-based adjustment, discussed in Sec. III C, the MaxAE for 17O in linear (gly)3 is reduced to 0.42 ppm, consistent with the slightly oscillating trend.
C. Accuracy: Dipole polarizabilities
As evident from Sec. II, the same general procedure can be used to calculate both magnetic shielding and dipole polarizability tensors. It is worth noting that the polarizability is an extensive global property of the system and as such is less useful to obtain for large molecules, which are the main targets of local correlation methods. Nonetheless, we have included the capability in our implementation and will evaluate the accuracy of the DLPNO-MP2 polarizabilities in this section. It is most convenient to use the relative error (RE) in the isotropic polarizability,
and its absolute value (ARE). This should be compared to the inherent error of MP2 vs more accurate methods such as CCSD(T), which is around 2% for closed-shell species.73 DHDFT performs slightly better but still over 1%. Therefore, we consider a DLPNO error of up to 0.5% to be reasonable.
Figure 11 shows the dependence of the error on TCutPNO for several molecules. The DLPNO-MP2 polarizabilities correctly converge toward the RI-MP2 reference values and are accurate enough for TCutPNO ≤ 10−7. The influence of PNO truncation for core electrons is almost negligible, which is to be expected, since they do not contribute significantly to the polarizability, especially since we did not include core-polarization basis functions in these calculations. We also note a number of kinks in the curves on the left—one for coronene at TCutPNO = 10−12 and two for ATP4− at 10−7 and 10−10 (the latter also leads to the slightly larger error in the right subplot). We associate these discrepancies with numerical instabilities due to two factors. The first one is near linear dependencies in the AO basis: even though TS was set to 10−5 for these calculations, the analogous threshold for the AO basis was 10−8 and the lowest overlap eigenvalues for caffeine, penicillin, coronene, and ATP4− were 2.7 × 10−7, 4.2 × 10−7, 2.2 × 10−8, and 4.6 × 10−7, respectively. The second factor is small differences in the occupation numbers of PNOs and CPNOs, as discussed in Sec. V B. Indeed, either using a cutoff of 10−6 for the MO orthonormalization or applying the level shift, introduced in Sec. III C, both eliminate these discrepancies (see the supplementary material).
The effect of domain truncation on the polarizability is shown in Fig. 12. As for the shieldings, a low value of TS results in erratic behavior and huge errors, especially for coronene, which already has linear dependency issues in the AO basis. At TS = 10−5, however, convergence toward the RI-MP2 reference is much smoother, and the errors are already sufficiently small at TCutDO = 3 × 10−2.
The influence of FCut on the polarizability, shown in Fig. 13, is rather small. Apparently, the default value of 10−5 is more than conservative enough.
Convergence with the pair prescreening threshold TCutPre, as well as the significance of the terms derived from ΔEPre, is shown in Fig. 14. With the default parameters, only 18% of pairs are skipped for the largest system, ATP4−; however with the high value of TCutDOij used here, between 35% and 67% of pairs are screened out already at TCutPre = 10−5, and the errors in the polarizability are below 0.2%. As for the shieldings, the dipole-based correction terms only make a difference for very high values of the threshold. The somewhat larger errors for coronene at TCutPre ≥ 10−4 are probably due to the strong delocalization in the system, together with the large value of TCutDOij.
The combined effect of all DLPNO approximations with the default NormalPNO parameters (and TS = 10−5, εscale = 0) is shown in Fig. 15. The signed error is plotted to illustrate that DLPNO-MP2 consistently underestimates the RI-MP2 polarizabilities. In turn, MP2 tends toward overestimation with respect to CCSD(T),73 so some fortuitous error compensation can be expected with DLPNO-MP2. In any case, the DLPNO-MP2 values only deviate from the RI-MP2 reference by 0.1%–0.2%, which is well within our 0.5% tolerance.
D. Computational efficiency
To compare the efficiency of our RI-MP2 and DLPNO-MP2 implementations, we performed calculations on a series of polyglycine chains in both linear and α-helical conformations. The latter are much denser and expected to be a bigger challenge for the local approximations. The total wall-clock times for the shielding calculations, including SCF and CPSCF solutions, are presented in Fig. 16. The steep scaling of the RI-MP2 calculations is apparent, and there is no difference between the two conformations, which is expected, as sparsity is not exploited apart from basic and rather conservative pre-screening. In contrast, both series of DLPNO calculations achieve quadratic effective scaling, although the prefactor for the helical conformation is about 4.5 times higher. The crossover point between RI- and DLPNO-MP2 is thus surpassed at (gly)5 (38 atoms, 160 electrons, and 994 basis functions) and (gly)9 (66 atoms, 280 electrons, and 1738 basis functions) in the linear and helical series, respectively. In other words, calculations that take more than a day or two (on eight cores) with RI-MP2 can be performed more efficiently with DLPNO-MP2. For the largest system reported here—(gly)15 (108 atoms, 460 electrons, and 3040 basis functions)—the difference is substantial: 1–3 days with DLPNO-MP2 (depending on the conformation) vs 36 days with RI-MP2. Before the crossover point, the DLPNO-MP2 calculation is at most two times slower [for the helical (gly)5] than the RI-MP2 one. Another worthwhile comparison is given in Fig. 17, namely, how much more expensive a shielding calculation is at the MP2-level vs the SCF-level (HF or hybrid DFT). Apparently, with NormalPNO and RIJCOSX, the DLPNO-MP2 shieldings are obtained at only 5–20 times the cost of the preceding HF shielding calculation (the latter is included in the DLPNO-MP2 timings), and the ratio is still going down at (gly)15. This is in stark contrast to RI-MP2, where the calculations with 400 correlated electrons are already more than 100 times more expensive than HF or hybrid DFT and hardly feasible beyond that. Based on the calculations performed for Fig. 9, the TightPNO (LoosePNO) settings increase (decrease) the total computation time by a factor of 1.3–2.2, compared to NormalPNO.
Detailed information on the most time-consuming steps for an even larger system—the vancomycin molecule—is given in Table V. In this case, a frozen-core calculation was performed with a DHDFT functional (DSD-PBEP86), also including implicit solvent DMSO effects using CPCM. Some key takeaways from the data are as follows: (1) The total time for the DLPNO-DHDFT calculation is only eight times longer than an equivalent hybrid DFT calculation (SCF plus SCF-level NMR timings). (2) The cost to calculate the DLPNO-MP2 response densities is roughly six times higher than the DLPNO-MP2 energy and density calculation, as for RI-MP2.72 (3) Reconstructing the SC amplitudes, NPAO space, and the perturbed (C)PNO coefficients twice takes about a third of the DLPNO-MP2 response time. (4) The equations related to orbital localization make up less than 1% of the computational time. (5) Timings for the evaluation of DFT functional, CPCM, and dipole-based prescreening terms are almost negligible and therefore not shown.
Atoms . | 176 . |
---|---|
Basis size (pcSseg-2) | 4 700 |
AuxJ size (def2/JK) | 9 097 |
AuxC size (cc-pwCVTZ/C) | 13 591 |
Core electrons | 218 |
Valence electrons | 542 |
Kept orbital pairs | 22 265 |
Mean PAO domain size | 1 013 |
Mean PNO domain size | 22 |
Mean AuxC domain size | 1 852 |
Total time | 4 791 |
SCF (15 cycles) | 355 |
SCF-level NMR | 242 |
→CPSCF right-hand size | 56 |
→CPSCF (8 cycles) | 175 |
DLPNO-MP2 energy + density | 579 |
transformation | 23 |
and dij construction | 88 |
→PNO amplitude iterations | 35 |
→NPAO space reconstruction | 163 |
→PNO-specific terms | 70 |
contributions | 84 |
→Fock response terms | 10 |
→Z-CPSCF (10 cycles) | 56 |
→ Localization + Z-CPL | 5 |
DLPNO-MP2 response density | 3 603 |
transformation | 276 |
and d″ij,λ construction | 809 |
→Tij,λ iterations | 160 |
→NPAO space reconstruction | 270 |
→PNO-specific terms | 671 |
and 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 |
transformation | 23 |
and dij construction | 88 |
→PNO amplitude iterations | 35 |
→NPAO space reconstruction | 163 |
→PNO-specific terms | 70 |
contributions | 84 |
→Fock response terms | 10 |
→Z-CPSCF (10 cycles) | 56 |
→ Localization + Z-CPL | 5 |
DLPNO-MP2 response density | 3 603 |
transformation | 276 |
and d″ij,λ construction | 809 |
→Tij,λ iterations | 160 |
→NPAO space reconstruction | 270 |
→PNO-specific terms | 671 |
and contributions | 748 |
→Fock response terms | 127 |
→Perturbed Z-CPSCF (8 cycles) | 176 |
→CPL + perturbed Z-CPL | 27 |
A final test is shown in Table VI. Here, DLPNO-MP2 polarizability calculations were carried out for a series of [n]helicene molecules. The results are compared to those obtained by Friese et al. with their Laplace-based RI-MP2 implementation.59 The calculated isotropic polarizabilities agree to within 0.3%, as expected. It is not fair to compare the timings directly due to the different computer architectures, but we can look at the trends and effective scaling of both implementations with the system size. Based on the calculations for n = 7–11, the RI-MP2 algorithm scales as , roughly equivalent to the 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 .
. | DLPNO-MP2 . | RI-MP259 . | |||
---|---|---|---|---|---|
N . | Nbas . | αiso . | ta . | αiso . | tb . |
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-MP259 . | |||
---|---|---|---|---|---|
N . | Nbas . | αiso . | ta . | αiso . | tb . |
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 |
On 12 Intel Xeon E5-2687Wv3@3.10GHz CPU cores with 8 GB of RAM per core.
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 13C and 1H 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)
pair natural orbitals (PNOs)
discarded/“complementary” pair natural orbitals (CPNOs)
union of PNOs and CPNOs
orbital/pair domain of orthonormal pseudo-canonical non-redundant projected atomic orbitals (NPAOs)
normalized redundant projected atomic orbitals (PAOs)
APPENDIX: DERIVATION OF THE PNO ROTATION STATIONARITY CONDITIONS
In the present work, the MO and PNO response is parameterized using full, complex-valued matrices, constrained with an explicit orthonormality condition, in contrast to the exponential parameterization used previously for DLPNO-MP2.69 This has the consequence that the final working equations for electric/magnetic perturbations, which are purely real/imaginary, must be obtained as linear combinations of several Lagrangian stationarity conditions. In this section, we show a detailed derivation of the equations related to PNO rotations and briefly discuss the analogous treatment of MO rotations.
The PNO rotation stationarity conditions have the following form:
The separate terms are as follows:
where quantities in the full PNO + CPNO basis are expanded as
The index of the intermediate spans the PAO domains of all pairs ik and jk included in the sum, i.e., those for which or , respectively, is greater than FCut. Note also that the amplitudes Tij are implicitly non-zero only in the PNO block. To obtain equations for vij, we take the following linear combination, where the contributions from CSC and CPNOO drop out:
Only rotations between PNOs and CPNOs result in the non-zero equations shown in Eq. (58). In addition, we use Eq. (55) (at λ = 0) to express Eq. (A5) as
and after expanding the semi-canonical residuals in Eq. (A7), we obtain a different expression for Eq. (A1),
where, padding Dij with zeros up to the full PNO + CPNO basis, we have substituted
If we define , we obtain a solution for xij,
Note that we obtain a different expression from the Hermitian conjugate,
but both equations are equivalent, provided that Eq. (A12) is fulfilled.
To obtain the perturbed multipliers vij,λ, we first express the perturbed PNO/CPNO stationarity conditions as
in which xij,λ cancels out to obtain Eq. (107). A solution for xij,λ can also be obtained from Eq. (A18),
This is equal to the Hermitian conjugate, obtained from Eq. (A19), provided that Eq. (A20) is fulfilled.
The derivation of the MO rotation stationarity conditions in analogous—the contributions from CMOO cancel in the linear combination,
which gives the Z-CPL, Z-CV, and Z-CPSCF equations (59), (67), and (69), respectively, while solutions for x can be obtained from either equation separately,
which are equivalent, provided that Eq. (A22) is fulfilled. The second derivatives are also treated analogously—we first express the stationarity conditions as
and then obain the perturbed Z-CPL, Z-CV, and Z-CPSCF equations (109), (112), and (113), respectively, from the linear combination in which xλ cancels,