Exercises#

In this exercise we will look at how substitution affects the vibrational frequencies, IR, and Raman spectra.

# import section
import veloxchem as vlx
import py3Dmol as p3d
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".
# define molecules and basis sets
ethene_xyz = """6
 ethene
 C          0.000000    -0.663984   0.000000
 C          0.000000     0.663984   0.000000
 H          0.919796    -1.223061   0.000000
 H         -0.919796    -1.223061   0.000000
 H          0.919796     1.223061   0.000000
 H         -0.919796     1.223061   0.000000
"""
ethene = vlx.Molecule.from_xyz_string(ethene_xyz)
ethene_basis = vlx.MolecularBasis.read(ethene, basis_set_label)
fluoroethene_xyz = """6
fluoroethene
 C          0.000000    -0.663984   0.000000
 C          0.000000     0.663984   0.000000
 F         1.519796    -1.223061   0.000000
 H         -0.919796    -1.223061   0.000000
 H          0.919796     1.223061   0.000000
 H         -0.919796     1.223061   0.000000
 """
fluoroethene = vlx.Molecule.from_xyz_string(fluoroethene_xyz)
fluoroethene_basis = vlx.MolecularBasis.read(fluoroethene, basis_set_label)
chloroethene_xyz = """6
chloroethene
 C          0.000000    -0.663984   0.000000
 C          0.000000     0.663984   0.000000
 Cl         1.519796    -1.223061   0.000000
 H         -0.919796    -1.223061   0.000000
 H          0.919796     1.223061   0.000000
 H         -0.919796     1.223061   0.000000
 """
chloroethene = vlx.Molecule.from_xyz_string(chloroethene_xyz)
chloroethene_basis = vlx.MolecularBasis.read(chloroethene, basis_set_label)
# to view the molecules
view = p3d.view(linked=True, viewergrid=(1,3),width=600,height=200)
view.addModel(ethene_xyz, 'xyz', viewer=(0,0))
view.addModel(fluoroethene_xyz, 'xyz', viewer=(0,1))
view.addModel(chloroethene_xyz, 'xyz', viewer=(0,2))
view.setStyle({'stick': {}})
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

SCF geometry optimization#

Before we can calculate the vibrational spectra, we first must optimize the geometries.

# Settings for SCF and gradient drivers
scf_settings = {'conv_thresh':1e-6}
method_settings = {}
# Run SCF for ethene
ethene_scfdrv = ...
...

# Run SCF for fluoroethene
fluoroethene_scfdrv = ...
...

# Run SCF for chloroethene
chloroethene_scfdrv = ...
...
# Run SCF for ethene
ethene_scfdrv = vlx.ScfRestrictedDriver()
ethene_scfdrv.update_settings(scf_settings, method_settings)
ethene_scfdrv.ostream.state = False 
ethene_scfdrv.compute(ethene, ethene_basis)

# Run SCF for fluoroethene
fluoroethene_scfdrv = vlx.ScfRestrictedDriver()
fluoroethene_scfdrv.update_settings(scf_settings, method_settings)
fluoroethene_scfdrv.ostream.state = False 
fluoroethene_scfdrv.compute(fluoroethene, fluoroethene_basis)

# Run SCF for chloroethene
chloroethene_scfdrv = vlx.ScfRestrictedDriver()
chloroethene_scfdrv.update_settings(scf_settings, method_settings)
chloroethene_scfdrv.ostream.state = False 
chloroethene_scfdrv.compute(chloroethene, chloroethene_basis)
# Set up the gradient and optimization dirvers:
ethene_grad_drv = ...

fluoroethene_grad_drv = ...
...

chloroethene_grad_drv = ...
...
# Set up the gradient and optimization dirvers:
ethene_grad_drv = vlx.ScfGradientDriver(ethene_scfdrv)
ethene_opt_drv = vlx.OptimizationDriver(ethene_grad_drv)

fluoroethene_grad_drv = vlx.ScfGradientDriver(fluoroethene_scfdrv)
fluoroethene_opt_drv = vlx.OptimizationDriver(fluoroethene_grad_drv)

chloroethene_grad_drv = vlx.ScfGradientDriver(chloroethene_scfdrv)
chloroethene_opt_drv = vlx.OptimizationDriver(chloroethene_grad_drv)
# Optimize the geometries
opt_ethene = ...
opt_fluoroethene = ...
opt_chloroethene = ...
# Either optimize the geometries
# Uncomment below
#ethene_opt_drv.ostream.state = True # to enable printout
#opt_ethene = ethene_opt_drv.compute(ethene, ethene_basis)
#fluoroethene_opt_drv.ostream.state = True # to enable printout
#opt_fluoroethene = fluoroethene_opt_drv.compute(fluoroethene, fluoroethene_basis)
#chloroethene_opt_drv.ostream.state = True
#opt_chloroethene = chloroethene_opt_drv.compute(chloroethene, chloroethene_basis)

# Or read them from file if already calculated
opt_ethene_xyz = open("../../data/ir_raman/opt_ethene.xyz", "r").read() # read optimized geometry from file
opt_fluoroethene_xyz = open("../../data/ir_raman/opt_fluoroethene.xyz", "r").read() # read optimized geometry from file
opt_chloroethene_xyz = open("../../data/ir_raman/opt_chloroethene.xyz", "r").read() # read optimized geometry from file

opt_ethene = vlx.Molecule.from_xyz_string(opt_ethene_xyz)
ethene_basis = vlx.MolecularBasis.read(opt_ethene, basis_set_label)

opt_fluoroethene = vlx.Molecule.from_xyz_string(opt_fluoroethene_xyz)
fluoroethene_basis = vlx.MolecularBasis.read(opt_fluoroethene, basis_set_label)

opt_chloroethene = vlx.Molecule.from_xyz_string(opt_chloroethene_xyz)
chloroethene_basis = vlx.MolecularBasis.read(opt_chloroethene, basis_set_label)
# Get optimized coordinates as xyz string
def get_xyz(molecule):
    natm = molecule.number_of_atoms()
    elements = molecule.get_labels()
    coords = molecule.get_coordinates() * bohr_in_angstroms()
    txt = "%d\n\n" % natm
    for i in range(natm):
        txt += elements[i] + " %15.7f %15.7f %15.7f\n" % (coords[i,0], coords[i,1], coords[i,2])
    return txt
# Visualize the optimized structures
...
# Compute SCF with the optimized geometries of all molecules
...
ethene_scfdrv.compute(opt_ethene, ethene_basis)
fluoroethene_scfdrv.compute(opt_fluoroethene, fluoroethene_basis)
chloroethene_scfdrv.compute(opt_chloroethene, chloroethene_basis)

Hessians#

Now, we can calculate the IR and Raman spectra for these optimized geometries.

# Settings for Hessian calculation
hessian_settings = {'do_raman': 'yes', 'print_depolarization_ratio':'no'}
# Create Hessian driver and update settings
ethene_hessian_drv = ...
...

fluoroethene_hessian_drv = ...
...

chloroethene_hessian_drv = ...
...
# Compute the Hessians:
...
...
...
# Create Hessian driver and update settings
ethene_hessian_drv = vlx.scfhessiandriver.ScfHessianDriver(ethene_scfdrv)
ethene_hessian_drv.update_settings(method_settings, hessian_settings)

fluoroethene_hessian_drv = vlx.scfhessiandriver.ScfHessianDriver(fluoroethene_scfdrv)
fluoroethene_hessian_drv.update_settings(method_settings, hessian_settings)

chloroethene_hessian_drv = vlx.scfhessiandriver.ScfHessianDriver(chloroethene_scfdrv)
chloroethene_hessian_drv.update_settings(method_settings, hessian_settings)

# Calculate
#ethene_hessian_drv.compute(ethene, ethene_basis)
#fluoroethene_hessian_drv.compute(fluoroethene, fluoroethene_basis)
#chloroethene_hessian_drv.compute(chloroethene, chloroethene_basis)

# Or read from file:
import h5py

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

labels = ['ethene', 'fluoroethene', 'chloroethene']

i = 0   
for driver in [ethene_hessian_drv, fluoroethene_hessian_drv, chloroethene_hessian_drv]:
    label = labels[i]

    driver.hessian = np.array(hf.get(label + '_hessian')) 
    driver.dipole_gradient = np.array( hf.get(label + '_dipolegrad'))
    driver.polarizability_gradient = np.array(hf.get(label + '_polgrad')) 
    
    i += 1

hf.close()
# Broadening function
def add_broadening(list_ex_energy, list_osci_strength, line_profile='Lorentzian', line_param=10, step=10):

    ...
    
    return x, y
# Broadening function
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
# To animate the normal mode we will need both the geometry and the displacements 
def get_normal_mode(molecule, normal_mode):

    ...
    
    return vib_xyz
# 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

Vibrational analysis#

To get a summary of the vibrational analysis, one simply needs to run the following command:

# Run for all molecules
ethene_hessian_drv.vibrational_analysis(opt_ethene)
...
...
# Run for all molecules
ethene_hessian_drv.vibrational_analysis(opt_ethene)
fluoroethene_hessian_drv.vibrational_analysis(opt_fluoroethene)
chloroethene_hessian_drv.vibrational_analysis(opt_chloroethene)
                                                   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),                          
                                         and Cartesian normal mode displacements.                                         
                                                                                                                          
                                                                                                                          
  Index:                        1                              2                              3               
  Frequency:                 927.96                         1120.00                        1154.18            
  Force constant:            0.5295                         0.8579                         1.1920             
  Reduced mass:              1.0436                         1.1608                         1.5188             
  IR intensity:              1.1152                        140.0931                        0.0003             
  Raman activ.:              0.0000                         0.0000                         5.6658             
  Normal mode:       X         Y         Z     |    X         Y         Z     |    X         Y         Z     |
  1 C             -0.0403    0.0000   -0.0000  | -0.0000   -0.0000   -0.0832  | -0.0000   -0.0000   -0.1525  |
  2 C             -0.0402    0.0000    0.0000  |  0.0000   -0.0000   -0.0835  |  0.0000    0.0000    0.1522  |
  3 H              0.2398    0.4378   -0.0001  |  0.0000    0.0001    0.4958  | -0.0000   -0.0000    0.4803  |
  4 H              0.2398   -0.4379    0.0002  | -0.0000    0.0001    0.4961  |  0.0000   -0.0001    0.4976  |
  5 H              0.2399   -0.4379    0.0001  | -0.0000   -0.0000    0.4973  | -0.0000    0.0001   -0.4788  |
  6 H              0.2397    0.4378   -0.0002  |  0.0000    0.0000    0.4969  | -0.0001   -0.0002   -0.4960  |
                                                                                                                          
  Index:                        4                              5                              6               
  Frequency:                 1157.33                        1378.66                        1512.66            
  Force constant:            0.7955                         1.7072                         1.7014             
  Reduced mass:              1.0080                         1.5245                         1.2620             
  IR intensity:              0.0000                         0.0000                         0.0000             
  Raman activ.:              0.0011                         1.2235                         50.5279            
  Normal mode:       X         Y         Z     |    X         Y         Z     |    X         Y         Z     |
  1 C              0.0000    0.0000    0.0018  |  0.1532    0.0000   -0.0000  |  0.0000   -0.1074   -0.0000  |
  2 C              0.0000   -0.0000   -0.0017  | -0.1532    0.0000    0.0000  | -0.0000    0.1075   -0.0000  |
  3 H             -0.0001   -0.0001   -0.5057  | -0.1454   -0.4660    0.0001  | -0.1897   -0.4570   -0.0000  |
  4 H             -0.0001    0.0001    0.4942  | -0.1452    0.4659    0.0001  |  0.1897   -0.4570    0.0001  |
  5 H             -0.0001    0.0001    0.5054  |  0.1455   -0.4662   -0.0001  | -0.1890    0.4559    0.0001  |
  6 H             -0.0000   -0.0001   -0.4945  |  0.1452    0.4658   -0.0001  |  0.1890    0.4561    0.0000  |
                                                                                                                          
  Index:                        7                              8                              9               
  Frequency:                 1631.37                        1846.11                        3312.91            
  Force constant:            1.7426                         5.7467                         6.7803             
  Reduced mass:              1.1113                         2.8619                         1.0485             
  IR intensity:              8.4043                         0.0000                         20.9734            
  Raman activ.:              0.0001                         16.3748                        0.0000             
  Normal mode:       X         Y         Z     |    X         Y         Z     |    X         Y         Z     |
  1 C             -0.0000    0.0687    0.0000  | -0.0000   -0.2902    0.0000  | -0.0000   -0.0429    0.0000  |
  2 C              0.0000    0.0683   -0.0000  |  0.0000    0.2903   -0.0000  | -0.0000   -0.0429   -0.0000  |
  3 H             -0.2845   -0.4077    0.0000  |  0.3822    0.2482   -0.0001  | -0.4284    0.2557   -0.0000  |
  4 H              0.2844   -0.4079   -0.0000  | -0.3821    0.2482   -0.0000  |  0.4284    0.2556    0.0000  |
  5 H              0.2846   -0.4085    0.0001  |  0.3825   -0.2486    0.0000  |  0.4289    0.2560    0.0000  |
  6 H             -0.2846   -0.4089    0.0001  | -0.3824   -0.2488   -0.0000  | -0.4286    0.2557    0.0000  |
                                                                                                                          
  Index:                       10                             11                             12               
  Frequency:                 3337.47                        3389.59                        3419.12            
  Force constant:            7.0900                         7.5393                         7.6928             
  Reduced mass:              1.0803                         1.1138                         1.1169             
  IR intensity:              0.0000                         0.0000                         43.0914            
  Raman activ.:             188.7189                       129.4466                        0.0001             
  Normal mode:       X         Y         Z     |    X         Y         Z     |    X         Y         Z     |
  1 C             -0.0001    0.0574   -0.0000  |  0.0694    0.0000    0.0000  | -0.0703   -0.0000    0.0000  |
  2 C             -0.0000   -0.0573   -0.0000  | -0.0693   -0.0000    0.0000  | -0.0704    0.0000   -0.0000  |
  3 H              0.4273   -0.2579   -0.0000  | -0.4198    0.2668    0.0000  |  0.4187   -0.2679    0.0000  |
  4 H             -0.4264   -0.2573    0.0000  | -0.4206   -0.2671   -0.0000  |  0.4192    0.2680   -0.0000  |
  5 H              0.4266    0.2574    0.0000  |  0.4199    0.2668   -0.0000  |  0.4189    0.2680    0.0000  |
  6 H             -0.4266    0.2573   -0.0000  |  0.4197   -0.2666    0.0000  |  0.4197   -0.2684    0.0000  |
                                                                                                                          
                                                                                                                          
                                                   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),                          
                                         and Cartesian normal mode displacements.                                         
                                                                                                                          
                                                                                                                          
  Index:                        1                              2                              3               
  Frequency:                 512.97                         784.04                         1005.42            
  Force constant:            0.4032                         0.5229                         1.5365             
  Reduced mass:              2.6004                         1.4438                         2.5798             
  IR intensity:              6.3843                         0.1889                         57.3584            
  Raman activ.:              2.5116                         10.6177                        7.0448             
  Normal mode:       X         Y         Z     |    X         Y         Z     |    X         Y         Z     |
  1 C              0.1705    0.1375    0.0000  | -0.0000   -0.0000    0.1791  | -0.1539    0.0459    0.0000  |
  2 C             -0.1372    0.1298   -0.0000  |  0.0000    0.0000   -0.0295  | -0.1648    0.1215   -0.0000  |
  3 F              0.0184   -0.1924   -0.0000  | -0.0000    0.0000   -0.0639  |  0.1793   -0.1175   -0.0000  |
  4 H              0.1047    0.2063   -0.0000  |  0.0000   -0.0000   -0.2629  |  0.0300   -0.1591   -0.0000  |
  5 H             -0.4757    0.5827    0.0000  |  0.0000   -0.0000    0.4926  |  0.2408   -0.4046    0.0000  |
  6 H             -0.3722   -0.3480   -0.0000  | -0.0000   -0.0000   -0.8070  |  0.1476    0.7853    0.0000  |
                                                                                                                          
  Index:                        4                              5                              6               
  Frequency:                 1090.60                        1109.86                        1252.17            
  Force constant:            0.7784                         1.0061                         2.0623             
  Reduced mass:              1.1107                         1.3863                         2.2324             
  IR intensity:             114.0030                        11.9584                        76.4288            
  Raman activ.:              0.0109                         7.4209                         3.1899             
  Normal mode:       X         Y         Z     |    X         Y         Z     |    X         Y         Z     |
  1 C             -0.0000    0.0000    0.0452  | -0.0000    0.0000   -0.1187  | -0.2563   -0.0283   -0.0000  |
  2 C             -0.0000    0.0000    0.0854  | -0.0000    0.0000    0.1424  |  0.1309    0.0303    0.0000  |
  3 F              0.0000   -0.0000   -0.0001  |  0.0000   -0.0000    0.0031  |  0.1267   -0.0175   -0.0000  |
  4 H              0.0000   -0.0000   -0.6233  | -0.0000    0.0000    0.6197  | -0.5223    0.2431    0.0000  |
  5 H              0.0000   -0.0000   -0.7558  |  0.0000   -0.0000   -0.2348  | -0.2515    0.5254    0.0000  |
  6 H              0.0000    0.0000   -0.1760  |  0.0000    0.0000   -0.7256  | -0.1202   -0.4621   -0.0000  |
                                                                                                                          
  Index:                        7                              8                              9               
  Frequency:                 1458.95                        1569.48                        1879.33            
  Force constant:            1.4750                         1.7011                         9.4714             
  Reduced mass:              1.1762                         1.1721                         4.5515             
  IR intensity:              0.3677                         3.5929                         79.5500            
  Raman activ.:              18.4529                        9.4538                         26.5674            
  Normal mode:       X         Y         Z     |    X         Y         Z     |    X         Y         Z     |
  1 C             -0.0576   -0.0450   -0.0000  | -0.0537    0.1020   -0.0000  |  0.1088   -0.3973   -0.0000  |
  2 C              0.0593    0.0739    0.0000  |  0.0139    0.0320    0.0000  | -0.0399    0.3857    0.0000  |
  3 F             -0.0238   -0.0043    0.0000  |  0.0081   -0.0137    0.0000  | -0.0239    0.0261    0.0000  |
  4 H              0.5432   -0.7157    0.0000  |  0.2438   -0.2107    0.0000  | -0.4357    0.0594    0.0000  |
  5 H             -0.1594    0.3814   -0.0000  |  0.3754   -0.4771   -0.0000  |  0.4774   -0.1785   -0.0000  |
  6 H              0.0452    0.0702   -0.0000  | -0.2980   -0.6512    0.0000  | -0.4123   -0.2358   -0.0000  |
                                                                                                                          
  Index:                       10                             11                             12               
  Frequency:                 3363.62                        3448.47                        3466.15            
  Force constant:            7.0856                         7.7686                         7.8601             
  Reduced mass:              1.0629                         1.1088                         1.1104             
  IR intensity:              0.1141                         1.2178                         8.5554             
  Raman activ.:              71.3931                       104.7873                        30.4417            
  Normal mode:       X         Y         Z     |    X         Y         Z     |    X         Y         Z     |
  1 C             -0.0099    0.0024   -0.0000  |  0.0450    0.0515   -0.0000  |  0.0412    0.0436   -0.0000  |
  2 C              0.0079   -0.0695    0.0000  | -0.0637   -0.0207    0.0000  |  0.0754   -0.0041    0.0000  |
  3 F              0.0001    0.0000    0.0000  |  0.0002   -0.0005   -0.0000  | -0.0006    0.0002    0.0000  |
  4 H              0.1018    0.0905    0.0000  | -0.5519   -0.5081    0.0000  | -0.4739   -0.4411    0.0000  |
  5 H              0.5590    0.4109   -0.0000  |  0.4220    0.3176    0.0000  | -0.3781   -0.2912   -0.0000  |
  6 H             -0.6377    0.2974   -0.0000  |  0.3480   -0.1671   -0.0000  | -0.5260    0.2588   -0.0000  |
                                                                                                                          
                                                                                                                          
                                                   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),                          
                                         and Cartesian normal mode displacements.                                         
                                                                                                                          
                                                                                                                          
  Index:                        1                              2                              3               
  Frequency:                 421.74                         691.77                         700.54             
  Force constant:            0.3291                         1.3261                         0.3966             
  Reduced mass:              3.1400                         4.7034                         1.3715             
  IR intensity:              0.7508                         55.9808                        16.1679            
  Raman activ.:              8.4334                         17.1578                        7.8631             
  Normal mode:       X         Y         Z     |    X         Y         Z     |    X         Y         Z     |
  1 C             -0.1491   -0.1752   -0.0000  |  0.4327   -0.0957   -0.0000  |  0.0000   -0.0000    0.1717  |
  2 C              0.2191   -0.1774    0.0000  |  0.0775   -0.1359    0.0000  |  0.0000   -0.0000   -0.0353  |
  3 Cl            -0.0448    0.1326    0.0000  | -0.1706    0.0873    0.0000  | -0.0000    0.0000   -0.0272  |
  4 H             -0.1848   -0.1269    0.0000  |  0.3495    0.0163    0.0000  |  0.0000    0.0000   -0.3394  |
  5 H              0.4837   -0.5633   -0.0000  | -0.2375    0.3227   -0.0000  | -0.0000    0.0000    0.4668  |
  6 H              0.4411    0.2278    0.0000  | -0.1926   -0.6499    0.0000  | -0.0000   -0.0000   -0.7971  |
                                                                                                                          
  Index:                        4                              5                              6               
  Frequency:                 1100.39                        1143.81                        1154.49            
  Force constant:            0.7758                         1.0717                         1.0512             
  Reduced mass:              1.0875                         1.3903                         1.3386             
  IR intensity:              67.9295                        45.1802                        15.9610            
  Raman activ.:              1.6158                         4.9232                         5.5759             
  Normal mode:       X         Y         Z     |    X         Y         Z     |    X         Y         Z     |
  1 C              0.0000    0.0000   -0.0847  | -0.0000   -0.0000   -0.0728  | -0.1105   -0.0089    0.0000  |
  2 C             -0.0000    0.0000   -0.0079  |  0.0000   -0.0000    0.1716  |  0.1277   -0.0212   -0.0000  |
  3 Cl            -0.0000   -0.0000   -0.0004  |  0.0000    0.0000    0.0008  |  0.0176    0.0032   -0.0000  |
  4 H              0.0000   -0.0000    0.7988  | -0.0000    0.0000    0.1611  | -0.3999    0.3532    0.0000  |
  5 H              0.0000   -0.0000    0.5488  | -0.0000    0.0000   -0.6441  | -0.2366    0.4901    0.0000  |
  6 H              0.0000    0.0000   -0.2313  | -0.0000   -0.0000   -0.7242  | -0.1888   -0.5951    0.0000  |
                                                                                                                          
  Index:                        7                              8                              9               
  Frequency:                 1419.98                        1558.93                        1836.24            
  Force constant:            1.4193                         1.7025                         8.2324             
  Reduced mass:              1.1947                         1.1890                         4.1440             
  IR intensity:              33.9957                        9.0960                         43.5328            
  Raman activ.:              12.1840                        23.4378                        21.4675            
  Normal mode:       X         Y         Z     |    X         Y         Z     |    X         Y         Z     |
  1 C              0.0841    0.0038   -0.0000  | -0.0198    0.1260    0.0000  | -0.0467    0.3587    0.0000  |
  2 C             -0.0804   -0.0558    0.0000  |  0.0005   -0.0132   -0.0000  |  0.0161   -0.3922    0.0000  |
  3 Cl             0.0098    0.0004    0.0000  |  0.0009   -0.0024   -0.0000  |  0.0026   -0.0049   -0.0000  |
  4 H             -0.5137    0.7707    0.0000  |  0.1379   -0.0498    0.0000  |  0.3274    0.0208   -0.0000  |
  5 H              0.0910   -0.3090    0.0000  |  0.3675   -0.5855    0.0000  | -0.5031    0.2471   -0.0000  |
  6 H              0.0343    0.1427   -0.0000  | -0.3054   -0.6252   -0.0000  |  0.4493    0.3022   -0.0000  |
                                                                                                                          
  Index:                       10                             11                             12               
  Frequency:                 3346.74                        3439.19                        3460.90            
  Force constant:            7.0153                         7.7706                         7.7719             
  Reduced mass:              1.0631                         1.1150                         1.1013             
  IR intensity:              0.4596                         0.6629                         7.5225             
  Raman activ.:             101.1155                        91.7141                        71.1723            
  Normal mode:       X         Y         Z     |    X         Y         Z     |    X         Y         Z     |
  1 C              0.0076   -0.0037   -0.0000  |  0.0211    0.0237   -0.0000  | -0.0621   -0.0573   -0.0000  |
  2 C             -0.0156    0.0685   -0.0000  | -0.0913   -0.0197    0.0000  | -0.0352    0.0100    0.0000  |
  3 Cl             0.0000   -0.0001    0.0000  |  0.0002   -0.0002    0.0000  |  0.0003   -0.0000    0.0000  |
  4 H             -0.0877   -0.0679    0.0000  | -0.2703   -0.2121    0.0000  |  0.7289    0.5784   -0.0000  |
  5 H             -0.5091   -0.3390    0.0000  |  0.6232    0.4334   -0.0000  |  0.1725    0.1238   -0.0000  |
  6 H              0.6908   -0.3624    0.0000  |  0.4768   -0.2621    0.0000  |  0.2480   -0.1381    0.0000  |
                                                                                                                          
                                                                                                                          

Think about the dipole moment of these molecules and consider the Hydrogen stretching modes. How do you expect the IR spectra to look like? Which molecule do you expect will have more intense IR-peaks?

# Plot the IR spectra
plt.figure(figsize=(7,4))

eth_x, eth_ir = ethene_hessian_drv.frequencies, ethene_hessian_drv.ir_intensities
flo_x, flo_ir = fluoroethene_hessian_drv.frequencies, fluoroethene_hessian_drv.ir_intensities
chl_x, chl_ir = chloroethene_hessian_drv.frequencies, chloroethene_hessian_drv.ir_intensities

eth_xl, eth_irl = add_broadening(eth_x, eth_ir, line_profile='Lorentzian', line_param=20, step=2)
flo_xl, flo_irl = add_broadening(flo_x, flo_ir, line_profile='Lorentzian', line_param=20, step=2)
chl_xl, chl_irl = add_broadening(chl_x, chl_ir, line_profile='Lorentzian', line_param=20, step=2)

plt.plot(eth_xl, eth_irl, label='Ethene')
plt.plot(flo_xl, flo_irl, label='Fluoroethene')
plt.plot(chl_xl, chl_irl, label='Chloroethene')

plt.xlabel('Wavenumber (cm**-1)')
plt.axis(xmin=3200, xmax=3500)
plt.axis(ymin=-0.2, ymax=1.5)
plt.ylabel('IR intensity (km/mol)')
plt.title("Calculated IR sepctra, H-stretching region")
plt.legend()
plt.tight_layout(); plt.show()
../../_images/vib_ir_exercises_31_0.png
# Plot the Raman spectra
plt.figure(figsize=(7,4))

eth_x, eth_raman = ethene_hessian_drv.frequencies, ethene_hessian_drv.raman_intensities
flo_x, flo_raman = fluoroethene_hessian_drv.frequencies, fluoroethene_hessian_drv.raman_intensities
chl_x, chl_raman = chloroethene_hessian_drv.frequencies, chloroethene_hessian_drv.raman_intensities

eth_xl, eth_ramanl = add_broadening(eth_x, eth_raman, line_profile='Lorentzian', line_param=20, step=2)
flo_xl, flo_ramanl = add_broadening(flo_x, flo_raman, line_profile='Lorentzian', line_param=20, step=2)
chl_xl, chl_ramanl = add_broadening(chl_x, chl_raman, line_profile='Lorentzian', line_param=20, step=2)

plt.plot(eth_xl, eth_ramanl, label='Ethene')
plt.plot(flo_xl, flo_ramanl, label='Fluoroethene')
plt.plot(chl_xl, chl_ramanl, label='Chloroethene')

plt.xlabel('Wavenumber (cm**-1)')
plt.axis(xmin=3200, xmax=3500)
plt.axis(ymin=-0.2, ymax=7)
plt.ylabel('Raman activity (A**4/amu)')
plt.title("Calculated Raman sepctra, H-stretching region")
plt.legend()
plt.tight_layout(); plt.show()
../../_images/vib_ir_exercises_32_0.png
# Get the displacements of the normal mode
ethene_h1 = get_normal_mode(ethene, ethene_hessian_drv.normal_modes[-1])
fluoroethene_h1 = get_normal_mode(fluoroethene, fluoroethene_hessian_drv.normal_modes[-1])
chloroethene_h1 = get_normal_mode(chloroethene, chloroethene_hessian_drv.normal_modes[-1])
# Animate the vibration
view = p3d.view(viewergrid=(1,3), width=600, height=200, linked=True)
view.addModel(ethene_h1, "xyz", {'vibrate': {'frames':10,'amplitude':0.75}}, viewer=(0,0))
view.addModel(fluoroethene_h1, "xyz", {'vibrate': {'frames':10,'amplitude':0.75}}, viewer=(0,1))
view.addModel(chloroethene_h1, "xyz", {'vibrate': {'frames':10,'amplitude':0.75}}, viewer=(0,2))
view.setStyle({'stick':{}})
view.animate({'loop': 'backAndForth'})
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

IR intensities and Raman activities#

To rationalize why the IR and Raman spectra look as they do, calculate how the dipole moment and polarizability change during particular vibrational motions. Look at the Hydrogen stretching modes and select a mode which is IR-active in ethene, but is suppressed in fluoroethene or chloroethene. What is the dipole moment in the optimized molecule? How does the dipole moment change during the vibration?

# Calculate the dipole moment of the optimized molecules
# For this we will use the FirstOrderProperties class from veloxchem
ethene_prop = vlx.firstorderprop.FirstOrderProperties()
ethene_prop.compute_scf_prop(opt_ethene, ethene_basis, ethene_scfdrv.scf_tensors)
ethene_dipole_moment = ethene_prop.get_property('dipole moment')
...

# Select normal mode and get the array of atomic displacements
...

# Use the atomic displacement array to construct several new molecular configurations
# along the vibrational mode, e.g. -0.75, -0.5, -0.25, 0.25, 0.5, 0.75 displacement 
...

# Calculate the dipole moment for the new configurations
...

# Plot as a function of displacement
...
ethene_prop = vlx.firstorderprop.FirstOrderProperties()
ethene_prop.compute_scf_prop(opt_ethene, ethene_basis, ethene_scfdrv.scf_tensors)
ethene_dipole_moment = ethene_prop.get_property('dipole moment')

fluoroethene_prop = vlx.firstorderprop.FirstOrderProperties()
fluoroethene_prop.compute_scf_prop(opt_fluoroethene, fluoroethene_basis, fluoroethene_scfdrv.scf_tensors)
fluoroethene_dipole_moment = fluoroethene_prop.get_property('dipole moment')

chloroethene_prop = vlx.firstorderprop.FirstOrderProperties()
chloroethene_prop.compute_scf_prop(opt_chloroethene, chloroethene_basis, chloroethene_scfdrv.scf_tensors)
chloroethene_dipole_moment = chloroethene_prop.get_property('dipole moment')

print("Ground state dipole moments:   x             y            z       \n")
print("Ethene                     : %5.2f a.u.    %5.2f a.u.   %5.2f a.u." % (ethene_dipole_moment[0],
                                                           ethene_dipole_moment[1],
                                                           ethene_dipole_moment[2] ))
print("Fluoroethene               : %5.2f a.u.    %5.2f a.u.   %5.2f a.u." % (fluoroethene_dipole_moment[0],
                                                           fluoroethene_dipole_moment[1],
                                                           fluoroethene_dipole_moment[2] ))
print("Chloroethene               : %5.2f a.u.    %5.2f a.u.   %5.2f a.u." % (chloroethene_dipole_moment[0],
                                                           chloroethene_dipole_moment[1],
                                                           chloroethene_dipole_moment[2] ))
print()

ethene_moments = []
fluoroethene_moments = []

natm = opt_ethene.number_of_atoms()

new_ethene_scfdrv = vlx.scfrestdriver.ScfRestrictedDriver()
new_ethene_scfdrv.ostream.state=False
ethene_displacements = ethene_hessian_drv.normal_modes[-1].reshape(natm, 3)
ethene_labels = opt_ethene.get_labels()

new_fluoroethene_scfdrv = vlx.scfrestdriver.ScfRestrictedDriver()
new_fluoroethene_scfdrv.ostream.state=False
fluoroethene_displacements = fluoroethene_hessian_drv.normal_modes[-1].reshape(natm, 3)
fluoroethene_labels = opt_fluoroethene.get_labels()

x_list = np.arange(-0.8,0.81,0.2)

# Calculate the dipole moment for different displacements
for x in x_list:
    ethene_coords = opt_ethene.get_coordinates() * bohr_in_angstroms()
    ethene_coords += x * ethene_displacements
    new_ethene = vlx.molecule.Molecule(ethene_labels, ethene_coords, units='angstrom')
    new_ethene_scfdrv.compute(new_ethene, ethene_basis)
    ethene_prop.compute_scf_prop(new_ethene, ethene_basis, new_ethene_scfdrv.scf_tensors)
    ethene_dipole_moment = ethene_prop.get_property('dipole moment')
    ethene_moments.append(ethene_dipole_moment)

for x in x_list:
    fluoro_coords = opt_fluoroethene.get_coordinates() * bohr_in_angstroms()
    fluoro_coords += x * fluoroethene_displacements
    new_fluoroethene = vlx.molecule.Molecule(fluoroethene_labels, fluoro_coords, units='angstrom')
    new_fluoroethene_scfdrv.compute(new_fluoroethene, fluoroethene_basis)
    ethene_prop.compute_scf_prop(new_fluoroethene, fluoroethene_basis, new_fluoroethene_scfdrv.scf_tensors)
    fluoroethene_dipole_moment = fluoroethene_prop.get_property('dipole moment')
    fluoroethene_moments.append(fluoroethene_dipole_moment)
    

ethene_moments_array = np.array(ethene_moments)
fluoroethene_moments_array = np.array(fluoroethene_moments)

# Plot
figure = plt.figure(figsize=(12,4))

plt1 = figure.add_subplot(1, 2, 1)
plt.plot(x_list, ethene_moments_array[:,0], '-o', label='x')
plt.plot(x_list, ethene_moments_array[:,1], '-o', label='y')
plt.plot(x_list, ethene_moments_array[:,2], '--o', label='z')
plt.axis(xmin=-0.85, xmax=0.85)
plt.axis(ymin=-0.85, ymax=0.85)
plt.xlabel('Displacement (%)')
plt.ylabel('Dipole moment (a.u.)')
plt.title("Ethene Dipole Moment along the %.2f cm-1 mode" % ethene_hessian_drv.frequencies[-1])
plt.legend()

plt2 = figure.add_subplot(1, 2, 2)
plt.plot(x_list, fluoroethene_moments_array[:,0], '-o', label='x')
plt.plot(x_list, fluoroethene_moments_array[:,1], '-o', label='y')
plt.plot(x_list, fluoroethene_moments_array[:,2], '--o', label='z')
plt.axis(xmin=-0.85, xmax=0.85)
plt.axis(ymin=-0.85, ymax=0.85)
plt.xlabel('Displacement (%)')
plt.ylabel('Dipole moment (a.u.)')
plt.title("Fluoroethene dipole moment along the %.2f cm-1 mode" % fluoroethene_hessian_drv.frequencies[-1])
plt.legend()

plt.tight_layout(); plt.show()
Ground state dipole moments:   x             y            z       

Ethene                     : -0.00 a.u.    -0.00 a.u.   -0.00 a.u.
Fluoroethene               : -0.79 a.u.     0.39 a.u.    0.00 a.u.
Chloroethene               : -0.77 a.u.     0.59 a.u.   -0.00 a.u.
../../_images/vib_ir_exercises_37_1.png

Now select a mode which is Raman-active in ethene, but is suppressed in fluoroethene or chloroethene. Calculate the polarizability of the optimized molecule. How does the polarizability change during the vibration?

# Calculate the polarizability of the optimized molecules
# For this, we need to run a linear response calculation
ethene_lrdrv = vlx.lrsolver.LinearResponseSolver()
ethene_pol_dict = ethene_lrdrv.compute(opt_ethene, ethene_basis, ethene_scfdrv.scf_tensors)
ethene_polarizability = ethene_pol_dict['response_functions']
...

# Select normal mode and get the array of atomic displacements
...

# Use the atomic displacement array to construct several new molecular configurations
# along the vibrational mode, e.g. -0.75, -0.5, -0.25, 0.25, 0.5, 0.75 displacement 
...

# Calculate the polarizability for the new configurations
...

# Plot component or norm as a function of displacement
...
# Calculate the polarizability of the optimized molecules
# For this, we need to run a linear response calculation
ethene_lrdrv = vlx.lrsolver.LinearResponseSolver()
ethene_lrdrv.ostream.state = False
ethene_pol_dict = ethene_lrdrv.compute(opt_ethene, ethene_basis, ethene_scfdrv.scf_tensors)
ethene_polarizability = ethene_pol_dict['response_functions']

fluoroethene_lrdrv = vlx.lrsolver.LinearResponseSolver()
fluoroethene_lrdrv.ostream.state = False
fluoroethene_pol_dict = fluoroethene_lrdrv.compute(opt_fluoroethene, fluoroethene_basis, fluoroethene_scfdrv.scf_tensors)
fluoroethene_polarizability = fluoroethene_pol_dict['response_functions']

chloroethene_lrdrv = vlx.lrsolver.LinearResponseSolver()
chloroethene_lrdrv.ostream.state = False
chloroethene_pol_dict = chloroethene_lrdrv.compute(opt_chloroethene, chloroethene_basis, chloroethene_scfdrv.scf_tensors)
chloroethene_polarizability = chloroethene_pol_dict['response_functions']

print("Ground state polarizabilities:    xx             yy           zz\n")
print("Ethene                       : %5.2f a.u.    %5.2f a.u.   %5.2f a.u."  % (ethene_polarizability[('x','x',0)],
                                                           ethene_polarizability[('y','y',0)],
                                                           ethene_polarizability[('z','z',0)]))
print("Fluoroethene                 : %5.2f a.u.    %5.2f a.u.   %5.2f a.u." % (fluoroethene_polarizability[('x','x',0)],
                                                           fluoroethene_polarizability[('y','y',0)],
                                                           fluoroethene_polarizability[('z','z',0)]))
print("Chloroethene                 : %5.2f a.u.    %5.2f a.u.   %5.2f a.u." % (chloroethene_polarizability[('x','x',0)],
                                                           chloroethene_polarizability[('y','y',0)],
                                                           chloroethene_polarizability[('z','z',0)]))
print()

# Save polarizability dictionary keys
keys = ethene_pol_dict['response_functions'].keys()

ethene_pol = []
fluoroethene_pol = []

# Create new linear response solvers to re-calculate the polarizabilities
new_ethene_lrdrv = vlx.lrsolver.LinearResponseSolver()
new_ethene_lrdrv.ostream.state = False

new_fluoroethene_lrdrv = vlx.lrsolver.LinearResponseSolver()
new_fluoroethene_lrdrv.ostream.state = False

# Calculate the polarizability for different displacements
for x in x_list:
    ethene_coords = opt_ethene.get_coordinates() * bohr_in_angstroms()
    ethene_coords += x * ethene_displacements
    new_ethene = vlx.molecule.Molecule(ethene_labels, ethene_coords, units='angstrom')
    new_ethene_scfdrv.compute(new_ethene, ethene_basis)

    new_ethene_lrdrv.is_converged = False
    new_ethene_pol_dict = new_ethene_lrdrv.compute(new_ethene, ethene_basis, new_ethene_scfdrv.scf_tensors)
    
    ethene_polarizability = np.zeros((9))
    i = 0
    for key in keys:
        ethene_polarizability[i] = new_ethene_pol_dict['response_functions'][key]
        i += 1
    ethene_pol.append(ethene_polarizability.reshape(3,3))

for x in x_list:
    fluoro_coords = opt_fluoroethene.get_coordinates() * bohr_in_angstroms()
    fluoro_coords += x * fluoroethene_displacements
    new_fluoroethene = vlx.molecule.Molecule(fluoroethene_labels, fluoro_coords, units='angstrom')
    new_fluoroethene_scfdrv.compute(new_fluoroethene, fluoroethene_basis)

    new_fluoroethene_lrdrv.is_converged = False
    new_fluoroethene_pol_dict = new_fluoroethene_lrdrv.compute(new_fluoroethene, fluoroethene_basis, new_fluoroethene_scfdrv.scf_tensors)
    
    fluoroethene_polarizability = np.zeros((9))
    i = 0
    for key in keys:
        fluoroethene_polarizability[i] = new_fluoroethene_pol_dict['response_functions'][key]
        i += 1
    fluoroethene_pol.append(fluoroethene_polarizability.reshape(3,3))

ethene_polarizability_array = np.array(ethene_pol)
fluoroethene_polarizability_array = np.array(fluoroethene_pol)

figure = plt.figure(figsize=(12,4))

# Plot
plt1 = figure.add_subplot(1, 2, 1)
plt.plot(x_list, ethene_polarizability_array[:,0,0], '-o', label='xx')
plt.plot(x_list, ethene_polarizability_array[:,1,1], '-o', label='yy')
plt.plot(x_list, ethene_polarizability_array[:,2,2], '--o', label='zz')
plt.axis(xmin=-0.85, xmax=0.85)
plt.axis(ymin=-32.0, ymax=2.0)
plt.xlabel('Displacement (%)')
plt.ylabel('Polarizability (a.u.)')
plt.title("Ethene polarizability along the %.2f cm-1 mode " % (ethene_hessian_drv.frequencies[-1]))
plt.legend()

plt2 = figure.add_subplot(1, 2, 2)
plt.plot(x_list, fluoroethene_polarizability_array[:,0,0], '-o', label='xx')
plt.plot(x_list, fluoroethene_polarizability_array[:,1,1], '-o', label='yy')
plt.plot(x_list, fluoroethene_polarizability_array[:,2,2], '--o', label='zz')
plt.axis(xmin=-0.85, xmax=0.85)
plt.axis(ymin=-32.0, ymax=2.0)
plt.xlabel('Displacement (%)')
plt.ylabel('Polarizability (a.u.)')
plt.title("Fluoroethene polarizability along the %.2f cm-1 mode " % (fluoroethene_hessian_drv.frequencies[-1]))
plt.legend()

plt.tight_layout(); plt.show()
Ground state polarizabilities:    xx             yy           zz

Ethene                       : -19.19 a.u.    -32.45 a.u.   -7.15 a.u.
Fluoroethene                 : -18.86 a.u.    -30.62 a.u.   -7.81 a.u.
Chloroethene                 : -31.94 a.u.    -41.05 a.u.   -10.12 a.u.
../../_images/vib_ir_exercises_40_1.png