[CHNOSZ-commits] r592 - in pkg/CHNOSZ: . inst man tests/testthat vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jul 28 14:06:42 CEST 2020
Author: jedick
Date: 2020-07-28 14:06:41 +0200 (Tue, 28 Jul 2020)
New Revision: 592
Added:
pkg/CHNOSZ/tests/testthat/test-objective.R
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/inst/NEWS.Rd
pkg/CHNOSZ/man/objective.Rd
pkg/CHNOSZ/man/subcrt.Rd
pkg/CHNOSZ/man/swap.basis.Rd
pkg/CHNOSZ/man/util.array.Rd
pkg/CHNOSZ/man/util.formula.Rd
pkg/CHNOSZ/man/util.misc.Rd
pkg/CHNOSZ/tests/testthat/test-swap.basis.R
pkg/CHNOSZ/vignettes/multi-metal.Rmd
Log:
Remove stopifnot() from examples, round 1
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2020-07-28 08:39:16 UTC (rev 591)
+++ pkg/CHNOSZ/DESCRIPTION 2020-07-28 12:06:41 UTC (rev 592)
@@ -1,6 +1,6 @@
Date: 2020-07-28
Package: CHNOSZ
-Version: 1.3.6-65
+Version: 1.3.6-66
Title: Thermodynamic Calculations and Diagrams for Geochemistry
Authors at R: c(
person("Jeffrey", "Dick", , "j3ffdick at gmail.com", role = c("aut", "cre"),
Modified: pkg/CHNOSZ/inst/NEWS.Rd
===================================================================
--- pkg/CHNOSZ/inst/NEWS.Rd 2020-07-28 08:39:16 UTC (rev 591)
+++ pkg/CHNOSZ/inst/NEWS.Rd 2020-07-28 12:06:41 UTC (rev 592)
@@ -241,6 +241,11 @@
\item TODO: comment out data modifications in demos/mosaic.R
+ \item TODO: move H2O92D.f and R wrapper to a separate package (so people
+ don't have to compile anything to install CHNOSZ updates).
+
+ \item TODO: get retrieve() into more examples/demos
+
}
}
Modified: pkg/CHNOSZ/man/objective.Rd
===================================================================
--- pkg/CHNOSZ/man/objective.Rd 2020-07-28 08:39:16 UTC (rev 591)
+++ pkg/CHNOSZ/man/objective.Rd 2020-07-28 12:06:41 UTC (rev 592)
@@ -78,29 +78,16 @@
\examples{
\dontshow{reset()}
-## a made-up system: 4 species, 1 condition
-loga1 <- t(-4:-1)
-loga2 <- loga1 + 1
-stopifnot(qqr(loga1) < 1)
-stopifnot(RMSD(loga1, loga1) == 0)
-stopifnot(RMSD(loga1, loga2) == 1)
-stopifnot(CVRMSD(loga1, loga2) == -0.4) # 1 / mean(-4:-1)
-stopifnot(spearman(loga1, loga2) == 1)
-stopifnot(spearman(loga1, rev(loga2)) == -1)
-# less statistical, more thermodynamical...
-stopifnot(all.equal(DGmix(loga1), -0.1234)) # as expected for decimal logarithms
-stopifnot(all.equal(DDGmix(loga1, loga2), 0.0004))
-
-## transforming an equilibrium assemblage of n-alkanes
+## Transforming an equilibrium assemblage of n-alkanes
basis(c("methane", "hydrogen"))
species(c("methane", "ethane", "propane", "butane"), "liq")
-# calculate equilibrium assemblages over a range of logaH2
-a <- affinity(H2=c(-10, -5, 101), exceed.Ttr=TRUE)
+# Calculate equilibrium assemblages over a range of logaH2
+a <- affinity(H2 = c(-10, -5, 101), exceed.Ttr = TRUE)
e <- equilibrate(a)
-# take a reference equilibrium distribution at logfH2 = -7.5
+# Take a reference equilibrium distribution at logfH2 = -7.5
loga1 <- list2array(e$loga.equil)[51, ]
Astar <- list2array(e$Astar)[51, ]
-# equilibrium at any other logfH2 is not equilibrium at logfH2 = -7.5
+# Equilibrium at any other logfH2 is not equilibrium at logfH2 = -7.5
DGtr.out <- DDGmix.out <- numeric()
for(i in 1:length(a$vals[[1]])) {
loga2 <- list2array(e$loga.equil)[i, ]
@@ -107,16 +94,12 @@
DGtr.out <- c(DGtr.out, DGtr(t(loga1), loga2, t(Astar)))
DDGmix.out <- c(DDGmix.out, DDGmix(t(loga1), loga2))
}
-# all(DGtr >= 0) is TRUE
-stopifnot(all(DGtr.out >= 0))
-# all(DDGmix >= 0) is FALSE
-stopifnot(!all(DDGmix.out >= 0))
-# a plot is also nice
-thermo.plot.new(xlim=range(a$vals[[1]]), xlab=axis.label("H2"),
- ylim=range(DDGmix.out, DGtr.out), ylab="energy")
-abline(h=0, lty=2)
-abline(v=-7.5, lty=2)
-text(-7.6, 2, "reference condition", srt=90)
+# Plot the results
+thermo.plot.new(xlim = range(a$vals[[1]]), xlab = axis.label("H2"),
+ ylim = range(DDGmix.out, DGtr.out), ylab = "energy")
+abline(h = 0, lty = 2)
+abline(v = -7.5, lty = 2)
+text(-7.6, 2, "reference condition", srt = 90)
lines(a$vals[[1]], DDGmix.out)
lines(a$vals[[1]], DGtr.out)
text(-6, 5.5, expr.property("DDGmix/2.303RT"))
@@ -123,8 +106,9 @@
text(-6, 2.3, expr.property("DGtr/2.303RT"))
title(main=paste("Transformation between metastable equilibrium\n",
"assemblages of n-alkanes"))
-# take-home message: use DGtr to measure distance from equilibrium in
-# open-system transformations (constant T, P, chemical potentials of basis species)
+# DGtr but not DDGmix is positive (or zero) everywhere
+stopifnot(all(DGtr.out >= 0))
+stopifnot(!all(DDGmix.out >= 0))
}
\references{
Modified: pkg/CHNOSZ/man/subcrt.Rd
===================================================================
--- pkg/CHNOSZ/man/subcrt.Rd 2020-07-28 08:39:16 UTC (rev 591)
+++ pkg/CHNOSZ/man/subcrt.Rd 2020-07-28 12:06:41 UTC (rev 592)
@@ -118,61 +118,45 @@
\examples{
\dontshow{reset()}
-## properties of species
+## Properties of species
subcrt("water")
-# calculating at different temperatures
-subcrt("water", T=seq(0, 100, 10))
-# calculating at even increments
-subcrt("water", T=seq(500, 1000, length.out=10),
- P=seq(5000, 10000, length.out=10))
-# calculating on a temperature-pressure grid
-subcrt("water", T=c(500, 1000), P=c(5000, 10000), grid="P")
-# to calculate only selected properties
-subcrt("water", property=c("G", "E"))
-# the properties of multiple species
-subcrt(c("glucose", "ethanol", "CO2"))
+# Change temperature
+subcrt("water", T = seq(0, 100, 20))
+# Change temperature and pressure
+T <- seq(500, 1000, 100)
+P <- seq(5000, 10000, 1000)
+subcrt("water", T = T, P = P)
+# Temperature-pressure grid
+subcrt("water", T = c(500, 1000), P = c(5000, 10000), grid = "P")
-## properties of reactions
-subcrt(c("H2O", "H+", "K-feldspar", "kaolinite", "K+", "SiO2"),
- c(-1, -2, -2, 1, 2, 4))
-subcrt(c("glucose", "ethanol", "CO2"), c(-1, 2, 2))
-# to specify the states
-subcrt(c("glucose", "ethanol", "CO2"), c(-1, 2, 2), c("aq", "aq", "gas"))
+## Properties of reactions
+subcrt(c("glucose", "ethanol", "CO2"), c(-1, 2, 2), T = 25)
+# Use CO2(gas) (or just change "CO2" to "carbon dioxide")
+subcrt(c("glucose", "ethanol", "CO2"), c(-1, 2, 2), c("aq", "aq", "gas"), T = 25)
-## auto balancing reactions
-# the basis species must first be defined
+## Automatically balance reactions
+# First define the basis species
basis(c("CO2", "H2O", "NH3", "H2S", "O2"))
-subcrt(c("glucose", "ethanol"), c(-1, 3))
-# a bug in CHNOSZ <0.9 caused the following
-# to initiate an infinite loop
+# Auto-balance adds the required amount of H2O and O2
+subcrt(c("ethanol", "glucose"), c(-3, 1), T = 37)
+# An example with H+
basis(c("H2O", "H2S", "O2", "H+"))
-subcrt(c("HS-", "SO4-2"), c(-1, 1))
-# because O2,aq is in the basis, this is a non-reaction
-# (O2,aq to O2,aq)
-subcrt("O2", 1, "aq")
-# but this one auto-balances into a reaction
-# (O2,aq to O2,gas)
-subcrt("O2", 1, "gas")
-# properties of a species and a formation
-# reaction for that species
-subcrt("C2H5OH") # species
-basis("CHNOS")
-subcrt("C2H5OH", 1) # reaction
+subcrt(c("HS-", "SO4-2"), c(-1, 1), T = 100)
-## mineral polymorphs
-# properties of the stable polymorph
+## Mineral polymorphs
+# Properties of the stable polymorph
subcrt("pyrrhotite")
-# properties of just the high-T phase
-subcrt(c("pyrrhotite"), state="cr2")
-# polymorphic transitions in a reaction
+# Properties of one of the polymorphs (metastable at other temperatures)
+subcrt(c("pyrrhotite"), state = "cr2")
+# Reactions automatically use stable polymorph
subcrt(c("pyrite", "pyrrhotite", "H2O", "H2S", "O2"), c(-1, 1, -1, 1, 0.5))
-## these produce messages about problems with the calculation
-# above the T, P limits for the H2O equations of state
-subcrt("alanine", T=c(2250, 2251), P=c(30000, 30001), grid="T")
-# Psat is not defined above the critical point
-## this also gives a warning
-suppressWarnings(subcrt("alanine", T=seq(0, 5000, by=1000)))
+## Messages about problems with the calculation
+# Above the T, P limits for the H2O equations of state
+subcrt("alanine", T = c(2250, 2251), P = c(30000, 30001), grid = "T")
+# Psat is not defined above the critical point;
+# this also gives a warning that is suppressed to keep the examples clean
+suppressWarnings(subcrt("alanine", T = seq(0, 5000, by = 1000)))
## minerals with phase transitions
# compare calculated values of heat capacity of iron with
@@ -194,69 +178,6 @@
T.units("C")
E.units("cal")
-## Skarn example after Johnson et al., 1992
-P <- seq(500, 5000, 500)
-# this is like the temperature specification used
-# in the example by Johnson et al., 1992
-T <- seq(0, 1000, 100)
-s <- suppressWarnings(subcrt(
- c("andradite", "carbon dioxide", "H2S", "Cu+", "quartz",
- "calcite", "chalcopyrite", "H+", "H2O"),
- coeff=c(-1, -3, -4, -2, 3, 3, 2, 2, 3),
- state=c("cr", "g", "aq", "aq", "cr", "cr", "cr", "aq", "liq"),
- P=P, T=T, grid="P"))
-# The results are not identical to SUPCRT92, as CHNOSZ has updated
-# parameters for species e.g. Cu+ from Shock et al., 1997.
-# Check the calculated phase transitions for chalcopyrite
-stopifnot(all.equal(s$polymorphs$chalcopyrite[1:11],
- c(1, 1, 1, 1, 1, 1, 2, 3, 3, 3, 3)))
-
-## Standard Gibbs energy of reactions with HCN and
-## formaldehyde, after Schulte and Shock, 1995 Fig. 1
-rxn1 <- subcrt(c("formaldehyde","HCN","H2O","glycolic acid","NH3"),
- c(-1,-1,-2,1,1),P=300)
-rxn2 <- subcrt(c("formaldehyde","HCN","H2O","glycine"),
- c(-1,-1,-1,1),P=300)
-plot(x=rxn1$out$T,rxn1$out$G/1000,type="l",ylim=c(-40,-10),
- xlab=axis.label("T"),ylab=axis.label("DG0r","k"))
-lines(rxn1$out$T,rxn2$out$G/1000)
-# write the reactions on the plot
-text(150, -14, describe.reaction(rxn1$reaction, iname=c(1,2,4)))
-text(200, -35, describe.reaction(rxn2$reaction, iname=c(1,2)))
-title(main=paste("Standard Gibbs energy of reactions",
- "after Schulte and Shock, 1995",sep="\n"))
-
-## Calculation of chemical affinities
-# after LaRowe and Helgeson, 2007, Fig. 3 (a): reduction of nicotinamide
-# adenine dinucleotide (NAD) coupled to oxidation of glucose
-# list the available NAD species
-info("NAD ")
-T <- seq(0, 120, 10)
-# oxidation of glucose (C6H12O6)
-basis(c("glucose", "H2O", "NH3", "CO2", "H+"), c(-3, 0, 999, -3, -7))
-s <- subcrt(c("NAD(ox)-", "NAD(red)-2"), c(-12, 12), logact=c(0, 0), T=T)
-# LH07's diagrams are shown per mole of electron (24 e- per 12 NAD)
-A <- s$out$A/24/1000
-plot(x=T, y=A, xlim=range(T), ylim=c(1.4, 5.4),
- xlab=axis.label("T"), ylab=axis.label("A", prefix="k"), type="l")
-text("NAD(ox)-/NAD(red)-2 = 1", x=53, y=median(A), srt=21)
-# different activity ratio
-s <- subcrt(c("NAD(ox)-","NAD(red)-2"), c(-12, 12), logact=c(1, 0), T=T)
-A <- s$out$A/24/1000
-lines(x=T, y=A)
-text("NAD(ox)-/NAD(red)-2 = 10", x=55, y=median(A), srt=24)
-# different activity ratio
-s <- subcrt(c("NAD(ox)-","NAD(red)-2"), c(-12, 12), logact=c(0, 1), T=T)
-A <- s$out$A/24/1000
-lines(x=T, y=A)
-text("NAD(ox)-/NAD(red)-2 = 0.1", x=52, y=median(A), srt=18)
-# print the reaction and chemical conditions on the plot
-text(0, 5.3, describe.reaction(s$reaction, iname=c(1, 2)), adj=0)
-text(0, 5.1, describe.basis(oneline=TRUE, ibasis=c(1, 2, 4, 5)), adj=0)
-# label the plot
-title(main=paste("Reduction of NAD coupled to oxidation of glucose",
- "after LaRowe and Helgeson, 2007", sep="\n"))
-
## Subzero (degrees C) calculations
# uncomment the following to try IAPWS95 instead of SUPCRT92
#water("IAPWS95")
@@ -265,7 +186,7 @@
sb <- subcrt(c("H2O", "Na+"), T=seq(-30, 10), P=1)$out
# start plot with extra room on right
opar <- par(mar=c(5, 4, 4, 4))
-# plot G
+# plot Delta G
plot(sb$water$T, sb$water$G, ylim=c(-63000, -56000), xlab=axis.label("T"),
ylab=axis.label("DG0"))
points(sb$`Na+`$T, sb$`Na+`$G, pch=2)
@@ -290,21 +211,16 @@
# at high temperatures (also, the gas is metastable at low
# temperatures, but that doesn't produce NA in the output)
sout2 <- subcrt(rep("octane", 2), c("liq", "gas"),
- c(-1, 1), T=seq(-50, 300, 0.1), P=2, exceed.Ttr=TRUE)$out
+ c(-1, 1), T = seq(-50, 300, 0.1), P = 2, exceed.Ttr = TRUE)$out
sout20 <- subcrt(rep("octane", 2), c("liq", "gas"),
- c(-1, 1), T=seq(-50, 300, 0.1), P=20, exceed.Ttr=TRUE)$out
+ c(-1, 1), T = seq(-50, 300, 0.1), P = 20, exceed.Ttr = TRUE)$out
# find T with the Gibbs energy of reaction that is closest to zero
Tvap2 <- sout2$T[which.min(abs(sout2$G))]
Tvap20 <- sout20$T[which.min(abs(sout20$G))]
-# the boiling point increases with pressure
-stopifnot(Tvap20 > Tvap2)
-# more precisely, the calculated boiling points should be near the
-# empirical values (digitized from Fig. 1 of Helgeson et al., 1998)
-Tvap_2bar <- 156
-Tvap_20bar <- 276
-stopifnot(abs(Tvap2 - Tvap_2bar) < 6)
-stopifnot(abs(Tvap20 - Tvap_20bar) < 25)
-# those comparisons would fail if varP were FALSE (the default)
+# Compare these with experimental values (Fig. 1 of Helgeson et al., 1998)
+Tvap2.exp <- 156
+Tvap20.exp <- 276
+# Reset varP to FALSE (the default)
thermo("opt$varP" = FALSE)
}
@@ -317,8 +233,6 @@
LaRowe, D. E. and Helgeson, H. C. (2007) Quantifying the energetics of metabolic reactions in diverse biogeochemical systems: electron flow and ATP synthesis. \emph{Geobiology} \bold{5}, 153--168. \url{https://doi.org/10.1111/j.1472-4669.2007.00099.x}
- Schulte, M. D. and Shock, E. L. (1995) Thermodynamics of Strecker synthesis in hydrothermal systems. \emph{Orig. Life Evol. Biosph.} \bold{25}, 161--173. \url{https://doi.org/10.1007/BF01581580}
-
Shock, E. L., Oelkers, E. H., Johnson, J. W., Sverjensky, D. A. and Helgeson, H. C. (1992) Calculation of the thermodynamic properties of aqueous species at high pressures and temperatures: Effective electrostatic radii, dissociation constants and standard partial molal properties to 1000 \degC and 5 kbar. \emph{J. Chem. Soc. Faraday Trans.} \bold{88}, 803--826. \url{https://doi.org/10.1039/FT9928800803}
}
Modified: pkg/CHNOSZ/man/swap.basis.Rd
===================================================================
--- pkg/CHNOSZ/man/swap.basis.Rd 2020-07-28 08:39:16 UTC (rev 591)
+++ pkg/CHNOSZ/man/swap.basis.Rd 2020-07-28 12:06:41 UTC (rev 592)
@@ -39,54 +39,22 @@
\examples{
\dontshow{reset()}
-## swapping basis species
+## Swapping basis species
# start with a preset basis definition
b1 <- basis("CHNOS+")
# swap H2(aq) for O2(gas)
-(b2 <- swap.basis("O2", "H2"))
-# the logarithm of activity calculated for H2
-# is equal to the one calculated from the equilibrium constant
-# for H2O = H2 + 0.5O2
-logK <- subcrt(c("oxygen","H2","H2O"), c(-0.5,-1,1), T=25)$out$logK
-# the equilibrium value of logaH2
-# (for logaH2O = 0 and logfO2 = -80)
-(logaH2 <- -logK + 40)
-stopifnot(all.equal(logaH2, b2$logact[5]))
-# put O2 back in
+b2 <- swap.basis("O2", "H2")
+# put oxygen back in
b3 <- swap.basis("H2", "oxygen")
-# we have returned to starting point
-stopifnot(all.equal(b1$logact, b3$logact))
-# demonstrating the interconvertibility between
-# chemical potentials of elements and logarithms
-# of activities of basis species at high temperature
+# Interconversion of chemical potentials of elements and
+# logarithms of activities of basis species at high temperature
basis("CHNOS+")
bl1 <- basis()$logact
-emu <- element.mu(T=100)
-bl2 <- basis.logact(emu, T=100)
-# note that basis.logact produces a named array
-stopifnot(all.equal(bl1, as.numeric(bl2)))
-
-## swapping basis species while species are defined
-## and using numeric species indices
-basis("MgCHNOPS+")
-# load some Mg-ATP species
-species(c("MgATP-2", "MgHATP-", "MgH2ATP", "Mg2ATP"))
-# swap in CO2(g) for CO2(aq)
-swap.basis("CO2", "carbon dioxide")
-a1 <- affinity()
-# swap in CH4(g) for CO2(g)
-swap.basis("carbon dioxide", "methane")
-a2 <- affinity()
-# the equilibrium fugacity of CH4 is *very* low
-# swap in CO2(aq) for CH4(g)
-swap.basis("methane", "CO2")
-a3 <- affinity()
-# swapping the basis species didn't affect the affinities
-# of the formation reactions of the species, since
-# the chemical potentials of the elements were unchanged
-stopifnot(all.equal(a1$values, a2$values))
-stopifnot(all.equal(a1$values, a3$values))
+emu <- element.mu(T = 100)
+bl2 <- basis.logact(emu, T = 100)
+# There's no difference
+round(bl2 - bl1, 10)
}
\concept{Extended workflow}
Modified: pkg/CHNOSZ/man/util.array.Rd
===================================================================
--- pkg/CHNOSZ/man/util.array.Rd 2020-07-28 08:39:16 UTC (rev 591)
+++ pkg/CHNOSZ/man/util.array.Rd 2020-07-28 12:06:41 UTC (rev 592)
@@ -38,53 +38,50 @@
\examples{
# start with a matrix
-x <- matrix(1:12,ncol=3)
+mat <- matrix(1:12, ncol = 3)
# pay attention to the following when
# writing examples that test for identity!
-identical(1*x,x) # FALSE
+identical(1 * mat, mat) # FALSE
# create two matrices that are multiples of the first
-a <- 1*x
-b <- 2*a
+a <- 1 * mat
+b <- 2 * mat
# these both have two dimensions of lengths 4 and 3
dim(a) # 4 3
# combine them to make an array with three dimensions
-c <- list2array(list(a,b))
+x <- list2array(list(a, b))
# the third dimension has length 2
-dim(c) # 4 3 2
-# the first slice of the third dimension == a
-stopifnot(identical( slice(c,3), a ))
-# the second slice of the third dimension == b
-stopifnot(identical( slice(c,3,2), b ))
+dim(x) # 4 3 2
+# the first slice of the third dimension
+slice(x, 3) # a
+# the second slice of the third dimension
+slice(x, 3, 2) # b
# 'slice' works just like the bracket operator
-c11 <- slice(c,1)
-c12 <- slice(c,1,2)
-c21 <- slice(c,2,1)
-c212 <- slice(c,2,1:2)
-stopifnot(identical( c11, c[1,,] ))
-stopifnot(identical( c12, c[2,,] ))
-stopifnot(identical( c21, c[,1,] ))
-stopifnot(identical( c212, c[,1:2,] ))
+slice(x, 1) # x[1, , ]
+slice(x, 1, 2) # x[2, , ]
+slice(x, 2, 1) # x[, 1, ]
+slice(x, 2, 1:2) # x[, 1:2, ]
+
# let us replace part of the array
-d <- slice(c,3,2,value=a)
+y <- slice(x, 3, 2, value = a)
# now the second slice of the third dimension == a
-stopifnot(identical( slice(d,3,2), a ))
+slice(y, 3, 2) # a
# and the sum across the third dimension == b
-stopifnot(identical( dimSums(d,3), b ))
+dimSums(y, 3) # b
# taking the sum removes that dimension
-dim(d) # 4 3 2
-dim(dimSums(d,1)) # 3 2
-dim(dimSums(d,2)) # 4 2
-dim(dimSums(d,3)) # 4 3
+dim(y) # 4 3 2
+dim(dimSums(y, 1)) # 3 2
+dim(dimSums(y, 2)) # 4 2
+dim(dimSums(y, 3)) # 4 3
# working with an 'affinity' object
\dontshow{reset()}
basis("CHNOS+")
species("alanine")
-a1 <- affinity(O2=c(-80,-60)) # at pH=7
-a2 <- affinity(O2=c(-80,-60),pH=c(0,14,7))
+a1 <- affinity(O2 = c(-80, -60)) # i.e. pH=7
+a2 <- affinity(O2 = c(-80, -60), pH = c(0, 14, 7))
# in the 2nd dimension (pH) get the 4th slice (pH=7)
-a3 <- slice.affinity(a2,2,4)
-stopifnot(all.equal(a1$values,a3$values))
+a3 <- slice.affinity(a2, 2, 4)
+stopifnot(all.equal(a1$values, a3$values))
}
\concept{Utility functions}
Modified: pkg/CHNOSZ/man/util.formula.Rd
===================================================================
--- pkg/CHNOSZ/man/util.formula.Rd 2020-07-28 08:39:16 UTC (rev 591)
+++ pkg/CHNOSZ/man/util.formula.Rd 2020-07-28 12:06:41 UTC (rev 592)
@@ -93,13 +93,12 @@
# that not quite equal to the value from water("G");
# probably using different entropies of the elements
-## average oxidation states of carbon
-stopifnot(ZC("CO2") == 4)
-stopifnot(ZC("CH4") == -4)
-stopifnot(ZC("CHNOSZ") == 7)
-si <- info(info("LYSC_CHICK"))
-stopifnot(si$formula == "C613H959N193O185S10")
-stopifnot(all.equal(ZC(si$formula), 0.0163132137031))
+## Average oxidation states of carbon
+ZC(c("CO2", "CH4", "CHNOSZ")) # 4, -4, 7
+si <- info("LYSC_CHICK")
+# Can use species index or formula
+ZC(si)
+ZC(info(si)$formula)
## calculate the chemical formulas, then
## ZC of all of the proteins in CHNOSZ' database
Modified: pkg/CHNOSZ/man/util.misc.Rd
===================================================================
--- pkg/CHNOSZ/man/util.misc.Rd 2020-07-28 08:39:16 UTC (rev 591)
+++ pkg/CHNOSZ/man/util.misc.Rd 2020-07-28 12:06:41 UTC (rev 592)
@@ -75,27 +75,26 @@
diff(substuff$out$iron$H)
})
-## scale logarithms of activity
+## Scale logarithms of activity
# suppose we have two proteins whose lengths are 100 and
# 200; what are the logarithms of activity of the proteins
# that are equal to each other and that give a total
# activity of residues equal to unity?
-logact <- c(-3,-3) # could be any two equal numbers
-length <- c(100,200)
+logact <- c(-3, -3) # could be any two equal numbers
+length <- c(100, 200)
logact.tot <- 0
-loga <- unitize(logact,length,logact.tot)
-# the proteins have equal activity
-stopifnot(identical(loga[1],loga[2]))
-# the sum of activity of the residues is unity
-stopifnot(isTRUE(all.equal(sum(10^loga * length),1)))
-## now, what if the activity of protein 2 is ten
-## times that of protein 1?
-logact <- c(-3,-2)
-loga <- unitize(logact,length,logact.tot)
-# the proteins have unequal activity
-stopifnot(isTRUE(all.equal(loga[2]-loga[1],1)))
+loga <- unitize(logact, length, logact.tot)
+# The proteins have equal activity
+loga[1] == loga[2]
+# The sum of activity of the residues is unity
+all.equal(sum(10^loga * length), 1)
+## What if the activity of protein 2 is ten times that of protein 1?
+logact <- c(-3, -2)
+loga <- unitize(logact, length, logact.tot)
+# The proteins have unequal activity,
# but the activities of residues still add up to one
-stopifnot(isTRUE(all.equal(sum(10^loga * length),1)))
+all.equal(loga[2] - loga[1], 1)
+all.equal(sum(10^loga * length), 1)
}
\concept{Thermodynamic calculations}
Added: pkg/CHNOSZ/tests/testthat/test-objective.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-objective.R (rev 0)
+++ pkg/CHNOSZ/tests/testthat/test-objective.R 2020-07-28 12:06:41 UTC (rev 592)
@@ -0,0 +1,16 @@
+context("objective")
+
+# 20200728 moved from objective.Rd
+test_that("calculations give expected results", {
+ loga1 <- t(-4:-1)
+ loga2 <- loga1 + 1
+ expect_lt(qqr(loga1), 1)
+ expect_equal(RMSD(loga1, loga1), 0)
+ expect_equal(RMSD(loga1, loga2), 1)
+ expect_equal(CVRMSD(loga1, loga2), -0.4) # 1 / mean(-4:-1)
+ expect_equal(spearman(loga1, loga2), 1)
+ expect_equal(spearman(loga1, rev(loga2)), -1)
+ # less statistical, more thermodynamical...
+ expect_equal(DGmix(loga1), -0.1234) # as expected for decimal logarithms
+ expect_equal(DDGmix(loga1, loga2), 0.0004)
+})
Modified: pkg/CHNOSZ/tests/testthat/test-swap.basis.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-swap.basis.R 2020-07-28 08:39:16 UTC (rev 591)
+++ pkg/CHNOSZ/tests/testthat/test-swap.basis.R 2020-07-28 12:06:41 UTC (rev 592)
@@ -29,3 +29,27 @@
# note: logact includes "PPM" for O2 (old) and H2 (new)
expect_identical(oldb$logact, newb$logact)
})
+
+# 20200728 moved from swap-basis.Rd
+test_that("swapping doesn't affect affinities of formation reactions of species", {
+ ## swapping basis species while species are defined
+ ## and using numeric species indices
+ basis("MgCHNOPS+")
+ # load some Mg-ATP species
+ species(c("MgATP-2", "MgHATP-", "MgH2ATP", "Mg2ATP"))
+ # swap in CO2(g) for CO2(aq)
+ swap.basis("CO2", "carbon dioxide")
+ a1 <- affinity()
+ # swap in CH4(g) for CO2(g)
+ swap.basis("carbon dioxide", "methane")
+ a2 <- affinity()
+ # the equilibrium fugacity of CH4 is *very* low
+ # swap in CO2(aq) for CH4(g)
+ swap.basis("methane", "CO2")
+ a3 <- affinity()
+ # swapping the basis species didn't affect the affinities
+ # of the formation reactions of the species, since
+ # the chemical potentials of the elements were unchanged
+ expect_equal(a1$values, a2$values)
+ expect_equal(a1$values, a3$values)
+})
Modified: pkg/CHNOSZ/vignettes/multi-metal.Rmd
===================================================================
--- pkg/CHNOSZ/vignettes/multi-metal.Rmd 2020-07-28 08:39:16 UTC (rev 591)
+++ pkg/CHNOSZ/vignettes/multi-metal.Rmd 2020-07-28 12:06:41 UTC (rev 592)
@@ -84,7 +84,7 @@
Flattening or simple overlay refers to independent calculations for two different systems that are displayed on the same diagram.
This example starts with a log*f*~O<sub>2</sub>~--pH base diagram for the C-O-H system then overlays a diagram for S-O-H.
-The second call to `affinity()` uses the argument recall feature, where the arguments are taken from the previous command.
+The second call to `affinity()` uses the argument recall feature, where the arguments after the first are taken from the previous command.
This allows calculations to be run at the same conditions for a different system.
This feature is also used in other examples in this vignette.
@@ -514,7 +514,7 @@
For this example, the speciation of aqueous sulfur is not considered; instead, the fugacity of S~2~ is a plotting variable.
The stable Fe-bearing minerals (Fe-S-O-H) are used as the changing basis species to make the diagram for Cu-bearing minerals (Cu-Fe-S-O-H).
Then, the two diagrams are flattened to show all minerals in a single diagram.
-Greener colors are used to indicate minerals with lower S~2~ and higher O~2~ in their formation reactions.
+Greener colors are used to indicate minerals with less S~2~ and more O~2~ in their formation reactions.
<button id="B-stack2" onclick="ToggleDiv('stack2')">Show code</button>
<div id="D-stack2" style="display: none">
More information about the CHNOSZ-commits
mailing list