ARTICLE pubs.acs.org/JPCA
Density Functional Calculation of Intermolecular Potentials Carl Nyeland* Institute of Chemistry, University of Copenhagen, Universitetsparken 5, DK 2100 Copenhagen Ø, Denmark ABSTRACT: Calculations of intermolecular potentials following the density functional theory (DFT) turn out to be very complicated without using some appropriate approximations. Most often the following three approximations have been considered. In one approximation the disturbed charge distributions during collisions are reduced to sums of undisturbed charge distributions from the colliding species. In another approximation, the so-called local density approximation (LDA), one neglects the fact that the intermolecular potentials that depend on charge densities also depend on gradients in the densities. In a third approximation one assumes that the intermolecular potential can be considered as a sum of two terms: a term for the longrange geometry and a term for the short-range geometry. In this Article the three approximations mentioned will be discussed for numerical accuracy for calculations of potentials between inert gas atoms and for calculations of potentials between surfaces and inert gas atoms. In the discussion a few other approximations will be mentioned too.
I. INTRODUCTION To be able to work out reasonably accurate calculations of intermolecular potentials following the density functional theory (DFT),1 one needs detailed descriptions of the deformations of charge distributions during the molecular interactions considered. This information could be obtained by appropriate DFT calculations, for instance from variation calculations, but is considered to be obtained from more simple approaches in the present Article. Anyway, such data are only scarcely available from the literature. Among others some results for interactions of two helium atoms from a HartreeFock calculation due to Gordon et al.2 will be considered. As mentioned from DFT calculations, structure dependent charge distributions could in principle be obtained by the use of the two well-known theorems of Hohenberg and Kohn,3 but a detailed method for numerical calculations of intermolecular potentials does not seem available at present. The two theorems will be mentioned briefly to make the approximations used below more acceptable: From the first HohenbergKohn theorem it follows that the energy of the ground state of a molecule or a van der Waals complex, for instance described as a BornOppenheimer state, is given as a general functional of the configuration-dependent charge distribution Fo(r). Formally this theorem can be written as a functional dependence E½Fo
ð1Þ
The second HohenbergKohn theorem provides the energy variation principle which can be written E½F > E½Fo
ð2Þ
where Fo(r) means the correct structure dependent charge distribution for the ground state; F(r) is the test function for the charge distribution that is different from the correct one, but it is preferable to have the correct symmetry. r 2011 American Chemical Society
It is important to mention that the functional E[F] is not given from the HK theorems or from the basic theory of Hohenberg and Kohn. It has to be obtained from other sources, in principle from the Schr€odinger equation or from analyses of other types of results, experimental or theoretical. A simple approximation that has often been considered for the energy funtional following the DFT is the so-called local density approximation (LDA). Here the energy functional is simplified to be a function of the charge density distribution, if possible the correct one, but not a functional of this distribution as considered in the proper DFT procedure. By using the ThomasFermiDirac (TFD)1 form without introducing further approximations, one gets in the LDA limit E½F T TF ðFÞ þ V coul ðFÞ þ V exch ðFÞ þ V corr ðFÞ
ð3Þ
where a general correlation term Vcorr by definition is added to the TFD form considered as the first three terms of eq 3 to obtain the proper LDA result. The general LDA correlation function was discussed in ref 1, especially for spinless systems. The basic theories used here are from refs 46. In ref 4 calculations for unpolarised electron gases as functions of the density F were presented, which in ref 5 were used to formulate a general function for the correlation energy. The radius of the WignerSeitz sphere rs, defined by F1 = 4/3πrs3, was used most often in the tabulations and formulas in ref 1. Here this function has been used for the correlation energy in the LDA limit. The higher-order results concerning the gradient expansion have been considered seriously during the last years.6 Some results for intermolecular potentials are given in Appendix 2 in Special Issue: J. Peter Toennies Festschrift Received: December 10, 2010 Revised: April 28, 2011 Published: May 11, 2011 6888
dx.doi.org/10.1021/jp111757u | J. Phys. Chem. A 2011, 115, 6888–6891
The Journal of Physical Chemistry A
ARTICLE
this Article. Equation 3 using the TFD functions was also the basis of the so-called electron-gas method suggested and used by Gordon and Kim in 1972.7 In the years since, the correlation term has been more well-defined and corrections for selfexchange and self- correlation for systems having a finite numbers of electrons have been added to the calculation procedure.8 When no useful information about disturbed charge distributions is available, a simple approximation to the disturbed charge distribution FAB for a van der Waals complex AB has often been considered. In this approximation the simple use of a summation of undisturbed charge distributions is used as discussed for instance in refs 911 FAB FA þ FB
Table 1. Potential Parameters for the Laterally Averaged of the Inert Gas AtomSurface Potential Given in eq 9 system
UMAa
aMAa
CvdWb
zvdWb
HeCu(001)
1.486
1.234
0.0562
0.32
NeCu(001) ArCu(001)
1.923 2.704
1.220 1.137
0.112 0.382
0.31 0.38
KrCu(001)
1.773
1.057
0.538
0.40
a
The results shown are obtained from exponential plots to the LDA results following eq 6. b vdW-data taken from ref 20.
ð4Þ
In Appendix 1 in the present Article some calculations on a few systems of interacting inert gas atoms based on disturbed distributions are considered to obtain more accurate DFT results in the LDA limit for intermolecular potentials. These results are compared with results obtained earlier8 by using the approximation of a sum of undisturbed charge distributions (eq 4). The approximation of additivity of charge distributions was discussed by Jorgensen and Allan,9 by Bone and Bader,10 and by Bentley.11,12 Gradient corrections to the LDA should be taken into account as discussed, for instance in ref 13. The numerical meaning of these corrections for interacting inert gas atoms will be considered in Appendix 2. In the main part of this Article a LDA method for inert gas atommetal surface potentials will be considered. Atomic units are used: 1 au length = 0.529 17 Å and 1 au energy = 27.21 eV.
II. CALCULATIONS OF INERT GAS ATOMMETAL SURFACE POTENTIALS FOLLOWING A DFT/LDA METHOD It seems useful to consider results of inert gas atommetal surface potentials from a calculation following a DFT method,1,14 but even in the LDA limit of this method one encounters further complications if one wants to consider all intermolecular distances, short-range distances as well as long-range distances, in a “seamless” calculation. Discussions of “seamless” methods for intermolecular potentials have been considered in the literature recently1517 mainly for atomatom potentials. The methods turn out to be very complicated to use. The approximation of using a summation of undisturbed charges in the short-range potentials was shown lately to give reasonable results,18 but as has often been discussed, the long-range part of the potentials is given from polarizations of the interacting charges and not even approximatively from the undisturbed charge distributions. As mentioned, an easy approximative method for calculations of short-range atomsurface potentials was considered, but to obtain full range intermolecular potentials, one needs to add results from short-range calculations to results from long-range calculations, possibly corrected for damping, as for instance discussed in ref 19. First some results from calculations of the inert gas atom metal surface short-range potentials following the method of ref 18 are given here. New numerical results for NeCu(001), ArCu(001), and KrCu(001) are presented together with the earlier published results for HeCu(001). Definitions of the quantities involved are given in the earlier publication18 concerning the LDA approximation. The short-range repulsive
Figure 1. Solid surfaceatom interaction model. ZA is the distance between the atom A and the topmost plane of surface atoms. R is the distance between the atom A and a specific surface atom and zA between the atom A and the jellium edge in the atomjellium interaction model. zb is the distance between the surface plane and the jellium edge. It is seen that ZA = zb þ zA.
potential V is given by the following sum of density functionals V ½F ¼ V kin ½F þ V coul ½F þ hexch V exch ½F þ hcorr V corr ½F ð5Þ where Vkin is the contribution to the intermolecular potential from the ThomasFermi kinetic energy, Vcoul is the classical electrostatic energy, Vexch is the Dirac exchange energy contribution, and Vcorr is the contribution from the local density (LDA) correlation energy. These results are based on charge distributions for the separate systems as discussed above and for which numerical values are given in refs 8 and 18. Approximate treatments of the correction factors for exchange and correlation are given in ref 18. The calculation procedure used for the shortrange atomatom potentials was exactly the same as in the earlier work. From the undisturbed charge distributions the DFT/LDA terms of eq 5 are calculated for interacting surface systems with inert gas atoms and for undisturbed systems from which the short-range potentials are obtained as differences. Also for the systems considered here we have realized that the resulting short-range atomatom potentials have approximately simple exponential forms aMA R V MA rep ðRÞ U MA e
ð6Þ
where R is the interatomic distance. The best fit parameters UMA and aMA were obtained from the method discussed in ref 18 and listed in Table 1 of this Article for atomatom potentials corresponding to the systems mentioned. As discussed in ref 18, Celli21 has shown that an exponential potential can be 6889
dx.doi.org/10.1021/jp111757u |J. Phys. Chem. A 2011, 115, 6888–6891
The Journal of Physical Chemistry A
ARTICLE
Table 2. Intermolecular Potential Results of the Minimum Data Zm and Vm for Different Systems, Using ZA as Distance Parametera Zm (au)
Table 3. Corrections to LDA Results Due to Disturbtions in Charge Distributionsa system
Vm (meV)
HeHe
eq 9 b ref 25, expc ref 22, theod eq 9 b ref 25, expc ref 22, theod HeCu(001) 8.6
8.6
7.3
3.25
6.3
5.9
NeCu(001) 8.0
8.1
6.6
7.97
12.4
15.8
ArCu(001) 7.2
7.8
6.3
33.8
54.6
61.2
KrCu(001) 6.8
7.7
6.0
52.9
78.5
96.6
a
The full potential V(ZA) for the inert gas atommetal surface is defined in eq 9. The jellium edge zb = 1.704 au. b Results obtained by use of data from Table 1. c Results from ref 25, mostly experimental. d Theoretical results from ref 22.
ArAr
ΔF/FHF
ΔV/Vsum
3.000
0.060
0.099d
7.098
0.0724
R
b
0.144
c
Results for some inert gasinert gas atom potentials are shown. ΔF = FHF Fsum, and ΔV = VLDA(FHF) VLDA(Fsum). “HF” means “HartreeFock”, and “sum” means summation of undisturbed charge distributions from the colliding atoms. b Data for ΔF and FHF from ref 2. c Data for ΔF and FHF from ref 11. d Data obtained by the LDA method following ref 8 using the charge distributions shown. a
Table 4. Intermolecular Potential LDA Results Due to Zeroand Second-Order Gradient Expansiona
summed in a straightforward, but approximate, way for the atoms in a single monolayer and the interacting atom to give the following analytical expression for the laterally averaged potential between the solid surface and the interacting atom (Figure 1),
R
V0 (DF), zero orderb
V0 (DF), second order includedc ArAr
3
0.113(1)
ð7Þ
4
0.125(0)
0.101(0)
5
0.249(1)
0.208(1)
Here ZA is the distance between the inert gas atom and a plane passing through the monolayer surface considered, and AC is the surface area per surface atom. Now that the repulsive average potential has been calculated, it is possible approximately to predict the full potential curve, by adding a term for the long-range dispersion potential for instance as discussed in ref 22.
6
0.410(2)
0.555(2)
3
0.180(1)
4
0.196(0)
0.163(0)
5
0.488(1)
0.400(1)
6
0.106(1)
0.107(1)
3
0.385(1)
0.372(1)
4
0.107(1)
0.101(1)
5 6
0.108(0) 0.302(1)
0.871(1) 0.259(1)
V SA rep ðZA Þ
2π ¼ ð1 þ aMA ZA ÞU MA eaMA ZA aMA 2 AC
f ðbðzA zvdW ÞÞ
CvdW ðzA zvdW Þ3
2π ð1 þ aMA ZA ÞU MA eaMA ZA aMA 2 AC CvdW f ðbðZA zb zvdW ÞÞ ðZA zb zvdW Þ3
0.171(1)
XeXe
a The systems ArAr, KrKr, and XeXe are considered. b V0 (DF) from ref 8. c V0 (DF) obtained by use of gradient correction data from Tables 2 and 3 of ref 13.
where b is set to the short-range value aMA. For different atommetal surface systems calculated here results for the locations of the approximate potential minima, the distance Zm and the potential Vm, are given in Table 2 for the potential model (eq 9). The results are obtained by use of potential parameters given in Table 1. The “area of sites” AC for the (001) surface for the face-centered cubic cell of Cu were obtained from values of the lattice constant,24 AC = 0.5d2 = 23.2 au for d = 6.818 au. For comparison, some mostly experimental results from ref 25 and some theoretical results from ref 22 are also given in Table 2. The results obtained here for the HeCu(001) systems do not seem very good, which might follow from the neglect of higher-order terms due to gradient expansions as discussed in Appendix 2, but results for other systems obtained here seem reasonable.
ð9Þ
The damping function f(x), which is not very well-known for atomsolid interactions, was taken as the one considered in ref 23 based on the better known damping function for atomatom interactions. Here it is approximated as f ðxÞ ¼ 1 ð1 þ x þ 0:5x2 Þex
KrKr
ð8Þ
Here f(x) is a damping function, for instance the one considered in ref 23, b is a parameter most often taken as the exponential parameter from the short-range potential considered, while CvdW and zvdW are the well-known van der Waals dispersion potential constants.20 The modified dispersion potential eq 8 is traditionally given as a function of zA, the distance to the electron surface. If the location of the jellium edge zb is estimated from the lattice constant d as zb = 0.25d due to the face-centered cubic unit cell for Cu(001), a smooth noncorrugated potential V(ZA), where both the long-range attractive part (eq 8) and the short-range repulsive part (eq 7) can be obtained as functions of ZA, the distance between the gas atom and the top-layer of the solid. By using the relationship for the definition of the jellium edge, zb = ZA zA (Figure 1), one can consider the long-range potential of eq 8 to be a function of ZA too. Finally, the following full potential result is obtained using ZA as distance parameter for both potential terms V ðZA Þ ¼
0.106(1)
ð10Þ
III. CONCLUSION In conclusion it should be mentioned that the approximation in using a sum of undisturbed charge distributions gives errors around 10% for inert gas atominert gas atom interactions. Similar numerical errors are obtained for neglecting derivations in the LDA method due to gradient expansions. Numerical 6890
dx.doi.org/10.1021/jp111757u |J. Phys. Chem. A 2011, 115, 6888–6891
The Journal of Physical Chemistry A estimates were not obtained for approximations assuming summations of parts of intermolecular potentials. Here a so-called “seamless” method would be the correct procedure, but no useful procedure seems available in the literature presently. Finally, it should be mentioned that the three approximations discussed above much reduce the calculation time for intermolecular potentials when used as sketched even though they contribute to some minor errors as estimated.
’ APPENDIX 1. INERT GAS ATOMINERT GAS ATOM POTENTIAL First results based on a detailed calculation of the disturbed charge distribution for HeHe from ref 2 will be discussed. In ref 2,Gordon et al.show that the distortion at the center of interatomic geometry is approximately ΔF ≈ 0.0025 au for an intermolecular distance of R = 3 au, compared with results obtained by using the method of summation of undisturbed charge distributions, eq 4. Following the method in ref 8 it is easy from this result to calculate that the short-range potential in the LDA limit will be reduced by approximately 10% due to the distribution just mentioned. A similar result is obtained for the system ArAr from distribution data from Bentley11 (Table 3). Only a few further results can be taken from the literature. Jorgensen and Allan have obtained ΔF ≈ 0.020 au for HeHe (R = 2 au).9 Bentley has realized a distortion of 2.5% in the HeHF system around the so-called BCP configuration.12 In summary, the approximation of using undisturbed charges for interaction of inert gas atoms seems acceptable in detailed calculations of intermolecular collisions, but more complicated methods might be necessary for more complicated systems.
ARTICLE
(4) Ceperley, D. M.; Alder, B. J. Phys. Rev. Lett. 1980, 45, 566. (5) Vosko, S. H.; Wilk, L.; Nusair, M. Can. J. Phys. 1980, 58, 1200. (6) Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.; Singh, D. J.; Fiolhais, C. Phys. Rev. B 1992, 46, 6671. Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865. See also:Tao, J.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Phy. Rev. Lett. 2003, 91, 146401. (7) Gordon, R. G.; Kim, Y. S. J. Chem. Phys. 1972, 56, 3122. (8) Nyeland, C.; Toennies, J. P. Chem. Phys. 1994, 188, 205. (9) Jorgensen, W. L.; Allen, L. C. J. Am. Chem. Soc. 1971, 93, 567. (10) Bone, R. G. A.; Bader, R. F. W. J. Phys. Chem. 1996, 100, 10892. (11) Bentley, J. J. Phys. Chem. 1998, A102, 6043. (12) Bentley, J. J. Phys. Chem. 2000, A104, 9630. (13) Nyeland, C. Chem. Phys. Lett. 2003, 370, 353. (14) Kohn, W.; Becke, A. D.; Parr, R. G. J. Phys. Chem. 1996, 100, 12974. (15) Kohn, W.; Meir, Y.; Makarov, D. E. Phys. Rev. Lett. 1998, 80, 4153. (16) Dion, M.; Rydberg, H.; Schr€oder, E.; Langreth, D. C.; Lundqvist, B. I. Phys. Rev. Lett. 2004, 92, 246401. (17) Furche, F.; Van Voorhis, T. J. Chem. Phys. 2005, 122, 164106. (18) Nyeland, C.; Toennies, J. P. Chem. Phys. 2006, 321, 285. (19) Patil, S. H.; Tang, K. T.; Toennies, J. P. J. Chem. Phys. 2002, 116, 8118. (20) Chizmeshya, A.; Zaremba, E. Surf. Sci. 1989, 220, 443. (21) Celli, V. In Helium Atom Scattering from Surfaces; Hulpke, E., Ed.; Springer Series in Surface Science, Vol.27; Springer: Berlin, 1992. (22) Chizmeshya, A.; Zaremba, E. Surf. Sci. 1992, 268, 432. (23) Tang, K. T.; Toennies, J. P. Surf. Sci. Lett. 1992, 279, L203. (24) See for instance: Ziman, J. M. Electrons and Phonons; Clarendon Press: Oxford, U.K., 1960. (25) Vidali, G.; Ihm, G.; Kim, H.-Y.; Cole, M. W. Surf. Sci. Rep. 1991, 12, 133.
’ APPENDIX 2. GRADIENT EXPANSIONS IN DENSITY FUNCTIONAL CALCULATIONS OF INTERMOLECULAR POTENTIALS A discussion of the topic shown for this Appendix was given to second-order in ref 13 for interactions between atoms of inert gases. A few results for the gradient expansions connected to the LDA method are given in Table 4, considered as estimations to LDA results due to second-order density gradients for the kinetic energy and the exchange energy. Also results for the LDA to zeroorder are presented. The higher-order results obtained for HeHe and NeNe using only up to second-order seems misleading for calculations of intermolecular potentials. Higher-order terms should be taken into account, but theories for these terms are not avaible in the literature. Results for the heavier systems are given to second-order too. For these systems it seems clear that the second-order corrections reduce the results for the LDA by about 1020%. ’ AUTHOR INFORMATION Corresponding Author
*Address: Madvigs Alle 16, DK-1829 Frederiksberg C, Denmark. E-mail:
[email protected]. Fax: þ45 35 32 03 22.
’ REFERENCES (1) Parr, R. G.; Yang, W. Density-Functional Theory of Atoms and Molecules; Oxford University Press: New York, 1989. (2) Gordon, M. D.; Secrest, D.; Llaguno, C. J. Chem. Phys. 1971, 55, 1046. (3) Hohenberg, P.; Kohn, W. Phys. Rev. 1964, 136, B864. 6891
dx.doi.org/10.1021/jp111757u |J. Phys. Chem. A 2011, 115, 6888–6891
