ccu.thermo package¶

Convenience functions for calculating thermodynamic quantities.

For a limited number of molecular species, ccu.thermo.chempot provides vibrational zero-point energies and a DFT-parametrized interpolation for chemical potentials between 0 and 1100 K. The parametrization is of the form:

(1)¶\[ \Delta \mu_{0 \rightarrow T} = b + m T + k_b T \log P\]

where \(T\) and \(P\) are temperature (in K) and pressure (in bar), respectively, \(k_b\) is the Boltzmann constant, and \(b(T)\) and \(m(T)\) are parametrized functions of \(T\) that are unique to each molecule. These two datasets allow one to calculate free energies at any given temperature between 1 and 1000K. That is, given a chemical potential \(\Delta \mu\) and vibrational zero-point energy \(ZPE\), the free energy \(G\) is given by

(2)¶\[ G = \Delta \mu_{0 \rightarrow T} + ZPE\]

For all other species, ccu.thermo.vibration provides a standardized, DFT-code-agnostic workflow that can be used to approximate the Hessian in the ground state, ccu.thermo.vibration.run_vibration(). This function is a thin wrapper around the ase thermochemistry utilities that conveniently logs all relevant thermochemistry data as well as other information (e.g., forces, frozen atoms, frequencies) to a file.

After running ccu.thermo.vibration.run_vibration(), one can then run ccu.thermo.gibbs.calculate_free_energy() to obtain the entropic correction (\(-TS\)) and the vibrational zero-point energy (\(ZPE\)). The Gibbs free energy is then calculated as:

(3)¶\[ G = E - TS + ZPE\]

where \(E\) is the DFT-calculated energy from ccu.thermo.vibration.run_vibration().

In both of the above cases, one has the option to use the CLI (ccu thermo chempot and ccu thermo gibbs) or the Python API (ccu.thermo.chempot.calculate() or ccu.thermo.vibration.run_vibration() with ccu.thermo.gibbs.calculate_free_energy()).

Submodules¶

ccu.thermo.chempot module¶

Calculate thermodynamic properties via a DFT-parametrized calculator.

This module defines classes for the calculation of thermodynamic data via DFT-parametrized models. In this way, zero-point energies, chemical potentials, and Gibbs free energies can be calculated without the need of expensive vibrational calculations. That hard work has already been done!

In particular, this module defines two main classes:

ChemPotDatabase: a database for parametrization data.
ChemPotCalculator: a calculator for computing thermodynamic data from parametrization data in a ChemPotDatabase.

Additionally, ccu ships with a pre-defined dataset for a number of molecules. The dataset was produced from a set of DFT-D3 calculations and includes zero-point energies calculated under the rigid rotator and harmonic oscillator approximations. Calculated values are valid below 1100K. Chemical potentials are calculated via calculate_chemical_potential() using equation (1).

Examples

List all molecules included in the pre-defined dataset.

>>> from ccu.thermo.chempot import ChemPotDatabase
>>> from ccu.thermo.chempot import load_zpe_data
>>> zpe_data = load_zpe_data()
>>> for molecule in zpe_data:
...     print(molecule)
CH4
...

Retrieve parametrization data for CO₂ valid at 350 K.

>>> from ccu.thermo.chempot import ChemPotDatabase
>>> database = ChemPotDatabase()
>>> query = {"molecule": "CO2", "temperature": 350}
>>> database.parameter_data.get(**query)[0]
0.130...

Retrieve the vibrational zero-point energy of N₂.

>>> from ccu.thermo.chempot import ChemPotDatabase
>>> database = ChemPotDatabase()
>>> database.zpe_data["N2"]
0.15...

Calculate the chemical potential of H₂O at IUPAC STP (273.15 K, 1 bar).

>>> from ccu.thermo.chempot import ChemPotCalculator
>>> calc = ChemPotCalculator()
>>> chem_pot, _ = calc.calculate("H2O")
>>> chem_pot
0.114...

or simply,

>>> from ccu.thermo.chempot import calculate
>>> calculate("H2O")[0]
0.114...

Calculate the gibbs free energy of CH₄ at 900 K and 2.4 bar.

>>> from ccu.thermo.chempot import calculate
>>> sum(calculate("CH4", 900, 2.4))
-0.759...

Use a custom database to perform the above calculations.

>>> from ccu.thermo.chempot import ChemPotCalculator
>>> from ccu.thermo.chempot import ChemPotDatabase
>>> from ccu.thermo.chempot import load_zpe_data
>>> # specify the paths to your database files
>>> parameter_data = load_parameter_data(..., external=True)
>>> zpe_data = load_zpe_data(...)
>>> # load the database into the calculator
>>> database = ChemPotDatabase(
...     parameter_data=parameter_data, zpe_data=zpe_data
... )
>>> calc = ChemPotCalculator(database)
>>> # and then as above...

class ccu.thermo.chempot.ChemPotCalculator(database: ChemPotDatabase | None = None)[source]¶

Bases: object

A chemical potential calculator.

Variables:: database – The underlying database used for calculations.

Example

>>> from ccu.thermo.chempot import ChemPotCalculator
>>>
>>> calc = ChemPotCalculator()
>>> calc.calculate("CO2", 298.15, 1.0)
-0.581856

Initialize the calculator from a database.

Parameters:: database – The database to use for the calculator. Defaults to None, in which case the pre-defined database is used.

calculate(molecule: str, temperature: float, pressure: float) → tuple[float, float][source]¶

Calculate the chemical potential of a molecule.

Parameters:

molecule – The molecule for which to calculate the chemical potential.
temperature – The temperature at which to calculate the chemical potential.
pressure – The pressure at which to calculated the chemical potential.

Returns:

A 2-tuple (chem_pot, zpe) whose elements represent the calculated chemical potential and retrieved vibrational zero-point energy, respectively.

property max_temperature: float¶

The maximum temperature supported by the underlying database.

Parameters:: tol – The amount by which to decrement the maximum reported upper limit in the underlying database.

Note

This property represents the maximum temperature of any data point irrespective of the molecule or parameter (“b” or “m”). Further, since the upper limit of an instance of ChemPotDataPoint represents an exclusive upper limit, the maximum reported upper limit is decremented by a tol to ensure non-empty queries with this value.

class ccu.thermo.chempot.ChemPotDataPoint(lower: float, molecule: str, param: Literal['b', 'm'], upper: float, value: float)[source]¶

Bases: NamedTuple

A thermodynamic parameter, valid in a certain temperature range.

Variables:

lower (float) – The lower temperature bound of the range (inclusive) in which the data point is valid.
molecule (str) – The molecule to which that data point applies.
term – The kind of term represented by the ChemPotDataPoint. If "b", the data point corresonds to a y-intercept. If "m", then the data point corresponds to a slope-like term.
upper (float) – The upper temperature bound of the range (exclusive) in which the data point is valid.
value (float) – The value of the parameter.

Create new instance of ChemPotDataPoint(lower, molecule, param, upper, value)

lower: float¶: Alias for field number 0

molecule: str¶: Alias for field number 1

param: Literal['b', 'm']¶: Alias for field number 2

upper: float¶: Alias for field number 3

value: float¶: Alias for field number 4

class ccu.thermo.chempot.ChemPotDatabase(parameter_data: list[~ccu.thermo.chempot.ChemPotDataPoint] = <factory>, zpe_data: dict[str, float] = <factory>)[source]¶

Bases: object

A database of parametrized thermochemical data.

Variables:

parameter_data (list[ccu.thermo.chempot.ChemPotDataPoint]) – A list of ChemPotDataPoint representing thermodynamic data for a parametrization.
zpe_data (dict[str, float]) – A dictionary mapping molecule names to vibrational zero-point energies.

Example

Retrieve b parameter data for CO₂ that is valid at 350 K.

>>> from ccu.thermo.chempot import ChemPotDatabase
>>> database = ChemPotDatabase()
>>> # Note that the following returns a list
>>> database.get(molecule="CO2", temperature=350, param="b")
[ChemPotDataPoint(param='b', lower=350.0, upper=400.0, molecule='CO...

Example

Retrieve ZPE data for CO₂.

>>> from ccu.thermo.chempot import ChemPotDatabase
>>> database = ChemPotDatabase()
>>> database.zpe_data["CO2"]
0.306

get(temperature: float | None = None, **query) → list[ChemPotDataPoint][source]¶

Retrieve values from the parametrization database.

Parameters:

temperature – The temperature to use as a search key. Defaults to None. Note, however, that specifying this instead results in all database entries with lower and upper values which bracket temperature to be used. If values with a specific lower or upper range are desired, use the keys "lower" and "upper", respectively.
**query – Specify filtering criteria with which to search for entries. Any keys corresponding to fields in a ChemPotDataPoint are valid.

Raises:

LookupError – You specified a key that is not a field of ChemPotDataPoint.

Returns:

A list of every ChemDataPoint with fields matching your query.

Note

Although it doesn’t affect the results returned, the order that keyword arguments are specied affects the speed in which results are returned. Except for temperature, filtering criteria is applied in the order it is supplied in the function call. If supplied, the temperature criteria is always applied first.

parameter_data: list[ChemPotDataPoint]¶

zpe_data: dict[str, float]¶

ccu.thermo.chempot.calculate(molecule: str, temperature: float = 273.15, pressure: float = 1.0) → tuple[float, float][source]¶

Calculate molecular thermodynamic properties with the default data.

Parameters:

molecule – The name of the molecule for which the calculation will be performed. Must be one of ccu.thermo.chempot.MOLECULES.
temperature – The temperature (in Kelvin) at which to perform the calculation. Defaults to 298.15.
pressure – The pressure (in bar) at which to perform the calculation. Defaults to 1.0.

Returns:

A 2-tuple (chem_pot, zpe) corresponding to the calculated chemical potential and vibrational zero-point energy, respectively.

Note

The model is only parametrized up to 1000K. Although values up to 1100K may still be reliable, at temperatures higher than 1100K, results become unreliable.

ccu.thermo.chempot.calculate_chemical_potential(b: float, m: float, temperature: float, pressure: float) → float[source]¶

Calculate the chemical potential from parametrized model.

This equation is consistent with (1).

Parameters:

b – The y-intercept.
m – The “slope”-like term in the interpolation.
temperature – The temperature at which to calculate the chemical potential.
pressure – The pressure at which to calculate the chemical potential.

Returns:

The chemical potential at the conditions specified.

ccu.thermo.chempot.load_parameter_data(*, filename: str | Path = 'param_data.csv', external: bool = False) → list[ChemPotDataPoint][source]¶

Load parametrized data from a file of comma-separated values.

Parameters:

filename – The file from which to load the data. Defaults to “param_data.csv”.
external – Whether to load the file from in external source. Defaults to False.

Returns:

A dictionary mapping molecule names to instances of ChemPotDataset.

Note

The fields of the CSV file must include those of a ChemPotDataPoint, and the values for the keys “lower”, “upper”, and “value” must be able to be converted into floats.

ccu.thermo.chempot.load_zpe_data(*, filename: str | Path = 'zpe_data.csv', external: bool = False) → dict[str, float][source]¶

Load vibrational zero-point energy data from a CSV file.

Parameters:

filename – The file from which to load the data. Defaults to “zpe_data.csv”.
external – Whether to load the file from an external source. Defaults to False.

Returns:

A dictionary mapping molecule names to their vibrational zero-point energies.

Note

The fields of the CSV file must include “molecule” and “zpe”, and the values corresponding to “zpe” must be able to be converted into floats.

ccu.thermo.chempot.print_instructions() → None[source]¶: Print instructions for using chempot interactively.

ccu.thermo.cli module¶

CLI utilities for thermochemistry.

ccu.thermo.cli.print_molecules(ctx, value, param: Parameter) → None[source]¶: Name all species for which chempot is parameterized.

ccu.thermo.cli.report_state(ctx: Context) → None[source]¶: Report the value and source of CLI parameters.

ccu.thermo.gibbs module¶

Functions for calculating Gibbs free energies from vibrational data.

ccu.thermo.gibbs.calculate_free_energy(log_file: TextIO | None = None, vib_file: TextIO | None = None, approximation: Literal['IDEAL_GAS', 'HARMONIC'] = 'HARMONIC', symmetry: int = 1, geometry: Literal['linear', 'nonlinear'] = 'nonlinear', transition_state: bool = False, frequency_threshold: float = 12, temperature: float = 273.15, pressure: float = 1.0, spin: int = 0, atoms_file: str = 'in.traj') → tuple[float, float, list[float]][source]¶

Calculate the thermodynamic corrections in the specified approximation.

Parameters:

log_file – A file in which to log the results.
vib_file – A file from which to read the frequencies.
approximation – Which approximation to use. One of "IDEAL_GAS" or "HARMONIC". Defaults to "HARMONIC".
symmetry – The symmetry number of the system. Defaults to 1.
geometry – The geometry of the system. One of "linear" or "nonlinear" Defaults to "nonlinear".
transition_state – Whether to treat the system as a transition state. Defaults to False.
frequency_threshold – All frequencies below this value with be normalized to this value. Defaults to 12.
temperature – The temperature for the system (in Kelvin). Note that this is only relevant if approximation="IDEAL_GAS". Defaults to 273.15.
pressure – The pressure for the system (in bar). Note that this is only relevant if approximation="IDEAL_GAS". Defaults to 1.0.
spin – The spin for the system. Note that this is only relevant if approximation="IDEAL_GAS". Defaults to 0.
atoms_file – A string pointing to the structure to which the vibrational calculations correspond. Defaults to "in.traj".

Returns:

A 3-tuple (ts, zpe, freq) where ts is the entropic correction (\(-TS\)), zpe is the vibrational zero-point energy (\(ZPE\)), and freq is a list of floats representing the frequencies used to calculate the zero-point energy.

Note

The free energy of the system can be calculated using (3)

ccu.thermo.gibbs.select_frequencies(vib_file: TextIO | None = None, geometry: Literal['linear', 'nonlinear'] = 'nonlinear', approximation: Literal['IDEAL_GAS', 'HARMONIC'] = 'HARMONIC', transition_state: bool = False, frequency_threshold: float = 12) → list[float][source]¶

Discard frequencies according to the geometry and approximations.

Parameters:

vib_file – An opened text file in which to save the vibrational data.
geometry – A string indicating the geometry. One of "linear" or "nonlinear". Defaults to "nonlinear".
approximation – A string indicating the approximation used in the treatment of the free energy. Defaults to "HARMONIC".
transition_state – Whether the frequencies correspond to a transition state. Defaults to False.
frequency_threshold – All frequencies below this value will be normalized to this value. Defaults to 12 cm^-1.

Returns:

The appropriate frequencies for the system.

ccu.thermo.vibration module¶

Standardized functions for vibrational calculations.

ccu.thermo.vibration.run_vibration(atoms: Atoms, *, label='vib') → None[source]¶

Run a vibrational calculation and log the output.

The equilibrium structure is saved as final.traj

Parameters:

atoms – An atoms.Atoms object with an attached calculator with which to run the vibration calculation.
label – A string used to name the summary and trajectory files. Defaults to "vib". If the default value is used, the output will be logged in vib.txt.

Raises:

RuntimeError – Unable to perform vibrational calculation due to missing forces.