Chemical reactions and chemical species

This language extension implements the domain concept of chemical reactions. The language enables defining chemical species and reactions, and computing thermochemical and electrochemical properties for the defined species and reactions.

Chemical species

The syntax for a species definition is:

Species <name> [, composition: <composition>] [table with parameters and properties]

The word Species is a keyword and the species’ name is mandatory. The specifications of composition and the properties table are optional. Here an example for water:

water = Species H2O

The composition is either a quoted string like 'H2O' or a reference to any parameter of string type. Here an example for water with specified composition:

w = Species water, composition: 'H2O'

The composition string can be any brute chemical formula formatting, e.g. ‘OH2’ a ‘HHO’ are accepted.

Currently, the optional table can include the following columns: energy, enthalpy, entropy, free_energy, zpe, and temperature. For example:

water = Species H2O (
    (energy: -14.22696999) [electron_volt],
    (zpe: 0.558) [eV],
    (entropy: 2.1669e-3) [eV/K],
    (temperature: 298.15) [K]
)

The provided input allows to evaluate the free energy at the given temperature. The free energy can be retrieved with water.free_energy. It is noted that water.free_energy is a series with the length of the series provided in the table in the species definition. To get the first element, water.free_energy[0] has to be used. Here an example with two temperatures:

water = Species H2O (
    (energy: -14.22696999, -14.22696999) [electron_volt / particle],
    (zpe: 0.558, 0.558) [eV / particle],
    (entropy: 206.5, 213.1) [J/K/mol],
    (temperature: 500., 600.) [K]
)
print(water.free_energy)
print('temperature, free energy:', water.temperature[0], water.free_energy[0])

In the latter example, water.free_energy includes two elements, for 500 and 600 K, respectively:

program output: >>>
(free_energy: -14.739080832009071, -14.994145508249682) [electron_volt / particle]
'temperature, free energy:', 500.0 [kelvin], -14.739080832009071 [electron_volt / particle]
<<<

The interpreter evaluates all thermochemistry quantities that can be determined by using the available data. For example, if the free energy is given then the energy and the enthalpy will be evaluated:

water = Species H2O, composition: 'H2O' (
    (free_energy: -14.739080832009071, -14.994145508249682) [electron_volt / particle],
    (zpe: 0.558, 0.558) [eV / particle],
    (entropy: 206.5, 213.1) [J/K/mol],
    (temperature: 500., 600.) [K]
)
print(water.energy, water.enthalpy)

Output:

program output: >>>
(energy: -14.22696999, -14.22696999) [electron_volt / particle] (enthalpy: -13.66896999, -13.66896999) [electron_volt / particle]
<<<

Chemical reactions

An equation of a chemical reaction is defined with the syntax:

Reaction <reactants> [=|->] <products> [: table with parameters and properties]

The word Reaction is a keyword. The individual reactants and products must be separated by a + sign and must be references to the variables defining the relevant species. Every reactant (or product) has a name, that must be the name of a Species, and a stoichiometric coefficient that is a real number. If no coefficient is specified, then 1.0 is assumed. Example:

h = Species H2; o2 = Species O2; h2o = Species H2O
r = Reaction 2 H2 + O2 = 2 H2O: ((free_energy: -4.916) [eV])

Thermochemical quantities for a reaction will be evaluated if the necessary data about the species is available. These quantities can be retrieved in the same way as for species. Example:

react = Reaction 2 H2 + O2 = 2 H2O
h2o = Species H2O (
    (free_energy: -2.458) [eV],
    (entropy: 2.1669e-3) [eV/K],
    (zpe: 0.558) [eV],
    (temperature: 298.15) [K]
)
h2 = Species H2 (
    (free_energy: 0.0) [eV],
    (zpe: 0.270) [eV],
    (entropy: 1.3613e-3) [eV/K],
    (temperature: 298.15) [K]
)
o2 = Species O2 (
    (free_energy: 0.0) [eV],
    (entropy: 2.1370e-3) [eV/K],
    (zpe: 0.098) [eV],
    (temperature: 298.15) [K]
)
print(react.free_energy)
print(react.enthalpy)
print(react.entropy)

program output: >>>
(free_energy: -4.916) [electron_volt]
(enthalpy: -5.07276727) [electron_volt]
(entropy: -0.0005258000000000007) [electron_volt / kelvin]
<<<

The calculation of thermochemical quantities assumes that the reaction equation is balanced. If all terms (species) in the equation have compositions specified then a check is performed and if the equation is not balanced the evaluation is interrupted with an error. Therefore, it is recommended to specify species compositions to enable these checks. If some compositions are not available for the check, a warning about this is issued.

Retrieve all properties of Species or Reaction

All properties of a Species or a Reaction parameter can be retrieved as Table using the properties keyword, for example:

print(H2O.properties)
print(H2O.composition)

program output: >>>
((free_energy: -2.458) [electron_volt], (entropy: 2.1669e-3) [electron_volt / kelvin], (zpe: 0.558) [electron_volt], (temperature: 298.15) [kelvin])
'H2O'
<<<

Using Species and Reaction as iterables

The Species and Reaction parameters are table-like iterables and can be used in any operations that are defined for the Table type. These operations are, e.g., subscripting, slicing, filter function and expression, map and reduce, and base on the table with properties.

Examples:

H2O = Species H2O (
    (free_energy: -14.739080832009071, -14.994145508249682) [eV],
    (temperature: 500., 600.) [K]
)

H2O_500K = H2O select free_energy where temperature == 500. [K]
H2O_500K_alt = filter((x: x.temperature == 500.0 [K]), H2O)
H2O_500K_free_table = map((x: {free_energy: x.free_energy}), H2O_500K)
H2O_500K_free_series = H2O_500K.free_energy

MO = Species MO
M = Species M
H2 = Species H2

react = Reaction MO + H2 = M + H2O:
    ((free_energy: -1., -2.) [eV],
     (temperature: 500., 600.) [K])

max(x, y) = if(y > x, y, x)
max_alt(x, y, ax, ay) = if(y > x, ay, ax)

react_550K = react where temperature > 550 [kelvin]
react_550K_alt = filter((x: x.temperature > 550 [K]), react)
react_max_temp = reduce((x, y: {temperature: max(x.temperature, y.temperature)}), react)
free_max_temp = reduce((x, y: {free_energy: max_alt(x.temperature, y.temperature, x.free_energy, y.free_energy)}), react);

Visualize reaction energy profiles

Waterfall type of reaction energy profiles can be visualized with the following statement:

view waterfall (<series of reactions>[, <'orr' | 'oer'>])

The first parameter is a series of reactions. The second optional string parameter indicates which reaction will be displayed - oxygen reduction reaction ('orr'), which is the default, or the oxygen evolution reaction ('oer'). A full example is provided in the script TextM-oxycat_electrochem.vm in the VRE Language software repository.