Surface Studies With the Embedded Atom Method

Published in Interface. Vol. 2, No. 1, Feb 1994

Dr. Daryush Ila, Alabama A&M University

Introduction. The embedded atom method[1-3] (EAM) has proven useful for modeling a variety of phenomena occurring at the surfaces of solids. Typically, the EAM is used to calculate the energy E for each of a large number of possible configurations of atoms. Since each calculation of E involves many atoms, practical applications of the method usually require large amounts of supercomputer time. In a practical application, one might test many atomic configurations to determine which has the lowest energy and is thus the most stable, or one might try to establish a "downhill" energy path in configuration space as the basis for a molecular dynamics simulation.

The EAM calculates E for a given configuration of atoms as

Fi (ri), the energy of embedding atom i into the crystal, depends only on the type of atom i and on

the superposition of atomic electron densities of all other atoms evaluated at atom i. The electrostatic interaction energy of atoms i and j, øij, depends only on their separation rij. Cutting off the atomic interaction beyond some maximum separation can serve to reduce the number of computations.

Hydrogen Adsorption. One EAM study was motivated by the observation, in low-energy electron-diffraction experiments on palladium (Pd), of diffraction spots associated with the adsorption of hydrogen (H). [4] The surface layer is thought to have the configuration designated c2x2 and shown in Figure 1. This surface layer corresponds to the 100 plane of a face-centered cubic (FCC) Pd crystal with H atoms at octahedral sites. Nearest neighbor repulsive and next nearest neighbor attractive H-H interactions of the indirect H-Pd-H "through-bonding" type suffice to order the overlayer.

An EAM-based simulation of the Pd-H system was run. Semi-empirical fittings to the bulk properties of Pd and to the heat of solution of H were used to obtain the required embedding and interaction energies and electron densities. At each step of the simulation, the code selects an atom at random, perturbs its position, calculates the new EAM energy, and uses a Metropolis algorithm [5] to decide whether to keep or reject the new configuration. About 10 [4] steps per mobile atom are required for meaningful results. Our model consists of 128 mobile Pd atoms and between 30 and 80 H atoms (all mobile). Thus, ~ 2x10[6] simulation steps are required for each run.

Every 200 steps, the code evaluates a lattice structure factor, a type of Fourier transform, to determine when a phase boundary (with respect to temperature) is crossed. These simulations agree with experiments in that when 63 H atoms are used, the c2x2 pattern is stable up to a higher temperature (Tc) relative to the other H:Pd ratios. However, the absolute Tc values for all ratios are lower than those observed experimentally.

Surface Melting. Another EAM study sought to verify the observation of surface melting up to 150K below the bulk melting temperature of lead (Pb). [6] In this case, the EAM energy is differentiated with respect to position for each atom to determine its acceleration, and a true time-dependent simulation of atomic motion is obtained. The model consists of 20(110) layers of 70 Pb atoms each. Five layers are pinned to their lattice points by enforcing zero velocity throughout the simulation. The positions of the atoms in the top 15 layers evolve in accord with the dynamic equations of motion integrated by the Nordsiek algorithm.[7] Again, the structure factor is used as an indicator of the transition between order and disorder. Figure 2 shows the decrease of the time-averaged square modulus of the structure factor (a decrease in order) as the temperature increases. Temperature control is achieved by subjecting the atoms to simulated collisions with the atoms of a surrounding heat bath.[8] The results of the simulation are consistent with experimental observations of a surface roughening transition near 400 K and more extensive surface disordering at about 550K. No estimate of the bulk melting temperature for Pb was obtained; the simulated bulk lattice layers remained solid at 650K, whereas the experimentally determined bulk melting point is 600K. Work in progress will calculate layer diffusivities to verify that disordering above 550K is surface melting and will examine more closely atomic configurations near 400K to verify surface roughing.

Conclusions. We see that the EAM was somewhat more successful in predicting the absolute temperature dependence of surface properties when used as the basis of a detailed simulation of atomic motion, as in the Pb study, rather than as a means of calculating the Metropolis "cost" function, as in the Pd-H study. Because the Metropolis algorithm is employed first as a minimization and then as a sampling scheme rather than as a means of following atomic motions in detail, H jumps between sites do not necessarily take kinetically realistic paths. Kinetically realistic jumps between the octahedral sites of FCC lattices require passage through an intermediate tetrahedral site to avoid cutting across a close packed direction, i.e., passing between Pd atoms at their distance of closest approach. By allowing direct jumps between octahedral sites, the simulation loses kinetic information but effects a much more rapid equilibration of the model and allows sampling of many more configurations.

Acknowledgments. We wish to thank Sandia National Laboratory for providing the DYNAMO molecular dynamics code and the Alabama Supercomputer Authority for Cray time.

References

1. M.S. Daw, M.I. Baskes, Phys. Rev. B 29, 6443

(1984).

2. S.M. Foiles, M.I. Baskes, M.S. Daw, Phys. Rev. Lett.

50,1285(1983).

3. D.J. Oh, R.A. Johnson, J. Mater, Res. 3,471 (1988).

4. P. Tibbits, M. Karimi, D. Ila, I. Dalins, and G.

Vidali, Mat. Res. Soc. Symp. Proc. 193,289 (1990).

5. N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth,

A.H. Teller, E. Teller, J. of Chem. Phys. 21,1087

(1953).

6. P. Tibbits, M. Karimi, D. Ila, I. Dalins, and G.

Vidali, J. Vac. Sci. Technol. A9(3),1937 (1991).

7. A. Nordseick, Math Comp. 16, 22 (1962).

8. H.C. Andersen, J. Chem. Phys. 72, 2384 (1980).