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")
## Activity of glycine as a function of those of
## HCN and formaldehyde (200 C, 300 bar)
## After Schulte and Shock, 1995, Fig. 5
# we can define the basis as this:
basis(c("formaldehyde","H2O","HCN","O2"))
species("glycine")
a <- affinity(HCHO=c(-10,-2,9),HCN=c(-18,-2,9),T=200,P=300)
# that gave us *affinities* (dimensionless) for logact(glycine)=-3
# (the default). we can now find the *activities* that
# correspond to affinity=0
logact.glycine <- species()$logact + a$values[[1]]
# note transposition of the z-value matrix in the following command
contour(x=-10:-2,y=seq(-18,-2,by=2),z=t(logact.glycine),
xlab=axis.label("HCHO"),ylab=axis.label("HCN"),
labcex=1,xaxs="i",yaxs="i")
title(main=paste("log activity glycine, after Schulte and Shock, 1995",
"200 deg C, 300 bar, logaH2O = 1",sep="\n"))
## 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(c("alanine","argininate","asparagine","aspartate","cysteine",
"glutamate","glutamine","glycine","histidine","isoleucine",
"leucine","lysinium","methionine","phenylalanine","proline",
"serine","threonine","tryptophan","tyrosine","valine"),
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))
Tc <- 18
T <- convert(Tc,"K")
# converting A (dimensionless) to G of reaction (cal/mol)
# is like converting log K to standard G of reaction
AS98.18 <-
convert(convert(as.numeric(affinity(T=Tc)$values),"G",T=T),"J")/1000
# the 100 degree C case
Tc <- 100
T <- convert(Tc,"K")
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)))
AS98.100 <-
convert(convert(as.numeric(affinity(T=Tc)$values),"G",T=T),"J")/1000
# the nominal carbon oxidation state
Z.C <- ZC(as.character(thermo$obigt$formula[thermo$species$ispecies]))
# put them together
print(data.frame(T100=AS98.100,T18=AS98.18,Z.C=Z.C))
# values not exactly reproducing AS98 - different amino acid parameters
# forget species to run next example
species(delete=TRUE)
## affinities of metabolic reactions
## after Amend and Shock, 2001, Fig. 7
basis(c("CO2","H2","NH3","O2","H2S","H+"))
basis(c("O2","H2"),"aq") # O2 and H2 were gas
species("H2O")
doplot <- function(T) {
res <- 20
a <- affinity(H2=c(-10,0,res),O2=c(-10,0,res),T=T)
T.K <- convert(T,"K") # temperature in Kelvin
a <- convert(a$values[[1]],"G",T.K) # affinities (cal/mol)
a <- convert(a,"J") # affinities (Joule)
contour(x=seq(-10,0,length.out=res),
y=seq(-10,0,length.out=res),z=t(a/1000),
labcex=1,xlab=axis.label("H2"),ylab=axis.label("O2"))
}
layout(matrix(c(1,1,2,3,4,5),ncol=2,byrow=TRUE),heights=c(1,4,4))
T <- c(25,55,100,150)
par(mar=c(0,0,0,0))
plot.new()
text(0.5,0.1,paste(c("H2(aq) + 0.5O2(aq) = H2O(liq)\n\n",
"after Amend and Shock, 2001")),cex=2)
par(mar=c(3,3,0.5,0.5),cex=1.3,mgp=c(2,1,0))
for(i in 1:length(T)) doplot(T[i])
# this is so the plots in the next examples show up OK
layout(matrix(1))
## continuation of last example, now in three dimensions
print(affinity(H2=c(-10,0,3),O2=c(-10,0,3),T=c(25,150,4))$values)
## calculations on a transect
# suppose that temperature and oxygen fugacity both change in space
# (say from 1 to 6 meters), and we have six values for each but want to
# interpolate them to make a plot with smooth curves
T <- splinefun(1:6,c(0,25,30,40,55,75))(seq(1,5,length.out=100))
O2 <- splinefun(1:6,c(-90,-83,-78,-73,-68,-63))(seq(1,5,length.out=100))
# what system could this be?
basis("CHNOS+")
species(paste("CSG",c("METBU","METVO","METTL","METJA"),sep="_"))
# now pass T and O2 to affinity, which because their lengths
# are greater than three, treats them like coordinates for a
# transect through chemical potential space rather than
# the definition of a 2-dimensional grid
a <- affinity(T=T,O2=O2)
diagram(a,ylim=c(-4,-2))
title(main=paste("Computed abundances of surface-layer proteins",
"as a function of T and logfO2",sep="\n"))Run the code above in your browser using DataLab