diagram: properties

diagrm>   ## Don't show:
diagrm> data(thermo)
thermo: loaded 1997 aqueous, 3089 total species to thermo$obigt
thermo: loaded 5264 proteins to thermo$ECO
thermo: loaded 6717 proteins to thermo$SGD
thermo: loaded 4155 localizations and 3570 abundances to thermo$yeastgfp

diagrm> ## End Don't show
diagrm>   ### 0-D property and logarithm of activity plots
diagrm>   ## properties of amino acids
diagrm>   basis("CHNOS")
    C H N O S ispecies logact state
CO2 1 0 0 2 0       69     -3    aq
H2O 0 2 0 1 0        1      0   liq
NH3 0 3 1 0 0       68     -4    aq
H2S 0 2 0 0 1       70     -7    aq
O2  0 0 0 2 0     2852    -80   gas

diagrm>   basis("H2S",-25)
    C H N O S ispecies logact state
CO2 1 0 0 2 0       69     -3    aq
H2O 0 2 0 1 0        1      0   liq
NH3 0 3 1 0 0       68     -4    aq
H2S 0 2 0 0 1       70    -25    aq
O2  0 0 0 2 0     2852    -80   gas

diagrm>   aa <- c("glutamic acid","methionine","isoleucine",
diagrm+     "leucine","tyrosine")

diagrm>   species(aa)

diagrm>   # Gibbs energy of formation reactions per CO2
diagrm>   a <- affinity(property="G")
affinity: temperature is 25 C
energy.args: pressure is Psat
subcrt: 10 species at 298.15 K and 1 bar (wet)

diagrm>   diagram(a,balance=1)
diagram: conservation coefficients are 1 1 1 1 1

diagrm>   # affinity of reactions per CO2
diagrm>   a <- affinity()
affinity: temperature is 25 C
energy.args: pressure is Psat
subcrt: 10 species at 298.15 K and 1 bar (wet)

diagrm>   diagram(a,property="A",do.plot=FALSE)
diagram: immobile component is CO2
diagram: conservation coefficients are 5 5 6 6 9

diagrm>   ## logarithms of activities of amino acids
diagrm>   # balanced on amino acid backbone, using abundance.old
diagrm>   d1 <- diagram(a,balance=1,do.plot=FALSE)
diagram: conservation coefficients are 1 1 1 1 1
diagram: log total activity of species (from species) is -2.30103

diagrm>   # balanced on amino acid backbone, using abundance.new
diagrm>   d2 <- diagram(a,balance=1,residue=TRUE,do.plot=FALSE)
diagram: conservation coefficients are 1 1 1 1 1
diagram: using residue equivalents
diagram: log total activity of species (from species) is -2.30103

diagrm>   # the last two give the same results
diagrm>   stopifnot(all.equal(as.numeric(d2$logact),
diagrm+     as.numeric(d1$logact)))

diagrm>   ## 1-D plots
diagrm>   a <- affinity(O2=c(-80,-70))
affinity: temperature is 25 C
energy.args: pressure is Psat
energy.args: variable 1 is O2 at 128 increments from -80 to -70
subcrt: 10 species at 298.15 K and 1 bar (wet)

diagrm>   # affinities of formation reactions
diagrm>   d <- diagram(a,property="A",balance=1,ylim=c(-80,30))
diagram: conservation coefficients are 1 1 1 1 1

diagrm>   ### 1-D property plot
diagrm>   ## for hydrous cordierite = cordierite + H2O
diagrm>   ## after Helgeson et al., 1978
diagrm>   basis(c("cordierite,hydrous","Mg+2","SiO2","H2O","O2","H+"))
basis: changed basis to Mg2Al3(AlSi5)O18*H2O Mg+2 SiO2 H2O O2 H+.
                     Al H Mg  O Si Z ispecies logact state
Mg2Al3(AlSi5)O18*H2O  4 2  2 19  5 0     2071      0    cr
Mg+2                  0 0  1  0  0 2        9      0    aq
SiO2                  0 0  0  2  1 0       72      0    aq
H2O                   0 2  0  1  0 0        1      0   liq
O2                    0 0  0  2  0 0     2852      0   gas
H+                    0 1  0  0  0 1        3      0    aq

diagrm>   species("cordierite")

diagrm>   # water.SUPCRT92 can only get us up to 5000 bar
diagrm>   # (7000 and 10000 bar shown in the original)
diagrm>   P <- c(1,2,3,5)*1000

diagrm>   col <- rainbow(length(P))

diagrm>   for(i in 1:length(P)) {
diagrm+     a <- affinity(property="logK",T=c(20,800),P=P[i])
diagrm+     diagram(a,add=(i!=1),ylim=c(-4,2),legend.x=NULL,
diagrm+       color=col[i],title="")
diagrm+   }
affinity: pressure is 1000 bar
energy.args: variable 1 is T at 128 increments from 293.15 to 1073.15
subcrt: 7 species at 128 values of T and P (wet)
diagram: immobile component is Mg2Al3(AlSi5)O18*H2O
diagram: conservation coefficients are 1
affinity: pressure is 2000 bar
energy.args: variable 1 is T at 128 increments from 293.15 to 1073.15
subcrt: 7 species at 128 values of T and P (wet)
diagram: immobile component is Mg2Al3(AlSi5)O18*H2O
diagram: conservation coefficients are 1
affinity: pressure is 3000 bar
energy.args: variable 1 is T at 128 increments from 293.15 to 1073.15
subcrt: 7 species at 128 values of T and P (wet)
diagram: immobile component is Mg2Al3(AlSi5)O18*H2O
diagram: conservation coefficients are 1
affinity: pressure is 5000 bar
energy.args: variable 1 is T at 128 increments from 293.15 to 1073.15
subcrt: 7 species at 128 values of T and P (wet)
diagram: immobile component is Mg2Al3(AlSi5)O18*H2O
diagram: conservation coefficients are 1

diagrm>   legend("topright",lty=1,col=col,legend=paste(P,"bar"))

diagrm>   title(main=paste("hydrous cordierite = cordierite + H2O",
diagrm+     "After Helgeson et al., 1978",sep="\n"),cex.main=0.9)

diagrm>   # this calculation could also be done through the
diagrm>   # subcrt function, without having to set up basis species
diagrm>
diagrm>