Benchmarking van der Waals functionals with noncontact RPA calculations on graphene-Ag(111)

We have benchmarked long range behavior of seven different van der Waals functionals comparing them with our ACF-RPA correlation calculations for grap

The isolation of graphene in 2004 by Novoselov et al. [1,2] initiated widespread research of this unique material. Although graphene appeared in laboratories decades earlier, these experiments revealed graphene's most prominent feature, its electronic structure. It is so special and has so many possible applications because of π and π * bands that touch in the so-called Dirac point exactly at the Fermi level where in addition, those bands show linear dispersion. After Novoselov's ground-breaking experiments graphene has been grown on different substrates, metallic and insulating, employing different techniques. Graphene grown on metallic substrates is interesting from several points of view. Graphene-metal contact is essential for future graphene nanoelectronic devices; for example, for measuring graphene transport properties. Also from the practical point of view it is relatively easy to grow good quality graphene on some metallic substrates such as copper, and then transfer it to the insulating or polymer support, leading to interesting applications [3]. Characterization of graphene adsorbed on metallic surfaces showed that some surfaces such as Pd(111), Co(0001), and Ni(111) significantly influence the electronic structure of graphene. On the other hand, adsorption of graphene on (111) surfaces of Cu, Au, Pt, Ir, and Ag has almost no effect on its band structure [4,5]. Understanding these phenomena in terms of density functional theory (DFT), the standard tool in theoretical solid state physics, is not easy. Although the simplest local density approximation (LDA) correctly predicted strong interaction of Pd(111), Co(0001), and Ni(111) surfaces with graphene and weak binding of graphene to (111) surfaces of Cu, Au, Pt, and Ag [6,7], more accurate generalized gradient approximation (GGA) leads to no binding at all [8]. However, there is no apparent reason why LDA would give more accurate results than GGA, especially in inhomogeneous systems like this one. The reason for the failure of local LDA and semilocal GGA approximations in describing graphene-metal systems is their inability to describe nonlocal dispersive forces (particularly van der Waals force) which are important in weakly bonded systems. * ivor.loncaric@gmail.com † vito@phy.hr To overcome these limitations of LDA and GGA approaches, the van der Waals density functional (vdW-DF) has been developed [9,10]. Applying vdW-DF to a graphene-metal system turned out to be partially successful [11,12]. Although interaction at larger distances is improved significantly, the main problems were the inability of vdW-DF to reproduce strong binding of graphene on (111) Ni and Pd surfaces and overestimation of equilibrium distances of weakly adsorbed graphene. Overestimation of the equilibrium distances, and overbinding in the intermediate separations also occurred in other weakly bounded systems. A possible solution to the first problem is to use less repulsive exchange functional than originally proposed revPBE [13]. There are several propositions in this direction. First Klimeš et al. [14] introduced a slightly modified PBE [15] and B88 [16] exchange functionals in 2010 (abbr. optPBE and optB88) and optB86b [17] in 2011. Secondly Cooper [18] introduced a new exchange functional (abbr. C09x) constructed by matching to different s = |∇n|/2k F n limits (n is charge density). In the limit of small s (slow variation/high density) the C09x functional is matched to gradient expansion approximation, and in the limit of high s it is matched to the revPBE functional. Both propositions improve equilibrium distances. The same group that developed vdW-DF tried to resolve both problems by introducing a second version of van der Waals functional, vdW-DF2 [19]. They changed originally used revPBE exchange with slightly modified PW86 [20,21] exchange and used large-N asymptote gradient correction to construct the vdW kernel. Most recently Hamada [22] developed another exchange functional to be used with the vdW-DF2 functional. There are also other vdW functionals inspired by the vdW-DF idea, e.g., by Vydrov and Van Voorhis [23][24][25]. Physical correctness of these functionals was subject to debate [26,27]. Faster implementation of the Vydrov and Van Voorhis VV10 functional is due to Sabatini et al. [28] (abbr. rVV10).
With so many options there is still a need for benchmarking these functionals. Here we will benchmark their long range behavior comparing them to our noncontact random phase approximation (RPA) calculations on the example of graphene on Ag(111). Graphene on Ag(111) was grown very recently and shown to be the weakest bounded so far [29] making it a good candidate for our purpose. To show that our conclusions are not a coincidence, that only applies to the Ag(111), we also performed the calculations with some of the functionals for Cu(111) and Al(111) surfaces.
Neutral objects such as atoms, molecules, or metallic surfaces in their LDA ground state produce the short range electrical field whose range, if they are not polar in their ground state, corresponds to the range of its ground state electronic density. Consequently, those objects separated such that their ground state electronic densities do not overlap will not interact. However, because of spontaneous long range fluctuations in the considered bodies, they get polarized and produce an electrical field with a range much longer than the range of its ground state electronic densities. This causes them to start to interact, i.e., to attract each other by van der Waals force. Specifically, this is the systematic failure of local and semilocal DFT-based methods in the effort to describe sparse matter's ground state electronic structure. For large separations the long range fluctuations give the dominant contribution to vdW interaction and the local dielectric response is valid approximation. Then vdW interaction can be obtained by using the Lifshitz formula which is based on the subtraction of zero point fluctuation energies of the collective modes in the coupled and decoupled objects [30]. However, at smaller separations the microscopic structure becomes visible and dispersive (nonlocal) quantum mechanical effects start influencing local theory results. In this case a common starting point is the expression for the exact exchange-correlation energy, the so-called adiabatic connection formula (ACF) [31], that uses nonlocal response function of the interacting systems. This causes direct use of the ACF to become computationally time consuming and various approximations have to be done. The DFT approximations are based on forcing the nonlocal into local response in order to construct computationally more efficient density functionals. However, most problems arise in constructing the DFT functionals that properly connect nondispersive (long range) and dispersive (short range) correlations and the intersystem exchange. The accuracy of such connection is of crucial importance for the determination of intersystem binding energy and equilibrium separation. For noncontact separations direct use of ACF becomes computationally less demanding. Namely, ACF then can be constructed from the response functions of the separated systems which can substantially reduce the size of the supercells used in calculation. In order to check the accuracy of the above-mentioned vdW functionals in the noncontact region we provide a calculation of correlation energy by using ACF adapted to the noncontact systems, developed in Ref. [30].
We shall use these results to benchmark the long range behavior of seven different van der Waals density functional options: vdW-DF [9], vdW-DF2 [19], C09x+vdW-DF(vdW-DF2) [18], optB86b+vdW-DF [17], rev-vdW-DF2 [22], and rVV10 [28]. Hamada and Otani [32] already used both vdW-DF and vdW-DF2 with and without employing C09x (in the non-self-consistent manner) to simulate graphene interaction with metal surfaces whose properties we will reproduce here by self-consistent calculation. Even though we are benchmarking long range behavior of vdW functionals we will also make some comparisons with the previous RPA calculations [33] and to some extent with experimental measurements of Kiraly et al. [29] in the range of adsorption distances.
For our DFT van der Waals calculations we used plane wave basis set with ultrasoft pseudopotentials as implemented in QUANTUM ESPRESSO [34]. Calculations involving van der Waals functionals were done self-consistently [35]. Due to factorization of the vdW integration kernel by Román-Pérez and Soler [36] time overhead in comparison with standard DFT calculation was negligible. We used the most stable graphenesurface geometries as given in Refs. [6,7]. Silver surface was modeled by a six-layer slab. Graphene lattice (2 × 2) was matched to √ 3 × √ 3 Ag(111) lattice to avoid a large mismatch between graphene and substrate, as shown in the upper part of Fig. 1. We used the value of nominal bulk silver lattice constant a Ag = 4.09Å and defined √ 3 × √ 3 Ag(111) lattice constant as a = a Ag √ 3 √ 2/2 = 5.01Å, thus expanding nominal graphene lattice by 2%. Separation of silver layers is then c = 2.36Å with neglecting marginal relaxation of Ag(111) surface. To avoid interaction between slabs we used a supercell with the vacuum separation equivalent to 16 layers of silver atoms, as shown in the lower part of Fig. 1 (supercell parameter in the z direction is z 0 = 50Å). Calculations involving Cu(111) and Al(111) surfaces are performed using similar geometries [37]. Kinetic energy cutoff for plane waves was 475 eV for wave functions, and 5700 eV for the density. We applied a 12 × 12 × 1 k-point grid to sample the Brillouin zone. Ultrasoft pseudopotentials were generated by atomic code that is part of QUANTUM ESPRESSO [34] (from the library's suggested and tested inputs).
For the calculation of noncontact RPA correlation energies we use the prescription from Ref. [30]: Here d represents the separation between the topmost metal layer and the graphene carbon layer, as shown in Fig. 2 which represents the geometry of the interacting systems. The metal surface is treated in the jellium model where we used 50-Å-thick positive background whose edge exceeds the topmost metal layer for c/2 as sketched by orange in Fig. 2 (c is the distance between metal layers). We modeled the silver metal surface using the electronic density parameter, i.e., the Wigner-Seitz radius, of r s = 3.00 a.u. [38]. Similarly, we used r s = 2.67 a.u. for the copper surface and r s = 2.99 a.u. for the aluminum surface. The metal is translationally invariant (homogeneous) in the x-y direction; however, Kohn-Sham (KS) orbitals in the z direction are obtained self-consistently [39] by using the Wigner LDA exchange-correlation functional [40]. In (1) D m (Q,ω) represents the metallic surface excitation propagator whose detailed calculation within RPA is presented in Ref. [41]. Moreover, D g (Q,ω) represents a graphene electronic excitation propagator whose detailed calculation is presented in Refs. [42,43]. However, the calculation of D g (Q,ω) is quite delicate and has to be briefly discussed here.
In the first stage of the calculation we determine the isolated graphene KS orbitals. Electronic structure calculations for graphene are performed using the plane-wave self-consistent code QUANTUM ESPRESSO [34]. For the graphene unit cell parameter we used a = 2.46Å, which is the experimental lattice constant. For the unit cell in the z direction (separation between periodically repeated graphene layers) we take L = 5a = 12.3Å. We employ the Perdew-Zunger LDA for the exchange-correlation (xc) potential [44]. An electronic temperature of k B T ≈ 0.1 eV is used to achieve convergence of the Kohn-Sham wave functions, with all energies extrapolated to 0 K. The electronic density is calculated using a 12 × 12 × 1 Monkhorst-Pack special k-point mesh, i.e., we use 19 special points in the irreducible Brillouin zone (BZ). We use the normconserving LDA-based pseudopotentials for carbon [45]. We find the energy spectrum converges for a 680 eV plane-wave cutoff.
Once the ground state Kohn-Sham wave functions have been calculated, we may calculate the two-dimensional (2D) density-density response function within RPA. The matrix of the noninteracting density-density response function χ 0 for a quasi-2D system may be expressed as Here S is the surface area of the unit cell, η is the peak broadening, and f n (K) and ε n (K) are the Fermi-Dirac occupation and energy, respectively, of the nth band at K. The summation over K runs through 40 603 points in the BZ or 3468 in the irreducible BZ which corresponds to 201 × 201 × 1 Monkhorst-Pack special k-point mesh. The n,m summation is performed over 50 bands, which is sufficient for a proper description of the high energy σ + π plasmon in graphene. The matrix elements M have the form where Q and G are the momentum transfer and reciprocal lattice vectors, respectively, in the (x,y) plane of the graphene surface. The integration is performed over the unit-cell surface S, and the plane-wave expansion of the Kohn-Sham wave functions is where V = S × L is the unit-cell volume, G = (G ,G z ) is a three-dimensional (3D) reciprocal lattice vector, r = (ρ,z) is a 3D position vector with z normal to the surface, and the coefficients C nK are obtained by solving the Kohn-Sham equations. Integration over the perpendicular z coordinate in expression (3) is not yet performed, so the matrix elements remain z dependent. The screened response function may be obtained from Eq. (2) by solving the combined integral-matrix equation where we have suppressed the Q and ω dependence to simplify the expression. Note that the z 1 and z 2 integrations from − L 2 to L 2 ensure that interactions between repeated images are removed. In this way we obtain a correct description for an isolated graphene monolayer. Here is the matrix of the bare Coulomb potential. Equation (5) can be solved by applying a Fourier transform in the z direction when it becomes a pure matrix equation [42]. A graphene electronic excitation propagator can be obtained from the screened response function (5) as For realistic crystal structures D g (Q,ω) is an anisotropic function of Q. This means that the intensity and frequency of the electronic modes depend on the propagation direction. For example, the π plasmon dispersion in graphene and carbon nanotubes splits if Q is in the → M direction but does not split if it is in the → K direction [43]. However, we approximated D g by an isotropic function by averaging over the high symmetry − M and → K directions Moreover, in (5) we put reciprocal lattice vectors G = 0 which means that we excluded crystal local field effects in the parallel direction. Factor e −2Qd in (1) defines the cutoff wave vector Q C in Q integration. Therefore, for the noncontact separations d 5Å investigated here, the cutoff wave vector can be estimated as Q C 1/d = 0.2Å −1 , so in the calculation we have chosen Q C = 0.3Å −1 . Because of smooth Q-dependent functions, Q integration is performed using the Simpson's method over 51 points. The ω integration is also performed using a Simpson's method with 2001 points up to 50 eV. In this way we ensure the inclusion of the entire spectral weight of the high energy σ + π plasmon. For the noncontact separation the main contribution to van der Waals interaction comes from interaction between collective excitations in the long wavelength limit, such as the surface plasmon in silver and the π plasmon in graphene. Moreover, because the silver surface plasmon ω S = 3 2r 3 s ≈ 6.4 eV is close to resonance with linearly dispersed graphene π plasmon ω π ≈ 4-7 eV [43], the vdW attraction possesses a long range tail. Similar behavior also occurs for the Cu(111) and Al(111) surfaces.
From the viewpoint of vdW functionals, the metallic long range correlations (surface plasmon) can be simply implemented in ACF. First metallic surface boundary (effective image plane) should be carefully determined and then the local Drude dielectric function (ω) can directly enter ACF. The inclusion of graphene π plasmon is not so straightforward. Namely, π plasmon is a Landau damped interband plasmon and there is no simple Drude-like parametrization of (ω). Therefore constructing the vdW functional with graphene as one of the interacting components is not feasible and the direct use of ACF is recommendable.
Binding curves for different vdW functionals are shown in Fig. 3. Binding energies E b are calculated as where E vdW is the energy of the graphene-Ag(111) system, E Ag(111) is the energy of Ag(111) slab, and E graphene is the energy of a graphene layer. All the energies are per the same unit cell from Fig. 1. To get binding energy per carbon atom their difference is divided by 8, as there are eight carbon atoms in the unit cell. The black line in Fig. 3 shows the correlation interaction energy obtained by using ACF-RPA, i.e., directly by using Eq. (1). In exploring Fig. 3 we will discuss three different regions: long range d 6Å, intermediate 4.5 d 6Å, and adsorption region d 4.5Å. In the long range region exchange energy is so small that only the correlation part is important, which is apparent from the  fact that curves with the same correlation part (vdW-DF or vdW-DF2) join in the same curve. Our noncontact RPA curve coincides with the supercell RPA result [33] at 6Å showing equivalence of the two approaches in the noncontact region. The vdW-DF2 curves are almost on top of the model RPA correlation curve, suggesting good long range behavior of the vdW-DF2 functional. On the other hand, vdW-DF as it is and combined with other exchange functionals (including the rVV10 functional) overestimates the ACF-RPA result for about 40% for almost all noncontact separations. In the intermediate region we start to see the consequences of using different exchange functionals. As our ACF-RPA result does not include exchange we cannot compare it to the vdW-DF results in the whole intermediate region where exchange is contributing. However, in the higher intermediate region, d 5Å, the exchange is already negligible and noncontact RPA results obtained here agree remarkably well with the RPA results from Ref. [33] that include exchange energy by evaluating exact exchange (Hartree-Fock) energy on the PBE orbitals. On this basis we can use any similarity with the ACF-RPA result to estimate how does particular exchange decay and which correlation functional is the most appropriate. We can see that the combined vdW-DF2(C09x) agrees the best with the ACF-RPA result in the higher intermediate region, followed by rev-vdW-DF2 and then by vdW-DF2.
Even though presented noncontact RPA calculations for Ag(111) agree well with the RPA calculations of Ref. [33], we investigated what are the effects of replacing metal surface with jellium slab and whether our analysis is valid for other metal surfaces. As a test of the importance of the atomic structure of metal versus the flat jellium model we calculated the change in binding energy per carbon atom given by the vdW functional when graphene is shifted to the right by 1/4 of the cell shown in Fig. 1. Differences are smaller than 10 −3 meV for distances larger than d 4.5Å, and smaller than 0.1 meV for distances larger than et al. [33] are shown for comparison. It can be noted that for this metallic surface both RPA calculations agree in the region where exchange is negligible (d 5Å). Again vdW-DF2 correlation closely follows the RPA results in the long range region.
Let us focus again on the Ag(111) surface and look at the behavior of the vdW functionals in the adsorption region where noncontact RPA is of little significance. To understand the results in the adsorption region we listed equilibrium distances and adsorption energies for all used vdW functionals in Table I, together with the RPA result [33] and the experimental equilibrium distance [29]. Experimental distance is measured by STM topography and because of differences in tunneling conductance between the two materials the authors of Ref. [29] argue that it could be too low. In agreement with expectations, vdW-DF and vdW-DF2 in their original form give too large adsorption distances. Other combinations of exchange functionals with the vdW-DF(2) correlation functional give equilibrium distances that agree well with the previous RPA result from Ref. [33]. Even though this RPA calculation gives adsorption energy that agrees best with the vdW-DF(C09x) result, those calculations are hard to converge and there is the possibility that this energy is slightly lower, in the range of energies obtained by using the rev-vdW-DF2 functional. We base this on the fact that for Ni(111) and Cu(111) surfaces binding energies given by the rev-vdW-DF2 functional are in very good agreement (up to 5 meV) with the same RPA calculations from Ref. [33] as shown in Ref. [22]. Another argument in this direction is that the experiments of Kiraly et al. [29] suggest weaker binding of graphene on Ag(111) than on Cu(111), while RPA results [33] predict 10 meV per carbon atom higher binding energy for Ag(111). Comparing our result with the non-self-consistent results of some of these functionals for our system, as given by Hamada and Otani [32], we note that the results almost do not differ (differences in adsorption energy are up to 4 meV), confirming that self-consistency is of little importance. To further investigate this system we used the vdW-DF2(C09x) combination of functionals to calculate the band structure of graphene on Ag(111) at the equilibrium distance d = 3.27Å as shown in Fig. 5. Color function C nk for every band n and k point in Fig. 5 was calculated as where z M is the plane in the middle between the graphene and the topmost silver layer and z 0 is supercell parameter in the z direction. In other words it shows probability that the electron with given energy and k belongs to the graphene layer. As one can see from Fig. 5, graphene's Dirac cone (in the K point) is preserved, as expected because of weak interaction. However, we can notice slight electronic doping of graphene, i.e., the Dirac cone moves −0.46 eV below the Fermi level. This result slightly overestimates previously reported values of −0.32 eV by LDA [7] and −0.40 eV by vdW-DF [11]. Moreover, the experimentally observed slight blue-shifted G band in the Raman spectrum [29] also indicates light electronic doping. The strength of the long range vdW tail is not of such physical relevance as the vdW energy in the shorter range, which defines the intersystem binding energy or equilibrium distance. However, the behavior of particular functionals in the noncontact region reflects in their behavior in the contact region, i.e., in its competition with the repulsive exchange energy. In this paper we showed that vdW-DF2(c09x) or rev-vdW-DF2 functionals apart from giving a good binding energy or equilibrium distance [18,22], also behave satisfactorily well in the noncontact region. They agree within 5 meV with the ACF-RPA result up to very close separations of d = 4.5Å (where the exchange turns on). Good agreement between our model noncontact RPA calculations and previous supercell RPA results [33] indicates that the details of the surface structure might be irrelevant for proper description of vdW force in the noncontact region. In conclusion, due to its precision in both long and short range interactions and its fast convergence the second van der Waals functional together with the appropriate exchange functional [vdW-DF2(c09x) or rev-vdW-DF2] looks like the most adequate functional in describing the similar systems in which dispersive forces are significant.
We thank M.Šunjić for useful suggestions. V.D. is grateful to the Donostia International Physics Center (DIPC) and Pedro M. Echenique for their hospitality during various stages of this research. Computational resources were provided by the DIPC computing center.