[CHNOSZ-commits] r726 - in pkg/CHNOSZ: . R demo inst man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Jun 4 07:29:26 CEST 2022
Author: jedick
Date: 2022-06-04 07:29:26 +0200 (Sat, 04 Jun 2022)
New Revision: 726
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/diagram.R
pkg/CHNOSZ/demo/contour.R
pkg/CHNOSZ/inst/NEWS.Rd
pkg/CHNOSZ/man/affinity_rank.Rd
pkg/CHNOSZ/vignettes/multi-metal.Rmd
Log:
Improve handling of line style for different species in predominance diagrams
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2022-05-06 10:19:39 UTC (rev 725)
+++ pkg/CHNOSZ/DESCRIPTION 2022-06-04 05:29:26 UTC (rev 726)
@@ -1,6 +1,6 @@
-Date: 2022-05-06
+Date: 2022-06-04
Package: CHNOSZ
-Version: 1.9.9-18
+Version: 1.9.9-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/R/diagram.R
===================================================================
--- pkg/CHNOSZ/R/diagram.R 2022-05-06 10:19:39 UTC (rev 725)
+++ pkg/CHNOSZ/R/diagram.R 2022-06-04 05:29:26 UTC (rev 726)
@@ -124,9 +124,12 @@
})
names(plotvals) <- names(eout$values)
plotvar <- eout$property
- # we change 'A' to 'A/(2.303RT)' so the axis label is made correctly
- # 20171027 use parentheses to avoid ambiguity about order of operations
- if(plotvar=="A") {
+ if(efun == "affinity_rank") {
+ plotvar <- "affinity_rank"
+ message(paste("diagram: plotting average affinity ranking for", length(plotvals), "groups"))
+ } else if(plotvar=="A") {
+ # we change 'A' to 'A/(2.303RT)' so the axis label is made correctly
+ # 20171027 use parentheses to avoid ambiguity about order of operations
plotvar <- "A/(2.303RT)"
if(nd==2 & type=="auto") message("diagram: using maximum affinity method for 2-D diagram")
else if(nd==2 & type=="saturation") message("diagram: plotting saturation lines for 2-D diagram")
@@ -200,7 +203,7 @@
## identify predominant species
predominant <- NA
H2O.predominant <- NULL
- if(plotvar %in% c("loga.equil", "alpha", "A/(2.303RT)") & type!="saturation") {
+ if(plotvar %in% c("loga.equil", "alpha", "A/(2.303RT)", "affinity_rank") & type!="saturation") {
pv <- plotvals
# some additional steps for affinity values, but not for equilibrated activities
if(eout.is.aout) {
@@ -349,6 +352,7 @@
if(missing(xlab)) xlab <- axis.label(eout$vars[1], basis=eout$basis, molality=molality)
if(missing(ylab)) {
ylab <- axis.label(plotvar, units="", molality=molality)
+ if(plotvar == "affinity_rank") ylab <- "Average affinity ranking"
# use ppb, ppm, ppt (or log ppb etc.) for converted values of solubility 20190526
if(grepl("solubility.", eout$fun, fixed=TRUE)) {
ylab <- strsplit(eout$fun, ".", fixed=TRUE)[[1]][2]
@@ -558,12 +562,9 @@
# DEFAULT method: loop over species
for(i in 1:(length(zvals)-1)) {
# get the "z" values
- z <- predominant
- # assign values to get one contour line between this species and all others
- i0 <- z==zvals[i]
- i1 <- z!=zvals[i]
- z[i0] <- 0
- z[i1] <- 1
+ # Use + 0 trick to convert T/F to 1/0 20220524
+ z <- (predominant == zvals[i]) + 0
+ z[is.na(z)] <- 0
# use contourLines() instead of contour() in order to get line coordinates 20181029
cLines <- contourLines(xs, ys, z, levels = 0.5)
if(length(cLines) > 0) {
@@ -570,11 +571,11 @@
# loop in case contourLines returns multiple lines
for(k in 1:length(cLines)) {
# draw the lines
- lines(cLines[[k]][2:3], lty = lty[i], col = col[i], lwd = lwd[i])
+ lines(cLines[[k]][2:3], lty = lty[zvals[i]], col = col[zvals[i]], lwd = lwd[zvals[i]])
}
}
# mask species to prevent double-plotting contour lines
- predominant[i0] <- NA
+ predominant[z == zvals[i]] <- NA
}
} else {
@@ -595,7 +596,7 @@
# loop in case contourLines returns multiple lines
for(k in 1:length(cLines)) {
# draw the lines
- mylty <- lty[i]
+ mylty <- lty[zvals[i]]
if(!is.null(lty.cr)) {
# use lty.cr for cr-cr boundaries 20190530
if(all(grepl("cr", eout$species$state[c(zvals[i], zvals[j])]))) mylty <- lty.cr
@@ -604,7 +605,7 @@
# use lty.aq for aq-aq boundaries 20190531
if(all(grepl("aq", eout$species$state[c(zvals[i], zvals[j])]))) mylty <- lty.aq
}
- lines(cLines[[k]][2:3], lty = mylty, col = col[i], lwd = lwd[i])
+ lines(cLines[[k]][2:3], lty = mylty, col = col[zvals[i]], lwd = lwd[zvals[i]])
}
}
}
Modified: pkg/CHNOSZ/demo/contour.R
===================================================================
--- pkg/CHNOSZ/demo/contour.R 2022-05-06 10:19:39 UTC (rev 725)
+++ pkg/CHNOSZ/demo/contour.R 2022-06-04 05:29:26 UTC (rev 726)
@@ -17,7 +17,7 @@
basis(c("Au", "Cl-", "H2S", "H2O", "oxygen", "H+"))
iaq <- info(c("Au(HS)2-", "AuHS", "AuOH", "AuCl2-"))
species(iaq)
-# This get us close to total S = 0.01 m
+# This gets us close to total S = 0.01 m
basis("H2S", -2)
# Calculate solution composition for 1 mol/kg NaCl
NaCl <- NaCl(T = T, P = P, m_tot=1)
Modified: pkg/CHNOSZ/inst/NEWS.Rd
===================================================================
--- pkg/CHNOSZ/inst/NEWS.Rd 2022-05-06 10:19:39 UTC (rev 725)
+++ pkg/CHNOSZ/inst/NEWS.Rd 2022-06-04 05:29:26 UTC (rev 726)
@@ -10,7 +10,7 @@
\newcommand{\s}{\ifelse{latex}{\eqn{_{#1}}}{\ifelse{html}{\out{<sub>#1</sub>}}{#1}}}
\newcommand{\S}{\ifelse{latex}{\eqn{^{#1}}}{\ifelse{html}{\out{<sup>#1</sup>}}{^#1}}}
-\section{Changes in CHNOSZ version 1.9.9-17 (2022-04-18)}{
+\section{Changes in CHNOSZ version 1.9.9-19 (2022-06-04)}{
\subsection{MAJOR CHANGE}{
\itemize{
@@ -85,6 +85,8 @@
\item Change default for \code{affinity()} to \code{loga.protein = 0}
(was -3).
+ \item Change license from GPL (>= 2) to GPL-3.
+
}
}
@@ -111,6 +113,17 @@
\item In \strong{multi-metal.Rmd}, work around absence of
\code{hcl.colors()} in \R < 3.6.0.
+ \item \code{diagram()}: Improve handling of length > 2 \code{col},
+ \code{lty} and \code{lwd} arguments for predominance diagrams. Because
+ field boundaries are drawn for the first species, then removing that
+ species before drawing boundaries for the second species, etc., species
+ for which different-colored or -styled lines are desired should be placed
+ at the top of the species list. This change makes possible one less call
+ to \code{diagram()} in the Mosaic Stacking sections in
+ \strong{multi-metal.Rmd}, and now the entire chalcopyrite field in Mosaic
+ Stacking 2 is bounded by a thick orange line, instead of just the
+ border shared with bornite.
+
}
}
Modified: pkg/CHNOSZ/man/affinity_rank.Rd
===================================================================
--- pkg/CHNOSZ/man/affinity_rank.Rd 2022-05-06 10:19:39 UTC (rev 725)
+++ pkg/CHNOSZ/man/affinity_rank.Rd 2022-06-04 05:29:26 UTC (rev 726)
@@ -46,7 +46,7 @@
nspecies <- sapply(groups, sum)
names <- paste0(names(groups), " (", nspecies, ")")
diagram(arank, fill = "terrain", font = 2, names = names, format.names = FALSE)
-title("Chemical affinity ranking of Rubisco proteins")
+title("Average affinity ranking of Rubisco proteins")
}
\concept{Extended workflow}
Modified: pkg/CHNOSZ/vignettes/multi-metal.Rmd
===================================================================
--- pkg/CHNOSZ/vignettes/multi-metal.Rmd 2022-05-06 10:19:39 UTC (rev 725)
+++ pkg/CHNOSZ/vignettes/multi-metal.Rmd 2022-06-04 05:29:26 UTC (rev 726)
@@ -69,7 +69,7 @@
## Resolution settings
# Change this to TRUE to make high-resolution plots
# (default is FALSE to save time in CRAN checks)
-hires <- FALSE
+hires <- TRUE
res1.lo <- 150
res1.hi <- 256
res1 <- ifelse(hires, res1.hi, res1.lo)
@@ -385,7 +385,7 @@
# To get the limits, convert range of affinities to eV/atom
arange <- rev(range(aFeVO4_vs_stable))
# This gets us to J/mol
-Jrange <- convert(convert(arange, "G"), "J")
+Jrange <- convert(arange, "G")
# And to eV/atom
eVrange <- Jrange / 1.602176634e-19 / 6.02214076e23 / 6
ylim <- formatC(eVrange, digits = 3, format = "f")
@@ -487,9 +487,10 @@
species(c(FeCu.cr, Cu.cr))
mFeCu <- mosaic(list(S.aq, Fe.cr), pH = pH, O2 = O2,
T = T, stable = list(NULL, dFe$predominant))
-diagram(mFeCu$A.species, add = TRUE, col = 2, col.names = 2, bold = TRUE, names = FeCu.abbrv, dy = c(0, 0, 0, 0, 0, 1, 0))
-col <- c("#FF8C00", rep(NA, 6))
-diagram(mFeCu$A.species, add = TRUE, col = col, lwd = 2, col.names = col, bold = TRUE, names = FeCu.abbrv)
+col <- c("#FF8C00", rep(2, 6))
+lwd <- c(2, 1, 1, 1, 1, 1, 1)
+dy = c(0, 0, 0, 0, 0, 1, 0)
+diagram(mFeCu$A.species, add = TRUE, col = col, lwd = lwd, col.names = col, bold = TRUE, names = FeCu.abbrv, dy = dy)
TPS <- c(describe.property(c("T", "P"), c(T, "Psat")), expression(sum(S) == 0.01*m))
legend("topright", TPS, bty = "n")
title("Cu-Fe-S-O-H (minerals only)", font.main = 1)
@@ -506,7 +507,7 @@
After that we add the legend and title.
-```{r stack1_2, echo=5:16, results = "hide", message = FALSE, fig.width = 6, fig.height = 5, out.width = "75%", fig.align = "center", pngquant = pngquant}
+```{r stack1_2, echo=5:17, results = "hide", message = FALSE, fig.width = 6, fig.height = 5, out.width = "75%", fig.align = "center", pngquant = pngquant}
```
This diagram has a distinctive chalcopyrite "hook" surrounded by a thin bornite field.
@@ -547,11 +548,6 @@
NaClexpr <- as.expression(bquote(NaCl == .(m_NaCl)*m))
aqexpr <- as.expression(bquote("("*aq*")"[italic(i)] == 10^.(logm_aq)*m))
-
-
-
-
-
# Setup basis species
basis(c("Cu+", "pyrite", "H2S", "oxygen", "H2O", "H+", "Cl-"))
basis("H2S", logmS)
@@ -570,22 +566,21 @@
# Add aqueous species 20210220
species(iCu.aq, logm_aq, add = TRUE)
-
-
-
-
-
-# TODO: limitation in mosaic() when using both solid and aqueous basis species (i.e. c(Fe.cr, Fe.aq)):
+# NOTE: limitation in mosaic() when using both solid and aqueous basis species (i.e. c(Fe.cr, Fe.aq)):
# Only the activity of the first-defined basis species is used for all basis species.
# The first-defined basis species is pyrite (logact = 0), but logact < 0 for the aq Fe species.
# Workaround: Adjust standard Gibbs energies of aq Fe species to virtually change their activities
-DG <- convert(-logm_aq, "G", T = convert(T, "K"))
+DG_J <- convert(-logm_aq, "G", T = convert(T, "K"))
+# We should use calories here because the database values are in calories 20220604
+stopifnot(all(info(iFe.aq)$E_units == "cal"))
+DG_cal <- convert(DG_J, "cal")
G.orig <- info(iFe.aq)$G
-G.new <- G.orig + DG
+G.new <- G.orig + DG_cal
mod.OBIGT(iFe.aq, G = G.new)
+## Mosaic with all Fe species as basis species
#mFeCu <- mosaic(list(S.aq, c(Fe.cr, Fe.aq)), pH = pH, O2 = O2, T = T, stable = list(NULL, dFe$predominant))
-# Use only predominant species as basis species (to speed up calculation) 20210224
+# Use only predominant Fe species as basis species (to speed up calculation) 20210224
predom <- dFe$predominant
ipredom <- sort(unique(as.numeric(predom)))
for(i in seq_along(ipredom)) predom[dFe$predominant == ipredom[i]] <- i
@@ -610,14 +605,15 @@
dy[names == "CuCl2-"] <- 2
cex[names == "Bn"] <- 0.8
srt[names == "Bn"] <- 85
-diagram(mFeCu$A.species, add = TRUE, col = 2, col.names = 2, names = names, bold = bold, dx = dx, dy = dy, cex.names = cex, srt = srt)
+# Highlight Ccp field
+col.names <- col <- rep(2, nrow(mFeCu$A.species$species))
+col[1] <- "#FF8C00"
+col.names[1] <- "#FF8C00"
+lwd <- rep(1, nrow(mFeCu$A.species$species))
+lwd[1] <- 2
+diagram(mFeCu$A.species, add = TRUE, lwd = lwd, col = col, col.names = col.names, names = names, bold = bold, dx = dx, dy = dy, cex.names = cex, srt = srt)
# Add second Cu label
text(12.3, -47, "Cu", col = 2, font = 2)
-# Highlight Bn-Ccp reaction
-col.names <- col <- rep(NA, nrow(mFeCu$A.species$species))
-col[3] <- "#FF8C00"
-col.names[1] <- "#FF8C00"
-diagram(mFeCu$A.species, add = TRUE, col = col, lwd = 2, col.names = col.names, bold = TRUE, names = names, fill = NA)
# Plot the Fe-system lines and names "on top" so they are not covered by fill colors
diagram(mFe$A.bases, add = TRUE, lty = 2, col = 4, names = FALSE, fill = NA)
@@ -843,16 +839,21 @@
# Sverjensky (1987) used Helgeson (1969) value, which is ca. -5.2
dlogK <- logK - -5.2
# Calculate the difference in ΔG° corresponding to this logK difference
-dG <- convert(dlogK, "G", T = convert(T, "K"))
+dG_J <- convert(dlogK, "G", T = convert(T, "K"))
+# We should use calories here because the database values are in calories 20220604
+stopifnot(info(info("CuCl2-"))$E_units == "cal")
+dG_cal <- convert(dG_J, "cal")
# Apply this difference to the CuCl2- entry in OBIGT
-newG <- info(info("CuCl2-"))$G + dG
+newG <- info(info("CuCl2-"))$G + dG_cal
mod.OBIGT("CuCl2-", G = newG)
# Do the same thing for CuCl3-2
logK <- subcrt(c("CuCl3-2", "Cu+", "Cl-"), c(-1, 1, 3), T = T)$out$logK
dlogK <- logK - -5.6
-dG <- convert(dlogK, "G", T = convert(T, "K"))
-newG <- info(info("CuCl3-2"))$G + dG
+dG_J <- convert(dlogK, "G", T = convert(T, "K"))
+stopifnot(info(info("CuCl3-2"))$E_units == "cal")
+dG_cal <- convert(dG_J, "cal")
+newG <- info(info("CuCl3-2"))$G + dG_cal
mod.OBIGT("CuCl3-2", G = newG)
# DIAGRAM 2
More information about the CHNOSZ-commits
mailing list