Local MP2 periodic study of rare-gas crystals

Chemical Physics Letters 467 (2009) 294–298

Contents lists available at ScienceDirect

Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett

Local MP2 periodic study of rare-gas crystals Migen Halo, Silvia Casassa, Lorenzo Maschio, Cesare Pisani * Dipartimento di Chimica IFM and Centre of Excellence NIS (Nanostructured Interfaces and Surfaces), Università di Torino, via P. Giuria 5, I-10125 Torino, Italy

a r t i c l e

i n f o

Article history: Received 27 September 2008 In final form 13 November 2008 Available online 21 November 2008

a b s t r a c t The CRYSCOR code, a computational tool which solves the second-order Møller–Plesset local equations for periodic systems, is applied to fcc rare-gas crystals Ne through Xe. Cohesive energies, equilibrium lattice parameters and bulk moduli are calculated; DFT calculations have been performed for purposes of comparison. Ó 2008 Elsevier B.V. All rights reserved.

1. Introduction Rare-gas crystals have aroused lively interest in recent years, as it is shown by numerous studies either in the frame of a cluster [1–3] or a periodic approach [4–8]. An accurate treatment of electron correlation is mandatory in order to correctly describe these prototypical Lennard–Jones systems: they thus represent an excellent test for post-Hartree–Fock (HF) methods, since both HF and density functional theory (DFT) are expected to fail in this case, the former because it does not account for electron correlation, the latter because of its poor description of weak intermolecular interactions, despite recent improvements [5]. In this Letter we present an ab-initio periodic treatment of rare gases Ne through Xe in their face-centered-cubic crystalline structure, by adopting a local-correlation approach as recently implemented in the CRYSCOR code [9]. CRYSCOR results from the confluence of knowledge and technology of two reference codes, CRYSTAL [10] and MOLPRO [11]. The local-correlation scheme, originally proposed for molecules by Pulay and Saebø [12,13] and implemented in MOLPRO, is reformulated in CRYSCOR for crystalline systems, which implies, among other features, the full exploitation of point and translational symmetry; the local approach, on the other hand, allows us to take advantage of the short-range nature of dynamical electron correlation. The reference HF solution is provided by the well-known periodic LCAO CRYSTAL code, to which CRYSCOR is strongly linked so as to represent in a sense a post-HF option of it. All these codes adopt a basis set of atomic orbitals (AO) which are in turn expressed as a linear combination of Gaussian type orbitals (GTO). Currently, only the second-order Møller–Plesset (MP2) treatment is implemented, therefore, only non-conducting systems can be treated. Since MP2 is the lowest level of post-HF perturbation theory, the limitations of the method are clear. The aim of this work is indeed not to provide data in good agreement with the experiment, but rather to explore the * Corresponding author. E-mail address: [email protected] (C. Pisani). 0009-2614/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2008.11.043

capabilities and limits of the program, in the framework of a robust and well-defined approach. The present study exploits the accurate exploration of the role of the various computational parameters performed in our recent papers concerning crystalline Argon [7] and Neon [8], respectively; the comprehensive treatment of the four rare gases permits us here to critically analyze their different behavior as concerns the adequacy of the MP2 approximation. DFT calculations using the PWGGA [14] and B3LYP [15] functionals have also been performed for purposes of comparison. The Letter is organized as follows: In Section 2, the main computational aspects are briefly illustrated; in Section 3, the calculated cohesive energies, equilibrium lattice parameters and bulk moduli are reported and commented on. 2. Computational details In this section we briefly comment on the main computational choices: basis set, correction of the basis-set-superposition-error (BSSE) and local-correlation parameters. We will not deal with the theoretical background of the method, since it has been thoroughly described elsewhere [16–19], nor justify in detail the computational settings adopted, as an exhaustive analysis of their role has been carried out in previous papers [7,8]. 2.1. Basis set The basis set (BS) is a delicate issue because the MP2 solution is particularly sensitive to it. In fact, to correctly describe the virtual space, the use of diffuse high-angular-momentum functions, not needed for the HF solution, becomes mandatory. Our reference BSs are the valence ones with a semi-empirical effective-core relativistic pseudopotential optimized by Nicklass et al. [21] and used as a starting point by Roscis´zewski and coworkers [1]. These sets, referred to in the following as A, are ð6s6p3d1f Þ=½4s4p3d1f  for Ar, Kr, Xe and ð7s7p3d1f Þ=½4s4p3d1f  for Ne. We could not further extend BS A to yield BSs B,C employed in Ref. [1] because they

295

M. Halo et al. / Chemical Physics Letters 467 (2009) 294–298

2.2. Basis-set-superposition-error BSSE is small but not completely negligible at the HF level, and is quite relevant at the MP2 level, therefore all interaction energies have been BSSE corrected. The standard counterpoise (CP) technique [23] has been used to correct the HF, DFT and hybrid-DFT cohesive energies, whereas for MP2 an alternative approach can be adopted which exploits the intrinsic additivity of the MP2 contributions. Practically, the energy due to intra-atomic correlation is subtracted from the total MP2 correlation energy, thus leaving only inter-atomic contributions. This simple technique, referred to as partition method (PM), has been shown to provide an MP2 energy which is almost BSSE-free, at least for such simple cases as the present ones [8]. In the following, we will essentially refer to this technique; however, a comparison between CP and PM results will be carried out in the next section. 2.3. Local-correlation parameters The use of the locality Ansatz is controlled essentially by three computational parameters:

7 Ar Kr Ne Xe Limit

6

5

100*[E(D1)-Elim]/Ecoh

contain very diffuse s and p functions which are not only harmful in the periodic case (quasi-linear dependence) but even of scarce utility, since crystalline systems are much more compact objects than molecules. As a matter of fact, the most diffuse s and p function had to be eliminated from BS A of Xe to achieve convergence in the self-consistent-field procedure. For the MP2 case, we tried the usefulness of enriching the BS with additional d- and f-type functions. The former ones were found of practically no use. As for the latter, we added to the A BS (which includes a single f-type GTO for Ne, Ar, Kr and Xe, with exponent 2.58, 0.931, 0.705, 0.516, respectively), a second and a third f-type GTO (with exponents 0.30, 0.20, 0.15, 0.15; 0.60, 0.55, 0.40, 0.40, respectively, optimized so as to recover the maximum MP2 correlation energy), so obtaining the two BSs A0 and A00 : their performance is analyzed in the next section.

4

3

2

1

0

-1 1.4

1.6

1.8

2

2.2

2.4

2.6

D1/a0 Fig. 1. Long-range dispersive contributions for rare gases Ne through Xe. Filled symbols represent the difference between MP2 energies including all contributions up to the truncation distance between WF pairs, E(D1), and the limiting value, Elim ¼ Eð1Þ; empty symbols the same quantity, including the ‘Lennard–Jones’ extrapolated contributions to energy. For each gas, energy differences are reported as a percentage of cohesive energy, D1 values in units of the respective experimental lattice parameter.

a relative distance D2 generally less than D1; for the remaining pairs within D1 a multipolar approximation is used [19]: the D2 values here adopted correspond to treating as close-by WF pairs those within 4th-neighbor distance. 2.5. Computational times

(a) The threshold for the cutting of the tails of Wannier functions (WF) and of projected atomic orbitals (PAO). WFs and PAOs describe the occupied and the virtual HF manifold, respectively; they are well localized, but their tails may extend in principle to infinite distance. The coefficients in their AO expansion are set to zero if they are smaller than the input parameter T: the value here adopted (T ¼ 0:0001) is quite tight. (b) The size of the domain associated to each WF, i.e. the number of atoms whose PAOs define the local excitation space: here, they consist of 13 atoms, i.e. they include all first neighbors of the atom on which the WF is centered. (c) The distance up to which the correlation between two electrons is explicitly treated; it is controlled by the input D1 parameter, which sets the maximum distance allowed between two WFs involved in a pair excitation: here it is set to 12 Å (10 Å for Ne). The correlation contribution from pairs farther than D1 is extrapolated by exploiting the ar 6 long-range behavior of dispersive interactions, with a coefficients derived from a best-fit analysis of the terms explicitly calculated up to D1 [19]. Fig. 1 demonstrates the excellent performance of this technique for the four rare gases. 2.4. Two-electron-integral calculation DFP, i.e. the density fitting technique adapted to the periodic case [17–20], has been used for the calculation of all 2-electron repulsion integrals involving pairs of close-by WFs that is, up to

With the settings reported, a single-point calculation for Ne 0 00 takes 16 000, 28 000 and 41 000 CPU s using BSs A, A and A ; the times for Xe are about the same. Data are referred to a single-processor Opteron PC, 2.2 GHz, 4 GB of memory. 3. Results and discussion The four panels of Fig. 2 report binding-energy-versus-distance curves for the four systems with different Hamiltonians and BSs. In each case, energies are referred to the atomic energy calculated in the same approximation; the MP2 energies are BSSE corrected according to the PM technique explained in the previous section, and include the extrapolated long-range contribution. As a reference, the experimental data concerning binding energy and equilibrium lattice parameter (corrected for the zero-point vibrational contribution as in Ref. [2]) are also shown. Consider first the results obtained with the one-electron Hamiltonians; among the BSs considered, only A was used for them, since additional f functions were seen to be not influential in those cases. HF is seen to provide no binding at all. This is not what happens with PWGGA: a clear minimum is found in all cases at finite distances; however, the calculated cohesive energies seem to have no relation with the experimental ones, and the equilibrium lattice parameter is badly overestimated. B3LYP is repulsive at all distances; interestingly enough, the B3LYP curves intersect the HF ones at distances between a0 and 2a0 : this anomalous description of the long-range inter-atomic interaction is due to the incorrect

296

M. Halo et al. / Chemical Physics Letters 467 (2009) 294–298

2000

8000

Ne

1000

4000

0

0

-1000

-4000

-2000 4

4.4

4.8

5.2

8000

5.6

-8000 4.8

4000

0

0

-4000

-4000

5.6

6

6.4

5.2

5.6

6

8000

Kr

4000

-8000 5.2

Ar

6.8

-8000

6.4

Xe

6

6.4

6.8

7.2

Fig. 2. Cohesive energy (lHa) as a function of the lattice parameter (Å) for the four rare gases obtained with different methods and BSs. HF: plus; B3LYP: stars; PWGGA: diamonds; HF + MP2: squares for BS A, circles for BS A0 , triangles for BS A00 . All MP2 energies have been corrected for BSSE using the PM technique. The experimental value is reported as a black cross. Note that for Ne the y-axis range is different.

behavior at those very small electron densities of the local-exchange functional employed in the B3LYP Hamiltonian. Generally speaking, MP2 nicely corrects the HF curves, yielding a crude but reasonable estimate of cohesive energy and equilibrium lattice parameter. When comparing MP2 results for different BSs, the importance of the first additional f function (set A0 ) is evident, while the second one (set A00 ) does not change the results appreciably (Krypton is the most affected). This fact alone does not mean that we are close to the MP2 limit. It has been proved however in the case of the dimer that the MP2 limit is above the experimental value for Neon and below it for Argon [22], in accordance with the present results. We think therefore that our estimates of the MP2 binding energies for the rare gases cannot be changed very much owing to a better description of the virtual manifold. The results for cohesive energy, equilibrium lattice parameter and bulk modulus of the four systems presently obtained with different basis sets are reported in Table 1 along with the corresponding experimental values. A few theoretical results recently reported in the literature and obtained with a variety of

approaches are also shown for the sake of comparison. Those by Roscis´zewski et al. [2] have already been commented on and represent in a sense the most advanced treatment presently feasible: in the frame of an incremental approach, these Authors were able to perform CCSD(T) calculations (namely, coupled-cluster with single and double excitations, and perturbative treatment of triples) with very rich BSs (B,D). A very different approach was followed by Civalleri et al. [5], who combined the exact Hartree–Fock exchange with the semi-empirical Wilson–Levy correlation functional (HF–WL), which aims at describing long-range dispersive interactions; again, a very good BS (QVZP) was used. The results by Harl and Kresse [6] were obtained in the frame of the adiabatic connection fluctuation–dissipation theory (ACFDT), either with LDA or PBE, with a plane-wave (PW) BS. In the last two lines of the table, dimer data (calculated and experimental [34]) are reported. The following can be noted for the crystals: MP2 results with basis set A0 present a mean error of 23% in cohesive energy values, which is rather large, though expected due to the method itself; on the other hand, the predictive power of MP2 for lattice equilibrium

297

M. Halo et al. / Chemical Physics Letters 467 (2009) 294–298 Table 1 Calculated and experimental equilibrium lattice parameter (Å), cohesive energy (lHa) and bulk modulus (kbar) for rare-gas crystals Ne–Xe. Ne

Ar

Method

BS

a

Crystal HF + MP2 HF + MP2 HF + MP2

A A0 A00

4.45 4.35 4.35

B D QVZP PW PW

4.388 4.314 4.15 4.7 4.5 4.35

A0

3.10 3.10

CCSD(T) CCSD(T) HF–WL ACFDT–LDA ACFDT–PBE Experimental Dimer HF + MP2 Experimental

Ecoh 568.8 741.0 754.3 800.7 971.6 641.7 404.3 624.8 1002.0 102.5 133.8

Kr

Xe

B

a

Ecoh

B

a

Ecoh

B

a

Ecoh

B

4.7 6.0 6.1

5.30 5.20 5.20

3015.3 3850.0 3884.7

18.6 26.7 27.8

5.70 5.50 5.50

4015.7 5293.0 5623.0

23.3 30.0 32.5

6.20 6.00 6.10

5743.9 7694.2 7712.4

30.7 38.8 40.9

14.3 – – – – 10.9

5.354 5.284 5.08 5.4 5.3 5.23

2716.4 3043.4 2438.0 2168.3 3050.3 3268.0

26.9 – – – – 23.8

5. 763 5.670 – 5.8 5.7 5. 61

3559.0 4203. 8 – 3234.0 4116.0 4502.0

28.1 – – – – 36.1

6.239 6.137 – – – 6. 10

5372. 2 6060.0 – – – 6239.0

33.8 – – – – 36.4

– –

3.75 3.58

458.6 453.5

– –

4.05 4.01

634.7 637.1

– –

4.40 4.36

882.9 893.9

– –

The first three lines are referred to the present MP2 results; the subsequent lines to the CCCSD(T) results by Roscis´zewski et al. [2]; to the HF–WL results by Civalleri et al. [5]; to the ACFDT–LDA and ACFDT–PBE results by Harl and Kresse [6]. The references for the experimental data for crystals are as follows: [24–29] for Ne and Ar, [28,30,31] for Kr and [28,32,33] for Xe. The corrections for the zero-point vibrational energy (ZPE) are taken from Ref. [2]. In the last two lines, the corresponding data for dimers, both calculated and experimental [34], are reported. All MP2 energies have been corrected for BSSE using the PM technique.

parameter is quite satisfactory; results by means of HF–WL functional are poor for both Ne and Ar, whereas ACFDT ones, like ours, show a maximum bias for Ne. This pattern is roughly repeated for the binding energy of the dimers, in the sense that Ne2 is less accurately described with respect to Ar2 , Kr2 and Xe2 . As concerns bulk moduli, it is seen that calculated values in the MP2 approximation increase with increasing basis set quality, consistently with the fact that the corresponding cohesive energies become larger, the lattice constants shorter: with BS A00 , the calculated bulk modulus is larger than the experiment for Ar and Xe, while for Kr and Ne it is still below the experimental value. Again, the result for Ne is quite far from the reference, which seems to confirm the inadequacy of the MP2 approximation for that case. Note however that even the very accurate CCSD(T) results are closer to the experiment as concerns binding energy and lattice constant than bulk modulus.

0

0

A

0.2

E/Ecoh

All the results commented on so far were obtained through the PM correction to MP2 energies. Fig. 3 compares PM and CP results for the two different basis sets A and A0 and the four rare-gas crystals. Cohesive energies and lattice parameters are there shown as a fraction of the respective experimental values, thus allowing to collect in one figure per each BS all the curves and to appreciate at a glance all their features. These data represent a check for our PM method, although a detailed analysis of the performance of the two techniques is to be found in a previous paper [8]. Generally speaking, it is seen that the PM–CP difference in cohesive energies is relevant, varying between 60% and 12% along the considered lattice parameters, only for Ne crystal, whereas it is much smaller, between 14% and 5%, for the other gases. Since in general the CP technique quite probably overestimates the MP2 BSSE whereas PM underestimates it [8], thus providing underbinding and overbinding energy curves, respectively,

0.4

0.4

0.6

0.6

0.8

0.8

1

1 Ne P M Ne C P Ar P M Ar C P KrPM KrCP Xe P M Xe C P

1.2

1.4

0.9

A’

0.2

0.95

1

D1/a0

1.05

1.1

Ne P M Ne C P Ar P M Ar C P KrPM KrCP Xe P M Xe C P

1.2

1.4

0.9

0.95

1

1.05

1.1

D1/a0

Fig. 3. HF + MP2 cohesive energy as a function of the lattice parameter calculated by using basis set A (left panel) or A0 (right panel), and by correcting the MP2 contribution for BSSE either using the PM (filled symbols) or the CP (empty symbols) method. For each gas, energies and lattice parameters are reported as a fraction of the respective experimental values.

298

M. Halo et al. / Chemical Physics Letters 467 (2009) 294–298

it can be said that the error due to this wrong estimate is significant only for Ne, whereas PM data for the other rare gases are reliable. It is also seen from the figure that PM–CP difference is rather dependent on the geometry, being larger for contracted lattice parameters and viceversa. This pattern is indeed an obvious consequence of the underestimation of intra-pair correlation energy by the PM technique – or else of the overestimation of it by CP – which occurs primarily for contracted lattice parameters. Also it has to be noted that the difference is scarcely dependent on the basis set quality, for the same geometry; this could permit extrapolation procedures between the two methods and between BSs of different quality. Such an application, however, would require further assessment of the generality and the limits of the PM technique. 4. Conclusions A robust periodic local MP2 technique has been applied to the study of fcc rare-gas crystals Ne to Xe. A strong point of our approach is its simplicity. Only few computational parameters have to be optimized before running a calculation, and in fact general criteria can be used to set reasonable values for them. In spite of the known limitations of the MP2 approximation, the results are rather close to the experiment, the case of Ne being the most critical one (and for this reason analyzed in depth in a previous paper [8]). Since long-range dispersive interactions are treated exactly in the present approximation, our method appears particularly suitable for application to those molecular crystals whose structure and stability largely depends on such forces which cannot be accounted for accurately with DFT local-exchange-correlation techniques. Work in this direction is getting on. Acknowledgement We are grateful to Denis Usvyat for useful discussions and suggestions.

References [1] K. Roscis´zewski, B. Paulus, P. Fulde, H. Stoll, Phys. Rev. B 60 (1999) 7905. [2] K. Roscis´zewski, B. Paulus, P. Fulde, H. Stoll, Phys. Rev. B 62 (2000) 5482. [3] P. Schwerdtfeger, N. Gaston, R.P. Krawczyk, R. Tonner, G.E. Moyano, Phys. Rev. B 73 (2006). [4] V.F. Lotrich, K. Szalewicz, Phys. Rev. Lett. 79 (1997) 1301. [5] B. Civalleri, D.S. Middlemiss, R. Orlando, C.C. Wilson, P. Ugliengo, Chem. Phys. Lett. 451 (2008) 287. [6] J. Harl, G. Kresse, Phys. Rev. B 77 (2008) 045136. [7] S. Casassa, M. Halo, L. Maschio, J. Phys.: Conf. Ser. 117 (2008) 012007. [8] M. Halo, S. Casassa, L. Maschio, C. Pisani, Phys. Chem. Chem. Phys. (2008), doi:10.1039/B812870G. [9] C. Pisani, S. Casassa, L. Maschio, M. Halo, D. Usvyat, M. Schütz, CRYSCOR User’s Manual, b-version, University of Turin, Torino, Italy, 2008. . [10] R. Dovesi et al., CRYSTAL06 User’s Manual, Università di Torino, Torino, 2006. . [11] MOLPRO is a package of ab initio programs written by H.-J. Werner et al., version 2006.1. . [12] P. Pulay, Chem. Phys. Lett. 100 (1983) 151. [13] S. Saebø, P. Pulay, Chem. Phys. Lett. 113 (1985) 13. [14] J.P. Perdew, in: P. Ziesche, H. Eschrig (Eds.), Electronic Structure of Solids, Academie Verlag, Berlin, 1991. [15] A.D. Becke, J. Chem. Phys. 98 (1993) 5648. [16] C. Pisani et al., J. Chem. Phys. 122 (2005) 094133. [17] L. Maschio, D. Usvyat, F. Manby, S. Casassa, C. Pisani, M. Schütz, Phys. Rev. B 76 (2007) 075101. [18] D. Usvyat, L. Maschio, F. Manby, S. Casassa, M. Schütz, C. Pisani, Phys. Rev. B 76 (2007) 075102. [19] C. Pisani, L. Maschio, S. Casassa, M. Halo, M. Schütz, D. Usvyat, J. Comput. Chem. 29 (2008) 2113. [20] L. Maschio, D. Usvyat, Phys. Rev. B 78 (2008) 073102. [21] A. Nicklass, M. Dolg, H. Stoll, H. Preuss, J. Chem. Phys. 102 (1995) 8943. [22] T. Van Mourick, A.K. Wilson, T.H. Dunning Jr., Mol. Phys. 96 (1999) 529. [23] S.F. Boys, F. Bernardi, Mol. Phys. 19 (1970) 553. [24] D.N. Batchelder, D.L. Losee, R.O. Simmons, Phys. Rev. 162 (1967) 767. [25] G.L. McConville, J. Chem. Phys. 60 (1974) 4093. [26] Y. Endoh, G. Shirane, J. Skalyo Jr., Phys. Rev. B 11 (1975) 1681. [27] O.G. Peterson, D.N. Batchelder, R.O. Simmons, Phys. Rev. 150 (1976) 703. [28] L.A. Schwalbe, R.K. Crawford, H.H. Chen, R.A. Aziz, J. Chem. Phys. 66 (1977) 4493. [29] S. Gewurtz, P. Stoicheff, Phys. Rev. B 10 (1974) 3487. [30] D.L. Losee, R.O. Simmons, Phys. Rev. 172 (1968) 944. [31] J. Skalyo Jr., Y. Endoh, G. Shirane, Phys. Rev. B 9 (1974) 1797. [32] D.R. Sears, H.P. Klug, J. Chem. Phys. 37 (1967) 3002. [33] N.A. Lurie, G. Shirane, J. Skalyo Jr., Phys. Rev. B 6 (1974) 2661. [34] J.F. Olgivie, F.Y. Wang, J. Mol. Struct. 273 (1992) 277.