1.2.4 Calculating Infra-Red (IR) Spectrum
It is possible to compute IR spectrum from a normal mode finite displacement calculation. Traditionally, quantum chemistry codes compute the IR intensities employing double-harmonic approximation, where the transition probability (oscillator strength) is proportional to the square of the first derivative of the total dipole moment of the molecule with respect to the normal coordinates. Here we will use the same apprach, where the normal mode derivative of the dipole moment would be computed using finite displacements, as implemented in PyEPFD. In this example, we will use CO2 molecule, however, the approach is applicable to periodic solids; see the following paper: https://onlinelibrary.wiley.com/doi/10.1002/anie.202420680 and section 17 of the associated supporting information for a recent example.
Here, the steps are very similar to the previous exercises: 1.2.1 and 1.2.2. There are two differences. First we have to instruct Qbox code to compute the total dipole monent, which can be done using the “set polarization MLWF_REF” command. Secondly, current implementation of Qbox does not compute correct dipole moment if we keep empty orbitals, therefore we have to use default nempty = 0. We are going to use the xyz file that we created in the exercise 1.2.1 but create a new input deck for Qbox calculations. Lets do it.
[2]:
%run ../../../utils/qbox_utils/xyz2qb.py enmfd_phonon.xyz
File enmfd_phonon.xyz has 3 number of atoms in each frame and a total 19 number of frames.
No of species: 2
Present species are:
O
C
The file enmfd_phonon.xyz is an i-PI output with pos_unit: angstrom and cell_unit: angstrom
Enter start frame number: (Default = 1)
Enter last frame number: (Default = 19)
Enter step frame: (Default = 1)
Enter qbox command for 1st iteration: (Default = ' randomize_wf, run -atomic_density 0 100 10')
Hint: more than one commands should be separated by a comma.
set polarization MLWF_REF, randomize_wf, run -atomic_density 0 100 10
Enter qbox run command: (Default = ' run 0 60 10')
Hint: more than one commands should be separated by a comma.
run 0 100 10
Enter qbox plot command: (Default = ' None')
Hints:
(1) Example: plot options filename (without .cube extension)
(2) More than one commands should be separated by a comma.
Enter qbox spectrum command: (Default = ' None')
Hints:
(1) Example: spectrum options filename (without .dat extension)
Do you want to save sample for each configurations (y/n)? (Default = n)
Enter filename of qbox input: (Default = qbox.i) ir_spectra-1.i
Enter xc: (Default = PBE)B3LYP
Enter wf_dyn: (Default = JD) PSDA
Enter ecut (Ry): (Default = 50.0)
Enter scf_tol (Ry): (Default = 1e-08) 1e-12
Enter nempty: (Default = 100) 0
Enter nspin (1/2): (Default = 1)
Enter delta_spin: (Default = None)
Enter net_charge: (Default = 0)
Enter pseudopotential: (Default = ONCV_PBE-1.0)
Warning!
Files:
O_ONCV_PBE-1.0.xml
C_ONCV_PBE-1.0.xml
must be present in the following pseudo-potential directory.
../pseudos/
As usual, lets have a look at the generated qbox input deck.
[3]:
%%bash
cat ir_spectra.i
species O ../pseudos/O_ONCV_PBE-1.0.xml
species C ../pseudos/C_ONCV_PBE-1.0.xml
# Frame 1/19 from enmfd_phonon.xyz#
set cell 20.000000 0.000000 0.000000 0.000000 20.000000 0.000000 0.000000 0.000000 20.000000
atom O1 O 0.000000 0.000000 2.186665
atom O2 O 0.000000 0.000000 -2.186662
atom C3 C 0.000000 0.000000 -0.000003
set ecut 50.00
set xc B3LYP
set wf_dyn PSDA
set scf_tol 1.00e-12
set polarization MLWF_REF
randomize_wf
run -atomic_density 0 100 10
# Frame 2/19 from enmfd_phonon.xyz#
set cell 20.000000 0.000000 0.000000 0.000000 20.000000 0.000000 0.000000 0.000000 20.000000
move O1 to -0.000000 -0.000000 2.475340
move O2 to -0.000000 -0.000000 -1.897987
move C3 to -0.000000 -0.000000 0.288672
run 0 100 10
# Frame 3/19 from enmfd_phonon.xyz#
set cell 20.000000 0.000000 0.000000 0.000000 20.000000 0.000000 0.000000 0.000000 20.000000
move O1 to 0.000000 0.000000 1.897990
move O2 to 0.000000 0.000000 -2.475337
move C3 to 0.000000 0.000000 -0.288678
run 0 100 10
# Frame 4/19 from enmfd_phonon.xyz#
set cell 20.000000 0.000000 0.000000 0.000000 20.000000 0.000000 0.000000 0.000000 20.000000
move O1 to 0.000000 0.448620 2.186665
move O2 to -0.000001 -0.172131 -2.186662
move C3 to -0.000000 0.138244 -0.000003
run 0 100 10
# Frame 5/19 from enmfd_phonon.xyz#
set cell 20.000000 0.000000 0.000000 0.000000 20.000000 0.000000 0.000000 0.000000 20.000000
move O1 to -0.000000 -0.448620 2.186665
move O2 to 0.000001 0.172131 -2.186662
move C3 to 0.000000 -0.138244 -0.000003
run 0 100 10
# Frame 6/19 from enmfd_phonon.xyz#
set cell 20.000000 0.000000 0.000000 0.000000 20.000000 0.000000 0.000000 0.000000 20.000000
move O1 to 0.352008 -0.000001 2.186665
move O2 to -0.355089 0.000000 -2.186662
move C3 to -0.001542 -0.000000 -0.000003
run 0 100 10
# Frame 7/19 from enmfd_phonon.xyz#
set cell 20.000000 0.000000 0.000000 0.000000 20.000000 0.000000 0.000000 0.000000 20.000000
move O1 to -0.352008 0.000001 2.186665
move O2 to 0.355089 -0.000000 -2.186662
move C3 to 0.001542 0.000000 -0.000003
run 0 100 10
# Frame 8/19 from enmfd_phonon.xyz#
set cell 20.000000 0.000000 0.000000 0.000000 20.000000 0.000000 0.000000 0.000000 20.000000
move O1 to -0.290402 -0.000001 2.186665
move O2 to -0.286941 0.000001 -2.186662
move C3 to -0.288672 -0.000000 -0.000003
run 0 100 10
# Frame 9/19 from enmfd_phonon.xyz#
set cell 20.000000 0.000000 0.000000 0.000000 20.000000 0.000000 0.000000 0.000000 20.000000
move O1 to 0.290402 0.000001 2.186665
move O2 to 0.286941 -0.000001 -2.186662
move C3 to 0.288672 0.000000 -0.000003
run 0 100 10
# Frame 10/19 from enmfd_phonon.xyz#
set cell 20.000000 0.000000 0.000000 0.000000 20.000000 0.000000 0.000000 0.000000 20.000000
move O1 to -0.000000 -0.100021 2.186665
move O2 to 0.000000 -0.416340 -2.186662
move C3 to -0.000000 -0.258181 -0.000003
run 0 100 10
# Frame 11/19 from enmfd_phonon.xyz#
set cell 20.000000 0.000000 0.000000 0.000000 20.000000 0.000000 0.000000 0.000000 20.000000
move O1 to 0.000000 0.100021 2.186665
move O2 to -0.000000 0.416340 -2.186662
move C3 to 0.000000 0.258181 -0.000003
run 0 100 10
# Frame 12/19 from enmfd_phonon.xyz#
set cell 20.000000 0.000000 0.000000 0.000000 20.000000 0.000000 0.000000 0.000000 20.000000
move O1 to 0.031448 0.010314 2.186665
move O2 to 0.031448 0.010314 -2.186662
move C3 to -0.083784 -0.027478 -0.000003
run 0 100 10
# Frame 13/19 from enmfd_phonon.xyz#
set cell 20.000000 0.000000 0.000000 0.000000 20.000000 0.000000 0.000000 0.000000 20.000000
move O1 to -0.031448 -0.010314 2.186665
move O2 to -0.031448 -0.010314 -2.186662
move C3 to 0.083784 0.027478 -0.000003
run 0 100 10
# Frame 14/19 from enmfd_phonon.xyz#
set cell 20.000000 0.000000 0.000000 0.000000 20.000000 0.000000 0.000000 0.000000 20.000000
move O1 to 0.010314 -0.031448 2.186665
move O2 to 0.010314 -0.031448 -2.186662
move C3 to -0.027478 0.083784 -0.000003
run 0 100 10
# Frame 15/19 from enmfd_phonon.xyz#
set cell 20.000000 0.000000 0.000000 0.000000 20.000000 0.000000 0.000000 0.000000 20.000000
move O1 to -0.010314 0.031448 2.186665
move O2 to -0.010314 0.031448 -2.186662
move C3 to 0.027478 -0.083784 -0.000003
run 0 100 10
# Frame 16/19 from enmfd_phonon.xyz#
set cell 20.000000 0.000000 0.000000 0.000000 20.000000 0.000000 0.000000 0.000000 20.000000
move O1 to 0.000000 0.000000 2.216857
move O2 to 0.000000 0.000000 -2.216853
move C3 to -0.000000 -0.000000 -0.000004
run 0 100 10
# Frame 17/19 from enmfd_phonon.xyz#
set cell 20.000000 0.000000 0.000000 0.000000 20.000000 0.000000 0.000000 0.000000 20.000000
move O1 to -0.000000 -0.000000 2.156473
move O2 to -0.000000 -0.000000 -2.156471
move C3 to 0.000000 0.000000 -0.000002
run 0 100 10
# Frame 18/19 from enmfd_phonon.xyz#
set cell 20.000000 0.000000 0.000000 0.000000 20.000000 0.000000 0.000000 0.000000 20.000000
move O1 to -0.000000 0.000000 2.195700
move O2 to -0.000000 0.000000 -2.177627
move C3 to 0.000000 -0.000000 -0.024074
run 0 100 10
# Frame 19/19 from enmfd_phonon.xyz#
set cell 20.000000 0.000000 0.000000 0.000000 20.000000 0.000000 0.000000 0.000000 20.000000
move O1 to 0.000000 -0.000000 2.177630
move O2 to 0.000000 -0.000000 -2.195697
move C3 to -0.000000 0.000000 0.024068
run 0 100 10
The next step is to submit those jobs and wait until all calculations finished properly. Once the calculations are finished, the obox outputs must be post-processed to compute the frequency and oscillator strengths and compute the spectra using these informations. Lets do that using the qbox output ir_spectra.r. We also need the output when we performed the normal mode displacement, i.e., the output when we ran ionic_mover class in exercise 1.2.1. We have copied it into the file enmfd.log.
Lets begin post-process all the information and compute the spectra.
[1]:
from pyepfd.vibrational_spectra import *
# Initiating the qbox_ir class
co2 = qbox_ir(file_path = 'enmfd.log',
qb_output_path = './',
qb_prefix = 'ir_spectra',
asr = 'lin',
save_path = 'co2-freq-osc-str.npz')
# calculating the frequency and oscillator strength
# using normal mode finite-displacement
co2.calc_freq_osc_str()
# Printing the frequency and oscillator strengths
co2.print_freq_osc_str()
# Saving the freq, and oscillator strengths at
# co2-freq-osc-str.npz
co2.save()
███████████
░░███░░░░░███
░███ ░███ █████ ████
░██████████ ░░███ ░███
░███░░░░░░ ░███ ░███
░███ ░███ ░███
█████ ░░███████
░░░░░ ░░░░░███
███ ░███
░░██████
░░░░░░
██████████ ███████████ ███████████ ██████████
░░███░░░░░█░░███░░░░░███░░███░░░░░░█░░███░░░░███
░███ █ ░ ░███ ░███ ░███ █ ░ ░███ ░░███
░██████ ░██████████ ░███████ ░███ ░███
░███░░█ ░███░░░░░░ ░███░░░█ ░███ ░███
░███ ░ █ ░███ ░███ ░ ░███ ███
██████████ █████ █████ ██████████
░░░░░░░░░░ ░░░░░ ░░░░░ ░░░░░░░░░░
PyEPFD version : 1.1
Author : Arpan Kundu
Author Email : arpan.kundu@gmail.com
Today : 2025-01-02 18:30:13.015826
*************************************************
CITATIONS
=================================================
Please cite the following 3 references:
(1) A. Kundu et al, Phys. Rev. Mater (2021), 5,
L070801,
(2) A. Kundu and G Galli,
J. Chem. Theory. Comput. (2023), 19, 4011
(3) https://pyepfd.readthedocs.io/en/latest/
*************************************************
Time spent on qbox class: 0.020743608474731445 s.
# Normal Mode Frequency and Dipole Oscillator Strength(IR)
#-------------------------------------------------------------------------------------------------
# Mode No Freq(cm-1) Oscillator Strength (Central/Forward/Backword) |
#------------------------------------------------------------------------------------------------
1 nan 3.0393204e-13 3.0820056e-13 3.0978919e-13
2 79.30 4.747106e-13 1.5042384e-12 1.5545001e-12
3 106.59 3.2775993e-17 3.8025423e-15 2.8191894e-15
4 nan 4.1862717e-14 7.3225319e-14 1.9810694e-14
5 nan 8.2335805e-14 1.4745235e-13 1.2574617e-13
6 676.09 1.6850272e-05 1.6848689e-05 1.6851855e-05
7 676.09 1.6850181e-05 1.6850613e-05 1.6849748e-05
8 1345.45 2.265674e-13 1.4833742e-12 1.5599035e-12
9 2361.85 0.00038054978 0.00038054979 0.00038054977
/home/arpan/.local/lib/python3.8/site-packages/pyepfd/vibrational_spectra.py:265: RuntimeWarning: invalid value encountered in sqrt
omega = np.sqrt(hessian) * ha2unit[freq_unit]
The print method would print all the modes that we have computed. As usual, the first 5 modes are not meaningful as they represent rigid body rotation and translations of the CO2 molecule. The save command would save only the last 4 modes as we have used asr = ‘lin’. Lets have a look at the data contained in co2-freq-osc-str.npz.
[2]:
import numpy as np
freq_osc_str = np.load('co2-freq-osc-str.npz')
print(freq_osc_str['omega'])
print(freq_osc_str['osc_str'])
[ 676.08784227 676.08784227 1345.44790558 2361.84521762]
[1.68502722e-05 1.68501808e-05 2.26567404e-13 3.80549776e-04]
We know that degenerate bending modes (676.0 cm-1) and asymetric stretching modes (2361.8 cm-1) are IR active as they involves change of total dipole moment, while the symmetric stretching mode (1345.4 cm-1) is not IR active. Indeed we see that the symmetric stretching mode has very low intensity. Lets use this data to compute a spectra.
The co2 class has these saved data, therefore we can directly use a method named calc_spectra() available within the co2 class. Lets do it and plot the spectra using matplotlib. We will do it using Lorentzinan broadening using three different broadening parameter of 100 cm-1.Note, if the co2-freq-osc-str.npz file is already present, then those data can be read and compute the spectra using the add_broadening() method, see the documentation for the vibrational_spectra module for more informations.
[3]:
import matplotlib.pyplot as plt
co2.calc_spectra(line_profile='Lorentzian',
line_param = 100, step = 1,
save_path = 'co2_ir.npz')
plt.xlabel("Wavenumbers (cm$^{-1}$)")
plt.ylabel("Intensities (a.u)")
plt.yticks([]) # Showing no y-ticks as units are arbitrary
plt.plot(co2.spectra[0],co2.spectra[1])
plt.show()
This would also create a co2_ir.npz, which could be later used to redraw the spectra as shown below.
[4]:
co2_spectra = np.load("co2_ir.npz")['spectra']
plt.xlabel("Wavenumbers (cm$^{-1}$)")
plt.ylabel("Intensities (a.u)")
plt.yticks([]) # Showing no y-ticks as units are arbitrary
plt.plot(co2_spectra[0],co2_spectra[1])
plt.show()
