Calculate equilibrium chemical activities of species from the affinities of formation of the species at unit activity.
1 2 3 4 5 6
equilibrate(aout, balance = NULL, loga.balance = NULL, ispecies = !grepl("cr", aout$species$state), normalize = FALSE, as.residue = FALSE, method = c("boltzmann", "reaction"), tol = .Machine$double.eps^0.25) equil.boltzmann(Astar, n.balance, loga.balance) equil.reaction(Astar, n.balance, loga.balance, tol = .Machine$double.eps^0.25) moles(eout)
list, output from
character or numeric, how to balance the transformations
numeric, which species to include
logical, normalize the molar formulas of species by the balancing coefficients?
logical, report results for the normalized formulas?
numeric, affinities of formation reactions excluding species contribution
numeric, number of moles of balancing component in the formation reactions of the species of interest
numeric (single value or vector), logarithm of total activity of balanced quantity
character, equilibration method to use
numeric, convergence tolerance for
list, output from
equilibrate calculates the chemical activities of species in metastable equilibrium, for constant temperature, pressure and chemical activities of basis species, using specified balancing constraints on reactions between species.
It takes as input
aout, the output from
affinity, giving the chemical affinities of formation reaction of each species, which may be calculated on a multidimensional grid of conditions.
aout can be the output from
mosaic, in which case the equilibrium activities of the formed species are calculated and combined with those of the changing basis species to make an object that can be plotted with
The equilibrium chemical activities of species are calculated using either the
equil.boltzmann functions, the latter only if the balance is on one mole of species.
equilibrate needs to be provided constraints on how to balance the reactions representing transformations between the species.
balance indicates the balancing component, according to the following scheme:
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 use the coefficients on the basis species that is present (i.e. with non-zero coefficients) in all formation reactions, 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 equivalent to ‘activity’ here.
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).
loga.balance is missing, its value is taken from the activities of species listed in
aout; this default is usually the desired operation.
The supplied value of
loga.balance may also be a vector of values, with length corresponding to the number of conditions in the calculations of
normalize if TRUE indicates to normalize the molar formulas of species by the balance coefficients.
This operation is intended for systems of 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.
After normalization and equilibration, the equilibrium activities are then re-scaled (for the original formulas of the species), unless
as.residue is TRUE.
equil.boltzmann is used to calculate the equilibrium activities if
balance is 1 (or when
as.residue is TRUE), otherwise
equil.reaction is called.
The default behavior can be overriden by specifying either boltzmann or reaction in
equil.reaction may be needed for systems with huge (negative or positive) affinities, where
equil.boltzmann produces a NaN result.
ispecies can be supplied to identify a subset of the species to include in the equilibrium calculation.
By default, this is all species except solids (species with cr state).
However, the stability regions of solids are still calculated (by a call to
diagram without plotting).
At all points outside of their stability region, the logarithms of activities of solids are set to -999.
Likewise, where any solid species is calculated to be stable, the logarithms of activities of all aqueous species are set to -999.
moles simply calculates the total number of moles of elements corresponding to the activities of formed species in the output from
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.equil (the calculated equilibrium activities of the species).
The input values to
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 base-10 logarithm activities 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 balancing component is equal to (
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.
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.
The solution is found using
uniroot with a flexible method for generating initial guesses.
equil.boltzmann, 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.
equil.boltzmann is called with
balance other than 1, it stops with an error.
Despite its name, this function does not generally produce a complete equilibrium.
It returns activities of species such that the affinities of formation reactions are equal to each other (and transformations between species have zero affinity); this is a type of metastable equilibrium.
Although they are equal to each other, the affinities are not necessarily equal to zero.
solubility to find complete equilibrium, where the affinities of the formation reactions become zero.
Dick, J. M. (2008) Calculation of the relative metastabilities of proteins using the CHNOSZ software package. Geochem. Trans. 9:10. https://doi.org/10.1186/1467-4866-9-10
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.boltzmann to parallelize intensive parts of the calculations.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
## 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 balancing component (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), rep(0.001, 100)) # probably close enough for most purposes stopifnot(all.equal(loga.boltz, loga.react)) # the first ionization constant (pKa) ipKa <- which.min(abs(loga.boltz[] - loga.boltz[])) 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)) # equilibrate with mosaic: balancing on two elements (N and C) # thanks to Kirt Robinson for the feature request and test system loga_N <- -4 loga_C <- -3 basis(c("CO2", "NH3", "O2", "H2O", "H+")) basis("NH3", loga_N) species(c("acetamide", "acetic acid", "acetate")) # this calculates equilibrium activities of NH3 and NH4+ for given loga_N # and calculates the corresponding affinities of the formed species m <- mosaic(c("NH3", "NH4+"), pH = c(0, 14)) # this calculates equilibrium activities of the formed species for given loga_C # and combines them with the activities of the changing basis species (NH3 and NH4+) eqc <- equilibrate(m, loga.balance = loga_C) diagram(eqc, ylim = c(-10, -2)) title(main = paste("log(total N in basis species) =", loga_N, "\nlog(total C in formed species) =", loga_C), font.main = 1)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.