VibrationalSpectroscopy - Maple Help

Home : Support : Online Help : Toolboxes : Quantum Chemistry : Courses : Lessons : VibrationalSpectroscopy

Vibrational Spectroscopy: Calculating the Rovibrational Spectrum of a Diatomic Molecule

 Overview In this exercise, you will calculate the rovibrational spectrum of a diatomic molecule of interest (e.g.HCl) using different levels of theory: a normal mode analysis and ab initio methods.     Normal Mode Analysis:     As a first approximation, you will use the harmonic oscillator / rigid rotor approximation to calculate the fundamental vibrational frequency, ${v}_{0},$ and rotational energies. In this approximation,  the total rovibrational energy is given by Eq. (1):        (1)     where n = 0, 1, 2, ... is the vibrational quantum number, l = 0, 1, 2, ... is the rotational quantum number, ${R}_{0}$is the equilibrium bond distance, and μ is the reduced mass of the diatomic.   Using the Bohr correspondence principle, one can calculate a given transition between any two rovibrational states.  However, selection rules for a rotating/vibrating heteronuclear diatomic are Δl = ±1 and Δn = ±1.  Therefore, the (0,0) --> (1,0) fundamental transition would not be observed. What would be observed are the P-branch (corresponding to all transitions with Δl = -1) and the R-branch (all transitions corresponding to Δl = +1):   P-branch:    = (2) R-branch:    = (3) Transition energies are typically expressed in wavenumbers (${\mathrm{cm}}^{-1}$):   .        (4)    You will use the GeometryOptimization function to calculate ${R}_{0}$ for the diatomic of interest and then use the VibrationalModes function to calculate the fundamental frequency, $\stackrel{~}{v}$You will then calculate the rotational energies given by the second term in Eq. (1) for l = 0, 1, ..., ${l}_{\mathrm{max}}$.  You will then calculate P- and R-branches according to Eqs. (2) and (3) and compare to experimental values.       Ab Initio Calculation   While the above approximation is simple, more sophisticated methods can be used to capture the inherent anharmonicity in the system. In the second portion of this activity, you will calculate the potential energy surface (PES), V(R), using density functional theory (DFT) and the 6-31g basis set for a range of internuclear distances, R.   Once V(R) has been constructed, the nuclear Schrodinger equation is given by   $\stackrel{\wedge }{H}$ (R) =  (R).           (5)   where is the corresponding rovibrational wavefunction, and the Hamiltonian, $\stackrel{\wedge }{H}$, is given by $\stackrel{\wedge }{H}$ = -(6)   While the details of how one solves Eq. (5) are beyond the scope of this exercise, the variational principle is used:  represent Eq. (5) in an underlying basis representation for Ψ and diagonalize the resulting matrix to get energy eigenvalues, ${E}_{n,l}$.  (For more on matrix representations and the variational principle, see Variational Theorem exercise).  Once ${E}_{n,l}$ values have been calculated, the transition frequencies can again be calculated using the appropriate selection rules and Eqs. (2) and (3), and the results can be compared to experiment.

Initialize

 >

Input Parameters

We begin by specifying the atoms in our diatomic.  Note that the energy levels and spacings will depend on the isotopes we specify!  The variables elementA and elementB are the element symbols for the atoms in our diatomic, and MassNumberA and MassNumberB are their respective mass numbers (integers).    The variables Rmin and Rmax define the range over which the PES is computed and determine the level of accuracy of the vibrational calculation. The smaller the range, the quicker the calculation but the less accurate the results.  Start with a reasonable guess and proceed through each step to calculate the fundamental frequency.  Vary Rmin and Rmax until some acceptable level of convergence has been reached.  Once this has been done, Rmin and Rmax remain fixed for the different electronic structure methods and basis sets.

Finally, we specify the electronic structure method and basis set here.  Begin with ESmethod = HartreeFock and AObasis = "sto-3g" and proceed.  Once the fundamental has been calculated, repeat each step with a new basis and/or method.

 > $\mathrm{elementA}≔H;$
 ${\mathrm{elementA}}{≔}{H}$ (3.1)
 > $\mathrm{MassNumberA}≔1;$
 ${\mathrm{MassNumberA}}{≔}{1}$ (3.2)
 > $\mathrm{elementB}≔\mathrm{Cl};$
 ${\mathrm{elementB}}{≔}{\mathrm{Cl}}$ (3.3)
 >
 ${\mathrm{MassNumberB}}{≔}{35}$ (3.4)
 >
 ${\mathrm{wavenum_observed}}{≔}{2885.80000000}$ (3.5)
 >
 ${\mathrm{Rmin}}{≔}{0.80000000}$ (3.6)
 >
 ${\mathrm{Rmax}}{≔}{3.00000000}$ (3.7)
 >
 ${\mathrm{ESmethod}}{≔}{\mathrm{DensityFunctional}}$ (3.8)
 > 
 ${\mathrm{ESbasis}}{≔}{"6-31g"}$ (3.9)

Normal Mode Analysis

In this section, calculate the fundamental mode of the diatomic using a normal mode analysis, which approximates the vibration of the diatomic as a harmonic oscillator.  First define the molecule with a bond length of 1 Angstrom and then use the GeometryOptimization function to redefine the molecule with the optimum bond length for the desired method and basis. Then use the VibrationalModes function to determine the energy (in wavenumbers) of the fundamental.  Calculate the percent error.

 > $\mathrm{molec}≔\left[\left[\mathrm{convert}\left(\mathrm{elementA},\mathrm{string}\right),0,0,0\right],\left[\mathrm{convert}\left(\mathrm{elementB},\mathrm{string}\right),0,0,1\right]\right];$
 ${\mathrm{molec}}{≔}\left[\left[{"H"}{,}{0}{,}{0}{,}{0}\right]{,}\left[{"Cl"}{,}{0}{,}{0}{,}{1}\right]\right]$ (4.1)
 > $\mathrm{molec},\mathrm{data}≔\mathrm{GeometryOptimization}\left(\mathrm{molec},\mathrm{ESmethod},\mathrm{basis}=\mathrm{ESbasis}\right);$
 ${\mathrm{molec}}{,}{\mathrm{data}}{≔}\left[\left[{"H"}{,}{0}{,}{0}{,}{0}\right]{,}\left[{"Cl"}{,}{3.31526629}{}{{10}}^{{-10}}{,}{-9.45757021}{}{{10}}^{{-10}}{,}{1.32081084}\right]\right]{,}{table}{}\left(\left[{\mathrm{mo_occ}}{=}\begin{array}{c}\left[\begin{array}{c}{2.00000000}\\ {2.00000000}\\ {2.00000000}\\ {2.00000000}\\ {2.00000000}\\ {2.00000000}\\ {⋮}\end{array}\right]\\ \hfill {\text{15 element Vector[column]}}\end{array}{,}{\mathrm{mo_symmetry}}{=}\begin{array}{c}\left[\begin{array}{c}{"A"}\\ {"A"}\\ {"A"}\\ {"A"}\\ {"A"}\\ {"A"}\\ {⋮}\end{array}\right]\\ \hfill {\text{15 element Vector[column]}}\end{array}{,}{\mathrm{converged}}{=}{1}{,}{\mathrm{dipole}}{=}\left[\begin{array}{c}{-4.45052230}{}{{10}}^{{-10}}\\ {1.26961595}{}{{10}}^{{-9}}\\ {-1.77310614}\end{array}\right]{,}{\mathrm{charges}}{=}\left[\begin{array}{c}{0.17963699}\\ {-0.17963699}\end{array}\right]{,}{\mathrm{aolabels}}{=}\begin{array}{c}\left[\begin{array}{c}{"0 H 1s"}\\ {"0 H 2s"}\\ {"1 Cl 1s"}\\ {"1 Cl 2s"}\\ {"1 Cl 3s"}\\ {"1 Cl 4s"}\\ {⋮}\end{array}\right]\\ \hfill {\text{15 element Vector[column]}}\end{array}{,}{\mathrm{mo_coeff}}{=}\begin{array}{c}\left[\begin{array}{ccccccc}{0.00025025}& {-0.00140120}& {0.00198866}& {0.}& {0.}& {0.16112321}& {\dots }\\ {-0.00087918}& {0.00517389}& {-0.00488115}& {0.}& {0.}& {0.05078387}& {\dots }\\ {0.99649842}& {-0.28577239}& {-0.00373959}& {0.}& {0.}& {0.08022513}& {\dots }\\ {0.01273518}& {1.02746865}& {0.01354787}& {0.}& {0.}& {-0.35999984}& {\dots }\\ {-0.00431431}& {0.03797009}& {-0.00168279}& {0.}& {0.}& {0.71942061}& {\dots }\\ {0.00266304}& {-0.01652813}& {0.00348467}& {0.}& {0.}& {0.26688961}& {\dots }\\ {⋮}& {⋮}& {⋮}& {⋮}& {⋮}& {⋮}& {}\end{array}\right]\\ \hfill {\text{15 Ã— 15 Matrix}}\end{array}{,}{\mathrm{group}}{=}{"C1"}{,}{\mathrm{mo_energy}}{=}\begin{array}{c}\left[\begin{array}{c}{-101.53709945}\\ {-9.47156815}\\ {-7.23564276}\\ {-7.22348732}\\ {-7.22348732}\\ {-0.85207881}\\ {⋮}\end{array}\right]\\ \hfill {\text{15 element Vector[column]}}\end{array}{,}{\mathrm{rdm1}}{=}\begin{array}{c}\left[\begin{array}{ccccccc}{0.22589453}& {0.17552213}& {0.00971109}& {-0.03851761}& {0.07259390}& {-0.06596280}& {\dots }\\ {0.17552213}& {0.15093664}& {-0.01247091}& {0.04740469}& {-0.07212500}& {-0.11222339}& {\dots }\\ {0.00971109}& {-0.01247091}& {2.16399561}& {-0.62777132}& {0.10108317}& {0.07277930}& {\dots }\\ {-0.03851761}& {0.04740469}& {-0.62777132}& {2.40835055}& {-0.51357131}& {-0.29614731}& {\dots }\\ {0.07259390}& {-0.07212500}& {0.10108317}& {-0.51357131}& {1.18360618}& {0.52178177}& {\dots }\\ {-0.06596280}& {-0.11222339}& {0.07277930}& {-0.29614731}& {0.52178177}& {0.27590606}& {\dots }\\ {⋮}& {⋮}& {⋮}& {⋮}& {⋮}& {⋮}& {}\end{array}\right]\\ \hfill {\text{15 Ã— 15 Matrix}}\end{array}{,}{\mathrm{populations}}{=}\begin{array}{c}\left[\begin{array}{c}{0.49412706}\\ {0.32623595}\\ {1.99971072}\\ {1.99676793}\\ {1.42627968}\\ {0.54423925}\\ {⋮}\end{array}\right]\\ \hfill {\text{15 element Vector[column]}}\end{array}{,}{\mathrm{e_tot}}{=}{-460.70734509}\right]\right)$ (4.2)
 > $\mathrm{emodes},\mathrm{vmodes}≔\mathrm{VibrationalModes}\left(\mathrm{molec},\mathrm{ESmethod},\mathrm{basis}=\mathrm{ESbasis}\right);$
 ${\mathrm{emodes}}{,}{\mathrm{vmodes}}{≔}\left[\begin{array}{c}{2713.23473969}\end{array}\right]{,}\left[\begin{array}{c}{-8.97667711}{}{{10}}^{{-7}}\\ {2.30719256}{}{{10}}^{{-6}}\\ {-0.98608093}\\ {3.91583629}{}{{10}}^{{-8}}\\ {-2.80189220}{}{{10}}^{{-7}}\\ {0.16626606}\end{array}\right]$ (4.3)
 > $\mathrm{PercentError}≔\frac{\left(\mathrm{emodes}\left[1\right]-\mathrm{wavenum_observed}\right)}{\mathrm{wavenum_observed}}\cdot 100;$
 ${\mathrm{PercentError}}{≔}{-5.97980665}$ (4.4)

Now, consider the rovibrational transitions allowed for a heteronuclear diatomic AB.   The vibrational selection rules are Δv = Δ±1,  Δl = ±1. Therefore, for transitions between v = 0 --> v = 1, there is a manifold of Δl = ±1 lines.

 > $\mathrm{mA}≔\mathrm{evalf}\left(\mathrm{Element}\left(\mathrm{elementA}\left[\mathrm{MassNumberA}\right],\mathrm{atomicmass},\mathrm{system}=\mathrm{SI}\right)\right);\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\mathrm{mB}≔\mathrm{evalf}\left(\mathrm{Element}\left(\mathrm{elementB}\left[\mathrm{MassNumberB}\right],\mathrm{atomicmass},\mathrm{system}=\mathrm{SI}\right)\right);$
 ${\mathrm{mA}}{≔}{1.67353398}{}{{10}}^{{-27}}$
 ${\mathrm{mB}}{≔}{5.80671857}{}{{10}}^{{-26}}$ (4.5)
 > $\mathrm{lmax}≔5;$
 ${\mathrm{lmax}}{≔}{5}$ (4.6)
 >
 ${\mathrm{R0}}{≔}{1.32081083}{}{{10}}^{{-10}}$ (4.7)
 > $\mathrm{reduced_mass}≔\frac{\mathrm{mA}\cdot \mathrm{mB}}{\mathrm{mA}+\mathrm{mB}};$
 ${\mathrm{reduced_mass}}{≔}{1.62665279}{}{{10}}^{{-27}}$ (4.8)
 > $\mathrm{moment_inertia}≔\mathrm{reduced_mass}\cdot {\mathrm{R0}}^{2};$
 ${\mathrm{moment_inertia}}{≔}{2.83776291}{}{{10}}^{{-47}}$ (4.9)
 >
 ${{\mathrm{Erot}}}_{{0}}{≔}{0.}$
 ${{\mathrm{Erot}}}_{{1}}{≔}{19.72874569}$
 ${{\mathrm{Erot}}}_{{2}}{≔}{59.18623706}$
 ${{\mathrm{Erot}}}_{{3}}{≔}{118.37247413}$
 ${{\mathrm{Erot}}}_{{4}}{≔}{197.28745688}$
 ${{\mathrm{Erot}}}_{{5}}{≔}{295.93118532}$ (4.10)
 > 
 > $#\mathrm{Calculate P- and R-branches for l = 0, 1, 2, ..., lmax}$

P-branch: Dl = -1

 >
 ${P}{≔}{2693.50599401}$
 ${P}{≔}{2673.77724832}$
 ${P}{≔}{2654.04850263}$
 ${P}{≔}{2634.31975694}$
 ${P}{≔}{2614.59101125}$ (4.11)

R-branch: Dl = +1

 >
 ${R}{≔}{2732.96348538}$
 ${R}{≔}{2752.69223107}$
 ${R}{≔}{2772.42097676}$
 ${R}{≔}{2792.14972244}$
 ${R}{≔}{2811.87846813}$ (4.12)

What do you notice about the spacing between lines in the P-branch?  in the R-branch? Are these observations consistent with an experimental IR spectrum of the diatomic AB?

 Answer There is equal spacing between lines in both P and R branches. This is consistent with a rigid rotor approximation!

Ab initio Calculation: Density Functional Theory

Construct Potential Energy Surface

In this section, use DFT to solve for the potential energy surface of diatomic AB.  Begin by using the ScientificConstants package to get the mass of each of the elements A and B. Use the QuantumChemistry package to calculate the PES for a range of separation values between Rmin and Rmax.     (Note, DFT may not converge for at larger values of R. This should not affect the final results.)

 > $\mathrm{mA}≔\mathrm{evalf}\left(\mathrm{Element}\left(\mathrm{elementA}\left[\mathrm{MassNumberA}\right],\mathrm{atomicmass},\mathrm{system}=\mathrm{Atomic}\right)\right);\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\mathrm{mB}≔\mathrm{evalf}\left(\mathrm{Element}\left(\mathrm{elementB}\left[\mathrm{MassNumberB}\right],\mathrm{atomicmass},\mathrm{system}=\mathrm{Atomic}\right)\right);$
 ${\mathrm{mA}}{≔}{1837.15393006}$
 ${\mathrm{mB}}{≔}{63744.36349535}$ (5.1.1)
 > $#\mathrm{mA≔evalf}\left(\mathrm{Element}\left(\mathrm{elementA},\mathrm{atomicweight},\mathrm{system}=\mathrm{Atomic}\right)\right);\mathrm{mB}≔\mathrm{evalf}\left(\mathrm{Element}\left(\mathrm{elementB},\mathrm{atomicweight},\mathrm{system}=\mathrm{Atomic}\right)\right);$
 >

Calculate Vibrational Energy Levels

Once you have a PES, you can solve the vibrational Schrodinger equation.  Calculate the energies for a range of angular momentum values and to calculate the P- and R-branches.

 > $n≔\mathrm{trunc}\left(20\cdot \left(\mathrm{Rmax}-\mathrm{Rmin}\right)\right):\mathrm{lmax}≔5:$
 >
 >
 >

Results

Fundamental Transition Frequency

We can now calculate the fundamental frequency ( in wavenumbers, 1/cm). This number can be compared to the experimental value for the AB diatomic of interest.  Enter the results into Table 1 below!

 > $\mathrm{wavenum}≔\left(\mathrm{Evl}\left[2,1\right]-\mathrm{Evl}\left[1,1\right]\right)\cdot 219476.329478$
 ${\mathrm{wavenum}}{≔}{2621.54229484}$ (6.1.1)
 > $\mathrm{PercentError}≔\frac{\left(\mathrm{wavenum}-\mathrm{wavenum_observed}\right)}{\mathrm{wavenum_observed}}\cdot 100$
 ${\mathrm{PercentError}}{≔}{-9.15717323}$ (6.1.2)

Table 1: Fundamental  Frequencies for Diatomic AB   (v0 = 2872.3 cm-1 for HCl.)

 ESMethod AObasis Fundamental Frequency v0 (cm-1) %Error Density Functional Theory 6-31G

P-Branch and R-Branch

As a final exercise, we consider the rovibrational transitions allowed for a heteronuclear diatomic AB using the method/basis that gave the best results in Table 1.   If necessary, recompute the rovibrational eigenvalues for this combination before proceeding.

P-branch: DJ = -1

 >
 ${P}{≔}{2552.23818551}$
 ${P}{≔}{2480.89533062}$
 ${P}{≔}{2407.67319968}$
 ${P}{≔}{2332.73352161}$
 ${P}{≔}{2256.23985068}$ (6.2.1)

R-branch: DJ = +1

 >
 ${R}{≔}{2688.65090282}$
 ${R}{≔}{2753.41043599}$
 ${R}{≔}{2815.67098378}$
 ${R}{≔}{2875.28678878}$
 ${R}{≔}{2932.11674112}$ (6.2.2)

What do you notice about the spacing between lines in the P-branch?  in the R-branch? Are these observations consistent with an experimental IR spectrum of the diatomic AB?

 Answer P-branch: spacing increases as J increases. R-branch: spacing decreases as J increases.

Appendix

 References J. C.Light, I. P.Hamilton, and J. V. Vill, J. Chem. Phys. 82 ,1400 (1985).