A systematic density functional theory study of the electronic structure of bulk and „001... surface of transition-metals carbides

A systematic study of the bulk and surface geometrical and electronic properties of a series of transition-metal carbidesTMC with TM=Ti, V, Zr, Nb, Mo, Hf, Ta, and Wd by first-principles methods is presented. It is shown that in these materials the chemical bonding is strongly covalent, the cohesive energies being directly related to the bonding-antibonding gap although the shift of the center of the Cs2sd band related peak in the density of states with respect to diamond indicates that some metal to carbon charge transfer does also take place. The s001d face of these metal carbides exhibits a noticeable surface rumpling which grows along the series. It is shown that neglecting surface relaxation results in very large errors on the surface energy and work function. The surface formation induces a significant shift of electronic energy levels with respect to the corresponding values in the bulk. The extent and nature of the shift can be understood from simple bonding-antibonding arguments and is enhanced by the structural rippling of this surface. © 2005 American Institute of Physics . fDOI: 10.1063/1.1888370 g


I. INTRODUCTION
Transition-metal carbides exhibit broad and amazing physical and chemical properties. 1 Their properties may be viewed as resulting from a combination of those of covalent solids, ionic crystals, and transition metals. Thus, they exhibit extreme hardness and brittleness as diamond and other covalent solids, very high melting points and simple crystal structures such as many ionic compounds, and electrical and thermal conductivities typical of metals. The curious behavior of metal carbides makes these materials very attractive from the fundamental and technological point of views. 1 For example, transition-metal carbides ͑TMC͒ are being used increasingly in heterogeneous catalysis because in many aspects they display a chemical behavior which is reminiscent of platinum and other transition metals such as Pd, Ru, or Rh but in addition exhibit important advantages over these bulk transition metals in activity, selectivity, and resistance to poisoning, especially by sulfur. 2 In some reactions, the TMC can also equal or surpass the catalytic activity of the precious metals. They can catalyze reactions of hydrogenation and dehydrogenation, hydrolysis isomerization, hydrodesulfurization, and hydrodenitrogenation. 3 The electronic structure and chemical bonding in bulk transition-metal carbides have been studied by experimental techniques, mainly X-ray photoemission 4 and more recently by near edge X-ray absorption fine structure ͑NEXAFS͒ spectroscopies 5 and also by theoretical methods. However, the amount of information is quite limited and sometimes contradictory. A key issue of these studies concerns the degree of charge transfer between metal and carbon atoms. 6,7 All studies indicate a certain amount of charge transfer from the metal to the carbon which is evidenced by clear shifts on the relevant bands. However, it is difficult to relate these shifts with the actual extent of the charge transfer. 8 Gruzalski and Zehner have shown that in bulk TaC and HfC, the direction of the metal core-levels shifts is different. 9 Surface corelevel shifts with respect to the bulk have been observed on most metal surfaces, [10][11][12][13] although there is almost no information for the corresponding carbides except for TiC. For TiC, density functional theory ͑DFT͒ calculations based on the local density approximation ͑LDA͒ suggest that there is no bulk-surface core-level shift for the C͑1s͒ level and that for the Ti͑2s͒ level it is of ϳ0.05 eV only. 14 However, the predicted lack of a shift for the C͑1s͒ core level is not consistent with subsequent photoemission experiments where a small surface-to-bulk core-level shift has been measured from the deconvolution of the experimental peak. 15 Careful theoretical analyses reveal that the final surface-to-bulk corelevel shift arises from the cancellation of various effects and, hence, is quite often of ϳ0.5 eV or even smaller. 8,16 Therefore, more precise experimental data would require the use of high-resolution core-level spectroscopy exploiting the properties of synchroton radiation.
Concerning the nature of the chemical bond, different bonding mechanisms have been discussed by various authors. In his comprehensive review, Chen 3 comments that the ionicity decreases along the series-ZrC more ionic than MoC-and that the bonding involves three main contributions. A metallic one arising from the rearrangement of the metal-metal bonds, a covalent part due to the formation of typical chemical bonds between metal and nonmetal atoms, and a third contribution from ionic bonding arising from the metal-to-carbon charge transfer. However, it is very difficult to extract information concerning the extent of each of these contributions by experimental techniques. This kind of information can better be obtained by electronic structure calculations. Based on full-potential linear muffin-tin orbital ͑FP-LMTO͒ DFT calculations within the generalized gradient approximation ͑GGA͒, Djellouli and Aourag 17 decompose the density of states ͑DOS͒ of TiC into three well separated regions, one clearly dominated by the C͑2s͒ orbitals, one bonding region ͑hybridization of d-metal orbitals with p-carbon orbitals͒, and the corresponding antibonding region just above the Fermi level, dominated by the metal d states but with small degree of hybridized p-carbon states. These results are also similar to those reported earlier by Häglund et al. [18][19][20][21][22][23][24] and more recently by Hugosson et al., 25,26 the later focusing as well on the electronic structure of metal carbide surfaces. These studies are also based on the use of the LMTO basis but within the LDA method. Furthermore, Djellouli and Aourag 17 propose a model for TiC where the chemical bond is different depending on the crystal direction: mainly covalent in the ͓001͔ but ionic along the ͓111͔ directions with significant charge transfer from the metal to the carbon. This is difficult to rationalize due to the symmetric rock-salt structure of TiC. However, this interpretation arises from the analysis of the electron density maps, a property which provides important qualitative information but which otherwise is difficult to use in a more quantitative way. Zhukov and Gubanov 27 use a similar approach and decompose the bulk modulus B 0 in contributions of the main features of the band structure. Zhang et al. have reported a study of the decomposition of cohesive energies of 3d carbides on different bonding contributions. 28 Finally, for metal carbides of cubic structure, a direct correlation between core-level shifts and B 0 has been reported. 29 Transition-metal carbide surfaces have been characterized experimentally through various surface science techniques ͑LEED, XPS, STM, etc.͒. The stability of the three low-index surfaces has been studied in detail for TiC 30 and it was concluded that the ͑001͒ surface is the most stable one, the ͑111͒ surface is metal terminated and that the ͑110͒ surface does not exist because of the formation of microscopic ͑310͒ oriented facets. Early impact collision ion scattering spectroscopy ͑ICISS͒ experiments carried out for TiC͑001͒ did not reveal significant rippling of the surfaces 31 although later experiments on TaC͑001͒ 32 show a clear displacement of the carbon atoms outwards to the vacuum accompanied by an inwards relaxation of the Ta atoms. Recent LEED experiments on TiC͑001͒ have also found evidence of a significant relaxation 33 and the same trend is observed in highresolution medium-energy ion scattering spectroscopy carried out for both TiC͑001͒ and TaC͑001͒. 34 This surface rum-pling has later been confirmed by LDA calculations using a FP-LMTO basis set. 35,36 Further experimental efforts have focused on determining if rippling takes place on other carbides such as HfC 37 although, apart from the recent work of Hugosson et al. 26 on the unrelaxed ͑001͒ surface of a wide series of transition-metal carbides, there is no experimental or theoretical information for the ͑001͒ surface of other transition-metal carbide surfaces.
In the present work we describe the crystal and electronic structures of a series of TMC having a 1:1 stoichiometry. In order to investigate trends in a systematic way, three different groups of transition-metal elements have been chosen. These are group 4 ͑Ti, Zr, Hf͒, group 5 ͑V, Nb, Ta͒, and group 6 ͑Cr, Mo, W͒. For most of the carbides under consideration ͑TiC, VC, ZrC, NbC, HfC, and TaC͒, the most stable phase corresponds to a fcc crystal packing; for MoC various stable phases exist and the ␦-MoC, although not being the most stable phase, has been chosen because it displays also a fcc structure. The tungsten carbide ͑WC͒ has been included in spite of exhibiting a distorted hexagonal close packed structure. Finally, chromium carbide has been excluded because it does not exhibit any phase with CrC stoichiometry. In a first step, bulk properties such as lattice parameter, bulk modulus, and electronic structure are calculated, whereas in a second step the electronic and geometric structure of the ͑001͒ surface is investigated in detail including surface formation energies and work functions. In particular, the influence of surface relaxation on these properties is stressed.

II. METHOD OF CALCULATION
The calculations were carried out in the framework of DFT within the GGA. All electron calculations were carried out with the DMol 3 code, 38 while frozen core calculations were performed with a parallel version of VASP4.5 code. [39][40][41] Two different implementations of the GGA were used: the exchange-correlation functional proposed by Perdew et al., 42 hereafter referred to as PW91, and the revised version of the Perdew-Burke-Ernzerhof functional, 43 usually referred to as RPBE. 44 In the PW91 calculations with VASP, a plane-wave basis set was used to span the valence electronic states and the core electrons were represented by the projected augmented wave ͑PAW͒ method of Blöchl. 45 This can be regarded as an all-electron frozen core method, which combines the accuracy of an all electron description with the computational simplicity of the pseudopotential approach, especially in the implementation of Kresse and Joubert. 46 This representation of the core states allows one to obtain converged results with a cutoff kinetic energy of 415 eV for the plane-wave basis set. A Monkhorst-Pack grid has been used to select the special k points necessary to carry out numerical integrations in the reciprocal space. The number of k points in each direction has been adapted for each system, until total energies computed with increasing sets of k points differ by less than 0.01 eV. The meshes used vary from 11ϫ 11ϫ 11 to 17ϫ 17ϫ 17 for bulk carbides and from 9 ϫ 9 ϫ 1 to 13ϫ 13ϫ 1 for the different ͑001͒ surfaces. The choice of odd numbers in the k-point grids is just for computational convenience. In this way, the high symmetry, points coincide with some of the k points on the Monkhorst-Pack grid. No point group symmetry constraints are imposed in the calculations.
A second set of periodic calculations has been carried out for selected systems using the RPBE functionals, which, in principle, is thought to provide more accurate adsorption energies. The RPBE calculations have been carried out using a numerical basis of double-quality plus polarization functions ͑DND in our notation͒ reported by Delley. 47,48 These calculations have been carried out using the DMol 3 computational code and a Brillouin unit cell of 21 k points along the most symmetric directions. However, even if the RPBE functional seems to provide adsorption energies closer to experiment than those obtained by the PW91 functional, 44 it has also been pointed out that it tends to lead to worse results for bulk properties such as lattice parameters and bulk moduli. 49 Moreover, there are examples where the RPBE functional may even overcorrect adsorption energies, thus predicting incorrect adsorption energies. 50 Therefore, a comparison of results obtained using both PW91 and RPBE seems very convenient.
For the bulk, structural optimization was performed using a conjugated gradient technique in which the iterative relaxation of atomic positions was stopped when the change in the total energy between successive steps was less than 0.001 eV. With this criterion, forces on the atoms are generally smaller than 0.1 eV/ Å. Once the lattice parameter was determined for each carbide, a model for the ͑001͒ surface has been constructed using a slab model approach. In this approach one uses a unit cell which is repeated periodically in two dimensions while it has a finite extent in the third one. In order to use the plane-wave basis set which is periodic in nature the corresponding slab is repeated in the third direction with the slabs separated by a sufficiently wide vacuum region. The vacuum width was progressively increased until energy variations were not significant. A vacuum width of 10 Å proved to be good enough for the present purposes. Slab models having four ͑Slab-4͒ or eight ͑Slab-8͒ atomic layers were tested and the three or four outermost atomic layers, respectively, were fully relaxed using the optimization procedure outlined above. Hence, the remaining atomic layers-one for Slab-4 and four for Slab-8-were fixed at the bulk geometry. The Slab-4 still exhibits some minor edge effects in the third direction, so that most of the discussion will focus on the results obtained with Slab-8. Nevertheless, the Slab-4 has also been used for comparison purposes, mainly to investigate the differences between the two different computational methods.

III. GEOMETRY AND ELECTRONIC STRUCTURE OF BULK METAL CARBIDES
The fcc crystal packing of most of the metal carbides studied in the present work-TiC, ZrC, HfC, VC, NbC, TaC, and ␦-MoC-results from an ABC layer stacking. This work focus on the ͑001͒ surface, which has the same amount of metal and carbon atoms. However, WC presents a distorted hcp packing since the WC structure is made of alternating ͑001͒ layers of metal and carbon atoms.
In a first step, the geometry optimization of each transition-metal carbide has been carried out; both the coordinates of the atoms in the unit cell and the unit cell parameter were allowed to relax. From the final geometric data the bulk modulus B 0 was computed as in Eq. ͑1͒, making use of one of the features of VASP which, for a given unit cell volume, provides the equivalent external pressure from the corresponding analytical energy gradient. The lattice parameter a 0 and B 0 values are reported in Table I and compared with available experimental data. Experimental cell parameters have been taken from the inorganic crystal structure database ͑ICSD͒. 51 In all cases, we found that the calculated a 0 values, both PW91 and RPBE, are in excellent agreement with experiment, the largest deviations being always less than 1% except for the ␦-MoC where a somehow larger deviation is found. Interestingly enough, both functionals provide similar results with small, albeit erratic, deviations from experiment.
Due to the uncontrolled amount of vacancies in these refractory materials, experimental B 0 values are scattered in a rather broad range, 52,53 and may also be affected by the temperature. This does not allow a direct comparison with the calculated data ͑Table I͒ although it appears that all cal-TABLE I. Calculated and experimental values for the lattice parameter ͑a 0 and a 0 expt.͒, calculated and range of experimental values of the bulk modulus ͑B 0 and B 0 expt., respectively͒ and the relative errors ͑%a 0 ͒ given in percent. Distances are in angstrom and bulk modulus in GPa. Notice that except for WC the metal-carbon distance is just a 0 / 2 whereas for the former is 1.34a 0 .  54 In order to investigate the electronic structure of the metal carbides studied in the present work we make use of the corresponding DOS, scaled to the Fermi level. The DOS of the different carbides exhibit similar features, the corresponding plot for TiC is given in Fig. 1. The projection of the DOS on the atomic components permits to decompose the main contributions to the chemical bonding. Here it is important to remark that because of the use of a plane-wave basis set the projection into local atomic states is not straightforward. In the present work the projected DOS is obtained using harmonic spherical functions in a given atomic radius and hence, although qualitatively correct, it is somehow arbitrary. Here, the radii of the atomic spheres used to carry out the integration of the projected DOS were the VASP standard internal values. An alternative representation can be obtained by means of the numerical basis calculations although this does also have a certain degree of arbitrariness. However, the DOS is mainly used for interpretation and it is unlikely that the qualitative description will change.
For TiC, a sharp peak at ϳ−33.5 eV corresponding to the Ti͑3p͒ orbitals appears, which is not shown in the corresponding figures. Clearly, these orbitals do not contribute to the bonding. The next peak appears at ϳ−9.5 eV, it is rather sharp and mainly of C͑2s͒ character although with some Ti͑3d͒ character; this peak is sufficiently far from the Fermi level so that it can also be considered as essentially nonbonding. Just below the Fermi level a very broad structure with two well defined main peaks shows up corresponding to a mixing of Ti͑3d͒ and C͑2p͒ orbitals. These peaks correspond to the bonding states, and their antibonding counterparts appear just above the Fermi level. These bonding and antibonding peaks exhibit a clear mixing of metal 3d and carbon 2p states although the latter has a stronger Ti͑3d͒ character.
Thus, the present results are in agreement with the previous DFT results based on a LMTO basis set. 17,18 The DOS for the other metal carbides studied in the present work are very similar; almost in all the cases, there are significant amounts of states on the Fermi level in agreement with the conducting nature of these materials. The main difference in the metalcarbide DOS appears in the relative position of the Fermi level between the bonding and antibonding states and three different patterns can be distinguished. In group 4 TMCs ͑TiC, ZrC, and HfC͒ and also in WC, the Fermi level is placed almost at the middle of the gap between the bonding and antibonding state ͑Fig. 1͒. Note that the behavior of WC is due to its different crystal structure which results in a DOS which is reminiscent of that of the group 4 carbides. In VC and TaC the extra electron shifts upwards the Fermi level making them better conducting materials ͑Fig. 2͒. The same happens in NbC but in this case the metal 5s band is well within the Fermi level; this is also the case for ␦-MoC ͑Fig.
3͒. This feature is not observed in the FP-LAPW of Amriou et al. since only projection onto the Nb͑4d͒ and C͑2s ,2p͒ is reported. Nevertheless, these authors claim that the Nb͑5s͒ band is vey high in energy due to the repulsion with the C͑2s͒ state. However, one must realize that the corresponding states are partially hybridized and not atomically pure. Moreover, this feature is as expected from the Nb 5s-shell partial occupancy of atomic Mo and Nb.
The analysis of the DOS provides useful information about the covalent contributions to the chemical bond but does not permit one to evaluate the importance of ionic contributions. A possibility consists in using some kind of population analysis but this is rather difficult when using a planewave basis set and always arbitrary when using a localized basis set due to the subjective partition of the overlap populations. A possible, although quite qualitative, alternative is to investigate the relative position of the center of a particular band, deep enough to be considered as nearly atomic, and its shift with respect to a given reference. Here, it is necessary to point out that energy level shifts are caused by a Although the latter is in many cases the dominant effect, the remaining contributions cannot be neglected. Hence, it is very difficult to extract charges from core-level shifts. Nevertheless, comparison of core-level shifts along a series of similar compounds in similar conditions can provide a reasonable qualitative trend for charge redistribution. To this end, we find it appropriate to compare the center of the C͑2s͒ band of the different carbides relative to the Fermi level with the corresponding value for bulk diamond which can be considered as a prototype of a pure covalent solid ͑Fig. 4͒ and for bulk graphite which exhibits a metallic like, delocalized electron density. There are various important features which emerge from the comparison of the C͑2s͒ levels. First, the C͑2s͒ peak for the carbides is much narrower than that of diamond and graphite indicating that the electrons are more localized. Second, in the carbides the C͑2s͒ is shifted towards lower binding energies, this suggests a more nega-tively charged C atoms. These two features are consistent with a non-negligible ionic contribution to the chemical bond. Note that TiC, ZrC, and HfC exhibit the largest shift with respect to diamond. According to the interpretation above, the ionicity of these materials decreases along the series in agreement with charge-transfer arguments estimated from XPS and NEXAFS core-level shifts. Interestingly enough, this interpretation is also supported by the net charges on the metal atoms estimated from the Mulliken population analysis carried out using the all electron DND localized basis set and the RPBE method ͑Table II͒.
The next important point involves the cohesive energies of the different metal carbides. The cohesive energy E coh is defined as the difference between the energy of the condensed system and that of the atoms in the vacuum. In this case there are only two atoms per unit cell ͑one C and one metal͒ and therefore E coh is simply the difference between the energy per unit cell and the energy of the carbon and of the metal atom; note that within this definition E coh is negative. The summary of E coh values is reported in Table II together with experimental values and those previously reported from LMTO/LDA calculations. 19,21 The present GGA values, either PW91 or RPBE, are in good agreement with experiment with deviations which, in general, do not exceed 0.5 eV. The improvement on the E coh values over previous LMTO/LDA calculations-with deviations with respect to experiment that in some cases are larger than 3 eV-is really clear. Nevertheless, one has to be aware that the atomic energies have been calculated placing the atom in a box surrounded by a sufficiently large vacuum width. Hence, the total angular space and spin moments are not well defined, the multiplet structure is lost and the only information about the atomic state is the electronic configuration. Previous work has shown that the GGA predicted ground state electronic configuration of the different transition-metal atoms is in agreement with experiment except for Ni and V, for further details see Refs. 55 and 56. The E coh values, either calculated or experimental, increase along the group but do not show any particular trend along the series. In order to try to understand the origin of the different stability and its rela-tionship with the electronic structure various correlations between E coh and electronic structure parameters have been investigated. No relation between the cohesive energy and either the peak positions, the cell parameter, or the centre of the metal d band has been found. However, for fcc TMC, it is found that E coh varies almost linearly with the energy separation between the centers of the bonding and antibonding peaks in the DOS ͑Fig. 5͒. The center of each peak is defined as the position at which the integral of the DOS is half that of the total value for the peak. Therefore, the relative stability of the metal carbides is dominated by the bondingantibonding splitting. This would seem to indicate that covalent effects play a dominant role. However, one must realize that the bonding-antibonding splitting is not symmetric and hence accounts as well for ionic contributions. Nevertheless, the linear behavior in Fig. 5 would seem to indicate that in spite of the fact that the chemical bond in these carbides has metallic, ionic, and covalent contributions, the latter grossly explains the relative order of stabilities. Moreover, one cannot claim that it is dominated by the covalent contribution since the C͑2s͒ energy level shifts indicate the presence of non-negligible charge transfer. This is in agreement with a  significant charge on the metal atom as predicted from the RPBE calculations carried out in a local basis set although the difference between compounds is rather small ͑see Table  II͒.

IV. GEOMETRY AND ELECTRONIC STRUCTURE OF THE METAL CARBIDES "001… SURFACE
Two different slab models have been used to simulate the ͑001͒ surface of the metal carbides studied in the present work, these are the Slab-4 and Slab-8 models described in Sec. II. Both models exhibit a surface made of mixed squares of metal and carbon atoms, and hence the surface composition respects the bulk stoichiometry. In the Slab-4 model the three outermost atomic layers have been allowed to relax in the ͓001͔ direction while keeping the fourth one fixed to mimic the bulk. In the Slab-8 models the same strategy is used but now four atomic layers are relaxed and four are fixed. The interlayer distances for the relaxed structures are given in Table III for the Slab-4 model. Since the metal and carbon atoms of a given atomic layer can relax in a slightly different way-rumpling-we use n and nЈ to define the carbon and metal atoms of the n atomic layer. Therefore, d nn Ј corresponds to the rumpling of the n atomic layer whereas d n Ј ͑n+1͒ corresponds to the minimum distance between atoms in different layers ͑see Fig. 6͒. For the Slab-4 model, the GGA calculated values obtained using either the all electron or the PAW implementations, and using the RPBE or PW91 exchange-correlation potentials, respectively, are in very good agreement except for ␦-MoC and TaC where the latter approach predicts a slightly larger rumpling effect in the sec- FIG. 5. ͑Color online͒ Calculated cohesive energy vs the bondingantibonding splitting ͑see E coh and ⌬ in Table III͒.   TABLE III. Intralayer and interlayer atomic perpendicular distances ͑in angstrom͒ obtained with the Slab-4 model and at the PW91 and RPBE levels of theory. n and nЈ define the carbon and metal atoms of the n atomic layer, d nn Ј corresponds to the rumpling of the n atomic layer, whereas d n Ј ͑n+1͒ corresponds to the minimum distance between atoms in different layers. The percent rumpling is also indicated for comparison. ond atomic layer. In any case, it is worth to point out that only the outermost atomic layers exhibit a significant rumpling. Hence, values for the rumpling of the innermost atomic layers are not reported in Table III. The results discussed below evidence that surface relaxation is important and should not be neglected. Nevertheless, values for the unrelaxed surfaces, even if somehow unphysical, are reported precisely to highlight the importance of surface relaxation effects. The agreement between the two different approaches for the Slab-4 model permits to use the computationally more efficient plane-wave basis with the PAW description of the core electrons to study the influence of the slab thickness on the surface final geometry. The results for the rumpling effect are presented in Table IV using the same definitions as for  Table III. The relaxation and rumpling values for the first layer are very similar for the Slab-4 and Slab-8 models. In contrast, the relaxation and rumpling values for the second and third atomic layers predicted by the Slab-8 models are noticeably smaller than the values predicted by the Slab-4 model. This result is important and implies that predictions from rather thin metal-carbide layers have to be handled with caution. In fact, the origin of these differences is the lack of enough bulk material below the relaxed atomic layers. Re-sults from Table IV show that the rumpling effect grows from group 4 to group 6 whereas along the groups no clear trends exist. Hence, the extent of intralayer relaxation seems to be driven by the number of metal d electrons. Unfortunately, comparison to experiment is only possible for TiC͑001͒ and TaC͑001͒. From a qualitative point of view the present calculations are in agreement with experimental evidence showing an outward relaxation of C atoms and an inward relaxation of the metal atoms. [32][33][34]37 Moreover, the present results for TiC͑001͒ and TaC͑001͒ are in excellent agreement with the recent quantitative low-energy electron diffraction ͑LEED͒ analysis of Tagawa et al. 33 for TiC͑001͒ ͑experimental value of 6.5% compared to the present calculated value of 5.1%͒ and Gruzalski et al. 37 for TaC͑001͒ ͑experimental value of 9% compared to the present calculated value of 8.9%͒ and represent a significant improvement over the LDA calculations of Kobayashi. 57 This author studied TiC͑001͒, TaC͑001͒, HfC͑001͒, NbC͑001͒, and ZrC͑001͒ and found that these surfaces exhibit a significant rippling although the calculated values appear to be 2% larger than the present ones.
The surface energy formation per surface atom E sur along the series has been calculated as in Eqs. ͑2͒ and ͑3͒ for the unrelaxed and relaxed surfaces, respectively, FIG. 6. ͑Color online͒ Schematic representation of the interlayer and intralayer ͑rumpling͒ atomic relaxation of the metal carbides ͑001͒ surface. TABLE IV. Intralayer and interlayer atomic perpendicular distances ͑in angstrom͒ obtained with the Slab-8 model at the PW91 level of theory. n and nЈ define the carbon and metal atoms of the n atomic layer, d nn Ј corresponds to the rumpling of the n atomic layer, whereas d n Ј ͑n+1͒ corresponds to the minimum distance between atoms in different layers. The percent rumpling is also indicated in parenthesis for comparison.
In these equations, n at-sup stands for the number of surface atoms in the slab unit cell; 4 in the present calculations. In the unrelaxed case, the factor 2 in the denominator comes from the fact that two surfaces are formed. The surface energy for the relaxed surface is obtained from the relaxed slab calculation, where one surface is kept fixed, subtracting the surface energy for the unrelaxed surface. This procedure ensures that the number of atoms in the bulk and in the slab is the same but it depends moderately on the number of atomic layers on the slab model. Table V reports the E sur values for the Slab-4 and Slab-8 models at the bulk and relaxed geometries, the small differences between the two models are indicative that in this case the results for the thicker model can be considered as converged. In principle this is a small difference although the relative error is larger than 50%. These differences seem to have two different origins: on one hand the use of the atomic sphere approximation ͑ASA͒ in the LMTO calculations, and on the other hand the use of the LDA exchange-correlation potential. For some carbides the use of the ASA introduced differences of ϳ0.15 eV in the surface energy as compared with the full potential LMTO results. 26 The remaining difference can be ascribed to the well-known tendency of LDA to overbind; this will also result in too large surface energies. Next, we comment on the work function values calculated for the Slab-8 model at the PW91 level ͑Table VI͒. The work function is computed as = V − E f , where V is the electrostatic potential in the vacuum and E f the corresponding Fermi level energy. Overall, the calculated values are within 0.5 eV from experiment which represents an average error of 10%-15%. It is important to remark that ignoring surface relaxation results in additional errors of about the same order of magnitude. In general, the inclusion of surface relaxation effects results in larger values of the calculated work function. The largest error is found for TaC and HfC which correspond to the 5d transition metals. It is likely that the error arises from the relativistic effects not completely included in the PAW potential. Overall, the present results represent a considerable improvement over those reported recently using the LMTO/LDA approach. However, while in most cases the present results have a lower error, it is also clear that the LMTO/LDA performs better for TaC and HfC. Since the present results are obtained using a more flexible representation of the valence and core densities and also by considering a complete geometry optimization for the surface, it is likely that the better performance of LMTO/LDA in these two systems is due to a cancellation error. To conclude this section we analyze the DOS of the Slab-8 model and compare to the corresponding bulk values. The most salient general features are a broadening of some peaks ͑the metal outermost valence p level, for example͒, a destabilization of the bonding peaks, and a stabilization of the antibonding ones. The first one is caused by symmetry breaking with respect to that of the bulk due to the appearance of the surface. The surface metal atoms loose the cubic symmetry and, hence, the otherwise threefold degenerate levels split into two components ͑a 1 -and e-like͒. The differential shift of bonding and antibonding levels contributes to decrease the gap with respect to that of the bulk. This effect is somehow enhanced when surface relaxation is properly taken into account. The reason behind this further effect is the relatively complex relaxation involving rippling of the atomic layers ͑Fig. 6͒. Each C atom moves upwards and the M atom moves downwards so that for each M atom there are four larger and one shorter M -C distances than in the bulk, the larger distances destabilizes the corresponding orbitals and since there is a 4:1 ratio, the final effect is a further destabilization of bonding contributions with a concomitant stabilization of the antibonding part.

V. CONCLUSIONS
In this work, the bulk and surface geometrical and electronic properties of a series of transition-metal carbides have been studied by DFT methods. For the bulk systems, the lattice parameters and bulk modulus are in good agreement with experimental data. The chemical bonding has been studied by analyzing the DOS plots and the shifts of the C͑2s͒ peak with respect to diamond. It is concluded that in these materials the chemical bonding is strongly covalent, the cohesive energies being directly related to the bondingantibonding bands. Nevertheless, the C͑2s͒ shifts with respect to diamond indicated that some metal to carbon charge transfer occurs. Group 4 metal carbides can be viewed as small band-gap semiconductors, with a weak overlap between the valence band and the conduction bands. Beyond group 4, metal carbides show a metallic character due to a partially filled d band, in agreement with the thermal and electric conductivity of these materials. However, in some cases ͑NbC and ␦-Moc͒, the metallic character is due to a partially occupied metal s band.
The ͑001͒ face of these metal carbides exhibits a noticeable surface rumpling which grows along the series. The extent of the surface rumpling predicted by the two levels of theoretical calculations ͑all electron with the RPBE functional and the PW91 one with the PAW description of the atomic cores͒ for the first layer are in good agreement with each other as well as with available experimental values. Nevertheless, accurate values for the second and third layers require including a few atomic layers to represent the bulk. Surface energies and work functions have also been predicted and compared to available theoretical data and experiment. One of the important conclusions is that neglecting surface relaxation results in very large errors on the surface energy and work function although in the latter the calculated values are slightly larger than experiment and further increase upon surface relaxation. Finally, the formation of a surface results in a noticeable shift of the most prominent features of the total DOS with respect to the corresponding values in the bulk. The extent and nature of the shift can be understood from simple bonding-antibonding arguments and is enhanced by the structural rippling of this surface.