Eigenvalues and Eigenvectors of a Linear Chain in 1D

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:

  (ii) a 5 atom chain:

  In[3]:=
  k= 1;
  
  In[4]:=
  mass = 1;
  
  In[5]:=
  vpot[3]
  
  Out[5]=
                     2                      2
  (-re - x[1] + x[2])    (-re - x[2] + x[3])
  -------------------- + --------------------
           2                      2
  
  In[6]:=
  MatrixForm[fmatrix[3]]
  
  Out[6]=
  1    -1   0
  
  -1   2    -1
  
  0    -1   1
  
  In[7]:=
  MatrixForm[eigenvalue[3]]
  
  Out[7]=
  3.
  
  
  1.
  
           -19
  1.0842 10
  
  In[8]:=
  MatrixForm[frequency[3]]
  
  Out[8]=
  0.275664
  
  
  0.159155
  
            -11
  5.24053 10
  
  In[9]:=
  MatrixForm[eigenvector[3]]
  
  Out[9]=
  -0.408248        0.816497         -0.408248
  
                              -20
  0.707107         -5.10372 10      -0.707107
  
  
  -0.57735         -0.57735         -0.57735

  (ii) 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, -0.511667},
   
                                     -20
   {-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}}