An Integral-Direct Linear-Scaling Second-Order ... - ACS Publications

Subscriber access provided by UNIV OF CALIFORNIA SAN DIEGO LIBRARIES

Article

An integral-direct linear-scaling second-order Møller–Plesset approach Péter R. Nagy, Gyula Samu, and Mihaly Kallay J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b00732 • Publication Date (Web): 12 Sep 2016 Downloaded from http://pubs.acs.org on September 18, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

An integral-direct linear-scaling second-order Møller–Plesset approach Péter R. Nagy,∗ Gyula Samu, and Mihály Kállay∗ MTA-BME Lendület Quantum Chemistry Research Group, Department of Physical Chemistry and Materials Science, Budapest University of Technology and Economics, H-1521 Budapest, P.O.Box 91, Hungary E-mail: [email protected]; [email protected]

Abstract An integral-direct, iteration-free, linear-scaling, local second-order Møller–Plesset (MP2) approach is presented, which is also useful for spin-scaled MP2 calculations as well as for the efficient evaluation of the perturbative terms of double-hybrid density functionals. The method is based on a fragmentation approximation: the correlation contributions of the individual electron pairs are evaluated in domains constructed for the corresponding localized orbitals, and the correlation energies of distant electron pairs are computed with multipole expansions. The required electron repulsion integrals are calculated directly invoking the density fitting approximation, the storage of integrals and intermediates is avoided. The approach also utilizes natural auxiliary functions to reduce the size of the auxiliary basis of the domains and thereby the operation count and memory requirement. Our test calculations show that the approach recovers 99.9 % of the canonical MP2 correlation energy and reproduces reaction energies with an average (maximum) error below 1 kJ/mol (4 kJ/mol). Our benchmark ∗

To whom correspondence should be addressed

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

calculations demonstrate that the new method enables MP2 calculations for molecules with more than 2300 atoms and 26000 basis functions on a single processor.

1

Introduction

The second-order Møller–Plesset approach (MP2) 1 is one of the simplest ab initio electron correlation methods. Since its development in the thirties of the previous century, the MP2 method has become a standard tool of quantum chemistry. Its relatively high accuracy and low cost make MP2 applicable to considerably bigger systems than what can be treated with more advanced correlation models, such as the configuration interaction or coupledcluster (CC) methods. Even though MP2 is a well-established theory, the development of improved MP2 variants is still an active field of quantum chemistry. In recent years several spin-component-scaled MP2 (SCS-MP2) schemes have been invented, 2–4 which supply more accurate results at the same or even lower cost. It is also notable that the so-called doublehybrid density functionals, 5 which have gained wide popularity over the last decade, also include an MP2-like post-Kohn–Sham correction. Besides the theoretical improvements there have been numerous attempts over the years to reduce the computational expenses of MP2, which scale as the fifth-power of the system size. It has been known for a long time that the costs of MP2 calculations can be significantly decreased by the resolution-of-identity or density-fitting (DF) techniques, 6 though these approaches do not alter the formal scaling of the method. Reduced or even linear scaling can be achieved by the atomic orbital (AO) or domain-based local MP2 approximations. Most of the AO-based MP2 algorithms rely on the idea of Almlöf, who realized that the energy denominators appearing in various correlation methods can be got rid of by a Laplace transform (LT) technique. 7,8 Using this idea Ayala and Scuseria obtained MP2 and CC expressions in the AO basis and showed that linear scaling with molecular size can be achieved exploiting the decay properties of various contributions. 9,10 The AO-based MP2

2


Page 2 of 58

Page 3 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


technique was later refined by Ochsenfeld and co-workers, 11,12 who combined the AO-MP2 algorithm with rigorous integral estimates 13 to fully utilize the distance decay of charge interactions within the AO basis. Further significant improvements were shown when using Cholesky-decomposed pseudodensity matrices and DF approximation for the electron repulsion integrals (ERIs). 14,15 A related massively parallel MP2 algorithm was published by VandeVondele et al., 16 which employs a hybrid Gaussian and plane wave basis and is particularly suitable for condensed-phase systems. The key feature of this scheme is that mixed AO-molecular orbital (MO) ERIs are used, which are computed by direct integration between the products of AOs and the electrostatic potential of the given occupied-virtual pair density. In addition to the simple AO basis low-order scaling correlated methods can also be developed using non-orthogonal orbitals. Maslen and Head-Gordon published an iteration-free local MP2 (LMP2) model using non-orthogonal occupied and virtual orbitals. 17 More recently the scaled opposite spin MP2 method was also implemented invoking similar techniques. 18 A general theory for the application of biorthogonal projected AO bases was developed by Weijo and co-workers, 19 who also demonstrated the advantages of their formalism by implementing a reduced-scaling MP2 approach. Here we should also mention a promising new approach, the stochastic MP2 algorithm of Neuhauser and co-workers, 20 which was demonstrated to yield accurate results for systems of thousands of electrons with semiempirical approximations. The common feature of the domain-based local MP2 methods is that occupied localized MOs (LMOs) are used, domains of atoms, AOs, auxiliary functions, and other quantities are constructed, and the indices of particular quantities, such as ERIs, wave function parameters, or other intermediates, are restricted to be in the domains. Usually the truncation thresholds for the assembly of the domains are less rigorous than those used in AO-based algorithms, but if they are tightened, the energy also converges to that obtained by the parent method. These approaches can be further classified as the methods which do not divide the system into small fragments and the ones which do.

3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The foundations for the methods in the first category were laid down in the 1980s by Pulay 21 and Kapuy et al., 22 who pointed out that the short-range nature of the electron correlation can be exploited if LMOs are used. Important steps towards the practical utilization of these ideas were made by Pulay and Saebø, who reformulated the Møller–Plesset perturbation theory in an orbital-invariant form permitting the use of arbitrary orbitals, and also reported the implementation of the approach. 23,24 They assigned local correlation domains to the occupied LMOs consisting of AOs projected onto the virtual space. These non-orthogonal redundant basis functions, the so-called projected atomic orbitals (PAOs), also attained great popularity among the developers of local correlation methods. Following the pioneering work of Pulay local correlation theories were further developed by Werner and Schütz and their co-workers. 25–32 These authors constructed orbital domains for each electron pair and classified the electron pairs according to distance criteria as strong, weak, distant, and very distant. Neglecting very distant pairs, treating distant pairs by multipole approximations, and using efficient integral prescreening algorithms linear scaling for LMP2 was achieved. 26,27 An important milestone along this path was the combination of local correlation and DF. 28 It not only significantly accelerated the LMP2 approach but also allowed for the assessment of local approximations for extended molecules, which also guided the way for the developers of local correlation models. Another notable step was the introduction of explicit correlation into LMP2 theory. 29 The resulting explicitly-correlated LMP2 method, which also makes use of the DF approximation, permitted MP2 calculations for medium-sized molecules with near basis set limit accuracy. The usefulness of Almlöf’s LT technique in the Pulay-type local methods was also investigated, and a LT DF-LMP2 method was developed and compared to previous LMP2 schemes. 30 It was concluded that the performance of the LT and the original approaches are comparable. The scope of the LMP2 methods developed by the Stuttgart group was also extended to the calculation of molecular properties by implementing analytical energy gradients 31 and nuclear magnetic resonance chemical shifts, 33 whose evaluation was later considerably speeded up by DF techniques. 32,34

4


Page 4 of 58

Page 5 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Recently the research of local correlation theories gained momentum after Neese revitalized the pair natural orbital (PNO) method. 35 The PNOs proposed by Edmiston and Krauss 36 and applied extensively by Meyer 37 are special MP2 natural orbitals (NOs) constructed for each electron pair separately allowing for a compact local representation of the correlated wave function. The success of PNOs also initiated the development of other NOlike correlation orbitals, which include the orbital-specific virtual (OSV) orbitals of Manby and co-workers 38 and the local NOs (LNOs) proposed by us. 39 Both PNOs 40 and the NOlike orbitals 38,41,42 were used to implement low-order scaling local MP2 algorithms. A cubicscaling hybrid OSV-PNO MP2 approach was also proposed by Hättig et al., 43,44 while a linear-scaling parallel LMP2 algorithm was developed by Werner and co-workers 45,46 combining PAOs, OSVs, and PNOs. The most recent advance in this direction is the PNO-based LMP2 implementation of Neese et al. 47 utilizing a general framework for the sparse representation of tensors. The development of local correlation approaches in the second category, the so-called fragmentation-based methods, also dates back to the eighties, when Förner and co-workers 48 developed a local CC doubles method. In that scheme domains of localized orbitals were defined, and the CC correlation energy was estimated as the sum of contributions evaluated within these domains and for their pairs. In the last three decades a series of similar methods was published the common feature of which is that the system is split up into several fragments, and contribution of the latter to the correlation energy is computed with the aid of the algorithms developed for the corresponding canonical correlation methods. Most of these schemes were also implemented at the MP2 level. In the incremental method proposed by Stoll 49 the fragment correlation energies are corrected by increments accounting for the correlation of pairs, triples, etc. of fragments until the desired accuracy is reached. A fully automated, parallel implementation of the incremental scheme for MP2 was reported by Friedrich. 50 The closely-related fragment molecular orbital method (FMO) developed by Kitaura and co-workers 51 was also implemented for MP2. 52

5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The divide-and-conquer local MP2 approaches proposed by Li and Li 53 and Kobayashi and Nakai 54 also decompose a large system into various fragments, and the correlation energy of the whole system is obtained as the sum of the contributions of the fragments and their surroundings. In the divide-expand-consolidate (DEC) family of approaches advocated by Jørgensen et al. 55,56 the fragments are defined by the individual atoms and their environment, and the correlation energy is decomposed into the contribution of the fragments and their pairs. At the MP2 level not only energies but analytical gradients are also available for the approach. 57 In the cluster-in-molecule (CIM) approach of Li and co-workers a domain is assembled for each occupied LMO, and the energy expression of correlation methods is expressed as the sum of the correlation contributions of occupied LMOs evaluated in the domain. 58,59 Several MP2 implementations of this model were reported, 42,60–62 and multilevel extensions of the theory are also available where the active center of a big molecule is described with a CC method, while its environment is treated at the MP2 level. 63 The CIM approach was recently combined with the FMO method by Gordon et al., who also reported a massively parallel MP2 implementation of the scheme. 64 In our previous papers 39,42 we successfully employed LNOs within the CIM framework to decrease the number of MOs in the fragments, and later we also proposed a modified ansatz reminiscent of the incremental scheme in which, to correctly account for long-range correlation, the pair correlation energies of distant electrons pairs are also calculated. 65 In this paper we report the development of an iteration-free, integral-direct, linear-scaling MP2 approach, which is a significantly improved version of our previous algorithm. 65 We have carefully re-designed the domain construction algorithms, completely avoided the storage of integrals or intermediates, and utilized natural auxiliary functions (NAFs). The new approach is by about an order of magnitude faster than the original one making it highly competitive with other low-order scaling MP2 algorithms, which will be demonstrated by benchmark calculations.

6


Page 6 of 58

Page 7 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2

Theory and implementation

2.1

The local MP2 Ansatz

The present local MP2 Ansatz is based on our recent linear-scaling scheme developed for direct random phase approximation (dRPA), MP2, and related methods. 65 Here we employ ideas from fragmentation approaches, such as the incremental expansion 49,66 up to orbital pairs, while the contribution of the individual fragments are evaluated independently in a manner motivated by the CIM 39,42,58,60,67 model. The Ansatz is introduced assuming that spin-restricted Hartree–Fock (HF) orbitals are available for a closed-shell system. Indices i, j , . . . ( a, b, . . . ) will refer to the spatial part of the occupied (virtual) canonical HF orbitals. Let us introduce the basic notation starting from the conventional MP2 correlation energy expression: c EMP2

=

Xh abij

i

2(ai|bj) − (aj|bi) tai,bj ,

(1)

where tai,bj ’s are the MP2 amplitudes. The four-center ERIs, which are denoted by Kai,bj = (ai|bj) using the Mulliken convention, will be evaluated by relying on the DF approximation, which reads for the ERIs as: K = IV−1 I = JJT ,

(2)

J = I(L−1 )T

(3)

where

with I containing the three-center two-electron integrals, that is, Iai,P = (ai|P ) and P referring to auxiliary basis functions. Additionally, (L−1 )T denotes the transpose of L−1 . Several alternatives have been proposed for the construction of matrix J. Here we choose to subject the two-center two-electron integral matrix of the auxiliary basis functions, with elements VP Q = (P |Q), to Cholesky-decomposition (CD) as V = LLT , where L is the lower 7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 58

triangular Cholesky-matrix. From that it is straightforward to find the V−1 = (L−1 )T L−1 expression and consequently the form of the fitting step in Eq. (3). Let us now rewrite Eq. (1) as the sum of contributions from individual orbitals as X

δEi

(4)

i 2(ai|bj) − (aj|bi) tai,bj .

(5)

c EMP2 =

i

with δEi =

Xh abj

The occupied and virtual orbitals can be assumed to be localized in Eq. (5) since the expression is invariant to separate unitary transformations of the occupied and virtual MOs. We can also suppose that PAOs 21 are constructed, which span the virtual space. The domain approximation to evaluate δEi —i.e., a truncation in the occupied and virtual subspaces of Eq. (5)—can be simply formulated in terms of LMOs and PAOs. First, those localized orbitals (j) are collected that strongly interact with LMO i, then PAOs are selected to describe the correlation of orbital i with the selected j’s. The domain of these LMOs and PAOs will be referred to as the extended domain (ED) of i and will be denoted by Ei . The occupied orbitals of Ei are selected based on the pair correlation energies of LMO i with all other LMOs. To that end, first, a smaller domain, hereafter referred to as primary domain (PD), is assembled for each occupied LMO. Then an approximate opposite-spin (OS) MP2 pair correlation energy is evaluated for the i–j LMO pair in the pair domain Pij including the PDs of i and j whose construction will be discussed in detail in Sects. 2.2 and 2.3. If this approximate pair correlation energy, hereafter denoted by δEij (Pij ), is greater than a threshold, the pair i–j is considered a strong pair and j is added to the ED of i and vice versa. On the other hand, if δEij (Pij ) is below the threshold, the pair i–j is considered a distant pair, and their contribution to the MP2 correlation energy is only taken into account as a second-order (pair) increment. Accordingly, the LMP2 correlation energy expression

8


Page 9 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


has the following form: c ELMP2 =

X

δEi (Ei ) +

i

X′

δEij (Pij ),

(6)

i Itol . Moreover, the transformation of Eq. (14) is only carried out for those LMOs, for which a tighter criterion is also met, that is Ea˜ν,P max Cνi > Itol . This second estimate is tighter than the one in ν shell

Eq. (17) because it corresponds to only a single LMO, the maximum is not taken for the i’s of the ED. We have found that Itol = 10−7 Eh can be chosen as a conservative default threshold which causes error on the microhartree scale in the LMP2 correlation energies. 19



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 58

In the following step, the occupied index of (˜ a i|P ) is transformed to the pseudo-canonical basis of the ED yielding the (˜ a ˜i|P ) set of integrals. The number of auxiliary functions in (˜ a ˜i|P ) can be roughly halved half by introducing NAFs 65,82 for the fitting of the pseudocanonical LMO-PAO charge densities of the ED. This leads to significant savings during the calculation of the δEi (Ei ) fragment energies for which both the operation count and the memory requirement is proportional to the size of the fitting basis. First, the WP′ Q =

X

(˜ a ˜i|P ) (˜ a ˜i|Q)

(18)

˜i a ˜

matrix is composed and fitted as W = (L−1 )T W′ L−1 . For that purpose, integrals (P |Q) are computed for the atoms of the LDF domain of all EDs from scratch to avoid the storage of (P |Q) for the whole molecule, a CD is performed on matrix VP Q = (P |Q), and the resulting lower triangle Cholesky-matrix is inverted leading to L−1 . The eigenvectors of W, i.e., the NAFs with eigenvalue smaller than a threshold, εNAF , are disregarded, while the remaining functions (P˜ ) will constitute the auxiliary basis of the ED, and their coefficients are collected into matrix N. Finally, the auxiliary index of (˜ a ˜i|P ) is fitted and transformed to the NAF basis in a single step with transformation matrix (L−1 )T N, arriving at the Ja˜ ˜i,P˜ integrals. Fragment energies in the extended domains Fragment energies of the EDs are evaluated similarly to the scheme outlined in Ref. 65. The main differences are that δEi (Ei ) is given in terms of canonicalized PAOs (instead of LNOs) and a NAF basis is constructed already in the EDs. In order to minimize the inaccuracy originating from the truncation of the NAF basis, 65,82 the ERIs of Ei , required for the energy calculation, are assembled using the full auxiliary basis of Ei as follows: (˜ a i|˜b ˜j) =

X PQ

−1 T

(˜ a i|P ) (L ) L

20

−1

PQ

(˜b ˜j|Q) ,


(19)

Page 21 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where the second index refers to only the single LMO, i. Next, MP2 amplitudes are evaluated, where the energy-denominators are factorized with either CD 83 (which is the second type of CD, hence it is denoted as CD2) or LT: 7,8 ˜

tia˜˜jb = −

XX ω

Ja˜ωi,P˜ J˜bω˜j,P˜ ,

(20)

P˜

with the sum over ω running for all Cholesky-vectors of CD2 or the integration quadrature points (QP) used in the LT. Integrals J˜bω˜j,P˜ are constructed as J˜bω˜j,P˜ = J˜b˜j,P˜ (eω )˜b˜j ,

(21)

using the corresponding elements of the Cholesky-vectors of CD2, (eω )˜b˜j or in the case of LT with (eω )˜b˜j =

√

wω e−(ε˜b −ε˜j )tω ,

(22)

where tω is the QP with wω as the corresponding weight. The remaining Ja˜ωi,P˜ integrals are obtained by transforming index ˜i of the Ja˜ω˜i,P˜ list back to LMO i. Finally, with the MP2 amplitudes of Eq. (20) and the ERIs of Eq. (19) at hand, the MP2 correlation energy contribution of the ED is computed as

δEi (Ei ) =

Xh a ˜˜b˜ j

i ˜ 2(˜ ai|˜b˜j) − (˜ a˜j|˜bi) tia˜˜jb ,

(23)

which concludes the present LMP2 algorithm. Scaling of the algorithm The formal scaling of all the operations with the various dimensions of the domains is presented in detail in our previous report. 65 We would like to merely repeat here that all manipulations in the EDs are asymptotically independent of the system size. For the sake

21



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

of completeness, we add that there are still higher than linear scaling steps in the LMP2 algorithm, which, however, usually take less than 2-3 % of the whole computation time. These steps are the PAO construction (Eq. 12) and the OS-MP2 pair energy computation for the LMO pairs. Additionally, the pre-requisite SCF calculation and Boys localization are at least cubic scaling, the former being significantly more demanding than LMP2 most of the time. Considering the memory requirement, again, the size of all the arrays within an ED is asymptotically linear scaling. For the sake of simplicity, we still have a few arrays stored in memory throughout the entire LMP2 computation, whose size is linear or quadratical scaling. However, the storage of such arrays, e.g., for the Fock matrix, untruncated PAO coefficients, and certain truncation maps, have not caused any problem so far for cases with up to 26000 basis functions.

2.4

Local fitting domains for HF calculations

The solution of the HF SCF equations for systems of several hundreds or thousands of atoms is as challenging as the subsequent correlation calculation. If DF is used, the Coulomb term is relatively cheap, and the evaluation of the exchange contribution is the rate-determining operation. 84 To break down the unfavorable, formal fourth-power scaling of the latter, several algorithms have been developed, 85–87 of which we previously adopted 42 the one proposed by Polly et al. 85 In this approach local fitting domains including a small subset of the original auxiliary basis are constructed for each occupied MO, and the generalized densities of the products of the MO and the AOs are fitted in this restricted set of functions. To select the atoms whose auxiliary functions are included in the domain, in each SCF iteration step, the MOs are localized, and for each LMO Löwdin atomic charges are computed. All atoms are selected which have a charge greater than 0.05, and all further atoms are included whose distance from the aforementioned atoms is smaller than a threshold, by default 5 bohrs. At the end of the SCF procedure an extra iteration is performed with a bigger threshold to get 22


Page 22 of 58

Page 23 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


an accurate HF energy. We have found that the accuracy of the final HF energy obtained with this algorithm is strongly basis set dependent. Obviously for basis sets including diffuse functions larger fitting domains would be required for the same accuracy than for other bases. However, the above scheme based on Löwdin charges and distance criteria is not sensitive enough to the properties of the basis set and may yield poor energies in particular cases. To obviate this problem, inspired by the success of the energy-based cutoffs in the local correlation treatment, we replaced the distance threshold by another one related to the ERIs over the AOs residing on the corresponding atoms. At the evaluation of the exchange contribution of occupied LMO i (µi|P )-type integrals are required, where µ stands for an arbitrary AO, and P is an auxiliary function of the local fitting domain. These integrals can be expanded in terms of AO integrals (µν|P ), which can be estimated by the notorious Schwarz inequality as 88 |(µν|P )| ≤

p p (µν|µν) (P |P ).

(24)

Consequently, the contribution of basis functions on atomic centers A and B to the exchange term is roughly proportional to the

XAB = max

µ∈A,ν∈B

p p (µν|µν) max (P |P ) P ∈A,B

(25)

quantity, where µ ∈ A refers to the AOs placed on atom A. We now suppose that the bigger the contribution of an atomic center to the exchange term is, the more important the fitting functions on the center are for the accurate fitting of the corresponding charge distributions. Thus, atom A is included in the local fitting domain of occupied LMO i if XAB is greater than a threshold for any atom B assigned to i based on the Löwdin charges. In practical calculation we have found that this scheme considerably reduces the basis set dependence of the results though does not eliminate entirely. Our experience shows that a threshold of 1 Eh is suitable in the regular SCF iterations. With this choice, for 23



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

basis sets excluding diffuse functions fitting domains of similar size are assembled as with the original distance-based selection algorithm and the criterion of 5 bohrs, that is, usually the first neighbors of the atoms assigned to the occupied LMO are included. With diffuse basis sets the local fitting domains are typically larger using the new algorithm. With a threshold of 10−3 Eh in the extra SCF iteration, an accuracy of a couple of µEh /atom can be achieved with respect to the conventional DF-HF energy if the basis set does not contain diffuse functions. For diffuse basis sets a more rigorous threshold, 10−5 Eh , is recommended. We also note that we replaced the conventional PM localization performed in the SCF iterations by a modified one where Löwdin charges are used instead of Mulliken charges. The advantage of the modified localization scheme lies in its increased speed and the slightly smaller fitting domains.

3

Results

In this section the convergence of the correlation energy with the various truncation thresholds are monitored, appropriate default values for the thresholds are proposed, the performance of our LMP2 method is benchmarked for correlation energies, reaction energies, and weak interactions, and finally representative timings are presented. Here we confine ourselves to the discussion of the results obtained with MP2 and note that the performance of the spin-scaled MP2 variants is virtually identical (see Tables 4 and 5 of the Supporting Information). Concerning double-hybrid functionals the performance of our algorithm was found to be similar in terms of the relative errors in the correlation energy for the LMP2 contribution to energies obtained by the B2PLYP 5 functional, the latter being presented in the Supporting Information. Note that absolute errors in the energy differences are significantly better for B2PLYP simply because the MP2-like term is scaled by 0.27 in the B2PLYP energy expression.

24


Page 24 of 58

Page 25 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3.1

Computational details

The LMP2 approach proposed in this paper has been implemented in the Mrcc suite of quantum chemical programs and will be available in the next release of the package. 89 In the calculations presented here Dunning’s (augmented) correlation-consistent polarized valence X-tuple-ζ basis sets [(aug-)cc-pVXZ, X = D, T, and Q], 90–92 the correlationconsistent polarized valence double- and triple-ζ bases for explicitly correlated wave functions (cc-pVXZ-F12, X = D and T) developed by Peterson et al., 93 as well as the split valence and triple-ζ valence basis sets including polarization functions (def2-SVP, def2-TZVP, and def2-TZVPP) bases optimized by Weigend and Ahlrichs 94 were used. For the gold atom the (augmented) correlation consistent valence double-ζ and triple-ζ pseudopotential (augcc-pVDZ-PP, and (aug-)cc-pVTZ-PP) basis sets of Peterson 95 were employed together with the corresponding effective core potential for the inner 60 electrons of the atom. 96 For all the AO basis sets the corresponding auxiliary bases of Weigend and co-workers were applied. 97–99 In the case of the cc-pVXZ-F12 AO basis sets the auxiliary bases corresponding to the aug-cc-pVXZ basis sets were applied. In all the HF and the reference canonical MP2 calculations the DF approximation was invoked. For systems with more than 5000 AOs the exchange contribution in the HF calculations was computed with local fitting domains as described in Sect. 2.4. The core electrons were kept frozen in all the correlation calculations. The Boys localization 73 was used throughout for the construction of occupied LMOs. The energy denominators within the EDs were removed by Cholesky decomposition. 83 The number of Cholesky vectors were automatically determined to make the diagonal elements of the residual matrix less than 10−4 (see Ref. 83 for details). For the performance analysis of the LMP2 Ansatz the following statistical error measures are utilized: maximum absolute error (MAX), mean absolute error (MAE), the root mean square (RMS) deviation, and standard deviation of the absolute errors (STD). Note that the latter measures the consistency of the errors. 25



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The reported computation times are wall-clock times determined on a machine with 64 GB of main memory and a 6-core 3.5 GHz Intel Xeon E5-1650 processor.

3.2

Benchmark sets and test systems

The accuracy of the LMP2 Ansatz is assessed with various benchmark sets. Reaction energies were evaluated for the test set of Neese, Wennmohs, and Hansen (NWH) assembled from 23 reactions including also molecules of 36 atoms 35 and the ISOL reaction energy set of Grimme and co-workers, 100 which is built from 48 organic molecules of 24 to 81 atoms. Weak interaction energies are tested for dimers collected in the S66 benchmark set of Řezáč and co-workers. 101 Additionally, we gathered a couple of weakly bound complexes from the literature which consists of much more extended molecules with up to 144 atoms. Four, mainly dispersion dominated complexes were obtained from the S12L set of Grimme, 102 namely 2b, 3b, 5a, and 6b, and two other H-bonded systems (GCGC and PHE) were taken from the L7 set of Hobza and co-workers. 103 We will refer to the set of the resulting six complexes here as S12L-L7. The performance of our LMP2 scheme was also assessed on the correlation energies of large molecules collected from the literature. The list includes beta-carotene, angiotensin, an amylose helix of 8 glucose units (amylose8 ), and DNA fragments of two and four adeninethymine base pairs (DNA2 and DNA4 ) from Ref. 11, sildenafil and vancomycin of Ref. 47, and the dalargin 42 hexapeptide (Tyr-D-Ala-Gly-Phe-Leu-Arg). Two additional systems are taken from Ref. 45 that are the reactant state of the amido-thiourea catalyzed enantioselective imine hydrocyanation (ATU for short) and a gold(I)-aminonitrene complex (AuA for short, with +1 charge). Moreover, some even larger systems were selected to illustrate the efficiency of our method, namely the crambin protein (with +2 charge) of Ref. 47, the 16-base-pair analogue of the above DNA fragments, DNA16 from Ref. 11, a water cluster with 569 molecules [(H2 O)569 ] of Ref. 13, and the core domain of the HIV-1 integrase enzyme (integrase for short, 26


Page 26 of 58

Page 27 of 58


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 58

where the appearance of Ei in the notation emphasizes that due to the truncations δEij (Ei ) depends on the parameters of the ED in which ED it is evaluated. Note that an additional factor of two appears to include the energies of both i–j and j–i pairs in the definition. The cc-pVTZ basis set was employed throughout this test, which allowed us to compute the above LMP2 pair energies with EDs including all atoms of the molecule. The resulting pair energies do not depend on the ED in this case, hence they will be simply referred to as δEij , but still not exact since the LMOs and the NAF basis had to be truncated according to the default settings, TEDo = 0.9998 and εNAF = 2 · 10−3 , respectively. In order to estimate the error originating from the pair approximations we introduce two error measures. The first one, E1 (εw ), gives the magnitude of the relative error in the pair c energy contribution to ELMP2 , the second term of Eq. (6). Explicitly, E1 (εw ) is the ratio of

the sum of the δEij (Pij ) pair energies which are lower than εw , i.e., which are included in Eq. (6), and the same sum computed with the reference δEij ’s, and it reads as

E1 (εw ) =

δEij (P ij )

An Integral-Direct Linear-Scaling Second-Order ... - ACS Publications

Recommend Documents