Eigenvalues and Eigenvectors of a Linear Chain in 1D

In one of the previous Fortran programming exercises, the aim was to determine the normal mode eigenvalues and eigenvectors of a linear chain of N atoms in one-dimension connected only by nearest neighbour harmonic bonds. The cartesian positions of each atom were specified by and the bondlengths by for . The atoms have mass m and the bonds have force constant k. The potential energy can be written as:

The force constant matrix can thus be written as:

and is zero elsewhere.

The above problem can be solved using the inbuilt commands of Mathematic to construct a number of functions by the following approach:

• Firstly, we can construct the potential energy (vpot) by doing a loop over the atoms.

  vpot[n_] := Sum[(k/2)*(x[i+1]-x[i]-re)^2, {i,n-1}]

• The second derivative matrix (d2vdx2) is generated from the potential energy (vpot).

  d2vdx2[n_] := Table[D[vpot[n], x[i], x[j]], {i,n}, {j,n}]

• We can generate the mass weighted force constant matrix, known as the F matrix (fmatrix), from the second derivative matrix (d2vdx2) by scaling by the mass m.

  fmatrix[n_] := d2vdx2[n]/mass

• The (real symmetric) F matrix can then be diagonalised to yield the eigenvalues and/or the eigenvectors. It is trivial to check the latter for orthonormality.

  eigenvalue[n_] := Eigenvalues[N[fmatrix[n]]]
  
  eigenvector[n_] := Eigenvectors[fmatrix[n]]

• The eigenvalues can be easily converted into the corresponding angular frequencies and the vibrational frequencies .

  omega[n_] := Sqrt[eigenvalue[n]]
  
  frequency[n_] := N[omega[n]/(2Pi)]

• Calculations can be performed by setting k and m to particular values and calling the functions with a particular chain length e.g. a 5 atom chain:

  In[10]:=
  k= 39.129396
  
  Out[10]=
  39.129396
  
  In[11]:=
  mass = 1
  
  Out[11]=
  1
  
  In[12]:=
  frequency[5]
  
  Out[12]=
                                                       -10
  {1.89369, 1.61087, 1.17036, 0.615296, 0. + 3.39013 10    I}
  
  In[13]:=
  eigenvalue[5]
  
  Out[13]=
                                                  -18
  {141.571, 102.442, 54.0755, 14.9461, -4.53724 10   }
  
  In[14]:=
  eigenvector[5]
  
  Out[14]=
  {{0.19544, -0.511667, 0.632456, -0.511667, 0.19544}, 
   
                                     -19
    {0.371748, -0.601501, -1.92959 10   , 0.601501, 
   
     -0.371748}, {-0.511667, 0.19544, 0.632456, 0.19544, 
   
                                                   -20
     -0.511667}, {-0.601501, -0.371748, -4.55515 10   , 
   
     0.371748, 0.601501}, {0.447214, 0.447214, 0.447214, 
   
     0.447214, 0.447214}}
  
  In[15]:=
  fmatrix[5]
  
  Out[15]=
  {{39.1294, -39.1294, 0, 0, 0}, {-39.1294, 78.2588, -39.1294, 0, 0}, 
   
    {0, -39.1294, 78.2588, -39.1294, 0}, {0, 0, -39.1294, 78.2588, -39.1294}, 
   
    {0, 0, 0, -39.1294, 39.1294}}

In the following section I give examples of applications for which we might have instead used a Fortran library routine.