Sanliang Linga, Matthew B. Watkinsa and Alexander L. Shlugerab
aDepartment of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK. E-mail: S.Ling@ucl.ac.uk
bLondon Centre for Nanotechnology, University College London, 17-19 Gordon Street, London WC1H 0AH, UK
First published on 26th September 2013
We evaluate the performance of different van der Waals (vdW) corrected density functional theory (DFT) methods in predicting the structure of perfect interfaces between the LiF(001), MgO(001), NiO(001) films on the Ag(001) surface and the resulting work function shift of Ag(001). The results demonstrate that including the van der Waals interaction is important for obtaining accurate interface structures and the metal work function shift. The work function shift results from a subtle interplay of several effects strongly affected by even small changes in the interface geometry. This makes the accuracy of theoretical methods insufficient for predicting the shift values better than within 0.2 eV. Most of the existing van der Waals corrected functionals are not particularly suited for studying metal/insulator interfaces. The lack of accurate experimental data on the interface geometries and surface rumpling of insulators hampers the calibration of existing and novel density functionals.
The metal/insulator interface distance is one of the most important factors that affect the local Δϕ. It is determined by the balance between the short-range (e.g. ionic or covalent bonding) and long-range (e.g. image and van der Waals) interactions.16 As pointed out by Jupille et al., in the limit of a noble metal interacting with a wide band gap insulator, the adhesion energy is dominated by long-range interactions.16 The theoretical description of such interfaces based on conventional DFT functionals may therefore fail to predict the correct interface geometry and adhesion energy as they do not capture long-range instantaneous electron correlation effects responsible for the van der Waals interactions.17,18 For example, Levchenko et al.19 recently used a hierarchy of DFT methods including PBE, HSE06, and the state-of-the-art exact exchange plus correlation in the random phase approximation (EX-cRPA/cRPA+) to calculate the adhesion energies of AuN clusters (N > 1) supported on MgO(100) surfaces and found that the van der Waals interaction, missing in PBE and HSE06 methods but accounted for by the EX-cRPA/cRPA+ method, makes a significant contribution to the adhesion energy of AuN clusters. It is therefore reasonable to assume that the van der Waals interaction may also play an important role in describing the properties of insulator thin films supported on metal substrates.
However, using accurate methods, such as EX-cRPA/cRPA+, to account for the van der Waals interaction at metal/insulator interfaces is extremely expensive and calculating forces with such methods is still a challenge. It is desirable to have a method that is relatively cheap and yet can correctly describe the van der Waals and other interactions at metal/insulator interfaces in a self-consistent manner. The recently developed van der Waals corrected DFT methods20–22 in principle provide such a possibility, but their accuracy has not yet been tested on real metal/insulator systems. They must describe the subtle interplay of three main effects that determine metal Δϕ due to adsorption of an insulator thin film:23,24 (i) Coulomb repulsion of metal electron density by insulator anions, which reduces the intrinsic metal surface dipole originating from electron density spilling out of the metal surface, and this decreases the effective work function of the interface system; (ii) interfacial charge transfer across the interface, which creates an interfacial dipole, either pointing towards the metal (electron transfer from insulator to metal) and decreasing the effective work function of the interface system, or pointing towards the insulator (electron transfer from metal to insulator) and increasing the effective work function of the interface system; (iii) insulator rumpling, which creates a dipole due to the finite separation between anionic and cationic atomic layers in an insulator, and this can either decrease the effective work function, if the anion is sitting closer to the metal substrate, or increase the effective work function, if the cation is sitting closer to the metal substrate. These effects are schematically presented in Fig. 1. We note that the effects of dipoles due to interfacial charge transfer and insulator rumpling on metal work function can be discussed in terms of the parallel plate capacitor model,25 which has been applied successfully in theoretical studies of self-assembled organic monolayers on metal substrates.26,27
Fig. 1 Various mechanisms that affect the shift of metal work function in contact with an insulator thin film. (a) Electron density spilling out of the metal surface naturally into the vacuum, creating a surface dipole Dm pointing outwards from the metal surface; (b) upon insulator deposition, due to Coulomb repulsion of the metal electron density by insulator anions, the metal surface dipole is decreased, resulting in a smaller surface dipole (denoted as Dcomp), and this decreases the effective work function; (c) electrons transfer from insulator to metal, resulting in a “positive” insulator and a “negative” metal, which creates an interfacial dipole (denoted as DCT) pointing towards the metal substrate and decreases the effective work function; (d) electrons transfer from metal to insulator, resulting in a “negative” insulator and a “positive” metal, which creates an interfacial dipole (denoted as DCT) pointing towards the insulator and increases the effective work function; (e) a rumpled insulator layer at the interface, i.e. cations sitting closer to the metal substrate, which creates an intrinsic dipole inside the insulator (denoted as Drump) pointing towards the insulator and increases the effective work function; (f) a rumpled insulator layer at the interface, i.e. anions sitting closer to the metal substrate, which creates an intrinsic dipole inside the insulator (denoted as Drump) pointing towards the metal substrate and decreases the effective work function. |
In this paper, we evaluate different van der Waals corrected DFT methods implemented in the VASP code using three model systems: LiF, MgO and NiO layers on the Ag(001) substrate. We optimize the interface geometries to compare the resulting interface structures and then discuss the effects on the shift of metal work function due to variations in the interface distance and insulator rumpling, as predicted by different methods. Thin MgO films on the Ag(001) surface have been studied extensively both experimentally14,15,28–38 and theoretically10,23,39–43 and are known to be rough on the nanoscale, forming some islands.30,37,38 Therefore Δϕ depends strongly on the film quality and morphology (see Table 1). However, not much is known about Δϕ at LiF/Ag and NiO/Ag interfaces. The effect of the film roughness has been discussed in our previous work.43 Here we demonstrate how different geometric factors at the interface of a perfectly matched flat insulator thin film affect the work function of the metallic support.
In addition to the conventional PBE method, the following van der Waals corrected DFT methods are considered: PBE + D2,44 PBE + TS,45 PBE + TS + SCS,46 optPBE,47 and vdW-DF2.48 For lattice parameters of isolated bulk materials, we also considered a meta-GGA, M06L.49 Depending on how the van der Waals interactions are included, these methods can be divided into three groups, i.e. a posteriori methods, including PBE + D2, PBE + TS, which add on top of the conventional PBE energy an additional van der Waals interaction term (represented by an additive pairwise C6/r6 expression), the non-local correlation functionals, including optPBE and vdW-DF2, which add non-local corrections to local or semi-local correlation functionals, and calculate the electronic structure in a self-consistent way, and PBE + TS + SCS where many-body screening effects are included to some degree.
For the first group of methods, different approaches exist to obtain C6 dispersion coefficients. For example, PBE + D2 calculates C6 coefficients as a function of the ionization potential and static polarizability of isolated atoms. While this approach exhibits good performance on the S22 database,50 there are two problems that might affect its performance for the metal/insulator interface systems: the lack of dependence of the C6 coefficients on different chemical environments and using the isolated atom as a reference can be inadequate when this atom forms a strong ionic bond with another atom. These problems are somewhat alleviated if one uses the PBE + TS and PBE + TS + SCS methods developed by Tkatchenko and co-workers. In the PBE + TS method,45 the atomic C6 coefficients are dependent upon their chemical environments by employing a Hirshfeld partitioning of the electron density to determine the effective volume for an atom inside a molecule/solid. This partition is then used as a scaling factor to calculate the effective atomic polarizability and the effective atomic C6 coefficients. The PBE + TS + SCS method goes a step further by solving the self-consistent screening equation of the frequency-dependent polarizability, which is then used to calculate the effective C6 coefficients. Such an approach enables one to take into account the effect of the dynamic electric field produced by the surrounding polarizable atoms on the effective C6 coefficient of a particular atom.
For the second group of methods, the conventional DFT functionals are supplemented with an extra non-local correlation energy term, which is represented by a double space integral of the electron density, with the integral kernel expressed by a pairwise formula.22 Different non-local correlation functionals have been proposed, and also different exchange functionals were employed. Many of these functionals have been reported to show qualitatively different results for molecular complexes47 and adsorption of organic molecules on surfaces of inorganic solids,51,52 compared with conventional GGA functionals. In this study we assess the performance of two typical non-local correlation functionals, optPBE and vdW-DF2.
While the first and third groups of methods do not modify the electronic structure of the interface other than indirectly through changing the interface geometry, the second group of methods solve the Kohn–Sham equation self-consistently and thus modify both the electronic structure as well as the geometry of an interface. Note that the performance of the second group of methods on the electronic structure of solid-state materials has not been fully tested and validated against high-level theoretical or experimental results, and it is not clear whether they provide correct band offsets at metal/insulator interfaces. Since performing such a benchmarking is beyond the scope of this paper, we use these methods only to obtain interface geometries.
Describing electron transfer between a metal and an insulator requires accurately calculating the band offset at the interface. Conventional GGA functionals, e.g. PBE and PW91, systematically underestimate band gaps of insulators.53 Therefore, we have performed electronic structure calculations for all the three metal/insulator interfaces with a screened hybrid functional (HSE06),54 which predicts reasonable insulator band gaps and is expected to better describe band offsets at the interfaces as well as the interfacial charge transfer. However, due to the computational demands of such hybrid functional calculations on mixed metal/insulator systems, the HSE06 calculations are single point electronic structure calculations at interface geometries optimized at PBE or PBE + D2 level.
We assume that a thin film of insulator is grown on top of the metal substrate and follows pseudomorphic growth. In other words, due to the relatively small (less than 5%) lattice mismatch between the two materials, the insulator is forced to adopt the lattice parameter of the metal substrate. Therefore, for all the three interfaces, we used Ag lattice parameters (as determined using the PBE functional) for constructing unit cells of interface slabs along X and Y directions. The most stable interface configuration was considered, with the anion located directly above the Ag atom and the cation above hollow sites. The smallest asymmetric unit cell included four layers of Ag and three layers of insulator (with a non-polar (001) surface). For NiO, we considered an antiferromagnetic (AFM2) magnetic structure of the Ni atoms. Therefore, each Ag layer has one (for LiF/Ag and MgO/Ag) or two (for NiO/Ag) Ag atoms, and each insulator layer has two (for LiF and MgO) or four (for NiO) atoms.
All calculations are performed using the Vienna Ab initio Simulation Package (VASP 5.3.3).55 A plane wave basis set with a 400 eV energy cutoff was employed to represent the wavefunctions, and the projector augmented wave (PAW) method56 was employed. For LiF and MgO a spin-unpolarized DFT method and for NiO a spin-polarized DFT + U method (with or without van der Waals corrections) were used. Following the previous calculations,24 a U value of 4.0 eV was used for d-electrons of Ni atoms in the interface layer of NiO to take into account the metal screening in the insulator, and a U value of 5.3 eV was used for d-electrons of Ni atoms in the surface and middle layers. Geometry optimizations of interface slabs were considered converged if the maximum force on relaxed atoms falls below 0.01 eV Å−1. Monkhorst–Pack k-point grids57 of (9 × 9 × 9) and (12 × 12 × 1) were used for isolated bulk materials and heterogeneous interface slabs, respectively. A vacuum layer exceeding 40 Å was used throughout the current study to reduce the artificial interaction between periodic images of interface slabs and to converge the electrostatic potential in the vacuum. Dipole corrections have been applied throughout the current study to eliminate interactions between total dipole moments of repeated interface slabs along the Z direction. All charge analyses were performed using a Bader method58 based on the tools developed by Henkelman and co-workers.59
PBE | PBE + D2 | PBE + TS | PBE + TS + SCS | optPBE | vdW-DF2 | M06L | Expt. | |
---|---|---|---|---|---|---|---|---|
Ag | 4.15 | 4.15 | 4.07 | 4.12 | 4.18 | 4.26 | 4.19 | 4.069 |
LiF | 4.07 | 4.07 | 3.67 | 3.98 | 4.07 | 4.07 | 3.99 | 4.010 |
MgO | 4.25 | 4.19 | 4.19 | 4.22 | 4.25 | 4.27 | 4.19 | 4.207 |
NiO | 4.19 | 4.19 | 4.16 | 4.19 | 4.19 | 4.31 | 4.23 | 4.195 |
One can see that all methods (except PBE + TS + SCS) overestimate the lattice constant of Ag. Notably vdW-DF2 overestimates by 4.7%, compared with 2.0% of PBE. This behaviour has also been observed in several previous studies.62,63 A similar overestimation of the Ag lattice parameter was also observed for optPBE, which can be regarded as belonging to the same family as the vdW-DF2 method, as both methods employ a non-local correlation functional to account for the van der Waals interaction and both methods were validated against the S22 data set consisting of molecular complexes. It is therefore important that metallic systems are included in fitting balanced vdW functionals for these types of metal/insulator interface systems.
We note that both schemes developed by Tkatchenko et al. (PBE + TS and PBE + TS + SCS) predict smaller lattice parameters for the four materials compared with PBE. This effect is most significant in the case of LiF, where the PBE + TS method gives a lattice parameter of 3.67 Å, almost 10% smaller than PBE. Similar underestimation of the lattice parameters by PBE + TS has also been reported for NaCl and KI,64 indicating that the PBE + TS method strongly overestimates the strength of van der Waals interactions in alkali halides. This is due to the fact that the Hirshfeld partitioning employed by the PBE + TS method fails for ionic solids with strong charge transfer character,65 such as alkali halides. These strongly ionic systems are too far from the atomic reference for which the C6 coefficients are calculated for the scaling based on atomic volumes to provide a good approximation. This problem is partially solved in the PBE + TS + SCS method, in which the underestimation of lattice parameters was reduced due to damping of the ionic polarizabilities by solving the self-consistent screening equation of the frequency-dependent polarizability,46e.g. the lattice parameter of LiF with PBE + TS + SCS is underestimated only by 2.2% compared with PBE.
In the current study, the interface distance is defined between the Ag atoms and anions at the interface (denoted as dint, see Fig. 2b). The definition of the insulator rumpling is more complicated, as we are dealing with a three-layer insulator thin film, and rumpling may exist in all three layers. As the first step, we define the insulator rumpling of the interface layer (denoted as drump_int, see Fig. 2b) as the difference in the Z coordinates of cations and anions at the metal/insulator interface: if cations are closer to the metal substrate, this is defined as positive rumpling, and if anions are closer to the metal substrate, negative rumpling. Note that rumpling also exists in the top two layers of the insulator (denoted as drump_mid and drump_surf, respectively, see Fig. 2b), due to propagation of insulator rumpling at the interface layer and intrinsic rumpling due to different polarizabilities of anions and cations at the surface layer.67
Fig. 2 Schematic illustrating the definitions of the interface distance and insulator rumpling in (a) free-standing and (b) metal-supported insulator thin film used in this paper. |
Rumpling in each insulator layer produces different dipole moments along the direction which is given by the relative positions of cations and anions, and the magnitude is determined by the charges on ions as well as the separation of the cationic and anionic atomic layers (difference in Z coordinates). As the first approximation, we consider three dipole moments due to rumpling in each insulator layer and neglect the dipole–dipole interactions, as well as different separations between each insulator layer given by different methods. Then the final dipole moment intrinsic to the whole system is a vector sum of the three dipole moments given by each insulator layer. If we neglect the difference in the charges on anions and cations in each insulator layer, the dipole moment in each insulator layer can be directly correlated to the insulator rumpling in each layer. In this way, the total dipole moment intrinsic to the insulator slab can be described by only one parameter, i.e. the total insulator rumpling (denoted as drump_tot, see Table 3), which is a vector sum of the insulator rumpling present in each insulator layer. Below we discuss these geometric components of the interface in more detail.
PBE | PBE + D2 | PBE + TS | |||||||
---|---|---|---|---|---|---|---|---|---|
dint | drump_int | drump_tot | dint | drump_int | drump_tot | dint | drump_int | drump_tot | |
LiF | 2.98 | −0.035 | 0.025 | 2.56 | 0.045 | 0.093 | 2.59 | −0.015 | 0.073 |
MgO | 2.70 | 0.012 | 0.045 | 2.51 | 0.034 | 0.073 | 2.52 | 0.027 | 0.073 |
NiO | 2.53 | 0.053 | 0.073 | 2.52 | 0.059 | 0.069 | 2.41 | 0.071 | 0.098 |
PBE + TS + SCS | optPBE | vdW-DF2 | |||||||
---|---|---|---|---|---|---|---|---|---|
dint | drump_int | drump_tot | dint | drump_int | drump_tot | dint | drump_int | drump_tot | |
LiF | — | — | — | 2.94 | −0.026 | 0.029 | 2.92 | −0.039 | 0.003 |
MgO | 2.37 | 0.061 | 0.108 | 2.69 | 0.016 | 0.044 | 3.04 | −0.006 | 0.022 |
NiO | 2.31 | 0.097 | 0.133 | 2.50 | 0.059 | 0.069 | 2.73 | 0.031 | 0.034 |
One can see that the general trend of the interface distance follows LiF/Ag > MgO/Ag > NiO/Ag, which reflects different strengths of the interaction between Ag and the three different insulators. As a simple rule of thumb, this ordering follows the ordering of insulator band gaps: LiF (14.2 eV) > MgO (7.8 eV) > NiO (4.3 eV). More strictly speaking, the adhesion between Ag and an insulator is determined by the relative position of the Ag Fermi level with respect to the top of the valence band and the bottom of the conduction band of the insulator.
As shown in Table 3, as a result of the van der Waals corrections, the interface distances for all three model systems are further reduced, and for LiF/Ag, this reduction can be as big as 0.4 Å using the PBE + D2 and PBE + TS methods. We note that this value might be overestimated due to the fact that the dispersion coefficients were predetermined and kept fixed during the calculations within the PBE + D2 method. PBE + TS also tends to overestimate the strength of van der Waals interactions in alkali halides. On the other hand, optPBE gives a similar interface distance to PBE. The vdW-DF2 method predicts even larger interface distances of 3.04 Å for MgO/Ag and 2.73 Å for NiO/Ag, which are by 0.3 Å and 0.2 Å larger than those of PBE, respectively, and appear to be too large compared to available experimental data (see Table 1). This indicates that the performance of these functionals needs to be further tested and validated against data sets which include high-level experimental or theoretical data on metal/insulator interfaces.
Comparing the reductions in interface distances due to the inclusion of various van der Waals corrections, LiF/Ag (0.4 Å) shows the biggest reduction, followed by MgO/Ag (0.2 Å with PBE + D2 and PBE + TS methods) and NiO/Ag (0.1 Å with PBE + TS). This reflects the interplay of the van der Waals and covalent forces in the three interface systems. For the weakly interacting LiF/Ag interface with essentially no interfacial charge transfer, the van der Waals interaction dominates the binding and thus larger reduction in the interface distance is expected if the van der Waals correction is included. For MgO/Ag and NiO/Ag, where there is some interfacial charge transfer, covalent forces start to compete with the van der Waals interaction, or even dominate. Therefore the reduction in interface distance is smaller (for MgO/Ag) and can even be neglected (for NiO/Ag).
Fig. 3 Electron density differences between (a) LiF/Ag, (b) MgO/Ag, and (c) NiO/Ag interface systems (with PBE geometries) and summations of electron densities of isolated Ag and insulator (LiF/MgO/NiO) slabs. Wireframes and solid hypersurfaces indicate electron density accumulation and reduction in the interface system, respectively. |
The change of sign of the insulator rumpling at the interface layer of LiF/Ag from PBE geometry to PBE + D2 geometry is probably an artefact owing to the fact that the PBE + D2 method used the atomic polarizability of a lithium atom (which is very polarizable) to calculate the dispersion coefficient, but actually the lithium cation in LiF is less polarizable than the fluorine anion. As there is almost no interfacial charge transfer, and the electron density accumulation at the metal hollow site due to the electrostatic compression effect is small in LiF/Ag with PBE geometry (see Fig. 3), insulator rumpling at the interface should be very similar to that of the isolated LiF slab, and thus the anion should sit closer to the metal (negative rumpling).
For MgO/Ag, our calculations using the PBE method show that no matter the thicknesses of the MgO slab (up to seven layers of MgO supported on Ag were considered), the total rumpling of the whole MgO slab is almost the same, i.e. around 0.04–0.06 Å, which explains why the shift of Ag work function due to MgO thin film saturates at the MgO thickness of three layers. We also found that with the PBE geometry, the total rumpling of the whole MgO slab is mainly determined by the surface layer (i.e. rumplings in the middle layers and interface layer almost cancelled each other out), which is due to different polarizabilities of the anions and cations at the surface.
dint | CT | Δϕ | |
---|---|---|---|
LiF/Ag | |||
PBE | 2.98 | 0.0018 | 0.65 |
PBE + D2 | 2.56 | −0.0010 | 0.87 |
PBE + TS | 2.59 | 0.0021 | 0.82 |
HSE06 (PBE) | 2.98 | 0.0065 | 0.77 |
HSE06 (PBE + D2) | 2.56 | −0.0029 | 0.91 |
MgO/Ag | |||
PBE | 2.70 | −0.0620 | 1.36 |
PBE + D2 | 2.51 | −0.0437 | 1.33 |
PBE + TS | 2.52 | −0.0806 | 1.18 |
PBE + TS + SCS | 2.37 | −0.0486 | 0.98 |
HSE06 (PBE) | 2.70 | −0.0405 | 1.39 |
HSE06 (PBE + D2) | 2.51 | −0.0507 | 1.45 |
NiO/Ag | |||
PBE | 2.53 | 0.0658 | 0.47 |
PBE + D2 | 2.52 | 0.0702 | 0.52 |
PBE + TS | 2.41 | 0.0743 | 0.42 |
PBE + TS + SCS | 2.31 | 0.0485 | 0.31 |
HSE06 (PBE) | 2.53 | 0.0561 | 0.60 |
HSE06 (PBE + D2) | 2.52 | 0.0603 | 0.65 |
dint | CT | Δϕ | |
---|---|---|---|
LiF/Ag | |||
PBE (flat) | 2.98 | 0.0011 | 0.92 |
PBE + D2 (flat) | 2.56 | −0.0065 | 2.04 |
MgO/Ag | |||
PBE (flat) | 2.70 | −0.0475 | 1.58 |
PBE + D2 (flat) | 2.51 | −0.0549 | 2.08 |
We also note that PBE predicts bigger interfacial charge transfer in the PBE geometry, while HSE06 predicts bigger interfacial charge transfer in the PBE + D2 geometry. The difference between PBE and HSE06 results again shows the importance of using a hybrid functional to correctly describe the band offset at the interface, which determines the interfacial charge transfer and thus the calculated shift of the metal work function.
The fact that relatively minor changes in the interface structure, as characterized by interface distance and insulator rumpling, can induce such a big difference (see Tables 4 and 5) in the calculated Δϕ values demonstrates that these interfaces are very delicate systems for testing DFT methods. Accurate experimental data are vital for developing both theoretical methods and technological applications, but Δϕ depends strongly on the film quality and interface structure (see Table 1). With such a degree of uncertainty in both the insulator thin film quality and experimental measurements, any rigorous comparisons between theoretical calculations and relevant experiments should be considered as qualitative rather than quantitative.
Experimentally, Δϕ can be measured in different ways: three of the most popular methods are Kelvin probe force microscopy, scanning tunnelling microscopy (STM) in dI/dV mode, and field emission resonance (FER).14 These methods average Δϕ over an area of several hundred square nanometres, or even more. For example, the Kelvin probe measurements of a MgO/Ag interface average over a surface area with a radius of 15–30 nm.14 This is larger than the MgO island size, and thus the bare metal substrate was most likely included in the measured work function shift. In this paper we ignored the character of film growth and focused on the idealized model of full coverage. However, for metal/insulator interfaces with large lattice mismatch, e.g. CsBr/Cu(100)8 and NaCl/Cu(110),71 the insulator films are likely to be rough and may include dislocations, corrugations and grain boundaries, and thus the interface distance and the insulator rumpling can change significantly from one small area to another affecting the shift.
To summarize, the results demonstrate that the accuracy of both experimental and theoretical methods is only sufficient for predicting Δϕ values within 0.2 eV at best. Most of the existing van der Waals corrected functionals are not particularly suited for studying metal/insulator interfaces (see ref. 21 and 77 and references therein for similar results on layered solids and organic/inorganic interfaces). The lack of accurate experimental data on the interface geometries and even surface rumpling of insulators hampers the calibration of existing and novel density functionals.
This journal is © the Owner Societies 2013 |