Lab 6: Vibrational Frequencies and Normal Modes

Back to Course Overview

As we saw in the previous lab, we get easy access to the first derivatives of the total energy with respect to various parameters in a basic self-consistent calculation of the density of a system via the Hellmann-Feynman theorem. However, if we want to access second derivatives, so that we can calculate, for example, vibrational frequencies, we need to do a little more.

Vibrational Modes of Methane

We’ll start by calculating the vibrational modes of a molecule. This calculation is set up in 2 parts:

  1. A self consistent calculation of the density and wavefunction. You might notice this is often the first step in many of the calculations we do.
  2. A DFPT calculation of the dynamical matrix and normal modes. As we have a molecule, we only need look at the gamma point here.

This calculation is set up in the 01_CH4 directory.

  • First take a look at the input file for the scf calculation. This doesn’t have any inputs you haven’t seen before. There are three things you should take note of though:
  1. We’ve specified the gamma point only by asking for a 1x1x1 unshifted grid, whereas previously we’ve usually explicitly asked for the gamma point. This latter option uses optimizations that can make the calculation a bit faster if we only use the gamma point, however it does not generate output that is compatible with the subsequent DFPT calculation of the frequencies, so we instead specify the gamma point in the manner done in the file.

  2. We’ve re-oriented the molecule relative to what you’ve seen previously, and the atomic positions have also been optimized. By specifying the atoms in this manner, where the components of the hydrogens along each direction have the same magnitude but different signs, the code is much better able to detect the symmetry of the molecule. This is very important for calculations of the vibrations. If the code understands that all the hydrogen atoms are equivalent by symmetry, it need only calculate the derivatives of the energy with respect to one of the hydrogen positions, and then it can populate the full dynamical matrix based on the symmetry of the system.

  3. We’ve increased the default value of ecutrho which is usually four times the value of ecutwfc. This can help alleviate some issues with acoustic mode frequencies at gamma mentioned below.

  • Use this input file to run a calculation with pw.x, and save the output to a file. Check it worked as expected.

  • Now take a look at the input file for the DFPT calculation. This is as follows:

phonons of CH4 (gamma only)
  tr2_ph = 1.0d-15
  asr = .true.
0.0 0.0 0.0
  • This input file is for the ph.x code which uses the output of a SCF calculation with pw.x to do a DFPT calculation of the dynamical matrix. Many more input options are available than we’ve used here. These are described in the INPUT_PH.txt file of the quantum espresso documentation.

    • Here the first line is a comment. This can be whatever description you like of your input file
    • This is followed by a single section: INPUTPH. We’ve specified two variables here.
      • tr2_ph which defines the threshold for self-consistency, with a default of 1.0d-12.
      • asr which turns on the use of the acoustic sum rule which sets the three translational modes at gamma to zero (or at least near zero).
    • Following this section, and depending on what input variables have been defined, we can give a list of wavevectors at which we want to calculate the vibrational frequencies. As with the electronic problem, we are modelling a molecule, so only the gamma point makes sense.
  • Use ph.x to run this calculation and save the output to a file.

  • You’ll see along with the output file, several files and folders have been generated:

    • matdyn (the name of this file can be changed in the input file) - this stores the dynamical matrix, which can then be passed to several other codes for further analysis. Take a look at it now.
      • The top of the file has the atoms, their masses (in internal units) and positions.
      • This is followed by the dynamical matrix (3x3 complex numbers for each possible pair of atoms - all the imaginary parts are zero at gamma).
      • Finally the frequencies of each calculated vibrational mode, both in THz and cm-1 are given along with the corresponding wavevector. A molecule with 5 atoms will have 15 possible vibrational modes.
    • _ph0 - contains intermediate and restart data for the ph.x calculation.
  • Now let’s look through the output of the ph.x calculation.

    • Near the top you’ll see an analysis of the number of calculations that need to be done given the symmetry of the system. You’ll note that only 3 of the 15 total representations need to be explicitly calculated. This is why it’s so important to have the code detect the symmetry of your system correctly. Potentially the calculation could have been 5 times slower.
    • This is followed by a self consistent calculation for each of the representations that need to be calculated.
    • Then we have the frequencies obtained by diagonalizing the dynamical matrix.
    • Finally we have a symmetry analysis of the resulting modes, giving the mode symmetry and whether they should be detectable by Raman or infrared spectroscropy or both.
  • You’ll likely have noticed we have 3 quite negative frequencies. These are the rotational modes of the molecule. It’s difficult to get these to come out to be zero in the DPFT calculation. In principle this could be enforced by a suitable constraint in the code, but this is not available in ph.x. It’s possible in some other codes however.

  • The subsequent three modes near zero are the three translational modes, and these have been enforced to be close to zero by the asr setting in the input file.

  • Note as the frequency is calculated using the square root of the curvature of the energy, a negative frequency is in fact a convention used to indicate that the frequency is imaginary, and the energy curvature is negative. This can indicate some instability in the system, so it is good to ensure you are starting from an optimized structure as you learned how to do last week.


  • Try running the calculation with a higher energy cut-off to see how well converged the vibrational frequencies are.

Phonon Bandstructure of Carbon Diamond

Now that you’ve seen how to calculate the frequencies at a single wavevector, the next step is to calculate the phonon band dispersion of a crystal. The directory 02_CarbonDiamond contains a set of input files for carbon diamond.

This calculation has five steps:

  1. Perform self-consistent calculation of the density and wavefunction.
  2. Calculate the dynamical matrix on a set of wavevectors. I’ll call these q-points from here on, and use k-points to refer to electronic wavevectors.
  3. Fourier transform our grid of dynamical matrices to obtain a set of real space force constants.
  4. Perform an inverse Fourier transform of the real space force constants to obtain the dynamical matrix at an arbitrary q-point. This allows us calculate mode frequencies for a fairly dense set of points along lines between high-symmetry points in the same manner as for the electronic band structure.
  5. Generate the plot.
  • First examine the scf input file in This is similar to those we’ve seen before. This time we’ve specified prefix. This can be useful if you’ve got different calculations in the same directory. Again, we’ve increased ecutrho from its default. The rest is as before.
  • Run pw.x using this input file and save the output.

    • The other output files generated will be in since we set the prefix.
    • Check the output to make sure it completed all right.
  • Now examine the ph.x input file in

    • This is structured as before, but now:
      • We have specified the same prefix as in the scf input file
      • We’ve set ldisp = .true. which says we’re going to calculate a grid of q-points.
      • nq1, nq2 and nq3 define our q-point grid then.
  • Run ph.x using this input file and save the output. This will take maybe 5 or 6 minutes to complete, so feel free to read ahead while it’s running.

  • Once complete you’ll see several files have been generated in the calculation directory.

    • matdyn0 has the q-point grid used followed by a list of the q-points at which the dynamical matrix has been calculated.
    • matdyn1, matdyn2, etc is the dynamical matrix at each calculated q-point as before.
  • Now take a look through the main output file.

    • You’ll see it’s quite similar to the methane case, but now contains a section for each calculated q-point where it figured how how many represenations need to be calculated based on the symmetry, an scf cycle for each, and the frequencies and mode symmetries obtained by diagonalizing the dynamical matrix at that q-point.
    • It’s good to look through this file and ensure you’re getting reasonable numbers for the frequencies before proceeding.
  • Next we’ll be generating the real space force constants from our grid of dynamical matrices using the q2r.x code.

  • q2r.x doesn’t have a help file like the other codes. You need to inspect the source file in PHonon/PH/q2r.f90 if you download a copy of the quantum espresso source code. The inputs are all described in a comment at the top of this file. On the mt-student server you can see this file at /opt/build/quantum-espresso/q-e-qe-6.3/PHonon/PH/q2r.f90.

  • Take a look at the q2r.x input file ` The contents are as follows:

    fildyn = 'matdyn'
    zasr = 'simple'
    flfrc = 'CD444.fc'
  • Here we’ve specified:

    • fildyn, which is the name of the dynamical matrix files from ph.x, where we left it at the default value of matdyn.
    • zasr, which is the approach used to enforce the acoustic sum rule and make the acoustic modes go to zero at the gamma point.
    • filefrc is the name of the file in which to output the force constants in real space.
  • Run q2r.x with this input file now and save the output. This will run almost instantly.

  • The output isn’t particularly interesting. The CD444.fc file has the force constants for each pair of atoms in each direction for each of the 4x4x4 unit cells in real space.

  • Finally we want to use this to generate our normal mode dispersion. We’ll be doing this with the matdyn.x code. As with q2r.x this doesn’t have a doc file describing its input variables. But you can check the comments at the top of its source file to get their details. On the mt-student server this is at /opt/build/quantum-espresso/q-e-qe-6.3/PHonon/PH/matdyn.f90.

  • Take a look at our input file The content is as follows:

    asr = 'simple'
    flfrc = 'CD444.fc'
    flfrq = 'CD-bands.freq'
    0.000 0.000 0.000 30
    0.375 0.375 0.750 10
    0.500 0.500 1.000 30
    1.000 1.000 1.000 30
    0.500 0.500 0.500 30
    0.000 0.500 0.500 30
    0.250 0.500 0.750 30
    0.500 0.500 0.500 0
  • Here we’re setting:

    • asr to again tell it to try to force the acoustic modes at gamma to go to zero.
    • flfrc to give it the name of the file with the real space force constants from the q2r.x calculation.
    • flfrq to give it the name of the output file to store its calculated frequencies.
    • dos to tell it we’re not calculating a density of states
    • q_in_band_form to tell it we want to calculate bands between high symmetry points.
    • We finally give it a number and list of high symmetry points with the number of points to calculate along each line, in the same way as we did for the electronic band structure.
  • Run matdyn.x using this input file now and save the output. Again this is very fast.

  • There’s very little actual output from the code itself, but it will generate the files CD-bands.freq and Both of which contain the frequencies along the lines we requested. Now we need to turn it into a plot.

  • Finally we want to generate a plot of these frequencies. We could do that directly with the previous output: is a multicolumn space separated list of the frequencies in cm-1.

    • It can be easier to use the plotband.x tool to help generate a plot.
    • Call this with plotband.x CD-bands.freq. This will then ask you for an Emin and Emax value - you should pick values equal to or below and above the numbers it suggests. Then it will ask you for an output file name. Pick “CD-bands-gpl.dat” here. You can then cancel further running of this code with ctrl+c. Note it has output the location of the high-symmetry points along its x-axis.
  • Once you’ve done this, a gnuplot script plotbands.gplt has been provided that you can use to generate the band structure plot. This is very similar to the one used for the electronic band structure, but the location of the high-symmetry points along the x-axis has changed as has the name of the data file we’re plotting. Run this with gnuplot plotbands.gplt.

Phonon Density of States of Carbon Diamond

As you might have noticed, we’ve already done most of the work needed to also calculate a phonon density of states. This can be done following the q2r.x calculation by choosing some different input options for matdyn.x. Take a look at the input file The contents are as follows:


This is quite similar to the band plot, but now we’re setting dos to true, choosing a dense q-point grid on which to recalculate our frequencies, and specifying an output for the density of states.

Run matdyn.x now with this input file. It’ll take a bit longer than the band calculation as it uses the tetrahedron method to integrate the density of states. The file it generates: CD-q20.dos can be plotted in whatever plotting software you like.


In this lab we have seen

  • For a molecule:
    • How to use the pw.x code to calculate a converged density and wavefunction and then use these with the ph.x code to calculate the vibrational modes using DFPT.
  • For a crystal:
    • How to use the pw.x code to calculate a converged density and wavefunction.
    • How to use these with the ph.x code to calculate the phonon modes on a grid of wavevectors.
    • How to transform these to obtain real-space force constants using the q2r.x code.
    • How to use the real-space force constants to obtain a phonon mode bandstructure along lines between high-symmetry points in the Brillouin zone.
    • How to plot the bandstructure with gnuplot in a similar way to how we plot the electronic bandstructure in a previous lab.

Extra: Task

  • Calculate the phonon band structure and density of states of silicon.