Diatomic: An Open-Source Excel Application to Calculate Ther-Modynamic Properties for Diatomic Molecule

1. Introduction

Diatomic molecules are important in many areas of modern science, such as atmospheric chemistry [1,2], astrobiology [3] and ultracold physics [4]. Despite their simplicity, diatomic molecules are valuable models for understanding the relationship between structure and thermodynamic properties in more complex systems. For this purpose, it is important to develop high-quality theoretical methods to calculate these properties. Two main strategies have usually been adopted to improve the accuracy of these methods: (i) improving the description of the electronic structure [5,6,7,8] and (ii) implementing better models to describe the (rotational and vibrational) nuclear motions [9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24].

The electronic structure refinements include [5,6]: (i) a high-level quantum method, (ii) a robust basis set, (iii) the extrapolation to the complete basis set (CBS) limit, (iv) the inclusion of relativistic and spin-orbit effects and (v) a very good description of the electronic correlation.

The rotational structure is usually described by a classical approach, since the temperatures of interest are normally much higher than the characteristic rotational temperatures. However, at low temperatures it is important to adopt a quantum description of this structure to obtain accurate values [9].

At low and moderate temperatures, the vibrational structure can be reasonably described by the harmonic oscillator model [10]. However, this approach is known to overestimate the vibrational frequencies [10], leading to significant errors at high temperatures. To overcome this limitation, several correction methods have been applied [11,12,13,14,15,16,17,18,19,20,21,22,23,24]. Most of them can be classified into two groups: (i) methods that use empirical scaling factors to correct the harmonic frequencies [13,14] and (ii) methods that use anharmonic potentials [15,16,17,18,19,20,21,22,23,24].

The thermodynamic properties of a molecular system can be calculated using general quantum chemistry packages [25,26,27]. However, the required computational resources and level of expertise are considerable.

In this work, I have developed an open-source Excel application to calculate thermodynamic properties for diatomic ideal gases. This application, called Diatomic, is very easy to use and requires only some basic knowledge of theoretical chemistry. In fact, an user of Diatomic only has to input a limited number of molecular constants (the electronic ground state energy and the corresponding spin multiplicity, the atomic masses, the rotational symmetry number, the equilibrium bond length, the harmonic vibrational wave number, and the well depth of the Morse potential) to perform a calculation of this type. This information is freely available from online resources. The pressure and a set of fourteen temperatures must also to be inputted.

Diatomic provides many features, often not available in general quantum chemistry packages, to perform and interpret this type of calculations. These include: (i) two models (classical and quantum) describing the rotational structure, (ii) two models (harmonic and Morse) describing the vibrational structure, (iii) energies and Boltzmann populations for the quantum rotational model, (iv) energies and Boltzmann populations for both vibrational models implemented in the application, (v) molar thermodynamic properties as a function of temperature available in very easy to interpret and process tables, (vi) availability of these tables organized by different criteria.

Most of the refinements to the electronic structure mentioned above can be easily incorporated into Diatomic because they only affect the electronic energy. For this purpose it is necessary to search the literature for these corrected values. However, my goal in this work is not to obtain very accurate values for the molar thermodynamic properties. I try to show here that users can get good results for these properties using Diatomic and molecular data freely available online. Since Diatomic is an open-source collaborative project, anyone can add new vibrational potentials to this application.

Diatomic was used here to calculate these properties for fifteen ideal gases as a function of temperature and constant pressure. The chemical species studied correspond to all diatomic molecules that can be obtained by combining the five elements: H, N, O, F and Cl (see Table 1).

These molecules are classified according to their symmetry (homonuclear or heteronuclear) and spin multiplicity (singlet, duplet or triplet) (see Table 2).

2. Materials and Methods

2.1. General formalism

According to the well-established formalism of statistical mechanics [28,29,30], all the most relevant thermodynamic properties of an ideal gas can be calculated from its partition function. For a molar canonical ensemble1 and within the Born-Oppenheimer approach, this partition function (Q(N_A, V, T)) can be calculated from its molecular components:

$$Q(N_A,V,T)={\displaystyle \frac{Z^{N_A}}{N_A!}}={\left(Z_{el}\right)}^{N_A}\times {\displaystyle \frac{{\left(Z_{transl}\right)}^{N_A}}{N_A!}}\times {\left(Z_{rot}\right)}^{N_A}\times {\left(Z_{vib}\right)}^{N_A}$$

(1)

In Equation (1), N_A is the Avogadro number and Z is the molecular partition function. In the same equation, Z_el, Z_transl, Z_rot and Z_vib are the (electronic, translational, rotational and vibrational) components of the partition function. Within this formalism, the molar thermodynamic properties can be decomposed in a similar way:

$$A_m=A_m^{el}+A_m^{transl}+A_m^{rot}+A_m^{vib}\hspace{0.17em}\left(\mathit{molar}\mathit{Helmholtz}\mathit{energy}\right)$$

(2)

$$G_m=A_m+RT\iff G_m=A_m^{el}+A_m^{transl}+A_m^{rot}+A_m^{vib}+RT \left(\mathit{molar}\mathit{Gibbs}\mathit{energy}\right)$$

(3)

$${〈E〉}_m={〈E^{el}〉}_m\hspace{0.33em}\hspace{0.17em}+{〈E^{transl}〉}_m\hspace{0.33em}+{〈E^{rot}〉}_m\hspace{0.33em}+{〈E^{vib}〉}_m\hspace{0.33em}\left(\mathit{average}\mathit{value}\mathit{of}\mathit{molar}\mathit{energy}\right)$$

(4)

$${〈H〉}_m={〈H〉}_m+RT\iff {〈H〉}_m={〈E^{el}〉}_m\hspace{0.33em}\hspace{0.17em}+{〈E^{transl}〉}_m\hspace{0.33em}+{〈E^{rot}〉}_m\hspace{0.33em}+{〈E^{vib}〉}_m+RT\hspace{0.33em}\left(\mathit{average}\mathit{value}\mathit{of}\mathit{molar}\mathit{enthalpy}\right)$$

(5)

$$S_m=S_m^{el}+S_m^{transl}+S_m^{rot}+S_m^{vib}\hspace{0.17em}\left(\mathit{molar}\mathit{entropy}\right)$$

(6)

$$C_{v,m}=C_{v,m}^{el}+C_{v,m}^{transl}+C_{v,m}^{rot}+C_{v,m}^{vib} \left(\mathit{molar}\mathit{heat}\mathit{capacity}\mathit{at}\mathit{constant}\mathit{volume}\right)$$

(7)

$$C_{p,m}=C_{v,m}+R\iff C_{p,m}=C_{v,m}^{el}+C_{v,m}^{transl}+C_{v,m}^{rot}+C_{v,m}^{vib}+R \left(\mathit{molar}\mathit{heat}\mathit{capacity}\mathit{at}\mathit{constant}\mathit{pressure}\right)$$

(8)

In some of the previous equations ((3), (5) and (8)) an additional non-specific component (RT or R) must also be included. The components of the molar thermodynamic properties can be calculated as

$$A_m^{comp}=-RT\ln Z_{comp},\hspace{0.33em}comp=el, rot\hspace{0.33em}\mathrm{and} vib;A_m^{transl}=-RT\ln {\displaystyle \frac{Z_{transl}}{N_A}}-RT$$

(9)

$${〈E^{comp}〉}_m=R{T}^2{\left({\displaystyle \frac{\partial \ln Z_{comp})}{\partial T}}\right)}_V,\hspace{0.33em}comp=el,\hspace{0.33em}transl,\hspace{0.33em}rot \mathrm{and} vib\iff$$

(10)

or

$${〈E^{comp}〉}_m=N_A{\displaystyle \sum _i{P}_i\hspace{0.33em}{E}_i^{comp}};\hspace{0.33em}comp=el,\hspace{0.33em}transl,\hspace{0.33em}rot \mathrm{and}\hspace{0.33em}vib$$

(11)

Equation (10) should be used for partition functions (Z_comp) calculated analytically and Equation (11) for partition functions calculated numerically. In the last equation, $E_i^{comp}$ is the energy of the quantum state i and P_i is the corresponding Boltzmann population given by the following formula,

$$P_i={\displaystyle \frac{\exp \left(-\beta \hspace{0.17em}{E}_i^{comp}\right)}{Z_{comp}}};\mathrm{with} \beta =k_B T$$

(12)

In which k_b is the Boltzmann constant.

$$S_m^{comp}=R\ln Z_{comp}+RT{\left({\displaystyle \frac{\partial \left(\ln Z_{comp}\right)}{\partial T}}\right)}_V;\hspace{0.33em}comp=el,\hspace{0.33em}transl,\hspace{0.33em}rot \mathrm{and}\hspace{0.33em}vib$$

(13)

$$S_m^{transl}=R\ln {\displaystyle \frac{Z_{transl}}{N_A}}+R+RT{\left({\displaystyle \frac{\partial \left(\ln Z_{transl}\right)}{\partial T}}\right)}_V$$

(14)

$$C_{v,m}^{comp}={\left({\displaystyle \frac{\partial {〈E^{comp}〉}_m}{\partial T}}\right)}_V=N_A{\displaystyle \frac{{\sigma }_{E^{comp}}^2}{k_B{T}^2}};\hspace{0.33em}comp=el,\hspace{0.33em}transl,\hspace{0.33em}rot \mathrm{and}\hspace{0.33em}vib$$

(15)

In Equation (15), ${\sigma }_{E^{comp}}^2$ is the variance associated with the energy component (E^comp) calculated as:

$${\sigma }_{E^{comp}}^2=〈E^{com{p}^2}〉-{〈E^{comp}〉}^2;\hspace{0.33em}comp=el,\hspace{0.33em}transl,\hspace{0.33em}rot \mathrm{and}\hspace{0.33em}vib$$

(16)

In Equation (16), $〈E^{comp}〉$ and $〈E^{com{p}^2}〉$ are respectively the average values of the energy component ($E^{comp}$) and of its square (${E^{comp}}^2$) given by:

$$〈E^{com\mathrm{p}}〉={\displaystyle \sum _i{P}_i\hspace{0.33em}{E}_i^{comp}};\hspace{0.33em}comp=el,\hspace{0.33em}transl,\hspace{0.33em}rot \mathrm{and}\hspace{0.33em}vib$$

(17)

$${〈{E^{comp}}^2〉}_m={\displaystyle \sum _i{P}_i\hspace{0.33em}{E}_i^{com{p}^2}};\hspace{0.33em}comp=el,\hspace{0.33em}transl,\hspace{0.33em}rot \mathrm{and}\hspace{0.33em}vib$$

(18)

2.2. Classification of the molar thermodynamic properties

The molar thermodynamic properties are classified according to the criteria presented in Table 3.

Some properties of this type depend on the zero point energy (ZPE), which corresponds to the molar energy at 0 K. This quantity can be calculated as the sum of two components: the zero point electronic energy (ZPEE) and the zero point vibrational energy (ZPVE):

$$ZPE=ZPEE+ZPVE$$

(19)

$$ZPEE=N_A\hspace{0.17em}{E}_0$$

(20)

$$ZPVE=N_A\hspace{0.17em}{E}_0^{vib}$$

(21)

In the last two equations, $E_0$ and $E_0^{vib}$ are the energies of the electronic and vibrational ground states, respectively.

The Type 1 molar thermodynamic properties correspond to the absolute values of the thermodynamic properties that depend on the zero-point energy. A list of these properties is given in Table 4.

These properties are functions of temperature (T) and a limited number of molecular constants, which will be presented in section 2.10. The molar entropy (S_m) also depends on the pressure (P) due to its translational component ($S_m^{transl}$).

The type 2 molar thermodynamic properties correspond to the thermal values of the thermodynamic properties that depend on zero-point energy. A list of these properties is presented in Table 5.

A generic type 2 molar thermodynamic property ($X_{m, therm}$) can be calculated as,

$$X_{m, therm}=X_m-X_{m, ref};\hspace{0.33em}X=A, G, 〈\mathrm{E}〉\hspace{0.33em} \mathrm{or} 〈H〉$$

(22)

in which $X_m$ is the molar value of the thermodynamic property (X) at a generic temperature (T) and $X_{m, r\mathrm{ef}}$ is the corresponding molar value at a reference temperature (T_ref). In this work, two different reference temperatures were used. The original molar thermal values, calculated from the Diatomic application, use the reference temperature of 0 K. Some thermal values, presented in section 3, use the reference temperature of 298.15 K for consistency with the corresponding experimental values.

The type 3 molar thermodynamic properties are those that are independent of zero-point energy. A list of these properties is presented in Table 6.

2.3. Classical limit

The general formalism, presented in section 2.1, is based in a quantum formulation. In this case, the partition function associated with a component comp is calculated as a summation over the respective quantum levels (i),

$$Z_{comp}={\displaystyle \sum _i{g}_i}\hspace{0.17em}{E}_i^{comp};\hspace{0.33em}comp=el,\hspace{0.33em}transl,\hspace{0.33em}rot \mathrm{and}\hspace{0.33em}vib,$$

(23)

In which g_i is the degeneracy of the quantum level i and $E_i^{comp}$ is the corresponding energy.

However, when the temperature is high enough that the energy difference between two successive levels is much smaller than k_BT, the summation (23) can be replaced by an appropriate integral (classical limit). In this case, it is valid the so-called equipartition theorem, according to which each classical quadratic energy term contributes with: R/2 for the molar heat capacity at constant volume (C_v,m) and RT/2 for the average value of the molar energy (<E>_m).

Each energy component is associated with a characteristic temperature (Θ_comp), which determines the applicability of the classical limit and the equipartition (T >> Θ_comp). These characteristic temperatures vary considerably depending on the nature of the component. Typically, the electronic temperatures (Θ_el) are very high, the vibrational temperatures (Θ_vib) are moderate or high, the rotational temperatures (Θ_rot) are low and translational temperatures (Θ_transl) are very low.

This theorem is rarely applicable to the electronic component, because the energy differences between two successive electronic levels are usually too large. For diatomic molecules, when this theorem is fully applied to the nuclear motion components (translational, rotational, and vibrational), these molar thermodynamic properties can be estimated as:

• The translational component, associated with three normal modes, contributes with 3R/2 to C_v,m and with 3RT/2 to <E>_m.
• The rotational component, associated with two normal modes, contributes with R to C_v,m and with RT to <E>_m.
• The vibrational component2, associated with one normal mode, contributes with R to C_v,m and with RT to <E>_m.

A translational or rotational normal mode contributes with a single value to the referred molar thermodynamic properties (R/2 to C_v,m and RT/2 to <E>_m), because it is associated with a single quadratic kinetic energy term in classical mechanics. The vibrational normal mode gives double contributions to the mentioned molar thermodynamic properties (R to C_v,m and with RT to <E>_m), due to its association with two quadratic energy terms (one kinetic and one potential) in classical mechanics.

For a temperature interval [T_ref, T], where the classical limit and the equipartition are valid for the nuclear motion components but not for the electronic component, the following approximate equations are valid:

$${〈E〉}_{m, therm}\approx {\displaystyle \frac{7}{2}}R\hspace{0.17em}\hspace{0.17em}\left(T\hspace{0.33em}-T_{ref}\right)$$

(24)

$${〈H〉}_{m, therm}\approx {\displaystyle \frac{9}{2}}\hspace{0.17em}\hspace{0.17em}R\hspace{0.33em}\left(T\hspace{0.33em}-T_{ref}\right)$$

(25)

$$C_{v,m}\approx {\displaystyle \frac{7}{2}}R$$

(26)

$$C_{p,m}\approx {\displaystyle \frac{9}{2}}R$$

(27)

$$S_m=S_{m, ref}+{\displaystyle \underset{T_{ref}}{\overset{T}{\int }}\frac{C_{p,m}}{T}}\hspace{0.17em}dT\iff S_m\approx S_{m, ref}+{\displaystyle \frac{9}{2}}R\hspace{0.33em}\ln \left({\displaystyle \frac{T}{T_{ref}}}\right)$$

(28)

In equation (28), $S_{m, ref}$ is the molar entropy at a reference temperature T_ref.

2.4. Theoretical models

In the Diatomic application, the translational structure is always treated with the classical approach (see section 2.6), whereas the electronic structure is treated in a distinctive manner (see section 2.5). Two models (classical and quantum) have been implemented for the rotational structure (see section 2.7). For the vibrational structure, two models (harmonic oscillator and Morse oscillator) were also available (see section 2.8). In addition, an almost fully classical model (classical rotation, classical vibration) (see section 2.3) is also used in this work. These theoretical models are summarized in Table 7.

2.5. Electronic structure

For the molecules studied, it was assumed that the energy differences between the electronic ground state and the electronic excited states are very large. Under these circumstances, only the first state is occupied for the temperatures studied. Furthermore, in the calculation of the electronic partition function (Z_el), it is usually assumed that the energy of the electronic ground state (E₀) is the zero of the energy scale,

$$Z_{el}=g_0$$

(29)

in which g₀ is the spin multiplicity of the electronic ground state.

This implies that some corrections have to be introduced to the general formalism presented in sections 2.1 and 2.3:

$$A_m^{el}=-RT\ln g_0+N_A\hspace{0.17em}{E}_0$$

(30)

$${〈E^{el}〉}_m=N_A{E}_0=ZPEE$$

(31)

$$C_{v,m}^{el}={\left({\displaystyle \frac{\partial {〈E^{el}〉}_m}{\partial T}}\right)}_V=N_A{\displaystyle \frac{{\sigma }_{E^{el}}^2}{k_B{T}^2}}=0$$

(32)

$$S_m^{el}=R\ln g_0$$

(33)

2.6. Translational structure

The translational structure is fully treated by the classical approach. Thus, the translational partition constant and the translational components of some relevant molar thermodynamic properties can be calculated, respectively, as follows:

$$Z_{transl}={\displaystyle \frac{V}{{\mathsf{\Lambda }}^3}}$$

(34)

$$A_m^{transl}=-RT\ln {\displaystyle \frac{k_{B}T}{P{\mathsf{\Lambda }}^3}}-RT$$

(35)

$${〈E^{transl}〉}_m={\displaystyle \frac{3}{2}}RT$$

(36)

$$C_{v,m}^{transl}={\displaystyle \frac{3}{2}}R$$

(37)

$$S_m^{transl}=R\ln {\displaystyle \frac{k_{B}T}{P{\mathsf{\Lambda }}^3}}+{\displaystyle \frac{5}{2}}R$$

(38)

In equations (34), (35) and (38), Λ is the thermal de Broglie wavelength given by,

$$\mathsf{\Lambda }={\left({\displaystyle \frac{h^2}{2\pi \hspace{0.17em}m\hspace{0.17em}{k}_{B}T}}\right)}^{1/2},$$

(39)

in which m is the molecular mass and h is the Planck constant.

2.7. Rotational structure

The rotational structure is treated using two alternative approaches (quantum and classical). Within the quantum formalism, the rotational partition constant can be calculated as:

$$\begin{array}{l}{Z}_{rot}={\displaystyle \frac{1}{{\sigma }_{rot}}}{\displaystyle \sum _{J=0}^{\infty }{g}_J\exp \left(-\frac{E_J^{rot}}{k_B\hspace{0.17em}T}\right)}, \mathrm{with} g_J=2J+1\iff Z_{rot}={\displaystyle \frac{1}{{\sigma }_{rot}}}{\displaystyle \sum _{J=0}^{\infty }{g}_J}\exp \left({\displaystyle \frac{-B\hspace{0.17em}J\left(J+1\right)}{k_{B}T}}\right)\iff \\ Z_{rot}={\displaystyle \frac{1}{{\sigma }_{rot}}}{\displaystyle \sum _{J=0}^{\infty }(2J+1)}\exp \left({\displaystyle \frac{-{\mathsf{\Theta }}_{rot}J\left(J+1\right)}{T}}\right)\end{array}$$

(40)

In the Diatomic application, a finite number of terms (J_max = 2000) are used in the summations included in equation (40). However, this was enough to obtain convergence in the rotational partition for all the calculations performed here.

In this equation, σ_rot is the rotational symmetry number (σ_rot = 2 for diatomic homonuclear molecules and σ_rot = 1 for diatomic heteronuclear molecules),$E_J^{rot}$ is the energy of the rotational quantum level J, $g_J$ is the respective degeneracy, B is the rotational constant and ${\mathsf{\Theta }}_{rot}$ is the rotational temperature. The last two quantities can be calculated using the following equations:

$$B={\displaystyle \frac{{\hslash }^2}{2I}}, \mathrm{with} \hslash ={\displaystyle \frac{h}{2\pi }}$$

(41)

$${\mathsf{\Theta }}_{rot}={\displaystyle \frac{B}{k_B}}$$

(42)

In equation (41), I is the moment of inertia given by,

$$I=\mu \hspace{0.17em}{r}_{eq}^2,\hspace{0.17em}\mu ={\displaystyle \frac{m_1\times m_2}{m}},$$

(43)

in which μ is the reduced mass of the diatomic molecule, m₁ the mass of its first atom, m₂ the mass of its second atom and r_eq its equilibrium bond length.

Within this formalism, the rotational components of some relevant molar thermodynamic properties can be calculated as:

$$A_m^{rot}=-RT\ln Z_{rot}\iff A_m^{rot}=-RT\ln \left[{\displaystyle \frac{1}{{\sigma }_{rot}}}{\displaystyle \sum _{J=0}^{J_{\max }}(2J+1)}\exp \left({\displaystyle \frac{-{\mathsf{\Theta }}_{rot}J\left(J+1\right)}{T}}\right)\right]; \mathrm{with} J_{\max }=2000$$

(44)

$${〈E^{\mathrm{rot}}〉}_{m,T}=N_A{\displaystyle \sum _{J=0}^{J_{\max }}{P}_{\mathrm{J}}\hspace{0.17em}{E}_J^{rot}; \mathrm{with} P_J=}(2J+1){\displaystyle \frac{\exp \left(-\beta \hspace{0.17em}{E}_J^{rot}\right)}{{\sigma }_{rot}\hspace{0.17em}{Z}_{rot}}}$$

(45)

$$C_{v,m}^{rot}=N_A{\displaystyle \frac{{\sigma }_{E^{rot}}^2}{k_B{T}^2}}=N_A{\displaystyle \frac{\left(〈E^{rot{ }^2}〉-{〈E^{ro\mathrm{t}}〉}^2\right)}{k_B{T}^2}}$$

(46)

In Equation (46), $〈E^{rot}〉$ and $〈E^{ro{t}^2}〉$ are the average values of the rotational energy ($E^{rot}$) and of its square (${E^{rot}}^2$), respectively, given by:

$$〈E^{rot}〉={\displaystyle \sum _{J=0}^{J_{\max }}{P}_J\hspace{0.33em}{E}_J^{rot}}$$

(47)

$$〈{E^{rot}}^2〉={\displaystyle \sum _{J=0}^{J_{\max }}{P}_J\hspace{0.33em}{E}_J^{ro{t}^2}}$$

(48)

$$S_m^{rot}={\displaystyle \frac{{〈E^{rot}〉}_m-A_m^{rot}}{T}}$$

(49)

Within the classical formalism, equations (40) to (49) must be reformulated as:

$$Z_{rot}={\displaystyle \frac{T}{{\sigma }_{rot}\hspace{0.17em}{\mathsf{\Theta }}_{rot}}}$$

(50)

$$A_m^{rot}=-RT\ln Z_{rot}\iff A_m^{rot}=-RT\ln \left({\displaystyle \frac{T}{{\sigma }_{rot}\hspace{0.17em}{\mathsf{\Theta }}_{rot}}}\right)$$

(51)

$${〈E^{rot}〉}_m=RT$$

(52)

$$C_{v,m}^{rot}=R$$

(53)

$$S_m^{rot}=R\hspace{0.33em}\left(\ln \left({\displaystyle \frac{T}{{\sigma }_{rot}\hspace{0.17em}{\mathsf{\Theta }}_{rot}}}\right)+1\right)$$

(54)

2.8. Vibrational structure

The vibrational structure is treated using two alternative approaches (harmonic oscillator and Morse oscillator). Within the harmonic oscillator formalism, the vibrational partition constant can be calculated as,

$$Z_{vib}={\displaystyle \frac{\exp \left[-\beta \hspace{0.17em}h{\nu }_{vib}/2\right]}{1-\exp \left[-\beta \hspace{0.17em}h{\nu }_{vib}\right]}}\iff Z_{vib}={\displaystyle \frac{\exp \left[-h{\nu }_{vib}/\left(2k_{B}T\right)\right]}{1-\exp \left[-h{\nu }_{vib}/\left(k_{B}T\right)\right]}}\iff Z_{vib}={\displaystyle \frac{\exp \left[-{\mathsf{\Theta }}_{vib}/\left(2T\right)\right]}{1-\exp \left[-{\mathsf{\Theta }}_{vib}/T\right]}},$$

(55)

where ν_vib is the vibrational frequency and Θ_vib the vibrational temperature, given respectively by,

$${\nu }_{vib}={\displaystyle \frac{1}{2\pi }}{\left({\displaystyle \frac{K}{\mu }}\right)}^{1/2}$$

(56)

in which Κ is the force constant of the diatomic molecule,

$${\mathsf{\Theta }}_{vib}={\displaystyle \frac{h{\nu }_{vib}}{k_B}}$$

(57)

Within this formalism, the vibrational components of some relevant molar thermodynamic properties can be calculated as:

$$A_m^{vib}=-RT\ln Z_{vib}\iff A_m^{vib}=-RT\ln \left({\displaystyle \frac{\exp \left[-{\mathsf{\Theta }}_{vib}/\left(2T\right)\right]}{1-\exp \left[-{\mathsf{\Theta }}_{vib}/T\right]}}\right)$$

(58)

$${〈E^{vib}〉}_m=R{\displaystyle \frac{{\mathsf{\Theta }}_{vib}}{2}}+R{\displaystyle \frac{{\mathsf{\Theta }}_{vib}}{\exp \left({\mathsf{\Theta }}_{vib}/T\right)-1}}\iff {〈E^{vib}〉}_m={\displaystyle \frac{N_A}{2}}h{\nu }_{vib}+{\displaystyle \frac{N_{A}h{\nu }_{vib}}{\exp \left[h{\nu }_{vib}/\left(k_{B}T\right)\right]-1}}$$

(59)

$$C_{v,m}^{vib}=N_A{\displaystyle \frac{{\sigma }_{E^{vib}}^2}{k_B{T}^2}}=R\hspace{0.33em}\left[\exp \left({\mathsf{\Theta }}_{vib}/T\right) {\left({\displaystyle \frac{{\mathsf{\Theta }}_{vib}/T}{\exp \left({\mathsf{\Theta }}_{vib}/T\right)-1}}\right)}^2\right]$$

(60)

$$S_m^{vib}=R\hspace{0.33em}\left({\displaystyle \frac{{\mathsf{\Theta }}_{vib}/T}{\exp \left({\mathsf{\Theta }}_{vib}/T\right)-1}}-\ln \left[1-\exp \left(-{\mathsf{\Theta }}_{vib}/T\right)\right]\right)$$

(61)

Within the Morse oscillator formalism, equations (55) to (61) must be reformulated as,

$$Z_{vib}={\displaystyle \sum _{v=0}^{v_{\max }}\exp \left(-\beta \hspace{0.17em}{E}_v^{vib}\right)},$$

(62)

where $E_v^{vib}$ is the vibrational energy of a diatomic molecule according to this formalism given by the equation (63):

$$E_v^{vib}=\left(v+{\displaystyle \frac{1}{2}}\right)\hspace{0.17em}\hspace{0.33em}h{\nu }_{vib}-{\left(v+{\displaystyle \frac{1}{2}}\right)}^2\hspace{0.17em}{\chi }_e\hspace{0.17em}h{\nu }_{vib}; v=0, 1, 2, 3, \dots , {v_{\max }}^3$$

(63)

In Equation (63), v is the vibrational quantum number, v_max+1 is the number of bound states associated with the Morse potential and χ_e is the anharmonicity constant. The last two quantities can be calculated by equations (64) and (65), respectively:

$$v_{\max }={\displaystyle \frac{1-{\chi }_e}{2{\chi }_e}}$$

(65)

$${\chi }_{\mathrm{e}}={\displaystyle \frac{h{\nu }_{vib}}{4D_e}}$$

(66)

In Equation (66), D_e is the well depth of the Morse potential.

$$A_m^{vib}=-RT\ln Z_{vib}\iff A_m^{vib}=-RT\ln \left[{\displaystyle \sum _{v=0}^{v_{\max }}\exp \left(-\beta \hspace{0.17em}{E}_{\mathrm{v}}^{\mathrm{vib}}\right)}\right]$$

(67)

$${〈E^{vib}〉}_m=N_A{\displaystyle \sum _{v=0}^{v_{\max }}{P}_v\hspace{0.17em}{E}_v^{vib}; \mathrm{with} P_v=}{\displaystyle \frac{\exp \left(-\beta \hspace{0.17em}{E}_v^{vib}\right)}{\hspace{0.17em}{Z}_{vib}}}$$

(68)

$$C_{v,m}^{vib}=N_A{\displaystyle \frac{\left(〈E^{v{\mathrm{ib} }^2}〉-{〈E^{vib}〉}^2\right)}{k_B{T}^2}}$$

(69)

In Equation (69), $〈E^{vib}〉$ and $〈E^{vi{b}^2}〉$ are the average values of the vibrational energy ($E^{vib}$) and of its square (${E^{vib}}^2$), respectively, given by:

$$〈E^{vib}〉={\displaystyle \sum _{v=0}^{v_{\max }}{P}_v\hspace{0.33em}{E}_v^{vib}}$$

(70)

$$〈{E^{vib}}^2〉={\displaystyle \sum _{v=0}^{v_{\max }}{P}_v\hspace{0.33em}{E}_v^{vi{b}^2}}$$

(71)

$$S_m^{vib}={\displaystyle \frac{{〈E^{vib}〉}_m-A_m^{vib}}{T}}$$

(72)

2.9. Standard molar enthalpy of formation

The standard molar enthalpy of formation of a diatomic heteronuclear gas AB (Δ_fH^o(AB, g)) is defined as the enthalpy variation, under standard conditions, associated with the following chemical equation:

½ A₂(g) + ½ B₂(g) ⇌ AB(g)

(73)

$${\mathsf{\Delta }}_f{H}^{\mathrm{o}}(\mathrm{AB}, \mathrm{g})={〈H^{\mathrm{o}}(\mathrm{AB}, \mathrm{g})〉}_m-{\displaystyle \frac{1}{2}}\left({〈H^{\mathrm{o}}({\mathrm{A}}_2, \mathrm{g})〉}_m+〈H^{\mathrm{o}}({\mathrm{B}}_2, \mathrm{g})〉\right)$$

(74)

In equation (74), the standard molar enthalpy of a diatomic gas X (<H^o (X, g)>_m; X= AB, A₂ or B₂) is the average value of the standard molar enthalpy of this chemical species at a pressure 1 bar and a specific temperature (usually 298.15 K).

2.10. Molecular constants

In the theoretical models considered, the molar thermodynamic properties depend only on a limited number of molecular constants (primary molecular constants). However, in many cases, the formalism previously presented in this section was expressed using other molecular constants (secondary molecular constants) that can be obtained from the primary ones. A list of molecular constants, used in this work, is presented in Table 8.

2.11. Diatomic application

A flowchart of the Diatomic application is shown in Figure 1.

This application is organized in three blocks of spreadsheets: block of data, block of preliminary calculations and block of results.

Yhe block of data contains the following information:

• Universal constants and conversion factors (provided by the application in the Universal_Constants spreadsheet);
• Pressure value used in the calculations (inputted by the user in the Universal_Constants spreadsheet);
• Primary molecular constants and 14 temperature values used in calculations (inputted by the user in the Molecular_Constants_1 spreadsheet4);
• Secondary molecular constants (calculated by the application in the Molecular_Constants_2 spreadsheet).

The block of preliminary calculations performs the following calculations:

• Rotational partition function using the quantum formulation (performed by the Quantum_Rotational_Part_Funct spreadsheet);
• Boltzmann populations and energies, for the rotational levels, using the quantum formulation (performed by the Quant_Rotational_Boltzmann_Pop spreadsheet);
• Vibrational partition function using the Morse oscillator formulation (calculated by the Morse spreadsheet);
• Vibrational partition function using the harmonic oscillator formulation (calculated by the Vib_Populations_Energies spreadsheet);
• Boltzmann populations and energies, for the vibrational states, using both the harmonic and Morse oscillators formulations (calculated by the Vib_Populations_Energies spreadsheet);
• Average vibrational molar energy using the Morse oscillator formulation (performed by the Vib_Populations_Energies spreadsheet);
• Vibrational molar heat capacity at constant volume using the Morse oscillator formulation (performed by the Vib_Populations_Energies spreadsheet);
• Canonical ensemble total partition function and its components, using the various theoretical models available (calculated by the Partition_Functions spreadsheet).

In the block of results, the results are presented/calculated using different criteria:

• In a first level, the results are presented/calculated by type as defined in Table 3 (Therm_Prop_1, Therm_Prop_2 and Therm_Prop_3 spreadsheets);
• At a second level, the results are presented/calculated by type as defined in Table 3 and by theoretical model as defined in Table 7 (type 1 molar thermodynamic properties with different theoretical models: Harm_Class_Therm_Prop_1, Morse_Class_Therm_Prop_1, Harm_Quantum_Therm_Prop_1 and Morse_Quantum_Therm_Prop_1 spreadsheets; type 2 molar thermodynamic properties with different theoretical models: Harm_Class_Therm_Prop_2, Morse_Class_Therm_Prop_2, Harm_Quantum_Therm_Prop_2 and Morse_Quantum_Therm_Prop_2 spreadsheets; type 3 molar thermodynamic properties with different theoretical models: Harm_Class_Therm_Prop_3, Morse_Class_Therm_Prop_3, Harm_Quantum_Therm_Prop_3 and Morse_Quantum_Therm_Prop_3 spreadsheets);
• In a third level, the results are presented/calculated by theoretical model as defined in Table 7 (harm, class: Harm_Class_All_Therm_Prop spreadsheet, Morse, class: Morse_Class_All_Therm_Prop spreadsheet, harm, quant: Harm_Quant_All_Therm_Prop spreadsheet, Morse, quant: Morse_Quant_All_Therm_Prop spreadsheet.

3. Results

In this study, the Diatomic application was used to calculate the molar thermodynamic properties presented in Table 3, Table 4, Table 5 and Table 6 for the molecules listed in Table 1 and Table 2. The first four models presented in Table 7 (harm, class; Morse, class; harm, quant and Morse, quant) were employed in these calculations. The molecular constants, calculated at a CCSD(T)/cc-pVQZ quantum level, were obtained from the Computational Chemistry Comparison and Benchmark DataBase [31]. These constants are presented in Table 9 and Table 10. These properties were calculated for 14 different temperatures (298.15 K, 373.15 K, 473.15 K, 573.15 K, 673.15 K, 773.15 K, 873.15 K, 973.15 K, 1073.15 K, 1173.15 K, 1273.15 K, 1373.15 K, 1473.15 K and 1573.15 K) and at a constant pressure of 1 bar. Thus, these are standard molar thermodynamic properties at different temperatures. These results are available as supplementary material (see Table S1 and files referred therein).

For comparison, these properties were also calculated using Gaussian 09 [25]. For each molecule studied, a geometry optimization with harmonic frequency calculation was performed. The maximum deviations between the results obtained by these two applications are 0.01, which is in the order of the precision of the thermal results provided by Gaussian 09 (two decimal places).

The standard molar enthalpies of formation at 298.15 K, for the diatomic heteronuclear gases studied, were further calculated using the Diatomic application and compared with the corresponding experimental values. The results obtained are presented in Table 11.

Subsequently, three additional molar thermodynamic properties (molar heat capacity at constant pressure (C_p,m), average value of the standard molar thermal enthalpy (<H^o>_{m, therm}) and standard molar entropy (S^o_m)) were calculated using the Diatomic application for all the molecules studied and compared with corresponding experimental values [32,33]. A reference temperature of 298.15 K was chosen for both the experimental and calculated results. The five theoretical models, presented in Table 7, were tested in these calculations. The temperature and pressure constraints were the same as those used in the global calculations. A reference temperature of 298.15 K was selected for both experimental and calculated results. The relative errors for the calculated values were determined, using the experimental values as reference. No experimental values are available for the nitrogen monochloride (NCl) molecule. The results obtained are presented in Tables 12 through 23.

The convergence of the molar heat capacity at constant pressure (C_p,m), calculated using the (harm, class) theoretical model, to the classical limit as a function of temperature is presented in Figure 2 for further analysis.

4. Discussion

The calculated and experimental values of standard molar enthalpy of formation, for the diatomic heteronuclear gas under study, are in reasonable agreement (see table 11). In addition, the following general tendencies can be deduced from the results:

• There is no significant difference in the performance of the 4 theoretical models used in these calculations.
• The relative errors obtained increase with the spin multiplicity of the molecular systems studied.

For the other three molar thermodynamic properties studied (molar heat capacity at pressure (C_p,m), average value of the standard molar thermal enthalpy (<H^o>_{m, therm}) and standard molar entropy (S^o_m)), the agreement between the experimental and calculated results is generally much better (see Table 12, Table 13, Table 14, Table 15, Table 16, Table 17, Table 18, Table 19, Table 20, Table 21, Table 22 and Table 23). Furthermore, the following general tendencies can be deduced from the results:

• The theoretical model, utilized to describe the rotational structure (class or quant), has an almost negligible effect on the results obtained. Only for hydrogen (which exhibits the highest rotational temperature in the investigated series) some residual effects occur at the lowest temperatures.
• For the lowest temperatures, the Morse oscillator and harmonic oscillator vibrational models have similar performances (similar relative errors).
• For moderate and high temperatures, the Morse oscillator model generally gives better results (smaller relative errors) than the harmonic oscillator model. This performance difference usually increases with temperature.
• The almost fully classical (class, class) model, which is based on the equipartition principle, usually gives the worst results (largest relative errors) in the calculation of the molar heat capacity at pressure and the average value of the standard molar thermal enthalpy.
• However, for the calculation of the standard molar entropy, this model gives similar results to the others. In fact, even for some open-shell systems (ClO, OF, NF, NO, and NH), it is usually associated with the best results (smallest relative errors).

The variations of the molar heat capacity at constant pressure (C_p,m), with temperature and the molecular structure, can be rationalized as follows:

• The general tendencies can be deduced from the (harm, class) model, with some anharmonic corrections at high temperatures.
• The quantum rotational effects are negligible, because the temperatures under study (T ≥ 298.15 K) are much higher than the typical rotational temperatures. The highest one of these is 87.44 K, associated with the hydrogen molecule.
• The equipartition principle is fully applied to the translational and rotational components of C_p,m. For all the molecules studied, the “classical” (translational + rotational + nonspecific) component is (3R/2 + R + R = 7R/2).
• Within the model adopted for the electronic structure, the electronic component ($C_{p,m}^{el}$) of this molar thermodynamic property is null (see equation 32).
• The corresponding vibrational component ($C_{p,m}^{vib}$) increases with temperature, because the variance of the vibrational energy (${\sigma }_{E^{vib}}^2$) follows the same tendency. In fact, the Boltzmann population of the vibrational ground state decreases and the Boltzmann populations of the vibrational excited states increase with temperature.
• These tendencies are favored by the decrease of the vibrational frequency (ν_vib), because this decreases the energy difference between two successive vibrational states.
• According to the harmonic model, the equipartition principle is also valid for the vibrational components of C_p,m at high or very high temperatures. Thus, this molar thermodynamic property tends to (3R/2 + R + R + R = 9R/2).
• This convergence to the classical limit is accelerated by the decrease of the vibrational frequency (ν_vib).
• These tendencies are illustrated by results presented Figure 2.
• At high temperatures, however, the anharmonic effects cannot be neglected. Therefore, the (Morse, class) model is important to understand the deviations from the harmonic behavior inherent in this situation.

The variations on the average value of the standard molar thermal enthalpy (<H^o>_{m, therm}), with temperature and the molecular structure, follow similar tendencies already described for the molar heat capacity at constant pressure. The main reasons for this behavior are:

• The two properties are closely correlated.
• The equipartition principle also applies here, despite the classical limit of <H^o>_{m, therm} is a function (9R (T-T_ref)/2) and not a constant (9R/2).

The factors, influencing the variations on C_p,m and <H^o>_m with temperature and the molecular structure, are also crucial to rationalize the corresponding variations on the standard molar entropy (S^o_m). However, this thermodynamic property is much more complex and new factors have to be included in the analysis for a full understanding of its behavior. For this purpose, the following issues must be considered:

• The equipartition principle is not directly applied to molar entropy. Even the approximate equation (28), derived for an interval [T_ref, T] where this principle is fully applied to the molar heat capacity at constant pressure and to the average value of molar thermal enthalpy, is not universal. In fact the reference molar entropy ($S_{m, ref}$) is a specific parameter of each molecule.
• The translational molar entropy increases with molecular mass m (see equations 38 and 39).
• The rotational molar entropy increases with the reduced mass μ, the equilibrium bond length r_eq and decreases with the rotational symmetry number σ_rot (see equations 41, 42, 43 and 54).
• The electronic molar entropy increases with the spin multiplicity of the electronic ground state g₀ (see equation 33).
• Consequently, the three components of the molar entropy above mentioned are also specific to each molecule.

5. Conclusions

In this paper, I presented Diatomic as an alternative to quantum chemistry packages for calculating thermodynamic properties of diatomic molecules. Diatomic is an open-source Excel application that requires only a limited number of molecular constants to perform these calculations. In addition, these constants can be easily obtained from freely available online databases.

In this new application, two theoretical models (Morse and harmonic) are available to describe the vibrational structure, and two theoretical models (classical and quantum) are available to describe the rotational structure. By combining the rotational and vibrational structures, four theoretical models can be used within Diatomic (harm-class, Morse-class, harm-quant, and Morse-quant). A fifth model (class-class), based on the equipartition principle, is also used here.

Diatomic was used to calculate four standard molar thermodynamic properties (standard molar enthalpy of formation, molar heat capacity at constant pressure, average value of the standard molar thermal enthalpy, and standard molar entropy) for a set of fifteen diatomic molecules. The results obtained with this application were compared with experimental data.

The theoretical model (classical or quantum) used to describe the rotational structure has a very limited effect on the accuracy of the calculated results. For moderate and high temperatures, the Morse model has better performance (smaller relative errors) than the harmonic model. At low temperatures, both models give similar results.

The results obtained with the molar heat capacity at constant pressure and the average value of the standard molar thermal enthalpy can be rationalized from a harmonic point of view (emphasizing the relevance of the vibrational frequency and the temperature parameters in the analysis), with some anharmonic corrections to the highest temperatures. However, to interpret the standard molar entropy results, some additional factors have to be included in the analysis. These include molecular constants (spin multiplicity of the electronic ground state, molecular mass, reduced mass, equilibrium bond length and rotational symmetry number) that affect the electronic, translational and rotational components of this molar thermodynamic property.