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
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.
We begin with initializing the QuantumChemistry package, along with the ScientificConstants and LinearAlgebra packages.
restart:withQuantumChemistry:withScientificConstants:withLinearAlgebra:Digits≔15:
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
# 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!
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.87396598108
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.4184641110−85.4184601510−8−0.707106931.2643075910−7−5.4184639210−80.70710663
emodes≔1842.27792298
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.9826676410−46
I_C≔1.9826676410−46
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.11466749106
Electronic partition function:
q_e≔2⋅molec_spin+1;
q_e≔1
Rotational temperatures and partition function:
theta_A≔evalfConstanth28⋅π2⋅ I_A⋅Constantk; # Kelvintheta_B≔evalfConstanth28⋅π2⋅ I_B⋅Constantk; # Kelvintheta_C≔evalfConstanth28⋅π2⋅ I_C⋅Constantk ; # Kelvin if I_A<1e−48 then I_A≔0; q_r≔1molec_symmetry_number⋅temptheta_B⋅theta_C12; num_modes≔3⋅num_atoms−5; else q_r≔Pi12molec_symmetry_number⋅temp32theta_A⋅theta_B⋅theta_C12; num_modes≔3⋅num_atoms−6; end if;
theta_A≔Float∞
theta_B≔2.03137108
theta_C≔2.03137108
I_A≔0
q_r≔73.38639462
num_modes≔1
Vibrational temperatures and partition function:
q_v≔1: for i from 1 by 1 to num_modes do theta_vi≔evalfConstanth⋅Reemodesi⋅Constantc⋅100Constantk; q_v≔q_v⋅exp−1.0⋅theta_vi2⋅temp1−exp−1.0⋅theta_vitemp; end do:q_v≔q_v;
q_v≔0.01173726
ZPE correction in J/mol:
zpe_correction≔0.0:for i from 1 by 1 to num_modes do zpe_correction≔zpe_correction+Reemodesi2⋅evalfConstantc⋅100⋅evalfConstanth⋅evalfConstantAvogadro_constant end do: zpe_correction≔zpe_correction;
zpe_correction≔11019.26903875
Now we can calculate total S, E, H, and G according to the equations provided in the Overview:
# Translational ContributionsS_t≔evalfConstantR⋅lnq_t+1+32;E_t≔evalf32⋅ConstantR⋅temp; # Electronic Contributions S_e≔evalfConstantR⋅lnq_e; # Rotational Contributions if I_A = 0 then S_r≔evalfConstantR⋅lnq_r+1 ; E_r≔evalfConstantR⋅temp; else S_r≔evalfConstantR⋅lnq_r+32; E_r≔evalf32⋅ ConstantR⋅temp; end if: # Vibrational Contributions S_v≔0:E_v≔0:for i from 1 by 1 to num_modes do S_v≔evalfS_v+ConstantR⋅theta_vitempexptheta_vitemp−1−ln1−exp−theta_vitemp; E_v≔evalfE_v+ConstantR⋅theta_vi⋅12+1exptheta_vitemp−1; end do:S_v≔S_v;E_v≔E_v; #Total Entropy S0molec_index≔S_e+S_t+S_v+S_r; # Total Energy E0molec_index≔molec_energy+E_t+E_v+E_r+zpe_correction; # Enthalpy H0molec_index≔evalfE0molec_index+ConstantR⋅temp; # Gibbs Free Energy G0molec_index≔evalfH0molec_index−temp⋅S0molec_index;
S_t≔151.96894459
E_t≔3718.43428406
S_e≔0.
S_v≔0.01132692
E_v≔11027.51974508
S02≔196.01147687
E02≔−3.87368354108
H02≔−3.87365875108
G02≔−3.87424316108
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)
ΔS≔ S03− S01−12⋅ S02;
ΔS≔S03−S01−98.00573844
ΔG≔G03− G01−12⋅ G021000.0;
ΔG≔0.00100000G03−0.00100000G01+193712.15788334
ΔH≔ H03− H01−12⋅ H021000.0;
ΔH≔0.00100000H03−0.00100000H01+193682.93747242
NIST lists the following values for enthalpy of formation and free-energy of formation:
ΔHf (kJ/mol)
Δ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
205.15
Using Hess's Law, calculate ΔH, ΔG, and ΔS for Rxn (1) and compare to your calculated results.
ΔSHess≔243.0 − 219.3−12⋅205.15;
ΔSHess≔−78.87500000
ΔHHess≔−52.6−52.4;
ΔHHess≔−105.00000000
ΔGHess≔−51.7−61.4;
ΔGHess≔−113.10000000
S02;
196.01147687
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
Download Help Document
