Day 9-11: First Principles: Phonons and Spectroscopy Tutorial¶

If you prefer you may do any of the tutorials on alternative types of calculations which you can find on http://www.castep.org

The aims of this tutorial session is to introduce you to running larger-scale CASTEP jobs on supercomputer clusters, creation of input files for and setup of CASTEP jobs, and analysis of the results using standalone tools. In order to achieve results without waiting too long for jobs to complete the initial runs will be small ones, but the aim by the end of this afternoon should be to set up larger and more significant runs. Today’s session will entirely comprise insulators or semiconductors.

You may need to consult the CASTEP phonon users guide, which may be accessed at http://www.tcm.phy.cam.ac.uk/castep/Phonons_Guide/Castep_Phonons.html

These calculations should be run on the dedicated high-performance cluster system SCARF. Please see the separate instructions for logging on to these systems and transferring files to and from the PC/Laptop system.

A. \(\Gamma\)-Point phonon in h-BN¶

The first exercise will take you through the construction of a \(\Gamma\) point phonon calculation and the generation of a simple model infrared spectrum of a semiconductor. It will also introduce you to the use of some additional analysis and visualisation tools.

BN is one of a family of Nitride semiconductors, which occurs in cubic zincblende (c-BN), hexagonal wurtzite (w-BN) and graphite-like hexagonal (w-BN) polymorphs. We will calculate the phonons, infrared spectrum and raman spectrum of h-BN.

Your starting point will be the structure which is provided in a castep .cell file named h-BN.cell which you can copy from

cp /work4/training/ccp5/abinitio/h-BN/h-BN.cell .
cp /work4/training/ccp5/abinitio/h-BN/h-BN.param .

The Phonons manual (http://www.tcm.phy.cam.ac.uk/castep/Phonons_Guide/Castep_Phonons.html) contains a very similar example for the Wurtzite-structure polymorph of BN.

Note in the cell file the %block species_pot which reads

%block species_pot

NCP

%endblock species_pot

to select the norm-conserving pseudopotentials library. Density functional perturbation calculations requre norm-conserving peeudopotentials (the computational savings of ultra-soft potentials is lost with all of the addtional terms these require in perturbation theory).

To specify a gamma-point phonon wavevector

phonon_kpoint_mp_grid 1 1 1

Test your configuration using

castep.mpi –dryrun <seedname>

Once you have in place the <seed>.cell and <seed>.param files you are ready to submit the CASTEP job. When it has finished, you can examine the output file h-BN.castep and find the frequencies. What you see is explained further in the Phonons user guide on the web link provided above. There is also a machine-readable file h-BN.phonon which contains the frequences, but also the eigenvectors which we will analyse.

Analysis of h-BN phonon output¶

We will use the academic tools to visualise the modes using the free Jmol visualiser.

Use secure file transfer to copy the .phonon file back to the cloud account (where Jmol is installed) or your own device (where you can install Jmol). Start jmol and bring up a console window from the right-mouse menu. Then type the Jmol command to load your .castep file

load h-BN.phonon {3 3 2} PACKED

(You may need to supply a full pathname for Jmol to find the .phonon file). Note that while you can read in the .phonon file from Jmol’s File menu, you need the additional opitions of the command line to display additional periodic repeats.

You can then use the “Tools-Vibrate” menu to turn on mode animation, and navigate the modes. Can you see from the mode eigenvectors which modes are ir active and which are raman active? Do you agree with the ir and raman activity printed in the .castep file?

Generation of ir spectrum¶

The easiest way to generate a simple model ir spectrum is to use the “dos.pl” tool. To run this in the most effective way and automatically display a plot, you will need to be running an X windows server. In that case the command

module load xmgrace

dos.pl -ir -xg h-BN.phonon

will generate a plot script and use the “xmgrace” plotting program to display it. You can create a GNUPLOT script instead of xmgrace by changing the “-xg” flag to “-gp”. Alternatively you can generate a GNUPLOT script without plotting by

dos.pl -ir -gp -np h-BN.phonon > h-BN-phonon.plt

Raman spectrum¶

The calculation of a raman intensities is fairly expensive compared to infrared matrix elements and it is therefore not turned on by default. To enable this, add the keyword

CALCULATE_RAMAN : TRUE

RAMAN_METHOD: DFPT

to the h-BN.param file, and resubmit the job. You can then use the “-raman” flag of dos.pl to generate and plot a raman spectrum.

You can then copy the “h-BN-phonon.plt” back to the PC and read this into GNUPLOT.

B. Molecular modes in benzene¶

The next part of this practical is to compute the modes and spectrum of a molecule and compare the result with a calculation of a molecular crystal. Our example is benzene. .cell and .param files describing a benzene molecule and also of a structure describing a high-pressure crystalline polymorph, phase III. These are in both in the course directory /work4/training/ccp5/abinitio and called benzene-..

First run the isolated molecule calculation. There are a few other considerations to take into account for an isolated molecule calculation.

The size of the simulation cell governs the interactions between periodic copies of the molecule and should be large enough that these are negligible.
The shape of the simulation cell governs the crystallographic point groups allowed in the handling of the symmetry. It should be chosen to be commensurate (as far as possible) with the molecular point group to maximise the use of symmetry. In the case of benzene it should obviously be hexagonal, and recommend a box 8A by 8A by 4A.
Recall that there is no electronic dispersion for a molecule, so only a single electronic k-point is needed. The general rule is that all “molecule in a box” calculations should use the \(\Gamma\) point only as CASTEP uses special performance optimisations in this case.

For simplicity use the local density approximation, the standard NCP library pseudopotentials.

In benzene not all atoms are on special symmetry sites so you should first perform a geometry optimisation. Then set up a follow-on calculation to compute the \(\Gamma\) point phonons. The Phonons User Guide on the WWW should help you fill in the details, but please ask if you are stuck.

Benzene Phase III: The next stage is to compute the gamma point phonon modes of the molecular crystal of benzene in the high pressure polymorph, Phase III. You should use a 2x2x2 grid of electronic k-points, as dispersion is non-zero in this molecular crystal. Make sure that symmetry is detected and enabled. Use the same .param file as for the molecular case to ensure the settings are the same.

Once these calculations have completed you should generate a phonon DOS and ir spectra as in the previous practical and compare the molecule with the molecular crystal. You can also use Jmol to identify the modes.

Are all your frequencies positive? If not, can you suggest why not? Think about periodic boundary conditions.