Find the quantities of chemical species, subject to constant elemental bulk composition of the system, that minimize the Gibbs energy of the system.
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)
matrix, chemical formulas of the species (elements on columns)
numeric, the Gibbs energies / RT, at a single temperature (length equal to number of species)
numeric, initial solution, a positive set of values (numbers of moles, length equal to number of species)
numeric, pressure in atmospheres
numeric, number of values of fractional distance change (lambda) tested at each step.
numeric, maximum number of iterations
numeric, Gibbs energy change of system, as a fraction of total system energy in previous step, below which iterations will stop
list, output from
logical, make a plot?
numeric, which elements for which to make plots
numeric, maximum difference in chemical potentials that counts as equilibrium
logical, don't output messages?
numeric, numbers of moles of the elements
character, method used for generating an initial solution
numeric, minimum mole number for 'central' method
numeric, which guess to return
numeric, which combination(s) of variable species to use (NULL for all)
numeric, species indices (rownumbers of
numeric, temperature in °C
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.guess are provided to make access to the thermodynamic database in CHNOSZ easier.
The default values of
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.
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 (
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
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
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
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 °C.
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
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
The second method for
guess ‘stoich’ (and the default for
run.wjd) is to test successive combinations of species whose elemental compositions are linearly independent.
The linearly independent combinations tested are all those from
ic is NULL, or only those identified by
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.
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
exceed.Ttr = TRUE) to generate the values of
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
The function gives an error if none of the guesses in
Y is near equilibrium, within the tolerance set by
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 (
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.
wjd returns a list with the problem definition and results: elements
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). The latter result is provided to assist in checking mass balance (mostly for debugging more than theoretical reasons).
White, W. B., Johnson, S. M. and Dantzig, G. B. (1958) Chemical equilibrium in complex mixtures. J. Chem. Phys. 28, 751–755. http://dx.doi.org/10.1063/1.1744264
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
## 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