Communication: Multipole approximations of distant pair energies in local correlation methods with pair natural orbitals

Werner, Hans-Joachim

doi:10.1063/1.4968595

The accuracy of multipole approximations for distant pair energies in local second-order Møller-Plesset perturbation theory (LMP2) as introduced by Hetzer et al. [Chem. Phys. Lett. 290, 143 (1998)] is investigated for three chemical reactions involving molecules with up to 92 atoms. Various iterative and non-iterative approaches are compared, using different energy thresholds for distant pair selection. It is demonstrated that the simple non-iterative dipole-dipole approximation, which has been used in several recent pair natural orbitals (PNO)-LMP2 and PNO-LCCSD (local coupled-cluster with singles and doubles) methods, may underestimate the distant pair energies by up to 50% and can lead to significant errors in relative energies, unless very tight thresholds are used. The accuracy can be much improved by including higher multipole orders and by optimizing the distant pair amplitudes iteratively along with all other amplitudes. A new approach is presented in which very small special PNO domains for distant pairs are used in the iterative approach. This reduces the number of distant pair amplitudes by 3 orders of magnitude and keeps the additional computational effort for the iterative optimization of distant pair amplitudes minimal.

The efficiency of local electron correlation methods depends crucially on a simplified treatment of the correlation of electrons in spatially distant localized molecular orbitals (LMOs) $ϕ_{i}$ and $ϕ_{j}$ ⁠. It is well known that the leading term of the distant pair correlation energies E_ij decays with $R_{i j}^{- 6}$ ⁠, where R_ij is the distance between the charge centers of the two LMOs. Typically, for $R > 12 a_{0}$ the distant pair energies become individually very small (⁠ $< 1 μ E_{h}$ ⁠). However, since the number of distant pairs increases quadratically with the number of correlated LMOs, the total energy contribution of distant pairs in large molecules can be substantial. In this communication, we describe various approximations to compute distant pair energies and investigate how these affect correlation and reaction energies of some larger systems.

In closed-shell local second-order Møller-Plesset perturbation theory (LMP2), the pair correlation energies for an orbital pair $i \geq j$ are given by

E_{i j} = (2 - δ_{i j}) \sum_{a, b \in [i j]} T_{a b}^{i j} [2 K_{a b}^{i j} - K_{b a}^{i j}],

(1)

where the indices a, b refer to suitable local virtual orbitals [e.g., projected atomic orbitals¹ (PAOs), orbital specific virtuals² (OSVs), or pair natural orbitals^3,4 (PNOs)]. $T_{a b}^{i j}$ are the amplitudes, and $K_{a b}^{i j} = (a i | b j)$ the two-electron repulsion integrals. The summation is limited to a subspace (pair domain) [ij] of virtual orbitals a, b that are spatially close to the LMOs i or j. We will assume that the orbitals in each of these domains are “pseudo-canonical,” i.e., that they diagonalize the Fock matrix within the domain, $F_{a b}^{i j} = ⟨ a^{i j} | \hat{f} | b^{i j} ⟩ = δ_{a b} \in_{a}^{i j}$ for $a^{i j}, b^{i j} \in [i j]$ ⁠. These orbitals are pair-specific, and within a domain they are orthonormal. However, the orbitals of different pair domains [ij] and [kl] are mutually non-orthogonal, with overlap matrix ${[𝐒^{[i j, k l]}]}_{a b} = ⟨ a^{i j} | b^{k l} ⟩$ ⁠. In the following, we will omit the superscripts of the orbital labels, since it should always be obvious to which pairs they belong.

The amplitudes are determined by solving first-order LMP2 equations, $R_{a b}^{i j} = 0$ for $a, b \in [i j]$ ⁠, with⁵

R_{a b}^{i j} = K_{a b}^{i j} + (ϵ_{a}^{i j} + ϵ_{b}^{i j} - f_{i i} - f_{j j}) T_{a b}^{i j} + G_{a b}^{i j} + G_{b a}^{j i}, G_{a b}^{i j} = - \sum_{k \neq i} f_{i k} {[𝐒^{[i j, k j]} 𝐓^{k j} 𝐒^{[k j, i j]}]}_{a b} .

(2)

Approximate amplitudes can be obtained by neglecting the coupling terms $G_{a b}^{i j}$

T_{a b}^{i j} \approx - \frac{K_{a b}^{i j}}{ϵ_{a}^{i j} + ϵ_{b}^{i j} - f_{i i} - f_{j j}} .

(3)

This is denoted as “semi-canonical” (SC) amplitude approximation, since the Fock matrix is not diagonal in the occupied space.

If i, j, a, b are local, the integrals (ai|bj) decay exponentially with the distances R_ai and R_bj. Thus, for distant pairs, one can limit the virtual labels to $a \in [i]$ and $b \in [j]$ ⁠, where [i] and [j] denote domains of virtual orbitals that are spatially close to the LMOs i and j, respectively. This means that distant pair energies are approximated as

E_{i j} \approx 4 \sum_{a \in [i], b \in [j]} T_{a b}^{i j} K_{a b}^{i j} .

(4)

Charge transfer contributions $a, b \in [i]$ or $a, b \in [j]$ ⁠, as well as exchange contributions $a \in [j]$ and $b \in [i]$ ⁠, are neglected. The impact of this “NOEX” approximation will be demonstrated later on. The amplitudes $T_{a b}^{i j}$ can either be approximated by Eq. (3) (“semi-canonical distant pair approximation”) or determined iteratively by solving the amplitude equations $R_{a b}^{i j} = 0$ for $a \in [i]$ and $b \in [j]$ (“iterative distant pair approximation’’).

The integrals $K_{a b}^{i j}$ describe the Coulomb interaction of effective charge distributions $ρ_{a i} (𝐫_{1})$ and $ρ_{b j} (𝐫_{2})$ ⁠,

K_{a b}^{i j} = \int d 𝐫_{1} \int d 𝐫_{2} ρ_{a i} (𝐫_{1}) \frac{1}{r_{12}} ρ_{a i} (𝐫_{2}),

(5)

with $ρ_{a i} (𝐫) = ϕ_{a} (𝐫) ϕ_{i} (𝐫)$ ⁠. Due to the orthonormality of the virtual and occupied orbitals, these charge distributions carry no charge, i.e., $\int ρ_{a i} (𝐫) d 𝐫 = 0$ ⁠. Thus, the lowest-order contribution at long-range is the dipole-dipole interaction, which decays with $R_{i j}^{- 3}$ ⁠. Consequently, the pair energies E_ij decay with $R_{i j}^{- 6}$ ⁠.

For non-overlapping charge distributions, these integrals can be approximated by multipole approximations, as first proposed in the LMP2 context by Hetzer et al.⁶ It is assumed that the charge centers of the charge distributions of $ρ_{a i} (𝐫)$ and $ρ_{b j} (𝐫)$ coincide with those of the LMOs, $𝐑_{i} = ⟨ i | 𝐫 | i ⟩$ and $𝐑_{j} = ⟨ j | 𝐫 | j ⟩$ ⁠, respectively. Let R = R_j $-$ R_i, R = |R|, and $\hat{𝐑} = 𝐑 / R$ ⁠. The integrals are then approximated as

K_{a b}^{i j} \approx {(a i | b j)}^{(p)} = \sum_{n = 2}^{p + 1} \frac{1}{R^{n + 1}} \sum_{m = 1}^{n - 1} \sum_{μ ν} A_{μ ν}^{(m, n - m)} (\hat{𝐑}) \times M_{a i}^{(m, μ)} (R_{i}) M_{b j}^{(n - m, ν)} (R_{j}) .

(6)

In this expansion all terms that decay with up to R^p+2 are included, and p is the highest order of the multipole moments involved. $A_{μ ν}^{(m, n - m)} (\hat{𝐑})$ is the “interaction matrix,” which only depends on $\hat{𝐑}$ ⁠. The index $μ$ runs over the (m + 1) (m + 2)/2 unique Cartesian multipole moment components of order m, and $M_{a i}^{(m, μ)} (𝐑_{i})$ are the integrals over these operator components, evaluated at the origin R_i. The sum over ν is restricted accordingly. Note that we use Cartesian multipole moment operators rather than traceless multipole operators as in other formulations.^7,8 More details of the derivation and explicit expressions up to p = 3 are given in the supplementary material.

The lowest-order approximation (p = 1) is the dipole-dipole approximation. Increasing the order p improves the accuracy at long range but leads to faster divergence at short range. Hetzer et al.⁶ therefore recommended to truncate the expansion at p = 3, which gives the best compromise of stability, accuracy, and efficiency.

Using the multipole approximations described above, the semi-canonical LMP2 pair energies $E_{i j}^{sc}$ [cf. Eqs. (3) and (4)] can be evaluated on-the-fly without storing the integrals. Though not linear scaling, this is very fast, even for large molecules. It is therefore convenient to compute the pair energies for all pairs for which $R_{i j} > R_{min}$ (using, e.g., $R_{min} = 5 a_{0}$ ⁠) and then select those for which $E_{i j}^{sc} < T_{dist}$ as distant pairs. All remaining pairs are fully included in the LMP2 calculation. The approximated distant pair energies are summed up to yield the total distant correlation energy $E_{dist}^{sc}$ ⁠, which is added to the final LMP2 energy. This approach has been used in the PNO-LMP2 and PNO-LCCSD (local coupled-cluster with singles and doubles) methods of Riplinger et al.^9,10 and ours,^11–13 as well as in recent PNO-NEVPT2¹⁴ and PNO-CASPT2¹⁵ methods. In all these works, only the semi-canonical dipole-dipole approximation was employed. Either PAO or OSV domains were used in these methods to compute the distant pair energy.

Alternatively, a better approximation is obtained by storing the approximate integrals (ai|bj)^(p) (⁠ $a \in [i]$ ⁠, $b \in [j]$ ⁠) and optimizing the corresponding amplitudes and pair energies iteratively along with all other pairs.

As will be shown below, the iterative approach is significantly more accurate than the semi-canonical one. However, it is also more expensive, since the semi-canonical PAO or OSV orbital domains [i] and [j] may be rather large. For accurate calculations using triple-ζ basis sets, the PAO and OSV orbital domains typically include 300 and 150 virtual orbitals, respectively. Since the number of distant pairs is very large, a straightforward use of these domains causes very high storage and central processing unit (CPU) demands.

This problem can be avoided by defining distant pair PNO domains [i]_dPNO and [j]_dPNO for each pair ij. In contrast to the PAO or OSV orbital domains, these are specific for each distant pair, and very few PNOs per pair domains are sufficient to optimally describe the correlation of the distant pair electrons. This strongly reduces the computational effort for optimizing the distant pair amplitudes. In general, PNOs are obtained by diagonalizing external pair density matrices (⁠ $\tilde{𝑻}$ ^ij = $2 𝐓^{i j} - 𝐓^{i j †})$

𝐃^{i j} = {\tilde{𝑻}}^{i j †} 𝐓^{i j} + {\tilde{𝑻}}^{i j} 𝐓^{i j †} (i \neq j)

(7)

in the basis of pseudo-canonical OSV pair domains [ij]_OSV. For distant pairs, the amplitudes $T_{a b}^{i j}$ are only non-zero for $a \in {[i]}_{OSV}, b \in {[j]}_{OSV}$ ⁠, and obtained using Eq. (3) and the multipole approximation of the integrals $K_{a b}^{i j}$ ⁠. Here [i]_OSV and [j]_OSV denote pseudo-canonical OSV orbital domains. The pair density then becomes block-diagonal, with the two blocks

{[𝐃^{i j}]}_{a b} = {[2 𝐓^{i j †} 𝐓^{i j}]}_{a b}, a, b \in {[j]}_{OSV},

(8)

{[𝐃^{i j}]}_{a b} = {[2 𝐓^{i j} 𝐓^{i j †}]}_{a b}, a, b \in {[i]}_{OSV} .

(9)

These blocks can be diagonalized separately, and the eigenvalues can be used to select distant pair dPNO domains [i]_dPNO and [j]_dPNO, respectively, using a PNO occupation number threshold $T_{PNO}^{occ}$ ⁠. Alternatively, a completeness criterion $T_{PNO}^{en}$ can be employed, and then as many PNOs are included as needed to recover for each pair at least $E_{i j}^{PNO} \geq T_{PNO}^{en} \times E_{i j}^{OSV}$ ⁠, where the semi-canonical pair energies are used (see the supplementary material for more details of this procedure). Since the occupation numbers of distant pairs are very small, we found the latter criterion more suitable and used $T_{PNO}^{en} = 0.997$ for all calculations presented in this paper. Test calculations with a stricter threshold of 0.999 yielded nearly the same results (see the supplementary material).

Finally, the new domains are made pseudo-canonical, and the integrals are transformed accordingly. The dPNO domains are very much smaller than the original OSV domains. Often only 3-6 functions are sufficient to recover more than 99.5% of the OSV pair energies, and therefore both the memory demand and the CPU time needed to optimize the distant pair amplitudes are dramatically reduced. Nevertheless, as will be demonstrated in the following, the results are almost unchanged. The transformation to the dPNO domains can be done on-the-fly when generating the integrals by the multipole approximations, so that only integrals for the very small domains need to be stored.

We demonstrate the effect of the various approximations for the dissociation of a gold(I)-aminonitrene complex (AuC₄₁H₄₅N₄P, for simplicity denoted Auamin, see Fig. 1), which has been used in previous benchmarks of PNO-LCCSD methods.^12,16 This is an excellent and difficult benchmark case since the correlation effects on the reaction energy are very large (Hartree-Fock, LMP2-F12, and LCCSD-F12 yield 92, 251, and 190 kJ mol⁻¹, respectively, using the VTZ-F12 basis set¹⁶). The results can be directly compared to an experimental gas-phase value¹⁷ of 196.5 $\pm 11.2$ kJ/mol, which has been obtained by subtracting the PW91/cc-pVTZ-pp zero-point correction¹⁷−8.2 kJ/mol from the measured value. The distant pair contribution to the reaction energy is the largest of all systems we have studied so far. Additional results for two other reactions are presented in the supplementary material.

FIG. 1.

View large Download slide

The Auamin reaction.

For all atoms except gold, we use the VTZ-F12 basis set,¹⁸ and for gold the aug-cc-pVTZ-PP basis¹⁹ with the ECP60MDF effective core potential.²⁰ The PAO, OSV, and PNO domains were selected as described in Ref. 11, using well established default thresholds (see the supplementary material for details).

The contributions $Δ E_{dist}$ of distant pairs to PNO-LMP2 reaction energies are presented in Table I for p = 1, 2, 3 and three distant pair selection thresholds of T_dist. Both the semi-canonical (non-iterative) and iterative multipole approximations have been used. For the ease of comparison, the total distant pair contributions obtained without multipole approximations are given in column $Δ E_{dist}$ ⁠; the remaining columns show the errors relative to these values that are introduced by the various approximations. In all cases, the distant pairs have been selected using the semi-canonical p = 1 pair energies, i.e., the number of distant pairs is independent of the approximation applied. In the cases where the distant pair amplitudes are optimized iteratively, all distant pairs were included. For comparison, also reaction energies computed without multipole approximations, but still neglecting the charge-transfer and exchange type integrals are shown (NOEX approximation).

TABLE I.

PNO-LMP2 distant pair energy contributions $Δ E_{dist}$ to the reaction energies without multipole approximations and errors relative to these values caused by various multipole approximations. All energies in kJ mol⁻¹, thresholds T_dist in $μ E_{h}$ ⁠. The numbers over the columns correspond to the multipole order applied. NOEX refers to calculations without multipole approximation but neglected exchange terms. N_dist is the number of distant pairs; the total number of pairs is 7503. $N_{OSV}^{av}$ and $N_{dPNO}^{av}$ are average orbital domain sizes for distant pairs. The computed total PNO-LMP2 reaction energy (without distant pair approximations) amounts to 239.7 kJ mol⁻¹.

					Errors of reaction energies
					Semi-canonical			Iterative			NOEX
T_dist	N_dist	$N_{OSV}^{av}$	$N_{dPNO}^{av}$	$Δ E_{dist}$	1	2	3	1	2	3	PNO	OSV
10	4812	143	3.8	22.02	−10.20	−7.13	−6.22	−7.91	−4.12	−2.89	−2.92	−2.86
3	3609	143	3.3	6.93	−2.71	−1.91	−1.77	−1.82	−0.79	−0.58	−0.55	−0.50
1	2355	143	3.1	1.98	−0.64	−0.46	−0.44	−0.46	−0.13	−0.10	−0.10	−0.08

					Errors of reaction energies
					Semi-canonical			Iterative			NOEX
T_dist	N_dist	$N_{OSV}^{av}$	$N_{dPNO}^{av}$	$Δ E_{dist}$	1	2	3	1	2	3	PNO	OSV
10	4812	143	3.8	22.02	−10.20	−7.13	−6.22	−7.91	−4.12	−2.89	−2.92	−2.86
3	3609	143	3.3	6.93	−2.71	−1.91	−1.77	−1.82	−0.79	−0.58	−0.55	−0.50
1	2355	143	3.1	1.98	−0.64	−0.46	−0.44	−0.46	−0.13	−0.10	−0.10	−0.08

The table demonstrates that the contribution of the distant pairs to the reaction energy is significant (up to 22 kJ mol⁻¹ for $T_{dist} = 10 μ E_{h}$ ⁠). Using the simple semi-canonical dipole-dipole approximation, this contribution is underestimated by 10 kJ mol⁻¹, while with the iterative scheme and p = 3 the error is reduced to 2.9 kJ mol⁻¹. The remaining error is almost exclusively due to the NOEX approximation.

We also compared calculations with the full OSV orbital domains to those with the small dPNO domains. As can be seen in Table I, the reduction of the domain sizes by the dPNO treatment is dramatic. Despite this, more than 99.5% of the distant pair energy is recovered with (on the average) only 3-4 dPNOs per pair, and the reaction energies change negligibly by only 0.02–0.06 kJ mol⁻¹.

The computed absolute distant pair energies of Auamin and some other molecules can be found in the supplementary material. With a large threshold of $10 μ E_{h}$ the approximate distant pair energies are in all cases significantly underestimated. There are three sources of errors: (i) the neglect of charge-transfer and exchange-type excitations for distant pairs (NOEX); (ii) the non-iterative semi-canonical (SC) approximation; and (iii) the multipole approximation for the integrals (MLTP). The largest error (about $- 15 %$ ⁠) arises from the NOEX approximation; the error caused by the SC approximation is also significant (⁠ $- 10 %$ to $- 15 %$ ⁠), while the error of the multipole approximation itself decreases quickly with p and is in all cases below $5 %$ for p = 3. For example, for the Auamin molecule, the errors of the distant pair correlation energy caused by the iterative MLTP approximation (relative to NOEX) amount to $- 25.7 %$ ⁠, $- 7.9 %$ ⁠, and $0 %$ for p = 1, 2, 3, respectively (T_dist = 10^{$-$ 5} E_h). The semi-canonical dipole-dipole approximation (p = 1) recovers only $55 %$ of the exact PNO-LMP2 distant pair energy for Auamin.

For the small threshold of $1 μ E_{h}$ ⁠, the errors caused by the NOEX approximation are reduced but still amount to about $- 7 %$ for the distant pair correlation energy of Auamin. The relative errors of the other approximations are quite similar as for the larger threshold, but the absolute errors are reduced by one order of magnitude.

Even though only 3 reactions have been investigated in the current work, we believe that these results are representative and sufficient to conclude that the numerous, but individually very small energy contributions of distant pairs in large systems are by no means negligible when summed up. Naturally, the effects of distant pair approximations grow with the molecular size. Furthermore, distant pairs have particularly large effects in systems with many aromatic groups, as demonstrated here by the Auamin reaction. It is therefore important to treat distant pairs as accurately as possible in a most cost effective way. The simple dipole-dipole approximation, which has been applied in previous works, is of limited accuracy. Higher multipole orders and the iterative optimization of distant pair amplitudes lead to significant improvements. It appears that $T_{dist} = 3 \times 10^{- 6} E_{h}$ is a safe choice. Very distant pairs with energies smaller than $T_{vdist} = 3 \times 10^{- 7} E_{h}$ can be treated by the non-iterative semi-canonical approximation without loss of accuracy (see the supplementary material). This restores linear scaling for the optimization of the remaining pairs in the iterative PNO-LMP2.

Finally we note that weak and distant pair energies computed by LMP2 are often significantly (typically 30%) overestimated relative to accurate coupled-cluster values. This effect can be reduced by using spin-component scaled²¹ (SCS-LMP2) distant pair energies. The use of LMP2 or SCS-LMP2 approximations for weak pairs (⁠ $10^{- 3} > E_{i j} > T_{dist}$ ⁠) in local coupled-cluster calculations can lead to large errors and is no longer recommended.¹² Significant computational savings for close and weak pair classes are possible, however, using approximate coupled-cluster treatments, as described in Ref. 12. Furthermore, basis set incompleteness errors may be large, even when triple-ζ basis sets are used. For the Auamin PNO-LMP2/VTZ-F12 reaction energy, the basis set error amounts to about 12 kJ mol⁻¹. This error can be very effectively reduced by explicit correlation (F12) treatments; at the same time, the F12 terms minimize errors caused by the domain approximations.^{11–13,22,23} This makes it possible to compute highly accurate LCCSD(T)-F12 energies for molecules with more than 100 atoms.

See supplementary material for explicit expressions for the multipole approximations up to p = 3 which are provided, and benchmark results for three reactions are presented.