edited by
JAMES P. BIRK Arizona Stale University
computer series, 152 A Simple Undergraduate Computer Modelling Experiment for PC Computer Systems Michael W. Jezercak University of Central Oklahoma Edmond, OK 73034 In recent vears. the increase in c o m ~ u t i n eDower and availability has led to rapidly growing a&& in the area of theoretical chemistrv. The abilitv to carrv out lendhv calculations a t higher speeds have illowed iesearche& tb look deeply into the energies and mechanisms involved in chemical reactions and scattering. With the lower cost and increased power of small personal computer systems, it is now possible for the undergraduate laboratories to introduce students to this growing field of chemistry. Presented is a theoretical one-dimensional atomic scattering experiment developed for the undergraduate physical chemistry laboratory, that is small enough for students to carry out and run on PC computer systems. The experiment is intended to model the scattering or trapping and desorption of a gas atom from a n attractive surface. Students gain pract&al experience in computer programming, elementary mathematics, and numerical analysis. More importantly, students are able to investigate the sensitivity of a variety of factors on scattering and absorption including vibrational frequencies, dissociation energies, force constants, and reduced mass. Instructors also will be able to include discussions on Maxwell-Boltzmann statistics. activation energies, and crystal growth.
The Model The subiect chosen for studv is a one-dimensional scattering of a single gas atom from a "surface" consisting of one atom fixed in place by a harmonic spring as illustrated in Figure 1. This system was chosen for the following reasons (1)Because the system is held to one dimension, there is a dramatic reduction in the complexities of the problem that can overshadow desired leamine ~0intswhile allowing students to complete the probl& in a reasonable amount of time. (2) This model provides a simulation of atomic trapping and crystal growth. Students may investigate the conditions required to form molecules and crystals. ~~
~
Tempe, A2 85281
(3) This syskm lends itself to variation of parameters such as mass, dissociation energies, ete. that allow students to adjust the system parameters to determine the optimum conditions for molecule formation. (4) This system allows students to work in Cartesian coordinates, with which they are more comfortable, urithout having to work with cross terms that can be discouraging mathematically, ( 5 ) The program is computationally small enough to run in a reasonable amount of time on ~ersonalconmuter svstems allowing departments without access to a large mainframe computer to carry out computer simulations. For this Droblem. students are introduced to Hamilton's equations of motion after a suitable discussion of classical mechanics including Newton's equations of motion (I1. The use of Hamilton's equations doubles the number of differential equations to be solved but is computationally easier in that only one integration is required for each. The total e n e r w for this system is given bv the Hamiltonian and expressed in terms of coordinates and momenta a s ~=~?i2rn + ,~;/2rn, + 0.5hz.2 + V&&
(1)
The first and second terms of eq 1are the kinetic energy expressions for particle s and particle g representing the surface atom and gas atom, respectively The third term is the harmonic potential energy allowing the surface atom the harto vibrate a t a freauencv selected bv " adiustine " monic force constaLt, k. itudents are required adjust k to fit desired vibrational frequencies. The fourth term is the gas-surface interaction potential. Asimple Morse function was selected to describe this interaction and is eiven as
-
V,(r,z,)
= D d l - e x p 1 4 ( x ~ - x ~ b ~ -De ell)~
(2)
This potential energy function is chosen because the three adjustable parameters, De, a, and Re, are related to dissociation energy, vibrational frequency, and equilibrium bond length, respectively. De is subtracted from the potential function to e v e -De as the eauilibrium bond enerw. This potential isuseful for the m&elling of open shell FAteractions and is shown in Figure 2. Hamilton's equations of motion are used to calculate the positions and momenta of the w s and surface atoms a t all times and are e v e n a s follows.
-
-
HI*
= dddt
(3)
and H I & = -dpldt
(4)
Taking derivatives of the system Hamiltonian with respect to the coordinates and momenta of each particle yields four first-order differential equations that are then solved by standard techniques.
gas atom Figure I . One-dimensional scattering model.
Procedure The following four-step procedure for carrying out this calculation is presented to the students. 1. Establish initial conditions and input parameters. 2. Calculate the derivatives (Hamilton's equations). 3. Integrate the equations of motion. Volume 70 Number 8 August 1993
637
4. Calculate new positions and momenta and evaluate trap-
ping and scattering criteria. Initial conditions and parameters are given to the student with a set of exercises. Due to the order of magnitude of numbers in the atomic realm. a rather arbitraw system of atomic units used. In the system chosen, a mass ;nit is selected to be 1amu. an enerw unit to be 1eV, and a distance unit to be 1 ~electing?hese three parameters necessarily defines the time unit. The units are, thus, defined as follows:
A.
1 mass unit = 1 amu 1 energy unit = 1ey 1distance unit = 1 A 1 time unit, tu, = 1.018 x 10-l4 s
This unit system provides numbers with orders of magnitude that &e within the capabilities of small computer svstems. Based on this unit system, a set of suggested parameters corresponding to a hydrogen-hydrogin collision are given in Table 1(2). Great latitude in the selection of these initial conditions is allowed for simulation of many different systems. For example, by placing the gas atom within the potential energy well, vibrational excitation of an atom on the surface may be examined. The surface atom initial position, x, is set to zero. This is the equilibrium position of the harmonic oscillator as given by the potential in eq. 1.At this initial position, all of the energy is in the form of kinetic energy. By assigning the oscillator a given amount of kinetic energy, the student mav then PO on to monitor enerw exchange between the suGace and gas atoms via colli&n. ~ e c a u s ethe atom is bound to the surface, and the surface is taken to be fixed in space, all the available energy for the surface atom is contained in vibrational motion. Because the energy of a n individual atom is 3kT,, by the equipartition theorem, the 1-D bound atom is assigned an energy of kT, reflecting a statistical average (3).This selection certainly is arbitrary as the energy of any single atom at temperature, T, is found within a rather broad distribution; however, this provides a reasonable value for the energy of the oscillator. From the "temperature" of the surface, the momentum of the surface atom is then given as
where rn is a random number selected from a uniform normalized distribution, and k is the surface force constant. This prevents the gas atom from striking the surface atom at the same point in its vibrational period. By averaging over phase angles, variances among different trajectories will be observed. To provide a statistical distribution, students are required to run 50-100 trajectories. These numbers are high enough to see significant effects of varying system parameters while allowing convenient computational t&e. For smaller computer systems, this number may be reduced cautiously. The differential equation integrator used is a fqurthorder Runge-Kutta integrator (5).This usually is given to the students for use in their oromam: but.. the instructor may have the students write the integration routine themselves or use an integrator provided in standard references. Predictor-corrector methods may be used to speed up the calculations once the RungeKutta integrator has started the trajectory, but programming complexity increases. The intemator must be able to find the derivatives given in eqs (3)and (4). The derivatives ofthe Hamiltonian usually are completed analytically and programmed into a subroutine that is accessed by the integrator, which then integrates the equations of motion, to determine the new positions and momenta of the gas and surface atoms. Table 1. Hydrogen-Hydrogen Parameters for OneDimensional Scattering Calculation Surface Force Constant Surface'Temperature" Gas Atom "Temperature" Time Step Morse Potential Parameters: Dissociation Energy Frequency Factor Equilibrium Bond Length Mass of Gas Atom Mass of Surface Atom
De = 4.4755 eV
a = 2.0 A-' Re = 0.7416 A
q=I.Oamu ms= 1.0 amu
The selection of the initial position of the gas atom is arbitrarv. I t usuallv is selected to be large enough to be outside t6e potentiai energy well of the surface atom, but not so large as to use up computational time waiting for the atom to reach the surface. For the potential energy parameters suggested, about 3.0 A usually is sufficient but may need to be increased if the well is extended by the Morse parameter, a. The gas atom initial momentum is selected from the average one-dimensionalkinetic energy of 0.5 kT. at the selected gas temperature. Like the surface atom, this selection is rather arbitrary, but, produces reasonable numbers that also reflect a statistical average. The resulting momentum is thus:
Pg= i m & T g ) In
(6)
If the instructor wishes, the student may be required to ~ a u s s i a nor select t h e initial momentum from Boltzmann distribution, although it is computationally more difficult for the undergradGte student invert such a distribution. If such a distribution is desired, a Monte Carlo method is well suited to this (4). It is critically imoortant to average the collisions of the gas and surface atoms over vibrational phase of the surface atom. This is done easily by adding to the initial gas atom position a random numbeFscaled by an attenuation factor.
a
638
Journal of Chemical Education
0
1
2
3
4
5
6
7
6
Gas Atom-SurfaceAtom Distance (A) Figure 2. Morse potential for gas atom-surface atom interaction with De = 4.4755 eV, a = 2.0 1/A, and Re= 0.7416 A.
Insertion of initial conditions into the integrator begins the program and allows the gas particle to move into and strike the surface atom. Programs can be checked for programming ermrs and numerical accuracy by calculating the energy of the system by eq. 1. This is a conservative system and the total system energy should remain wnstant; however, numerical error can produce variances in the calculated energy. The time step selected gives variances in total energy ofless than one part per 10 thousand. An additional check for program accuracy is a back integration. In this check, a single trajectory is allowed to run for about 500 time steps. The program then changes the sign on the time step, ts,used in the integrator and run for a n additional 500 time steps. The calculated coordinates and momenta of this '%a& integration" should be the reverse of the forward calculation. Finally, criteria must be included into the program to test whether or not t h e atom has been scattered or trapped. Criteria for trapping is somewhat arbitrary. Trapping can be considered to have occurred if the atom undergoes six or more turning points on the surface, that is, if the momentum of the particle has changed direction more than six times. It is imvortant to recoenize that this is a small system. The energy lost to the h;~rmonicosciilator likelv will be transferred back to the
