next up previous
Next: findit - see demos Up: CHNOSZ examples Previous: [1] objective

[3] revisit

## Don't show: data(thermo)
thermo$obigt: 1809 aqueous, 3368 total species
## End(Don't show) ## example of defining a new objective function # count the species with logarithms of activity greater than loga2 count <- function(loga1, loga2) rowSums(loga1 > loga2) # set the attribute indicating the type of optimum attr(count, "optimum") <- "maximal" # equilibrate a system of amino acids 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 3095 -80 gas
species(aminoacids(""))
CO2 H2O NH3 H2S O2 ispecies logact state name 1 3 2 1 0 -3.0 1504 -3 aq alanine 2 3 1 1 1 -2.5 1511 -3 aq cysteine 3 4 2 1 0 -3.0 1509 -3 aq aspartic acid 4 5 3 1 0 -4.5 1514 -3 aq glutamic acid 5 9 4 1 0 -10.0 1526 -3 aq phenylalanine 6 2 1 1 0 -1.5 1516 -3 aq glycine 7 6 0 3 0 -5.0 1518 -3 aq histidine 8 6 5 1 0 -7.5 1520 -3 aq isoleucine 9 6 4 2 0 -7.0 1522 -3 aq lysine 10 6 5 1 0 -7.5 1521 -3 aq leucine 11 5 3 1 1 -5.5 1525 -3 aq methionine 12 4 1 2 0 -3.0 1508 -3 aq asparagine 13 5 3 1 0 -5.5 1527 -3 aq proline 14 5 2 2 0 -4.5 1513 -3 aq glutamine 15 6 1 4 0 -5.5 1505 -3 aq arginine 16 3 2 1 0 -2.5 1528 -3 aq serine 17 4 3 1 0 -4.0 1529 -3 aq threonine 18 5 4 1 0 -6.0 1533 -3 aq valine 19 11 3 2 0 -11.5 1530 -3 aq tryptophan 20 9 4 1 0 -9.5 1531 -3 aq tyrosine
a <- affinity(O2=c(-80, -60))
energy.args: temperature is 25 C energy.args: pressure is Psat energy.args: variable 1 is log_f(O2) at 128 values from -80 to -60 subcrt: 25 species at 298.15 K and 1 bar (wet)
e <- equilibrate(a)
balance: coefficients are moles of CO2 in formation reactions equilibrate: balancing coefficients are 3 3 4 5 9 2 6 6 6 6 5 4 5 5 6 3 4 5 11 9 equilibrate: logarithm of total moles of CO2 is -0.97061622231479
# make a plot r <- revisit(e, "count", -5)
revisit: calculating count in 1 dimensions
title(main="amino acids with metastable log activities > -5")

Image revisit1

 

# can also make a 2-D plot a <- affinity(O2=c(-74, -60, 64), H2O=c(-3, 3, 64))
energy.args: temperature is 25 C energy.args: pressure is Psat energy.args: variable 1 is log_f(O2) at 64 values from -74 to -60 energy.args: variable 2 is log_a(H2O) at 64 values from -3 to 3 subcrt: 25 species at 298.15 K and 1 bar (wet)
e <- equilibrate(a)
balance: coefficients are moles of CO2 in formation reactions equilibrate: balancing coefficients are 3 3 4 5 9 2 6 6 6 6 5 4 5 5 6 3 4 5 11 9 equilibrate: logarithm of total moles of CO2 is -0.97061622231479
r <- revisit(e, "count", -5, style.2D="image", plot.optval=FALSE)
revisit: calculating count in 2 dimensions
title(main="amino acids with metastable log activities > -5")

Image revisit2

 

## 'revisit' calculations for amino acids opar <- par(mfrow=c(2, 2)) 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 3095 -80 gas H+ 0 1 0 0 0 1 3 -7 aq
species(aminoacids(""))
CO2 H2O NH3 H2S O2 H+ ispecies logact state name 1 3 2 1 0 -3.0 0 1504 -3 aq alanine 2 3 1 1 1 -2.5 0 1511 -3 aq cysteine 3 4 2 1 0 -3.0 0 1509 -3 aq aspartic acid 4 5 3 1 0 -4.5 0 1514 -3 aq glutamic acid 5 9 4 1 0 -10.0 0 1526 -3 aq phenylalanine 6 2 1 1 0 -1.5 0 1516 -3 aq glycine 7 6 0 3 0 -5.0 0 1518 -3 aq histidine 8 6 5 1 0 -7.5 0 1520 -3 aq isoleucine 9 6 4 2 0 -7.0 0 1522 -3 aq lysine 10 6 5 1 0 -7.5 0 1521 -3 aq leucine 11 5 3 1 1 -5.5 0 1525 -3 aq methionine 12 4 1 2 0 -3.0 0 1508 -3 aq asparagine 13 5 3 1 0 -5.5 0 1527 -3 aq proline 14 5 2 2 0 -4.5 0 1513 -3 aq glutamine 15 6 1 4 0 -5.5 0 1505 -3 aq arginine 16 3 2 1 0 -2.5 0 1528 -3 aq serine 17 4 3 1 0 -4.0 0 1529 -3 aq threonine 18 5 4 1 0 -6.0 0 1533 -3 aq valine 19 11 3 2 0 -11.5 0 1530 -3 aq tryptophan 20 9 4 1 0 -9.5 0 1531 -3 aq tyrosine
# chemical affinities as a function of logarithm of oxygen fugacity a <- affinity(O2=c(-85, -60))
energy.args: temperature is 25 C energy.args: pressure is Psat energy.args: variable 1 is log_f(O2) at 128 values from -85 to -60 subcrt: 26 species at 298.15 K and 1 bar (wet)
# shows the equilibrium abundances of the amino acids e <- equilibrate(a)
balance: coefficients are moles of CO2 in formation reactions equilibrate: balancing coefficients are 3 3 4 5 9 2 6 6 6 6 5 4 5 5 6 3 4 5 11 9 equilibrate: logarithm of total moles of CO2 is -0.97061622231479
diagram(e) mtitle(c("20 amino acids", "balanced on CO2")) # show a legend with input constraints db <- describe.basis(ibasis=3) dp <- describe.property("T", 25) legend("bottomright", c(dp, db)) # default is to plot coefficient of variation r <- revisit(e)
revisit: calculating CV in 1 dimensions
# show a title with the optimal conditions mincv <- format(r$optimum, digits=3) t1 <- paste("minimum coeff of variation,", mincv, "at:") # the logfO2 that minimized the C.V. basis("O2", r$x)
C H N O S Z ispecies logact state CO2 1 0 0 2 0 0 69 -3.00000 aq H2O 0 2 0 1 0 0 1 0.00000 liq NH3 0 3 1 0 0 0 68 -4.00000 aq H2S 0 2 0 0 1 0 70 -7.00000 aq O2 0 0 0 2 0 0 3095 -67.67717 gas H+ 0 1 0 0 0 1 3 -7.00000 aq
t2 <- describe.basis(ibasis=5) mtitle(c(t1, as.expression(t2))) # chemical affinities as a function of two other variables a <- affinity(NH3=c(-10, 10, 40), T=c(0, 80, 40))
energy.args: pressure is Psat energy.args: variable 1 is log_a(NH3) at 40 values from -10 to 10 energy.args: variable 2 is T at 40 values from 273.15 to 353.15 K subcrt: 26 species at 40 values of T and P (wet)
diagram(a, fill="heat")
balance: coefficients are moles of CO2 in formation reactions diagram: plotting A/2.303RT from affinity(), divided by balancing coefficients
# show a legend with input constraints db <- describe.basis(ibasis=5) legend("bottomright", as.expression(db)) # contour plot of the CV e <- equilibrate(a)
balance: coefficients are moles of CO2 in formation reactions equilibrate: balancing coefficients are 3 3 4 5 9 2 6 6 6 6 5 4 5 5 6 3 4 5 11 9 equilibrate: logarithm of total moles of CO2 is -0.97061622231479
r <- revisit(e)
revisit: calculating CV in 2 dimensions
# show a title with the optimal conditions mincv <- format(r$optimum, digits=3) t1 <- paste("minimum coeff of variation,", mincv, "at:") # the logaNH3 and T that minimized the C.V. basis("NH3", r$x)
C H N O S Z ispecies logact state CO2 1 0 0 2 0 0 69 -3.000000 aq H2O 0 2 0 1 0 0 1 0.000000 liq NH3 0 3 1 0 0 0 68 2.820513 aq H2S 0 2 0 0 1 0 70 -7.000000 aq O2 0 0 0 2 0 0 3095 -67.677165 gas H+ 0 1 0 0 0 1 3 -7.000000 aq
db <- describe.basis(ibasis=3) dp <- describe.property("T", r$y) t2 <- substitute(list(dp, db), list(dp=dp[[1]], db=db[[1]])) mtitle(c(t1, as.expression(t2))) par(opar)

Image revisit3

 


next up previous
Next: findit - see demos Up: CHNOSZ examples Previous: [1] objective