# 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](scl.md#string-type) like `'H2O'` or a reference to any parameter of [string type](scl.md#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`](#chemical-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](#chemical-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](scl.md#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](scl.md#series) of [reactions](#chemical-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.