Interaction potentials from periodic density-functional theory calculations: Molecular-dynamics simulations of Au clusters deposited on the TiN (001) surface

J. Chem. Phys. 123, 244706 (2005); https://doi.org/10.1063/1.2140702 123, 244706 © 2005 American Institute of Physics. Interaction potentials from periodic densityfunctional theory calculations: Moleculardynamics simulations of Au clusters deposited on the TiN (001) surface Cite as: J. Chem. Phys. 123, 244706 (2005); https://doi.org/10.1063/1.2140702 Submitted: 21 June 2005 . Accepted: 01 November 2005 . Published Online: 28 December 2005 N. Cruz Hernández, and Javier Fdez. Sanz ARTICLES YOU MAY BE INTERESTED IN A climbing image nudged elastic band method for finding saddle points and minimum energy paths The Journal of Chemical Physics 113, 9901 (2000); https://doi.org/10.1063/1.1329672 TiC and TiN supported platinum monolayer as high-performance catalysts for CO oxidation: A DFT study The Journal of Chemical Physics 149, 054705 (2018); https://doi.org/10.1063/1.5038857 Interaction potentials from periodic density-functional theory calculations: Molecular-dynamics simulations of Au clusters deposited on the TiN „001…surface N. Cruz Hernández and Javier Fdez. Sanza兲 Departamento de Química Física, Facultad de Química, Universidad de Sevilla, E-41012 Sevilla, Spain 共Received 21 June 2005; accepted 1 November 2005; published online 28 December 2005兲 Molecular-dynamics simulations of gold particles deposited on a TiN 共001兲surface have been accounted for through classical pair potentials describing the atom force field. The interaction between Ti–N, Ti–Ti, N–N, Au–Au, Au–Ti, and Au–N pairs was estimated by following a procedure in which the interaction energy between two sets of atoms is estimated from density-functional calculations performed with periodic boundary conditions using plane waves as basis set. The pair potentials were expressed as the sum of two contributions: long range in a Coulomb form and a short-range term, which included the rest of the energy contributions. Simulations of the TiN 共001兲 isolated surface reproduced the already described surface relaxation, with a rippling parameter in agreement with that found from a purely first-principles approach. Simulations of gold deposition on such surfaces showed the formation of metal clusters with well-defined fcc structure and epitaxially grown. © 2005 American Institute of Physics.关DOI: 10.1063/1.2140702兴 I. INTRODUCTION Density-functional theory,1using plane-wave functions and periodic boundary conditions 共p-DFT兲, constitutes one of the most powerful tools nowadays in computer material science 共see Ref. 2 for a review兲. Within the p-DFT framework, static properties of medium-sized systems 共about 100 atoms兲are currently being simulated. However, in spite of the brilliant elegance of this methodology, its application to large systems is restricted by the high computational cost. Therefore, the study of problems with relatively large models must be undertaken using other kinds of methods. An additional aspect to consider concerns the study of the dynamic behavior of a system. For instance, ab initio moleculardynamics 共aiMD兲simulations3based on the DFT calculations constitute an extraordinarily powerful tool in this field. Recently, computational power has increased to the point where the aiMD methods based on the p-DFT calculations can be used to describe many processes of chemical interest and may be a suitable option if the interest is in systems of relatively small dimensions. With this method it is possible to take into account the changes of the electronic structure of the system during the simulation. However, morphology of three-dimensional metal particles on the surface can be obtained only if we use models of larger dimensions. Classical molecular-dynamics 共cMD兲simulations4con- stitute one of the most effective methodological alternatives for the study of problems that require models with large numbers of particles. With this approach it is possible to use models with thousands of particles. However, in the cMD methods the correct election of the force field that defines the interaction between atoms has long been an important bottleneck in the theoretical modeling of many types of physical systems. Intermolecular and interatomic interaction potentials 共force fields兲are used in a vast area of modern science. In the simplest cases, i.e., for interactions between atoms and small linear molecules, the potentials are usually computed on a first-principles basis or fitted to experimental data. Assuming a global functional for the system, a full potential-energy surface 共PES兲can be obtained by fitting the data or by interpolating in some manner between the known data points. Therefore, once a set of potential energies is computed on a grid of points, an analytic function has to be fitted to these data in order to be used in various types of simulations. Many forms of such functions have been proposed in literature but a consensus as to which of them is most effective is yet to be determined. Along with this, the fact that currently available accurate analytical fits of finite discrete sets of ab initio potential-energy data typically involve large numbers of adjustable parameters. To avoid any requirement of chemical intuition as well as the possibility of unphysical features being introduced, an alternative procedure makes use of interpolation methods that do not need a mathematical expression to represent the interaction between atoms. In this case the potential surface is just a set of data on a dense grid and the actual values of the potential are estimated from, for instance, spline interpolation. One of the most common ways to represent the interaction in a group of atoms is by using pairwise interaction potentials. According to this, the sum of all pair interactions results in a global PES formulation. Such interaction potentials may take varied forms and usefully represented as a sum of a long-range Coulomb interaction and short-range interaction. In most cases the short-range contribution includes parameters that are fit in order to agree with the experimental a兲Author to whom correspondence should be addressed. Electronic mail: [email protected] THE JOURNAL OF CHEMICAL PHYSICS 123, 244706 共2005兲 0021-9606/2005/123共24兲/244706/6/$22.50 © 2005 American Institute of Physics123, 244706-1 results or ab initio calculations. However, independent of the expression that has been used to describe the short-range interaction potential, a more practical procedure to perform the simulations is to build up a table of values of the potential and forces with respect to the square distance. Then, the contributions to the energy and forces are obtained by numerical interpolation. This alternative reduces the computational cost of the simulations but can also be used to set up interaction potential that does not match a given analytical expression. The interaction potentials used in this work follow this approach: instead of looking for an expression to analytically fit the potential, the interaction between two species is estimated in a range 共dense兲of interatomic distances and directly used in the simulations. With this aim, the short-range potential contribution is estimated in this work from periodic DFT calculations using a recently outlined5mathematical procedure that will be fully described in Sec. II. In this procedure the interaction energy between a given set of particles is represented by means of a collection of vectors and coefficients obtained by solving a system of linear equations built up using the DFT energies obtained in a close dense range of interatomic distances. A simplified variant to this method that used a diagonal approach to the system of equations was also employed to simulate the behavior of Pd64 clusters deposited on a clean ␣ -Al2O3共0001兲surface.6Both the morphology of the transition-metal particle and the stability with respect to the temperature were found to agree with the experimental data. The procedure to obtain interaction pair potentials, Sec. III, will be subsequently used to analyze the structural properties of Au clusters deposited on titanium nitride 共001兲sur- face by means of the MD simulations. Titanium nitride, TiN, is a material of high interesting technological and catalytic applications because of its exceptional physical and chemical properties, which make it attractive both from fundamental and experimental points of view.7,8 Transition-metal nitrides are commonly referred to as refractory metals as they have extreme hardness and a high melting point, properties comparable to those of covalent crystals, as well as good corrosion resistance and metallic conductivity close to the reported values of pure transition metals. TiN has also attracted interest as a support for the deposition of transition-metal clusters, both for engineering purposes 共sensors, electronic devices, etc.,兲and for catalysis.9,10 In spite of this, only a few theoretical works11–14 have been devoted to study its surface properties 共for a review on TiN and other nitrides see Ref. 15兲. In particular, we have recently reported an investigation on the relaxation of the 共001兲surface for several transitionmetal nitrides as well as the interaction of Au atoms with those surfaces using p-DFT calculations.16 In the present paper, as a part of our continuous effort to unravel the properties of such interfaces, we report on the classical MD simulations of the TiN 共001兲surface and gold clusters deposited on it 共Sec. IV兲. II. THEORETICAL ASPECTS Let us assume two groups of atoms 共A and B兲inside a super cell. They are placed so that repeating the cell in the three directions two films of A and B separated by d共Å兲are created 共Fig. 1兲. In addition, the model has an enough empty space on top of group A in order to minimize the interaction energy with other films obtained by repeating the cell in z direction. The idea of the method is to generate a set of pair potentials between atoms of groups A and B, so that when calculating the interaction energy between A and B films at distance d共using cMD interaction potentials兲,EMD int共A,B,d兲,we obtain the same result as using the p-DFT calculations Ep-DFT int共A,B,d兲, Ep-DFT int共A,B,d兲⬇EMD int共A,B,d兲.共1兲 The interaction energy Ep-DFT int共A,B,d兲is computed by Ep-DFT int共A,B,d兲=Ep-DFT A,B,d−Ep-DFT A−Ep-DFT B,共2兲 where the first term of the right member is the energy of the cell including A and B atoms at distance d, while the next two terms correspond to the energies of isolated A and B groups. When the distance between the films is large enough 共d⬁兲the interaction between the two groups is mainly of long-range type. On the other hand, the right member of 共1兲can be written in a similar way to Eq. 共2兲, EMD int共A,B,d兲=EMD A,B,d−EMD A−EMD B.共3兲 FIG. 1. Scheme of the supercell containing two groups of atoms 共A and B兲 used to calculate the interaction energy from the DFT calculations. The vacuum space over group A is large enough to neglect the interaction with outer cell film. 244706-2 N. Cruz Hernández and J. Fdez. Sanz J. Chem. Phys. 123, 244706 共2005兲 Taking into account that the interaction potential for MD is a sum of long-range and short-range terms, the interaction energy is EMD int共A,B,d兲=共EMD long共A,B,d兲−EMD long共A兲−EMD long共B兲兲+EMD ⬘short共A,B,d兲, 共4兲 where EMD ⬘short共A,B,d兲=EMD short共A,B,d兲−EMD short共A兲−EMD short共B兲共5兲 is the contribution to the interaction energy only due to the short-range potentials between A and B atoms. Increasing the value of dcauses EMD ⬘short共A,B,d兲to become zero rapidly for any rational set of potentials chosen. The distance for which EMD ⬘short共A,B,d兲⬇0 determines the value of the parameter d⬁, mentioned above. The actual value of d⬁ depends on the configuration adopted for atoms in both groups and on the cutoff distance 共rc兲proposed for the shortrange potentials. At the same time, rcis indirectly related to the form of the long-range potential, i.e., the charges assigned to the atoms. The key to this procedure is the correct selection of such charges in the groups 共q兵A其,q兵B其兲 and the cutoff radius rcfor each pair of atoms. In principle, expression 共1兲is always fulfilled for any set of charges, q兵A其,q兵B其, and rc, but not always in all cases is a reasonable short distance d⬁between the groups found so that the short-range interaction vanishes. Replacing Eqs. 共4兲and 共2兲in 共1兲we obtain Ep-DFT A,B,d−Ep-DFT A−Ep-DFT B−共EMD long共A,B,d兲−EMD long共A兲−EMD long共B兲兲 =EMD ⬘short共A,B,d兲共6兲 and a similar equation for d⬁, Ep-DFT A,B,d⬁−Ep-DFT A−Ep-DFT B−共EMD long共A,B,d⬁兲−EMD long共A兲−EMD long共B兲兲 =EMD ⬘short共A,B,d⬁兲⬇0. 共7兲 With the aim of avoiding the calculation of Ep-DFT A, Ep-DFT B,EMD long共A兲, and EMD long共B兲, Eq. 共7兲is subtracted from the set of equations 共6兲. Following this, we also assure that the set of the pairwise potentials vanishes at rc. Finally, we obtain the set of equations 共as many as points dichosen兲: Ep-DFT A,B,d−Ep-DFT A,B,d⬁−共EMD long共A,B,d兲−EMD long共A,B,d⬁兲兲=EMD ⬘short共A,B,d兲. 共8兲 In Eq. 共8兲the terms related to the long-range interaction are easily computed using the Ewald17 sum procedure once the atomic charges of the A and B groups have been assigned. Let us assume that groups A and B are composed of different species of atoms. Then there are as many interaction potentials to obtain as the number of pair species jk that can be formed. We will call these the pair potential jk or Pjk. In addition to the cutoff radius for a given jk pair, rc jk, let us define a minimum radius 共rmin兲, so that two species jk cannot be closer than the distance rmin jk for any type of calculation 共p-DFT or cMD兲. In principle, our task is to obtain a set of functions Pjk共r兲 with rmin jk ⬍r⬍rc jk. Because Pjk共r兲are continuous, we do not need necessarily to know the mathematical form of Pjk共r兲, but the values of these functions at some points ri jk, where rmin jk ⬍ri jk⬍rc jk, i.e., to change to the discrete space on r. Therefore, our task is to obtain i=1,2,…,Njk components of the vector Pjk for each jk pair of species. The right-hand side of Eq. 共8兲can then be written as a linear combination of Pjk i as EMD ⬘short共A,B,di兲=兺 j苸A k苸B 兺 t=1 t=Njk ␣ i,jktPjk t.共9兲 To solve this system we need as many equations 共calculated configurations兲as unknown short-range vector components, N共with N=兺j苸A k苸BNjk兲, and, in addition, the equations must be linearly independent. After solving the system of equations 共9兲, the short-range contribution to the pairwise interaction potential is obtained. In Eq. 共9兲, the coefficients ␣ i,jkt represent the number of times that the term Pjk tappears in the right-hand side of Eq. 共8兲. They are calculated by forming all the possible pairs of atoms between A and B groups, fulfilling the condition rmin jk ⬍r⬍rc jk, taking into account not only the computational cell but also the neighbor cells. III. DETERMINATION OF INTERACTION POTENTIALS The first step of the simulations was to determine the force field governing the system behavior that, according to our scheme, and considering first the support, TiN, involves determination of three parameters PTiTi,PTiN, and PNN for a set values of rin the rmin,rcrange 共here set to 1 and 5 Å, respectively兲. With this aim we used an initial supercell in which groups A and B 共two subcells兲are placed at a distance d, which fulfills the rmin jk ⬍r⬍rc jk condition. Because of periodic boundary conditions, groups A and B actually represent two slabs of TiN, and there is a vacuum in the supercell to prevent intercell interaction between slabs. As shown in Fig. 1, groups A and B consisted of two Ti atoms and two N atoms arranged in slabs two-layer thick. In order to get a dense mesh the distance interval should be divided in, say, 1000 steps, which obviously would involve quite cumbersome computations. Instead of that we preferred to compute the energy at slab distances dgoing from 1 to 5 Å, step of 0.1 Å 共41 geometries兲, and to perform then an interpolation. The energy of these 41 structures was estimated from the DFT calculations using the VASP 4.5 共Refs. 18 and 19兲code and the generalized gradient approach 共GGA兲.20 Full details on these calculations can be found in Ref. 13. After cubic spline interpolation a set of 1000 linearly independent equations of the form 共9兲was obtained. Since we have only a third of the equations required, we repeat the process for two more initial different configurations, obtaining a final set of 3000 equations that after solving lead to the 1000 components of vectors PTiTi,PTiN, and PNN.5 To simulate the Au/TiN interface we now need to estimate the interaction concerning Au atoms: Au–Au, Au–Ti, and Au–N 共assuming transferability of the PTiTi,PTiN, and PNN vectors兲. As already reported, isolated gold atoms prefer to bind the TiN 共001兲surface on top of Ti atoms, with a small charge transfer towards the surface. However, when 244706-3 Simulation of Au clusters deposited on TiN 共001兲surface J. Chem. Phys. 123, 244706 共2005兲 Au2dimers are supported, no such a transfer was observed. Similar behavior was also observed in the Pd/ ␣ -Al2O3inter- face where adsorption of the metal cluster involved only a small polarization.6,21 These findings enabled us to limit the interaction potential of Au to the short-range contribution since it is obvious that in the interface higher coverage does occur. In a similar way to that exposed for TiN, three initial configurations were chosen to determine the pair interaction potentials Au–Ti, Au–N, and Au–Au 共Fig. 2兲. Here group A consists of two Au atoms in all cases and group B was either a slab of TiN 共001兲five-layer thick or two gold atoms deposited on such a slab. Gold atoms were placed on top of Ti and on top of N 关Figs. 2共a兲and 2共b兲兴. In Fig. 2共c兲two gold atoms in group A interact with two gold atoms of group B, allowing to determine an interaction potential of the gold-adsorbed gold type. Following these configurations, the complete set of EMD ⬘short共A,B,d兲values were obtained by variation of the distance dso that 1.0 Å⬍d⬍5.0 Å. The total pair interaction potential curves are shown in Fig. 3. IV. MD SIMULATIONS OF Au CLUSTERS DEPOSITED ON TiN „001…SURFACE The MD simulations were undertaken in the microcanonical ensemble using the DLគPOLY computer code.22 The computational box consisted of a cube 34.032 Å side with a TiN slab six-layer thick and an enough vacuum gap in the 具001典direction. In the simulations, the two lowest layers were kept frozen and 1024 atoms were moved simultaneously. The initial geometry of particles in the slab was set coherent with the cell parameter obtained from the DFT calculations: a=4.254 Å 关experimental value 4.238 Å 共Ref. 23兲兴. On this surface, gold clusters of different sizes 共200– 350 atoms兲were deposited. The simulation proceeds by first assigning initial velocities to all atoms according to a Maxwell-Boltzmann distribution that depends on the temperature of the systems 共300 K兲. From this starting point, Newton’s equations of motion were solved by numeric integration using the leap-frog algorithm with a time step of 0.1 fs.24 The simulation schedule consisted of a thermalization run of 10 ps followed by a simulation of 10 ps without temperature control to ensure that the system reached equilibrium. After this stage, a simulation of 15 ps was carried out to collect data for statistics. The conservation of energy was better than 1 in 10 000. The starting step of the simulations was to investigate the behavior of the TiN 共001兲isolated surface. The surface shows to be stable at that temperature, and the major features can be deduced from the radial distribution function 共RDF兲 reproduced in Fig. 4. For the Ti–N pair the first band is large, with a peak at 2.1 Å, in agreement with the Ti–N interatomic distance of 2.127 Å obtained from the DFT calculations. According to the fcc structure of TiN, this peak is followed by another at about 共a/2兲冑3 Å due to the second-neighbor shell. FIG. 2. The three initial configurations used to determine the pair potentials describing the interaction of Au with the surface. FIG. 3. Interaction potentials for Au–Au, Au–Ti, and Au–N pairs. FIG. 4. Radial distribution functions 共RDFs兲for Ti–N, N–N, and Ti–Ti pairs obtained from the MD simulations of TiN 共001兲surface. 244706-4 N. Cruz Hernández and J. Fdez. Sanz J. Chem. Phys. 123, 244706 共2005兲 For Ti–Ti and N–N pairs, the most outstanding features are broad peaks slightly falling below 3 Å, i.e., somewhat lower than the theoretical distance of 共a/2兲冑2 Å, followed by bands about 4 Å, i.e., a共Å兲. Compared to the RDF obtained from bulk simulations,5two main differences may be observed. Primarily, the first band of Ti–Ti and N–N pairs appear to be slightly shifted to shorter interatomic distances, suggesting that there is some contraction of the surface. Secondly, for these pairs, a shoulder in the second band is clearly observed. These shoulders are a consequence of the rippling relaxation of the surface, in which the outermost layers of nitrogen and titanium atoms are not in the same plane. The difference between the positions of these shoulders gives a rippling of 0.12 Å, in quite good agreement with that found from p-DFT calculations 共0.17 Å兲.13 This result constitutes a test of the effectiveness of the pairwise interaction potentials obtained in the preceding sections. In order to simulate the Au/TiN interface two gold clusters consisting of 216 and 343 atoms were added to the surface. 共Since the results are practically identical, only the 343 atoms deposit will be reported.兲A snapshot of the system after the simulations is reproduced in Fig. 5. The gold atoms closest to the surface are found on top of Ti atoms, in agreement with the largest Au–Ti interaction. In the second layer, gold atoms are placed on the vertical above the outermost nitrogen atoms, bound to four first-layer gold atoms. The next layer is identical to the first one, according to an almost layer-by-layer ABABAB pattern. In fact, the structure of the Au atoms is like that of the titanium lattice in the nitride, corresponding to a compact packing fcc in agreement with the natural form that is determined experimentally for the bulk. This 共001兲TiN 储 共001兲Au epitaxial relation is similar to that observed for gold clusters deposited on the MgO 共001兲 surface from both high-resolution electron microscopy experiments and MD simulations.25–28 On the other hand, it is worth to note that there is a misfit between the Au and TiN lattices, being the Au cell parameter somewhat lower 共a =4.078 and 4.254 Å, respectively兲. This misfit mainly affects the first layer of gold atoms, which is found to be slightly stressed, however, beyond this layer; the Au cluster gradually contracts to recover the metallic structure. Such an accommodation of the gold lattice has already been observed in the deposition of Au clusters on MgO 共001兲surfaces.28 The analysis of the Au cluster shows a morphology corresponding to an irregular polyhedron truncated along specific directions, exposing 共111兲and 共001兲facets. The misfit between the TiN and Au lattices is also responsible for the wideness of the first band observed in the Au–Au RDF reported in Fig. 6. The main peak approximately appears at 2.9 Å, which is close to the experimental distance of 2.88 Å.29 There is also another effect due to the undercoordination of Au atoms lying on the outermost layers of the cluster. These atoms tend to shorten the bond distances and give rise to the shoulder observed at lower distances. With respect to the other two pairs, the first peak of Au–Ti RDF appears at 2.7 Å and corresponds to Au–Ti pairs of atoms directly bonded, i.e., the first layer of gold and outermost layer of the support. The second band, much more intense, also arises from these two layers, but involve Au–Ti atoms not directly bound 共one Au by each four Ti and vice versa兲. Finally, the most prominent feature of the Au–N RDF is the peak at 3.4 Å originated from the first-layer Au atoms and outermost N atoms. V. SUMMARY AND CONCLUSIONS In this paper we report on a procedure that allows us to determine interaction pairwise potentials specifically tailored for the system under study. These ad hoc potentials are expressed as a sum of long- and short-range contributions and obtained, in our case, from the DFT calculations carried out under periodic boundary conditions. After solving a set of equations linearly independent, the potentials are obtained and directly used without further fitting to any analytical expression. Although the procedure is in principle of general purpose, our intention is to use it to simulate the surface properties of materials, and we report on an application of FIG. 5. Snapshot of the Au343 cluster deposited on the TiN 共001兲surface. FIG. 6. Radial distribution functions 共RDFs兲for Au–Au, Au–Ti, and Au–N pairs obtained from the MD simulations of Au cluster deposited on the TiN 共001兲surface. 244706-5 Simulation of Au clusters deposited on TiN 共001兲surface J. Chem. Phys. 123, 244706 共2005兲 the complete procedure to TiN. The MD simulations of TiN 共001兲surface shows noticeable rumpling of titanium and nitrogen outermost layers, with a rippling relaxation parameter estimated from the radial distribution functions of 0.12 Å. This value is in reasonable agreement with that obtained using a full DFT approach, and supports the reliability of the potentials. Simulations of a gold particle deposited on the TiN 共001兲 surface shows that the transition-metal particle epitaxially grows on the surface, in such a way that the Au atoms are on top of surface Ti atoms. The gold particle shows fcc structure and mainly exposes the 共111兲and 共001兲facets. Because of the misfit of the Au and TiN lattices, the layer of gold atoms closest to the substrate is dilated, and gradually contracts towards the external region of the cluster to recover the metallic interatomic distance. ACKNOWLEDGMENTS The project was financially supported by the Spanish Ministerio de Ciencia y Tecnologia—Project No. MAT2005- 01872 and the Ramon y Cajal program 共awarded to N.C.H.兲—and the Junta de Andalucía 共Project No. FQM- 132兲. 1R. G. Parr and W. Yang, Density-Functional Theory of Atoms and Molecules 共Oxford University Press, New York, 1989兲. 2M. C. Payne, M. P. Teter, D. C. Allan, T. A. Arias, and J. D. Joannopoulos, Rev. Mod. Phys. 64,1045共1992兲. 3R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 共1985兲. 4M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids 共Claredon, Oxford 1987兲. 5N. C. Hernández and J. F. Sanz, Comput. Mater. Sci. 共in press兲. 6N. C. Hernández and J. F. Sanz, J. Phys. Chem. B 105, 12111 共2001兲. 7M. Wittmer, J. Vac. Sci. Technol. A 3, 1797 共1985兲. 8L. E. Toth, Transition Metal Carbides and Nitrides 共Academic, New York, 1971兲. 9S. Tsubota, M. Haruta, T. Kobayashi, A. Ueda, Y. Nakahara, and G. Poncelet Preparation of Catalysts V 共Elsevier, New York, 1991兲,p.695. 10 M. A. Centeno and M. J. Lloret, Mater. Sci. Forum 383, 111 共2002兲. 11 M. Marlo and V. Milman, Phys. Rev. B 62, 2899 共2000兲. 12 J.-Q. Li, Y.-F. Zhang, S.-C. Xiang, and Y.-N. Chiu, J. Mol. Struct.: THEOCHEM 530, 209 共2000兲. 13 N. C. Hernández, J. Graciani, and J. F. Sanz, Surf. Sci. 541,217共2003兲. 14 C. V. Ciobanu, D. T. Tambe, and V. B. Shenoy, Surf. Sci. 582,145 共2005兲. 15 L. I. Johansson, Surf. Sci. Rep. 21, 177 共1995兲. 16 N. C. Hernández and J. F. Sanz, Mater. Sci. Forum 383,177共2002兲. 17 P. Ewald, Ann. Phys. 共N.Y.兲64, 253 共1921兲. 18 G. Kresse and J. Hafner, Phys. Rev. B 47, R558 共1993兲. 19 G. Kresse and J. Furthmuller, Comput. Mater. Sci. 6,15共1996兲. 20 J. Perdew, J. Chevary, S. Vosko, K. Jackson, M. Pederson, D. Singh, and D. Fiolhais, Phys. Rev. B 46, 6671 共1992兲. 21 R. B. Gomes, F. Illas, N. C. Hernández, A. Márquez, and J. F. Sanz, Phys. Rev. B 65, 125414 共2002兲. 22 W. Smith and T. R. Forester, J. Mol. Graphics 14, 136 共1996兲. 23 N. Schoenberg, Acta Chem. Scand. 8, 213 共1954兲. 24 L. Verlet, Phys. Rev. 159,98共1967兲. 25 S. Giorgio, C. Chapon, C. R. Henry, G. Nihoul, and J. M. Penisson, Philos. Mag. A 64,87共1991兲. 26 T. Kizuka, T. Kachi, and N. Tanaka, Z. Phys. D: At., Mol. Clusters 26, S58 共1993兲. 27 T. Kizuka and N. Tanaka, Phys. Rev. B 56, R10079 共1997兲. 28 B. Pauwels, G. Van Tendeloo, W. Bouwen, L. Theil Kuhn, P. Lievens, H. Lei, and M. Hou, Phys. Rev. B 62, 10383 共2000兲. 29 A. Maeland and T. B. Flanagan, Can. J. Phys. 42, 2364 共1964兲. 244706-6 N. Cruz Hernández and J. Fdez. Sanz J. Chem. Phys. 123, 244706 共2005兲