[CHNOSZ-commits] r665 - in pkg/CHNOSZ: . R inst man tests/testthat

Wed Apr 7 08:14:08 CEST 2021

Author: jedick
Date: 2021-04-07 08:14:08 +0200 (Wed, 07 Apr 2021)
New Revision: 665

Modified:
   pkg/CHNOSZ/DESCRIPTION
   pkg/CHNOSZ/R/ionize.aa.R
   pkg/CHNOSZ/inst/NEWS.Rd
   pkg/CHNOSZ/inst/TODO
   pkg/CHNOSZ/man/subcrt.Rd
   pkg/CHNOSZ/tests/testthat/test-ionize.aa.R
Log:
Fix bug where ionize.aa() was affected by E.units()


Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================

--- pkg/CHNOSZ/DESCRIPTION	2021-04-07 04:19:46 UTC (rev 664)
+++ pkg/CHNOSZ/DESCRIPTION	2021-04-07 06:14:08 UTC (rev 665)
@@ -1,6 +1,6 @@
 Date: 2021-04-07
 Package: CHNOSZ
-Version: 1.4.0-34
+Version: 1.4.0-35
 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/R/ionize.aa.R
===================================================================
--- pkg/CHNOSZ/R/ionize.aa.R	2021-04-07 04:19:46 UTC (rev 664)
+++ pkg/CHNOSZ/R/ionize.aa.R	2021-04-07 06:14:08 UTC (rev 665)
@@ -28,7 +28,10 @@
   # what property are we after
   sprop <- c("G", property)
   if(property %in%  c("A", "Z")) sprop <- "G"
-  sout <- subcrt(c(ineutral, icharged), T=T[!dupPT], P=P[!dupPT], property=sprop)$out
+  # Use convert=FALSE so we get results in calories 20210407
+  # (means we need to supply temperature in Kelvin)
+  TK <- convert(T, "K")
+  sout <- subcrt(c(ineutral, icharged), T=TK[!dupPT], P=P[!dupPT], property=sprop, convert = FALSE)$out
   # the G-values
   Gs <- sapply(sout, function(x) x$G)
   # keep it as a matrix even if we have only one unique T, P-combo
@@ -40,7 +43,7 @@
   DG <- t(sapply(pTP, function(x) DG[match(x, uPT), , drop=FALSE]))
   # the pK values (-logK) 
   DG <- DG * charges
-  pK <- apply(DG, 2, function(x) convert(x, "logK", T=convert(T, "K")))
+  pK <- apply(DG, 2, function(x) convert(x, "logK", T=TK))
   # keep it as a matrix even if we have only one T, P-combo
   if(lmax==1) pK <- t(pK)
   if(identical(ret.val, "pK")) {

Modified: pkg/CHNOSZ/inst/NEWS.Rd
===================================================================
--- pkg/CHNOSZ/inst/NEWS.Rd	2021-04-07 04:19:46 UTC (rev 664)
+++ pkg/CHNOSZ/inst/NEWS.Rd	2021-04-07 06:14:08 UTC (rev 665)
@@ -10,7 +10,7 @@
 \newcommand{\s}{\ifelse{latex}{\eqn{_{#1}}}{\ifelse{html}{\out{<sub>#1</sub>}}{#1}}}
 \newcommand{\S}{\ifelse{latex}{\eqn{^{#1}}}{\ifelse{html}{\out{<sup>#1</sup>}}{^#1}}}
 
-\section{Changes in CHNOSZ version 1.4.0-31 (2021-03-22)}{
+\section{Changes in CHNOSZ version 1.4.0-35 (2021-04-07)}{
 
   \subsection{NEW FEATURES}{
     \itemize{
@@ -95,6 +95,9 @@
       \item Add option \code{thermo()$opt$ionize.aa} to control calculation of
       properties of ionized proteins in \code{affinity()}.
 
+      \item Fix bug where protein ionization calculations were incorrected
+      affected by the setting of \code{E.units()}.
+
     }
   }
 

Modified: pkg/CHNOSZ/inst/TODO
===================================================================
--- pkg/CHNOSZ/inst/TODO	2021-04-07 04:19:46 UTC (rev 664)
+++ pkg/CHNOSZ/inst/TODO	2021-04-07 06:14:08 UTC (rev 665)
@@ -28,15 +28,6 @@
 
 - make swap.basis() work when DEW model is active
 
-- Check that protein ionization calculations are not affected by E.units()
-library(testthat)
-reset()
-basis(c("CO2", "H2", "NH4+", "H2O", "H2S", "H+"))
-a1 <- affinity(iprotein = 1)
-E.units("J")
-a2 <- affinity(iprotein = 1)
-expect_equal(a1$values, a2$values)
-
 [20210212]
 
 - Fix formatting of messages from palply() (seen with read.fasta())
@@ -43,8 +34,6 @@
 
 [20210219]
 
-- Don't use data(thermo) in demo Shiny app (citrate)
-
 - Check all mineral abbreviations in Berman_cr.csv.
 
 [20210301]

Modified: pkg/CHNOSZ/man/subcrt.Rd
===================================================================
--- pkg/CHNOSZ/man/subcrt.Rd	2021-04-07 04:19:46 UTC (rev 664)
+++ pkg/CHNOSZ/man/subcrt.Rd	2021-04-07 06:14:08 UTC (rev 665)
@@ -25,7 +25,7 @@
   \item{exceed.Ttr}{logical, calculate Gibbs energies of mineral phases and other species beyond their transition temperatures?}
   \item{exceed.rhomin}{logical, return properties of species in the HKF model below 0.35 g cm\S{-3}?}
   \item{logact}{numeric, logarithms of activities of species in reaction}
-  \item{convert}{logical, are input and output units of T and P those of the user (\code{TRUE}) (see \code{\link{T.units}}), or are they Kelvin and bar (\code{FALSE})?}
+  \item{convert}{logical, are units of T, P, and energy settable by the user (default) (see \code{\link{T.units}})?}
   \item{autobalance}{logical, attempt to automatically balance reaction with basis species?}
   \item{IS}{numeric, ionic strength(s) at which to calculate adjusted molal properties, mol kg\eqn{^{-1}}{^-1}}
 }
@@ -55,8 +55,8 @@
 On the other hand, if the \code{water} option is \samp{SUPCRT} (the default), \code{E} and \code{kT} can be calculated for water but not for aqueous species.
 (This is not an inherent limitation in either formulation, but it is just a matter of implementation.)
 
-Depending on the units currently defined (\code{\link{E.units}}) the values of \code{G}, \code{H}, \code{S} and \code{Cp} are returned using calories or Joules as the unit of energy, but only if \code{convert} is \code{TRUE}.
-Likewise, the input values of \code{T} and \code{P} are interpreted to have the units specified through \code{\link{T.units}} and \code{\link{P.units}}, but setting \code{convert} to \code{FALSE} forces \code{subcrt} to treat them as Kelvin and bar, respectively.
+If \code{convert} is \code{TRUE} (the default), the input values of \code{T} and \code{P} are interpreted to have the units given by \code{\link{T.units}} and \code{\link{P.units}} (default: \degC and bar), and the output values of \code{G}, \code{H}, \code{S} and \code{Cp} are based on the units given in \code{\link{E.units}} (default: calories).
+If \code{convert} is \code{FALSE}, the user units (\code{\link{T.units}}, \code{\link{P.units}}, and \code{\link{E.units}}) are ignored, and \code{T} and \code{P} are taken to be in Kelvin and bar, and the returned values of \code{G}, \code{H}, \code{S} and \code{Cp} are in calories.
 
 A chemical reaction is defined if \code{coeff} is given.
 In this mode the standard molal properties of species are summed according to the stoichiometric \code{coeff}icients, where negative values denote reactants.

Modified: pkg/CHNOSZ/tests/testthat/test-ionize.aa.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-ionize.aa.R	2021-04-07 04:19:46 UTC (rev 664)
+++ pkg/CHNOSZ/tests/testthat/test-ionize.aa.R	2021-04-07 06:14:08 UTC (rev 665)
@@ -104,6 +104,17 @@
   expect_equal(ionize.aa(aa, property="A")[1, ], A.protein, check.attributes=FALSE)
 })
 
+# Test added 20210407
+test_that("protein ionization calculations are not affected by E.units()", {
+  E.units("cal")
+  basis(c("CO2", "H2", "NH4+", "H2O", "H2S", "H+"))
+  a1 <- affinity(iprotein = 1)
+  E.units("J")
+  a2 <- affinity(iprotein = 1)
+  expect_equal(a1$values, a2$values)
+  reset()
+})
+
 # references
 
 # Dick, J. M., LaRowe, D. E. and Helgeson, H. C. (2006) 

