HarmonicOscillator - Maple Help
For the best experience, we recommend viewing online help using Google Chrome or Microsoft Edge.

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;

Digits≔15

(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;

A≔H

B≔F

(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;

dataA≔table⁡names=hydrogen,symbol=H,meltingpoint=13.81000000⁢K,atomicnumber=1,electronegativity=2.10000000,ionizationenergy=13.59840000⁢eV,atomicweight=1.00794000⁢amu,name=hydrogen,boilingpoint=20.28000000⁢K

(2.2.1.2)
> 

dataB ≔ AtomicDataB;

dataB≔table⁡names=fluorine,symbol=F,meltingpoint=53.53000000⁢K,atomicnumber=9,electronegativity=3.98000000,ionizationenergy=17.42280000⁢eV,atomicweight=18.99840320⁢amu,name=fluorine,electronaffinity=3.40119000⁢eV,boilingpoint=85.03000000⁢K

(2.2.1.3)

Set the masses

> 

mA ≔ dataAatomicweight;

mA≔1.00794000⁢amu

(2.2.1.4)
> 

mB ≔ dataBatomicweight;

mB≔18.99840320⁢amu

(2.2.1.5)

Compute the reduced mass in amu

> 

mu0 ≔ mA⋅mBmA+mB;

μ0≔0.95715895⁢amu

(2.2.1.6)

Convert the reduced mass from amu to kg

> 

mu ≔ convertmu0,units,'kg';

μ≔1.58940092⁢10−27⁢kg

(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;

bonds≔0.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;

molecules≔H,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;

pes≔−2.84449898⁢R5+17.02510742⁢R4−41.28045339⁢R3+50.73972578⁢R2−31.33797062⁢R−92.31177987

(2.2.2.4)

The potential energy surface (curve) can be plotted

> 

plotpesR, R=bonds1..bonds−1, 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;

eq≔−14.22249490⁢R4+68.10042970⁢R3−123.84136016⁢R2+101.47945156⁢R−31.33797062=0

(2.2.2.5)

Solving the resulting equation yields the equilibrium bond length

> 

R_eq ≔ fsolveeq, R=bonds1..bonds−1;

R_eq≔0.90161288

(2.2.2.6)

Differentiating the potential energy curve yields

> 

d2pes ≔ diffpes,R$2;

d2pes≔−56.88997962⁢R3+204.30128909⁢R2−247.68272033⁢R+101.47945156

(2.2.2.7)

We define the atomic units of the spring constant

> 

unit_factor ≔ UnitsUnit'hartree'UnitsUnit'angstrom'2;

unit_factor≔E0Å2

(2.2.2.8)

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

> 

k0 ≔ subsR=R_eq, d2pes⋅unit_factor;

k0≔2.54705679Å2⁢E0

(2.2.2.9)

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

> 

k ≔ convertk0,units,'Jm2';

k≔1110.45171985⁢Jm2

(2.2.2.10)

Angular Frequency

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

> 

omega0 ≔ sqrtkmu;

ω0≔8.35859168⁢1014⁢Jm2kg

(2.2.3.1)
> 

omega ≔ simplifyomega0;

ω≔8.35859168⁢1014⁢1s

(2.2.3.2)

Energy Spacing

We define Z from Maple's knowledge of scientific constants

> 

hbar ≔ evalfScientificConstantsConstant'hbar', units;

ℏ≔1.05457180⁢10−34⁢m2⁢kgs

(2.2.4.1)

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

> 

dE0 ≔ hbar⋅omega;

dE0≔8.81473508⁢10−20⁢m2⁢kgs⁢1s

(2.2.4.2)

The units can be combined with the simplify command

> 

dE ≔ simplifydE0;

dE≔8.81473508⁢10−20⁢J

(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..bonds−1, axes=boxed, labels='R','E', color=blue, thickness=3:p_ho≔ plotW0+k0⋅R−R_eq2, R=bonds1..bonds−1, 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).


Download Help Document