# Equilibrium Chemical Activities of Species

### Description

Calculate equilibrium chemical activities of species from the affinities of formation of the species at unit activity.

### Usage

1 2 3 4 5 6 |

### Arguments

`aout` |
list, output from |

`balance` |
character or numeric, how to balance the transformations |

`ispecies` |
numeric, which species to include |

`normalize` |
logical, normalize the molar formulas of species by the balancing coefficients? |

`as.residue` |
logical, report results for the normalized formulas? |

`Astar` |
numeric, affinities of formation reactions excluding species contribution |

`n.balance` |
numeric, number of moles of conserved component in the formation reactions of the species of interest |

`loga.balance` |
numeric, logarithm of total activity of balanced quantity |

`method` |
character, equilibration method to use |

### Details

`equilibrate`

provides an interface to calculate the chemical activities of species in metastable equilibrium, in an open system at constant temperature and pressure and chemical activities of basis species, and with linear balancing constraints on transformations.

It takes as input `aout`

, the output from `affinity`

, which may be calculated from a multidimensional grid of conditions.
The equilibrium chemical activities of species are calculated using either the `equil.reaction`

or `equil.boltzmann`

functions, the latter only if the balance is on one mole of species.

As `aout`

contains the chemical affinities of formation reactions of each species of interest, `equilibrate`

in order to function needs to be provided constraints on how to balance the reactions representing transformations between the species.
`balance`

returns the balancing coefficients, where `balance`

indicates the balancing constraints, according to the following scheme:

NULL | autoselect using `which.balance` |

name of basis species | balance on this basis species |

length | balance on length of proteins |

1 | balance on one mole of species |

numeric vector | user-defined constraints |

The default value of NULL for `balance`

indicates to select the first shared basis species in all formation reactions identified using `which.balance`

, or if that fails, to set the balance to 1.
However, if all the species (as listed in code `aout$species`

) are proteins (have an underscore character in their names), the default value of NULL for `balance`

indicates to use length as the balance.

NOTE: the summation of activities assumes an ideal system, so ‘molality’ is implied by ‘activity’ in the following.
`loga.balance`

gives the logarithm of the total activity of `balance`

(which is total activity of species for 1 or total activity of amino acid residue-equivalents for length).
If `loga.balance`

is missing, its value is taken from the activities of species listed in `aout`

; this default is usually the desired operation.

`normalize`

if TRUE indicates to normalize the molar formulas of species by the balance coefficients.
This operation is intended for systems of polymers, such as proteins, whose conventional formulas are much larger than the basis speices.
The normalization also applies to the balancing coefficients, which as a result consist of 1s.
`normalize`

has the same effect as did `diagram(..., residue=TRUE)`

in versions of CHNOSZ before 0.9-9.
After normalization and equilibration, the equilibrium activities are then re-scaled (for the original formulas of the species), unless `as.residue`

is TRUE.

`ispecies`

can be supplied to identify a subset of the species to include in the calculation.

`equil.boltzmann`

is used to calculate the equilibrium activities if `balance`

is 1 (or when `normalize`

or `as.residue`

is TRUE), otherwise `equil.reaction`

is called.
The default behavior can be overriden by specifying either boltzmann or reaction in `method`

.
Using `equil.reaction`

may be needed for systems with huge (negative or positive) affinities, where `equil.boltzmann`

produces a NaN result.

### Value

`equil.reaction`

and `equil.boltzmann`

each return a list with dimensions and length equal to those of `Astar`

, giving the `log10`

of the equilibrium activities of the species of interest.
`equilibrate`

returns a list, containing first the values in `aout`

, to which are appended `m.balance`

(the balancing coefficients if `normalize`

is TRUE, a vector of 1s otherwise), `n.balance`

(the balancing coefficients if `normalize`

is FALSE, a vector of 1s otherwise), `loga.balance`

, `Astar`

, and `loga.equil`

(the calculated equilibrium activities of the species).
`balance`

returns a list containing the balancing coefficients (`n`

) and a textual description (`description`

).

### Algorithms

The input values to `equil.reaction`

and `equil.boltzmann`

are in a list, `Astar`

, all elements of the list having the same dimensions; they can be vectors, matrices, or higher-dimensionsal arrays.
Each list element contains the chemical affinities of the formation reactions of one of the species of interest (in dimensionless base-10 units, i.e. A/2.303RT), calculated at unit activity of the species of interest.
The equilibrium activities (in base-10 log units) of the species of interest returned by either function satisfy the constraints that 1) the final chemical affinities of the formation reactions of the species are all equal and 2) the total activity of the conserved component is equal to (`loga.balance`

).
The first constraint does *not* impose a complete equilibrium, where the affinities of the formation reactions are all equal to zero, but allows for a metastable equilibrium, where the affinities of the formation reactions are equal to each other.

In `equil.reaction`

(the algorithm described in Dick, 2008 and the only one available prior to CHNOSZ-0.8), the calculations of relative abundances of species are based on a solving a system of equations representing the two constraints stated above.
A close approximation of the solution is found using `uniroot`

.
Prior to CHNOSZ_0.9-9, the values in the `Astar`

were used to generate initial guesses of the logarithms of activities of species; values of `loga.balance`

that were hugely different from these guesses could generate errors from `uniroot`

such as "f() values at end points not of opposite sign".
Now (from version 0.9-9), a more flexible method for generating guesses is in place.

In `equil.boltzmann`

(algorithm available beginning with CHNOSZ-0.8), the chemical activities of species are calculated using the Boltzmann distribution.
This calculation is faster than the algorithm of `equil.reaction`

, but is limited to systems where the transformations are all balanced on one mole of species.
If `equil.boltzmann`

is called with `balance`

other than 1, it stops with an error.

### References

Dick, J. M. (2008) Calculation of the relative metastabilities of proteins using the CHNOSZ software package. *Geochem. Trans.* **9**:10. http://dx.doi.org/10.1186/1467-4866-9-10

### See Also

`diagram`

has examples of using `equilibrate`

to make equilibrium activity diagrams. `revisit`

can be used to perform further analysis of the equilibrium activities.
`palply`

is used by both `equil.reaction`

and `equil.boltzmann`

to parallelize intensive parts of the calculations if parallel is loaded.

### Examples

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 | ```
## equilibrium in a simple system:
## ionization of carbonic acid
basis("CHNOS+")
species(c("CO2", "HCO3-", "CO3-2"))
# set unit activity of the species (0 = log10(1))
species(1:3, 0)
# calculate Astar (for unit activity)
res <- 100
Astar <- affinity(pH=c(0, 14, res))$values
# the logarithms of activity for a total activity
# of the balanced quantity (C or CO2) equal to 0.001
loga.boltz <- equil.boltzmann(Astar, c(1, 1, 1), 0.001)
# calculated another way
loga.react <- equil.reaction(Astar, c(1, 1, 1), 0.001)
# probably close enough for most purposes
stopifnot(all.equal(loga.boltz, loga.react))
# the first ionization constant (pKa)
ipKa <- which.min(abs(loga.boltz[[1]] - loga.boltz[[2]]))
pKa.equil <- seq(0, 14, length.out=res)[ipKa]
# calculate logK directly
logK <- subcrt(c("CO2","H2O","HCO3-","H+"), c(-1, -1, 1, 1), T=25)$out$logK
# we could decrease tolerance here by increasing res
stopifnot(all.equal(pKa.equil, -logK, tolerance=1e-2))
``` |