Molecular Vibrations and Group Theory.

The author describes a package for obtaining the vibrational modes of non-linear polyatomic molecules. The computed normal modes are visualized through Mathematica—generated animations.

by Gustavo Demarco

This paper describes a purely classical discussion of the vibrations of molecules (regarded as systems of several interacting nuclei). Such a description is necessary prior to any quantum mechanical investigation of this problem.

Applying the methods of group theory enables us to carry out a classification of the normal modes of vibration from the knowledge of the symmetry group of the molecule with the nuclei in its equilibrium positions. Then, the utilization of projection operators [1] allow us to solve the more complex problem of determining the normal coordinates or, what is the same, the displacements of the nuclei in each normal mode. This allows a beautiful visualization of the vibrating molecule with Mathematica generated animations. The functions necessary for doing these computations are bundled in the package Molecules.m. 

In section 1, I summarize the group theoretic considerations necessary for the problem of molecular vibrations. The needed information about several point groups consists of a set of data files (one for each group) generated with previous packages charactr.m [2] and Projectors.m [3]. In section 2, I explain the use of these data files. Section 3 shows in more detail the utilization of the package Molecules.m solving the problem of vibrations in the  (Subscript[NH, 3]) molecule. Section 4 contains the strictly necessary steps to solve any molecule vibration problem using (Subscript[UF, 6]) as example. Section 5 contains a discussion for the case when there are involved complex irreducible representations and section 6 gives a brief discussion about the force matrix. Section 7 contains the conclusions.

The largest molecule we considered was the fullerene (Subscript[C, 60]). In this case, although the computations can be performed in symbolic form, we modified some of the functions of the package to do the calculations with floating-point arithmetic to avoid time consuming operations. The symbolic results are such large expressions that no insight is gained doing the computations in symbolic form. Also, the visualization of the vibrating (Subscript[C, 60]) molecule is done as a polyhedra with faces and not as a collection of points joined by bonds (as we have done in the other examples) so we have also adapted the visualization function for this case.

In the electronic distribution I included the three files (the notebook vib_C60num.nb and the data files C60.dat.num and bases.dat.num) necessary to perform the animation of the (Subscript[C, 60]) fullerene (buckyball). Also included in the electronic distribution are the data files for the 32 point groups and a notebook examples.nb with several examples:  (Subscript[C, 2]) (Subscript[H, 4]), (Subscript[C, 2]) (Subscript[H, 6]), (Subscript[C, 6]) (Subscript[H, 6]), CH(Subscript[Cl, 3]), (Subscript[H, 2]) O and (Subscript[OsF, 8]).

As a final remark, I should mention that throughout the paper, I start a fresh Mathematica session in each section.

2.  Group Theoretic Considerations.

In general, a system of N nuclei have 3N degrees of freedom. Of these, three correspond to a translational motion and three (two for a linear molecule) to a rotational motion of the system as a whole. This leaves 3N – 6 pure vibrational degrees of freedom (3N-5 for a linear molecule).

The energy of the system executing small vibrations is a quadratic form in the velocities and in the cartesian coordinates of the particles. We can always find a set of coordinates (Subscript[Q, (Subscript[α, i])]) (normal coordinates) such that the energy takes a diagonal form,

\(E\  = \ 
      1\/2\ \(∑\+\(i, α\)\ Q\& . \_\(α\_i\)\%2\)\  + \ 
        1\/2\ \ \(∑\+α\ ω\_α\%2\ \ \(∑\+i\ \ Q\_\(α\_i\)\%2\)\)\)

where (Subscript[ω, α]) are the frequencies of vibration. The sum over i =1,2,...,(Subscript[f, α]) means that (Subscript[ω, α]) has associated (Subscript[f, α]) normal coordinates.

From symmetry considerations [4], it can be shown that the (Subscript[f, α]) normal coordinates belonging to (Subscript[ω, α]) give some irreducible representation (IR) of the symmetry group of the molecule. To find the IR's present in Equation 1, we need the reducible representation (RR) given by all the coordinates together (we call this RR the total representation TR). Decomposing it into irreducible parts, we determine the symmetry properties of the vibrations. We are interested only in the purely vibrational degrees of freedom. Thus we must exclude from the TR the six degrees of freedom corresponding to pure translations and rotations of the molecule as a whole. These six modes have null frequencies. We call the RR obtained by excluding these six modes the total vibrational representation (TVR).

To find the IR's contained in the TVR we must compute the character of this representation. Then we project it onto the characters of the IR's of the symmetry group G of the molecule using the orthogonality relation between these [5].
  
The character of the TVR for a symmetry operation T can be written in the following form,

\(\(χ\^TVR\) \((T)\)\  = \ \(\(U\_M\) \((T)\)\(\ \)\(×\)\(\ \)\(U\_G\) \((T)\)\(\ \)\)\)

where (Subscript[U, M]) (T) depends on T and on the distribution of the atoms in the molecule and  (Subscript[U, G]) (T) depends only on the operation T.
  The symmetry transformation T can only be one of the following: pure rotation R, rotorreflexion S and inversion I. It can be shown [4] that  (Subscript[U, M]) (T)  can take the following values,

\(\(\(\(\(\(\(a\)\()\)\)\ \ \ \(U\_M\) \((R)\)\  = N\_R - 2

                b\)\()\)\)\ \ \ \(U\_M\) \((S)\)\  = \ N\_S

          c\)\()\)\)\ \ \ \(U\_M\) \((I)\)\  = \ N\_I\)

where  (Subscript[N, i]) is the number of atoms that do not change their position under the symmetry operation i. It can also be shown [4] that,

\(\(\(\(\(\(\(a\)\()\)\)\ \ \ \(U\_G\) \((R)\)\  = \ 
                  1 + 2\ Cos[Φ]

                b\)\()\)\)\ \ \ \(U\_G\) \((S)\)\  = \ \(-1\) + 
              2\ Cos[Φ]

          c\)\()\)\)\ \ \ \(U\_G\) \((I)\)\  = \ \(-3\)\)

where Φ is the angle of rotation in the operation R or S (remember that any rotorreflexion S can be expressed as S = I·R(Φ) for a suitable chossen pure rotation R(Φ)).
  The numbers  (Subscript[U, G]) (T) depend only on G and are included in the data files for the point groups with the constant name vibrchar.
  Now we consider the problem of determining the displacements of the nuclei in each normal mode (a detailed discussion of this point can be viewed from reference [1]). First we should obtain the TR. This is a 3N dimensional representation acting on vectors of the form,

v = {(Subscript[x, 1]), (Subscript[y, 1]), (Subscript[z, 1]), (Subscript[x, 2]), (Subscript[y, 2]), (Subscript[z, 2]),  ... , (Subscript[x, N]), (Subscript[y, N]), 
        (Subscript[z, N])}

where (Subscript[x, i]), (Subscript[y, i]), (Subscript[z, i]) are the cartesian coordinates of the i nuclei. This TR of G is simply the direct product representation obtained from the three dimensional representation of G (that is, the representation of rotation matrices in (ℝ^3)) and the representation obtained by permuting the nuclei for each group operation, that is, an N dimensional permutation of G. For a given symmetry operation T we denote by (Γ^TR)(T) the matrix that represents this operation acting on the vectors of Equation (5), that is, (Γ^TR)(T) are the matrices in TR. Now suppose that a given IR (Γ^((ν))) of G is contained (Subscript[a, ν]) times in TR. Thus, in the space (ℝ^(3  N)) on which TR acts, there are (Subscript[a, ν]) linearly independent vectors  (Subscript[v, j]) (j=1,...,(Subscript[a, ν])) that transform as the first row of (Γ^((ν))). To obtain these (Subscript[a, ν]) vectors we construct the following projector matrices [1],

(Subscript[℘, 11]^((ν)))=  ((Subscript[n, ν])/g)(∑_(T ∈ G))\(\([\(Γ\^\((ν)\)\) \((T)\)]\)\_11\%*\)   (Γ^TR)(T)  .

where (Subscript[n, ν])is the dimension of the matrices (Γ^((ν)))(T) in the IR (Γ^((ν))). g is the number of elements in G.
  Thus, (Subscript[℘, 11]^((ν))) is a 3N×3N matrix. This matrix has  (Subscript[a, ν])  linearly independent rows which can be taken as the (Subscript[a, ν]) linearly independent vectors (Subscript[v, j]). Each of these vectors transforms as the first row of (Γ^((ν))), and each of them allow us to generate the (Subscript[a, ν]) copies of (Γ^((ν))) contained in (Γ^TR). However, since we are interested in the pure vibrational modes, we have to eliminate from the set {(Subscript[v, j])}_(j = 1,  ... , (Subscript[a, ν])) those vectors that correspond to the zero frequency modes (pure rotations and translations). We start eliminating the translations. For each vector (Subscript[v, j]) we extract the components (Subscript[u, a])={(Subscript[x, a]),(Subscript[y, a]),(Subscript[z, a])} for each atom a (see Equation 5). In Mathematica notation,

{(Subscript[u, 1]),  ... , (Subscript[u, N])} = Partition[(Subscript[v, j]), 3];

These (Subscript[u, a]), when multiplied globally by an infinitesimal constant, are the infinitesimal displacements that transform as the first row of one of the modes ⋁.
  The displacement of the centre of mass of the molecule corresponding to the displacement vectors (Subscript[v, j]) is ,

\(u\_CM = \(1\/M\) \(∑\+a m\_a\ u\_a\)\ , \ \ \ \ \ \ \ \ M = ∑\+a m\_a\)

Now we substract (Subscript[u, CM]) to each (Subscript[u, a]) ,and reconstruct from them the modified vectors (Subscript[v, j]) that do not contain translations,

\(\(u\_a\%'\) = u\_a - u\_CM\)

\(\(\(v\_j\%'\) = 
        Flatten[{\(u\_1\%'\),  ... , \(u\_N\%'\)}];\)\)

Some of these new vectors {\(\(v\_j\%'\)\)}_(j = 1,  ... , (Subscript[a, ν])) will be null vectors and correspond to pure translations. We discard them. 
  Now we proceed with the global rotations. For simplicity, let us suppose that the centre of mass coincides with the origin of coordinates, that is, (Subscript[R, CM])= 0. The "angular momentum" (strictly speaking, an infinitesimal rotation due to the fact that the (Subscript[u, a]) are infinitesimal displacements, not velocities) of the molecule to first order in (Subscript[u, a]) is,

\(L = ∑\+a m\_a\ r\_a\%0×\ \(u\_a\%'\)\)

where (Subscript[r, a]^0) are the equilibrium positions of the nuclei.
  If we perform a global infinitesimal rotation Ω (an infinitesimal rotation can be caracterized by vector with direction given by the axis of rotation and modulus equal to the angle of rotation), L changes to \(\(L\^'\)\),

\(\(L\^'\) = ∑\+a m\_a\ r\_a\%0×\ \(u\_a\%''\)\)

where \(\(\(\ \)\(u\_a\%''\)\)\)=  \(\(u\_a\%'\)\)- (Subscript[h, a])  and  (Subscript[h, a]) = (Subscript[r, a]^0)× Ω .
  We have to find Ω such that,

\(\(L\^'\) = 0\)

The answer is, in component notation,

(Subscript[Ω, j]) = (Subscript[(((U^(-1)))), ji]) (Subscript[L, i])

where

\(U\_im = ∑\+a m\_a[\(\((\ 
                    r\_a\%0\ )\)\_m\) \((\ r\_a\%0\ )\)\_i - \(δ\_\(\(im\)\(\ \)\)\) \((\ 
                r\_a\%0\ ·r\_a\%0)\)]\)

Now, from the vectors \(\(\(\ \)\(u\_a\%''\)\)\) we obtain the vectors \(\(\(\ \)\(v\_j\%''\)\)\),

\(\(\(v\_j\%''\) = 
        Flatten[{\(u\_1\%''\),  ... , \(u\_N\%''\)}];\)\)

From this set, the non-null vectors are the pure vibrational modes that transform as the first row of the IR (Γ^((ν))). Let  (Subscript[b, ν]) be the number of non-null vectors and rename them as  (Subscript[v, j]). Now we construct the following (Subscript[n, ν])-1 projection matrices [1],

(Subscript[℘, 1  m]^((ν)))=  ((Subscript[n, ν])/g)(∑_(T ∈ G))\(\([\(Γ\^\((ν)\)\) \((T)\)]\)\_\(1  
                  m\)\%*\)   (Γ^TR)(T)  .

From each vector (Subscript[v, j]) we obtain the following  (Subscript[n, ν])-1 vectors,

\(v\_j\%α = \(v\_j·℘\_\(1  α\)\%\((ν)\)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ α = 2\),  ... , n\_ν\)

where (Subscript[v, j])≡ (Subscript[v, j]^1)  is taken as a row matrix. 
  The (Subscript[n, ν]) vectors {(Subscript[v, j]^α)}_(α = 1,  ... , (Subscript[n, ν])) transform as the IR (Γ^((ν))). Thus, from the (Subscript[b, ν]) vectors {(Subscript[v, j])}_(j = 1,  ... , (Subscript[b, ν])) we obtain the (Subscript[b, ν]) copies of (Γ^((ν))) contained in (Γ^TVR), that is, we have (Subscript[b, ν]) linearly independent basis {(Subscript[v, j]^1),...,(Subscript[v, j]^((Subscript[n, ν])))}_(j = 1,  ... , (Subscript[b, ν])) for the IR (Γ^((ν))). However, in general, these do not correspond to the pure vibrational normal modes {(Subscript[Q, (Subscript[α, 1])]),...,(Subscript[Q, (Subscript[α, (Subscript[n, ν])])])}_(α = 1,  ... , (Subscript[b, ν])) (see Equation 1) but to a linear combination of them. This problem cannot be solved from purely group theoretic considerations. Instead, we have to take into account the force matrix, that is, an explicit knowledge from the dynamics is necessary. A few considerations about these points are relegated to section 7.

3.  The Point Groups Database.

With the help of our previous packages charactr.m [2] and Projectors.m [3], it is easy to obtain all the necessary information needed to treat molecular vibration problems. For several  point groups, we generated a file <group>.dat, where <group> stands for a standard name of the group. When one of these files is loaded in Mathematica, it defines the following constants,

  comments   
  		this is a string that tells what rotation operations have been used to generate the group. This is necessary 
  		to give the location for the atoms that generate the molecule from the group operations.
  group		
  		the representation of the group by rotation matrices in three diemensions.
  pos
  		a list whose elements are the classes of the group given by the position of their elements within the 
  		group.
  g
  		the number of elements within each class.
  ch
  		the character table.
  vibrchar
  		we have already defined this quantity in the previous section (see Equation 4 and the comments 
  		thereafter).
  symbols
  		the names for the IR's following Landau notation [4].
  reality
  		each element from this list gives the reality properties of each IR. Its elements can take the following 
  		values:
  		1	 the IR is equivalent to as real representation (real representation).
  		0	the IR is complex and not equivalent to the complex conjugate representation (complex 
  			representation).
  		-1	the IR is equivalent to the complex conjugate representation but it is not equivalent to a real one 
  			(pseudo real representation).
  irreps
  		the matrices for the IR's. All the representations but the five dimensional representations for the 
  		icosahedral groups are unitary representations. The five dimensional representations for the 
  		icosahedral groups are equivalent to unitary representations, but the computation is very time 
  		consuming (at least in our old 486).

4. Vibrations of the ammonium molecule ((Subscript[NH, 3])).

This molecule is a regular tetrahedron, with the N atom at the vertex and the H atoms at the corners of the base. The symmetry group of this molecule is (Subscript[C, 3  v]). For any molecule with symmetry group G, the atoms can be grouped in disjoints sets (orbits). All the atoms within an orbit have positions that can be generated from the position of any one of them by applying to it all the symmetry operations of the group G. Thus, each orbit can be generated from one atom and the complete molecule can be generated from a minimal set of atoms (generators), one from each orbit. In the present case, the complete molecule can be obtained from the N atom and one of the H atoms. The orbit of the N atom contains only this atom and the orbit of the H atom contains the three H atoms. The comments in the data file C3v.dat helps us to choose an appropiate location for these generating atoms.

<< "Molecules.m";
<< "PointGroups/C3v.dat";
comments

We put the N atom (atomic mass = 14) at {0,0,1} and an H generating atom (atomic mass = 1) at {0,1,0},

locations = {{0, 0, 1}, {0, 1, 0}};
masses = {14, 1};
cls = {GrayLevel[0], RGBColor[0, 1, 0]};

where cls are the colors we assign to the generating atoms  in the animations.
  Now we generate the coordinates, masses and colors of the remaining atoms. All the atoms in each orbit have the same mass and we assign the same color to all the atoms within each orbit,

mol = GenerateMolecule[group, locations]

masses = GenerateMasses[group, mol, locations, masses]

cls = GenerateMasses[group, mol, locations, cls]

Before classifying the nomal modes and obtaining the displacement vector for each normal mode, we will verify that we have correctly generated the molecule with a 3D plot. First we generate bonds joining the atoms from two generating bonds given by the atoms at the end of the bonds,

bg = {{1, 2}, {2, 3}};

bonds = GenerateBonds[group, bg, mol]

bnscoor = Map[mol[[#]] &, bonds];

\(\(Show[Graphics3D[Map[Line, bnscoor]], 
        Graphics3D[{PointSize[0.05], 
			Map[Point, mol]}]];\)\)

Now we compute the vibrational modes,

modes = NormalVibrations[group, pos, ch, vibrchar, mol]

We can see a standard notation for the modes,

symbols . modes

The "\!\(A\_1\)" are symmetrical modes (non-degenerate) and the E modes are double-degenerate modes.
  Now we should obtain (Γ^(TR ))(totalrep). As we have said in section 2, this representation is the direct product representation of the representation obtained by permuting the atoms with the symmetry operations (permrep) and the representation given by the rotation matrices in (ℝ^3) (group). In fact, permrep is a direct product of other permutation representations, one for each orbit,

Map[
      MatrixForm, (permrep = 
          Map[AtomsPermutations[mol, #] &, group])]

totalrep = MapThread[DirectProduct, {permrep, group}];

This is a 12×12 representation. 
     Now we proceed to compute the set of (Subscript[a, ν]) linearly indepent rows of the projector  (Subscript[℘, 11]^((ν))), or what is the same, the vectors that transform as the first row (see Equation 6) for each of the IR's present in  (Γ^(TVR ))(see the variable modes before).  We start with (Γ^((1))),

vecs1 = FirstRowVectors[irreps[[1]], totalrep]

We can see that a linear combination of these three vectors should correspond to a pure translation or rotation since there are only two pure vibrational modes for (Γ^((1))) (modes[[1]]=2). To eliminate this zero frequency mode we start fixing the centre of mass of the molecule, that is, we compute the vectors \(\(v\_j\%'\)\) (see Equation 10),

vecs1 = FixTranslations[vecs1, masses]

There is no rotational mode remaining in vecs1. Since the (Subscript[A, 1]) modes are one-dimensional, there are no additional vectors in these modes (the matrices in (Γ^((1))) have only one row, that is they are numbers).
  Now we proceed with (Γ^((3))), the other pure vibrational modes,

vecs3 = FirstRowVectors[irreps[[3]], totalrep]

vecs3 = FixTranslations[vecs3, masses]

Thus, a linear combination of these vectors corresponds to a pure rotation, and we should eliminate it,

vecs3 = FixRotations[vecs3, masses, mol]

These are the vectors that transform as the first row of (Γ^((3))) in (Γ^TVR), but (Γ^((3)))is two dimensional, thus we should compute the correspondig vectors that transform as the second row of  (Γ^((3))) (each E vibrational mode has a two-fold degenaration), that is, we should compute the vectors (Subscript[v, j]^E), j=2 (see Equation 18).

base31 = 
      GenerateBasisVectors[vecs3[[1]], irreps[[3]], 
        totalrep]

base32 = 
      GenerateBasisVectors[vecs3[[2]], irreps[[3]], 
        totalrep]

Now from each of the vectors in vecs1, base31 and base32 we can make an animation for the vibrating molecule that corresponds to each of the vibrational modes. Remember however that these are not the pure vibrational modes but a linear combination of them. To obtain them we should give values to the force matrix and find the linear combination of modes that diagonalizes the force matrix (see Section 7).
  Let us see, for example, the molecule vibrating in the E mode that corresponds to base32[[2]].

ShowVibratingMolecule[mol, base32[[1]], bonds, NumPic -> 2, 
        Colors -> cls];

The option NumPic is the number of distincts pictures generated. The animation contains a series of NumPic pictures and  the reversed sequence to give continuity to the movement. You can also control the amplitude of the movement and the viewpoint.

5. Uranium Hexafluoride.

In the (Subscript[UF, 6]) molecule, the U atom is at the centre of an octahedron with the F atoms at the vertices. The symmetry group of the molecule is (Subscript[O, h]). We will solve the problem of molecular vibrations with the minimal set of necessary functions from the package Molecules.m. For simplicity we will assign atomic mass 2 to the U atom and 1 to the F atoms. (The U atom is much heavier than the F atoms thus the displacement for the U atom in each normal mode is much smaller than the displacements for the F atoms. Using the real atomic masses complicates the visualization of the animations.) 
  The generators for the molecule are the U atom and one of the F atoms. We put the U atom at the origin and an F atom at {0,0,1}.

<< "Molecules.m";
<< "PointGroups/Oh.dat";
comments

locations = {{0, 0, 0}, {0, 0, 1}};
masses = {2, 1};
cls = {GrayLevel[0], RGBColor[1, 0, 0]};
mol = GenerateMolecule[group, locations]

masses = GenerateMasses[group, mol, locations, masses]

cls = GenerateMasses[group, mol, locations, cls];

The normal modes are,

modes = NormalVibrations[group, pos, ch, vibrchar, mol]

modes . symbols

For visualization purposes we put bonds only between the F atoms,

bg = {{5, 6}};
bonds = GenerateBonds[group, bg, mol]

Now we compute (Γ^TR),

totalrep = VibrationalRepresentation[group, mol];

Now we obtain all the basis {(Subscript[v, j]^1),...,(Subscript[v, j]^((Subscript[n, ν])))}_(j = 1,  ... , (Subscript[b, ν])) (see Equation 18) for the pure vibrational modes for each (Γ^((ν))),

bases = 
        GenerateVibrationalBasis[irreps, modes, totalrep, mol, 
          masses];

bases is a list with length equal to the number of IR's of the symmetry group G. The element bases[[ν]] is a list with (Subscript[b, ν]) elements, where (Subscript[b, ν]) is the number of independent pure vibrational modes that transorm as (Γ^((ν))). Finally, bases[[ν, j]] is the list of the (Subscript[n, ν]) vectors within each independent vibrational mode in correspondence with the rows of (Γ^((ν))) and with the same energy. For example, there are two modes (Subscript[((Subscript[F, 1])), u]) (symbols[[10]]=(Subscript[((Subscript[F, 1])), u])), each of them containing three vibrations with the same energy ((Subscript[n, (Subscript[((Subscript[F, 1])), u])])= 3). Thus, the animation for the vector that transforms as the second row in bases[[10, 2]] should be obtained by executing the following comand,

ShowVibratingMolecule[mol, bases[[10, 2, 2]], bonds, 
        Colors -> cls];

The modes (Subscript[((Subscript[A, 1])), g]), (Subscript[E, g]), (Subscript[((Subscript[F, 2])), g]) and (Subscript[((Subscript[F, 2])), u]) appear only one time in (Γ^TVR). Thus, the basis vectors we have obtained for these modes are the pure vibrational modes; that is, we do not need to consider the force matrix to obtain them. This is not the case for the mode (Subscript[((Subscript[F, 1])), u]) that appears two times in (Γ^TVR).

6. Complex Representations.

The symmetry groups of the two previous examples have only real representations. However, there are several point groups that have complex irreducible representations. We will not treat the case of pseudo-real representations, because none of the 32 point groups have this type of representation, although the treatment would be very similar to the case of complex representations. In this case, there are some displacement vectors obtained by the projector method that are complex. Obviously displacement vectors should be real. One must take into account another symmetry present in vibrating molecules, that is, the time reversal symmetry. Here,  group theory tells us that there are extra degeneracies when the symmetry group has complex or pseudo-real representations. In the case of a normal mode transforming as a complex IR of dimension (Subscript[n, ν]), there is another normal mode transforming as the non-equivalent conjugate IR with the same energy. Thus, there are 2(Subscript[n, ν]) real, linearly independent combinations of these two normal modes that must be regarded as forming a single "physically irreducible" (real) [4] representation of twice the dimension (Subscript[n, ν]). In the case of pseudo real representation, the IR and the complex conjugate one are equivalent but they are not equivalent to a real one. In this case, if there is a normal mode with the symmetry of this IR, there should be another normal mode transforming with the same IR and with the same energy such that the two normal modes form a physically irreducible representation.
 As an example we will consider a ficticious (Subscript[H, 3]) molecule with the H atoms at the vertex of an equilateral triangle. The symmetry group of this molecule is (Subscript[D, 3  h]) and has two normal modes, one one-dimensional with symmetry \(\(A\_1\%'\)\) and the other, bidimensional, with symmetry "\!\(\(E\^'\)\)". (Subscript[D, 3  h]) has only real representations. So we will forget some of the symmetries of the molecule and consider (Subscript[C, 3]) as the symmetry group. This group has two one-dimensional complex IR's. As we will see, time-reversal symmetry forces the same degeneracy as in the case of (Subscript[D, 3  h]). It is possible to consider molecules that have only the (Subscript[C, 3]) symmetry, but they have at least six atoms and so the displacement vectors are eighteen dimensional. For pedagogical purpuses, it is sufficient to consider the  (Subscript[H, 3]) molecule with (Subscript[C, 3]) symmetry,

<< "Molecules.m";
<< "PointGroups/C3.dat";

points = {{0, 1, 0}};
mol = GenerateMolecule[group, points]

masses = GenerateMasses[group, mol, points, {1}];
cls = GenerateMasses[group, mol, points, {GrayLevel[0]}];
bonds = GenerateBonds[group, {{1, 2}}, mol]

modes = NormalVibrations[group, pos, ch, vibrchar, mol]

symbols

The types of the IR's are given by the quantity reality,

reality

Thus \(\(E\_+\)\) and "\!\(\(E\_-\)\)" are complex conjugate representations and the modes transforming as these IR's have the same energy, according to time-reversal symmetry. Now we compute the displacement vectors for each normal mode,

vibrep = VibrationalRepresentation[group, mol];
basis = 
        GenerateVibrationalBasis[irreps, modes, vibrep, mol, 
          masses];
basis = ComplexExpand[basis]

As we can see, the basis vectors for  \(\(E\_+\)\) and "\!\(\(E\_-\)\)" are complex and conjugate to each other and together form a real, two dimensional physically irreducible representation.
  From a complex normalized basis vector (Subscript[v, j]^α) (where α indicates an IR, j the row and we do not indicate the index that tells to which of the (Subscript[b, α]) copies of (Γ^((α))) present in (Γ^TVR)belongs this vector) like basis[[2,1,1]] and the complex conjugate ((Subscript[v, j]^α)\(TraditionalForm\`\()\^*\)\) (basis[[3,1,1]]) we can construct two real normalized independent vectors [1],

(Subscript[a, j]^α)=(1/(Sqrt[2])){(Subscript[v, j]^α)+((Subscript[v, j]^α)\()\^*\)}
(Subscript[d, j]^α)=(1/(i(Sqrt[2]))){(Subscript[v, j]^α)-((Subscript[v, j]^α)\()\^*\)}

Thus (Subscript[a, j]^α) and (Subscript[d, j]^α) do not transform as the j row of  (Γ^((α))), but rather they are basis vectors of the physically irreducible representation,

(({{(2^(-1/2))1, -((i2)^(-1/2))1}, {(2^(-1/2))1, ((i2)^(-1/2))1}}))(({{(Γ^((α))), 0}, {0, \(\(Γ\^\((α)\)\)\^*\)}}))(({{(2^(-1/2))1, -((i2)^(-1/2))1}, {(2^(-1/2))1, ((i2)^(-1/2))1}}))

which is a representation contained in (Γ^TVR) and where 1 is the identity matrix of dimension (Subscript[n, ν]).
  The function RealVibrationalBasis allows us to compute the real basis vectors. The second argument is a list whose elements are lists of two integers that give the positions in the output from GenerateBasisVectors of the complex representations and their conjugate representations. In the present case, the second and thirth IR's are conjugate to each other, so the second argument should be {{2,3}}. The output gives the (Subscript[a, j]^α) vectors in the position formerly ocupied by the second IR, and the (Subscript[d, j]^α) in the position occupied by the third IR,

realbasis = RealVibrationalBasis[basis, {{2, 3}}]

Now we could perform the animations with the following commands,

ShowVibratingMolecule[mol, realbasis[[1, 1, 1]], bonds, 
        Colors -> cls];

ShowVibratingMolecule[mol, realbasis[[2, 1, 1]], bonds, 
        Colors -> cls];

ShowVibratingMolecule[mol, realbasis[[3, 1, 1]], bonds, 
        Colors -> cls];

7. The Force Matrix.

The potential energy of a molecule executing small vibrations in terms of the cartesian coordinates of the individual atoms is of the form,

V = (v^T) F v

where v is the vector of Equation 5 considered as a column matrix and F is the force matrix. F is a real, symmetrical matrix with positive (pure vibrations) and zero (translations and rotations) eigenvalues. The number of independent parameters in F is determined by the symmetry properties of the molecule. F acquires a particularly (although in general not real if there are complex IR's) simple form when considered in a basis of displacement vectors that transform as the rows of the IR's that we have calculated in the previous sections. (Remember that these may not be the pure normal modes but a linear combination of them). If a given IR (Γ^((ν))) of dimension (Subscript[n, ν]) appears (Subscript[b, ν]) times in (Γ^TVR), we will denote the basis vectors for the  m copy (m=1,...,(Subscript[b, ν])) with (Subscript[v, mj]^ν) where j indicates the row of (Γ^((ν))) according to which (Subscript[v, mj]^ν) transforms.
  Group theory tells us that the matrix elements of F in this basis are of the form,

\(⟨v\_mj\%ν | F | 
          v\_nk\%μ⟩ = δ\_νμ\ δ\_jk\ F\_mn\%\((ν)\)\)

that is, F has non-null matrix elements only between vectors that transform with the same row in the same IR, and the value does not depend on the row.
  To obtain the energy of vibrations (eigenvalues of F) and the normal modes (eigenvectors of F) we have to diagonalize the force matrix. Thus, when the given IR (Γ^((ν))) appears in (Γ^TVR), we have to consider only the submatrix of F whose matrix elements are,

(Subscript[f, mn]^((ν))) = ⟨(Subscript[v, m1]^ν) | F | 
            (Subscript[v, n1]^ν)⟩ =  
        (Subscript[F, mn]^((ν)))                      m, 
    n = 1,  ... , (Subscript[b, ν])

This  (Subscript[b, ν])× (Subscript[b, ν]) matrix appears (Subscript[n, ν])times in the diagonal of F, due to Equation 22. Thus the eigenvalues of (f^((ν)))have multiplicity (Subscript[n, ν]) in F. The eigenvectors determine the coefficients for the linear combinations of vectors (Subscript[v, m1]^ν) that are the pure normal modes of the molecule that transform as the first row of (Γ^((ν))) (if they are complex we should make the unitary transformation of the previous section). These coefficients are the same for the combinations of the other rows (Subscript[v, mj]^ν). The number of independent coefficients in (f^((ν))) depends on the reality of (Γ^((ν))). If (Γ^((ν))) is real, (f^((ν))) is a real symmetrical matrix. If (Γ^((ν))) is complex, (f^((ν))) is hermition and the IR \(TraditionalForm\`\(\(Γ\^\(\((ν)\)\(\ \)\)\)\^*\)\), due to time reversal symmetry, has \(\(\(f\^\((ν)\)\)\^*\)\) (the conjugate of (f^((ν)))) as the force submatrix.
  As a simple example, consider a fictitious molecule like the one of the previous section but with an additional H atom at the centre of the triangle. Again we consider the symmetry of the molecule as being (Subscript[C, 3]),

<< "Molecules.m";
<< "PointGroups/C3.dat";

points = {{0, 1, 0}, {0, 0, 0}};
masses = {1, 1};
mol = GenerateMolecule[group, points]

masses = GenerateMasses[group, mol, points, masses];
modes = NormalVibrations[group, pos, ch, vibrchar, mol];
modes . symbols

Thus, the F matrix has three 2×2 (f^((ν))) submatrices, one for each IR (in this case, since all the IR's are one dimensional, the submatrices (f^((ν)))appears only one time on the diagonal, that is, there are only (Subscript[v, m1]^ν) vectors because (Subscript[n, ν])=1 for all the IR's).  We call these three submatrices as (f^A), \(f\^\(E\_-\)\)and \(f\^\(E\_+\)\). Due to time reversal, \(f\^\(E\_+\)\) is the conjugate of \(f\^\(E\_-\)\) and in general they are complex. (f^A) is real because A is a real IR. 
  We will find the pure normal modes given the coefficients of the (f^((ν))) matrices. We should first compute the displacement vectors as in the previous sections,

vibrep = VibrationalRepresentation[group, mol];
basis = 
        GenerateVibrationalBasis[irreps, modes, vibrep, mol, 
          masses];
basis = ComplexExpand[basis]

We start with the A representation. The vectors (Subscript[v, 11]^A) and (Subscript[v, 21]^A) are,

vA11 = basis[[1, 1, 1]]

vA21 = basis[[1, 2, 1]]

The most general form for the (f^A)matrix is,

MatrixForm[(FA = {{V1, c1}, {c1, V2}})]

The eigenvalues of FA are the coefficients for the linear combinations of the vectors vA11  and vA21 that give the true normal A modes,

{coeff1, coeff2} = Eigenvectors[FA]

Thus, the pure normal A modes are,

purevA11 = coeff1 . {vA11, vA21}

purevA21 = coeff2 . {vA11, vA21}

Giving values to the coefficients V1, V2 and c1, we could compute the numerical values for these vectors and perform the animations. We will not do this here.
  Now we can start with the complex representations \(\(E\_+\)\) and \(\(E\_-\)\). Since they are conjugate, we only need to obtain the pure complex displacement vectors for \(\(E\_+\)\) and then we can apply the procedure of the previous section to obtain the pure real displacement vectors.
  The vectors \(v\_11\%\(E\_+\)\) and \(v\_21\%\(E\_+\)\) are,

vEplus11 = basis[[2, 1, 1]]

vEplus21 = basis[[2, 2, 1]]

The most general form for the \(f\^\(E\_+\)\)matrix is,

MatrixForm[(FEplus = {{W1, b1 + I d1}, {b1 - I d1, 
              W2}})]

The eigenvalues of FEplus are the coefficients for the linear combinations of the vectors vEplus11  and vEplus21 that give the pure  normal \(\(E\_+\)\) modes,

{coeff1, coeff2} = Eigenvectors[FEplus];

Finally, the pure normal \(\(E\_+\)\) modes are,

purevEplus11 = coeff1 . {vEplus11, vEplus21};

purevEplus21 = coeff2 . {vEplus11, vEplus21};

These are complex vectors. Now we can perform the unitary transformation of Equations 19 and 20 to obtain the pure real displacement vectors that correspond to the two copies of the physically irreducible representation  E = \(TraditionalForm\`\(E\_+\)\) ⊕ \(TraditionalForm\`\(E\_-\)\) that appear in (Γ^TVR).
  Finally, we will consider the same molecule but with the full symmetry group, that is (Subscript[D, 3  h]). In this case there is a two dimensional IR,

<< "PointGroups/D3h.dat";
Off[General::spell1];

points = {{0, 1, 0}, {0, 0, 0}};
masses = {1, 1};
mol = GenerateMolecule[group, points]

masses = GenerateMasses[group, mol, points, masses];
modes = NormalVibrations[group, pos, ch, vibrchar, mol];
modes . symbols

vibrep = VibrationalRepresentation[group, mol];
basis = 
        GenerateVibrationalBasis[irreps, modes, vibrep, mol, 
          masses];

vE11 = basis[[5, 1, 1]];
vE21 = basis[[5, 2, 1]];

The most general form for the (f^(E')) matrix is the same as the previous FA,

FE = FA;
{coeff1, coeff2} = Eigenvectors[FE];

Thus, the pure vectors (Subscript[v, 11]^(E')) and (Subscript[v, 21]^(E')) are,

purevE11 = coeff1 . {vE11, vE21};
purevE21 = coeff2 . {vE11, vE21};

Now we have to compute the pure vectors corresponding to  (Subscript[v, 12]^(E')) and (Subscript[v, 22]^(E')), that is the second row partners for purevE11 and purevE21. But as we have seen in Equation 22 , the coefficients for the linear combinations are coeff1 and coeff2,

vE12 = basis[[5, 1, 2]];
vE22 = basis[[5, 2, 2]];
purevE12 = coeff1 . {vE12, vE22};
purevE22 = coeff2 . {vE12, vE22};

8. Conclusion.

We have implemented a Mathematica package to solve classical molecular vibrations problems using the projector method. In an introductory course of Group Theory for Physics students, this method is only marginally mentioned. This is due to the fact that to solve small molecules, there are usually some tricks that one can exploit to avoid the use of projectors. For example, in reference [6], in the computation of the displacement vectors for the (Subscript[NH, 3]) molecule, three methods are explained, each with more algebra than the previous one, that can be used only in certain limited problems. The projector method is usually deferred as the last recourse if the tricks do not work. Although it is certainly true that the work needed to perform a calculation with the projector method by hand is tedious, it is the only with general applicability and it is conceptually simple. The use of Mathematica allows one to concentrate on the general aspects of group theory that have general applicability to this kind of problem. Also these principles appear in other physical problems. For example, the projector method is applicable to the molecular orbitals (MO)  method using the linear combination of atomic orbitals (LCAO) method for solving electronic structure problems in molecules. Basically, one has only to replace the first argument in the function VibrationalRepresentation  by the atomic orbital representation under study, that is the restriction to the symmetry group of the molecule of one of the irreducible representations (D^((j))) of the rotation group. In fact, the package Molecules.m is a subset of a more general package we developed that includes tools to treat MO problems by the LCAO method.

REFERENCES

[1] Cornwell, J. F., Group Theory in Physics, Academic Press, 1984.

[2] Sosnovsky, A. and Demarco, G., Math in Ed. and Res. 6 (1997) 5.

[3] Sosnovsky, A.  and Demarco, G., submmited to Math. in Ed. and Res.  (1997) .

[4] Landau, L. D. and Lifshitz E.M., Quantum Mechanics, Pergamon, 1989

[5] Hamermesh M., Group Theory and Its Application to Physical Problems, Addison-Wesley, 1962.

[6] Lax M., Symmetry Principles in Solid State and Molecular Physics, Wiley, 1974.

ABOUT THE AUTHOR

Gustavo Demarco obtained a Ph.D. in Physics from Instituto Balseiro in 1992. Since that time he has been at Centro Atomico Bariloche working in the Isotope Separation Group and now in the Advanced Designs and Economic Evaluation Division. He is also an Assistant Researcher of CONICET and an Assistant at Instituto Balseiro. His main  interests are research in pulse propagation in resonant and near resonant media, and the utiliziation of Mathematica as a pedagogical tool in Physics courses.

Gustavo Demarco
Centro Atomico Bariloche
8400-Bariloche
R.N.
Argentina
demarco@cab.cnea.edu.ar

ELECTRONIC SUBSCRIPTIONS

Included in the distribution for each electronic subscription is the file molvib.nb, containing Mathematica code for the material described in this article.