Update value of gas constant (8.314463 J K-1 mol-1)

git-svn-id: svn://scm.r-forge.r-project.org/svnroot/chnosz/pkg/CHNOSZ@796 edb9625f-4e0d-4859-8d74-9fd3b1da38cb

DESCRIPTION

@@ -1,6 +1,6 @@
-Date: 2023-06-24
+Date: 2023-07-01
 Package: CHNOSZ
-Version: 2.0.0-15
+Version: 2.0.0-16
 Title: Thermodynamic Calculations and Diagrams for Geochemistry
 Authors@R: c(
     person("Jeffrey", "Dick", , "[email protected]", role = c("aut", "cre"),

R/AD.R

@@ -8,7 +8,8 @@ AD <- function(property = NULL, parameters = NULL, T = 298.15, P = 1, isPsat = T
   # Some constants (from Akinfiev and Diamond, 2004 doi:10.1016/j.fluid.2004.06.010)
   MW <- 18.0153 # g mol-1
   NW <- 1000 / MW # mol kg-1
-  R <- 8.31441 # J K-1 mol-1
+  #R <- 8.31441 # J K-1 mol-1  20190219
+  R <- 8.314463  # https://physics.nist.gov/cgi-bin/cuu/Value?r 20230630
   # R expressed in volume units
   RV <- 10 * R # cm3 bar K-1 mol-1
 

R/nonideal.R

@@ -34,8 +34,10 @@ nonideal <- function(species, speciesprops, IS, T, P, A_DH, B_DH, m_star = NULL,
     else stop("invalid method (", thermo$opt$nonideal, ")")
   }
 
-  #R <- 1.9872  # gas constant, cal K^-1 mol^-1
-  R <- 8.314445  # gas constant, J K^-1 mol^-1  20220325
+  # Gas constant
+  #R <- 1.9872  # cal K^-1 mol^-1  Value used in SUPCRT92
+  #R <- 8.314445  # = 1.9872 * 4.184 J K^-1 mol^-1  20220325
+  R <- 8.314463  # https://physics.nist.gov/cgi-bin/cuu/Value?r 20230630
 
   # Function to calculate extended Debye-Huckel equation and derivatives using Alberty's parameters
   Alberty <- function(prop = "loggamma", Z, I, T) {

R/protein.info.R

@@ -206,8 +206,11 @@ protein.equil <- function(protein, T = 25, loga.protein = 0, digits = 4) {
     G0protform <- G0protform + G0ionization
   }
   # Standard Gibbs energy of formation reaction of non/ionized residue equivalents, dimensionless
-  if(E_units == "cal") R <- 1.9872  # gas constant, cal K^-1 mol^-1
-  if(E_units == "J") R <- 8.314445  # gas constant, J K^-1 mol^-1  20220325
+  # Gas constant
+  #if(E_units == "cal") R <- 1.9872  # gas constant, cal K^-1 mol^-1
+  #if(E_units == "J") R <- 8.314445  # = 1.9872 * 4.184 J K^-1 mol^-1  20220325
+  if(E_units == "J") R <- 8.314463  # https://physics.nist.gov/cgi-bin/cuu/Value?r 20230630
+  if(E_units == "cal") R <- 8.314463 / 4.184
   G0res.RT <- G0protform/R/TK/plength
   mymessage("protein.equil [1]: per residue, reaction to form ", iword, " protein from basis species has G0/RT of ", signif(G0res.RT[1], digits))
   # Coefficients of basis species in formation reactions of residues

R/util.units.R

@@ -100,8 +100,10 @@ convert <- function(value, units, T = 298.15,
     if(units == 'cal') value <- value / Jcal
   }
   else if(units %in% c('g', 'logk')) {
-    #R <- 1.9872  # Gas constant, cal K^-1 mol^-1
-    R <- 8.314445  # Gas constant, J K^-1 mol^-1  20220325
+    # Gas constant
+    #R <- 1.9872  # cal K^-1 mol^-1  Value used in SUPCRT92
+    #R <- 8.314445  # = 1.9872 * 4.184 J K^-1 mol^-1  20220325
+    R <- 8.314463  # https://physics.nist.gov/cgi-bin/cuu/Value?r 20230630
     if(units == 'logk') value <- value / (-log(10) * R * T)
     if(units == 'g') value <- value * (-log(10) * R * T)
   }

inst/tinytest/test-DEW.R

@@ -1,3 +1,6 @@
+# Tests for DEW implementation in CHNOSZ
+# First version 20170925
+
 # Load default settings for CHNOSZ
 reset()
 
@@ -119,7 +122,8 @@ T <- c(300, 1100)
 P <- c(10000, 60000)
 # Adjust calculated logKs for different values of the
 # gas constant used in the DEW spreadsheet and CHNOSZ
-RoverR <- 1.9858775 / 1.9872
+#RoverR <- 1.9858775 / 1.9872 # 20170925
+RoverR <- 1.9858775 / (8.314463 / 4.184) # 20230630
 # Calculate logK for each reaction
 logK1 <- subcrt(c("H2O", "CO2", "H2CO3"), c(-1, -1, 1), T = T, P = P)$out$logK / RoverR
 logK2 <- subcrt(c("AlO2(SiO2)-", "H+", "Al+3", "H2O", "SiO2"), c(-1, -4, 1, 2, 1), T = T, P = P)$out$logK / RoverR

inst/tinytest/test-affinity.R

@@ -122,9 +122,15 @@ expect_equal(A.2303RT_max, A.2303RT_ref, tolerance = 1e-3, info = info)
 # TODO: add comparison with results from loading proteins via species()
 
 info <- "affinity() for proteins (with/without 'iprotein') returns same value as in previous package versions"
-# These values were calculated using versions 0.6, 0.8 and 0.9-7 (25 degrees C, 1 bar, basis species "CHNOS" or "CHNOS+")
-A.2303RT.nonionized <- -3795.297
-A.2303RT.ionized <- -3075.222
+## These values were calculated using versions 0.6, 0.8 and 0.9-7 (25 degrees C, 1 bar, basis species "CHNOS" or "CHNOS+")
+#A.2303RT.nonionized <- -3795.297
+#A.2303RT.ionized <- -3075.222
+## Calculated with version 2.0.0 (util.units() has R = 8.314445)
+#A.2303RT.nonionized <- -3795.291
+#A.2303RT.ionized <- -3075.215
+# Calculated with version 2.0.0-16 (util.units() has R = 8.314463)
+A.2303RT.nonionized <- -3794.69
+A.2303RT.ionized <- -3074.613
 # First for nonionized protein
 basis("CHNOS")
 # Try it with iprotein

inst/tinytest/test-logmolality.R

@@ -37,8 +37,8 @@ expect_equal(out1[[3]]$G - out0[[3]]$G, -convert(loggam_CO2, "G"), info = info)
 # What is the equilibrium constant for the reaction CO2 + H2O = H+ + HCO3-?
 # (this is the standard state property at IS = 0)
 logK <- subcrt(c("CO2", "H2O", "H+", "HCO3-"), c(-1, -1, 1, 1), T = 25)$out$logK
-# We get logK = -6.344694 (rounded)
-expect_true(maxdiff(logK, -6.344694) < 1e-6, info = info)
+# We get logK = -6.34468 (rounded)
+expect_true(maxdiff(logK, -6.34468) < 1e-6, info = info)
 
 ### What is the affinity of the reaction at pH=7 and molalities of HCO3- and CO2 = 10^-3?
 ## Case 1: ionic strength = 0, so gamma = 0 and activity = molality

inst/tinytest/test-protein.info.R

@@ -19,9 +19,22 @@ mod.OBIGT("[UPBB]", G = -21436, H = -45220, S = 1.62, E_units = "cal")
 basis("CHNOS+")
 suppressMessages(swap.basis("O2", "H2"))
 pequil <- protein.equil(protein, loga.protein=-3)
+
+# Before 20230630:
+# With R = 8.314445 (= 1.9872 * 4.184) in util.units(), ionize.aa(pinfo(protein)) gives:
+#            20        18
+#[1,] -56.06509 -55.87025
 # Astar/RT in the paragraph following Eq. 23, p. 6 of DS11
 # (truncated because of rounding)
-expect_true(any(grepl(c("0\\.435.*1\\.36"), pequil)), info = info)
+# expect_true(any(grepl(c("0\\.435.*1\\.36"), pequil)), info = info)
+
+# After 20230630:
+# With R = 8.314463 in util.units(), ionize.aa(pinfo(protein)) gives:
+#            20        18
+#[1,] -56.06511 -55.87027
+# The differences propagate up, so the test was changed on 20230630:
+expect_true(any(grepl(c("0\\.437.*1\\.36"), pequil)), info = info)
+
 # log10 activities of the proteins in the left-hand column of the same page
 expect_true(any(grepl(c("-3\\.256.*-2\\.834"), pequil)), info = info)