Raman spectroscopy#

The Raman scattering process involves the inelasic scattering of light in the visible range, leaving the molecule in a different vibrational state. The associated scattering amplitude is related to the transition polarizability tensor \(\boldsymbol{\alpha}\), as discussed in more details here. In the case of normal Raman scattering described using the harmonic approximation for molecular vibrations, the transition polarizability associated to a normal mode \(k\) depends on the derivative of the electronic polarizability with respect to the normal coordinate \(Q_k\). As a consequence, only modes along which the polarizability changes will be Raman active.

Raman spectra of H2O, NH3 and CH4#

Here, we will calculate the Raman spectra of H2O, NH3 and CH4 and determine how going from O over N to C influences the vibrational frequencies of the Hydrogen stretching modes.

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

from matplotlib import pyplot as plt
import numpy as np

basis_set_label = '6-31G'
* Warning * Environment variable OMP_NUM_THREADS not set.
* Warning * Setting OMP_NUM_THREADS to 8.
* Warning * Setting MKL_THREADING_LAYER to "GNU".

To compute a meaningful Raman spectrum, the molecular geometries must be optimized at the same level of theory and using the same basis set. Let us assume that we have optimized the geometries at the HF/6-31G theory level and we have these geometries stored in xyz files. Otherwise, we would start by performing the necessary geometry optimizations.

# Prepare molecules and basis sets:
h2o_xyz = open("../../data/ir_raman/opt_h2o.xyz", "r").read() # read optimized geometry from file
nh3_xyz = open("../../data/ir_raman/opt_nh3.xyz", "r").read() # read optimized geometry from file
ch4_xyz = open("../../data/ir_raman/opt_ch4.xyz", "r").read() # read optimized geometry from file

h2o = vlx.Molecule.from_xyz_string(h2o_xyz)
h2o_basis = vlx.MolecularBasis.read(h2o, basis_set_label)

nh3 = vlx.Molecule.from_xyz_string(nh3_xyz)
nh3_basis = vlx.MolecularBasis.read(nh3, basis_set_label)

ch4 = vlx.Molecule.from_xyz_string(ch4_xyz)
ch4_basis = vlx.MolecularBasis.read(ch4, basis_set_label)

print("\n       (a) H2O                   (b) NH3                  (c) CH4")
viewer = p3d.view(viewergrid=(1,3),width=600,height=250, linked=True)
viewer.addModel(h2o_xyz, 'xyz', viewer=(0,0))
viewer.addModel(nh3_xyz, 'xyz', viewer=(0,1))
viewer.addModel(ch4_xyz, 'xyz', viewer=(0,2))
viewer.setStyle({'stick': {}})
viewer.show()
       (a) H2O                   (b) NH3                  (c) CH4

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

To compute the Raman spectrum, we need the vibrational frequencies and Raman intensities. The frequencies are eigenvalues of the molecular Hessian, while the Raman intensities are proportional to the gradient of the polarizability. The first step is to determine the Hessian matrix.

# Settings:
scf_settings = {'conv_thresh': 1.0e-5}
method_settings = {} # HF;  use 'xcfun': 'b3lyp', 'grid_level': 4, etc. for a DFT calculation 
# Numerical Hessian
hessian_settings = {'do_raman': 'yes', 'print_depolarization_ratio':'yes'}

# We first need to run SCF calculations:
h2o_scfdrv = vlx.ScfRestrictedDriver()
h2o_scfdrv.update_settings(scf_settings, method_settings)
h2o_scfdrv.ostream.state = False # To remove printout
h2o_scfdrv.compute(h2o, h2o_basis)

nh3_scfdrv = vlx.ScfRestrictedDriver()
nh3_scfdrv.update_settings(scf_settings, method_settings)
nh3_scfdrv.ostream.state = False
nh3_scfdrv.compute(nh3, nh3_basis)

ch4_scfdrv = vlx.ScfRestrictedDriver()
ch4_scfdrv.update_settings(scf_settings, method_settings)
ch4_scfdrv.ostream.state = False
ch4_scfdrv.compute(ch4, ch4_basis)
# Create a Hessian driver object, update settings, and compute 
h2o_hessian_drv = vlx.scfhessiandriver.ScfHessianDriver(h2o_scfdrv)
h2o_hessian_drv.update_settings(method_settings, hessian_settings)

nh3_hessian_drv = vlx.scfhessiandriver.ScfHessianDriver(nh3_scfdrv)
nh3_hessian_drv.update_settings(method_settings, hessian_settings)

ch4_hessian_drv = vlx.scfhessiandriver.ScfHessianDriver(ch4_scfdrv)
ch4_hessian_drv.update_settings(method_settings, hessian_settings)
# Calculate hessians from scratch
h2o_hessian_drv.compute(h2o, h2o_basis)
nh3_hessian_drv.compute(nh3, nh3_basis)
ch4_hessian_drv.compute(ch4, ch4_basis)
# Or read from checkpoint file:

import h5py

fname = '../../data/ir_raman/raman.h5'
hf = h5py.File(fname, "r")

labels = ['h2o', 'nh3', 'ch4']
i = 0
    
for driver in [h2o_hessian_drv, nh3_hessian_drv, ch4_hessian_drv]:
    label = labels[i]
    
    driver.hessian = np.array(hf.get(label+'_hessian')) 
    driver.polarizability_gradient = np.array(hf.get(label+'_polgrad')) 
    driver.dipole_gradient = np.array( hf.get(label+'_dipolegrad'))
    
    i += 1

hf.close()
h2o_hessian_drv.vibrational_analysis(h2o)
                                                   Vibrational Analysis                                                   
                                                  ======================                                                  
                                                                                                                          
                 Harmonic frequencies (in cm**-1), force constants (in mdyne/A), reduced masses (in amu),                 
                          IR intensities (in km/mol), Raman scattering activities (in A**4/amu),                          
                      parallel and perpendicular Raman scattering activities, depolarization ratios,                      
                                         and Cartesian normal mode displacements.                                         
                                                                                                                          
                                                                                                                          
  Index:                        1                              2                              3               
  Frequency:                 1736.87                        3988.20                        4145.10            
  Force constant:            1.9404                         9.7200                         11.0231            
  Reduced mass:              1.0917                         1.0372                         1.0889             
  IR intensity:             123.0358                        2.9544                         54.2936            
  Raman activ.:              10.6274                        90.1326                        40.0871            
  Parallel Raman:            97.0717                       943.7930                       292.0906            
  Perp. Raman:               38.4408                       205.5058                       219.0680            
  Depol. ratio:              0.3960                         0.2177                         0.7500             
  Normal mode:       X         Y         Z     |    X         Y         Z     |    X         Y         Z     |
  1 O             -0.0554   -0.0410    0.0290  |  0.0327    0.0242   -0.0171  |  0.0158   -0.0544   -0.0467  |
  2 H              0.5217    0.0427   -0.4725  | -0.1273   -0.6466   -0.2544  |  0.1681    0.6495    0.2172  |
  3 H              0.3571    0.6079    0.0128  | -0.3919    0.2622    0.5260  | -0.4197    0.2144    0.5246  |
                                                                                                                          
                                                                                                                          
nh3_hessian_drv.vibrational_analysis(nh3)
                                                   Vibrational Analysis                                                   
                                                  ======================                                                  
                                                                                                                          
                 Harmonic frequencies (in cm**-1), force constants (in mdyne/A), reduced masses (in amu),                 
                          IR intensities (in km/mol), Raman scattering activities (in A**4/amu),                          
                      parallel and perpendicular Raman scattering activities, depolarization ratios,                      
                                         and Cartesian normal mode displacements.                                         
                                                                                                                          
                                                                                                                          
  Index:                        1                              2                              3               
  Frequency:                 597.47                         1814.46                        1814.47            
  Force constant:            0.2525                         2.1054                         2.1054             
  Reduced mass:              1.2003                         1.0854                         1.0854             
  IR intensity:             644.8631                        46.9030                        46.9030            
  Raman activ.:              4.6612                         6.2677                         6.2677             
  Parallel Raman:            58.6274                        45.6689                        45.6690            
  Perp. Raman:               0.8086                         34.2517                        34.2517            
  Depol. ratio:              0.0138                         0.7500                         0.7500             
  Normal mode:       X         Y         Z     |    X         Y         Z     |    X         Y         Z     |
  1 N              0.0794    0.0113    0.0915  |  0.0081   -0.0767    0.0024  | -0.0579   -0.0045    0.0508  |
  2 H             -0.4452   -0.0243   -0.3601  |  0.1611    0.7186   -0.1654  | -0.1089    0.2583   -0.1202  |
  3 H             -0.3156    0.0200   -0.4779  | -0.3086    0.4013    0.3445  |  0.3386   -0.4046   -0.0974  |
  4 H             -0.3420   -0.1531   -0.4336  |  0.0343   -0.0540   -0.2130  |  0.5753    0.2093   -0.4884  |
                                                                                                                          
  Index:                        4                              5                              6               
  Frequency:                 3779.47                        3983.60                        3983.60            
  Force constant:            8.5217                         10.2768                        10.2768            
  Reduced mass:              1.0125                         1.0991                         1.0991             
  IR intensity:              0.5571                         16.1431                        16.1430            
  Raman activ.:             109.6648                        39.5442                        39.5443            
  Parallel Raman:           1223.8453                      288.1351                       288.1354            
  Perp. Raman:              174.5125                       216.1014                       216.1015            
  Depol. ratio:              0.1426                         0.7500                         0.7500             
  Normal mode:       X         Y         Z     |    X         Y         Z     |    X         Y         Z     |
  1 N              0.0122    0.0017    0.0141  | -0.0518    0.0536    0.0383  |  0.0367    0.0639   -0.0398  |
  2 H              0.3677   -0.1623   -0.4143  |  0.4485   -0.2013   -0.5554  | -0.1983    0.1137    0.2440  |
  3 H             -0.3412   -0.4046    0.2305  |  0.0379    0.0729   -0.0208  | -0.4994   -0.5646    0.2943  |
  4 H             -0.1966    0.5426   -0.0122  |  0.2330   -0.6159    0.0444  |  0.1873   -0.4371    0.0143  |
                                                                                                                          
                                                                                                                          
ch4_hessian_drv.vibrational_analysis(ch4)
                                                   Vibrational Analysis                                                   
                                                  ======================                                                  
                                                                                                                          
                 Harmonic frequencies (in cm**-1), force constants (in mdyne/A), reduced masses (in amu),                 
                          IR intensities (in km/mol), Raman scattering activities (in A**4/amu),                          
                      parallel and perpendicular Raman scattering activities, depolarization ratios,                      
                                         and Cartesian normal mode displacements.                                         
                                                                                                                          
                                                                                                                          
  Index:                        1                              2                              3               
  Frequency:                 1516.76                        1516.77                        1516.77            
  Force constant:            1.6021                         1.6021                         1.6021             
  Reduced mass:              1.1820                         1.1820                         1.1820             
  IR intensity:              21.4495                        21.4494                        21.4495            
  Raman activ.:              2.9988                         2.9987                         2.9988             
  Parallel Raman:            21.8503                        21.8500                        21.8504            
  Perp. Raman:               16.3878                        16.3875                        16.3878            
  Depol. ratio:              0.7500                         0.7500                         0.7500             
  Normal mode:       X         Y         Z     |    X         Y         Z     |    X         Y         Z     |
  1 C              0.1208    0.0145    0.0316  |  0.0074    0.1011   -0.0744  |  0.0340   -0.0734   -0.0963  |
  2 H             -0.3408    0.2769    0.0631  |  0.0868   -0.3096    0.4673  | -0.3867    0.0305    0.2628  |
  3 H             -0.3837   -0.3044   -0.2252  | -0.2187   -0.2678    0.4258  |  0.1738   -0.0487    0.3378  |
  4 H             -0.3008    0.0648   -0.3929  |  0.2984   -0.3257    0.0707  |  0.0638    0.4618    0.2676  |
  5 H             -0.4146   -0.2098    0.1782  | -0.2543   -0.3013   -0.0769  | -0.2559    0.4307    0.2793  |
                                                                                                                          
  Index:                        4                              5                              6               
  Frequency:                 1708.70                        1708.70                        3181.65            
  Force constant:            1.7339                         1.7339                         6.0118             
  Reduced mass:              1.0080                         1.0080                         1.0080             
  IR intensity:              0.0000                         0.0000                         0.0000             
  Raman activ.:              39.5165                        39.5165                       145.7008            
  Parallel Raman:           287.9329                       287.9329                       1857.8604           
  Perp. Raman:              215.9497                       215.9497                        0.0000             
  Depol. ratio:              0.7500                         0.7500                         0.0000             
  Normal mode:       X         Y         Z     |    X         Y         Z     |    X         Y         Z     |
  1 C              0.0000    0.0000   -0.0000  | -0.0000   -0.0000   -0.0000  | -0.0000   -0.0000   -0.0000  |
  2 H             -0.4146    0.2792    0.0130  |  0.1175    0.1950   -0.4452  |  0.2537    0.3660    0.2273  |
  3 H             -0.3848   -0.2703   -0.1699  | -0.0637   -0.1956    0.4557  |  0.3128   -0.3724   -0.1161  |
  4 H              0.3506   -0.1577    0.3197  |  0.0127   -0.4426   -0.2322  | -0.3562   -0.1710    0.3063  |
  5 H              0.4487    0.1487   -0.1628  | -0.0664    0.4432    0.2217  | -0.2103    0.1773   -0.4175  |
                                                                                                                          
  Index:                        7                              8                              9               
  Frequency:                 3296.17                        3296.17                        3296.18            
  Force constant:            7.0367                         7.0367                         7.0367             
  Reduced mass:              1.0993                         1.0993                         1.0993             
  IR intensity:              37.8001                        37.8000                        37.8000            
  Raman activ.:              66.7339                        66.7338                        66.7338            
  Parallel Raman:           486.2497                       486.2491                       486.2492            
  Perp. Raman:              364.6872                       364.6868                       364.6869            
  Depol. ratio:              0.7500                         0.7500                         0.7500             
  Normal mode:       X         Y         Z     |    X         Y         Z     |    X         Y         Z     |
  1 C             -0.0340   -0.0610    0.0584  | -0.0033   -0.0620   -0.0666  | -0.0844    0.0270   -0.0209  |
  2 H              0.1655    0.2357    0.1716  |  0.3808    0.5344    0.3245  |  0.1386    0.2390    0.1385  |
  3 H             -0.0733    0.0606    0.0392  | -0.3633    0.4152    0.1170  |  0.3920   -0.4859   -0.1592  |
  4 H              0.5517    0.2531   -0.4667  | -0.1206   -0.0737    0.0854  |  0.2417    0.1337   -0.2324  |
  5 H             -0.2385    0.1776   -0.4405  |  0.1425   -0.1372    0.2672  |  0.2337   -0.2087    0.5026  |
                                                                                                                          
                                                                                                                          
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)))

        # 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 Raman spectra
plt.figure(figsize=(7,4))

h2o_x,h2o_y = h2o_hessian_drv.frequencies, h2o_hessian_drv.raman_intensities
nh3_x,nh3_y = nh3_hessian_drv.frequencies, nh3_hessian_drv.raman_intensities
ch4_x,ch4_y = ch4_hessian_drv.frequencies, ch4_hessian_drv.raman_intensities

h2o_xl, h2o_yl = add_broadening(h2o_x, h2o_y, line_profile='Lorentzian', line_param=10, step=10)
nh3_xl, nh3_yl = add_broadening(nh3_x, nh3_y, line_profile='Lorentzian', line_param=10, step=10)
ch4_xl, ch4_yl = add_broadening(ch4_x, ch4_y, line_profile='Lorentzian', line_param=10, step=10)

plt.plot(h2o_xl, h2o_yl, label='H2O')
plt.plot(nh3_xl, nh3_yl, label='NH3')
plt.plot(ch4_xl, ch4_yl, label='CH4')

plt.xlabel('Wavenumber (cm**-1)')
plt.ylabel('Raman activity (A**4/amu)')
plt.title("Calculated Raman sepctra of H2O, NH3, CH4")
plt.legend()
plt.tight_layout(); plt.show()
../../_images/vib_ir_raman_14_0.png

Let’s compare the H-stretching modes:

# 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):
        # add coordinates:
        vib_xyz += elements[i] + " %15.7f %15.7f %15.7f " % (coords[i,0], coords[i,1], coords[i,2])
        # add displacements:
        vib_xyz += "%15.7f %15.7f %15.7f\n" % (nm[i,0], nm[i,1], nm[i,2])
    return vib_xyz
h2o_h1 = get_normal_mode(h2o, h2o_hessian_drv.normal_modes[-2])
h2o_h2 = get_normal_mode(h2o, h2o_hessian_drv.normal_modes[-1])

nh3_h1 = get_normal_mode(nh3, nh3_hessian_drv.normal_modes[-3])
nh3_h2 = get_normal_mode(nh3, nh3_hessian_drv.normal_modes[-2])
nh3_h3 = get_normal_mode(nh3, nh3_hessian_drv.normal_modes[-1])

ch4_h1 = get_normal_mode(ch4, ch4_hessian_drv.normal_modes[-4])
ch4_h2 = get_normal_mode(ch4, ch4_hessian_drv.normal_modes[-3])
ch4_h3 = get_normal_mode(ch4, ch4_hessian_drv.normal_modes[-2])
ch4_h4 = get_normal_mode(ch4, ch4_hessian_drv.normal_modes[-1])
print("These are the two H stretching modes of H2O.")
view = p3d.view(viewergrid=(1,2), width=300, height=200, linked=True)
view.addModel(h2o_h1, "xyz", {'vibrate': {'frames':10,'amplitude':0.75}}, viewer=(0,0))
view.addModel(h2o_h2, "xyz", {'vibrate': {'frames':10,'amplitude':0.75}}, viewer=(0,1))
view.setStyle({'stick':{}})
view.animate({'loop': 'backAndForth'})
view.rotate(120, "x")
view.zoomTo()
view.show()
These are the two H stretching modes of H2O.

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("These are the three H stretching modes of NH3.")
view = p3d.view(viewergrid=(1,3), width=600, height=200, linked=True)
view.addModel(nh3_h1, "xyz", {'vibrate': {'frames':10,'amplitude':0.75}}, viewer=(0,0))
view.addModel(nh3_h2, "xyz", {'vibrate': {'frames':10,'amplitude':0.75}}, viewer=(0,1))
view.addModel(nh3_h3, "xyz", {'vibrate': {'frames':10,'amplitude':0.75}}, viewer=(0,2))
view.setStyle({'stick':{}})
view.animate({'loop': 'backAndForth'})
view.rotate(90, "x")
view.zoomTo()
view.show()
These are the three H stretching modes of NH3.

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("These are the four H stretching modes of CH4.")
view = p3d.view(viewergrid=(1,4), width=600, height=200, linked=True)
view.addModel(ch4_h1, "xyz", {'vibrate': {'frames':10,'amplitude':0.75}}, viewer=(0,0))
view.addModel(ch4_h2, "xyz", {'vibrate': {'frames':10,'amplitude':0.75}}, viewer=(0,1))
view.addModel(ch4_h3, "xyz", {'vibrate': {'frames':10,'amplitude':0.75}}, viewer=(0,2))
view.addModel(ch4_h4, "xyz", {'vibrate': {'frames':10,'amplitude':0.75}}, viewer=(0,3))
view.setStyle({'stick':{}})
view.animate({'loop': 'backAndForth'})
view.zoomTo()
view.show()
These are the four H stretching modes of CH4.

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