StatisticalThermodynamics - Maple Help
For the best experience, we recommend viewing online help using Google Chrome or Mozilla Firefox.

Online Help

All Products    Maple    MapleSim


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

Statistical Thermodynamics: Calculating Thermodynamic Functions DH, DG, and DS

Copyright (c) RDMCHEM LLC 2019

 

 

Overview

Initialize

Specify Molecule

Calculate Thermodynamics

Calculate ΔH, ΔG, and ΔS (after completing above for all reactants and products)

Appendix

Overview

 

In this activity, you will apply statistical thermodynamics to calculate the enthalpy, free energy, and entropy of reaction for the combustion of ethene to form ethylene oxide at 298 K:

 

 

Rxn (1)

 

To simplify the calculations, you will treat each species as an ideal gas. The thermal entropy (S) (J/mol K) and the internal energy (E) (J/mol) are given by:

 

(1)

 

and

(2)

 

where q is the molecular partition function qe=qt qe qr qv , where the subscript t, e, r, and v refer to the translational, electronic, rotational, and vibrational contributions, respectively.  To calculate q, we further assume each species is in the ground electronic state (with degeneracy 2 spin + 1) with rovibrational energies given by harmonic oscillator and rigid rotor approximations. These assumptions lead to the following expressions for q, S, and E:

 

Translational:

(3)

(4)

(5)

where P = 1 atm = 100.325 kPa.

 

Electronic:

(6)

(7)

 

The thermal contribution to internal energy is 0 since the qe does not depend on temperature.

 

Rotational:

Linear Molecule:

 

(8)

(9)

(10)

where σr is the rotational symmetry number and Θr is the rotational temperature.  Each molecule we will consider in this activity is characterized by σr = 2.  

Nonlinear Polyatomic Molecule:

 

(11)

(12)

(13)

 

where ΘA, ΘB, ΘC are the rotational temperatures corresponding to the principal moments of inertia.

 

Vibrational:

(14)

(15)

(16)

 

where s = 3N - 5 normal modes for a linear molecule and s = 3N - 6 for a nonlinear polyatomic.

 

Once the translational, electronic, rotational, and vibrational thermal contributions to the total entropy and internal energy have been calculated, total entropy (S), internal energy (E), enthalpy (H), and free-energy (G) can be calculated as follows:

 

(17)

(18)

(19)

(20)

 

Note we have corrected for the vibrational zero-point-energy for the internal energy.  We will calculate S, H, and G for each of the reactants and product in Rxn (1) and calculate ΔH, ΔG, and ΔS as

 

ΔHrxn = ∑nHprod−∑nHreact

(21)

etc.

 

Initialize

 

We begin with initializing the QuantumChemistry package, along with the ScientificConstants and LinearAlgebra packages.

 

> 

restart:withQuantumChemistry:withScientificConstants:withLinearAlgebra:Digits≔15:

 

Specify Molecule

 

Here we specify the coordinates of the molecule for which we want to calculate thermodynamics as well as some related calculate parameters, such as symmetry number, charge, AObasis, etc.  Note that the geometries provided in .xyz data files are at the HartreeFock / STO-3G level.   If you wish to use a different method/basis combination, begin with the provided .xyz coordinates for each molecule and uncomment the geometry optimization line below.

 

> 

#molec≔ReadXYZethylene.xyz;

> 

molec≔ReadXYZO2.xyz;

molec≔O,0.,0.,0.,O,0.,0.,1.22169673

(3.1)
> 

# molec≔ReadXYZethyleneoxide.xyz;

> 

molec_symmetry_number≔2 :    # symmetry number 

> 

molec_spin≔0:         # molec_spin = 0 for singlet, molec_spin = 1 for triplet

> 

molec_charge≔0:

> 

freqmethod≔'HartreeFock':

> 

freqbasis≔sto-3g:

> 

energymethod≔'HartreeFock':

> 

energybasis≔sto-3g:

> 

freq_scaling≔0.8953:   # see Ref. 2 for scaling factors (HF = 0.8953, DFT = 0.9613

> 

temp≔298.15:

> 

molec_index≔2:     # identifies which reactant or product. Doesn't matter which is 1, 2, or 3, just keep track!

> 

Calculate Thermodynamics

 

If  not using Hartree-Fock and sto-3g, uncomment the following line:

 

> 

#molec,molec_data≔GeometryOptimizationmolec,energymethod,basis=energybasis,spin=2⋅ molec_spin,charge=molec_charge;

 

Calculate electronic energy in J/mol using AOmethod, AObasis, spin, and charge defined above: 

 

> 

molec_energy≔Energymolec,method=energymethod,basis=energybasis,spin=2⋅molec_spin,charge=molec_charge; au2Jmol≔4.359744650ⅇ−18⋅evalfConstantAvogadro_constant: molec_energy≔molec_energy⋅au2Jmol

molec_energy≔−147.55157168

molec_energy≔−3.87396598⁢108

(4.1)

 

 

To calculate the vibrational energies, perform a vibrational analysis:

 

> 

emodes,vmodes ≔ VibrationalModesmolec,freqmethod,spin=2⋅molec_spin,charge=molec_charge, basis=freqbasis; emodes≔emodes⋅freq_scaling;

emodes,vmodes≔2057.72134813,5.41846411⁢10−85.41846015⁢10−8−0.707106931.26430759⁢10−7−5.41846392⁢10−80.70710663

emodes≔1842.27792298

(4.2)

 

 

Calculate principal moments of inertia:

 

> 

num_atoms≔nopsmolec:mass_tot≔0.0:for i from 1 by 1 to num_atoms do      atomi≔parsemoleci1;      xposi≔moleci2⋅1e−10;      yposi≔moleci3⋅1e−10;      zposi≔moleci4⋅1e−10;     massi≔evalfElementatomi,atomicweight,system=SI;     mass_tot≔mass_tot+massi;end do:xpos_cm≔0.0:ypos_cm≔0.0:zpos_cm≔0.0:for i from 1 by 1 to num_atoms do     xpos_cm≔xpos_cm+massimass_tot⋅xposi;       ypos_cm≔ypos_cm+massimass_tot⋅yposi;      zpos_cm≔zpos_cm+massimass_tot⋅zposi;end do:for i from 1 by 1 to num_atoms do              #places origin at center of mass    xposi≔xposi−xpos_cm;    yposi≔yposi−ypos_cm;    zposi≔zposi−zpos_cm;end do:moment_inertia≔Matrix3,3,0:for i  from 1 by 1 to num_atoms do    moment_inertia1,1≔moment_inertia1,1+massi⋅yposi2+zposi2;   moment_inertia1,2≔moment_inertia1,2−massi⋅xposi⋅yposi;  moment_inertia1,3≔moment_inertia1,3−massi⋅xposi⋅zposi;  moment_inertia2,2≔moment_inertia2,2+massi⋅xposi2+zposi2;  moment_inertia2,3≔moment_inertia2,3−massi⋅yposi⋅zposi;  moment_inertia3,3≔moment_inertia3,3+massi⋅xposi2+yposi2;end do:moment_inertia2,1=moment_inertia1,2: moment_inertia3,1=moment_inertia1,3:moment_inertia3,2=moment_inertia2,3:diag≔sortReEigenvaluesmoment_inertia:I_A≔diag1;I_B≔diag2;I_C≔diag3;

I_A≔0.

I_B≔1.98266764⁢10−46

I_C≔1.98266764⁢10−46

(4.3)

 

 

Now we are ready to calculate partition functions!

 

Translational partition function:

 

> 

q_t≔evalf2⋅Pi⋅mass_tot⋅Constantk⋅tempConstanth232⋅Constantk⋅temp101325.0 ;

q_t≔7.11466749⁢106

(4.4)

 

Electronic partition function:

 

> 

q_e≔2⋅molec_spin+1; 

q_e≔1

(4.5)

 

 

Rotational temperatures and partition function:

 

> 

theta_A≔evalfConstanth28⋅π2⋅ I_A⋅Constantk&semi;   &num; Kelvintheta_B≔evalfConstanth28⋅π2⋅ I_B⋅Constantk&semi;     &num; Kelvintheta_C≔evalfConstanth28⋅π2⋅ I_C⋅Constantk &semi;     &num; Kelvin if I_A<1e−48 then           I_A≔0&semi;           q_r≔1molec_symmetry_number⋅temptheta_B⋅theta_C12&semi;            num_modes≔3⋅num_atoms−5&semi; else           q_r≔Pi12molec_symmetry_number⋅temp32theta_A⋅theta_B⋅theta_C12&semi;           num_modes≔3⋅num_atoms−6&semi; end if&semi;

theta_A≔Float⁡∞

theta_B≔2.03137108

theta_C≔2.03137108

I_A≔0

q_r≔73.38639462

num_modes≔1

(4.6)

 

Vibrational temperatures and partition function:

 

> 

 q_v≔1&colon; for i from 1 by 1 to num_modes do             theta_vi≔evalfConstanth⋅Reemodesi⋅Constantc⋅100Constantk&semi;              q_v≔q_v⋅exp−1.0⋅theta_vi2⋅temp1−exp−1.0⋅theta_vitemp&semi;  end do&colon;q_v≔q_v&semi;

q_v≔0.01173726

(4.7)

 

ZPE correction in J/mol:

 

> 

 zpe_correction≔0.0&colon;for i from 1 by 1 to num_modes do zpe_correction≔zpe_correction&plus;Reemodesi2⋅evalfConstantc⋅100⋅evalfConstanth⋅evalfConstantAvogadro_constant end do&colon; zpe_correction≔zpe_correction&semi; 

zpe_correction≔11019.26903875

(4.8)

 

 

Now we can calculate total S, E, H, and G according to the equations provided in the Overview:

 

> 

 &num; Translational ContributionsS_t≔evalfConstantR⋅lnq_t&plus;1&plus;32&semi;E_t≔evalf32⋅ConstantR⋅temp&semi;  &num; Electronic Contributions S_e≔evalfConstantR⋅lnq_e&semi;  &num; Rotational Contributions if I_A &equals; 0 then         S_r≔evalfConstantR⋅lnq_r&plus;1 &semi;        E_r≔evalfConstantR⋅temp&semi; else        S_r≔evalfConstantR⋅lnq_r&plus;32&semi;       E_r≔evalf32⋅ ConstantR⋅temp&semi; end if&colon;  &num; Vibrational Contributions S_v≔0&colon;E_v≔0&colon;for i from 1 by 1 to num_modes do     S_v≔evalfS_v&plus;ConstantR⋅theta_vitempexptheta_vitemp−1−ln1−exp−theta_vitemp&semi;    E_v≔evalfE_v&plus;ConstantR⋅theta_vi⋅12&plus;1exptheta_vitemp−1&semi; end do&colon;S_v≔S_v&semi;E_v≔E_v&semi;  &num;Total Entropy  S0molec_index≔S_e&plus;S_t&plus;S_v&plus;S_r&semi;  &num; Total Energy  E0molec_index≔molec_energy&plus;E_t&plus;E_v&plus;E_r&plus;zpe_correction&semi;   &num; Enthalpy  H0molec_index≔evalfE0molec_index&plus;ConstantR⋅temp&semi;  &num; Gibbs Free Energy G0molec_index≔evalfH0molec_index−temp⋅S0molec_index&semi;         

S_t≔151.96894459

E_t≔3718.43428406

S_e≔0.

S_v≔0.01132692

E_v≔11027.51974508

S02≔196.01147687

E02≔−3.87368354⁢108

H02≔−3.87365875⁢108

G02≔−3.87424316⁢108

(4.9)

 

Repeat the above calculations for each reactant and product in Rxn (1). Don't forget to change the molec_index!

> 

Calculate DH, DG, and DS (after completing above for all reactants and products)

 

 

> 

&Delta;S≔ S03− S01−12⋅ S02&semi;

&Delta;S≔S03−S01−98.00573844

(5.1)
> 

&Delta;G≔G03− G01−12⋅ G021000.0&semi;

&Delta;G≔0.00100000⁢G03−0.00100000⁢G01+193712.15788334

(5.2)
> 

&Delta;H≔ H03− H01−12⋅ H021000.0&semi;

&Delta;H≔0.00100000⁢H03−0.00100000⁢H01+193682.93747242

(5.3)
> 

 

 

NIST lists the following values for enthalpy of formation and free-energy of formation:

 

 

&Delta;Hf (kJ/mol)

&Delta;Gf (kJ/mol) (calc.)

S (1 bar) (J/molK)

ethylene

52.4

61.4

219.3

ethylene oxide

-52.6

-51.7

243.0

oxygen

0

0

205.15

   

Using Hess's Law, calculate ΔH, ΔG, and ΔS for Rxn (1) and compare to your calculated results.

 

 

> 

&Delta;SHess≔243.0 − 219.3−12⋅205.15&semi;

&Delta;SHess≔−78.87500000

(5.4)
> 

&Delta;HHess≔−52.6−52.4&semi;

&Delta;HHess≔−105.00000000

(5.5)
> 

&Delta;GHess≔−51.7−61.4&semi;

&Delta;GHess≔−113.10000000

(5.6)
> 

S02&semi;

196.01147687

(5.7)
> 

 

Appendix

References

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

2) JPCA scaling factors 2007

3) https://webbook.nist.gov/cgi/inchi/InChI=1S/C2H4/c1-2/h1-2H2

4) https://webbook.nist.gov/cgi/inchi/InChI=1S/C2H4O/c1-2-3-1/h1-2H2