Simulating the Raman Spectra of a Benzene crystal using DFPT
This tutorial shares many similarities with the benzene molecule tutorial, albeit we are now investigating a benzene crystal system (pictured above). Our calculations thus have to be adjusted to deal with a system in periodic boundary conditions.
Workflow
We are going to carry out:
- Geometry Optimization: Perform a geometry optimization to obtain the equilibrium structure of the benzene crystal.
- 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.
- Diagonalization: Diagonalize the Hessian matrix to obtain the vibrational frequencies and normal modes of the crystal.
- 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.
- 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.
- Raman Intensity Calculation: Use the polarizability tensor to compute the Raman intensities for each vibrational mode.
- 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 will group the calculations for this part of the tutorial and move into it (N.B. if you are currently in the directory of part 1, you will need to first cd .. out of it):
mkdir 2_benzene_crystal
cd 2_benzene_crystal
Key differences when considering periodic systems
- The keyword,
output dipole, does not work for periodic systems, thus we must move to a Berry phase formalism. We can do this by employing the keywordoutput polarization. - For periodic systems we must use
DFPT dielectricto obtain the polarizability, rather thanDFPT polarizability. - The keyword
output polarizationcannot currently be used withDFPT dielectric, thus IR and Raman spectra require separate calculations. - The keyword
output polarizationcannot currently be used with ScaLAPACK parallelism, so we must use LAPACK parallelism and specifyKS_method serialincontrol.in. This will limit the scalability of our calculations.
Geometry optimisation
As with the molecular calculation, first we need to obtain an optimized geometry upon which we can compute the IR and Raman spectra. Let us create a directory in which to perform the relaxation:
mkdir relaxation
cd relaxation
geometry.in file using the provided geometry below:
lattice_vector 7.36014186 -0.02880954 0.03884039
lattice_vector 0.03886960 9.36606084 -0.02633631
lattice_vector -0.03503709 0.01893797 6.73427189
atom 6.52419029 2.34236335 6.69077530 H
atom 4.48475498 7.03340714 6.66648672 H
atom 6.54243911 2.31747328 3.32326563 H
atom 4.50288053 7.00839124 3.29985039 H
atom 0.83979120 7.01381239 0.05599872 H
atom 2.87923853 2.32276933 0.08027794 H
atom 0.82151326 7.03872010 3.42347877 H
atom 2.86110502 2.34778327 3.44692256 H
atom 5.53910512 0.70810271 1.56500429 H
atom 5.50972230 5.39151547 1.55188061 H
atom 5.53592525 3.94936044 4.92360220 H
atom 5.50636507 8.63292218 4.90953211 H
atom 1.82487386 8.64807263 5.18176485 H
atom 1.85426069 3.96466245 5.19489997 H
atom 1.82806526 5.40681643 1.82317377 H
atom 1.85762818 0.72325578 1.83723586 H
atom 0.96680077 1.61385864 5.11748017 H
atom 2.69274094 6.28993960 5.11278062 H
atom 0.99033149 3.08112362 1.74524429 H
atom 2.71630921 7.75692850 1.74232411 H
atom 6.39717547 7.74231633 1.62928776 H
atom 4.67123882 3.06624155 1.63400081 H
atom 6.37365326 6.27505683 5.00152781 H
atom 4.64768503 1.59924522 5.00444652 H
atom 6.87379544 1.31090478 6.72804724 C
atom 4.12664179 6.00474223 6.70014708 C
atom 6.90026041 3.34630231 3.35493715 C
atom 4.15302252 8.04006835 3.32750960 C
atom 0.49018479 8.04527038 0.01872335 C
atom 3.23734213 3.35143807 0.04663197 C
atom 0.46372772 6.00987768 3.39183670 C
atom 3.21096729 1.31610766 3.41926188 C
atom 6.34013793 0.39203726 0.89728421 C
atom 4.71347042 5.08162828 0.87560063 C
atom 6.33915086 4.25549953 4.25389522 C
atom 4.71238003 8.94517341 4.23167573 C
atom 1.02384154 8.96413750 5.84948585 C
atom 2.65051125 4.27455339 5.87117963 C
atom 1.02484225 5.10067922 2.49288475 C
atom 2.65161088 0.41100261 2.51509406 C
atom 0.53048056 0.91723433 5.83093697 C
atom 3.11583889 5.59007719 5.83101687 C
atom 0.55987497 3.78495388 2.45518196 C
atom 3.14524972 8.45764016 2.45625296 C
atom 6.83349853 8.43894063 0.91583277 C
atom 4.24814175 3.76610384 0.91576398 C
atom 6.80411653 5.57122406 4.29159676 C
atom 4.21874052 0.89853616 4.29051750 C
and then create the control.in file with the following keywords:
xc pbe
spin none
sc_accuracy_rho 1E-6
relax_geometry bfgs 1E-4
charge_mix_param 0.4
k_grid 6 6 6
This is again similar to the molecular case, with the addition of the following keywords:
k_grid: specifies a \(6\times6\times6\) k-grid upon which to perform the calculation. A user should ensure the convergence of this parameter, here we have used the rule of thumb that \(n_i \times a_i > 40\)Å, for the lattice vector length \(a_i\) and the number of k-points \(n_i\) along the k-space direction \(i\).charge_mix_param: specifies a mixing parameter of 0.4 for the Pular mixer during the SCF cycle. Similar todfpt_mixingparameter introduced in the molecular secton, this reduces the number of SCF steps required.
We will then append light species defaults to this file
clims-prepare-run --species light
Then, you can run the calculation by e.g. submitting to HPC resources
sbatch fhiaims_job.sh
Generation of IR spectrum
Next we will generated displaced structures once more. As we need to perform separate calculations for the IR and Rama spectra for the crystal, let us create a directory for the IR first:
cd ..
mkdir IR_spectrum
cd IR_spectrum
control.in file from the relaxation folder
cp ../relaxation/geometry.in.next_step geometry.in
cp ../relaxation/control.in .
control.in file by removing the relax_geometry bfgs 1E-4 keyword, and adding compute_forces .true. and final_forces_cleaned .true. like in the molecular case, plus the following keywords:
output polarization: this calculates the electronic and ionic contributions to the polarization in periodic systems via the Berry-phase formalism along a specific reciprocal lattice vector using an \(n_1 \times n_2 \times n_3\) k-point grid. To compute the 3D polarization we therefore need three of these keywords.KS_method serial: specifies a calculation with LAPACK parallelisation should be used, asoutput polarizationis currently not compatible with ScaLAPACK.
The keywords section of the control.in file is then:
xc pbe
spin none
sc_accuracy_rho 1E-6
charge_mix_param 0.4
k_grid 6 6 6
output polarization 1 20 5 5
output polarization 2 5 20 5
output polarization 3 5 5 20
KS_method serial
compute_forces .true.
final_forces_cleaned .true.
with the three output polarization keywords computing the polarization in each lattice vector direction, with 20 k-points in the direction of the lattice vector, and 5 in the other directions. This k-grid is computed via Fourier interpolation of the SCF k-grid. Light species defaults should once more be appended below the keywords.
We will then generate the displaced structures
python get_vibrations.py crystalIR 0
fhiaims_job_loop.sh script from the first part, accounting for the change in the number of atoms and the naming of the folders, as well as the increased resources that the calculations will require due to the larger size of the system:
Example modified submission script fhiaims_job_loop.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
#SBATCH --tasks-per-node=192
#
# Wall clock limit:
#SBATCH --time=06:00:00
module purge
module load fhi-aims
# Create arrays of all possible values
atoms=({0..47})
coords=(0 1 2)
displacements=(-0.0025 0.0025)
# Loop through all combinations
for atom in "${atoms[@]}"; do
for coord in "${coords[@]}"; do
for displacement in "${displacements[@]}"; do
dir_name="crystalIR.i_atom_${atom}.i_coord_${coord}.displ_${displacement}"
echo "Processing $dir_name"
cd $dir_name || exit 1
# Run aims
date
mpirun aims.x > aims.out
date
cd ..
done
done
done
Once all calculations have finished we can process and produce an IR spectrum.
python get_vibrations.py -Ip --xmin 500 --xmax 3400 crystalIR 2
-I argument tells the script to calculate only IR intensities, and the other arguments are the same as the molecular case.
This produces the same types of output files as the molecular calculation. If we open the file crystalIR.ir we will notice that there are three frequencies at approximately 0 cm\(^{-1}\). These are known as acoustic modes, which correspond to translations of the lattice.
If we open the IR spectrum, we should see the following:
We can again visualize the modes using Jmol (see benzene molecule tutorial), for example, here is the symmetric CH stretch for the crystalline system:
Generation of Raman spectrum
In a separate directory, we will now regenerate finite-difference displaced structures, with DFPT enabled, to simulate the Raman spectrum.
cd ..
mkdir Raman_spectrum
cd Raman_spectrum
control.in and geometry.in file from the IR folder
cp ../IR_spectrum/control.in .
cp ../IR_spectrum/geometry.in .
output polarization and KS_method serial keywords, and add the dfpt_mixing and dfpt_sc_accuracy_dm keywords like in the molecular example. as well as:
DFPT dielectric: calculates the dielectric constant for extended systems using DFPT.
The keywords section of the control.in file will then look like the following:
xc pbe
spin none
sc_accuracy_rho 1E-6
charge_mix_param 0.4
k_grid 6 6 6
dfpt_mixing 0.4
DFPT dielectric
dfpt_sc_accuracy_dm 1E-5
compute_forces .true.
final_forces_cleaned .true.
The dfpt_sc_accuracy_dm parameter is reduced slightly with respect to the molecular case to speed up the calculation.
Once more we will generate the finite difference displaced structures:
python get_vibrations.py crystalRaman 0
fhiaims_job_loop.sh script, remember to modify the crystalIR name to crystalRaman, as well as approximately doubling the time requested to compute.
Then, once more, to process the results we use the get_vibrations.py script:
python get_vibrations.py -Rp --xmin 500 --xmax 3400 crystalRaman 2
-R argument tells the script to calculate only Raman intensities.
If we open the Raman spectrum, we should see the following:
Solutions
You find the solutions to all the above exercises by clicking on the button below.