Appendix 1: An introduction of thermodynamic data calculation

In order to make users understand how Shermo works, an overview of calculation of thermochemistry quantities is given in this section.

A.1 Basic concept

The thermochemistry quantities of common interest include internal energy (U), enthalpy (H), entropy (S), Gibbs free energy (G), constant volume heat capacity (C_V) and constant pressure heat capacity (C_P). The partition function (q), which plays a key role in statistical thermodynamics, is also frequently investigated.

Below equations can be found in any physical chemistry textbook

\[ \begin{aligned} H(T) = U(T) + P \times V \\ G(T) = H(T) - T \times S(T) \end{aligned} \]

where \(T\) is temperature, \(P\) is pressure, \(V\) is volume.

In Shermo, all calculations are based on ideal gas assumption and are reported for 1 mole substance. In this case, PV=RT and C_P=C_V+R, where R= 8.3144648 J·(mol·K)⁻¹ is ideal gas constant.

The U, H and G can be written as sum of electronic energy (\(\varepsilon_{\rm ele}\)) and corresponding thermal corrections:

\[ \begin{aligned} & U(0) = H(0) = G(0) = \varepsilon_{\rm ele} + {\rm ZPE} \\ & U(T) = \varepsilon_{\rm ele} + U_{\rm corr}(T) \\ & H(T) = \varepsilon_{\rm ele} + H_{\rm corr}(T) \\ & G(T) = \varepsilon_{\rm ele} + G_{\rm corr}(T) \end{aligned} \]

where ZPE is acronym of zero-point energy, which is sometimes also referred to as ZPVE (zero-point vibrational energy).

At finite temperature, the thermal corrections can also be written as sum of ZPE and contribution by heating the system from 0 K to current temperature

\[ \begin{aligned} & U_{\rm corr}(T) = {\rm ZPE} + U_{0\to T} = {\rm ZPE} + [U(T)-U(0)] = {\rm ZPE} + \int C_V(T)\,{\rm d}T \\ & H_{\rm corr}(T) = {\rm ZPE} + H_{0\to T} = {\rm ZPE} + [H(T)-H(0)] = {\rm ZPE} + \int C_P(T)\,{\rm d}T \\ & G_{\rm corr}(T) = {\rm ZPE} + G_{0\to T} = {\rm ZPE} + [G(T)-G(0)] \end{aligned} \]

The thermochemistry quantities can be approximated as sum of contributions from translation (trans), rotation (rot), vibration (vib) and electron (ele) by ignoring their coupling

\[ \begin{aligned} & U_{\rm corr} = U_{\rm trans} + U_{\rm rot} + U_{\rm vib} + U_{\rm ele} \\ & S = S_{\rm trans} + S_{\rm rot} + S_{\rm vib} + S_{\rm ele} \\ & C_V = C_{V,\rm trans} + C_{V,\rm rot} + C_{V,\rm vib} + C_{V,\rm ele} \end{aligned} \]

Similarly for H_corr(T) and G_corr(T).

The H_corr and U_corr, as well as C_P and C_V, only differ by their translation parts

\[ \begin{aligned} & H_{\rm trans} = U_{\rm trans} + RT \\ & C_{P,\rm trans} = C_{V,\rm trans} + R \end{aligned} \]

Partition function (q) is a dimensionless quantity playing crucial role in statistical thermodynamic theory, it represents average number of thermal accessible quantum states. Formally, all aforementioned quantities can be derived based on q of a molecule, for example

\[ \begin{aligned} & U = kT^2\left(\frac{\partial \ln Q}{\partial T}\right)_V \\ & S = k\ln Q + kT\left(\frac{\partial \ln Q}{\partial T}\right)_V \end{aligned} \]

where \(Q\) is ensemble partition function. Under ideal gas assumption, for an ensemble consisting of \(N\) molecules, the \(Q\) relates to molecular \(q\) with

\[ Q = \frac{q^N}{N!} \]

When coupling between translation, rotation, vibrations and electron is ignored, q can be expressed as

\[ \begin{aligned} & q_{\rm tot} = q_{\rm trans}q_{\rm rot}q_{\rm vib}q_{\rm ele} \\ & = \sum_a^{\rm trans} \exp\left(\frac{-\varepsilon_{\rm trans}^a}{kT}\right) \sum_b^{\rm rot} \exp\left(\frac{-\varepsilon_{\rm rot}^b}{kT}\right) \sum_c^{\rm vib} \exp\left(\frac{-\varepsilon_{\rm vib}^c}{kT}\right) \sum_d^{\rm ele} \exp\left(\frac{-\varepsilon_{\rm ele}^d}{kT}\right) \end{aligned} \]

where \(\varepsilon\) denotes energy at a molecular level. For example, \(\varepsilon_{\rm rot}^b\) denotes molecular rotation state \(b\). The \(k\) is Boltzmann constant \(1.3806503 \times 10^{-23}\) J·K⁻¹.

Next, let us see how the aforementioned quantities are calculated, all working equations are given. Unless otherwise specified, unit of U, H and G is J·mol⁻¹, unit of S, C_V and C_P is J·(mol·K)⁻¹.

A.2 Translation contribution

Contribution of molecular overall translation motion to various thermodynamic data can be well calculated as follows if temperature is not extremely low (<5 K for most molecules)

\[ \begin{aligned} & q_{\rm trans} = \left(\frac{2\pi mkT}{h^2}\right)^{3/2}\frac{RT}{P} \\ & U_{\rm trans}(T) = (3/2)RT \\ & H_{\rm trans}(T) = U_{\rm trans}(T) + RT \\ & C_{V,\rm trans} = (3/2)R \\ & C_{P,\rm trans} = C_{V,\rm trans} + R \\ & S_{\rm trans}(T) = R[\ln(q_{\rm trans}/N_{\rm A}) + 5/2] \end{aligned} \]

where \(h\) is Planck constant, \(6.6260696 \times 10^{-34}\) J·s, \(N_{\rm A}\) is Avogadro constant \(6.02214179 \times 10^{23}\) , \(m\) is mass of current molecule. It is noteworthy that in all thermochemistry quantities, pressure only affects translation contribution to entropy (S_trans).

A.3 Rotation contribution

Molecular rotation also has contribution to various thermodynamic data. Of course, a single atom does not have rotation contribution.

Linear molecules:

\[ \begin{aligned} & q_{\rm rot} = \frac{8\pi^2 IkT}{\sigma h^2} \\ & U_{\rm rot}(T) = RT \\ & C_{V,\rm rot}(T) = R \\ & S_{\rm rot}(T) = R(\ln q_{\rm rot} + 1) \end{aligned} \]

Nonlinear molecules:

\[ \begin{aligned} & q_{\rm rot} = \frac{8\pi^2}{\sigma h^3}(2\pi kT)^{3/2}\sqrt{I_AI_BI_C} \\ & U_{\rm rot}(T) = (3/2)RT \\ & C_{V,\rm rot}(T) = (3/2)R \\ & S_{\rm rot}(T) = R(\ln q_{\rm rot} + 3/2) \end{aligned} \]

where \(\sigma\) is rotational symmetry number, \(\{I\}\) are principal moments of inertia.

The \(\sigma\) is determined according to point group:

Point group	\(\sigma\)	Point group	\(\sigma\)
C1, Ci, Cs, Cv	1	Cn, Cnv, Cnh	n
Dh	2	Dn, Dnh, Dnd	2n
Ih	60	Sn	n/2
Oh	24	T, Td	12

It is worth to note that in Shermo, the point group is detected based on current geometry based on my adapted SYVA code. The original SYVA code was described in Comput. Phys. Commun., 215, 156 (2017).

The principal moments of inertia are the three eigenvalues of the moments of inertia tensor:

\[ \mathbf{I} = \left[\begin{matrix} \sum_i m_i(y_i^2+z_i^2) & -\sum_i m_ix_iy_i & -\sum_i m_ix_iz_i \\ -\sum_i m_iy_ix_i & \sum_i m_i(x_i^2+z_i^2) & -\sum_i m_iy_iz_i \\ -\sum_i m_iz_ix_i & -\sum_i m_iz_iy_i & \sum_i m_i(x_i^2+y_i^2) \end{matrix}\right] \]

where \(x_i\), \(y_i\), \(z_i\) are Cartesian coordinates of atom \(i\), the \(m_i\) is mass of atom \(i\).

For a linear molecule, two of its principal moments of inertia are identical, and the remaining one is zero. Shermo determines whether or not a molecule is linear by checking its lowest principal moment of inertia. If it is less than 0.001 k·gm², then the molecule will be regarded as linear.

The rotation constants are also outputted by Shermo, they are closely related to principal moments of inertia:

\[ A = \frac{h}{8\pi^2I_A} \qquad B = \frac{h}{8\pi^2I_B} \qquad C = \frac{h}{8\pi^2I_C} \]

If the unit of \(I\) is kg·m², then the unit of the corresponding rotation constants will be Hz.

Shermo also outputs rotation temperatures:

\[ T_A^{\rm rot} = \frac{hA}{k} \qquad T_B^{\rm rot} = \frac{hB}{k} \qquad T_C^{\rm rot} = \frac{hC}{k} \]

A.4 Vibration contribution

Every real frequency has contribution to partition function and various thermodynamic quantities, while imaginary frequency modes are ignored during calculation. Harmonic approximation is assumed below; this is also known as rigid-rotor harmonic oscillator (RRHO) model.

Total vibration partition function is product of contributions of individual vibrational modes:

\[ \begin{aligned} q_{\rm vib}^{(\rm BOT)} = \prod_i \frac{e^{-h\nu_i/(2kT)}}{1-e^{-h\nu_i/(kT)}} \\ q_{\rm vib}^{(V=0)} = \prod_i \frac{1}{1-e^{-h\nu_i/(kT)}} \end{aligned} \]

where \(\nu_i\) is frequency of vibration mode \(i\) in Hz. The \(q_{\rm vib}^{(\rm BOT)}\) is the partition function calculated by taking bottom of potential energy surface as zero point, while \(q_{\rm vib}^{(V=0)}\) is the partition function calculated by viewing vibrational ground state as zero point, and it reflects thermal excitation of vibrational states. \(q_{\rm vib}^{(V=0)}\) is more commonly used since its magnitude is not so small as \(q_{\rm vib}^{(\rm BOT)}\).

Thermodynamic quantities due to vibrations can be written as sum of individual contribution of each vibration mode:

\[ \begin{aligned} & U(0)_{\rm vib} = {\rm ZPE} = \frac{1}{2}\sum_i h\nu_i \\ & U_{0\to T}^{\rm vib} \equiv U(T)_{\rm vib} - U(0)_{\rm vib} = RT\sum_i \left(\frac{h\nu_i}{kT}\right)\frac{e^{-h\nu_i/(kT)}}{1-e^{-h\nu_i/(kT)}} \\ & C_{V,\rm vib}(T) = R\sum_i \left(\frac{h\nu_i}{kT}\right)^2 \frac{e^{-h\nu_i/(kT)}}{[1-e^{-h\nu_i/(kT)}]^2} \\ & S_{\rm vib}(T) = R\sum_i \left\{\frac{h\nu_i}{kT}\frac{e^{-h\nu_i/(kT)}}{1-e^{-h\nu_i/(kT)}} - \ln[1-e^{-h\nu_i/(kT)}]\right\} \end{aligned} \]

A vibration mode with higher frequency has larger contribution to ZPE, while a mode with lower frequency has larger contribution to the other three terms.

Frequency scale factor

The error of thermodynamic quantities calculated in above way has various sources:

(1) Representation of potential energy surface is inaccurate. This is inevitable for any practical calculation level.

(2) Anharmonic effect is ignored

(3) Torsional motions of some groups or fragments should be treated as free or hindered rotor rather than (an)harmonic oscillator.

(4) Some torsional motions may be strongly coupled with others

The errors (1) and (2) are most prominent and relatively systematical, many studies found that they can be largely eliminated by applying scale factors for vibrational frequencies. Specifically, before calculating a thermodynamic quantity, the raw harmonic frequencies are multiplied by a pre-fitted factor (\(\lambda\)). The \(\lambda\) is dependent of both theoretical method and basis set, in other words, different calculation level has different optimal \(\lambda\). Many papers fitted \(\lambda\) for ZPE by taking experimental data as reference, such as J. Chem. Theory Comput., 6, 2872 (2010), while some other papers also present \(\lambda\) for \(S\) and \(U_{0\to T}^{\rm vib}\), such as J. Phys. Chem., 100, 16502 (1996). Different \(\lambda\) should be applied to frequencies prior to calculating different thermodynamic quantities, thus Shermo allows user to simultaneously specify different \(\lambda\). Usually, the \(\lambda\) for ZPE deviates 1.0 notably (except for pure DFT functionals), therefore it should always be considered if accuracy of the result is important. In contrast, the \(\lambda\) for \(S\) and \(U_{0\to T}^{\rm vib}\) is usually close to 1.0 and thus often be ignored.

To maximally eliminate the aforementioned error source (3), hindered rotor treatment should be employed, such as using the freq=hinderedrotor keyword in Gaussian, however this is beyond scope of Shermo. Using the black box quasi-RRHO treatment as shown below can also partially relieve this error.

The way to deal with error source (4) is more involved and sophisticated method is needed, such as the MS-T method, see J. Chem. Theory Comput., 9, 1356 (2013).

Quasi-RRHO treatment for low frequencies

The breakdown of harmonic approximation of low-frequency modes is well known, especially for Gibbs free energy. The \(S\) and heating contribution to \(U\) due to mode \(i\) can be written as

\[ \begin{aligned} & S^i(T) = R\left\{\frac{h\nu_i}{kT}\frac{1}{e^{h\nu_i/(kT)}-1} - \ln[1-e^{-h\nu_i/(kT)}]\right\} \\ & U_{0\to T}^i(T) = \left(\frac{h\nu_i}{kT}\right)\frac{RT}{e^{h\nu_i/(kT)}-1} \end{aligned} \]

Clearly, as the vibration frequency \(\nu_i\) approaches zero, the above two quantities diverge. There are some methods to avoid this issue and they are usually referred to as quasi-RRHO (QRRHO) or modified RRHO (mRRHO), and usually frequency scaled factor is simultaneously used, in this case they are known as modified scaled RRHO (msRRHO).

Truhlar’s QRRHO model: The first QRRHO treatment is the one suggested by Truhlar, see J. Phys. Chem. B, 115, 14556 (2011) for example. This method artificially raises all (unscaled) frequencies that lower than 100 cm⁻¹ to 100 cm⁻¹, this threshold was chosen empirically by the authors. In Shermo, this treatment is enabled in the calculation of \(S\), \(C_{\rm V}/C_{\rm P}\) and \(U/H/G_{0 \to T}\) when "ilowfreq" parameter is set to 1, and the threshold can be adjusted by "ravib" parameter. In J. Phys. Chem. A, 123, 7430 (2019), it was shown that for small molecules, the accuracy of result is insensitive to the choice of the threshold, while for molecular complexes, the threshold of 100 cm⁻¹ performs best.

Grimme’s QRRHO model: Another very popular QRRHO treatment is the one proposed by Grimme in Chem. Eur. J., 18, 9955 (2012), it only applies to entropy. Essentially, this method performs interpolation between the vibrational entropy calculated under harmonic approximation (S_RRHO) and the one of free-rotor model (S_FR):

\[ \begin{aligned} S = w(v)S_{\rm RRHO} + [1-w(v)]S_{\rm FR} \\ w(v) = \frac{1}{1+(v_0/v)^\alpha} \qquad \alpha = 4 \end{aligned} \]

where \(\nu\) is the unscaled frequency of the mode under consideration, while \(\nu_0\) is usually taken as 100 cm⁻¹ (determined by intpvib in settings.ini), other values within 50~150 cm⁻¹ can also give similar result. The S_FR is calculated as follows

\[ S_{\rm FR} = R\left(\frac{1}{2} + \ln\sqrt{\frac{8\pi^3\mu'kT}{h^2}}\right) \]

with

\[ \mu' = \frac{\mu B_{\rm av}}{\mu + B_{\rm av}} \qquad \mu = \frac{h}{8\pi^2\nu} \]

where \(\mu\) is moment of inertia of a free-rotor with the same frequency. \(\mu'\) is effective moment of inertia. The \(B_{\rm av}\) corresponds to the average molecular moment of inertia and was empirically chosen to be \(10^{-44}\) kg·m².

The Grimme's QRRHO can be enabled by setting "ilowfreq=2". It is evidently more elegant than artificially raising low frequencies and thus should be more accurate and robust. Therefore, I always recommend setting "ilowfreq" to 2 if the system is quite flexible (e.g. a large molecule with many rotatable bonds or molecular complexes), since it is always accompanied by lots of low-frequency modes. It was found that for a flexible organic molecule consisting of several hundreds of atoms, the difference between the G computed based on RRHO and QRRHO can reach as high as several kcal/mol, clearly in this case employing QRRHO is compulsory if quantitative accuracy of the result is important. For rigid molecules, the two QRRHO methods do not bring additional benefits since there is no low frequency.

Minenkov’s QRRHO model: In J. Comput. Chem., 44, 1807 (2023), Otlyotov and Minenkov showed that if the interpolation between harmonic approximation and free-rotor model is not only applied to entropy as above, but also applied to vibrational contribution to internal energy, then the free energy will be slightly more accurate, and quality of internal energy will also be improved. This consideration also makes the Grimme’s QRRHO model more consistent. This model corresponds to "ilowfreq=3". The interpolation equation for the internal energy part is as follows

\[ U_{\rm vib} = w(v)U_{\rm RRHO} + [1-w(v)]U_{\rm FR} \]

where definition of \(w\) is exactly the same as that in Grimme’s QRRHO model, \(U_{\rm FR} = RT/2\) is free-rotor contribution to internal energy, and \(U_{\rm RRHO}\) is vibrational internal energy estimated by harmonic approximation. Notice that when using this model, ZPE and \(U_{0\to T}\) cannot be separated, and thus they are not outputted by Shermo.

Note that the \(S_{\rm RRHO}\) and \(U_{\rm RRHO}\) in above formula are calculated based on scaled frequencies.

Consideration of small imaginary frequencies

Note

See my blog article “On the calculation of thermodynamic quantities with small imaginary frequencies” (http://sobereva.com/699, in Chinese) for more comprehensive discussion about this topic.

For systems such as flexible macromolecules and weakly interacting complexes, where certain dimensions of potential energy surface are very smooth, after geometric optimization and vibration analysis, it is often found that there are small imaginary frequencies, such as 9.8 and 27.3 cm⁻¹, which are generally less than 50-100 cm⁻¹. This is usually caused by the issue that geometric optimization was not converged accurately enough, or the integration grid employed is not fine enough (this problem is particularly significant for functionals that are sensitive to the integration grid such as M06-2X). In this case, the imaginary frequency modes usually correspond to the motion such as methyl rotation, flexible skeleton deformation, and intermolecular relative translation or rotation. Note that the imaginary frequency (usually above 100 cm⁻¹) that exists due to the transition state structure or the symmetry of the initial structure is higher than the actual one does not belong to the situation discussed here.

When encountering such small imaginary frequencies, the most ideal approach is to use tighter geometry optimization convergence threshold or use better integration grid to redo the optimization and vibration analysis so that the imaginary frequencies are completely eliminated before calculating thermodynamic quantities. However, for large systems, eliminating all imaginary frequencies often requires very high computational time and thus is difficult to achieve. In fact, it is not unacceptable to calculate thermodynamic quantities in the presence of such small imaginary frequencies. These frequencies actually do not cause obvious problems in the calculation accuracy of electronic energy, because gradients of potential energy surface in the direction of normal vectors of these imaginary frequencies are fairly low, so perturbating structure in these directions does not affect electronic energy notably. In addition, according to the equation of ZPE, the contribution of these small imaginary frequencies to ZPE is also negligible. Therefore, the existence of small imaginary frequencies will not cause evident inaccuracy in the calculation of electronic energy and U(0).

When calculating thermodynamic corrections to free energy, internal energy and enthalpy at finite temperature, these small imaginary frequencies need special consideration, because they have a significant contribution to entropy, heat capacity and the change in internal energy caused by heating from 0 K, so they cannot be simply ignored. In Section 3.2 of Angew. Chem. Int. Ed., 2022, e202205735, Grimme et al. suggest directly treating these small imaginary frequencies as real frequencies when using the aforementioned QRRHO model. Although this approach is obviously not strict, it is indeed reasonable to a certain extent. As shown in Figure 2 of Chem. Eur. J., 18, 9955 (2012), contributions of small frequencies to entropy are close. After eliminating the small imaginary frequencies by refining structure, one will definitely get same number of real frequencies with small wave numbers, so directly replacing the small imaginary frequencies with real frequencies of same magnitude may be an acceptable approximation, perhaps much better than completely ignoring their contributions. In addition, as shown in Figure 1 of J. Phys. Chem., 100, 16502 (1996), in the low-frequency range, the change in internal energy caused by heating from 0 K to expected temperature is insensitive to frequency, which also somewhat justified the treatment of setting the small imaginary frequencies as real ones. Note that this treatment is never suitable for the case when using standard RRHO model (ilowfreq = 0), because entropy in low-frequency range is very sensitive to frequency value in this case.

In Shermo, imaginary frequencies with magnitude smaller than “imagfreq” parameter (cm⁻¹) in settings.ini will be treated as real frequencies. For example, if it is set to 50, then 34.9i cm⁻¹ will be automatically replaced with 34.9 cm⁻¹ after loading data from input file. Obviously, if “imagfreq” is set to 0, then no imaginary frequency will be set to real.

Unfortunately, according to my brief test in http://sobereva.com/699 (a molecular complex with imaginary frequency of -6.4i cm⁻¹) with ilowfreq = 2 or 3, the treatment of setting the small imaginary frequencies to real ones does not reduce the difference in G_corr(298.15 K) between the case with and without the imaginary frequency, but increases notably the difference! So, setting small imaginary frequencies to real ones can never be regarded as a robust treatment. It should be aware that the aforementioned justification of this treatment is evidently too optimistic, actually during refining the structure from that with imaginary frequencies to exact minimum, not only the imaginary frequencies change, but also values of many real low frequencies will change, which makes situation fairly complicated. However, according to theoretical analysis and my test, setting the small imaginary frequencies to real ones is beneficial for calculating U_corr when there are small imaginary frequencies. See http://sobereva.com/699 for detailed discussion.

A.5 Electron contribution

Contribution to thermodynamic data from electrons comes from two parts, the electronic energy and electron excitation. The contribution due to the latter one is given below, where \(\varepsilon_i\) is electronic energy level \(i\) with respect to ground state, and \(g_i\) is degeneracy of level \(i\).

\[ \begin{aligned} & q_{\rm ele} = \sum_{i=0}^{\infty} g_ie^{-\varepsilon_i/(kT)} \\ & U(T)_{\rm ele} - U(0)_{\rm ele} = RT\frac{\sum_i g_i[\varepsilon_i/(kT)]e^{-\varepsilon_i/(kT)}}{\sum_i g_ie^{-\varepsilon_i/(kT)}} \\ & S_{\rm ele} = R\ln\left(\sum_i g_ie^{-\varepsilon_i/(kT)}\right) + R\frac{\sum_i g_i[\varepsilon_i/(kT)]e^{-\varepsilon_i/(kT)}}{\sum_i g_ie^{-\varepsilon_i/(kT)}} \\ & C_{V,\rm ele} = R\frac{\sum_i g_i[\varepsilon_i/(kT)]^2e^{-\varepsilon_i/(kT)}}{\sum_i g_ie^{-\varepsilon_i/(kT)}} - R\left\{\frac{\sum_i g_i[\varepsilon_i/(kT)]e^{-\varepsilon_i/(kT)}}{\sum_i g_ie^{-\varepsilon_i/(kT)}}\right\}^2 \end{aligned} \]

At room temperature, only the states with excitation energy less than a few tenths of an eV may have notable contribution to thermodynamic data, therefore, contribution of electronic excitations is usually not explicitly taken into account since the excitation energies are commonly one order of magnitude higher (except for the cases that energy splitting is caused by spin-orbit coupling, the magnitude of lowest excitation energy in this situation is usually about 0.01 eV).

If ground state is degenerate, electronic contribution to entropy will be significant and never be ignored. If the degeneracy of ground state is g₀, then the above expressions can be reduced to

\[ \begin{aligned} & q_{\rm ele} = g_0 \\ & S_{\rm ele} = R\ln g_0 \end{aligned} \]

As an instance, a doublet radical has spin multiplicity of 2, thus electronic contribution to entropy should be Rln2 = 5.76 J·(mol·K)⁻¹. Note that ground state degeneracy does not contribute to U and C_V.

A.6 Conformation weighting

Flexible molecules usually have many thermal accessible conformations (i.e.their relative free energies with respect to the lowest one are about the order of magnitude of RT, which corresponds to 2.48 kJ/mol at standard condition). Similarly, molecular complexes usually have many thermal accessible configurations. Appropriate consideration of the conformers/configurations is crucial for obtaining reliable and meaningful thermodynamic quantities for these kinds of systems.

The conformation (or configuration, similarly hereinafter) weighted U, H and C_V/C_P can be easily obtained via the following way

\[ U = \sum_i^{\rm conf} p_iU_i \qquad H = \sum_i^{\rm conf} p_iH_i \qquad C = \sum_i^{\rm conf} p_iC_i \]

where \(p_i\) stands for Boltzmann weight of conformer \(i\).

The way of calculating conformation weighted entropy is different to other thermodynamic quantities, because an additional term, the conformation entropy (\(S^{\rm conf}\)), must be taken into account

\[ \begin{aligned} S = S^{\rm conf} + \sum_i^{\rm conf} p_iS_i \\ S^{\rm conf} = -R\sum_i^{\rm conf} p_i\ln p_i \end{aligned} \]

The \(S^{\rm conf}\) comes from the interconversion between various conformations.

The conformation weighted Gibbs free energy can be obtained in usual way based on conformation weighted H and S, namely G=H-TS.

The Boltzmann weights can be calculated straightforwardly via relative Gibbs free energy

\[ \begin{aligned} p_i = \frac{e^{-\Delta G_i/RT}}{\sum_j e^{-\Delta G_j/RT}} \times 100\% \\ \Delta G_i = G_i - G_{\rm lowest} \end{aligned} \]

where G_lowest is the lowest Gibbs free energy of all considered conformers.