data(thermo)## set up a system and calculate
## chemical affinities of formation reactions
basis(c("SiO2","MgO","H2O","O2"),c(-5,-5,0,999))
species(c("quartz","enstatite","forsterite"))
# chemical affinities (A/2.303RT) at 25 deg C and 1 bar
affinity()
# at higher temperature and pressure
affinity(T=500,P=2000)
# ten different temperatures at one pressure
affinity(T=c(500,1000,10),P=2000)
# at 25 temperatures and pressures
affinity(T=c(500,1000,5),P=c(1000,5000,5))
# as a function of logarithm of activity of MgO
affinity(MgO=c(-10,-5,10))
## equilibrium constants of formation reactions
affinity(property="logK")
# Standard molal Gibbs energies of species,
# user units (default: cal/mol)
affinity(property="G.species")
# Standard molal Gibbs energies of reactions
affinity(property="G")
## amino acid synthesis at low and high temperatures
## after Amend and Shock, 1998
# select the basis species and species of interest
# and set their activities, first for the 18 degree C case
basis(c("H2O", "CO2", "NH4+", "H2", "H+", "H2S"),
log10(c(1, 1e-4, 5e-8, 2e-9, 5e-9, 1e-15)))
species(sort(aminoacids("Z")),
log10(c(3.9, 0.7, 1.1, 3.3, 0.5, 3.8, 1.0, 5.8, 1.2, 0.7,
0.8, 1.0, 2.8, 0.5, 0.5, 4.6, 5.8, 0.6, 0.9, 2.8)/1e9))
T <- 18
TK <- convert(T, "K")
# calculate A/2.303RT (dimensionless), convert to G of reaction (cal/mol)
a <- affinity(T=T)
G.18.cal <- convert(unlist(a$values), "G", T=TK)
# covvert to kJ/mol
G.18.kJ <- convert(G.18.cal, "J")/1000
# the 100 degree C case
basis(c("H2O", "CO2", "NH4+", "H2", "H+", "H2S"),
log10(c(1, 2.2e-3, 2.9e-6, 3.4e-4, 1.9e-6, 1.6e-3)))
species(1:20, log10(c(2.8e-9, 5.0e-10, 7.9e-10, 2.4e-9, 3.6e-10,
2.7e-9, 7.2e-10, 4.2e-9, 8.6e-10, 5.0e-10, 5.7e-10, 7.2e-10, 2.0e-9,
3.6e-10,3.6e-10, 3.3e-9, 4.2e-9, 4.3e-10, 6.5e-10, 2.0e-9)))
T <- 100
TK <- convert(T, "K")
a <- affinity(T=T)
G.100.cal <- convert(unlist(a$values), "G", T=TK)
G.100.kJ <- convert(G.100.cal, "J")/1000
# the average oxidation states of carbon
Z.C <- ZC(thermo$obigt$formula[thermo$species$ispecies])
# put everything together a la Table 3 in the paper
print(out <- data.frame(G.18=G.18.kJ, G.100=G.100.kJ, Z.C=Z.C))
# make a plot; set units to get correct label
E.units("J")
plot(out$Z.C, out$G.18, pch=20, xlim=c(-1.1, 1.1), ylim=c(-200, 500),
xlab=axis.label("ZC"), ylab=axis.label("DGr"))
points(out$Z.C, out$G.100, col="red", pch=20)
legend("topleft", pch=c(20, 20), col=c("black", "red"),
legend=describe.property(c("T", "T"), c(18, 100)))
title(main="Amino acid synthesis, after Amend and Shock, 1998")
# 9 amino acids have negative delta Gr under hydrothermal conditions
# (cf. AS98 with 11; we are using more recent thermodynamic data)
stopifnot(sum(out$G.100 < 0)==9)
# reset units and species to run next examples
E.units("cal")
species(delete=TRUE)
## affinities of metabolic reactions
## after Amend and Shock, 2001, Fig. 7
# use aq state for all basis species (including O2)
basis(c("CO2", "H2", "NH3", "O2", "H2S", "H+"), "aq")
# we're going to make H2O
species("H2O")
# a function to create the plots
doplot <- function(T) {
res <- 20
# calculate affinity/2.303RT as a function of loga(H2) and loga(O2)
a <- affinity(H2=c(-10, 0, res), O2=c(-10, 0, res), T=T)
T.K <- convert(T, "K") # temperature in Kelvin
acal <- convert(a$values[[1]], "G", T.K) # affinity (cal/mol)
akJ <- convert(acal, "J")/1000 # affinity (kJ/mol)
# now contour the values
xyvals <- seq(-10, 0, length.out=res)
contour(x=xyvals, y=xyvals, z=t(akJ), levels=seq(-150, -250, -20),
labcex=1, xlab=axis.label("H2"), ylab=axis.label("O2"))
# show the temperature
legend("topleft", bg="white", cex=1,
legend=describe.property("T", T, digits=0, ret.val=TRUE) )
}
# plot layout with space for title at top
layout(matrix(c(1, 1, 2, 3, 4, 5), ncol=2, byrow=TRUE), heights=c(1, 4, 4))
par(mar=c(0, 0, 0, 0))
plot.new()
# we use subcrt() to generate a reaction for titling the plot
rxnexpr <- describe.reaction(subcrt("H2O", 1)$reaction, states="all")
# also in the title is the property with its units
E.units("J")
Gexpr <- axis.label("DGr", prefix="k")[[1]]
text(0.5, 0.6, substitute(paste(G~~"for"~~r), list(G=Gexpr, r=rxnexpr)), cex=2)
text(0.5, 0.2, "after Amend and Shock, 2001 Figure 7", cex=2)
# now make the plots
par(mar=c(3, 3, 0.5, 0.5), cex=1.3, mgp=c(2, 1, 0))
sapply(c(25, 55, 100, 150), doplot)
# affinity() can handle the three dimensions simultaneously
print(affinity(H2=c(-10, 0, 3), O2=c(-10, 0, 3), T=c(25, 150, 4))$values)
# this is so the plots in the next examples show up OK
E.units("cal")
layout(matrix(1))
## calculations along a transect: methanogenesis and biosynthetic
## reactions in hydrothermal systems, after Shock and Canovas, 2010
# this file has their mixing path results for Rainbow hydrothermal field
file <- system.file("extdata/cpetc/SC10_Rainbow.csv", package="CHNOSZ")
rb <- read.csv(file, check.names=FALSE)
# write all synthesis reactions in terms of these basis species
# it's okay not to set the activities of the basis species now
# because they'll be changing along with temperature
basis(c("CO2", "H2", "NH4+", "H2O", "H2S", "H+"))
# now a selection of the species from SC10, with activities equal to 1e-6
species(c("CH4", "formaldehyde", "ethylene", "glycolic acid",
"n-nonanoic acid", "leucine", "aspartic acid", "tryptophan", "deoxyribose",
"adenine", "cytosine"), -6)
# the exception is methane; unlike SC10 we use a constant activity 1e-3
# (accounting for variable activities of the species of interest here
# is possible but would require longer code ....)
species("CH4", -3)
# synchronized change of temperature and five basis activities
a <- affinity(T=rb$T, CO2=rb$CO2, H2=rb$H2, `NH4+`=rb$`NH4+`, H2S=rb$H2S, pH=rb$pH)
# the tricky part: affinity() uses dimensionless values (A/2.303RT)
# but we want to show the values in cal/mol
a$values <- lapply(a$values, function(val) {
-convert(val, "G", T=convert(a$vals[[1]], "K")) })
# if we didn't have balance=1 here the values would be
# divided by the number of moles of CO2 in the reactions ...
diagram(a, balance=1, ylim=c(-100000, 100000), ylab=axis.label("A"),
col=topo.colors(4), lwd=2)
# add a zero-affinity line and a title
abline(h=0, lty=2, lwd=2)
title(main="Affinities of organic synthesis, after Shock and Canovas, 2010")Run the code above in your browser using DataLab