Simulating the IR and Raman Spectra of Benzene using DFPT

We are going to simulate the infra-red (IR) and Raman spectra of a benzene molecule (see above) using a combination of finite differences and DFPT.

Workflow

We are going to carry out:

1. Geometry Optimization: Perform a geometry optimization to obtain the equilibrium structure of the benzene.
2. Hessian Calculation: Compute the Hessian matrix, which contains second derivatives of the energy with respect to atomic displacements. This can be done using finite difference methods, where small displacements are applied to each atom, and the numerical derivative of the forces are used to approximate the Hessian.
3. Diagonalization: Diagonalize the Hessian matrix to obtain the vibrational frequencies and normal modes of the molecule.
4. IR Intensity Calculation: Calculate the IR intensities for each vibrational mode. This involves evaluating the transition dipole moment between the ground state and excited state for each mode.
5. Polarizability Calculation: Compute the polarizability tensor using DFPT. This involves calculating the second derivatives of the total energy with respect to an external electric field.
6. Raman Intensity Calculation: Use the polarizability tensor to compute the Raman intensities for each vibrational mode.
7. Spectrum Generation: Finally, generate the IR/Raman spectrum by plotting the IR/Raman intensities as a function of frequencies.

Let us first create a directory in which we will group the calculations for this part of the tutorial and move into it:

mkdir 1_benzene
cd 1_benzene

Geometry optimisation

As with any FHI-aims calculation, (see the Basics of Running FHI-aims tutorial) we require a geometry.in file, and a control.in file, and we follow the rule of "one directory per calculation". Therefore, for the geometry optimization, let us make a directory and move to it:

mkdir relaxation
cd relaxation

To obtain a geometry.in file, one can use a simple online structure builder, or a molecular visualiser such as Jmol to create an .xyz file, and then edit the formatting to match that of the FHI-aims geometry.in style. For convenience, we provide a geometry.in file for benzene below:

geometry.in

atom       1.22178135     -0.16537864      2.16412218 H
atom       0.68642597     -0.09302867      1.21582627 C
atom      -0.71157626     -0.03547226      1.20415364 C
atom      -1.26664300     -0.06315013      2.14343935 H
atom      -1.39799604      0.05762918     -0.01158666 C
atom      -2.48842294      0.10240052     -0.02047005 H
atom      -0.68642026      0.09302085     -1.21585173 C
atom      -1.22181116      0.16536325     -2.16412612 H
atom       0.71156432      0.03543109     -1.20417122 C
atom       1.26657727      0.06309837     -2.14349204 H
atom       1.39799643     -0.05768697      0.01157649 C
atom       2.48842432     -0.10242659      0.02047991 H

It is important that we relax this structure into a stable ground state, with well converged forces. Particularly because we are employing finite differences to calculate the Hessian, which takes a difference between 2 force values (that are small since the molecule is close to the ground-state).

Thus, for this system we are going to use a couple of control.in keywords to ensure the accuracy of our results.

1. The relaxation criterion (relax_geometry bfgs) is set to 1E-4, ensuring efficient convergence without excessive computational costs.
2. We set the accuracy of the charge density (sc_accuracy_rho) to 1E-6. By setting an accurately converged density, we ensure that the forces will also be computed accurately.

The keywords section of the control.in file will then look like the following:

control.in

xc                 pbe
spin               none
sc_accuracy_rho    1E-6
relax_geometry bfgs 1E-4

We have also set a PBE exchange-correlation funtional with no spin. We will also append lightdense species defaults to this control.in file using clims

clims-prepare-run --species lightdense

lightdense species defaults

lightdense species defaults are a recent augmentation to the "standard" FHI-aims species defaults of "light", "intermediate", and "tight". They utilise the same basis functions as the light species defaults, but with denser integration grids, which can provide more accurate forces.
In principle, one should use tight species defaults to ensure convergence with respect to basis functions, as well as having even denser integration grids for more precise forces. However, in the interest of reducing computational time for this tutorial, lightdense species defaults were found to give qualitatively correct results for this example.

We can now run the FHI-aims calculation. For webinar users, you can copy the provided SLURM run script fhiaims_job.sh from your home directory to your current directory

cp ~/fhiaims_job.sh .

and then submit to the HPC cluster

sbatch fhiaims_job.sh

Example submission script fhiaims_job.sh

fhiaims_job.sh

#!/bin/bash -l
# Standard output and error:
#SBATCH -o ./slurm.stdout
#SBATCH -e ./slurm.stderr
# Initial working directory:
#SBATCH -D ./
# Job Name:
#SBATCH -J FHI-aims
# Number of nodes and MPI tasks per node:
#SBATCH --nodes=1
# HPC7a
#SBATCH --tasks-per-node=192
#
# Wall clock limit:
#SBATCH --time=00:15:00

module purge
module load fhi-aims aims

date
mpirun aims.x > aims.out

Running calculations on other resources

If you are running the calculations on a different machine, make sure to adjust the submission script specific to your system and submission system. You can also alternatively run it locally on your own PC, as this is not too expensive a calculation.

Once the calculation has finished, you should have the following output files in your directory:

Output files

filename	description
aims.out	Standard aims output file
geometry.in.next_step	Relaxed geometry which we shall use to calculate the IR/Raman spectrum
hessian.aims	Estimate of the Hessian by the geometry optimisation algorithm

Generation of finite-difference displaced structures

Next, we need to generate the finite-difference displaced structures, and then for each of them:

• Compute the dipole moment in order to obtain the IR intensities.
• Compute the polarizability tensor with DFPT to obtain the Raman intensities.

Let us create and move to a new directory for this set of calculations

cd ..
mkdir displacements
cd displacements

Inside this directory let us copy over the relaxed benzene geometry and rename it from geometry.in.next_step to geometry.in, as well as the control.in file from the relaxation

cp ../relaxation/geometry.in.next_step geometry.in
cp ../relaxation/control.in .

We need to modify the control.in file by removing the relax_geometry bfgs 1E-4 line, and adding the following keywords:

• output dipole : will calculate the dipole moment that will be used for the IR intensities.
• compute_forces .true. : will compute the forces, which are needed to obtain the second derivative of the energy from the numerical change on the forces arising from the small finite displacements.
• final_forces_cleaned .true. : removes spurious unitary force components (rotation and translation of the whole structure due to residual numerical noise) which is necessary for the finite-difference calculation of vibrational frequencies.
• DFPT polarizability : will calculate the polarizability tensor with DFPT that will be used for the Raman intensities.
• dfpt_sc_accuracy_dm 1E-6 : increases the accuracy threshold for the first order density matrix in the DFPT cycle to yield a more precise polarazability tensor.
• dfpt_mixing 0.5 : increases the mixing parameter of the Pulay mixer for the DFPT steps. This rather aggressive value speeds up the convergence, although one should take care when setting this value. For insulators with a large HOMO-LUMO gap, typically a large value such as this can help speed up convergence. However, for metals with a small or no HOMO-LUMO gap, a more conservative value such as 0.02 should be chosen.

The keywords section of the control.in file should then look like the following:

control.in

xc                 pbe
spin               none
sc_accuracy_rho    1E-6
output dipole 
compute_forces .true.
final_forces_cleaned .true.
DFPT polarizability
dfpt_sc_accuracy_dm 1E-6
dfpt_mixing 0.5

And we should once more have lightdense species defaults appended below. Now, let us create the displaced structures. We will use a python script that is available in the utilities/ folder of your FHI-aims installation, named get_vibrations.py. Webinar users should find this in their home directory, and can copy it to their current working folder:

cp ~/get_vibrations.py .

It is also available in the solutions folder of this tutorial.

get_vibrations.py script

This script can compute basic things, such as structure displacements, and IR and Raman spectra of molecules and periodic structures (which will be covered in Part 3 of the tutorial). However, note that for periodic structures it does not include non-analytic corrections^{1 2} and it cannot do more advanced things such as calculating vibrational cross-sections. Treatment of vibrational spectra is ongoing work in the FHI-aims code, and users that have specific functionality requests are invited to contact
aims-coordinators@ms1p.org.

There is also a perl script numerical_vibrations.pl available inside src/vibrations of the FHI-aims code, which is capable of doing similar things, plus additionally computing free energy contributions. For the purposes of this tutorial, we prefer to work with the Python script. Interested readers can consult the FHI-aims manual for further details.

The vibrations script provides a help string which shows the possible modes of operation and options for this script:

python get_vibrations.py -h

The script can be called as such:

python get_vibrations.py <name> <mode>

where <name> is just chosen to name the subdirectories and output files, whilst <mode> is replaced with 0, 1 or 2 depending on the chosen run mode. The run modes are as follows:

get_vibrations.py run modes

run mode	description
0	Generate finite-difference displaced structures
1	Calculate vibrational modes
2	Calculate vibrational modes and IR/Raman intensities

We shall begin by calling the script to generate displaced structures of the geometry.in in subdirectories of the current directory, as well as placing copies of control.in in these subdirectories.

python get_vibrations.py benz 0

This should generate 72 subdirectories, each with a control.in and a geometry.in corresponding to a displaced structure. There are 72 subdirectories as each of the 12 atoms are displaced in the forward and negative direction, for each of the 3 Cartesian axes (12 x 2 x 3 = 72). A file containing the masses of the elements will also be created in the current directory.

Output files

file or dir name	description
<name>.i_atom_*.i_coord_*.displ_*	Subdirectories are generated with this naming structure, each containing a control.in and a geometry.in.
masses.<name>.dat	Geometry file with masses of elements. Can be ignored.

DFPT calculation submission on displaced structures

We will now run each of the 72 DFPT calculations with FHI-aims. We can do this via a submission script that contains a loop to go through each of the subdirectories just created, and run the respective FHI-aims calculation.

Example submission script fhiaims_job_loop.sh

fhiaims_job_loop.sh

#!/bin/bash -l
#SBATCH -o ./slurm.stdout
#SBATCH -e ./slurm.stderr
#SBATCH -D ./
#SBATCH -J FHI-aims
#SBATCH --nodes=1
#SBATCH --tasks-per-node=192
#SBATCH --time=00:15:00
#SBATCH --exclusive

module purge
module load fhi-aims

# Function to run a single aims calculation
run_aims() {
    local dir_name=$1
    cd $dir_name || exit 1
    # Check whether a calculation has already been completed here
    if grep -q "Have a nice day" aims.out; then
        echo "Calculation in $dir_name already complete"
    else
        echo "Starting calculation in $dir_name"
        date
        mpirun -np 8 aims.x > aims.out &
    fi
    cd ..
}


# Variables for the folder names
atoms=(0 1 2 3 4 5 6 7 8 9 10 11)
coords=(0 1 2)
displacements=(-0.0025 0.0025)

# Loop over all the subdirectories
for coord in "${coords[@]}"; do
    for displacement in "${displacements[@]}"; do
        for atom in "${atoms[@]}"; do
            dir_name="benz.i_atom_${atom}.i_coord_${coord}.displ_${displacement}"
            run_aims "$dir_name"
        done
    done
    # Wait for the batch of calculations to finish
    wait
done

For webinar users, this script is once more available in your home folder and can be copied to the current directory

cp ~/fhiaims_job_loop.sh .

and submitted

sbatch fhiaims_job_loop.sh

Once the calculations are finished, each of the 72 subdirectories should contain an aims.out file with "Have a nice day" printed at the end, if the calculation finished successfully. You can check this with the following command:

grep -ri "Have a nice day" benz* | wc -l

which should print 72 to the screen.

Processing results and simulation of IR and Raman spectra

We can now use the get_vibrations.py script in run mode 2 to process the results. This involves computing the Hessian matrix via finite differences of the forces, diagonalizing to obtain the vibrational frequencies and normal modes, computing the IR and Raman intensities, and hence generating the IR and Raman spectra:

python get_vibrations.py -Mp --xmin 500 --xmax 3400 benz 2

The -M argument tells the script to calculate both IR and Raman intensities. The -p argument tells the script to plot the IR and Raman spectra. The --xmin argument provides a lower x limit of the spectra in cm\(^{-1}\). The --xmax argument provides a upper x limit of the spectra in cm\(^{-1}\).

This prints some information to the screen, including the vibrational modes and intensities. There should be six frequencies close to 0 cm\(^{-1}\), which correspond to three translations and three rotations of the structure. The remaining \(3N-6\) frequencies are the vibrational normal modes. The script also produces a number of output files, important ones are listed below:

Output files

file or dir name	description
<name>.ir	List of 3N modes, their frequencies, zero-point energies and IR intensities.
<name>.ir	List of 3N modes, their frequencies, zero-point energies and Raman intensities.
<name>.xyz	XYZ file of vibrational modes. Can visualize the modes using Jmol (see below).
<name>_IR_spectrum.pdf	Image of IR spectrum
<name>_Raman_spectrum.pdf	Image of Raman spectrum
<name>_IR_spectrum.csv	Text file of simulated spectrum, same units as plot.
<name>_Raman_spectrum.csv	Text file of simulated spectrum, same units as plot.

Visualizing the vibrational modes

The vibrational modes can be visualised with Jmol, which can be downloaded for free.

Opening the .xyz file in Jmol, and selecting Tools --> Vibrate --> Start vibration begins an animation of the vibrational modes. Arrows on the far right of the toolbar allow switching between the different vibrational modes.

The last vibrational mode of our benzene molecule should be a symmetric stretch of all CH bonds:

The CH symmetric stretch mode involves a symmetric stretching of the CH bonds. Since the dipole moment does not change during this symmetric vibration, this mode should be IR inactive. There is a change in polarisability under this vibration, so this mode should be Raman active.

Infrared Spectrum

If we now open <name>_IR_spectrum.pdf, we should see a spectrum with four peaks.

We can make a quick comparison to experiment, by plotting our simulated IR spectrum with the IR absorbance spectrum of benzene (arbitrarily scaled) measured in the gas phase³:

Our simulated peak positions match those observed in experiment, however, there are additional peaks observed in experiment, as well as fine structure missing in the simulation.

Our simulation does not include vibrational overtones, combination bands, nor any anharmonic effects. In particular, overtone and combination bands are responsible for the peaks in the 1700-2000 cm\(^{-1}\) region ⁴.

Raman Spectrum

If we now open <name>_Raman_spectrum.pdf, we should see the following spectrum:

By comparing the visualization of the vibrational modes, you should be able to assign the IR and Raman activity of each mode and assign the peaks in each spectrum.

Solutions

You find the solutions to all the above exercises by clicking on the button below.

Show solutions

References

---

1. R. M. Pick, M. H. Cohen, R. M. Martin, Microscopic Theory of Force Constants in the Adiabatic Approximation, 1, 910-920, (1970) ↩

2. S. Baroni, S. de Gironcoli, A. Dal Corso, P. Giannozzi, Phonons and related crystal properties from density-functional perturbation theory, 73, 515-562, (2001) ↩

3. NIST Standard Reference Database 69: NIST Chemistry WebBook ↩

4. https://www.spectroscopyonline.com/view/benzene-fingers-part-i-overtone-and-combination-bands ↩