# IR spectroscopy#

The infrared (IR) absorption process involves the absorption of light with low photon energy (from approx. 800 nm to approx. 1000 $$\mu$$m) which promotes the excitation of molecular vibrations. Of the entire IR photon energy range, the photons which excite fundamental vibrations of covalent bonds have energies in the range 2.5-25 $$\mu$$m [Cen15, NRS18]. The absorbed photon energies (i.e., IR absorption peak positions) correspond to the molecular vibrational frequencies, while the IR intensities are related to the IR linear absorption cross-section [NRS18]:

$\begin{equation*} \sigma_k(\omega)= \frac{\pi}{3\epsilon_0c}\left(\frac{\partial \mu}{\partial Q_k}\right)^2\frac{\gamma_{k}}{(\omega_{k0}-\omega)^2+\gamma_{k}^2}\,. \end{equation*}$

Here, $$\mu$$ is the molecular dipole moment of the electronic ground state, $$Q_k$$ is the normal coordinate corresponding to the vibrational mode $$k$$, $$\omega$$ is the angular frequency of the absorbed electromagnetic radiation, $$\omega_{k0}$$ is the angular frequency difference between the excited vibrational state $$|k\rangle$$ and the ground state, and $$\gamma_{k}$$ is the half-width broadening associated with the inverse lifetime of $$|k\rangle$$. As can be seen from the equation, the IR absorption cross-section depends on the derivative of the dipole moment with respect to the normal coordinate $$Q_k$$. This means that the only modes which have non-zero IR absorption cross-sections are those associated with a change in the dipole moment [SH07].

These molecular vibrations constitute fingerprints of different functional groups present in the system and thus are important for chemical characterization and molecular identification.

The first step in an IR spectrum calculation is to determine the normal modes and frequencies. These are obtained from the eigenstates and eigenvalues of the Hessian matrix which collects all the second order derivatives of the energy $$E$$ with respect to the nuclear coordinates. For the ground state, the Hessian can be determined analytically, or numerically based on the analytical gradient. The dipole moment gradient, required to compute IR absorption cross-sections, can also be determined analytically.

# Import section
import veloxchem as vlx
import numpy as np
import sys
import py3Dmol as p3d
from veloxchem.veloxchemlib import bohr_in_angstroms

# Define the molecule and basis set
molecule_string = """3

O 0.000 0.000  0.000
H 0.000 0.000  0.950
H 0.896 0.000 -0.317"""
basis_set_label = 'sto-3g'

molecule = vlx.Molecule.from_xyz_string(molecule_string)
basis = vlx.MolecularBasis.read(molecule, basis_set_label)

# to view the molecules
view = p3d.view(linked=True, viewergrid=(1,1),width=300,height=200)
view.setStyle({'stick': {}})
view.rotate(-90, "x")
view.zoomTo()
view.show()


You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

# Run the SCF calculation
scf_settings = {}
method_settings = {} #{'xcfun':'b3lyp'}
scfdrv = vlx.ScfRestrictedDriver()
scfdrv.update_settings(scf_settings, method_settings)
scfdrv.compute(molecule, basis)


Self Consistent Field Driver Setup
====================================

Wave Function Model             : Spin-Restricted Hartree-Fock
Initial Guess Model             : Superposition of Atomic Densities
Convergence Accelerator         : Two Level Direct Inversion of Iterative Subspace
Max. Number of Iterations       : 50
Max. Number of Error Vectors    : 10
Convergence Threshold           : 1.0e-06
ERI Screening Scheme            : Cauchy Schwarz + Density
ERI Screening Mode              : Dynamic
ERI Screening Threshold         : 1.0e-12
Linear Dependence Threshold     : 1.0e-06

* Info * Nuclear repulsion energy: 9.2514793147 a.u.

* Info * Overlap matrix computed in 0.00 sec.

* Info * Kinetic energy matrix computed in 0.00 sec.

* Info * Nuclear potential matrix computed in 0.00 sec.


* Info * Orthogonalization matrix computed in 0.06 sec.


* Info * SAD initial guess computed in 0.00 sec.


* Info * Starting Reduced Basis SCF calculation...
* Info * ...done. SCF energy in reduced basis set: -74.959319307971 a.u. Time: 0.08 sec.


* Info * Overlap matrix computed in 0.00 sec.

* Info * Kinetic energy matrix computed in 0.00 sec.

* Info * Nuclear potential matrix computed in 0.00 sec.


* Info * Orthogonalization matrix computed in 0.00 sec.



Iter. | Hartree-Fock Energy | Energy Change | Gradient Norm | Max. Gradient | Density Change
--------------------------------------------------------------------------------------------
1       -74.959319310077    0.0000000000      0.00002657      0.00000643      0.00000000

                  2       -74.959319310295   -0.0000000002      0.00000889      0.00000203      0.00002226

                  3       -74.959319310327   -0.0000000000      0.00000049      0.00000011      0.00001225


*** SCF converged in 3 iterations. Time: 0.01 sec.


               Spin-Restricted Hartree-Fock:
-----------------------------
Total Energy                       :      -74.9593193103 a.u.
Electronic Energy                  :      -84.2107986251 a.u.
Nuclear Repulsion Energy           :        9.2514793147 a.u.
------------------------------------
Gradient Norm                      :        0.0000004869 a.u.

Ground State Information
------------------------
Charge of Molecule            :  0.0
Multiplicity (2S+1)           :  1.0
Magnetic Quantum Number (M_S) :  0.0

Spin Restricted Orbitals
------------------------

Molecular Orbital No.   1:
--------------------------
Occupation: 2.000 Energy:  -20.23340 a.u.
(   1 O   1s  :     0.99)

Molecular Orbital No.   2:
--------------------------
Occupation: 2.000 Energy:   -1.26571 a.u.
(   1 O   1s  :    -0.23) (   1 O   2s  :     0.83) (   2 H   1s  :     0.16)
(   3 H   1s  :     0.16)

Molecular Orbital No.   3:
--------------------------
Occupation: 2.000 Energy:   -0.62927 a.u.
(   1 O   1p+1:     0.35) (   1 O   1p0 :    -0.49) (   2 H   1s  :    -0.44)
(   3 H   1s  :     0.44)

Molecular Orbital No.   4:
--------------------------
Occupation: 2.000 Energy:   -0.44167 a.u.
(   1 O   2s  :    -0.52) (   1 O   1p+1:     0.65) (   1 O   1p0 :     0.46)
(   2 H   1s  :     0.27) (   3 H   1s  :     0.27)

Molecular Orbital No.   5:
--------------------------
Occupation: 2.000 Energy:   -0.38765 a.u.
(   1 O   1p-1:     1.00)

Molecular Orbital No.   6:
--------------------------
Occupation: 0.000 Energy:    0.60284 a.u.
(   1 O   2s  :     0.91) (   1 O   1p+1:     0.58) (   1 O   1p0 :     0.41)
(   2 H   1s  :    -0.81) (   3 H   1s  :    -0.81)

Molecular Orbital No.   7:
--------------------------
Occupation: 0.000 Energy:    0.76592 a.u.
(   1 O   1p+1:    -0.58) (   1 O   1p0 :     0.82) (   2 H   1s  :    -0.84)
(   3 H   1s  :     0.84)



Now we can calculate the Hessian

# Create the Hessian Driver and use it to compute an analytical SCF Hessian
hessian_settings = {}
hessian_drv = vlx.scfhessiandriver.ScfHessianDriver(scfdrv)
hessian_drv.update_settings(method_settings, hessian_settings)
hessian_drv.compute(molecule, basis)


SCF Hessian Driver
====================


                                   *** Time spent in Hessian calculation: 12.81 sec ***


hessian_drv.print_hessian(molecule)

                                           Numerical Hessian (Hartree/Bohr**2)
-------------------------------------

Coord.       1 O(x)          1 O(y)          1 O(z)          2 H(x)          2 H(y)          2 H(z)

1 O(x)       0.77341797     -0.00000043     -0.22069621     -0.05793220     -0.00000004      0.05915246
1 O(y)      -0.00000043     -0.08462381     -0.00000045     -0.00000004      0.04250770     -0.00000040
1 O(z)      -0.22069621     -0.00000045      0.93182676     -0.07097811     -0.00000050     -0.79549537
2 H(x)      -0.05793220     -0.00000004     -0.07097811      0.05352700     -0.00000006     -0.01075576
2 H(y)      -0.00000004      0.04250770     -0.00000050     -0.00000006     -0.02940283     -0.00000045
2 H(z)       0.05915246     -0.00000040     -0.79549537     -0.01075576     -0.00000045      0.80288069
3 H(x)      -0.71548574     -0.00000038      0.29167301      0.00440526     -0.00000020     -0.04839772
3 H(y)      -0.00000052      0.04211599     -0.00000015     -0.00000014     -0.01310488      0.00000001
3 H(z)       0.16154235     -0.00000009     -0.13633144      0.08173279     -0.00000013     -0.00738539

Coord.       3 H(x)          3 H(y)          3 H(z)

1 O(x)      -0.71548574     -0.00000052      0.16154235
1 O(y)      -0.00000038      0.04211599     -0.00000009
1 O(z)       0.29167301     -0.00000015     -0.13633144
2 H(x)       0.00440526     -0.00000014      0.08173279
2 H(y)      -0.00000020     -0.01310488     -0.00000013
2 H(z)      -0.04839772      0.00000001     -0.00738539
3 H(x)       0.71108043     -0.00000044     -0.24327534
3 H(y)      -0.00000044     -0.02901132     -0.00000002
3 H(z)      -0.24327534     -0.00000002      0.14371682



The normal modes and frequencies are then determined by diagonalizing the Hessian. First, the translation and rotational modes must be projected out, leaving $$3N-6$$ (or $$3N-5$$ for a linear molecule) vibrational modes for a system with $$N$$ atoms. This is done under the hood with the help of geomeTRIC.

The IR intensities are determined as discussed in more detail here.

hessian_drv.vibrational_analysis(molecule)

                                                   Vibrational Analysis
======================

Harmonic frequencies (in cm**-1), force constants (in mdyne/A), reduced masses (in amu),
IR intensities (in km/mol),
and Cartesian normal mode displacements.

Index:                        1                              2                              3
Frequency:                 1927.08                        4545.78                        4896.44
Force constant:            2.3750                         12.6969                        15.3533
Reduced mass:              1.0855                         1.0429                         1.0869
IR intensity:              15.0373                        30.3498                        12.8846
Normal mode:       X         Y         Z     |    X         Y         Z     |    X         Y         Z     |
1 O              0.0587    0.0000    0.0415  | -0.0397    0.0000   -0.0274  |  0.0416    0.0000   -0.0595  |
2 H             -0.7054    0.0000    0.0092  | -0.0299    0.0000    0.7000  | -0.0002    0.0000    0.7111  |
3 H             -0.2266   -0.0000   -0.6677  |  0.6605   -0.0000   -0.2654  | -0.6594    0.0000    0.2331  |


def add_broadening(list_ex_energy, list_osci_strength, line_profile='Lorentzian', line_param=10, step=10):
x_min = np.amin(list_ex_energy) - 50
x_max = np.amax(list_ex_energy) + 50
x = np.arange(x_min, x_max, step)
y = np.zeros((len(x)))
#print(x)
#print(y)

# go through the frames and calculate the spectrum for each frame
for xp in range(len(x)):
for e, f in zip(list_ex_energy, list_osci_strength):
if line_profile == 'Gaussian':
y[xp] += f * np.exp(-(
(e - x[xp]) / line_param)**2)
elif line_profile == 'Lorentzian':
y[xp] += 0.5 * line_param * f / (np.pi * (
(x[xp] - e)**2 + 0.25 * line_param**2))
return x, y

# plot the IR spectrum
from matplotlib import pyplot as plt

plt.figure(figsize=(7,4))
# Spectrum from -CHF
x1,y1 = hessian_drv.frequencies, hessian_drv.ir_intensities
plt.plot(x1i,y1i)
#plt.plot(x1,y1,'x')
plt.xlabel('wavenumber (cm**-1)')
plt.ylabel('IR intensity (km/mol)')
plt.title("Calculated IR sepctrum of the water molecule")
plt.tight_layout(); plt.show()


Let’s now visualize the three vibrational normal modes. The two high frequency modes are OH stretching modes (one symmetryc and one asymmetric), while the lower energy mode is a bond angle bending.

# To animate the normal mode we will need both the geometry and the displacements
def get_normal_mode(molecule, normal_mode):
elements = molecule.get_labels()
coords = molecule.get_coordinates() * bohr_in_angstroms() # To transform from au to A
natm = molecule.number_of_atoms()
vib_xyz = "%d\n\n" % natm
nm = normal_mode.reshape(natm, 3)
for i in range(natm):
vib_xyz += elements[i] + " %15.7f %15.7f %15.7f " % (coords[i,0], coords[i,1], coords[i,2])
vib_xyz += "%15.7f %15.7f %15.7f\n" % (nm[i,0], nm[i,1], nm[i,2])
return vib_xyz

vib_1 = get_normal_mode(molecule, hessian_drv.normal_modes[0])
vib_2 = get_normal_mode(molecule, hessian_drv.normal_modes[1])
vib_3 = get_normal_mode(molecule, hessian_drv.normal_modes[2])

print("This is the bending mode at 1927.09 cm-1.")
view = p3d.view(width=300, height=300)
view.addModel(vib_1, "xyz", {'vibrate': {'frames':10,'amplitude':0.75}})
view.setStyle({'stick':{}})
view.animate({'loop': 'backAndForth'})
view.rotate(-90, "x")
view.zoomTo()
view.show()

This is the bending mode at 1927.09 cm-1.


You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

print("This is the symmetric stretching mode at 4545.78 cm-1.")
view = p3d.view(width=300, height=300)
view.addModel(vib_2, "xyz", {'vibrate': {'frames':10,'amplitude':0.75}})
view.setStyle({'stick':{}})
view.animate({'loop': 'backAndForth'})
view.rotate(-90, "x")
view.zoomTo()
view.show()

This is the symmetric stretching mode at 4545.78 cm-1.


You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

print("This is the asymmetric stretching mode at 4896.44 cm-1.")
view = p3d.view(width=300, height=300)
view.addModel(vib_3, "xyz", {'vibrate': {'frames':10,'amplitude':0.75}})
view.setStyle({'stick':{}})
view.animate({'loop': 'backAndForth'})
view.rotate(-90, "x")
view.zoomTo()
view.show()

This is the asymmetric stretching mode at 4896.44 cm-1.


You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol