Nothing
inst/doc/anintro.R
In CHNOSZ: Thermodynamic Calculations and Diagrams for Geochemistry
## ----options, include=FALSE---------------------------------------------------
options(width = 80)
options(digits = 6)
## ----HTML, include=FALSE------------------------------------------------------
## Some frequently used HTML expressions
logfO2 <- "log<i>f</i><sub>O<sub>2</sub></sub>"
# Use lowercase here because these tend to be variable names in the examples
zc <- "<i>Z</i><sub>C</sub>"
o2 <- "O<sub>2</sub>"
h2o <- "H<sub>2</sub>O"
sio2 <- "SiO<sub>2</sub>"
ch4 <- "CH<sub>4</sub>"
## ----setup, include=FALSE-----------------------------------------------------
library(knitr)
## From "Tufte Handout" example dated 2016-12-27
# Invalidate cache when the tufte version changes
opts_chunk$set(tidy = FALSE, cache.extra = packageVersion('tufte'))
options(htmltools.dir.version = FALSE)
## Adjust plot margins
## First one from https://yihui.name/knitr/hooks/
knit_hooks$set(small.mar = function(before, options, envir) {
if (before) par(mar = c(4.2, 4.2, .1, .1)) # smaller margin on top and right
})
knit_hooks$set(tiny.mar = function(before, options, envir) {
if (before) par(mar = c(.1, .1, .1, .1)) # tiny margin all around
})
knit_hooks$set(smallish.mar = function(before, options, envir) {
if (before) par(mar = c(4.2, 4.2, 0.9, 0.9)) # smallish margins on top and right
})
## Use pngquant to optimize PNG images
knit_hooks$set(pngquant = hook_pngquant)
pngquant <- "--speed=1 --quality=0-25"
# pngquant isn't available on R-Forge ...
if (!nzchar(Sys.which("pngquant"))) pngquant <- NULL
# Set dpi 20231129
knitr::opts_chunk$set(
dpi = if(nzchar(Sys.getenv("CHNOSZ_BUILD_LARGE_VIGNETTES"))) 72 else 50
)
hidpi = if(nzchar(Sys.getenv("CHNOSZ_BUILD_LARGE_VIGNETTES"))) 100 else 50
## http://stackoverflow.com/questions/23852753/knitr-with-gridsvg
## Set up a chunk hook for manually saved plots.
knit_hooks$set(custom.plot = hook_plot_custom)
## Hook to change <img /> to <embed /> -- required for interactive SVG
hook_plot <- knit_hooks$get("plot")
knit_hooks$set(plot = function(x, options) {
x <- hook_plot(x, options)
if (!is.null(options$embed.tag) && options$embed.tag) x <- gsub("<img ", "<embed ", x)
x
})
## http://stackoverflow.com/questions/30530008/hook-to-time-knitr-chunks
now = Sys.time()
knit_hooks$set(timeit = function(before) {
if (before) { now <<- Sys.time() }
else {
paste("%", sprintf("Chunk rendering time: %s seconds.\n", round(Sys.time() - now, digits = 3)))
}
})
timeit <- NULL
## Colorize messages 20171031
## Adapted from https://gist.github.com/yihui/2629886#file-knitr-color-msg-rnw
color_block = function(color) {
function(x, options) sprintf('<pre style="color:%s">%s</pre>', color, x)
}
knit_hooks$set(warning = color_block('magenta'), error = color_block('red'), message = color_block('blue'))
## ----install_CHNOSZ, eval=FALSE-----------------------------------------------
# install.packages("CHNOSZ")
## ----library_CHNOSZ-----------------------------------------------------------
library(CHNOSZ)
## ----reset--------------------------------------------------------------------
reset()
## ----pseudocode, eval=FALSE---------------------------------------------------
# basis(...)
# species(...)
# a <- affinity(...)
# e <- equilibrate(a) ## optional
# diagram(e) ## or diagram(a)
# reset() ## clear settings for next calculation
## ----info_CH4-----------------------------------------------------------------
info("CH4")
## ----info_CH4_gas-------------------------------------------------------------
info("CH4", "gas")
## ----info_names_gas-----------------------------------------------------------
info("methane")
info("oxygen")
info("carbon dioxide")
## ----info_S_S2----------------------------------------------------------------
info("S")
info("S2")
## ----iCH4, message=FALSE------------------------------------------------------
iCH4 <- info("CH4")
info(iCH4)
## ----info_info_water----------------------------------------------------------
info(info("water"))
## ----width180, include=FALSE------------------------------------------------------------------------------------------------------------------------------------------------------
options(width = 180)
## ----info_acid--------------------------------------------------------------------------------------------------------------------------------------------------------------------
info("acid")
## ----width80, include=FALSE---------------------------------------------------
options(width = 80)
## ----info_ribose--------------------------------------------------------------
info(" ribose")
## ----info_CH4_formula, message=FALSE------------------------------------------
info(iCH4)$formula
## ----makeup_iCH4--------------------------------------------------------------
makeup(iCH4)
as.chemical.formula(makeup(iCH4))
## ----ZC_iCH4, message=FALSE---------------------------------------------------
ZC(iCH4)
ZC(info(iCH4)$formula)
ZC(makeup(iCH4))
## ----subcrt_water-------------------------------------------------------------
subcrt("water")
## ----subcrt_water_grid--------------------------------------------------------
subcrt("water", T = c(400, 500, 600), P = c(200, 400, 600), grid = "P")$out$water
## ----subcrt_water_plot, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, out.width="100%", echo=FALSE, message=FALSE, fig.cap="Isothermal contours of density (g cm<sup>-3</sup>) and pressure (bar) of water.", cache=TRUE, pngquant=pngquant, timeit=timeit----
sres <- subcrt("water", T=seq(0,1000,100), P=c(NA, seq(1,500,1)), grid="T")
water <- sres$out$water
plot(water$P, water$rho, type = "l")
## ----subcrt_water_plot, eval=FALSE--------------------------------------------
# sres <- subcrt("water", T=seq(0,1000,100), P=c(NA, seq(1,500,1)), grid="T")
# water <- sres$out$water
# plot(water$P, water$rho, type = "l")
## ----units_CH4----------------------------------------------------------------
T.units("K")
P.units("MPa")
subcrt("CH4", T = 298.15, P = 0.1)$out$CH4$G
E.units("cal")
subcrt("CH4", T = 298.15, P = 0.1)$out$CH4$G
## ----convert_G, message=FALSE-------------------------------------------------
(CH4dat <- info(info("CH4")))
convert(CH4dat$G, "J")
## ----reset--------------------------------------------------------------------
reset()
## ----subcrt_CO2---------------------------------------------------------------
subcrt(c("CO2", "CO2"), c("gas", "aq"), c(-1, 1), T = seq(0, 250, 50))
## ----CO2_logK, echo=FALSE, message=FALSE--------------------------------------
T <- seq(0, 350, 10)
CO2 <- subcrt(c("CO2", "CO2"), c("gas", "aq"), c(-1, 1), T = T)$out$logK
CO <- subcrt(c("CO", "CO"), c("gas", "aq"), c(-1, 1), T = T)$out$logK
CH4 <- subcrt(c("CH4", "CH4"), c("gas", "aq"), c(-1, 1), T = T)$out$logK
logK <- data.frame(T, CO2, CO, CH4)
## ----CO2_plot, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, out.width="100%", echo=FALSE, message=FALSE, fig.cap="Calculated equilibrium constants for dissolution of CO<sub>2</sub>, CO, and CH<sub>4</sub>.", cache=TRUE, pngquant=pngquant, timeit=timeit----
matplot(logK[, 1], logK[, -1], type = "l", col = 1, lty = 1,
xlab = axis.label("T"), ylab = axis.label("logK"))
text(80, -1.7, expr.species("CO2"))
text(240, -2.37, expr.species("CO"))
text(300, -2.57, expr.species("CH4"))
## ----CO2_logK, eval=FALSE-----------------------------------------------------
# T <- seq(0, 350, 10)
# CO2 <- subcrt(c("CO2", "CO2"), c("gas", "aq"), c(-1, 1), T = T)$out$logK
# CO <- subcrt(c("CO", "CO"), c("gas", "aq"), c(-1, 1), T = T)$out$logK
# CH4 <- subcrt(c("CH4", "CH4"), c("gas", "aq"), c(-1, 1), T = T)$out$logK
# logK <- data.frame(T, CO2, CO, CH4)
## ----CO2_plot, eval=FALSE-----------------------------------------------------
# matplot(logK[, 1], logK[, -1], type = "l", col = 1, lty = 1,
# xlab = axis.label("T"), ylab = axis.label("logK"))
# text(80, -1.7, expr.species("CO2"))
# text(240, -2.37, expr.species("CO"))
# text(300, -2.57, expr.species("CH4"))
## ----subcrt_unbalanced, results="hide"----------------------------------------
subcrt(c("CO2", "CH4"), c(-1, 1))
## ----basis_singular, error=TRUE-----------------------------------------------
basis(c("CO2", "H2", "H2CO2"))
## ----basis_CHO----------------------------------------------------------------
basis(c("CO2", "H2", "H2O"))
## ----basis_CHOZ---------------------------------------------------------------
basis(c("CO2", "H2", "H2O", "H+"))
## ----subcrt_acetoclastic, message=FALSE---------------------------------------
subcrt(c("acetate", "CH4"), c(-1, 1))$reaction
## ----subcrt_methanogenesis, message=FALSE-------------------------------------
acetate_oxidation <- subcrt("acetate", -1)
hydrogenotrophic <- subcrt("CH4", 1)
acetoclastic <- subcrt(c("acetate", "CH4"), c(-1, 1))
## ----describe_reaction_plot, fig.margin=TRUE, fig.width=3.5, fig.height=1.8, tiny.mar=TRUE, out.width="100%", pngquant=pngquant, timeit=timeit----
plot(0, 0, type = "n", axes = FALSE, ann=FALSE, xlim=c(0, 5), ylim=c(5.2, -0.2))
text(0, 0, "acetoclastic methanogenesis", adj = 0)
text(5, 1, describe.reaction(acetoclastic$reaction), adj = 1)
text(0, 2, "acetate oxidation", adj = 0)
text(5, 3, describe.reaction(acetate_oxidation$reaction), adj = 1)
text(0, 4, "hydrogenotrophic methanogenesis", adj = 0)
text(5, 5, describe.reaction(hydrogenotrophic$reaction), adj = 1)
## ----basis_mayumi, message=FALSE, results="hide"------------------------------
basis(c("CO2", "H2", "H2O", "H+"))
basis(c("CO2", "H2"), "gas")
basis(c("H2", "pH"), c(-3.92, 7.3))
## ----affinity_acetoclastic, message=FALSE-------------------------------------
subcrt(c("acetate", "CH4"), c(-1, 1),
c("aq", "gas"), logact = c(-3.4, -0.18), T = 55, P = 50)$out
## ----affinity_hydrogenotrophic, message=FALSE---------------------------------
subcrt("CH4", 1, "gas", logact = -0.18, T = 55, P = 50)$out
## ----rxnfun, message=FALSE----------------------------------------------------
rxnfun <- function(coeffs) {
subcrt(c("acetate", "CH4"), coeffs,
c("aq", "gas"), logact = c(-3.4, -0.18), T = 55, P = 50)$out
}
## ----methanogenesis_plot, fig.margin=TRUE, fig.width=4.1, fig.height=4.1, small.mar=TRUE, out.width="100%", echo=FALSE, message=FALSE, fig.cap="Gibbs energies of acetate oxidation and methanogenesis (after Mayumi et al., 2013).", cache=TRUE, pngquant=pngquant, timeit=timeit----
Adat <- lapply(c(-3, 3), function(logfCO2) {
basis("CO2", logfCO2)
data.frame(logfCO2,
rxnfun(c(0, 0))$A,
rxnfun(c(-1, 0))$A,
rxnfun(c(-1, 1))$A,
rxnfun(c(0, 1))$A
)
})
Adat <- do.call(rbind, Adat)
matplot(Adat[, 1], -Adat[, -1]/1000, type = "l", lty = 1, lwd = 2,
xlab = axis.label("CO2"), ylab = axis.label("DG", prefix = "k"))
legend("topleft", c("acetate oxidation", "acetoclastic methanogenesis",
"hydrogenotrophic methanogenesis"), lty = 1, col = 2:4)
## ----methanogenesis_plot, eval=FALSE------------------------------------------
# Adat <- lapply(c(-3, 3), function(logfCO2) {
# basis("CO2", logfCO2)
# data.frame(logfCO2,
# rxnfun(c(0, 0))$A,
# rxnfun(c(-1, 0))$A,
# rxnfun(c(-1, 1))$A,
# rxnfun(c(0, 1))$A
# )
# })
# Adat <- do.call(rbind, Adat)
# matplot(Adat[, 1], -Adat[, -1]/1000, type = "l", lty = 1, lwd = 2,
# xlab = axis.label("CO2"), ylab = axis.label("DG", prefix = "k"))
# legend("topleft", c("acetate oxidation", "acetoclastic methanogenesis",
# "hydrogenotrophic methanogenesis"), lty = 1, col = 2:4)
## ----reset, message=FALSE-----------------------------------------------------
reset()
## ----basis_CHNOSZ, results="hide"---------------------------------------------
basis("CHNOSe")
## ----species_sulfur-----------------------------------------------------------
species(c("H2S", "HS-", "HSO4-", "SO4-2"))
## ----affinity-----------------------------------------------------------------
unlist(affinity()$values)
## ----EhpH_plot, fig.margin=TRUE, fig.width=4, fig.height=4, out.width="100%", echo = FALSE, message=FALSE, cache=TRUE, fig.cap="Aqueous sulfur species at 25 °C.", pngquant=pngquant, timeit=timeit----
a <- affinity(pH = c(0, 12), Eh = c(-0.5, 1))
diagram(a, limit.water = TRUE)
## ----EhpH_plot, echo=TRUE, eval=FALSE-----------------------------------------
# a <- affinity(pH = c(0, 12), Eh = c(-0.5, 1))
# diagram(a, limit.water = TRUE)
## ----EhpH_plot_color, fig.margin=TRUE, fig.width=4, fig.height=4, smallish.mar=TRUE, out.width="100%", echo=FALSE, message=FALSE, cache=TRUE, fig.cap="The same plot, with different colors and labels.", pngquant=pngquant, timeit=timeit----
diagram(a, fill = "terrain", lwd = 2, lty = 3,
names = c("hydrogen sulfide", "bisulfide", "bisulfate", "sulfate"),
las = 0)
water.lines(a, col = 6, lwd = 2)
## ----EhpH_plot_color, echo=TRUE, eval=FALSE-----------------------------------
# diagram(a, fill = "terrain", lwd = 2, lty = 3,
# names = c("hydrogen sulfide", "bisulfide", "bisulfate", "sulfate"),
# las = 0)
# water.lines(a, col = 6, lwd = 2)
## ----retrieve-----------------------------------------------------------------
retrieve("Mn", c("O", "H"), "aq")
retrieve("Mn", c("O", "H"), "cr")
## ----retrieve_diagram, fig.margin=TRUE, fig.width=5, fig.height=5, out.width="100%", message=FALSE, results = "hide", cache=TRUE, fig.cap="Eh-pH diagram for the Mn-O-H system.", pngquant=pngquant, timeit=timeit----
# Set decimal logarithm of activity of aqueous species,
# temperature and plot resolution
logact <- -4
T <- 100
res <- 400
# Start with the aqueous species
basis(c("Mn+2", "H2O", "H+", "e-"))
iaq <- retrieve("Mn", c("O", "H"), "aq")
species(iaq, logact)
aaq <- affinity(pH = c(4, 16, res), Eh = c(-1.5, 1.5, res), T = T)
# Show names for only the metastable species here
names <- names(iaq)
names[!names(iaq) %in% c("MnOH+", "MnO", "HMnO2-")] <- ""
diagram(aaq, lty = 2, col = "#4169E188", names = names, col.names = 4)
# Overlay mineral stability fields
icr <- retrieve("Mn", c("O", "H"), "cr")
species(icr, add = TRUE)
# Supply the previous result from affinity() to use
# argument recall (for plotted variables and T)
acr <- affinity(aaq)
diagram(acr, add = TRUE, bold = acr$species$state=="cr", limit.water = FALSE)
# Add legend
legend <- c(
bquote(log * italic(a)["Mn(aq)"] == .(logact)),
bquote(italic(T) == .(T) ~ degree*C)
)
legend("topright", legend = as.expression(legend), bty = "n")
## ----info_CuCl, results="hide"------------------------------------------------
info(" CuCl")
## ----copper_setup, echo=TRUE, results="hide"----------------------------------
basis(c("Cu", "H2S", "Cl-", "H2O", "H+", "e-"))
basis("H2S", -6)
basis("Cl-", -0.7)
species(c("CuCl", "CuCl2-", "CuCl3-2", "CuCl+", "CuCl2", "CuCl3-", "CuCl4-2"))
species(c("chalcocite", "tenorite", "cuprite", "copper"), add = TRUE)
## ----info_chalcocite, message=FALSE-------------------------------------------
info(info("chalcocite", c("cr", "cr2", "cr3")))$T
## ----copper_mosaic, fig.margin=TRUE, fig.width=4, fig.height=4, out.width="100%", message=FALSE, cache=TRUE, fig.cap="Copper minerals and aqueous complexes with chloride, 200 °C.", pngquant=pngquant, timeit=timeit----
T <- 200
res <- 200
bases <- c("H2S", "HS-", "HSO4-", "SO4-2")
m1 <- mosaic(bases, pH = c(0, 12, res), Eh=c(-1.2, 0.75, res), T=T)
diagram(m1$A.species, lwd = 2)
diagram(m1$A.bases, add = TRUE, col = 4, col.names = 4, lty = 3,
italic = TRUE)
water.lines(m1$A.species, col = "blue1")
## ----rainbow_data-------------------------------------------------------------
file <- system.file("extdata/cpetc/SC10_Rainbow.csv", package = "CHNOSZ")
rb <- read.csv(file, check.names = FALSE)
## ----rainbow_species, results="hide"------------------------------------------
basis(c("CO2", "H2", "NH4+", "H2O", "H2S", "H+"))
species("CH4", -3)
species(c("adenine", "cytosine", "aspartic acid", "deoxyribose",
"CH4", "leucine", "tryptophan", "n-nonanoic acid"), -6)
## ----rainbow_affinity, message=FALSE------------------------------------------
a <- affinity(T = rb$T, CO2 = rb$CO2, H2 = rb$H2,
`NH4+` = rb$`NH4+`, H2S = rb$H2S, pH = rb$pH)
T <- convert(a$vals[[1]], "K")
a$values <- lapply(a$values, convert, "G", T)
a$values <- lapply(a$values, convert, "cal")
a$values <- lapply(a$values, `*`, -0.001)
## ----rainbow_diagram, fig.margin=TRUE, fig.width=4, fig.height=4, out.width="100%", echo=FALSE, message=FALSE, cache=TRUE, fig.cap="Affinities of organic synthesis in a hydrothermal system, after Shock and Canovas (2010).", pngquant=pngquant, timeit=timeit----
diagram(a, balance = 1, ylim = c(-100, 100), ylab = quote(italic(A)*", kcal/mol"),
col = rainbow(8), lwd = 2, bg = "slategray3")
abline(h = 0, lty = 2, lwd = 2)
## ----rainbow_diagram, eval=FALSE----------------------------------------------
# diagram(a, balance = 1, ylim = c(-100, 100), ylab = quote(italic(A)*", kcal/mol"),
# col = rainbow(8), lwd = 2, bg = "slategray3")
# abline(h = 0, lty = 2, lwd = 2)
## ----PPM_basis, results="hide", message=FALSE---------------------------------
basis(c("FeS2", "H2S", "O2", "H2O"))
species(c("pyrite", "magnetite"))
species("pyrrhotite", "cr2", add = TRUE)
## ----PPM_affinity, message=FALSE----------------------------------------------
unlist(affinity(T = 300, P = 100)$values)
## ----PPM_setup, results="hide"------------------------------------------------
mod.buffer("PPM", "pyrrhotite", "cr2")
basis(c("H2S", "O2"), c("PPM", "PPM"))
## ----PPM_activities, message=FALSE--------------------------------------------
unlist(affinity(T = 300, P = 100, return.buffer = TRUE)[1:3])
## ----demo_buffer_noecho, fig.margin=TRUE, fig.width=4, fig.height=4, out.width="100%", message=FALSE, echo=FALSE, cache=TRUE, fig.cap="Values of log<i>f</i><sub>H<sub>2</sub></sub> corresponding to mineral buffers or to given activities of aqueous species.", pngquant=pngquant, timeit=timeit----
demo(buffer, echo = FALSE)
## ----PPM_basis, echo = FALSE, results = "hide"--------------------------------
basis(c("FeS2", "H2S", "O2", "H2O"))
species(c("pyrite", "magnetite"))
species("pyrrhotite", "cr2", add = TRUE)
## ----PPM_setup, echo = FALSE, results = "hide", message = FALSE---------------
mod.buffer("PPM", "pyrrhotite", "cr2")
basis(c("H2S", "O2"), c("PPM", "PPM"))
## ----PPM_affinity, message = FALSE--------------------------------------------
unlist(affinity(T = 300, P = 100)$values)
## ----demo_buffer, eval=FALSE--------------------------------------------------
# demo(buffer)
## ----bjerrum_diagram, fig.margin=TRUE, fig.width=3, fig.height=6, out.width="100%", echo=FALSE, results="hide", message=FALSE, cache=TRUE, fig.cap="Three views of carbonate speciation: affinity, activity, degree of formation.", pngquant=pngquant, timeit=timeit----
par(mfrow = c(3, 1))
basis("CHNOS+")
species(c("CO2", "HCO3-", "CO3-2"))
a25 <- affinity(pH = c(4, 13))
a150 <- affinity(pH = c(4, 13), T = 150)
diagram(a25, dy = 0.4)
diagram(a150, add = TRUE, names = FALSE, col = "red")
e25 <- equilibrate(a25, loga.balance = -3)
e150 <- equilibrate(a150, loga.balance = -3)
diagram(e25, ylim = c(-6, 0), dy = 0.15)
diagram(e150, add = TRUE, names = FALSE, col = "red")
diagram(e25, alpha = TRUE, dy = -0.25)
diagram(e150, alpha = TRUE, add = TRUE, names = FALSE, col = "red")
## ----bjerrum_diagram, echo=1:7, eval=FALSE------------------------------------
# par(mfrow = c(3, 1))
# basis("CHNOS+")
# species(c("CO2", "HCO3-", "CO3-2"))
# a25 <- affinity(pH = c(4, 13))
# a150 <- affinity(pH = c(4, 13), T = 150)
# diagram(a25, dy = 0.4)
# diagram(a150, add = TRUE, names = FALSE, col = "red")
# e25 <- equilibrate(a25, loga.balance = -3)
# e150 <- equilibrate(a150, loga.balance = -3)
# diagram(e25, ylim = c(-6, 0), dy = 0.15)
# diagram(e150, add = TRUE, names = FALSE, col = "red")
# diagram(e25, alpha = TRUE, dy = -0.25)
# diagram(e150, alpha = TRUE, add = TRUE, names = FALSE, col = "red")
## ----bjerrum_diagram, echo=8:11, eval=FALSE-----------------------------------
# par(mfrow = c(3, 1))
# basis("CHNOS+")
# species(c("CO2", "HCO3-", "CO3-2"))
# a25 <- affinity(pH = c(4, 13))
# a150 <- affinity(pH = c(4, 13), T = 150)
# diagram(a25, dy = 0.4)
# diagram(a150, add = TRUE, names = FALSE, col = "red")
# e25 <- equilibrate(a25, loga.balance = -3)
# e150 <- equilibrate(a150, loga.balance = -3)
# diagram(e25, ylim = c(-6, 0), dy = 0.15)
# diagram(e150, add = TRUE, names = FALSE, col = "red")
# diagram(e25, alpha = TRUE, dy = -0.25)
# diagram(e150, alpha = TRUE, add = TRUE, names = FALSE, col = "red")
## ----bjerrum_diagram, echo=12:13, eval=FALSE----------------------------------
# par(mfrow = c(3, 1))
# basis("CHNOS+")
# species(c("CO2", "HCO3-", "CO3-2"))
# a25 <- affinity(pH = c(4, 13))
# a150 <- affinity(pH = c(4, 13), T = 150)
# diagram(a25, dy = 0.4)
# diagram(a150, add = TRUE, names = FALSE, col = "red")
# e25 <- equilibrate(a25, loga.balance = -3)
# e150 <- equilibrate(a150, loga.balance = -3)
# diagram(e25, ylim = c(-6, 0), dy = 0.15)
# diagram(e150, add = TRUE, names = FALSE, col = "red")
# diagram(e25, alpha = TRUE, dy = -0.25)
# diagram(e150, alpha = TRUE, add = TRUE, names = FALSE, col = "red")
## ----corundum, fig.margin=TRUE, fig.width=4, fig.height=4, out.width="100%", results="hide", message=FALSE, cache=TRUE, fig.cap="Solubility of corundum (green line) and equilibrium concentrations of aqueous species (black lines).", pngquant=pngquant, timeit=timeit----
add.OBIGT("SLOP98")
basis(c("Al+3", "H2O", "H+", "O2"))
species("corundum")
iaq <- c("Al+3", "AlO2-", "AlOH+2", "AlO+", "HAlO2")
s <- solubility(iaq, pH = c(0, 10), IS = 0, in.terms.of = "Al+3")
diagram(s, type = "loga.balance", ylim = c(-10, 0), lwd = 4, col = "green3")
diagram(s, add = TRUE, adj = c(0, 1, 2.1, -0.2, -1.5), dy = c(0, 0, 4, -0.3, 0.1))
legend("topright", c("25 °C", "1 bar"), text.font = 2, bty = "n")
reset()
## ----Alberty------------------------------------------------------------------
oldnon <- nonideal("Alberty")
## ----subcrt_IS----------------------------------------------------------------
subcrt(c("MgATP-2", "MgHATP-", "MgH2ATP"),
T = c(25, 100), IS = c(0, 0.25), property = "G")$out
## ----info_ATP, results="hide"-------------------------------------------------
info(" ATP")
## ----T_100--------------------------------------------------------------------
T <- 100
## ----ATP, eval=FALSE, echo=2:6------------------------------------------------
# par(mfrow = c(1, 4), mar = c(3.1, 3.6, 2.1, 1.6), mgp = c(1.8, 0.5, 0))
# basis("MgCHNOPS+")
# species(c("ATP-4", "HATP-3", "H2ATP-2", "H3ATP-", "H4ATP"))
# a <- affinity(pH = c(3, 9), T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, tplot = FALSE)
# title(main = describe.property("T", T))
# alphas <- do.call(rbind, d$plotvals)
# nH <- alphas * 0:4
# Hlab <- substitute(italic(N)[H^`+`])
# plot(a$vals[[1]], colSums(nH), type = "l", xlab = "pH", ylab=Hlab, lty=2, col=2)
# a <- affinity(pH = c(3, 9), IS = 0.25, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# nH <- alphas * 0:4
# lines(a$vals[[1]], colSums(nH))
# legend("topright", legend = c("I = 0 M", "I = 0.25 M"), lty = 2:1, col = 2:1, cex = 0.8)
# ATP.H <- substitute("ATP and H"^`+`)
# title(main = ATP.H)
# species(c("MgATP-2", "MgHATP-", "MgH2ATP", "Mg2ATP"), add = TRUE)
# Hplot <- function(pMg, IS = 0.25) {
# basis("Mg+2", -pMg)
# a <- affinity(pH = c(3, 9), IS = IS, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# NH <- alphas * c(0:4, 0, 1, 2, 0)
# lines(a$vals[[1]], colSums(NH), lty = 7 - pMg, col = 7 - pMg)
# }
# plot(c(3, 9), c(0, 2), type = "n", xlab = "pH", ylab = Hlab)
# lapply(2:6, Hplot)
# legend("topright", legend = paste("pMg = ", 2:6), lty = 5:1, col = 5:1, cex = 0.8)
# ATP.H.Mg <- substitute("ATP and H"^`+`~"and Mg"^`+2`)
# title(main = ATP.H.Mg)
# Mgplot <- function(pH, IS = 0.25) {
# basis("pH", pH)
# a <- affinity(`Mg+2` = c(-2, -7), IS = IS, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# NMg <- alphas * species()$`Mg+`
# lines(-a$vals[[1]], colSums(NMg), lty = 10 - pH, col = 10 - pH)
# }
# Mglab <- substitute(italic(N)[Mg^`+2`])
# plot(c(2, 7), c(0, 1.2), type = "n", xlab = "pMg", ylab = Mglab)
# lapply(3:9, Mgplot)
# legend("topright", legend = paste("pH = ", 3:9), lty = 7:1, col = 7:1, cex = 0.8)
# title(main = ATP.H.Mg)
## ----ATP, eval=FALSE, echo=8:11-----------------------------------------------
# par(mfrow = c(1, 4), mar = c(3.1, 3.6, 2.1, 1.6), mgp = c(1.8, 0.5, 0))
# basis("MgCHNOPS+")
# species(c("ATP-4", "HATP-3", "H2ATP-2", "H3ATP-", "H4ATP"))
# a <- affinity(pH = c(3, 9), T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, tplot = FALSE)
# title(main = describe.property("T", T))
# alphas <- do.call(rbind, d$plotvals)
# nH <- alphas * 0:4
# Hlab <- substitute(italic(N)[H^`+`])
# plot(a$vals[[1]], colSums(nH), type = "l", xlab = "pH", ylab=Hlab, lty=2, col=2)
# a <- affinity(pH = c(3, 9), IS = 0.25, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# nH <- alphas * 0:4
# lines(a$vals[[1]], colSums(nH))
# legend("topright", legend = c("I = 0 M", "I = 0.25 M"), lty = 2:1, col = 2:1, cex = 0.8)
# ATP.H <- substitute("ATP and H"^`+`)
# title(main = ATP.H)
# species(c("MgATP-2", "MgHATP-", "MgH2ATP", "Mg2ATP"), add = TRUE)
# Hplot <- function(pMg, IS = 0.25) {
# basis("Mg+2", -pMg)
# a <- affinity(pH = c(3, 9), IS = IS, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# NH <- alphas * c(0:4, 0, 1, 2, 0)
# lines(a$vals[[1]], colSums(NH), lty = 7 - pMg, col = 7 - pMg)
# }
# plot(c(3, 9), c(0, 2), type = "n", xlab = "pH", ylab = Hlab)
# lapply(2:6, Hplot)
# legend("topright", legend = paste("pMg = ", 2:6), lty = 5:1, col = 5:1, cex = 0.8)
# ATP.H.Mg <- substitute("ATP and H"^`+`~"and Mg"^`+2`)
# title(main = ATP.H.Mg)
# Mgplot <- function(pH, IS = 0.25) {
# basis("pH", pH)
# a <- affinity(`Mg+2` = c(-2, -7), IS = IS, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# NMg <- alphas * species()$`Mg+`
# lines(-a$vals[[1]], colSums(NMg), lty = 10 - pH, col = 10 - pH)
# }
# Mglab <- substitute(italic(N)[Mg^`+2`])
# plot(c(2, 7), c(0, 1.2), type = "n", xlab = "pMg", ylab = Mglab)
# lapply(3:9, Mgplot)
# legend("topright", legend = paste("pH = ", 3:9), lty = 7:1, col = 7:1, cex = 0.8)
# title(main = ATP.H.Mg)
## ----ATP, eval=FALSE, echo=12:17----------------------------------------------
# par(mfrow = c(1, 4), mar = c(3.1, 3.6, 2.1, 1.6), mgp = c(1.8, 0.5, 0))
# basis("MgCHNOPS+")
# species(c("ATP-4", "HATP-3", "H2ATP-2", "H3ATP-", "H4ATP"))
# a <- affinity(pH = c(3, 9), T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, tplot = FALSE)
# title(main = describe.property("T", T))
# alphas <- do.call(rbind, d$plotvals)
# nH <- alphas * 0:4
# Hlab <- substitute(italic(N)[H^`+`])
# plot(a$vals[[1]], colSums(nH), type = "l", xlab = "pH", ylab=Hlab, lty=2, col=2)
# a <- affinity(pH = c(3, 9), IS = 0.25, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# nH <- alphas * 0:4
# lines(a$vals[[1]], colSums(nH))
# legend("topright", legend = c("I = 0 M", "I = 0.25 M"), lty = 2:1, col = 2:1, cex = 0.8)
# ATP.H <- substitute("ATP and H"^`+`)
# title(main = ATP.H)
# species(c("MgATP-2", "MgHATP-", "MgH2ATP", "Mg2ATP"), add = TRUE)
# Hplot <- function(pMg, IS = 0.25) {
# basis("Mg+2", -pMg)
# a <- affinity(pH = c(3, 9), IS = IS, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# NH <- alphas * c(0:4, 0, 1, 2, 0)
# lines(a$vals[[1]], colSums(NH), lty = 7 - pMg, col = 7 - pMg)
# }
# plot(c(3, 9), c(0, 2), type = "n", xlab = "pH", ylab = Hlab)
# lapply(2:6, Hplot)
# legend("topright", legend = paste("pMg = ", 2:6), lty = 5:1, col = 5:1, cex = 0.8)
# ATP.H.Mg <- substitute("ATP and H"^`+`~"and Mg"^`+2`)
# title(main = ATP.H.Mg)
# Mgplot <- function(pH, IS = 0.25) {
# basis("pH", pH)
# a <- affinity(`Mg+2` = c(-2, -7), IS = IS, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# NMg <- alphas * species()$`Mg+`
# lines(-a$vals[[1]], colSums(NMg), lty = 10 - pH, col = 10 - pH)
# }
# Mglab <- substitute(italic(N)[Mg^`+2`])
# plot(c(2, 7), c(0, 1.2), type = "n", xlab = "pMg", ylab = Mglab)
# lapply(3:9, Mgplot)
# legend("topright", legend = paste("pH = ", 3:9), lty = 7:1, col = 7:1, cex = 0.8)
# title(main = ATP.H.Mg)
## ----ATP, eval=FALSE, echo=21-------------------------------------------------
# par(mfrow = c(1, 4), mar = c(3.1, 3.6, 2.1, 1.6), mgp = c(1.8, 0.5, 0))
# basis("MgCHNOPS+")
# species(c("ATP-4", "HATP-3", "H2ATP-2", "H3ATP-", "H4ATP"))
# a <- affinity(pH = c(3, 9), T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, tplot = FALSE)
# title(main = describe.property("T", T))
# alphas <- do.call(rbind, d$plotvals)
# nH <- alphas * 0:4
# Hlab <- substitute(italic(N)[H^`+`])
# plot(a$vals[[1]], colSums(nH), type = "l", xlab = "pH", ylab=Hlab, lty=2, col=2)
# a <- affinity(pH = c(3, 9), IS = 0.25, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# nH <- alphas * 0:4
# lines(a$vals[[1]], colSums(nH))
# legend("topright", legend = c("I = 0 M", "I = 0.25 M"), lty = 2:1, col = 2:1, cex = 0.8)
# ATP.H <- substitute("ATP and H"^`+`)
# title(main = ATP.H)
# species(c("MgATP-2", "MgHATP-", "MgH2ATP", "Mg2ATP"), add = TRUE)
# Hplot <- function(pMg, IS = 0.25) {
# basis("Mg+2", -pMg)
# a <- affinity(pH = c(3, 9), IS = IS, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# NH <- alphas * c(0:4, 0, 1, 2, 0)
# lines(a$vals[[1]], colSums(NH), lty = 7 - pMg, col = 7 - pMg)
# }
# plot(c(3, 9), c(0, 2), type = "n", xlab = "pH", ylab = Hlab)
# lapply(2:6, Hplot)
# legend("topright", legend = paste("pMg = ", 2:6), lty = 5:1, col = 5:1, cex = 0.8)
# ATP.H.Mg <- substitute("ATP and H"^`+`~"and Mg"^`+2`)
# title(main = ATP.H.Mg)
# Mgplot <- function(pH, IS = 0.25) {
# basis("pH", pH)
# a <- affinity(`Mg+2` = c(-2, -7), IS = IS, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# NMg <- alphas * species()$`Mg+`
# lines(-a$vals[[1]], colSums(NMg), lty = 10 - pH, col = 10 - pH)
# }
# Mglab <- substitute(italic(N)[Mg^`+2`])
# plot(c(2, 7), c(0, 1.2), type = "n", xlab = "pMg", ylab = Mglab)
# lapply(3:9, Mgplot)
# legend("topright", legend = paste("pH = ", 3:9), lty = 7:1, col = 7:1, cex = 0.8)
# title(main = ATP.H.Mg)
## ----ATP, eval=FALSE, echo=22:30----------------------------------------------
# par(mfrow = c(1, 4), mar = c(3.1, 3.6, 2.1, 1.6), mgp = c(1.8, 0.5, 0))
# basis("MgCHNOPS+")
# species(c("ATP-4", "HATP-3", "H2ATP-2", "H3ATP-", "H4ATP"))
# a <- affinity(pH = c(3, 9), T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, tplot = FALSE)
# title(main = describe.property("T", T))
# alphas <- do.call(rbind, d$plotvals)
# nH <- alphas * 0:4
# Hlab <- substitute(italic(N)[H^`+`])
# plot(a$vals[[1]], colSums(nH), type = "l", xlab = "pH", ylab=Hlab, lty=2, col=2)
# a <- affinity(pH = c(3, 9), IS = 0.25, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# nH <- alphas * 0:4
# lines(a$vals[[1]], colSums(nH))
# legend("topright", legend = c("I = 0 M", "I = 0.25 M"), lty = 2:1, col = 2:1, cex = 0.8)
# ATP.H <- substitute("ATP and H"^`+`)
# title(main = ATP.H)
# species(c("MgATP-2", "MgHATP-", "MgH2ATP", "Mg2ATP"), add = TRUE)
# Hplot <- function(pMg, IS = 0.25) {
# basis("Mg+2", -pMg)
# a <- affinity(pH = c(3, 9), IS = IS, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# NH <- alphas * c(0:4, 0, 1, 2, 0)
# lines(a$vals[[1]], colSums(NH), lty = 7 - pMg, col = 7 - pMg)
# }
# plot(c(3, 9), c(0, 2), type = "n", xlab = "pH", ylab = Hlab)
# lapply(2:6, Hplot)
# legend("topright", legend = paste("pMg = ", 2:6), lty = 5:1, col = 5:1, cex = 0.8)
# ATP.H.Mg <- substitute("ATP and H"^`+`~"and Mg"^`+2`)
# title(main = ATP.H.Mg)
# Mgplot <- function(pH, IS = 0.25) {
# basis("pH", pH)
# a <- affinity(`Mg+2` = c(-2, -7), IS = IS, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# NMg <- alphas * species()$`Mg+`
# lines(-a$vals[[1]], colSums(NMg), lty = 10 - pH, col = 10 - pH)
# }
# Mglab <- substitute(italic(N)[Mg^`+2`])
# plot(c(2, 7), c(0, 1.2), type = "n", xlab = "pMg", ylab = Mglab)
# lapply(3:9, Mgplot)
# legend("topright", legend = paste("pH = ", 3:9), lty = 7:1, col = 7:1, cex = 0.8)
# title(main = ATP.H.Mg)
## ----ATP, eval=FALSE, echo=31:32----------------------------------------------
# par(mfrow = c(1, 4), mar = c(3.1, 3.6, 2.1, 1.6), mgp = c(1.8, 0.5, 0))
# basis("MgCHNOPS+")
# species(c("ATP-4", "HATP-3", "H2ATP-2", "H3ATP-", "H4ATP"))
# a <- affinity(pH = c(3, 9), T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, tplot = FALSE)
# title(main = describe.property("T", T))
# alphas <- do.call(rbind, d$plotvals)
# nH <- alphas * 0:4
# Hlab <- substitute(italic(N)[H^`+`])
# plot(a$vals[[1]], colSums(nH), type = "l", xlab = "pH", ylab=Hlab, lty=2, col=2)
# a <- affinity(pH = c(3, 9), IS = 0.25, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# nH <- alphas * 0:4
# lines(a$vals[[1]], colSums(nH))
# legend("topright", legend = c("I = 0 M", "I = 0.25 M"), lty = 2:1, col = 2:1, cex = 0.8)
# ATP.H <- substitute("ATP and H"^`+`)
# title(main = ATP.H)
# species(c("MgATP-2", "MgHATP-", "MgH2ATP", "Mg2ATP"), add = TRUE)
# Hplot <- function(pMg, IS = 0.25) {
# basis("Mg+2", -pMg)
# a <- affinity(pH = c(3, 9), IS = IS, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# NH <- alphas * c(0:4, 0, 1, 2, 0)
# lines(a$vals[[1]], colSums(NH), lty = 7 - pMg, col = 7 - pMg)
# }
# plot(c(3, 9), c(0, 2), type = "n", xlab = "pH", ylab = Hlab)
# lapply(2:6, Hplot)
# legend("topright", legend = paste("pMg = ", 2:6), lty = 5:1, col = 5:1, cex = 0.8)
# ATP.H.Mg <- substitute("ATP and H"^`+`~"and Mg"^`+2`)
# title(main = ATP.H.Mg)
# Mgplot <- function(pH, IS = 0.25) {
# basis("pH", pH)
# a <- affinity(`Mg+2` = c(-2, -7), IS = IS, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# NMg <- alphas * species()$`Mg+`
# lines(-a$vals[[1]], colSums(NMg), lty = 10 - pH, col = 10 - pH)
# }
# Mglab <- substitute(italic(N)[Mg^`+2`])
# plot(c(2, 7), c(0, 1.2), type = "n", xlab = "pMg", ylab = Mglab)
# lapply(3:9, Mgplot)
# legend("topright", legend = paste("pH = ", 3:9), lty = 7:1, col = 7:1, cex = 0.8)
# title(main = ATP.H.Mg)
## ----ATP, eval=FALSE, echo=36:44----------------------------------------------
# par(mfrow = c(1, 4), mar = c(3.1, 3.6, 2.1, 1.6), mgp = c(1.8, 0.5, 0))
# basis("MgCHNOPS+")
# species(c("ATP-4", "HATP-3", "H2ATP-2", "H3ATP-", "H4ATP"))
# a <- affinity(pH = c(3, 9), T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, tplot = FALSE)
# title(main = describe.property("T", T))
# alphas <- do.call(rbind, d$plotvals)
# nH <- alphas * 0:4
# Hlab <- substitute(italic(N)[H^`+`])
# plot(a$vals[[1]], colSums(nH), type = "l", xlab = "pH", ylab=Hlab, lty=2, col=2)
# a <- affinity(pH = c(3, 9), IS = 0.25, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# nH <- alphas * 0:4
# lines(a$vals[[1]], colSums(nH))
# legend("topright", legend = c("I = 0 M", "I = 0.25 M"), lty = 2:1, col = 2:1, cex = 0.8)
# ATP.H <- substitute("ATP and H"^`+`)
# title(main = ATP.H)
# species(c("MgATP-2", "MgHATP-", "MgH2ATP", "Mg2ATP"), add = TRUE)
# Hplot <- function(pMg, IS = 0.25) {
# basis("Mg+2", -pMg)
# a <- affinity(pH = c(3, 9), IS = IS, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# NH <- alphas * c(0:4, 0, 1, 2, 0)
# lines(a$vals[[1]], colSums(NH), lty = 7 - pMg, col = 7 - pMg)
# }
# plot(c(3, 9), c(0, 2), type = "n", xlab = "pH", ylab = Hlab)
# lapply(2:6, Hplot)
# legend("topright", legend = paste("pMg = ", 2:6), lty = 5:1, col = 5:1, cex = 0.8)
# ATP.H.Mg <- substitute("ATP and H"^`+`~"and Mg"^`+2`)
# title(main = ATP.H.Mg)
# Mgplot <- function(pH, IS = 0.25) {
# basis("pH", pH)
# a <- affinity(`Mg+2` = c(-2, -7), IS = IS, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# NMg <- alphas * species()$`Mg+`
# lines(-a$vals[[1]], colSums(NMg), lty = 10 - pH, col = 10 - pH)
# }
# Mglab <- substitute(italic(N)[Mg^`+2`])
# plot(c(2, 7), c(0, 1.2), type = "n", xlab = "pMg", ylab = Mglab)
# lapply(3:9, Mgplot)
# legend("topright", legend = paste("pH = ", 3:9), lty = 7:1, col = 7:1, cex = 0.8)
# title(main = ATP.H.Mg)
## ----ATP, eval=FALSE, echo=45:47----------------------------------------------
# par(mfrow = c(1, 4), mar = c(3.1, 3.6, 2.1, 1.6), mgp = c(1.8, 0.5, 0))
# basis("MgCHNOPS+")
# species(c("ATP-4", "HATP-3", "H2ATP-2", "H3ATP-", "H4ATP"))
# a <- affinity(pH = c(3, 9), T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, tplot = FALSE)
# title(main = describe.property("T", T))
# alphas <- do.call(rbind, d$plotvals)
# nH <- alphas * 0:4
# Hlab <- substitute(italic(N)[H^`+`])
# plot(a$vals[[1]], colSums(nH), type = "l", xlab = "pH", ylab=Hlab, lty=2, col=2)
# a <- affinity(pH = c(3, 9), IS = 0.25, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# nH <- alphas * 0:4
# lines(a$vals[[1]], colSums(nH))
# legend("topright", legend = c("I = 0 M", "I = 0.25 M"), lty = 2:1, col = 2:1, cex = 0.8)
# ATP.H <- substitute("ATP and H"^`+`)
# title(main = ATP.H)
# species(c("MgATP-2", "MgHATP-", "MgH2ATP", "Mg2ATP"), add = TRUE)
# Hplot <- function(pMg, IS = 0.25) {
# basis("Mg+2", -pMg)
# a <- affinity(pH = c(3, 9), IS = IS, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# NH <- alphas * c(0:4, 0, 1, 2, 0)
# lines(a$vals[[1]], colSums(NH), lty = 7 - pMg, col = 7 - pMg)
# }
# plot(c(3, 9), c(0, 2), type = "n", xlab = "pH", ylab = Hlab)
# lapply(2:6, Hplot)
# legend("topright", legend = paste("pMg = ", 2:6), lty = 5:1, col = 5:1, cex = 0.8)
# ATP.H.Mg <- substitute("ATP and H"^`+`~"and Mg"^`+2`)
# title(main = ATP.H.Mg)
# Mgplot <- function(pH, IS = 0.25) {
# basis("pH", pH)
# a <- affinity(`Mg+2` = c(-2, -7), IS = IS, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# NMg <- alphas * species()$`Mg+`
# lines(-a$vals[[1]], colSums(NMg), lty = 10 - pH, col = 10 - pH)
# }
# Mglab <- substitute(italic(N)[Mg^`+2`])
# plot(c(2, 7), c(0, 1.2), type = "n", xlab = "pMg", ylab = Mglab)
# lapply(3:9, Mgplot)
# legend("topright", legend = paste("pH = ", 3:9), lty = 7:1, col = 7:1, cex = 0.8)
# title(main = ATP.H.Mg)
## ----ATP, fig.fullwidth=TRUE, fig.width=10, fig.height=2.5, dpi=hidpi, out.width="100%", echo=FALSE, message=FALSE, results="hide", fig.cap="Binding of H<sup>+</sup> and Mg<sup>+2</sup> to ATP at 100 °C and *I* = 0 M (first plot) or *I* = 0.25 M (third and fourth plots).", cache=TRUE, pngquant=pngquant, timeit=timeit----
par(mfrow = c(1, 4), mar = c(3.1, 3.6, 2.1, 1.6), mgp = c(1.8, 0.5, 0))
basis("MgCHNOPS+")
species(c("ATP-4", "HATP-3", "H2ATP-2", "H3ATP-", "H4ATP"))
a <- affinity(pH = c(3, 9), T = T)
e <- equilibrate(a)
d <- diagram(e, alpha = TRUE, tplot = FALSE)
title(main = describe.property("T", T))
alphas <- do.call(rbind, d$plotvals)
nH <- alphas * 0:4
Hlab <- substitute(italic(N)[H^`+`])
plot(a$vals[[1]], colSums(nH), type = "l", xlab = "pH", ylab=Hlab, lty=2, col=2)
a <- affinity(pH = c(3, 9), IS = 0.25, T = T)
e <- equilibrate(a)
d <- diagram(e, alpha = TRUE, plot.it = FALSE)
alphas <- do.call(rbind, d$plotvals)
nH <- alphas * 0:4
lines(a$vals[[1]], colSums(nH))
legend("topright", legend = c("I = 0 M", "I = 0.25 M"), lty = 2:1, col = 2:1, cex = 0.8)
ATP.H <- substitute("ATP and H"^`+`)
title(main = ATP.H)
species(c("MgATP-2", "MgHATP-", "MgH2ATP", "Mg2ATP"), add = TRUE)
Hplot <- function(pMg, IS = 0.25) {
basis("Mg+2", -pMg)
a <- affinity(pH = c(3, 9), IS = IS, T = T)
e <- equilibrate(a)
d <- diagram(e, alpha = TRUE, plot.it = FALSE)
alphas <- do.call(rbind, d$plotvals)
NH <- alphas * c(0:4, 0, 1, 2, 0)
lines(a$vals[[1]], colSums(NH), lty = 7 - pMg, col = 7 - pMg)
}
plot(c(3, 9), c(0, 2), type = "n", xlab = "pH", ylab = Hlab)
lapply(2:6, Hplot)
legend("topright", legend = paste("pMg = ", 2:6), lty = 5:1, col = 5:1, cex = 0.8)
ATP.H.Mg <- substitute("ATP and H"^`+`~"and Mg"^`+2`)
title(main = ATP.H.Mg)
Mgplot <- function(pH, IS = 0.25) {
basis("pH", pH)
a <- affinity(`Mg+2` = c(-2, -7), IS = IS, T = T)
e <- equilibrate(a)
d <- diagram(e, alpha = TRUE, plot.it = FALSE)
alphas <- do.call(rbind, d$plotvals)
NMg <- alphas * species()$`Mg+`
lines(-a$vals[[1]], colSums(NMg), lty = 10 - pH, col = 10 - pH)
}
Mglab <- substitute(italic(N)[Mg^`+2`])
plot(c(2, 7), c(0, 1.2), type = "n", xlab = "pMg", ylab = Mglab)
lapply(3:9, Mgplot)
legend("topright", legend = paste("pH = ", 3:9), lty = 7:1, col = 7:1, cex = 0.8)
title(main = ATP.H.Mg)
## ----oldnon-------------------------------------------------------------------
nonideal(oldnon)
## ----pinfo_LYSC_CHICK---------------------------------------------------------
p1 <- pinfo("LYSC_CHICK")
p2 <- pinfo(c("SHH", "OLIG2"), "HUMAN")
pinfo(c(p1, p2))
## ----formula_LYSC_CHICK-------------------------------------------------------
pl <- protein.length("LYSC_CHICK")
pf <- protein.formula("LYSC_CHICK")
list(length = pl, protein = pf, residue = pf / pl,
ZC_protein = ZC(pf), ZC_residue = ZC(pf / pl))
## ----subcrt_LYSC_CHICK, message=FALSE-----------------------------------------
subcrt("LYSC_CHICK")$out[[1]][1:6, ]
## ----protein_Cp, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, out.width="100%", echo=FALSE, message=FALSE, fig.cap='The heat capacity calculated by group additivity closely approximates experimental values for aqueous proteins. For a related figure showing the effects of ionization in the calculations, see <span style="color:blue">`?ionize.aa`</span>.', cache=TRUE, pngquant=pngquant, timeit=timeit----
PM90 <- read.csv(system.file("extdata/cpetc/PM90.csv", package = "CHNOSZ"))
plength <- protein.length(colnames(PM90)[2:5])
Cp_expt <- t(t(PM90[, 2:5]) / plength)
matplot(PM90[, 1], Cp_expt, type = "p", pch = 19,
xlab = axis.label("T"), ylab = axis.label("Cp0"), ylim = c(110, 280))
for(i in 1:4) {
pname <- colnames(Cp_expt)[i]
aq <- subcrt(pname, "aq", T = seq(0, 150))$out[[1]]
cr <- subcrt(pname, "cr", T = seq(0, 150))$out[[1]]
lines(aq$T, aq$Cp / plength[i], col = i)
lines(cr$T, cr$Cp / plength[i], col = i, lty = 2)
}
legend("right", legend = colnames(Cp_expt),
col = 1:4, pch = 19, lty = 1, bty = "n", cex = 0.9)
legend("bottomright", legend = c("experimental", "calculated (aq)",
"calculated (cr)"), lty = c(NA, 1, 2), pch = c(19, NA, NA), bty = "n")
## ----protein_Cp, eval=FALSE, echo=1:5-----------------------------------------
# PM90 <- read.csv(system.file("extdata/cpetc/PM90.csv", package = "CHNOSZ"))
# plength <- protein.length(colnames(PM90)[2:5])
# Cp_expt <- t(t(PM90[, 2:5]) / plength)
# matplot(PM90[, 1], Cp_expt, type = "p", pch = 19,
# xlab = axis.label("T"), ylab = axis.label("Cp0"), ylim = c(110, 280))
# for(i in 1:4) {
# pname <- colnames(Cp_expt)[i]
# aq <- subcrt(pname, "aq", T = seq(0, 150))$out[[1]]
# cr <- subcrt(pname, "cr", T = seq(0, 150))$out[[1]]
# lines(aq$T, aq$Cp / plength[i], col = i)
# lines(cr$T, cr$Cp / plength[i], col = i, lty = 2)
# }
# legend("right", legend = colnames(Cp_expt),
# col = 1:4, pch = 19, lty = 1, bty = "n", cex = 0.9)
# legend("bottomright", legend = c("experimental", "calculated (aq)",
# "calculated (cr)"), lty = c(NA, 1, 2), pch = c(19, NA, NA), bty = "n")
## ----protein_Cp, eval=FALSE, echo=-(1:5)--------------------------------------
# PM90 <- read.csv(system.file("extdata/cpetc/PM90.csv", package = "CHNOSZ"))
# plength <- protein.length(colnames(PM90)[2:5])
# Cp_expt <- t(t(PM90[, 2:5]) / plength)
# matplot(PM90[, 1], Cp_expt, type = "p", pch = 19,
# xlab = axis.label("T"), ylab = axis.label("Cp0"), ylim = c(110, 280))
# for(i in 1:4) {
# pname <- colnames(Cp_expt)[i]
# aq <- subcrt(pname, "aq", T = seq(0, 150))$out[[1]]
# cr <- subcrt(pname, "cr", T = seq(0, 150))$out[[1]]
# lines(aq$T, aq$Cp / plength[i], col = i)
# lines(cr$T, cr$Cp / plength[i], col = i, lty = 2)
# }
# legend("right", legend = colnames(Cp_expt),
# col = 1:4, pch = 19, lty = 1, bty = "n", cex = 0.9)
# legend("bottomright", legend = c("experimental", "calculated (aq)",
# "calculated (cr)"), lty = c(NA, 1, 2), pch = c(19, NA, NA), bty = "n")
## ----protein_ionization, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, out.width="100%", echo=FALSE, results="hide", message=FALSE, fig.cap='Affinity of ionization of proteins. See [<span style="color:blue">demo(ionize)</span>](../demo) for ionization properties calculated as a function of temperature and pH.', cache=TRUE, pngquant=pngquant, timeit=timeit----
ip <- pinfo(c("CYC_BOVIN", "LYSC_CHICK", "MYG_PHYCA", "RNAS1_BOVIN"))
basis("CHNOS+")
a_ion <- affinity(pH = c(0, 14), iprotein = ip)
basis("CHNOS")
a_nonion <- affinity(iprotein = ip)
plot(c(0, 14), c(50, 300), xlab = "pH", ylab = quote(italic(A/2.303*RT)), type="n")
for(i in 1:4) {
A_ion <- as.numeric(a_ion$values[[i]])
A_nonion <- as.numeric(a_nonion$values[[i]])
lines(a_ion$vals[[1]], A_ion - A_nonion, col=i)
}
legend("topright", legend = a_ion$species$name,
col = 1:4, lty = 1, bty = "n", cex = 0.9)
## ----protein_ionization, eval=FALSE-------------------------------------------
# ip <- pinfo(c("CYC_BOVIN", "LYSC_CHICK", "MYG_PHYCA", "RNAS1_BOVIN"))
# basis("CHNOS+")
# a_ion <- affinity(pH = c(0, 14), iprotein = ip)
# basis("CHNOS")
# a_nonion <- affinity(iprotein = ip)
# plot(c(0, 14), c(50, 300), xlab = "pH", ylab = quote(italic(A/2.303*RT)), type="n")
# for(i in 1:4) {
# A_ion <- as.numeric(a_ion$values[[i]])
# A_nonion <- as.numeric(a_nonion$values[[i]])
# lines(a_ion$vals[[1]], A_ion - A_nonion, col=i)
# }
# legend("topright", legend = a_ion$species$name,
# col = 1:4, lty = 1, bty = "n", cex = 0.9)
## ----basis_CHNOS, echo=FALSE, results="hide"----------------------------------
basis("CHNOS")
## ----species_protein, results="hide", message=FALSE, echo=1:2-----------------
species(c("CYC_BOVIN", "LYSC_CHICK", "MYG_PHYCA", "RNAS1_BOVIN"))
a_nonion_species <- affinity()
## ----nonion_species_values----------------------------------------------------
unlist(a_nonion_species$values)
## ----nonion_values------------------------------------------------------------
unlist(a_nonion$values)
## ----rubisco_svg, echo=FALSE, results="hide", fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, out.width="100%", fig.ext='svg', custom.plot=TRUE, embed.tag=TRUE, fig.cap='Average oxidation state of carbon in Rubisco compared with optimal growth temperature of organisms. **This is an interactive image.** Move the mouse over the points to show the names of the organisms, and click to open a reference in a new window. (Made with [**RSVGTipsDevice**](https://cran.r-project.org/package=RSVGTipsDevice) using code that can be found in the source of this document.)'----
## copies the premade SVG image to the knitr figure path
file.copy("rubisco.svg", fig_path(".svg"))
## the code for making the SVG image -- not "live" in the vignette because RSVGTipsDevice isn't available on Windows
#if(require(RSVGTipsDevice)) {
# datfile <- system.file("extdata/cpetc/rubisco.csv", package = "CHNOSZ")
# fastafile <- system.file("extdata/protein/rubisco.fasta", package = "CHNOSZ")
# dat <- read.csv(datfile)
# aa <- read.fasta(fastafile)
# Topt <- (dat$T1 + dat$T2) / 2
# idat <- match(dat$ID, substr(aa$protein, 4, 9))
# aa <- aa[idat, ]
# ZC <- ZC(protein.formula(aa))
# pch <- match(dat$domain, c("E", "B", "A")) - 1
# col <- match(dat$domain, c("A", "B", "E")) + 1
# # because the tooltip titles in the SVG file are shown by recent browsers,
# # we do not need to draw the tooltips explicitly, so set toolTipMode=0
# devSVGTips("rubisco.svg", toolTipMode=0, title="Rubisco")
# par(cex=1.4)
# # unfortunately, plotmath can't be used with devSVGTips,
# # so axis labels here don't contain italics.
# plot(Topt, ZC, type="n", xlab="T, °C", ylab="ZC")
# n <- rep(1:9, 3)
# for(i in seq_along(Topt)) {
# # adjust cex to make the symbols look the same size
# cex <- ifelse(pch[i]==1, 2.5, 3.5)
# points(Topt[i], ZC[i], pch=pch[i], cex=cex, col=col[i])
# URL <- dat$URL[i]
# setSVGShapeURL(URL, target="_blank")
# setSVGShapeContents(paste0("<title>", dat$species[i], "</title>"))
# text(Topt[i], ZC[i], n[i], cex = 1.2)
# }
# abline(v = c(36, 63), lty = 2, col = "grey")
# legend("topright", legend = c("Archaea", "Bacteria", "Eukaryota"),
# pch = c(2, 1, 0), col = 2:4, cex=1.5, pt.cex = c(3, 2.3, 3), bty="n")
# dev.off()
#}
## ----rubisco_ZC, fig.keep="none", message=FALSE-------------------------------
datfile <- system.file("extdata/cpetc/rubisco.csv", package = "CHNOSZ")
fastafile <- system.file("extdata/protein/rubisco.fasta", package = "CHNOSZ")
dat <- read.csv(datfile)
aa <- read.fasta(fastafile)
Topt <- (dat$T1 + dat$T2) / 2
idat <- match(dat$ID, substr(aa$protein, 4, 9))
aa <- aa[idat, ]
ZC <- ZC(protein.formula(aa))
pch <- match(dat$domain, c("E", "B", "A")) - 1
col <- match(dat$domain, c("A", "B", "E")) + 1
plot(Topt, ZC, pch = pch, cex = 2, col = col,
xlab = expression(list(italic(T)[opt], degree*C)),
ylab = expression(italic(Z)[C]))
text(Topt, ZC, rep(1:9, 3), cex = 0.8)
abline(v = c(36, 63), lty = 2, col = "grey")
legend("topright", legend = c("Archaea", "Bacteria", "Eukaryota"),
pch = c(2, 1, 0), col = 2:4, pt.cex = 2)
## ----rubisco_O2, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, out.width="100%", echo=FALSE, results="hide", message=FALSE, fig.cap="Elemental compositions of proteins projected into different sets of basis species.", cache=TRUE, pngquant=pngquant, timeit=timeit----
layout(matrix(1:4, nrow = 2))
par(mgp = c(1.8, 0.5, 0))
pl <- protein.length(aa)
ZClab <- axis.label("ZC")
nO2lab <- expression(bar(italic(n))[O[2]])
nH2Olab <- expression(bar(italic(n))[H[2]*O])
lapply(c("CHNOS", "QEC"), function(thisbasis) {
basis(thisbasis)
pb <- protein.basis(aa)
nO2 <- pb[, "O2"] / pl
plot(ZC, nO2, pch = pch, col = col, xlab = ZClab, ylab = nO2lab)
nH2O <- pb[, "H2O"] / pl
plot(ZC, nH2O, pch = pch, col = col, xlab = ZClab, ylab = nH2Olab)
mtext(thisbasis, font = 2)
})
## ----rubisco_O2, eval=FALSE---------------------------------------------------
# layout(matrix(1:4, nrow = 2))
# par(mgp = c(1.8, 0.5, 0))
# pl <- protein.length(aa)
# ZClab <- axis.label("ZC")
# nO2lab <- expression(bar(italic(n))[O[2]])
# nH2Olab <- expression(bar(italic(n))[H[2]*O])
# lapply(c("CHNOS", "QEC"), function(thisbasis) {
# basis(thisbasis)
# pb <- protein.basis(aa)
# nO2 <- pb[, "O2"] / pl
# plot(ZC, nO2, pch = pch, col = col, xlab = ZClab, ylab = nO2lab)
# nH2O <- pb[, "H2O"] / pl
# plot(ZC, nH2O, pch = pch, col = col, xlab = ZClab, ylab = nH2Olab)
# mtext(thisbasis, font = 2)
# })
## ----yeastgfp-----------------------------------------------------------------
protein <- c("YDL195W", "YHR098C", "YIL109C", "YLR208W", "YNL049C", "YPL085W")
abundance <- c(1840, 12200, NA, 21400, 1720, 358)
ina <- is.na(abundance)
## ----add_protein_yeast, message=FALSE-----------------------------------------
ip <- match(protein[!ina], thermo()$protein$protein)
## ----unitize------------------------------------------------------------------
pl <- protein.length(ip)
logact <- unitize(numeric(5), pl)
logabundance <- unitize(log10(abundance[!ina]), pl)
## ----yeastplot, eval=FALSE, echo=1:5------------------------------------------
# par(mfrow = c(1, 3))
# basis("CHNOS+")
# a <- affinity(O2 = c(-80, -73), iprotein = ip, loga.protein = logact)
# e <- equilibrate(a)
# diagram(e, ylim = c(-5, -2), col = 1:5, lwd = 2)
# e <- equilibrate(a, normalize = TRUE)
# diagram(e, ylim = c(-5, -2.5), col = 1:5, lwd = 2)
# abline(h = logabundance, lty = 1:5, col = 1:5)
## ----yeastplot, eval=FALSE, echo=6:8------------------------------------------
# par(mfrow = c(1, 3))
# basis("CHNOS+")
# a <- affinity(O2 = c(-80, -73), iprotein = ip, loga.protein = logact)
# e <- equilibrate(a)
# diagram(e, ylim = c(-5, -2), col = 1:5, lwd = 2)
# e <- equilibrate(a, normalize = TRUE)
# diagram(e, ylim = c(-5, -2.5), col = 1:5, lwd = 2)
# abline(h = logabundance, lty = 1:5, col = 1:5)
## ----yeastplot, fig.width=6, fig.height=2.5, dpi=hidpi, out.width="100%", echo=FALSE, message=FALSE, results="hide", cache=TRUE, fig.cap="ER-to-Golgi proteins: calculations without and with length normalization.", pngquant=pngquant, timeit=timeit----
par(mfrow = c(1, 3))
basis("CHNOS+")
a <- affinity(O2 = c(-80, -73), iprotein = ip, loga.protein = logact)
e <- equilibrate(a)
diagram(e, ylim = c(-5, -2), col = 1:5, lwd = 2)
e <- equilibrate(a, normalize = TRUE)
diagram(e, ylim = c(-5, -2.5), col = 1:5, lwd = 2)
abline(h = logabundance, lty = 1:5, col = 1:5)
## ----read_csv-----------------------------------------------------------------
file <- system.file("extdata/protein/POLG.csv", package = "CHNOSZ")
aa_POLG <- read.csv(file, as.is = TRUE, nrows = 5)
## ----read_fasta, message=FALSE------------------------------------------------
file <- system.file("extdata/protein/EF-Tu.aln", package = "CHNOSZ")
aa_Ef <- read.fasta(file, iseq = 1:2)
## ----seq2aa-------------------------------------------------------------------
aa_PRIO <- seq2aa("
MANLGCWMLVLFVATWSDLGLCKKRPKPGGWNTGGSRYPGQGSPGGNRYPPQGGGGWGQP
HGGGWGQPHGGGWGQPHGGGWGQPHGGGWGQGGGTHSQWNKPSKPKTNMKHMAGAAAAGA
VVGGLGGYMLGSAMSRPIIHFGSDYEDRYYRENMHRYPNQVYYRPMDEYSNQNNFVHDCV
NITIKQHTVTTTTKGENFTETDVKMMERVVEQMCITQYERESQAYYQRGSSMVLFSSPPV
ILLISFLIFLIVG
", "PRIO_HUMAN")
## ----protein_length-----------------------------------------------------------
myaa <- rbind(aa_Ef, aa_PRIO)
protein.length(myaa)
## ----thermo_refs_table, eval=FALSE--------------------------------------------
# thermo.refs() ## shows table in a browser
## ----width180, include=FALSE------------------------------------------------------------------------------------------------------------------------------------------------------
options(width = 180)
## ----thermo_refs_numeric----------------------------------------------------------------------------------------------------------------------------------------------------------
iATP <- info("ATP-4")
iMgATP <- info("MgATP-2")
thermo.refs(c(iATP, iMgATP))
## ----thermo_refs_character--------------------------------------------------------------------------------------------------------------------------------------------------------
thermo.refs(c("HDNB78", "MGN03"))
## ----thermo_refs_subcrt, message=FALSE--------------------------------------------------------------------------------------------------------------------------------------------
sres <- subcrt(c("C2H5OH", "O2", "CO2", "H2O"), c(-1, -3, 2, 3))
thermo.refs(sres)
## ----width80, include=FALSE---------------------------------------------------
options(width = 80)
## ----thermo_refs_browse, eval=FALSE-------------------------------------------
# iFo <- info("forsterite")
# ref <- thermo.refs(iFo)
# browseURL(ref$URL) ## opens a link to worldcat.org
## ----info_S3, results="hide"--------------------------------------------------
info(info("S3-"))
## ----check_OBIGT--------------------------------------------------------------
file <- system.file("extdata/adds/OBIGT_check.csv", package = "CHNOSZ")
dat <- read.csv(file, as.is = TRUE)
nrow(dat)
## ----citation_CHNOSZ, results="asis", echo = FALSE----------------------------
cref <- citation("CHNOSZ")
print(cref[1], style = "html")
## ----citation_multimetal, results="asis", echo = FALSE------------------------
print(cref[2], style = "html")
## ----maintainer_CHNOSZ--------------------------------------------------------
maintainer("CHNOSZ")
## ----file_edit_anintro, eval=FALSE--------------------------------------------
# file.edit(system.file("doc/anintro.Rmd", package = "CHNOSZ"))
## ----the_end------------------------------------------------------------------
###### ## ## ## ## ###### ##### #####
## ##===## ## \\## ## ## \\ //
###### ## ## ## ## ###### ##### #####
## ----sessionInfo--------------------------------------------------------------
sessionInfo()
Try the CHNOSZ package in your browser
Any scripts or data that you put into this service are public.
CHNOSZ documentation built on May 29, 2024, 3:30 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
## ----options, include=FALSE---------------------------------------------------
options(width = 80)
options(digits = 6)
## ----HTML, include=FALSE------------------------------------------------------
## Some frequently used HTML expressions
logfO2 <- "log<i>f</i><sub>O<sub>2</sub></sub>"
# Use lowercase here because these tend to be variable names in the examples
zc <- "<i>Z</i><sub>C</sub>"
o2 <- "O<sub>2</sub>"
h2o <- "H<sub>2</sub>O"
sio2 <- "SiO<sub>2</sub>"
ch4 <- "CH<sub>4</sub>"
## ----setup, include=FALSE-----------------------------------------------------
library(knitr)
## From "Tufte Handout" example dated 2016-12-27
# Invalidate cache when the tufte version changes
opts_chunk$set(tidy = FALSE, cache.extra = packageVersion('tufte'))
options(htmltools.dir.version = FALSE)
## Adjust plot margins
## First one from https://yihui.name/knitr/hooks/
knit_hooks$set(small.mar = function(before, options, envir) {
if (before) par(mar = c(4.2, 4.2, .1, .1)) # smaller margin on top and right
})
knit_hooks$set(tiny.mar = function(before, options, envir) {
if (before) par(mar = c(.1, .1, .1, .1)) # tiny margin all around
})
knit_hooks$set(smallish.mar = function(before, options, envir) {
if (before) par(mar = c(4.2, 4.2, 0.9, 0.9)) # smallish margins on top and right
})
## Use pngquant to optimize PNG images
knit_hooks$set(pngquant = hook_pngquant)
pngquant <- "--speed=1 --quality=0-25"
# pngquant isn't available on R-Forge ...
if (!nzchar(Sys.which("pngquant"))) pngquant <- NULL
# Set dpi 20231129
knitr::opts_chunk$set(
dpi = if(nzchar(Sys.getenv("CHNOSZ_BUILD_LARGE_VIGNETTES"))) 72 else 50
)
hidpi = if(nzchar(Sys.getenv("CHNOSZ_BUILD_LARGE_VIGNETTES"))) 100 else 50
## http://stackoverflow.com/questions/23852753/knitr-with-gridsvg
## Set up a chunk hook for manually saved plots.
knit_hooks$set(custom.plot = hook_plot_custom)
## Hook to change <img /> to <embed /> -- required for interactive SVG
hook_plot <- knit_hooks$get("plot")
knit_hooks$set(plot = function(x, options) {
x <- hook_plot(x, options)
if (!is.null(options$embed.tag) && options$embed.tag) x <- gsub("<img ", "<embed ", x)
x
})
## http://stackoverflow.com/questions/30530008/hook-to-time-knitr-chunks
now = Sys.time()
knit_hooks$set(timeit = function(before) {
if (before) { now <<- Sys.time() }
else {
paste("%", sprintf("Chunk rendering time: %s seconds.\n", round(Sys.time() - now, digits = 3)))
}
})
timeit <- NULL
## Colorize messages 20171031
## Adapted from https://gist.github.com/yihui/2629886#file-knitr-color-msg-rnw
color_block = function(color) {
function(x, options) sprintf('<pre style="color:%s">%s</pre>', color, x)
}
knit_hooks$set(warning = color_block('magenta'), error = color_block('red'), message = color_block('blue'))
## ----install_CHNOSZ, eval=FALSE-----------------------------------------------
# install.packages("CHNOSZ")
## ----library_CHNOSZ-----------------------------------------------------------
library(CHNOSZ)
## ----reset--------------------------------------------------------------------
reset()
## ----pseudocode, eval=FALSE---------------------------------------------------
# basis(...)
# species(...)
# a <- affinity(...)
# e <- equilibrate(a) ## optional
# diagram(e) ## or diagram(a)
# reset() ## clear settings for next calculation
## ----info_CH4-----------------------------------------------------------------
info("CH4")
## ----info_CH4_gas-------------------------------------------------------------
info("CH4", "gas")
## ----info_names_gas-----------------------------------------------------------
info("methane")
info("oxygen")
info("carbon dioxide")
## ----info_S_S2----------------------------------------------------------------
info("S")
info("S2")
## ----iCH4, message=FALSE------------------------------------------------------
iCH4 <- info("CH4")
info(iCH4)
## ----info_info_water----------------------------------------------------------
info(info("water"))
## ----width180, include=FALSE------------------------------------------------------------------------------------------------------------------------------------------------------
options(width = 180)
## ----info_acid--------------------------------------------------------------------------------------------------------------------------------------------------------------------
info("acid")
## ----width80, include=FALSE---------------------------------------------------
options(width = 80)
## ----info_ribose--------------------------------------------------------------
info(" ribose")
## ----info_CH4_formula, message=FALSE------------------------------------------
info(iCH4)$formula
## ----makeup_iCH4--------------------------------------------------------------
makeup(iCH4)
as.chemical.formula(makeup(iCH4))
## ----ZC_iCH4, message=FALSE---------------------------------------------------
ZC(iCH4)
ZC(info(iCH4)$formula)
ZC(makeup(iCH4))
## ----subcrt_water-------------------------------------------------------------
subcrt("water")
## ----subcrt_water_grid--------------------------------------------------------
subcrt("water", T = c(400, 500, 600), P = c(200, 400, 600), grid = "P")$out$water
## ----subcrt_water_plot, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, out.width="100%", echo=FALSE, message=FALSE, fig.cap="Isothermal contours of density (g cm<sup>-3</sup>) and pressure (bar) of water.", cache=TRUE, pngquant=pngquant, timeit=timeit----
sres <- subcrt("water", T=seq(0,1000,100), P=c(NA, seq(1,500,1)), grid="T")
water <- sres$out$water
plot(water$P, water$rho, type = "l")
## ----subcrt_water_plot, eval=FALSE--------------------------------------------
# sres <- subcrt("water", T=seq(0,1000,100), P=c(NA, seq(1,500,1)), grid="T")
# water <- sres$out$water
# plot(water$P, water$rho, type = "l")
## ----units_CH4----------------------------------------------------------------
T.units("K")
P.units("MPa")
subcrt("CH4", T = 298.15, P = 0.1)$out$CH4$G
E.units("cal")
subcrt("CH4", T = 298.15, P = 0.1)$out$CH4$G
## ----convert_G, message=FALSE-------------------------------------------------
(CH4dat <- info(info("CH4")))
convert(CH4dat$G, "J")
## ----reset--------------------------------------------------------------------
reset()
## ----subcrt_CO2---------------------------------------------------------------
subcrt(c("CO2", "CO2"), c("gas", "aq"), c(-1, 1), T = seq(0, 250, 50))
## ----CO2_logK, echo=FALSE, message=FALSE--------------------------------------
T <- seq(0, 350, 10)
CO2 <- subcrt(c("CO2", "CO2"), c("gas", "aq"), c(-1, 1), T = T)$out$logK
CO <- subcrt(c("CO", "CO"), c("gas", "aq"), c(-1, 1), T = T)$out$logK
CH4 <- subcrt(c("CH4", "CH4"), c("gas", "aq"), c(-1, 1), T = T)$out$logK
logK <- data.frame(T, CO2, CO, CH4)
## ----CO2_plot, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, out.width="100%", echo=FALSE, message=FALSE, fig.cap="Calculated equilibrium constants for dissolution of CO<sub>2</sub>, CO, and CH<sub>4</sub>.", cache=TRUE, pngquant=pngquant, timeit=timeit----
matplot(logK[, 1], logK[, -1], type = "l", col = 1, lty = 1,
xlab = axis.label("T"), ylab = axis.label("logK"))
text(80, -1.7, expr.species("CO2"))
text(240, -2.37, expr.species("CO"))
text(300, -2.57, expr.species("CH4"))
## ----CO2_logK, eval=FALSE-----------------------------------------------------
# T <- seq(0, 350, 10)
# CO2 <- subcrt(c("CO2", "CO2"), c("gas", "aq"), c(-1, 1), T = T)$out$logK
# CO <- subcrt(c("CO", "CO"), c("gas", "aq"), c(-1, 1), T = T)$out$logK
# CH4 <- subcrt(c("CH4", "CH4"), c("gas", "aq"), c(-1, 1), T = T)$out$logK
# logK <- data.frame(T, CO2, CO, CH4)
## ----CO2_plot, eval=FALSE-----------------------------------------------------
# matplot(logK[, 1], logK[, -1], type = "l", col = 1, lty = 1,
# xlab = axis.label("T"), ylab = axis.label("logK"))
# text(80, -1.7, expr.species("CO2"))
# text(240, -2.37, expr.species("CO"))
# text(300, -2.57, expr.species("CH4"))
## ----subcrt_unbalanced, results="hide"----------------------------------------
subcrt(c("CO2", "CH4"), c(-1, 1))
## ----basis_singular, error=TRUE-----------------------------------------------
basis(c("CO2", "H2", "H2CO2"))
## ----basis_CHO----------------------------------------------------------------
basis(c("CO2", "H2", "H2O"))
## ----basis_CHOZ---------------------------------------------------------------
basis(c("CO2", "H2", "H2O", "H+"))
## ----subcrt_acetoclastic, message=FALSE---------------------------------------
subcrt(c("acetate", "CH4"), c(-1, 1))$reaction
## ----subcrt_methanogenesis, message=FALSE-------------------------------------
acetate_oxidation <- subcrt("acetate", -1)
hydrogenotrophic <- subcrt("CH4", 1)
acetoclastic <- subcrt(c("acetate", "CH4"), c(-1, 1))
## ----describe_reaction_plot, fig.margin=TRUE, fig.width=3.5, fig.height=1.8, tiny.mar=TRUE, out.width="100%", pngquant=pngquant, timeit=timeit----
plot(0, 0, type = "n", axes = FALSE, ann=FALSE, xlim=c(0, 5), ylim=c(5.2, -0.2))
text(0, 0, "acetoclastic methanogenesis", adj = 0)
text(5, 1, describe.reaction(acetoclastic$reaction), adj = 1)
text(0, 2, "acetate oxidation", adj = 0)
text(5, 3, describe.reaction(acetate_oxidation$reaction), adj = 1)
text(0, 4, "hydrogenotrophic methanogenesis", adj = 0)
text(5, 5, describe.reaction(hydrogenotrophic$reaction), adj = 1)
## ----basis_mayumi, message=FALSE, results="hide"------------------------------
basis(c("CO2", "H2", "H2O", "H+"))
basis(c("CO2", "H2"), "gas")
basis(c("H2", "pH"), c(-3.92, 7.3))
## ----affinity_acetoclastic, message=FALSE-------------------------------------
subcrt(c("acetate", "CH4"), c(-1, 1),
c("aq", "gas"), logact = c(-3.4, -0.18), T = 55, P = 50)$out
## ----affinity_hydrogenotrophic, message=FALSE---------------------------------
subcrt("CH4", 1, "gas", logact = -0.18, T = 55, P = 50)$out
## ----rxnfun, message=FALSE----------------------------------------------------
rxnfun <- function(coeffs) {
subcrt(c("acetate", "CH4"), coeffs,
c("aq", "gas"), logact = c(-3.4, -0.18), T = 55, P = 50)$out
}
## ----methanogenesis_plot, fig.margin=TRUE, fig.width=4.1, fig.height=4.1, small.mar=TRUE, out.width="100%", echo=FALSE, message=FALSE, fig.cap="Gibbs energies of acetate oxidation and methanogenesis (after Mayumi et al., 2013).", cache=TRUE, pngquant=pngquant, timeit=timeit----
Adat <- lapply(c(-3, 3), function(logfCO2) {
basis("CO2", logfCO2)
data.frame(logfCO2,
rxnfun(c(0, 0))$A,
rxnfun(c(-1, 0))$A,
rxnfun(c(-1, 1))$A,
rxnfun(c(0, 1))$A
)
})
Adat <- do.call(rbind, Adat)
matplot(Adat[, 1], -Adat[, -1]/1000, type = "l", lty = 1, lwd = 2,
xlab = axis.label("CO2"), ylab = axis.label("DG", prefix = "k"))
legend("topleft", c("acetate oxidation", "acetoclastic methanogenesis",
"hydrogenotrophic methanogenesis"), lty = 1, col = 2:4)
## ----methanogenesis_plot, eval=FALSE------------------------------------------
# Adat <- lapply(c(-3, 3), function(logfCO2) {
# basis("CO2", logfCO2)
# data.frame(logfCO2,
# rxnfun(c(0, 0))$A,
# rxnfun(c(-1, 0))$A,
# rxnfun(c(-1, 1))$A,
# rxnfun(c(0, 1))$A
# )
# })
# Adat <- do.call(rbind, Adat)
# matplot(Adat[, 1], -Adat[, -1]/1000, type = "l", lty = 1, lwd = 2,
# xlab = axis.label("CO2"), ylab = axis.label("DG", prefix = "k"))
# legend("topleft", c("acetate oxidation", "acetoclastic methanogenesis",
# "hydrogenotrophic methanogenesis"), lty = 1, col = 2:4)
## ----reset, message=FALSE-----------------------------------------------------
reset()
## ----basis_CHNOSZ, results="hide"---------------------------------------------
basis("CHNOSe")
## ----species_sulfur-----------------------------------------------------------
species(c("H2S", "HS-", "HSO4-", "SO4-2"))
## ----affinity-----------------------------------------------------------------
unlist(affinity()$values)
## ----EhpH_plot, fig.margin=TRUE, fig.width=4, fig.height=4, out.width="100%", echo = FALSE, message=FALSE, cache=TRUE, fig.cap="Aqueous sulfur species at 25 °C.", pngquant=pngquant, timeit=timeit----
a <- affinity(pH = c(0, 12), Eh = c(-0.5, 1))
diagram(a, limit.water = TRUE)
## ----EhpH_plot, echo=TRUE, eval=FALSE-----------------------------------------
# a <- affinity(pH = c(0, 12), Eh = c(-0.5, 1))
# diagram(a, limit.water = TRUE)
## ----EhpH_plot_color, fig.margin=TRUE, fig.width=4, fig.height=4, smallish.mar=TRUE, out.width="100%", echo=FALSE, message=FALSE, cache=TRUE, fig.cap="The same plot, with different colors and labels.", pngquant=pngquant, timeit=timeit----
diagram(a, fill = "terrain", lwd = 2, lty = 3,
names = c("hydrogen sulfide", "bisulfide", "bisulfate", "sulfate"),
las = 0)
water.lines(a, col = 6, lwd = 2)
## ----EhpH_plot_color, echo=TRUE, eval=FALSE-----------------------------------
# diagram(a, fill = "terrain", lwd = 2, lty = 3,
# names = c("hydrogen sulfide", "bisulfide", "bisulfate", "sulfate"),
# las = 0)
# water.lines(a, col = 6, lwd = 2)
## ----retrieve-----------------------------------------------------------------
retrieve("Mn", c("O", "H"), "aq")
retrieve("Mn", c("O", "H"), "cr")
## ----retrieve_diagram, fig.margin=TRUE, fig.width=5, fig.height=5, out.width="100%", message=FALSE, results = "hide", cache=TRUE, fig.cap="Eh-pH diagram for the Mn-O-H system.", pngquant=pngquant, timeit=timeit----
# Set decimal logarithm of activity of aqueous species,
# temperature and plot resolution
logact <- -4
T <- 100
res <- 400
# Start with the aqueous species
basis(c("Mn+2", "H2O", "H+", "e-"))
iaq <- retrieve("Mn", c("O", "H"), "aq")
species(iaq, logact)
aaq <- affinity(pH = c(4, 16, res), Eh = c(-1.5, 1.5, res), T = T)
# Show names for only the metastable species here
names <- names(iaq)
names[!names(iaq) %in% c("MnOH+", "MnO", "HMnO2-")] <- ""
diagram(aaq, lty = 2, col = "#4169E188", names = names, col.names = 4)
# Overlay mineral stability fields
icr <- retrieve("Mn", c("O", "H"), "cr")
species(icr, add = TRUE)
# Supply the previous result from affinity() to use
# argument recall (for plotted variables and T)
acr <- affinity(aaq)
diagram(acr, add = TRUE, bold = acr$species$state=="cr", limit.water = FALSE)
# Add legend
legend <- c(
bquote(log * italic(a)["Mn(aq)"] == .(logact)),
bquote(italic(T) == .(T) ~ degree*C)
)
legend("topright", legend = as.expression(legend), bty = "n")
## ----info_CuCl, results="hide"------------------------------------------------
info(" CuCl")
## ----copper_setup, echo=TRUE, results="hide"----------------------------------
basis(c("Cu", "H2S", "Cl-", "H2O", "H+", "e-"))
basis("H2S", -6)
basis("Cl-", -0.7)
species(c("CuCl", "CuCl2-", "CuCl3-2", "CuCl+", "CuCl2", "CuCl3-", "CuCl4-2"))
species(c("chalcocite", "tenorite", "cuprite", "copper"), add = TRUE)
## ----info_chalcocite, message=FALSE-------------------------------------------
info(info("chalcocite", c("cr", "cr2", "cr3")))$T
## ----copper_mosaic, fig.margin=TRUE, fig.width=4, fig.height=4, out.width="100%", message=FALSE, cache=TRUE, fig.cap="Copper minerals and aqueous complexes with chloride, 200 °C.", pngquant=pngquant, timeit=timeit----
T <- 200
res <- 200
bases <- c("H2S", "HS-", "HSO4-", "SO4-2")
m1 <- mosaic(bases, pH = c(0, 12, res), Eh=c(-1.2, 0.75, res), T=T)
diagram(m1$A.species, lwd = 2)
diagram(m1$A.bases, add = TRUE, col = 4, col.names = 4, lty = 3,
italic = TRUE)
water.lines(m1$A.species, col = "blue1")
## ----rainbow_data-------------------------------------------------------------
file <- system.file("extdata/cpetc/SC10_Rainbow.csv", package = "CHNOSZ")
rb <- read.csv(file, check.names = FALSE)
## ----rainbow_species, results="hide"------------------------------------------
basis(c("CO2", "H2", "NH4+", "H2O", "H2S", "H+"))
species("CH4", -3)
species(c("adenine", "cytosine", "aspartic acid", "deoxyribose",
"CH4", "leucine", "tryptophan", "n-nonanoic acid"), -6)
## ----rainbow_affinity, message=FALSE------------------------------------------
a <- affinity(T = rb$T, CO2 = rb$CO2, H2 = rb$H2,
`NH4+` = rb$`NH4+`, H2S = rb$H2S, pH = rb$pH)
T <- convert(a$vals[[1]], "K")
a$values <- lapply(a$values, convert, "G", T)
a$values <- lapply(a$values, convert, "cal")
a$values <- lapply(a$values, `*`, -0.001)
## ----rainbow_diagram, fig.margin=TRUE, fig.width=4, fig.height=4, out.width="100%", echo=FALSE, message=FALSE, cache=TRUE, fig.cap="Affinities of organic synthesis in a hydrothermal system, after Shock and Canovas (2010).", pngquant=pngquant, timeit=timeit----
diagram(a, balance = 1, ylim = c(-100, 100), ylab = quote(italic(A)*", kcal/mol"),
col = rainbow(8), lwd = 2, bg = "slategray3")
abline(h = 0, lty = 2, lwd = 2)
## ----rainbow_diagram, eval=FALSE----------------------------------------------
# diagram(a, balance = 1, ylim = c(-100, 100), ylab = quote(italic(A)*", kcal/mol"),
# col = rainbow(8), lwd = 2, bg = "slategray3")
# abline(h = 0, lty = 2, lwd = 2)
## ----PPM_basis, results="hide", message=FALSE---------------------------------
basis(c("FeS2", "H2S", "O2", "H2O"))
species(c("pyrite", "magnetite"))
species("pyrrhotite", "cr2", add = TRUE)
## ----PPM_affinity, message=FALSE----------------------------------------------
unlist(affinity(T = 300, P = 100)$values)
## ----PPM_setup, results="hide"------------------------------------------------
mod.buffer("PPM", "pyrrhotite", "cr2")
basis(c("H2S", "O2"), c("PPM", "PPM"))
## ----PPM_activities, message=FALSE--------------------------------------------
unlist(affinity(T = 300, P = 100, return.buffer = TRUE)[1:3])
## ----demo_buffer_noecho, fig.margin=TRUE, fig.width=4, fig.height=4, out.width="100%", message=FALSE, echo=FALSE, cache=TRUE, fig.cap="Values of log<i>f</i><sub>H<sub>2</sub></sub> corresponding to mineral buffers or to given activities of aqueous species.", pngquant=pngquant, timeit=timeit----
demo(buffer, echo = FALSE)
## ----PPM_basis, echo = FALSE, results = "hide"--------------------------------
basis(c("FeS2", "H2S", "O2", "H2O"))
species(c("pyrite", "magnetite"))
species("pyrrhotite", "cr2", add = TRUE)
## ----PPM_setup, echo = FALSE, results = "hide", message = FALSE---------------
mod.buffer("PPM", "pyrrhotite", "cr2")
basis(c("H2S", "O2"), c("PPM", "PPM"))
## ----PPM_affinity, message = FALSE--------------------------------------------
unlist(affinity(T = 300, P = 100)$values)
## ----demo_buffer, eval=FALSE--------------------------------------------------
# demo(buffer)
## ----bjerrum_diagram, fig.margin=TRUE, fig.width=3, fig.height=6, out.width="100%", echo=FALSE, results="hide", message=FALSE, cache=TRUE, fig.cap="Three views of carbonate speciation: affinity, activity, degree of formation.", pngquant=pngquant, timeit=timeit----
par(mfrow = c(3, 1))
basis("CHNOS+")
species(c("CO2", "HCO3-", "CO3-2"))
a25 <- affinity(pH = c(4, 13))
a150 <- affinity(pH = c(4, 13), T = 150)
diagram(a25, dy = 0.4)
diagram(a150, add = TRUE, names = FALSE, col = "red")
e25 <- equilibrate(a25, loga.balance = -3)
e150 <- equilibrate(a150, loga.balance = -3)
diagram(e25, ylim = c(-6, 0), dy = 0.15)
diagram(e150, add = TRUE, names = FALSE, col = "red")
diagram(e25, alpha = TRUE, dy = -0.25)
diagram(e150, alpha = TRUE, add = TRUE, names = FALSE, col = "red")
## ----bjerrum_diagram, echo=1:7, eval=FALSE------------------------------------
# par(mfrow = c(3, 1))
# basis("CHNOS+")
# species(c("CO2", "HCO3-", "CO3-2"))
# a25 <- affinity(pH = c(4, 13))
# a150 <- affinity(pH = c(4, 13), T = 150)
# diagram(a25, dy = 0.4)
# diagram(a150, add = TRUE, names = FALSE, col = "red")
# e25 <- equilibrate(a25, loga.balance = -3)
# e150 <- equilibrate(a150, loga.balance = -3)
# diagram(e25, ylim = c(-6, 0), dy = 0.15)
# diagram(e150, add = TRUE, names = FALSE, col = "red")
# diagram(e25, alpha = TRUE, dy = -0.25)
# diagram(e150, alpha = TRUE, add = TRUE, names = FALSE, col = "red")
## ----bjerrum_diagram, echo=8:11, eval=FALSE-----------------------------------
# par(mfrow = c(3, 1))
# basis("CHNOS+")
# species(c("CO2", "HCO3-", "CO3-2"))
# a25 <- affinity(pH = c(4, 13))
# a150 <- affinity(pH = c(4, 13), T = 150)
# diagram(a25, dy = 0.4)
# diagram(a150, add = TRUE, names = FALSE, col = "red")
# e25 <- equilibrate(a25, loga.balance = -3)
# e150 <- equilibrate(a150, loga.balance = -3)
# diagram(e25, ylim = c(-6, 0), dy = 0.15)
# diagram(e150, add = TRUE, names = FALSE, col = "red")
# diagram(e25, alpha = TRUE, dy = -0.25)
# diagram(e150, alpha = TRUE, add = TRUE, names = FALSE, col = "red")
## ----bjerrum_diagram, echo=12:13, eval=FALSE----------------------------------
# par(mfrow = c(3, 1))
# basis("CHNOS+")
# species(c("CO2", "HCO3-", "CO3-2"))
# a25 <- affinity(pH = c(4, 13))
# a150 <- affinity(pH = c(4, 13), T = 150)
# diagram(a25, dy = 0.4)
# diagram(a150, add = TRUE, names = FALSE, col = "red")
# e25 <- equilibrate(a25, loga.balance = -3)
# e150 <- equilibrate(a150, loga.balance = -3)
# diagram(e25, ylim = c(-6, 0), dy = 0.15)
# diagram(e150, add = TRUE, names = FALSE, col = "red")
# diagram(e25, alpha = TRUE, dy = -0.25)
# diagram(e150, alpha = TRUE, add = TRUE, names = FALSE, col = "red")
## ----corundum, fig.margin=TRUE, fig.width=4, fig.height=4, out.width="100%", results="hide", message=FALSE, cache=TRUE, fig.cap="Solubility of corundum (green line) and equilibrium concentrations of aqueous species (black lines).", pngquant=pngquant, timeit=timeit----
add.OBIGT("SLOP98")
basis(c("Al+3", "H2O", "H+", "O2"))
species("corundum")
iaq <- c("Al+3", "AlO2-", "AlOH+2", "AlO+", "HAlO2")
s <- solubility(iaq, pH = c(0, 10), IS = 0, in.terms.of = "Al+3")
diagram(s, type = "loga.balance", ylim = c(-10, 0), lwd = 4, col = "green3")
diagram(s, add = TRUE, adj = c(0, 1, 2.1, -0.2, -1.5), dy = c(0, 0, 4, -0.3, 0.1))
legend("topright", c("25 °C", "1 bar"), text.font = 2, bty = "n")
reset()
## ----Alberty------------------------------------------------------------------
oldnon <- nonideal("Alberty")
## ----subcrt_IS----------------------------------------------------------------
subcrt(c("MgATP-2", "MgHATP-", "MgH2ATP"),
T = c(25, 100), IS = c(0, 0.25), property = "G")$out
## ----info_ATP, results="hide"-------------------------------------------------
info(" ATP")
## ----T_100--------------------------------------------------------------------
T <- 100
## ----ATP, eval=FALSE, echo=2:6------------------------------------------------
# par(mfrow = c(1, 4), mar = c(3.1, 3.6, 2.1, 1.6), mgp = c(1.8, 0.5, 0))
# basis("MgCHNOPS+")
# species(c("ATP-4", "HATP-3", "H2ATP-2", "H3ATP-", "H4ATP"))
# a <- affinity(pH = c(3, 9), T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, tplot = FALSE)
# title(main = describe.property("T", T))
# alphas <- do.call(rbind, d$plotvals)
# nH <- alphas * 0:4
# Hlab <- substitute(italic(N)[H^`+`])
# plot(a$vals[[1]], colSums(nH), type = "l", xlab = "pH", ylab=Hlab, lty=2, col=2)
# a <- affinity(pH = c(3, 9), IS = 0.25, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# nH <- alphas * 0:4
# lines(a$vals[[1]], colSums(nH))
# legend("topright", legend = c("I = 0 M", "I = 0.25 M"), lty = 2:1, col = 2:1, cex = 0.8)
# ATP.H <- substitute("ATP and H"^`+`)
# title(main = ATP.H)
# species(c("MgATP-2", "MgHATP-", "MgH2ATP", "Mg2ATP"), add = TRUE)
# Hplot <- function(pMg, IS = 0.25) {
# basis("Mg+2", -pMg)
# a <- affinity(pH = c(3, 9), IS = IS, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# NH <- alphas * c(0:4, 0, 1, 2, 0)
# lines(a$vals[[1]], colSums(NH), lty = 7 - pMg, col = 7 - pMg)
# }
# plot(c(3, 9), c(0, 2), type = "n", xlab = "pH", ylab = Hlab)
# lapply(2:6, Hplot)
# legend("topright", legend = paste("pMg = ", 2:6), lty = 5:1, col = 5:1, cex = 0.8)
# ATP.H.Mg <- substitute("ATP and H"^`+`~"and Mg"^`+2`)
# title(main = ATP.H.Mg)
# Mgplot <- function(pH, IS = 0.25) {
# basis("pH", pH)
# a <- affinity(`Mg+2` = c(-2, -7), IS = IS, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# NMg <- alphas * species()$`Mg+`
# lines(-a$vals[[1]], colSums(NMg), lty = 10 - pH, col = 10 - pH)
# }
# Mglab <- substitute(italic(N)[Mg^`+2`])
# plot(c(2, 7), c(0, 1.2), type = "n", xlab = "pMg", ylab = Mglab)
# lapply(3:9, Mgplot)
# legend("topright", legend = paste("pH = ", 3:9), lty = 7:1, col = 7:1, cex = 0.8)
# title(main = ATP.H.Mg)
## ----ATP, eval=FALSE, echo=8:11-----------------------------------------------
# par(mfrow = c(1, 4), mar = c(3.1, 3.6, 2.1, 1.6), mgp = c(1.8, 0.5, 0))
# basis("MgCHNOPS+")
# species(c("ATP-4", "HATP-3", "H2ATP-2", "H3ATP-", "H4ATP"))
# a <- affinity(pH = c(3, 9), T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, tplot = FALSE)
# title(main = describe.property("T", T))
# alphas <- do.call(rbind, d$plotvals)
# nH <- alphas * 0:4
# Hlab <- substitute(italic(N)[H^`+`])
# plot(a$vals[[1]], colSums(nH), type = "l", xlab = "pH", ylab=Hlab, lty=2, col=2)
# a <- affinity(pH = c(3, 9), IS = 0.25, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# nH <- alphas * 0:4
# lines(a$vals[[1]], colSums(nH))
# legend("topright", legend = c("I = 0 M", "I = 0.25 M"), lty = 2:1, col = 2:1, cex = 0.8)
# ATP.H <- substitute("ATP and H"^`+`)
# title(main = ATP.H)
# species(c("MgATP-2", "MgHATP-", "MgH2ATP", "Mg2ATP"), add = TRUE)
# Hplot <- function(pMg, IS = 0.25) {
# basis("Mg+2", -pMg)
# a <- affinity(pH = c(3, 9), IS = IS, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# NH <- alphas * c(0:4, 0, 1, 2, 0)
# lines(a$vals[[1]], colSums(NH), lty = 7 - pMg, col = 7 - pMg)
# }
# plot(c(3, 9), c(0, 2), type = "n", xlab = "pH", ylab = Hlab)
# lapply(2:6, Hplot)
# legend("topright", legend = paste("pMg = ", 2:6), lty = 5:1, col = 5:1, cex = 0.8)
# ATP.H.Mg <- substitute("ATP and H"^`+`~"and Mg"^`+2`)
# title(main = ATP.H.Mg)
# Mgplot <- function(pH, IS = 0.25) {
# basis("pH", pH)
# a <- affinity(`Mg+2` = c(-2, -7), IS = IS, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# NMg <- alphas * species()$`Mg+`
# lines(-a$vals[[1]], colSums(NMg), lty = 10 - pH, col = 10 - pH)
# }
# Mglab <- substitute(italic(N)[Mg^`+2`])
# plot(c(2, 7), c(0, 1.2), type = "n", xlab = "pMg", ylab = Mglab)
# lapply(3:9, Mgplot)
# legend("topright", legend = paste("pH = ", 3:9), lty = 7:1, col = 7:1, cex = 0.8)
# title(main = ATP.H.Mg)
## ----ATP, eval=FALSE, echo=12:17----------------------------------------------
# par(mfrow = c(1, 4), mar = c(3.1, 3.6, 2.1, 1.6), mgp = c(1.8, 0.5, 0))
# basis("MgCHNOPS+")
# species(c("ATP-4", "HATP-3", "H2ATP-2", "H3ATP-", "H4ATP"))
# a <- affinity(pH = c(3, 9), T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, tplot = FALSE)
# title(main = describe.property("T", T))
# alphas <- do.call(rbind, d$plotvals)
# nH <- alphas * 0:4
# Hlab <- substitute(italic(N)[H^`+`])
# plot(a$vals[[1]], colSums(nH), type = "l", xlab = "pH", ylab=Hlab, lty=2, col=2)
# a <- affinity(pH = c(3, 9), IS = 0.25, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# nH <- alphas * 0:4
# lines(a$vals[[1]], colSums(nH))
# legend("topright", legend = c("I = 0 M", "I = 0.25 M"), lty = 2:1, col = 2:1, cex = 0.8)
# ATP.H <- substitute("ATP and H"^`+`)
# title(main = ATP.H)
# species(c("MgATP-2", "MgHATP-", "MgH2ATP", "Mg2ATP"), add = TRUE)
# Hplot <- function(pMg, IS = 0.25) {
# basis("Mg+2", -pMg)
# a <- affinity(pH = c(3, 9), IS = IS, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# NH <- alphas * c(0:4, 0, 1, 2, 0)
# lines(a$vals[[1]], colSums(NH), lty = 7 - pMg, col = 7 - pMg)
# }
# plot(c(3, 9), c(0, 2), type = "n", xlab = "pH", ylab = Hlab)
# lapply(2:6, Hplot)
# legend("topright", legend = paste("pMg = ", 2:6), lty = 5:1, col = 5:1, cex = 0.8)
# ATP.H.Mg <- substitute("ATP and H"^`+`~"and Mg"^`+2`)
# title(main = ATP.H.Mg)
# Mgplot <- function(pH, IS = 0.25) {
# basis("pH", pH)
# a <- affinity(`Mg+2` = c(-2, -7), IS = IS, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# NMg <- alphas * species()$`Mg+`
# lines(-a$vals[[1]], colSums(NMg), lty = 10 - pH, col = 10 - pH)
# }
# Mglab <- substitute(italic(N)[Mg^`+2`])
# plot(c(2, 7), c(0, 1.2), type = "n", xlab = "pMg", ylab = Mglab)
# lapply(3:9, Mgplot)
# legend("topright", legend = paste("pH = ", 3:9), lty = 7:1, col = 7:1, cex = 0.8)
# title(main = ATP.H.Mg)
## ----ATP, eval=FALSE, echo=21-------------------------------------------------
# par(mfrow = c(1, 4), mar = c(3.1, 3.6, 2.1, 1.6), mgp = c(1.8, 0.5, 0))
# basis("MgCHNOPS+")
# species(c("ATP-4", "HATP-3", "H2ATP-2", "H3ATP-", "H4ATP"))
# a <- affinity(pH = c(3, 9), T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, tplot = FALSE)
# title(main = describe.property("T", T))
# alphas <- do.call(rbind, d$plotvals)
# nH <- alphas * 0:4
# Hlab <- substitute(italic(N)[H^`+`])
# plot(a$vals[[1]], colSums(nH), type = "l", xlab = "pH", ylab=Hlab, lty=2, col=2)
# a <- affinity(pH = c(3, 9), IS = 0.25, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# nH <- alphas * 0:4
# lines(a$vals[[1]], colSums(nH))
# legend("topright", legend = c("I = 0 M", "I = 0.25 M"), lty = 2:1, col = 2:1, cex = 0.8)
# ATP.H <- substitute("ATP and H"^`+`)
# title(main = ATP.H)
# species(c("MgATP-2", "MgHATP-", "MgH2ATP", "Mg2ATP"), add = TRUE)
# Hplot <- function(pMg, IS = 0.25) {
# basis("Mg+2", -pMg)
# a <- affinity(pH = c(3, 9), IS = IS, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# NH <- alphas * c(0:4, 0, 1, 2, 0)
# lines(a$vals[[1]], colSums(NH), lty = 7 - pMg, col = 7 - pMg)
# }
# plot(c(3, 9), c(0, 2), type = "n", xlab = "pH", ylab = Hlab)
# lapply(2:6, Hplot)
# legend("topright", legend = paste("pMg = ", 2:6), lty = 5:1, col = 5:1, cex = 0.8)
# ATP.H.Mg <- substitute("ATP and H"^`+`~"and Mg"^`+2`)
# title(main = ATP.H.Mg)
# Mgplot <- function(pH, IS = 0.25) {
# basis("pH", pH)
# a <- affinity(`Mg+2` = c(-2, -7), IS = IS, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# NMg <- alphas * species()$`Mg+`
# lines(-a$vals[[1]], colSums(NMg), lty = 10 - pH, col = 10 - pH)
# }
# Mglab <- substitute(italic(N)[Mg^`+2`])
# plot(c(2, 7), c(0, 1.2), type = "n", xlab = "pMg", ylab = Mglab)
# lapply(3:9, Mgplot)
# legend("topright", legend = paste("pH = ", 3:9), lty = 7:1, col = 7:1, cex = 0.8)
# title(main = ATP.H.Mg)
## ----ATP, eval=FALSE, echo=22:30----------------------------------------------
# par(mfrow = c(1, 4), mar = c(3.1, 3.6, 2.1, 1.6), mgp = c(1.8, 0.5, 0))
# basis("MgCHNOPS+")
# species(c("ATP-4", "HATP-3", "H2ATP-2", "H3ATP-", "H4ATP"))
# a <- affinity(pH = c(3, 9), T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, tplot = FALSE)
# title(main = describe.property("T", T))
# alphas <- do.call(rbind, d$plotvals)
# nH <- alphas * 0:4
# Hlab <- substitute(italic(N)[H^`+`])
# plot(a$vals[[1]], colSums(nH), type = "l", xlab = "pH", ylab=Hlab, lty=2, col=2)
# a <- affinity(pH = c(3, 9), IS = 0.25, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# nH <- alphas * 0:4
# lines(a$vals[[1]], colSums(nH))
# legend("topright", legend = c("I = 0 M", "I = 0.25 M"), lty = 2:1, col = 2:1, cex = 0.8)
# ATP.H <- substitute("ATP and H"^`+`)
# title(main = ATP.H)
# species(c("MgATP-2", "MgHATP-", "MgH2ATP", "Mg2ATP"), add = TRUE)
# Hplot <- function(pMg, IS = 0.25) {
# basis("Mg+2", -pMg)
# a <- affinity(pH = c(3, 9), IS = IS, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# NH <- alphas * c(0:4, 0, 1, 2, 0)
# lines(a$vals[[1]], colSums(NH), lty = 7 - pMg, col = 7 - pMg)
# }
# plot(c(3, 9), c(0, 2), type = "n", xlab = "pH", ylab = Hlab)
# lapply(2:6, Hplot)
# legend("topright", legend = paste("pMg = ", 2:6), lty = 5:1, col = 5:1, cex = 0.8)
# ATP.H.Mg <- substitute("ATP and H"^`+`~"and Mg"^`+2`)
# title(main = ATP.H.Mg)
# Mgplot <- function(pH, IS = 0.25) {
# basis("pH", pH)
# a <- affinity(`Mg+2` = c(-2, -7), IS = IS, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# NMg <- alphas * species()$`Mg+`
# lines(-a$vals[[1]], colSums(NMg), lty = 10 - pH, col = 10 - pH)
# }
# Mglab <- substitute(italic(N)[Mg^`+2`])
# plot(c(2, 7), c(0, 1.2), type = "n", xlab = "pMg", ylab = Mglab)
# lapply(3:9, Mgplot)
# legend("topright", legend = paste("pH = ", 3:9), lty = 7:1, col = 7:1, cex = 0.8)
# title(main = ATP.H.Mg)
## ----ATP, eval=FALSE, echo=31:32----------------------------------------------
# par(mfrow = c(1, 4), mar = c(3.1, 3.6, 2.1, 1.6), mgp = c(1.8, 0.5, 0))
# basis("MgCHNOPS+")
# species(c("ATP-4", "HATP-3", "H2ATP-2", "H3ATP-", "H4ATP"))
# a <- affinity(pH = c(3, 9), T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, tplot = FALSE)
# title(main = describe.property("T", T))
# alphas <- do.call(rbind, d$plotvals)
# nH <- alphas * 0:4
# Hlab <- substitute(italic(N)[H^`+`])
# plot(a$vals[[1]], colSums(nH), type = "l", xlab = "pH", ylab=Hlab, lty=2, col=2)
# a <- affinity(pH = c(3, 9), IS = 0.25, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# nH <- alphas * 0:4
# lines(a$vals[[1]], colSums(nH))
# legend("topright", legend = c("I = 0 M", "I = 0.25 M"), lty = 2:1, col = 2:1, cex = 0.8)
# ATP.H <- substitute("ATP and H"^`+`)
# title(main = ATP.H)
# species(c("MgATP-2", "MgHATP-", "MgH2ATP", "Mg2ATP"), add = TRUE)
# Hplot <- function(pMg, IS = 0.25) {
# basis("Mg+2", -pMg)
# a <- affinity(pH = c(3, 9), IS = IS, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# NH <- alphas * c(0:4, 0, 1, 2, 0)
# lines(a$vals[[1]], colSums(NH), lty = 7 - pMg, col = 7 - pMg)
# }
# plot(c(3, 9), c(0, 2), type = "n", xlab = "pH", ylab = Hlab)
# lapply(2:6, Hplot)
# legend("topright", legend = paste("pMg = ", 2:6), lty = 5:1, col = 5:1, cex = 0.8)
# ATP.H.Mg <- substitute("ATP and H"^`+`~"and Mg"^`+2`)
# title(main = ATP.H.Mg)
# Mgplot <- function(pH, IS = 0.25) {
# basis("pH", pH)
# a <- affinity(`Mg+2` = c(-2, -7), IS = IS, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# NMg <- alphas * species()$`Mg+`
# lines(-a$vals[[1]], colSums(NMg), lty = 10 - pH, col = 10 - pH)
# }
# Mglab <- substitute(italic(N)[Mg^`+2`])
# plot(c(2, 7), c(0, 1.2), type = "n", xlab = "pMg", ylab = Mglab)
# lapply(3:9, Mgplot)
# legend("topright", legend = paste("pH = ", 3:9), lty = 7:1, col = 7:1, cex = 0.8)
# title(main = ATP.H.Mg)
## ----ATP, eval=FALSE, echo=36:44----------------------------------------------
# par(mfrow = c(1, 4), mar = c(3.1, 3.6, 2.1, 1.6), mgp = c(1.8, 0.5, 0))
# basis("MgCHNOPS+")
# species(c("ATP-4", "HATP-3", "H2ATP-2", "H3ATP-", "H4ATP"))
# a <- affinity(pH = c(3, 9), T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, tplot = FALSE)
# title(main = describe.property("T", T))
# alphas <- do.call(rbind, d$plotvals)
# nH <- alphas * 0:4
# Hlab <- substitute(italic(N)[H^`+`])
# plot(a$vals[[1]], colSums(nH), type = "l", xlab = "pH", ylab=Hlab, lty=2, col=2)
# a <- affinity(pH = c(3, 9), IS = 0.25, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# nH <- alphas * 0:4
# lines(a$vals[[1]], colSums(nH))
# legend("topright", legend = c("I = 0 M", "I = 0.25 M"), lty = 2:1, col = 2:1, cex = 0.8)
# ATP.H <- substitute("ATP and H"^`+`)
# title(main = ATP.H)
# species(c("MgATP-2", "MgHATP-", "MgH2ATP", "Mg2ATP"), add = TRUE)
# Hplot <- function(pMg, IS = 0.25) {
# basis("Mg+2", -pMg)
# a <- affinity(pH = c(3, 9), IS = IS, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# NH <- alphas * c(0:4, 0, 1, 2, 0)
# lines(a$vals[[1]], colSums(NH), lty = 7 - pMg, col = 7 - pMg)
# }
# plot(c(3, 9), c(0, 2), type = "n", xlab = "pH", ylab = Hlab)
# lapply(2:6, Hplot)
# legend("topright", legend = paste("pMg = ", 2:6), lty = 5:1, col = 5:1, cex = 0.8)
# ATP.H.Mg <- substitute("ATP and H"^`+`~"and Mg"^`+2`)
# title(main = ATP.H.Mg)
# Mgplot <- function(pH, IS = 0.25) {
# basis("pH", pH)
# a <- affinity(`Mg+2` = c(-2, -7), IS = IS, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# NMg <- alphas * species()$`Mg+`
# lines(-a$vals[[1]], colSums(NMg), lty = 10 - pH, col = 10 - pH)
# }
# Mglab <- substitute(italic(N)[Mg^`+2`])
# plot(c(2, 7), c(0, 1.2), type = "n", xlab = "pMg", ylab = Mglab)
# lapply(3:9, Mgplot)
# legend("topright", legend = paste("pH = ", 3:9), lty = 7:1, col = 7:1, cex = 0.8)
# title(main = ATP.H.Mg)
## ----ATP, eval=FALSE, echo=45:47----------------------------------------------
# par(mfrow = c(1, 4), mar = c(3.1, 3.6, 2.1, 1.6), mgp = c(1.8, 0.5, 0))
# basis("MgCHNOPS+")
# species(c("ATP-4", "HATP-3", "H2ATP-2", "H3ATP-", "H4ATP"))
# a <- affinity(pH = c(3, 9), T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, tplot = FALSE)
# title(main = describe.property("T", T))
# alphas <- do.call(rbind, d$plotvals)
# nH <- alphas * 0:4
# Hlab <- substitute(italic(N)[H^`+`])
# plot(a$vals[[1]], colSums(nH), type = "l", xlab = "pH", ylab=Hlab, lty=2, col=2)
# a <- affinity(pH = c(3, 9), IS = 0.25, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# nH <- alphas * 0:4
# lines(a$vals[[1]], colSums(nH))
# legend("topright", legend = c("I = 0 M", "I = 0.25 M"), lty = 2:1, col = 2:1, cex = 0.8)
# ATP.H <- substitute("ATP and H"^`+`)
# title(main = ATP.H)
# species(c("MgATP-2", "MgHATP-", "MgH2ATP", "Mg2ATP"), add = TRUE)
# Hplot <- function(pMg, IS = 0.25) {
# basis("Mg+2", -pMg)
# a <- affinity(pH = c(3, 9), IS = IS, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# NH <- alphas * c(0:4, 0, 1, 2, 0)
# lines(a$vals[[1]], colSums(NH), lty = 7 - pMg, col = 7 - pMg)
# }
# plot(c(3, 9), c(0, 2), type = "n", xlab = "pH", ylab = Hlab)
# lapply(2:6, Hplot)
# legend("topright", legend = paste("pMg = ", 2:6), lty = 5:1, col = 5:1, cex = 0.8)
# ATP.H.Mg <- substitute("ATP and H"^`+`~"and Mg"^`+2`)
# title(main = ATP.H.Mg)
# Mgplot <- function(pH, IS = 0.25) {
# basis("pH", pH)
# a <- affinity(`Mg+2` = c(-2, -7), IS = IS, T = T)
# e <- equilibrate(a)
# d <- diagram(e, alpha = TRUE, plot.it = FALSE)
# alphas <- do.call(rbind, d$plotvals)
# NMg <- alphas * species()$`Mg+`
# lines(-a$vals[[1]], colSums(NMg), lty = 10 - pH, col = 10 - pH)
# }
# Mglab <- substitute(italic(N)[Mg^`+2`])
# plot(c(2, 7), c(0, 1.2), type = "n", xlab = "pMg", ylab = Mglab)
# lapply(3:9, Mgplot)
# legend("topright", legend = paste("pH = ", 3:9), lty = 7:1, col = 7:1, cex = 0.8)
# title(main = ATP.H.Mg)
## ----ATP, fig.fullwidth=TRUE, fig.width=10, fig.height=2.5, dpi=hidpi, out.width="100%", echo=FALSE, message=FALSE, results="hide", fig.cap="Binding of H<sup>+</sup> and Mg<sup>+2</sup> to ATP at 100 °C and *I* = 0 M (first plot) or *I* = 0.25 M (third and fourth plots).", cache=TRUE, pngquant=pngquant, timeit=timeit----
par(mfrow = c(1, 4), mar = c(3.1, 3.6, 2.1, 1.6), mgp = c(1.8, 0.5, 0))
basis("MgCHNOPS+")
species(c("ATP-4", "HATP-3", "H2ATP-2", "H3ATP-", "H4ATP"))
a <- affinity(pH = c(3, 9), T = T)
e <- equilibrate(a)
d <- diagram(e, alpha = TRUE, tplot = FALSE)
title(main = describe.property("T", T))
alphas <- do.call(rbind, d$plotvals)
nH <- alphas * 0:4
Hlab <- substitute(italic(N)[H^`+`])
plot(a$vals[[1]], colSums(nH), type = "l", xlab = "pH", ylab=Hlab, lty=2, col=2)
a <- affinity(pH = c(3, 9), IS = 0.25, T = T)
e <- equilibrate(a)
d <- diagram(e, alpha = TRUE, plot.it = FALSE)
alphas <- do.call(rbind, d$plotvals)
nH <- alphas * 0:4
lines(a$vals[[1]], colSums(nH))
legend("topright", legend = c("I = 0 M", "I = 0.25 M"), lty = 2:1, col = 2:1, cex = 0.8)
ATP.H <- substitute("ATP and H"^`+`)
title(main = ATP.H)
species(c("MgATP-2", "MgHATP-", "MgH2ATP", "Mg2ATP"), add = TRUE)
Hplot <- function(pMg, IS = 0.25) {
basis("Mg+2", -pMg)
a <- affinity(pH = c(3, 9), IS = IS, T = T)
e <- equilibrate(a)
d <- diagram(e, alpha = TRUE, plot.it = FALSE)
alphas <- do.call(rbind, d$plotvals)
NH <- alphas * c(0:4, 0, 1, 2, 0)
lines(a$vals[[1]], colSums(NH), lty = 7 - pMg, col = 7 - pMg)
}
plot(c(3, 9), c(0, 2), type = "n", xlab = "pH", ylab = Hlab)
lapply(2:6, Hplot)
legend("topright", legend = paste("pMg = ", 2:6), lty = 5:1, col = 5:1, cex = 0.8)
ATP.H.Mg <- substitute("ATP and H"^`+`~"and Mg"^`+2`)
title(main = ATP.H.Mg)
Mgplot <- function(pH, IS = 0.25) {
basis("pH", pH)
a <- affinity(`Mg+2` = c(-2, -7), IS = IS, T = T)
e <- equilibrate(a)
d <- diagram(e, alpha = TRUE, plot.it = FALSE)
alphas <- do.call(rbind, d$plotvals)
NMg <- alphas * species()$`Mg+`
lines(-a$vals[[1]], colSums(NMg), lty = 10 - pH, col = 10 - pH)
}
Mglab <- substitute(italic(N)[Mg^`+2`])
plot(c(2, 7), c(0, 1.2), type = "n", xlab = "pMg", ylab = Mglab)
lapply(3:9, Mgplot)
legend("topright", legend = paste("pH = ", 3:9), lty = 7:1, col = 7:1, cex = 0.8)
title(main = ATP.H.Mg)
## ----oldnon-------------------------------------------------------------------
nonideal(oldnon)
## ----pinfo_LYSC_CHICK---------------------------------------------------------
p1 <- pinfo("LYSC_CHICK")
p2 <- pinfo(c("SHH", "OLIG2"), "HUMAN")
pinfo(c(p1, p2))
## ----formula_LYSC_CHICK-------------------------------------------------------
pl <- protein.length("LYSC_CHICK")
pf <- protein.formula("LYSC_CHICK")
list(length = pl, protein = pf, residue = pf / pl,
ZC_protein = ZC(pf), ZC_residue = ZC(pf / pl))
## ----subcrt_LYSC_CHICK, message=FALSE-----------------------------------------
subcrt("LYSC_CHICK")$out[[1]][1:6, ]
## ----protein_Cp, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, out.width="100%", echo=FALSE, message=FALSE, fig.cap='The heat capacity calculated by group additivity closely approximates experimental values for aqueous proteins. For a related figure showing the effects of ionization in the calculations, see <span style="color:blue">`?ionize.aa`</span>.', cache=TRUE, pngquant=pngquant, timeit=timeit----
PM90 <- read.csv(system.file("extdata/cpetc/PM90.csv", package = "CHNOSZ"))
plength <- protein.length(colnames(PM90)[2:5])
Cp_expt <- t(t(PM90[, 2:5]) / plength)
matplot(PM90[, 1], Cp_expt, type = "p", pch = 19,
xlab = axis.label("T"), ylab = axis.label("Cp0"), ylim = c(110, 280))
for(i in 1:4) {
pname <- colnames(Cp_expt)[i]
aq <- subcrt(pname, "aq", T = seq(0, 150))$out[[1]]
cr <- subcrt(pname, "cr", T = seq(0, 150))$out[[1]]
lines(aq$T, aq$Cp / plength[i], col = i)
lines(cr$T, cr$Cp / plength[i], col = i, lty = 2)
}
legend("right", legend = colnames(Cp_expt),
col = 1:4, pch = 19, lty = 1, bty = "n", cex = 0.9)
legend("bottomright", legend = c("experimental", "calculated (aq)",
"calculated (cr)"), lty = c(NA, 1, 2), pch = c(19, NA, NA), bty = "n")
## ----protein_Cp, eval=FALSE, echo=1:5-----------------------------------------
# PM90 <- read.csv(system.file("extdata/cpetc/PM90.csv", package = "CHNOSZ"))
# plength <- protein.length(colnames(PM90)[2:5])
# Cp_expt <- t(t(PM90[, 2:5]) / plength)
# matplot(PM90[, 1], Cp_expt, type = "p", pch = 19,
# xlab = axis.label("T"), ylab = axis.label("Cp0"), ylim = c(110, 280))
# for(i in 1:4) {
# pname <- colnames(Cp_expt)[i]
# aq <- subcrt(pname, "aq", T = seq(0, 150))$out[[1]]
# cr <- subcrt(pname, "cr", T = seq(0, 150))$out[[1]]
# lines(aq$T, aq$Cp / plength[i], col = i)
# lines(cr$T, cr$Cp / plength[i], col = i, lty = 2)
# }
# legend("right", legend = colnames(Cp_expt),
# col = 1:4, pch = 19, lty = 1, bty = "n", cex = 0.9)
# legend("bottomright", legend = c("experimental", "calculated (aq)",
# "calculated (cr)"), lty = c(NA, 1, 2), pch = c(19, NA, NA), bty = "n")
## ----protein_Cp, eval=FALSE, echo=-(1:5)--------------------------------------
# PM90 <- read.csv(system.file("extdata/cpetc/PM90.csv", package = "CHNOSZ"))
# plength <- protein.length(colnames(PM90)[2:5])
# Cp_expt <- t(t(PM90[, 2:5]) / plength)
# matplot(PM90[, 1], Cp_expt, type = "p", pch = 19,
# xlab = axis.label("T"), ylab = axis.label("Cp0"), ylim = c(110, 280))
# for(i in 1:4) {
# pname <- colnames(Cp_expt)[i]
# aq <- subcrt(pname, "aq", T = seq(0, 150))$out[[1]]
# cr <- subcrt(pname, "cr", T = seq(0, 150))$out[[1]]
# lines(aq$T, aq$Cp / plength[i], col = i)
# lines(cr$T, cr$Cp / plength[i], col = i, lty = 2)
# }
# legend("right", legend = colnames(Cp_expt),
# col = 1:4, pch = 19, lty = 1, bty = "n", cex = 0.9)
# legend("bottomright", legend = c("experimental", "calculated (aq)",
# "calculated (cr)"), lty = c(NA, 1, 2), pch = c(19, NA, NA), bty = "n")
## ----protein_ionization, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, out.width="100%", echo=FALSE, results="hide", message=FALSE, fig.cap='Affinity of ionization of proteins. See [<span style="color:blue">demo(ionize)</span>](../demo) for ionization properties calculated as a function of temperature and pH.', cache=TRUE, pngquant=pngquant, timeit=timeit----
ip <- pinfo(c("CYC_BOVIN", "LYSC_CHICK", "MYG_PHYCA", "RNAS1_BOVIN"))
basis("CHNOS+")
a_ion <- affinity(pH = c(0, 14), iprotein = ip)
basis("CHNOS")
a_nonion <- affinity(iprotein = ip)
plot(c(0, 14), c(50, 300), xlab = "pH", ylab = quote(italic(A/2.303*RT)), type="n")
for(i in 1:4) {
A_ion <- as.numeric(a_ion$values[[i]])
A_nonion <- as.numeric(a_nonion$values[[i]])
lines(a_ion$vals[[1]], A_ion - A_nonion, col=i)
}
legend("topright", legend = a_ion$species$name,
col = 1:4, lty = 1, bty = "n", cex = 0.9)
## ----protein_ionization, eval=FALSE-------------------------------------------
# ip <- pinfo(c("CYC_BOVIN", "LYSC_CHICK", "MYG_PHYCA", "RNAS1_BOVIN"))
# basis("CHNOS+")
# a_ion <- affinity(pH = c(0, 14), iprotein = ip)
# basis("CHNOS")
# a_nonion <- affinity(iprotein = ip)
# plot(c(0, 14), c(50, 300), xlab = "pH", ylab = quote(italic(A/2.303*RT)), type="n")
# for(i in 1:4) {
# A_ion <- as.numeric(a_ion$values[[i]])
# A_nonion <- as.numeric(a_nonion$values[[i]])
# lines(a_ion$vals[[1]], A_ion - A_nonion, col=i)
# }
# legend("topright", legend = a_ion$species$name,
# col = 1:4, lty = 1, bty = "n", cex = 0.9)
## ----basis_CHNOS, echo=FALSE, results="hide"----------------------------------
basis("CHNOS")
## ----species_protein, results="hide", message=FALSE, echo=1:2-----------------
species(c("CYC_BOVIN", "LYSC_CHICK", "MYG_PHYCA", "RNAS1_BOVIN"))
a_nonion_species <- affinity()
## ----nonion_species_values----------------------------------------------------
unlist(a_nonion_species$values)
## ----nonion_values------------------------------------------------------------
unlist(a_nonion$values)
## ----rubisco_svg, echo=FALSE, results="hide", fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, out.width="100%", fig.ext='svg', custom.plot=TRUE, embed.tag=TRUE, fig.cap='Average oxidation state of carbon in Rubisco compared with optimal growth temperature of organisms. **This is an interactive image.** Move the mouse over the points to show the names of the organisms, and click to open a reference in a new window. (Made with [**RSVGTipsDevice**](https://cran.r-project.org/package=RSVGTipsDevice) using code that can be found in the source of this document.)'----
## copies the premade SVG image to the knitr figure path
file.copy("rubisco.svg", fig_path(".svg"))
## the code for making the SVG image -- not "live" in the vignette because RSVGTipsDevice isn't available on Windows
#if(require(RSVGTipsDevice)) {
# datfile <- system.file("extdata/cpetc/rubisco.csv", package = "CHNOSZ")
# fastafile <- system.file("extdata/protein/rubisco.fasta", package = "CHNOSZ")
# dat <- read.csv(datfile)
# aa <- read.fasta(fastafile)
# Topt <- (dat$T1 + dat$T2) / 2
# idat <- match(dat$ID, substr(aa$protein, 4, 9))
# aa <- aa[idat, ]
# ZC <- ZC(protein.formula(aa))
# pch <- match(dat$domain, c("E", "B", "A")) - 1
# col <- match(dat$domain, c("A", "B", "E")) + 1
# # because the tooltip titles in the SVG file are shown by recent browsers,
# # we do not need to draw the tooltips explicitly, so set toolTipMode=0
# devSVGTips("rubisco.svg", toolTipMode=0, title="Rubisco")
# par(cex=1.4)
# # unfortunately, plotmath can't be used with devSVGTips,
# # so axis labels here don't contain italics.
# plot(Topt, ZC, type="n", xlab="T, °C", ylab="ZC")
# n <- rep(1:9, 3)
# for(i in seq_along(Topt)) {
# # adjust cex to make the symbols look the same size
# cex <- ifelse(pch[i]==1, 2.5, 3.5)
# points(Topt[i], ZC[i], pch=pch[i], cex=cex, col=col[i])
# URL <- dat$URL[i]
# setSVGShapeURL(URL, target="_blank")
# setSVGShapeContents(paste0("<title>", dat$species[i], "</title>"))
# text(Topt[i], ZC[i], n[i], cex = 1.2)
# }
# abline(v = c(36, 63), lty = 2, col = "grey")
# legend("topright", legend = c("Archaea", "Bacteria", "Eukaryota"),
# pch = c(2, 1, 0), col = 2:4, cex=1.5, pt.cex = c(3, 2.3, 3), bty="n")
# dev.off()
#}
## ----rubisco_ZC, fig.keep="none", message=FALSE-------------------------------
datfile <- system.file("extdata/cpetc/rubisco.csv", package = "CHNOSZ")
fastafile <- system.file("extdata/protein/rubisco.fasta", package = "CHNOSZ")
dat <- read.csv(datfile)
aa <- read.fasta(fastafile)
Topt <- (dat$T1 + dat$T2) / 2
idat <- match(dat$ID, substr(aa$protein, 4, 9))
aa <- aa[idat, ]
ZC <- ZC(protein.formula(aa))
pch <- match(dat$domain, c("E", "B", "A")) - 1
col <- match(dat$domain, c("A", "B", "E")) + 1
plot(Topt, ZC, pch = pch, cex = 2, col = col,
xlab = expression(list(italic(T)[opt], degree*C)),
ylab = expression(italic(Z)[C]))
text(Topt, ZC, rep(1:9, 3), cex = 0.8)
abline(v = c(36, 63), lty = 2, col = "grey")
legend("topright", legend = c("Archaea", "Bacteria", "Eukaryota"),
pch = c(2, 1, 0), col = 2:4, pt.cex = 2)
## ----rubisco_O2, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, out.width="100%", echo=FALSE, results="hide", message=FALSE, fig.cap="Elemental compositions of proteins projected into different sets of basis species.", cache=TRUE, pngquant=pngquant, timeit=timeit----
layout(matrix(1:4, nrow = 2))
par(mgp = c(1.8, 0.5, 0))
pl <- protein.length(aa)
ZClab <- axis.label("ZC")
nO2lab <- expression(bar(italic(n))[O[2]])
nH2Olab <- expression(bar(italic(n))[H[2]*O])
lapply(c("CHNOS", "QEC"), function(thisbasis) {
basis(thisbasis)
pb <- protein.basis(aa)
nO2 <- pb[, "O2"] / pl
plot(ZC, nO2, pch = pch, col = col, xlab = ZClab, ylab = nO2lab)
nH2O <- pb[, "H2O"] / pl
plot(ZC, nH2O, pch = pch, col = col, xlab = ZClab, ylab = nH2Olab)
mtext(thisbasis, font = 2)
})
## ----rubisco_O2, eval=FALSE---------------------------------------------------
# layout(matrix(1:4, nrow = 2))
# par(mgp = c(1.8, 0.5, 0))
# pl <- protein.length(aa)
# ZClab <- axis.label("ZC")
# nO2lab <- expression(bar(italic(n))[O[2]])
# nH2Olab <- expression(bar(italic(n))[H[2]*O])
# lapply(c("CHNOS", "QEC"), function(thisbasis) {
# basis(thisbasis)
# pb <- protein.basis(aa)
# nO2 <- pb[, "O2"] / pl
# plot(ZC, nO2, pch = pch, col = col, xlab = ZClab, ylab = nO2lab)
# nH2O <- pb[, "H2O"] / pl
# plot(ZC, nH2O, pch = pch, col = col, xlab = ZClab, ylab = nH2Olab)
# mtext(thisbasis, font = 2)
# })
## ----yeastgfp-----------------------------------------------------------------
protein <- c("YDL195W", "YHR098C", "YIL109C", "YLR208W", "YNL049C", "YPL085W")
abundance <- c(1840, 12200, NA, 21400, 1720, 358)
ina <- is.na(abundance)
## ----add_protein_yeast, message=FALSE-----------------------------------------
ip <- match(protein[!ina], thermo()$protein$protein)
## ----unitize------------------------------------------------------------------
pl <- protein.length(ip)
logact <- unitize(numeric(5), pl)
logabundance <- unitize(log10(abundance[!ina]), pl)
## ----yeastplot, eval=FALSE, echo=1:5------------------------------------------
# par(mfrow = c(1, 3))
# basis("CHNOS+")
# a <- affinity(O2 = c(-80, -73), iprotein = ip, loga.protein = logact)
# e <- equilibrate(a)
# diagram(e, ylim = c(-5, -2), col = 1:5, lwd = 2)
# e <- equilibrate(a, normalize = TRUE)
# diagram(e, ylim = c(-5, -2.5), col = 1:5, lwd = 2)
# abline(h = logabundance, lty = 1:5, col = 1:5)
## ----yeastplot, eval=FALSE, echo=6:8------------------------------------------
# par(mfrow = c(1, 3))
# basis("CHNOS+")
# a <- affinity(O2 = c(-80, -73), iprotein = ip, loga.protein = logact)
# e <- equilibrate(a)
# diagram(e, ylim = c(-5, -2), col = 1:5, lwd = 2)
# e <- equilibrate(a, normalize = TRUE)
# diagram(e, ylim = c(-5, -2.5), col = 1:5, lwd = 2)
# abline(h = logabundance, lty = 1:5, col = 1:5)
## ----yeastplot, fig.width=6, fig.height=2.5, dpi=hidpi, out.width="100%", echo=FALSE, message=FALSE, results="hide", cache=TRUE, fig.cap="ER-to-Golgi proteins: calculations without and with length normalization.", pngquant=pngquant, timeit=timeit----
par(mfrow = c(1, 3))
basis("CHNOS+")
a <- affinity(O2 = c(-80, -73), iprotein = ip, loga.protein = logact)
e <- equilibrate(a)
diagram(e, ylim = c(-5, -2), col = 1:5, lwd = 2)
e <- equilibrate(a, normalize = TRUE)
diagram(e, ylim = c(-5, -2.5), col = 1:5, lwd = 2)
abline(h = logabundance, lty = 1:5, col = 1:5)
## ----read_csv-----------------------------------------------------------------
file <- system.file("extdata/protein/POLG.csv", package = "CHNOSZ")
aa_POLG <- read.csv(file, as.is = TRUE, nrows = 5)
## ----read_fasta, message=FALSE------------------------------------------------
file <- system.file("extdata/protein/EF-Tu.aln", package = "CHNOSZ")
aa_Ef <- read.fasta(file, iseq = 1:2)
## ----seq2aa-------------------------------------------------------------------
aa_PRIO <- seq2aa("
MANLGCWMLVLFVATWSDLGLCKKRPKPGGWNTGGSRYPGQGSPGGNRYPPQGGGGWGQP
HGGGWGQPHGGGWGQPHGGGWGQPHGGGWGQGGGTHSQWNKPSKPKTNMKHMAGAAAAGA
VVGGLGGYMLGSAMSRPIIHFGSDYEDRYYRENMHRYPNQVYYRPMDEYSNQNNFVHDCV
NITIKQHTVTTTTKGENFTETDVKMMERVVEQMCITQYERESQAYYQRGSSMVLFSSPPV
ILLISFLIFLIVG
", "PRIO_HUMAN")
## ----protein_length-----------------------------------------------------------
myaa <- rbind(aa_Ef, aa_PRIO)
protein.length(myaa)
## ----thermo_refs_table, eval=FALSE--------------------------------------------
# thermo.refs() ## shows table in a browser
## ----width180, include=FALSE------------------------------------------------------------------------------------------------------------------------------------------------------
options(width = 180)
## ----thermo_refs_numeric----------------------------------------------------------------------------------------------------------------------------------------------------------
iATP <- info("ATP-4")
iMgATP <- info("MgATP-2")
thermo.refs(c(iATP, iMgATP))
## ----thermo_refs_character--------------------------------------------------------------------------------------------------------------------------------------------------------
thermo.refs(c("HDNB78", "MGN03"))
## ----thermo_refs_subcrt, message=FALSE--------------------------------------------------------------------------------------------------------------------------------------------
sres <- subcrt(c("C2H5OH", "O2", "CO2", "H2O"), c(-1, -3, 2, 3))
thermo.refs(sres)
## ----width80, include=FALSE---------------------------------------------------
options(width = 80)
## ----thermo_refs_browse, eval=FALSE-------------------------------------------
# iFo <- info("forsterite")
# ref <- thermo.refs(iFo)
# browseURL(ref$URL) ## opens a link to worldcat.org
## ----info_S3, results="hide"--------------------------------------------------
info(info("S3-"))
## ----check_OBIGT--------------------------------------------------------------
file <- system.file("extdata/adds/OBIGT_check.csv", package = "CHNOSZ")
dat <- read.csv(file, as.is = TRUE)
nrow(dat)
## ----citation_CHNOSZ, results="asis", echo = FALSE----------------------------
cref <- citation("CHNOSZ")
print(cref[1], style = "html")
## ----citation_multimetal, results="asis", echo = FALSE------------------------
print(cref[2], style = "html")
## ----maintainer_CHNOSZ--------------------------------------------------------
maintainer("CHNOSZ")
## ----file_edit_anintro, eval=FALSE--------------------------------------------
# file.edit(system.file("doc/anintro.Rmd", package = "CHNOSZ"))
## ----the_end------------------------------------------------------------------
###### ## ## ## ## ###### ##### #####
## ##===## ## \\## ## ## \\ //
###### ## ## ## ## ###### ##### #####
## ----sessionInfo--------------------------------------------------------------
sessionInfo()
Try the CHNOSZ package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.