Communication: Non-additivity of van der Waals interactions between nanostructures Jianmin Tao and John P. Perdew Citation: The Journal of Chemical Physics 141, 141101 (2014); doi: 10.1063/1.4897957 View online: http://dx.doi.org/10.1063/1.4897957 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/141/14?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Toward transferable interatomic van der Waals interactions without electrons: The role of multipole electrostatics and many-body dispersion J. Chem. Phys. 141, 034101 (2014); 10.1063/1.4885339 Non-additivity of polarizabilities and van der Waals C6 coefficients of fullerenes J. Chem. Phys. 138, 114107 (2013); 10.1063/1.4795158 The van der Waals coefficients between carbon nanostructures and small molecules: A time-dependent density functional theory study J. Chem. Phys. 131, 164708 (2009); 10.1063/1.3256238 TimeDependent Density Functional Theory Calculation of Van der Waals Coefficient of Metal Clusters AIP Conf. Proc. 1108, 114 (2009); 10.1063/1.3117119 The van der Waals interaction between protein molecules in an electrolyte solution J. Chem. Phys. 120, 2005 (2004); 10.1063/1.1634955
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.160.4.77 On: Sat, 20 Dec 2014 07:43:25
THE JOURNAL OF CHEMICAL PHYSICS 141, 141101 (2014)
Communication: Non-additivity of van der Waals interactions between nanostructures Jianmin Tao1,a) and John P. Perdew2 1 2
Department of Chemistry, University of Pennsylvania, Philadelphia, Pennsylvania 19104-6323, USA Department of Physics, Temple University, Philadelphia, Pennsylvania 19122, USA
(Received 22 August 2014; accepted 1 October 2014; published online 10 October 2014) Due to size-dependent non-additivity, the van der Waals interaction (vdW) between nanostructures remains elusive. Here we first develop a model dynamic multipole polarizability for an inhomogeneous system that allows for a cavity. The model recovers the exact zero- and high-frequency limits and respects the paradigms of condensed matter physics (slowly varying density) and quantum chemistry (one- and two-electron densities). We find that the model can generate accurate vdW coefficients for both spherical and non-spherical clusters, with an overall mean absolute relative error of 4%, without any fitting. Based on this model, we study the non-additivity of vdW interactions. We find that there is strong non-additivity of vdW interactions between nanostructures, arising from electron delocalization, inequivalent contributions of atoms, and non-additive many-body interactions. Furthermore, we find that the non-additivity can have increasing size dependence as well as decreasing size dependence with cluster size. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4897957] In recent years, there has been an explosive interest in the study of organic nanomaterials. Part of the interest is motivated by microelectronic devices, with the capability to outperform traditional technology. However, due to the strong van der Waals (vdW) interactions between nanostructures,1 simulation of nanostructured materials presents a challenge to density functional theory (DFT). Therefore, there is a great demand for knowledge of van der Waals interactions between large molecules. The van der Waals interaction is an important long-range intermolecular interaction arising from instantaneous collective charge fluctuations. It is much weaker than the chemical bond at equilibrium, but it may become dominant in otherwise non-bonded situations. It is ubiquitous, and many phenomena arising from it (e.g., sublimation of solids, high interlayer mobility of layered materials, and physisorption) may be “observed” in our daily life. It determines higher-order structures of biomolecular chains (e.g., DNA double helix structure2 ) that define biological activities. It affects many properties of molecules and π -π stacking in nucleic bases.3 In the study of vdW interactions, the central task is to calculate the vdW coefficients, which, by second-order perturbation theory, can be expressed in terms of the dynamic multipole polarizability α l (iu).4–6 Since the polarizability is highly sensitive to basis size,7 reliable calculations require a large basis set. This significantly increases computational cost. Development of a simple and yet accurate model dynamic polarizability provides an alternative solution. Many models have been proposed8–11 to simulate α l (iu) and the vdW coefficients. In particular, Becke and Johnson8 proposed a model for vdW coefficients, based on the assumption that the vdW energy is due to the dipole interaction crea) Email:
[email protected]. ~jianmint/.
URL:
0021-9606/2014/141(14)/141101/4/$30.00
http://www.sas.upenn.edu/
ated by an electron and the associated exchange hole. This model can generate C6 with a mean absolute relative error (MARE) of 3% for free atom pairs. Based on this model, Steinmann and Corminboeuf12 have proposed a densitydependent dispersion correction, which seems promising. Tkatchenko and Scheffler10 proposed a method for intermolecular vdW coefficients, with good accuracy of less than 6%. Recently, Tkatchenko et al.11 developed an accurate α 1 (iu) to include non-additive many-body (i.e., 3-body, 4body, · · · ) interactions, improving the performance of the Tkatchenko-Scheffler method for nanosystems. It was shown11, 13, 14 that pairwise-based models can make a large error when applied to large nanostructures, due to the non-additive nature of vdW interactions. To overcome this problem, we developed a simple model15 for α l (iu), which extends the original formula16 for the dipole polarizability to the multipole one. Since our model treats the vdW interaction between two whole molecules, not the sum of pair-wise interactions between atomic fragments of the two molecules, it includes all non-additivity effects via the input static polarizability, as discussed below. It was shown15 that this model yields for free atom pairs C6 , C8 , and C10 in excellent agreement (3%) with accurate reference values. However, the solidsphere model may not be suitable to cage molecules such as fullerenes.14 The reason is that the cavity affects α l (iu) via both the electron density and the frequency dependence. The solid-sphere model can capture the first effect, but not the second one. The aim of this work is to study the non-additivity of vdW interactions between nanostructures, such as fullerenes and nanoclusters. For this purpose, we first develop a model α l (iu) for an inhomogeneous system that allows for a cavity. Based on this model, the non-additivity of vdW interactions between nanoclusters is investigated. Let us consider a shell of density, which decays rapidly outside the surface formed
141, 141101-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.160.4.77 On: Sat, 20 Dec 2014 07:43:25
141101-2
J. Tao and J. P. Perdew
by the nuclear framework. For such a system, we define an l-dependent effective shell thickness by Rl = Rn + tl /2, with Rl being the l-dependent outer radius of the shell and Rn the average radius of the nuclear framework. Then we generalize the exact dynamic polarizability of a conducting spherical shell of uniform density14, 17, 18 to a shell of inhomogeneous density. Our generalization is made by imposing the following constraints: (i) recovery of the exact static and high-frequency limits, and (ii) respect for two paradigms of condensed matter physics (slowly varying density) and quantum chemistry (one- and two-electron densities). Constraints (i) ensure the correct transfer of physics from one order to another, while constraints (ii) guarantee the consistent accuracy for different densities. A model that satisfies all four constraints is assumed to take the form 1 2l + 1 Rl 3 r 2l−2 dl4 ωl2 d r 4 2 , (1) αl (iu) = 2 4π dl Rl −tl dl ωl + u 1 − βl ρl where βl (r) = ωl2 (r)ω˜ l2 (r)/[(ωl2 (r) + u2 )(ω˜ l2 (r) + u2 )] describes the coupling of the sphere and cavity plasmon oscillations, and ρ l = (1√− tl /Rl )2l + 1 describes the shape of the shell. ωl (r) = ωp (r) l/(2l + 1) is the generalized local √ sphere plasmon frequency, ω˜ l (r) = ωp (r) (l + 1)/(2l + 1) √ is the local cavity plasmon frequency, and ωp (r) = 4π n(r). (Atomic units are used.) Rl and dl can be fixed by the zero- and high-frequency limits (see the supplementary material).15, 19 The vdW coefficients C2k are integrals over u that can be evaluated analytically (see the supplementary material).15, 19 From Eq. (1), we see that, in the tl → Rl (solid sphere) limit, ρ l = 0 and our hollow-sphere model correctly reduces to the solid-sphere model15 2l−2 4 2 d ω r 2l + 1 (2) αl (iu) = d 3 r 4 2 l l2 . 4π dl dl ωl + u In the high-frequency limit, our hollow-sphere model also reduces to the solid-sphere model of Eq. (2). However, in the static limit, the local polarizability model of Eq. (2) predicts αl (0) = [Rl2l+1 − (Rl − tl )2l+1 ]/dl . For C60 (dl ≈ 1), this yields α 1 (0) ≈ 380, which is too small, compared to the TDDFT value 534.7 The error increases with fullerene size. In the large-size limit, the local polarizability of Eq. (2) yields α l (0) = (2l + 1)tl R2l /dl , which is significantly smaller than the correct value αl (0) = Rl2l+1 /dl . Our model may be regarded as an interpolation of α l (iu) between exact zero- and high-frequency limits. As such, the accurate static multipole polarizability α l (0), which can be found from the literature (e.g., Ref. 7), is needed as input. (For large systems, α l (0) can be estimated from α 1 (0); see discussion below.) Since the model respects the two paradigms, it is expected to be accurate for systems of any size. This expected size-independent performance can be understood from the energy-gap viewpoint. A remarkable feature of the present model is that it is valid for any value of tl , including the tl /Rl → 1 (solid-sphere) and tl /Rl → 0 (large-size) limits. The general validity and nearly the same accuracy for all orders allow us to study on the same footing the vdW interactions for systems that may or may not have a cavity. The model may
J. Chem. Phys. 141, 141101 (2014)
be simplified with the single-frequency approximation (SFA), which assumes that (i) only valence electrons are polarizable, and (ii) the density is uniform within the shell and zero otherwise. This simplified SFA takes the same expression as the full hollow-sphere model [Eq. (1)], but with n(r) replaced by the average valence electron density n¯ = N/Vl . Here N is the number of valence electrons, and Vl is the l-dependent shell volume. Note that our SFA is different from the classical shell model (CSM).18 In SFA, the cutoff radius is determined by Rl = [α l (0)]1/(2l+1) , rather than by the sharp single boundary Rl = R1 . While less accurate than the full hollow-sphere expression (1), it requires no detailed knowledge of the density and no numerical integration (see the supplementary material).19 It is worth pointing out that our model is not limited to spherical systems. The non-spherical effect can enter the model via the exact static and high-frequency limits. For nonspherical systems, we use the isotropic α¯ l (0). In addition, since the energy gap is related to system size and since the energy gap can be also accounted for via α l (0), the model should be valid for any system size and the accuracy should be size-independent. The general validity of the model is the most appealing feature and has been demonstrated by applying it to non-spherical H-terminated silicon cluster HM SiN . These and other results are reported in detail in the supplementary material.19 Now we apply SFA to NaN -NaN (no cavity) and C60 NaN cluster pairs. We find that the mean absolute relative deviations of C6 from the TD Hartree-Fock (TDHF) values20 are only 6% and 7%, respectively, for N = 2−20. To make a comparison, we repeat calculations for C60 -NaN pairs with the solid-sphere model within SFA. We find that MARE increases to 22%, which is significantly larger than that for the hollow-sphere model-SFA by a factor of 3. For a C60 dimer, the solid-sphere model predicts C6 = 86.6 × 103 and C8 = 365.5 × 105 , which are 32% and 27% smaller than the TDDFT values (C6 = 126.5 × 103 , C8 = 498.7 × 105 ).21 Next we compare the full hollow-sphere expression with SFA for C60 -Li pair. Our calculation shows that the full hollowsphere model (C6 = 8.31 × 103 ) is indeed closer to the experimental value 8.05 × 103 (Ref. 22) than SFA (9.14 × 103 ). Finally, we apply SFA to H-terminated Si clusters from SiH4 to Si172 H120 . Calculation shows that our C6 are in excellent agreement (2%) with more expensive TDDFT values.23 The overall performance of the model is nearly size-independent for the clusters considered here (see Tables S5, S6, and S8 of the supplementary material19 ), with a mean absolute relative deviation of only 4%. (TDDFT has an error of 5%,23 while the error of TDHF is smaller. So they are reliable and can be used as standard for comparison.) As an important application, we study the vdW coefficients for fullerene (CN ) pairs. For small fullerenes, α 1 (0) can be calculated accurately from TDDFT,7 while α 2 (0) and α 3 (0) may be estimated from CSM. For a C60 dimer, SFA yields C6 = 136.7 × 103 and C8 = 545.2 × 105 , which agree well (about 8%) with the TDDFT values. However, there is a large discrepancy of C10 between the SFA (1495 × 107 ) and the literature value (446.0 × 107 ),21 which arises from the error in estimating α 3 (0). This is similar to sodium clusters.24 For instance, for Na20 -Na20 , SFA predicts
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.160.4.77 On: Sat, 20 Dec 2014 07:43:25
J. Tao and J. P. Perdew
J. Chem. Phys. 141, 141101 (2014)
Na cluster pair
0.4 2
10 100 N (cluster size)
0.4 (c) Si cluster pair
3
0.3 0.2
2 5 10 15 20 25 30 N
Fullerene pair
20
100
1000 5000 N
0.4 0.3
(d) Ge cluster pair
C6
(a)
1.2
C8
1.1
C10
1 0
0.2 0
C2k(t)/C2k(t=2.8)
(b)
0.4
1.3
2
4 6 t (bohr)
8
1.3 1.2
C6
(c)
C8 C10
1.1 1 2
4
6
t (bohr)
8
8 C10
6
(b)
C8 C6
4 2 0 0
C2k(M)/C2k(M=60)
0.6
0.6
ing trend with M. The failure to follow the normal N × M size dependence is more serious for NaM -NaN , CM -NaN , and CN -CN than for GeN -GeN . In the bulk (N → ∞) limit, in which the energy gap in CN (graphene) and NaN vanishes, the (valence) electrons are nonlocal and the non-additivity of vdW interactions is so strong that the pairwise picture completely fails.14, 34 Even with finite energy gap, it may produce a large error, as observed in H-terminated Si clusters11 and other systems.35–37 From this study, we see that non-additivity of vdW interactions may arise from electron delocalization (or energy gap) and/or from nonadditive many-body interactions,11, 13 and/or from inequivalent contributions of each atom pair,14 in which atoms further from the center contribute more, and this inequivalence increases from the leading-order to higher-order coefficients. While these effects have strong size-dependence and are included in our model via the input α l (0), they are missing in the standard atom pairwise picture. Consequently, the error of this picture grows with system size and may affect the performance of other pairwise-based models for nanostructures.11 In order to get a fuller understanding of the shape and size effects, we imagine modifying C60 in each of four ways that change the vdW coefficients (or the vdW coefficients per atom pair) between two modified fullerenes. (a) First, we fix the number of valence electrons N and the outer radius Rl of the shell, while varying tl from 0 to Rl . This corresponds to compressing or expanding the material in the shell without changing the outer radius. The result is plotted in Fig. 2(a). It is observed from Fig. 2(a) that, relative to the tl = 0 case, C6 (or C6 per atom pair) monotonically increases with tl and reaches a maximum as the shell evolves into a solid sphere, suggesting that the non-additivity increases with shell thickness. However, C8 and C10 display a local maximum and then drop to a minimum, which is slightly larger than the tl = 0 value. (b) Then we hold N and the valence electron density fixed, while varying tl and Rl . When tl is equal to Rl , the shell collapses into a solid sphere. This corresponds to changing C2k(t)/C2k(t=0)
(a)
0.8
C6 (X103)/atom pair
1.1 1
C6 (X10 )/atom pair
C6 (X103)/atom pair
C6 (X103)/atom pair
C8 = 239 × 106 and C10 = 157 × 109 , while CSM (with t = R) yields larger C8 = 246 × 106 and C10 = 169 × 109 . For large fullerenes (N ≥ 110), evaluation of Rn is more reliable25, 26 than actual calculation of α 1 (0), due to the basis size issue. Here instead of taking α 1 (0) from the literature,26, 27 we estimate it from α 1 (0) = (Rn + t1 /2)3 , in which t1 is obtained from t1 = 3.24 + 0.16(N − 110)/(N + 1), an interpolation between the thickness of C110 (3.24) and the large-N limit (3.4).27 We apply this formula to calculate α 1 (0) for C60 (Ih ), C70 (D5h ), C78 (C2v ), and C84 (D2 ), using the Rn estimated from Ref. 25. The results are 525, 638, 733, and 811, respectively, which are in good agreement (about 2%) with TDHF28 (537, 659, 748, 806). This simple method provides the best estimates of α 1 (0) for large fullerenes. They are tabulated as input in Table S1 of the supplementary material.19 An error in the estimate of α 1 (0) may carry over twice that error to C6 and slightly more to C8 . The atom pairwise picture29, 30 is perhaps the most popular method for the evaluation of C2k between large systems. In this picture, the C2k are obtained as a sum of atom pair interactions, where each pair involves one atom from each of the two systems. Because of simplicity, it has been widely used8, 29, 31, 32 in chemistry and physics. However, examination of its validity has been focused on small molecules. Here we will examine it for a variety of nanocluster pairs. For a cluster pair with M atoms in cluster A and N in cluster B, the vdW coefficients per atom pair is33 C2k /MN . If all the atoms of a cluster make equal contributions (as for C6 in the simplest version of the atom pairwise picture), C2k per atom pair should be nearly size-independent.14 To our surprise, this picture only works fairly well with C6 for SiN -SiN (no cavity), but not for other cluster pairs, as shown in Fig. 1. From Figs. 1(a) and 1(b) we see that C6 per atom pair for CN -CN exhibits very strong non-additivity, leading to a rapidly increasing trend with N ,13, 14, 28 while the trend is opposite for NaN -NaN . Figure 1(d) shows an oscillating increasing trend for GeN GeN , suggesting a very possible failure of this picture, as N increases. Our calculation also shows that, for C60 -NaN , C6 per atom pair displays a decreasing trend with N , like NaN -NaN , while for a fixed N , CM -NaN shows an increas-
C2k(t)/C2k(t=2.8)
141101-3
2
1 0.8 0.6 0.4 0.2 0
4 6 t (bohr)
8
(d) C6, C8, C10 0
20
40
60
M (Number of atoms)
0.2 2 5 10 15 20 25 30 N
FIG. 1. C6 per atom pair in SFA and CSM vs. cluster size for (a) NaN -NaN , (b) CN -CN , (c) SiN -SiN , and (d) GeN -GeN pairs. C6 /N2 tends to increase or decrease with N according as α 1 (0)/N does.
FIG. 2. Effects on the vdW coefficients in CSM due to identical modifications of each half of a C60 -C60 pair. Panel (a): N and R are held fixed. At t = 0, C6 = 0.033 × 103 , C8 = 0.123 × 105 , C10 = 0.382 × 107 . Panel (b): N and the average valence electron density n are held fixed. Panel (c): R and n are held fixed. Panel (d): R and t are held fixed. The parameters fixed at their C60 values are N = 240, t = 2.8, R = 8.1, n = 0.150.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.160.4.77 On: Sat, 20 Dec 2014 07:43:25
141101-4
J. Tao and J. P. Perdew
the way in which the uncompressed and unexpanded material is shaped. The vdW coefficients reach a minimum at tl = Rl , because the radius of the shell is smallest in this case and so is the static polarizability. When tl tends to zero, Rl 1/2 → ∞ as (1/tl ). This makes the polarizability grow rapidly and so do the vdW coefficients, as shown in Fig. 2(b). This leads to the decrease of non-additivity. (c) Next, we hold Rl and the valence electron density at our C60 value, but vary tl from its C60 value to Rl . This shows the effect on C6 , C8 , and C10 of gradually filling in the cavity with about 23 additional atoms. From Fig. 2(c) we see that all the coefficients increase with tl , but with the strongest effect on C6 , followed by C8 and C10 . In particular we observe that atoms further from the center contribute more to the vdW coefficients, and especially to the higher-order ones. This results in the increase of nonadditivity. (d) Finally we fix Rl and tl at our C60 values, but vary the atom number M from 0 to 60. This shows the effect of filling in vacancies in the single-walled fullerene. We observe from Fig. 2(d) that the ratio of the vdW coefficients to their M = 60 values grows with M rapidly, leading to the increase of non-additivity as well. In conclusion, we have developed a model α l (iu) for inhomogeneous systems that allow for a cavity. The vdW coefficients generated from this model are consistently accurate for both small and large molecules (i.e., no size dependence). The needed inputs are the accurate static multipole polarizability α l (0), Rn (for systems with a cavity) and the electron density (the average valence electron density in SFA). Although computing α l (0) needs some effort, it is much cheaper, compared to the full frequency calculation of TDDFT and other ab initio methods. α l (0) only involves a single-point calculation, which may be performed even with ground-state DFT,38 while for α l (iu), we have to repeat calculations at many frequency points.39 Based on the model, we study the nonadditivity of vdW interactions. We find that there is strong size-dependent non-additivity between nanostructures. The knowledge of non-addivity and the long-range nature of vdW is important for our understanding of ubiquitous intermolecular forces and also for developing better DFT methods.40 J.T. thanks Professor Aparna Chakrabarti for valuable communications. J.P.P. thanks Professor K. A. Jackson for helpful discussion and for sharing the work of Ref. 38 prior publication. This work was supported by two NSF grants under Grant Nos. DMR-1305135 (JPP) and CHE 1261918 (JT), and by the Air Force Office of Scientific Research, Air Force Materiel Command, USAF, under Grant No. FA9550-10-10248 (via Professor Andrew M. Rappe). Computational support was provided by the HPCMO and the NERSC. 1 R.
Loutfy and E. Wexler, Perspectives of Fullerene Nanotechnology edited ¯ by E. Osawa (Kluwer Academic Publishers, New York, 2002), Part VII, pp. 275–280.
J. Chem. Phys. 141, 141101 (2014) 2 R.
A. DiStasio, O. A. von Lilienfeld, and A. Tkatchenko, Proc. Natl. Acad. Sci. U.S.A. 109, 14791 (2012). 3 C. A. Hunter and J. K. M. Sanders, J. Am. Chem. Soc. 112, 5525 (1990). 4 S. H. Patil and K. T. Tang, J. Chem. Phys. 106, 2298 (1997). 5 J. Mitroy and M. W. J. Bromley, Phys. Rev. A 68, 052714 (2003) (Erratum); 71, 019902 (2005) (Erratum). 6 S. J. A. van Gisbergen, J. G. Snijders, and E. J. Baerends, J. Chem. Phys. 103, 9347 (1995). 7 M. van Faassen, L. Jensen, J. A. Berger, and P. L. de Boeij, Chem. Phys. Lett. 395, 274 (2004). 8 A. D. Becke and E. R. Johnson, J. Chem. Phys. 127, 154108 (2007). 9 O. A. Vydrov and T. Van Voorhis, Phys. Rev. A 81, 062708 (2010). 10 A. Tkatchenko and M. Scheffler, Phys. Rev. Lett. 102, 073005 (2009). 11 A. Tkatchenko, R. A. DiStasio, R. Car, and M. Scheffler, Phys. Rev. Lett. 108, 236402 (2012). 12 S. N. Steinmann and C. Corminboeuf, J. Chem. Theory Comput. 7, 3567 (2011). 13 V. V. Gobre and A. Tkatchenko, Nat. Commun. 4, 2341 (2013). 14 A. Ruzsinszky, J. P. Perdew, J. Tao, G. I. Csonka, and J. M. Pitarke, Phys. Rev. Lett. 109, 233203 (2012). 15 J. Tao, J. P. Perdew, and A. Ruzsinszky, Proc. Natl. Acad. Sci. U.S.A. 109, 18 (2012). 16 J. Tao, J. P. Perdew, and A. Ruzsinszky, Phys. Rev. B 81, 233102 (2010). 17 A. A. Lucas, L. Henrard, and Ph. Lambin, Phys. Rev. B 49, 2888 (1994). 18 J. P. Perdew, J. Tao, P. Hao, A. Ruzsinszky, G. I. Csonka, and J. M. Pitarke, J. Phys.: Condens. Matter 24, 424207 (2012). 19 See supplementary material at http://dx.doi.org/10.1063/1.4897957 for detailed derivation. 20 A. Jiemchooroj, P. Norman, and B. E. Sernelius, J. Chem. Phys. 125, 124306 (2006). 21 J. M. Pacheco and J. P. P. Ramalho, Phys. Rev. Lett. 79, 3873 (1997). 22 A. Banerjee, J. Autschbach, and A. Chakrabaru, Phys. Rev. A 78, 032704 (2008). 23 S. Botti, A. Castro, X. Andrade, A. Rubio, and M. A. L. Marques, Phys. Rev. B 78, 035333 (2008). 24 L. Serra, F. Garcias, M. Barranco, J. Navarro, L. C. Balbas, A. Rubio, and A. Mananes, J. Phys.: Condens. Matter 1, 10391 (1989). 25 G. B. Adams, M. O’Keeffe, and R. S. Ruoff, J. Phys. Chem. 98, 9465 (1994). 26 R. R. Zope, T. Baruah, M. R. Pederson, and B. I. Dunlap, Phys. Rev. B 77, 115452 (2008). 27 G. K. Gueorguiev, J. M. Pacheco, and D. Tománek, Phys. Rev. Lett. 92, 215501 (2004). 28 J. Kauczor, P. Norman, and W. A. Saidi, J. Chem. Phys. 138, 114107 (2013). 29 V. A. Parsegian, Van der Waals Forces (Cambridge University Press, Cambridge, 2005). 30 H. C. Hamaker, Physica 4, 1058 (1937). 31 J. F. Dobson and T. Gould, J. Phys.: Condens. Matter 24, 073201 (2012). 32 S. Grimme, Comput. Mol. Sci. 1, 211 (2011). 33 There is a small additional contribution for C and C given by 21C MA 8 10 6 i=1 M M M j =1B δij2 /(MA MB ) and i=1A j =1B [126C6 δij4 + 36C8 δij2 ]/(MA MB ), re-
spectively. Here δij = |R − rij | and rij is the distance of atom pairs. F. Dobson, A. White, and A. Rubio, Phys. Rev. Lett. 96, 073201 (2006). 35 R.-F. Liu, J. G. Ángyán, and J. F. Dobson, J. Chem. Phys. 134, 114106 (2011). 36 J. F. Dobson, Int. J. Quantum Chem. 114, 1157 (2014). 37 A. J. Misquitta, J. Spencer, A. J. Stone, and A. Alavi, Phys. Rev. B 82, 075312 (2010). 38 L. Ma, K. A. Jackson, J. Wang, M. Horoi, and J. Jellinek, Phys. Rev. B 89, 035429 (2014). 39 A. Derevianko, S. G. Porsev, and J. F. Babb, At. Data. Nucl. Data Tables 96, 323 (2010). 40 R. Peverati and D. G. Truhlar, Philos. Trans. R. Soc. A 372, 20120476 (2014). 34 J.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.160.4.77 On: Sat, 20 Dec 2014 07:43:25
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.160.4.77 On: Sat, 20 Dec 2014 07:43:25
THE JOURNAL OF CHEMICAL PHYSICS 141, 141101 (2014)
Communication: Non-additivity of van der Waals interactions between nanostructures Jianmin Tao1,a) and John P. Perdew2 1 2
Department of Chemistry, University of Pennsylvania, Philadelphia, Pennsylvania 19104-6323, USA Department of Physics, Temple University, Philadelphia, Pennsylvania 19122, USA
(Received 22 August 2014; accepted 1 October 2014; published online 10 October 2014) Due to size-dependent non-additivity, the van der Waals interaction (vdW) between nanostructures remains elusive. Here we first develop a model dynamic multipole polarizability for an inhomogeneous system that allows for a cavity. The model recovers the exact zero- and high-frequency limits and respects the paradigms of condensed matter physics (slowly varying density) and quantum chemistry (one- and two-electron densities). We find that the model can generate accurate vdW coefficients for both spherical and non-spherical clusters, with an overall mean absolute relative error of 4%, without any fitting. Based on this model, we study the non-additivity of vdW interactions. We find that there is strong non-additivity of vdW interactions between nanostructures, arising from electron delocalization, inequivalent contributions of atoms, and non-additive many-body interactions. Furthermore, we find that the non-additivity can have increasing size dependence as well as decreasing size dependence with cluster size. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4897957] In recent years, there has been an explosive interest in the study of organic nanomaterials. Part of the interest is motivated by microelectronic devices, with the capability to outperform traditional technology. However, due to the strong van der Waals (vdW) interactions between nanostructures,1 simulation of nanostructured materials presents a challenge to density functional theory (DFT). Therefore, there is a great demand for knowledge of van der Waals interactions between large molecules. The van der Waals interaction is an important long-range intermolecular interaction arising from instantaneous collective charge fluctuations. It is much weaker than the chemical bond at equilibrium, but it may become dominant in otherwise non-bonded situations. It is ubiquitous, and many phenomena arising from it (e.g., sublimation of solids, high interlayer mobility of layered materials, and physisorption) may be “observed” in our daily life. It determines higher-order structures of biomolecular chains (e.g., DNA double helix structure2 ) that define biological activities. It affects many properties of molecules and π -π stacking in nucleic bases.3 In the study of vdW interactions, the central task is to calculate the vdW coefficients, which, by second-order perturbation theory, can be expressed in terms of the dynamic multipole polarizability α l (iu).4–6 Since the polarizability is highly sensitive to basis size,7 reliable calculations require a large basis set. This significantly increases computational cost. Development of a simple and yet accurate model dynamic polarizability provides an alternative solution. Many models have been proposed8–11 to simulate α l (iu) and the vdW coefficients. In particular, Becke and Johnson8 proposed a model for vdW coefficients, based on the assumption that the vdW energy is due to the dipole interaction crea) Email:
[email protected]. ~jianmint/.
URL:
0021-9606/2014/141(14)/141101/4/$30.00
http://www.sas.upenn.edu/
ated by an electron and the associated exchange hole. This model can generate C6 with a mean absolute relative error (MARE) of 3% for free atom pairs. Based on this model, Steinmann and Corminboeuf12 have proposed a densitydependent dispersion correction, which seems promising. Tkatchenko and Scheffler10 proposed a method for intermolecular vdW coefficients, with good accuracy of less than 6%. Recently, Tkatchenko et al.11 developed an accurate α 1 (iu) to include non-additive many-body (i.e., 3-body, 4body, · · · ) interactions, improving the performance of the Tkatchenko-Scheffler method for nanosystems. It was shown11, 13, 14 that pairwise-based models can make a large error when applied to large nanostructures, due to the non-additive nature of vdW interactions. To overcome this problem, we developed a simple model15 for α l (iu), which extends the original formula16 for the dipole polarizability to the multipole one. Since our model treats the vdW interaction between two whole molecules, not the sum of pair-wise interactions between atomic fragments of the two molecules, it includes all non-additivity effects via the input static polarizability, as discussed below. It was shown15 that this model yields for free atom pairs C6 , C8 , and C10 in excellent agreement (3%) with accurate reference values. However, the solidsphere model may not be suitable to cage molecules such as fullerenes.14 The reason is that the cavity affects α l (iu) via both the electron density and the frequency dependence. The solid-sphere model can capture the first effect, but not the second one. The aim of this work is to study the non-additivity of vdW interactions between nanostructures, such as fullerenes and nanoclusters. For this purpose, we first develop a model α l (iu) for an inhomogeneous system that allows for a cavity. Based on this model, the non-additivity of vdW interactions between nanoclusters is investigated. Let us consider a shell of density, which decays rapidly outside the surface formed
141, 141101-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.160.4.77 On: Sat, 20 Dec 2014 07:43:25
141101-2
J. Tao and J. P. Perdew
by the nuclear framework. For such a system, we define an l-dependent effective shell thickness by Rl = Rn + tl /2, with Rl being the l-dependent outer radius of the shell and Rn the average radius of the nuclear framework. Then we generalize the exact dynamic polarizability of a conducting spherical shell of uniform density14, 17, 18 to a shell of inhomogeneous density. Our generalization is made by imposing the following constraints: (i) recovery of the exact static and high-frequency limits, and (ii) respect for two paradigms of condensed matter physics (slowly varying density) and quantum chemistry (one- and two-electron densities). Constraints (i) ensure the correct transfer of physics from one order to another, while constraints (ii) guarantee the consistent accuracy for different densities. A model that satisfies all four constraints is assumed to take the form 1 2l + 1 Rl 3 r 2l−2 dl4 ωl2 d r 4 2 , (1) αl (iu) = 2 4π dl Rl −tl dl ωl + u 1 − βl ρl where βl (r) = ωl2 (r)ω˜ l2 (r)/[(ωl2 (r) + u2 )(ω˜ l2 (r) + u2 )] describes the coupling of the sphere and cavity plasmon oscillations, and ρ l = (1√− tl /Rl )2l + 1 describes the shape of the shell. ωl (r) = ωp (r) l/(2l + 1) is the generalized local √ sphere plasmon frequency, ω˜ l (r) = ωp (r) (l + 1)/(2l + 1) √ is the local cavity plasmon frequency, and ωp (r) = 4π n(r). (Atomic units are used.) Rl and dl can be fixed by the zero- and high-frequency limits (see the supplementary material).15, 19 The vdW coefficients C2k are integrals over u that can be evaluated analytically (see the supplementary material).15, 19 From Eq. (1), we see that, in the tl → Rl (solid sphere) limit, ρ l = 0 and our hollow-sphere model correctly reduces to the solid-sphere model15 2l−2 4 2 d ω r 2l + 1 (2) αl (iu) = d 3 r 4 2 l l2 . 4π dl dl ωl + u In the high-frequency limit, our hollow-sphere model also reduces to the solid-sphere model of Eq. (2). However, in the static limit, the local polarizability model of Eq. (2) predicts αl (0) = [Rl2l+1 − (Rl − tl )2l+1 ]/dl . For C60 (dl ≈ 1), this yields α 1 (0) ≈ 380, which is too small, compared to the TDDFT value 534.7 The error increases with fullerene size. In the large-size limit, the local polarizability of Eq. (2) yields α l (0) = (2l + 1)tl R2l /dl , which is significantly smaller than the correct value αl (0) = Rl2l+1 /dl . Our model may be regarded as an interpolation of α l (iu) between exact zero- and high-frequency limits. As such, the accurate static multipole polarizability α l (0), which can be found from the literature (e.g., Ref. 7), is needed as input. (For large systems, α l (0) can be estimated from α 1 (0); see discussion below.) Since the model respects the two paradigms, it is expected to be accurate for systems of any size. This expected size-independent performance can be understood from the energy-gap viewpoint. A remarkable feature of the present model is that it is valid for any value of tl , including the tl /Rl → 1 (solid-sphere) and tl /Rl → 0 (large-size) limits. The general validity and nearly the same accuracy for all orders allow us to study on the same footing the vdW interactions for systems that may or may not have a cavity. The model may
J. Chem. Phys. 141, 141101 (2014)
be simplified with the single-frequency approximation (SFA), which assumes that (i) only valence electrons are polarizable, and (ii) the density is uniform within the shell and zero otherwise. This simplified SFA takes the same expression as the full hollow-sphere model [Eq. (1)], but with n(r) replaced by the average valence electron density n¯ = N/Vl . Here N is the number of valence electrons, and Vl is the l-dependent shell volume. Note that our SFA is different from the classical shell model (CSM).18 In SFA, the cutoff radius is determined by Rl = [α l (0)]1/(2l+1) , rather than by the sharp single boundary Rl = R1 . While less accurate than the full hollow-sphere expression (1), it requires no detailed knowledge of the density and no numerical integration (see the supplementary material).19 It is worth pointing out that our model is not limited to spherical systems. The non-spherical effect can enter the model via the exact static and high-frequency limits. For nonspherical systems, we use the isotropic α¯ l (0). In addition, since the energy gap is related to system size and since the energy gap can be also accounted for via α l (0), the model should be valid for any system size and the accuracy should be size-independent. The general validity of the model is the most appealing feature and has been demonstrated by applying it to non-spherical H-terminated silicon cluster HM SiN . These and other results are reported in detail in the supplementary material.19 Now we apply SFA to NaN -NaN (no cavity) and C60 NaN cluster pairs. We find that the mean absolute relative deviations of C6 from the TD Hartree-Fock (TDHF) values20 are only 6% and 7%, respectively, for N = 2−20. To make a comparison, we repeat calculations for C60 -NaN pairs with the solid-sphere model within SFA. We find that MARE increases to 22%, which is significantly larger than that for the hollow-sphere model-SFA by a factor of 3. For a C60 dimer, the solid-sphere model predicts C6 = 86.6 × 103 and C8 = 365.5 × 105 , which are 32% and 27% smaller than the TDDFT values (C6 = 126.5 × 103 , C8 = 498.7 × 105 ).21 Next we compare the full hollow-sphere expression with SFA for C60 -Li pair. Our calculation shows that the full hollowsphere model (C6 = 8.31 × 103 ) is indeed closer to the experimental value 8.05 × 103 (Ref. 22) than SFA (9.14 × 103 ). Finally, we apply SFA to H-terminated Si clusters from SiH4 to Si172 H120 . Calculation shows that our C6 are in excellent agreement (2%) with more expensive TDDFT values.23 The overall performance of the model is nearly size-independent for the clusters considered here (see Tables S5, S6, and S8 of the supplementary material19 ), with a mean absolute relative deviation of only 4%. (TDDFT has an error of 5%,23 while the error of TDHF is smaller. So they are reliable and can be used as standard for comparison.) As an important application, we study the vdW coefficients for fullerene (CN ) pairs. For small fullerenes, α 1 (0) can be calculated accurately from TDDFT,7 while α 2 (0) and α 3 (0) may be estimated from CSM. For a C60 dimer, SFA yields C6 = 136.7 × 103 and C8 = 545.2 × 105 , which agree well (about 8%) with the TDDFT values. However, there is a large discrepancy of C10 between the SFA (1495 × 107 ) and the literature value (446.0 × 107 ),21 which arises from the error in estimating α 3 (0). This is similar to sodium clusters.24 For instance, for Na20 -Na20 , SFA predicts
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.160.4.77 On: Sat, 20 Dec 2014 07:43:25
J. Tao and J. P. Perdew
J. Chem. Phys. 141, 141101 (2014)
Na cluster pair
0.4 2
10 100 N (cluster size)
0.4 (c) Si cluster pair
3
0.3 0.2
2 5 10 15 20 25 30 N
Fullerene pair
20
100
1000 5000 N
0.4 0.3
(d) Ge cluster pair
C6
(a)
1.2
C8
1.1
C10
1 0
0.2 0
C2k(t)/C2k(t=2.8)
(b)
0.4
1.3
2
4 6 t (bohr)
8
1.3 1.2
C6
(c)
C8 C10
1.1 1 2
4
6
t (bohr)
8
8 C10
6
(b)
C8 C6
4 2 0 0
C2k(M)/C2k(M=60)
0.6
0.6
ing trend with M. The failure to follow the normal N × M size dependence is more serious for NaM -NaN , CM -NaN , and CN -CN than for GeN -GeN . In the bulk (N → ∞) limit, in which the energy gap in CN (graphene) and NaN vanishes, the (valence) electrons are nonlocal and the non-additivity of vdW interactions is so strong that the pairwise picture completely fails.14, 34 Even with finite energy gap, it may produce a large error, as observed in H-terminated Si clusters11 and other systems.35–37 From this study, we see that non-additivity of vdW interactions may arise from electron delocalization (or energy gap) and/or from nonadditive many-body interactions,11, 13 and/or from inequivalent contributions of each atom pair,14 in which atoms further from the center contribute more, and this inequivalence increases from the leading-order to higher-order coefficients. While these effects have strong size-dependence and are included in our model via the input α l (0), they are missing in the standard atom pairwise picture. Consequently, the error of this picture grows with system size and may affect the performance of other pairwise-based models for nanostructures.11 In order to get a fuller understanding of the shape and size effects, we imagine modifying C60 in each of four ways that change the vdW coefficients (or the vdW coefficients per atom pair) between two modified fullerenes. (a) First, we fix the number of valence electrons N and the outer radius Rl of the shell, while varying tl from 0 to Rl . This corresponds to compressing or expanding the material in the shell without changing the outer radius. The result is plotted in Fig. 2(a). It is observed from Fig. 2(a) that, relative to the tl = 0 case, C6 (or C6 per atom pair) monotonically increases with tl and reaches a maximum as the shell evolves into a solid sphere, suggesting that the non-additivity increases with shell thickness. However, C8 and C10 display a local maximum and then drop to a minimum, which is slightly larger than the tl = 0 value. (b) Then we hold N and the valence electron density fixed, while varying tl and Rl . When tl is equal to Rl , the shell collapses into a solid sphere. This corresponds to changing C2k(t)/C2k(t=0)
(a)
0.8
C6 (X103)/atom pair
1.1 1
C6 (X10 )/atom pair
C6 (X103)/atom pair
C6 (X103)/atom pair
C8 = 239 × 106 and C10 = 157 × 109 , while CSM (with t = R) yields larger C8 = 246 × 106 and C10 = 169 × 109 . For large fullerenes (N ≥ 110), evaluation of Rn is more reliable25, 26 than actual calculation of α 1 (0), due to the basis size issue. Here instead of taking α 1 (0) from the literature,26, 27 we estimate it from α 1 (0) = (Rn + t1 /2)3 , in which t1 is obtained from t1 = 3.24 + 0.16(N − 110)/(N + 1), an interpolation between the thickness of C110 (3.24) and the large-N limit (3.4).27 We apply this formula to calculate α 1 (0) for C60 (Ih ), C70 (D5h ), C78 (C2v ), and C84 (D2 ), using the Rn estimated from Ref. 25. The results are 525, 638, 733, and 811, respectively, which are in good agreement (about 2%) with TDHF28 (537, 659, 748, 806). This simple method provides the best estimates of α 1 (0) for large fullerenes. They are tabulated as input in Table S1 of the supplementary material.19 An error in the estimate of α 1 (0) may carry over twice that error to C6 and slightly more to C8 . The atom pairwise picture29, 30 is perhaps the most popular method for the evaluation of C2k between large systems. In this picture, the C2k are obtained as a sum of atom pair interactions, where each pair involves one atom from each of the two systems. Because of simplicity, it has been widely used8, 29, 31, 32 in chemistry and physics. However, examination of its validity has been focused on small molecules. Here we will examine it for a variety of nanocluster pairs. For a cluster pair with M atoms in cluster A and N in cluster B, the vdW coefficients per atom pair is33 C2k /MN . If all the atoms of a cluster make equal contributions (as for C6 in the simplest version of the atom pairwise picture), C2k per atom pair should be nearly size-independent.14 To our surprise, this picture only works fairly well with C6 for SiN -SiN (no cavity), but not for other cluster pairs, as shown in Fig. 1. From Figs. 1(a) and 1(b) we see that C6 per atom pair for CN -CN exhibits very strong non-additivity, leading to a rapidly increasing trend with N ,13, 14, 28 while the trend is opposite for NaN -NaN . Figure 1(d) shows an oscillating increasing trend for GeN GeN , suggesting a very possible failure of this picture, as N increases. Our calculation also shows that, for C60 -NaN , C6 per atom pair displays a decreasing trend with N , like NaN -NaN , while for a fixed N , CM -NaN shows an increas-
C2k(t)/C2k(t=2.8)
141101-3
2
1 0.8 0.6 0.4 0.2 0
4 6 t (bohr)
8
(d) C6, C8, C10 0
20
40
60
M (Number of atoms)
0.2 2 5 10 15 20 25 30 N
FIG. 1. C6 per atom pair in SFA and CSM vs. cluster size for (a) NaN -NaN , (b) CN -CN , (c) SiN -SiN , and (d) GeN -GeN pairs. C6 /N2 tends to increase or decrease with N according as α 1 (0)/N does.
FIG. 2. Effects on the vdW coefficients in CSM due to identical modifications of each half of a C60 -C60 pair. Panel (a): N and R are held fixed. At t = 0, C6 = 0.033 × 103 , C8 = 0.123 × 105 , C10 = 0.382 × 107 . Panel (b): N and the average valence electron density n are held fixed. Panel (c): R and n are held fixed. Panel (d): R and t are held fixed. The parameters fixed at their C60 values are N = 240, t = 2.8, R = 8.1, n = 0.150.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.160.4.77 On: Sat, 20 Dec 2014 07:43:25
141101-4
J. Tao and J. P. Perdew
the way in which the uncompressed and unexpanded material is shaped. The vdW coefficients reach a minimum at tl = Rl , because the radius of the shell is smallest in this case and so is the static polarizability. When tl tends to zero, Rl 1/2 → ∞ as (1/tl ). This makes the polarizability grow rapidly and so do the vdW coefficients, as shown in Fig. 2(b). This leads to the decrease of non-additivity. (c) Next, we hold Rl and the valence electron density at our C60 value, but vary tl from its C60 value to Rl . This shows the effect on C6 , C8 , and C10 of gradually filling in the cavity with about 23 additional atoms. From Fig. 2(c) we see that all the coefficients increase with tl , but with the strongest effect on C6 , followed by C8 and C10 . In particular we observe that atoms further from the center contribute more to the vdW coefficients, and especially to the higher-order ones. This results in the increase of nonadditivity. (d) Finally we fix Rl and tl at our C60 values, but vary the atom number M from 0 to 60. This shows the effect of filling in vacancies in the single-walled fullerene. We observe from Fig. 2(d) that the ratio of the vdW coefficients to their M = 60 values grows with M rapidly, leading to the increase of non-additivity as well. In conclusion, we have developed a model α l (iu) for inhomogeneous systems that allow for a cavity. The vdW coefficients generated from this model are consistently accurate for both small and large molecules (i.e., no size dependence). The needed inputs are the accurate static multipole polarizability α l (0), Rn (for systems with a cavity) and the electron density (the average valence electron density in SFA). Although computing α l (0) needs some effort, it is much cheaper, compared to the full frequency calculation of TDDFT and other ab initio methods. α l (0) only involves a single-point calculation, which may be performed even with ground-state DFT,38 while for α l (iu), we have to repeat calculations at many frequency points.39 Based on the model, we study the nonadditivity of vdW interactions. We find that there is strong size-dependent non-additivity between nanostructures. The knowledge of non-addivity and the long-range nature of vdW is important for our understanding of ubiquitous intermolecular forces and also for developing better DFT methods.40 J.T. thanks Professor Aparna Chakrabarti for valuable communications. J.P.P. thanks Professor K. A. Jackson for helpful discussion and for sharing the work of Ref. 38 prior publication. This work was supported by two NSF grants under Grant Nos. DMR-1305135 (JPP) and CHE 1261918 (JT), and by the Air Force Office of Scientific Research, Air Force Materiel Command, USAF, under Grant No. FA9550-10-10248 (via Professor Andrew M. Rappe). Computational support was provided by the HPCMO and the NERSC. 1 R.
Loutfy and E. Wexler, Perspectives of Fullerene Nanotechnology edited ¯ by E. Osawa (Kluwer Academic Publishers, New York, 2002), Part VII, pp. 275–280.
J. Chem. Phys. 141, 141101 (2014) 2 R.
A. DiStasio, O. A. von Lilienfeld, and A. Tkatchenko, Proc. Natl. Acad. Sci. U.S.A. 109, 14791 (2012). 3 C. A. Hunter and J. K. M. Sanders, J. Am. Chem. Soc. 112, 5525 (1990). 4 S. H. Patil and K. T. Tang, J. Chem. Phys. 106, 2298 (1997). 5 J. Mitroy and M. W. J. Bromley, Phys. Rev. A 68, 052714 (2003) (Erratum); 71, 019902 (2005) (Erratum). 6 S. J. A. van Gisbergen, J. G. Snijders, and E. J. Baerends, J. Chem. Phys. 103, 9347 (1995). 7 M. van Faassen, L. Jensen, J. A. Berger, and P. L. de Boeij, Chem. Phys. Lett. 395, 274 (2004). 8 A. D. Becke and E. R. Johnson, J. Chem. Phys. 127, 154108 (2007). 9 O. A. Vydrov and T. Van Voorhis, Phys. Rev. A 81, 062708 (2010). 10 A. Tkatchenko and M. Scheffler, Phys. Rev. Lett. 102, 073005 (2009). 11 A. Tkatchenko, R. A. DiStasio, R. Car, and M. Scheffler, Phys. Rev. Lett. 108, 236402 (2012). 12 S. N. Steinmann and C. Corminboeuf, J. Chem. Theory Comput. 7, 3567 (2011). 13 V. V. Gobre and A. Tkatchenko, Nat. Commun. 4, 2341 (2013). 14 A. Ruzsinszky, J. P. Perdew, J. Tao, G. I. Csonka, and J. M. Pitarke, Phys. Rev. Lett. 109, 233203 (2012). 15 J. Tao, J. P. Perdew, and A. Ruzsinszky, Proc. Natl. Acad. Sci. U.S.A. 109, 18 (2012). 16 J. Tao, J. P. Perdew, and A. Ruzsinszky, Phys. Rev. B 81, 233102 (2010). 17 A. A. Lucas, L. Henrard, and Ph. Lambin, Phys. Rev. B 49, 2888 (1994). 18 J. P. Perdew, J. Tao, P. Hao, A. Ruzsinszky, G. I. Csonka, and J. M. Pitarke, J. Phys.: Condens. Matter 24, 424207 (2012). 19 See supplementary material at http://dx.doi.org/10.1063/1.4897957 for detailed derivation. 20 A. Jiemchooroj, P. Norman, and B. E. Sernelius, J. Chem. Phys. 125, 124306 (2006). 21 J. M. Pacheco and J. P. P. Ramalho, Phys. Rev. Lett. 79, 3873 (1997). 22 A. Banerjee, J. Autschbach, and A. Chakrabaru, Phys. Rev. A 78, 032704 (2008). 23 S. Botti, A. Castro, X. Andrade, A. Rubio, and M. A. L. Marques, Phys. Rev. B 78, 035333 (2008). 24 L. Serra, F. Garcias, M. Barranco, J. Navarro, L. C. Balbas, A. Rubio, and A. Mananes, J. Phys.: Condens. Matter 1, 10391 (1989). 25 G. B. Adams, M. O’Keeffe, and R. S. Ruoff, J. Phys. Chem. 98, 9465 (1994). 26 R. R. Zope, T. Baruah, M. R. Pederson, and B. I. Dunlap, Phys. Rev. B 77, 115452 (2008). 27 G. K. Gueorguiev, J. M. Pacheco, and D. Tománek, Phys. Rev. Lett. 92, 215501 (2004). 28 J. Kauczor, P. Norman, and W. A. Saidi, J. Chem. Phys. 138, 114107 (2013). 29 V. A. Parsegian, Van der Waals Forces (Cambridge University Press, Cambridge, 2005). 30 H. C. Hamaker, Physica 4, 1058 (1937). 31 J. F. Dobson and T. Gould, J. Phys.: Condens. Matter 24, 073201 (2012). 32 S. Grimme, Comput. Mol. Sci. 1, 211 (2011). 33 There is a small additional contribution for C and C given by 21C MA 8 10 6 i=1 M M M j =1B δij2 /(MA MB ) and i=1A j =1B [126C6 δij4 + 36C8 δij2 ]/(MA MB ), re-
spectively. Here δij = |R − rij | and rij is the distance of atom pairs. F. Dobson, A. White, and A. Rubio, Phys. Rev. Lett. 96, 073201 (2006). 35 R.-F. Liu, J. G. Ángyán, and J. F. Dobson, J. Chem. Phys. 134, 114106 (2011). 36 J. F. Dobson, Int. J. Quantum Chem. 114, 1157 (2014). 37 A. J. Misquitta, J. Spencer, A. J. Stone, and A. Alavi, Phys. Rev. B 82, 075312 (2010). 38 L. Ma, K. A. Jackson, J. Wang, M. Horoi, and J. Jellinek, Phys. Rev. B 89, 035429 (2014). 39 A. Derevianko, S. G. Porsev, and J. F. Babb, At. Data. Nucl. Data Tables 96, 323 (2010). 40 R. Peverati and D. G. Truhlar, Philos. Trans. R. Soc. A 372, 20120476 (2014). 34 J.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.160.4.77 On: Sat, 20 Dec 2014 07:43:25
