[CHNOSZ-commits] r967 - in pkg/CHNOSZ: . R demo inst inst/extdata/Berman inst/extdata/OBIGT inst/extdata/OBIGT/testing inst/extdata/misc inst/extdata/thermo inst/tinytest man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Jun 21 10:52:52 CEST 2026
Author: jedick
Date: 2026-06-21 10:52:51 +0200 (Sun, 21 Jun 2026)
New Revision: 967
Added:
pkg/CHNOSZ/R/JANAF.to.OBIGT.R
pkg/CHNOSZ/demo/chalcocite.R
pkg/CHNOSZ/demo/gold1.R
pkg/CHNOSZ/demo/gold2.R
pkg/CHNOSZ/demo/neodymium.R
pkg/CHNOSZ/demo/sulfur.R
pkg/CHNOSZ/inst/extdata/OBIGT/SLOP98-a.csv
pkg/CHNOSZ/inst/extdata/OBIGT/SLOP98-b.csv
pkg/CHNOSZ/inst/extdata/OBIGT/inorganic_liq.csv
pkg/CHNOSZ/inst/extdata/misc/C-127.txt
pkg/CHNOSZ/inst/extdata/misc/Fe-001.txt
pkg/CHNOSZ/inst/extdata/misc/S-004.txt
pkg/CHNOSZ/inst/tinytest/test-JANAF.to.OBIGT.R
pkg/CHNOSZ/inst/tinytest/test-S_liq.R
pkg/CHNOSZ/man/JANAF.to.OBIGT.Rd
Removed:
pkg/CHNOSZ/demo/gold.R
pkg/CHNOSZ/demo/sum_S.R
pkg/CHNOSZ/inst/TODO
pkg/CHNOSZ/inst/extdata/OBIGT/SLOP98.csv
pkg/CHNOSZ/inst/extdata/OBIGT/testing/IGEM.csv
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/NAMESPACE
pkg/CHNOSZ/R/add.OBIGT.R
pkg/CHNOSZ/R/affinity.R
pkg/CHNOSZ/R/diagram.R
pkg/CHNOSZ/R/equilibrate.R
pkg/CHNOSZ/R/info.R
pkg/CHNOSZ/R/phosphorylate.R
pkg/CHNOSZ/R/solubility.R
pkg/CHNOSZ/R/subcrt.R
pkg/CHNOSZ/R/thermo.R
pkg/CHNOSZ/R/util.affinity.R
pkg/CHNOSZ/R/util.data.R
pkg/CHNOSZ/R/util.legend.R
pkg/CHNOSZ/R/util.misc.R
pkg/CHNOSZ/demo/00Index
pkg/CHNOSZ/demo/Pourbaix.R
pkg/CHNOSZ/demo/aluminum.R
pkg/CHNOSZ/demo/arsenic.R
pkg/CHNOSZ/demo/buffer.R
pkg/CHNOSZ/demo/contour.R
pkg/CHNOSZ/demo/demos.R
pkg/CHNOSZ/demo/demos_png.R
pkg/CHNOSZ/demo/glycinate.R
pkg/CHNOSZ/demo/minsol.R
pkg/CHNOSZ/demo/phosphorylate.R
pkg/CHNOSZ/demo/potassium.R
pkg/CHNOSZ/demo/references.R
pkg/CHNOSZ/demo/saturation.R
pkg/CHNOSZ/demo/uranyl.R
pkg/CHNOSZ/inst/NEWS.Rd
pkg/CHNOSZ/inst/extdata/Berman/sympy.R
pkg/CHNOSZ/inst/extdata/OBIGT/inorganic_aq.csv
pkg/CHNOSZ/inst/extdata/OBIGT/inorganic_cr.csv
pkg/CHNOSZ/inst/extdata/OBIGT/organic_aq.csv
pkg/CHNOSZ/inst/extdata/OBIGT/organic_gas.csv
pkg/CHNOSZ/inst/extdata/misc/OBIGT_check.csv
pkg/CHNOSZ/inst/extdata/thermo/refs.csv
pkg/CHNOSZ/inst/extdata/thermo/stoich.csv.xz
pkg/CHNOSZ/inst/tinytest/test-AD.R
pkg/CHNOSZ/inst/tinytest/test-affinity.R
pkg/CHNOSZ/inst/tinytest/test-diagram.R
pkg/CHNOSZ/inst/tinytest/test-equilibrate.R
pkg/CHNOSZ/inst/tinytest/test-mix.R
pkg/CHNOSZ/inst/tinytest/test-nonideal.R
pkg/CHNOSZ/inst/tinytest/test-phosphorylate.R
pkg/CHNOSZ/inst/tinytest/test-rank.affinity.R
pkg/CHNOSZ/inst/tinytest/test-subcrt.R
pkg/CHNOSZ/inst/tinytest/test-util.legend.R
pkg/CHNOSZ/man/Berman.Rd
pkg/CHNOSZ/man/affinity.Rd
pkg/CHNOSZ/man/diagram.Rd
pkg/CHNOSZ/man/examples.Rd
pkg/CHNOSZ/man/extdata.Rd
pkg/CHNOSZ/man/mod.buffer.Rd
pkg/CHNOSZ/man/mosaic.Rd
pkg/CHNOSZ/man/nonideal.Rd
pkg/CHNOSZ/man/phosphorylate.Rd
pkg/CHNOSZ/man/solubility.Rd
pkg/CHNOSZ/man/subcrt.Rd
pkg/CHNOSZ/man/thermo.Rd
pkg/CHNOSZ/man/water.Rd
pkg/CHNOSZ/vignettes/FAQ.Rmd
pkg/CHNOSZ/vignettes/OBIGT.Rmd
pkg/CHNOSZ/vignettes/OBIGT.bib
pkg/CHNOSZ/vignettes/anintro.Rmd
pkg/CHNOSZ/vignettes/multi-metal.Rmd
pkg/CHNOSZ/vignettes/vig.bib
Log:
Merge 23 GitHub commits (e231bb2 to 14c98a0)
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2026-06-02 12:33:23 UTC (rev 966)
+++ pkg/CHNOSZ/DESCRIPTION 2026-06-21 08:52:51 UTC (rev 967)
@@ -1,6 +1,6 @@
-Date: 2026-06-02
+Date: 2026-06-21
Package: CHNOSZ
-Version: 2.2.0-33
+Version: 2.2.0-55
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/NAMESPACE
===================================================================
--- pkg/CHNOSZ/NAMESPACE 2026-06-02 12:33:23 UTC (rev 966)
+++ pkg/CHNOSZ/NAMESPACE 2026-06-21 08:52:51 UTC (rev 967)
@@ -62,7 +62,9 @@
# added 20220620
"stack_mosaic",
# added 20251207
- "phosphorylate", "phospho_plot", "hyphen.in.pdf"
+ "phosphorylate", "phospho.plot", "hyphen.in.pdf",
+# added 20260605
+ "JANAF.to.OBIGT"
)
# Load shared objects
Added: pkg/CHNOSZ/R/JANAF.to.OBIGT.R
===================================================================
--- pkg/CHNOSZ/R/JANAF.to.OBIGT.R (rev 0)
+++ pkg/CHNOSZ/R/JANAF.to.OBIGT.R 2026-06-21 08:52:51 UTC (rev 967)
@@ -0,0 +1,122 @@
+# JANAF.to.OBIGT.R
+# Convert JANAF data to OBIGT format
+# 20260604 jmd
+
+# The first four lines of https://janaf.nist.gov/tables/C-127.txt
+# Ethyne (C2H2) C2H2(g)
+# T(K) Cp S -[G-H(Tr)]/T H-H(Tr) delta-f H delta-f G log Kf
+# 0 0.000 0.000 INFINITE -10.012 227.288 227.288 INFINITE
+# 100 29.347 163.294 234.338 -7.104 227.078 221.011 -115.444
+
+JANAF.to.OBIGT <- function(file, abbrv = NULL, T_max = 1500, MAE_max = 1, plot_Cp = FALSE) {
+
+ # Preprocess file: replace strings to get equal numbers of columns for read.table()
+ lines <- readLines(file)
+ # For S1(l) and C2H2(g)
+ lines <- gsub("Cp LAMBDA MAXIMUM", "NA\tNA\tNA", lines)
+ lines <- gsub("TRANSITION", "NA\tNA\tNA", lines)
+ # For S1(l)
+ lines <- gsub("\t\t\t\t\t\t\t\t\t\t\t\t\t\t", "\tNA\tNA\tNA\tNA\tNA\tNA\tNA", lines)
+ lines <- gsub("BETA <--> LIQUID", "NA\tNA\tNA", lines)
+ lines <- gsub("FUGACITY = 1 bar", "NA\tNA\tNA", lines)
+ tmpfile <- tempfile(fileext=".txt")
+ writeLines(lines, tmpfile)
+
+ # Read the table, skipping the first three lines
+ dat <- read.table(tmpfile, skip = 3)
+
+ # Process abbrv argument
+ abbrv = ifelse(is.null(abbrv), NA_character_, abbrv)
+
+ # Add column names
+ colnames(dat) <- c("T", "Cp", "S", "DG", "DH", "H", "G", "logK")
+ # Read the first line to get name, formula, and state
+ first <- readLines(tmpfile, n = 1)
+ name <- tolower(strsplit(first, " ")[[1]][1])
+ # Get formula with state, e.g. C2H2(g)
+ formula_state <- strsplit(first, "\t")[[1]][2]
+ # Get formula by itself, e.g. C2H2
+ formula <- strsplit(formula_state, "\\(")[[1]][1]
+ # Get state (inside parentheses)
+ m <- regexec("\\(([^)]*)\\)", formula_state)
+ state_orig <- regmatches(formula_state, m)[[1]][2]
+ # Convert g to gas
+ # TODO: work with other states (aq, cr)
+ state <- switch(state_orig, g = "gas", l = "liq", NA)
+ if(is.na(state)) stop(paste("unrecognized state:", state_orig))
+ # Format today's date for ISO 8601, e.g. 2026-06-04
+ date <- format(Sys.time(), "%Y-%m-%d")
+ # Find row for 298.15 K
+ i298 <- dat$T==298.15
+
+ # Find rows to use for Cp values (100 to 1500 K)
+ iCp <- dat$T >= 100 & dat$T <= T_max
+ Cp <- dat$Cp[iCp]
+ # Remove NA Cp values
+ iNA <- is.na(Cp)
+ Cp <- Cp[!iNA]
+ iCp[iNA] <- FALSE
+ # Use linear model to fit Cp equation
+ # Cp = a + b*T + c*T^-2 + d*T^-0.5 + e*T^2
+ T <- dat$T[iCp]
+ T_2 <- T ^ -2
+ T_0.5 <- T ^ -0.5
+ T2 <- T ^ 2
+ Cp_lm <- lm(Cp ~ T + T_2 + T_0.5 + T2)
+ Cp_predicted <- predict(Cp_lm)
+
+ if(plot_Cp) {
+ # Plot actual values
+ plot(T, Cp, xlab = quote(italic(T)~"(K)"), ylab = axis.label("Cp"))
+ # Predict Cp at closely-spaced T values
+ T_pred <- min(T):max(T)
+ newdata <- data.frame(
+ T = T_pred,
+ T_2 = T_pred ^ -2,
+ T_0.5 = T_pred ^ -0.5,
+ T2 = T_pred ^ 2
+ )
+ Cp_pred <- predict(Cp_lm, newdata)
+ lines(T_pred, Cp_pred)
+ legend("topleft", c("actual", "predicted"), pch = c(1, NA), lty = c(0, 1))
+ title(formula_state, font.main = 1)
+ }
+
+ # Calculate MAE
+ MAE <- mean(abs(Cp_predicted - Cp))
+ message(paste0("JANAF.to.OBIGT: MAE for Cp of ", name, " ", formula_state, ": ", round(MAE, 2)))
+ # Error if MAE > max MAE
+ stopifnot(MAE <= MAE_max)
+
+ # Put together parameters
+ PAR <- list(
+ name = name,
+ abbrv = abbrv,
+ formula = formula,
+ state = state,
+ ref1 = "JANAF98",
+ ref2 = "OBIGT26",
+ date = date,
+ model = "CGL",
+ E_units = "J",
+ G = dat$G[i298] * 1000,
+ H = dat$H[i298] * 1000,
+ S = dat$S[i298],
+ Cp = dat$Cp[i298],
+ V = 0,
+ a = signif(Cp_lm$coefficients[1], 5),
+ b = signif(Cp_lm$coefficients[2], 5),
+ c = signif(Cp_lm$coefficients[3], 5),
+ d = signif(Cp_lm$coefficients[4], 5),
+ e = signif(Cp_lm$coefficients[5], 5),
+ f = 0,
+ lambda = 0,
+ T = 1500,
+ zap = TRUE
+ )
+
+ # Add data to OBIGT
+ inew <- do.call(mod.OBIGT, PAR)
+ return(inew)
+
+}
Modified: pkg/CHNOSZ/R/add.OBIGT.R
===================================================================
--- pkg/CHNOSZ/R/add.OBIGT.R 2026-06-02 12:33:23 UTC (rev 966)
+++ pkg/CHNOSZ/R/add.OBIGT.R 2026-06-21 08:52:51 UTC (rev 967)
@@ -113,10 +113,12 @@
# Add/replace entries in thermo$OBIGT from values saved in a file
# Only replace if force == TRUE
thermo <- get("thermo", CHNOSZ)
- to1 <- thermo$OBIGT
- id1 <- paste(to1$name,to1$state)
+ olddat <- thermo$OBIGT
+ oldid <- paste(olddat$name, olddat$state)
- # `file` should be the path to a CSV file, or the name (without path or .csv) of a file in the package's OBIGT or OBIGT/testing
+ # `file` should be the path to a CSV file, or the base name (before hyphen,
+ # without path or .csv) of a file in the package's OBIGT or OBIGT/testing
+ # Example: file = "SLOP98" reads both SLOP98-a.csv and SLOP98-b.csv 20260617
if(!file.exists(file)) {
# List all files in OBIGT and OBIGT/testing
@@ -125,28 +127,42 @@
all_files <- c(OBIGT_files, testing_files)
# Match argument to file names without path or .csv suffix
all_names <- gsub(".csv", "", basename(all_files))
- ifile <- match(file, all_names)
- if(is.na(ifile)) stop(paste(file, "is not a file and doesn't match any files in the OBIGT database"))
- file <- all_files[ifile]
+ # Split file names at hyphen 20260617
+ all_names <- sapply(strsplit(all_names, "-"), "[", 1)
+ # Get all matches
+ matches <- which(all_names %in% file)
+ if(length(matches) == 0) {
+ stop(paste(file, "is not a file and doesn't match any files in the OBIGT database"))
+ } else if(length(matches) == 1) {
+ file <- all_files[matches]
+ } else {
+ # Read all matching files, combine them,
+ # and put into a single temporary file for further processing 20260617
+ file_contents <- lapply(all_files[matches], read.csv)
+ combined_contents <- do.call(rbind, file_contents)
+ file <- tempfile(fileext = ".csv")
+ write.csv(combined_contents, file, row.names = FALSE)
+ }
}
# Read data from the file
- to2 <- read.csv(file, as.is = TRUE)
- Etxt <- paste(unique(to2$E_units), collapse = " and ")
+ newdat <- read.csv(file, as.is = TRUE)
+ Etxt <- paste(unique(newdat$E_units), collapse = " and ")
# Load only selected species if requested
if(!is.null(species)) {
- idat <- match(species, to2$name)
+ idat <- match(species, newdat$name)
ina <- is.na(idat)
- if(!any(ina)) to2 <- to2[idat, ]
+ if(!any(ina)) newdat <- newdat[idat, ]
else stop(paste("file", file, "doesn't have", paste(species[ina], collapse = ", ")))
}
- id2 <- paste(to2$name,to2$state)
+ newid <- paste(newdat$name, newdat$state)
# Check if the data table is compatible with thermo$OBIGT
- if(!identical(colnames(to1), colnames(to2))) stop(paste(file, "does not have same column names as thermo$OBIGT data frame."))
+ if(!identical(colnames(olddat), colnames(newdat)))
+ stop(paste(file, "does not have same column names as thermo$OBIGT data frame."))
# Match the new species to existing ones
- does.exist <- id2 %in% id1
- ispecies.exist <- na.omit(match(id2, id1))
+ does.exist <- newid %in% oldid
+ ispecies.exist <- na.omit(match(newid, oldid))
nexist <- sum(does.exist)
# Keep track of the species we've added
inew <- numeric()
@@ -153,27 +169,28 @@
if(force) {
# Replace existing entries
if(nexist > 0) {
- to1[ispecies.exist, ] <- to2[does.exist, ]
- to2 <- to2[!does.exist, ]
+ olddat[ispecies.exist, ] <- newdat[does.exist, ]
+ newdat <- newdat[!does.exist, ]
inew <- c(inew, ispecies.exist)
}
} else {
# Ignore any new entries that already exist
- to2 <- to2[!does.exist, ]
+ newdat <- newdat[!does.exist, ]
nexist <- 0
}
# Add new entries
- if(nrow(to2) > 0) {
- to1 <- rbind(to1, to2)
- inew <- c(inew, (length(id1)+1):nrow(to1))
+ alldat <- olddat
+ if(nrow(newdat) > 0) {
+ alldat <- rbind(olddat, newdat)
+ inew <- c(inew, (length(oldid)+1):nrow(alldat))
}
# Commit the change
- thermo$OBIGT <- to1
+ thermo$OBIGT <- alldat
rownames(thermo$OBIGT) <- 1:nrow(thermo$OBIGT)
assign("thermo", thermo, CHNOSZ)
# Give the user a message
message("add.OBIGT: read ", length(does.exist), " rows; made ",
- nexist, " replacements, ", nrow(to2), " additions [energy units: ", Etxt, "]")
+ nexist, " replacements, ", nrow(newdat), " additions [energy units: ", Etxt, "]")
#message("add.OBIGT: use OBIGT() or reset() to restore default database")
return(invisible(inew))
}
Modified: pkg/CHNOSZ/R/affinity.R
===================================================================
--- pkg/CHNOSZ/R/affinity.R 2026-06-02 12:33:23 UTC (rev 966)
+++ pkg/CHNOSZ/R/affinity.R 2026-06-21 08:52:51 UTC (rev 967)
@@ -7,15 +7,16 @@
#source("util.character.R")
#source("util.list.R")
#source("subcrt.R")
-#source("buffer.R")
+#source("mod.buffer.R")
#source("util.args.R")
#source("util.data.R")
+#source("util.misc.R")
#source("species.R")
#source("info.R")
#source("hkf.R")
#source("cgl.R")
-affinity <- function(..., property = NULL, sout = NULL, exceed.Ttr = FALSE, exceed.rhomin = FALSE,
+affinity <- function(..., property = NULL, sout = NULL, warn.Ttr = TRUE, exceed.rhomin = FALSE,
return.buffer = FALSE, return.sout = FALSE, balance = "PBB", iprotein = NULL, loga.protein = 0, transect = NULL) {
# ...: variables over which to calculate
# property: what type of energy
@@ -22,8 +23,8 @@
# (G.basis, G.species, logact.basis, logK, logQ, A)
# return.buffer: return buffered activities
# balance: balance protein buffers on PBB
- # exceed.Ttr: extrapolate Gibbs energies
- # of minerals beyond their T-limits?
+ # warn.Ttr: warn if Gibbs energies are extrapolated
+ # above transition temperature for minerals?
# sout: provide a previously calculated output from subcrt
# iprotein: build these proteins from residues (speed optimization)
@@ -52,7 +53,7 @@
# The argument list
args <- energy.args(args.orig, transect = transect)
- args <- c(args, list(sout = sout, exceed.Ttr = exceed.Ttr, exceed.rhomin = exceed.rhomin))
+ args <- c(args, list(sout = sout, warn.Ttr = warn.Ttr, exceed.rhomin = exceed.rhomin))
# The user-defined species (including basis species, formed species, and proteins)
thermo <- get("thermo", CHNOSZ)
@@ -252,7 +253,7 @@
# Content of return value depends on buffer request
if(return.buffer) return(c(tb, list(vars = vars, vals = vals)))
# For argument recall, include all arguments (except sout) in output 20190117
- allargs <- c(args.orig, list(property = property, exceed.Ttr = exceed.Ttr, exceed.rhomin = exceed.rhomin,
+ allargs <- c(args.orig, list(property = property, warn.Ttr = warn.Ttr, exceed.rhomin = exceed.rhomin,
return.buffer = return.buffer, balance = balance, iprotein = iprotein, loga.protein = loga.protein))
# Add IS value only if it given as an argument 20171101
# (even if its value is 0, the presence of IS will trigger diagram() to use "m" instead of "a" in axis labels)
Modified: pkg/CHNOSZ/R/diagram.R
===================================================================
--- pkg/CHNOSZ/R/diagram.R 2026-06-02 12:33:23 UTC (rev 966)
+++ pkg/CHNOSZ/R/diagram.R 2026-06-21 08:52:51 UTC (rev 967)
@@ -48,7 +48,9 @@
if(!(efun %in% c("affinity", "rank.affinity", "equilibrate") | grepl("solubilit", efun)))
stop("'eout' is not the output from one of these functions: affinity, rank.affinity, equilibrate, or solubility")
# For solubility(), default type is loga.balance 20210303
- if(grepl("solubilit", efun) & missing(type)) type <- "loga.balance"
+ if(grepl("solubility", efun) & missing(type)) type <- "loga.balance"
+ # Use loga.equil to plot stabilities of multiple substrates 20260614
+ if(grepl("solubilities", efun) & missing(type)) type <- "loga.equil"
# Check balance argument for rank.affinity() 20220416
if(efun == "rank.affinity") {
if(!identical(balance, 1)) {
@@ -266,7 +268,7 @@
}
else names <- as.character(eout$species$name)
# Remove non-unique organism or protein names
- if(all(grepl("_", names))) {
+ if(all(is.protein(names))) {
is.pname <- TRUE
# Everything before the underscore (the protein)
pname <- gsub("_.*$", "", names)
@@ -290,7 +292,7 @@
}
## Apply formatting to chemical formulas 20170204
- if(all(grepl("_", names))) is.pname <- TRUE
+ if(all(is.protein(names))) is.pname <- TRUE
if(format.names & !is.pname) {
# Check if names are a deparsed expression (used in mix()) 20200718
parsed <- FALSE
Modified: pkg/CHNOSZ/R/equilibrate.R
===================================================================
--- pkg/CHNOSZ/R/equilibrate.R 2026-06-02 12:33:23 UTC (rev 966)
+++ pkg/CHNOSZ/R/equilibrate.R 2026-06-21 08:52:51 UTC (rev 967)
@@ -85,7 +85,7 @@
}
## Normalize the molar formula by the balance coefficients
m.balance <- n.balance
- isprotein <- grepl("_", as.character(aout$species$name))
+ isprotein <- is.protein(as.character(aout$species$name))
if(any(normalize) | as.residue) {
if(any(n.balance < 0)) stop("one or more negative balancing coefficients prohibit using normalized molar formulas")
n.balance[normalize|as.residue] <- 1
@@ -363,7 +363,7 @@
# The index of the basis species that might be balanced
ibalance <- numeric()
# Deal with proteins
- isprotein <- grepl("_", as.character(aout$species$name))
+ isprotein <- is.protein(as.character(aout$species$name))
if(is.null(balance) & all(isprotein)) balance <- "length"
# Try to automatically find a balance
if(is.null(balance)) {
Modified: pkg/CHNOSZ/R/info.R
===================================================================
--- pkg/CHNOSZ/R/info.R 2026-06-02 12:33:23 UTC (rev 966)
+++ pkg/CHNOSZ/R/info.R 2026-06-21 08:52:51 UTC (rev 967)
@@ -8,6 +8,7 @@
## If this file is interactively sourced, the following are also needed to provide unexported functions:
#source("util.data.R")
+#source("util.misc.R")
info <- function(species = NULL, state = NULL, check.it = TRUE) {
@@ -56,7 +57,7 @@
# First look for exact match
ispecies <- info.character(species[i], state[i])
# If no exact match and it's not a protein, show approximate matches (side effect of info.approx)
- if(identical(ispecies, NA) & !grepl("_", species[i])) ispecies.notused <- info.approx(species[i], state[i])
+ if(identical(ispecies, NA) & !is.protein(species[i])) ispecies.notused <- info.approx(species[i], state[i])
# Do not accept multiple matches
if(length(ispecies) > 1) ispecies <- NA
return(ispecies)
@@ -286,10 +287,11 @@
}
# If we got here there were no approximate matches
# 20190127 look for the species in optional data files
- for(opt in c("SLOP98", "SUPCRT92", "AD")) {
+ for(opt in c("SLOP98-a", "SLOP98-b", "SUPCRT92", "AD")) {
optdat <- read.csv(system.file(paste0("extdata/OBIGT/", opt, ".csv"), package = "CHNOSZ"), as.is = TRUE)
if(species %in% optdat$name) {
- message('info.approx: ', species, ' is in an optional database; use add.OBIGT("', opt, '", "', species, '") to load it')
+ message('info.approx: ', species, ' is in an optional database; use add.OBIGT("',
+ strsplit(opt, "-")[[1]][1], '", "', species, '") to load it')
return(NA)
}
}
Modified: pkg/CHNOSZ/R/phosphorylate.R
===================================================================
--- pkg/CHNOSZ/R/phosphorylate.R 2026-06-02 12:33:23 UTC (rev 966)
+++ pkg/CHNOSZ/R/phosphorylate.R 2026-06-21 08:52:51 UTC (rev 967)
@@ -1,6 +1,6 @@
# Calculate affinity of phosphorylation reactions taking account of speciation
# 20251206 first version (extracted from sugars paper script) jmd
-# 20251208 add phospho_plot()
+# 20251208 add phospho.plot()
# 20251224 add Mg species
phosphorylate <- function(reactant, P_source, loga_reactant = 0, loga_product = 0, loga_P_source = 0, loga_P_remainder = 0, const_pH = 7, loga_Mg = -999, ...) {
@@ -220,7 +220,7 @@
} # end of phosphorylate()
# Define a function to make the plots for a given reaction
-phospho_plot <- function(reactant, P_source, loga_Mg = -999, res = 50) {
+phospho.plot <- function(reactant, P_source, loga_Mg = -999, res = 50) {
# Reaction-independent settings (activities of species)
# The product (phosphorylated species)
@@ -410,5 +410,5 @@
# Return the calculated value
G.kJ
-} # end of phospho_plot()
+} # end of phospho.plot()
Modified: pkg/CHNOSZ/R/solubility.R
===================================================================
--- pkg/CHNOSZ/R/solubility.R 2026-06-02 12:33:23 UTC (rev 966)
+++ pkg/CHNOSZ/R/solubility.R 2026-06-21 08:52:51 UTC (rev 967)
@@ -2,7 +2,7 @@
# 20181031 first version jmd
# 20181106 work on the output from affinity(); no "equilibrate()" needed!
# 20190117 add find.IS and test for dissociation reaction
-# 20210319 use vector of aqueous species as main argument (with back-compatibility for affinity output) and handle multiple minerals
+# 20210319 use vector of aqueous species as main argument (with back-compatibility for affinity output) and handle multiple substrates
## If this file is interactively sourced, the following are also needed to provide unexported functions:
#source("equilibrate.R")
@@ -13,9 +13,9 @@
#source("mosaic.R")
#source("basis.R")
-# Function to calculate solubilities of multiple minerals 20210303
-# species() should be used first to load the minerals (all bearing the same metal)
-# 'iaq' lists aqueous species that can be produced by dissolution of the minerals
+# Function to calculate solubilities of multiple substrates 20210303
+# species() should be used first to load the substrates (all bearing the same metal)
+# 'iaq' lists aqueous species that can be produced by dissolution of the substrates
# '...' contains arguments for affinity() or mosaic() (i.e. plotting variables)
solubility <- function(iaq, ..., in.terms.of = NULL, dissociate = FALSE, find.IS = FALSE) {
@@ -36,13 +36,13 @@
# Use current basis species as a template for the solubility calculations
ispecies <- basis()$ispecies
logact <- basis()$logact
- # The current formed species are the minerals to be dissolved
- mineral <- species()
- if(is.null(mineral)) stop("please load minerals or gases with species()")
+ # The current formed species are the substrates to be dissolved
+ substrates <- species()
+ if(is.null(substrates)) stop("please load substrates (e.g. minerals or gases) with species()")
if(!find.IS) {
- # Get subcrt() output for all aqueous species and minerals 20210322
- # Add aqueous species here - the minerals are already present
+ # Get subcrt() output for all aqueous species and substrates 20210322
+ # Add aqueous species here - the substrates are already present
lapply(iaq, species, add = TRUE)
# Also add basis species for mosaic()!
if(is.mosaic) lapply(unlist(ddd$bases), species, add = TRUE)
@@ -49,15 +49,16 @@
sout <- suppressMessages(do.call(affinity, c(affargs, return.sout = TRUE)))
}
- # Make a list to store the calculated solubilities for each mineral
- slist <- list()
- # Loop over minerals
- for(i in seq_along(mineral$ispecies)) {
+ # Make lists to store the calculated solubilities (loga.balance) and speciation (loga.equil) for each substrate
+ balance_list <- list()
+ equil_list <- list()
+ # Loop over substrates
+ for(i in seq_along(substrates$ispecies)) {
# Print message
- message(paste("solubility: calculating for", mineral$name[i]))
- # Define basis species with the mineral first (so it will be dissolved)
- ispecies[1] <- mineral$ispecies[i]
- logact[1] <- mineral$logact[i]
+ message(paste("solubility: calculating for", substrates$name[i]))
+ # Define basis species with the substrate first (so it will be dissolved)
+ ispecies[1] <- substrates$ispecies[i]
+ logact[1] <- substrates$logact[i]
# Use numeric values first and put in buffer names second (needed for demo/gold.R)
loga.numeric <- suppressWarnings(as.numeric(logact))
basis(ispecies, loga.numeric)
@@ -70,33 +71,40 @@
} else {
if(is.mosaic) aout <- suppressMessages(mosaic(..., sout = sout)) else aout <- suppressMessages(affinity(..., sout = sout))
}
- # Calculate solubility of this mineral
+ # Calculate solubility of this substrate
scalc <- solubility_calc(aout, in.terms.of = in.terms.of, dissociate = dissociate, find.IS = find.IS)
- # Store the solubilities in the list
- slist[[i]] <- scalc$loga.balance
+ # Store the results
+ balance_list[[i]] <- scalc$loga.balance
+ equil_list[[i]] <- scalc$loga.equil
}
# Restore the original thermodynamic system settings
assign("thermo", thermo, CHNOSZ)
- if(length(mineral$ispecies) == 1) {
- # For one mineral, return the results of the solubility calculation
+ if(length(substrates$ispecies) == 1) {
+ # For one substrate, return the results of the solubility calculation
scalc
} else {
- # For multiple minerals, the overall solubility is the *minimum* among all the minerals
- smin <- do.call(pmin, slist)
- # Put this into the last-computed 'solubility' object
- scalc$loga.balance <- smin
- scalc$loga.equil <- slist
- scalc$species <- mineral
- # Change the function name stored in the object so diagram() plots loga.balance automatically
- scalc$fun <- "solubilities"
- # Return the object
- scalc
+ # Create two outputs (substrate and aqueous) using the last-computed 'solubility' object as a template
+ substrate <- aqueous <- scalc
+ # For multiple substrates, the overall solubility is the *minimum* among all the substrates
+ imin <- which.pmax(balance_list, maximum = FALSE)
+ substrate$loga.balance <- parallel_index(balance_list, imin)
+ substrate$loga.equil <- balance_list
+ substrate$species <- substrates
+ # Assign the function name so diagram() uses minimum solubility to identify the stable substrate
+ substrate$fun <- "solubilities"
+ # Now get the activities of aqueous species
+ aqueous$loga.balance <- parallel_index(balance_list, imin)
+ for(i in 1:length(aqueous$loga.equil)) {
+ this_equil_list <- sapply(equil_list, "[", i)
+ aqueous$loga.equil[[i]] <- parallel_index(this_equil_list, imin)
+ }
+ # Return the objects
+ return(list(substrate = substrate, aqueous = aqueous))
}
}
-
# The "nuts and bolts" of solubility calculations
# Moved from solubility() to solubility_calc() 20210318
solubility_calc <- function(aout, dissociate = NULL, find.IS = FALSE, in.terms.of = NULL) {
Modified: pkg/CHNOSZ/R/subcrt.R
===================================================================
--- pkg/CHNOSZ/R/subcrt.R 2026-06-02 12:33:23 UTC (rev 966)
+++ pkg/CHNOSZ/R/subcrt.R 2026-06-21 08:52:51 UTC (rev 967)
@@ -3,6 +3,7 @@
# 20060817 jmd
## If this file is interactively sourced, the following are also needed to provide unexported functions:
+#source("util.misc.R")
#source("util.args.R")
#source("util.character.R")
#source("info.R")
@@ -15,7 +16,7 @@
#source("cgl.R")
subcrt <- function(species, coeff = 1, state = NULL, property = c("logK", "G", "H", "S", "V", "Cp"),
- T = seq(273.15, 623.15, 25), P = "Psat", grid = NULL, convert = TRUE, exceed.Ttr = FALSE,
+ T = seq(273.15, 623.15, 25), P = "Psat", grid = NULL, convert = TRUE, warn.Ttr = TRUE,
exceed.rhomin = FALSE, logact = NULL, autobalance = TRUE, use.polymorphs = TRUE, IS = 0) {
# Revise the call if the states are the second argument
@@ -25,7 +26,7 @@
# This is missing coeff and T in order that missing values are correctly detected further below 20230621
newargs <- list(species = species, state = newstate,
property = property, P = P, grid = grid, convert = convert,
- exceed.Ttr = exceed.Ttr, exceed.rhomin = exceed.rhomin, logact = logact,
+ warn.Ttr = warn.Ttr, exceed.rhomin = exceed.rhomin, logact = logact,
autobalance = autobalance, use.polymorphs = use.polymorphs, IS = IS)
if(!missing(state)) {
if(is.numeric(state[1])) newcoeff <- state else stop("If they are both given, one of arguments 2 and 3 should be numeric reaction coefficients")
@@ -73,7 +74,7 @@
if(identical(P, "Psat") & any(T > 647.067)) {
nover <- sum(T > 647.067)
if(nover==1) vtext <- "value" else vtext <- "values"
- warnings <- c(warnings, paste0("P = 'Psat' undefined for T > Tcritical (", nover, " T ", vtext, ")"))
+ warnings <- c(warnings, warning("P = 'Psat' undefined for T > Tcritical (", nover, " T ", vtext, ")"))
}
# Are we gridding?
@@ -252,8 +253,8 @@
if(!is.null(logact)) {
ila <- match(colnames(bc), rownames(thermo$basis))
nla <- !(can.be.numeric(thermo$basis$logact[ila]))
- if(any(nla)) warning("subcrt: logact values of basis species",
- c2s(rownames(thermo$basis)[ila]), "are NA.")
+ if(any(nla)) warnings <- c(warnings, warning("subcrt: logact values of basis species",
+ c2s(rownames(thermo$basis)[ila]), "are NA."))
logact <- c(logact, thermo$basis$logact[ila])
}
# Warn user and do it!
@@ -267,9 +268,9 @@
newstate <- c(state, b.state)
return(subcrt(species = newspecies, coeff = newcoeff, state = newstate,
property = property, T = outvert(T, "K"), P = P, grid = grid, convert = convert, logact = logact,
- exceed.Ttr = exceed.Ttr, exceed.rhomin = exceed.rhomin, IS = IS))
- } else warnings <- c(warnings, paste("reaction among", paste(species, collapse = ","), "was unbalanced, missing", as.chemical.formula(miss)))
- } else warnings <- c(warnings, paste("reaction among", paste(species, collapse = ","), "was unbalanced, missing", as.chemical.formula(miss)))
+ warn.Ttr = warn.Ttr, exceed.rhomin = exceed.rhomin, IS = IS))
+ } else warnings <- c(warnings, warning("reaction among ", paste(species, collapse = ","), " was unbalanced, missing ", as.chemical.formula(miss)))
+ } else warnings <- c(warnings, warning("reaction among ", paste(species, collapse = ","), " was unbalanced, missing ", as.chemical.formula(miss)))
}
}
@@ -313,7 +314,7 @@
if(any(ilowrho)) {
for(i in 1:length(p.aq)) p.aq[[i]][ilowrho, ] <- NA
if(sum(ilowrho) == 1) ptext <- "pair" else ptext <- "pairs"
- warnings <- c(warnings, paste0("below minimum density for applicability of revised HKF equations (", sum(ilowrho), " T,P ", ptext, ")"))
+ warnings <- c(warnings, warning("below minimum density for applicability of revised HKF equations (", sum(ilowrho), " T,P ", ptext, ")"))
}
}
# Calculate properties using Akinfiev-Diamond model 20190219
@@ -343,7 +344,7 @@
# phases are beyond their temperature range
if("G" %in% eosprop) {
# 20080304 This code is weird and hard to read - needs a lot of cleanup!
- # 20120219 Cleaned up somewhat; using exceed.Ttr and NA instead of do.phases and 999999
+ # 20120219 Cleaned up somewhat; using warn.Ttr and NA instead of do.phases and 999999
# The numbers of the cgl species (becomes 0 for any that aren't cgl)
ncgl <- iscgl
ncgl[iscgl] <- 1:nrow(param)
@@ -361,13 +362,12 @@
Ttr <- Ttr(iphases[i-1], iphases[i], P = P, dPdT = dPdTtr(iphases[i-1], iphases[i]))
if(all(is.na(Ttr))) next
if(any(T <= Ttr)) {
- status.Ttr <- "(extrapolating G)"
- if(!exceed.Ttr) {
- # put NA into the value of G
- p.cgl[[ncgl[i]]]$G[T <= Ttr] <- NA
- status.Ttr <- "(using NA for G)"
- }
- #message(paste("subcrt: some points below transition temperature for", myname, mystate, status.Ttr))
+ # Put NA into the value of G
+ # NOTE: this is the key to getting the right polymorphs (mask the low-T ones at high T) 20260615
+ p.cgl[[ncgl[i]]]$G[T <= Ttr] <- NA
+ status.Ttr <- "(using NA for G)"
+ # This message comes up for any polymorphic transition - comment it but leave it for debugging 20260615
+ #message(paste("subcrt: some points below transition temperature for", myname, mystate))
}
}
}
@@ -387,23 +387,15 @@
}
if(!is.polymorphic.transition) {
- # Check if we're above the T limit for a Cp equation or a phase change (e.g. melting or vaporization)
+ # Check if we're above the transition temperature or a T limit for a Cp equation
if(all(is.na(Ttr))) next
if(all(Ttr == 0)) next
- is.phase_change <- any(thermo$OBIGT$model[iphases[i]] == "CGL_Ttr")
if(any(T > Ttr)) {
- if(is.phase_change) {
- if(exceed.Ttr) {
- message(paste0("subcrt: showing G for ", myname, "(", mystate, ") above its stability limit of ", Ttr, " K (use exceed.Ttr = FALSE to prevent this)"))
- } else {
- message(paste0("subcrt: setting G to NA for ", myname, "(", mystate, ") above its stability limit of ", Ttr, " K (use exceed.Ttr = TRUE to output G)"))
- p.cgl[[ncgl[i]]]$G[T > Ttr] <- NA
- }
- } else {
- if(! exceed.Ttr) {
- # Warn if we're above a Cp limit (but don't change the output)
- warning(paste0("above T limit of ", Ttr, " K for the Cp equation for ", myname, "(", mystate, ")"))
- }
+ # Message if we're above a T limit for Cp equation
+ message(paste0("subcrt: above T limit of ", Ttr, " K for the Cp equation for ", myname, "(", mystate, ")"))
+ if(warn.Ttr & any(thermo$OBIGT$model[iphases[i]] == "CGL_Ttr")) {
+ # Add a warning for transition temperature if warn.Ttr is TRUE
+ warnings <- c(warnings, warning("above transition temperature of ", Ttr, " K for ", myname, "(", mystate, ")"))
}
}
}
@@ -460,9 +452,9 @@
# Assemble the Gibbs energies for each species
for(j in 1:length(are.polymorphs)) {
G.this <- outprops[[are.polymorphs[j]]]$G
-# if(sum(is.na(G.this)) > 0 & exceed.Ttr) warning(paste("subcrt: NAs found for G of ",
+# if(sum(is.na(G.this)) > 0) warnings <- c(warnings, warning(paste("subcrt: NAs found for G of ",
# reaction$name[are.polymorphs[j]], " ", reaction$state[are.polymorphs[j]], " at T-P point(s) ",
-# c2s(which(is.na(G.this)), sep = " "), sep = ""), call. = FALSE)
+# c2s(which(is.na(G.this)), sep = " "), sep = ""), call. = FALSE))
if(j == 1) G <- as.data.frame(G.this)
else G <- cbind(G, as.data.frame(G.this))
}
@@ -477,8 +469,8 @@
#ps <- 1
# - above temperature limit for the highest-T phase (subcrt.Rd skarn example) --> use highest-T phase 20171110
ps <- ncol(G)
- if(exceed.Ttr) warning("subcrt: stable polymorph for ", reaction$name[are.polymorphs[ps]], " at T-P point ", j,
- " undetermined (using ", reaction$state[are.polymorphs[ps]], ")", call. = FALSE)
+ warnings <- c(warnings, warning("subcrt: stable polymorph for ", reaction$name[are.polymorphs[ps]], " at T-P point ", j,
+ " undetermined (using ", reaction$state[are.polymorphs[ps]], ")", call. = FALSE))
}
stable.polymorph <- c(stable.polymorph, ps)
out.new.entry[j, ] <- outprops[[ are.polymorphs[ps] ]][j, ]
@@ -618,9 +610,6 @@
}
}
# Add warnings to output 20180922
- if(length(warnings) > 0) {
- OUT <- c(OUT, list(warnings = warnings))
- for(warn in warnings) warning(warn)
- }
+ if(length(warnings) > 0) OUT <- c(OUT, list(warnings = warnings))
return(OUT)
}
Modified: pkg/CHNOSZ/R/thermo.R
===================================================================
--- pkg/CHNOSZ/R/thermo.R 2026-06-02 12:33:23 UTC (rev 966)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 967
More information about the CHNOSZ-commits
mailing list