[CHNOSZ-commits] r693 - in pkg/CHNOSZ: . R demo inst inst/extdata/Berman/testing inst/tinytest man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Feb 3 13:13:19 CET 2022
Author: jedick
Date: 2022-02-03 13:13:18 +0100 (Thu, 03 Feb 2022)
New Revision: 693
Added:
pkg/CHNOSZ/R/Berman.R
pkg/CHNOSZ/inst/extdata/Berman/testing/BA96_Berman.csv
pkg/CHNOSZ/inst/tinytest/test-Berman.R
pkg/CHNOSZ/man/Berman.Rd
Removed:
pkg/CHNOSZ/R/berman.R
pkg/CHNOSZ/inst/extdata/Berman/testing/BA96_berman.csv
pkg/CHNOSZ/inst/tinytest/test-berman.R
pkg/CHNOSZ/man/berman.Rd
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/NAMESPACE
pkg/CHNOSZ/R/cgl.R
pkg/CHNOSZ/R/examples.R
pkg/CHNOSZ/R/info.R
pkg/CHNOSZ/demo/lambda.R
pkg/CHNOSZ/inst/NEWS.Rd
pkg/CHNOSZ/man/CHNOSZ-package.Rd
pkg/CHNOSZ/man/diagram.Rd
pkg/CHNOSZ/man/extdata.Rd
pkg/CHNOSZ/man/thermo.Rd
pkg/CHNOSZ/vignettes/OBIGT.Rmd
Log:
Change "berman" to "Berman"
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2022-02-03 11:57:50 UTC (rev 692)
+++ pkg/CHNOSZ/DESCRIPTION 2022-02-03 12:13:18 UTC (rev 693)
@@ -1,6 +1,6 @@
Date: 2022-02-03
Package: CHNOSZ
-Version: 1.4.1-18
+Version: 1.4.1-19
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 2022-02-03 11:57:50 UTC (rev 692)
+++ pkg/CHNOSZ/NAMESPACE 2022-02-03 12:13:18 UTC (rev 693)
@@ -50,7 +50,7 @@
"nonideal", "uniprot.aa",
# added 20170301 or later
"GHS_Tr", "calculateDensity", "calculateGibbsOfWater",
- "calculateEpsilon", "calculateQ", "water.DEW", "berman",
+ "calculateEpsilon", "calculateQ", "water.DEW", "Berman",
"bgamma",
# added 20171121 or later
"dumpdata", "thermo.axis", "solubility", "NaCl",
Copied: pkg/CHNOSZ/R/Berman.R (from rev 692, pkg/CHNOSZ/R/berman.R)
===================================================================
--- pkg/CHNOSZ/R/Berman.R (rev 0)
+++ pkg/CHNOSZ/R/Berman.R 2022-02-03 12:13:18 UTC (rev 693)
@@ -0,0 +1,203 @@
+# CHNOSZ/Berman.R 20170930
+# 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
+ Pr <- 1
+ Tr <- 298.15
+ # 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 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)) {
+ if(userfile!="") {
+ if(file.exists(userfile)) {
+ userfileexists <- TRUE
+ BDat_user <- read.csv(userfile, as.is=TRUE)
+ dat <- rbind(BDat_user, dat)
+ } 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)
+ dat <- dat[!duplicated(dat$name), ]
+ # 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) 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,
+ # 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[, 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)
+ if(check.G) {
+ GfPrTr_calc <- HfPrTr - Tr * (SPrTr - SPrTr_elements)
+ Gdiff <- GfPrTr_calc - GfPrTr
+ #if(is.na(GfPrTr)) warning(paste0(name, ": GfPrTr(table) is NA"), call.=FALSE)
+ if(!is.na(GfPrTr)) if(abs(Gdiff) >= 1000) warning(paste0(name, ": GfPrTr(calc) - GfPrTr(table) is too big! == ",
+ round(GfPrTr_calc - GfPrTr), " J/mol"), call.=FALSE)
+ # (the tabulated GfPrTr is unused below)
+ }
+
+ ### 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)
+ #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) ) -
+ # 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)
+ 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)
+ 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))
+ Ga <- Ha - T * S
+
+ ### Polymorphic transition properties ***
+ if(!is.na(Tlambda) & !is.na(Tref) & any(T > Tref) & calc.transition) {
+ # 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
+ # Eq. 9: Tlambda at P
+ Tlambda_P <- Tlambda + dTdP * (P - 1)
+ # Eq. 8a: Cp at P
+ Td <- Tlambda - Tlambda_P
+ Tprime <- T + Td
+ # 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
+ 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
+ iTtr <- T > Tref
+ Ttr <- T[iTtr]
+ Tlambda_P <- Tlambda_P[iTtr]
+ Td <- Td[iTtr]
+ # Handle NA values 20180925
+ Tlambda_P[is.na(Tlambda_P)] <- Inf
+ # The upper integration limit is Tlambda_P
+ Ttr[Ttr >= Tlambda_P] <- Tlambda_P[Ttr >= Tlambda_P]
+ # 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
+ x3 <- 2 * l1 * l2 + 3 * l2^2 * Td
+ x4 <- l2 ^ 2
+ # Eqs. 10, 11, 12
+ 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
+ Ga <- Ga + Gtr
+ Ha <- Ha + Htr
+ S <- S + Str
+ Cp <- Cp + Cptr
+ }
+
+ ### Disorder thermodynamic properties ###
+ if(!is.na(Tmin) & !is.na(Tmax) & any(T > Tmin) & calc.disorder) {
+ # Starting disorder contributions are 0
+ Cpds <- Hds <- Sds <- Vds <- Gds <- numeric(ncond)
+ # The lower integration limit is Tmin
+ iTds <- T > Tmin
+ Tds <- T[iTds]
+ # 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
+ Hds[iTds] <- d0*(Tds - Tmin) + d1*(Tds^0.5 - Tmin^0.5)/0.5 +
+ d2*(Tds^-1 - Tmin^-1)/-1 + d3*(Tds^2 - Tmin^2)/2 + d4*(Tds^3 - Tmin^3)/3
+ Sds[iTds] <- d0*(log(Tds) - log(Tmin)) + d1*(Tds^-0.5 - Tmin^-0.5)/-0.5 +
+ d2*(Tds^-2 - Tmin^-2)/-2 + d3*(Tds - Tmin) + d4*(Tds^2 - Tmin^2)/2
+ # "d5 is a constant computed in such as way as to scale the disordring enthalpy to the volume of disorder" (Berman, 1988)
+ # 20180331: however, having a "d5" that isn't a coefficient in the same equation as d0, d1, d2, d3, d4 is confusing notation
+ # therefore, CHNOSZ now uses "Vad" for this variable, following the notation in the Theriak-Domino manual
+ # Eq. 18; we can't do this if Vad == 0 (dolomite and gehlenite)
+ 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
+ Hds <- Hds + Vds * (P - Pr)
+ # 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
+ #Gds[ihigh] <- Gds[ihigh] - (T[ihigh] - Tmax) * Sds[ihigh]
+ # 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
+ # equations affect the computed values only for Hds, not Gds
+ Gds <- Hds - T * Sds
+ # Apply the disorder contributions
+ Ga <- Ga + Gds
+ Ha <- Ha + Hds
+ S <- S + Sds
+ V <- V + Vds
+ Cp <- Cp + Cpds
+ }
+
+ ### (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)
+ Gf <- Ga + Tr * SPrTr_elements
+ # 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)
+}
Deleted: pkg/CHNOSZ/R/berman.R
===================================================================
--- pkg/CHNOSZ/R/berman.R 2022-02-03 11:57:50 UTC (rev 692)
+++ pkg/CHNOSZ/R/berman.R 2022-02-03 12:13:18 UTC (rev 693)
@@ -1,203 +0,0 @@
-# CHNOSZ/berman.R 20170930
-# 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
- Pr <- 1
- Tr <- 298.15
- # 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 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)) {
- if(userfile!="") {
- if(file.exists(userfile)) {
- userfileexists <- TRUE
- BDat_user <- read.csv(userfile, as.is=TRUE)
- dat <- rbind(BDat_user, dat)
- } 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)
- dat <- dat[!duplicated(dat$name), ]
- # 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) 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,
- # 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[, 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)
- if(check.G) {
- GfPrTr_calc <- HfPrTr - Tr * (SPrTr - SPrTr_elements)
- Gdiff <- GfPrTr_calc - GfPrTr
- #if(is.na(GfPrTr)) warning(paste0(name, ": GfPrTr(table) is NA"), call.=FALSE)
- if(!is.na(GfPrTr)) if(abs(Gdiff) >= 1000) warning(paste0(name, ": GfPrTr(calc) - GfPrTr(table) is too big! == ",
- round(GfPrTr_calc - GfPrTr), " J/mol"), call.=FALSE)
- # (the tabulated GfPrTr is unused below)
- }
-
- ### 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)
- #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) ) -
- # 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)
- 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)
- 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))
- Ga <- Ha - T * S
-
- ### Polymorphic transition properties ***
- if(!is.na(Tlambda) & !is.na(Tref) & any(T > Tref) & calc.transition) {
- # 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
- # Eq. 9: Tlambda at P
- Tlambda_P <- Tlambda + dTdP * (P - 1)
- # Eq. 8a: Cp at P
- Td <- Tlambda - Tlambda_P
- Tprime <- T + Td
- # 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
- 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
- iTtr <- T > Tref
- Ttr <- T[iTtr]
- Tlambda_P <- Tlambda_P[iTtr]
- Td <- Td[iTtr]
- # Handle NA values 20180925
- Tlambda_P[is.na(Tlambda_P)] <- Inf
- # The upper integration limit is Tlambda_P
- Ttr[Ttr >= Tlambda_P] <- Tlambda_P[Ttr >= Tlambda_P]
- # 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
- x3 <- 2 * l1 * l2 + 3 * l2^2 * Td
- x4 <- l2 ^ 2
- # Eqs. 10, 11, 12
- 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
- Ga <- Ga + Gtr
- Ha <- Ha + Htr
- S <- S + Str
- Cp <- Cp + Cptr
- }
-
- ### Disorder thermodynamic properties ###
- if(!is.na(Tmin) & !is.na(Tmax) & any(T > Tmin) & calc.disorder) {
- # Starting disorder contributions are 0
- Cpds <- Hds <- Sds <- Vds <- Gds <- numeric(ncond)
- # The lower integration limit is Tmin
- iTds <- T > Tmin
- Tds <- T[iTds]
- # 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
- Hds[iTds] <- d0*(Tds - Tmin) + d1*(Tds^0.5 - Tmin^0.5)/0.5 +
- d2*(Tds^-1 - Tmin^-1)/-1 + d3*(Tds^2 - Tmin^2)/2 + d4*(Tds^3 - Tmin^3)/3
- Sds[iTds] <- d0*(log(Tds) - log(Tmin)) + d1*(Tds^-0.5 - Tmin^-0.5)/-0.5 +
- d2*(Tds^-2 - Tmin^-2)/-2 + d3*(Tds - Tmin) + d4*(Tds^2 - Tmin^2)/2
- # "d5 is a constant computed in such as way as to scale the disordring enthalpy to the volume of disorder" (Berman, 1988)
- # 20180331: however, having a "d5" that isn't a coefficient in the same equation as d0, d1, d2, d3, d4 is confusing notation
- # therefore, CHNOSZ now uses "Vad" for this variable, following the notation in the Theriak-Domino manual
- # Eq. 18; we can't do this if Vad == 0 (dolomite and gehlenite)
- 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
- Hds <- Hds + Vds * (P - Pr)
- # 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
- #Gds[ihigh] <- Gds[ihigh] - (T[ihigh] - Tmax) * Sds[ihigh]
- # 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
- # equations affect the computed values only for Hds, not Gds
- Gds <- Hds - T * Sds
- # Apply the disorder contributions
- Ga <- Ga + Gds
- Ha <- Ha + Hds
- S <- S + Sds
- V <- V + Vds
- Cp <- Cp + Cpds
- }
-
- ### (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)
- Gf <- Ga + Tr * SPrTr_elements
- # 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)
-}
Modified: pkg/CHNOSZ/R/cgl.R
===================================================================
--- pkg/CHNOSZ/R/cgl.R 2022-02-03 11:57:50 UTC (rev 692)
+++ pkg/CHNOSZ/R/cgl.R 2022-02-03 12:13:18 UTC (rev 693)
@@ -16,7 +16,7 @@
PAR <- parameters[k, ]
if(all(is.na(PAR[9:21]))) {
# use Berman equations (parameters not in thermo()$OBIGT)
- properties <- berman(PAR$name, T = T, P = P)
+ properties <- Berman(PAR$name, T = T, P = P)
iprop <- match(property, colnames(properties))
values <- properties[, iprop, drop=FALSE]
} else {
Modified: pkg/CHNOSZ/R/examples.R
===================================================================
--- pkg/CHNOSZ/R/examples.R 2022-02-03 11:57:50 UTC (rev 692)
+++ pkg/CHNOSZ/R/examples.R 2022-02-03 12:13:18 UTC (rev 693)
@@ -9,7 +9,7 @@
"util.array", "util.blast", "util.data", "util.expression", "util.legend", "util.plot",
"util.fasta", "util.formula", "util.misc", "util.seq", "util.units",
"util.water", "taxonomy", "info", "retrieve", "add.OBIGT", "protein.info",
- "hkf", "water", "IAPWS95", "subcrt", "berman",
+ "hkf", "water", "IAPWS95", "subcrt", "Berman",
"makeup", "basis", "swap.basis", "species", "affinity",
"solubility", "equilibrate",
"diagram", "mosaic", "mix",
Modified: pkg/CHNOSZ/R/info.R
===================================================================
--- pkg/CHNOSZ/R/info.R 2022-02-03 11:57:50 UTC (rev 692)
+++ pkg/CHNOSZ/R/info.R 2022-02-03 12:13:18 UTC (rev 693)
@@ -197,9 +197,9 @@
if(all(is.na(this[, 9:21])) & this$name != "water") {
# 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)
+ 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/demo/lambda.R
===================================================================
--- pkg/CHNOSZ/demo/lambda.R 2022-02-03 11:57:50 UTC (rev 692)
+++ pkg/CHNOSZ/demo/lambda.R 2022-02-03 12:13:18 UTC (rev 693)
@@ -21,8 +21,8 @@
Tlab <- axis.label("T")
# calculate properties at 1 kbar with and without transition
-Qz_1bar <- berman("quartz", T=T, units="J")
-Qz_1bar_notrans <- berman("quartz", T=T, calc.transition=FALSE, units="J")
+Qz_1bar <- Berman("quartz", T=T, units="J")
+Qz_1bar_notrans <- Berman("quartz", T=T, calc.transition=FALSE, units="J")
## Fig. 1a: volume
plot(TC, Qz_1bar$V, type="l", xlab=Tlab, ylab=Vlab, ylim=c(22.5, 24))
@@ -40,8 +40,8 @@
## calculate finite difference derivative of dGlambda/dP = Vlambda between 1 and 1.001 bar
Glambda_1bar <- Qz_1bar_notrans$G - Qz_1bar$G
-Qz_1.001bar <- berman("quartz", T=T, P=1.001, units="J")
-Qz_1.001bar_notrans <- berman("quartz", T=T, P=1.001, calc.transition=FALSE, units="J")
+Qz_1.001bar <- Berman("quartz", T=T, P=1.001, units="J")
+Qz_1.001bar_notrans <- Berman("quartz", T=T, P=1.001, calc.transition=FALSE, units="J")
Glambda_1.001bar <- Qz_1.001bar_notrans$G - Qz_1.001bar$G
dGlambdadP <- (Glambda_1.001bar - Glambda_1bar) / 0.001
# we're using Joules, so multiply by ten to get cm^3
@@ -48,7 +48,7 @@
Vlambda <- -10 * dGlambdadP
VQz <- Qz_1bar$V + Vlambda
# above 848K, we have beta quartz
-Qz_beta <- berman("quartz,beta", T=T, P=1, units="J")
+Qz_beta <- Berman("quartz,beta", T=T, P=1, units="J")
VQz[T >= 848.15] <- Qz_beta$V[T >= 848.14]
lines(TC, VQz, lty=2)
legend("topleft", legend="1 bar", bty="n")
@@ -61,8 +61,8 @@
labplot("b")
## calculate properties at 10 kbar with and without transition
-Qz_10bar <- berman("quartz", T=T, P=10000, units="J")
-Qz_10bar_notrans <- berman("quartz", T=T, P=10000, calc.transition=FALSE, units="J")
+Qz_10bar <- Berman("quartz", T=T, P=10000, units="J")
+Qz_10bar_notrans <- Berman("quartz", T=T, P=10000, calc.transition=FALSE, units="J")
## Fig. 1c: heat capacity
plot(TC, Qz_10bar$Cp, type="l", xlab=Tlab, ylab=Cplab)
lines(TC, Qz_10bar_notrans$Cp, lty=3)
@@ -75,8 +75,8 @@
Pkb <- seq(1, 50, 1)
P <- 1000 * Pkb
T <- Tlambda + dTdP * (P - 1)
-Qz_withtrans <- berman("quartz", T=T, P=P, units="J")
-Qz_notrans <- berman("quartz", T=T, P=P, calc.transition=FALSE, units="J")
+Qz_withtrans <- Berman("quartz", T=T, P=P, units="J")
+Qz_notrans <- Berman("quartz", T=T, P=P, calc.transition=FALSE, units="J")
Qz_lambda <- Qz_withtrans - Qz_notrans
Plab <- expression(list(italic(P), "kb"))
plot(Pkb, Qz_lambda$G, type="l", ylim=c(-300, -50), ylab=axis.label("DlG"), xlab=Plab)
Modified: pkg/CHNOSZ/inst/NEWS.Rd
===================================================================
--- pkg/CHNOSZ/inst/NEWS.Rd 2022-02-03 11:57:50 UTC (rev 692)
+++ pkg/CHNOSZ/inst/NEWS.Rd 2022-02-03 12:13:18 UTC (rev 693)
@@ -52,6 +52,9 @@
\item \code{info()} with a numeric argument now includes values of
\code{G}, \code{H}, \code{S}, and \code{V} at 25 \degC and 1 bar for
minerals with thermodynamic parameters in the Berman equations.
+
+ \item All function, variables, and data file names now use capitalized
+ \samp{Berman}, not \samp{berman}.
}
}
Copied: pkg/CHNOSZ/inst/extdata/Berman/testing/BA96_Berman.csv (from rev 692, pkg/CHNOSZ/inst/extdata/Berman/testing/BA96_berman.csv)
===================================================================
--- pkg/CHNOSZ/inst/extdata/Berman/testing/BA96_Berman.csv (rev 0)
+++ pkg/CHNOSZ/inst/extdata/Berman/testing/BA96_Berman.csv 2022-02-03 12:13:18 UTC (rev 693)
@@ -0,0 +1,21 @@
+name,GfPrTr,HfPrTr,SPrTr,VPrTr,k0,k1,k2,k3,k4,k5,k6,v1,v2,v3,v4,Tlambda,Tref,dTdP,l1,l2,DtH,Tmax,Tmin,d0,d1,d2,d3,d4,Vad
+almandine,,-5265440,341.51,11.529,621.43,-3287.9,-15081000,2211870000,0,0,0,1.8599,0.00074711,-0.057,4.34E-05,,,,,,,,,,,,,,
+anthophyllite,,-12074470,537.4,26.33,1233.79,-7134,-22163800,2333940000,0,0,0,2.8105,0.00062894,-0.1139,0,,,,,,,,,,,,,,
+cordierite,,-9161480,416.23,23.311,954.39,-7962.3,-2317300,-370210000,0,0,0,0.3003,0.00018017,-0.1158,0,,,,,,,,,,,,,,
+Fe-cordierite,,-8430550,482.18,23.706,983.48,-8403.7,-1870300,-85680000,0,0,0,0.4265,0,-0.17,0,,,,,,,,,,,,,,
+fayalite,,-1477170,151.73,4.639,252,-2013.7,0,-62190000,0,0,0,2.621,0.00084233,-0.0822,0.0001944,,,,,,,,,,,,,,
+ferrosilite,,-1192910,96.47,3.295,174.2,-1393,-454400,-37710000,0,0,0,3.1406,0.000804,-0.0911,3.03E-05,,,,,,,,,,,,,,
+forsterite,,-2174420,94.18,4.36,233.18,-1801.6,0,-267940000,0,0,0,2.9464,0.00088633,-0.0791,0.0001351,,,,,,,,,,,,,,
+geikelite,,-1570520,74.41,3.077,146.2,-416,-3999800,402330000,0,0,0,2.3314,0.00088328,-0.0529,0,,,,,,,,,,,,,,
+hematite,,-826740,87.36,3.027,146.86,0,-5576800,525630000,0,0,0,3.831,1.65E-05,-0.0479,3.04E-05,955,298,0,-0.07403,0.00027921,1287,,,,,,,,
+hercynite,,-1945360,123.13,4.08,251.77,-2044.4,-1348300,131500000,0,0,0,1.5819,0.00096276,-0.051,0,,,,,,,,,,,,,,
+ilmenite,,-1233320,108.5,3.17,150,-441.6,-3323700,348150000,0,0,0,2.3314,0.00088328,-0.0529,0,,,,,,,,,,,,,,
+iron-a,,0,27.45,0.709,51.87,-379.4,-2543000,506800000,0,0,0,4.5071,0.00014104,-0.0602,0,1042,298,-0.00057,-0.22208,0.00034823,1000,,,,,,,,
+iron-g,,7720,35.64,0.691,66.24,-1237.1,6373300,-1060400000,0,0,0,4.7153,0.00014451,-0.0467,-0.0001445,,,,,,,,,,,,,,
+magnetite,,-1116960,146.04,4.452,207.93,0,-7243300,664360000,0,0,0,3.0291,0.0013847,-0.0582,0.0001751,848,298,0,-0.19502,0.00061037,1565,,,,,,,,
+enstatite,,-1546040,66.18,3.133,166.58,-1200.6,-2270600,279150000,0,0,0,2.4656,0.0007467,-0.0749,4.47E-05,,,,,,,,,,,,,,
+orthocorundum,,-1634950,33.93,3.123,119.38,774.8,-6509100,422880000,0,0,0,2.4656,0.0007467,-0.0749,4.47E-05,,,,,,,,,,,,,,
+pyrope,,-6284740,268.8,11.311,590.9,-2827,-13320800,1260330000,0,0,0,2.2519,0.00037044,-0.0576,4.42E-05,,,,,,,,,,,,,,
+rutile,,-944750,50.88,1.883,77.84,0,-3367800,402940000,0,0,0,2.5716,0.00015409,-0.0454,5.85E-05,,,,,,,,,,,,,,
+spinel,,-2302160,80.37,3.978,244.67,-2004,0,0,0,0,0,2.1691,0.00050528,-0.0489,0,,,,,,,,,,,,,,
+talc,,-5898960,261.21,13.61,664.11,-5187.2,-2147200,-327370000,0,0,0,2.5616,0,-0.1847,0.0005878,,,,,,,,,,,,,,
Deleted: pkg/CHNOSZ/inst/extdata/Berman/testing/BA96_berman.csv
===================================================================
--- pkg/CHNOSZ/inst/extdata/Berman/testing/BA96_berman.csv 2022-02-03 11:57:50 UTC (rev 692)
+++ pkg/CHNOSZ/inst/extdata/Berman/testing/BA96_berman.csv 2022-02-03 12:13:18 UTC (rev 693)
@@ -1,21 +0,0 @@
-name,GfPrTr,HfPrTr,SPrTr,VPrTr,k0,k1,k2,k3,k4,k5,k6,v1,v2,v3,v4,Tlambda,Tref,dTdP,l1,l2,DtH,Tmax,Tmin,d0,d1,d2,d3,d4,Vad
-almandine,,-5265440,341.51,11.529,621.43,-3287.9,-15081000,2211870000,0,0,0,1.8599,0.00074711,-0.057,4.34E-05,,,,,,,,,,,,,,
-anthophyllite,,-12074470,537.4,26.33,1233.79,-7134,-22163800,2333940000,0,0,0,2.8105,0.00062894,-0.1139,0,,,,,,,,,,,,,,
-cordierite,,-9161480,416.23,23.311,954.39,-7962.3,-2317300,-370210000,0,0,0,0.3003,0.00018017,-0.1158,0,,,,,,,,,,,,,,
-Fe-cordierite,,-8430550,482.18,23.706,983.48,-8403.7,-1870300,-85680000,0,0,0,0.4265,0,-0.17,0,,,,,,,,,,,,,,
-fayalite,,-1477170,151.73,4.639,252,-2013.7,0,-62190000,0,0,0,2.621,0.00084233,-0.0822,0.0001944,,,,,,,,,,,,,,
-ferrosilite,,-1192910,96.47,3.295,174.2,-1393,-454400,-37710000,0,0,0,3.1406,0.000804,-0.0911,3.03E-05,,,,,,,,,,,,,,
-forsterite,,-2174420,94.18,4.36,233.18,-1801.6,0,-267940000,0,0,0,2.9464,0.00088633,-0.0791,0.0001351,,,,,,,,,,,,,,
-geikelite,,-1570520,74.41,3.077,146.2,-416,-3999800,402330000,0,0,0,2.3314,0.00088328,-0.0529,0,,,,,,,,,,,,,,
-hematite,,-826740,87.36,3.027,146.86,0,-5576800,525630000,0,0,0,3.831,1.65E-05,-0.0479,3.04E-05,955,298,0,-0.07403,0.00027921,1287,,,,,,,,
-hercynite,,-1945360,123.13,4.08,251.77,-2044.4,-1348300,131500000,0,0,0,1.5819,0.00096276,-0.051,0,,,,,,,,,,,,,,
-ilmenite,,-1233320,108.5,3.17,150,-441.6,-3323700,348150000,0,0,0,2.3314,0.00088328,-0.0529,0,,,,,,,,,,,,,,
-iron-a,,0,27.45,0.709,51.87,-379.4,-2543000,506800000,0,0,0,4.5071,0.00014104,-0.0602,0,1042,298,-0.00057,-0.22208,0.00034823,1000,,,,,,,,
-iron-g,,7720,35.64,0.691,66.24,-1237.1,6373300,-1060400000,0,0,0,4.7153,0.00014451,-0.0467,-0.0001445,,,,,,,,,,,,,,
-magnetite,,-1116960,146.04,4.452,207.93,0,-7243300,664360000,0,0,0,3.0291,0.0013847,-0.0582,0.0001751,848,298,0,-0.19502,0.00061037,1565,,,,,,,,
-enstatite,,-1546040,66.18,3.133,166.58,-1200.6,-2270600,279150000,0,0,0,2.4656,0.0007467,-0.0749,4.47E-05,,,,,,,,,,,,,,
-orthocorundum,,-1634950,33.93,3.123,119.38,774.8,-6509100,422880000,0,0,0,2.4656,0.0007467,-0.0749,4.47E-05,,,,,,,,,,,,,,
-pyrope,,-6284740,268.8,11.311,590.9,-2827,-13320800,1260330000,0,0,0,2.2519,0.00037044,-0.0576,4.42E-05,,,,,,,,,,,,,,
-rutile,,-944750,50.88,1.883,77.84,0,-3367800,402940000,0,0,0,2.5716,0.00015409,-0.0454,5.85E-05,,,,,,,,,,,,,,
-spinel,,-2302160,80.37,3.978,244.67,-2004,0,0,0,0,0,2.1691,0.00050528,-0.0489,0,,,,,,,,,,,,,,
-talc,,-5898960,261.21,13.61,664.11,-5187.2,-2147200,-327370000,0,0,0,2.5616,0,-0.1847,0.0005878,,,,,,,,,,,,,,
Copied: pkg/CHNOSZ/inst/tinytest/test-Berman.R (from rev 692, pkg/CHNOSZ/inst/tinytest/test-berman.R)
===================================================================
--- pkg/CHNOSZ/inst/tinytest/test-Berman.R (rev 0)
+++ pkg/CHNOSZ/inst/tinytest/test-Berman.R 2022-02-03 12:13:18 UTC (rev 693)
@@ -0,0 +1,118 @@
+# test-Berman.R 20171001
+
+# Load default settings for CHNOSZ
+reset()
+
+# Make sure all Berman minerals are listed with units of J in OBIGT 20220203
+info <- "Berman minerals are listed with units of J in OBIGT"
+file <- system.file("extdata/OBIGT/Berman_cr.csv", package = "CHNOSZ")
+dat <- read.csv(file)
+expect_true(all(dat$E_units == "J"), info = info)
+
+# The maximum absolute pairwise difference between x and y
+maxdiff <- function(x, y) max(abs(y - x))
+
+info <- "high-T,P calculated properties are similar to precalculated ones"
+# Reference values for G were taken from the spreadsheet Berman_Gibbs_Free_Energies.xlsx
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 693
More information about the CHNOSZ-commits
mailing list