常见分子吉布斯自由能计算
另一篇文章讲表面吸附分子的振动自由能计算，见叠加态：表面吸附分子的振动自由能计算。表面吸附分子和理想气体分子的振动自由能计算都可以用vaspkit实现，对应task为501和502.
分子自由能的实验值可参考叠加态：几种气体热力学量
电化学理论计算需要求分子体系的吉布斯自由能，其各项的计算如图一所示。其中求旋转熵S_rot涉及到分子对称数σ，常见分子的对称性可在wikipedia查询如图二，对称性对应的对称数见图三。具体计算可使用python的ASE库，代码示例附在后面。代码中的分子能量及振动频率由第一性原理计算得到。
原理及公式
图一：https://wiki.fysik.dtu.dk/ase/ase/thermochemistry/thermochemistry.html
分子对称性
https://en.wikipedia.org/wiki/Molecular_symmetry
查表得到对称性对应的转动对称数
C.J. Cramer. Essentials of Computational Chemistry, Second Edition. Wiley, 2004.
利用pymatgen直接给出分子的对称数
from ase.io import read,write
import numpy as np
from pymatgen.symmetry.analyzer import PointGroupAnalyzer
from pymatgen.io.ase import AseAtomsAdaptor
trans=AseAtomsAdaptor()
atoms=read("POSCAR")
mol=trans.get_molecule(atoms)
pga=PointGroupAnalyzer(mol)
sym_N=pga.get_rotational_symmetry_number()
#print(E,vib_e,sym_N)ASE计算气体自由能(基于理想气体近似)
from ase.build import molecule
from ase.calculators.emt import EMT
from ase.optimize import QuasiNewton
from ase.vibrations import Vibrations
from ase.thermochemistry import IdealGasThermo
import numpy as np
H2 = molecule('H2')
CO = molecule('CO')
CO2= molecule('CO2')
H2O= molecule('H2O')
E_H2 = -6.6876
E_CO = -14.792
E_CO2 = -22.9776
E_H2O = -14.1996
vib_H2 = np.array([0.5305489])
vib_CO = np.array([0.2636537])
vib_CO2 = np.array([0.293035,0.16338,0.07855,0.07827])
vib_H2O = np.array([0.449186,0.432445,0.197995])
T_H2 = IdealGasThermo(vib_energies=vib_H2,
potentialenergy=E_H2,
atoms=H2,
geometry='linear',
symmetrynumber=2, spin=0)
G_H2 = T_H2.get_gibbs_energy(temperature=298.15, pressure=101325.)
T_CO = IdealGasThermo(vib_energies=vib_CO,
potentialenergy=E_CO,
atoms=CO,
geometry='linear',
symmetrynumber=1, spin=0)
G_CO = T_CO.get_gibbs_energy(temperature=298.15, pressure=101325.)
T_CO2 = IdealGasThermo(vib_energies=vib_CO2,
potentialenergy=E_CO2,
atoms=CO2,
geometry='linear',
symmetrynumber=2, spin=0)
G_CO2 = T_CO2.get_gibbs_energy(temperature=298.15, pressure=101325.)
T_H2O = IdealGasThermo(vib_energies=vib_H2O,
potentialenergy=E_H2O,
atoms=H2O,
geometry='nonlinear',
symmetrynumber=2, spin=0)
pressure_H2O=101325.*0.035
G_H2O = T_H2O.get_gibbs_energy(temperature=298.15, pressure=pressure_H2O)Table A: DFT energies, Zero point energy, enthalpic temperature correction and entropy correction at 25 C. The entropies of H2(g), CO(g) and CO2(g) are calculated at 1 atm, while the entropy of H2O(g=l) is calculated at 0.035 atm, which corresponds to the vapor pressure of liquid water. Units are in eV.
| EDFT | EZPE | 积分项 | -T*S | G | |
| H2 | -6.688 | 0.265 | 0.091 | -0.402 | -6.735 |
| CO | -14.792 | 0.132 | 0.091 | -0.611 | -15.182 |
| CO2 | -22.978 | 0.307 | 0.099 | -0.662 | -23.235 |
| H2O | -14.200 | 0.540 | 0.104 | -0.670 | -14.227 |
