Molecular Dynamics and the Diatomic Molecule as a Harmonic Oscilla= tor 2

Molecular Dynamics of Water Mono= mer[1]

C. W. David

Department of Chemistry

University of Connecticut

Storrs, Connecticut 06269-3060

Carl.David@uconn.edu

1/26/2006

Abstract<= /p>

&= nbsp; Simulation of the classical molecular dynamics of a water molecule can be useful in explaining nomral modes of motion, Fourier Transforms, and fundamental frequencies of vibration, as illustrated herein.

Introduction<= /b>

When learning about molecular vibrations, students can become mystified about the relationship = of normal modes to actual atomic coordinates and their time variations. They a= lso need help in understanding why the potential energy function, usually expre= ssed in simple terms vis-à-vis bond length extensions, bond angle distort= ions and dihedral angle deformations, needs to be re-constituted into forms appropriate to those normal modes. Finally, the static diagrams used in illustrating normal modes do not necessarily convey the compete information needed to aquire mastery of this particular subject.

It is possible to c= arry out a molecular dynamics simulation of a simple molecule, water in our case, which allows students to not only “see” each aspect of the moti= on (and of the simulation itself) but manipulate the parameters in their poten= tial energy model to more fully appreciate molecular vibrations (and rotations a= nd translations, vide ante). In ea= rlier work the HCℓ molecule was used to illustrate how molecular dynamics c= ould be carried out. Since the simulations presented here are done in a symbolic algebraic system, students remain close to the analytical formulations they= are used to without becoming distracted by high level coding in Fortran or othe= r (more efficient) programming language. Finally, an appreciation of the Fourier Transform which allows “picking out” the vibrational frequencies from the molecular dynamics run’s output data can be enhanced with a hands on activity which transcends the use of them as applied herein, and extends the understanding to other spheres in chemistry where this transfor= m is used.

The Molecular Mechanics Model

            We are concerned here with a single molecule of water (larger systems might so= on exhaust the resources of table top computers). The simplest forcefield for = the molecule H₁-O-H₂ is

        &= nbsp;           &nbs= p;    (1)

where each OH distance is d= istorted relative to the equilibrium value (~0.95 Å), and the H₁-O-= H₂ bond angle, = J, is distorted relative to its equil= ibrium value of about 104.5^o. H₁ and H₂ are labels which allow us to distinguish the two protons when and if we choose to treat HDO or HTO or DTO, the three isotopomers with unsymmetrical masses, or D2O or T₂O.

            The problem for programmers in dealing with the vibrational motions of this molecule is evaluating the force on each nucleus. These forces consists of three components on each atom(nucleus), x-, y- and z- each, and corresponds vectorially to the equationwhere i =3D1 an= d 3, say, refer to protons, and i=3D2 refers to the oxygen. The choice is arbitr= ary. Since , it is a purely mechanical task to obtain expressions for the nine partial derivatives required to obtain the three force vectors. We= use Maple to do this, since we want those expressions for instantaneous evaluat= ion when needed during a single step of the molecular dynamics simulation. Then, since , we can calcul= ate the acceleration each nucleus experiences, given the force acting on it, and from there calculate the cha= nge in velocity (and the change therefore of position), i.e., integrate Newton’s equations for all nine “coordinates” pertinant to the (water) molecule.

      &= nbsp;     Once we have the simulation running (debugged), we can change the parameters of = the potential energy model and predict the values of the frequencies which we expect to see in the IR (i.e., predict the normal modes’ frequencies = of vibration). Working in 9 coordinates, while needing only three to describe = the vibrations means that we have six “unused” coordinates, three f= or translation and three for rotation, which we could activate by suitable cho= ices of initial positions and velocities of the 3 atoms (9 initial conditions in position, 9 initial conditions in velocity). Whether translating as an enti= ty (or not) and whether rotating (or not) the vibrations (when small) remain essen= tially the same.

      &= nbsp;     Finally, to “understand” the normal modes, we can set the water molecule’s initial nuclear coordinates such that one of the three nor= mal modes is (approximately) dominant, and watch for the dominant frequency of vibration which emerges, thus assigning (at least initial displacements) a normal mode to its associated frequency of vibration. Tinkering with the fo= rce constants while exercising a particular normal mode allows us to obtain the value of that force constant which generates a frequency which agrees best = with the experimental one.

Implementing The Molecular Mech= anics Model in Maple

            All Maple programs begin with some prologue material. In our case, we need to m= ake the Fourier Transform available, hence “inttrans” in the following code fragment:

= re= start;

= wi= th(inttrans):

= de= lta_r :=3D 0.1e-8;

= de= lta_theta :=3D 0.1;

= h = :=3D 1.0e-16;#(this is the time step in seconds)GOOD VALUE for m=3D10=

= kO= H :=3D 7.76e5;#dynes/cm, pg 218 Barrow

= ka= lpha :=3D 0.699e5*r_e^2;#dynes/rad

= m = :=3D 10;# FFT power of 2

We choose to use 0.1 Angstrom as the amount of displacement an OH bond will be subject to and 0.1 degrees as the amount that the H-O-H angle will be disto= rted from its equilibrium position. Further, we pick a time step (in seconds) corresponding to 10 picoseconds.

        &= nbsp;   Lastly, we choose the numerical values of the force constants we are going to use in modeling the motion of this molecule, as well as the power of two (2^{m<=
/sup>)
setting for the number of Fourier Transform points that we are going to col=
lect.}

=             N= ext, we define a magnitude function:

= ma= g :=3D proc (RealPart,ImaginaryPart)

= sq= rt(RealPart^2+ImaginaryPart^2);

= en= d proc;

= pr= intlevel :=3D 0;

an= d set the printlevel to a non-debugging value of 0.

        &= nbsp;   We will need the distance of separation between two nuclei, and define a funct= ion which will obtain that value:

= di= st :=3D proc(r1,r2) local t;

= t = :=3D sqrt((r1[1]-r2[1])^2+(r1[2]-r2[2])^2+(r1[3]-r2[3])^2);

= re= turn(t);

= end proc;

We are ready to compute the forces as functions. We define one of the H atoms coordinates as x₁₁₌, x_12,and x_13,i.e., , while we define the oxygen’s coordinates simply b= y O₁, O₂, and O₃. This means that the difference vector which we he= re designate as is given by= :

= t1= :=3D [x11,x12,x13]-[O1,O2,O3]:

= t2= :=3D [x21,x22,x23]-[O1,O2,O3]:

We= need the magnitudes of these vectors to calculate the forces, so we use our previously defined mag function:

= t1= _mag :=3D sqrt(t1[1]^2+t1[2]^2+t1[3]^2):

= t2= _mag :=3D sqrt(t2[1]^2+t2[2]^2+t2[3]^2):

La= stly, we need to know the magnitude of the instantaneous value of the H-O-H bond angle which we obtain using one of the vector dot product formulations:

= th= eta_mag :=3D arccos((t1[1]*t2[1]+t1[2]*t2[2]+t1[3]*t2[3])/(t1_mag*t2_mag)):

i.e., . Now, we have all the functions necessary to obtain a functional description of the potential energy function for an isolated wat= er molecule, subject to the simplest model possible. We have for the simplest potential energy model (translated into Maple):

= t3= :=3D (kOH/2)*((t1_mag-r_e)^2+(t2_mag-r_e)^2)+(kalpha/2)*(theta_mag-theta_e)^2:

We= note that r_e is the designated “equilibrium” bond length,= the value of r for which the potential energy curve has a minimum. Now we can t= ake the nine partial derivatives required to get the three components of force = for each of the three nuclei:

= d1= 1 :=3D diff(t3,x11):

= d1= 2 :=3D diff(t3,x12):

= d1= 3 :=3D diff(t3,x13):

= d2= 1 :=3D diff(t3,x21):

= d2= 2 :=3D diff(t3,x22):

= d2= 3 :=3D diff(t3,x23):

= dO= 1 :=3D diff(t3,O1):

= dO= 2 :=3D diff(t3,O2):

= dO= 3 :=3D diff(t3,O3):

For future use, we define a potential energy function which will be useful later (the last line is “returned” by the subroutine as its “value”):

= V = :=3D proc(r_H_1,r_H_2,r_O) local t1,t2,t1_mag,t2_mag,theta_mag;

= t1= :=3D r_H_1-r_O;

= t2= :=3D r_H_2-r_O;

= t1= _mag :=3D sqrt(t1[1]^2+t1[2]^2+t1[3]^2);

= t2= _mag :=3D sqrt(t2[1]^2+t2[2]^2+t2[3]^2);

= th= eta_mag :=3D arccos((t1[1]*t2[1]+t1[2]*t2[2]+t1[3]*t2[3])/(t1_mag*t2_mag)):

= (k= OH/2)*((t1_mag-r_e)^2+(t2_mag-r_e)^2)+(kalpha/2)*(theta_mag-theta_e)^2:

= end proc;

We need some water-specific constants defined:

= r_= e :=3D 0.9584e-8;#cm

= th= eta_e :=3D evalf((104.5)*Pi/180);

= th= eta_e_delta :=3D evalf(delta_theta*Pi/180);

=

= m_= H_1 :=3D 1.0078/(6.023e23):#grams/atom

= m_= H :=3D m_H_1;#choose your isotope

=

= m_= O :=3D 16/(6.023e23);#choose your isotope

(s= ubstituting 2.014102 amu for 1.0078 would allow us to treat D₂O, vide infra)

It is time to set up the normal mode indexing that will be used. Our indexing system defines “norm_mode” =3D 0 to me= an we will set the initial coordinate= s, while “normal mode” =3D1, 2, or 3 means we are attempting the textbook standard displacements according to which normal mode we are interested in,= as indicated in the code:

= no= rm_mode :=3D 1;

= r_= O :=3D [0,0,0];#set this in common, but change later

= if norm_mode =3D 0 then

= r_H_1 :=3D [0.9e-8,0.5e-8,0]:=

= r_H_2 :=3D [-0.7e-8,0.6e-8,0]:

= el= if norm_mode =3D 2 then

= #W= AGGING

= pr= int(` theta varying normal mode`);

= x_H_1 :=3D r_e*cos(theta_e/2+theta= _e_delta);

= y_H_1 :=3D r_e*sin(theta_e/2+theta_e_delta);

= pr= int (`x,y =3D `,x_H_1, y_H_1);

= x_H_2 :=3D x_H_1;

= y_H_2 :=3D - y_H_1;

= r_H_1 :=3D [x_H_1,y_H_1,0];

= r_H_2 :=3D [x_H_2,y_H_2,0];

= pr= int (`h-h dist =3D `,dist(r_H_1,r_H_2));

= el= if norm_mode =3D 1 then

= #S= YMMETRIC STRETCH

= pr= int(` symmetric stretch normal mode`);

=

= x_H_1 :=3D (r_e+delta_r)*cos(theta= _e/2);

= y_H_1 :=3D (r_e+delta_r)*sin(theta= _e/2);

= x_H_2 :=3D x_H_1;

= y_H_2 :=3D -y_H_1;

= r_H_1 :=3D [x_H_1,y_H_1,0];

= r_H_2 :=3D [x_H_2,y_H_2,0];

= pr= int (`h-h dist =3D `,dist(r_H_1,r_H_2));

=

= el= if norm_mode =3D 3 then

= #A= NTISYMMETRIC STRETCH

= pr= int(` ANTI symmetric stretch normal mode`);

=

= x_H_1 :=3D (r_e+0.2e-8)*cos(theta_= e/2);

= y_H_1 :=3D (r_e+0.2e-8)*sin(theta_= e/2);

= x_H_2 :=3D x_H_1;

= y_H_2 :=3D -y_H_1;

= r_H_1 :=3D [x_H_1,y_H_1,0];

= r_H_2 :=3D [x_H_2,y_H_2,0];

= r_= O :=3D [0.1e-8,0,0];

=

= pr= int (`h-h dist =3D `,dist(r_H_1,r_H_2));

=

= en= d if;

= pr= int (`starting H1 coords =3D `,r_H_1);

= pr= int (`starting H2 coords =3D `,r_H_2);

= pr= int (`starting O coords =3D `,r_O);

We= have included some explanatory printing to help understand what is going on.

        &= nbsp;   Arbitrarily, we assign all the nuclei zero velocity (initially) although experimentation here is certainly encouraged. We also zero the forces!

= v_= H_1 :=3D [0,0,0]:

= v_= H_2 :=3D [0,0,0]:

= v_= H_1_p :=3D [0,0,0]:

= v_= H_2_p :=3D [0,0,0]:

= v_= O :=3D [0,0,0]:

= v_= O_p :=3D [0,0,0]:

=

= F_= H_1 :=3D [0,0,0]:

= F_= H_2 :=3D [0,0,0]:

= F_= O :=3D [0,0,0]:

= F_= H_1_p :=3D [0,0,0]:

= F_= H_2_p :=3D [0,0,0]:

= F_= O_p :=3D [0,0,0]:

Next, we compute the potential energy at the start of the motion:

= V_= start :=3DV(r_H_1,r_H_2,r_O):

= pr= int (`starting potential energy =3D `,V_start);

He= re, we use “m” defined at the outset, and set up some arrays which wil= l be used to gather information:

= n_= stop :=3D 2^m;

= l = :=3D array(1..n_stop);

= y = :=3D array(1..n_stop);

        &= nbsp;   All the preparatory matters have been taken care of and we are ready to start s= imulating, one time step at a (computer) cycle. We need a control statement:

For I from 1 by 1 to n_stop do

And then, as the saying went, “away we go”.

Inside each time step, preparatory to calculating the inst= antaneous force on each nucleus, we substitute the now current values of all the coordinates into the appropriate partial derivative terms, and evaluate them numerically. The minus sign is used to connect the force to “minus= 221; the gradient, as noted above. We have :

= t5= :=3D subs(x11=3Dr_H_1[1],x12=3Dr_H_1[2],x13=3Dr_H_1[3], x21=3Dr_H_2[1],x22=3Dr_H_2[2],x23=3Dr_H_2[3],

= O1= =3Dr_O[1],O2=3Dr_O[2],O3=3Dr_O[3],

= d1= 1):

= F_= H_1[1] :=3D -evalf(t5):

=

= t5= :=3D subs(x11=3Dr_H_1[1],x12=3Dr_H_1[2],x13=3Dr_H_1[3], x21=3Dr_H_2[1],x22=3Dr_H_2[2],x23=3Dr_H_2[3],

= O1= =3Dr_O[1],O2=3Dr_O[2],O3=3Dr_O[3],

= d1= 2):

= F_= H_1[2] :=3D -evalf(t5):

=

= t5= :=3D subs(x11=3Dr_H_1[1],x12=3Dr_H_1[2],x13=3Dr_H_1[3], x21=3Dr_H_2[1],x22=3Dr_H_2[2],x23=3Dr_H_2[3],

= O1= =3Dr_O[1],O2=3Dr_O[2],O3=3Dr_O[3],

= d1= 3):

= F_= H_1[3] :=3D -evalf(t5):

=

= t5= :=3D subs(x11=3Dr_H_1[1],x12=3Dr_H_1[2],x13=3Dr_H_1[3], x21=3Dr_H_2[1],x22=3Dr_H_2[2],x23=3Dr_H_2[3],

= O1= =3Dr_O[1],O2=3Dr_O[2],O3=3Dr_O[3],

= d2= 1):

= F_= H_2[1] :=3D -evalf(t5):

=

= t5= :=3D subs(x11=3Dr_H_1[1],x12=3Dr_H_1[2],x13=3Dr_H_1[3], x21=3Dr_H_2[1],x22=3Dr_H_2[2],x23=3Dr_H_2[3],

= O1= =3Dr_O[1],O2=3Dr_O[2],O3=3Dr_O[3],

= d2= 2):

= F_= H_2[2] :=3D -evalf(t5):

=

= t5= :=3D subs(x11=3Dr_H_1[1],x12=3Dr_H_1[2],x13=3Dr_H_1[3], x21=3Dr_H_2[1],x22=3Dr_H_2[2],x23=3Dr_H_2[3],

= O1= =3Dr_O[1],O2=3Dr_O[2],O3=3Dr_O[3],

= d2= 3):

= F_= H_2[3] :=3D -evalf(t5):

=

= t5= :=3D subs(x11=3Dr_H_1[1],x12=3Dr_H_1[2],x13=3Dr_H_1[3], x21=3Dr_H_2[1],x22=3Dr_H_2[2],x23=3Dr_H_2[3],

= O1= =3Dr_O[1],O2=3Dr_O[2],O3=3Dr_O[3],

= dO= 1):

= F_= O[1] :=3D -evalf(t5):

= t5= :=3D subs(x11=3Dr_H_1[1],x12=3Dr_H_1[2],x13=3Dr_H_1[3], x21=3Dr_H_2[1],x22=3Dr_H_2[2],x23=3Dr_H_2[3],

= O1= =3Dr_O[1],O2=3Dr_O[2],O3=3Dr_O[3],

= dO= 2):

= F_= O[2] :=3D -evalf(t5):

=

= t5= :=3D subs(x11=3Dr_H_1[1],x12=3Dr_H_1[2],x13=3Dr_H_1[3], x21=3Dr_H_2[1],x22=3Dr_H_2[2],x23=3Dr_H_2[3],

= O1= =3Dr_O[1],O2=3Dr_O[2],O3=3Dr_O[3],

= dO= 3):

= F_= O[3] :=3D -evalf(t5):

Ha= ving the instantaneous forces, we now can obtain the (primed) new positions of t= he nuclei as predicted from the old ones, the velocities of the nuclei and the forces on the nuclei, according to Verlet’s algorithm.

= r_H_1_p :=3D r_H_1 + h*v_H_1 + ((h^2)/(2*m_H))*F_H_1;

= r_H_2_p :=3D r_H_2 + h*v_H_2 + ((h^2)/(2*m_H))*F_H_2;

=

= r_O_p :=3D r_O + h*v_O + ((h^2)/(2*m_O))*F_O;

(W= ere we to allow the two masses of the H nuclei to be different, e.g., allowing one= to be about double of the other (for example), then we could treat HDO as well. This would require slight coding changes in the center line (above). Of cou= rse, we could add code for Tritium in the mass definitions, and treat all the isotopomers.)

Once we’ve “moved” the nuclei, we need to recalculate the forces, which is indicated below (some repetitive code has = been removed).

= t5= :=3D subs(x11=3Dr_H_1_p[1],x12=3Dr_H_1_p[2],x13=3Dr_H_1_p[3],x21=3Dr_H_2_p[1],x2= 2=3Dr_H_2_p[2],x23=3Dr_H_2_p[3],

= O1= =3Dr_O_p[1],O2=3Dr_O_p[2],O3=3Dr_O_p[3],

= d1= 1):

= F_= H_1_p[1] :=3D -evalf(t5):

=

= t5= :=3D subs(x11=3Dr_H_1_p[1],x12=3Dr_H_1_p[2],x13=3Dr_H_1_p[3],x21=3Dr_H_2_p[1],x2= 2=3Dr_H_2_p[2],x23=3Dr_H_2_p[3],

= O1= =3Dr_O_p[1],O2=3Dr_O_p[2],O3=3Dr_O_p[3],

= d1= 2):

= F_= H_1_p[2] :=3D -evalf(t5):

= #r= epetitive code omitted here

= t5= :=3D subs(x11=3Dr_H_1_p[1],x12=3Dr_H_1_p[2],x13=3Dr_H_1_p[3],x21=3Dr_H_2_p[1],x2= 2=3Dr_H_2_p[2],x23=3Dr_H_2_p[3],

= O1= =3Dr_O_p[1],O2=3Dr_O_p[2],O3=3Dr_O_p[3],

= dO= 3):

= F_= O_p[3] :=3D -evalf(t5):

=

= We= also need to recompute the velocities (primed) to their new values, based on the= new forces (primed)
v_H_1_p :=3D v_H_1 + (h/(2*m_H))*(F_H_1_p + F_H_1):

= v_= H_2_p :=3D v_H_2 + (h/(2*m_H))*(F_H_2_p + F_H_2):

= v_= O_p :=3D v_O + (h/(2*m_O))*(F_O_p + F_O):

=

= r_= H_1 :=3D r_H_1_p:#update coordinates

= r_= H_2 :=3D r_H_2_p:

= r_= O :=3D r_O_p:

(W= e are again forced to change the coding slightly if we are interested in treating water with two unequal hydrogen isotopes.) Next, we need to update the coordinates and velocities from the new values in this cycle to the “= to be old values” for the next cycle, i.e.,

= v_= H_1 :=3D v_H_1_p:#update velocities

= v_= H_2 :=3D v_H_2_p:

= v_= O :=3D v_O_p:

We’ve achieved what was requi= red, one cycle of the simulation. Now we need to do some auxiliary computations, which will allow analysis of our results. First, we compute the magnitudes = of the velocities of the three nuclei, precursors of the computation of the to= tal kinetic energy of the molecule:

= v_= H_1_mag :=3D sqrt(v_H_1[1]^2+v_H_1[2]^2+v_H_1[3]^2);#for kinetic energy<= /span>

= v_= H_2_mag :=3D sqrt(v_H_2[1]^2+v_H_2[2]^2+v_H_2[3]^2);

= v_= O_mag :=3D sqrt(v_O[1]^2+v_O[2]^2+v_O[3]^2);

Ne= xt, we zero the imaginary part of the Fourier transform output (not really necessa= ry) and then compute the kinetic energy, the potential energy, and the total en= ergy at this time step.We have:

= y[i] :=3D 0;#imaginary part

=

= KE(i) :=3D 1/2*(m_H*v_H_1_mag^2+m_H*v_H_2_mag^2+m_O*v_O_mag^2):#

= PE(i) :=3D V(r_H_1,r_H_2,r_O):

= E_total(i) :=3D evalf(KE(i)+PE(i))= :

= E_totaloverV[i] :=3D (E_total(i)/V_start)*100;

= E_totaloverV_plot(i) :=3D (E_total(i)/V_start)*100;

=

Lastly, but most importantly from t= he perspective of visualization, we compute the H₁-H₂ distance (instantaneously) for two purposes. First, we wish to plot this resultant value set as a function to “time” (hence r_plot) and second, we will carry out our Fourier transform on it also (hence l[i]),

= l[i] :=3D dist(r_H_1,r_H_2)*1e8:#co= nvert to angstrom

= r_plot(i) :=3D l[i]:

= yy[i] :=3D 0;#imaginary part<= /o:p>

and then we’re done, with a single time step. The loop proceeds n_stop times.

= en= d do;

Next, we carry out the Fourier Transform on l[i],y[i]

= i = :=3D 'i':

= FF= T(m,l,y);

and plot the internuclear (proton-proton) separation as a function of “time”:<= /span>

= pr= int(` r_plot(i):`);

= PL= OT(POINTS(seq([i,r_plot(i)],i=3D1..n_stop), SYMBOL(DIAMOND),LEGEND(`r versus t`)));

Then, we plot the Fourier transf= orm of that data:

= i = :=3D 'i':

= Fr= eqSpectrum :=3D [seq([(i-1),(2*mag(l[i],y[i])/(2^m))],i=3D1..floor(((2^m)/2)))]:<= /o:p>

= #p= lot([seq(FreqSpectrum[j],j=3D2..floor(((2^m)/2)))],title=3D"Fourier Transform");

= pl= ot([seq(FreqSpectrum[j],j=3D2..40)],title=3D"Fourier Transform");

No= te that the commented out line (above) would plot the entire Fourier Transform, but lowering the upper limit of the sequence allows the reader to better see where the maximum is (or maxima are).

        &= nbsp;   Next, we compute some reportable molecularly interesting results of these computations:

= fr= equency :=3D number_of_cycles/(h*n_stop);

= pr= int(`sec^-1 =3D `,frequency,` times # of peaks`);

= pr= int(`cm^-1 =3D `,evalf(frequency/3e10),` times # of peaks`);

=

=

= pr= int(` KE(i):`);

= PL= OT(POINTS(seq([i,KE(i)],i=3D1..n_stop), SYMBOL(CROSS),LEGEND(`kinetic energy versus t`)));

= pr= int(` PE(i):`);

=

= PL= OT(POINTS(seq([i,PE(i)],i=3D1..n_stop), SYMBOL(CIRCLE),LEGEND(`potential energy versus t`)));

= pr= int(` E_total(i):`);

= l_= plot :=3D [[ n, E_totaloverV(n)] $n=3D1..n_stop]:

= plot(l_plot, n=3D1..n_stop, style=3Dline,symbol=3Dcircle,labels=3D[`time

= step`,`% deviation of total energy`],labeldirections=3D[HORIZONTAL,

= VERTICAL]);

And we’re done. The last bit of code is included for completeness, i.e., it’s optional. However, it is really interesting to see that the kine= tic and potential energies are 180^o of out phase, and that when added together give almost a perfectly constant total energy (one has to check the ordinate of the energy plot to see that the variations in total energy are incredibly small, but real). These last plots are not included in this manuscript.

Using the Code (1)

Choosing force constants, (7.76x10⁵dynes/cm for= k_OH and 0.699x10⁵*(0.9584)² dynes/rad) from older literat= ure allows us to test the basic concept, i.e., run the code using the three nor= mal mode choices, to check that we get suitable harmonic motion. For the bond a= ngle stretch mode, we obtain nice, in fact beautiful, harmonic motion (see Figure 1).

The Fourier Transform of this motion is shown in Figure 2, where one sees that there are about 9 cycles during the total sample.

For the anti-symmetric mode on the other hand, we find, indeed, that the motion is anything but harmonic, which is quite suggestive= of the need to improve the initial coordinate guesses that we made. Figure 3 s= hows this complex motion (as reflected in the r_HH(t) data). Perhaps it might be wise to tinker with the initial velocities as well as the initial positions in this case.

Using the Code (2)

Th= ere is no question that seeing the motion vividly displayed this way, accompanied = by the Fourier Transform which informs us of the frequencies, is gratifying. B= ut one can do more with the simulation.

Fi= rst, one can obtain the best force constants (inside the model definition) conso= nant with the experimental data. Some experimental data is shown in Table 1.

We= can adjust the time step (h) in the code, to get as close as possible to an integral number of cycles inside the 2^m datum, and then calculate the resultant frequency. Tinkering with the force constants would then allow getting the closest fit of simulation frequencies to experimental frequenci= es. Then, with just a slight alteration in the coding, one could, for instance, predict the spectrum of HDO.   Consider that the frequency (u_2<=
/span>) = is given primitively by the expression:

which translates to <= /span>

Tinkering with the two force constants and getting a closer approximation to an integral number of cycles/experiment (varying h) yields frequencies which can be made to come = as close as one desires to the experimental value.

By changing the velocity components= , one can see that they are, when not overly large, ignorable. This means that the separation of vibration from rotation and translation is warranted under th= ose assumptions.

By changing the initial coordinates= to reflect enormous distortions, one can see the motion deteriorate rapidly in= to something ostensibly chaotic.

Using the Code (3)

Another application which might be of interest consists of altering the potential energy model employed. Clearly, one could add cross terms to the simplified function (see Equation 1). Completely different mod= els might also be worthy of exploration. Finally, extending the code to include more than 3 nuclei, or changing the nature of the nuclei considered in the 3 nucleus case, is straightforward.

Discussion

It is easy to hand wave arguments in elementary physical chemistry which lead students to a glib form of understanding of topics whi= ch, upon closer investigation, are really completely confusing to students. The example used here allows a form of concrete realization of vibrational ener= gy concepts which elude understanding under normal teaching conditions. Its inclusion as an exercise, especially for students who will not continue in physical chemistry, seem warranted. The hard thing is to understand the V=3Df(nine coordinates) and therefore that therefore are also functions of nine coordinates, and is non-trivial to derive and/or evaluate.
Tables

=

= H₂O(6)

= D₂O⁽⁵⁾

= H₂O(alternate)(7)

n_=
1=

= 3832.17

= 2763.80

= 3657

n_=
2=

= 1648.47

= 1206.39

= 1594.7

n_=
3=

= 3942.53

= 2888.78

= 3755.7

Ta= ble1: Experimental IR frequencies for water.

Fi= gure 1: Bond angle stretching normal mode, showing beautiful harmonic motion.

Figure 2. The Fourier Transform of the oscillation shown in Figure 1. There are about 9 cycles in the total motion displayed.

= <= span style=3D'mso-no-proof:no'>
Figure 3: The r_HH distance as a function of time, showing a comp= lex motion indicative that a single normal mode is not being used.

Literature Cited

(1)  &n= bsp;     Rempe, S. B.; Jonson, H Chem. Educator= , 1998, 3, 1

(2)  &n= bsp;     David, C. W. Journal of Chemical Education= , submitted

(3)  &n= bsp;     Verlet, L Phys. Rev. 1967, 98, 159<= /o:p>

(4)  &n= bsp;     7.684 Millidyne/Angstrom and 0.707*0.954 Millidyne/Angstrom resppectively, Bernat= h, P. F. Spectra of Atoms and Molecules, Oxford Univ. Press, 1995, page= 223

(5)  &n= bsp;     Barrow, G. M. Introduction to Molecular Spectroscopy, 1962, page 218.

(6)  &n= bsp;     Kuchitsu, K; Bartell, L. S. J. Chem. Phys= ., 1962, 36, 2460

(7)  &n= bsp;     Herzberg, G Molecular Spectra and Molecular Structure, D. van Nostrand Co., Inc. 1966, page 585

=	= H₂O(6)	= D₂O⁽⁵⁾	= H₂O(alternate)(7)
n_= 1=	= 3832.17	= 2763.80	= 3657
n_= 2=	= 1648.47	= 1206.39	= 1594.7
n_= 3=	= 3942.53	= 2888.78	= 3755.7