- Email: [email protected]
A novel method for simulating atoms and molecules interacting with ionic surfaces
Surface Science North-Holland
286 (1993) 182-189
‘surface science
A novel method for simulating atoms and molecules interacting with ionic surfaces Zhi-Hong
Huang
and
Hua
Guo
Department of Chemistry, University of Toledo, Toledo, OH 43606, USA Received
13 October
1992; accepted
for publication
29 December
1992
A new method for computer simulations of atoms and molecules interacting with a moving ionic surface is introduced. This approach is capable of evaluating the long-range Coulomb interaction energy for an ion on the surface and at any position inside a crystal. The main idea is to break down the total Coulomb energy to contributions from different layers in the crystal. The interaction energy of an ion with a layer to which it does not belong can be calculated with a two-dimensional Fourier method. For the interaction of the ion with other ions in the same layer, a quasi-two-dimensional Ewald sum is used to determine the interaction energy. The method is shown to be computationally efficient in evaluating the long-range Coulomb interaction energy of a system consisting of moving surface ions.
1. Introduction Interaction of atoms or molecules with solid surfaces plays a central role in many surface processes, such as molecular beam scattering from a surface [ll, photodissociation of adsorbed species [2] and chemical reactions on surfaces [3]. The dynamics of adspecies depends heavily upon the topographic features of the interaction potential. Hence, an accurate description of the interaction potential is especially important for computer simulations, either quantum or classical, of surface processes. Despite of expeditious accumulation of experimental data and rapid advances in theories, our knowledge of the gas/surface interaction is still far from complete. Such a lack of understanding is largely due to the vast number of particles involved in a gas/surface complex. Information on the gas/surface interaction may be obtained experimentally. For example, the adsorption energy of an adspecies can be acquired from temperature programmed desorption (TPD) [4] and information concerning the orientation of an adsorbate can be extracted from 0039-6028/93/$06.00
0 1993 - Elsevier
Science
Publishers
electron-stimulated desorption ion angular distribution (ESDIAD) [5]. However, a complete experimental determination of the potential surface is almost impossible because of the formidable dimensionality of the problem. Empirical, semiempirical and ab initio methods have been suggested to calculate the interaction potential for gas/surface complexes and many of them yield very satisfactory results [6-81. Ionic surfaces present a unique example for gas/surface interaction. LiF and MgO surfaces have been used in experiments to study scattering and other surface processes [l]. Most ionic solids are chemically inert and interact with adspecies through electrostatic forces which can be calculated quite accurately even without the knowledge of electron charge distribution. One added advantage of these solids for experimentalists is the fact that highly perfect crystal faces can be prepared with relative ease. If the surface motion is not important to the outcome of the chemical process, a rigid surface model would be sufficient. Under such circumstances, the interaction between an adsorbed molecule and the stationary surface can be characterized fairly well by a two-
B.V. All rights
reserved
Z.-H. Huang, H. Guo / Simulation of gas /surface
dimensional Fourier transform method [6]. In this approach, the total electrostatic potential at the surface due to all the positive and negative ions in the solid can be expanded into a Fourier series due to the periodicity of the crystal surface. This method has been used in many surface simulations [9] and proven to be very efficient because the Fourier expansion quickly converges in the reciprocal space. However, surface motion is important for many processes that involve large energy transfer between the surface and molecule [ll. The inclusion of surface motion, even in a finite region, would break the periodicity of the surface and therefore preclude the application of the Fourier method. The attractive force between oppositely charged ions in an ionic solid is essentially the Coulomb interaction, while a short-range repulsive potential due to Pauli repulsion prevents the crystal from collapsing. Both the long- and shortrange interactions can be adequately represented by pairwise potentials. One simplistic approach is to model the surface as a cluster [lO,lll. Such a model allows a direct pairwise evaluation of the gas/surface and surface/surface interactions and has been used to simulate photodissociation dynamics of adsorbed molecules [ll-151. However, one needs to be cautious in using such a model to represent an infinite solid in which long-range interaction plays a dominant role. Because of the long-range nature of the Coulomb forces between ions in the solid, a large number of surface atoms may have to be included in the modeling of the surface. Even that cannot guarantee an accurate description of the potential because it is well known that the sum of the long-range Coulomb interaction is only conditionally convergent and in principle should include an infinite number of ions 1161. For crystals, a properly chosen cluster model may be acceptable because of the regularity in structure, which may result in cancellations of long-range forces. As shown below, the approach is a reasonable approximation in the evaluation the Coulomb interaction, though errors may arise due to the truncation of the solid. Another possible alternative is to employ the three-dimensional Ewald sum [17] to a stack of slabs by constructing a supercell with part of
interaction
183
space filled and the rest empty. However, a very large unit cell may be required in order to eliminate the unphysical interaction between two consecutive slabs. Consequently, one may have an inefficient computational procedure. In this paper, we introduce a new method for calculating the Coulomb sum at ionic surfaces, which is comparable with the cluster model in using computing resources and yet includes an infinite number of ions in the substrate. The basic idea is to separate the crystal into two regions, a “primary zone” consisting of atoms which are close to the adsorbed atom/molecule and the other region including all the other atoms in the substrate at their equilibrium positions. The electrostatic potential outside the solid or in the primary zone is evaluated by taking the sum over all the ions at their equilibrium positions in the solid and then adding the difference due to the movements of ions in the primary zone. As shown below, the most difficult part is the calculation of the interaction energy between a substrate ion and the other ions in the same layer. We, in this paper, derive a quasi-two-dimensional Ewald sum to compute this interaction energy. The paper is organized as follows: in section 2, the model for an atom or molecule interacting with an ionic surface and the formulae for the evaluation of the Coulomb interaction are presented. Comparisons of this method to the cluster model and the two-dimensional Fourier transform method are made in section 3. A summary is given in section 4.
2. Model 2.1. General considerations
The major objective of this work is to develop an accurate yet efficient computational method to handle the molecular dynamics description of chemical reactions occurring at ionic surfaces. The simple model adopted in this work consists of a molecule (or an atom) sitting on an ionic surface. Since we are interested in molecule/surface energy exchange, we define a three-dimensional region in the substrate, consisting of a
184
Z.-H. Huang, H. Guo / Simulation of gas /surface
number of movable ions (e.g., 6 x 6 x 3). The adsorbed molecule will be placed on top of the primary zone. All the substrate ions outside the primary zone are assumed to be stationary and at their ideal lattice sites. The surface ions are treated as point charges and the pairwise interaction between two ions contains the long-range Coulomb interaction and short-range repulsion, for which only the nearest and the second nearest neighbor interactions are considered. It has been shown that such simple pairwise functions give satisfactory lattice structures [11,131. On the border between the primary zone and rest of the crystal, the movable ions inside are effectively connected with the static ions outside. An improvement of the current model will be to treat the border ions with the generalized Langevin equations (GLE) [lo], which will allow energy exchange between the primary zone and the rest of the solid. The total energy of the system is the sum of the intramolecular, molecule-substrate, and substrate-substrate interaction potential energies. If there is no chemical bonding between the molecule and the surface, the gas phase intramolecular potential should still be a good representation of the adsorbate intramolecular interaction. All the interactions except electrostatic terms are short-range and can be evaluated readily with a direct sum. In what follows, we concentrate on the long-range electrostatic interactions among surface ions and adsorbate. To calculate the electrostatic part of the molecule-surface potential, one may model the molecule as a collection of point charges, dipoles, and even quadrupoles. The interaction energy of those charges and multipoles, with an ionic surface at a point r above the primary zone, is determined by the Coulomb potential generated by all the ions in the solid. For a surface ion far away from the adsorbate, the vibrational motion of the ion near its equilibrium position makes little difference to the overall Coulomb potential near the adsorbate. Therefore, they can be treated as if they were fixed at their equilibrium positions. The vibrational motion of the ions inside the primary zone cannot be ignored and has to be taken into account. The overall Coulomb poten-
interaction
tial at r = (x, y, 2) can hence be written follows: Q(r)
as
= Qo(r) + A@(r) =
where qi, ri, and roi are the charge, actual and ideal positions of the ith ion in the crystal, respectively. In eq. (11, we have identified two types of ions: the ions in the primary zone (j E P> and the rest which are fixed at their ideal lattice positions. Thus, the electrostatic field generated by the semi-infinite crystal can be approximated as the sum of the contributions from all the ions at their ideal positions, which is represented by the first sum Q&r), and the correction due to the movement of the ions in the primary zone given by the second sum A@(r). The second sum can be carried out directly because only a finite number of ions are involved. We devote the next subsection to the calculation of the first sum, i.e., the long-range Coulomb potential at a point r, from the contribution of the substrate ions in their ideal lattice positions. 2.2. Coulomb potential of a semi-infinite crystal In this work, we present only the results for simple sodium chloride type surface (facecentered cubic lattice). Extension to other types of lattice structure should be straightforward. Assuming that the ions in the unit cell have the following arrangement: +q located at (0, 0, 01, (a/2, a/2, 01, (a/2, 0, a/2), and (0, a/2, a/2), and -4 at (a/2, 0, 01, (0, a/2, 01, (0, 0, a/2), and (a/2, a/2, a/2). The potential at an arbitrary point r = (x, y, z) due to ions in the substrate is
@o(r)=9
c c ,Na+:._r,
jEunit
cell
n
I
i
1
INa+r,-r+dl
I’
(2)
Z.-H. Huang, H. Guo / Simulation of gas /surface
where N= nxex + nYeY+ n,e,, n, and nY are integers from --m to +w, excluding (O,O), II, from - ~0to 0; e,, ey , and e, are the unit vectors along the X-, y- and z-axis, respectively, and d = (a/2, 0, 0) is the vector from the origin to the site of a negative ion. Note that all the ions take their ideal positions in eq. (2). Let us first discuss the Coulomb field above the surface for the adsorbate. Methods exist for summing all the electrostatic contributions from a perfect semi-infinite crystal [6]. The sum can be carried out in a layer-by-layer fashion, taking advantage of the periodicity in the X- and y-coordinates. The Coulomb potential generated by one layer can be readily expressed in terms of a two-dimensional Fourier series and the overall potential at the surface is therefore an infinite sum of all the layers in the crystal. It can be shown I63 that the lowest order approximation to the overall potential from all the layers below r is given as
interaction
185
sional Fourier transform method. It is easy to show [6] that the electrostatic potential due to a perfect two-dimensional mesh of ions can be written as follows, provided that z i:s of same order of magnitude as, or larger than, a,
Iz-zlI/a) xcos(24 where zl is “ + ” sign is positive and layers below
cos(27r$i,
the position of the f-layer and the taken if the ion at (X = 0, y = 0) is “-” otherwise. The sum of all the z is therefore [6]
xcos(273 4273 xcos
(
2lry
a )’
where a is the lattice constant. Note that the field decreases exponentially with the z-coordinate. As shown by Steele [6], this approximation gives an error of approximately 1% if 2 is comparable to the lattice constant a, which is satisfied in most gas/surface processes. Let us now turn to the Coulomb potential inside the crystal. Since ions in the primary zone are allowed to move, one needs to know the potential field at the positions of those ions. Our strategy is similar to the above method, i.e., treating the overall potential as the sum of contributions from individual layers. The Coulomb field at any point inside the crystal can thus be broken into three terms: @ = @D&w+
@above
+
@Ewald
9
(4)
i.e., the sum of contributions from layers below and above, and ions in the same layer, if r happens to lie in one of the layers (Q,,,,,). The first two terms can be evaluated using the two-dimen-
(5)
(6)
where zI is the position of the layer just below z. Note that eq. (6) is basically the same as eq. (3) except for a possible change of sign. For Qjabove, there are a finite number of layers, each of which has a contribution given by eq. (5). Problems arise if the ion is too close to a layer or in the layer. Unfortunately, the Fourier transformation method, upon which the derivation of eq. (5) is based, fails when I z - zI I is small compared to a [6]. Alternative approaches have to be developed. In the following discussion, we derive a quasi-two-dimensional Ewald sum to handle the Coulomb field very close or in the plane of an infinite two-dimensional mesh. We first consider the sum over the positive ions. The sum of negative ions can be treated in the same manner. Using the definition of the complimentary error function [18], we have 12m -=& I
exp( -r*k*)
= &kGexp(
-r2k2)
dk
dk + ierfc(Gr).
(7)
186
Z.-H. Huang, H. Guo / Simulation of gas /surface
We define R =xe,
the following
quantities,
+ye,,
interaction
Eq. (11) becomes
(8)
r=R+Aze
Z)
(9)
n = n,e,
+ n,,eY,
( 10)
2.6 = 12
Ccig.R
where e,, eY, and e, are unit vectors along the three axes, and II, and IZ~ are integers. The sum over a layer of positive ions can be written as
c In-f-1l
erfc(GIna +c
n
exp(-Ina-R12k2-Az2k2)
dk
erfc(G)na-rJ) +c n
(m-r1
(11)
’
where G is the convergence factor which controls the convergence of the two sums. Let F(R,
k)=&Eexp(--lnali’k’). n
= cf< g
s,k)
= G/F(R,k)
e’g’R,
(13)
e-ig’R dR,
where A is the total surface area. Substituting (12) into eq. (14) and integrating, we have
f(s,k) = $exp
, i-& 1
which leads to the expression
F(R,k)
2vG = -pp
(17)
.
The integral in the first sum has an analytical form [18], which leads to the following expression, c 1 n Ina-rl =- ,:
c y[eK” g
erfc(&
+ GAzj
erfc(GIna-rI)
where g’s are the two-dimensional reciprocal tice vectors, i.e., g. nu = 27~ X integers. Fourier transform f(g,k) is given as f(g,k)
Ina-rl
(12)
By virtue of its periodicity in (x, y), F(R,k) can be expanded into a Fourier series in the reciprocal space: F( R,k)
n
-rl)
latThe
(14) eq.
(15)
for the F(R,k),
(16)
+c n
Ina-rI
’
( 18)
Eq. (18) is equivalent to the formula stated by Parry [19], who derived it as a limiting case of the 3D Ewald sum by assuming the interlayer spacing in the z-direction to be infinite. The sum in eq. (18) includes a term with g = 0, which is divergent. If the layer is electrically neutral such as a (001) layer of the NaCl crystal, this term can be excluded since the contributions from the positive and negative ions cancel each other. In cases where net charge is not zero, such as in the (111) layer of NaCl, the divergence can also be removed by including multiple layer contributions [20,21]. The g = 0 term has also been found to be important in non-centrosymmetric crystal structures [22]. The formula has been extended and expressed in terms of charge speading functions of arbitrary analytical forms [21]. To evaluate the potential energy as one pulls out an ion from its equilibrium position, one needs to subtract the field generated by the ion itself. Assuming the ion being selected is the one
Z-H. Huang, H. Guo / S~~~t~n
of gas /surface
187
interaction
at r0 = 0, one should replace the n = 0 term in the last sum by the following, lim
r+O
erfc( Gr) r
-2G/v?,
(19)
The intralayer interaction can then be obtained from eqs. (18) and (19). The method introduced above has been applied to the study of adsorption and photodissociation of CH,I adsorbed on a LiF(OO1) surface. We have implemented this method on a Cray YMP and our simulations showed that this method is numerically stable and reasonably efficient. The details of implementation of the method and simulations of photodissociation of methyl iodide on an LiF surface are published elsewhere [231.
3. Comparison
_,4
L-L-LL 0.0
0.5
/ ! .o
1
15
!
/
I
2.0
2.5
3.0
_L_
._..
3.5
i
4.0
zc&
Fig. I. Plot of electrostatic potential energy of a central surface ion as it is pulled away in the direction of surface normal. Other ions in the solid are assumed to be stationary and at their ideal lattice positions. The solid line corresponds to the quasi-two-dimensional Ewald sum method, while the dotted and the dashed lines are for the two-dimensional Fourier transform method and the cluster model, respectively. Variable .z is measured from the lattice positions of the ion.
with other methods
In this section, we will examine the difference between the quasi-ho-dimensional Ewald sum model and two other approximate methods. We choose a LiF(001) surface since we have used both the cluster model and the quasi-two-dimensional Ewald sum method to study photodissociation of several molecules on this surface. For simplicity, we assume that both Li and F carry a charge of “ +e” and “ -e”, respectively, though it may be more reasonable to use partial charges to represent the Li and F ions [24]. We first examine the Coulomb interaction energy q@(r) as one pulls out a central ion in the top layer along the surface normal. Fig. 1 is the plot of the energy calculated using three different methods, the solid line for quasi-two-dimensional Ewald sum, the dotted line for the two-dimensional Fourier method, and the dashed line for the cluster model. There are 6 x 6 X 3 ions in the cIuster model. As shown in fig. 1, our quasi-twodimensional Ewald sum method gives an excellent agreement with the cluster model at small z. In the caIculation of the Ewald sum, the number of g and R points are the same and set to be 36 while the convergence factor G is 0.214 A-‘. The electrostatic potentiat calculated using the
ho-dimensional Fourier transform method, as we pointed out earlier, diverges for small z. Fig. 1 shows that the Fourier method is acceptable only when z > 1 A. Nevertheless, it is a good representation of the field for interlayer interaction. It should be pointed out that eq. (18) is exact and the accuracy can be improved by including more g and n points. The results shown in fig. 1 justify our approach presented in section 1. I.e., for intralayer interaction, the quasi-do-dimensional Ewald sum method is a very good appro~mation, while for intertayer interaction, the Fourier method is quite satisfactory. According to eq. (3), the electrostatic potential on an ionic surface decays exponentially into the vacuum, while for a cluster model, the potential is expected to decrease only as a certain power of r. In fig. 2, we plot the potential energy of a test charge over a central positive ion as a function of z. Again, the cluster model agrees quite well with the Fourier method, with an error in the potential energy only in tens of meV. It is surprising to note the accuracy of the cluster model, even with a relatively small size of 6 X 6 X 3 atoms, is so accurate in the determination of the electrostatic potential above the sur-
Z.-H. Huang, H. Guo / Simulation of gas /surface
188
- - - Fourier
I) 05
1
0.m
i_
’
-n.lls 1
2
3
4
’ 5
6
7
’
’
R
9
1 10
z ca,
Fig. 2. The potential energy of a test charge as a function of z over a central positive ion. The solid and dashed lines are, respectively, for the cluster model and the Fourier transform method.
face. This is probably due to the perfectly regular arrangement of the ions in the cluster, which leads to a cancellation at the adsorption site in the center of the cluster surface. This conclusion offers some justification for the cluster models that have been adopted in many surface related simulations [ll-141. However, one should bear in mind that the cluster model is a finite-size model and caution must be exercised in order to avoid possible errors. In fig. 3, we replot the potential energy using the three methods for an ion at the corner of the cluster. It is clear that due to the finite size, the cluster model is no longer accurate, while the other two methods are not affected. It is expected that Q, calculated using the
interaction
cluster model is position dependent. One may attach the cluster to some fixed ions to avoid the problem, which provides not only a more accurate potential field above the surface but also the forces needed to stabilize the ions at the boundaries of the cluster. However, the number of fixed ions could be significantly larger than the moving ions. For a 6 X 6 X 3 cluster, if two fixed layers of ions are added, we have 232 fixed ions. Undoubtedly, the evaluation of potential energies and forces is now more costly since for each ion, an extra sum over the fixed ions needs to be performed. We have used both a cluster model and the quasi-two-dimensional method in computer simulations of photodissociation of small molecules adsorbed on a LiF(001) surface, and found that the efficiency is comparable for the two methods.
4. Summary In this paper, we have introduced a new method for treating the long-range Coulomb interaction in a system consisting of moving ionic surfaces. The calculation of the Coulomb interaction energy is broken down into sums over layers. For an ion in one layer interacting with another layer, a two-dimensional Fourier method is used to evaluate the energy. For an ion interacting with other ions in the same layer, a quasi-two-dimensional Ewald sum is performed. The method is computationally efficient and yet able to evaluate the Coulomb potential to arbitrary precision. It is intended to be implemented in Monte Carlo or molecular dynamics simulations.
Acknowledgement
This work was supported by the National Science Foundation (Grant No. CHE-9116501).
References [I] R.B. Gerber,
Fig. 3. Same as in fig. 1 except for a corner
ion.
Chem. Rev. 87 (1987) 29. [2] J.C. Polanyi and H. Rieley, Photochemistry in the Adsorbed State, Eds. C.T. Rettner and M.N.R. Ashfold (The Royal Society of Chemistry, Cambridge, 1991);
Z.-H. Huang, H. Guo / Simulation of gas /surface
[3] [4] [5] [6] [7] [8] [9] [lo] [Ill [12]
X.-L. Zhou, X.-Y. Zhu and J.M. White, Surf. Sci. Rep. 13 (1991) 73. G.A. Somorjai, Chemistry in Two-dimensions: Surfaces (Cornell University Press, Ithaca, NY, 1981). D.A. King, Surf. Sci. 47 (1975) 384. R.D. Ramsier and J.T. Yates, Jr., Surf. Sci. Rep. 12 (1991) 243. W.A. Steele, The Interaction of Gases with Solid Surfaces (Pergamon, New York, 1974). A.M. Stonehem and J.H. Harding, Annu. Rev. Phys. Chem. 37 (1986) 53. J.C. Tully, Annu. Rev. Phys. Chem. 31 (1980) 319. J.C. Polanyi, R.J. Williams and S.F. O’Shea, J. Chem. Phys. 94 (1991) 978. R.R. Lucchese and J.C. Tully, J. Chem. Phys. 80 (1984) 3451. MI. McCarthy and R.B. Gerber, J. Chem. Phys. 93 (19901 887. MI. McCarthy, R.B. Gerber and M. Shapiro, J. Chem. Phys. 92 (1990) 7780.
interaction
189
(131 Z.-H. Huang and H. Guo, J. Chem. Phys. 96 (1992) 8564. [14] Z.-H. Huang and H. Guo, J. Chem. Phys. 97 (1992) 2110. [15] J.M. Watson, I. NoorBatcha and R.R. Lucchese, J. Chem. Phys. 96 (1992) 7771. [16] Simulation of Liquids and Solids, Eds. G. Ciccotti, D. Frenkel and I.R. McDonald (North-Holland, Arnsterdam, 1987) ch. 5. 1171 P. Ewald, Ann. Phys. 64 (1921) 253. [18] Handbook of Mathematical Functions, Eds. M. Abramowitz and I.A. Stegun (Dover, New York, 1970) ch. 7. [19] D.E. Parry, Surf. Sci. 49 (1975) 433. [20] D.E. Parry, Surf. Sci. 54 (1976) 195. [21] D.M. Heyes and F. van Swol, J. Chem. Phys. 75 (1981) 5051. [22] M.W. Deem, J.M. Newsam and S.K. Sinha, J. Phys. Chem. 94 (1990) 8356. [23] Z.-H. Huang and H. Guo, J. Chem. Phys., in press. [24] G. Dolling, H.G. Smith, R.M. Nicklow, P.R. Vijayaraghavan and M.K. Wilkinson, Phys. Rev. 168 (1968) 970.
286 (1993) 182-189
‘surface science
A novel method for simulating atoms and molecules interacting with ionic surfaces Zhi-Hong
Huang
and
Hua
Guo
Department of Chemistry, University of Toledo, Toledo, OH 43606, USA Received
13 October
1992; accepted
for publication
29 December
1992
A new method for computer simulations of atoms and molecules interacting with a moving ionic surface is introduced. This approach is capable of evaluating the long-range Coulomb interaction energy for an ion on the surface and at any position inside a crystal. The main idea is to break down the total Coulomb energy to contributions from different layers in the crystal. The interaction energy of an ion with a layer to which it does not belong can be calculated with a two-dimensional Fourier method. For the interaction of the ion with other ions in the same layer, a quasi-two-dimensional Ewald sum is used to determine the interaction energy. The method is shown to be computationally efficient in evaluating the long-range Coulomb interaction energy of a system consisting of moving surface ions.
1. Introduction Interaction of atoms or molecules with solid surfaces plays a central role in many surface processes, such as molecular beam scattering from a surface [ll, photodissociation of adsorbed species [2] and chemical reactions on surfaces [3]. The dynamics of adspecies depends heavily upon the topographic features of the interaction potential. Hence, an accurate description of the interaction potential is especially important for computer simulations, either quantum or classical, of surface processes. Despite of expeditious accumulation of experimental data and rapid advances in theories, our knowledge of the gas/surface interaction is still far from complete. Such a lack of understanding is largely due to the vast number of particles involved in a gas/surface complex. Information on the gas/surface interaction may be obtained experimentally. For example, the adsorption energy of an adspecies can be acquired from temperature programmed desorption (TPD) [4] and information concerning the orientation of an adsorbate can be extracted from 0039-6028/93/$06.00
0 1993 - Elsevier
Science
Publishers
electron-stimulated desorption ion angular distribution (ESDIAD) [5]. However, a complete experimental determination of the potential surface is almost impossible because of the formidable dimensionality of the problem. Empirical, semiempirical and ab initio methods have been suggested to calculate the interaction potential for gas/surface complexes and many of them yield very satisfactory results [6-81. Ionic surfaces present a unique example for gas/surface interaction. LiF and MgO surfaces have been used in experiments to study scattering and other surface processes [l]. Most ionic solids are chemically inert and interact with adspecies through electrostatic forces which can be calculated quite accurately even without the knowledge of electron charge distribution. One added advantage of these solids for experimentalists is the fact that highly perfect crystal faces can be prepared with relative ease. If the surface motion is not important to the outcome of the chemical process, a rigid surface model would be sufficient. Under such circumstances, the interaction between an adsorbed molecule and the stationary surface can be characterized fairly well by a two-
B.V. All rights
reserved
Z.-H. Huang, H. Guo / Simulation of gas /surface
dimensional Fourier transform method [6]. In this approach, the total electrostatic potential at the surface due to all the positive and negative ions in the solid can be expanded into a Fourier series due to the periodicity of the crystal surface. This method has been used in many surface simulations [9] and proven to be very efficient because the Fourier expansion quickly converges in the reciprocal space. However, surface motion is important for many processes that involve large energy transfer between the surface and molecule [ll. The inclusion of surface motion, even in a finite region, would break the periodicity of the surface and therefore preclude the application of the Fourier method. The attractive force between oppositely charged ions in an ionic solid is essentially the Coulomb interaction, while a short-range repulsive potential due to Pauli repulsion prevents the crystal from collapsing. Both the long- and shortrange interactions can be adequately represented by pairwise potentials. One simplistic approach is to model the surface as a cluster [lO,lll. Such a model allows a direct pairwise evaluation of the gas/surface and surface/surface interactions and has been used to simulate photodissociation dynamics of adsorbed molecules [ll-151. However, one needs to be cautious in using such a model to represent an infinite solid in which long-range interaction plays a dominant role. Because of the long-range nature of the Coulomb forces between ions in the solid, a large number of surface atoms may have to be included in the modeling of the surface. Even that cannot guarantee an accurate description of the potential because it is well known that the sum of the long-range Coulomb interaction is only conditionally convergent and in principle should include an infinite number of ions 1161. For crystals, a properly chosen cluster model may be acceptable because of the regularity in structure, which may result in cancellations of long-range forces. As shown below, the approach is a reasonable approximation in the evaluation the Coulomb interaction, though errors may arise due to the truncation of the solid. Another possible alternative is to employ the three-dimensional Ewald sum [17] to a stack of slabs by constructing a supercell with part of
interaction
183
space filled and the rest empty. However, a very large unit cell may be required in order to eliminate the unphysical interaction between two consecutive slabs. Consequently, one may have an inefficient computational procedure. In this paper, we introduce a new method for calculating the Coulomb sum at ionic surfaces, which is comparable with the cluster model in using computing resources and yet includes an infinite number of ions in the substrate. The basic idea is to separate the crystal into two regions, a “primary zone” consisting of atoms which are close to the adsorbed atom/molecule and the other region including all the other atoms in the substrate at their equilibrium positions. The electrostatic potential outside the solid or in the primary zone is evaluated by taking the sum over all the ions at their equilibrium positions in the solid and then adding the difference due to the movements of ions in the primary zone. As shown below, the most difficult part is the calculation of the interaction energy between a substrate ion and the other ions in the same layer. We, in this paper, derive a quasi-two-dimensional Ewald sum to compute this interaction energy. The paper is organized as follows: in section 2, the model for an atom or molecule interacting with an ionic surface and the formulae for the evaluation of the Coulomb interaction are presented. Comparisons of this method to the cluster model and the two-dimensional Fourier transform method are made in section 3. A summary is given in section 4.
2. Model 2.1. General considerations
The major objective of this work is to develop an accurate yet efficient computational method to handle the molecular dynamics description of chemical reactions occurring at ionic surfaces. The simple model adopted in this work consists of a molecule (or an atom) sitting on an ionic surface. Since we are interested in molecule/surface energy exchange, we define a three-dimensional region in the substrate, consisting of a
184
Z.-H. Huang, H. Guo / Simulation of gas /surface
number of movable ions (e.g., 6 x 6 x 3). The adsorbed molecule will be placed on top of the primary zone. All the substrate ions outside the primary zone are assumed to be stationary and at their ideal lattice sites. The surface ions are treated as point charges and the pairwise interaction between two ions contains the long-range Coulomb interaction and short-range repulsion, for which only the nearest and the second nearest neighbor interactions are considered. It has been shown that such simple pairwise functions give satisfactory lattice structures [11,131. On the border between the primary zone and rest of the crystal, the movable ions inside are effectively connected with the static ions outside. An improvement of the current model will be to treat the border ions with the generalized Langevin equations (GLE) [lo], which will allow energy exchange between the primary zone and the rest of the solid. The total energy of the system is the sum of the intramolecular, molecule-substrate, and substrate-substrate interaction potential energies. If there is no chemical bonding between the molecule and the surface, the gas phase intramolecular potential should still be a good representation of the adsorbate intramolecular interaction. All the interactions except electrostatic terms are short-range and can be evaluated readily with a direct sum. In what follows, we concentrate on the long-range electrostatic interactions among surface ions and adsorbate. To calculate the electrostatic part of the molecule-surface potential, one may model the molecule as a collection of point charges, dipoles, and even quadrupoles. The interaction energy of those charges and multipoles, with an ionic surface at a point r above the primary zone, is determined by the Coulomb potential generated by all the ions in the solid. For a surface ion far away from the adsorbate, the vibrational motion of the ion near its equilibrium position makes little difference to the overall Coulomb potential near the adsorbate. Therefore, they can be treated as if they were fixed at their equilibrium positions. The vibrational motion of the ions inside the primary zone cannot be ignored and has to be taken into account. The overall Coulomb poten-
interaction
tial at r = (x, y, 2) can hence be written follows: Q(r)
as
= Qo(r) + A@(r) =
where qi, ri, and roi are the charge, actual and ideal positions of the ith ion in the crystal, respectively. In eq. (11, we have identified two types of ions: the ions in the primary zone (j E P> and the rest which are fixed at their ideal lattice positions. Thus, the electrostatic field generated by the semi-infinite crystal can be approximated as the sum of the contributions from all the ions at their ideal positions, which is represented by the first sum Q&r), and the correction due to the movement of the ions in the primary zone given by the second sum A@(r). The second sum can be carried out directly because only a finite number of ions are involved. We devote the next subsection to the calculation of the first sum, i.e., the long-range Coulomb potential at a point r, from the contribution of the substrate ions in their ideal lattice positions. 2.2. Coulomb potential of a semi-infinite crystal In this work, we present only the results for simple sodium chloride type surface (facecentered cubic lattice). Extension to other types of lattice structure should be straightforward. Assuming that the ions in the unit cell have the following arrangement: +q located at (0, 0, 01, (a/2, a/2, 01, (a/2, 0, a/2), and (0, a/2, a/2), and -4 at (a/2, 0, 01, (0, a/2, 01, (0, 0, a/2), and (a/2, a/2, a/2). The potential at an arbitrary point r = (x, y, z) due to ions in the substrate is
@o(r)=9
c c ,Na+:._r,
jEunit
cell
n
I
i
1
INa+r,-r+dl
I’
(2)
Z.-H. Huang, H. Guo / Simulation of gas /surface
where N= nxex + nYeY+ n,e,, n, and nY are integers from --m to +w, excluding (O,O), II, from - ~0to 0; e,, ey , and e, are the unit vectors along the X-, y- and z-axis, respectively, and d = (a/2, 0, 0) is the vector from the origin to the site of a negative ion. Note that all the ions take their ideal positions in eq. (2). Let us first discuss the Coulomb field above the surface for the adsorbate. Methods exist for summing all the electrostatic contributions from a perfect semi-infinite crystal [6]. The sum can be carried out in a layer-by-layer fashion, taking advantage of the periodicity in the X- and y-coordinates. The Coulomb potential generated by one layer can be readily expressed in terms of a two-dimensional Fourier series and the overall potential at the surface is therefore an infinite sum of all the layers in the crystal. It can be shown I63 that the lowest order approximation to the overall potential from all the layers below r is given as
interaction
185
sional Fourier transform method. It is easy to show [6] that the electrostatic potential due to a perfect two-dimensional mesh of ions can be written as follows, provided that z i:s of same order of magnitude as, or larger than, a,
Iz-zlI/a) xcos(24 where zl is “ + ” sign is positive and layers below
cos(27r$i,
the position of the f-layer and the taken if the ion at (X = 0, y = 0) is “-” otherwise. The sum of all the z is therefore [6]
xcos(273 4273 xcos
(
2lry
a )’
where a is the lattice constant. Note that the field decreases exponentially with the z-coordinate. As shown by Steele [6], this approximation gives an error of approximately 1% if 2 is comparable to the lattice constant a, which is satisfied in most gas/surface processes. Let us now turn to the Coulomb potential inside the crystal. Since ions in the primary zone are allowed to move, one needs to know the potential field at the positions of those ions. Our strategy is similar to the above method, i.e., treating the overall potential as the sum of contributions from individual layers. The Coulomb field at any point inside the crystal can thus be broken into three terms: @ = @D&w+
@above
+
@Ewald
9
(4)
i.e., the sum of contributions from layers below and above, and ions in the same layer, if r happens to lie in one of the layers (Q,,,,,). The first two terms can be evaluated using the two-dimen-
(5)
(6)
where zI is the position of the layer just below z. Note that eq. (6) is basically the same as eq. (3) except for a possible change of sign. For Qjabove, there are a finite number of layers, each of which has a contribution given by eq. (5). Problems arise if the ion is too close to a layer or in the layer. Unfortunately, the Fourier transformation method, upon which the derivation of eq. (5) is based, fails when I z - zI I is small compared to a [6]. Alternative approaches have to be developed. In the following discussion, we derive a quasi-two-dimensional Ewald sum to handle the Coulomb field very close or in the plane of an infinite two-dimensional mesh. We first consider the sum over the positive ions. The sum of negative ions can be treated in the same manner. Using the definition of the complimentary error function [18], we have 12m -=& I
exp( -r*k*)
= &kGexp(
-r2k2)
dk
dk + ierfc(Gr).
(7)
186
Z.-H. Huang, H. Guo / Simulation of gas /surface
We define R =xe,
the following
quantities,
+ye,,
interaction
Eq. (11) becomes
(8)
r=R+Aze
Z)
(9)
n = n,e,
+ n,,eY,
( 10)
2.6 = 12
Ccig.R
where e,, eY, and e, are unit vectors along the three axes, and II, and IZ~ are integers. The sum over a layer of positive ions can be written as
c In-f-1l
erfc(GIna +c
n
exp(-Ina-R12k2-Az2k2)
dk
erfc(G)na-rJ) +c n
(m-r1
(11)
’
where G is the convergence factor which controls the convergence of the two sums. Let F(R,
k)=&Eexp(--lnali’k’). n
= cf< g
s,k)
= G/F(R,k)
e’g’R,
(13)
e-ig’R dR,
where A is the total surface area. Substituting (12) into eq. (14) and integrating, we have
f(s,k) = $exp
, i-& 1
which leads to the expression
F(R,k)
2vG = -pp
(17)
.
The integral in the first sum has an analytical form [18], which leads to the following expression, c 1 n Ina-rl =- ,:
c y[eK” g
erfc(&
+ GAzj
erfc(GIna-rI)
where g’s are the two-dimensional reciprocal tice vectors, i.e., g. nu = 27~ X integers. Fourier transform f(g,k) is given as f(g,k)
Ina-rl
(12)
By virtue of its periodicity in (x, y), F(R,k) can be expanded into a Fourier series in the reciprocal space: F( R,k)
n
-rl)
latThe
(14) eq.
(15)
for the F(R,k),
(16)
+c n
Ina-rI
’
( 18)
Eq. (18) is equivalent to the formula stated by Parry [19], who derived it as a limiting case of the 3D Ewald sum by assuming the interlayer spacing in the z-direction to be infinite. The sum in eq. (18) includes a term with g = 0, which is divergent. If the layer is electrically neutral such as a (001) layer of the NaCl crystal, this term can be excluded since the contributions from the positive and negative ions cancel each other. In cases where net charge is not zero, such as in the (111) layer of NaCl, the divergence can also be removed by including multiple layer contributions [20,21]. The g = 0 term has also been found to be important in non-centrosymmetric crystal structures [22]. The formula has been extended and expressed in terms of charge speading functions of arbitrary analytical forms [21]. To evaluate the potential energy as one pulls out an ion from its equilibrium position, one needs to subtract the field generated by the ion itself. Assuming the ion being selected is the one
Z-H. Huang, H. Guo / S~~~t~n
of gas /surface
187
interaction
at r0 = 0, one should replace the n = 0 term in the last sum by the following, lim
r+O
erfc( Gr) r
-2G/v?,
(19)
The intralayer interaction can then be obtained from eqs. (18) and (19). The method introduced above has been applied to the study of adsorption and photodissociation of CH,I adsorbed on a LiF(OO1) surface. We have implemented this method on a Cray YMP and our simulations showed that this method is numerically stable and reasonably efficient. The details of implementation of the method and simulations of photodissociation of methyl iodide on an LiF surface are published elsewhere [231.
3. Comparison
_,4
L-L-LL 0.0
0.5
/ ! .o
1
15
!
/
I
2.0
2.5
3.0
_L_
._..
3.5
i
4.0
zc&
Fig. I. Plot of electrostatic potential energy of a central surface ion as it is pulled away in the direction of surface normal. Other ions in the solid are assumed to be stationary and at their ideal lattice positions. The solid line corresponds to the quasi-two-dimensional Ewald sum method, while the dotted and the dashed lines are for the two-dimensional Fourier transform method and the cluster model, respectively. Variable .z is measured from the lattice positions of the ion.
with other methods
In this section, we will examine the difference between the quasi-ho-dimensional Ewald sum model and two other approximate methods. We choose a LiF(001) surface since we have used both the cluster model and the quasi-two-dimensional Ewald sum method to study photodissociation of several molecules on this surface. For simplicity, we assume that both Li and F carry a charge of “ +e” and “ -e”, respectively, though it may be more reasonable to use partial charges to represent the Li and F ions [24]. We first examine the Coulomb interaction energy q@(r) as one pulls out a central ion in the top layer along the surface normal. Fig. 1 is the plot of the energy calculated using three different methods, the solid line for quasi-two-dimensional Ewald sum, the dotted line for the two-dimensional Fourier method, and the dashed line for the cluster model. There are 6 x 6 X 3 ions in the cIuster model. As shown in fig. 1, our quasi-twodimensional Ewald sum method gives an excellent agreement with the cluster model at small z. In the caIculation of the Ewald sum, the number of g and R points are the same and set to be 36 while the convergence factor G is 0.214 A-‘. The electrostatic potentiat calculated using the
ho-dimensional Fourier transform method, as we pointed out earlier, diverges for small z. Fig. 1 shows that the Fourier method is acceptable only when z > 1 A. Nevertheless, it is a good representation of the field for interlayer interaction. It should be pointed out that eq. (18) is exact and the accuracy can be improved by including more g and n points. The results shown in fig. 1 justify our approach presented in section 1. I.e., for intralayer interaction, the quasi-do-dimensional Ewald sum method is a very good appro~mation, while for intertayer interaction, the Fourier method is quite satisfactory. According to eq. (3), the electrostatic potential on an ionic surface decays exponentially into the vacuum, while for a cluster model, the potential is expected to decrease only as a certain power of r. In fig. 2, we plot the potential energy of a test charge over a central positive ion as a function of z. Again, the cluster model agrees quite well with the Fourier method, with an error in the potential energy only in tens of meV. It is surprising to note the accuracy of the cluster model, even with a relatively small size of 6 X 6 X 3 atoms, is so accurate in the determination of the electrostatic potential above the sur-
Z.-H. Huang, H. Guo / Simulation of gas /surface
188
- - - Fourier
I) 05
1
0.m
i_
’
-n.lls 1
2
3
4
’ 5
6
7
’
’
R
9
1 10
z ca,
Fig. 2. The potential energy of a test charge as a function of z over a central positive ion. The solid and dashed lines are, respectively, for the cluster model and the Fourier transform method.
face. This is probably due to the perfectly regular arrangement of the ions in the cluster, which leads to a cancellation at the adsorption site in the center of the cluster surface. This conclusion offers some justification for the cluster models that have been adopted in many surface related simulations [ll-141. However, one should bear in mind that the cluster model is a finite-size model and caution must be exercised in order to avoid possible errors. In fig. 3, we replot the potential energy using the three methods for an ion at the corner of the cluster. It is clear that due to the finite size, the cluster model is no longer accurate, while the other two methods are not affected. It is expected that Q, calculated using the
interaction
cluster model is position dependent. One may attach the cluster to some fixed ions to avoid the problem, which provides not only a more accurate potential field above the surface but also the forces needed to stabilize the ions at the boundaries of the cluster. However, the number of fixed ions could be significantly larger than the moving ions. For a 6 X 6 X 3 cluster, if two fixed layers of ions are added, we have 232 fixed ions. Undoubtedly, the evaluation of potential energies and forces is now more costly since for each ion, an extra sum over the fixed ions needs to be performed. We have used both a cluster model and the quasi-two-dimensional method in computer simulations of photodissociation of small molecules adsorbed on a LiF(001) surface, and found that the efficiency is comparable for the two methods.
4. Summary In this paper, we have introduced a new method for treating the long-range Coulomb interaction in a system consisting of moving ionic surfaces. The calculation of the Coulomb interaction energy is broken down into sums over layers. For an ion in one layer interacting with another layer, a two-dimensional Fourier method is used to evaluate the energy. For an ion interacting with other ions in the same layer, a quasi-two-dimensional Ewald sum is performed. The method is computationally efficient and yet able to evaluate the Coulomb potential to arbitrary precision. It is intended to be implemented in Monte Carlo or molecular dynamics simulations.
Acknowledgement
This work was supported by the National Science Foundation (Grant No. CHE-9116501).
References [I] R.B. Gerber,
Fig. 3. Same as in fig. 1 except for a corner
ion.
Chem. Rev. 87 (1987) 29. [2] J.C. Polanyi and H. Rieley, Photochemistry in the Adsorbed State, Eds. C.T. Rettner and M.N.R. Ashfold (The Royal Society of Chemistry, Cambridge, 1991);
Z.-H. Huang, H. Guo / Simulation of gas /surface
[3] [4] [5] [6] [7] [8] [9] [lo] [Ill [12]
X.-L. Zhou, X.-Y. Zhu and J.M. White, Surf. Sci. Rep. 13 (1991) 73. G.A. Somorjai, Chemistry in Two-dimensions: Surfaces (Cornell University Press, Ithaca, NY, 1981). D.A. King, Surf. Sci. 47 (1975) 384. R.D. Ramsier and J.T. Yates, Jr., Surf. Sci. Rep. 12 (1991) 243. W.A. Steele, The Interaction of Gases with Solid Surfaces (Pergamon, New York, 1974). A.M. Stonehem and J.H. Harding, Annu. Rev. Phys. Chem. 37 (1986) 53. J.C. Tully, Annu. Rev. Phys. Chem. 31 (1980) 319. J.C. Polanyi, R.J. Williams and S.F. O’Shea, J. Chem. Phys. 94 (1991) 978. R.R. Lucchese and J.C. Tully, J. Chem. Phys. 80 (1984) 3451. MI. McCarthy and R.B. Gerber, J. Chem. Phys. 93 (19901 887. MI. McCarthy, R.B. Gerber and M. Shapiro, J. Chem. Phys. 92 (1990) 7780.
interaction
189
(131 Z.-H. Huang and H. Guo, J. Chem. Phys. 96 (1992) 8564. [14] Z.-H. Huang and H. Guo, J. Chem. Phys. 97 (1992) 2110. [15] J.M. Watson, I. NoorBatcha and R.R. Lucchese, J. Chem. Phys. 96 (1992) 7771. [16] Simulation of Liquids and Solids, Eds. G. Ciccotti, D. Frenkel and I.R. McDonald (North-Holland, Arnsterdam, 1987) ch. 5. 1171 P. Ewald, Ann. Phys. 64 (1921) 253. [18] Handbook of Mathematical Functions, Eds. M. Abramowitz and I.A. Stegun (Dover, New York, 1970) ch. 7. [19] D.E. Parry, Surf. Sci. 49 (1975) 433. [20] D.E. Parry, Surf. Sci. 54 (1976) 195. [21] D.M. Heyes and F. van Swol, J. Chem. Phys. 75 (1981) 5051. [22] M.W. Deem, J.M. Newsam and S.K. Sinha, J. Phys. Chem. 94 (1990) 8356. [23] Z.-H. Huang and H. Guo, J. Chem. Phys., in press. [24] G. Dolling, H.G. Smith, R.M. Nicklow, P.R. Vijayaraghavan and M.K. Wilkinson, Phys. Rev. 168 (1968) 970.