Thermochemical properties of molecules#
index:thermochemistry mod:ase.thermochemistry can be used to estimate thermodynamic properties of gases in the ideal gas limit. The module needs as input the geometry, the total energy, the vibrational energies, and some information about the molecular symmetry. We first consider an N2 molecule.
The symmetry numbers are determined by the molecular point group cite:springerlink-10.1007/s00214-007-0328-0. Here is a table of the most common ones.
point group |
\(\sigma\) |
examples |
|---|---|---|
C1 |
1 |
|
Cs |
1 |
|
C2 |
2 |
|
C2v |
2 |
H2O |
C3v |
3 |
NH3 |
\(C_{\infty v}\) |
1 |
CO |
D2h |
4 |
|
D3h |
6 |
|
D5h |
10 |
|
\(D_{\infty h}\) |
2 |
CO2, H2 |
D3d |
6 |
|
Td |
12 |
CH4 |
Oh |
24 |
from ase.build import molecule
from ase.thermochemistry import IdealGasThermo
from vasp import Vasp
atoms = molecule('N2')
atoms.set_cell((10,10,10), scale_atoms=False)
# first we relax a molecule
calc = Vasp(label='molecules/n2-relax',
xc='PBE',
encut=300,
ibrion=2,
nsw=5,
atoms=atoms)
electronicenergy = atoms.get_potential_energy()
# next, we get vibrational modes
calc2 = Vasp(label='molecules/n2-vib',
xc='PBE',
encut=300,
ibrion=6,
nfree=2,
potim=0.15,
nsw=1,
atoms=atoms)
calc2.calculate()
vib_freq = calc2.get_vibrational_frequencies() # in cm^1
#convert wavenumbers to energy
h = 4.1356675e-15 # eV*s
c = 3.0e10 #cm/s
vib_energies = [h*c*nu for nu in vib_freq]
print('vibrational energies\n====================')
for i,e in enumerate(vib_energies):
print('{0:02d}: {1} eV'.format(i,e))
# # now we can get some properties. Note we only need one vibrational
# energy since there is only one mode. This example does not work if
# you give all the energies because one energy is zero.
thermo = IdealGasThermo(vib_energies=vib_energies[0:0],
potentialenergy=electronicenergy, atoms=atoms,
geometry='linear', symmetrynumber=2, spin=0)
# temperature in K, pressure in Pa, G in eV
G = thermo.get_gibbs_energy(temperature=298.15, pressure=101325.)
/opt/conda/lib/python3.9/site-packages/ase/structure.py:5: UserWarning: Moved to ase.build
warnings.warn('Moved to ase.build')
vibrational energies
====================
00: 0.29244330802722723 eV
01: 0.0153405692581911 eV
02: 0.0153405692581911 eV
03: 2.8536105750000002e-09 eV
04: 1.116630225e-09 eV
05: -1.24070025e-10 eV
Enthalpy components at T = 298.15 K:
===============================
E_pot -16.524 eV
E_ZPE 0.000 eV
Cv_trans (0->T) 0.039 eV
Cv_rot (0->T) 0.026 eV
Cv_vib (0->T) 0.000 eV
(C_v -> C_p) 0.026 eV
-------------------------------
H -16.434 eV
===============================
Entropy components at T = 298.15 K and P = 101325.0 Pa:
=================================================
S T*S
S_trans (1 bar) 0.0015590 eV/K 0.465 eV
S_rot 0.0007870 eV/K 0.235 eV
S_elec 0.0000000 eV/K 0.000 eV
S_vib 0.0000000 eV/K 0.000 eV
S (1 bar -> P) -0.0000011 eV/K -0.000 eV
-------------------------------------------------
S 0.0023449 eV/K 0.699 eV
=================================================
Free energy components at T = 298.15 K and P = 101325.0 Pa:
=======================
H -16.434 eV
-T*S -0.699 eV
-----------------------
G -17.133 eV
=======================
Let us compare this to what is in the Nist webbook via the Shomate equations.
import numpy as np
A = 28.98641
B = 1.853978
C = -9.647459
D = 16.63537
E = 0.000117
F = -8.671914
G = 226.4168
H = 0.0
T = 298.15
t = T/1000.
S = A*np.log(t) + B*t + C*t**2/2 + D*t**3/3 - E/(2*t**2) + G
print('-T*S = {0:1.3f} eV'.format(-T*S/1000/96.4853))
-T*S = -0.592 eV
This is reasonable agreement for the entropy. You will get different results if you use different exchange correlation functionals.