User Tools

Site Tools


Programming Project #1: Molecular Geometry Analysis

The purpose of this project is to introduce you to fundamental C-language programming techniques in the context of a scientific problem, viz. the calculation of the internal coordinates (bond lengths, bond angles, dihedral angles), moments of inertia, and rotational constants of a polyatomic molecule. A concise set of instructions for this project may be found here.

We thank Dr. Yukio Yamaguchi of the University of Georgia for the original version of this project.

Step 1: Read the Coordinate Data from Input

The input to the program is the set of Cartesian coordinates of the atoms (in bohr) and their associated atomic numbers. A sample molecule (acetaldehyde) to use as input to the program is:

6  0.000000000000     0.000000000000     0.000000000000
6  0.000000000000     0.000000000000     2.845112131228
8  1.899115961744     0.000000000000     4.139062527233
1 -1.894048308506     0.000000000000     3.747688672216
1  1.942500819960     0.000000000000    -0.701145981971
1 -1.007295466862    -1.669971842687    -0.705916966833
1 -1.007295466862     1.669971842687    -0.705916966833
(You may cut and paste the above or download the file from here.) The first line above is the number of atoms (an integer), while the remaining lines contain the z-values and x-, y-, and z-coordinates of each atom (one integer followed by three double-precision floating-point numbers).

After downloading the file to your computer (to a file called “geom.dat”, for example), you must open the file, read the data from each line into appropriate variables, and finally close the file.

  • Hint 1: Opening and closing the file stream
  • Hint 2: Reading the number of atoms
  • Hint 3: Storing the z-values and the coordinates

Step 2: Bond Lengths

Calculate the interatomic distances using the expression:

R_{ij}=\sqrt{(x_i-x_j)^2 + (y_i - y_j)^2 + (z_i - z_j)^2}

where x, y, and z are Cartesian coordinates and i and j denote atomic indices.

Step 3: Bond Angles

Calculate all possible bond angles. For example, the angle, \phi_{ijk}, between atoms i-j-k, where j is the central atom is given by:

\cos {\phi_{ijk}} = {\bf \vec{e}_{ji} } \cdot {\bf \vec{e}_{jk}}

where the {\bf \vec{e}_{ij}} are unit vectors between the atoms, e.g.,

        e_{ij}^x = - \left(x_i - x_j \right)/R_{ij},\ \ \ \ \ 
        e_{ij}^y = - \left(y_i - y_j \right)/R_{ij},\ \ \ \ \ 
        e_{ij}^z = - \left(z_i - z_j \right)/R_{ij}

  • Hint 1: Memory allocation for the unit vectors
  • Hint 2: Avoiding a divide-by-zero
  • Hint 3: Memory allocation for the bond angles
  • Hint 4: Smart printing
  • Hint 5: Trigonometric functions

Step 4: Out-of-Plane Angles

Calculate all possible out-of-plane angles. For example, the angle \theta_{ijkl} for atom i out of the plane containing atoms j-k-l (with k as the central atom, connected to i) is given by:

       \sin {\theta_{ijkl}} =  \frac{{\bf \vec{e}_{kj} \times \vec{e}_{kl} }}{ \sin {\phi_{jkl}}} \cdot  {\bf \vec{e}_{ki} }

Step 5: Torsion/Dihedral Angles

Calculate all possible torsional angles. For example, the torsional angle 
for the atom connectivity i-j-k-l is given by:

\cos {\tau_{ijkl}} = \frac{ ({\bf \vec{e}_{ij} \times \vec{e}_{jk} }) \cdot ( {\bf \vec{e}_{jk} \times \vec{e}_{kl} } ) }{  \sin {\phi_{ijk}} \  \ \sin{\phi_{jkl}}}

Can you also determine the sign of the torsional angle?

Step 6: Center-of-Mass Translation

Find the center of mass of the molecule:

        X_{c.m.} = \frac{ \sum_i m_i x_i}{\sum_i m_i},\ \ \ \ 
        Y_{c.m.} = \frac{ \sum_i m_i y_i}{\sum_i m_i},\ \ \ \ 
        Z_{c.m.} = \frac{ \sum_i m_i z_i}{\sum_i m_i},

where mi is the mass of atom i and the summation runs over all atoms in the molecule.

Translate the input coordinates of the molecule to the center-of-mass.

Step 7: Principal Moments of Inertia

Calculate elements of the moment of inertia tensor.


      I_{xx} = \sum_i m_i \left( y_i^2 + z_i^2 \right),\ \ \ \ 
      I_{yy} = \sum_i m_i \left( x_i^2 + z_i^2 \right),\ \ \ \
      I_{zz} = \sum_i m_i \left( x_i^2 + y_i^2 \right)


      I_{xy} = \sum_i m_i  x_i y_i,\ \ \ \ 
      I_{xz} = \sum_i m_i  x_i z_i,\ \ \ \ 
      I_{yz} = \sum_i m_i  y_i z_i.

Diagonalize the inertia tensor to obtain the principal moments of inertia:

      I_a \leq I_b \leq I_c

Report the moments of inertia in amu bohr2, amu Å2, and g cm2.

Based on the relative values of the principal moments, determine the molecular rotor type: linear, oblate, prolate, asymmetric.

Step 8: Rotational Constants

Compute the rotational constants in cm-1 and MHz:

\[ A = \frac{h}{8{\pi^2} c I_a} \] \[ B = \frac{h}{8\pi^2 c I_b} \]
\[ C = \frac{h}{8\pi^2 c I_c} \]

A \geq B \geq C

Test Cases


E.B. Wilson, J.C. Decius, and P.C. Cross, Molecular Vibrations, McGraw-Hill, 1955.

crawdad/programming/project1.txt · Last modified: 2014/08/29 10:10 by crawdad