diagram: amino acid example

diagrm> ## Don't show: 
diagrm> data(thermo)
thermo$obigt has 1800 aqueous, 2925 total species

diagrm> ## End Don't show
diagrm> ## chemical affinities of formation (for equal activities)
diagrm> ## and equilibrium activities (for equal affinities) of amino acids
diagrm> opar <- par(mfrow=c(2,2))

diagrm> mycol <- rev(rainbow(5))

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

diagrm> # comment out next line to find predominance of methionine
diagrm> basis("H2S",-25)
    C H N O S Z ispecies logact state
CO2 1 0 0 2 0 0       69     -3    aq
H2O 0 2 0 1 0 0        1      0   liq
NH3 0 3 1 0 0 0       68     -4    aq
H2S 0 2 0 0 1 0       70    -25    aq
O2  0 0 0 2 0 0     2691    -80   gas
H+  0 1 0 0 0 1        3     -7    aq

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

diagrm> species(aa)
  CO2 H2O NH3 H2S   O2 H+ ispecies logact state          name
1   5   3   1   0 -4.5  0     1514     -3    aq glutamic acid
2   5   3   1   1 -5.5  0     1525     -3    aq    methionine
3   6   5   1   0 -7.5  0     1520     -3    aq    isoleucine
4   6   5   1   0 -7.5  0     1521     -3    aq       leucine
5   9   4   1   0 -9.5  0     1531     -3    aq      tyrosine

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

diagrm> diagram(a,what="A",col=mycol,yline=3,
diagrm+   main=paste("affinity of formation reactions per CO2\n",
diagrm+   "logfO2 = -80, logaH2O = 0"))
diagram: immobile component is CO2 
diagram: conservation coefficients are 5 5 6 6 9
Hit <Return> to see next plot: 

diagrm> # affinities of reactions as a function of oxygen fugacity
diagrm> a <- affinity(O2=c(-80,-66))
affinity: temperature is 25 C
energy.args: pressure is Psat
energy.args: variable 1 is O2 at 128 increments from -80 to -66 
subcrt: 11 species at 298.15 K and 1 bar (wet) 

diagrm> diagram(a,what="A",ylim=c(-15,5),col=mycol,lwd=2)
diagram: immobile component is CO2 
diagram: conservation coefficients are 5 5 6 6 9

diagrm> title("affinity per CO2; logaH2O = 0")

diagrm> ## logarithms of activities of amino acids
diagrm> # reactions balanced on CO2
diagrm> diagram(a,legend.x="bottomright",col=mycol,lwd=2)
diagram: immobile component is CO2 
diagram: conservation coefficients are 5 5 6 6 9
diagram: log total activity of CO2 (from species) is -1.508638

diagrm> title("equilibrium activities; logaH2O = 0")

diagrm> # in two dimensions
diagrm> a <- affinity(O2=c(-80,-66),H2O=c(-8,4))
affinity: temperature is 25 C
energy.args: pressure is Psat
energy.args: variable 1 is O2 at 128 increments from -80 to -66 
energy.args: variable 2 is H2O at 128 increments from -8 to 4 
subcrt: 11 species at 298.15 K and 1 bar (wet) 

diagrm> diagram(a,color=mycol)
diagram: immobile component is CO2 
diagram: conservation coefficients are 5 5 6 6 9

diagrm> title("highest equilibrium activities")

diagrm> # reset the plot device
diagrm> par(opar)