Thermodynamic functions calculator

  Chemcraft can compute the entropy and the Gibbs energy of a molecule by its molecular parameters, including the geometry (moments of inertia) and vibrational frequencies. This computation is made in "rigid rotor - harmonic oscillator" approximation (RRHO). The main aim of this utility is the possibility to cut low frequencies, which produce big errors because of ahharmonicity. You can specify a threshold value, e.g. 100 cm**-1, and the values of the frequencies below this threshold will be increased to this threshold value. The formulas for the computation were taken from the following refs.: 


Donald A. McQuarrie, John D. Simon. Molecular Thermodynamics. University Science Books, 1999. 


Thermochemistry in Gaussian. Joseph W. Ochterski, Ph.D. June 2, 2000. https://gaussian.com/thermo/

ASE thermochemistry computation. https://wiki.fysik.dtu.dk/ase/ase/thermochemistry/thermochemistry.html 


  Here are the formulas implemented in Chemcraft:

 


 
  Other functions:

 

  You can download the Delphi source code which computes the thermodynamic functions in Chemcraft here: chemdata_u_b646.zip. 

Troubleshooting 
  Currently this utility in Chemcraft does not always determine correctly the point groups of highly-symmetrical molecules (point groups T, Td, Th, O, Oh, I, Ih), so the Srot values can be computed improperly correctly for such molecules. If your molecule is symmetrical, you should check the point group in the printout. Before applying the frequencies threshold, we recommend to compare the S, C, G values printed by Chemcraft (without using the frequencies threshold) with the values printed by your QC package. We have performed such a comparison for a number of output files produced by different packages. For Gaussian, the results for Gibbs energies and their components are almost equal, except some molecules with point groups T, Td, Th, O, Oh, I, Ih. For Orca, the values of most thermodynamic functions are almost equal, but the Svib values are slightly different; this is caused by the fact that the version of Orca we used (5.0.4) uses the Svib computation algorithm other than the standard RRHO (in Orca, the QRRHO approach is implemented). For GAMESS-US (two files were checked), a good match was found all components of energies (TRANS., ROT., VIB., TOTAL), heat capacities (all too), entropies (all too). For Priroda package, a rather good match was obtained for most components except rotational S, because we were able to check only a file for symmetrical molecule. For Psi, a good match was obtained for all components of G except the electronic S printed in file. For QChem files we checked, a good match was obtained for all enthalpy and entropy components. For Spartan files, a good match was found for all enthalpy and entropy components.

 

  Experimental verification of the Gibbs energies computation


  We have performed a computation of pKa values (pKa=∆Ga/2.303RT) of some carboxylic acids and phenols, which can be compared with experimental ones from literature. The set of compounds studies included 21 carbon and benzoic acids, and 14 phenols. The pKa values do not directly match the experimental ones from literature, but the correlations are clearly seen (as usual for the quantum chemistry). The molecules were computed via Orca 5.0.4 package with method PCM B3LYP/6-31++G(D,P), with explicit water molecules. 

  Our models contain explicit water molecules which produce very small frequencies (like 10 cm**-1), and we found that the anharmonic origin of these frequencies produces very big errors.

  The mean absolute deviations (mean differences between the experimental pKa values and the values calculated by linear fit parameters) are in the range of 0.08 - 0.17 pKa units. For the computations of Gibbs energies and pKa values, the MADs for the carboxylic acids obtained are 0.170 pKa units for Orca QRRHO computations without the low frequency threshold, 0.126 pKa units for Orca QRRHO computation with 100 cm**-1 threshold, 0.214 pKa units for Chemcraft RRHO computation without the threshold, 0.149 for Chemcraft RRHO computation with the threshold 40 cm**-1, 0.144 with 100 cm**-1 threshold, 0.124 with 300 cm**-1 threshold, and 0.126 with 800 cm**-1 threshold. For RRHO computations of phenols, the MADs are 0.104 pKa units for Orca QRRHO computations without the low frequency threshold, 0.111 pKa units for Orca QRRHO computation with 100 cm**-1 threshold, 0.161 pKa units for Chemcraft RRHO computation without the threshold, 0.118 for Chemcraft RRHO computation with the threshold 40 cm**-1, 0.090 with 100 cm**-1 threshold, 0.082 with 300 cm**-1 threshold, and 0.081 with 800 cm**-1 threshold. So, our study confirms that applying the threshold of 100-800 cm**-1 is a good approach which improves the agreement of computed Gibbs energies with the experiment. 

  At these pictures you can see the comparison of correlations between the computed and experimental pKa values obtained in our study with different approaches:

  We have published our results in Molecular Physics (DOI 10.1080/00268976.2024.2396535):

https://simplecompchem.blogspot.com/2024/09/qspr-prediction-of-acidities-of.html

  Some people can find somewhat strange, that the best agreement with experiment in our study was achieved with 300 cm**-1 and 800 cm**-1 thresholds, while it is usually recommended to use the value of 100 cm**-1 threshold. Possibly this is a statistical deviation caused by the fact that the number of compounds is not big enough, and, more importantly, there are much more significant sources of errors here (the main source of errors in an incomplete taking into account of the specific solvatation). We found that with the 800 cm-1 thresholds, the agreement with the experiment (MADS for linear fits) is almost the same as if no entropy computation is performed (when the Kohn-Sham energies of deprotonation are correlated with the pKa values, not the computed pKa values). This means that the role of entropy computations is almost neglible, and still if these computations are performed without the thresholds - an additional error appears due to the neglect of the anharmonicity of small modes.

 

Citation

We suppose that a proper citation can be like this:

"The Gibbs energies were computed with Chemcraft program [1], in which the commonly used formulas for RRHO approach [2] have been implemented. An approach of raising low vibrational frequencies to a specified threshold value, firstly suggested by Truhlar et all [3], was used."

[1] See the Citation page

[2] Donald A. McQuarrie, John D. Simon. Molecular Thermodynamics. University Science Books, 1999. 

[3] Ribeiro, R. F.; Marenich, A. V.; Cramer, C. J. & Truhlar, D. G. "Use of Solution-Phase Vibrational Frequencies in Continuum Models for the Free Energy of Solvation”. J. Phys. Chem. B 2011, 115, 49, 14556–14562

BACK