# wjd: Gibbs Energy Minimization by Steepest Descent In CHNOSZ: Thermodynamic Calculations for Geobiochemistry

## Description

Find the quantities of chemical species, subject to constant elemental bulk composition of the system, that minimize the Gibbs energy of the system.

## Usage

 ``` 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``` ``` wjd( A = matrix(c( 1,2,2,0,0,1,0,0,0,1, 0,0,0,1,2,1,1,0,0,0, 0,0,1,0,0,0,1,1,2,1),ncol=3, dimnames=list(NULL,c("H","N","O"))), G0.RT = c( -10.021,-21.096,-37.986,-9.846,-28.653, -18.918,-28.032,-14.640,-30.594,-26.111), Y = c(0.1,0.35,0.5,0.1,0.35,0.1,0.1,0.1,0.1,0.1), P = 51, nlambda = 101, imax = 10, Gfrac = 1e-7 ) element.potentials(w, plot.it=FALSE, iplot=1:ncol(w\$A)) is.near.equil(w, tol=0.01, quiet=FALSE) guess( A = matrix(c( 1,2,2,0,0,1,0,0,0,1, 0,0,0,1,2,1,1,0,0,0, 0,0,1,0,0,0,1,1,2,1),ncol=3, dimnames=list(NULL,c("H","N","O"))), B = c(2,1,1), method="stoich", minX = 0.001, iguess = 1, ic = NULL ) run.wjd(ispecies, B = NULL, method = "stoich", Y = run.guess(ispecies, B, method), P=1, T=25, nlambda=101, imax = 10, Gfrac = 1e-7, tol = 0.01) run.guess(ispecies, B = NULL, method = "stoich", iguess = NULL) equil.potentials(w, tol=0.01, T=25) ```

## Arguments

 `A` matrix, chemical formulas of the species (elements on columns) `G0.RT` numeric, the Gibbs energies / RT, at a single temperature (length equal to number of species) `Y` numeric, initial solution, a positive set of values (numbers of moles, length equal to number of species) `P` numeric, pressure in atmospheres `nlambda` numeric, number of values of fractional distance change (lambda) tested at each step. `imax` numeric, maximum number of iterations `Gfrac` numeric, Gibbs energy change of system, as a fraction of total system energy in previous step, below which iterations will stop `w` list, output from `wjd` `plot.it` logical, make a plot? `iplot` numeric, which elements for which to make plots `tol` numeric, maximum difference in chemical potentials that counts as equilibrium `quiet` logical, don't output messages? `B` numeric, numbers of moles of the elements `method` character, method used for generating an initial solution `minX` numeric, minimum mole number for 'central' method `iguess` numeric, which guess to return `ic` numeric, which combination(s) of variable species to use (NULL for all) `ispecies` numeric, species indices (rownumbers of `thermo`\$obigt) `T` numeric, temperature in \degC

## Details

`wjd` implements the steepest descent algorithm for Gibbs energy minimization in a closed system described by White et al., 1958. “Gibbs energy” (G) referred to here is the same as the “free energy” (F) used by White et al., 1958. `wjd` itself is independent of other functions or datasets in CHNOSZ, but `run.wjd` and `run.guess` are provided to make access to the thermodynamic database in CHNOSZ easier.

The default values of `A`, `G0.RT`, `Y` and `P` correspond to the example problem described by White et al., 1958 for gases in the H, N, O system at 3500 K. Note that for this, and for any other equilibrium problem that can be simulated using this function, the mole quantities in `Y` must all be positive numbers. Operationally, this vector defines not only the “initial solution” but also the bulk composition of the system; it is not possible to define the bulk composition using mole numbers of elements alone. The `dimnames` attribute in the default value for `A` gives the names of the elements; this attribute is not necessary for the function to operate, but is used in the examples below to help label the plots.

White et al., 1958 describe in detail the computation of the direction of steepest descent by means of Lagrange multipliers. They propose an iterative solution to the energy minimization problem, where at each step the mole numbers of species are recomputed and a new steepest descent direction calculated from there. However, the authors only give general guidelines for computing the value of lambda, that is, the fraction of the total distance the system actually moves in the direction of steepest descent from the current point at each iteration. The two constraints given for determining the value of lambda are that all mole numbers of species are positive, and the Gibbs energy of the system actually decreases (the minimum point is not passed). In the code described here, at each iteration the minimum value of lambda, not exceeding unity, that violates the first condition is computed directly (a value of one is assigned if the mole numbers remain positive through the entire range). With the default setting of `nlambda`, 101 values of lambda at even intervals from 0 to this maximum permissible value are tested for the second condition at each iteration, and the highest conforming value is selected. If a value of 0 occurs, it means that the algorithm has reached an endpoint independently of the iteration and convergence tests (`rho` and `Gfrac`; see below). If this occurs, the value of `nlambda` might have to be increased depending on the user's needs.

The number of iterations is controlled by `imax` and `Gfrac`. The maximum number of iterations is set by `imax`; it can even be zero, though such a setting might only be useful in combination with `element.potentials` to characterize the initial state of a minimization problem. Within the limit of `imax`, the iterations continue until the magnitutde of the change in total Gibbs energy of the system, as a fraction of the system's energy in the previous iteration, is lower than the value specified in `Gfrac`. For the first example below, the default setting of `Gfrac` causes the algorithm to stop after six iterations.

Using the output of `wjd`, provided in the argument `w`, `element.potentials` calculates the chemical potentials of the elements in the system. It does so by combining the values of `G0.RT` of species with the inverses of stoichiometric matrices of combinations of species whose elemental compositions are linearly independent from each other. These possible combinations are constructed using the function `invertible.combs`. The value returned by `element.potentials` is a matrix, with each column corresponding to a different element and each row to a different combination of species. The entries in the matrix are the chemical potentials of the elements divided by RT. If `plot.it` is set to TRUE, the chemical potentials of the elements are plotted as a function of species combination number, with as many plots as elements, unless `iplot` is set to another value (e.g. c(1,3) for only elements 1 and 3). In the first example below, the number of unique combinations of species is 120, but only 76 of these combinations provide stoichiometric independence.

There is no guarantee that the algorithm will converge on a global (or even be close to a local) minimum. However, some tests are available to help assess the likelihood that a solution is close to equilibrium. A necessary condition of equilibrium is that the chemical potentials of the elements be independent of the particular combination of species used to compute them. `is.near.equil` compares the chemical potentials for each element, computed using `element.potentials`, with the value of `tol`. If, for each element, the range of potentials/RT (difference between minimum and maximum) is less than `tol`, the result is TRUE, otherwise the function returns FALSE, and prints a message unless `quiet` is TRUE. The default value of `tol` corresponds to an energy of 0.01 * 1.9872 * 298.15 = ca. 6 cal/mol at 25 \degC.

One of the constraints of the algorithm coded in `wjd` is that the initial solution, and all iterations, require positive (non-zero) numbers of moles of every species. Often, when investigating an equilibrium problem, the stoichiometric constraints are expressed most readily in terms of bulk composition – numbers of moles of each element. `guess` is a function to make initial guesses about the numbers of moles of all species in the system subject to the positivity constraints. Its system-specific arguments are the stoichiometric matrix `A` (as defined above for `wjd`) and the bulk composition vector `B`, giving the number of moles of each element, in the same order that the elements appear in `A`. The first `method` available in `guess` generates the ‘central’ solution of the system of linear equations using the `xranges` function from limSolve. The central solution is the mean of ranges of unknowns. The inequality constraint, or minimum number of moles of any species, is given by `minX`.

The second method for `guess` ‘stoich’ (and the default for `run.guess` and `run.wjd`) is to test successive combinations of species whose elemental compositions are linearly independent. The linearly independent combinations tested are all those from `invertible.combs` if `ic` is NULL, or only those identified by `ic`. Each combination is referred to as the ‘variable’ species; the moles of all ‘other’ species are set to a single value. This value is determined by the constraint that the greatest proportion, relative to the bulk composition in `B`, of any element contributed by all the ‘other’ species is equal to a value in `max.frac` (see code). The function tests nine hard-coded values of `max.frac` from 0.01 to 0.99, at each one solving for the moles of the ‘variable’ species that make up the difference in numbers of moles of elements. If the numbers of moles of all the ‘variable’ species is positive, the guess is accepted. The first accepted guess is returned if `iguess` is 1 (the default); other values of `iguess` indicate which guess to return. If `iguess` is NULL, all results are returned in a list, with non-successful guesses indicated by NA. In the first example below, of the 76 combinations of species whose elemental compositions are linearly independent, 32 yield guesses following this method.

`run.wjd` is a wrapper function to run `wjd`, provided the `ispecies` in the thermodynamic database (`thermo\$obigt`), the chemical composition of the system in `B`, and the pressure `P` and temperature `T`; the latter are passed to `subcrt` (with `exceed.Ttr = TRUE`) to generate the values of `G0.RT` for `wjd`. Alternatively to `B`, the initial guess of numbers of moles of species can be provided in `Y`; otherwise as many combinations of `Y` as returned from `run.guess` are tested until one is found that `is.near.equil`. The function gives an error if none of the guesses in `Y` is near equilibrium, within the tolerance set by `tol`.

`run.guess` is a wrapper function to call `guess` using the stoichiometric matrix `A` built from the `ispecies` indices in the thermodynamic database.

`equil.potentials` returns the average (`colMeans`) of `element.potentials(w)`, or NULL if `is.near.equil(w, tol=tol)` is FALSE. The output of this function can be used as the `emu` argument for `basis.logact` to calculate the corresponding activities of the basis species.

## Value

`wjd` returns a list with the problem definition and results: elements `A`, `G0.RT`, `Y`, and `P` are as supplied in the arguments; the results are in `X` (final mole numbers of species), `F.Y` (Gibbs energy of the system at initial conditions and after each iteration), `lambda` (value used for lambda at each iteration), and `elements` (matrix of moles of elements at initial conditions and after each iteration; iterations on the columns and elements on the rows).

## References

White, W. B., Johnson, S. M. and Dantzig, G. B. (1958) Chemical equilibrium in complex mixtures. J. Chem. Phys. 28, 751–755. https://doi.org/10.1063/1.1744264

`demo("wjd")` for a longer example, and `invertible.combs`, used by `element.potentials` to find combinations of species that are compositionally independent.
 ``` 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 41 42 43 44``` ```## run the function with default settings to reproduce ## the example problem in White et al., 1958 w <- wjd() # the mole fractions are very close to those shown in the # last column of Table III in the paper print(w\$X) # the Gibbs energy of the system decreased, # and by a smaller amount, at each iteration print(diff(w\$F.Y)) # there are 76 unique combinations of species that can be # used to calculate the chemical potentials of the elements stopifnot(nrow(invertible.combs(w\$A))==76) # what the scatter in those chemical potentials looks like ep <- element.potentials(w, plot.it=TRUE) # the differences in chemical potentials / RT are all less than 0.01 is.near.equil(w) # TRUE ## run the algorithm for only one iteration w <- wjd(imax=1) # the scatter in chemical potentials is much greater ep <- element.potentials(w, plot.it=TRUE) # and we're pretty far from equilibrium is.near.equil(w) # FALSE ## test all of the guesses of inititial mole quantities ## provided by guess() using default bulk composition (H2NO) # 9 of them are not is.near.equil with the tolerance lowered to 0.0001 sapply( 1:32, function(i) is.near.equil(wjd(Y=guess(method = "stoich", iguess=i)), tol=0.0001) ) ## using run.wjd(): 20 crystalline amino acids iaa <- info(aminoacids(""), "cr") # starting with one mole of each amino acid w <- run.wjd(iaa, Y=rep(1, 20), imax=20) # the following is TRUE (FALSE if tol is left at default) is.near.equil(w, tol=0.012) # in this assemblage, what are the amino acids # in order of increasing abundance? aminoacids()[order(w\$X)] # because the elements are redistributed among the species, # the total number of moles of species does not remain constant sum(w\$X) # <20 ```