[CHNOSZ-commits] r692 - in pkg/CHNOSZ: . R inst/tinytest man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Feb 3 12:57:50 CET 2022
Author: jedick
Date: 2022-02-03 12:57:50 +0100 (Thu, 03 Feb 2022)
New Revision: 692
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/berman.R
pkg/CHNOSZ/R/info.R
pkg/CHNOSZ/R/thermo.R
pkg/CHNOSZ/inst/tinytest/test-berman.R
pkg/CHNOSZ/man/thermo.Rd
Log:
Move read.csv for Berman data files to thermo()
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2022-02-03 09:57:37 UTC (rev 691)
+++ pkg/CHNOSZ/DESCRIPTION 2022-02-03 11:57:50 UTC (rev 692)
@@ -1,6 +1,6 @@
Date: 2022-02-03
Package: CHNOSZ
-Version: 1.4.1-17
+Version: 1.4.1-18
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/berman.R
===================================================================
--- pkg/CHNOSZ/R/berman.R 2022-02-03 09:57:37 UTC (rev 691)
+++ pkg/CHNOSZ/R/berman.R 2022-02-03 11:57:50 UTC (rev 692)
@@ -1,28 +1,22 @@
# CHNOSZ/berman.R 20170930
-# calculate thermodynamic properties of minerals using equations from:
+# Calculate thermodynamic properties of minerals using equations from:
# Berman, R. G. (1988) Internally-consistent thermodynamic data for minerals
# 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") {
- # reference temperature and pressure
+ # Reference temperature and pressure
Pr <- 1
Tr <- 298.15
- # get T and P to be same length
+ # Make T and P the same length
ncond <- max(length(T), length(P))
T <- rep(T, length.out=ncond)
P <- rep(P, length.out=ncond)
- # get thermodynamic parameters from data files
- path <- system.file("extdata/Berman/", package="CHNOSZ")
- files <- dir(path, "\\.csv$")
- # put files in reverse chronological order (youngest first)
- files <- rev(files[order(sapply(strsplit(files, "_"), "[", 2))])
- # read the parameters from each file
- dat <- list()
- for(i in 1:length(files)) dat[[i]] <- read.csv(file.path(path, files[i]), as.is = TRUE)
- # assemble the parameters in a single data frame
- dat <- do.call(rbind, dat)
- # is there a user-supplied data file?
+
+ # Get parameters in the Berman equations
+ # Start with thermodynamic parameters provided with CHNOSZ
+ dat <- thermo()$Berman
+ # Is there a user-supplied data file?
userfile <- get("thermo", CHNOSZ)$opt$Berman
userfileexists <- FALSE
if(!is.na(userfile)) {
@@ -34,33 +28,35 @@
} else stop("the file named in thermo()$opt$Berman (", userfile, ") does not exist")
}
}
- # remove duplicates (only the first, i.e. most recent entry is kept)
+ # Remove duplicates (only the first, i.e. most recent entry is kept)
dat <- dat[!duplicated(dat$name), ]
- # remove the multipliers on volume parameters
+ # Remove the multipliers on volume parameters
vcols <- 13:16 # columns with v1, v2, v3, v4
multexp <- c(5, 5, 5, 8)
dat[, vcols] <- t(t(dat[, vcols]) / 10^multexp)
- # if name is missing, return the entire data frame (used in test-berman.R)
- if(missing(name)) return(dat)
- # which row has data for this mineral?
- irow <- which(dat$name == name)
- if(length(irow)==0) {
- if(userfileexists) stop("Data for ", name, " not available. Please add it to ", userfile)
- if(!userfileexists) stop("Data for ", name, " not available. Please add it to your_data_file.csv and run thermo('opt$Berman' = 'path/to/your_data_file.csv')")
+ # If name is missing, return the entire data frame (used in test-berman.R)
+ if(missing(name)) return(dat) else {
+ # Which row has data for this mineral?
+ irow <- which(dat$name == name)
+ if(length(irow)==0) {
+ if(userfileexists) stop("Data for ", name, " not available. Please add it to ", userfile)
+ if(!userfileexists) stop("Data for ", name, " not available. Please add it to your_data_file.csv and run thermo('opt$Berman' = 'path/to/your_data_file.csv')")
+ }
+ dat <- dat[irow, ]
}
- # the function works fine with just the following assign() call,
+
+ # The function works fine with just the following assign() call,
# but an explicit dummy assignment here is used to avoid "Undefined global functions or variables" in R CMD check
GfPrTr <- HfPrTr <- SPrTr <- Tlambda <- Tmax <- Tmin <- Tref <- VPrTr <-
d0 <- d1 <- d2 <- d3 <- d4 <- Vad <- dTdP <- k0 <- k1 <- k2 <- k3 <-
k4 <- k5 <- k6 <- l1 <- l2 <- v1 <- v2 <- v3 <- v4 <- NA
- # assign values to the variables used below
- for(i in 1:ncol(dat)) assign(colnames(dat)[i], dat[irow, i])
- # get the entropy of the elements using the chemical formula in thermo()$OBIGT
- file <- system.file("extdata/OBIGT/Berman_cr.csv", package = "CHNOSZ")
- dat <- read.csv(file)
- formula <- dat$formula[match(name, dat$name)]
+ # Assign values to the variables used below
+ for(i in 1:ncol(dat)) assign(colnames(dat)[i], dat[, i])
+ # 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")
- # check that G in data file is the G of formation from the elements --> Benson-Helgeson convention (DG = DH - T*DS)
+ # 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)
Gdiff <- GfPrTr_calc - GfPrTr
@@ -70,14 +66,14 @@
# (the tabulated GfPrTr is unused below)
}
- ### thermodynamic properties ###
- # calculate Cp and V (Berman, 1988 Eqs. 4 and 5)
+ ### Thermodynamic properties ###
+ # Calculate Cp and V (Berman, 1988 Eqs. 4 and 5)
# k4, k5, k6 terms from winTWQ documentation (doi:10.4095/223425)
Cp <- k0 + k1 * T^-0.5 + k2 * T^-2 + k3 * T^-3 + k4 * T^-1 + k5 * T + k6 * T^2
P_Pr <- P - Pr
T_Tr <- T - Tr
V <- VPrTr * (1 + v1 * T_Tr + v2 * T_Tr^2 + v3 * P_Pr + v4 * P_Pr^2)
- ## calculate Ga (Ber88 Eq. 6) (superseded 20180328 as it does not include k4, k5, k6)
+ ## Calculate Ga (Ber88 Eq. 6) (superseded 20180328 as it does not include k4, k5, k6)
#Ga <- HfPrTr - T * SPrTr + k0 * ( (T - Tr) - T * (log(T) - log(Tr)) ) +
# 2 * k1 * ( (T^0.5 - Tr^0.5) + T*(T^-0.5 - Tr^-0.5) ) -
# k2 * ( (T^-1 - Tr^-1) - T / 2 * (T^-2 - Tr^-2) ) -
@@ -84,23 +80,23 @@
# k3 * ( (T^-2 - Tr^-2) / 2 - T / 3 * (T^-3 - Tr^-3) ) +
# VPrTr * ( (v3 / 2 - v4) * (P^2 - Pr^2) + v4 / 3 * (P^3 - Pr^3) +
# (1 - v3 + v4 + v1 * (T - Tr) + v2 * (T - Tr)^2) * (P - Pr) )
- # calculate Ha (symbolically integrated using sympy - expressions not simplified)
+ # Calculate Ha (symbolically integrated using sympy - expressions not simplified)
intCp <- T*k0 - Tr*k0 + k2/Tr - k2/T + k3/(2*Tr^2) - k3/(2*T^2) + 2.0*k1*T^0.5 - 2.0*k1*Tr^0.5 +
k4*log(T) - k4*log(Tr) + k5*T^2/2 - k5*Tr^2/2 - k6*Tr^3/3 + k6*T^3/3
intVminusTdVdT <- -VPrTr + P*(VPrTr + VPrTr*v4 - VPrTr*v3 - Tr*VPrTr*v1 + VPrTr*v2*Tr^2 - VPrTr*v2*T^2) +
P^2*(VPrTr*v3/2 - VPrTr*v4) + VPrTr*v3/2 - VPrTr*v4/3 + Tr*VPrTr*v1 + VPrTr*v2*T^2 - VPrTr*v2*Tr^2 + VPrTr*v4*P^3/3
Ha <- HfPrTr + intCp + intVminusTdVdT
- # calculate S (also symbolically integrated)
+ # Calculate S (also symbolically integrated)
intCpoverT <- k0*log(T) - k0*log(Tr) - k3/(3*T^3) + k3/(3*Tr^3) + k2/(2*Tr^2) - k2/(2*T^2) + 2.0*k1*Tr^-0.5 - 2.0*k1*T^-0.5 +
k4/Tr - k4/T + T*k5 - Tr*k5 + k6*T**2/2 - k6*Tr**2/2
intdVdT <- -VPrTr*(v1 + v2*(-2*Tr + 2*T)) + P*VPrTr*(v1 + v2*(-2*Tr + 2*T))
S <- SPrTr + intCpoverT - intdVdT
- # calculate Ga --> Berman-Brown convention (DG = DH - T*S, no S(element))
+ # Calculate Ga --> Berman-Brown convention (DG = DH - T*S, no S(element))
Ga <- Ha - T * S
- ### polymorphic transition properties ***
+ ### Polymorphic transition properties ***
if(!is.na(Tlambda) & !is.na(Tref) & any(T > Tref) & calc.transition) {
- # starting transition contributions are 0
+ # Starting transition contributions are 0
Cptr <- Htr <- Str <- Gtr <- numeric(ncond)
## Ber88 Eq. 8: Cp at 1 bar
#Cplambda_1bar <- T * (l1 + l2 * T)^2
@@ -109,23 +105,23 @@
# Eq. 8a: Cp at P
Td <- Tlambda - Tlambda_P
Tprime <- T + Td
- # with the condition that Tref < Tprime < Tlambda(1bar)
+ # With the condition that Tref < Tprime < Tlambda(1bar)
iTprime <- Tref < Tprime & Tprime < Tlambda
- # handle NA values (arising from NA in input P values e.g. Psat above Tcritical) 20180925
+ # Handle NA values (arising from NA in input P values e.g. Psat above Tcritical) 20180925
iTprime[is.na(iTprime)] <- FALSE
Tprime <- Tprime[iTprime]
Cptr[iTprime] <- Tprime * (l1 + l2 * Tprime)^2
- # we got Cp, now calculate the integrations for H and S
- # the lower integration limit is Tref
+ # We got Cp, now calculate the integrations for H and S
+ # The lower integration limit is Tref
iTtr <- T > Tref
Ttr <- T[iTtr]
Tlambda_P <- Tlambda_P[iTtr]
Td <- Td[iTtr]
- # handle NA values 20180925
+ # Handle NA values 20180925
Tlambda_P[is.na(Tlambda_P)] <- Inf
- # the upper integration limit is Tlambda_P
+ # The upper integration limit is Tlambda_P
Ttr[Ttr >= Tlambda_P] <- Tlambda_P[Ttr >= Tlambda_P]
- # derived variables
+ # Derived variables
tref <- Tref - Td
x1 <- l1^2 * Td + 2 * l1 * l2 * Td^2 + l2^2 * Td^3
x2 <- l1^2 + 4 * l1 * l2 * Td + 3 * l2^2 * Td^2
@@ -135,7 +131,7 @@
Htr[iTtr] <- x1 * (Ttr - tref) + x2 / 2 * (Ttr^2 - tref^2) + x3 / 3 * (Ttr^3 - tref^3) + x4 / 4 * (Ttr^4 - tref^4)
Str[iTtr] <- x1 * (log(Ttr) - log(tref)) + x2 * (Ttr - tref) + x3 / 2 * (Ttr^2 - tref^2) + x4 / 3 * (Ttr^3 - tref^3)
Gtr <- Htr - T * Str
- # apply the transition contributions
+ # Apply the transition contributions
Ga <- Ga + Gtr
Ha <- Ha + Htr
S <- S + Str
@@ -142,14 +138,14 @@
Cp <- Cp + Cptr
}
- ### disorder thermodynamic properties ###
+ ### Disorder thermodynamic properties ###
if(!is.na(Tmin) & !is.na(Tmax) & any(T > Tmin) & calc.disorder) {
- # starting disorder contributions are 0
+ # Starting disorder contributions are 0
Cpds <- Hds <- Sds <- Vds <- Gds <- numeric(ncond)
- # the lower integration limit is Tmin
+ # The lower integration limit is Tmin
iTds <- T > Tmin
Tds <- T[iTds]
- # the upper integration limit is Tmax
+ # The upper integration limit is Tmax
Tds[Tds > Tmax] <- Tmax
# Ber88 Eqs. 15, 16, 17
Cpds[iTds] <- d0 + d1*Tds^-0.5 + d2*Tds^-2 + d3*Tds + d4*Tds^2
@@ -164,18 +160,18 @@
if(Vad != 0) Vds <- Hds / Vad
# Berman puts the Vds term directly into Eq. 19 (commented below), but that necessarily makes Gds != Hds - T * Sds
#Gds <- Hds - T * Sds + Vds * (P - Pr)
- # instead, we include the Vds term with Hds
+ # Instead, we include the Vds term with Hds
Hds <- Hds + Vds * (P - Pr)
- # disordering properties above Tmax (Eq. 20)
+ # Disordering properties above Tmax (Eq. 20)
ihigh <- T > Tmax
- # again, Berman put the Sds term (for T > Tmax) into Eq. 20 for Gds (commented below), which would also make Gds != Hds - T * Sds
+ # Again, Berman put the Sds term (for T > Tmax) into Eq. 20 for Gds (commented below), which would also make Gds != Hds - T * Sds
#Gds[ihigh] <- Gds[ihigh] - (T[ihigh] - Tmax) * Sds[ihigh]
- # instead, we add the Sds[ihigh] term to Hds
+ # Instead, we add the Sds[ihigh] term to Hds
Hds[ihigh] <- Hds[ihigh] - (T[ihigh] - Tmax) * Sds[ihigh]
- # by writing Gds = Hds - T * Sds, the above two changes w.r.t. Berman's
+ # By writing Gds = Hds - T * Sds, the above two changes w.r.t. Berman's
# equations affect the computed values only for Hds, not Gds
Gds <- Hds - T * Sds
- # apply the disorder contributions
+ # Apply the disorder contributions
Ga <- Ga + Gds
Ha <- Ha + Hds
S <- S + Sds
@@ -183,17 +179,17 @@
Cp <- Cp + Cpds
}
- ### (for testing) use G = H - TS to check that integrals for H and S are written correctly
+ ### (for testing) Use G = H - TS to check that integrals for H and S are written correctly
Ga_fromHminusTS <- Ha - T * S
if(!isTRUE(all.equal(Ga_fromHminusTS, Ga))) stop(paste0(name, ": incorrect integrals detected using DG = DH - T*S"))
- ### thermodynamic and unit conventions used in SUPCRT ###
- # use entropy of the elements in calculation of G --> Benson-Helgeson convention (DG = DH - T*DS)
+ ### Thermodynamic and unit conventions used in SUPCRT ###
+ # Use entropy of the elements in calculation of G --> Benson-Helgeson convention (DG = DH - T*DS)
Gf <- Ga + Tr * SPrTr_elements
- # the output will just have "G" and "H"
+ # The output will just have "G" and "H"
G <- Gf
H <- Ha
- # convert J to cal
+ # Convert J to cal
if(grepl("cal", units)) {
G <- convert(Gf, "cal")
H <- convert(Ha, "cal")
@@ -200,7 +196,7 @@
S <- convert(S, "cal")
Cp <- convert(Cp, "cal")
}
- # convert J/bar to cm^3/mol
+ # 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)
Modified: pkg/CHNOSZ/R/info.R
===================================================================
--- pkg/CHNOSZ/R/info.R 2022-02-03 09:57:37 UTC (rev 691)
+++ pkg/CHNOSZ/R/info.R 2022-02-03 11:57:50 UTC (rev 692)
@@ -196,9 +196,10 @@
this <- OBIGT2eos(this, this$state)
if(all(is.na(this[, 9:21])) & this$name != "water") {
- # Get G, H, S, and V from datasets using Berman equations 20220203
- properties <- subset(berman(), name == this$name)
- this[, c("G", "H", "S", "V")] <- properties[, c("GfPrTr", "HfPrTr", "SPrTr", "VPrTr")] * c(1, 1, 1, 10)
+ # Get G, H, S, and V for minerals with Berman parameters 20220203
+ bermandat <- berman()
+ bermandat <- bermandat[bermandat$name == this$name, ]
+ this[, c("G", "H", "S", "V")] <- bermandat[, c("GfPrTr", "HfPrTr", "SPrTr", "VPrTr")] * c(1, 1, 1, 10)
}
# Identify any missing GHS values
Modified: pkg/CHNOSZ/R/thermo.R
===================================================================
--- pkg/CHNOSZ/R/thermo.R 2022-02-03 09:57:37 UTC (rev 691)
+++ pkg/CHNOSZ/R/thermo.R 2022-02-03 11:57:50 UTC (rev 692)
@@ -14,6 +14,7 @@
element = read.csv(file.path(thermodir, "element.csv"), as.is=1:3),
OBIGT = NULL,
refs = NULL,
+ Berman = NULL,
buffer = read.csv(file.path(thermodir, "buffer.csv"), as.is=1:3),
protein = read.csv(file.path(thermodir, "protein.csv"), as.is=1:4),
groups = read.csv(file.path(thermodir, "groups.csv"), row.names=1, check.names=FALSE),
@@ -22,10 +23,23 @@
species = NULL,
opar = NULL
)
+
# store stoich as matrix (with non-unique row names), not data frame
formula <- thermo$stoich[, 1]
thermo$stoich <- as.matrix(thermo$stoich[, 2:ncol(thermo$stoich)])
rownames(thermo$stoich) <- formula
+
+ # Get parameters in Berman equations from data files 20220203
+ path <- system.file("extdata/Berman/", package = "CHNOSZ")
+ files <- dir(path, "\\.csv$")
+ # Put files in reverse chronological order (youngest first)
+ files <- rev(files[order(sapply(strsplit(files, "_"), "[", 2))])
+ # Read the parameters from each file
+ Berman <- list()
+ for(i in 1:length(files)) Berman[[i]] <- read.csv(file.path(path, files[i]), as.is = TRUE)
+ # Assemble the parameters in a single data frame
+ thermo$Berman <- do.call(rbind, Berman)
+
# give a summary of what we are doing
if(!"thermo" %in% ls(CHNOSZ)) message("reset: creating \"thermo\" object")
else message("reset: resetting \"thermo\" object")
Modified: pkg/CHNOSZ/inst/tinytest/test-berman.R
===================================================================
--- pkg/CHNOSZ/inst/tinytest/test-berman.R 2022-02-03 09:57:37 UTC (rev 691)
+++ pkg/CHNOSZ/inst/tinytest/test-berman.R 2022-02-03 11:57:50 UTC (rev 692)
@@ -59,11 +59,6 @@
# 20191116
expect_false(any(is.na(subcrt("K-feldspar", P = 1, T = seq(273.15, 303.15, 5), convert = FALSE)$out[[1]]$G)), info = info)
-
-
-# The next set of tests are long ... only run them "at home" 20220129
-if(!at_home()) exit_file("Skipping long test")
-
# Get parameters for all available minerals
dat <- berman()
mineral <- unique(dat$name)
Modified: pkg/CHNOSZ/man/thermo.Rd
===================================================================
--- pkg/CHNOSZ/man/thermo.Rd 2022-02-03 09:57:37 UTC (rev 691)
+++ pkg/CHNOSZ/man/thermo.Rd 2022-02-03 11:57:50 UTC (rev 692)
@@ -238,6 +238,8 @@
\code{...} \tab numeric \tab Stoichiometry, one column for each element present in any species
}
+ \item \code{thermo()$Berman}
+ A data frame with thermodynamic parameters for minerals in the Berman equations, assembled from files in \samp{extdata/Berman} and used in \code{\link{berman}}.
} % end of itemize with long descriptions
More information about the CHNOSZ-commits
mailing list