We report an implementation of the molecular gradient using the divide-expand-consolidate resolution of the identity second-order Møller-Plesset perturbation theory (DEC-RI-MP2). The new DEC-RI-MP2 gradient method combines the precision control as well as the linear-scaling and massively parallel features of the DEC scheme with efficient evaluations of the gradient contributions using the RI approximation. We further demonstrate that the DEC-RI-MP2 gradient method is capable of calculating molecular gradients for very large molecular systems. A test set of supramolecular complexes containing up to 158 atoms and 1960 contracted basis functions has been employed to demonstrate the general applicability of the DEC-RI-MP2 method and to analyze the errors of the DEC approximation. Moreover, the test set contains molecules of complicated electronic structures and is thus deliberately chosen to stress test the DEC-RI-MP2 gradient implementation. Additionally, as a showcase example the full molecular gradient for insulin (787 atoms and 7604 contracted basis functions) has been evaluated.

Electronic energy derivatives with respect to nuclear coordinates play a central role in contemporary computational chemistry as they open a possibility of calculating a variety of properties.1 For instance, the first energy derivatives are used in the determination of equilibrium and saddle point geometries of molecules, and the second derivatives are used for calculating force constants, harmonic vibrational frequencies, and infrared (IR) and Raman intensities. Higher derivatives are the basis for calculating more subtle effects revealing themselves in spectroscopy as vibration-rotation constants, anharmonic constants, sextic distortion constants, etc.2 

The theory of self-consistent field (SCF) energy derivatives has its roots in the work of Bratoz.3 More modern developments with major applications were brought forward by the work of Pulay.4 In his pioneering work, Pulay demonstrated that properly implemented analytic SCF gradients are not only more accurate but also more efficient than numerical derivatives. This work laid the foundation for the correlated wave function methods, which then followed for a variety of models—Møller-Plesset perturbation theory (MP2,5 MP3,6 MP47), configuration interaction (CI),8–11 and coupled-cluster (CC)12–15 gradients.

In this work we focus on the MP2 molecular gradient. Formulations of the MP2 molecular gradient in the canonical molecular orbital (MO) basis5,16–27 suffer from a high-order computational scaling with system size and thus cannot be used for large molecular systems. Efficient implementations extend the applicability range,28 but for very large molecular systems MP2 gradients are still out of reach using the traditional MO formulation.

The computational cost of the conventional MP2 gradient implementation can be expressed as kN5, where N is the number of the basis functions in the calculation and k is a prefactor. Obviously, the computational cost may be reduced using two different strategies: (i) reducing the prefactor k and (ii) reducing the N5 scaling to a lower power, ideally to a linear dependence on N. Either of these strategies requires that approximations are introduced in the conventional MP2 algorithm.

One of the most common strategies to reduce the prefactor k is to use the resolution of the identity (RI) approximation for the evaluation of two-electron integrals. The RI method has its roots in the work of Baerends et al.29 It was further developed by Dunlap30,31 and by Almløf32 and Ahlrichs.33,34 The RI method is now routinely used for the evaluation of gradients35–39 and has also been extended to Hessian evaluation.40–42 

To reduce the N5 scaling of the MP2 and other correlated models it is necessary to exploit the local nature of the electron correlation. The idea of utilizing the wave function locality originates in the work on local CC methods by Pulay43 and Saebø and Pulay,44 which has been further developed and matured by the works of Hampel and Werner45 and Schütz and Werner.46–48 Many other methods for obtaining reduced scaling have been proposed, including atomic-orbital-based CC,49–51 the natural-linear-scaling approach,52 the cluster-in-molecule approach,53–55 the smooth local CC approach,56 the divide-and-conquer approach,57 the fragment-molecular-orbital approach,58 the incremental scheme,59,60 Laplace-MP2,61–64 the dynamically screened local correlation method,65 local pair natural orbital approximations,66,67 and orbital-specific virtual orbital approximations.68,69

Our group has recently proposed an alternative local scheme—the divide-expand-consolidate (DEC) method. The method allows the CC energy to be evaluated in a linear-scaling and massively parallel fashion.70–72 In the DEC scheme the standard CC calculation is divided into a set of independent CC fragment calculations each using their own subset of the total orbital space. This is possible because we use localized occupied and localized unoccupied HF orbitals. A crucial feature of the DEC method is the precision control obtained by systematically expanding the local orbital spaces in a black-box manner during the calculation to ensure that the fragment energies are determined to a preset threshold denoted the fragment optimization threshold (FOT).73 This, in turn, defines the total CC correlation energy as a sum of the fragment energies to a preset threshold compared to a standard molecular calculation. Moreover, the DEC scheme can be extended to the geometric derivatives without changing the underlying structure of the method.74 

Several other local correlation methods have been extended to enable the evaluation of the MP2 molecular gradient. For example, MP2 molecular gradients were implemented in the context of the local CC method by Werner and coworkers,75 and this approach was also extended to employ density-fitting techniques.76 The MP2 molecular gradient has also been presented for the fragment-molecular-orbital method.77 Furthermore, Schweizer et al.78 reformulated the MP2 gradient equations in the atomic orbital (AO) basis using Laplace transformations of the energy denominator.79 

In the present article, we report an implementation of the DEC-RI-MP2 gradient as an extension of the DEC-MP2 molecular gradient developed in Ref. 74. The DEC-RI-MP2 inherits the precision control as well as the linear-scaling and massively parallel attributes of the DEC-MP2 gradient. Thus, the RI approximation does not change the formal linear scaling of the DEC method, which comes from the division of the total orbital space into many small local orbital spaces. Rather, the use of the RI approximation lowers the cost of the gradient contribution evaluations for each individual fragment. The combination of the DEC and RI approximations results in a method capable of calculating molecular gradients for very large molecular systems. The errors of the DEC-RI-MP2 gradients are analyzed using a test set of supramolecular complexes, and the results demonstrate the general applicability of the DEC method. The test set contains molecules of complicated electronic structures and is thus deliberately chosen to stress test our DEC-RI-MP2 implementation. Furthermore, as a showcase example the full molecular gradient for insulin is evaluated.

The paper is structured as follows. Section II describes the general DEC-MP2 energy and gradient theory followed by a discussion of the modifications needed for the evaluation of the RI approximated terms. Section III contains numerical results on a test set of very large molecular systems, including a molecular gradient calculation on insulin. In Section IV we give some concluding remarks and future perspectives.

In Sections II A and II B we summarize the main equations defining the DEC-MP2 energy and DEC-MP2 molecular gradient,74 respectively. In Section II C we apply RI approximations to the DEC-MP2 gradient equations to obtain expressions for the DEC-RI-MP2 gradient. Throughout the paper we use the following index conventions.

  • α, β, μ, ν: atomic orbital (AO) indices.

  • i, j, k, l: occupied molecular orbital (MO) indices.

  • a, b, c, d: virtual MO indices.

  • p, q, r, s: MO indices of unspecified occupation.

  • Υ, Λ, Φ: auxiliary AO indices.

The MP2 correlation energy for a closed-shell system may be written as

E corr = i j a b t i j a b ( 2 g a i b j g b i a j ) ,
(1)

where gpqrs is a two-electron repulsion integral (ERI) using the Mulliken notation, gpqrs = (pq|rs). The MP2 amplitudes t i j a b are determined by solving the amplitude equations,

0 = g a i b j + c ( t i j c b F a c + t i j a c F b c ) k ( t k j a b F k i + t i k a b F k j ) ,
(2)

where F is the Fock matrix. Equations (1) and (2) are valid for any choice of optimized HF basis. In the following we assume that the HF orbitals have been localized.

In the DEC scheme, the localized HF orbitals are assigned to atomic sites P, Q, R, … (e.g., based on Löwdin atomic charges73), and the calculation of the correlation energy in Eq. (1) is partitioned into many small and independent fragment calculations. The occupied and virtual MOs assigned to atomic site P are denoted P ¯ and P ¯ , respectively. In this work, we consider two different partitioning schemes. In the occupied partitioning scheme, the correlation energy E corr o is written as

E corr o = P E P o + P > Q Δ E P Q o ,
(3)

where the atomic fragment energy E P o and pair fragment energy Δ E P Q o are given by

E P o = i j P ¯ a b [ P ¯ ] t i j a b ( 2 g a i b j g b i a j ) ,
(4)
Δ E P Q o = ( i P ¯ j Q ¯ + i Q ¯ j P ¯ ) a b [ P ¯ ] [ Q ¯ ] t i j a b ( 2 g a i b j g b i a j )
(5)

and [ P ¯ ] is the set of localized virtual HF orbitals spatially close to atomic site P. The virtual summation restriction is possible due to the locality of the two-electron integrals and MP2 amplitudes when expressed in the localized HF basis.71 

Alternatively, the correlation energy may be evaluated using the virtual partitioning scheme,72 

E corr v = P E P v + P > Q Δ E P Q v ,
(6)
E P v = a b P ¯ i j [ P ¯ ] t i j a b ( 2 g a i b j g b i a j ) ,
(7)
Δ E P Q v = ( a P ¯ b Q ¯ + a Q ¯ b P ¯ ) i j [ P ¯ ] [ Q ¯ ] t i j a b ( 2 g a i b j g b i a j ) ,
(8)

where [ P ¯ ] is the set of localized occupied HF orbitals spatially close to atomic site P.

The amplitudes entering Eqs. (4) and (7) are determined by solving the MP2 amplitude equations in Eq. (2) in a local orbital space where the occupied and virtual orbitals are restricted to the [ P ¯ ] and [ P ¯ ] spaces, respectively. Similarly, the MP2 amplitudes for pair fragments are determined by solving the MP2 amplitude equations in union of spaces from atomic fragment calculations. In practice, [ P ¯ ] and [ P ¯ ] are determined in a black box manner to ensure that the error in the atomic fragment energy E P o is below the fragment optimization threshold (FOT) as described in Ref. 73. The occupied and virtual partitioning schemes in Eqs. (3) and (6), respectively, are equally valid strategies for partitioning the correlation energy into a sum of atomic fragment and pair fragment contributions. Since the local orbital spaces are determined to the FOT precision in a black-box manner, both approaches yield the total correlation energy to FOT precision.72 For the energy evaluation it is thus sufficient to consider either Eq. (3) or Eq. (6). However, for the gradient evaluation, the locality analysis in Ref. 74 shows that some terms are most easily evaluated using the occupied partitioning scheme, while other terms require the virtual partitioning scheme to exploit the locality of the associated MO integrals.

For atomic fragment P, the local HF orbitals are also expanded in a restricted set of atomic orbitals, the so-called atomic extent (AE) {P}, which is a subset of the full set of atomic orbitals in the vicinity of atomic site P. The AE is determined such that the approximate MO ϕ ̃ r P assigned to atomic site P,

ϕ ̃ r P = μ { P } χ μ C ̃ μ r P ,
(9)

resembles the true MO ϕ r P as much as possible in a least squares sense.73 

In principle, the number of fragments to consider scales quadratically with system size. However, the pair fragment energies for distant pairs describe dispersion effects, which decay with the inverse pair distance to the sixth power.72 Pair fragments for distant atoms may therefore be omitted without affecting the precision of the total correlation energy. When pair screening is employed, the DEC-MP2 scheme becomes linear-scaling with system size. In the simplest implementation a distance cutoff can be used to screen away pairs with negligible energy contributions. A more sophisticated pair screening procedure will be presented in a forthcoming publication.

The DEC-MP2 molecular gradient expressions were derived in detail in Ref. 71. Here we summarize the main equations that define the DEC-MP2 molecular gradient to provide the necessary framework for discussing the RI approximations in Section II C.

The DEC-MP2 molecular gradient is obtained by differentiating the MP2 Lagrangian LMP2 with respect to the nuclear displacements x for the 3 Cartesian coordinates for each atom and subsequently applying locality approximations to simplify the evaluation of the gradient intermediates. The DEC-MP2 molecular gradient for a nuclear displacement x, L MP2 ( x ) , may be written as

L MP2 ( x ) = h nuc ( x ) + L Θ ( x ) + L 1-el ( x ) + L cou-ex ( x ) + L reort ( x ) .
(10)

In Eq. (10) and throughout the paper the notation Z(x) refers to a partial derivative, i.e., the AO integrals contained in Z are differentiated with respect to nuclear displacement x, while the wave function parameters (MO coefficients and MP2 amplitudes/multipliers) contained in Z are not. The nuclear repulsion contribution to the energy is denoted h nuc ( x ) , and L Θ ( x ) contains the terms where MP2 multipliers t ̄ i j a b are contracted with differentiated two-electron integrals,

L Θ ( x ) = P L Θ , P ( x ) + P > Q Δ L Θ , PQ ( x ) ,
(11)
L Θ , P ( x ) = i j P ¯ a b [ P ¯ ] t ̄ i j a b g a i b j ( x ) ,
(12)
Δ L Θ , PQ ( x ) = ( i P ¯ j Q ¯ + i Q ¯ j P ¯ ) a b [ P ¯ ] [ Q ¯ ] t ̄ i j a b g a i b j ( x ) ,
(13)

where the MP2 Lagrange multipliers can be calculated directly from the MP2 amplitudes according to

t ̄ i j a b = 4 t i j a b 2 t i j b a .
(14)

The one-electron L 1-el ( x ) , Coulomb/exchange L cou-ex ( x ) , and reorthonormalization L reort ( x ) contributions to the gradient can be expressed as

L 1-el ( x ) = μ ν h μ ν ( x ) ( 2 D + ρ AO ) μ ν ,
(15)
L cou-ex ( x ) = 2 μ ν G μ ν ( x ) ( D ) ( D + ρ AO ) μ ν ,
(16)
L reort ( x ) = μ ν S μ ν ( x ) W μ ν ,
(17)

where S μ ν ( x ) is the differentiated AO overlap matrix and h μ ν ( x ) is a differentiated one-electron integral in the AO basis,

S μ ν ( x ) = μ | ν ( x ) ,
(18)
h μ ν ( x ) = μ | 1 2 2 K Z K r K 1 | ν ( x ) ,
(19)

with ZK and rK representing the nuclear charge and the nuclear-electronic distance for nucleus K, respectively. The G(A) transformation on a general matrix A is given as the sum of a Coulomb contribution J(A) and an exchange contribution K(A),

G μ ν ( A ) = 2 J μ ν ( A ) K μ ν ( A ) ,
(20)
J μ ν ( A ) = α β g μ ν α β A α β ,
(21)
K μ ν ( A ) = α β g μ β α ν A α β ,
(22)

and the corresponding quantities with differentiated two-electron AO integrals (not differentiated matrix A) are defined as

G μ ν ( x ) ( A ) = 2 J μ ν ( x ) ( A ) K μ ν ( x ) ( A ) ,
(23)
J μ ν ( x ) ( A ) = α β g μ ν α β ( x ) A α β ,
(24)
K μ ν ( x ) ( A ) = α β g μ β α ν ( x ) A α β .
(25)

D is the HF density matrix,

D μ ν = i C μ i C ν i ,
(26)

and ρAO is the MP2 correlation density matrix in the AO basis, which is connected to the corresponding MO matrix ρMO via the following transformation:

ρ AO = C ρ MO C T .
(27)

ρMO is partitioned into occupied/virtual blocks,

ρ MO = X i j κ ̄ i a T κ ̄ a i Y a b .
(28)

The locality analysis in Ref. 74 shows that the X intermediate is most efficiently evaluated using the virtual partitioning scheme in the following manner:

X i j = P ( X P ) i j + P > Q ( Δ X PQ ) i j ,
(29)
( X P ) i j = a b P ¯ k [ P ¯ ] t k i b a t ̄ k j b a , ( i j [ P ¯ ] ) ,
(30)
( Δ X PQ ) i j = ( a P ¯ b Q ¯ + a Q ¯ b P ¯ ) k [ P ¯ ] [ Q ¯ ] t k i b a t ̄ k j b a ( i j [ P ¯ ] [ Q ¯ ] ) .
(31)

Note that (XP)ij is zero if i and/or j is outside the [ P ¯ ] space, and similarly for (ΔXPQ)ij. Similarly, using the locality arguments of Ref. 74, the Y intermediate may be evaluated using the occupied partitioning scheme,

Y a b = P ( Y P ) a b + P > Q ( Δ Y PQ ) a b ,
(32)
( Y P ) a b = i j P ¯ c [ P ¯ ] t j i c a t ̄ j i c b ( a b [ P ¯ ] ) ,
(33)
( Δ Y PQ ) a b = ( i P ¯ j Q ¯ + i Q ¯ j P ¯ ) c [ P ¯ ] [ Q ¯ ] t j i c a t ̄ j i c b ( a b [ P ¯ ] [ Q ¯ ] ) .
(34)

The κ ̄ a i multipliers are determined by solving coupled-perturbed HF equations of the form

2 a κ ̄ a l F a d 2 i κ ̄ d i F l i + 2 G d l ( κ ̄ s ) = Φ d l vo Φ l d ov + 2 G d l ( M ) ,
(35)

where κ ̄ s has been symmetrized and transformed to the AO basis,

κ ̄ μ ν s = a i ( κ ̄ a i C μ a C ν i + κ ̄ a i C ν a C μ i )
(36)

and M is determined from the X and Y matrices,

M μ ν = i k C μ i C ν k X i k a c C μ a C ν c Y a c .
(37)

G d l ( κ ̄ s ) and Gdl(M) are determined using Eq. (20) where the first and the second index are transformed to the virtual and the occupied MO basis, respectively, e.g.,

G d l ( κ ̄ s ) = α β C α d C β l G α β ( κ ̄ s ) .
(38)

The Φ matrix is written in block form,

Φ = Φ oo Φ ov Φ vo Φ vv ,
(39)

where o/v denotes occupied/virtual orbital indices, and the different blocks are given by

Φ a k vo = P ( Φ P vo ) a k + P > Q ( Δ Φ PQ vo ) a k ,
(40)
( Φ P vo ) a k = 2 i j P ¯ c [ P ¯ ] t ̄ j i c a g c j k i ( a [ P ¯ ] , k [ P ¯ ] ) ,
(41)
( Δ Φ PQ vo ) a k = 2 ( i P ¯ j Q ¯ + i Q ¯ j P ¯ ) c [ P ¯ ] [ Q ¯ ] t ̄ j i c a g c j k i ( a [ P ¯ ] [ Q ¯ ] , k [ P ¯ ] [ Q ¯ ] ) ,
(42)
Φ i c ov = P ( Φ P ov ) i c + P > Q ( Δ Φ PQ ov ) i c ,
(43)
( Φ P ov ) i c = 2 a b P ¯ k [ P ¯ ] t ̄ k i b a g b k a c ( i [ P ¯ ] , c [ P ¯ ] ) ,
(44)
( Δ Φ PQ ov ) i c = 2 ( a P ¯ b Q ¯ + a Q ¯ b P ¯ ) k [ P ¯ ] [ Q ¯ ] t ̄ k i b a g b k a c ( i [ P ¯ ] [ Q ¯ ] , c [ P ¯ ] [ Q ¯ ] ) ,
(45)
Φ i j oo = P ( Φ P oo ) i j + P > Q ( Δ Φ PQ oo ) i j ,
(46)
( Φ P oo ) i j = 2 a b P ¯ k [ P ¯ ] t ̄ k i b a g b k a j , ( i j [ P ¯ ] ) ,
(47)
( Δ Φ PQ oo ) i j = 2 ( a P ¯ b Q ¯ + a Q ¯ b P ¯ ) k [ P ¯ ] [ Q ¯ ] t ̄ k i b a g b k a j ( i j [ P ¯ ] [ Q ¯ ] ) ,
(48)
Φ a b vv = P ( Φ P vv ) a b + P > Q ( Δ Φ PQ vv ) a b ,
(49)
( Φ P vv ) a b = 2 i j P ¯ c [ P ¯ ] t ̄ j i c a g c j b i , ( a b [ P ¯ ] ) ,
(50)
( Δ Φ PQ vv ) a b = 2 ( i P ¯ j Q ¯ + i Q ¯ j P ¯ ) c [ P ¯ ] [ Q ¯ ] t ̄ j i c a g c j b i ( a b [ P ¯ ] [ Q ¯ ] ) .
(51)

The reorthonormalization contribution in Eq. (17) is written in terms of the intermediate matrix W defined as

W = 2 DFD + 1 2 Φ + ρ AO F D ̄ + D G ( ρ AO ) D ̄ ,
(52)

where all matrices are given in the AO basis, and D ̄ is an effective density matrix, which is obtained by allowing the indices in Eq. (26) to run over occupied as well as virtual orbitals,

D ̄ μ ν = p C μ p C ν p .
(53)

The DEC-MP2 gradient scheme may thus be summarized as follows.

  1. Perform HF calculation and localize occupied and virtual HF orbitals.

  2. For each atomic fragment P, perform the following.

    • Calculate MP2 amplitudes t i j a b for the occupied partitioning scheme ( i j P ¯ and a b [ P ¯ ] ) and the virtual partitioning scheme ( a b P ¯ and i j [ P ¯ ] ).

    • Calculate fragment intermediates L Θ , P ( x ) , XP, YP, and ΦP using Eqs. (12), (30), (33), (44), (41), (47), and (50).

    • Update full intermediates L Θ ( x ) , X, Y, and Φ with calculated fragment intermediates using Eqs. (11), (29), (32), (43), (40), (46), and (49).

  3. Repeat step 2 for the important pair fragments (the number of pair fragments to consider scales linearly with fragment size when pair screening is employed).

  4. Solve the coupled-perturbed HF equations in Eq. (35) to determine the κ ̄ matrix.

  5. Determine ρAO and W using Eqs. (27), (28), and (52).

  6. Calculate the molecular gradient using Eq. (10).

Having summarized the working equations for the DEC-MP2 molecular gradient, we now apply the RI approximation using the symmetric decomposition32,35 to write the gaibj integral as

g a i b j g a i b j = Λ Υ ( a i | Λ ) ( Λ | Υ ) 1 ( Υ | b j ) = Λ Υ Φ ( a i | Λ ) ( Λ | Φ ) 1 / 2 ( Φ | Υ ) 1 / 2 ( Υ | b j ) = Φ C a i Φ C b j Φ ,
(54)

where g denotes the RI approximated two-electron repulsion integrals (ERIs), (ai|Λ) is a 3-center ERI, (Λ|Φ) is a 2-center ERI, C a i Φ = Λ ( a i | Λ ) ( Λ | Φ ) 1 / 2 defines the symmetrically decomposed fitting coefficients, and capital greek letters {Υ, Λ, Φ} refer to auxiliary AO indices.

Calculating the gaibj integrals in Section II A according to Eq. (54) with auxiliary AOs assigned to atomic sites in the AE (Υ ∈ {P} for atomic fragment P) defines the DEC-RI-MP2 model for the evaluation of the correlation energy.80 The current DEC-RI-MP2 implementation does not use RI approximations for the initial HF calculation, and the Fock matrix has thus been calculated without employing the RI approximation. Consequently, in the MP2 amplitude equations in Eq. (2), the RI approximation in applied for the gaibj integrals, but not for the Fock matrix elements.

For the gradient evaluation it is convenient to use the non-symmetric formulation where integrals are expressed in the following manner:

g a i b j g a i b j = Υ ( a i | Υ ) D b j Υ ,
(55)

with the non-symmetric fitting coefficients D a i Υ = Λ ( Υ | Λ ) 1 ( Λ | a i ) . The formulation using the symmetrically decomposed fitting coefficients is convenient for the energy evaluation as the storage requirements are reduced by a factor of 2, since the storage of the 3-center ERIs (ai|Υ) is not required. However, in the case of the integral derivatives, the symmetric formulation does not lead to reduced memory requirements, and the formulation using the non-symmetric fitting coefficients leads to simpler expressions.

As described above, the initial HF calculation does not employ the RI approximation, while the correlated part of the calculation does. For this reason, only some of the integrals of the DEC-RI-MP2 gradient are evaluated using the RI approximation. In summary, the DEC-RI-MP2 molecular gradient may be evaluated using the procedure summarized at the end of Section II B with the following modifications for each fragment calculation.

  • Determine RI-MP2 amplitudes t i j a b (and hence the associated multipliers t ̄ i j a b = 4 t i j a b 2 t i j b a ) by solving the amplitude equation in Eq. (2) using RI-approximated ERIs80 calculated using the symmetric formulation (Eq. (54)).

  • Evaluate all ERIs entering the Φ matrix in Eqs. (40) to (51) using the symmetric RI formulation.

  • Replace the derivative ERIs g a i b j ( x ) by g a i b j ( x ) calculated using the non-symmetric RI formulation, i.e.,
    g a i b j ( x ) = Υ ( a i | Υ ) ( x ) D b j Υ Λ Υ ( a i | Λ ) ( Λ | Υ ) ( x ) D b j Υ + Υ D a i Υ ( Υ | b j ) ( x ) .
    (56)

The last point results in the modified L Θ ( x ) terms

L Θ ( x ) = i j a b t ̄ i j a b g a i b j ( x ) P L Θ , P ( x ) + P > Q Δ L Θ , PQ ( x ) ,
(57)
L Θ , P ( x ) = i j P ¯ a b [ P ¯ ] t ̄ i j a b g a i b j ( x ) ,
(58)
Δ L Θ , PQ ( x ) = ( i P ¯ j Q ¯ + i Q ¯ j P ¯ ) a b [ P ¯ ] [ Q ¯ ] t ̄ i j a b g a i b j ( x ) .
(59)

Note that the derivative integrals in Eq. (23) are not to be replaced according to Eq. (56) as these derivative integrals originate from differentiation of Fock matrix elements, which have been calculated without employing the RI approximation (since the HF calculation does not employ the RI approximation).

Before discussing the computational results we review a few aspects of the implementation with emphasis on the RI modified L Θ ( x ) fragment terms in Eq. (57). Exploiting permutational symmetry and rearranging the terms L Θ ( x ) (without involving any approximations) may be written as

L Θ ( x ) = 2 i j a b t ̄ i j a b Υ ( a i | Υ ) ( x ) D b j Υ i j a b t ̄ i j a b Λ Υ ( a i | Λ ) ( Λ | Υ ) ( x ) D b j Υ
(60)
= 2 i a Υ ( a i | Υ ) ( x ) ( j b t ̄ i j a b D b j Υ ) Λ Υ ( Λ | Υ ) ( x ) ( i j a b t ̄ i j a b ( a i | Λ ) D b j Υ ) .
(61)

Defining the two intermediates

Θ i a Υ = 2 j b t ̄ i j a b D b j Υ ,
(62)
G Λ Υ = 1 2 i a ( a i | Λ ) Θ i a Υ ,
(63)

we may write L Θ ( x ) in the simple form

L Θ ( x ) = i a Υ ( a i | Υ ) ( x ) Θ i a Υ Λ Υ ( Λ | Υ ) ( x ) G Λ Υ .
(64)

A similar expression may be obtained for the atomic fragments and pair fragment contributions in Eqs. (58) and (59),

L Θ , P ( x ) = i P ¯ a [ P ¯ ] Υ { P } ( a i | Υ ) ( x ) Θ i a Υ P Λ Υ { P } ( Λ | Υ ) ( x ) G Λ Υ P ,
(65)
Δ L Θ , PQ ( x ) = k P Q a [ P ¯ ] [ Q ¯ ] Υ { P } { Q } ( a k | Υ ) ( x ) Θ k a Υ P Q Λ Υ { P } { Q } ( Λ | Υ ) ( x ) G Λ Υ P Q ,
(66)

where we have introduced,

Θ i a Υ P = 2 j P b [ P ¯ ] t ̄ i j a b D b j Υ ,
(67)
G Λ Υ P = 1 2 i P a [ P ¯ ] ( a i | Λ ) Θ i a Υ P ,
(68)
Θ k a Υ P Q = { 2 j Q b [ P ¯ ] [ Q ¯ ] t ̄ i j a b D b j Υ ( k P ) 2 i P b [ P ¯ ] [ Q ¯ ] t ̄ i j a b D b j Υ ( k Q ) ,
(69)
G Λ Υ P Q = 1 2 a [ P ¯ ] [ Q ¯ ] k P Q ( a k | Λ ) Θ k a Υ P Q .
(70)

Table I shows the storage requirements and the computational cost for constructing the intermediates of the L Θ , P ( x ) contribution to the DEC-RI-MP2 gradient calculation. The computational cost is determined as follows:

  • The occupied MOs assigned to P ( P ¯ ) denoted the occupied energy orbital space (EOS) of dimension OEOS.

  • The occupied/virtual MOs that are local to P ( [ P ¯ ] / [ P ¯ ] ) denoted the amplitude orbital space (AOS) of dimension OAOS/VAOS.

  • The number of standard AOs NAO in the AE ({P}).

  • The number of auxiliary AOs Naux in the AE.

  • The number of atoms Natoms in the AE.

Note that these dimensions all refer to an atomic fragment P, not the full molecular system. For pair fragments we use unions of these spaces for the atoms involved (e.g., P and Q).

TABLE I.

Storage requirements and the computational cost for the dominant step for constructing the intermediates of the L Θ , P ( x ) contribution to the DEC-RI-MP2 gradient calculation.

Intermediate Storage Computational cost
t i j a b   O EOS 2 V AOS 2   N aux O AOS 2 V AOS 2  
t ̄ i j a b   O EOS 2 V AOS 2   O EOS 2 V AOS 2  
Θ i a Υ P   NauxOEOSVAOS   N aux O EOS 2 V AOS 2  
G Λ Υ P   N aux 2   N aux 2 O EOS V AOS  
L Θ , P ( x )   3Natoms  3 N atoms N aux O EOS V AOS + 3 N atoms N aux 2  
(Λ|Υ)(x)  3 N atoms N aux 2   3 N atoms N aux 2  
(ak|Υ)(x)  3NatomsNauxOEOSVAOS  3 N atoms N aux N AO 2 O EOS  
Intermediate Storage Computational cost
t i j a b   O EOS 2 V AOS 2   N aux O AOS 2 V AOS 2  
t ̄ i j a b   O EOS 2 V AOS 2   O EOS 2 V AOS 2  
Θ i a Υ P   NauxOEOSVAOS   N aux O EOS 2 V AOS 2  
G Λ Υ P   N aux 2   N aux 2 O EOS V AOS  
L Θ , P ( x )   3Natoms  3 N atoms N aux O EOS V AOS + 3 N atoms N aux 2  
(Λ|Υ)(x)  3 N atoms N aux 2   3 N atoms N aux 2  
(ak|Υ)(x)  3NatomsNauxOEOSVAOS  3 N atoms N aux N AO 2 O EOS  

Naturally the storage in memory of (ak|Υ)(x) and (Λ|Υ)(x) quickly becomes infeasible. However, using a batching of the auxiliary functions (and the AOs if needed) these integrals may be constructed in parts and immediately contracted with Θ k a Υ P and G Λ Υ P . The time critical part of the DEC-RI-MP2 molecular gradient scheme is the evaluation of the (ak|Υ)(x) derivative integral which depends critically on the size of the atomic extent of the fragment Natoms.

All DEC-RI-MP2 calculations were done using a development version of the LSDalton electronic structure program, which is a part of Dalton2016 suite.81,82 The reference numbers were obtained with ORCA suite of programs version 3.0.83,84 Correlation consistent basis set of double-ζ quality (cc-pVDZ) was used for all calculations,85 and the corresponding matching cc-pVDZ-RI basis set was used for the RI terms. The self-consistent-field equations were converged using the “Tight” keywords for both programs, which corresponds to the energy change of 10−8 hartree. The coupled-perturbed HF equations Eq. (35) converged to reach a residual norm of 10−8. The fragment optimization threshold (FOT) in the DEC calculations was chosen based on our previous investigation74 to target an accuracy of 10−3 hartree/angstrom in the final molecular gradient.

We have employed the S12L test set of Grimme et al.86 to estimate the errors of the DEC approximation. The test set consists of supramolecular complexes of up to 158 atoms (the cucurbit[7]uril (CB7) host inclusion complex with the bis(trimethylammoniomethyl)ferrocene (FECP) was excluded as it has an open-shell metal core). The test set contains molecules of complicated electronic structures and is thus deliberately chosen to demonstrate the general applicability of the DEC method.

The individual complexes of the test set are as follows (see Figure 1): two “tweezer” complexes with tetracyanoquinone (TCNQ) and 1,4-dicyanobenzene—TCNA@tweezer and DCB@tweezer; two “pincer” complexes of organic π-systems—3c@pincer and 3d@pincer; fullerenes C60 and C70 in a so-called “buckycatcher”—C60@catcher and C70@catcher; an amide macrocycle (mcycle) that binds benzoquinone (BQ) and glycine anhydride (GLH)—BQ@mcycle and GLH@mcycle; the cucurbit[6]uril (CB6) cation complexes of butylammonium (BuNH4) and propylammonium (PrNH4)—BuNH4@CB6 and PrNH4@CB6; and the cucurbit[7]uril (CB7) host with 1-hydroxyadamantane (ADOH)—ADOH@CB7. For the insulin molecule considered in Section III C, the molecular geometry was obtained as described in Ref. 87 where the original geometry was taken from Ref. 88. All molecular geometries used in this work are listed in the supplementary material.89 

FIG. 1.

The reduced S12L test set used in this study.

FIG. 1.

The reduced S12L test set used in this study.

Close modal

All DEC-RI-MP2 gradient calculations were run on the Titan supercomputer, a Cray XK7 system deployed by the US Department of Energy at the Oak Ridge National Laboratory. Each Titan node consists of a 16-core AMD Opteron 6274 processor with 32 GB of main memory. Titan’s nodes are connected according to a 3D torus topology with the Cray Gemini network interfaces. The application-level communication layer is provided by the Cray MPI library implementing most of the MPI 3.0 features, a derivative of MPICH-2 from the Argonne National Laboratory.

The DEC errors were estimated by comparing the results of conventional RI-MP2 calculations with those obtained using the DEC-RI-MP2 method. Direct comparison with the conventional implementation is the most efficient and simple way to conduct the error analysis provided that the reference analytical molecular gradient calculation can be carried out. The DEC method is aimed at calculations on very large molecules, and we therefore consider a test set of rather large molecular systems (Figure 1). In order to enable the determination of reference numbers for this test set, we have restricted ourselves to the correlation consistent basis set of double-ζ quality. Obtaining reference gradients without the RI approximation to test the RI performance benefits and errors is impossible for the considered test set. Moreover, the RI acceleration and errors are now very well understood and documented in the literature.35–39 In the DEC context, introducing the RI approximation means that each fragment calculation benefits from RI approximation in the very same way as any conventional RI-MP2 calculation. Thus, the RI performance benefits and errors are out of scope for the paper and the presented error analysis reflects solely errors coming from DEC approximation.

Tables II and III summarize the results of the test set calculations and associated DEC errors. As it is evident from the tables, the absolute maximum DEC errors are of the order 10−3 hartree/angstrom. This is in accordance with the results obtained previously74 for the DEC-MP2 gradient for smaller molecular systems. Also, the DEC-RI-MP2 gradient errors are in general independent of the system size or type of electronic structure, even for the highly delocalized C60@catcher and C70@catcher. The full molecular gradients are given in the supplementary material89 for future comparison with other methods.

TABLE II.

Summary of the DEC-RI-MP2 test set calculation for the S12L test set. The number of regular AOs and auxiliary AOs is given as well as the HF, correlation, total energies, and the gradient norms.

Complex Basis/AUX basis HF Correlation Total energy Gradient norm
TCNA@tweezer  982/3724  −2282.511 183 59  −7.7317  −2290.2429  0.1542 
DCB@tweezer  898/3388  −2022.198 430 64  −6.9032  −2029.1016  0.1174 
3c@pincer  1341/5082  −3804.668 221 82  −11.7084  −3816.3906  0.0802 
3d@pincer  1190/4500  −3697.985 014 30  −10.1780  −3708.1630  0.0976 
C60@catcher  1820/7112  −4560.549 432 22  −15.5399  −4576.0894  0.1226 
C70@catcher  1960/7672  −4939.295 517 56  −16.8315  −4956.1270  0.1147 
BQ@mcycle  1404/5208  −3272.631 842 13  −10.7488  −3283.3806  0.0709 
GLH@mcycle  1394/5180  −3238.224 661 99  −10.6522  −3248.8769  0.0743 
BuNH4@CB6  1318/4984  −3802.794 931 92  −11.5053  −3814.3002  0.0676 
PrNH4@CB6  1294/4900  −3763.755 541 33  −11.3561  −3775.1117  0.0669 
ADOH@CB7  1620/6132  −4651.225 865 58  −14.1529  −4665.3788  0.0734 
Complex Basis/AUX basis HF Correlation Total energy Gradient norm
TCNA@tweezer  982/3724  −2282.511 183 59  −7.7317  −2290.2429  0.1542 
DCB@tweezer  898/3388  −2022.198 430 64  −6.9032  −2029.1016  0.1174 
3c@pincer  1341/5082  −3804.668 221 82  −11.7084  −3816.3906  0.0802 
3d@pincer  1190/4500  −3697.985 014 30  −10.1780  −3708.1630  0.0976 
C60@catcher  1820/7112  −4560.549 432 22  −15.5399  −4576.0894  0.1226 
C70@catcher  1960/7672  −4939.295 517 56  −16.8315  −4956.1270  0.1147 
BQ@mcycle  1404/5208  −3272.631 842 13  −10.7488  −3283.3806  0.0709 
GLH@mcycle  1394/5180  −3238.224 661 99  −10.6522  −3248.8769  0.0743 
BuNH4@CB6  1318/4984  −3802.794 931 92  −11.5053  −3814.3002  0.0676 
PrNH4@CB6  1294/4900  −3763.755 541 33  −11.3561  −3775.1117  0.0669 
ADOH@CB7  1620/6132  −4651.225 865 58  −14.1529  −4665.3788  0.0734 
TABLE III.

Test set error analysis results from comparing the DEC-RI-MP2 gradient with the conventional RI-MP2 gradient. The ABS stands for absolute maximum deviation, MAD—mean absolute deviation, and RMSD—root mean square deviation.

Complex ABS × 10−3 MAD × 10−3 RMSD × 10−3
TCNA@tweezer  4.4  0.5  0.9 
DCB@tweezer  3.0  0.4  0.6 
3c@pincer  4.0  0.4  0.6 
3d@pincer  3.9  0.4  0.7 
C60@catcher  4.7  0.8  1.2 
C70@catcher  6.4  1.3  2.0 
BQ@mcycle  1.8  0.3  0.4 
GLH@mcycle  1.9  0.3  0.4 
BuNH4@CB6  1.2  0.2  0.3 
PrNH4@CB6  1.3  0.2  0.3 
ADOH@CB7  1.6  0.3  0.4 
Complex ABS × 10−3 MAD × 10−3 RMSD × 10−3
TCNA@tweezer  4.4  0.5  0.9 
DCB@tweezer  3.0  0.4  0.6 
3c@pincer  4.0  0.4  0.6 
3d@pincer  3.9  0.4  0.7 
C60@catcher  4.7  0.8  1.2 
C70@catcher  6.4  1.3  2.0 
BQ@mcycle  1.8  0.3  0.4 
GLH@mcycle  1.9  0.3  0.4 
BuNH4@CB6  1.2  0.2  0.3 
PrNH4@CB6  1.3  0.2  0.3 
ADOH@CB7  1.6  0.3  0.4 

As an example of a computationally demanding system, which is not accessible with a conventional implementation, we have chosen insulin (see Figure 2). This molecule contains 787 atoms and 7604 regular and 28 148 auxiliary basis functions. All computational thresholds were kept at the same level as for the test set examples. The DEC-RI-MP2 gradient was obtained in less than 10 h using one third of the Titan supercomputer (6000 nodes). The results of the insulin calculation are presented in Table IV.

TABLE IV.

Summary of the DEC-RI-MP2 insulin calculation.

Hartree-Fock energy  −21 645.424 375 24 
Correlation energy  −60.610 9 
Total RI-MP2 energy  −21 706.035 3 
DEC-RI-MP2 GRADIENT NORM  0.775 0 
DEC-RI-MP2 GRADIENT RMS  0.015 9 
Hartree-Fock energy  −21 645.424 375 24 
Correlation energy  −60.610 9 
Total RI-MP2 energy  −21 706.035 3 
DEC-RI-MP2 GRADIENT NORM  0.775 0 
DEC-RI-MP2 GRADIENT RMS  0.015 9 
FIG. 2.

The insulin molecule.

FIG. 2.

The insulin molecule.

Close modal

A direct comparison with a conventional implementation is impossible, but we might consider the trace of the correlation density matrix divided by the number of electrons, Tr(ρMO)/Nel, as a size-intensive error measure.74 This measure is zero if no DEC approximation is invoked (FOT = 0), while it is nonzero for practical FOT values. For the DEC-MP2 model we have previously found Tr(ρMO)/Nel to be −7.4 × 10−6 for insulin using the same geometry, basis, and FOT value as employed here. The corresponding number for DEC-RI-MP2 is 9.6 × 10−6. As expected, these numbers are thus of the same order of magnitude, since the same FOT has been used.

Finally, we note that the DEC-RI-MP2 molecular gradient for insulin is provided in the supplementary material89 to enable future comparisons to lower-level methods, such as the atomic forces determined from classical force field calculations. In general, the DEC-RI-MP2 gradient scheme might provide useful benchmark data for large molecular systems for semi-empirical methods which are often parameterized based on CC calculations for small molecular systems.

We have presented the DEC-RI-MP2 molecular gradient scheme, which combines the precision control as well as the linear-scaling and massively parallel features of the DEC scheme with efficient evaluations of the gradient contributions using the RI approximation. The DEC-RI-MP2 molecular gradient was successfully applied to a test set of large molecular systems, including complicated electronic structures, to demonstrate the general and wide applicability of the method. We have furthermore performed a large-scale DEC-RI-MP2 gradient calculation on the insulin molecule, which contains 787 atoms and 7604 regular and 28 148 auxiliary basis functions. This calculation was completed in less than 10 h using one third of the Titan supercomputer.

The overall goal of the DEC scheme is to provide a framework for treating large molecular systems at the CC level of theory, i.e., systems where the conventional implementation hits a scaling wall. Furthermore, the DEC scheme is linear-scaling with system size and designed to utilize thousands of computing cores on a large high performance computer (HPC) architecture. However, the DEC method is rather inefficient for small molecular systems. As a rule of thumb, if the canonical calculation can be performed within a reasonable time frame, then the canonical calculation will use fewer computational resources than the corresponding DEC calculation. However, if one is willing to use many compute nodes, then the time-to-solution provided by the DEC method will generally be much shorter than for a canonical implementation. Nonetheless, the main advantage of the DEC method is the ability to treat very large molecular systems where canonical implementations encounter a scaling wall. Such calculations are computationally demanding, but linear-scaling and massively parallel. The insulin calculation presented in this work represents such a case, and it was carried out in less than 10 h. It is one of the largest correlated gradient calculations known to date.

In many aspects the comparison of DEC and canonical algorithms mirrors the rise of direct SCF, which was initially only an advantage when disk space was insufficient. However, the development of linear-scaling algorithms, improved integral screening, and density-fitting technology made direct SCF competitive, even in cases of sufficient disk space. Direct SCF has also benefitted from HPC developments which favor computations rather than memory access. Similar developments might happen for the DEC scheme due the continued algorithmic and HPC architecture improvements.

Finally, an additional obvious advantage that should be emphasized is the relatively low effort required for extending the DEC scheme to the determination of molecular gradients and other molecular properties, while the linear-scaling and massively parallel features of the DEC scheme are retained.

This research used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the Department of Energy under Contract No. DE-AC05-00OR22725. The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme No. (FP/2007-2013)/ERC Grant Agreement No. 291371. D.B. acknowledges the Marie Curie Individual Fellowship funding for “DECOS,” Project No. 657514.

1.
H. F.
Schaefer
and
Y.
Yamaguchi
,
J. Mol. Struct.: THEOCHEM
28
,
369
(
1986
).
2.
P.
Jørgensen
and
J.
Simons
,
Geometrical Derivatives of Energy Surfaces and Molecular Properties
(
D. Reidel Publishing Company
,
Holland
,
1986
).
3.
S.
Bratoz
,
R.
Daudel
,
M.
Roux
, and
M.
Allavena
,
J. Mol. Struct.: THEOCHEM
32
,
412
(
1960
).
4.
P.
Pulay
,
Mol. Phys.
17
,
197
(
1969
).
5.
J. A.
Pople
,
R.
Krishnan
,
H. B.
Schlegel
, and
J. S.
Binkley
,
Int. J. Quantum Chem. Symp.
13
,
225
(
1979
).
6.
G.
Fitzgerald
,
R.
Harrison
,
W. D.
Laidig
, and
R. J.
Bartlett
,
J. Chem. Phys.
82
,
4379
(
1985
).
7.
J.
Gauss
and
D.
Cremer
,
Chem. Phys. Lett.
153
,
303
(
1988
).
8.
B. R.
Brooks
,
W. D.
Laidig
,
P.
Saxe
,
J. D.
Goddard
,
Y.
Yamaguchi
, and
H. F.
Schaefer
III
,
J. Chem. Phys.
72
,
4652
(
1980
).
9.
R.
Krishnan
,
H.
Schlegel
, and
J. A.
Pople
,
J. Chem. Phys.
72
,
4654
(
1980
).
10.
J.
Gauss
and
D.
Cremer
,
Chem. Phys. Lett.
150
,
280
(
1988
).
11.
J.
Gauss
and
D.
Cremer
,
Chem. Phys. Lett.
163
,
549
(
1989
).
12.
A.
Scheiner
,
G. E.
Scuseria
,
J. E.
Rice
,
T. J.
Lee
, and
H. F.
Schaefer
III
,
J. Chem. Phys.
87
,
5361
(
1987
).
13.
G. E.
Scuseria
,
J. Chem. Phys.
94
,
442
(
1991
).
14.
T. J.
Lee
and
A. P.
Rendell
,
J. Chem. Phys.
94
,
6229
(
1991
).
15.
J.
Gauss
and
J. F.
Stanton
,
Phys. Chem. Chem. Phys.
2
,
2047
(
2000
).
16.
N. C.
Handy
,
R. D.
Amos
,
J. F.
Gaw
,
J. E.
Rice
, and
E. D.
Simandiras
,
Chem. Phys. Lett.
120
,
151
(
1985
).
17.
E.
Rice
and
R. D.
Amos
,
Chem. Phys. Lett.
122
,
585
(
1985
).
18.
E. D.
Simandiras
,
R. D.
Amos
, and
N. C.
Handy
,
Chem. Phys.
114
,
9
(
1987
).
19.
P.
Jørgensen
and
T.
Helgaker
,
J. Chem. Phys.
89
,
1560
(
1988
).
20.
T.
Helgaker
,
P.
Jørgensen
, and
N. C.
Handy
,
Theor. Chim. Acta
76
,
227
(
1989
).
21.
M. J.
Frisch
,
M.
Head-Gordon
, and
J. A.
Pople
,
Chem. Phys. Lett.
166
,
275
(
1990
).
22.
M. J.
Frisch
,
M.
Head-Gordon
, and
J. A.
Pople
,
Chem. Phys. Lett.
166
,
281
(
1990
).
23.
F.
Haase
and
R.
Ahlrichs
,
J. Comput. Chem.
14
,
907
(
1983
).
24.
I. M. B.
Nielsen
,
Chem. Phys. Lett.
255
,
210
(
1996
).
25.
G. D.
Fletcher
,
A. P.
Rendell
, and
P.
Sherwood
,
Mol. Phys.
91
,
431
(
1997
).
26.
M.
Head-Gordon
,
Mol. Phys.
96
,
673
(
1999
).
27.
C. M.
Aikens
,
S. P.
Webb
,
R. L.
Bell
,
G. D.
Fletcher
,
M. W.
Schmidt
, and
M. S.
Gordon
,
Theor. Chim. Acta
110
,
233
(
2003
).
28.
K.
Ishimura
,
P.
Pulay
, and
S.
Nagase
,
J. Comput. Chem.
28
,
2034
(
2007
).
29.
E. J.
Baerends
,
D. E.
Ellis
, and
P.
Ros
,
Chem. Phys.
2
,
41
(
1973
).
30.
J. W.
Mintmire
and
B. I.
Dunlap
,
Phys. Rev. A
25
,
88
(
1982
).
31.
B. I.
Dunlap
,
J. W. D.
Connolly
, and
J. R.
Sabin
,
J. Chem. Phys.
71
,
3396
(
1979
).
32.
O.
Vahtras
,
J.
Almlöf
, and
M.
Feyereisen
,
Chem. Phys. Lett.
213
,
514
(
1993
).
33.
K.
Eichkorn
,
O.
Treutler
,
H.
Ohm
,
M.
Haser
, and
R.
Ahlrichs
,
Chem. Phys. Lett.
240
,
283
(
1995
).
34.
K.
Eichkorn
,
F.
Weigend
,
O.
Treutler
, and
R.
Ahlrichs
,
Theor. Chem. Acc.
97
,
119
(
1997
).
35.
F.
Weigend
and
M.
Häser
,
Theor. Chim. Acta
97
,
331
(
1997
).
36.
Y. M.
Rhee
,
R. A.
DiStasio
,
R. C.
Lochan
, and
M.
Head-Gordon
,
Chem. Phys. Lett.
426
,
197
(
2006
).
37.
R. A.
DiStasio
,
R. P.
Steele
,
Y. M.
Rhee
,
Y.
Shao
, and
M.
Head-Gordon
,
J. Comput. Chem.
28
,
839
(
2007
).
38.
S.
Kossmann
and
F.
Neese
,
J. Chem. Theory Comput.
6
,
2325
(
2010
).
39.
U.
Bozkaya
and
C. D.
Sherrill
,
J. Chem. Phys.
144
,
174103
(
2016
).
40.
P.
Deglmann
,
K.
May
,
F.
Furche
, and
R.
Ahlrichs
,
Chem. Phys. Lett.
384
,
103
(
2004
).
41.
D.
Bykov
,
T.
Petrenko
,
R.
Izsak
,
S.
Kossmann
,
U.
Becker
,
E.
Valeev
, and
F.
Neese
,
Mol. Phys.
113
,
1961
(
2015
).
42.
D. H.
Friese
,
C.
Hattig
, and
J.
Kossmann
,
J. Chem. Theory Comput.
9
,
1469
(
2013
).
43.
P.
Pulay
,
Chem. Phys. Lett.
100
,
151
(
1983
).
44.
S.
Saebø
and
P.
Pulay
,
Annu. Rev. Phys. Chem.
44
,
213
(
1993
).
45.
C.
Hampel
and
H.-J.
Werner
,
J. Chem. Phys.
104
,
6286
(
1996
).
46.
M.
Schütz
,
G.
Hetzer
, and
H.-J.
Werner
,
J. Chem. Phys.
111
,
5691
(
1999
).
47.
M.
Schütz
and
H.-J.
Werner
,
J. Chem. Phys.
114
,
661
(
2001
).
48.
H.-J.
Werner
and
M.
Schütz
,
J. Chem. Phys.
135
,
144116
(
2011
).
49.
G. E.
Scuseria
and
P. Y.
Ayala
,
J. Chem. Phys.
111
,
8330
(
1999
).
50.
O.
Christiansen
,
P.
Manninen
,
P.
Jørgensen
, and
J.
Olsen
,
J. Chem. Phys.
124
,
084103
(
2006
).
51.
V.
Weijo
,
P.
Manninen
,
P.
Jørgensen
,
O.
Christiansen
, and
J.
Olsen
,
J. Chem. Phys.
127
,
074106
(
2007
).
52.
N.
Flocke
and
R. J.
Bartlett
,
J. Chem. Phys.
121
,
10935
(
2004
).
53.
S.
Li
,
J.
Ma
, and
Y.
Jiang
,
J. Comput. Chem.
23
,
237
(
2002
).
54.
W.
Li
,
P.
Piecuch
,
J. R.
Gour
, and
S.
Li
,
J. Chem. Phys.
131
,
114109
(
2009
).
55.
W.
Li
and
P.
Piecuch
,
J. Phys. Chem. A
114
,
8644
(
2010
).
56.
J. E.
Subotnik
,
A.
Sodt
, and
M.
Head-Gordon
,
J. Chem. Phys.
125
,
074116
(
2006
).
57.
M.
Kobayashi
and
H.
Nakai
,
J. Chem. Phys.
129
,
044103
(
2008
).
58.
D. G.
Fedorov
and
K.
Kitaura
,
J. Chem. Phys.
123
,
134103
(
2005
).
59.
H.
Stoll
,
Chem. Phys. Lett.
191
,
548
(
1992
).
60.
J.
Friedrich
,
M.
Hanrath
, and
M.
Dolg
,
J. Chem. Phys.
126
,
154110
(
2007
).
61.
M.
Häser
,
Theor. Chim. Acta
87
,
147
(
1993
).
62.
P. Y.
Ayala
and
G. E.
Scuseria
,
J. Chem. Phys.
110
,
3660
(
1999
).
63.
D. S.
Lambrecht
,
B.
Doser
, and
C.
Ochsenfeld
,
J. Chem. Phys.
123
,
184102
(
2005
).
64.
B.
Doser
,
D. S.
Lambrecht
, and
C.
Ochsenfeld
,
Phys. Chem. Chem. Phys.
10
,
3335
(
2008
).
65.
A. A.
Auer
and
M.
Nooijen
,
J. Chem. Phys.
125
,
024104
(
2006
).
66.
F.
Neese
,
F.
Wennmohs
, and
A.
Hansen
,
J. Chem. Phys.
130
,
114108
(
2009
).
67.
F.
Neese
,
A.
Hansen
, and
D. G.
Liakos
,
J. Chem. Phys.
131
,
064103
(
2009
).
68.
J.
Yang
,
Y.
Kurashige
,
F. R.
Manby
, and
G. K. L.
Chan
,
J. Chem. Phys.
134
,
044123
(
2011
).
69.
J.
Yang
,
G. K. L.
Chan
,
F. R.
Manby
,
M.
Schütz
, and
H.-J.
Werner
,
J. Chem. Phys.
136
,
114105
(
2012
).
70.
M.
Ziółkowski
,
B.
Jansík
,
T.
Kjærgaard
, and
P.
Jørgensen
,
J. Chem. Phys.
133
,
014107
(
2010
).
71.
K.
Kristensen
,
M.
Ziółkowski
,
B.
Jansík
,
T.
Kjærgaard
, and
P.
Jørgensen
,
J. Chem. Theory Comput.
7
,
1677
(
2011
).
72.
I.-M.
Høyvik
,
K.
Kristensen
,
B.
Jansík
, and
P.
Jørgensen
,
J. Chem. Phys.
136
,
014105
(
2012
).
73.
P.
Ettenhuber
,
P.
Baudin
,
T.
Kjærgaard
,
P.
Jørgensen
, and
K.
Kristensen
,
J. Chem. Phys.
144
,
164116
(
2016
).
74.
K.
Kristensen
,
P.
Jørgensen
,
B.
Jansík
,
T.
Kjærgaard
, and
S.
Reine
,
J. Chem. Phys.
137
,
114102
(
2012
).
75.
A. E.
Azhary
,
G.
Rauhut
,
P.
Pulay
, and
H.-J.
Werner
,
J. Chem. Phys.
108
,
5185
(
1998
).
76.
M.
Schütz
,
H.-J.
Werner
,
R.
Lindh
, and
F. R.
Manby
,
J. Chem. Phys.
121
,
737
(
2004
).
77.
T.
Nagata
,
D. G.
Fedorov
,
K.
Ishimura
, and
K.
Kitaura
,
J. Chem. Phys.
135
,
044110
(
2011
).
78.
S.
Schweizer
,
B.
Doser
, and
C.
Ochsenfeld
,
J. Chem. Phys.
128
,
154101
(
2008
).
79.
J.
Almlöf
,
Chem. Phys. Lett.
181
,
319
(
1991
).
80.
P.
Baudin
,
P.
Ettenhuber
,
S.
Reine
,
K.
Kristensen
, and
T.
Kjærgaard
,
J. Chem. Phys.
144
,
054102
(
2016
).
81.
K.
Aidas
,
C.
Angeli
,
K. L.
Bak
,
V.
Bakken
,
R.
Bast
,
L.
Boman
,
O.
Christiansen
,
R.
Cimiraglia
,
S.
Coriani
,
P.
Dahle
,
E. K.
Dalskov
,
U.
Ekström
,
T.
Enevoldsen
,
J. J.
Eriksen
,
P.
Ettenhuber
,
B.
Fernández
,
L.
Ferrighi
,
H.
Fliegl
,
L.
Frediani
,
K.
Hald
,
A.
Halkier
,
C.
Hättig
,
H.
Heiberg
,
T.
Helgaker
,
A. C.
Hennum
,
H.
Hettema
,
E.
Hjertenæs
,
S.
Høst
,
I.-M.
Høyvik
,
M. F.
Iozzi
,
B.
Jansik
,
H. J. A.
Jensen
,
D.
Jonsson
,
P.
Jørgensen
,
J.
Kauczor
,
S.
Kirpekar
,
T.
Kjærgaard
,
W.
Klopper
,
S.
Knecht
,
R.
Kobayashi
,
H.
Koch
,
J.
Kongsted
,
A.
Krapp
,
K.
Kristensen
,
A.
Ligabue
,
O. B.
Lutnæs
,
J. I.
Melo
,
K. V.
Mikkelsen
,
R. H.
Myhre
,
C.
Neiss
,
C. B.
Nielsen
,
P.
Norman
,
J.
Olsen
,
J. M. H.
Olsen
,
A.
Osted
,
M. J.
Packer
,
F.
Pawlowski
,
T. B.
Pedersen
,
P. F.
Provasi
,
S.
Reine
,
Z.
Rinkevicius
,
T. A.
Ruden
,
K.
Ruud
,
V.
Rybkin
,
P.
Salek
,
C. C. M.
Samson
,
A. S.
de Merás
,
T.
Saue
,
S. P. A.
Sauer
,
B.
Schimmelpfennig
,
K.
Sneskov
,
A. H.
Steindal
,
K. O.
Sylvester-Hvid
,
P. R.
Taylor
,
A. M.
Teale
,
E. I.
Tellgren
,
D. P.
Tew
,
A. J.
Thorvaldsen
,
L.
Thøgersen
,
O.
Vahtras
,
M. A.
Watson
,
D. J. D.
Wilson
,
M.
Ziolkowski
, and
H.
Ågren
,
WIREs: Comput. Mol. Sci.
4
,
269
(
2013
).
82.
LSDalton, a linear-scaling molecular electronic structure program, Release Dalton2016, 2015, see http://daltonprogram.org.
83.
F.
Neese
,
F.
Wennmohs
,
U.
Becker
,
D.
Bykov
,
D.
Ganyushin
,
A.
Hansen
,
R.
Izsak
,
D.
Liakos
,
C.
Kollmar
,
S.
Kossmann
,
D. A.
Pantazis
,
T.
Petrenko
,
C.
Reimann
,
C.
Riplinger
,
M.
Roemelt
,
B.
Sandhofer
,
I.
Schapiro
,
K.
Sivalingam
, and
B.
Wezisla
, MPI CEC, Mlheim a.d. Ruhr, Germany (2015), ORCA, version 3.0.
84.
F.
Neese
,
WIREs: Comput. Mol. Sci.
2
,
73
(
2012
).
85.
T. H.
Dunning
, Jr.
,
J. Chem. Phys.
90
,
1007
(
1989
).
86.
S.
Grimme
,
Chem.–A Eur. J.
18
,
9955
(
2012
).
87.
K.
Kristensen
,
I.-M.
Høyvik
,
B.
Jansík
,
P.
Jørgensen
,
T.
Kjærgaard
,
S.
Reine
, and
J.
Jakowski
,
Phys. Chem. Chem. Phys.
14
,
15706
(
2012
).
88.
G. D.
Smith
,
W. A.
Pangborn
, and
R. H.
Blessing
,
Acta Crystallogr., Sect. D: Biol. Crystallogr.
59
,
474
(
2003
).
89.
See supplementary material at http://dx.doi.org/10.1063/1.4956454 for molecular geometries and DEC-RI-MP2 molecular gradients for the S12L test set and insulin.

Supplementary Material