[CHNOSZ-commits] r714 - in pkg/CHNOSZ: . R demo inst inst/extdata/OBIGT inst/tinytest man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Mar 25 12:13:39 CET 2022
Author: jedick
Date: 2022-03-25 12:13:38 +0100 (Fri, 25 Mar 2022)
New Revision: 714
Removed:
pkg/CHNOSZ/vignettes/eos-regress.Rmd
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/AD.R
pkg/CHNOSZ/R/Berman.R
pkg/CHNOSZ/R/EOSregress.R
pkg/CHNOSZ/R/cgl.R
pkg/CHNOSZ/R/hkf.R
pkg/CHNOSZ/R/info.R
pkg/CHNOSZ/R/ionize.aa.R
pkg/CHNOSZ/R/logB_to_OBIGT.R
pkg/CHNOSZ/R/nonideal.R
pkg/CHNOSZ/R/protein.info.R
pkg/CHNOSZ/R/subcrt.R
pkg/CHNOSZ/R/util.data.R
pkg/CHNOSZ/R/util.formula.R
pkg/CHNOSZ/R/util.units.R
pkg/CHNOSZ/R/water.R
pkg/CHNOSZ/demo/E_coli.R
pkg/CHNOSZ/demo/affinity.R
pkg/CHNOSZ/demo/comproportionation.R
pkg/CHNOSZ/demo/lambda.R
pkg/CHNOSZ/demo/protein.equil.R
pkg/CHNOSZ/inst/NEWS.Rd
pkg/CHNOSZ/inst/extdata/OBIGT/AD.csv
pkg/CHNOSZ/inst/tinytest/test-AD.R
pkg/CHNOSZ/inst/tinytest/test-Berman.R
pkg/CHNOSZ/inst/tinytest/test-DEW.R
pkg/CHNOSZ/inst/tinytest/test-EOSregress.R
pkg/CHNOSZ/inst/tinytest/test-IAPWS95.R
pkg/CHNOSZ/inst/tinytest/test-affinity.R
pkg/CHNOSZ/inst/tinytest/test-ionize.aa.R
pkg/CHNOSZ/inst/tinytest/test-logmolality.R
pkg/CHNOSZ/inst/tinytest/test-protein.info.R
pkg/CHNOSZ/inst/tinytest/test-subcrt.R
pkg/CHNOSZ/inst/tinytest/test-util.R
pkg/CHNOSZ/inst/tinytest/test-water.R
pkg/CHNOSZ/man/Berman.Rd
pkg/CHNOSZ/man/eos.Rd
pkg/CHNOSZ/man/util.units.Rd
pkg/CHNOSZ/man/water.Rd
pkg/CHNOSZ/vignettes/anintro.Rmd
pkg/CHNOSZ/vignettes/customizing.Rmd
Log:
Use Joules internally in most functions
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2022-03-25 04:26:50 UTC (rev 713)
+++ pkg/CHNOSZ/DESCRIPTION 2022-03-25 11:13:38 UTC (rev 714)
@@ -1,6 +1,6 @@
Date: 2022-03-25
Package: CHNOSZ
-Version: 1.4.3-6
+Version: 1.4.3-7
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 2022-03-25 04:26:50 UTC (rev 713)
+++ pkg/CHNOSZ/R/AD.R 2022-03-25 11:13:38 UTC (rev 714)
@@ -25,9 +25,9 @@
# Density (g cm-3)
rho1 <- waterTP$rho / 1000
# Entropy (dimensionless)
- S1 <- waterTP$S / 1.9872
+ S1 <- waterTP$S / R
# Heat capacity (dimensionless)
- Cp1 <- waterTP$Cp / 1.9872
+ Cp1 <- waterTP$Cp / R
# Volume (cm3 mol-1)
V1 <- waterTP$V
# Calculate properties of ideal H2O gas
@@ -54,13 +54,11 @@
if(property[[j]] == "G") {
# Get gas properties (J mol-1)
- G_gas <- convert(gasprops$G, "J", T = T)
+ G_gas <- gasprops$G
# Calculate G_hyd (J mol-1)
G_hyd <- R*T * ( -log(NW) + (1 - PAR$xi) * log(f1) + PAR$xi * log(RV * T * rho1 / MW) + rho1 * (PAR$a + PAR$b * (1000/T)^0.5) )
# Calculate the chemical potential (J mol-1)
G <- G_gas + G_hyd
- # Convert J to cal
- G <- convert(G, "cal", T = T)
# Insert into data frame of properties
myprops$G <- G
}
@@ -67,7 +65,7 @@
if(property[[j]] == "S") {
# Get S_gas
- S_gas <- convert(gasprops$S, "J", T = T)
+ S_gas <- gasprops$S
# Calculate S_hyd
S_hyd <- R * (
(1 - PAR$xi) * (S1 - S1_g)
@@ -79,13 +77,12 @@
)
)
S <- S_gas + S_hyd
- S <- convert(S, "cal", T = T)
myprops$S <- S
}
if(property[[j]] == "Cp") {
# Get Cp_gas
- Cp_gas <- convert(gasprops$Cp, "J", T = T)
+ Cp_gas <- gasprops$Cp
# Calculate Cp_hyd
Cp_hyd <- R * (
(1 - PAR$xi) * (Cp1 - Cp1_g)
@@ -95,7 +92,6 @@
+ PAR$b * (-0.25 * 10^1.5 * T^-1.5 * rho1 + 10^1.5 * T^-0.5 * drho1_dT + 10^1.5 * T^0.5 * d2rho1_dT2)
)
Cp <- Cp_gas + Cp_hyd
- Cp <- convert(Cp, "cal", T = T)
myprops$Cp <- Cp
}
@@ -109,7 +105,7 @@
}
}
- # Calculate enthalpy (NOTE: this is in calories) 20220206
+ # Calculate enthalpy 20220206
myprops$H <- myprops$G - 298.15 * entropy(PAR$formula) + T * myprops$S
out[[i]] <- myprops
}
@@ -122,7 +118,7 @@
# Get H2O fugacity (bar)
GH2O_P <- water("G", T = T, P = P)$G
GH2O_1 <- water("G", T = T, P = 1)$G
- f1 <- exp ( (GH2O_P - GH2O_1) / (1.9872 * T) )
+ f1 <- exp ( (GH2O_P - GH2O_1) / (8.31441 * T) )
# For Psat, calculate the real liquid-vapor curve (not 1 bar below 100 degC)
if(isPsat) {
P <- water("Psat", T = T, P = "Psat", P1 = FALSE)$Psat
Modified: pkg/CHNOSZ/R/Berman.R
===================================================================
--- pkg/CHNOSZ/R/Berman.R 2022-03-25 04:26:50 UTC (rev 713)
+++ pkg/CHNOSZ/R/Berman.R 2022-03-25 11:13:38 UTC (rev 714)
@@ -4,7 +4,7 @@
# in the system Na2O-K2O-CaO-MgO-FeO-Fe2O3-Al2O3-SiO2-TiO2-H2O-CO2.
# J. Petrol. 29, 445-522. https://doi.org/10.1093/petrology/29.2.445
-Berman <- function(name, T = 298.15, P = 1, check.G=FALSE, calc.transition=TRUE, calc.disorder=TRUE, units="cal") {
+Berman <- function(name, T = 298.15, P = 1, check.G=FALSE, calc.transition=TRUE, calc.disorder=TRUE) {
# Reference temperature and pressure
Pr <- 1
Tr <- 298.15
@@ -55,7 +55,7 @@
# Get the entropy of the elements using the chemical formula in thermo()$OBIGT
OBIGT <- thermo()$OBIGT
formula <- OBIGT$formula[match(name, OBIGT$name)]
- SPrTr_elements <- convert(entropy(formula), "J")
+ SPrTr_elements <- entropy(formula)
# Check that G in data file is the G of formation from the elements --> Benson-Helgeson convention (DG = DH - T*DS)
if(check.G) {
GfPrTr_calc <- HfPrTr - Tr * (SPrTr - SPrTr_elements)
@@ -189,15 +189,8 @@
# The output will just have "G" and "H"
G <- Gf
H <- Ha
- # Convert J to cal
- if(grepl("cal", units)) {
- G <- convert(Gf, "cal")
- H <- convert(Ha, "cal")
- S <- convert(S, "cal")
- Cp <- convert(Cp, "cal")
- }
# Convert J/bar to cm^3/mol
V <- V * 10
- data.frame(T=T, P=P, G=G, H=H, S=S, Cp=Cp, V=V)
+ data.frame(T = T, P = P, G = G, H = H, S = S, Cp = Cp, V = V)
}
Modified: pkg/CHNOSZ/R/EOSregress.R
===================================================================
--- pkg/CHNOSZ/R/EOSregress.R 2022-03-25 04:26:50 UTC (rev 713)
+++ pkg/CHNOSZ/R/EOSregress.R 2022-03-25 11:13:38 UTC (rev 714)
@@ -199,7 +199,7 @@
names(out) <- c("(Intercept)", "invTTheta2", "TXBorn")
} else if(property=="V") {
iis <- iis[,c("a1", "a2", "a3", "a4", "omega")]
- # calculate sigma and xi and convert to volumetric units: 1 cal = 41.84 cm^3 bar
+ # calculate sigma and xi and convert to volumetric units: 1 J = 10 cm^3 bar
sigma <- convert( iis$a1 + iis$a2 / (2600 + P), "cm3bar" )
xi <- convert( iis$a3 + iis$a4 / (2600 + P), "cm3bar" )
omega <- convert( iis$omega, "cm3bar" )
Modified: pkg/CHNOSZ/R/cgl.R
===================================================================
--- pkg/CHNOSZ/R/cgl.R 2022-03-25 04:26:50 UTC (rev 713)
+++ pkg/CHNOSZ/R/cgl.R 2022-03-25 11:13:38 UTC (rev 714)
@@ -1,36 +1,37 @@
# CHNOSZ/cgl.R
-# calculate standard thermodynamic properties of non-aqueous species
+# Calculate standard thermodynamic properties of non-aqueous species
# 20060729 jmd
cgl <- function(property = NULL, parameters = NULL, T = 298.15, P = 1) {
- # calculate properties of crystalline, liquid (except H2O) and gas species
+ # Calculate properties of crystalline, liquid (except H2O) and gas species
Tr <- 298.15
Pr <- 1
- # the number of T, P conditions
+ # The number of T, P conditions
ncond <- max(c(length(T), length(P)))
- # initialize output list
+ # Initialize output list
out <- list()
- # loop over each species
+ # Loop over each species
for(k in 1:nrow(parameters)) {
- # the parameters for *this* species
+ # The parameters for *this* species
PAR <- parameters[k, ]
if(all(is.na(PAR[9:21]))) {
- # use Berman equations (parameters not in thermo()$OBIGT)
+ # Use Berman equations (parameters not in thermo()$OBIGT)
properties <- Berman(PAR$name, T = T, P = P)
iprop <- match(property, colnames(properties))
values <- properties[, iprop, drop=FALSE]
} else {
- # in CHNOSZ, we have
+ # In CHNOSZ, we have
# 1 cm^3 bar --> convert(1, "calories") == 0.02390057 cal
# but REAC92D.F in SUPCRT92 uses
cm3bar_to_cal <- 0.023901488 # cal
- # start with NA values
+ cm3bar_to_J <- convert(cm3bar_to_cal, "J")
+ # Start with NA values
values <- data.frame(matrix(NA, ncol = length(property), nrow=ncond))
colnames(values) <- property
- # a test for availability of heat capacity coefficients (a, b, c, d, e, f)
+ # A test for availability of heat capacity coefficients (a, b, c, d, e, f)
# based on the column assignments in thermo()$OBIGT
if(any(!is.na(PAR[, 14:19]))) {
- # we have at least one of the heat capacity coefficients;
+ # We have at least one of the heat capacity coefficients;
# zero out any NA's in the rest (leave lambda and T of transition (columns 19-20) alone)
PAR[, 14:19][, is.na(PAR[, 14:19])] <- 0
# calculate the heat capacity and its integrals
@@ -37,7 +38,7 @@
Cp <- PAR$a + PAR$b*T + PAR$c*T^-2 + PAR$d*T^-0.5 + PAR$e*T^2 + PAR$f*T^PAR$lambda
intCpdT <- PAR$a*(T - Tr) + PAR$b*(T^2 - Tr^2)/2 + PAR$c*(1/T - 1/Tr)/-1 + PAR$d*(T^0.5 - Tr^0.5)/0.5 + PAR$e*(T^3-Tr^3)/3
intCpdlnT <- PAR$a*log(T / Tr) + PAR$b*(T - Tr) + PAR$c*(T^-2 - Tr^-2)/-2 + PAR$d*(T^-0.5 - Tr^-0.5)/-0.5 + PAR$e*(T^2 - Tr^2)/2
- # do we also have the lambda parameter (Cp term with adjustable exponent on T)?
+ # Do we also have the lambda parameter (Cp term with adjustable exponent on T)?
if(!is.na(PAR$lambda) & !identical(PAR$lambda, 0)) {
# equations for lambda adapted from Helgeson et al., 1998 (doi:10.1016/S0016-7037(97)00219-6)
if(PAR$lambda == -1) intCpdT <- intCpdT + PAR$f*log(T/Tr)
@@ -45,32 +46,32 @@
intCpdlnT <- intCpdlnT + PAR$f*(T^PAR$lambda - Tr^PAR$lambda) / PAR$lambda
}
} else {
- # use constant heat capacity if the coefficients are not available
+ # Use constant heat capacity if the coefficients are not available
Cp <- PAR$Cp
intCpdT <- PAR$Cp*(T - Tr)
intCpdlnT <- PAR$Cp*log(T / Tr)
- # in case Cp is listed as NA, set the integrals to 0 at Tr
+ # In case Cp is listed as NA, set the integrals to 0 at Tr
intCpdT[T==Tr] <- 0
intCpdlnT[T==Tr] <- 0
}
- # volume and its integrals
+ # Volume and its integrals
if(PAR$name %in% c("quartz", "coesite")) {
- # volume calculations for quartz and coesite
+ # Volume calculations for quartz and coesite
qtz <- quartz_coesite(PAR, T, P)
V <- qtz$V
intVdP <- qtz$intVdP
intdVdTdP <- qtz$intdVdTdP
} else {
- # for other minerals, volume is constant (Helgeson et al., 1978)
+ # For other minerals, volume is constant (Helgeson et al., 1978)
V <- rep(PAR$V, ncond)
- # if the volume is NA, set its integrals to zero
+ # If the volume is NA, set its integrals to zero
if(is.na(PAR$V)) intVdP <- intdVdTdP <- numeric(ncond)
else {
- intVdP <- PAR$V*(P - Pr) * cm3bar_to_cal
+ intVdP <- PAR$V*(P - Pr) * cm3bar_to_J
intdVdTdP <- 0
}
}
- # get the values of each of the requested thermodynamic properties
+ # Get the values of each of the requested thermodynamic properties
for(i in 1:length(property)) {
if(property[i] == "Cp") values[, i] <- Cp
if(property[i] == "V") values[, i] <- V
@@ -77,7 +78,7 @@
if(property[i] == "E") values[, i] <- rep(NA, ncond)
if(property[i] == "kT") values[, i] <- rep(NA, ncond)
if(property[i] == "G") {
- # calculate S * (T - Tr), but set it to 0 at Tr (in case S is NA)
+ # Calculate S * (T - Tr), but set it to 0 at Tr (in case S is NA)
Sterm <- PAR$S*(T - Tr)
Sterm[T==Tr] <- 0
values[, i] <- PAR$G - Sterm + intCpdT - T*intCpdlnT + intVdP
@@ -85,13 +86,14 @@
if(property[i] == "H") values[, i] <- PAR$H + intCpdT + intVdP - T*intdVdTdP
if(property[i] == "S") values[, i] <- PAR$S + intCpdlnT - intdVdTdP
}
- } # end calculations using parameters from thermo()$OBIGT
+ } # End calculations using parameters from thermo()$OBIGT
out[[k]] <- values
- } # end loop over species
+ } # End loop over species
+
return(out)
}
-### unexported function ###
+### Unexported function ###
# calculate GHS and V corrections for quartz and coesite 20170929
# (these are the only mineral phases for which SUPCRT92 uses an inconstant volume)
@@ -109,7 +111,7 @@
ca <- -0.4973e-4
VPtTta <- 23.348
VPrTtb <- 23.72
- Stran <- 0.342
+ Stran <- convert(0.342, "J")
# constants from REAC92D.f
VPrTra <- 22.688 # VPrTr(a-quartz)
Vdiff <- 2.047 # VPrTr(a-quartz) - VPrTr(coesite)
@@ -134,11 +136,12 @@
V <- V - Vdiff
}
cm3bar_to_cal <- 0.023901488
- GVterm <- cm3bar_to_cal * (VPrTra * (P - Pstar) + VPrTtb * (Pstar - Pr) -
+ cm3bar_to_J <- convert(cm3bar_to_cal, "J")
+ GVterm <- cm3bar_to_J * (VPrTra * (P - Pstar) + VPrTtb * (Pstar - Pr) -
0.5 * ca * (2 * Pr * (P - Pstar) - (P^2 - Pstar^2)) -
ca * k * (T - Tr) * (P - Pstar) +
k * (ba + aa * ca * k) * (T - Tr) * log((aa + P/k) / (aa + Pstar/k)))
- SVterm <- cm3bar_to_cal * (-k * (ba + aa * ca * k) *
+ SVterm <- cm3bar_to_J * (-k * (ba + aa * ca * k) *
log((aa + P/k) / (aa + Pstar/k)) + ca * k * (P - Pstar)) - Sstar
# note the minus sign on "SVterm" in order that intdVdTdP has the correct sign
list(intVdP=GVterm, intdVdTdP=-SVterm, V=V)
Modified: pkg/CHNOSZ/R/hkf.R
===================================================================
--- pkg/CHNOSZ/R/hkf.R 2022-03-25 04:26:50 UTC (rev 713)
+++ pkg/CHNOSZ/R/hkf.R 2022-03-25 11:13:38 UTC (rev 714)
@@ -1,43 +1,42 @@
# CHNOSZ/hkf.R
-# calculate thermodynamic properties using equations of state
+# Calculate thermodynamic properties using equations of state
# 11/17/03 jmd
-## if this file is interactively sourced, the following is also needed to provide unexported functions:
+## If this file is interactively sourced, the following is also needed to provide unexported functions:
#source("util.args.R")
hkf <- function(property = NULL, parameters = NULL, T = 298.15, P = 1,
contrib = c("n", "s", "o"), H2O.props="rho") {
- # calculate G, H, S, Cp, V, kT, and/or E using
- # the revised HKF equations of state
+ # Calculate G, H, S, Cp, V, kT, and/or E using the revised HKF equations of state
# H2O.props - H2O properties needed for subcrt() output
- # constants
+ # Constants
Tr <- 298.15 # K
Pr <- 1 # bar
Theta <- 228 # K
Psi <- 2600 # bar
- # make T and P equal length
+ # Make T and P equal length
if(!identical(P, "Psat")) {
if(length(P) < length(T)) P <- rep(P, length.out = length(T))
if(length(T) < length(P)) T <- rep(T, length.out = length(P))
}
- # nonsolvation, solvation, and origination contribution
+ # Nonsolvation, solvation, and origination contribution
notcontrib <- ! contrib %in% c("n", "s", "o")
if(TRUE %in% notcontrib) stop(paste("contrib must be in c('n', 's', 'o); got", c2s(contrib[notcontrib])))
- # get water properties
+ # Get water properties
# rho - for subcrt() output and g function
# Born functions and epsilon - for HKF calculations
H2O.props <- c(H2O.props, "QBorn", "XBorn", "YBorn", "epsilon")
thermo <- get("thermo", CHNOSZ)
if(grepl("SUPCRT", thermo$opt$water)) {
- # using H2O92D.f from SUPCRT92: alpha, daldT, beta - for partial derivatives of omega (g function)
+ # Using H2O92D.f from SUPCRT92: alpha, daldT, beta - for partial derivatives of omega (g function)
H2O.props <- c(H2O.props, "alpha", "daldT", "beta")
}
if(grepl("IAPWS", thermo$opt$water)) {
- # using IAPWS-95: NBorn, UBorn - for compressibility, expansibility
+ # Using IAPWS-95: NBorn, UBorn - for compressibility, expansibility
H2O.props <- c(H2O.props, "NBorn", "UBorn")
}
if(grepl("DEW", thermo$opt$water)) {
- # using DEW model: get beta to calculate dgdP
+ # Using DEW model: get beta to calculate dgdP
H2O.props <- c(H2O.props, "beta")
}
H2O <- water(H2O.props, T = c(Tr, T), P = c(Pr, P))
@@ -45,38 +44,39 @@
H2O.PT <- H2O[-1, ]
ZBorn <- -1 / H2O.PT$epsilon
ZBorn.PrTr <- -1 / H2O.PrTr$epsilon
- # a list to store the result
+ # A list to store the result
aq.out <- list()
nspecies <- nrow(parameters)
for(k in seq_len(nspecies)) {
- # loop over each species
+ # Loop over each species
PAR <- parameters[k, ]
- # substitute Cp and V for missing EoS parameters
- # here we assume that the parameters are in the same position as in thermo()$OBIGT
- # we don't need this if we're just looking at solvation properties (Cp_s_var, V_s_var)
+ # Substitute Cp and V for missing EoS parameters
+ # Here we assume that the parameters are in the same position as in thermo()$OBIGT
+ # We don't need this if we're just looking at solvation properties (Cp_s_var, V_s_var)
if("n" %in% contrib) {
- # put the heat capacity in for c1 if both c1 and c2 are missing
+ # Put the heat capacity in for c1 if both c1 and c2 are missing
if(all(is.na(PAR[, 18:19]))) PAR[, 18] <- PAR$Cp
- # put the volume in for a1 if a1, a2, a3 and a4 are missing
- if(all(is.na(PAR[, 14:17]))) PAR[, 14] <- convert(PAR$V, "calories")
- # test for availability of the EoS parameters
+ # Put the volume in for a1 if a1, a2, a3 and a4 are missing
+ if(all(is.na(PAR[, 14:17]))) PAR[, 14] <- convert(PAR$V, "joules")
+ # Test for availability of the EoS parameters
hasEOS <- any(!is.na(PAR[, 14:21]))
- # if at least one of the EoS parameters is available, zero out any NA's in the rest
+ # If at least one of the EoS parameters is available, zero out any NA's in the rest
if(hasEOS) PAR[, 14:21][, is.na(PAR[, 14:21])] <- 0
}
- # compute values of omega(P,T) from those of omega(Pr,Tr)
- # using g function etc. (Shock et al., 1992 and others)
+ # Compute values of omega(P,T) from those of omega(Pr,Tr)
+ # Using g function etc. (Shock et al., 1992 and others)
omega <- PAR$omega # omega.PrTr
- # its derivatives are zero unless the g function kicks in
+ # The derivatives are zero unless the g function kicks in
dwdP <- dwdT <- d2wdT2 <- numeric(length(T))
Z <- PAR$Z
omega.PT <- rep(PAR$omega, length(T))
if(!identical(Z, 0) & !is.na(Z) & !identical(PAR$name, "H+")) {
- # compute derivatives of omega: g and f functions (Shock et al., 1992; Johnson et al., 1992)
+ # Compute derivatives of omega: g and f functions (Shock et al., 1992; Johnson et al., 1992)
rhohat <- H2O.PT$rho/1000 # just converting kg/m3 to g/cm3
g <- gfun(rhohat, convert(T, "C"), P, H2O.PT$alpha, H2O.PT$daldT, H2O.PT$beta)
- # after SUPCRT92/reac92.f
- eta <- 1.66027E5
+ # After SUPCRT92/reac92.f
+ # eta needs to be converted to Joules! 20220325
+ eta <- convert(1.66027E5, "J")
reref <- Z^2 / (omega/eta + Z/(3.082 + 0))
re <- reref + abs(Z) * g$g
omega.PT <- eta * (Z^2/re - Z/(3.082 + g$g))
@@ -86,16 +86,16 @@
dwdT <- (-eta * Z3 * g$dgdT)
d2wdT2 <- (2 * eta * Z4 * g$dgdT^2 - eta * Z3 * g$d2gdT2)
}
- # loop over each property
+ # Loop over each property
w <- NULL
for(i in 1:length(property)) {
PROP <- property[i]
- # over nonsolvation, solvation, or origination contributions
+ # Over nonsolvation, solvation, or origination contributions
hkf.p <- numeric(length(T))
for(icontrib in contrib) {
- # various contributions to the properties
+ # Various contributions to the properties
if(icontrib == "n") {
- # nonsolvation ghs equations
+ # Nonsolvation ghs equations
if(PROP == "H") {
p.c <- PAR$c1*(T-Tr) - PAR$c2*(1/(T-Theta)-1/(Tr-Theta))
p.a <- PAR$a1*(P-Pr) + PAR$a2*log((Psi+P)/(Psi+Pr)) +
@@ -114,9 +114,9 @@
p.a <- PAR$a1*(P-Pr) + PAR$a2*log((Psi+P)/(Psi+Pr)) +
(PAR$a3*(P-Pr) + PAR$a4*log((Psi+P)/(Psi+Pr)))/(T-Theta)
p <- p.c + p.a
- # at Tr,Pr, if the origination contribution is not NA, ensure the solvation contribution is 0, not NA
+ # If the origination contribution is not NA at Tr,Pr, ensure the solvation contribution is 0, not NA
if(!is.na(PAR$G)) p[T==Tr & P==Pr] <- 0
- # nonsolvation cp v kt e equations
+ # Nonsolvation cp v kt e equations
} else if(PROP == "Cp") {
p <- PAR$c1 + PAR$c2 * ( T - Theta ) ^ (-2)
} else if(PROP == "V") {
@@ -127,15 +127,15 @@
p <- (convert(PAR$a2, "cm3bar") +
convert(PAR$a4, "cm3bar") / (T - Theta)) * (Psi + P) ^ (-2)
} else if(PROP == "E") {
- p <- convert( - (PAR$a3 + PAR$a4 / convert((Psi + P), "calories")) *
+ p <- convert( - (PAR$a3 + PAR$a4 / convert((Psi + P), "joules")) *
(T - Theta) ^ (-2), "cm3bar")
}
}
if(icontrib == "s") {
- # solvation ghs equations
+ # Solvation ghs equations
if(PROP == "G") {
p <- -omega.PT*(ZBorn+1) + omega*(ZBorn.PrTr+1) + omega*H2O.PrTr$YBorn*(T-Tr)
- # at Tr,Pr, if the origination contribution is not NA, ensure the solvation contribution is 0, not NA
+ # If the origination contribution is not NA at Tr,Pr, ensure the solvation contribution is 0, not NA
if(!is.na(PAR$G)) p[T==Tr & P==Pr] <- 0
}
if(PROP == "H")
@@ -143,7 +143,7 @@
omega*(ZBorn.PrTr+1) - omega*Tr*H2O.PrTr$YBorn
if(PROP == "S")
p <- omega.PT*H2O.PT$YBorn + (ZBorn+1)*dwdT - omega*H2O.PrTr$YBorn
- # solvation cp v kt e equations
+ # Solvation cp v kt e equations
if(PROP == "Cp") p <- omega.PT*T*H2O.PT$XBorn + 2*T*H2O.PT$YBorn*dwdT +
T*(ZBorn+1)*d2wdT2
if(PROP == "V") p <- -convert(omega.PT, "cm3bar") *
@@ -154,18 +154,18 @@
if(PROP == "E") p <- -convert(omega, "cm3bar") * H2O.PT$UBorn
}
if(icontrib == "o") {
- # origination ghs equations
+ # Origination ghs equations
if(PROP == "G") {
p <- PAR$G - PAR$S * (T-Tr)
- # don"t inherit NA from PAR$S at Tr
+ # Don't inherit NA from PAR$S at Tr
p[T==Tr] <- PAR$G
}
else if(PROP == "H") p <- PAR$H
else if(PROP == "S") p <- PAR$S
- # origination eos equations (Cp, V, kT, E): senseless
+ # Origination eos equations (Cp, V, kT, E): senseless
else p <- 0 * T
}
- # accumulate the contribution
+ # Accumulate the contribution
hkf.p <- hkf.p + p
}
wnew <- data.frame(hkf.p)
@@ -174,10 +174,11 @@
colnames(w) <- property
aq.out[[k]] <- w
}
- return(list(aq=aq.out, H2O=H2O.PT))
+
+ return(list(aq = aq.out, H2O = H2O.PT))
}
-### unexported functions ###
+### Unexported functions ###
gfun <- function(rhohat, Tc, P, alpha, daldT, beta) {
## g and f functions for describing effective electrostatic radii of ions
Modified: pkg/CHNOSZ/R/info.R
===================================================================
--- pkg/CHNOSZ/R/info.R 2022-03-25 04:26:50 UTC (rev 713)
+++ pkg/CHNOSZ/R/info.R 2022-03-25 11:13:38 UTC (rev 714)
@@ -228,7 +228,6 @@
} else if(isBerman) {
# Calculate Cp for Berman minerals 20220208
calcCp <- Berman(this$name)$Cp
- if(this$E_units == "J") calcCp <- convert(calcCp, "J")
this$Cp <- calcCp
}
# check tabulated volumes - only for aq (HKF equation)
Modified: pkg/CHNOSZ/R/ionize.aa.R
===================================================================
--- pkg/CHNOSZ/R/ionize.aa.R 2022-03-25 04:26:50 UTC (rev 713)
+++ pkg/CHNOSZ/R/ionize.aa.R 2022-03-25 11:13:38 UTC (rev 714)
@@ -2,7 +2,7 @@
# rewritten ionization function 20120526 jmd
ionize.aa <- function(aa, property="Z", T=25, P="Psat", pH=7, ret.val=NULL, suppress.Cys=FALSE) {
- # calculate the additive ionization property of proteins with amino acid
+ # Calculate the additive ionization property of proteins with amino acid
# composition in aa as a function of vectors of T, P and pH;
# property if NULL is the net charge, if not NULL is one of the subcrt() properties
# or "A" to calculate A/2.303RT for the ionization reaction
@@ -12,64 +12,64 @@
T <- rep(T, length.out=lmax)
P <- rep(P, length.out=lmax)
pH <- rep(pH, length.out=lmax)
- # turn pH into a matrix with as many columns as ionizable groups
+ # Turn pH into a matrix with as many columns as ionizable groups
pH <- matrix(rep(pH, 9), ncol=9)
- # turn charges into a matrix with as many rows as T,P,pH conditions
+ # Turn charges into a matrix with as many rows as T,P,pH conditions
charges <- c(-1, -1, -1, 1, 1, 1, -1, 1, -1)
charges <- matrix(rep(charges, lmax), nrow=lmax, byrow=TRUE)
- # the rownumbers of the ionizable groups in thermo()$OBIGT
+ # The rownumbers of the ionizable groups in thermo()$OBIGT
neutral <- c("[Cys]", "[Asp]", "[Glu]", "[His]", "[Lys]", "[Arg]", "[Tyr]", "[AABB]", "[AABB]")
charged <- c("[Cys-]","[Asp-]","[Glu-]","[His+]","[Lys+]","[Arg+]","[Tyr-]","[AABB+]","[AABB-]")
ineutral <- info(neutral, "aq")
icharged <- info(charged, "aq")
- # we'll only call subcrt() with the unique pressure/temperature combinations
+ # We'll only call subcrt() with the unique pressure/temperature combinations
pTP <- paste(T, P)
dupPT <- duplicated(pTP)
- # what property are we after
+ # What property are we after
sprop <- c("G", property)
if(property %in% c("A", "Z")) sprop <- "G"
- # Use convert=FALSE so we get results in calories 20210407
+ # Use convert=FALSE so we get results in Joules 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
+ # 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
+ # Keep it as a matrix even if we have only one unique T, P-combo
if(length(pTP[!dupPT])==1) Gs <- t(Gs)
- # now the Gibbs energy difference for each group
+ # Now the Gibbs energy difference for each group
DG <- Gs[, 10:18, drop=FALSE] - Gs[, 1:9, drop=FALSE]
- # build a matrix with one row for each of the (possibly duplicated) T, P values
+ # Build a matrix with one row for each of the (possibly duplicated) T, P values
uPT <- unique(pTP)
DG <- t(sapply(pTP, function(x) DG[match(x, uPT), , drop=FALSE]))
- # the pK values (-logK)
+ # The pK values (-logK)
DG <- DG * charges
- 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
+ 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")) {
colnames(pK) <- charged
return(pK)
}
- # now to calculate alpha! - degrees of formation of the charged groups
+ # Now to calculate alpha! - degrees of formation of the charged groups
alpha <- 1 / (1 + 10 ^ (charges * (pH - pK)))
- # suppress cysteine ionization if requested
+ # Suppress cysteine ionization if requested
if(suppress.Cys) alpha[, 1] <- 0
if(identical(ret.val, "alpha")) return(alpha)
- # now to calculate the properties of the ionizable groups - can be charges,
+ # Now to calculate the properties of the ionizable groups - can be charges,
# the chemical affinities of the ionization reactions,
# or another property from subcrt()
if(identical(property, "Z")) aavals <- charges
else if(identical(property, "A")) aavals <- - charges * (pH - pK)
else {
- # it's not charge, so compile it from the subcrt output
+ # It's not charge, so compile it from the subcrt output
# the property-values
icol <- match(property, colnames(sout[[1]]))
aavals <- sapply(sout, function(x) x[,icol])
- # keep it as a matrix even if we have only one unique T, P-combo
+ # Keep it as a matrix even if we have only one unique T, P-combo
if(length(pTP[!dupPT])==1) aavals <- t(aavals)
- # build a matrix with one row for each of the (possibly duplicated) T, P values
+ # Build a matrix with one row for each of the (possibly duplicated) T, P values
aavals <- t(sapply(pTP, function(x) aavals[match(x, uPT), , drop=FALSE], USE.NAMES=FALSE))
- # the property difference for each group
+ # The property difference for each group
aavals <- aavals[, 10:18, drop=FALSE] - aavals[, 1:9, drop=FALSE]
}
if(identical(ret.val, "aavals")) {
@@ -76,21 +76,21 @@
colnames(aavals) <- charged
return(aavals)
}
- # the contribution from each group to the ionization property of the protein
+ # The contribution from each group to the ionization property of the protein
aavals <- aavals * alpha
- # now work with 'aa'; the next line is so that a missing argument shows
- # the name of this function in the error message
+ # Now work with 'aa'; the next line is so that a missing argument shows
+ # The name of this function in the error message
aa <- aa
- # the columns where we find the counts of ionizable groups
+ # The columns where we find the counts of ionizable groups
iionize <- match(c("Cys", "Asp", "Glu", "His", "Lys", "Arg", "Tyr", "chains", "chains"), colnames(aa))
aa <- as.matrix(aa[, iionize])
- # add it all up
+ # Add it all up
out <- apply(aa, 1, function(x) {
aavals %*% x
})
- # keep it as a matrix even if we have only one T, P-combo
+ # Keep it as a matrix even if we have only one T, P-combo
if(lmax==1) out <- t(out)
rownames(out) <- rownames(aavals)
- # that's all folks!
+ # That's all folks!
return(out)
}
Modified: pkg/CHNOSZ/R/logB_to_OBIGT.R
===================================================================
--- pkg/CHNOSZ/R/logB_to_OBIGT.R 2022-03-25 04:26:50 UTC (rev 713)
+++ pkg/CHNOSZ/R/logB_to_OBIGT.R 2022-03-25 11:13:38 UTC (rev 714)
@@ -23,11 +23,6 @@
Gr <- convert(logB, "G", TK)
# logK0 gives values for ΔG°r of the reaction with ΔG°f = 0 for the formed species
Gr0 <- convert(logK0, "G", TK)
- # TODO: fix convert() so that it uses Joules
- if(E.units() == "J") {
- Gr <- convert(Gr, "J")
- Gr0 <- convert(Gr0, "J")
- }
# Calculate ΔG°f of the formed species
Gf <- Gr - Gr0
Modified: pkg/CHNOSZ/R/nonideal.R
===================================================================
--- pkg/CHNOSZ/R/nonideal.R 2022-03-25 04:26:50 UTC (rev 713)
+++ pkg/CHNOSZ/R/nonideal.R 2022-03-25 11:13:38 UTC (rev 714)
@@ -34,8 +34,8 @@
else stop("invalid method (", thermo$opt$nonideal, ")")
}
- R <- 1.9872 # gas constant, cal K^-1 mol^-1
- #R <- 8.31446261815324 # gas constant, J K^-1 mol^-1 20220325
+ #R <- 1.9872 # gas constant, cal K^-1 mol^-1
+ R <- 8.314445 # gas constant, J K^-1 mol^-1 20220325
# 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 2022-03-25 04:26:50 UTC (rev 713)
+++ pkg/CHNOSZ/R/protein.info.R 2022-03-25 11:13:38 UTC (rev 714)
@@ -149,8 +149,6 @@
}
protein.equil <- function(protein, T=25, loga.protein=0, digits=4) {
- # For now we have to use calories 20220325
- if(thermo()$opt$E.units != "cal") stop('please run E.units("cal") first')
out <- character()
mymessage <- function(...) {
message(...)
@@ -159,6 +157,9 @@
}
# show the individual steps in calculating metastable equilibrium among proteins
mymessage("protein.equil: temperature from argument is ", T, " degrees C")
+ # Display units
+ E_units <- E.units()
+ mymessage("protein.equil: energy units is ", E_units)
TK <- convert(T, "K")
# get the amino acid compositions of the proteins
aa <- pinfo(pinfo(protein))
@@ -189,18 +190,21 @@
G0basissum <- colSums(t(protbasis) * G0basis)
# standard Gibbs energies of nonionized proteins
G0prot <- unlist(suppressMessages(subcrt(pname, T=T, property="G")$out))
- # standard Gibbs energy of formation reaction of nonionized protein, cal/mol
+ # standard Gibbs energy of formation reaction of nonionized protein, E_units/mol
G0protform <- G0prot - G0basissum
- mymessage("protein.equil [1]: reaction to form nonionized protein from basis species has G0(cal/mol) of ", signif(G0protform[1], digits))
+ mymessage("protein.equil [1]: reaction to form nonionized protein from basis species has G0(", E_units, "/mol) of ", signif(G0protform[1], digits))
if(ionize.it) {
- # standard Gibbs energy of ionization of protein, cal/mol
+ # standard Gibbs energy of ionization of protein, J/mol
G0ionization <- suppressMessages(ionize.aa(aa, property="G", T=T, pH=pH))[1, ]
- mymessage("protein.equil [1]: ionization reaction of protein has G0(cal/mol) of ", signif(G0ionization[1], digits))
- # standard Gibbs energy of formation reaction of ionized protein, cal/mol
+ # standard Gibbs energy of ionization of protein, E_units/mol
+ if(E_units == "cal") G0ionization <- convert(G0ionization, "cal")
+ mymessage("protein.equil [1]: ionization reaction of protein has G0(", E_units, "/mol) of ", signif(G0ionization[1], digits))
+ # standard Gibbs energy of formation reaction of ionized protein, E_units/mol
G0protform <- G0protform + G0ionization
}
# standard Gibbs energy of formation reaction of non/ionized residue equivalents, dimensionless
- R <- 1.9872 # gas constant, cal K^-1 mol^-1
+ 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
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/subcrt.R
===================================================================
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 714
More information about the CHNOSZ-commits
mailing list