data(thermo)
### 0-D property and logarithm of activity plots
## properties of amino acids
basis("CHNOS")
basis("H2S",-25)
species(c("glutamic acid","methionine","isoleucine",
"leucine","tyrosine"))
# Gibbs energy of formation reactions per CO2
a <- affinity(property="G")
diagram(a,balance=1)
# affinity of reactions per CO2
a <- affinity()
diagram(a,property="A")
## logarithms of activities of amino acids
# balanced on CO2
diagram(a)
# balanced on amino acid backbone, using abundance.old
d1 <- diagram(a,balance=1)
# balanced on amino acid backbone, using abundance
d2 <- diagram(a,balance=1,residue=TRUE)
stopifnot(all(as.numeric(d2$logact) -
as.numeric(d1$logact) < thermo$opt$cutoff))
## 1-D plots
a <- affinity(O2=c(-80,-70))
# affinities of formation reactions
d <- diagram(a,property="A",balance=1,ylim=c(-80,30))
# logarithm of activity
d <- diagram(a,balance=1,ylim=c(-9,0))
title(main=paste("Relative abundances of amino acids"))
### 1-D property plot
## for hydrous cordierite = cordierite + H2O
## after Helgeson et al., 1978
basis(c("cordierite,hydrous","Mg+2","SiO2","H2O","O2","H+"))
species("cordierite")
# water.SUPCRT92 can only get us up to 5000 bar
# (7000 and 10000 bar shown in the original)
P <- c(1,2,3,5)*1000
col <- rainbow(length(P))
for(i in 1:length(P)) {
a <- affinity(property="logK",T=c(20,800),P=P[i])
diagram(a,add=(i!=1),ylim=c(-4,2),legend.x=NULL,
color=col[i],title="")
}
legend("topright",lty=1,col=col,legend=paste(P,"bar"))
title(main=paste("hydrous cordierite = cordierite + H2O",
"After Helgeson et al., 1978",sep=""),cex.main=0.9)
# this calculation could also be done through the
# subcrt function, without having to set up basis species
### 1-D logarithm of activity plots
### (speciation diagrams)
## Aqueous sulfur species (after Seewald, 1997)
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-'))
diagram(affinity(O2=c(-50,-15,50),T=325,P=350),
ylim=c(-30,0),logact=-2,legend.x="topleft")
title(main=paste('Aqueous sulfur speciation, 325 degC, 350 bar<n>',
'After Seewald, 1997'))
# try it with and without the logact argument (total activity of
# the balanced quantity, in this case H2S aka sulfur)
## Degrees of formation of ionized forms of glycine
## as a function of pH
## After Fig. 1 of Aksu and Doyle, 2001
basis('CHNOS+')
species(ispecies <- info(c('glycinium','glycine','glycinate')))
t <- affinity(pH=c(0,14))
diagram(t,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
# Mg+2 here is a perfectly mobile component
# cf. LH07, where bulk composition of Mg+2 is specified
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-'))
diagram(affinity(T=c(0,120,25)),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)
## speciation of phosphate as a function of ionic strength
basis('CHNOPS+')
T <- c(25,100)
species(c('HPO4-2','H2PO4-'))
diagram(affinity(IS=c(0,0.25),T=T[1]),ylim=c(-3.2,-2.8))
title(main=paste('Aqueous phosphate speciation, pH=7',
'25 and 100 deg C - black and red lines',sep='<n>'))
diagram(affinity(IS=c(0,0.25),T=T[2]),ylim=c(-3.2,-2.8),add=TRUE,color='red')
## phosphate predominance f(IS,pH)
diagram(affinity(pH=c(6.8,7.2),IS=c(0,0.25),T=T[1]))
title(main=paste('Aqueous phosphate predominance, 25 deg C',
'and 100 deg C (dotted overlay)',sep='<n>'))
diagram(affinity(pH=c(6.8,7.2),IS=c(0,0.25),T=T[2]),add=TRUE,dotted=2,names=NULL)
### predominance diagrams (equal activities of
### species as a function of two variables)
## T-P phase diagram for a unary system (SiO2)
## after Ernst, 1976, Fig. 4.4
# the system is SiO2 (one component) but
# the basis species require all the elements:
# basis(c('SiO2','O2'))
# that would identify SiO2(aq) which is okay
# but brings in calculations of properties of
# water which are relatively slow.
basis(c('quartz','O2')) # cr, gas
# browse the species of interest
info('SiO2 ')
# we'll take the crystalline ones
t <- info('SiO2','cr')
species(t)
# calculate chemical affinities
# the do.phases argument is necessary for
# dealing with the phase transitions of minerals
t <- affinity(T=c(0,2000),P=c(0,30000),do.phases=TRUE)
# a predominance diagram coincides with a
# phase diagram for a *unary* system
diagram(t)
title(main='Phases of SiO2<nafter>Ernst, 1976')
# NOTE: phase diagrams for systems of more
# than one component are not implemented
## Distribution of copper on Eh-pH diagram
## after Fig. 15 of Pourbaix, 1949
basis(c('Cu+2','Cl-','H2O','H+','e-'))
basis('Cl-',-2)
# two aqueous species, three solid ones
# (our CuCl is aq but it was cr in P's Fig. 15)
species(c('CuCl','Cu+2','copper','cuprite','tenorite'))
# (which is equivalent to ...)
# species(c('CuCl','Cu+2','Cu','Cu2O','CuO'),c('','','','','','cr'))
t <- affinity(pH=c(0,7),Eh=c(-0.1,1))
# this one corresponds to activity contours of
# aqueous species at 10^-3 (the default aq activity in CHNOSZ)
diagram(t,color=NULL)
# here we set activities to unity; aq-cr boundaries change
species(c('CuCl','Cu+2'),c(0,0))
t <- affinity(pH=c(0,7),Eh=c(-0.1,1))
diagram(t,add=TRUE,names=NULL,col="blue",color=NULL)
water.lines()
title(main=paste('H2O-Cl-(Cu); activities of 10^-3 (black)<n>',
'or 0 (blue); after Pourbaix, 1949'))
# we could do it at higher temperatures too
# t <- affinity(pH=c(0,7),Eh=c(-0.1,1),T=100)
# diagram(t,add=TRUE,dotted=2,col='red')
## Eh-pH diagrams for copper-water-glycine
## After Fig. 2 of Aksu and Doyle, 2001
## We need to add some species and change some Gibbs energies.
# copy two rows of the data object (actually one row twice)
# and change them as needed
i <- info(c('Cu(Gly)+','glycinate','e-','H+'))
n <- nrow(thermo$obigt <- rbind(thermo$obigt,thermo$obigt[rep(i[1],2),]))
thermo$obigt$name[n-1] <- 'Cu(Gly)2-'
thermo$obigt$name[n] <- 'HCu(Gly)+2'
thermo$obigt$formula[n-1] <- makeup(makeup(c(i[1],i[2],i[3]),''),'')
thermo$obigt$formula[n] <- makeup(makeup(c(i[1],i[4]),''),'')
# In Fig 2b, total log activities of Cu (Cu_T)
# and glycine (L_T) are -4 and -1
basis(c('Cu+2','H2O','H+','e-','glycine','CO2'),c(999,0,999,999,-1,999))
# solid species
species(c('copper','cuprite','tenorite'))
# aqueous species
species(c('glycinium','glycine','glycinate','Cu+','Cu+2','CuO2-2','HCuO2-',
'Cu(Gly)+','Cu(Gly)2','Cu(Gly)2-','HCu(Gly)+2'),-4)
ispecies <- species()$ispecies
# update the Gibbs energies using A&D's Table 1 and Table II
logK <- c(convert(convert(c(0,-146,-129.7,-384.061,-370.647,-314.833,
49.98,65.49,-183.6,-258.5,-298.2)*1000,'cal'),'logK'),15.64,10.1,2.92)
# do it in order so later species take account of prev. species' values
for(i in 1:length(logK)) {
G <- convert(logK[i],'G')
if(i==12) G <- G + thermo$obigt$G[ispecies[8]] +
2*thermo$obigt$G[ispecies[6]]
if(i==13) G <- G + thermo$obigt$G[ispecies[7]] +
2*thermo$obigt$G[ispecies[6]]
if(i==14) G <- G + thermo$obigt$G[ispecies[11]]
thermo$obigt$G[ispecies[i]] <- G
} # done with changing Gibbs free energies!
# we have to get some leftovers out of there or diagram() gets confused
species(c('glycinium','glycine','glycinate'),delete=TRUE)
# make a plot to see if it's working
ispecies <- ispecies[-(1:6)]
afun <- function(cu,gly) {
# from above: our fifth basis species is glycine(-ate,-ium)
basis(rownames(basis())[5],gly)
t <- match(ispecies,species()$ispecies)
species(t,cu)
affinity(pH=c(0,16),Eh=c(-0.6,1.0))
}
diagram(afun(-4,-1))
title(main=paste('Aqueous Copper + Glycine, 25 deg C, 1 bar',
'After Aksu and Doyle, 2001 Fig. 2b',sep='<n>'))
# What's missing? Try glycinate not glycine in reactions at ph > ~9.8
basis(c('Cu+2','H2O','H+','e-','glycinate','CO2'),
c(999,0,999,999,-2,999))
species(c('copper','cuprite','tenorite','Cu+','Cu+2','CuO2-2','HCuO2-',
'Cu(Gly)+','Cu(Gly)2','Cu(Gly)2-','HCu(Gly)+2'))
loga_Cu <- -4
loga_Gly <- -1
diagram(afun(loga_Cu,loga_Gly),color=NULL,col='blue',
names=species()$name,col.names='blue',add=TRUE)
water.lines()
# the glycine ionization constants could be calculated using
# subcrt, here they are taken from A&D Table II
abline(v=c(2.35,9.778),lty=3)
# now put glycinium (low pH) in the basis
basis(c('Cu+2','H2O','H+','e-','glycinium','CO2'),c(999,0,999,999,-2,999))
species(c('copper','cuprite','tenorite','Cu+','Cu+2','CuO2-2','HCuO2-',
'Cu(Gly)+','Cu(Gly)2','Cu(Gly)2-','HCu(Gly)+2'))
diagram(afun(loga_Cu,loga_Gly),color=NULL,col='green',
names=NULL,col.names='green',add=TRUE)
# let's forget our changes to 'thermo' so the example
# below that uses glycine will work as expected
data(thermo)
## a pe-pH diagram for hydrated iron sulfides,
## goethite and pyrite, After Majzlan et al., 2006
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),pe=c(-5,23))
diagram(a)
water.lines(yaxis="pe")
title(main=paste("Fe-S-O-H After Majzlan et al., 2006",
describe(thermo$basis[2:3,],digits=2),sep="<n>"))
## cysteine-cysteinate-cystine Eh-pH
## at 25 and 100 deg C
basis("CHNOSe")
species(c("cysteine","cysteinate","cystine"))
a <- affinity(pH=c(5,10),Eh=c(-0.5,0))
diagram(a,color=NULL)
water.lines()
a <- affinity(pH=c(5,10),Eh=c(-0.5,0),T=100)
diagram(a,col="red",add=TRUE,names=NULL)
water.lines(T=convert(100,"K"),col="red")
title(main=paste("Cysteine Cysteinate Cystine",
"25 and 100 deg C",sep="<n>"))
## Soil Organic Matter CO2-H2O, O2-H2O (after Tardy et al., 1997)
# NH3 is conserved, and H2O is on an axis of the
# diagrams, so their activities are nonsensical here.
# formulas for aqueous species, names for phases ...
basis(c('NH3','water','CO2','O2'),c(999,999,-2.5,-28))
# switch to gaseous CO2 (aq is the default)
basis('CO2','gas')
# load the species of interest
species(c('microflore','plantes','acide fulvique',
'acide humique','humine'))
# proceed with the diagrams
diagram(affinity(H2O=c(-0.6,0.1),CO2=c(-3,-1)))
title(main=paste('Relative stabilities of soil organic matter<n>',
'after Tardy et al., 1997'))
# this would be the O2-H2O diagram
# diagram(affinity(H2O=c(-1,0.5),O2=c(-28.5,-27.5)))
## Aqueous Aluminum Species F-/OH- (afer Tagirov and Schott, 2001)
# The coefficients on H+ and O2 in all the formation reactions
# are zero, so the number of basis species here is three. Al+3
# becomes the conservant, and F- and OH- are being plotted ...
# so their initial activities don't have to be set.
basis(c('Al+3','F-','OH-','H+','O2'),rep(999,5))
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'))
# Increase the x- and y- resolution from the default and calculate
# and plot predominance limits. Names of charged basis species,
# such as "H+", "e-" and the ones shown here, should be quoted
# when given as arguments to affinity(). The OH- values shown here
# correspond to pH=c(0,14) (at unit activity of water).
a <- affinity("OH-"=c(-14,0),"F-"=c(-1,-8),T=200)
diagram(a)
title(main=paste('Aqueous aluminium species, T=200 C, P=Psat<n>',
'after Tagirov and Schott, 2001'))
# We could do this to overlay boundaries for a different pressure
# a.P <- affinity("OH-"=c(-14,0),"F-"=c(-1,-8),T=200,P=5000)
# diagram(a.P,names=NULL,color=NULL,col='blue',add=TRUE)
## Mineral equilibrium activity diagram
## (After Bowers et al., 1984, p. 246)
# Consider the system denoted by BJH84 as
# HCl-H2O-CaO-CO2-MgO-(SiO2) at 300 deg C and 1000 bar
# For our purposes (to make diagrams as a function of the
# activities of Ca+2 and Mg+2), we use this basis.
# notes regarding 999's:
# HCl, CO2, O2 have zero stoichiometric coeffs in the species
# Ca+2, Mg+2 are going to be on the diagram
# SiO2 will be conserved
# (try changing any of the 999's, the diagram will be unaffected)
basis(c('HCl','H2O','Ca+2','CO2','Mg+2','SiO2','O2','H+'),
c(999,0,999,999,999,999,999,-7))
# we have to tell CHNOSZ which species to include; there are
# others that could be in this system
species(c('quartz','talc','forsterite','tremolite','diopside',
'wollastonite','monticellite','merwinite'))
# calculate the chemical affinities of formation reactions
#t <- affinity('Mg+2'=c(-14,-8),'Ca+2'=c(-12,-4),T=300,P=1000)
t <- affinity('Mg+2'=c(-12,-4),'Ca+2'=c(-8,0),T=300,P=1000)
diagram(t)
title(main=paste('HCl-H2O-CaO-CO2-MgO-(SiO2) at 300 deg C,<n>',
'1000 bar and pH=7. After Bowers et al., 1984'))
# note: BJH84 use a different method for representing
# the axes of the diagrams, similar to (a_Ca+2)/(a_H+)^2,
# so this is so far an approximate reproduction of their diagram.
## Fe-S-O at 200 deg C
## (After Helgeson, 1970)
basis(c('Fe','O2','S2'))
species(c('iron','ferrous-oxide','magnetite',
'hematite','pyrite','pyrrhotite'))
a <- affinity(S2=c(-50,0),O2=c(-90,-10),T=200)
diagram(a,color='heat')
title(main=paste('Fe-S-O, 200 degrees C, 1 bar',
'After Helgeson, 1970',sep='<n>'))
## Nucleobase - Amino Acid Interaction Eh-H2O
## association of amino acids with
## nucleotides in the genetic code.
# for this example we try a unique basis definition
basis(c('CO2','H2O','glutamine','e-','H+'),c(-3,0,-3,0,-7))
species(c('uracil','cytosine','adenine','guanine',
'phenylalanine','proline','lysine','glycine'),'aq')
# this loaded four nucleobases and four related amino acids
# (coded for by the homocodon triplets)
# check out the predominance diagrams
a.1 <- affinity(H2O=c(-5,0),Eh=c(-0.5,0))
diagram(a.1,color=NULL)
# overlay a different temperature
a.2 <- affinity(H2O=c(-5,0),Eh=c(-0.5,0),T=100)
diagram(a.2,col='red',add=TRUE,names=NULL)
# start make a title for the plot
tb <- thermo$basis # includes activities of basis species
# exclude those that are on the axes
tb <- tb[!((rownames(tb) %in% c('e-','H2O'))),]
title(main=paste('Nucleobases and amino acids,',
'T=25 and 100 C at Psat<n>',
describe(tb,T=NULL,P=NULL)),cex.main=0.9)
## Temperature-Pressure: kayanite-sillimanite-andalusite
# this is a system of a single component (Al2SiO5)
# but we are required to provide put more than one
# species in the basis definition
basis(c('kyanite','quartz','oxygen'))
species(c('kyanite','sillimanite','andalusite'))
diagram(affinity(T=c(0,1000,200),P=c(500,5000,200)),color=NULL)
title(main='Al2SiO5')</n>
<references>Aksu, S. and Doyle, F. M., 2001. Electrochemistry of copper in aqueous glycine solutions. <em>J. Electrochem. Soc.</em>, 148, B51-B57.
Bowers, T. S., Jackson, K. J. and Helgeson, H. C., 1984. <em>Equilibrium Activity Diagrams for Coexisting Minerals and Aqueous Solutions at Pressures and Temperatures to 5 kb and 600</em> $^{\circ}$<em>C</em>. Springer-Verlag, Heidelberg, 397 p.
Dick, J. M., 2008. Calculation of the relative metastabilities of proteins using the CHNOSZ software package. <em>Geochem. Trans.</em>, 9:10.
Ernst, W. G., 1976. <em>Petrologic Phase Equilibria</em>. W. H. Freeman, San Francisco, 333 p.
Helgeson, H. C., 1970. A chemical and thermodynamic model of ore deposition in hydrothermal systems. <em>Mineral. Soc. Amer. Spec. Pap.</em>, 3, 155-186.
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. <em>Am. J. Sci.</em>, 278-A, 1-229.
LaRowe, D. E. and Helgeson, H. C., 2007. Quantifying the energetics of metabolic reactions in diverse biogeochemical systems: electron flow and ATP synthesis. <em>Geobiology</em>, 5, 153-168.
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. <em>Eur. J. Mineral.</em>, 18, 175-186.
Pourbaix, M. J. N., 1949. <em>Thermodynamics of dilute aqueous solutions</em>. Edward Arnold & Co., London, 136 p.
Seewald, J. S., 1997. Mineral redox buffers and the stability of organic compounds under hydrothermal conditions. <em>Mat. Res. Soc. Symp. Proc.</em>, 432, 317-331.
Tagirov, B. and Schott, J., 2001. Aluminum speciation in crustal fluids revisited. <em>Geochim. Cosmochim. Acta</em>, 65, 3965-3992.
Tardy, Y., Schaul, R. and Duplay, J., 1997. Thermodynamic stability fields of humus, microflora and plants. <em>C. R. Acad. Sci. Paris</em>, 324, 969-976.</references>
<keyword>misc</keyword></n></n></n></n></n></n></n></n></nafter></n></n></n></n></n>Run the code above in your browser using DataLab