[CHNOSZ-commits] r762 - in pkg/CHNOSZ: . R demo inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Dec 13 13:52:21 CET 2022
Author: jedick
Date: 2022-12-13 13:52:21 +0100 (Tue, 13 Dec 2022)
New Revision: 762
Added:
pkg/CHNOSZ/demo/yttrium.R
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/NaCl.R
pkg/CHNOSZ/demo/00Index
pkg/CHNOSZ/demo/sphalerite.R
pkg/CHNOSZ/inst/NEWS.Rd
pkg/CHNOSZ/man/NaCl.Rd
pkg/CHNOSZ/man/examples.Rd
Log:
Add demo/yttrium.R
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2022-12-11 04:21:41 UTC (rev 761)
+++ pkg/CHNOSZ/DESCRIPTION 2022-12-13 12:52:21 UTC (rev 762)
@@ -1,6 +1,6 @@
-Date: 2022-12-11
+Date: 2022-12-13
Package: CHNOSZ
-Version: 1.9.9-53
+Version: 1.9.9-54
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/NaCl.R
===================================================================
--- pkg/CHNOSZ/R/NaCl.R 2022-12-11 04:21:41 UTC (rev 761)
+++ pkg/CHNOSZ/R/NaCl.R 2022-12-13 12:52:21 UTC (rev 762)
@@ -11,7 +11,7 @@
# Store existing thermo data frame
thermo <- get("thermo", CHNOSZ)
# Get length of longest variable
- nTP <- max(length(T), length(P))
+ nTP <- max(length(T), length(P), length(pH))
pH.arg <- pH
pH <- rep(pH, length.out = nTP)
Modified: pkg/CHNOSZ/demo/00Index
===================================================================
--- pkg/CHNOSZ/demo/00Index 2022-12-11 04:21:41 UTC (rev 761)
+++ pkg/CHNOSZ/demo/00Index 2022-12-13 12:52:21 UTC (rev 762)
@@ -30,3 +30,4 @@
Pourbaix Eh-pH diagram for Fe-O-H with equisolubility lines
E_coli Gibbs energy of biomass synthesis in E. coli
rank.affinity Affinity ranking for proteins in yeast nutrient limitation
+yttrium logB.to.OBIGT fits at 800 and 1000 bar and Y speciation in NaCl solution at varying pH
Modified: pkg/CHNOSZ/demo/sphalerite.R
===================================================================
--- pkg/CHNOSZ/demo/sphalerite.R 2022-12-11 04:21:41 UTC (rev 761)
+++ pkg/CHNOSZ/demo/sphalerite.R 2022-12-13 12:52:21 UTC (rev 762)
@@ -6,18 +6,20 @@
# Set up chemical system
basis(c("Zn+2", "Cl-", "H2S", "H2O", "O2", "H+"))
+basis("H2S", log10(0.05))
species("ZnS")
iaq <- retrieve("Zn", c("O", "H", "Cl", "S"), "aq")
-# a function to make a single plot
+# A function to make a single plot
plotfun <- function(T = 400, P = 500, m_tot = 0.1, pHmin = 4, logppmmax = 3) {
- # Calculate NaCl speciation from simplified model
- NaCl <- NaCl(m_tot = m_tot, T = T, P = P)
- basis("Cl-", log10(NaCl$m_Cl))
- basis("H2S", log10(0.05))
+ # Get pH values
+ res <- 100
+ pH <- seq(pHmin, 10, length.out = res)
+ # Calculate speciation in NaCl-H2O system at given pH
+ NaCl <- NaCl(m_tot = m_tot, T = T, P = P, pH = pH)
# Calculate solubility with mosaic (triggered by bases argument) to account for HS- and H2S speciation
- s <- solubility(iaq, bases = c("H2S", "HS-"), pH = c(pHmin, 10), T = T, P = P, IS = NaCl$IS)
+ s <- solubility(iaq, bases = c("H2S", "HS-"), pH = pH, "Cl-" = log10(NaCl$m_Cl), T = T, P = P, IS = NaCl$IS)
# Convert log activity to log ppm
sp <- convert(s, "logppm")
@@ -42,15 +44,15 @@
# A function to make a page of plots
pagefun <- function() {
- # set the values of temperature, pressure, and total NaCl
+ # Set the values of temperature, pressure, and total NaCl
T <- c(400, 400, 250, 250, 100, 100)
- # use a list to be able to mix numeric and character values for P
+ # Use a list to be able to mix numeric and character values for P
P <- list(500, 500, "Psat", "Psat", "Psat", "Psat")
m_tot <- c(0.1, 1, 0.1, 1, 0.1, 1)
- # the plots have differing limits
+ # The plots have differing limits
pHmin <- c(4, 4, 2, 2, 2, 2)
logppmmax <- c(3, 3, 2, 2, 0, 0)
- # make the plots
+ # Make the plots
par(mfrow = c(3, 2))
for(i in 1:6) plotfun(T = T[i], P = P[[i]], m_tot = m_tot[i], pHmin = pHmin[i], logppmmax = logppmmax[i])
}
@@ -59,7 +61,7 @@
pngfun <- function() {
png("sphalerite.png", width = 1000, height = 1200, pointsize = 24)
pagefun()
- # add an overall title
+ # Add an overall title
par(xpd = NA)
text(1, 14, "Solubility of sphalerite, after Akinfiev and Tagirov, 2014, Fig. 13", cex = 1.5)
par(xpd = FALSE)
Added: pkg/CHNOSZ/demo/yttrium.R
===================================================================
--- pkg/CHNOSZ/demo/yttrium.R (rev 0)
+++ pkg/CHNOSZ/demo/yttrium.R 2022-12-13 12:52:21 UTC (rev 762)
@@ -0,0 +1,133 @@
+# CHNOSZ/demo/yttrium.R
+# Make diagrams similar to Figure 11 of Guan et al. (2020)
+# https://doi.org/10.1016/j.gca.2020.04.015
+# 20221122 jmd version 1
+
+library(CHNOSZ)
+
+# Function to add Y(III)-Cl species using logB values from Table 9 of Guan et al. (2020)
+add.Y.species <- function(P, plot.it = FALSE) {
+ # Temperature
+ T <- c(25, seq(50, 500, 50))
+ if(P == 800) {
+ logB <- list(
+ c(1.04, 0.48, 0.13, 0.43, 1.12, 2.06, 3.21, 4.58, 6.26, 8.39, 11.14),
+ c(-9.14, -7.34, -4.14, -1.37, 1.12, 3.46, 5.78, 8.24, 11.05, 14.48, 18.77),
+ c(-14, -11.48, -7.06, -3.25, 0.14, 3.32, 6.45, 9.74, 13.47, 18.01, 23.65),
+ c(-15.94, -13.2, -8.39, -4.27, -0.61, 2.79, 6.15, 9.67, 13.63, 18.45, 24.41)
+ )
+ } else if(P == 1000) {
+ logB <- list(
+ c(1.13, 0.54, 0.16, 0.43, 1.09, 2, 3.1, 4.39, 5.9, 7.69, 9.85),
+ c(-9.33, -7.51, -4.3, -1.55, 0.9, 3.18, 5.4, 7.68, 10.15, 12.94, 16.16),
+ c(-14.24, -11.71, -7.27, -3.49, -0.14, 2.95, 5.95, 9.01, 12.3, 16, 20.25),
+ c(-16.19, -13.43, -8.62, -4.52, -0.91, 2.41, 5.62, 8.89, 12.4, 16.33, 20.84)
+ )
+ } else stop("logB values for P =", P, "are not available here")
+ # Define species and coefficients in formation reactions
+ species <- list(
+ c("Y+3", "Cl-", "YCl+2"),
+ c("Y+3", "Cl-", "YCl2+"),
+ c("Y+3", "Cl-", "YCl3"),
+ c("Y+3", "Cl-", "YCl4-")
+ )
+ coeffs <- list(
+ c(-1, -1, 1),
+ c(-1, -2, 1),
+ c(-1, -3, 1),
+ c(-1, -4, 1)
+ )
+ # Fit the formation constants to thermodynamic parameters and add them to OBIGT
+ for(i in 1:4) logB.to.OBIGT(logB[[i]], species[[i]], coeffs[[i]], T = T, P = P, tolerance = 0.6, npar = 5)
+ # Plot the given and fitted values
+ if(plot.it) {
+ par(mfrow = c(2, 2))
+ for(i in 1:4) {
+ sres <- subcrt(species[[i]], coeffs[[i]], T = T, P = P)
+ plot(T, sres$out$logK, type = "l", xlab = axis.label("T"), ylab = axis.label("logK"))
+ points(T, logB[[i]], pch = 19)
+ title(describe.reaction(sres$reaction))
+ if(i == 1) legend("topleft", c("Guan et al. (2020)", "logB.to.OBIGT()"), pch = c(19, NA), lty = c(NA, 1))
+ legend("bottomright", paste(P, "bar"), bty = "n")
+ }
+ }
+}
+
+# Function to plot distribution of Y(III) chloride species at T and P
+Y_Cl <- function() {
+
+ # Define total molality of NaCl
+ # Start at 0.1 because we can't use 0 in the logarithmic value needed for affinity()
+ mNaCl <- seq(0.1, 4.9, 0.2)
+
+ # Define T, P, and pH values
+ Ts <- c(200, 350, 500)
+ Ps <- c(800, 800, 1000)
+ pHs <- c(3, 0.3)
+
+ # Setup plot
+ par(mfrow = c(3, 2))
+ par(cex = 0.9)
+
+ # Loop over T and P
+ for(i in 1:3) {
+ T <- Ts[i]
+ P <- Ps[i]
+
+ # Add new species
+ add.Y.species(P)
+ # Setup chemical system
+ basis(c("Y+3", "Cl-", "e-"))
+ species(c("Y+3", "YCl+2", "YCl2+", "YCl3", "YCl4-"))
+
+ # Loop over pH
+ for(j in 1:2) {
+ pH <- pHs[j]
+
+ # Calculate molality of Cl- and ionic strength
+ NaCl_props <- suppressMessages(lapply(mNaCl, NaCl, T = T, P = P, pH = pH))
+ # Turn the list into a data frame
+ NaCl_props <- do.call(rbind, lapply(NaCl_props, as.data.frame))
+ # Calculate affinity of formation reactions
+ a <- affinity("Cl-" = log10(NaCl_props$m_Cl), IS = NaCl_props$IS, T = T, P = P)
+ # Calculate species distribution for total Y(III) equal to 0.01 m
+ m_Y <- 0.01
+ e <- equilibrate(a, loga.balance = log10(m_Y))
+ # Make x-axis show total m(NaCl) instead of logm(Cl-) 20221208
+ e$vals[[1]] <- mNaCl
+ # Set colors similar to Guan et al. (2020)
+ col <- 2:6
+ # Only label lines above 1/20 = 0.05 mol fraction
+ mol <- 10^do.call(cbind, lapply(e$loga.equil, as.data.frame))
+ molfrac <- mol / rowSums(mol)
+ ilab <- apply(molfrac > 0.05, 2, any)
+ names <- e$species$name
+ names[!ilab] <- ""
+ d <- diagram(e, alpha = TRUE, xlim = c(0, 5), ylim = c(0, 1),
+ xlab = expr.species("NaCl", molality = TRUE), ylab = "Fraction of Y-Cl species",
+ names = names, col = col, lty = 1, lwd = 2, mar = c(3, 3.5, 2.5, 3.5))
+ # Calculate and plot coordination number
+ CN <- 1 * d$plotvals[[2]] + 2 * d$plotvals[[3]] + 3 * d$plotvals[[4]] + 4 * d$plotvals[[5]]
+ # Rescale to y-axis limits [0, 1]
+ CN_scaled <- CN / 4
+ lines(mNaCl, CN_scaled, lty = 3, lwd = 3, col = "darkorange")
+ # Add ticks for CN
+ axis(4, seq(0, 1, 0.25), labels = 0:4, lwd.ticks = 2, tcl = -0.5, mgp = c(1.7, 0.8, 0))
+ mtext("Cl coordination number", 4, 1.7, las = 0, cex = par("cex"))
+ # Make title
+ lab <- lTP(T, P)
+ lab[[4]] <- bquote(pH == .(pH))
+ title(lab, cex.main = 1)
+
+ }
+
+ if(i==1) mtext("After Guan et al. (2020)", line = 0.8, adj = -2.7, cex = 0.9, font = 2)
+ }
+}
+
+# Run the functions to amke plots for the demo
+opar <- par(no.readonly = TRUE)
+add.Y.species(800, plot.it = TRUE)
+add.Y.species(1000, plot.it = TRUE)
+Y_Cl()
+par(opar)
Modified: pkg/CHNOSZ/inst/NEWS.Rd
===================================================================
--- pkg/CHNOSZ/inst/NEWS.Rd 2022-12-11 04:21:41 UTC (rev 761)
+++ pkg/CHNOSZ/inst/NEWS.Rd 2022-12-13 12:52:21 UTC (rev 762)
@@ -12,7 +12,7 @@
% links to vignettes 20220723
\newcommand{\viglink}{\ifelse{html}{\out{<a href="../CHNOSZ/doc/#1.html"><strong>#1.Rmd</strong></a>}}{\bold{#1.Rmd}}}
-\section{Changes in CHNOSZ version 1.9.9-53 (2022-12-11)}{
+\section{Changes in CHNOSZ version 1.9.9-54 (2022-12-13)}{
\subsection{MAJOR USER-VISIBLE CHANGES}{
\itemize{
@@ -129,6 +129,21 @@
}
}
+ \subsection{NEW FEATURES: OTHER}{
+ \itemize{
+
+ \item \code{NaCl()} has been rewritten to include pH dependence (i.e.
+ formation of HCl as well as NaCl) by using \code{affinity()} and
+ \code{equilibrate()} instead of algebraic equations.
+
+ \item Add demo \code{yttrium.R} to show speciation of Y-Cl complexes as a
+ function of NaCl concentration, pH, \emph{T}, and \emph{P}, after
+ \href{https://doi.org/10.1016/j.gca.2020.04.015}{Guan et al. (2020)}.
+
+ }
+ }
+
+
\subsection{REMOVED FEATURES}{
\itemize{
@@ -205,10 +220,6 @@
\item \code{seq2aa()} now has the \strong{sequence} argument first and a
default of NA for \strong{protein} (the protein name).
- \item \code{NaCl()} has been rewritten to include pH dependence (i.e.
- formation of HCl as well as NaCl) by using \code{affinity()} and
- \code{equilibrate()} instead of algebraic equations.
-
\item Remove unused \samp{cutoff} value in \code{thermo()$opt}.
}
Modified: pkg/CHNOSZ/man/NaCl.Rd
===================================================================
--- pkg/CHNOSZ/man/NaCl.Rd 2022-12-11 04:21:41 UTC (rev 761)
+++ pkg/CHNOSZ/man/NaCl.Rd 2022-12-13 12:52:21 UTC (rev 762)
@@ -47,6 +47,7 @@
\seealso{
This function is used in a few demos (\code{demo("contour")}, \code{demo("gold")}, \code{demo("minsol")}, \code{demo("sphalerite")}).
+\code{demo("yttrium")} uses the \code{pH} argument.
}
\examples{\dontshow{reset()}
@@ -53,14 +54,14 @@
# Ionic strength calculated with HCh version 3.7 (Shvarov and Bastrakov, 1999)
# at 1000 bar, 100, 200, and 300 degress C, and 1 to 6 molal NaCl
m.HCh <- 1:6
-IS.HCh <- list(`100`=c(0.992, 1.969, 2.926, 3.858, 4.758, 5.619),
- `300`=c(0.807, 1.499, 2.136, 2.739, 3.317, 3.875),
- `500`=c(0.311, 0.590, 0.861, 1.125, 1.385, 1.642))
+IS.HCh <- list(`100` = c(0.992, 1.969, 2.926, 3.858, 4.758, 5.619),
+ `300` = c(0.807, 1.499, 2.136, 2.739, 3.317, 3.875),
+ `500` = c(0.311, 0.590, 0.861, 1.125, 1.385, 1.642))
# Total molality in the calculation with NaCl()
m_tot <- seq(1, 6, 0.5)
N <- length(m_tot)
# Where we'll put the calculated values
-IS.calc <- data.frame(`100`=numeric(N), `300`=numeric(N), `500`=numeric(N))
+IS.calc <- data.frame(`100` = numeric(N), `300` = numeric(N), `500` = numeric(N))
# NaCl() is *not* vectorized over m_tot, so we use a loop here
for(i in 1:length(m_tot)) {
NaCl.out <- NaCl(m_tot[i], c(100, 300, 500), P = 1000)
@@ -68,17 +69,17 @@
}
# Plot ionic strength from HCh and NaCl() as points and lines
col <- c("black", "red", "orange")
-plot(c(1,6), c(0,6), xlab="NaCl (mol/kg)", ylab=axis.label("IS"), type="n")
+plot(c(1,6), c(0,6), xlab = "NaCl (mol/kg)", ylab = axis.label("IS"), type = "n")
for(i in 1:3) {
# NOTE: the differences are probably mostly due to different models
# for the properties of NaCl(aq) (HCh: B.Ryhzenko model;
# CHONSZ: revised HKF with parameters from Shock et al., 1997)
- points(m.HCh, IS.HCh[[i]], col=col[i])
- lines(m_tot, IS.calc[, i], col=col[i])
+ points(m.HCh, IS.HCh[[i]], col = col[i])
+ lines(m_tot, IS.calc[, i], col = col[i])
}
# Add legend and title
dprop <- describe.property(rep("T", 3), c(100, 300, 500))
-legend("topleft", dprop, lty=1, pch=1, col=col)
+legend("topleft", dprop, lty = 1, pch = 1, col = col)
title(main="H2O + NaCl; HCh (points) and 'NaCl()' (lines)")
}
Modified: pkg/CHNOSZ/man/examples.Rd
===================================================================
--- pkg/CHNOSZ/man/examples.Rd 2022-12-11 04:21:41 UTC (rev 761)
+++ pkg/CHNOSZ/man/examples.Rd 2022-12-13 12:52:21 UTC (rev 762)
@@ -62,6 +62,7 @@
\item{comproportionation}{Gibbs energy of sulfur comproportionation (Amend et al., 2020)}
\item{Pourbaix}{Eh-pH diagram for Fe-O-H with equisolubility lines (Pourbaix, 1974)}
\item{E_coli}{Gibbs energy of biomass synthesis in \emph{E. coli} (LaRowe and Amend, 2016)}
+ \item{yttrium}{\code{\link{logB.to.OBIGT}} fits at 800 and 1000 bar and Y speciation in \code{\link{NaCl}} solution at varying pH (Guan et al., 2020)}
}
For either function, if \code{save.png} is TRUE, the plots are saved in \code{\link{png}} files whose names begin with the names of the help topics or demos.
@@ -113,6 +114,8 @@
Garrels, R. M. and Christ, C. L. (1965) \emph{Solutions, Minerals, and Equilibria}, Harper & Row, New York, 450 p. \url{https://www.worldcat.org/oclc/517586}
+Guan, Q., Mei, Y., Etschmann, B., Testemale, D., Louvel, M. and Brugger, J. (2020) Yttrium complexation and hydration in chloride-rich hydrothermal fluids: A combined \emph{ab initio} molecular dynamics and \emph{in situ} X-ray absorption spectroscopy study. \emph{Geochim. Cosmochim. Acta} \bold{281}, 168--189. \doi{10.1016/j.gca.2020.04.015}
+
Johnson, J. W., Oelkers, E. H. and Helgeson, H. C. (1992) SUPCRT92: A software package for calculating the standard molal thermodynamic properties of minerals, gases, aqueous species, and reactions from 1 to 5000 bar and 0 to 1000\degC. \emph{Comp. Geosci.} \bold{18}, 899--947. \doi{10.1016/0098-3004(92)90029-Q}
LaRowe, D. E. and Amend, J. P. (2016) The energetics of anabolism in natural settings. \emph{ISME J.} \bold{10}, 1285--1295. \doi{10.1038/ismej.2015.227}
More information about the CHNOSZ-commits
mailing list