Equilibrium Chemical Activity Diagrams

Share:

Description

Plot equilibrium chemical activity (1-D speciation) or equal-activity (2-D predominance) diagrams as a function of chemical activities of basis speecies, temperature and/or pressure.

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
  diagram(eout, what = "loga.equil", alpha = FALSE, normalize = FALSE,
    as.residue = FALSE, balance=NULL, groups=as.list(1:length(eout$values)),
    xrange=NULL, mar=NULL, yline=par("mgp")[1]+0.7, side=1:4,
    ylog=TRUE, xlim=NULL, ylim=NULL, xlab=NULL, ylab=NULL,
    cex=par("cex"), cex.names=1, cex.axis=par("cex"),
    lty=NULL, lwd=par("lwd"), dotted=0,
    bg=par("bg"), col=par("col"), col.names=par("col"), fill=NULL,
    names=NULL, main=NULL, legend.x="topright", add=FALSE, plot.it=TRUE,
    tplot=TRUE, ...)
  strip(affinity, ispecies = NULL, col = NULL, ns = NULL, 
    xticks = NULL, ymin = -0.2, xpad = 1, cex.names = 0.7)
  find.tp(x)

Arguments

eout

list, object returned by equilibrate

or affinity

what

character, what property to calculate and plot

alpha

logical, for speciation diagrams, plot degree of formation instead of activities?

normalize

logical, divide chemical affinities by balance coefficients (rescale to whole formulas)?

as.residue

logical, divide chemical affinities by balance coefficients (no rescaling)?

balance

character, balancing constraint; see balance

groups

list of numeric, groups of species to consider as a single effective species

xrange

numeric, range of x-values between which predominance field boundaries are plotted

mar

numeric, margins of plot frame

yline

numeric, margin line on which to plot the y-axis name

side

numeric, which sides of plot to draw axes

xlim

numeric, limits of x-axis

ylim

numeric, limits of y-axis

xlab

character, label to use for x-axis

ylab

character, label to use for y-axis

ylog

logical, use a logarithmic y-axis (on 1D degree diagrams)?

cex

numeric, character expansion (scaling relative to current)

cex.names

numeric, character expansion factor to be used for names of species on plots

cex.axis

numeric, character expansion factor for names of axes

lty

numeric, line types to be used in plots

lwd

numeric, line width

dotted

numeric, how often to skip plotting points on predominance field boundaries (to gain the effect of dotted or dashed boundary lines)

bg

character, background color for legend

col

character, color of activity lines (1D diagram) or predominance field boundaries (2D diagram), or colors of bars in a strip diagram (strip)

col.names

character, colors for labels of species

fill

character, colors used to fill predominance fields

names

character, names of species for activity lines or predominance fields

main

character, a main title for the plot; NULL means to plot no title

legend.x

character, description of legend placement passed to legend

add

logical, add to current plot?

plot.it

logical, make a plot?

tplot

logical, set up plot with thermo.plot.new?

affinity

list, object returned by affinity

ispecies

numeric, which species to consider (default of NULL is to consider all species)

ns

numeric, numbers of species, used to make inset plots for strip diagrams

xticks

numeric, location of supplemental tick marks on x-axis

ymin

numeric, lower limit of y-axis

xpad

numeric, amount to extend x-axis on each side

x

matrix, value of the predominant list element from diagram

...

additional arguments passed to plot or barplot

Details

diagram takes as its primary input the results from equilibrate and displays diagrams representing the equilibrium chemical activities of the species. 0-D diagrams, at a single point, are shown as barplots. 1-D diagrams, for a single variable on the x-axis, are plotted as lines. 2-D diagrams, for two variables, are plotted as predominance fields. The allowed variables are any that affinity accepts: temperature, pressure, or the chemical activities of the basis species

If alpha is TRUE, the fractional degrees of formation (ratios of activities to total activity) are plotted. This setting is useful for visualizing the molalities (activities for ideal interactions) of monomer groups in a system of biopolymers (e.g. proteins). If groups is supplied, the activities of the species identified in each numeric vector of this list are summed together and subsequently treated as a single species; the names of the new species are taken from the names of this list.

For 1-D diagrams, the default setting for the y-axis is a logarithmic scale (unless alpha is TRUE) with limits corresponding to the range of logarithms of activities (or 0,1 if alpha is TRUE); these actions can be overridden by ylog and ylim. A legend is placed at the location identified by legend.x, or omitted if legend.x is FALSE. If legend.x is NA, instead of a legend, the lines are labeled with the names of the species near the maximum value. The line type and line width can be controlled with lty and lwd, respectively.

On 2-D diagrams, the fields represent the species with the highest equilibrium activity. fill determines the color of the predominance fields, col that of the boundary lines. By default, heat.colors are used to fill the predominance fields in diagrams on the screen plot device. The style of the boundary lines on 2-D diagrams can be altered by supplying one or more non-zero integers in dotted, which indicates the fraction of line segments to omit; a value of 1 or NULL for dotted has the effect of not drawing the boundary lines.

normalize and as.residue apply only to the 2-D diagrams, and only when eout is the output from affinity. With normalize, the activity boundaries are calculated as between the residues of the species (the species divided by the balance coefficients), then the activities are rescaled to the whole species formulas. With as.residue, the activity boundaries are calculated as between the residues of the species, and no rescaling is performed.

For all diagrams, the names of the species and their colors in col.names can be changed, as can cex, cex.names, and cex.axis to adjust the overall character expansion factors (see par) and those of the species names and axis labels. The x- and y-axis labels are automatically generated unless they are supplied in xlab and ylab. A new plot is started unless add is TRUE. If plot.it is FALSE, no plot will be generated but all the intermediate computations will be performed and the results returned.

diagram also accepts the output from affinity, for which three actions are possible: 1) plot a property of a reaction, such as the equilibrium constant (logK) (0-D or 1-D); 2) plot the equilibrium activity of a selected basis species in all of the formation reactions (0-D, 1-D or 2-D); 3) plot predominance fields, based on the relative magnitudes of the affinities of the formation reactions (2-D only).

Some of the arguments have different effects when the output from affinity is being used instead of the equilibrium activities from equilibrate. If what is missing, option (1) is selected if the number of dimensions is 0 or 1, and option (3) is selected if the number of dimensions is 2. The latter is referred to as the maximum affinity method. In cases where it applies (see Warning), the maximum affinity method is much faster than an equilibrium calculation. balance is the option, sent to balance, that determines the balancing coefficients (this argument has no effect on the calculations of equilibrium activities.) If what is the name of a basis species, it refers to option (2) above. A contour plot is made in the case of 2-D diagrams of the equilibrium activity of a basis species (see demos("CO2Ac"), and only the first species of interest is used in the calculation; a warning is produced if there is more than one.

A different incarnation of 1-D speciation diagrams is provided by strip. This function generates any number of strip diagrams in a single plot. The diagrams are made up of colors bars whose heights represent the relative abundances of species; the color bars are arranged in order of abundance and the total height of the stack of colors bars is constant. If ispecies is a list, the number of strip diagrams is equal to the number of elements of the list, and the elements of this list are numeric vectors that identify the species to consider for each diagram. The strips are labeled with the names of ispecies. If col is NULL, the colors of the bars are generated using rainbow. Supplemental ticks can be added to the x-axis at the locations specified in xtick; they are larger than the standard ticks and have colors corresponding to those of the color bars. ymin can be decreased in order to add more space at the bottom of the plot, and xpad can be changed in order to increase or decrease the size of the x-axis relative to the width of the strips. An inset dot-and-line plot is created below each strip if ns is given. This argument has the same format as ispecies, and can be used e.g. to display the relative numbers of species for comparison with the stability calculations.

find.tp finds the locations in a matrix of integers that are surrounded by the greatest number of different values. The function counts the unique values in a 3x3 grid around each point and returns a matrix of indices (similar to which(..., arr.ind = TRUE)) for the maximum count (ties result in more than one pair of indices). It can be used with the output from diagram for calculations in 2 dimensions to approximately locate the triple points on the diagram.

Value

For speciation diagrams, an invisible list of the chemical activities of the species, or their degrees of formation (if alpha is TRUE), at each point. For predominance diagrams, an invisible list with elements species, the dataframe describing the species, out, which species predominates at each grid point, and A, a list of the calculated values of the chemical affinity (per balanced quantity) (log10 dimensionless) at each point.

Warning

The maximum affinity method is premised on the calculation of the affinities of formation reactions at equal activities of the species of interest. Then, the species with the highest affinity of formation, after normalization by the balancing coefficients, corresponds to the predominant species in an equilibrium calculation. The examples below “work” because they are relatively simple - the balancing coefficients are unity or all the same value (aqueous aluminum example), or the species are solids with unit activities (the mineral examples). The examples shown for proteins elsewhere also take the balancing coefficients to unity, after normalizing by protein length. However, if aqueous species are present with different balancing coefficients, the maximum affinity method is not dependable, as shown in the TCA metabolite example below.

References

Aksu, S. and Doyle, F. M. (2001) Electrochemistry of copper in aqueous glycine solutions. J. Electrochem. Soc. 148, B51–B57. http://dx.doi.org/10.1149/1.1344532

Helgeson, H. C. (1970) A chemical and thermodynamic model of ore deposition in hydrothermal systems. Mineral. Soc. Amer. Spec. Pap. 3, 155–186. http://www.worldcat.org/oclc/583263

Helgeson, H. C., Delany, J. M., Nesbitt, H. W. and Bird, D. K. (1978) Summary and critique of the thermodynamic properties of rock-forming minerals. Am. J. Sci. 278-A, 1–229. http://www.worldcat.org/oclc/13594862

LaRowe, D. E. and Helgeson, H. C. (2007) Quantifying the energetics of metabolic reactions in diverse biogeochemical systems: electron flow and ATP synthesis. Geobiology 5, 153–168. http://dx.doi.org/10.1111/j.1472-4669.2007.00099.x

Majzlan, J., Navrotsky, A., McClesky, R. B. and Alpers, C. N. (2006) Thermodynamic properties and crystal structure refinement of ferricopiapite, coquimbite, rhomboclase, and Fe2(SO4)3(H2O)5. Eur. J. Mineral. 18, 175–186. http://dx.doi.org/10.1127/0935-1221/2006/0018-0175

Seewald, J. S. (1997) Mineral redox buffers and the stability of organic compounds under hydrothermal conditions. Mat. Res. Soc. Symp. Proc. 432, 317–331. http://dx.doi.org/10.1557/PROC-432-317

Seewald, J. S. (2001) Aqueous geochemistry of low molecular weight hydrocarbons at elevated temperatures and pressures: Constraints from mineral buffered laboratory experiments Geochim. Cosmochim. Acta 65, 1641–1664. http://dx.doi.org/10.1016/S0016-7037(01)00544-0

Tagirov, B. and Schott, J. (2001) Aluminum speciation in crustal fluids revisited. Geochim. Cosmochim. Acta 65, 3965–3992. http://dx.doi.org/10.1016/S0016-7037(01)00705-0

See Also

Other examples are present in the help for protein and buffer, and even more can be found in demos. See the vignette hotspring.Rnw for an example of the strip charts.

Examples

  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
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
### 1-D diagrams: logarithms of activities

## Aqueous sulfur species (after Seewald, 1997 and 2001)
basis("CHNOS+")
basis("pH", 5)
species(c("H2S", "S2-2", "S3-2", "S2O3-2", "S2O4-2", 
  "S3O6-2", "S5O6-2", "S2O6-2", "HSO3-", "SO2", "HSO4-"))
a <- affinity(O2=c(-50, -15), T=325, P=350)
e <- equilibrate(a, loga.balance=-2)
diagram(e, ylim=c(-30, 0), legend.x="topleft")
title(main=paste("Aqueous sulfur speciation, 325 degC, 350 bar\n",
  "After Seewald, 1997"))
# try it with and without the loga.balance argument (total activity of 
# the balanced quantity, in this case H2S aka sulfur) 

## Degrees of formation of ionized forms of glycine
## After Fig. 1 of Aksu and Doyle, 2001 
basis("CHNOS+")
species(ispecies <- info(c("glycinium", "glycine", "glycinate")))
a <- affinity(pH=c(0, 14))
e <- equilibrate(a)
diagram(e, alpha=TRUE, lwd=1)
title(main=paste("Degrees of formation of aqueous glycine species\n",
  "after Aksu and Doyle, 2001"))

## Degrees of formation of ATP species as a function of 
## temperature, after LaRowe and Helgeson, 2007, Fig. 10b
# to make a similar diagram, activity of Mg+2 here is set to
# 10^-4, which is different from LH07, who used 10^-3 total molality
basis(c("CO2", "NH3", "H2O", "H3PO4", "O2", "H+", "Mg+2"),
  c(999, 999, 999, 999, 999, -5, -4))
species(c("HATP-3", "H2ATP-2", "MgATP-2", "MgHATP-"))
a <- affinity(T=c(0, 120, 25))
e <- equilibrate(a)
diagram(e, alpha=TRUE)  
title(main=paste("Degrees of formation of ATP species,\n",
  "pH=5, log(aMg+2)=-3. After LaRowe and Helgeson, 2007"),
  cex.main=0.9)

### 2-D diagrams: predominance diagrams
### these use the maximum affinity method

## Fe-S-O at 200 deg C, after Helgeson, 1970
basis(c("Fe", "O2", "S2"))
species(c("iron", "ferrous-oxide", "magnetite",
  "hematite", "pyrite", "pyrrhotite"))
# include the high-temperature phase of pyrrhotite
species("pyrrhotite", "cr2")
a <- affinity(S2=c(-50, 0), O2=c(-90, -10), T=200)
diagram(a, fill="heat")
title(main=paste("Fe-S-O, 200 degrees C, 1 bar",
  "After Helgeson, 1970", sep="\n"))

## pe-pH diagram for hydrated iron sulfides,
## goethite and pyrite, after Majzlan et al., 2006
# add some of these species to the database
add.obigt()
basis(c("Fe+2", "SO4-2", "H2O", "H+", "e-"), 
  c(0, log10(3), log10(0.75), 999, 999))
species(c("rhomboclase", "ferricopiapite", "hydronium jarosite",
  "goethite", "melanterite", "pyrite"))
a <- affinity(pH=c(-1, 4, 256), pe=c(-5, 23, 256))
d <- diagram(a, main="Fe-S-O-H, after Majzlan et al., 2006")
# the first four species show up along the top of the diagram
stopifnot(all.equal(unique(t(d$predominant)[256,]), 1:4))
water.lines(yaxis="pe")
text(3, 22, describe.basis(thermo$basis[2:3,], digits=2, oneline=TRUE))
text(3, 21, describe.property(c("T", "P"), c(25, 1), oneline=TRUE))
# reset the database
data(thermo)

## Aqueous Aluminum Species F-/OH-, after Tagirov and Schott, 2001
# some of the species are not in the default databse
add.obigt()
# the 999s have no effect on the diagram:
# pH and log_a(F-) are plotting variables
# O2 is not in the formation reactions
# Al+3 is the balanced quantity
basis(c("Al+3", "F-", "H+", "O2", "H2O"), c(rep(999, 4), 0))
species(c("Al+3", "Al(OH)4-", "AlOH+2", "Al(OH)2+", "Al(OH)3",
  "AlF+2", "AlF2+", "AlF3", "AlF4-", "Al(OH)2F2-", "Al(OH)2F",
  "AlOHF2"), "aq")
a <- affinity(pH=c(0, 10), "F-"=c(-1, -9), T=200)
diagram(a, fill="heat")
title(main=paste("Aqueous aluminium species, T=200 C, P=Psat\n",
  "after Tagirov and Schott, 2001"))
# restore thermodynamic database to default
data(thermo)

## Temperature-Pressure: kayanite-sillimanite-andalusite
# cf. Fig. 49 of Helgeson et al., 1978
# this is a system of one component (Al2SiO5), but we need the same 
# number of basis species as elements; and avoid using H2O or other 
# aqueous species because of T/P limits of the water() calculations;
basis(c("corundum", "quartz", "oxygen"))
species(c("kyanite", "sillimanite", "andalusite"))
# database has transition temperatures of kyanite and andalusite
# at 1 bar only, so we permit calculation at higher temperatures
a <- affinity(T=c(200, 900, 99), P=c(0, 9000, 101), exceed.Ttr=TRUE)
d <- diagram(a, fill=NULL)
bexpr <- sapply(c("Al2O3", "SiO2", "H2O"), expr.species, simplify=FALSE)
btext <- substitute(Al2O3 - SiO2 - H2O, unlist(bexpr))
mtitle(c(as.expression(btext), "after Helgeson et al., 1978"))
# find the approximate position of the triple point
tp <- find.tp(d$predominant)
Ttp <- a$vals[[1]][tp[1, 2]]
Ptp <- rev(a$vals[[2]])[tp[1, 1]]
points(Ttp, Ptp, pch=10, cex=5)
# some testing of the overall geometry
stopifnot(species()$name[d$predominant[1, 1]]=="andalusite")
stopifnot(species()$name[d$predominant[1, 101]]=="kyanite")
stopifnot(species()$name[d$predominant[99, 101]]=="sillimanite")

## calculate the equilibrium logarithm of activity of a 
## basis species in different reactions
basis("CHNOS")
species(c("ethanol", "lactic acid", "deoxyribose", "ribose"))
a <- affinity(T=c(0, 150))
diagram(a, what="O2", legend.x="topleft", col=rev(rainbow(4)), lwd=2)
title(main="equilibrium logfO2 for 1e-3 mol/kg of CO2 and ... ")

Want to suggest features or report bugs for rdrr.io? Use the GitHub issue tracker.