First-principles study of the adsorption of atomic H on Ni (111), (100) and (110)

Surface Science 459 (2000) 287–302 www.elsevier.nl/locate/susc

First-principles study of the adsorption of atomic H on Ni (111), (100) and (110) G. Kresse *, J. Hafner Institut fu¨r Materialphysik and Center for Computational Material Science, Universita¨t Wien, Sensengasse 8/12, A-1090 Vienna, Austria Received 29 December 1999; accepted for publication 21 March 2000

Abstract The adsorption of atomic H on Ni (111), (100) and (110) surfaces is studied using spin-polarized gradientcorrected density-functional theory. Full H-coverage and 1/4 H-coverage are considered. Trends in the adsorption energy, repulsion between adsorbed hydrogen atoms, and diffusion barriers are discussed in detail. In addition, calculations are presented for the experimentally observed Ni(111) (2×2)-2H superstructure and for the Ni(110) (1×2)-3H pairing-row reconstruction. The results are compared with experiments, showing a good agreement with the structural models derived from experimental data. We use a simulated annealing approach to show the global stability of the pairing-row reconstruction. © 2000 Elsevier Science B.V. All rights reserved. Keywords: Chemisorption; Density functional calculations; Hydrogen atom; Low index single crystal surfaces; Nickel; Surface relaxation and reconstruction

1. Introduction The adsorption of H on low-index Ni surfaces has been investigated experimentally and theoretically for more than 30 years. Consequently, a wealth of experimental data has accumulated, and a clear picture of the features induced by the hydrogen adsorption has emerged. It is interesting, however, that only a few accurate density-functional calculations have been reported to date. The major motivation of the present work is to fill this gap. We start with a brief summary of the experimental data: On all low-index surfaces, H adsorbs dissociatively. The (111) and the (100) surfaces remain unreconstructed under H exposure. On * Corresponding author. Fax: +43-1-4277-9514/13. E-mail address: [email protected] (G. Kresse)

both surfaces, hydrogen is believed to adsorb in the hollow site, and the saturation coverage is around h=1 monolayer (ML). At a coverage of h=0.5 ML, the (111) surface exhibits a welldefined (2×2)-2H superstructure that is stable up to a temperature of 250 K. On the basis of early low-energy electron diffraction (LEED) measurements, it was suggested that the hydrogen atoms adsorb in the fcc and hcp hollows, forming a graphitic over-layer [1,2]. This conclusion was later confirmed by accurate tensor-LEED ( T-LEED) measurements [3], which also indicate that the Ni substrate exhibits a considerable buckling. Interestingly, on the (100) surface, either no ordering [4–8] or only a rather poorly ordered p(2×2) superstructure was found [9]. This indicates that the interaction between adsorbed H atoms differs at least quantitatively on both surfaces, and ab-initio calculations might help to explain the observed differences.

0039-6028/00/$ - see front matter © 2000 Elsevier Science B.V. All rights reserved. PII: S0 0 39 - 6 0 28 ( 00 ) 0 04 5 7 -X

288

G. Kresse, J. Hafner / Surface Science 459 (2000) 287–302

In contrast to the other two low-index surfaces, the more open (110) surface shows a rather complex behavior with varying reconstructions upon H exposure. At full coverage and low temperatures, a fairly well-ordered (2×1)-2H superstructure is observed. [10–12] When the coverage is increased to h=1.5, a first-order phase transition to a (1×2)-3H pairing-row reconstruction is found [12–17]. Irrespective of coverage, at temperatures above 200 K, a one-dimensionally ordered ‘streak’ phase is stable [15]. Little is known about the H adsorption energies on the (110) surface since the reconstruction to the streak phase makes an interpretation of temperature-programmed desorption experiments difficult. Therefore, first-principles calculations are the only reliable source for adsorption energies. For the (110) surface, we restrict our present study to the (2×1)-2H superstructure and the (1×2)-3H pairing-row reconstruction because rather precise models are available for both phases [12], whereas only little is know about the ‘streak’ phase. As we have already stressed, only a few theoretical first-principles calculations have been published. For the (111) surface, density-functional studies for the high-coverage case have been reported by Paul and Sautet [18], and CI calculations have been performed by Yang and Whitten [19]. For the (100) surface, calculations for the high-coverage case have been reported too [20]. However, many questions have yet to be addressed. A thorough comparison of the adsorption energies and adsorption geometries on the three low-index surfaces is missing, for instance. Another question concerns the repulsion between H adatoms adsorbed on the surface. Repulsion could be either purely electrostatic, or substrate-mediated. In the later case, the adsorbed species modifies the electronic states in a way that reduces the adsorption energy of neighboring H atoms. We will investigate which effect is more important for H atoms adsorbed on the Ni surface. Finally, we address the question of whether the magnetic properties of the surface have any influence on the adsorption energy by performing calculations for artificial non-magnetic Ni. The paper is organized as follows. First, we briefly discuss the computational method used in

this work. As a point of reference, we then present results for bulk Ni, the H dimer (Section 3.1), 2 and the clean low-index surfaces (Section 3.2). In Section 3.3, we discuss our results for H adsorption in high-symmetry sites on all three low-index surfaces at 1/4 and full coverage. The results of fully spin-polarized (magnetic) calculations are compared to non-spin-polarized calculations. Then, we address the repulsion between adsorbed H adatoms and the work function changes (Section 3.3.3). The (111) (2×2)-2H superstructure is studied in Section 3.4, and detailed results for the (110) pairing-row reconstruction are reported in Section 3.5. We finish with Discussions and Conclusions in Section 4.

2. Computational methods The first-principles calculations are based on density-functional theory (see, for example, Refs. [21,22]) employing a plane wave basis set [23,24]. To solve the Kohn Sham equations, we use the Vienna ab-initio simulation package ( VASP) [25– 27], which performs an iterative solution of Kohn Sham equations via an unconstrained band-byband minimization of the norm of the residual vector of each eigenstate and an optimized charge density mixing. The electron–ion interaction is described by the projector augmented wave (PAW ) method, as proposed by Blo¨chl [28]. The PAW method is a frozen core all-electron method and, in many aspects, resembles the ultrasoft pseudopotential method of Vanderbilt [29]. It has the advantage that the exact shape of the valence wavefunctions is taken into account, and this in turn improves the description of the magnetic materials considerably [30]. The PAW potential for Ni and H are similar to those proposed and tested in Ref. [30]. The energy cut-off was set to 270 eV. In all calculations, we use the generalized gradient corrections (GGA) of Perdew and Wang [31,32] commonly refered to as PW91. The slab supercell approach with periodic boundaries is employed to model the surface, and the Brillouin zone (BZ) sampling is based on the Monkhorst– Pack technique [33]. For the calculation of the fractional occupancies, a broadening approach

G. Kresse, J. Hafner / Surface Science 459 (2000) 287–302

proposed by Methfessel and Paxton [34] is used with N=1 and s=0.2 eV. In some cases, a tetrahedron method with quadratic corrections is applied to check the results obtained with the smearing scheme [35]. The relaxation of the ionic positions into the ground state geometry is performed with either a conjugate gradient algorithm or a quasiNewton scheme [36 ]. Most calculations presented here take into account spin-polarization. This means that bulk Ni and the Ni surface have a magnetic moment. However, a few calculations were performed nonspin-polarized, i.e. the magnetic moments were forced to be zero. The reason to do so is twofold. First, using such spin-restricted calculations it is possible to investigate whether the magnetic moment has a significant impact on the adsorption geometry and adsorption energy. Second, non spin-polarized calculations are a factor of two faster, and because we intend to perform further calculations on the dissociation of H on Ni sur2 faces, they can serve as a test to determine whether faster non-magnetic calculations are sufficiently accurate for the description of the H dissociation 2 on these surfaces.

3. Results 3.1. Bulk Ni and the H dimer 2 The results for bulk Ni have already been reported elsewhere [30]. For comprehensiveness, we summarize them again in Table 1 comparing them to previous calculations [37,38] and experiTable 1 Equilibrium lattice constant, a, cohesive energy, E , bulk modcoh ulus, B, and magnetic moment, M, of ferromagnetic fcc-Ni fcc-Ni

˚ 3) a (A

E (eV ) coh

B (GPa)

M (m ) B

Expa Present PAW US-PPb FLAPWc

3.52 3.52 3.53 3.52

4.44 4.78 4.76 –

1.86 1.94 1.95 2.0

0.61 0.61 0.62 0.60

a Ref. [39]. b US-PP, Ref. [37]. c FLAPW, Ref. [38]

289

ment [39]. In the calculation of the cohesive energy, we have taken into account that the GGA ground state of the Ni atom is magnetic and nonspherical [40] (we obtained, in accordance with Ref. [40], DE =0.12 eV ). The agreement aspherical with experiment is quite good. However, even with GGA and non-spherical corrections to the atom, the cohesive energy is somewhat overestimated. For the H dimer, we find a bond length of ˚ , a2 binding energy of E=4.58 eV and a =0.750 A 0 a vibrational frequency of v=4300 cm−1. These values are typical for gradient-corrected calculations of H , and they agree reasonably well 2 ˚ , E=4.75 eV, with experiment (a =0.741 A 0 v=4395 cm−1) [41,42]. We have also checked on the convergence with respect to the plane-wave cut-off and found that an increase of the energy cut-off to 400 eV yields essentially the same binding energies (DE<10 meV ), lattice constant, and bond ˚ ). length (Dr<0.005 A

3.2. Clean surface The clean Ni surfaces have recently been studied in detail by Mittendorfer, Eichler, and Hafner [43] employing the same plane-wave code and ultrasoft pseudopotentials. In our present calculations, we use the PAW method so that we have to check on the properties of the clean surfaces again. In addition, we have tested at this point how many layers are required to evaluate surface properties ( like the magnetic moment in the surface layer, the surface energy, and the relaxation of the top most layers) with reasonable accuracy. In order to do that, we considered slabs with four and six layers [five and six layers for the (110) surface]. The results are summarized in Table 2 together with the results of Ref. [43]. In our four-layer calculations, only the topmost layer was relaxed, whereas in the five- and six-layer calculations, two surface layers were allowed to move. A k-point grid of 8×8×1 Monkhorst Pack points was found to be sufficient for all evaluated properties, but to be on the safe side, the calculations were performed with a 12×12×1 Monkhorst Pack grid for the (111) and (100) surface, and a 12×8×1 grid for

290

G. Kresse, J. Hafner / Surface Science 459 (2000) 287–302

Table 2 Results for the clean low-index surfaces of fcc Ni for 4 layers [(110) — five layers] and six layersa Ni (111) Four layers Six layers US-PPa (100) Four layers Six layers US-PPa (110) Five layers Six layers US-PPa

Dd

Dd (%) 23

M (m ) S B

s (eV )

w (eV )

−1.3 −1.2 −0.9

– 0.1 0.0

0.65 0.66 0.68

0.63 0.63 0.65

5.07 5.11

−3.3 −3.6 −3.6

– 1.0 1.4

0.71 0.71 0.76

0.84 0.83 0.85

4.97 4.97

−11.1 −10.4 −10.3

3.2 1.8 3.2

0.76 0.76

1.20 1.21 1.24

4.55 4.60

12

(%)

a Dd and Dd are the changes of the lattice spacing between the first and second layer, and between the second and third layer, 12 23 respectively. M , s, and w are the local magnetic moment of the surface atoms, the surface energies and the work-function, respectively. S b Nine layers, ultrasoft-PP, PW91, Ref. [43].

the (110) surface1. With the 8×8×1 grid, the relaxation of the topmost layer is accurate to within 0.2%, and the 12×12×1 grid yields practically converged results. The only exceptions are the magnetic moments of the surface layer, and the reported magnetic moments have been calculated with more k-points [(111) — 16×16×1, (100) — 16×16×1, (110) — 12×8×1)]. Within numerical accuracy, our six-layer calculations yield essentially the same results as the calculations of Mittendorfer et al. with nine layers and ultrasoft pseudopotentials ( US-PPs). This shows that the discrepancies between US-PPs and the PAW method are small for Ni (as already discussed in Ref. [30]), but it also indicates that six-layer slabs are certainly sufficient for a very accurate representation of the real surface. The only notable difference between the present calculation and Ref. [43] is that the magnetic moment of the surface layer increases from the close-packed to the more open surfaces [(111) — 0.65m , B (100) — 0.71m , (110) — 0.76m ]. This is exactly B B what one would expect, but in Ref. [43], this trend was not quite reproduced. We attribute these discrepancies mainly to the fact that we have 1 To obtain the same accuracy for the (110) surface as for the other two surface, we can decrease the density of k-points in the [11: 0] direction, since the cell is longer in that direction.

significantly increased the k-point grids for the evaluation of the magnetic moments. The second important result of the present calculation is that four-layer calculations for the (111) and (100) surface and five-layer calculations for the (110) surface yield results that are very close to that for six layers — at least if the clean surface is considered. We have also tested, whether four and five layers, respectively, are sufficient to describe the absorption energy of hydrogen accurately. Our findings are that the absolute adsorption energies for four layers [and five layers for (110)] agree within 10 meV with those for six layers. The relative error between different sites is, in all inspected cases, only 5 meV. We also performed tests with respect to the k-point sampling, and found that 12×12×1 and 12×8×1 k-points are sufficient to describe the energies with an accuracy of approximately 2 meV. The tests for the (111) surface are summarized in Table 3. These tests justify the following parameters: the (111) and the (100) surface are modeled by four layers. For the (110) surface, we decided to use a six-layer slab, although our convergence tests show that five layers would be adequate. The H atoms are placed asymmetrically on one side of the slab. A 12×12×1 k-point grid [12×8×1 for (110)] is chosen for calculations with primitive surface cells.

291

G. Kresse, J. Hafner / Surface Science 459 (2000) 287–302

Table 3 Adsorption energy E =E −E −1/2E(H ) of H on Ni (111) at full coverage including surface relaxation ad H/Ni-slab Ni-slab 2 (in meV/atom)a Slab thickness Relaxed k-points

Four layers Two layers 8×8×1

Four layers Two layers 12×12×1

Four layers Two layers 16×16×1

Six layers Three layers 12×12 ×1

fcc hcp Bridge

−572 −558 −343

−567 −553 −352

−568 −552 −351

−573 −562 −360

a Results for different k-point sets and slabs are compared. For the four-layer slabs, the H atom was adsorbed on one side of the slab, and for the six-layer slab, H was adsorbed on both sides of the slab.

When the surface cell is laterally doubled, the k-point grid is reduced in the corresponding direction; e.g. for a (2×2) cell, a 6×6×1 k-point grid is used. With this set-up, we estimate that the errors in the absolute adsorption energies are 10 meV, and relative errors are expected to be around 5 meV. 3.3. H adsorption on the unreconstructed surfaces at full and 1/4 H coverage The high-symmetry sites considered in the present work are shown in Fig. 1, and the results for the adsorption energies at 1/4 and full H coverage are summarized in Tables 4 and 5, respectively. The most important geometrical parameters — the height of the adatom above the surface and the shortest HMNi distance — are shown in Table 6.

For all calculations, we used as initial geometry the relaxed geometry of the clean slab, and — in addition — allowed the two topmost Ni layers [in the case of the (110) surface, the three topmost Ni layers] to relax. The adsorption energy is defined as E =E −E −1/2E( H ), ad H/Ni-slab Ni-slab 2 where we used the theoretical binding energy of H . We report the results for magnetic calculations 2 and for calculations in which the magnetic moments were forced to be zero. The later calculations do not correspond to the reality and only serve as a test to determine whether the surface magnetism has a significant influence on the adsorption energies. Indeed, we find the effect to be quite considerable: on the (111) and (100) surface, the differences between the magnetic and non-magnetic results are 70–90 and 50–70 meV, respectively. On the more open (110) surface, the

Fig. 1. Representation of the high-symmetry sites considered in the present work: TP, top; BR, bridge; LB, long bridge; SB, short bridge; HL, hollow; and PT, pseudo-threefold.

292

G. Kresse, J. Hafner / Surface Science 459 (2000) 287–302

Table 4 Adsorption energy E =E −E −1/2E(H ) of H on ad H/Ni-slab Ni-slab 2 the low-index surfaces of Ni calculated for high-symmetry sites at 1/4 coverage including surface relaxation (in meV/atom)a (111) Four layers Spin-polarized Hollow (fcc) −598 (hcp) −587 Bridge (SB) −456 (LB) Top 1 Non-spin-polarized Hollow (fcc) −692 (hcp) −676 Bridge (SB) −542 (LB) Top −74

(−588) (−576) (−446) ( 7)

(100) Four layers

(110) Six layers

−567 (−530) −450 −226 −457 (−426) −391 −409 16 ( 30) 86 −626 −520 −57

(−426)b (−217)c (−364) (−386) ( 101)

−556b −315c −521 −505 −50

a Results from non spin-polarized and spin-polarized calculations are compared (values in parentheses are for an unrelaxed fixed substrate). b Pseudo-threefold-coordinated. c Hollow site. Table 5 Adsorption energy E =E −E −1/2E(H ) of H on ad H/Ni-slab Ni-slab 2 the low-index surfaces of Ni calculated for high-symmetry sites at full coverage including surface relaxation (in meV/atom)a (111) Four layers Spin-polarized Hollow (fcc) −567 (hcp) −553 Bridge (SB) −352 (LB) Top 195 Non-spin-polarized Hollow (fcc) −685 (hcp) −662 Bridge (SB) −467 (LB) Top 88

(−557) (−539) (−341) ( 207)

(100) Four layers

(110) Six layers

−584 (−544) −443 −185 −391 (−374) −423 −406 169 ( 188) 205 −628 −448 106

(−434)b (−163)c (−401) (−393) ( 223)

−548b −285c −536 −492 71

a Results from non spin-polarized and spin-polarized calculations are compared (values in parentheses are for an unrelaxed fixed substrate). b Pseudo-threefold-coordinated. c Hollow site.

effect of magnetism even amounts to 100– 150 meV. From a quantitative point of view, magnetism obviously must be taken into account, since

Table 6 Geometrical parameters of H on the Ni surface in high-symmetry sites (magnetic calculations, full coverage)

˚) Height, h (A 0 Hollow (fcc) (hcp) Bridge (SB) (LB) Top ˚) HMNi distance, d (A 0 Hollow (fcc) (hcp) Bridge (SB) (LB) Top

(111)

(100)

(110)

0.87 0.87 1.01

0.45

1.47

1.47

0.45a 0.46b 1.02 0.13 1.47

1.67 1.67 1.60

1.81

1.47

1.47

1.01

1.60

1.66a 1.58b 1.61 1.76 1.47

a Pseudo-threefold-coordinated. b Hollow site.

neglecting it results in an overestimation of the adsorption energies by up to 20%. However, the relative energetic ordering of the high-symmetry sites and the geometrical parameters are reproduced remarkably well by non-magnetic calculations. The HMNi bond length, for instance, changes by less than 1%. Before discussing several trends in detail, we would like to compare our results with those of other density-functional studies. For the (111) surface at full coverage, Paul and Sautet [18] found an adsorption energy of E =−620 and ad −440 meV for the hollow (HL) and the bridge (BR) sites, respectively2. Their calculations were non-magnetic, employed less k-points and did not include substrate relaxation. Given their approximations, their results agree within the error bars with our calculations. The agreement with the results of Mattsson et al. for the (100) surface is slightly worse (E =−470 and −350 meV for HL and BR sites) ad [20]. The precise reason for these discrepancies is difficult to disentangle, but we have reason to believe that our present calculations are more accurate. The strongest evidence is the good agreement of the adsorption energies and adsorp2 For consistency, we have calculated the adsorption energy from the tabulated binding energies of Ref. [18] using the theoretical PW91 energy of the H dimer, E =4.58. 2 H2

G. Kresse, J. Hafner / Surface Science 459 (2000) 287–302

tion geometries for different surfaces, as discussed in Section 3.3.1. 3.3.1. Common features of high-symmetry sites on the three surfaces We concentrate first on the low-coverage case, because this allows us to avoid any complications arising from the interaction between adsorbed H atoms, a topic we will address at a later point (Section 3.3.3). Let us consider first relaxation: for the (111) surface, substrate relaxation changes the energy by about 10 meV for all sites except TP. For the (100) and (110) surface, the energy gain due to relaxation increases to about 30 meV. The relative energetic ordering of the adsorption sites, however, remains unaffected by substrate relaxation. It is also interesting to note that the energy gain due to relaxation is roughly equal for low coverage and high coverage. The second remarkable result that follows from Tables 4 and 6 is that the adsorption energies and geometries are very similar for comparable adsorption sites. This similarity is most obvious for the bridge sites [short bridge (SB) site on the (110) surface] and the top ( TP) sites. In the non-magnetic case, the adsorption energies differ only by 20 meV. Somewhat larger differences are found for the (110) surface, if magnetism is included. Here, the TP and SB sites are some 70 meV less stable than the TP and BR sites on the other two surfaces. We attribute this to the increased magnetic moment of the surface atom, which in turn reduces the adsorption energy. The hollow sites differ geometrically on the three surfaces. On the (111) surface, the H atom is coordinated to three surface atoms, whereas on the (100) surface, it is bonded to four Ni atoms. As a consequence of bond-order conservation, the HMNi bond length is slightly larger on the (100) ˚ ) than on the (111) surface surface (1.81 A ˚ ). The adsorption energies are, nevertheless, (1.67 A very similar. On the more open (110) surface, the adsorption energy in the HL site is considerably smaller than for the other two surfaces, mainly because the approach of the hydrogen atom is hindered by the subsurface Ni atom. The distance ˚ , close to the NiMH to this Ni atom is only 1.58 A ˚, bond length for atop adsorption. With 2.2 A

293

however, the distance to the other four surface Ni atoms is too large to allow for a strong bonding. Therefore, the HL site is rather unfavorable on the (110) surface. If bond length and coordination are considered, the closest analogous to the (111) HL site is the pseudo-threefold (PT ) coordinated site between the HL and the SB site (see Fig. 1). ˚, The distance to the two surface atoms is 1.66 A which is somewhat shorter than the HMNi bond length for the HL site on the (111) surface; the distance to the subsurface Ni atom is elongated to ˚ . This site is the most stable site on the 1.75 A (110) surface, but the adsorption energy is smaller than adsorption in the HL site on the (111) surface, since the bond to the subsurface Ni atom is relatively weak. Finally, the long bridge (LB) (110) site differs from other bridge sites such that the H adatom is located almost in the plane of the surface layer. Simple geometrical arguments are again responsible for this behavior: here, the H atom has four nearest neighbors, two surface and two subsurface atoms, and the HMNi distances are 1.76 and ˚ , respectively. These values are close to the 1.80 A HMNi bond length for the fourfold hollow site on the (100) surface. 3.3.2. Trends in adsorption energies, and diffusion barriers Table 4 clearly shows that the adsorption energies in the most favorable sites are very similar for the (111) and (100) surface, at low coverage. This agrees with the experimental observation: adsorption energies of −470 meV/H adatom (21.5± 1.0 kcal/mol H ) [2] and −500 meV/H adatom 2 (23 kcal/mol H ) [4–6 ] have been reported for the 2 (111) surface and (100) surface, respectively. Our results agree reasonably well with the experimental values, although the theoretical values are slightly too large. Recent systematic calculations for the adsorption of small molecules on different transition-metal surfaces indicate that this is a deficiency of the PW91 functional (see, for example, Ref. [44]). For the (110) surface, no adsorption energies have been reported. Since the surface undergoes various surface reconstructions upon heating, it is not possible to derive them from temperature-

294

G. Kresse, J. Hafner / Surface Science 459 (2000) 287–302

programmed desorption experiments (micro-calorimetry data are also not available). Our results indicate that the adsorption energies on the unreconstructed (110) surface are significantly lower than on the other low-index surfaces. Diffusion of H atoms on the surface is another important issue that we would like to address briefly. Although, for a full understanding, a quantum mechanical treatment of the H atoms would be required, we can make at least a few tentative comments. On the (111) surface, the diffusion barrier is 150 meV (fcc BR); on the more open (100) surface, the barrier drops to 110 meV (HL BR). For the (110) surface, we have performed thorough scans of the full HMNi potential energy surface in the preparation of work on the dissociation of molecular hydrogen. These investigations show that the LB site is a local minimum. Diffusion from the PT to the LB site is accomplished almost on a straight line connecting both sites (compare also Ref. [45]). The diffusion barrier along this path is only 100 meV. With 60 meV, the diffusion barrier from the PT via the SB to the next PT site is even lower, which indicates that H atoms can move freely between troughs. We can conclude that all three low-index surfaces exhibit rather small H-diffusion barriers, but contrary to what one might expect, the most corrugated surface exhibits the lowest barriers. Before continuing, we want to make one further remark. At full coverage, our calculations yield quantitatively somewhat different results for the diffusion barriers. For the (111) and (100) surfaces, barriers of 210 and 190 meV, respectively, are obtained. These values are almost 50% higher than those at low coverage. This means that diffusion barriers evaluated at high coverage are not completely reliable. The reason is that the repulsion between H adatoms varies significantly between different adsorption sites, which increases the energy at the transition state (BR sites) compared to the HL sites. 3.3.3. Electrostatic repulsion between H adatoms and work function To facilitate a simple comparison between the low- and high-coverage cases, the differences in the adsorption energies are summarized in Table 7.

Table 7 Energy difference (in eV/atom) between high-coverage and lowcoverage cases (magnetic calculations, relaxed substrate)

Hollow (fcc) (hcp) Bridge (SB) (LB) Top

(111)

(100)

(110)

0.03 0.03 0.10

−0.02

0.20

0.15

0.01a 0.04b −0.03 0.00 0.12

0.07

a Pseudo-threefold-coordinated. b Hollow site.

Two trends become clear upon inspection of the results. First, the repulsion is strongest for the (111) and smallest for the (110) surface. The trend is particularly evident for the TP site where the repulsion decreases from 0.20 eV (111) to 0.15 eV (100) and finally to 0.12 eV (110). Second, the repulsion is strongest for those hydrogen atoms that are a large distance away from the surface. This behavior can be reconciled easily with an electrostatic model. Charge flows from the surface towards the H adatom, causing a dipole moment and, in turn, an electrostatic repulsion between the H atoms. The resulting repulsion is greatest for the (111) surface where the packing of H atoms is densest, and least for the (110) surface with the smallest packing ratio. Screening explains why the repulsion is stronger for the TP site than for the HL site: if the distance of the H adatom from the surface is large, the surface electrons cannot screen the resulting dipole moment, and the repulsion is large ( TP site). However, if the H atom is located close to the surface, the metal electrons screen the accumulated charge on the H adatom, reducing the repulsion between them. The repulsion is therefore particularly low for the hollow sites on the (100) surface and the LB site on the (110) surface. For the HL(100) site, we find even a small attraction between the H atoms. Interestingly, we also find no repulsion between the SB sites and between the PT sites on the (110) surface, although we would have expected repulsion on electrostatic grounds alone. The attraction between the atoms in the HL(100) sites and SB(110) sites can be caused only by band-structure-related energy gains.

G. Kresse, J. Hafner / Surface Science 459 (2000) 287–302

In an attempt to gain a more qualitative insight on the electrostatics, we have evaluated the difference electron densities Dr=r −r −r , H/Ni-slab Ni-slab H where r , r and r are the electron H/Ni-slab Ni-slab H densities for H adsorbed on the slab, for the clean slab and for an isolated H atom, respectively. For four selected configurations, iso-surfaces of the difference electron density are shown in Fig. 2. Our findings are that charge always accumulates on the H adatom. When the charge density difference is integrated for a sphere with a radius of ˚ around the H atom, a charge transfer of 0.9 A 0.15–0.20 electrons to the H adatom is found. The

295

charge transfer is generally largest for the top and bridge sites (0.19–0.20) and somewhat smaller for the hollow sites (0.15–0.17). These values fit rather well with the observed repulsive interaction between the H atoms adsorbed atop: assuming a charge transfer of 0.19 electrons from the surface Ni atom to the H atom and using a simple pointcharge model, repulsions of 0.22, 0.19 and 0.10 eV for the (111), (100) and (110) surface at full coverage are obtained. This is in almost perfect agreement with our exact first-principles results. For the bridge site, the agreement is more approximate; values of 0.11, 0.10 and 0.06 eV, are obtained, assuming a point charge on the H atom and a second charge below the adatom in the

Fig. 2. Iso-surfaces of the difference electron density [r −r −r ] for (a) H in the fcc (111) HL site, (b) and the (100) HL H/Ni-slab Ni-slab H site, (c) the (100) BR site, and (d) the (100) TP site. Views are from the side. Charge flows from the dark into the light regions. Small black balls indicate the position of the H and Ni atoms, and the H atom is always enclosed in the light region.

296

G. Kresse, J. Hafner / Surface Science 459 (2000) 287–302

surface plane. Although the simple point-charge model is certainly oversimplified, the essential features responsible for repulsion in the TP and BR site are obviously captured by this picture. The second quantity that is affected by the charge accumulation on the H adatom is the workfunction. On the basis of the simple point charge model, one would expect the change in workfunction to greatest for the TP site and smallest for the HL sites. This is indeed what we find. At first sight, one would also anticipate that due to packing, the work-function change is larger for the (111) surface than for the (100) surface. For atop adsorption at full coverage, the oversimplified picture now, however, fails. The calculated workfunction changes are 0.15, 0.54 and 1.40 eV for the (111), (100) and (110) surfaces contrary to our expectation, and in addition, the values are much smaller than expected on the basis of the point-charge model (#3–5 eV ). The reason is that the real charge distribution is significantly more complicated than had been assumed so far. Fig. 3 shows the difference electron densities, Dr(z), integrated over planes of constant height over the surface

P

Dr(z)= Dr(x, y, z) dx dy

(1)

as a function of the distance from the Ni surface layer. The z-axis is chosen to be perpendicular to the surface, positive z-values point away from the slab, and negative values correspond to positions inside the slab. The induced surface dipolemoment, p, is then defined as

P

p= Dr(z)z dz, and is essentially responsible for the change in work-function. We have considered three adsorption sites: (111) HL, for which the discussion is postponed to the Section 3.4, (111) TP and (110) TP. For the (111) TP and (110) TP sites, the charge redistribution looks qualitatively remarkably similar: below and above the surface atom, charge has been depleted as a result of the depletion of the Ni d orbital. In addition, a charge z ˚ above depletion is visible in the vacuum region 3 A

Fig. 3. Difference electron density [r −r −r ] integH/Ni-slab Ni-slab H rated for a plane parallel to the surface as a function of the distance from the topmost layer for (a) H in the (111) HL site, (b) (111) TP site, and (c) the (110) TP site.

the surface, which is due to the fact that the center of mass of the charge density on the H adatom is shifted towards the surface. The later effect is also clearly visible in Fig. 2d. Both effects — the charge transfer to the H adatom and the shift of the center of the charge density towards the surface — counterbalance each other. In turn, the workfunction change is generally small: a calculation of the dipole moment arising from the difference electron density shows that the cancellation is almost perfect for the TP(111) site, whereas for the TP(110) site, the mutual cancellation is not so strong, which is in agreement with the calculated work-function changes. As already pointed out before, the work-function change is smaller for the stable adsorption site than for the TP site. At full coverage, we found values of 0.10, 0.07 and 0.06 eV for the

G. Kresse, J. Hafner / Surface Science 459 (2000) 287–302

HL(111), HL(100) and PT(110) sites. These values turn out to be very sensitive to k-point sampling, slab thickness and vacuum width. Even though we used symmetric slabs with six layers ˚ of vacuum for the evaluation of and at least 10 A the work function, the errors might be up to 0.02 eV. Agreement with experiment is, however, reasonable. Extrapolation of the experimental values to full coverage yields 0.08 eV for the (111) surface [2], and 0.06 [7] and 0.11 [46 ] for the (100) surface. For the (111) surface, we will discuss the work-function changes in more detail in Section 3.4. 3.4. Ni(111): the (2×2)-2H superstructure The (2×2)-2H superstructure on Ni(111) is known to be stable up to 200 K at a coverage of 1/2 ML. The precise atomic arrangement of the H adatoms was unknown for some time and has finally been identified by accurate T-LEED measurements [2,3]. Interestingly, the later measurements show a considerable degree of surface buckling, which we would like to confirm with theoretical methods [3]. Our calculations were performed in a (2×2) supercell with four layers and a 6×6×1 Monkhorst–Pack k-point grid. When the substrate was allowed to relax, we found an adsorption energy of −620 meV/ H atom. Interestingly, this is 20 meV lower than the adsorption energy at 1/4 coverage and 50 meV lower than the adsorption energy at full coverage, which relates to the experimental observation that the (2×2)-2H superstructure is observed in a rather large temperature and concentration regime. It is interesting to compare this with the situation on the (100) surface. Here, we found that the adsorption energies were −567, −565 and −584 meV for 1/4, 1/2 [c(2×2)] and 1 ML H coverage, respectively. Clearly, the c(2×2) superstructure is not stable for this surface. For the (2×2)-2H/Ni(111) phase, our final ˚ , dB= structural parameters are dH=0.92 A ˚ , d =2.044 A ˚ , d =2.030 A ˚ , where dH is 0.08 A 12 23 the height of H atoms above the surface layer, and dB is the buckling in the first layer and d 12 and d are the inter-layer spacings (for details on 23

297

how the structural parameters are defined, we refer to Ref. [3]). In contrast to the clean surface, the first layer is now relaxed outward by approximately Dd /d =0.5%, and the second layer 12 0 remains almost unaffected by the hydrogen adsorption (Dd /d =0%). Our values agree well with 23 0 those determined by T-LEED [3] (dH=0.98± ˚ , dB=0.04±0.02 A ˚, ˚, 0.08 A d =2.03±0.03 A 12 ˚ d =2.05±0.05 A) and confirm that the surface is 23 buckled. In the theoretical calculations, the buckling parameter is even somewhat larger than that determined from experiment, and all other structural parameters agree within the error bars with experiment. One aspect that has gained a lot of experimental attention, is the behavior of the work-function. Experimentally, an almost linear increase is observed up to 0.5 ML coverage [2]. Then, the work function decreases. This observation was initially attributed to two different adsorbed species: the increase is due to the occupation of a positively polarized species and the decrease due to the occupation of a second negatively polarized state [2]. In our calculations, we were able to observe precisely the same dependency, and the work-function changes are 80, 160, and 100 meV at 1/4, 1/2 and full coverage, which is in almost perfect agreement with experiment [2]. However, in our case, only the hollow sites were occupied, making the non-linear behavior even more astonishing. For a better understanding, we turn again to the induced charge density Dr(z) [Eq. (1)] shown in Fig. 3. As for the TP sites, a significant accumulation of charge is found on the H adatom, and again, part of the charge stems from the surface and part from the vacuum. Per H adatom, the charge density flow is very similar for 1/4, 1/2 and full ML coverage, so that we show only results for full coverage. When the surface dipole moment is calculated from the induced charge density, its value first increases towards a maximum and then reduces by almost two orders of magnitude because of the charge depleted from the vacuum. This shows that the work-function is extremely sensitive to the exact shape of the induced charge density, and this sensitivity must be responsible for the observed non-linearity. However, the complex charge rearrangement unfortunately also makes it

298

G. Kresse, J. Hafner / Surface Science 459 (2000) 287–302

difficult for us to determine the precise reason for the observed behavior. 3.5. Ni(110) surface 3.5.1. (2×1)-2H phase, at h=1 ML The structural behavior of the (110) surface under H exposure was already briefly discussed in Section 1. At least three different stable structures have been observed experimentally: (a) a (2×1)-2H phase, in which the H atoms are believed to occupy the PT sites (Fig. 4a), (b) a (1×2)-3H pairing-row reconstruction (Fig. 4b and c), and (c) the streak phase. In order to investigate the relative stability of these phases, we have explored several structural models. The most

Fig. 4. Top view of the explored structures of H on Ni (110): (a) (2×1)-2H lattice gas structure (no reconstruction), (b), (c) and (d) three models for the (1×2)-3H pairing-row reconstruction, and (e) (1×2)-3H missing row reconstruction.

important models are depicted in Fig. 4. The first model (Fig. 4a) is the unreconstructed (2×1) lattice gas phase, which is stable only at rather low temperatures and coverage. Including substrate relaxation, we find an adsorption energy of −485 meV per H atom, which is lower than the (1×1) structure explored at full coverage in Section 3.3.3 (−450 meV ). 3.5.2. (1×2)-3H phase, at h=1.5 ML To study the reconstructed phase, we first explored the two models depicted in Fig. 4b and c, which have both been suggested on the basis of experimental measurements (see, for example, Ref. [12]). In the calculation, a (1×2) supercell with six layers was used, and the top four layers were allowed to relax3. The k-point sampling was performed with a 12×4 Monkhorst Pack grid (corresponding to 12×8 in the primitive cell ). We find that model (c) is unstable and relaxes to model (b). After relaxation, the substrate is reconstructed for model (b), and the H adsorption energy is found to be −444 meV/H. Compared to the 1 ML coverage, the adsorption energy has only slightly decreased, and the surface energy per Ni atom is approximately 200 meV lower than for the (2×1)-2H phase. This shows that the reconstruction is an exothermic process if the surface is exposed to H . We have also considered two 2 starting models similar to (b) and (c), but with one H adatom initially in the LB site instead of the HL site. After relaxation, these two models converge to the same geometry, but the final structure was approximately 50 meV per H-adatom less stable than model (b). Finally, model (d ), which has not been proposed before in the literature, was investigated. Here, two atoms sit in adjacent PT sites, and the third H atom is located in the SB site. This arrangement allows for an efficient screening between the H atoms and a rather larger mutual HMH distance. Due to the symmetry, this model shows no reconstruction. With an adsorption energy of −450 meV per H atom, our calculation predicts that this model is 3 To obtain reliable results for the reconstructed surface, we had to relax four layers. For unreconstructed structures, it is sufficient to include the relaxation of three layers.

G. Kresse, J. Hafner / Surface Science 459 (2000) 287–302

slightly more stable than the pairing-row reconstruction (b). Before discussing this unexpected result in more detail, we would like to compare the results for our pairing-row structure [model (b)] with experiment in more detail. The work-function for model (b) is 400 meV larger than for the clean surface, in good agreement with experiment (510 meV according to Ref. [15]). The final structural parameters that we obtain for ˚ , d =1.27 A ˚ , BU= model (b) are: d =1.26 A 12 ˚ 23 ˚ 0.26 A and LS=0.34 A, where d and d are the 12 23 layer spacing, BU is the buckling observed in the subsurface Ni layer, and LS is the lateral displacement of the paired rows with respect to their ideal bulk terminated positions (for a more thorough definition of the structural parameters, we refer to Ref. [16 ]). The layer spacings, d and d , corre12 23 spond to an outward relaxation of 1.0, 2.0% of the first and second layer with respect to the bulk. Our values agree very well with the structural ˚, parameters determined by T-LEED: d =1.27 A 12 ˚ ˚ , BU=0.25 A ˚ and LS=0.30 d =1.31 A A, which 23 corroborates that model (b) corresponds indeed to the structure observed experimentally. The final positions of the H atoms adatoms are particularly interesting since, experimentally, their determination is very difficult. In fact, neon diffraction experiments [12] have favored model (c). However, our calculations exclude with certainty that this structure is stable. In addition, taking into account the strong electrostatic repulsion between H adatoms, even common sense implies that the distance between the H adatoms is too small in the proposed model. In our final relaxed model (b), the positions of the H adatoms are rather similar to those on the unreconstructed surface: in the PT site, the distance between the H atom and the surface Ni atom is ˚ , and the bond length to the subsurface Ni 1.67 A ˚ . In the HL site, the H atom is now atom is 1.75 A ˚ above the subsurface Ni atom, and located 1.76 A the distance between the H and the surface Ni ˚. atoms is 1.90 A The stability of the pairing-row reconstruction is driven by two major factors. Firstly, electrostatics: as discussed in Section 3.3.3, the H atoms carry a charge of approximately 0.2 e−. In the troughs, there is no charge to screen the charge

299

density accumulated on the H atoms, so the two H atoms in the PT sites experience a significant repulsion. The pairing-row reconstruction allows the H adatoms to increase their distance by almost ˚ , in turn reducing the repulsive interaction. 0.4 A Also important is the fact that the adsorption energy in the HL site is increased by the reconstruction to −380 meV (to be compared with −220 meV for the unreconstructed surface). The increase is mainly due to the reduction of the bond length between the H adatom and the four Ni ˚ (unreconstructed sursurface atoms from 2.1 A ˚ face) to 1.9 A (reconstructed surface). In fact, ˚ is quite close to the HMNi bond length in 1.9 A the fourfold hollow site on the (100) surface. In addition to the pairing-row reconstruction, we have also explored the possibility of a missing row reconstruction, as shown in Fig. 4e. Since this reconstruction requires a considerable Ni mass transport, it does not compete directly with the phases discussed so far. Our findings are that the missing row structure is in fact even more stable than either model (b) or (d ); per H adatom, the adsorption energy increases from −450 to −500 meV. It is rather tempting to associate the missing row reconstruction with the experimentally observed ‘streak phase’, but more models — for instance, H in subsurface positions — or larger facets must be explored before a definite answer can be given. Although the structural parameters of the pairing-row reconstruction agree well with experiment, the results for the energetics remain controversial. Our zero temperature results clearly favor model (d ) over model (b). This could be attributed to technical errors, but increasing the layer thickness from six to eight layers or increasing the k-point mesh to 16×6×1 gives the same result within a few meV. A possible solution to this problem is that we might have overlooked other plausible models. To exclude this, we performed simulated annealing runs and finite temperature MDs, which have the advantage that no precise structural input is required. For transition-metal surfaces, such calculations are still very rare, since they are believed to be rather expensive. Our calculations indicate that modern first-principles methods are now sufficiently fast to perform such calculations

300

G. Kresse, J. Hafner / Surface Science 459 (2000) 287–302

Table 8 Parameters for the simulated annealing (cf. text) Run

Layers

Mobile

k-points

Spin-polarized

1 2 3 4

5 5 6 6

3 4 4 4

8×4×1 8×4×1 8×4×1 12×4×1

No No No Yes

routinely. We performed a total of four runs, for which the important parameters are summarized in Table 8. The k-point sampling was performed with 10 and 12 k-points that were derived from a 8×4×1 and 12×4×1 Monkhorst–Pack grid for the high-symmetry case. Our tests indicate that this approximation is well founded, typical errors due to the restricted k-point set being only around 10 and 5 meV, respectively. To save computer time, spin-polarization was taken into account only in one run. All simulations were started at an initial temperature of 800 K, and the temperature was gradually decreased to 300 K; a total of about 1500 molecular dynamic steps were performed in each case (for the spin-polarized calculation, even 2500 steps were performed ). At the end, the systems were quenched into the nearest local minimum. With the exception of the spin-polarized calculation, each calculation took less than 3 days on a workstation, which shows that such calculations are perfectly feasible. Not unexpectedly, the MD simulation, in which three layers were allowed to move, ended up in structure (d ). Static calculations with the same technical set-up (five layers, three layers relaxed ) show that model (d ) is 80 meV more stable than model (b) with these technical parameters. But interestingly, all other simulations ended in model (b). This can be explained only by assuming that an entropy contribution stabilizes the pairing-row reconstruction over the unreconstructed model since, with these technical parameters, model (d ) is more stable than the reconstructed model (b) at zero temperature (see above).

4. Discussions and conclusions We have presented extensive first-principles calculations for the adsorption of H on the three low-

index Ni surfaces. The calculations clearly show that surface magnetism is essential for an accurate quantitative description of the energetics. In nonspin-polarized calculations, the H adsorption energies are overestimated by 100 meV. However, the geometries and the relative ordering of the highsymmetry sites are reproduced quite well, even in non-spin-polarized calculations. Our results further show that the adsorption geometries and adsorption energies are remarkably similar on the two close-packed low-index surfaces. The most stable adsorption site in both cases is the HL site, and the theoretical adsorption energy is approximately 25 kcal/mol of H . On the more 2 open (110) surface, the adsorption energy is significantly smaller, and the most stable site is the pseudo-threefold site. The energy differences between the most stable sites and the bridge sites were used to estimate the diffusion barriers. Our findings are that the diffusion barriers are low on all three investigated surfaces, but contrary to intuition, the barriers are smallest on the corrugated (110) surface. Here, a maximum barrier of 100 meV is encountered for diffusion parallel to the troughs from the pseudothreefold site to the meta-stable long bridge site, and the barrier for diffusion between the troughs almost vanishes. For the other two surfaces, the diffusion barriers are 110 meV (100) and 140 meV (111). We can explain this trend by observing that the adsorption energies in the bridge sites are very similar for all three surfaces, whereas adsorption in the most stable site yields a larger energy gain on the close-packed surface than on the corrugated surface decreasing the diffusion barriers. We have investigated the repulsion between H adatoms in detail. Our results indicate that repulsion is mainly due to the electrostatic interaction between H adatoms, which carry a negative net charge of approximately 0.2 e− once they are adsorbed on the Ni surface. Only if the H adatom is close to the surface can the surface electrons screen the charge on the H adatoms, which in turn reduces the repulsion between the adatoms. For the most stable sites on the (100) and (110) surface, this is the case, i.e. the height of the H ˚ , and the adsorption adatom is less than 0.5 A energy is essentially independent of coverage.

G. Kresse, J. Hafner / Surface Science 459 (2000) 287–302

Since, on the (111) surface, the height of the H ˚ even in the HL site, a significant adatom is 0.9 A repulsion between adsorbed H adatoms is found, which decreases the adsorption energy at full coverage by approximately 50 meV per H adatom. In addition, we find that the H atoms are 20 meV more strongly bound in the (2×2)-2H reconstruction than in the (2×2)-1H structure, which correlates well with the experimental observation that the (2×2)-2H phase is stable over a wide temperature and pressure regime. For the (111) surface, we have investigated the work-function with varying H coverage. In agreement with experiment, the work-function increases by 150 meV up to 1/2 ML coverage and then decreases by 50 meV (full coverage). The precise reason for this behavior is difficult to determine, since a fairly significant charge redistribution is observed on the surface, and even very subtle changes of the charge density cause a significant change in the work-function. However, the essential point is that only hollow sites were occupied in our calculations. An initial increase and successive decrease of the work-function with coverage is therefore not necessarily related to the occupation of positively and negatively charged species, as often concluded from experiments. In the final section of the present work, we have concentrated on the pairing-row reconstruction of the (110) surface. The work-function change and the structural parameters of the final reconstructed model agree well with experiment, but our calculations go beyond experiment, since the H positions can be determined. In agreement with the pairingrow reconstruction of the Pd (110) surface [45,47], the H atoms are found to be located in the PT and HL sites. The pairing-row reconstruction allows the H atoms in the PT sites to maximize their distance and, at the same time, the reconstruction increases the adsorption energy of the H atom in the HL site. In addition to the static calculations, which require a reasonable structural input, we have applied simulated annealing and finite temperature molecular dynamics to the (1×2)-3H supercell. To the best of our knowledge, this is the first time that simulated annealing has been combined with first-principles calculations for a transition-metal

301

surface, and our results are encouraging in many respects. First, the calculations are fast: typically one simulated annealing calculation took only 1 week on a simple workstation. The construction of reasonable (1×2) trial models and their relaxation required approximately the same time but much more ‘human’ control. Since the simulated annealing approach is also not biased by the structural input, this route is even accessible if precise structural models are not available. In our case, the finite temperature simulations also gave a new physical insight. Although we found one structural model with a slightly lower zero temperature energy than the pairing-row reconstruction, molecular dynamics favor the pairing-row reconstruction at finite temperature. This indicates that vibrational effects — zero point vibrations at very low temperatures — might play a crucial role in the formation of the reconstruction.

Acknowledgements We thank A. Eichler for helpful discussions and comments. This work has been supported by the Austrian Science Funds within the Joint Research Project on ‘Gas–Surface Interaction’ Grant No. S8106-PHY.

References [1] J. Behm, K. Christmann, G. Ertl, Solid State Comm. 25 (1978) 763. [2] K. Christmann, R. J. Behm, G. Ertl, M.A. Van Hove, W.H. Weinberg, J. Chem. Phys. 70 (1979) 4168. [3] L. Hammer, H. Landskron, W. Nichtl-Pecher, A. Fricke, K. Heinz, K. Mu¨ller, Phys. Rev. B 47 (1993) 15969. [4] J. Lapujoulade, K.S. Neil, Surf. Sci. 35 (1973) 288. [5] K. Christmann, G. Ertl, O. Schober, Surf. Sci. 40 (1973) 61. [6 ] K. Christmann, O. Schober, G. Ertl, M. Neumann, J. Chem. Phys. 60 (1974) 4528. [7] K. Christmann, Z. NaturForsch. 34a (1979) 22. [8] K.H. Rieder, H. Wilsch, Surf. Sci. 131 (1983) 245. [9] S. Andersson, Chem. Phys. Lett. 55 (1978) 185. [10] V. Penka, K. Christmann, G. Ertl, Surf. Sci. 136 (1984) 307. [11] W. Reimer, V. Penka, M. Skottke, R.J. Behm, G. Ertl, Surf. Sci. 186 (1987) 45. [12] G. Parschau, B. Burg, K.H. Rieder, Surf. Sci. 293 (1993) L830.

302

G. Kresse, J. Hafner / Surface Science 459 (2000) 287–302

[13] C.W. Tucker, Surf. Sci. 26 (1971) 311. [14] K.H. Rieder, T. Engel, Phys. Rev. Lett. 43 (1979) 373. [15] K. Christmann, F. Chehab, V. Penka, G. Ertl, Surf. Sci. 152–153 (1985) 356. [16 ] G. Kleinle, V. Penka, R.J. Behm, G. Ertl, W. Moritz, Phys. Rev. Lett. 58 (1987) 148. [17] A. Grigo, D. Badt, H. Wengelnik, H. Neddermayer, Surf. Sci. 331–333 (1995) 1077. [18] J.-F. Paul, P. Sautet, Surf. Sci. 356 (1996) L403. [19] H. Yang, J.L. Whitten, J. Chem. Phys. 98 (1993) 5039. [20] T.R. Mattsson, G. Wahnstro¨m, L. Bengtsson, B. Hammer, Phys. Rev. B 56 (1997) 2258. [21] W. Kohn, L. Sham, Phys. Rev. A 140 (1965) 1133. [22] R.O. Jones, O. Gunnarsson, Rev. Mod. Phys. 61 (1989) 689. [23] R. Car, M. Parrinello, Phys. Rev. Lett. 55 (1985) 2471. [24] M.C. Payne, M.P. Teter, D.C. Allan, T.A. Arias, J.D. Joannopoulos, Rev. Mod. Phys. 64 (1992) 1045. [25] G. Kresse, J. Hafner, Phys. Rev. B 48 (1993) 13115. [26 ] G. Kresse, J. Furthmu¨ller, Comput. Mater. Sci. 6 (1996) 15. [27] G. Kresse, J. Furthmu¨ller, Phys. Rev. B 54 (1996) 11169. [28] P.E. Blo¨chl, Phys. Rev. B 50 (1994) 17953. [29] D. Vanderbilt, Phys. Rev. B 41 (1990) 7892. [30] G. Kresse, D. Joubert, Phys. Rev. B 59 (1998) 1758. [31] Y. Wang, J.P. Perdew, Phys. Rev. B 44 (1991) 13298. [32] J.P. Perdew, J.A. Chevary, S.H. Vosko, K.A. Jackson, M.R. Pederson, D.J. Singh, C. Fiolhais, Phys. Rev. B 46 (1992) 6671.

[33] H.J. Monkhorst, J.D. Pack, Phys. Rev. B 13 (1972) 5188. [34] M. Methfessel, A.T. Paxton, Phys. Rev. B 40 (1989) 3616. [35] P.E. Blo¨chl, O. Jepsen, O.K. Andersen, Phys. Rev. B 49 (1994) 16223. [36 ] W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling, Numerical Recipes, Cambridge University Press, New York, 1986. [37] E.G. Moroni, G. Kresse, J. Hafner, J. Furthmu¨ller, Phys. Rev. B 56 (1997) 15629. [38] H.C. Herper, E. Hoffmann, P. Entel, private communication. [39] C. Kittel, Introduction to Solid State Physics, sixth ed.Wiley, New York, 1986. [40] P.H.T. Philipsen, E.J. Baerends, Phys. Rev. 54 (1996) 5326. [41] K.P. Huber, G. Herzberg, Molecular Structure and Molecular Spectra IV: Constants of Diatomic Molecules, Van Rostrand-Reinhold, New York, 1979. [42] K.P. Huber, in: D.E. Gray ( Ed.), American Institute of Physics HandbookMcGraw-Hill, New York, 1979.. [43] F. Mittendorfer, A. Eichler, J. Hafner, Surf. Sci. 423 (1999) 1. [44] B. Hammer, L.B. Hansen, J.K. Norskov, Phys. Rev. B 59 (1999) 7413. [45] W. Dong, V. Ledentu, P. Sautet, G. Kresse, J. Hafner, Surf. Sci. 377–379 (1997) 56. [46 ] B.E. Koel, D.E. Peebles, J.M. White, Surf. Sci. 125 (1983) 709. [47] D. Tomanek, S. Wilke, M. Scheffler, Phys. Rev. Lett. 79 (1997) 1329.