Evaluation of New Density Functional with Account of van der Waals Forces by Use of Experimental H2 Physisorption Data on Cu(111) Kyuho Lee,1 André K. Kelkkanen,2, 3 Kristian Berland,3 Stig Andersson,4 David C. Langreth,1 Elsebeth Schröder,3 Bengt I. Lundqvist,2, 3, 5 and Per Hyldgaard3 arXiv:1109.0726v1 [cond-mat.mtrl-sci] 4 Sep 2011 1 Department of Physics and Astronomy, Rutgers University, Piscataway, NJ 08854-8019, USA 2 Center for Atomic-scale Materials Design, Department of Physics, Technical University of Denmark, DK-2800 Kongens Lyngby, Denmark 3 Department of Microtechnology and Nanoscience, MC2, Chalmers University of Technology, SE-41296 Göteborg, Sweden 4 Department of Physics, Göteborg University, SE-41296 Göteborg, Sweden 5 Department of Applied Physics, Chalmers University of Technology, SE-41296 Göteborg, Sweden (Dated: September 6, 2011) Detailed experimental data for physisorption potential-energy curves of H2 on low-indexed faces of Cu challenge theory. Recently, density-functional theory has been developed to also account for nonlocal correlation effects, including van der Waals forces. We show that one functional, denoted vdW-DF2, gives a potential-energy curve promisingly close to the experiment-derived physisorptionenergy curve. The comparison also gives indications for further improvements of the functionals. PACS numbers: 71.15.-m, 73.90.+f, 68.35.Np Density-functional theory (DFT) gives in principle exact descriptions of stability and structure of electron systems, but in practice approximations have to be made to describe electron exchange and correlation (XC) [1, 2]. Evaluation of XC functionals is commonly done by comparison with results from other accurate electronstructure theories or by comparing with relevant experiments, typically providing test numbers only for one or two measurables. This Letter illustrates the advantages of a third approach, which builds the experiment-theory calibration based on extensive experimental data, in this case a full physisorption potential derived from surfacephysics measurements. In the physisorption regime, resonant elastic backscattering-diffraction experiments provide a detailed quantitative knowledge. Here, data obtained for H2 and D2 on Cu surfaces [3–6] are used as a demanding benchmark for the performance of adsorbate potentialenergy curves (PECs) calculated with a nonempirical theory for extended systems. Density functionals that aspire to account for nonlocal electron-correlation effects, including van der Waals (vdW) forces, can then be assessed. In particular, we study the vdW-DF method [7–10] and show that calculations with versions of it provide a promising description of the physisorption potential for H2 on the Cu(111) surface and that the most recent one, vdW-DF2 [10], is more accurate than the first one [7–9] and other tested functionals. Sparse matter is abundant. Dense matter, also abundant, is since long successfully described by DFT. The recent extensions of DFT functionals to regions of low electron density, where the ubiquitous vdW forces are particularly relevant, render DFT useful also for sparse matter. In the vdW-DF functional the vdW interactions and correlations are expressed in the density n(r) as a truly nonlocal six-dimensional integral [7–9]. Its key in- gredients are (i) its origin in the adiabatic connection formula [11–13], (ii) an approximate coupling-constant integration, (iii) the use of an approximate dielectric function in a single-pole form, (iv) which is fully nonlocal and satisfies known limits, sum rules, and invariances, and (v) whose pole strength is determined by a sum rule and whose pole position is scaled to give the approximate gradient-corrected electron-gas ground-state energy locally. There are no empirical or fitted parameters, just references to general theoretical criteria. Like composite molecules, adsorption systems have electrons in separate molecule-like regions, with exponentially decaying tails in between. Then the slowly varying electron gas, used in the original vdW-DF method [7– 9, 14], might not be the most appropriate reference system for the gradient correction [15]. Although promising results have been obtained for a variety of systems, including adsorption [16, 17], there is room for improvements. For the recent vdW-DF2 functional, the gradient coefficient of the B88 exchange functional [18] is used for the determination of the internal functional [Eq. (12) of Ref. 7] within the nonlocal correlation functional. This is based on application of the large-N asymptote [19, 20] on appropriate molecular systems. Using this method, Elliott and Burke [21] have shown from first principles that the correct exchange gradient coefficient β for an isolated atom (monomer) is essentially identical to the B88 value, which had been previously determined empirically [18]. Thus in the internal functional, vdW-DF2 [10] replaces Zab in that equation with the value implied by the β of B88. This procedure defines the relationship between the kernels of vdW-DF and vdW-DF2 for the nonlocal correlation energy. Like vdW-DF, vdW-DF2 is a transferable functional based on physical principles and approximation and without empirical input. The vdW-DF method also needs to choose an over- 2 all exchange functional to obtain the exchange contribution to the interaction energy between two monomers (for example). The original vdW-DF uses the revPBE [22] exchange functional, good at separations in typical vdW complexes [7–9]. The latter choice can be improved on [23–27]. Recent studies suggest that the PW86 exchange functional [28] most closely duplicates Hartree-Fock interaction energies both for atoms [24] and molecules [25]. The vdW-DF2 functional [10] employs the PW86R functional [25], which more closely reproduces the PW86 integral form at lower densities than those considered by the original PW86 authors. Evaluation of XC functionals with respect to other theoretical results is often done systematically, e.g., by benchmarking against the S22 data sets [29–33]. The S22 sets have twenty-two prototypical small molecular duplexes for noncovalent interactions (hydrogen-bonded, dispersion-dominated, and mixed) in biological molecules and provides PECs at a very accurate level of wavefunction methods, in particular the CCSD(T) method. However, by necessity, the electron systems in such sets have finite size. The original vdW-DF performs well on the S22 dataset [29–33], except for hydrogen-bonded duplexes (underbinding by about 15% [10, 16]). The vdW-DF2 functional reduces the mean absolute deviations of binding energy and equilibrium separation significantly [10]. Shapes of PECs away from the equilibrium separation are greatly improved. The long-range part of the vdW interaction, particularly crucial for extended systems, has a weaker attraction in the vdW-DF2, thus reducing the error to 8 meV at separations 1 Å away from equilibrium [10]. Experimental information provides the ultimate basis for assessing functionals. The vdW-DF functional has been promising in applications to a variety of systems [16], but primarily vdW bonded ones, typically tested on binding-energy and/or bond-length values that happen to be available. The vdW-DF2 functional has also been successfully applied to some extended systems, like graphene and graphite [10], metal-organic-frameworks systems [34], molecular crystal systems [35], physisorption systems [36, 37], liquid water [38] and layered oxides [39]. However, those studies are of the common kind that focus on comparison against just a few accessible observations. The key step taken by the present work is to benchmark a full PEC in an extended system. Fortunately, for almost two decades, accurate experimental values for the eigenenergies of H2 and D2 molecules bound to Cu surfaces [3, 4] have been waiting for theoretical account and assessment. The rich data bank covers results for the whole shape of the physisorption potentials. The H2 -Cu system is particularly demanding for the vdW-DF and vdW-DF2 functionals and alike. On one hand, H2 is a small molecule with a large HOMO-LUMO gap, far from the low-frequency polarization modes as- sumed in derivation of these functionals [7–10]. On the other hand, Cu is a metal, which has created particular concerns [40]. Chemically inert atoms and molecules adsorb physically on cold metal surfaces [5], with characteristic low desorption temperatures ranging from a few K (He) to tens of K (say Ar and CH4 ). Adsorption energies, which take values from a few meV to around 100 meV, may be determined from thermal-desorption or isosteric-heatof-adsorption measurements. For light adsorbates, like He and H2 , gas-surface-scattering experiments, involving resonance structure of the elastic backscattering, provide a more direct and elegant method, with accurate and detailed measurements of bound-level sequences in the potential well. The availability of isotopes with widely different masses (3 He, 4 He, H2 , D2 ) permits a unique assignment of the levels and a determination of the well depth and ultimately a qualified test of model potentials [6]. The bound-level sequences of specific concern here were obtained using nozzle beams of para-H2 and normalD2 , i.e., the beams are predominantly composed of j = 0 molecules. This implies that the measured bound-state energies, ǫn (listed in Table I for H2 and D2 on Cu(111)), refer to an isotropic distribution of the molecular orientation. For this particular ordering, all ǫn values fall accurately on a common curve when plotted √ versus the mass-reduced level number η = (n + 1/2)/ m. This implies a level assignment that is compatible with a single gas-surface potential for the two hydrogen isotopes [4]. A third-order polynomial fit to the data yields for η = 0 a potential-well depth D = 29.5 meV. The experimental energy levels in the H2 -Cu(111) PEC (see Table I) may be analyzed [3–5] within the traditional theoretical picture [41, 42] of the interaction between inert adsorbates and metal surfaces: The PEC is then approximated as a superposition of a long-range vdW attraction, VvdW , and a short-range Pauli repulsion, VR , reflecting overlap between tails of metal Bloch electrons and the adparticle’s closed-shell electrons [5, 41]. This results in a laterally and angularly averaged potential Vo (z) = VR (z) + VvdW (z), where z is the normal distance of the H2 bond center from the jellium edge. The bound-level sequences in Table I can be accurately reproduced (< 0.3 meV) by such a physisorption potential [4, 5] (Fig. 1), having a well depth of 28.9 meV and a potential minimum located 3.50 Å outside the topmost layer of copper ion cores. From the measured intensities of the first-order diffraction beams, a very small lateral variation of the H2 -Cu(111) potential can be deduced, ∼ 0.5 meV at the potential-well minimum. A direct solution of the Schrödinger equation in Vo (z) reproduces the four low-energy eigenvalues to within 3% of the measured ones. It is therefore consistent [43] to benchmark against the very accurately constructed experimental physisorption curve in Fig. 1. 3 FIG. 1: Experimentally determined effective physisorption potential for H2 on Cu(111) [3], compared with potentialenergy curves for H2 on Cu(111), calculated for the atop site in GGA-revPBE, GGA-PBE, vdW-DF2, and vdW-DF. TABLE I: Sequences of bound-state energies for H2 and D2 on Cu(111). DFT eigenvalues are calculated with the vdWDF2 potential of Fig. 1, and experimental numbers are from [3, 4]. n ǫn [meV] H2 D2 DFT Exper. DFT Exper. 0 −32.6 −23.9 −34.4 1 −21.3 −15.5 −26.0 −19.0 2 −12.1 −8.7 −18.7 −12.9 3 −5.4 −5.0 −12.4 −8.9 4 −7.4 −5.6 5 −3.5 −3.3 Figure 1 shows our comparison of density functional PECs against the experimental physisorption potential for H2 in an orientation above an atop site on the Cu(111) surface. The PECs of Fig. 1 are calculated with the vdWDF [7] and vdW-DF2 functionals [10] as well as with two generalized gradient approximations (GGA-PBE and GGA-revPBE). We use an efficient vdW algorithm [44] adapted from SIESTA’s [45] vdW code within a modified version of the plane-wave code ABINIT [46]. The vdW interaction is treated fully self-consistently [9] (allowing also vdW forces to relax the adsorption geometry). The computational costs are the same with vdWDF and vdW-DF2. Our choice of Troullier-Martins-type norm-conserving pseudopotentials and a high cutoff energy (70 Ry) ensures excellent convergence; the GPAW code [47], in its default mode, gives similar but less accurate results [48]. We stress that there is neither a damping nor saturation function in vdW-DF and vdW-DF2 calculations. The need for an account of nonlocal corre- FIG. 2: Interaction potential for H2 on Cu(111), calculated self-consistently with the vdW-DF2 functional [10] in the bridge, hollow and atop sites. lations for the description of vdW forces is illustrated by the GGA curves giving inadequate PECs. The calculated well depth in vdW-DF, 53 meV, should be compared with the measured one, 29.5 meV [4], and the one calculated from Vo (z), suitably parametrized, 28.9 meV [3] (Fig. 1). We find that the vdW-DF2 PEC lies close to the experimental physisorption potential, both at the equilibrium position and at separations further away from the surface. Calculated PECs of H2 above bridge, atop, and center sites on the Cu(111) surface are shown in Fig. 2, their closeness illustrating the lack of corrugation on this surface, as in experimental findings [3, 4]. Similar results for bridge, atop, and center sites on Cu(100) show the PECs on the two surfaces to be very close to each other, just like for the experimental result [3, 4]. The vdW-DF2 equilibrium separation is about 3.5 Å, like the value that is deduced from the experimental data and for reasonable physical assumptions about the parameters [3, 4]. A further refined comparison is provided by the boundstate eigenvalues for the point-of-gravity motion of H2 on Cu(111), both experimental data and those from the vdW-DF2 potential (Table I). In addition to PEC shapes and eigenenergy values, there should be comparison with values for well depth and equilibrium separation, experimental ones being (29.5 meV; 3.5 Å) and those of the vdW-DF2 potential in Fig. 1 (37 meV; 3.5 Å). Viewing the facts that (i) vdW-DF2 is a first-principles method, where characteristic electron energies are typically in the eV range, and (ii) the test system and results are very demanding, as other popular methods deviate significantly more from the experimental curve (for instance, application of the DFT-D3(PBE) method [49], with atompairwise specific dispersion coefficients and cutoff radii computed from first principles, gives (−88 meV; 2.8 Å) 4 for the PEC minimum point). We judge this as very promising. So is the relative closeness of experimental and calculated eigenenergy values in Table I. The discrepancies between the eigenvalues signal that the vdW-DF2 PEC might not have the right shape for H2 on Cu(111). Reference 10 shows that vdW-DF2 benchmarks very well against the S22 data sets. It is possible that the metallic nature of the H2 /Cu(111) system causes modifications in the details of the electrodynamical response; H2 /Cu(111) physisorption constitutes a very strong challenge for the density functional. For a well established conclusion, a more accurate theory is called for. In summary, accurate and extensive experimental data for the physisorption PEC of H2 on Cu(111) are used to evaluate vdW-DF functionals and their adequacy for metal surfaces. The Cu(111) surface is chosen here, as its flatness gives clarity in the analysis and eliminates several side-issues that could have made interpretations fuzzier. More generally, there exists an accurate data bank of experimental physisorption information that challenges every density functional approximation to produce relevant PECs. We propose that such surface-related PEC benchmarking should find a broader usage, supplementing, for example, S22 comparisons as an accelerator in density-functional development. Several qualitative similarities are found for both vdWDF and vdW-DF2 functionals. The vdW-DF2 functional gives PECs in a useful qualitative and quantitative agreement with the experimental PEC, i.e. with respect to well depth, equilibrium separation, and curvature of PEC near the well bottom, and thus zero-point vibration frequency. This is very promising for applications of this nonlocal correlation functional at short and intermediate separations, relevant for the adsorption. However, the accuracy of experimental data is high enough to stimulate a more detailed analysis of all aspects of the theoretical description. This should be valuable for the further XC-functional development. The Swedish National Infrastructure for Computing (SNIC) is acknowledged for providing computer allocation and the Swedish Research Council for providing support to ES and PH. AK and BIL thank the Lundbeck foundation for sponsoring the center for Atomic-scale Materials Design and the Danish Center for Scientific Computing for providing computational resources. Work by KL and DCL is supported by NSF DMR-0801343. Professor David Langreth was active in all aspects of the research until his untimely death in May 2011. We would like to express our sense both of personal loss and of loss to our discipline. [1] P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964). [2] W. Kohn and L.J. Sham, Phys. Rev. 140, A1133 (1965). [3] S. Andersson and M. Persson, Phys. Rev. Lett. 70, 202 (1993). [4] S. Andersson and M. Persson, Phys. Rev. B 48, 5685 (1993). [5] See, e.g., M. Persson and S. Andersson, Chapter 4, “Physisorption Dynamics at Metal Surfaces”, in Handbook of Surface Science, Vol. 3 (Eds. E. Hasselbrink and B.I. Lundqvist), Elsevier, Amsterdam (2008), p. 95. [6] R. J. Le Roy, Surf. Sci. 59, 541 (1976). [7] M. Dion, H. Rydberg, E. Schröder, D.C. Langreth, and B.I. Lundqvist, Phys. Rev. Lett. 92, 246401 (2004) and 95, 109902(E) (2005). [8] D.C. Langreth, M. Dion, H. Rydberg, E. Schröder, P. Hyldgaard, and B.I. Lundqvist, Int. J. Quant. Chem. 101, 599 (2005). [9] T. Thonhauser, V.R. Cooper, S. Li, A. Puzder, P. Hyldgaard, and D.C. Langreth, Phys. Rev. B 76, 125112 (2007). [10] K. Lee, É.D. Murray, L. Kong, B.I. Lundqvist, and D.C. Langreth, Phys. Rev. B (RC) 82, 081101 (2010). [11] D.C. Langreth and J.P. Perdew, Solid State Commun. 17, 1425 (1975). [12] O. Gunnarsson and B.I. Lundqvist, Phys. Rev. B 13, 4274 (1976). [13] D.C. Langreth and J.P. Perdew, Phys. Rev. B 15, 2884 (1977). [14] H. Rydberg, B.I. Lundqvist, D.C. Langreth, and M. Dion, Phys. Rev. B 62, 6997 (2000). [15] D.C. Langreth and S.H. Vosko, in “Density Functional Theory of Many-Fermion Systems”, ed. S.B. Trickey, Academic Press, Orlando, 1990. [16] D.C. Langreth, B.I. Lundqvist, S.D. Chakarova-Käck, V.R. Cooper, M. Dion, P. Hyldgaard, A. Kelkkanen, J. Kleis, L. Kong, S. Li, P.G. Moses, E. Murray, A. Puzder, H. Rydberg, E. Schröder, and T. Thonhauser, J. Phys.: Cond. Mat. 21, 084203 (2009). [17] Y.N. Zhang, F. Hanke, V. Bortolani, M. Persson, and R. Q. Wu, Phys. Rev. Lett. 106, 236103 (2011). [18] A.D. Becke, Phys. Rev. A 38, 3098 (1988). [19] J. Schwinger, Phys. Rev. A 22, 1827 (1980). [20] J. Schwinger, Phys. Rev. A 24, 2353 (1981). [21] P. Elliott and K. Burke, Can. J. Chem. 87, 1485 (2009). [22] Y. Zhang and W. Yang, Phys. Rev. Lett. 80, 890 (1998). [23] A. Puzder, M. Dion, and D.C. Langreth, J. Chem. Phys. 126, 164105 (2006). [24] F. O. Kannemann and A. D. Becke, J. Chem. Theory Comput. 5, 719 (2009). [25] É.D. Murray, K. Lee, and D.C. Langreth, Jour. Chem. Theor. Comput. 5, 2754 (2009). [26] J. Klimeš, D.R. Bowler, and A. Michaelides, J. Phys.: Condens. Matter 22, 022201 (2010). [27] V.R. Cooper, Phys. Rev. B 81, 161104(R) (2010). [28] J.P. Perdew and Y. Wang, Phys. Rev. B 33, 8800(R) (1986). [29] P. Jurečka, J. Šponer, J. Černý, and P. Hobza, Phys. Chem. Chem. Phys. 8, 1985 (2006). [30] D. Sherrill, T. Takatani, and E. G. Hohenstein, J. Phys. 5 Chem. A 113, 10146 (2009). [31] L. F. Molnar, X. He, B. Wang, and K. M. Merz, J. Chem. Phys. 131, 065102 (2009). [32] T. Takatani, E.G. Hohenstein, M. Malagoli, M.S. Marshall, and C.D. Sherrill, J. Chem. Phys. 132, 144104 (2010). [33] R. Podeszwa, K. Patkowski, and K. Szalewicz, Phys. Chem. Chem. Phys. 12, 5974 (2010). [34] L.Z. Kong, G. Román-Pérez, J.M. Soler, and D.C. Langreth, Phys. Rev. Lett. 103, 096103 (2009). [35] K. Berland, Ø. Borck, and P. Hyldgaard, Comp. Phys. Commun. 182, 1800 (2011). [36] K. Lee, Y. Morikawa, and D.C. Langreth, Phys. Rev. B 82, 155461 (2010). [37] J. Wyrick, D.-H. Kim, D. Sun, Z .Cheng, W. Lu, Y. Zhu, K. Berland, Y.S. Kim, E. Rotenberg, M. Luo, P. Hyldgaard, T.L. Einstein, and L. Bartels, Nano Letters 11, 2944 (2011). [38] A. Møgelhøj, A. Kelkkanen, K.T. Wikfeldt, J. Schiøtz, J.J. Mortensen, L.G.M. Pettersson, B.I. Lundqvist, K.W. Jacobsen, A. Nilsson, and J.K. Nørskov, J. Phys. Chem. B (DOI: 10.1021/jp2040345), in print. [39] E. Londero and E. Schröder, Computer Phys. Commun. 182, 1805 (2011). [40] J.F. Dobson, Surf. Sci. 605, 1621 (2011) and references therein. [41] E. Zaremba and W. Kohn, Phys. Rev. B 15, 1769 (1977). [42] P. Nordlander and J. Harris, J. Phys. C 17, 1141 (1984). [43] We note that VvdW (z) in the potential Vo (z) involves a prefactor f (z), which introduces a saturation of the attraction at atomic-scale separations. The function f (z) [44] [45] [46] [47] [48] [49] [50] lacks a rigorous prescription and this results in a level of arbitrariness to VvdW (z). It should be noted, though, that the saturation function f (z) of Refs. 3–5, 42 has a formal resemblance to the damping functions introduced to adjust for double counting in semi-empirical DFT descriptions [49] of dispersive interactions. We emphasize that f (z) of VvdW (z) is likely to work across several length scales [50] because VvdW (z) includes the effects of image planes and is consistent with the Zaremba-Kohn formulation of physisorption [41]. G. Román-Pérez and J.M. Soler, Phys. Rev. Lett. 103, 096102 (2009). P. Ordejón, E. Artacho, and J.M. Soler, Phys. Rev. 53, 10441(R) (1996); J.M. Soler, E. Artacho, J.D. Gale, A. Garcı́a, J. Junquera, P. Ordejón, and D. Sánchez-Portal, J. Phys.: Condens. Matter 14, 2745 (2002). X. Gonze, J.-M. Beuken, R. Caracas, F. Detraux, M. Fuchs, G.-M. Rignanese, L. Sindic, M. Verstraete, G. Zerah, F. Jollet, M. Torrent, A. Roy, M. Mikami, Ph. Ghosez, J.-Y. Raty, and D.C. Allan, Comput. Mater. Sci. 25, 478 (2002). J.J. Mortensen, L.B. Hansen, and K.W. Jacobsen, Phys. Rev. B 71, 035109 (2005). A. Kelkkanen, Implementation of van der Waals Forces in DFT for Nanoelectronics and Water Structure, Thesis (CAMd, DTU, 2011). S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, J. Chem. Phys. 132, 154104 (2010). K. Berland and P. Hyldgaard, J. Chem. Phys. 132, 134705 (2010).