Modes of vibration of the benzene molecule

Have you set up the local environment?

If not, do that now before proceeding.

The goal of this exercise is to compute the vibrational frequency of a molecule (benzene).

Hint

Enter the Benzene directory.

Before running the calculation to compute the vibrational frequencies, the first step is to relax the geometrical structure of the system under study.

To start with, we run a conjugate gradient minimization to relax the atomic positions. The input file has been prepared for you in the file benzene.relax.fdf. See how the structure of the benzene molecule has been introduced in the Z-matrix format (especification of internal variables, such as distances, angles, and torsional angles). This allows the minimization including some constraints in the symmetry in a trivial way. We just allow the C-C and C-H distances to relax:

%block Zmatrix
molecule
2 0 0 0  xm1 ym1  zm1   0 0 0
2 1 0 0  CC  90.0 60.0  0 0 0
2 2 1 0  CC  CCC  90.0  0 0 0
2 3 2 1  CC  CCC  0.0   0 0 0
2 4 3 2  CC  CCC  0.0   0 0 0
2 5 4 3  CC  CCC  0.0   0 0 0
1 1 2 3  CH  CCH  180.0 0 0 0
1 2 1 7  CH  CCH  0.0   0 0 0
1 3 2 8  CH  CCH  0.0   0 0 0
1 4 3 9  CH  CCH  0.0   0 0 0
1 5 4 10 CH  CCH  0.0   0 0 0
1 6 5 11 CH  CCH  0.0   0 0 0
constants
  ym1 5.00
  zm1 0.00
  CCC 120.0
  CCH 120.0
variables
  CC 1.390
  CH 1.090
constraints
  xm1 CC -1.0 3.903229
%endblock Zmatrix

We run Siesta to carry out the relaxation:

siesta  benzene.relax.fdf > benzene.relax.out

Computing the force constants in real space

We have prepared an input file, called benzene.ifc.fdf to run Siesta and compute the interatomic force constants in real space. In this, we have already copied the relaxed coordinates and unit cell from the benzene.XV generated previously into the AtomicCoordinatesAndAtomicSpecies block. We have also included the atomic masses after the coordinates of each atom. This will be useful for the next step with the vibra code

Here are the relevant sections:

LatticeConstant     1.0 Bohr
%block LatticeVectors
 21.938124322       0.000000000       0.000000000
  0.000000000      20.556799916       0.000000000
  0.000000000       0.000000000      11.910755412
%endblock LatticeVectors

AtomicCoordinatesFormat NotScaledCartesianBohr
%block AtomicCoordinatesAndAtomicSpecies
   4.732644349   9.448630623   0.000000000  2   12.01070
   6.054339080  11.737873050   0.000000000  2   12.01070
   8.697728543  11.737873050   0.000000000  2   12.01070
  10.019423274   9.448630623   0.000000000  2   12.01070
   8.697728543   7.159388196   0.000000000  2   12.01070
   6.054339080   7.159388196   0.000000000  2   12.01070
   2.637688094   9.448630623   0.000000000  1    1.00794
   5.006860953  13.552158387   0.000000000  1    1.00794
   9.745206671  13.552158387   0.000000000  1    1.00794
  12.114379529   9.448630623   0.000000000  1    1.00794
   9.745206671   5.345102860   0.000000000  1    1.00794
   5.006860953   5.345102806   0.000000000  1    1.00794
%endblock AtomicCoordinatesAndAtomicSpecies

Note

If you want to copy your own coordinates from the XV file, you’ll need to do some rearrangements. Remember that for vibra, the AtomicCoordinatesAndAtomicSpecies block needs to look like this:

X    Y    Z    SpeciesIndex   Mass

To compute the interatomic force constant in real space, we have to run Siesta:

siesta < benzene.ifc.fdf > benzene.ifc.out

The interatomic force constant matrix in real space is stored in a file with name of the form SystemLabel.FC.

Computing the dynamical matrix at the Gamma point, and (phonon) modes

Once the interatomic force constants in real space have been computed, a discrete Fourier transform is performed to compute the dynamical matrix in reciprocal space. Then, the dynamical matrix is diagonalized and its eigenvalues and eigenvectors are computed. This is done using the VIBRA code.

In the case of a molecule, only the Gamma point is relevant. It is specified in the same way as to compute the electronic band structure, in the same file benzene.ifc.fdf:

Eigenvectors    .true.        # Compute both phonon eigenvalues and eigenvectors
BandLinesScale  pi/a
%block BandLines
1   0.0   0.0   0.0   \Gamma  # Only the Gamma point (enough for a molecule)
%endblock BandLines

To compute the vibrational frequencies:

vibra < benzene.ifc.fdf > vibra.out

The output of this code is:

  • SystemLabel.bands: with the different mode frequencies (in cm^-1). They are stored in the same way as the electronic band structure.

  • SystemLabel.vectors: with the eigenmodes at Gamma (the format is self-explained).

Note

Some of the modes might have negative frequencies. How could that be?

How to visualize the normal modes

After getting the .vectors file (calculated by vibra) and the .XV file (computed in Siesta), run the vib2xsf program.

You have to answer a few question on the fly, regarding the name of the files where the .vectors are stored, the units to be used to introduce the lattice vectors (Bohrs or Angstroms), the zero of coordinates, the unit cell lattice vectors, the first mode to visualize, the last mode to visualize, the amplitude of the modes to be visualized, and the number of steps in the movie.

You can play a little bit, but to save time we have prepared all the answers in the file vib2xsf.dat for you. Just run:

vib2xsf < vib2xsf.dat

This will produce two files per mode:

  • .XSF file: contains a static structures (as in .XV)

  • .AXSF file: contains the animation of a phonon, for a (user-chosen) amplitude and number of steps.

They can be visualized using OVITO:

$ ovito benzene.Mode_1.AXSF (Animation XCrySDen Structure File)

Or open ovito

$ ovito
File > Load File
Load benzen.Mode_*.AXSF

Note

It might be interesting to analyze those modes with negative or small frequencies.