M3.D4: Thermodynamic potentials for spontaneity and other thermodynamic properties#
Learning Objectives#
Calculate thermodynamic potentials using PySCF of several diatomic ideal gases
Ensembles and variables and Thermodynamic potentials for the 2nd law#
Thermodynamic functions can be calculated as a function of its partition function. We will just summarize the most relevant thermodynamic potentials down here.
Extensive property |
Corresponding intensive property |
|---|---|
(depends on the size of the system) |
(it can be measured at any point of the system) |
Volume (V) |
Pressure (p) |
Energy (E) |
Temperature (T) |
Moles (N) |
Chemical potential (\(\mu\)) |
Ensembles |
Potential measuring spontaneity |
|---|---|
Microcanonical (NVE) |
Entropy \(\Delta S_{N,V,E} > 0\) |
Canonical (NVT) |
Helmholtz Free Energy \(\Delta A_{N,V,T} < 0\) |
Isothermal-Isobaric (NPT) |
Gibbs Free Energy \(\Delta G_{N,P,T} < 0\) |
In the microcanonical ensemble, meaning, energy, volume and number of particles being constant, the entropy has a simple expression that we have seen before.
Where W implies the number of states (the partition function if all states were equally probable).
However, when instead of constant energy we keep temperature constant, that is, considering the NVT ensemble, keeping temperature, volume and number of particle constant not all states “j” are equally probable, so the average would be:
and using now the canonical partition function instead of \(p_j\) we get a potential that is not as elegant.
with \(\Delta S_{N,V,E}>0\) as the criteria for spontaneity
In the canonical ensemble, it is the free energy (Helmholtz free energy) that has a simple expression
with \(\Delta A_{V,T}<0\) as the criteria for spontaneity
In the so-called isothermic isobaric ensemble where the variables held constant are N,P,T the Gibbs free energy is the thermodynamic potential with the simple expression
with \(\Delta G_{P,T}<0\) as the criteria for spontaneity
Before we proceed with more mathematical expressions, it may be useful to look at real calculations of molecules and see what contributions from vibrational, translational, rotational, and electronic we should expect
Calculating thermodynamic potentials using pyscf#
In a separate notebook, for each of the following molecules: CO, N2, O2(use singlet!), F2, HF, HCl, CN(-), OH(-), H2. See below an example for CO.
Optimize the structure
Run a frequency analysis of the optimized structure. Save the output each in a separate file
Build a function that for a given output file it reports a list of thermodynamic data such as the electronic, translational, rotational, and vibrational contributions to thermal energy, enthalpy, and free energy.
Build a table with those data.
Identify the most important contributions: compare what contributes the most to the internal/thermal energy, is it translational, rotational, and vibrational contributions?
Compare among molecules.
#pip install geometric
from pyscf import gto
from pyscf.geomopt.geometric_solver import optimize
from pyscf.hessian import thermo
from pyscf import lib
# Convert from Bohr to Angstroms
bohr_to_angstrom = lib.parameters.BOHR
molCO = gto.M(
atom = 'O 0 0 0; C 0 0 1.3', # in Angstrom
basis = '631g',
symmetry = True,
charge=0,
spin=0,
verbose=3,
output='m3d2_co_fq.out'
)
# Perform the Hartree-Fock calculation
co_sp = molCO.RHF.run()
co_opt = optimize(co_sp).run()
print("Optimized CO structure:")
print(co_opt.atom_coords()*bohr_to_angstrom)
hessian = co_sp.Hessian().kernel()
# Frequency analysis
freq_info = thermo.harmonic_analysis(co_sp.mol, hessian)
# Thermochemistry analysis at 298.15 K and 1 atmospheric pressure
thermo_info = thermo.thermo(co_sp, freq_info['freq_au'], 298.15, 101325)
print('Rotation constant')
print(thermo_info['rot_const'])
print('Zero-point energy')
print(thermo_info['ZPE' ])
print('Internal energy at 0 K')
print(thermo_info['E_0K' ])
print('Internal energy at 298.15 K')
print(thermo_info['E_tot' ])
print('Enthalpy energy at 298.15 K')
print(thermo_info['H_tot' ])
print('Gibbs free energy at 298.15 K')
print(thermo_info['G_tot' ])
print('Heat capacity at 298.15 K')
print(thermo_info['Cv_tot'])
output file: m3d2_co_fq.out
geometric-optimize called with the following command line:
/Users/pratr001/miniconda3/envs/pchem/lib/python3.12/site-packages/ipykernel_launcher.py --f=/Users/pratr001/Library/Jupyter/runtime/kernel-v3248b0861d0d945a8ae00f635235d9ec3a3ebeebc.json
())))))))))))))))/
())))))))))))))))))))))))),
*)))))))))))))))))))))))))))))))))
#, ()))))))))/ .)))))))))),
#%%%%, ()))))) .))))))))*
*%%%%%%, )) .. ,))))))).
*%%%%%%, ***************/. .)))))))
#%%/ (%%%%%%, /*********************. )))))))
.%%%%%%# *%%%%%%, *******/, **********, .))))))
.%%%%%%/ *%%%%%%, ** ******** .))))))
## .%%%%%%/ (%%%%%%, ,****** /)))))
%%%%%% .%%%%%%# *%%%%%%, ,/////. ****** ))))))
#% %% .%%%%%%/ *%%%%%%, ////////, *****/ ,)))))
#%% %%% %%%# .%%%%%%/ (%%%%%%, ///////. /***** ))))).
#%%%%. %%%%%# /%%%%%%* #%%%%%% /////) ****** ))))),
#%%%%##% %%%# .%%%%%%/ (%%%%%%, ///////. /***** ))))).
## %%% .%%%%%%/ *%%%%%%, ////////. *****/ ,)))))
#%%%%# /%%%%%%/ (%%%%%% /)/)// ****** ))))))
## .%%%%%%/ (%%%%%%, ******* ))))))
.%%%%%%/ *%%%%%%, **. /******* .))))))
*%%%%%%/ (%%%%%% ********/*..,*/********* *))))))
#%%/ (%%%%%%, *********************/ )))))))
*%%%%%%, ,**************/ ,))))))/
(%%%%%% () ))))))))
#%%%%, ()))))) ,)))))))),
#, ()))))))))) ,)))))))))).
()))))))))))))))))))))))))))))))/
())))))))))))))))))))))))).
())))))))))))))),
-=# geomeTRIC started. Version: 1.1 #=-
Current date and time: 2025-04-17 09:39:39
#========================================================#
#| Arguments passed to driver run_optimizer(): |#
#========================================================#
customengine <pyscf.geomopt.geometric_solver.PySCFEngine object at 0x179ce6570>
input /var/folders/1b/pw4sp5tj3_16kjl12vvf8vt40000gq/T/tmpxbj0cqve/f3510fde-899e-48dd-b6a4-3c5e95daf0f8
logIni /Users/pratr001/miniconda3/envs/pchem/lib/python3.12/site-packages/pyscf/geomopt/log.ini
----------------------------------------------------------
Custom engine selected.
Bonds will be generated from interatomic distances less than 1.20 times sum of covalent radii
6 internal coordinates being used (instead of 6 Cartesians)
Internal coordinate system (atoms numbered from 1):
Distance 1-2
Translation-X 1-2
Translation-Y 1-2
Translation-Z 1-2
Rotation-A 1-2
Rotation-B 1-2
Rotation-C 1-2
<class 'geometric.internal.Distance'> : 1
<class 'geometric.internal.TranslationX'> : 1
<class 'geometric.internal.TranslationY'> : 1
<class 'geometric.internal.TranslationZ'> : 1
<class 'geometric.internal.RotationA'> : 1
<class 'geometric.internal.RotationB'> : 1
<class 'geometric.internal.RotationC'> : 1
> ===== Optimization Info: ====
> Job type: Energy minimization
> Maximum number of optimization cycles: 300
> Initial / maximum trust radius (Angstrom): 0.100 / 0.300
> Convergence Criteria:
> Will converge when all 5 criteria are reached:
> |Delta-E| < 1.00e-06
> RMS-Grad < 3.00e-04
> Max-Grad < 4.50e-04
> RMS-Disp < 1.20e-03
> Max-Disp < 1.80e-03
> === End Optimization Info ===
Step 0 : Gradient = 2.498e-01/2.498e-01 (rms/max) Energy = -112.6189499615
Hessian Eigenvalues: 5.00000e-02 5.00000e-02 5.00000e-02 ... 5.00000e-02 5.00000e-02 6.71931e-01
Step 1 : Displace = 9.837e-02/9.837e-02 (rms/max) Trust = 1.000e-01 (=) Grad = 7.714e-02/7.714e-02 (rms/max) E (change) = -112.6652860802 (-4.634e-02) Quality = 0.998
Hessian Eigenvalues: 5.00000e-02 5.00000e-02 5.00000e-02 ... 5.00000e-02 5.00000e-02 8.79414e-01
Step 2 : Displace = 2.321e-02/2.321e-02 (rms/max) Trust = 1.414e-01 (+) Grad = 4.568e-02/4.568e-02 (rms/max) E (change) = -112.6663834887 (-1.097e-03) Quality = 0.324
Hessian Eigenvalues: 5.00000e-02 5.00000e-02 5.00000e-02 ... 5.00000e-02 5.00000e-02 1.40021e+00
Step 3 : Displace = 8.632e-03/8.632e-03 (rms/max) Trust = 1.414e-01 (=) Grad = 4.454e-03/4.454e-03 (rms/max) E (change) = -112.6672145698 (-8.311e-04) Quality = 1.115
Hessian Eigenvalues: 5.00000e-02 5.00000e-02 5.00000e-02 ... 5.00000e-02 5.00000e-02 1.26368e+00
Step 4 : Displace = 9.326e-04/9.326e-04 (rms/max) Trust = 2.000e-01 (+) Grad = 3.002e-04/3.002e-04 (rms/max) E (change) = -112.6672219101 (-7.340e-06) Quality = 0.935
Hessian Eigenvalues: 5.00000e-02 5.00000e-02 5.00000e-02 ... 5.00000e-02 5.00000e-02 1.34887e+00
Step 5 : Displace = 5.889e-05/5.889e-05 (rms/max) Trust = 2.828e-01 (+) Grad = 3.235e-06/3.235e-06 (rms/max) E (change) = -112.6672219430 (-3.293e-08) Quality = 0.985
Hessian Eigenvalues: 5.00000e-02 5.00000e-02 5.00000e-02 ... 5.00000e-02 5.00000e-02 1.34887e+00
Converged! =D
#==========================================================================#
#| If this code has benefited your research, please support us by citing: |#
#| |#
#| Wang, L.-P.; Song, C.C. (2016) "Geometry optimization made simple with |#
#| translation and rotation coordinates", J. Chem, Phys. 144, 214108. |#
#| http://dx.doi.org/10.1063/1.4952956 |#
#==========================================================================#
Time elapsed since start of run_optimizer: 0.757 seconds
Optimized CO structure:
[[0. 0. 0.08466366]
[0. 0. 1.21533634]]
Rotation constant
(array([ inf, 43.58846729, 43.58846729]), 'GHz')
Zero-point energy
(np.float64(0.002706164641888301), 'Eh')
Internal energy at 0 K
(np.float64(-112.6162437968851), 'Eh')
Internal energy at 298.15 K
(np.float64(-112.61386574425926), 'Eh')
Enthalpy energy at 298.15 K
(np.float64(-112.61292155970314), 'Eh')
Gibbs free energy at 298.15 K
(np.float64(-112.63563803079069), 'Eh')
Heat capacity at 298.15 K
(np.float64(8.25633742976514e-06), 'Eh/K')
Questions#
For the Canonical NVT ensemble use PySCF to build a table calculating the entropy, the enthalpy, and the free energy at room temperature of the molecules you already calculated in the previous M3.D3 notebook.