Page 166 - Numerical methods for chemical engineering

152 3 Matrix eigenvalue analysis

angles at their natural values of θ a = 2π/3. Let θ αβγ be
(
r [αβ] · r [γβ] [αβ] [β] [α]
[α] [β] [γ ] r = r − r
θ αβγ r ,r ,r = acos [γβ] [β] [γ ] (3.275)
r
r
[αβ] [γβ] r = r − r
We resist bending of this angle away from its natural value θ a = 2π/3 by adding to the
potential energy an angle-bending contribution
[a] [α] [β] [γ ] 1 [α] [β] [γ ] 2
U r ,r ,r = K a θ αβγ r ,r ,r (3.276)
αβγ 2 − θ a
For each angle in the energy model for this 2-D molecule, we use
K a = 1 θ a = 2π/3 (3.277)

For a real molecule, we would include additional terms from van der Waal’s and electro-
static nonbonded interactions, but here for simplicity we consider only the bonded energy
contributions. Let the coordinate vector be
q = [x 1 y 1 x 2 y 2 x 3 y 3 x 4 y 4 x 5 y 5 x 6 y 6 ] T (3.278)

The potential energy of the molecule is
[b] [b] [b] [b] [b]
U(q) = U + U + U + U + U
12 13 14 25 26
[a] [a] [a] [a] [a] [a]
+ U + U + U + U + U + U (3.279)
213 214 314 125 126 526
and the kinetic energy is

m 1 ˙r m 1 ˙r m 2 ˙r m 2 ˙r m 2 ˙r m 2 ˙r
[2]
[5]
[6]
[3]
[4]
[1]
K(q, ˙q) = + + + + + (3.280)
2 2 2 2 2 2
where m 1 = 10 and m 2 = 1. We can then compute the mass matrix M such that
1 T
K(q, ˙q) = ˙q M ˙q (3.281)
2
Your task is to write a MATLAB program that computes the vibrational frequencies, through
the following steps:
(a) Write a routine that sets the atomic positions, with r [1] = 0, such that all bond lengths
and angles take their natural values. This is a minimum potential energy state. Store
the results as a state vector ˆ q.
(b) Write a routine that returns, for input q, the potential energy U(q).
(c) Using the routine from (b), write a routine to compute the elements of the Hessian
matrix
2
∂ U
[H(ˆq)] mn = [H(ˆq)] nm = (3.282)

∂q m q n ˆq
This can be done numerically by finite differences,

∂ ∂U (∂U/∂q n )| ˆq+εe [m] − (∂U/∂q n )| ˆq−εe [m]
[H(ˆq)] mn = ≈ (3.283)
∂q m ∂q n ˆ q 2ε