Thermochemistry

ASE contains a thermochemistry module that lets the user derive commonly desired thermodynamic quantities of molecules and crystalline solids from ASE output and some user-specified parameters. The following classes are currently available in this module:

• The ideal-gas limit (in which translational and rotational degrees of freedom are taken into account) and the harmonic limit (generally used for adsorbates, in which all degrees of freedom are treated harmonically). See IdealGasThermo and HarmonicThermo.

• The hindered translator / hindered rotor model (used for adsorbates, in which two degrees of freedom are translational, one is rotational, and the remaining 3N-3 are vibrational). See HinderedThermo.

• A crystalline solid model (in which a lattice of N atoms is treated as a system of 3N independent harmonic oscillators). See CrystalThermo.

All non-periodic classes rely on good vibrational energies being fed to the calculators, which can be calculated with the vibrations module. Likewise, the crystalline solid model depends on an accurate phonon density of states; this is readily calculated using the phonons module.

If you want to calculate the thermodynamic properties of a transition state, you can use the alternative constructor BaseThermoChem.from_transition_state() to create an instance and automatically removing a single imaginary frequency.

ASE uses an approach, in which each individual vibrational mode is represented by an instance of an AbstractMode class. This class is a base class for the different types of modes (e.g. HarmonicMode or RRHOMode). The Thermochemistry classes then create a list of these modes and use them to calculate the thermodynamic properties. The user can also create their own modes by creating their own list of modes. Note, that not all Thermochemistry classes make use of this yet.

Ideal-gas limit

The thermodynamic quantities of ideal gases are calculated by assuming that all spatial degrees of freedom are independent and separable into translational, rotational, and vibrational degrees of freedom. The IdealGasThermo class supports calculation of enthalpy (\(H\)), entropy (\(S\)), and Gibbs free energy (\(G\)), and has the interface listed below.

class ase.thermochemistry.IdealGasThermo(vib_energies: Sequence[complex], geometry: Literal['linear', 'nonlinear', 'monatomic'], potentialenergy: float = 0.0, atoms: Atoms | None = None, symmetrynumber: int | None = None, spin: float | None = None, natoms: int | None = None, vib_selection: Literal['exact', 'all', 'highest', 'abs_highest'] | None = 'highest', ignore_imag_modes: bool = False, modes: Sequence[AbstractMode] | None = None)

Class for calculating thermodynamic properties of a molecule based on statistical mechanical treatments in the ideal gas approximation.

Inputs for enthalpy calculations:

vib_energieslist

a list of the vibrational energies of the molecule (e.g., from ase.vibrations.Vibrations.get_energies). The number of expected vibrations is calculated by the geometry and the number of atoms. By default, the lowest energies are cut to enforce the expected length. By setting vib_selection, this behaviour can be changed. If either atoms and natoms is unspecified, the full list is used. Units are eV.

geometry‘monatomic’, ‘linear’, or ‘nonlinear’

geometry of the molecule

potentialenergyfloat

the potential energy in eV (e.g., from atoms.get_potential_energy) (if potentialenergy is unspecified, then the methods of this class can be interpreted as the energy corrections)

natomsinteger

the number of atoms, used along with ‘geometry’ to determine how many vibrations to use. (Not needed if an atoms object is supplied in ‘atoms’ or if the user desires the entire list of vibrations to be used.)

vib_selection‘exact’, ‘highest’ (default), ‘abs_highest’, ‘all’

selection of input vibrational energies considered to be true vibrations (excluding translations and rotations) implied by the geometry and number of atoms. Only applied if number of atoms is provided, either with natoms or atoms. - ‘exact’: no selection; the number of input vibrational energies must be equal to the number of true vibrations - ‘highest’: select vibrational energies whose square is the highest , i.e. large real energies followed by small real energies, small imaginary energies and large imaginary energies. Will not catch unrealistically large imaginary frequencies, resulting from e.g. unrelaxed molecules. This is the default. - ‘abs_highest’: select vibrational energies whose absolute number are highest, real or imaginary. Will fail intentionally for unrealistically large imaginary frequencies, but also unintentionally for small real energies, corresponding e.g. to frustrated rotations in a larger molecule. - ‘all’: Use all input vibrational energies without checking the number of them.

ignore_imag_modesbool

If ‘True’, any imaginary frequencies present after the 3N-3 cut will be removed in the calculation of the thermochemical properties. If ‘False’ (default), an error will be raised if imaginary frequencies are present after the 3N-3 cut.

Extra inputs needed for entropy / free energy calculations:

atomsan ASE atoms object

used to calculate rotational moments of inertia and molecular mass

symmetrynumberinteger

symmetry number of the molecule. See, for example, Table 10.1 and Appendix B of C. Cramer “Essentials of Computational Chemistry”, 2nd Ed.

spinfloat

the total electronic spin. (0 for molecules in which all electrons are paired, 0.5 for a free radical with a single unpaired electron, 1.0 for a triplet with two unpaired electrons, such as O_2.)

get_enthalpy(temperature: float, verbose: bool = True) → float

Returns the enthalpy, in eV, in the ideal gas approximation at a specified temperature (K).

get_entropy(temperature: float, pressure: float = 6.241509125883258e-07, verbose: bool = True) → float

Returns the entropy, in eV/K, in the ideal gas approximation at a specified temperature (K) and pressure (Pa).

get_gibbs_energy(temperature: float, pressure: float, verbose: bool = True) → float

Returns the Gibbs free energy, in eV, in the ideal gas approximation at a specified temperature (K) and pressure (Pa).

get_internal_energy(temperature: float, verbose: bool = True) → float

Returns the internal energy, in eV, in the ideal gas approximation at a specified temperature (K).

Example

The IdealGasThermo class would generally be called after an energy optimization and a vibrational analysis. The user needs to supply certain parameters if the entropy or free energy are desired, such as the geometry and symmetry number. An example on the nitrogen molecule is:

from ase.build import molecule
from ase.calculators.emt import EMT
from ase.optimize import QuasiNewton
from ase.thermochemistry import IdealGasThermo
from ase.vibrations import Vibrations

atoms = molecule('N2')
atoms.calc = EMT()
dyn = QuasiNewton(atoms)
dyn.run(fmax=0.01)
potentialenergy = atoms.get_potential_energy()

vib = Vibrations(atoms)
vib.run()
vib_energies = vib.get_energies()

thermo = IdealGasThermo(
    vib_energies=vib_energies,
    potentialenergy=potentialenergy,
    atoms=atoms,
    geometry='linear',
    symmetrynumber=2,
    spin=0,
    vib_selection='highest',
)
G = thermo.get_gibbs_energy(temperature=298.15, pressure=101325.0)

This will give the thermodynamic summary output:

Enthalpy components at T = 298.15 K:
===============================
Enthalpy components at T = 298.15 K:
===============================
E_pot                  0.263 eV
E_ZPE                  0.076 eV
Cv_trans (0->T)        0.039 eV
Cv_rot (0->T)          0.026 eV
Cv_vib (0->T)          0.000 eV
-------------------------------
U                      0.404 eV
===============================
(C_v -> C_p)           0.026 eV
-------------------------------
H                      0.429 eV
===============================

Entropy components at T = 298.15 K and P = 101325.0 Pa:
=================================================
{"":15s}{"S":13s}     {"T*S:13s}
S_trans (1 bar)    0.0015590 eV/K        0.465 eV
S_rot              0.0004101 eV/K        0.122 eV
S_elec             0.0000000 eV/K        0.000 eV
S_vib              0.0000016 eV/K        0.000 eV
S (1 bar -> P)    -0.0000011 eV/K       -0.000 eV
-------------------------------------------------
S                  0.0019695 eV/K        0.587 eV
=================================================

Free energy components at T = 298.15 K and P = 101325.0 Pa:
=======================
    H          0.429 eV
 -T*S         -0.587 eV
-----------------------
    G         -0.158 eV
=======================

Harmonic limit

In the harmonic limit, all degrees of freedom are treated harmonically. The HarmonicThermo class supports the calculation of internal energy, entropy, and free energy. This class returns the Helmholtz free energy; if the user assumes the pV term (in H = U + pV) is zero this can also be interpreted as the Gibbs free energy. This class uses all of the energies given to it in the vib_energies list; this is a list as can be generated with the .get_energies() method of ase.vibrations.Vibrations, but the user should take care that all of these energies are real (non-imaginary). If imaginary values are encountered, by default an error is raised. By setting ignore_imag_modes to either ‘True’ or ‘False’, the imaginary modes can be deleted or an error raised, respectively. Note, that the HarmonicThermo class does not calculate \(S_\text{trans}\), \(S_\text{rot}\) and \(S_\text{elec}\). If these are desired, the user should calculate them separately and add them to the entropy. Use the functions IdealGasThermo.get_ideal_entropy() for that purpose. The class HarmonicThermo has the interface described below.

class ase.thermochemistry.HarmonicThermo(vib_energies: Sequence[complex], potentialenergy: float = 0.0, ignore_imag_modes: bool = False, modes: Sequence[AbstractMode] | None = None)

Class for calculating thermodynamic properties in the approximation that all degrees of freedom are treated harmonically. Often used for adsorbates.

Note: This class does not include the translational and rotational contributions to the entropy by default. Use the get_ideal_entropy method for that and add them manually.

Inputs:

vib_energieslist

a list of the harmonic energies of the adsorbate (e.g., from ase.vibrations.Vibrations.get_energies). The number of energies should match the number of degrees of freedom of the adsorbate; i.e., 3*n, where n is the number of atoms. Note that this class does not check that the user has supplied the correct number of energies. Units of energies are eV. Note, that if ‘modes’ are given, these energies are not used.

potentialenergyfloat

the potential energy in eV (e.g., from atoms.get_potential_energy) (if potentialenergy is unspecified, then the methods of this class can be interpreted as the energy corrections)

ignore_imag_modesbool

If ‘True’, any imaginary frequencies will be removed in the calculation of the thermochemical properties. If ‘False’ (default), an error will be raised if any imaginary frequencies are present.

modeslist of AbstractMode

A list of mode objects. If not provided, HarmonicMode objects will be created from the vib_energies. This is useful if you want to replace individual modes with non-harmonic modes.

get_entropy(temperature: float, pressure: float = 6.241509125883258e-07, verbose: bool = True) → float

Returns the entropy, in eV/K at a specified temperature (K).

Note: This does not include the translational and rotational contributions to the entropy. Use the get_ideal_entropy method for that.

Parameters:

• temperature (float) – The temperature in Kelvin.

• pressure (float) – Not used, but kept for compatibility with other classes.

• verbose (bool) – If True, print the contributions to the entropy.

Return type:

float

get_helmholtz_energy(temperature: float, verbose: bool = True) → float

Returns the Helmholtz free energy, in eV, in the harmonic approximation at a specified temperature (K).

get_internal_energy(temperature: float, verbose: bool = True) → float

Returns the internal energy, in eV, in the harmonic approximation at a specified temperature (K).

Derivations

There are some derivations of the standard procedure for calculating the thermodynamic properties in the harmonic limit. Currently, ASE provides a class QuasiHarmonicThermo based on the quasi-Harmonic approximation by Cramer, Truhlar and coworkers (doi: 10.1021/jp205508z) and a MSRRHOThermo based on the modified rigid-rotor-harmonic-oscillator (msRRHO) approximation by Grimme et al. (doi: 10.1002/chem.201200497 and doi: 10.1039/D1SC00621E).

class ase.thermochemistry.QuasiHarmonicThermo(vib_energies: Sequence[complex], potentialenergy: float = 0.0, ignore_imag_modes: bool = False, modes: Sequence[AbstractMode] | None = None, raise_to: float = 0.01239841973964072)

Subclass of HarmonicThermo, including the quasi-harmonic approximation of Cramer, Truhlar and coworkers doi: 10.1021/jp205508z.

Note: This class not include the translational and rotational contributions to the entropy by default. Use the get_ideal_entropy method for that and add them manually.

Inputs:

vib_energieslist

a list of the energies of the vibrations (e.g., from ase.vibrations.Vibrations.get_energies). The number of energies should match the number of degrees of freedom of the adsorbate; i.e., 3*n, where n is the number of atoms. Note that this class does not check that the user has supplied the correct number of energies. Units of energies are eV.

potentialenergyfloat

the potential energy in eV (e.g., from atoms.get_potential_energy) (if potentialenergy is unspecified, then the methods of this class can be interpreted as the energy corrections)

ignore_imag_modesbool

If ‘True’, any imaginary frequencies will be removed in the calculation of the thermochemical properties. If ‘False’ (default), an error will be raised if any imaginary frequencies are present.

modeslist of AbstractMode

A list of mode objects. If not provided, HarmonicMode objects will be created from the raised vib_energies. This is useful if you want to replace individual modes with non-harmonic modes.

raise_tofloat

The value to which all frequencies smaller than this value will be raised. Unit is eV. Defaults to \(100cm^{-1} = 0.012398 eV\).

class ase.thermochemistry.MSRRHOThermo(vib_energies: Sequence[complex], atoms: Atoms, potentialenergy: float = 0.0, tau: float = 35.0, nu_scal: float = 1.0, treat_int_energy: bool = False, modes: Sequence[AbstractMode] | None = None)

Subclass of QuasiHarmonicThermo, including Grimme’s scaling method based on doi: 10.1002/chem.201200497 and doi: 10.1039/D1SC00621E.

Note: This class not include the translational and rotational contributions to the entropy by default. Use the get_ideal_entropy method for that and add them manually.

We enforce treating imaginary modes as Grimme suggests (converting them to real by multiplying them with \(-i\)). So make sure to check your input energies.

Inputs:

vib_energieslist

a list of the energies of the vibrations (e.g., from ase.vibrations.Vibrations.get_energies). The number of energies should match the number of degrees of freedom of the atoms object; i.e., 3*n, where n is the number of atoms. Note that this class does not check that the user has supplied the correct number of energies. Units of energies are eV. Note, that if ‘modes’ are given, these energies are not used.

atoms: an ASE atoms object

used to calculate rotational moments of inertia and molecular mass

taufloat

the vibrational energy threshold in \(cm^{-1}\), namcomplexed \(\tau\) in doi: 10.1039/D1SC00621E. Values close or equal to 0 will result in the standard harmonic approximation. Defaults to \(35cm^{-1}\).

nu_scalfloat

Linear scaling factor for the vibrational frequencies. Named \(\nu_{scal}\) in doi: 10.1039/D1SC00621E. Defaults to 1.0, check the Truhlar group database for values corresponding to your level of theory. Note that for \(\nu_{scal}=1.0\) this method is equivalent to the quasi-RRHO method in doi: 10.1002/chem.201200497.

treat_int_energybool

Extend the msRRHO treatement to the internal energy. If False, only the entropy contribution as in Grimmmes paper is considered. If true, the approach of Otlyotov and Minenkov doi: 10.1002/jcc.27129 is used.

modeslist of AbstractMode

A list of mode objects. If not provided, RRHOMode objects will be created from the raised vib_energies. This is useful if you want to replace individual modes with non-harmonic modes.

Hindered translator / hindered rotor model

The hindered translator / hindered rotor model bridges the gap between the 2D gas (i.e. free translator / free rotor) and the 2D lattice gas (i.e. harmonic oscillator). For an adsorbate containing N atoms, two degrees of freedom are treated as hindered translations in the two directions parallel to the surface, one degree of freedom is treated as a hindered rotation about the axis perpendicular to the surface, and the remaining 3N-3 degrees of freedom are treated as vibrations. The HinderedThermo class supports the calculation of internal energy, entropy, free energy, and zero point energy (included in the internal energy). All of the thermodynamic properties calculated here are at the standard state surface concentration (defined here such that a 2D ideal gas at that concentration has 2/3 the translational entropy of a 3D ideal gas at 1 bar pressure, so that \(\theta^0\) = 0.012 at 298 K for a surface with \(10^{15}\) sites/cm²). This class returns the Helmholtz free energy; if the user assumes that the pV term (in G = U + pV - TS) is zero then this free energy can also be interpreted as the Gibbs free energy. This class depends on the user defined translation barrier (trans_barrier_energy) and rotational barrier (rot_barrier_energy) for the adsorbate to move on the surface in order to calculate the translational and rotational degrees of freedom. To calculate the vibrational degrees of freedom, all 3N vibrational energies must be supplied in the vib_energies list and the 3N-3 largest vibrational energies are used to calculate the vibrational contribution; this is a list as can be generated with the .get_energies() method of ase.vibrations.Vibrations. The class HinderedThermo has the interface described below.

class ase.thermochemistry.HinderedThermo(vib_energies, trans_barrier_energy, rot_barrier_energy, sitedensity, rotationalminima, potentialenergy=0.0, mass=None, inertia=None, atoms=None, symmetrynumber=1, ignore_imag_modes: bool = False)

Class for calculating thermodynamic properties in the hindered translator and hindered rotor model where all but three degrees of freedom are treated as harmonic vibrations, two are treated as hindered translations, and one is treated as a hindered rotation.

Inputs:

vib_energieslist

a list of all the vibrational energies of the adsorbate (e.g., from ase.vibrations.Vibrations.get_energies). If atoms is not provided, then the number of energies must match the number of degrees of freedom of the adsorbate; i.e., 3*n, where n is the number of atoms. Note that this class does not check that the user has supplied the correct number of energies. Units of energies are eV.

trans_barrier_energyfloat

the translational energy barrier in eV. This is the barrier for an adsorbate to diffuse on the surface.

rot_barrier_energyfloat

the rotational energy barrier in eV. This is the barrier for an adsorbate to rotate about an axis perpendicular to the surface.

sitedensityfloat

density of surface sites in cm^-2

rotationalminimainteger

the number of equivalent minima for an adsorbate’s full rotation. For example, 6 for an adsorbate on an fcc(111) top site

potentialenergyfloat

the potential energy in eV (e.g., from atoms.get_potential_energy) (if potentialenergy is unspecified, then the methods of this class can be interpreted as the energy corrections)

massfloat

the mass of the adsorbate in amu (if mass is unspecified, then it will be calculated from the atoms class)

inertiafloat

the reduced moment of inertia of the adsorbate in amu*Ang^-2 (if inertia is unspecified, then it will be calculated from the atoms class)

atomsan ASE atoms object

used to calculate rotational moments of inertia and molecular mass

symmetrynumberinteger

symmetry number of the adsorbate. This is the number of symmetric arms of the adsorbate and depends upon how it is bound to the surface. For example, propane bound through its end carbon has a symmetry number of 1 but propane bound through its middle carbon has a symmetry number of 2. (if symmetrynumber is unspecified, then the default is 1)

ignore_imag_modesbool

If ‘True’, any imaginary frequencies present after the 3N-3 cut will be removed in the calculation of the thermochemical properties. If ‘False’ (default), an error will be raised if imaginary frequencies are present after the 3N-3 cut.

get_entropy(temperature, pressure=6.241509125883258e-07, verbose=True)

Returns the entropy, in eV/K, in the hindered translator and hindered rotor model at a specified temperature (K).

get_helmholtz_energy(temperature, verbose=True)

Returns the Helmholtz free energy, in eV, in the hindered translator and hindered rotor model at a specified temperature (K).

get_internal_energy(temperature, verbose=True)

Returns the internal energy (including the zero point energy), in eV, in the hindered translator and hindered rotor model at a specified temperature (K).

get_zero_point_energy(verbose=True)

Returns the zero point energy, in eV, in the hindered translator and hindered rotor model

Example

The HinderedThermo class would generally be called after an energy optimization and a vibrational analysis. The user needs to supply certain parameters, such as the vibrational energies, translational energy barrier, rotational energy barrier, surface site density, number of equivalent minima in a full rotation, and the number of symmetric arms of the adsorbate as it rotates on the surface. The user also needs to supply either the mass of the adsorbate and the reduced moment of inertia of the adsorbate as it rotates on the surface or the user can supply the atoms object from which the mass and an approximate reduced moment of inertia may be determined. An example for ethane on a platinum (111) surface is:

from numpy import array

from ase.thermochemistry import HinderedThermo

vibs = array(
    [
        3049.060670,
        3040.796863,
        3001.661338,
        2997.961647,
        2866.153162,
        2750.855460,
        1436.792655,
        1431.413595,
        1415.952186,
        1395.726300,
        1358.412432,
        1335.922737,
        1167.009954,
        1142.126116,
        1013.918680,
        803.400098,
        783.026031,
        310.448278,
        136.112935,
        112.939853,
        103.926392,
        77.262869,
        60.278004,
        25.825447,
    ]
)
vib_energies = vibs / 8065.54429  # convert to eV from cm^-1
trans_barrier_energy = 0.049313  # eV
rot_barrier_energy = 0.017675  # eV
sitedensity = 1.5e15  # cm^-2
rotationalminima = 6
symmetrynumber = 1
mass = 30.07  # amu
inertia = 73.149  # amu Ang^-2

thermo = HinderedThermo(
    vib_energies=vib_energies,
    trans_barrier_energy=trans_barrier_energy,
    rot_barrier_energy=rot_barrier_energy,
    sitedensity=sitedensity,
    rotationalminima=rotationalminima,
    symmetrynumber=symmetrynumber,
    mass=mass,
    inertia=inertia,
)

F = thermo.get_helmholtz_energy(temperature=298.15)

This will give the thermodynamic summary output:

Internal energy components at T = 298.15 K:
===============================
E_pot                  0.000 eV
E_trans                0.049 eV
E_rot                  0.018 eV
E_vib                  0.076 eV
E_ZPE                  1.969 eV
-------------------------------
U                      2.112 eV
===============================

Entropy components at T = 298.15 K:
=================================================
                           S               T*S
S_trans            0.0005074 eV/K        0.151 eV
S_rot              0.0002287 eV/K        0.068 eV
S_vib              0.0005004 eV/K        0.149 eV
S_con              0.0005044 eV/K        0.150 eV
-------------------------------------------------
S                  0.0017409 eV/K        0.519 eV
=================================================

Free energy components at T = 298.15 K:
=======================
    U          2.112 eV
 -T*S         -0.519 eV
-----------------------
    F          1.593 eV
=======================

Crystals

In this model a crystalline solid is treated as a periodic system of independent harmonic oscillators. The CrystalThermo class supports the calculation of internal energy (\(U\)), entropy (\(S\)) and Helmholtz free energy (\(F\)), and has the interface listed below.

class ase.thermochemistry.CrystalThermo(phonon_DOS, phonon_energies, formula_units=None, potentialenergy=0.0)

Class for calculating thermodynamic properties of a crystalline solid in the approximation that a lattice of N atoms behaves as a system of 3N independent harmonic oscillators.

Inputs:

phonon_DOSlist

a list of the phonon density of states, where each value represents the phonon DOS at the vibrational energy value of the corresponding index in phonon_energies.

phonon_energieslist

a list of the range of vibrational energies (hbar*omega) over which the phonon density of states has been evaluated. This list should be the same length as phonon_DOS and integrating phonon_DOS over phonon_energies should yield approximately 3N, where N is the number of atoms per unit cell. If the first element of this list is zero-valued it will be deleted along with the first element of phonon_DOS. Units of vibrational energies are eV.

potentialenergyfloat

the potential energy in eV (e.g., from atoms.get_potential_energy) (if potentialenergy is unspecified, then the methods of this class can be interpreted as the energy corrections)

formula_unitsint

the number of formula units per unit cell. If unspecified, the thermodynamic quantities calculated will be listed on a per-unit-cell basis.

classmethod from_transition_state(*args, **kwargs) → CrystalThermo

Not yet implemented for CrystalThermo but needed to overwrite BaseThermoChem method.

get_entropy(temperature, verbose=True)

Returns the entropy, in eV/K, of crystalline solid at a specified temperature (K).

get_helmholtz_energy(temperature, verbose=True)

Returns the Helmholtz free energy, in eV, of crystalline solid at a specified temperature (K).

get_internal_energy(temperature, verbose=True)

Returns the internal energy, in eV, of crystalline solid at a specified temperature (K).

Example

The CrystalThermo class will generally be called after an energy optimization and a phonon vibrational analysis of the crystal. An example for bulk gold is:

from ase.calculators.emt import EMT
from ase.optimize import QuasiNewton
from ase.phonons import Phonons
from ase.spacegroup import crystal
from ase.thermochemistry import CrystalThermo

# Set up gold bulk and attach EMT calculator
a = 4.078
atoms = crystal(
    'Au',
    (0.0, 0.0, 0.0),
    spacegroup=225,
    cellpar=[a, a, a, 90, 90, 90],
    pbc=(1, 1, 1),
)
calc = EMT()
atoms.calc = calc
qn = QuasiNewton(atoms)
qn.run(fmax=0.05)
potentialenergy = atoms.get_potential_energy()

# Phonon analysis
N = 5
ph = Phonons(atoms, calc, supercell=(N, N, N), delta=0.05)
try:
    ph.run()
    ph.read(acoustic=True)
finally:
    ph.clean()
dos = ph.get_dos(kpts=(40, 40, 40)).sample_grid(npts=3000, width=5e-4, xmin=0.0)
phonon_energies = dos.get_energies()
phonon_DOS = dos.get_weights()


# Calculate the Helmholtz free energy
thermo = CrystalThermo(
    phonon_energies=phonon_energies,
    phonon_DOS=phonon_DOS,
    potentialenergy=potentialenergy,
    formula_units=4,
)
F = thermo.get_helmholtz_energy(temperature=298.15)

This will give the thermodynamic summary output:

Internal energy components at T = 298.15 K,
on a per-formula-unit basis:
===============================
E_pot                 0.0022 eV
E_ZPE                 0.0138 eV
E_phonon              0.0642 eV
-------------------------------
U                     0.0801 eV
===============================

Entropy components at T = 298.15 K,
on a per-formula-unit basis:
=================================================
                           S               T*S
-------------------------------------------------
S                  0.0005397 eV/K       0.1609 eV
=================================================

Helmholtz free energy components at T = 298.15 K,
on a per-formula-unit basis:
=======================
    U         0.0801 eV
 -T*S        -0.1609 eV
-----------------------
    F        -0.0808 eV
=======================

Imaginary Frequencies

In all thermochemistry classes, imaginary frequencies can be handled in two ways. By default, an error is raised if any imaginary frequencies are found. However, by setting the ignore_imag_modes parameter to either ‘True’, the imaginary modes can be deleted. For transition states, the alternative constructor BaseThermoChem.from_transition_state() can be used to create an instance and automatically removing a single imaginary frequency.

However, often imaginary frequencies are artifacts of the numerical calculation and cannot be avoided. In the literature, several approaches have been proposed to deal with these frequencies. One approach is to invert the imaginary frequencies by multiplying them with \(-i\), this is the approach taken in the MSRRHOThermo class. Below you find a variety of one-liner code snippets showing how this could be achieved manually:

vib_energies = np.abs(vib_energies)
vib_energies = vib_energies.real + vib_energies.imag
vib_energies[vib_energies.imag > 0] = vib_energies[vib_energies.imag > 0].imag
vib_energies[np.abs(vib_energies.imag) > 0] = vib_energies[np.abs(vib_energies.imag) > 0].imag

Background

Ideal gas

The conversion of electronic structure calculations to thermodynamic properties in the ideal-gas limit is well documented; see, for example, Chapter 10 of Cramer, 2004. The key equations used in the IdealGasThermo class are summarized here.

C.J. Cramer. Essentials of Computational Chemistry, Second Edition. Wiley, 2004.

The ideal-gas enthalpy is calculated from extrapolation of the energy at 0 K to the relevant temperature (for an ideal gas, the enthalpy is not a function of pressure):

\[H(T) = E_\text{elec} + E_\text{ZPE} + \int_0^\text{T} C_P \, \text{d}T\]

where the first two terms are the electronic energy and the zero-point energy, and the integral is over the constant-pressure heat capacity. The heat capacity is separable into translational, rotational, vibrational, and electronic parts (plus a term of \(k_\text{B}\) to switch from constant-volume to constant-pressure):

\[C_P = k_\text{B} + C_{V\text{,trans}} + C_{V\text{,rot}} + C_{V\text{,vib}} + C_{V\text{,elec}}\]

The translational heat capacity is 3/2 \(k_\text{B}\) for a 3-dimensional gas. The rotational heat capacity is 0 for a monatomic species, \(k_\text{B}\) for a linear molecule, and 3/2 \(k_\text{B}\) for a nonlinear molecule. In this module, the electronic component of the heat capacity is assumed to be 0. The vibrational 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). The integrated form of the vibrational heat capacity is:

\[\int_0^T C_{V,\text{vib}} \text{d}T = \sum_i^\text{vib DOF} \frac{\epsilon_i}{e^{\epsilon_i / k_\text{B} T} - 1 }\]

where \(\epsilon_i\) are the energies associated with the vibrational frequencies, \(\epsilon_i = h \omega_i\).

The ideal gas entropy can be calculated as a function of temperature and pressure as:

\[\begin{split}S(T,P) &= S(T,P^\circ) - k_\text{B} \ln \frac{P}{P^\circ} \\ &= S_\text{trans} + S_\text{rot} + S_\text{elec} + S_\text{vib} - k_\text{B} \ln \frac{P}{P^\circ}\end{split}\]

where the translational, rotational, electronic, and vibrational components are calculated as below. (Note that the translational component also includes components from the Stirling approximation, and that the vibrational degrees of freedom are enumerated the same as in the above.)

\[S_\text{trans} = k_\text{B} \left\{ \ln \left[ \left( \frac{2 \pi M k_\text{B} T}{h^2} \right)^{3/2} \frac{k_\text{B} T}{P^\circ} \right] + \frac{5}{2} \right\}\]

\[\begin{split}S_\text{rot} = \left\{ \begin{array}{ll} 0 & \text{, if monatomic} \\ k_\text{B} \left[ \ln \left( \frac{8\pi^2 I k_\text{B}T}{\sigma h^2}\right) + 1 \right] & \text{, if linear} \\ k_\text{B} \left\{ \ln \left[ \frac{\sqrt{\pi I_\text{A} I_\text{B} I_\text{C}}}{\sigma} \left(\frac{8\pi^2 k_\text{B} T}{h^2}\right)^{3/2}\right] + \frac{3}{2} \right\} & \text{, if nonlinear} \\ \end{array} \right.\end{split}\]

\[S_\text{vib} = k_\text{B} \sum_i^\text{vib DOF} \left[ \frac{\epsilon_i}{k_\text{B}T\left(e^{\epsilon_i/k_\text{B}T}-1\right)} - \ln \left( 1 - e^{-\epsilon_i/k_\text{B}T} \right)\right]\]

\[S_\text{elec} = k_\text{B} \ln \left[ 2 \times \left(\text{total spin}\right) + 1\right]\]

\(I_\text{A}\) through \(I_\text{C}\) are the three principle moments of inertia for a non-linear molecule. \(I\) is the degenerate moment of inertia for a linear molecule. \(\sigma\) is the symmetry number of the molecule.

The ideal-gas Gibbs free energy is then just calculated from the combination of the enthalpy and entropy:

\[G(T,P) = H(T) - T\, S(T,P)\]

Note

Note, that in comparison to other implementations the following can play a role:

• The reference pressure used in ASE is 1.0e5 Pa (1 bar). Other codes might use 1.0 atm. This can lead to small differences in the calculated ideal entropy contributions. You can manually overwrite the value of the reference pressure in the BaseThermoChem instances and its subclasses by setting BaseThermoChem.referencepressure after initialization.

• ASE uses the masses of the atoms object. These might be slightly different from the masses used in other codes. You can overwrite them in the Atoms object.

• Some codes use an average inertia of \(B_\text{av}=10^{-44}\) kg m² as in doi: 10.1002/chem.201200497, ASE calculates the real mean inertia by calling get_moments_of_inertia().

Harmonic limit

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. This is the approach implemented in the HarmonicThermo class. Thus, the internal energy and entropy of the adsorbate are calculated as

\[U(T) = E_\text{elec} + E_\text{ZPE} + \sum_i^\text{harm DOF} \frac{\epsilon_i}{e^{\epsilon_i / k_\text{B} T} - 1 }\]

\[S = k_\text{B} \sum_i^\text{harm DOF} \left[ \frac{\epsilon_i}{k_\text{B}T\left(e^{\epsilon_i/k_\text{B}T}-1\right)} - \ln \left( 1 - e^{-\epsilon_i/k_\text{B}T} \right)\right]\]

and the Helmholtz free energy is calculated as

\[F(T) = U(T) - T\, S(T)\]

In this case, the number of harmonic energies (\(\epsilon_i\)) used in the summation is generally \(3N\), where \(N\) is the number of atoms in the adsorbate. If the user 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\).

Hindered translator / hindered rotor

The conversion of electronic structure calculations to thermodynamic properties in the hindered translator / hindered rotor model was developed for adsorbates on close packed surfaces and is documented by Sprowl, Campbell, and Arnadottir, 2016. The key equations used in the HinderedThermo class are summarized here.

L.H. Sprowl, C.T. Campbell, and L. Arnadottir. Hindered Translator and Hindered Rotor Models for Adsorbates: Partition Functions and Entropies. J. Phys. Chem. C, 2016, 120 (18), pp 9719-9731.

L.H. Sprowl, C.T. Campbell, and L. Arnadottir. Correction to “Hindered Translator and Hindered Rotor Models for Adsorbates: Partition Functions and Entropies”. J. Phys. Chem. C, 2017, 121 (17), pp 9655-9655.

C.T. Campbell, L.H. Sprowl, and L. Arnadottir. Equilibrium Constants and Rate Constants for Adsorbates: Two-Dimensional (2D) Ideal Gas, 2D Ideal Lattice Gas, and Ideal Hindered Translator Models. J. Phys. Chem. C, 2016, 120 (19), pp 10283-10297.

The \(3N-3\) largest vibrational frequencies are used to calculate the vibrational contributions to the internal energy and the entropy. The remaining three degrees of freedom are calculated from two translational contributions and one rotational contribution of the adsorbate. The energy barriers for the adsorbate to translate and rotate on a close packed surface are used to calculate the translational and rotational frequencies, respectively. From the translational and rotational frequencies, the translational and rotational contributions to the internal energy and the entropy of the adsorbate are determined. The calculation of the translational frequency is:

\[\nu_{trans} = \sqrt{\frac{W_{trans}}{2mA}}\]

where \(W_{trans}\) is the translational energy barrier, \(m\) is the mass of the adsorbate, and \(A\) is the area per surface atom, or the inverse of the surface site density. The rotational frequency is calculated as:

\[\nu_{rot} = \frac{1}{2\pi}\sqrt{\frac{n^2W_{rot}}{2I}}\]

where \(W_{rot}\) is the rotational energy barrier, \(n\) is the number of equivalent energy minima in a full rotation of the adsorbate, and \(I\) is the reduced moment of inertia of the adsorbate about its surface bond. Two variables are now introduced, a unitless temperature

\[T_i = \frac{kT}{h\nu_i}\]

and a unitless energy barrier

\[r_i = \frac{W_i}{h\nu_i}\]

to ease the internal energy and entropy calculations.

The internal energy of the adsorbate is calculated as:

\[U(T) = E_\text{elec} + E_\text{ZPE} + E_\text{trans} + E_\text{rot} + E_\text{vib}\]

where \(E_{trans}\) and \(E_{rot}\) are:

\[E_i = k_\text{B}T \left( \frac{1/T_i}{\exp\left[1/T_i\right]-1} -\frac{1}{2} - \frac{1}{\left(2+16r_i\right)T_i} + \frac{r_i}{2T_i} \left( 1 - \frac{\text{I}_1\left[r_i/2T_i\right]}{\text{I}_0\left[r_i/2T_i\right]}\right) \right)\]

where \(I_{n}\) is the nth-order modified Bessel function of the first kind. Similarly for the harmonic limit, \(E_{vib}\) is:

\[E_\text{vib} = k_\text{B}T \sum_i^\text{3N-3} \left( \frac{1/T_i}{\exp\left[1/T_i\right]-1} \right)\]

The entropy of the adsorbate is calculated as:

\[S = S_\text{trans} + S_\text{rot} + S_\text{vib} + S_\text{con}\]

where \(S_{trans}\) and \(S_{rot}\) are:

\[S_i = k_\text{B} \left( \frac{1/T_i}{\exp\left[1/T_i\right]-1} - \ln \left[ 1 - \exp\left[-\frac{1}{T_i}\right]\right] - \frac{1}{2} - \frac{r_i}{2T_i}\frac{\text{I}_1\left[r_i/2T_i\right]}{\text{I}_0\left[r_i/2T_i\right]} + \ln\left[\left(\frac{\pi r_i}{T_i}\right)^{1/2}\text{I}_0\left[\frac{r_i}{2T_i}\right]\right] \right)\]

and \(S_{vib}\) is:

\[S_\text{vib} = k_\text{B} \sum_i^\text{3N-3} \left( \frac{1/T_i}{\exp\left[1/T_i\right]-1} - \ln \left[ 1 - \exp\left[-\frac{1}{T_i}\right]\right] \right)\]

\(S_{con}\) is a concentration related entropy and is calculated as:

\[S_\text{con} = k_\text{B} \left( 1 - \ln\left[A\left(\frac{N}{A}\right)^0\right] \right)\]

where

\[\left(\frac{N}{A}\right)^0 = e^{1/3}\left(\frac{N_A \text{ 1 bar}}{RT}\right)\]

The Helmholtz free energy is calculated as:

\[F(T) = U(T) - T\, S(T)\]

If the user 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\).

Crystalline solid

The derivation of the partition function for a crystalline solid is fairly straight-forward and can be found, for example, in Chapter 11 of McQuarrie, 2000.

D.A. McQuarrie. Statistical Mechanics. University Science Books, 2000.

The treatment implemented in the CrystalThermo class depends on introducing normal coordinates to the entire crystal and treating each atom in the lattice as an independent harmonic oscillator. This yields the partition function

\[Z = \prod_{j=1}^\text{3N} \left( \frac{e^{-\frac{1}{2}\epsilon_j/k_\text{B}T}}{1 - e^{-\epsilon_j/k_\text{B}T}} \right) e^{-E_\text{elec} / k_\mathrm{B}T}\]

where \(\epsilon_j\) are the \(3N\) vibrational energy levels and \(E_\text{elec}\) is the electronic energy of the crystalline solid. Now, taking the logarithm of the partition function and replacing the resulting sum with an integral (assuming that the energy level spacing is essentially continuous) gives

\[-\ln Z = E_\text{elec}/k_\text{B}T + \int_0^\infty \left[ \ln \left( 1 - e^{-\epsilon/k_\text{B}T} \right) + \frac{\epsilon}{2 k_\text{B} T} \right]\sigma (\epsilon) \text{d}\epsilon\]

Here \(\sigma (\epsilon)\) represents the degeneracy or phonon density of states as a function of vibrational energy. Once this function has been determined (i.e. using the ase.phonons module), it is a simple matter to calculate the canonical ensemble thermodynamic quantities; namely the internal energy, the entropy and the Helmholtz free energy.

\[\begin{split}U(T) &= -\left( \frac{\partial \ln Z}{\partial \frac{1}{k_\text{B}T} } \right)_\text{N,V} \\ &= E_\text{elec} + \int_0^\infty \left[ \frac{\epsilon}{e^{\epsilon/k_\text{B}T} - 1} + \frac{\epsilon}{2} \right]\sigma (\epsilon) \text{d}\epsilon\end{split}\]

\[\begin{split}S(T) &= \frac{U}{T} + k_\text{B} \ln Z \\ &= \int_0^\infty \left[ \frac{\epsilon}{T} \frac{1}{e^{\epsilon/k_\text{B}T} - 1} - k_\text{B} \ln \left(1 - e^{-\epsilon/k_\text{B}T} \right) \right]\sigma (\epsilon) \text{d}\epsilon\end{split}\]

\[F(T) = U(T) - T\, S(T,P)\]

Abstract Base Classes

class ase.thermochemistry.BaseThermoChem(vib_energies: Sequence[float] | None = None, atoms: Atoms | None = None, modes: Sequence[AbstractMode] | None = None, spin: float | None = None)

Abstract base class containing common methods used in thermochemistry calculations.

static combine_contributions(contrib_dicts: Sequence[Dict[str, float]]) → Dict[str, float]

Combine the contributions from multiple modes.

classmethod from_transition_state(vib_energies: Sequence[complex], *args, **kwargs) → BaseThermoChem

Create a new instance for a transition state.

This will work just as the standard constructor, but will remove one imaginary frequency from the given vib_energies first. If there is more than one imaginary frequency, an error will be raised.

Return type:

BaseThermoChem instance

get_ZPE_correction() → float

Returns the zero-point vibrational energy correction in eV.

get_ideal_entropy(temperature: float, translation: bool = False, vibration: bool = False, rotation: bool = False, geometry: Literal['linear', 'nonlinear', 'monatomic'] | None = None, electronic: bool = False, pressure: float | None = None, symmetrynumber: int | None = None) → Tuple[float, Dict[str, float]]

Returns the entropy, in eV/K and a dict of the contributions.

Parameters:

• temperature (float) – The temperature in Kelvin.

• translation (bool) – Include translational entropy.

• vibration (bool) – Include vibrational entropy.

• rotation (bool) – Include rotational entropy.

• geometry (str) – The geometry of the molecule. Options are ‘nonlinear’, ‘linear’, and ‘monatomic’.

• electronic (bool) – Include electronic entropy.

• pressure (float) – The pressure in Pa. Only needed for the translational entropy.

• symmetrynumber (int) – The symmetry number of the molecule. Only needed for linear and nonlinear molecules.

Returns:

• Tuple of one float and one dict

• The float is the total entropy in eV/K.

• The dict contains the contributions to the entropy.

static get_ideal_rotational_energy(geometry: Literal['linear', 'nonlinear', 'monatomic'], temperature: float) → float

Returns the rotational heat capacity times T in eV.

Parameters:

• geometry (str) – The geometry of the molecule. Options are ‘nonlinear’, ‘linear’, and ‘monatomic’.

• temperature (float) – The temperature in Kelvin.

Return type:

float

get_ideal_trans_entropy(atoms: Atoms, temperature: float) → float

Returns the translational entropy in eV/K.

Parameters:

• atoms (ase.Atoms) – The atoms object.

• temperature (float) – The temperature in Kelvin.

Return type:

float

static get_ideal_translational_energy(temperature: float) → float

Returns the translational heat capacity times T in eV.

Parameters:

temperature (float) – The temperature in Kelvin.

Return type:

float

get_vib_energy_contribution(temperature: float) → float

Calculates the change in internal energy due to vibrations from 0K to the specified temperature for a set of vibrations given in eV and a temperature given in Kelvin. Returns the energy change in eV.

get_vib_entropy_contribution(temperature: float, return_list: bool = False) → float | Sequence[float]

Calculates the entropy due to vibrations for a set of vibrations given in eV and a temperature given in Kelvin. Returns the entropy in eV/K.

print_contributions(contributions: Dict[str, float], verbose: bool) → None

Print the contributions.

property vib_energies: ndarray

Get the vibrational energies in eV.

For backwards compatibility, this will return the modes’ energies if they are available. Otherwise, it will return the stored vib_energies.

class ase.thermochemistry.AbstractMode(energy: float)

Abstract base class for mode objects.

Input:

energy: float

Initialize the mode with its energy in eV. Note: This energy refers to the vibrational energy of the mode. Get it e.g. from ase.vibrations.Vibrations.get_energies.

Individual Mode Classes

class ase.thermochemistry.HarmonicMode(energy: float)

Class for a single harmonic mode.

get_ZPE_correction() → float

Returns the zero-point vibrational energy correction in eV.

get_entropy(temperature: float, contributions: bool) → float | Tuple[float, Dict[str, float]]

Returns the entropy in the harmonic approximation at a specified temperature (K).

Parameters:

• temperature (float) – The temperature in Kelvin.

• contributions (bool) – If True, return the contributions to the entropy as a dict too. Here it will only contain the vibrational entropy.

Return type:

float or Tuple[float, Dict[str, float]]

get_internal_energy(temperature: float, contributions: bool) → float | Tuple[float, Dict[str, float]]

Returns the internal energy in the harmonic approximation at a specified temperature (K).

get_vib_energy_contribution(temperature: float) → float

Calculates the change in internal energy due to vibrations from 0K to the specified temperature for a set of vibrations given in eV and a temperature given in Kelvin. Returns the energy change in eV.

get_vib_entropy_contribution(temperature: float) → float

Calculates the entropy due to vibrations given in eV and a temperature given in Kelvin. Returns the entropy in eV/K.

class ase.thermochemistry.RRHOMode(energy: float, mean_inertia: float, tau: float = 35.0, treat_int_energy: bool = False)

Class for a single RRHO mode, including Grimme’s scaling method based on doi: 10.1002/chem.201200497 and doi: 10.1039/D1SC00621E.

Inputs:

mean_inertiafloat

The mean moment of inertia in amu*angstrom^2. Use \(np.mean(ase.Atoms.get_moments_of_inertia())\) to calculate.

taufloat

the vibrational energy threshold in \(cm^{-1}\), named \(\tau\) in doi: 10.1039/D1SC00621E. Values close or equal to 0 will result in the standard harmonic approximation. Defaults to \(35cm^{-1}\).

treat_int_energybool

Extend the msRRHO treatement to the internal thermal energy. If False, only the entropy contribution as in Grimmmes paper is modified according to the RRHO scheme. If True, the approach of Otlyotov and Minenkov doi: 10.1002/jcc.27129 is used to modify the internal energy contribution.

get_RRHO_entropy_r(temperature: float) → float

Calculates the rotation of a rigid rotor for low frequency modes.

Equation numbering from doi: 10.1002/chem.201200497

Returns the entropy contribution in eV/K.

get_entropy(temperature: float, contributions: bool) → float | Tuple[float, Dict[str, float]]

Returns the entropy in the harmonic approximation at a specified temperature (K).

Parameters:

• temperature (float) – The temperature in Kelvin.

• contributions (bool) – If True, return the contributions to the entropy as a dict too. Here it will only contain the vibrational entropy.

Return type:

float or Tuple[float, Dict[str, float]]

get_internal_energy(temperature: float, contributions: bool) → float | Tuple[float, Dict[str, float]]

Returns the internal energy in the msRRHO approximation at a specified temperature (K).

If self.treat_int_energy is True, the approach of Otlyotov and Minenkov doi: 10.1002/jcc.27129 is used. Otherwise, the approach of Grimme doi: 10.1002/chem.201200497 is used.

static get_rrho_internal_energy_r_contribution(temperature: float) → float

Calculates the rotation of a rigid rotor contribution.

Equation numbering from doi: 10.1002/jcc.27129

Returns the internal energy contribution in eV.

get_rrho_internal_energy_v_contribution(temperature: float) → float

RRHO Vibrational Internal Energy Contribution from doi: 10.1002/jcc.27129.

Returns the internal energy contribution in eV.