wjd: Gibbs Energy Minimization by Steepest Descent

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

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=c("central", "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

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)

numeric, initial solution, a positive set of values (numbers of moles, length equal to number of species)

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

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?

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

numeric, which combination(s) of variable species to use (NULL for all)

ispecies

numeric, species indices (rownumbers of thermo$obigt)

numeric, temperature in $^{\circ}$C

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). The latter result is provided to assist in checking mass balance (mostly for debugging more than theoretical reasons).

encoding

UTF-8

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 $^{\circ}$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 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 in 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.

References

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

Examples

Run this code

opar <- par(no.readonly=TRUE)
## 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
par(thermo$opar)

## 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

## run.wjd with proteins: cell periphery of yeast
# get the proteins in the requested location
y <- yeastgfp("cell.periphery")
# get the amino acid compositions of the proteins
aa <- more.aa(y$protein, "Sce")
# don't use those with NA abundance or sequence
ina <- is.na(y$abundance) | is.na(aa$chains)
aa <- aa[!ina, ]
# let's try normalizing the proteins to single residues
# columns 6:25 are the actual amino acid counts
aa.625 <- aa[, 6:25]
aa[, 6:25] <- aa.625 / rowSums(aa.625)
# add proteins to thermo$protein
add.protein(aa)
# add proteins to thermo$obigt
iobigt <- info(paste(aa$protein, aa$organism, sep="_"))
# use equal initial abundances, with total equal to yeastGFP abundances
Y <- rep(mean(y$abundance[!ina]), length(y$abundance[!ina]))
# run the Gibbs energy minimization
w <- run.wjd(iobigt, Y=Y, imax=100)
# make a log-log plot
plot(log10(y$abundance[!ina]), log10(w$X), xlim=c(1.5, 5), ylim=c(1.5, 5))
# get the element potentials (tolerating "close enough" to equilibrium)
emu <- equil.potentials(w, tol=1e7)
# then the logarithms of activities of the basis species
basis("CHNOS")
bl <- basis.logact(emu)
# make a title and legend
title(main="calculated vs observed abundances: yeast cell periphery")
basis(names(bl), bl)
legend("topleft", describe.basis(digits=2))

Run the code above in your browser using DataLab