HarmonicOscillator - Maple Help

Online Help

All Products    Maple    MapleSim


Vibrational Motion and the Harmonic Oscillator

Copyright (c) RDMCHEM LLC 2021

 

 

Overview

Vibrational Motion

References

Overview

Vibrational motion of molecules can be modeled around the equilibrium geometry to a good approximation by a harmonic oscillator.  In this lesson we will model several diatomic molecules by harmonic oscillators.   The harmonic oscillator model consists of a particle in a parabolic potential:

Vx = 12kx2

where x is the displacement coordinate and k is a constant, often known as the spring constant.  From the mass m of the particle and the spring constant k the angular frequency ω of the oscillator can be computed from

    ω = km.

A diatomic molecule can be modeled by a harmonic oscillator through two steps.  (1) The diatomic molecule A-B with masses mA and mB of the atoms A and B separated by a distance R can be represented by a single particle with reduced mass μ oscillating about the equilibrium bond distance Req−the distance at which the minimum potential energy occurs.  Specifically, we can define the reduced mass μ as

μ = m__Am__Bm__A+m__B

and the displacement coordinate x as

x = R  R__eq .

(2) The potential energy of the diatomic molecule as a function of the displacement coordinate x can be approximated as a parabolic potential.  While the actual potential energy W(x) is not parabolic (see the blue curve in Fig. 1), it can be approximated in the vicinity of the equilibrium bond distance Req by a parabolic potential  (see the red curve in Fig.1)

Wx W__0+12kx2  

where the force constant k is the second derivative of the actual potential energy W(x) with respect to x (or by the chain rule R)

k = ⅆ2ⅆx2Wx=ⅆ2ⅆR2WR .

Figure 1: Exact (blue) and harmonic (red) potential energy curves of diatomic nitrogen

 

In this lesson for several diatomic molecules we will explore their approximations as harmonic oscillators by computing for each of them: a force constant k, a reduced mass μ, and the angular frequency ω.  From the angular frequency we can compute the energies of the oscillators

E__n=n+12 ℏω

and the energy differences

ΔE=ℏω .

The harmonic oscillator approximation can be extended to polyatomic molecules through coupled harmonic oscillators and their normal modes (which is discussed in a different lesson).

Vibrational Motion

We will model three diatomic molecules hydrogen fluoride HF, diatomic nitrogen N2, and carbon monoxide CO.  Computed properties will complete the following Table:

 

Table 1: Harmonic Oscillator Approximation of Diatomic Molecules

Molecule

Reduced Mass (μ)

Spring Constant (k)

Angular Frequency (ω)

Energy Spacing (ΔE)

HF

 

 

 

 

N2

 

 

 

 

CO

 

 

 

 

 

(a) Use the worksheet below to compute these properties for HF.

 

(b) Modify the worksheet below to complete the table for N2 and CO.

 

Quantum Chemistry

We set the number of Digits to be used in computations to 15 and load the Quantum Chemistry package using Maple's with command.

Digits  15;

Digits15

(2.1.1)

withQuantumChemistry;

AOLabels,ActiveSpaceCI,ActiveSpaceSCF,AtomicData,BondAngles,BondDistances,Charges,ChargesPlot,ContractedSchrodinger,CorrelationEnergy,CoupledCluster,DensityFunctional,DensityPlot3D,Dipole,DipolePlot,Energy,ExcitationEnergies,ExcitationSpectra,ExcitationSpectraPlot,ExcitedStateEnergies,ExcitedStateSpins,FullCI,GeometryOptimization,HartreeFock,Interactive,Isotopes,MOCoefficients,MODiagram,MOEnergies,MOIntegrals,MOOccupations,MOOccupationsPlot,MOSymmetries,MP2,MolecularData,MolecularDictionary,MolecularGeometry,NuclearEnergy,NuclearGradient,OscillatorStrengths,Parametric2RDM,PlotMolecule,Populations,Purify2RDM,RDM1,RDM2,RTM1,ReadXYZ,Restore,Save,SaveXYZ,SearchBasisSets,SearchFunctionals,SkeletalStructure,Thermodynamics,TransitionDipolePlot,TransitionDipoles,TransitionOrbitalPlot,TransitionOrbitals,Variational2RDM,VibrationalModeAnimation,VibrationalModes,Video

(2.1.2)

Diatomic Molecule A-B

Reduced Mass

Here we set the variables A and B to the strings of the atoms in the diatomic molecule (By changing these values, you can use the worksheet to treat other diatomic molecules!)

A  H;B  F;

AH

BF

(2.2.1.1)

We use the AtomicData command in the Quantum Chemistry package to obtain the masses as well as other data

dataA  AtomicDataA;

dataAtablenames=hydrogen,symbol=H,meltingpoint=13.81000000K,atomicnumber=1,electronegativity=2.10000000,ionizationenergy=13.59840000eV,atomicweight=1.00794000amu,name=hydrogen,boilingpoint=20.28000000K

(2.2.1.2)

dataB  AtomicDataB;

dataBtablenames=fluorine,symbol=F,meltingpoint=53.53000000K,atomicnumber=9,electronegativity=3.98000000,ionizationenergy=17.42280000eV,atomicweight=18.99840320amu,name=fluorine,electronaffinity=3.40119000eV,boilingpoint=85.03000000K

(2.2.1.3)

Set the masses

mA  dataAatomicweight;

mA1.00794000amu

(2.2.1.4)

mB  dataBatomicweight;

mB18.99840320amu

(2.2.1.5)

Compute the reduced mass in amu

mu0  mAmBmA+mB;

μ00.95715895amu

(2.2.1.6)

Convert the reduced mass from amu to kg

mu  convertmu0,units,'kg';

μ1.5894009210−27kg

(2.2.1.7)

Spring constant

To compute the equilibrium bond length and spring constant, we select a set of bond distances from the roots of the sixth-order Chebyshev polynomial that are suitable for interpolation

bonds  mapx  x/5+1.05, fsolveexpandChebyshevT6,x;

bonds0.85681483,0.90857864,0.99823619,1.10176381,1.19142136,1.24318517

(2.2.2.1)

We define a list of molecular geometries with each geometry corresponding to one of the bond distances

molecules mapR A,0,0,0,B,0,0,R, bonds;

moleculesH,0,0,0,F,0,0,0.85681483,H,0,0,0,F,0,0,0.90857864,H,0,0,0,F,0,0,0.99823619,H,0,0,0,F,0,0,1.10176381,H,0,0,0,F,0,0,1.19142136,H,0,0,0,F,0,0,1.24318517

(2.2.2.2)

The energies for each geometry may be then readily computed with the Energy command in the Quantum Chemistry package.

energies  mapEnergy,molecules,basis=cc-pVDZ;

energies−100.01686137,−100.01964383,−100.01018182,−99.98693789,−99.96201956,−99.94683591

(2.2.2.3)

We use polynomial interpolation to generate a polynomial in the bond distance R

pes  interpbonds,energies,R;

pes2.84449898R5+17.02510742R441.28045339R3+50.73972578R231.33797062R92.31177987

(2.2.2.4)

The potential energy surface (curve) can be plotted

plotpesR, R=bonds1..bonds1, axes=boxed, labels='R','E', color=blue, thickness=3;

Finally, we differential the potential energy curve with respect to R and set the derivative to zero.

eq  diffpes, R = 0;

eq14.22249490R4+68.10042970R3123.84136016R2+101.47945156R31.33797062=0

(2.2.2.5)

Solving the resulting equation yields the equilibrium bond length

R_eq  fsolveeq, R=bonds1..bonds1;

R_eq0.90161288

(2.2.2.6)

Differentiating the potential energy curve yields

d2pes  diffpes,R$2;

d2pes56.88997962R3+204.30128909R2247.68272033R+101.47945156

(2.2.2.7)

We define the atomic units of the spring constant

unit_factor  UnitsUnit'hartree'UnitsUnit'angstrom'2;

unit_factorE0Å2

(2.2.2.8)

Evaluating the second derivative at the equilibrium bond distance gives the spring constant

k0  subsR=R_eq, d2pesunit_factor;

k02.54705679Å2E0

(2.2.2.9)

We convert the atomic units to standard international (SI) units

k  convertk0,units,'Jm2';

k1110.45171985Jm2

(2.2.2.10)

Angular Frequency

We compute the angular frequency ω in terms of k and μ and then simplify the units

omega0  sqrtkmu;

ω08.358591681014Jm2kg

(2.2.3.1)

omega  simplifyomega0;

ω8.3585916810141s

(2.2.3.2)

Energy Spacing

We define Z from Maple's knowledge of scientific constants

hbar  evalfScientificConstantsConstant'hbar', units;

ℏ1.0545718010−34m2kgs

(2.2.4.1)

Using Z and ω, we compute the ΔE between energy levels (which is a constant)

dE0  hbaromega;

dE08.8147350810−20m2kgs1s

(2.2.4.2)

The units can be combined with the simplify command

dE  simplifydE0;

dE8.8147350810−20J

(2.2.4.3)

Comparison

The computed potential energy curve (blue) can be plotted against the harmonic potential energy curve (red) to assess the approximation

k0  convertk0,unit_free:W0  subsR=R_eq,pes:p_ex plotpes, R=bonds1..bonds1, axes=boxed, labels='R','E', color=blue, thickness=3:p_ho plotW0+k0RR_eq2, R=bonds1..bonds1, axes=boxed, labels='R','E', color=red, thickness=3:plotsdisplayp_ex,p_ho;

References

1. 

D. J. Griffiths and D. F. Schroeter, Introduction to Quantum Mechanics 3rd Edition (Cambridge University Press, 2018).

2. 

I. N. Levine, Quantum Chemistry 7th Edition (Pearson, New York, 2017).

3. 

J. J. Sakurai and J. Napolitano, Modern Quantum Mechanics 2nd Edition (Cambridge University Press, Cambridge, 2017).

4. 

J. P. Lowe, Quantum Chemistry Illustrated Edition (Academic Press, New York, 2012).

5. 

P. W. Atkins and R. S. Friedman, Molecular Quantum Mechanics 5th Edition (Oxford University Press, Oxford, 2010).

6. 

D. A. McQuarrie, Quantum Chemistry 2nd Edition (University Science, New York, 2007).

7. 

D. A. McQuarrie and J. D. Simon, Physical Chemistry: A Molecular Approach (University Science, New York, 1997).