Normal modes of a triangular molecule

Copyright (c) 1994, Wm. Martin McClain
wmm@chem.wayne.edu
Dept. Chemistry, Wayne State University,
Detroit, Michigan 48202

Intoduction

The normal modes of motion of a molecule with N atoms in three dimensions are determined by 3N second order coupled differential equations; in two dimensions, there are 2N equations. The main mathematical problem is to uncouple the equations; if the harmonic approximation is used this can be done exactly and the solution then follows in closed form by elementary methods. 
							
							In this script we work out a two dimensional problem using a "brute force" method that diagonalizes a 2N-by-2N matrix.  It is called "brute force" because it makes no use of any tricks to reduce the size of the matrix.  These tricks were important in pre-computer days when diagonalization took a major effort, but they are now arguably obsolete.  Indeed, the brute force method works for asymmetric molecules of up to 100 atoms or more,  such as small peptides.  This script shows the simplest possible example, a nice little triangular molecule in two dimensions.
							
							This molecule has bond stretching forces, but because it is triangular, it does not need specific bond angle forces.  After you understand this calculation, you will be able make your own program for molecules that need bond angle forces (though it is not trivial!).  Hints will be available in a separate script, if you need them.
							
							It was at one time thought that triatomic oxygen (ozone, O3) was an equilateral triangle, but it is not.  In fact, no element forms a triatomic molecule that is equilateral in its ground state.  But this calculation played a central role in determining that ozone is not equilateral, and you will see the result that proves this point.  The carbon atoms in cyclopropane, C3 H6, make an equilateral triangle, and its lower frequency vibrations almost obey this model.

Preliminaries needed by the computer

Off[General::spell,General::spell1]

<<LinearAlgebra`Orthogonalization`

<<Miscellaneous`PhysicalConstants`

nAvogadro = AvogadroConstant Mole

cInCmPerSec = SpeedOfLight (100 Second/Meter)

Step 1:  write the potential V in internal coordinates

Since three points determine a plane, all vibrational motions of a triangle are two dimensional.  Three internal coordinates determine completely the shape of the molecule:  We let dr12 be the extension from rest length of the bond between atoms 1 and 2, and similarly for dr23 and dr31. In rest position all extensions are zero;  as the molecule vibrates, they take both positive and negative values.  The vibrational potential energy, written in internal coordinates, is

Vint = kHooke (dr12^2 + dr23^2 +dr31^2);

Step 2:  transform V to Cartesian coordinates

The molecule has atom 1 at {x1,y1}, etc., and the bond lengths are taken to be 2s at rest. Therefore

internalToCartesian = 
{dr12->Sqrt[(x1-x2)^2 + (y1-y2)^2]-2s,
 dr23->Sqrt[(x2-x3)^2 + (y2-y3)^2]-2s,
 dr31->Sqrt[(x3-x1)^2 + (y3-y1)^2]-2s};

Now we use this replacement rule on the potential energy function Vint to give Vcar, the potential energy in terms of cartesian coordinates :

Vcar = Vint /. internalToCartesian

Vectorize and put in numerical constants

Rename all Cartesian coordinates as z[i], with i running 1 to 6 :

toZs = {x1->z[1], y1->z[2],
        x2->z[3], y2->z[4],
        x3->z[5], y3->z[6]};

The bond length parameter s has no influence on the modes or the frequencies, because the problem will be solved solely in terms of the partial derivatives of Vcar wrt the Cartesian coordinates, and s disappears as soon as a partial is taken.  It is therefore convenient to set s to 1 at this point.

We must now pick a value of the Hooke's Law constant for the three equivalent springs.  We take a value, learned by experience, which is not unreasonable for chemical bonds:

data = {kHooke->250000 (*dyne/cm*), s->1};

VcarN = Vcar /. toZs /. data

Step 3:  define a standard equilibrium position

We need a standard equilibrium position for the triangle, at which we will evaluate the partial derivatives.  We take it centered on the origin, peak pointing up.  We number the atoms:

     atom 1= peak,    atom 2 = left end of base,    atom 3 = right end of base.  
         
 A little geometry shows that in this rest position the atom coordinates are

rest = {
	x1->  0,	y1-> +s Sqrt[3.]/2,
	x2-> -s,	y2-> -s Sqrt[3.]/2,
	x3-> +s,	y3-> -s Sqrt[3.]/2}/.toZs/.data

As a geometry check, we calculate the rest values of the atom-atom distances:

{ Sqrt[(x1-x2)^2 + (y1-y2)^2],
  Sqrt[(x2-x3)^2 + (y2-y3)^2],
  Sqrt[(x3-x1)^2 + (y3-y1)^2] }/. 
  					toZs /.rest // PowerExpand

As a further visual geometry check, we draw a picture of the molecule at rest.

restpoints = {rp[1],rp[2],rp[3]}=
	{{x1,y1},{x2,y2},{x3,y3}} /. toZs /. rest

stdPlot={PlotRange->{{-1.4,1.4},{-1.4,1.4}},
		AspectRatio->1};

restPicture = Show[Graphics[
		Line[{rp[1],rp[2],rp[3],rp[1]}]],stdPlot];

Step 4:  expand V about standard position, using the harmonic approximation

If we expand Vcar in multidimensional power series we will have

(1)  We may pick any zero of energy Vo that we desire, and above we have taken the potential energy to be zero when all bond lengths are 2s.  Note that the value of s cannot affect any higher terms because it vanishes when any derivative of Vcar is taken.

(2)  The negative of the  first derivative of Vcar wrt any Cartesian coordinate is the force along that coordinate.  At equilibrium, this is zero.  

Because of these two facts, a multidimensional Taylor expansion of Vcar about an equilibrium configuration is particularly simple:  the zeroth order term vanishes because of (1), and all the first order terms vanish because of (2).  Thus the Taylor expansion starts with second order terms.  

The harmonic approximation is the neglect of cubic and higher terms.  In our simple Vcar the cubic terms vanish exactly, but in real molecules their effects may easily be seen in molecular behaviour. However, the quadratic terms by themselves present the most general problem that can be solved in closed form.  That is the real reason we neglect the higher terms.  And, or course, is true that the lower vibrational levels are often well approximated by the harmonic terms alone.

The numerical value of kMat, in erg / cm^2, or dyne / cm,  is

kMatN = 
Table[D[VcarN,z[i],z[j]]/.rest,{i,1,6},{j,1,6}]

kMatN ==Transpose[kMatN]

Yes, kMatN always is symmetric.  What is it in our construction method that causes this?  It is a little bit hidden, but a theorem that you saw in differential calculus is at the root of it.

Step 5: write the differential equations of motion

Generally, the force on Cartesian coordinate zeta sub i  is

The next step is not totally obvious, so you might want to do this explicitly in Mathematica  for yourself:  It may be shown that in harmonic approximation the general force formula reduces to

But the force is also given by Newton's third law as

Eliminating the force between these last two equations, we have the equations of motion.  There are six of them, as i runs from 1 to 6.

We could solve the six simultaneous equations by an eigenvalue technique if it were not for the mass factor on the left.  To get rid of it, we must make our first transform into an abstract space, the rootmass weighted coordinate space.  Don't try to envision it or think what it means physically; the effort will not yield anything worthwhile.  Trust the algebra.

Step 6:  from Cartesian to rootmass coordinates

We rewrite the six equations of motion in an eccentric, but fruitful, way:

Now dividing through by root m sub i, and regrouping cleverly, we get

Now we make the abstract transform to "rootmass weighted" q coordinates:

and rewrite the six equations in terms of these symbols (i running 1 to 6) :

In rootmass weighted coordinates, the coefficient of the second derivative by time (on the left) has disappeared, and we can solve the six simultaneous equations by the "well known" techniques of linear algebra.  We now return to Mathematica to calculate actual values for the c matrix, which we will call krm (for krootmass).

To show how to treat unequal masses in the general problem, we take atom 1 to be of mass 16.1 dalton (i.e., gram per mole), and atoms 2 and 3 of mass 16.0 dalton.  This small inequality in masses has another good effect on the problem that you will see later on.

We now define a redundant mass vector, in which the mass of each particle is repeated twice in a row (in a three dimensional problem, three times in a row).  In other words, each coordinate has its own mass.  In this formalism it is a miracle, but nevertheless it is always true, that the x, y, and z coordinates of each particle have the same mass.  That is just a strange property of all "particles".  This long mass vector makes the algebra go very smoothly, both symbolically and in computation.

massVec = 
{16.1,16.1,16.,16.,16.,16.}/nAvogadro; (*grams*)

Now construct the root-mass weighted force constant matrix:

krm = 
kMatN/Outer[Times,Sqrt[massVec],Sqrt[massVec]];

krm==Transpose[krm]

Why is krm always symmetric?  Why is it important that it be symmetric?

Step 7:  from rootmass to normal coordinates; solve the eqns of motion in normal coordinates

Symbolic introduction

Suppose  it were within your power to start with matrix c , and from it to calculate two matrices S  and D , such that S , c , and D  would obey

If you had such  S  and D matrices, what could you do with them? You might apply them to the equations of motion in the following way:

You would then be ready for a second abstract coordinate transform.  You could let

when DsubM is not zero.  If it is zero, the solution is (in Newton's words) "a state of rest, or of uniform motion in a straight line":

VoilÇ.  The differential equations are solved. Do not feel disappointed that you can not yet see how the molecule moves;  that will come, but it requires a lot of back transformation.  Before we get to that, however, let us carry out concretely the program outlined above.

Calculation of  Smat

It is in fact quite possible to start with c (or krm) and find S and D. It is done by the Mathematica command Eigensystem, plus some clearing up with the GramSchmidt operator.  We multiply matrix krm by 10^-28 for this step, so that Chop will work correctly :

{scaledVals,SmatRawT}=Eigensystem[10^-28 krm]//Chop

In scaledVals, note that two of the modes are nearly degenerate, differing only in the sixth decimal place.  If we had not made one atom slightly heavier than the the other two, this degeneracy would have been exact.  We do have an exact degeneracy in the last three eigenvectors, all of which have 0 as their eigenvalue.  The ratio of the eigenvalues have an interesting simplicity about them:

eValRatios = scaledVals/scaledVals[[3]]

This predicts that equilateral triangular molecules should have two fundamental frequencies, with ratio close to Sqrt[2].  If we had not added that extra bit of mass to one oxygen atom, the eigenvalue ratios would have been exact integers.

A digression into reality

Convert the eigenvalues to vibrational frequencies, in traditional units :

freqs = ws = 
Sqrt[10^28 scaledVals]/(2 Pi cInCmPerSec) // N  
(*reciprocal centimeters*)

The infrared spectrum of ozone can be nicely interpreted in terms of  three fundamental frequencies (Herzberg, IR and Raman, p. 286).  But only two fundamentals are provided by this model.  Therefore the model is wrong for ozone.

The two strongest Raman lines of cyclopropane (CH2)3 have frequencies of 1189 and 866 wavenumbers (Herzberg, ibid., p. 352).

(1189/866) // N

So their ratio is about 3% lower than the value Sqrt[2] predicted by a simple model in which the  CH2 groups move rigidly together like an atom.

The mode frequencies are the square roots of the eigenvalues.  We adopt a strange time unit equal to one cycle of the lower frequency vibration:

{w1,w2,w3,w4,w5,w6} = Sqrt[eValRatios]

We now write the solutions for the normal modes:

normalModes = Array[Q,6];

Qsolns = {
Q[1] -> A1 Cos[w1 t + p1],
Q[2] -> A2 Cos[w2 t + p2],
Q[3] -> A3 Cos[w3 t + p3],
Q[4] -> x0 + vx t,
Q[5] -> y0 + vy t,
Q[6] -> a0 + wz t};

The arbitrary constants of integration (initial conditions) that distinguish all possible motions from each other are the amplitudes (A1, A2, A3) and phases (p1, p2, p2) of the vibrational motions, and the initial positions (x0, y0, a0), and velocities (vx, vy, and wz) of the uniform motions.

Eigensystem does not automatically orthogonalize degenerate eigenvectors.  Test this out for yourself, multiplying SmatRawT by its own transpose (in both possible multiplication orders).  This non-orthogonality is contrary to one of our fundamental requirements on matrix S, so we must cure it with a Gram-Schmidt orthogonalization:

SmatT = GramSchmidt[SmatRawT]//Chop;

Smat = Transpose[SmatT];

Now we check that Smat has the properties we required of it.
(1) Is it unitary?

Smat.SmatT // Chop

Yes.  It is not necessary to test the other mutiplication order, but you may if you like.  What general theorem permits me to be confident that the other order will also produce a unit matrix?  

(2)  Does Smat diagonalize krm by a similarity transform?

SmatT.(10^-28 krm).Smat // Chop;
MatrixForm[%]

Yes.  All is well.

Step 8:  pleasing combinations of degenerate modes

There is nothing simple or meaningful about the degenerate orthogonal vectors created by  the Gram-Schmidt procedure.  If you want them in a form that is easy to visualize, you must recombine them by hand.

We illustrate this by recombining the three zero frequency modes in a meaningful way.  (Generally, the zero frequency modes are thrown out.  But it is good to see them dealt with explicitly in one simple example, so here they are.)  We lift out the last three eigenvectors, as returned by GramSchmidt :

basis0 = Take[SmatT,-3]

Legend has it that this little vector space should be spanned by an x-translation, a y-translation, and a z-rotation.  Let's see if the legend is true.

In Cartesian space, mode vectors representing x-translation amd y-translation would be

xTransCar = {1,0,1,0,1,0};
yTransCar = {0,1,0,1,0,1};

But our orthogonalized mode vectors are in rootmass weighted space, and in rootmass space the translation vectors are

xTransRM = Sqrt[massVec]*xTransCar // Normalize

yTransRM = Sqrt[massVec]*yTransCar // Normalize

These are obviously orthogonal.  But are they spanned by basis0?  To find out, dot them against all three eigenvectors:

xTCoeffs = Map[Dot[xTransRM,#]&,basis0]

Generally, it takes six basis vectors to span a vector with six components.  But if the legend is true, the three in basis0 should do the job:

xTCoeffs.basis0  // Chop

Yes!  Vector xTransRM is spanned by basis0.

yTCoeffs = Map[Dot[yTransRM,#]&,basis0]

yTCoeffs.basis0  // Chop

Yes again; vector yTransRM is also spanned.  Now all we have to do is find the vector in this space that is orthogonal to both xTransRM and yTransRM.  It should be the legendary rotational displacement mode.  

To get the third vector, we will use GramSchmidt on the two translation vectors, together with a third vector selected from basis0.  We will get the best accuracy by using one that is already as orthogonal as possible to both the translations.  Which one is this?

{{xTransRM.SmatT[[4]],yTransRM.SmatT[[4]]},
 {xTransRM.SmatT[[5]],yTransRM.SmatT[[5]]},
 {xTransRM.SmatT[[6]],yTransRM.SmatT[[6]]}}

It looks like SmatT[[4]] should be the best.

zeroModesNeat=
GramSchmidt[{xTransRM,yTransRM,SmatT[[4]]}] // Chop

The last vector could plausibly be parallel to an infinitesimal rotation about z; note that the peak atom (#1) is moving purely horizontally, which it does when the triangle undergoes an infinitesimal rotation about z.  We will plot all of these motions later.

Now join these neat zero frequency vectors onto the three vibration vectors, to make a basis set that spans all possible atomic displacements:

SmatNeatT = Join[Take[SmatT,3],zeroModesNeat]

SmatNeat = Transpose[SmatNeatT];

At this point the unitarity and the diagonalizing ability of SmatNeat should be checked.  We leave this as an exercise.

Step 9:  make the back-transforms to Cartesian

The definition of the normal modes Q is

so the inverse of the normal mode transform is

We already have the six Q(t) solutions; therefore we implement this as a replacement rule that will turn q's into expressions containing Q's:

qVec = Array[q,6];

q2Qgeneral = 
		qVec->SmatNeat.normalModes // Thread;
		Print[ColumnForm[%]]

We can discard the uniform modes by setting their amplitudes to zero:

q2Qvibr=
		q2Qgeneral/.{Q[4]->0,Q[5]->0,Q[6]->0}

Or we can do the opposite:

q2Quniform=
		q2Qgeneral/.{Q[1]->0,Q[2]->0,Q[3]->0}

Remember, the q's are still root-mass weighted.  Now we unweight the vibrational vectors to give a picture of nuclear motion in "real" x,y space.

Because the masses are so tiny (in grams), their inverses are quite large.  We therefore adjust the length scale by a factor of 10^-10, and again implement the transform as a replacement:

zVec = Array[z,6];
z2q = zVec-> 10^-10 qVec/Sqrt[massVec] // Thread;
ColumnForm[%]

In a molecule with different kinds of atoms in it, the numbers would of course vary more widely.

Step 10: the answer, in algebraic and graphic form

the answer in algebraic form

Now writing an explicit solution in Cartesian atomic motion coordinates is trivial.  We start with the symbols z that we want.  We replace each z in terms of the corresponding q; then q in terms of the Q's, then the Q's in terms of their harmonic (or uniform) solutions.  The most general vibrational motion is given by

Vibrational motions

zVibration = zVec/.z2q/.q2Qvibr/.Qsolns // Expand;
Print[%]

Uniform motions

The most general uniform motion is given by

zUniform = zVec/.z2q/.q2Quniform/.Qsolns // Expand;
Print[%]

the answer in graphic form

Vibrational motions

mode 1  (the breathing mode)

m1 = Partition[zVibration,2] /. {A2->0,A3->0}

m1SnapshotData =
restpoints + m1 /. {A1->.03,p1->0,t->0}

m1Snapshot = Show[restPicture,Graphics[
{PointSize[.04],Map[Point,m1SnapshotData]},stdPlot]];

mode 2  (asymmetric stretch of sides, degenerate with 3)

m2 = Partition[zVibration,2]/.{A1->0,A3->0}

m2data = restpoints + m2 /. {A2->.03,p2->0,t->0};
m2Snapshot = Show[restPicture,Graphics[
  {PointSize[.04],Map[Point,m2data]},stdPlot]];

mode 3  (bending mode for top angle, degenerate with 2)

m3 = Partition[zVibration,2]/.{A1->0,A2->0}

m3data = restpoints + m3 /. {A3->.03,p3->0,t->0};
m3Snapshot = Show[restPicture,Graphics[
  {PointSize[.04],Map[Point,m3data]},stdPlot]];

Show[GraphicsArray[
{{m1Snapshot,m2Snapshot,m3Snapshot}}]];

1260 wvnbr         891 wvnbr            891 wvnbr
                 mode 1															mode 2                 mode 3

Uniform motions

mode 4  (x-translation)

m4 = Partition[zUniform,2]/.
						{y0->0,vy->0,a0->0,wz->0}

m4data = restpoints + m4 /. {x0->0.03,t->0}

m4Snapshot = Show[restPicture, Graphics[
  {PointSize[.04],Map[Point,m4data]},stdPlot]];

mode 5  (y-translation)

m5 = Partition[zUniform,2]/.
						{x0->0,vx->0,a0->0,wz->0}

m5data = restpoints + m5 /. {y0->0.03,t->0}

m5Snapshot = Show[restPicture, Graphics[
  {PointSize[.04],Map[Point,m5data]},stdPlot]];

mode 6  (z-rotation)

m6 = Partition[zUniform,2] /. 
						{x0->0,y0->0,vx->0,vy->0}

m6data = restpoints + m6 /. {a0->0.03,t->0}

m6Snapshot = Show[restPicture, Graphics[
  {PointSize[.04],Map[Point,m6data]},stdPlot]];

Show[GraphicsArray[
{{m4Snapshot,m5Snapshot,m6Snapshot}}]];

x-translation          y-translation									z-rotation

Animations of these modes, and their linear combinations, are easily made in Mma .  See the notebook "Triangle mode movies" to see how it is done.