[CHNOSZ-commits] r760 - in pkg/CHNOSZ: . R man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Dec 9 07:17:33 CET 2022
Author: jedick
Date: 2022-12-09 07:17:33 +0100 (Fri, 09 Dec 2022)
New Revision: 760
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/logB.to.OBIGT.R
pkg/CHNOSZ/R/util.args.R
pkg/CHNOSZ/R/util.expression.R
pkg/CHNOSZ/man/logB.to.OBIGT.Rd
pkg/CHNOSZ/vignettes/customizing.Rmd
Log:
logB.to.OBIGT(): use 3 fitting parameters by default
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2022-12-08 09:00:48 UTC (rev 759)
+++ pkg/CHNOSZ/DESCRIPTION 2022-12-09 06:17:33 UTC (rev 760)
@@ -1,6 +1,6 @@
-Date: 2022-12-08
+Date: 2022-12-09
Package: CHNOSZ
-Version: 1.9.9-51
+Version: 1.9.9-52
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/logB.to.OBIGT.R
===================================================================
--- pkg/CHNOSZ/R/logB.to.OBIGT.R 2022-12-08 09:00:48 UTC (rev 759)
+++ pkg/CHNOSZ/R/logB.to.OBIGT.R 2022-12-09 06:17:33 UTC (rev 760)
@@ -3,7 +3,7 @@
# 20220324 v1 Regress three parameters (G, S, and Cp)
# 20221202 v2 Regress HKF parameters (assume constant pressure and optimize omega parameter for charged species)
-logB.to.OBIGT <- function(logB, species, coeffs, T, P, optimize.omega = TRUE, tolerance = 0.05) {
+logB.to.OBIGT <- function(logB, species, coeffs, T, P, npar = 3, optimize.omega = FALSE, tolerance = 0.05) {
# We need at least five temperature values to optimize omega
if(optimize.omega & length(unique(T)) < 5) {
@@ -10,6 +10,11 @@
message("logB.to.OBIGT: forcing optimize.omega = FALSE because there are < 5 T values")
optimize.omega <- FALSE
}
+ # We need five parameters to optimize omega 20221209
+ if(optimize.omega & npar != 5) {
+ message("logB.to.OBIGT: forcing optimize.omega = FALSE because npar != 5")
+ optimize.omega <- FALSE
+ }
## Get the formula of the formed species (used as the species name in OBIGT)
### NOTE: the formed species has to be *last* in this list
@@ -52,9 +57,38 @@
if(!optimize.omega) {
- # Perform linear regression (constant initial omega)
- Glm <- lm(Gf ~ S_var + c1_var + c2_var + omega_var)
- omega <- Glm$coefficients["omega_var"]
+ # Default values
+ omega <- NA
+ c2 <- NA
+ c1 <- NA
+ S <- NA
+ if(npar == 5) {
+ # Perform linear regression with constant omega
+ Glm <- lm(Gf ~ S_var + c1_var + c2_var + omega_var)
+ omega <- Glm$coefficients["omega_var"]
+ c2 <- Glm$coefficients["c2_var"]
+ c1 <- Glm$coefficients["c1_var"]
+ S <- Glm$coefficients["S_var"]
+ G <- Glm$coefficients["(Intercept)"]
+ } else if(npar == 4) {
+ # Now with fewer parameters
+ Glm <- lm(Gf ~ S_var + c1_var + c2_var)
+ c2 <- Glm$coefficients["c2_var"]
+ c1 <- Glm$coefficients["c1_var"]
+ S <- Glm$coefficients["S_var"]
+ G <- Glm$coefficients["(Intercept)"]
+ } else if(npar == 3) {
+ Glm <- lm(Gf ~ S_var + c1_var)
+ c1 <- Glm$coefficients["c1_var"]
+ S <- Glm$coefficients["S_var"]
+ G <- Glm$coefficients["(Intercept)"]
+ } else if(npar == 2) {
+ Glm <- lm(Gf ~ S_var)
+ S <- Glm$coefficients["S_var"]
+ G <- Glm$coefficients["(Intercept)"]
+ } else if(npar == 1) {
+ G <- mean(Gf)
+ } else stop("invalid value for npar (should be from 1 to 5)")
} else {
@@ -73,15 +107,15 @@
omega <- optimize(Sqfun, c(-1e10, 1e10))$minimum
# Construct the linear model with this omega
Glm <- Gfun(omega)
+ # Get coefficients other than omega
+ G <- Glm$coefficients["(Intercept)"]
+ S <- Glm$coefficients["S_var"]
+ c1 <- Glm$coefficients["c1_var"]
+ c2 <- Glm$coefficients["c2_var"]
}
- # Get coefficients
- G <- Glm$coefficients["(Intercept)"]
- S <- Glm$coefficients["S_var"]
- c1 <- Glm$coefficients["c1_var"]
- c2 <- Glm$coefficients["c2_var"]
- # NAs propagate as zero in the EOS
+ # NAs propagate as zero in the HKF equations
if(is.na(S)) S <- 0
if(is.na(c1)) c1 <- 0
if(is.na(c2)) c2 <- 0
Modified: pkg/CHNOSZ/R/util.args.R
===================================================================
--- pkg/CHNOSZ/R/util.args.R 2022-12-08 09:00:48 UTC (rev 759)
+++ pkg/CHNOSZ/R/util.args.R 2022-12-09 06:17:33 UTC (rev 760)
@@ -1,14 +1,14 @@
# CHNOSZ/util.args.R
-# functions to create argument lists and get name of calling function
+# Functions to create argument lists and get name of calling function
-### unexported functions ###
+### Unexported functions ###
-# force T and P to equal length
-# also looks for the keyword Psat in the value of P and substitutes calculated values of the saturation vapor pressure
-TP.args <- function(T=NULL, P=NULL) {
- # keep the [1] here because some functions (e.g. subcrt) will repeat "Psat"
+# Force T and P to equal length
+# Also looks for the keyword Psat in the value of P and substitutes calculated values of the saturation vapor pressure
+TP.args <- function(T = NULL, P = NULL) {
+ # Keep the [1] here because some functions (e.g. subcrt) will repeat "Psat"
if(identical(P[1], "Psat")) {
- P <- water("Psat", T, P="Psat")[, 1]
+ P <- water("Psat", T, P = "Psat")[, 1]
# water.SUPCRT92 issues its own warnings about
# exceeding Psat's temperature limit
if(get("thermo", CHNOSZ)$opt$water == "IAPWS95")
@@ -15,28 +15,28 @@
if(sum(is.na(P))>0)
warning('TP.args: NAs in Psat (likely T > Tc where Tc = 647.096 K)',call.=FALSE)
}
- if(length(P) < length(T) & !is.null(P)) P <- rep(P, length.out=length(T))
- else if(length(T) < length(P) & !is.null(T)) T <- rep(T, length.out=length(P))
- # something we do here so the SUPCRT water calculations work
- T[T==273.15] <- 273.16
- return(list(T=T, P=P))
+ if(length(P) < length(T) & !is.null(P)) P <- rep(P, length.out = length(T))
+ else if(length(T) < length(P) & !is.null(T)) T <- rep(T, length.out = length(P))
+ # Something we do here so the SUPCRT water calculations work
+ T[T == 273.15] <- 273.16
+ return(list(T = T, P = P))
}
-# make the argument lowercase, then transform a, c, g, and l to aq, gas, cr, and liq
-state.args <- function(state=NULL) {
+# Make the argument lowercase, then transform a, c, g, and l to aq, gas, cr, and liq
+state.args <- function(state = NULL) {
if(is.null(state) | is.numeric(state[1])) return(state)
- # normalize state arguments
+ # Normalize state arguments
for(i in 1:length(state)) {
- if(tolower(state[i])=='a') state[i] <- 'aq'
- if(tolower(state[i])=='c') state[i] <- 'cr'
- if(tolower(state[i])=='g') state[i] <- 'gas'
- if(tolower(state[i])=='l') state[i] <- 'liq'
+ if(tolower(state[i]) == 'a') state[i] <- 'aq'
+ if(tolower(state[i]) == 'c') state[i] <- 'cr'
+ if(tolower(state[i]) == 'g') state[i] <- 'gas'
+ if(tolower(state[i]) == 'l') state[i] <- 'liq'
}
return(state)
}
-caller.name <- function(n=2) {
- # returns the name of the calling function n frames up
+caller.name <- function(n = 2) {
+ # Returns the name of the calling function n frames up
# (n=2: the caller of the function that calls this one)
# or character() if called interactively
if(sys.nframe() < n) name <- character()
@@ -43,7 +43,7 @@
else {
sc <- sys.call(-n)[[1]]
name <- tryCatch(as.character(sc), error = function(e) character())
- # also return character() if the value from sys.call is
+ # Also return character() if the value from sys.call is
# the function itself (why does this sometimes happen,
# e.g. when called from affinity()?)
}
Modified: pkg/CHNOSZ/R/util.expression.R
===================================================================
--- pkg/CHNOSZ/R/util.expression.R 2022-12-08 09:00:48 UTC (rev 759)
+++ pkg/CHNOSZ/R/util.expression.R 2022-12-09 06:17:33 UTC (rev 760)
@@ -1,67 +1,67 @@
# CHNOSZ/util.expression.R
-# write descriptions of chemical species, properties, reactions, conditions
-# modified from describe(), axis.label() 20120121 jmd
+# Write descriptions of chemical species, properties, reactions, conditions
+# Modified from describe(), axis.label() 20120121 jmd
-## if this file is interactively sourced, the following are also needed to provide unexported functions:
+## If this file is interactively sourced, the following are also needed to provide unexported functions:
#source("util.character.R")
expr.species <- function(species, state = "aq", value = NULL, log = FALSE, molality = FALSE, use.state = FALSE, use.makeup = FALSE) {
- # make plotting expressions for chemical formulas
+ # Make plotting expressions for chemical formulas
# that include subscripts, superscripts (if charged)
# and optionally designations of states +/- loga or logf prefix
if(length(species) > 1) (stop("more than one species"))
- # convert to character so that "1", "2", etc. don't get converted to chemical formulas via makeup()
+ # Convert to character so that "1", "2", etc. don't get converted to chemical formulas via makeup()
species <- as.character(species)
if(use.makeup) {
- # the counts of elements in the species:
+ # The counts of elements in the species:
# here we don't care too much if an "element" is a real element
# (listed in thermo()$element), so we suppress warnings
elements <- suppressWarnings(try(makeup(species), TRUE))
} else elements <- split.formula(species)
- # if species can't be parsed as a chemical formula, we don't do the formula formatting
+ # If species can't be parsed as a chemical formula, we don't do the formula formatting
if(inherits(elements, "try-error") | !is.numeric(elements)) expr <- species
else {
- # where we'll put the expression
+ # Where we'll put the expression
expr <- ""
- # loop over elements
+ # Loop over elements
for(i in 1:length(elements)) {
if(names(elements)[i] != 'Z') {
- # append the elemental symbol
+ # Append the elemental symbol
expr <- substitute(paste(a, b), list(a=expr, b=names(elements)[i]))
- # recover the coefficient
+ # Recover the coefficient
coeff <- elements[i]
if(coeff!=1) {
- # append the coefficient
+ # Append the coefficient
expr <- substitute(a[b], list(a=expr, b=as.character(coeff)))
}
} else {
- # for charged species, don't show "Z" but do show e.g. "+2"
+ # For charged species, don't show "Z" but do show e.g. "+2"
coeff <- elements[i]
if(coeff==-1) coeff <- "-"
else if(coeff==1) coeff <- "+"
else if(coeff > 0) coeff <- paste("+", as.character(coeff), sep="")
else coeff <- as.character(coeff)
- # append the coefficient as a superscript
+ # Append the coefficient as a superscript
expr <- substitute(a^b, list(a=expr, b=coeff))
}
}
}
- # write the physical state
+ # Write the physical state
if(use.state) {
- # subscript it if we're not giving the value
+ # Subscript it if we're not giving the value
if(is.null(value)) expr <- substitute(a[group('(',italic(b),')')],list(a=expr, b=state))
else expr <- substitute(a*group('(',italic(b),')'),list(a=expr, b=state))
}
- # write a variable and value if given
+ # Write a variable and value if given
if(!is.null(value) | log | molality) {
- # write [logarithm of] activity or fugacity (or molality 20171101)
+ # Write [logarithm of] activity or fugacity (or molality 20171101)
var <- "a"
if(molality) var <- "m"
if(state %in% c("g", "gas")) var <- "f"
expr <- substitute(italic(a)[b], list(a = var, b = expr))
- # use the logarithm?
+ # Use the logarithm?
if(log) expr <- substitute(log ~ a, list(a = expr))
- # write the value if not NULL or NA
+ # Write the value if not NULL or NA
if(!is.null(value)) {
if(!is.na(value)) expr <- substitute(a == b, list(a = expr, b = value))
}
@@ -70,16 +70,16 @@
}
expr.property <- function(property, molality=FALSE) {
- # a way to make expressions for various properties
+ # A way to make expressions for various properties
# e.g. expr.property('DG0r') for standard molal Gibbs
- # energy change of reaction
+ # Energy change of reaction
propchar <- s2c(property)
expr <- ""
- # some special cases
+ # Some special cases
if(is.na(property)) return("")
if(property=="logK") return(quote(log~italic(K)))
if(property=="logB") return(quote(log~beta))
- # grepl here b/c diagram() uses "loga.equil" and "loga.basis"
+ # Use grepl here b/c diagram() uses "loga.equil" and "loga.basis"
if(grepl("loga", property)) {
if(molality) return(quote(log~italic(m)))
else return(quote(log~italic(a)))
@@ -90,12 +90,12 @@
if(property=="pe") return("pe")
if(property=="IS") return(quote(italic(I)))
if(property=="ZC") return(quote(italic(Z)[C]))
- # process each character in the property abbreviation
+ # Process each character in the property abbreviation
prevchar <- character()
for(i in 1:length(propchar)) {
if(i > 1) prevchar <- thischar
thischar <- propchar[i]
- # unless indicated below, uppercase letters are italicized
+ # Unless indicated below, uppercase letters are italicized
# and lowercase letters are italicized and subscripted
# (includes f for property of formation and r for property of reaction)
if(thischar %in% LETTERS) thisexpr <- substitute(italic(a), list(a=thischar))
@@ -111,7 +111,7 @@
if(thischar=='0' & !can.be.numeric(prevchar)) thisexpr <- substitute(degree)
if(thischar=='l') thisexpr <- substitute(a[lambda], list(a=""))
if(thischar=="'") thisexpr <- substitute(minute)
- # put it together
+ # Put it together
expr <- substitute(a*b, list(a=expr, b=thisexpr))
}
return(expr)
@@ -118,8 +118,8 @@
}
expr.units <- function(property, prefix="", per="mol") {
- # make an expression describing units
- # unless we have match below, there will be no units
+ # Make an expression describing units
+ # Unless we have match below, there will be no units
# (e.g., logK, pH, pe)
expr <- ""
# A, G, H - energy
@@ -138,7 +138,7 @@
# T - temperature
if(grepl("T", property)) {
expr <- substitute(a, list(a=T.units()))
- # add a degree sign for Celsius
+ # Add a degree sign for Celsius
if(T.units()=="C") expr <- substitute(degree*a, list(a=expr))
}
# Eh - electrical potential
@@ -146,9 +146,9 @@
# IS - ionic strength
if(grepl("IS", property)) expr <- substitute(a, list(a=mol~kg^-1))
if(!expr=="") {
- # add prefix if appropriate
+ # Add prefix if appropriate
if(!prefix=="") expr <- substitute(a*b, list(a=prefix, b=expr))
- # add mol^-1 if appropriate
+ # Add mol^-1 if appropriate
if(!any(sapply(c("P", "T", "Eh", "IS"), function(x) grepl(x, property))))
expr <- substitute(a~b^-1, list(a=expr, b=per))
}
@@ -155,39 +155,38 @@
return(expr)
}
-axis.label <- function(label, units=NULL, basis=thermo()$basis, prefix="", molality=FALSE) {
- # make a formatted axis label from a generic description
- # it can be a chemical property, condition, or chemical activity in the system
- # if the label matches one of the basis species
- # or if the state is specified, it's a chemical activity
- # 20090826: just return the argument if a comma is already present
+axis.label <- function(label, units = NULL, basis = thermo()$basis, prefix = "", molality = FALSE) {
+ # Make a formatted axis label from a generic description
+ # It can be a chemical property, condition, or chemical activity in the system;
+ # if the label matches one of the basis species or if the state is specified, it's a chemical activity
+ # 20090826: Just return the argument if a comma is already present
# (it's good for custom labels that shouldn't be italicized)
if(grepl(",", label)) return(label)
if(label %in% rownames(basis)) {
- # 20090215: the state this basis species is in
+ # 20090215: The state this basis species is in
state <- basis$state[match(label, rownames(basis))]
- # get the formatted label
+ # Get the formatted label
desc <- expr.species(label, state = state, log = TRUE, molality = molality)
} else {
- # the label is for a chemical property or condition
- # make the label by putting a comma between the property and the units
+ # The label is for a chemical property or condition
+ # Make the label by putting a comma between the property and the units
property <- expr.property(label, molality=molality)
if(is.null(units)) units <- expr.units(label, prefix=prefix)
- # no comma needed if there are no units
+ # No comma needed if there are no units
if(units=="") desc <- substitute(a, list(a=property))
else desc <- substitute(a~"("*b*")", list(a=property, b=units))
}
- # done!
+ # Done!
return(desc)
}
describe.basis <- function(basis = thermo()$basis, ibasis = 1:nrow(basis),
digits = 1, oneline = FALSE, molality = FALSE, use.pH = TRUE) {
- # make expressions for the chemical activities/fugacities of the basis species
+ # Make expressions for the chemical activities/fugacities of the basis species
propexpr <- valexpr <- character()
for(i in ibasis) {
if(can.be.numeric(basis$logact[i])) {
- # we have an as.numeric here in case the basis$logact is character
+ # We have an as.numeric here in case the basis$logact is character
# (by inclusion of a buffer for one of the other basis species)
if(rownames(basis)[i]=="H+" & use.pH) {
propexpr <- c(propexpr, "pH")
@@ -198,7 +197,7 @@
valexpr <- c(valexpr, format(round(as.numeric(basis$logact[i]), digits), nsmall=digits))
}
} else {
- # a non-numeric value is the name of a buffer
+ # A non-numeric value is the name of a buffer
valexpr <- c(valexpr, basis$logact[i])
# propexpr is pH, activity or fugacity
if(rownames(basis)[i]=="H+" & use.pH) propexpr <- c(propexpr, "pH")
@@ -205,7 +204,7 @@
else propexpr <- c(propexpr, expr.species(rownames(basis)[i], state=basis$state[i], value = NA, log = FALSE, molality = molality))
}
}
- # write an equals sign between the property and value
+ # Write an equals sign between the property and value
desc <- character()
for(i in seq_along(propexpr)) {
thisdesc <- substitute(a==b, list(a=propexpr[[i]], b=valexpr[[i]]))
@@ -219,7 +218,7 @@
}
describe.property <- function(property=NULL, value=NULL, digits=0, oneline=FALSE, ret.val=FALSE) {
- # make expressions for pressure, temperature, other conditions
+ # Make expressions for pressure, temperature, other conditions
if(is.null(property) | is.null(value)) stop("property or value is NULL")
propexpr <- valexpr <- character()
for(i in 1:length(property)) {
@@ -233,15 +232,15 @@
}
valexpr <- c(valexpr, as.expression(thisvalexpr))
}
- # with ret.val=TRUE, return just the value with the units (e.g. 55 degC)
+ # With ret.val=TRUE, return just the value with the units (e.g. 55 degC)
if(ret.val) return(valexpr)
- # write an equals sign between the property and value
+ # Write an equals sign between the property and value
desc <- character()
for(i in seq_along(propexpr)) {
if(is.na(value[i])) thisdesc <- propexpr[[i]]
else thisdesc <- substitute(a==b, list(a=propexpr[[i]], b=valexpr[[i]]))
if(oneline) {
- # put all the property/value equations on one line, separated by commas
+ # Put all the property/value equations on one line, separated by commas
if(i==1) desc <- substitute(a, list(a=thisdesc))
else desc <- substitute(list(a, b), list(a=desc, b=thisdesc))
} else desc <- c(desc, thisdesc)
@@ -250,89 +249,89 @@
}
describe.reaction <- function(reaction, iname=numeric(), states=NULL) {
- # make an expression describing the reaction that is
+ # Make an expression describing the reaction that is
# the 'reaction' part of subcrt() output
reactexpr <- prodexpr <- character()
- # loop over the species in the reaction
+ # Loop over the species in the reaction
for(i in 1:nrow(reaction)) {
- # get the name or the chemical formula of the species
+ # Get the name or the chemical formula of the species
if(i %in% iname) species <- reaction$name[i]
else {
- # should the chemical formula have a state?
+ # Should the chemical formula have a state?
if(identical(states,"all") | i %in% states) species <- expr.species(reaction$formula[i], state=reaction$state[i], use.state=TRUE)
else species <- expr.species(reaction$formula[i], state=reaction$state[i])
}
- # get the absolute value of the reaction coefficient
+ # Get the absolute value of the reaction coefficient
abscoeff <- abs(reaction$coeff[i])
- # put the coefficient in if it's not 1
+ # Put the coefficient in if it's not 1
if(abscoeff==1) coeffspec <- species
else {
- # we put in some space if the coefficient comes before a name
+ # We put in some space if the coefficient comes before a name
if(i %in% iname) coeffspec <- substitute(a~b, list(a=abscoeff, b=species))
else coeffspec <- substitute(a*b, list(a=abscoeff, b=species))
}
- # is it a reactant or product?
+ # Is it a reactant or product?
if(reaction$coeff[i] < 0) {
if(length(reactexpr)==0) reactexpr <- substitute(a, list(a=coeffspec))
else reactexpr <- substitute(a+b, list(a=reactexpr, b=coeffspec))
} else {
- # it's a product
+ # It's a product
if(length(prodexpr)==0) prodexpr <- substitute(a, list(a=coeffspec))
else prodexpr <- substitute(a+b, list(a=prodexpr, b=coeffspec))
}
}
- # put an equals sign between reactants and products
- # change this to unicode for the reaction double-arrow 20190218 \u21cc
- # go back to equals - the arrow only works out-of-the-box on linux 20190817
+ # Put an equals sign between reactants and products
+ # Change this to unicode for the reaction double-arrow 20190218 \u21cc
+ # Go back to equals - the arrow doesn't work out-of-the-box on all OSes 20190817
desc <- substitute(a ~ "=" ~ b, list(a=reactexpr, b=prodexpr))
return(desc)
}
-# make formatted text for activity ratio 20170217
-# allow changing the bottom ion 20200716
+# Make formatted text for activity ratio 20170217
+# Allow changing the bottom ion 20200716
ratlab <- function(top = "K+", bottom = "H+", molality = FALSE) {
- # the charges
+ # The charges
Ztop <- makeup(top)["Z"]
Zbottom <- makeup(bottom)["Z"]
- # the text for the exponents
+ # The text for the exponents
exp.bottom <- as.character(Ztop)
exp.top <- as.character(Zbottom)
if(exp.top=="1") exp.top <- ""
if(exp.bottom=="1") exp.bottom <- ""
- # the expression for the top and bottom
+ # The expression for the top and bottom
expr.top <- expr.species(top)
expr.bottom <- expr.species(bottom)
- # with molality, change a to m
+ # With molality, change a to m
a <- ifelse(molality, "m", "a")
- # the final expression
+ # The final expression
substitute(log~(italic(a)[expr.top]^exp.top / italic(a)[expr.bottom]^exp.bottom),
list(a = a, expr.top = expr.top, exp.top = exp.top, expr.bottom = expr.bottom, exp.bottom = exp.bottom))
}
-# make formatted text for thermodynamic system 20170217
+# Make formatted text for thermodynamic system 20170217
syslab <- function(system = c("K2O", "Al2O3", "SiO2", "H2O"), dash="-") {
for(i in seq_along(system)) {
expr <- expr.species(system[i])
- # use en dash here
+ # Use en dash here
if(i==1) lab <- expr else lab <- substitute(a*dash*b, list(a=lab, dash=dash, b=expr))
}
lab
}
-### unexported function ###
+### Unexported function ###
split.formula <- function(formula) {
- ## like makeup(), but split apart the formula based on
+ ## Like makeup(), but split apart the formula based on
## numbers (subscripts); don't scan for elemental symbols 20171018
- # if there are no numbers or charge, return the formula as-is
+ # If there are no numbers or charge, return the formula as-is
# Change [0-9]? to [\\.0-9]* (recognize decimal point in numbers, recognize charge longer than one digit) 20220608
if(! (grepl("[\\.0-9]", formula) | grepl("\\+[\\.0-9]*$", formula) | grepl("-[\\.0-9]*$", formula))) return(formula)
- # first split off charge
+ # First split off charge
# (assume that no subscripts are signed)
Z <- 0
hascharge <- grepl("\\+[\\.0-9]*$", formula) | grepl("-[\\.0-9]*$", formula)
if(hascharge) {
- # for charge, we match + or - followed by zero or more numbers at the end of the string
+ # For charge, we match + or - followed by zero or more numbers at the end of the string
if(grepl("\\+[\\.0-9]*$", formula)) {
fsplit <- strsplit(formula, "+", fixed=TRUE)[[1]]
if(is.na(fsplit[2])) Z <- 1 else Z <- as.numeric(fsplit[2])
@@ -339,7 +338,7 @@
}
if(grepl("-[\\.0-9]*$", formula)) {
fsplit <- strsplit(formula, "-")[[1]]
- # for formula=="H-citrate-2", unsplit H-citrate
+ # For formula=="H-citrate-2", unsplit H-citrate
if(length(fsplit) > 2) {
f2 <- tail(fsplit, 1)
f1 <- paste(head(fsplit, -1), collapse="-")
@@ -349,16 +348,16 @@
}
formula <- fsplit[1]
}
- # to get strings, replace all numbers with placeholder (#), then split on that symbol
- # the outer gsub is to replace multiple #'s with one
+ # To get strings, replace all numbers with placeholder (#), then split on that symbol
+ # The outer gsub is to replace multiple #'s with one
numhash <- gsub("#+", "#", gsub("[\\.0-9]", "#", formula))
strings <- strsplit(numhash, "#")[[1]]
- # to get coefficients, replace all characters (non-numbers) with placeholder, then split
+ # To get coefficients, replace all characters (non-numbers) with placeholder, then split
charhash <- gsub("#+", "#", gsub("[^\\.0-9]", "#", formula))
coeffs <- strsplit(charhash, "#")[[1]]
- # if the first coefficient is empty, remove it
+ # If the first coefficient is empty, remove it
if(coeffs[1]=="") coeffs <- tail(coeffs, -1) else {
- # if the first string is empty, treat the first coefficient as a leading string (e.g. in 2-octanone)
+ # If the first string is empty, treat the first coefficient as a leading string (e.g. in 2-octanone)
if(strings[1]=="") {
strings[2] <- paste0(coeffs[1], strings[2])
coeffs <- tail(coeffs, -1)
@@ -365,14 +364,14 @@
strings <- tail(strings, -1)
}
}
- # if we're left with no coefficients, just return the string
+ # If we're left with no coefficients, just return the string
if(length(coeffs)==0 & Z==0) return(strings)
- # if we're missing a coefficient, append one
+ # If we're missing a coefficient, append one
if(length(coeffs) < length(strings)) coeffs <- c(coeffs, 1)
- # use strings as names for the numeric coefficients
+ # Use strings as names for the numeric coefficients
coeffs <- as.numeric(coeffs)
names(coeffs) <- strings
- # include charge if it is not 0
+ # Include charge if it is not 0
if(Z!=0) coeffs <- c(coeffs, Z=Z)
return(coeffs)
}
Modified: pkg/CHNOSZ/man/logB.to.OBIGT.Rd
===================================================================
--- pkg/CHNOSZ/man/logB.to.OBIGT.Rd 2022-12-08 09:00:48 UTC (rev 759)
+++ pkg/CHNOSZ/man/logB.to.OBIGT.Rd 2022-12-09 06:17:33 UTC (rev 760)
@@ -7,7 +7,8 @@
}
\usage{
- logB.to.OBIGT(logB, species, coeffs, T, P, optimize.omega = TRUE, tolerance = 0.05)
+ logB.to.OBIGT(logB, species, coeffs, T, P, npar = 3,
+ optimize.omega = FALSE, tolerance = 0.05)
}
\arguments{
@@ -16,6 +17,7 @@
\item{coeffs}{Coefficients of the formation reaction}
\item{T}{Temperature (units of \code{\link{T.units}})}
\item{P}{Temperature (\samp{Psat} or numerical values with units of \code{\link{P.units}})}
+ \item{npar}{numeric, number of parameters to use in fitting}
\item{optimize.omega}{logical, optimize the omega parameter (relevant for charged species)?}
\item{tolerance}{Tolerance for checking equivalence of input and calculated \logB values}
}
@@ -22,27 +24,31 @@
\details{
-This function can be used to update the \code{\link{OBIGT}} thermodynamic database with parameters fit to formation constants of aqueous species as a function of temperature.
-The \code{logB} argument should have decimal logarithm of equilibrium constants (\logB) for a formation reaction.
+This function updates the \code{\link{OBIGT}} thermodynamic database with parameters fit to formation constants of aqueous species as a function of temperature.
+The \code{logB} argument should have decimal logarithm of formation constants for an aqueous complex (\logB).
The formation reaction is defined by \code{species} and \code{coeffs}.
-All species in the formation reaction must be available in OBIGT except for the \emph{last} species, which is the newly formed species.
+All species in the formation reaction must be present in OBIGT except for the \emph{last} species, which is the newly formed species.
The function works as follows.
First, the properties of the known species are combined with \logB to calculate the standard Gibbs energy (G[T]) of the formed species at each value of \code{T} and \code{P}.
-Then, the values of G[T] are used to retrieve the thermodynamic parameters \code{G}, \code{S}, \code{c1}, \code{c2}, and \code{omega}.
-The first two of these are standard-state values at 25 \degC and 1 bar, and the last three are parameters in the revised HKF equations (e.g. Eq. B25 of Shock et al., 1992).
-Volumetric parameters (\a1 to \a4) in the HKF equations currently aren't included, so the resulting fits are valid only at the pressure corresponding to the given values of \logB.
+Then, \code{\link{lm}} is used to fit one or more of the thermodynamic parameters \code{G}, \code{S}, \code{c1}, \code{c2}, and \code{omega} to the values of G[T].
+The first two of these parameters are standard-state values at 25 \degC and 1 bar, and the last three are parameters in the revised HKF equations (e.g. Eq. B25 of Shock et al., 1992).
The fitted parameters for the formed species are then added to OBIGT.
-\code{\link{all.equal}} is used to test for approximate equivalence of the input and calculated equilibrium constants; if the mean absolute difference exceeds \code{tolerance}, an error occurs.
+Finally, \code{\link{all.equal}} is used to test for approximate equivalence of the input values of \logB and calculated equilibrium constants; if the mean absolute difference exceeds \code{tolerance}, an error occurs.
-The number of parameters used in the fit is not more than the number of data values (\emph{n}), and the parameters are added in the order stated above.
+To avoid overfitting, only the first three parameters (code{G}, \code{S}, and \code{c1}) are used by default.
+More parameters (up to 5) or fewer (down to 1) can be selected by changing \code{npar}.
+Volumetric parameters (\a1 to \a4) in the HKF equations currently aren't included, so the resulting fits are valid only at the input pressure values.
+
+Independent of \code{npar}, the number of parameters used in the fit is not more than the number of data values (\emph{n}).
If \emph{n} is less than 5, then the values of the unused parameters are set to 0 for addition to OBIGT.
For instance, a single value of \logB would be fit only with \code{G}, implying that computed values of G[T] have no temperature dependence.
-For uncharged species, the value of \omega is a constant, but for charged species, it is a function of \T and \P as described by the \dQuote{g function} (Shock et al., 1992).
-By default an optimization step is used to refine the value of \code{omega} at 25 \degC and 1 bar for charged species.
-In case the data are not sufficient to constrain \code{omega} (see first example), this step can be skipped by setting \code{optimize.omega} to FALSE; then, it is assumed that \omega is constant during the fitting.
-In OBIGT, a constant omega for a charged species can be enforced by setting the \code{z} parameter to 0, but that is not done by this function.
+The value of \omega is a constant in the revised HKF equations for uncharged species, but for charged species, it is a function of \T and \P as described by the \dQuote{g function} (Shock et al., 1992).
+An optimization step is available to refine the value of \code{omega} at 25 \degC and 1 bar for charged species taking its temperature dependence into account for the fitting.
+However, representative datasets are not well represented by variable \code{omega} (see first example), so this step is skipped by default.
+Furthermore, \code{logB.to.OBIGT} by default also sets the \code{z} parameter in OBIGT to 0; this enforces a constant \omega for charged species in calculations with \code{\link{subcrt}} (note that this is a parameter for the HKF equations and does not affect reaction balancing).
+Set \code{optimize.omega} to TRUE to override the defaults and activate variable \omega for charged species; this takes effect only if \code{npar == 5} and \emph{n} > 5.
}
@@ -65,7 +71,7 @@
opar <- par(mfrow = c(2, 1))
for(optimize.omega in c(TRUE, FALSE)) {
# Fit the parameters with or without variable omega
- inew <- logB.to.OBIGT(logB, species, coeffs, T, P, optimize.omega = optimize.omega)
+ inew <- logB.to.OBIGT(logB, species, coeffs, T, P, npar = 5, optimize.omega = optimize.omega)
# Print the new database entry
info(inew)
# Plot experimental logB
@@ -76,10 +82,11 @@
sres <- subcrt(species, coeffs, T = Tfit)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 760
More information about the CHNOSZ-commits
mailing list