.. _thermo:
Thermodynamics
==============
If you have successfully done a :ref:`vibfreq` calculation, now you have access to the thermodynamic functions such as enthalpy (:math:`H`), entropy (:math:`S`) and the Gibbs free energy (:math:`G`)!
Let's follow an example and use this to predict the relative energy ordering of some `butadiene isomers `_.
.. image:: thermo/butene.png
:align: center
:width: 450
using the following keywords and the geminal isomer as an example::
!B3LYP DEF2-TZVP OPT FREQ
* XYZFILE 0 1 gem-butene.xyz
Thermodynamic functions
-----------------------
Assuming you had !FREQ on the input, immediately after the frequencies and normal modes, one will see that the thermochemistry section starts::
--------------------------
THERMOCHEMISTRY AT 298.15K
--------------------------
Temperature ... 298.15 K
Pressure ... 1.00 atm
Total Mass ... 56.11 AMU
.. note::
These properties will be computed for a given temperature and pressure, assuming a gas-phase molecule and using statistical mechanics. For thermodynamics in solution, check the :ref:`solvation` section.
Inner energy
............
The first function is the inner energy (:math:`U`), that is computed following the output::
------------
INNER ENERGY
------------
The inner energy is: U= E(el) + E(ZPE) + E(vib) + E(rot) + E(trans)
E(el) - is the total energy from the electronic structure calculation
= E(kin-el) + E(nuc-el) + E(el-el) + E(nuc-nuc)
E(ZPE) - the the zero temperature vibrational energy from the frequency calculation
E(vib) - the the finite temperature correction to E(ZPE) due to population
of excited vibrational states
E(rot) - is the rotational thermal energy
E(trans)- is the translational thermal energy
Summary of contributions to the inner energy U:
Electronic energy ... -157.17473665 Eh
Zero point energy ... 0.10741046 Eh 67.40 kcal/mol
Thermal vibrational correction ... 0.00247510 Eh 1.55 kcal/mol
Thermal rotational correction ... 0.00141627 Eh 0.89 kcal/mol
Thermal translational correction ... 0.00141627 Eh 0.89 kcal/mol
-----------------------------------------------------------------------
Total thermal energy -157.06201854 Eh
The *Electronic energy* is what you get from the `FINAL SINGLE POINT ENERGY`, then the *Zero-Point Energy (ZPE)* is calculated, together with the *thermal correction*, which has to do with the rotational, vibrational and translational degrees of freedom.
The *Total correction* is amount the of energy you have to add to your *FINAL SINGLE POINT* electronic energy to convert :math:`E_{el}` into :math:`U`.
Enthalpy
........
Then the calculation of the enthalpy. That is simply :math:`U + k_BT`::
--------
ENTHALPY
--------
The enthalpy is H = U + kB*T
kB is Boltzmann's constant
Total free energy ... -157.06201854 Eh
Thermal Enthalpy correction ... 0.00094421 Eh 0.59 kcal/mol
-----------------------------------------------------------------------
Total Enthalpy ... -157.06107433 Eh
Again, the *Thermal Enthalpy correction* is what one has to add to :math:`U` to compute :math:`H`.
Entropy
.......
The entropy term follows::
-------
ENTROPY
-------
The entropy contributions are T*S = T*(S(el)+S(vib)+S(rot)+S(trans))
S(el) - electronic entropy
S(vib) - vibrational entropy
S(rot) - rotational entropy
S(trans)- translational entropy
The entropies will be listed as mutliplied by the temperature to get
units of energy
Electronic entropy ... 0.00000000 Eh 0.00 kcal/mol
Vibrational entropy ... 0.00394490 Eh 2.48 kcal/mol
Rotational entropy ... 0.01155172 Eh 7.25 kcal/mol
Translational entropy ... 0.01805279 Eh 11.33 kcal/mol
-----------------------------------------------------------------------
Final entropy term ... 0.03354941 Eh 21.05 kcal/mol
with its different components and a *Final entropy term* that is already multiplied by the temperature, the :math:`TS` term that is latter used for the free energy.
.. _gibbs:
Gibbs free energy
.................
Finally, the free energy term is printed::
-------------------
GIBBS FREE ENERGY
-------------------
The Gibbs free energy is G = H - T*S
Total enthalpy ... -157.06107433 Eh
Total entropy correction ... -0.03354941 Eh -21.05 kcal/mol
-----------------------------------------------------------------------
Final Gibbs free energy ... -157.09462374 Eh
For completeness - the Gibbs free energy minus the electronic energy
G-E(el) ... 0.08011291 Eh 50.27 kcal/mol
The last term, *G-E(el)* is exactly what you have to add to your electronic energy to transform :math:`E_{el}` into :math:`G^o` at **a given temperature and pressure**.
.. note::
This *G-E(el)* only depends on the geometry and frequencies, and thus can be computed in a lower level method, such as DFT, and later combined with the :math:`E_{el}` of a higher level method, such as DLPNO-CCSD(T) to achieve a very good prediction of the :math:`G^o`!
Comparison to experiment
------------------------
By taking the Gibbs free energy values directly from the DFT results::
!B3LYP DEF2-TZVP OPT FREQ
* XYZFILE 0 1 butene.xyz
the predicted energy differences are :math:`1.02` :math:`kJmol^{-1}` and :math:`5.99` :math:`kJmol^{-1}`, which captures the differences between cis and trans, but is certainly not good with respect to the most stable geminal isomer.
This can be made better by computing the electronic energy :math:`E_{el}` at a high-quality DLPNO-CCSD(T) level, given by the `FINAL SINGLE POINT ENERGY`::
!DLPNO-CCSD(T) aug-cc-pVTZ AUTOAUX RIJCOSX
* XYZFILE 0 1 butene_optimized_DFT.xyz
and adding the *G-E(el)* term from DFT to form a better final Gibbs free energy, and the results then are:
.. table:: Comparison of calculated versus experimental energies for butene isomers in kJ/mol
:align: center
=============== ======== =========== ======
Compound DFT CCSD(T)+DFT Exp.
=============== ======== =========== ======
1 0 0 0
2 1.02 4.45 5.3
3 5.99 9.11 10.0
=============== ======== =========== ======
Thermodynamics at different temperatures
----------------------------------------
To compute these thermodynamic properties at different temperatures, or at several, one can use::
%FREQ TEMP 77, 298, 330, 450 END
and the thermodynamic functions will be separated by the headers::
--------------------------
THERMOCHEMISTRY AT 77K
--------------------------
[...]
--------------------------
THERMOCHEMISTRY AT 298K
--------------------------
[...]
and so on.
Changing the temperature after the calculation
..............................................
If you want to compute thermodynamic functions **after** you calculation has already finished, that can be done setting a simple input such as::
!PRINTTHERMOCHEM
%GEOM
INHESSNAME "basename.hess"
END
%FREQ TEMP 290, 295, 300
END
* XYZFILE 0 1 geometry.xyz
as one can save a lot of time from recalculations.
.. note::
One could use this to also include isotopic effects, by changing the masses in the `basename.hess` file and recomputing the thermodynamics (see :ref:`KIE`).
Starting structures
-------------------
gem-butene
..........
::
C -1.43336 0.99760 -0.11833
C -0.17392 0.62436 -0.39435
C 0.14450 -0.70205 -1.02166
C 1.00254 1.50776 -0.09468
H 0.65348 -0.55631 -1.97995
H -0.75651 -1.29570 -1.20886
H 0.79819 -1.28546 -0.36521
H -2.27870 0.35250 -0.33839
H -1.65494 1.95886 0.33562
H 1.55177 1.73378 -1.01439
H 0.69898 2.45896 0.35516
H 1.68283 1.00994 0.60389
trans-butene
............
::
C 1.30292 0.51517 0.05694
C -0.04707 -0.10202 0.22754
H 1.76922 0.66709 1.03514
H 1.25468 1.48282 -0.45237
H 1.94420 -0.14853 -0.53097
C -1.18840 0.46276 -0.19477
C -2.53838 -0.15441 -0.02414
H -3.00427 -0.30734 -1.00238
H -3.17998 0.50979 0.56286
H -2.49026 -1.12154 0.48617
H -0.07366 -1.06762 0.72803
H -1.16180 1.42834 -0.69529
cis-butene
..........
::
C 1.48629 -0.48153 0.16062
C 0.03276 -0.80512 0.28503
H 1.68200 0.56889 -0.06324
H 1.92979 -1.08298 -0.63914
H 1.99954 -0.72748 1.09552
C -1.01110 0.03009 0.16420
H -0.18277 -1.85160 0.49931
C -0.98494 1.49624 -0.12334
H -2.00868 -0.39072 0.28765
H -1.47672 2.03930 0.68986
H 0.02271 1.90131 -0.23449
H -1.53322 1.69981 -1.04849