- Email: [email protected]
Implementation of Surján’s density matrix formulae for calculating second-order Møller–Plesset energy
Chemical Physics Letters 420 (2006) 250–255 www.elsevier.com/locate/cplett
Implementation of Surja´n’s density matrix formulae for calculating second-order Møller–Plesset energy Masato Kobayashi, Hiromi Nakai
*
Department of Chemistry, School of Science and Engineering, Waseda University, 3-4-1 Okubo, Shinjuku-ku, Tokyo 169-8555, Japan Received 20 October 2005 Available online 20 January 2006
Abstract We numerically assess the method for obtaining second-order Møller–Plesset (MP2) energy from the Hartree–Fock density matrix (DM) recently proposed by Surja´n [Surja´n, Chem. Phys. Lett. 406 (2005) 318]. It is confirmed that Surja´n’s method, referred to as DM-Laplace MP2, can obtain MP2 energy accurately by means of appropriate integral quadrature and a matrix exponential evaluation scheme. Numerical tests reveal that the Euler–Maclaurin and the Romberg numerical integration schemes can achieve milli-hartree accuracy with small quadrature points. This Letter also indicates the possibility of the application of DM-Laplace MP2 to linear-scaling selfconsistent field techniques, which give approximate DM. Ó 2006 Elsevier B.V. All rights reserved.
1. Introduction Second-order Møller–Plesset (MP2) perturbation theory [1] is the simplest molecular orbital (MO) method that can deal with electron correlation. MP2 is also the most practical method enabling large-scale MO calculations including correlation effects. So far, many researchers have proposed various techniques for accelerating the evaluation of MP2 energy in terms of the localized orbital [2–4], the resolution of the identity (RI) technique for integral evaluations [5], and the Laplace transformation of the energy denominator [6–11]. The closed-shell MP2 correlation energy is given by using spatial occupied orbitals i, j, . . ., virtual orbitals a, b, . . ., andRRstandard electron repulsion integral notation of ðiajjbÞ ¼ iðr1 Þjðr2 Þr1 12 aðr1 Þbðr2 Þ dr1 dr2 as follows: ð2Þ Ecorr
occ X vir X ðiajjbÞ½2ðiajjbÞ ðibjjaÞ ¼ ; ea þ eb ei ej i;j a;b
ð1Þ
where ei and ea are the orbital energies of occupied and virtual orbitals, respectively. In the recent Letter by Surja´n [12], *
Corresponding author. Fax: +81 3 3205 2504. E-mail address: [email protected] (H. Nakai).
0009-2614/$ - see front matter Ó 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2005.12.088
the explicit MP2 correlation energy functional Eð2Þ corr ½P, where the argument function P is the Hartree–Fock (HF) density matrix (DM), was presented. Surja´n’s strategy is based on so-called atomic orbital (AO) Laplace MP2 form [8,11]. Eq. (1) can be rewritten by means of the Laplace transformation as Z 1 s X ð2Þ Ecorr ¼ e2 ðsÞ ds wa e2 ðsa Þ ð2Þ 0
a¼1
with s quadrature points {sa} and their weights {wa}. The integrand of Eq. (2) is given in AO basis by XX e2 ðsa Þ ¼ X alc Y amd X ajk Y aer ðcdjjeÞ½2ðlmjkrÞ ðlrjkmÞ cdje lmkr
¼
X
ðlmjk rÞ½2ðlmjkrÞ ðlrjkmÞ;
ð3Þ
lmkr
where Xa and Ya are the energy-weighted DMs of electron and hole, respectively, and are given by Xa ¼
occ X
eei sa Ci CTi ;
ð4Þ
eea sa Ca CTa .
ð5Þ
i
Ya ¼
vir X a
M. Kobayashi, H. Nakai / Chemical Physics Letters 420 (2006) 250–255
According to Surja´n [12], the energy-weighted DMs of Eqs. (4) and (5) can be represented by means of AO-based DMs instead of orbital energies as follows: X a ¼ e sa S a
Y ¼e
1
sa S
F
P;
1
F
Q;
ð6Þ ð7Þ
where S is an overlap matrix, F is a Fock matrix, and P and Q are the HF DMs of electron and hole in AO basis, P¼
occ X
Ci CTi ;
ð8Þ
Ca CTa ;
ð9Þ
i
Q¼
vir X a
respectively. Note that we only assume the Hartree–Fock– Roothaan equation FCi = eiSCi so far. Furthermore, an inverse overlap matrix can be eliminated by assuming the commutation relationship of FPS = SPF and idempotency of PS. Finally, the following working formulae are obtained: Xa ¼ esa PF P; a
Y ¼e
sa QF
Q.
ð10Þ ð11Þ
Therefore, MP2 energy can be obtained as the functional of DMs of electron P and hole Q. This method for evaluating MP2 energy is referred to as DM-Laplace MP2. It should be noted that these DMs satisfy P þ Q ¼ S1 .
ð12Þ
DM-Laplace MP2 method has prospects of enabling large-scale MP2 calculation. DM-Laplace MP2 requires only HF DM which can often be constructed by linearscaling SCF methods [13,14]. Furthermore, the method is based on the AO-Laplace MP2 scheme whose computational cost can also be reduced to scaling linearly depending on the system size and is well assessed by Ayala and Scuseria [11]. Although Surja´n suggested the applicability of the DM-Laplace MP2 method, numerical results have not been shown yet. In this Letter, the DM-Laplace MP2 calculation code was implemented to the GAMESS quantum chemistry package [15] and numerical assessments of the method were carried out. We indicate two problem areas in the DM-Laplace MP2 calculation and provide feasible solutions. We first attend to the numerical quadrature of Eq. (2) in Section 2. Next, a method of evaluating energy-weighted DMs X and Y is discussed in Section 3. In addition, in Section 4, we examine whether the DM-Laplace MP2 can apply to a method for constructing approximate DMs. Finally, concluding remarks are presented in Section 5.
251
gested the efficiency of a least-square fit of the 1/x function over Dijab ! s X X 1 wa expðDijab sa Þ ¼ min !; ð13Þ Dijab a¼1 ijab where Dijab ¼ ea þ eb ei ej .
ð14Þ
This scheme provides quite accurate results even with a small number of quadrature points. For instance, according to Ha¨ser and Almlo¨f [7], micro-hartree accuracy can be achieved with the number of quadrature points s = 8– 10. This scheme, however, is not suited for the DM-Laplace MP2 calculation because the DM-Laplace MP2 was proposed as a method for evaluating MP2 energy without the use of orbital energies, which appear in Eq. (14). Wilson and Almlo¨f [10] have reported the exponentially decaying behavior of the integrand of Eq. (2), i.e., e2(s) with the value of s. The Gauss–Laguerre quadrature scheme is suitable for the semi-infinite integral whose integrand can be approximated as the product of polynomials and an exponential function [16]. Fig. 1 shows the decay of e2(s) with respect to s in the Laplace-MP2 calculation of a benzene molecule with the 6-31G basis set [17]. In Fig. 1, the Gauss–Laguerre quadrature points corresponding to s = 10 are also shown by small squares. Since e2(s) is a quasi single exponential function, as seen from Fig. 1, it is expected that the Gauss–Laguerre quadrature is suited for the numerical integration of Eq. (2). On the other hand, Ayala and Scuseria [11] have suggested the efficiency of using Euler–Maclaurin summation [18,19] for this numerical quadrature, with a proper change of variable (CoV) from s to r: Z 1 Z 1 Z 1 ds ð2Þ Ecorr ¼ e2 ðsÞ ds ¼ e2 ðsÞ dr ¼ f2 ðrÞ dr dr 0 0 0 " # s X 1 k 1 f2 ð15Þ þ ðf2 ð0Þ þ f2 ð1ÞÞ . s þ 1 k¼1 sþ1 2
2. Numerical quadrature In the Laplace MP2 calculations, numerical quadrature points {sa} and weight parameters {wa} have to be determined in some fashion. Several authors [7,10,11] have sug-
Fig. 1. Decay behavior of integrand e2(s) for benzene (6-31G). (h) Gauss– Laguerre quadrature points for s = 10.
252
M. Kobayashi, H. Nakai / Chemical Physics Letters 420 (2006) 250–255
The integration error dEð2Þ corr of the above formula can be evaluated as 2k h 1 i X B2k 1 ð2k1Þ ð2k1Þ ð2Þ ¼ f2 ð1Þ f2 ð0Þ dEcorr ð2kÞ! s þ 1 k¼1 ¼
f20 ð0Þ 12ðs þ 1Þ
2
ð5Þ
f2000 ð0Þ
f2 ð0Þ þ 720ðs þ 1Þ4 30240ðs þ 1Þ6
ð7Þ
f2 ð0Þ 1209600ðs þ 1Þ
8
þ
ð16Þ
because e2(s) and its all higher derivatives are zero at s = 1, where B2k are Bernoulli numbers. Although arbitrary CoV choice is possible in this scheme, it is best to choose one for which its Jacobian ds/dr and as large a number of its derivatives as possible are zero for r = 0. Here, we apply two formulae for the CoV; type (A) is s¼
r2 ðr3 þ r4 þ r5 þ 0:9r6 Þ=4 2
ð17Þ
þ r2 tanðpr=2Þ.
ð18Þ
ð1 rÞ
where T nð0Þ is the approximate value of the integral by the trapezoid method with 2n division " # s X 1 k 1 ð0Þ Tn ¼ f2 ð20Þ þ ðf2 ð0Þ þ f2 ð1ÞÞ ; s þ 1 k¼1 sþ1 2 with s = 2n 1. The explicit formulae for n = 3 (s = 7) and n = 4 (s = 15) are 2 1 1 3 ð3Þ 256f 2 T3 ¼ þ 88f 2 þ 256f 2 2835 8 4 8 1 5 3 7 þ109f 2 þ 256f 2 þ 88f 2 þ 256f 2 2 8 4 8 ð21Þ and ð4Þ T4
and type (B) s¼
r3 0:9r4 ð1 rÞ
2
Their Jacobians are both zero, and the first derivative of type (B)’s Jacobian also becomes zero for r = 0. We additionally tested the Romberg quadrature method [16,18], which is the extrapolation of some results of the extended trapezoid method based on the idea of Richardson’s deferred approach to the limit, with the same CoV, ð2Þ i.e., Eqs. (17) and (18). Specifically, Ecorr can be approxiðnÞ mated to T n by means of the following recurrence relation: ðk1Þ
T ðkÞ n ¼
4k T ðk1Þ T n1 n ; k 4 1
ð19Þ
2 1 1 32768f 2 ¼ þ 11008f 2 722925 16 8 3 1 5 þ32768f 2 þ13864f 2 þ 32768f 2 16 4 16 3 7 1 þ 11008f 2 þ 32768f 2 þ 13779f 2 8 16 2 9 5 11 þ 32768f 2 þ 11008f 2 þ 32768f 2 16 8 16 3 13 þ 13864f 2 þ32768f 2 4 16 7 15 þ 11008f 2 þ 32768f 2 ; ð22Þ 8 16
respectively. We have implemented the AO-Laplace MP2 and the DM-Laplace MP2 calculation code in the GAMESS quantum chemistry package [15] based on Ha¨ser’s quadratic
Table 1 Dependence of MP2 correlation energy of benzene on quadrature method [hartree] Method
s
6-31G
6-31G*
Eð2Þ corr
(Diff.)
Eð2Þ corr
(Diff.)
Gauss–Laguerre
5 10 15 50
0.512724 0.519626 0.521514 0.523061
(+0.010345) (+0.003444) (+0.001556) (+0.000009)
0.733451 0.766420 0.774904 0.784609
(+0.051311) (+0.018342) (+0.009858) (+0.000152)
Euler–Maclaurin (A)
5 10 15 50
0.516278 0.521006 0.522094 0.522972
(+0.006791) (+0.002064) (+0.000976) (+0.000098)
0.770241 0.780406 0.782704 0.784557
(+0.014521) (+0.004355) (+0.002057) (+0.000205)
Euler–Maclaurin (B)
5 10 15 50
0.522936 0.523065 0.523069 0.523068
(+0.000133) (+0.000004) (+0.000001) (+0.000002)
0.784540 0.784761 0.784760 0.784758
(+0.000221) (+0.000000) (+0.000001) (+0.000003)
Romberg (A)
7 15
0.522629 0.523123
(+0.000441) (0.000053)
0.784643 0.784666
(+0.000118) (+0.000096)
Romberg (B)
7 15
0.523683 0.523227
(0.000613) (0.000157)
0.783803 0.784919
(+0.000958) (0.000157)
Canonical
–
0.523070
0.784761
M. Kobayashi, H. Nakai / Chemical Physics Letters 420 (2006) 250–255
scaling algorithm [8]. In this algorithm, six steps of integral cut-off using Schwarz’s inequality are adopted (see Refs. [8,11] for details). We applied Schwarz’s inequality threshold of h = 1011 a.u. in all the following calculations. Table 1 shows the all-electron MP2 correlation energy of a benzene molecule obtained from the DM-Laplace scheme by several quadrature methods. The Chebyshev (PF) scheme was used for the evaluation of energyweighted DMs, which is explained in detail in Section 3. The 6-31G [17] and 6-31G* [20] basis sets were adopted. ‘(A)’ and ‘(B)’ on the quadrature scheme mean the CoV. The discrepancy with respect to the energy obtained from conventional canonical MP2 calculation, which is expressed in the bottom row of the table, is also shown in parentheses. The correlation energy of each quadrature method approaches the canonical one as the number of quadrature points s increases. The energies on the Gauss–Laguerre scheme with s = 15 have errors of about 1.6 and 9.9 milli-hartrees with the 6-31G and 6-31G* basis, respectively. The Euler–Maclaurin and the Romberg schemes, however, achieve milli-hartree accuracy even with 5–7 quadrature points if a suitable CoV is adopted. In the Euler–Maclaurin quadrature, type (B) CoV gives a considerably more precise estimate of energy (attaining micro-hartree accuracy with s = 10) than type (A) because the Jacobian as well as its first derivative become zero for r = 0 and the first term of the r.h.s. of Eq. (16) vanishes in type (B) CoV. On the other hand, the choice of CoV makes little difference in the Romberg method since the method gives the extrapolation to the Euler–Maclaurin summation. The Euler–Maclaurin scheme gives the upper bound of the canonical MP2 energy. This should be true in any case in type (A) CoV since the first term of the r.h.s. of Eq. (16) remains as f20 ð0Þ ¼ 2e2 ð0Þ with e2(0) > 0. In type (B) CoV, however, this behavior does not necessarily satisfy and it depends on the shape of e2(s). 3. Evaluation of energy-weighted density matrix We have proposed two working formulae for the evaluation of energy-weighted DMs X and Y in the DM-Laplace MP2. Eqs. (6) and (7) referred to as S1F formulae require the inverse overlap matrix, while Eqs. (10) and (11) referred to as PF formulae do not. In principle, both formulae should give the same result if a Fock F and a DM P satisfy the commutation relationship of FPS = SPF, which is fulfilled in the self-consistent field (SCF) solution. Whichever formula is used for energy-weighted DMs, the calculation of a matrix exponential is necessary. A number of methods for evaluating a matrix exponential have been presented. We apply two rational approximation algorithms for exp(x) function, that is, (6,6)-degree Pade´ and (14,14)-degree Chebyshev expansion using the EXPOKIT software package [21]. In the practical calculation of matrix exponentials using the Pade´ or the Chebyshev approximation, two working equations, i.e., S1F formulae and PF formulae, may no longer provide the same
253
Table 2 Dependence of MP2 correlation energy of benzene on energy-weighted DM evaluating method [hartree] Method
6-31G
6-31G*
Pade´ (S1F) Pade´ (PF) Chebyshev (S1F) Chebyshev (PF)
–a 0.5230654543 0.5230654543 0.5230654543
–a 0.7847609764 0.7847609764 0.7847609764
AO-Laplace
0.5230654543
0.7847609764
a
Energy-weighted DM diverged for a large s value.
matrix due to the character of the algorithm and/or numerical stability. Table 2 shows the DM-Laplace MP2 correlation energy of a benzene molecule by using the Pade´ and the Chebyshev approximation for the evaluation of matrix exponentials with the 6-31G and 6-31G* basis sets. We used the Euler–Maclaurin (B) scheme with s = 10. The results in the AO-Laplace scheme, in which MOs and orbital energies are used to construct energy-weighted DMs by Eqs. (4) and (5), are also shown. Except for the Pade´ approximant case using S1F formulae, where the energy-weighted DM diverges for a large s value, correlation energies obtained from the DM-Laplace MP2 calculation completely accord with those from the AO-Laplace MP2 with 1010 hartree accuracy. It is found that the present DM-Laplace MP2 scheme is equivalent both in theory and in practice to the AO-Laplace MP2 as long as an appropriate method for evaluating the energy-weighted DM is chosen. 4. DM-Laplace MP2 calculation using approximate density matrix Surja´n [12] has suggested the applicability of the DMLaplace MP2 method to linear-scaling SCF techniques [13,14] constructing approximate DMs. In linear-scaling calculations, DMs will include unexpected errors. If the DM-Laplace MP2 is combined with these techniques, this noise may seriously affect the obtained correlation energy. We investigated the effect of the noise upon the correlation energy by performing the DM-Laplace MP2 calculation with DMs involving artificial random errors. Table 3 shows mean errors (MEs) and mean absolute errors (MAEs) of the HF energies EHF, the MP2 correlation energies Eð2Þ corr , and the total MP2 energies EMP2 of 50 samples of DMs of a benzene molecule including at maximum 10% random noise with the 6-31G basis set. The difference of EHF arises from the difference of random samples. Only the Chebyshev (PF) scheme is a robust method for giving correlation energies when the noisy DMs are used. Table 4 shows MEs and MAEs of the MP2 energies of 50 samples of DMs of a benzene molecule with the 631G basis set including 10%, 1%, 0.1% and 0.01% random noise at maximum, respectively. The Chebyshev scheme with PF formulae was adopted for matrix exponentials. 23 MAEs of the correlation energies Eð2Þ corr are about 10 times as small as those of HF energies EHF. This ratio
254
M. Kobayashi, H. Nakai / Chemical Physics Letters 420 (2006) 250–255
Table 3 Dependence of HF and MP2 energy errors with noisy DMs on energy-weighted DM evaluating methods (benzene, 6-31G) [hartree] Method
ME Pade´ (S1F) Pade´ (PF) Chebyshev (S1F) Chebyshev (PF) No discrepancy a
Eð2Þ corr
EHF +0.182331 0.059809 0.225462 +0.085511
EMP2
MAE
ME
MAE
ME
MAE
2.101904 1.810934 2.466067 2.089454
–a –a 11.976461 0.016030
–a –a 11.976461 0.016518
–a –a 11.233509 +0.069480
–a –a 11.233509 2.086687
230.623696
0.523123
231.146818
Energy-weighted DM diverged for a large s value.
Table 4 Effect of DM discrepancy on HF and MP2 energies of benzene (6-31G) [hartree]a DM discrepancy
ME 10% 1% 0.1% 0.01% No discrepancy a
Eð2Þ corr
EHF +0.085511 +0.040437 0.008195 0.000059
EMP2
MAE
ME
MAE
2.089454 0.207304 0.023166 0.001974
0.016030 0.000183 +0.000061 +0.000058
0.016518 0.000615 0.000083 0.000058
230.623696
0.523123
ME 0.069480 0.040254 0.008134 0.000002
MAE 2.086687 0.206971 0.023116 0.001973
231.146818
The Chebyshev (PF) scheme was adopted for the evaluation of energy-weighted DMs.
corresponds to that of actual HF energy and correlation energy. The noise included in the DM affects the HF energy much more strongly than the correlation energy. By combining with linear-scaling HF techniques providing accurate HF energy, the DM-Laplace MP2 scheme has a great potential for the efficient evaluation of correlation energy. 5. Concluding remarks In this Letter, we have assessed the numerical calculation of MP2 correlation energy by means of Surja´n’s DM-based method, referred to as DM-Laplace MP2. The DM-Laplace scheme can provide MP2 correlation energy with reasonable (milli-hartree or micro-hartree) accuracy by using an appropriate numerical integration algorithm. Since the Gaussian quadrature method by means of a least-square fitting of the 1/x function over all orbital energies should not be applied to the DM-Laplace MP2 method, we adopted three schemes: the Gauss–Laguerre, the Euler–Maclaurin, and the Romberg quadratures. Milli-hartree accuracy can be achieved by the Euler–Maclaurin and the Romberg scheme with small quadrature points. In the DM-Laplace MP2 calculation, the evaluation of matrix exponentials is required to construct the energy-weighted DMs. It was found that the Chebyshev expansion scheme using PF formulae is stable even if noise exists in the reference DM. The influence of DM errors upon the correlation energy is remarkably smaller than upon HF energy. The DM-Laplace MP2 opened up new possibility of large-scale correlated MO calculations by the combination of the linear-scaling HF [13,14] and linear-scaling AOLaplace MP2 [11] calculations. For that purpose, sufficient investigation of the behavior of approximate DMs
obtained from linear-scaling SCF techniques is important. In addition, linear-scaling calculation of matrix exponentials is required. Nevertheless, the latter problem should be simple to solve since a DM and a Fock are both sparse matrices. We are now investigating the former problem and will report on it in the following study. Furthermore, the derivation of the energy derivatives in the DM-Laplace MP2 is important for geometry optimizations and frequency calculations of massive systems. Extension of AO-Laplace MP2 derivatives derived by Ha¨ser [8] to DM-Laplace scheme is strongly expected. Acknowledgements We thank Mr. Michio Katouda for valuable discussions about the Laplace MP2. This work was supported by a Grant-in-Aid for Exploratory Research ‘KAKENHI 16655010’ from the Japanese Ministry of Education, Culture, Sports, Science and Technology (MEXT), by a NAREGI Nano Science Project from MEXT, and by the 21st Century Center Of Excellence (21COE) ‘Practical NanoChemistry’ from MEXT. The calculations were performed in part at the Research Center for Computational Science (RCCS) of National Institutes of Natural Sciences. References [1] C. Møller, M.S. Plesset, Phys. Rev. 46 (1934) 618. [2] P. Pulay, S. Saebø, Theor. Chim. Acta 69 (1986) 357. [3] M. Schu¨tz, G. Hetzer, H.-J. Werner, J. Chem. Phys. 111 (1999) 5691. [4] G. Hetzer, M. Schu¨tz, H. Stoll, H.-J. Werner, J. Chem. Phys. 113 (2000) 9443. [5] M. Feyereisen, G. Fitzgerald, A. Komornicki, Chem. Phys. Lett. 208 (1993) 359.
M. Kobayashi, H. Nakai / Chemical Physics Letters 420 (2006) 250–255 [6] [7] [8] [9] [10] [11] [12] [13]
J. Almlo¨f, Chem. Phys. Lett. 181 (1991) 319. M. Ha¨ser, J. Almlo¨f, J. Chem. Phys. 96 (1992) 489. M. Ha¨ser, Theor. Chim. Acta 87 (1993) 147. G. Rauhut, P. Pulay, Chem. Phys. Lett. 248 (1996) 223. A.K. Wilson, J. Almlo¨f, Theor. Chim. Acta 95 (1997) 49. P.Y. Ayala, G.E. Scuseria, J. Chem. Phys. 110 (1999) 3660. P.R. Surja´n, Chem. Phys. Lett. 406 (2005) 318. S. Goedecker, Rev. Mod. Phys. 71 (1999) 1085, and references therein. [14] S.Y. Wu, C.S. Jayanthi, Phys. Rep. 358 (2002) 1, and references therein. [15] M.W. Schmidt, K.K. Baldridge, J.A. Boatz, S.T. Elbert, M.S. Gordon, J.J. Jensen, S. Koseki, N. Matsunaga, K.A. Nguyen,
[16]
[17] [18] [19] [20] [21]
255
S. Su, T.L. Windus, M. Dupuis, J.A. Montgomery, J. Comput. Chem. 14 (1993) 1347. W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in FORTRAN 77, second edn., Cambridge University Press, Cambridge, England, 1992. W.J. Hehre, R. Ditchfield, J.A. Pople, J. Chem. Phys. 56 (1972) 2257. P.J. Davis, P. Rabinowitz, Methods of Numerical Integration, second edn., Academic Press, FL, USA, 1984. C.W. Murray, N.C. Handy, G.J. Laming, Mol. Phys. 78 (1993) 997. P.C. Hariharan, J.A. Pople, Theor. Chim. Acta 28 (1973) 213. R.B. Sidje, ACM Trans. Math. Software 24 (1998) 130.
Implementation of Surja´n’s density matrix formulae for calculating second-order Møller–Plesset energy Masato Kobayashi, Hiromi Nakai
*
Department of Chemistry, School of Science and Engineering, Waseda University, 3-4-1 Okubo, Shinjuku-ku, Tokyo 169-8555, Japan Received 20 October 2005 Available online 20 January 2006
Abstract We numerically assess the method for obtaining second-order Møller–Plesset (MP2) energy from the Hartree–Fock density matrix (DM) recently proposed by Surja´n [Surja´n, Chem. Phys. Lett. 406 (2005) 318]. It is confirmed that Surja´n’s method, referred to as DM-Laplace MP2, can obtain MP2 energy accurately by means of appropriate integral quadrature and a matrix exponential evaluation scheme. Numerical tests reveal that the Euler–Maclaurin and the Romberg numerical integration schemes can achieve milli-hartree accuracy with small quadrature points. This Letter also indicates the possibility of the application of DM-Laplace MP2 to linear-scaling selfconsistent field techniques, which give approximate DM. Ó 2006 Elsevier B.V. All rights reserved.
1. Introduction Second-order Møller–Plesset (MP2) perturbation theory [1] is the simplest molecular orbital (MO) method that can deal with electron correlation. MP2 is also the most practical method enabling large-scale MO calculations including correlation effects. So far, many researchers have proposed various techniques for accelerating the evaluation of MP2 energy in terms of the localized orbital [2–4], the resolution of the identity (RI) technique for integral evaluations [5], and the Laplace transformation of the energy denominator [6–11]. The closed-shell MP2 correlation energy is given by using spatial occupied orbitals i, j, . . ., virtual orbitals a, b, . . ., andRRstandard electron repulsion integral notation of ðiajjbÞ ¼ iðr1 Þjðr2 Þr1 12 aðr1 Þbðr2 Þ dr1 dr2 as follows: ð2Þ Ecorr
occ X vir X ðiajjbÞ½2ðiajjbÞ ðibjjaÞ ¼ ; ea þ eb ei ej i;j a;b
ð1Þ
where ei and ea are the orbital energies of occupied and virtual orbitals, respectively. In the recent Letter by Surja´n [12], *
Corresponding author. Fax: +81 3 3205 2504. E-mail address: [email protected] (H. Nakai).
0009-2614/$ - see front matter Ó 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2005.12.088
the explicit MP2 correlation energy functional Eð2Þ corr ½P, where the argument function P is the Hartree–Fock (HF) density matrix (DM), was presented. Surja´n’s strategy is based on so-called atomic orbital (AO) Laplace MP2 form [8,11]. Eq. (1) can be rewritten by means of the Laplace transformation as Z 1 s X ð2Þ Ecorr ¼ e2 ðsÞ ds wa e2 ðsa Þ ð2Þ 0
a¼1
with s quadrature points {sa} and their weights {wa}. The integrand of Eq. (2) is given in AO basis by XX e2 ðsa Þ ¼ X alc Y amd X ajk Y aer ðcdjjeÞ½2ðlmjkrÞ ðlrjkmÞ cdje lmkr
¼
X
ðlmjk rÞ½2ðlmjkrÞ ðlrjkmÞ;
ð3Þ
lmkr
where Xa and Ya are the energy-weighted DMs of electron and hole, respectively, and are given by Xa ¼
occ X
eei sa Ci CTi ;
ð4Þ
eea sa Ca CTa .
ð5Þ
i
Ya ¼
vir X a
M. Kobayashi, H. Nakai / Chemical Physics Letters 420 (2006) 250–255
According to Surja´n [12], the energy-weighted DMs of Eqs. (4) and (5) can be represented by means of AO-based DMs instead of orbital energies as follows: X a ¼ e sa S a
Y ¼e
1
sa S
F
P;
1
F
Q;
ð6Þ ð7Þ
where S is an overlap matrix, F is a Fock matrix, and P and Q are the HF DMs of electron and hole in AO basis, P¼
occ X
Ci CTi ;
ð8Þ
Ca CTa ;
ð9Þ
i
Q¼
vir X a
respectively. Note that we only assume the Hartree–Fock– Roothaan equation FCi = eiSCi so far. Furthermore, an inverse overlap matrix can be eliminated by assuming the commutation relationship of FPS = SPF and idempotency of PS. Finally, the following working formulae are obtained: Xa ¼ esa PF P; a
Y ¼e
sa QF
Q.
ð10Þ ð11Þ
Therefore, MP2 energy can be obtained as the functional of DMs of electron P and hole Q. This method for evaluating MP2 energy is referred to as DM-Laplace MP2. It should be noted that these DMs satisfy P þ Q ¼ S1 .
ð12Þ
DM-Laplace MP2 method has prospects of enabling large-scale MP2 calculation. DM-Laplace MP2 requires only HF DM which can often be constructed by linearscaling SCF methods [13,14]. Furthermore, the method is based on the AO-Laplace MP2 scheme whose computational cost can also be reduced to scaling linearly depending on the system size and is well assessed by Ayala and Scuseria [11]. Although Surja´n suggested the applicability of the DM-Laplace MP2 method, numerical results have not been shown yet. In this Letter, the DM-Laplace MP2 calculation code was implemented to the GAMESS quantum chemistry package [15] and numerical assessments of the method were carried out. We indicate two problem areas in the DM-Laplace MP2 calculation and provide feasible solutions. We first attend to the numerical quadrature of Eq. (2) in Section 2. Next, a method of evaluating energy-weighted DMs X and Y is discussed in Section 3. In addition, in Section 4, we examine whether the DM-Laplace MP2 can apply to a method for constructing approximate DMs. Finally, concluding remarks are presented in Section 5.
251
gested the efficiency of a least-square fit of the 1/x function over Dijab ! s X X 1 wa expðDijab sa Þ ¼ min !; ð13Þ Dijab a¼1 ijab where Dijab ¼ ea þ eb ei ej .
ð14Þ
This scheme provides quite accurate results even with a small number of quadrature points. For instance, according to Ha¨ser and Almlo¨f [7], micro-hartree accuracy can be achieved with the number of quadrature points s = 8– 10. This scheme, however, is not suited for the DM-Laplace MP2 calculation because the DM-Laplace MP2 was proposed as a method for evaluating MP2 energy without the use of orbital energies, which appear in Eq. (14). Wilson and Almlo¨f [10] have reported the exponentially decaying behavior of the integrand of Eq. (2), i.e., e2(s) with the value of s. The Gauss–Laguerre quadrature scheme is suitable for the semi-infinite integral whose integrand can be approximated as the product of polynomials and an exponential function [16]. Fig. 1 shows the decay of e2(s) with respect to s in the Laplace-MP2 calculation of a benzene molecule with the 6-31G basis set [17]. In Fig. 1, the Gauss–Laguerre quadrature points corresponding to s = 10 are also shown by small squares. Since e2(s) is a quasi single exponential function, as seen from Fig. 1, it is expected that the Gauss–Laguerre quadrature is suited for the numerical integration of Eq. (2). On the other hand, Ayala and Scuseria [11] have suggested the efficiency of using Euler–Maclaurin summation [18,19] for this numerical quadrature, with a proper change of variable (CoV) from s to r: Z 1 Z 1 Z 1 ds ð2Þ Ecorr ¼ e2 ðsÞ ds ¼ e2 ðsÞ dr ¼ f2 ðrÞ dr dr 0 0 0 " # s X 1 k 1 f2 ð15Þ þ ðf2 ð0Þ þ f2 ð1ÞÞ . s þ 1 k¼1 sþ1 2
2. Numerical quadrature In the Laplace MP2 calculations, numerical quadrature points {sa} and weight parameters {wa} have to be determined in some fashion. Several authors [7,10,11] have sug-
Fig. 1. Decay behavior of integrand e2(s) for benzene (6-31G). (h) Gauss– Laguerre quadrature points for s = 10.
252
M. Kobayashi, H. Nakai / Chemical Physics Letters 420 (2006) 250–255
The integration error dEð2Þ corr of the above formula can be evaluated as 2k h 1 i X B2k 1 ð2k1Þ ð2k1Þ ð2Þ ¼ f2 ð1Þ f2 ð0Þ dEcorr ð2kÞ! s þ 1 k¼1 ¼
f20 ð0Þ 12ðs þ 1Þ
2
ð5Þ
f2000 ð0Þ
f2 ð0Þ þ 720ðs þ 1Þ4 30240ðs þ 1Þ6
ð7Þ
f2 ð0Þ 1209600ðs þ 1Þ
8
þ
ð16Þ
because e2(s) and its all higher derivatives are zero at s = 1, where B2k are Bernoulli numbers. Although arbitrary CoV choice is possible in this scheme, it is best to choose one for which its Jacobian ds/dr and as large a number of its derivatives as possible are zero for r = 0. Here, we apply two formulae for the CoV; type (A) is s¼
r2 ðr3 þ r4 þ r5 þ 0:9r6 Þ=4 2
ð17Þ
þ r2 tanðpr=2Þ.
ð18Þ
ð1 rÞ
where T nð0Þ is the approximate value of the integral by the trapezoid method with 2n division " # s X 1 k 1 ð0Þ Tn ¼ f2 ð20Þ þ ðf2 ð0Þ þ f2 ð1ÞÞ ; s þ 1 k¼1 sþ1 2 with s = 2n 1. The explicit formulae for n = 3 (s = 7) and n = 4 (s = 15) are 2 1 1 3 ð3Þ 256f 2 T3 ¼ þ 88f 2 þ 256f 2 2835 8 4 8 1 5 3 7 þ109f 2 þ 256f 2 þ 88f 2 þ 256f 2 2 8 4 8 ð21Þ and ð4Þ T4
and type (B) s¼
r3 0:9r4 ð1 rÞ
2
Their Jacobians are both zero, and the first derivative of type (B)’s Jacobian also becomes zero for r = 0. We additionally tested the Romberg quadrature method [16,18], which is the extrapolation of some results of the extended trapezoid method based on the idea of Richardson’s deferred approach to the limit, with the same CoV, ð2Þ i.e., Eqs. (17) and (18). Specifically, Ecorr can be approxiðnÞ mated to T n by means of the following recurrence relation: ðk1Þ
T ðkÞ n ¼
4k T ðk1Þ T n1 n ; k 4 1
ð19Þ
2 1 1 32768f 2 ¼ þ 11008f 2 722925 16 8 3 1 5 þ32768f 2 þ13864f 2 þ 32768f 2 16 4 16 3 7 1 þ 11008f 2 þ 32768f 2 þ 13779f 2 8 16 2 9 5 11 þ 32768f 2 þ 11008f 2 þ 32768f 2 16 8 16 3 13 þ 13864f 2 þ32768f 2 4 16 7 15 þ 11008f 2 þ 32768f 2 ; ð22Þ 8 16
respectively. We have implemented the AO-Laplace MP2 and the DM-Laplace MP2 calculation code in the GAMESS quantum chemistry package [15] based on Ha¨ser’s quadratic
Table 1 Dependence of MP2 correlation energy of benzene on quadrature method [hartree] Method
s
6-31G
6-31G*
Eð2Þ corr
(Diff.)
Eð2Þ corr
(Diff.)
Gauss–Laguerre
5 10 15 50
0.512724 0.519626 0.521514 0.523061
(+0.010345) (+0.003444) (+0.001556) (+0.000009)
0.733451 0.766420 0.774904 0.784609
(+0.051311) (+0.018342) (+0.009858) (+0.000152)
Euler–Maclaurin (A)
5 10 15 50
0.516278 0.521006 0.522094 0.522972
(+0.006791) (+0.002064) (+0.000976) (+0.000098)
0.770241 0.780406 0.782704 0.784557
(+0.014521) (+0.004355) (+0.002057) (+0.000205)
Euler–Maclaurin (B)
5 10 15 50
0.522936 0.523065 0.523069 0.523068
(+0.000133) (+0.000004) (+0.000001) (+0.000002)
0.784540 0.784761 0.784760 0.784758
(+0.000221) (+0.000000) (+0.000001) (+0.000003)
Romberg (A)
7 15
0.522629 0.523123
(+0.000441) (0.000053)
0.784643 0.784666
(+0.000118) (+0.000096)
Romberg (B)
7 15
0.523683 0.523227
(0.000613) (0.000157)
0.783803 0.784919
(+0.000958) (0.000157)
Canonical
–
0.523070
0.784761
M. Kobayashi, H. Nakai / Chemical Physics Letters 420 (2006) 250–255
scaling algorithm [8]. In this algorithm, six steps of integral cut-off using Schwarz’s inequality are adopted (see Refs. [8,11] for details). We applied Schwarz’s inequality threshold of h = 1011 a.u. in all the following calculations. Table 1 shows the all-electron MP2 correlation energy of a benzene molecule obtained from the DM-Laplace scheme by several quadrature methods. The Chebyshev (PF) scheme was used for the evaluation of energyweighted DMs, which is explained in detail in Section 3. The 6-31G [17] and 6-31G* [20] basis sets were adopted. ‘(A)’ and ‘(B)’ on the quadrature scheme mean the CoV. The discrepancy with respect to the energy obtained from conventional canonical MP2 calculation, which is expressed in the bottom row of the table, is also shown in parentheses. The correlation energy of each quadrature method approaches the canonical one as the number of quadrature points s increases. The energies on the Gauss–Laguerre scheme with s = 15 have errors of about 1.6 and 9.9 milli-hartrees with the 6-31G and 6-31G* basis, respectively. The Euler–Maclaurin and the Romberg schemes, however, achieve milli-hartree accuracy even with 5–7 quadrature points if a suitable CoV is adopted. In the Euler–Maclaurin quadrature, type (B) CoV gives a considerably more precise estimate of energy (attaining micro-hartree accuracy with s = 10) than type (A) because the Jacobian as well as its first derivative become zero for r = 0 and the first term of the r.h.s. of Eq. (16) vanishes in type (B) CoV. On the other hand, the choice of CoV makes little difference in the Romberg method since the method gives the extrapolation to the Euler–Maclaurin summation. The Euler–Maclaurin scheme gives the upper bound of the canonical MP2 energy. This should be true in any case in type (A) CoV since the first term of the r.h.s. of Eq. (16) remains as f20 ð0Þ ¼ 2e2 ð0Þ with e2(0) > 0. In type (B) CoV, however, this behavior does not necessarily satisfy and it depends on the shape of e2(s). 3. Evaluation of energy-weighted density matrix We have proposed two working formulae for the evaluation of energy-weighted DMs X and Y in the DM-Laplace MP2. Eqs. (6) and (7) referred to as S1F formulae require the inverse overlap matrix, while Eqs. (10) and (11) referred to as PF formulae do not. In principle, both formulae should give the same result if a Fock F and a DM P satisfy the commutation relationship of FPS = SPF, which is fulfilled in the self-consistent field (SCF) solution. Whichever formula is used for energy-weighted DMs, the calculation of a matrix exponential is necessary. A number of methods for evaluating a matrix exponential have been presented. We apply two rational approximation algorithms for exp(x) function, that is, (6,6)-degree Pade´ and (14,14)-degree Chebyshev expansion using the EXPOKIT software package [21]. In the practical calculation of matrix exponentials using the Pade´ or the Chebyshev approximation, two working equations, i.e., S1F formulae and PF formulae, may no longer provide the same
253
Table 2 Dependence of MP2 correlation energy of benzene on energy-weighted DM evaluating method [hartree] Method
6-31G
6-31G*
Pade´ (S1F) Pade´ (PF) Chebyshev (S1F) Chebyshev (PF)
–a 0.5230654543 0.5230654543 0.5230654543
–a 0.7847609764 0.7847609764 0.7847609764
AO-Laplace
0.5230654543
0.7847609764
a
Energy-weighted DM diverged for a large s value.
matrix due to the character of the algorithm and/or numerical stability. Table 2 shows the DM-Laplace MP2 correlation energy of a benzene molecule by using the Pade´ and the Chebyshev approximation for the evaluation of matrix exponentials with the 6-31G and 6-31G* basis sets. We used the Euler–Maclaurin (B) scheme with s = 10. The results in the AO-Laplace scheme, in which MOs and orbital energies are used to construct energy-weighted DMs by Eqs. (4) and (5), are also shown. Except for the Pade´ approximant case using S1F formulae, where the energy-weighted DM diverges for a large s value, correlation energies obtained from the DM-Laplace MP2 calculation completely accord with those from the AO-Laplace MP2 with 1010 hartree accuracy. It is found that the present DM-Laplace MP2 scheme is equivalent both in theory and in practice to the AO-Laplace MP2 as long as an appropriate method for evaluating the energy-weighted DM is chosen. 4. DM-Laplace MP2 calculation using approximate density matrix Surja´n [12] has suggested the applicability of the DMLaplace MP2 method to linear-scaling SCF techniques [13,14] constructing approximate DMs. In linear-scaling calculations, DMs will include unexpected errors. If the DM-Laplace MP2 is combined with these techniques, this noise may seriously affect the obtained correlation energy. We investigated the effect of the noise upon the correlation energy by performing the DM-Laplace MP2 calculation with DMs involving artificial random errors. Table 3 shows mean errors (MEs) and mean absolute errors (MAEs) of the HF energies EHF, the MP2 correlation energies Eð2Þ corr , and the total MP2 energies EMP2 of 50 samples of DMs of a benzene molecule including at maximum 10% random noise with the 6-31G basis set. The difference of EHF arises from the difference of random samples. Only the Chebyshev (PF) scheme is a robust method for giving correlation energies when the noisy DMs are used. Table 4 shows MEs and MAEs of the MP2 energies of 50 samples of DMs of a benzene molecule with the 631G basis set including 10%, 1%, 0.1% and 0.01% random noise at maximum, respectively. The Chebyshev scheme with PF formulae was adopted for matrix exponentials. 23 MAEs of the correlation energies Eð2Þ corr are about 10 times as small as those of HF energies EHF. This ratio
254
M. Kobayashi, H. Nakai / Chemical Physics Letters 420 (2006) 250–255
Table 3 Dependence of HF and MP2 energy errors with noisy DMs on energy-weighted DM evaluating methods (benzene, 6-31G) [hartree] Method
ME Pade´ (S1F) Pade´ (PF) Chebyshev (S1F) Chebyshev (PF) No discrepancy a
Eð2Þ corr
EHF +0.182331 0.059809 0.225462 +0.085511
EMP2
MAE
ME
MAE
ME
MAE
2.101904 1.810934 2.466067 2.089454
–a –a 11.976461 0.016030
–a –a 11.976461 0.016518
–a –a 11.233509 +0.069480
–a –a 11.233509 2.086687
230.623696
0.523123
231.146818
Energy-weighted DM diverged for a large s value.
Table 4 Effect of DM discrepancy on HF and MP2 energies of benzene (6-31G) [hartree]a DM discrepancy
ME 10% 1% 0.1% 0.01% No discrepancy a
Eð2Þ corr
EHF +0.085511 +0.040437 0.008195 0.000059
EMP2
MAE
ME
MAE
2.089454 0.207304 0.023166 0.001974
0.016030 0.000183 +0.000061 +0.000058
0.016518 0.000615 0.000083 0.000058
230.623696
0.523123
ME 0.069480 0.040254 0.008134 0.000002
MAE 2.086687 0.206971 0.023116 0.001973
231.146818
The Chebyshev (PF) scheme was adopted for the evaluation of energy-weighted DMs.
corresponds to that of actual HF energy and correlation energy. The noise included in the DM affects the HF energy much more strongly than the correlation energy. By combining with linear-scaling HF techniques providing accurate HF energy, the DM-Laplace MP2 scheme has a great potential for the efficient evaluation of correlation energy. 5. Concluding remarks In this Letter, we have assessed the numerical calculation of MP2 correlation energy by means of Surja´n’s DM-based method, referred to as DM-Laplace MP2. The DM-Laplace scheme can provide MP2 correlation energy with reasonable (milli-hartree or micro-hartree) accuracy by using an appropriate numerical integration algorithm. Since the Gaussian quadrature method by means of a least-square fitting of the 1/x function over all orbital energies should not be applied to the DM-Laplace MP2 method, we adopted three schemes: the Gauss–Laguerre, the Euler–Maclaurin, and the Romberg quadratures. Milli-hartree accuracy can be achieved by the Euler–Maclaurin and the Romberg scheme with small quadrature points. In the DM-Laplace MP2 calculation, the evaluation of matrix exponentials is required to construct the energy-weighted DMs. It was found that the Chebyshev expansion scheme using PF formulae is stable even if noise exists in the reference DM. The influence of DM errors upon the correlation energy is remarkably smaller than upon HF energy. The DM-Laplace MP2 opened up new possibility of large-scale correlated MO calculations by the combination of the linear-scaling HF [13,14] and linear-scaling AOLaplace MP2 [11] calculations. For that purpose, sufficient investigation of the behavior of approximate DMs
obtained from linear-scaling SCF techniques is important. In addition, linear-scaling calculation of matrix exponentials is required. Nevertheless, the latter problem should be simple to solve since a DM and a Fock are both sparse matrices. We are now investigating the former problem and will report on it in the following study. Furthermore, the derivation of the energy derivatives in the DM-Laplace MP2 is important for geometry optimizations and frequency calculations of massive systems. Extension of AO-Laplace MP2 derivatives derived by Ha¨ser [8] to DM-Laplace scheme is strongly expected. Acknowledgements We thank Mr. Michio Katouda for valuable discussions about the Laplace MP2. This work was supported by a Grant-in-Aid for Exploratory Research ‘KAKENHI 16655010’ from the Japanese Ministry of Education, Culture, Sports, Science and Technology (MEXT), by a NAREGI Nano Science Project from MEXT, and by the 21st Century Center Of Excellence (21COE) ‘Practical NanoChemistry’ from MEXT. The calculations were performed in part at the Research Center for Computational Science (RCCS) of National Institutes of Natural Sciences. References [1] C. Møller, M.S. Plesset, Phys. Rev. 46 (1934) 618. [2] P. Pulay, S. Saebø, Theor. Chim. Acta 69 (1986) 357. [3] M. Schu¨tz, G. Hetzer, H.-J. Werner, J. Chem. Phys. 111 (1999) 5691. [4] G. Hetzer, M. Schu¨tz, H. Stoll, H.-J. Werner, J. Chem. Phys. 113 (2000) 9443. [5] M. Feyereisen, G. Fitzgerald, A. Komornicki, Chem. Phys. Lett. 208 (1993) 359.
M. Kobayashi, H. Nakai / Chemical Physics Letters 420 (2006) 250–255 [6] [7] [8] [9] [10] [11] [12] [13]
J. Almlo¨f, Chem. Phys. Lett. 181 (1991) 319. M. Ha¨ser, J. Almlo¨f, J. Chem. Phys. 96 (1992) 489. M. Ha¨ser, Theor. Chim. Acta 87 (1993) 147. G. Rauhut, P. Pulay, Chem. Phys. Lett. 248 (1996) 223. A.K. Wilson, J. Almlo¨f, Theor. Chim. Acta 95 (1997) 49. P.Y. Ayala, G.E. Scuseria, J. Chem. Phys. 110 (1999) 3660. P.R. Surja´n, Chem. Phys. Lett. 406 (2005) 318. S. Goedecker, Rev. Mod. Phys. 71 (1999) 1085, and references therein. [14] S.Y. Wu, C.S. Jayanthi, Phys. Rep. 358 (2002) 1, and references therein. [15] M.W. Schmidt, K.K. Baldridge, J.A. Boatz, S.T. Elbert, M.S. Gordon, J.J. Jensen, S. Koseki, N. Matsunaga, K.A. Nguyen,
[16]
[17] [18] [19] [20] [21]
255
S. Su, T.L. Windus, M. Dupuis, J.A. Montgomery, J. Comput. Chem. 14 (1993) 1347. W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in FORTRAN 77, second edn., Cambridge University Press, Cambridge, England, 1992. W.J. Hehre, R. Ditchfield, J.A. Pople, J. Chem. Phys. 56 (1972) 2257. P.J. Davis, P. Rabinowitz, Methods of Numerical Integration, second edn., Academic Press, FL, USA, 1984. C.W. Murray, N.C. Handy, G.J. Laming, Mol. Phys. 78 (1993) 997. P.C. Hariharan, J.A. Pople, Theor. Chim. Acta 28 (1973) 213. R.B. Sidje, ACM Trans. Math. Software 24 (1998) 130.