thermo

esta.general.thermo

this is thermo class for calculating the energy contribution coming from the translational, rotational , vibrational and electronic degrees of freedom towards entropy, internal energy, andspecific heat

In Nutshell: we need three components

enthalpy component H = E + PV

entropy components S

free energy components G = H -TS = E + PV - TS

Model under which the thermodynamical quantities can be calculed:

1. Ideal gas approximation

2. Harmonic approx

3. Hindered rotor and vibrator

4. crystalline materials case

Notes

For More information related to theory, see the following:

see gaussian thermo pdf file,

or chapter 10: essential of compuational chemisty,

or html file in the present folder

or 'A Note on Thermochemistry' by openmopac manual

or Introduction to computational chemistry by Jensen chapter 14.5

author: sk, email: sonukumar.physics@gmail.com

ThermoDyn

thermodyn class to find energies (eV or other units) for entropy, internal energy, and specific heat with respect to translational, rotational, vibrational, and electronic degrees of freedom

Note theat all energies are in eV or SI units

for theory: see gaussian thermo pdf file, or chapter 10: essential of compuational chemisty, or html file in the present folder

__init__(E_dft=None, vib_freq=None, E_zpe=None, S_vib=None, E_vib=None, temp=None, type_molecule=None, spin=None, symbols=None, pos=None, symmetry_number=None, pressure=None)

Parameters:

• E_dft (scalar, default: None ) –

  dft energy = electronic energy = E_ele in eV

• vib_freq (array, default: None ) –

  array of vibrational freq in eV

• E_zpe (scalar, default: None ) –

  zero point energy in eV

• S_vib (scalar, default: None ) –

  vibrational entropy in eV/K

• E_vib (scalar, default: None ) –

  internal (thermal) energy energy contribution due to vibrational motion in eV -- note when multipying exp factors etc... sk

• temp (scalar, default: None ) –

  temp in Kelvin

• type_molecule (str, default: None ) –

  linear or non-linear

Returns:

• it returns various terms present in the free energy and free energy –

get_zpe()

get quantum mechanical zpe energy in eV

get_S_vib()

get S_vib, vibrational entropy in eV/K

get_E_vib()

get vibration contribution to internal energy:

it contains vibtaration heat capacity contains 3N−6 degrees of freedom for nonlinear molecules and 3N−5 degrees of freedom for linear molecules (where N is the number of atoms)

get_moment_inertia()

get moment of inertia of the molecule or other possible system

formula: see A note on thermochemistry or else where!!

Notes

Taken from A note on thermochemisty:

IA = Σimi(RAi)2, where i runs over all atoms in the system, mi is the mass of the atom in amu, and RAi is the distance from the axis of rotation, A, to atom i in Ångstroms.

The axes of rotation are calculated as follows:

First, a 3 by 3 matrix, t, is constructed, with the elements of t being:

t1,1 = Y2 + Z2 = Σimi(yi2 + zi2) t1,2 = -X.Y = -Σimixiyi t2,2 = X2 + Z2 = Σimi(xi2 + zi2) t1,3 = -X.Z = -Σimixizi t2,3 = -Y.Z = -Σimiyizi t3,3 = X2 + Y2 = Σimi(xi2 + yi2)

where mi is the mass of the atom in amu, and xi, yi, and zi, are the Cartesian coordinates of the atoms, in Ångstroms.
Then t is diagonalized. The resulting eigenvalues, (amu Ångstrom2), are divided by N.A2, where N=Avogardo's number and A = number of Ångstroms in a centimeter, to give the moments of inertia in g.cm2. Because a useful unit is 10-40.g.cm2, the moments of inertia are multiplied by 1040 before being printed.
The eigenvectors associated with the eigenvalues are the axes of rotation, A, B, and C.

Useful conversion factors 1 g cm2 = 1.660540x10-40 (amu Ångstrom2) Rotational constants in cm-1: A = hN1016/(8π2c)/(amu Ångstrom2) A (in MHz) = 5.053791x105/(amu Ångstrom2) A (in cm-1) = 5.053791x105/c(amu Ångstrom2) = 16.85763/(amu Ångstrom2)

Parameters:

• coordinates (array) –

  cartesian coordinates of rank 2; [:,:]

• symbols (list) –

  list of atomic symbols (need to calculate the Molecular mass or total mass)

Returns:

• vector ( array ) –

  moment of inertia about axis of rotation A, B, and C;

• principle axis : array –

  array of rank2, principle axis of rotations, A, B, and C about which moment of inertia is calculated

ThermoIdealGas

Bases: ThermoDyn

calculate thermodynamic quantities within the ideal gas approximation

get_enthalpy(l_ignore=False)

calculate enthapy (H) for ideal gas; H (P,T) = H (T); pressure dependence is not there

if l_ignore is True, then enthalpy is set to zero .. 9/11/2024

Notes

H(P,T) = E + PV = E_electronic (E_DFT) + E_zpe + integrate over T (specific heat capacity at const P)

H(T) = Eelec + EZPE + ∫ C_P dT; integral from 0→T

Cp = kB + CV,trans + CV,rot + CV,vib + CV,elec

∫CV,vib dT = ∑ ϵi/exp ( ϵi/kBT−1) ; for linear mol = 3N-5 and 3N-6 nonlinear vibrations

get_entropy(l_ignore_S_trans=None)

calculate entropy S(T,P) = S(T,P_0) − kB ln P/P_0

Notes

S = S_trans + S_rot + S_elec + S_vib − kB ln (P/P_o)

1 Pa = 1 N / m2

get_Gibbs_free_energy()

calculate the free energy G at a given temp and pressue

Notes

(T,P ) = H(T) − TS (T,P) + E_DFT

E_DFT is energy from DFT calculation ror some other nonDFT calculations

ThermoHarmonic

Harmonic thermodynamics

NOTE:

Harmonic limit: taken from ASE website: The conversion of electronic structure calculation information into thermodynamic properties is less established for adsorbates. However, the simplest approach often taken is to treat all 3N degrees of freedom of the adsorbate harmonically since the adsorbate often has no real translational or rotational degrees of freedom.

U(T) = E_ ext{elec} + E_ ext{ZPE} + \sum_i^ ext{harm DOF} rac{\epsilon_i}{e^{\epsilon_i / k_ ext{B} T} - 1 }

S = k_ ext{B} \sum_i^ ext{harm DOF} \left[ rac{\epsilon_i}{k_ ext{B}T\left(e^{\epsilon_i/k_ ext{B}T}-1 ight)} - \ln \left( 1 - e^{-\epsilon_i/k_ ext{B}T} ight) ight]

i.e. U(T)=Eelec+EZPE+∑iharm DOF ϵi eϵi/kBT−1

S=kB ∑iharm DOF[ϵi kBT(e ϵi/kBT−1)−ln(1−e−ϵi/kBT)]

__init__(E_dft=None, E_zpe=None, vib_freq=None, temp=None)

Returns:

• see respecitve methods in the class –

get_zpe(ignore=False)

Returns the zero-point vibrational energy eV.

get_E_vib(l_ignore=False, set_zero=None)

find the in internal energy part due to vibrations from 0K to a given temperature in eV. All DOF are considered within the harmonic approximation

if l_ignore is true then internal energy contribution is set to zero set_zero also does the same thing as l_ignore, but in set_zero first inter. energy is calculated and then set to zero!

Output the vibration contribution to internal energy:

get_S_vib()

get vibrational entropy in eV/K (3N DOF) given in eV

output is entropy and entropy energy in eV/K and ev/k, respectively.

get_internal_energy()

find internal energy in eV at a given T in the harmonic approximation taking all DOF (3N) at a given temperature (K).

Notes

E(T) = U(T) = Eelec+EZPE+∑iharm DOF ϵi eϵi/kBT−1

get_Helmholtz_free_energy()

calculate the Helmholtz free energy

Notes

G(T ) = H(T) − TS (T) # + E_DFT # P depenedence negligble; PV term may be...

F(T) = U(T)−TS(T) -- Helmholtz free energy

assuming PV term to be negligible, we have G(T) equal to F(T)

If we assumes that the pV term in H=U+pV is negligible, then the Helmholtz free energy can be used to approximate the Gibbs free energy, as G=F+pV.

ThermoHindredRotVib

Hinder vibrational and rotation oscillator thermodynamics

ThermoLattice

Atomic lattice thermodynamics