[CHNOSZ-commits] r796 - in pkg/CHNOSZ: . R inst/tinytest
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Jul 1 09:38:27 CEST 2023
Author: jedick
Date: 2023-07-01 09:38:26 +0200 (Sat, 01 Jul 2023)
New Revision: 796
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/AD.R
pkg/CHNOSZ/R/nonideal.R
pkg/CHNOSZ/R/protein.info.R
pkg/CHNOSZ/R/util.units.R
pkg/CHNOSZ/inst/tinytest/test-DEW.R
pkg/CHNOSZ/inst/tinytest/test-affinity.R
pkg/CHNOSZ/inst/tinytest/test-logmolality.R
pkg/CHNOSZ/inst/tinytest/test-protein.info.R
Log:
Update value of gas constant (8.314463 J K-1 mol-1)
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2023-06-24 01:59:05 UTC (rev 795)
+++ pkg/CHNOSZ/DESCRIPTION 2023-07-01 07:38:26 UTC (rev 796)
@@ -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 at R: c(
person("Jeffrey", "Dick", , "j3ffdick at gmail.com", role = c("aut", "cre"),
Modified: pkg/CHNOSZ/R/AD.R
===================================================================
--- pkg/CHNOSZ/R/AD.R 2023-06-24 01:59:05 UTC (rev 795)
+++ pkg/CHNOSZ/R/AD.R 2023-07-01 07:38:26 UTC (rev 796)
@@ -8,7 +8,8 @@
# 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
Modified: pkg/CHNOSZ/R/nonideal.R
===================================================================
--- pkg/CHNOSZ/R/nonideal.R 2023-06-24 01:59:05 UTC (rev 795)
+++ pkg/CHNOSZ/R/nonideal.R 2023-07-01 07:38:26 UTC (rev 796)
@@ -34,8 +34,10 @@
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) {
Modified: pkg/CHNOSZ/R/protein.info.R
===================================================================
--- pkg/CHNOSZ/R/protein.info.R 2023-06-24 01:59:05 UTC (rev 795)
+++ pkg/CHNOSZ/R/protein.info.R 2023-07-01 07:38:26 UTC (rev 796)
@@ -206,8 +206,11 @@
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
Modified: pkg/CHNOSZ/R/util.units.R
===================================================================
--- pkg/CHNOSZ/R/util.units.R 2023-06-24 01:59:05 UTC (rev 795)
+++ pkg/CHNOSZ/R/util.units.R 2023-07-01 07:38:26 UTC (rev 796)
@@ -100,8 +100,10 @@
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)
}
Modified: pkg/CHNOSZ/inst/tinytest/test-DEW.R
===================================================================
--- pkg/CHNOSZ/inst/tinytest/test-DEW.R 2023-06-24 01:59:05 UTC (rev 795)
+++ pkg/CHNOSZ/inst/tinytest/test-DEW.R 2023-07-01 07:38:26 UTC (rev 796)
@@ -1,3 +1,6 @@
+# Tests for DEW implementation in CHNOSZ
+# First version 20170925
+
# Load default settings for CHNOSZ
reset()
@@ -119,7 +122,8 @@
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
Modified: pkg/CHNOSZ/inst/tinytest/test-affinity.R
===================================================================
--- pkg/CHNOSZ/inst/tinytest/test-affinity.R 2023-06-24 01:59:05 UTC (rev 795)
+++ pkg/CHNOSZ/inst/tinytest/test-affinity.R 2023-07-01 07:38:26 UTC (rev 796)
@@ -122,9 +122,15 @@
# 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
Modified: pkg/CHNOSZ/inst/tinytest/test-logmolality.R
===================================================================
--- pkg/CHNOSZ/inst/tinytest/test-logmolality.R 2023-06-24 01:59:05 UTC (rev 795)
+++ pkg/CHNOSZ/inst/tinytest/test-logmolality.R 2023-07-01 07:38:26 UTC (rev 796)
@@ -37,8 +37,8 @@
# 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
Modified: pkg/CHNOSZ/inst/tinytest/test-protein.info.R
===================================================================
--- pkg/CHNOSZ/inst/tinytest/test-protein.info.R 2023-06-24 01:59:05 UTC (rev 795)
+++ pkg/CHNOSZ/inst/tinytest/test-protein.info.R 2023-07-01 07:38:26 UTC (rev 796)
@@ -19,9 +19,22 @@
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)
More information about the CHNOSZ-commits
mailing list