Rename 'm_tot' argument in NaCl() to 'm_NaCl'

git-svn-id: svn://scm.r-forge.r-project.org/svnroot/chnosz/pkg/CHNOSZ@867 edb9625f-4e0d-4859-8d74-9fd3b1da38cb

DESCRIPTION

@@ -1,6 +1,6 @@
-Date: 2025-01-07
+Date: 2025-01-08
 Package: CHNOSZ
-Version: 2.1.0-38
+Version: 2.1.0-39
 Title: Thermodynamic Calculations and Diagrams for Geochemistry
 Authors@R: c(
     person("Jeffrey", "Dick", , "[email protected]", role = c("aut", "cre"),

R/NaCl.R

@@ -2,11 +2,11 @@
 # Calculate ionic strength and molalities of species given a total molality of NaCl
 # 20181102 First version (jmd)
 # 20181106 Use activity coefficients of Na+ and Nacl
-# 20221122 Make it work with m_tot = 0
+# 20221122 Make it work with m_NaCl = 0
 # 20221210 Rewritten to include pH dependence;
 #   uses affinity() and equilibrate() instead of algebraic equations
 
-NaCl <- function(m_tot = 1, T = 25, P = "Psat", pH = NA, attenuate = FALSE) {
+NaCl <- function(m_NaCl = 1, T = 25, P = "Psat", pH = NA, attenuate = FALSE) {
 
   # Store existing thermo data frame
   thermo <- get("thermo", CHNOSZ)
@@ -16,51 +16,51 @@ NaCl <- function(m_tot = 1, T = 25, P = "Psat", pH = NA, attenuate = FALSE) {
   pH <- rep(pH, length.out = nTP)
 
   # Start with complete dissociation into Na+ and Cl-,
-  # so ionic strength and molality of Na+ are equal to m_tot
-  m_Na <- IS <- m_tot
+  # so ionic strength and molality of Na+ are equal to m_NaCl
+  m_Naplus <- IS <- m_NaCl
   # Make them same length as T and P
   IS <- rep(IS, length.out = nTP)
-  m_Na <- rep(m_Na, length.out = nTP)
-  # Set tolerance for convergence to 1/100th of m_tot
-  tolerance <- m_tot / 100
+  m_Naplus <- rep(m_Naplus, length.out = nTP)
+  # Set tolerance for convergence to 1/100th of m_NaCl
+  tolerance <- m_NaCl / 100
 
-  # If m_tot is 0, return 0 for all variables 20221122
+  # If m_NaCl is 0, return 0 for all variables 20221122
   zeros <- rep(0, nTP)
-  if(m_tot == 0) return(list(IS = zeros, m_Na = zeros, m_Cl = zeros, m_NaCl = zeros, m_HCl = zeros))
+  if(m_NaCl == 0) return(list(IS = zeros, m_Naplus = zeros, m_Clminus = zeros, m_NaCl0 = zeros, m_HCl0 = zeros))
 
   maxiter <- 100
   for(i in 1:maxiter) {
     # Setup chemical system and calculate affinities
     if(identical(pH.arg, NA)) {
       basis(c("Na+", "Cl-", "e-"))
       species(c("Cl-", "NaCl"))
-      a <- suppressMessages(affinity(T = T, P = P, "Na+" = log10(m_Na), IS = IS, transect = TRUE))
+      a <- suppressMessages(affinity(T = T, P = P, "Na+" = log10(m_Naplus), IS = IS, transect = TRUE))
     } else {
       basis(c("Na+", "Cl-", "H+", "e-"))
       species(c("Cl-", "NaCl", "HCl"))
-      a <- suppressMessages(affinity(T = T, P = P, pH = pH, "Na+" = log10(m_Na), IS = IS, transect = TRUE))
+      a <- suppressMessages(affinity(T = T, P = P, pH = pH, "Na+" = log10(m_Naplus), IS = IS, transect = TRUE))
     }
     # Speciate Cl-
-    e <- suppressMessages(equilibrate(a, loga.balance = log10(m_tot)))
+    e <- suppressMessages(equilibrate(a, loga.balance = log10(m_NaCl)))
     # Get molality of each Cl-bearing species
-    m_Cl <- 10^e$loga.equil[[1]]
-    m_NaCl <- 10^e$loga.equil[[2]]
-    if(identical(pH.arg, NA)) m_HCl <- NA else m_HCl <- 10^e$loga.equil[[3]]
+    m_Clminus <- 10^e$loga.equil[[1]]
+    m_NaCl0 <- 10^e$loga.equil[[2]]
+    if(identical(pH.arg, NA)) m_HCl0 <- NA else m_HCl0 <- 10^e$loga.equil[[3]]
     # Store previous ionic strength and molality of Na+
     IS_prev <- IS
-    m_Na_prev <- m_Na
+    m_Naplus_prev <- m_Naplus
     # Calculate new molality of Na+ and deviation
-    m_Na <- m_tot - m_NaCl
+    m_Naplus <- m_NaCl - m_NaCl0
     # Only go halfway to avoid overshoot
-    if(attenuate) m_Na <- (m_Na_prev + m_Na) / 2
-    dm_Na <- m_Na - m_Na_prev
+    if(attenuate) m_Naplus <- (m_Naplus_prev + m_Naplus) / 2
+    dm_Naplus <- m_Naplus - m_Naplus_prev
     # Calculate ionic strength and deviation
-    IS <- (m_Na + m_Cl) / 2
+    IS <- (m_Naplus + m_Clminus) / 2
     dIS <- IS - IS_prev
     # Keep going until the deviations in ionic strength and molality of Na+ at all temperatures are less than tolerance
-    converged <- abs(dIS) < tolerance & abs(dm_Na) < tolerance
+    converged <- abs(dIS) < tolerance & abs(dm_Naplus) < tolerance
     if(all(converged)) {
-      # Add one step without attenuating the deviation of m_Na
+      # Add one step without attenuating the deviation of m_Naplus
       if(attenuate) attenuate <- FALSE else break
     }
   }
@@ -72,6 +72,6 @@ NaCl <- function(m_tot = 1, T = 25, P = "Psat", pH = NA, attenuate = FALSE) {
   # Restore thermo data frame
   assign("thermo", thermo, CHNOSZ)
   # Return the calculated values
-  list(IS = IS, m_Na = m_Na, m_Cl = m_Cl, m_NaCl = m_NaCl, m_HCl = m_HCl)
+  list(IS = IS, m_Naplus = m_Naplus, m_Clminus = m_Clminus, m_NaCl0 = m_NaCl0, m_HCl0 = m_HCl0)
 
 }

demo/contour.R

@@ -20,15 +20,18 @@ species(iaq)
 # This gets us close to total S = 0.01 m
 basis("H2S", -2)
 # Calculate solution composition for 1 mol/kg NaCl
-NaCl <- NaCl(m_tot = 1, T = T, P = P)
-basis("Cl-", log10(NaCl$m_Cl))
+NaCl <- NaCl(m_NaCl = 1, T = T, P = P)
+basis("Cl-", log10(NaCl$m_Clminus))
 # Calculate affinity with changing basis species
 bases <- c("H2S", "HS-", "HSO4-", "SO4-2")
 m <- mosaic(bases, pH = c(2, 10, res), O2 = c(-41, -29, res), T = T, P = P, IS = NaCl$IS, blend = blend)
 # Show predominance fields for S-species
-diagram(m$A.bases, col = 4, col.names = 4, lty = 2, italic = TRUE)
+diagram(m$A.bases, col = 8, col.names = 8, lty = 2, italic = TRUE)
 # Show predominance fields for Au-species
-diagram(m$A.species, add=TRUE, col = 2, col.names = 2, lwd = 2, bold = TRUE)
+# Plot multiple times to get deeper color
+diagram(m$A.species, add=TRUE, col = 7, col.names = 7, lwd = 2, bold = TRUE)
+diagram(m$A.species, add=TRUE, col = 7, col.names = 7, lwd = 2, bold = TRUE)
+diagram(m$A.species, add=TRUE, col = 7, col.names = 7, lwd = 2, bold = TRUE)
 
 # Calculate and plot solubility of Au (use named 'bases' argument to trigger mosaic calculation)
 species("Au")

demo/gold.R

@@ -102,12 +102,12 @@ Au_pH2 <- function() {
 # NaCl solution with total chloride equal to specified NaCl + KCl solution,
 # then estimate the molality of K+ in that solution 20181109
 chloride <- function(T, P, m_NaCl, m_KCl) {
-  NaCl <- NaCl(m_tot = m_NaCl + m_KCl, T = T, P = P)
+  NaCl <- NaCl(m_NaCl = m_NaCl + m_KCl, T = T, P = P)
   # Calculate logK of K+ + Cl- = KCl, adjusted for ionic strength
   logKadj <- subcrt(c("K+", "Cl-", "KCl"), c(-1, -1, 1), T = T, P = P, IS = NaCl$IS)$out$logK
   # What is the molality of K+ from 0.5 mol KCl in solution with 2 mol total Cl
-  m_K <- m_KCl / (10^logKadj * NaCl$m_Cl + 1)
-  list(IS = NaCl$IS, m_Cl = NaCl$m_Cl, m_K = m_K)
+  m_Kplus <- m_KCl / (10^logKadj * NaCl$m_Clminus + 1)
+  list(IS = NaCl$IS, m_Clminus = NaCl$m_Clminus, m_Kplus = m_Kplus)
 }
 
 # log(m_Au)-T diagram like Fig. 2B of Williams-Jones et al., 2009
@@ -123,7 +123,7 @@ Au_T1 <- function() {
   # Calculate solubility of gold
   species("Au")
   iaq <- info(c("Au(HS)2-", "AuHS", "AuOH", "AuCl2-"))
-  s <- solubility(iaq, T = seq(150, 550, 10), `Cl-` = log10(chl$m_Cl), `K+` = log10(chl$m_K), P = 1000, IS = chl$IS)
+  s <- solubility(iaq, T = seq(150, 550, 10), `Cl-` = log10(chl$m_Clminus), `K+` = log10(chl$m_Kplus), P = 1000, IS = chl$IS)
   # Make diagram and show total log molality
   diagram(s, type = "loga.equil", ylim = c(-10, -3), col = col, lwd = 2, lty = 1)
   diagram(s, add = TRUE, lwd = 3, lty = 2)
@@ -154,10 +154,10 @@ Au_T2 <- function() {
   # Calculate solubility of gold
   species("Au")
   iaq <- info(c("Au(HS)2-", "AuHS", "AuOH", "AuCl2-"))
-  s <- solubility(iaq, T = seq(150, 550, 10), `Cl-` = log10(chl$m_Cl), `K+` = log10(chl$m_K), P = 1000, IS = chl$IS)
+  s <- solubility(iaq, T = seq(150, 550, 10), `Cl-` = log10(chl$m_Clminus), `K+` = log10(chl$m_Kplus), P = 1000, IS = chl$IS)
 #  # Uncomment to calculate solubility considering speciation of sulfur
 #  bases <- c("H2S", "HS-", "SO4-2", "HSO4-")
-#  s <- solubility(iaq, bases = bases, T = seq(150, 550, 10), `Cl-` = log10(chl$m_Cl), `K+` = log10(chl$m_K), P = 1000, IS = chl$IS)
+#  s <- solubility(iaq, bases = bases, T = seq(150, 550, 10), `Cl-` = log10(chl$m_Clminus), `K+` = log10(chl$m_Kplus), P = 1000, IS = chl$IS)
   # Make diagram and show total log molality
   diagram(s, type = "loga.equil", ylim = c(-10, -3), col = col, lwd = 2, lty = 1)
   diagram(s, add = TRUE, lwd = 3, lty = 2)

demo/minsol.R

@@ -26,8 +26,8 @@ basis("H2S", log10(Stot))
 # Molality of NaCl
 mNaCl <- 1000 * wNaCl / (mass("NaCl") * (1 - wNaCl))
 # Estimate ionic strength and molality of Cl-
-sat <- NaCl(m_tot = mNaCl, T = T)
-basis("Cl-", log10(sat$m_Cl))
+NaCl <- NaCl(m_NaCl = mNaCl, T = T)
+basis("Cl-", log10(NaCl$m_Clminus))
 
 # Add minerals and aqueous species
 icr <- retrieve(metal, c("Cl", "S", "O"), state = "cr")
@@ -38,9 +38,9 @@ species(iaq, logm_metal, add = TRUE)
 
 # Calculate affinities and make diagram
 bases <- c("H2S", "HS-", "HSO4-", "SO4-2")
-m <- mosaic(bases, pH = pH, O2 = O2, T = T, P = P, IS = sat$IS)
+m <- mosaic(bases, pH = pH, O2 = O2, T = T, P = P, IS = NaCl$IS)
 d <- diagram(m$A.species, bold = TRUE)
-diagram(m$A.bases, add = TRUE, col = "slategray", lwd = 2, lty = 3, names = NA)
+diagram(m$A.bases, add = TRUE, col = 8, col.names = 8, lty = 3, italic = TRUE)
 title(bquote(log * italic(m)[.(metal)*"(aq) species"] == .(logm_metal)))
 label.figure("A")
 
@@ -62,14 +62,14 @@ par(xpd = FALSE)
 
 # Make diagram for minerals only 20201007
 species(icr)
-mcr <- mosaic(bases, pH = pH, O2 = O2, T = T, P = P, IS = sat$IS)
+mcr <- mosaic(bases, pH = pH, O2 = O2, T = T, P = P, IS = NaCl$IS)
 diagram(mcr$A.species, col = 2)
 label.figure("B")
 
 # Calculate *minimum* solubility among all the minerals 20201008
 # (i.e. saturation condition for the solution)
 # Use solubility() 20210303
-s <- solubility(iaq, bases = bases, pH = pH, O2 = O2, T = T, P = P, IS = sat$IS, in.terms.of = metal)
+s <- solubility(iaq, bases = bases, pH = pH, O2 = O2, T = T, P = P, IS = NaCl$IS, in.terms.of = metal)
 # Specify contour levels
 levels <- seq(-12, 9, 3)
 diagram(s, levels = levels, contour.method = "flattest")

demo/sphalerite.R

@@ -11,15 +11,15 @@ species("ZnS")
 iaq <- retrieve("Zn", c("O", "H", "Cl", "S"), "aq")
 
 # A function to make a single plot
-plotfun <- function(T = 400, P = 500, m_tot = 0.1, pHmin = 4, logppmmax = 3) {
+plotfun <- function(T = 400, P = 500, m_NaCl = 0.1, pHmin = 4, logppmmax = 3) {
   # 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)
+  NaCl <- NaCl(m_NaCl = m_NaCl, 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 = pH, "Cl-" = log10(NaCl$m_Cl), T = T, P = P, IS = NaCl$IS)
+  s <- solubility(iaq, bases = c("H2S", "HS-"), pH = pH, "Cl-" = log10(NaCl$m_Clminus), T = T, P = P, IS = NaCl$IS)
 
   # Convert log activity to log ppm
   sp <- convert(s, "logppm")
@@ -31,7 +31,7 @@ plotfun <- function(T = 400, P = 500, m_tot = 0.1, pHmin = 4, logppmmax = 3) {
   abline(v = pKw / 2, lty = 2, lwd = 2, col = "blue1")
 
   # Add legend
-  l <- lex(lNaCl(m_tot), lTP(T, P))
+  l <- lex(lNaCl(m_NaCl), lTP(T, P))
   legend("topright", legend = l, bty = "n")
 }
 
@@ -48,13 +48,13 @@ pagefun <- function() {
   T <- c(400, 400, 250, 250, 100, 100)
   # 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)
+  m_NaCl <- c(0.1, 1, 0.1, 1, 0.1, 1)
   # The plots have differing limits
   pHmin <- c(4, 4, 2, 2, 2, 2)
   logppmmax <- c(3, 3, 2, 2, 0, 0)
   # 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])
+  for(i in 1:6) plotfun(T = T[i], P = P[[i]], m_NaCl = m_NaCl[i], pHmin = pHmin[i], logppmmax = logppmmax[i])
 }
 
 # A function to make a png file with all the plots

demo/sum_S.R

@@ -17,14 +17,14 @@ P <- 1000
 # Calculate solution composition for 0.8 mol/kg NaCl
 # Based on activity of Cl- from Fig. 18 of Skirrow and Walsh (2002)
 m_NaCl = 0.8
-NaCl <- NaCl(m_tot = m_NaCl, T = T, P = P, pH = pH)
+NaCl <- NaCl(m_NaCl = m_NaCl, T = T, P = P, pH = pH)
 IS <- NaCl$IS
 
 # Setup chemical system
 # Use oxygen instead of O2 to get the gas
 basis(c("Fe", "SO4-2", "oxygen", "H+", "Cl-", "H2O"))
 basis("pH", pH)
-basis("Cl-", log10(NaCl$m_Cl))
+basis("Cl-", log10(NaCl$m_Clminus))
 # Add minerals as formed species
 species(c("pyrrhotite", "pyrite", "hematite", "magnetite"))
 
@@ -38,7 +38,7 @@ for(metal in c("Fe", "Au")) {
 
   # Plot diagram for Fe minerals and show sulfur species
   d <- diagram(m$A.species, bold = TRUE, lwd = 2)
-  diagram(m$A.bases[[1]], add = TRUE, col = 4, col.names = 4, lty = 4, dx = -2.5, dy = c(-4, 0, 0, 6))
+  diagram(m$A.bases[[1]], add = TRUE, col = 8, col.names = 8, lty = 4, dx = -2.5, dy = c(-4, 0, 0, 6), italic = TRUE)
   # Add water stability limit
   water.lines(d)
 

demo/uranyl.R

@@ -5,41 +5,41 @@ library(CHNOSZ)
 
 # Conditions
 logm_U <- log10(3.16e-5)
-m_tot <- 1  # mol NaCl / kg H2O
+m_NaCl <- 1  # mol NaCl / kg H2O
 T <- 200
 P <- "Psat"
 pH_lim <- c(2, 10)
 CS_lim <- c(-4, 1)
 res <- 500
 
 # Calculations for NaCl
-NaCl <- NaCl(m_tot = m_tot, T = T, P = P)
+NaCl <- NaCl(m_NaCl = m_NaCl, T = T, P = P)
 IS <- NaCl$IS
-logm_Na <- log10(NaCl$m_Na)
-logm_Cl <- log10(NaCl$m_Cl)
+logm_Naplus <- log10(NaCl$m_Naplus)
+logm_Clminus <- log10(NaCl$m_Clminus)
 
 # Total carbonate-pH
 iaq <- retrieve("U", ligands = c("C", "O", "H", "Cl", "Na"), state = "aq")
 icr <- retrieve("U", ligands = c("C", "O", "H", "Cl", "Na"), state = "cr")
 basis(c("UO2+2", "CO3-2", "Na+", "Cl-", "H+", "H2O", "O2"))
-basis(c("Na+", "Cl-"), c(logm_Na, logm_Cl))
+basis(c("Na+", "Cl-"), c(logm_Naplus, logm_Clminus))
 species(iaq, logm_U)
 species(icr, add = TRUE)
 bases <- c("CO3-2", "HCO3-", "CO2")
 m <- mosaic(bases, pH = c(pH_lim, res), "CO3-2" = c(CS_lim, res), T = T, P = P, IS = IS)
 diagram(m$A.species)
-diagram(m$A.bases, add = TRUE, col = 2, lty = 2, col.names = 2)
+diagram(m$A.bases, add = TRUE, col = 8, lty = 2, col.names = 8, italic = TRUE)
 title("Uranyl-carbonate complexation at 200 \u00b0C, after Migdisov et al., 2024", font.main = 1)
 
 # Total sulfate-pH
 iaq <- retrieve("U", ligands = c("S", "O", "H", "Cl", "Na"), state = "aq")
 icr <- retrieve("U", ligands = c("S", "O", "H", "Cl", "Na"), state = "cr")
 basis(c("UO2+2", "SO4-2", "Na+", "Cl-", "H+", "H2O", "O2"))
-basis(c("Na+", "Cl-"), c(logm_Na, logm_Cl))
+basis(c("Na+", "Cl-"), c(logm_Naplus, logm_Clminus))
 species(iaq, logm_U)
 species(icr, add = TRUE)
 bases <- c("SO4-2", "HSO4-", "HS-", "H2S")
 m <- mosaic(bases, pH = c(pH_lim, res), "SO4-2" = c(CS_lim, res), T = T, P = P, IS = IS)
 diagram(m$A.species)
-diagram(m$A.bases, add = TRUE, col = 2, lty = 2, col.names = 2)
+diagram(m$A.bases, add = TRUE, col = 8, lty = 2, col.names = 8, italic = TRUE)
 title("Uranyl-sulfate complexation at 200 \u00b0C, after Migdisov et al., 2024", font.main = 1)

inst/NEWS.Rd

@@ -15,7 +15,7 @@
 \newcommand{\Cp}{\ifelse{latex}{\eqn{C_P}}{\ifelse{html}{\out{<I>C<sub>P</sub></I>}}{Cp}}}
 \newcommand{\DG0}{\ifelse{latex}{\eqn{{\Delta}G^{\circ}}}{\ifelse{html}{\out{&Delta;<I>G</I>&deg;}}{ΔG°}}}
 
-\section{Changes in CHNOSZ version 2.1.0-38 (2025-01-07)}{
+\section{Changes in CHNOSZ version 2.1.0-39 (2025-01-08)}{
 
   \subsection{OBIGT DEFAULT DATA}{
     \itemize{
@@ -110,11 +110,16 @@
       \item Merge \file{extdata/adds} and \file{extdata/cpetc} into
       \file{extdata/misc}.
 
-      \item In \code{diagram}, the default for the \strong{type} argument when
-      using the output from \code{solubility} is now \code{loga.balance} (sum
-      of activities of aqueous species) rather than \code{loga.equil}
+      \item In \code{diagram()}, the default for the \strong{type} argument
+      when using the output from \code{solubility()} is now \code{loga.balance}
+      (sum of activities of aqueous species) rather than \code{loga.equil}
       (activities of individual aqueous species).
 
+      \item In \code{NaCl()}, rename \strong{m_tot} to \strong{m_NaCl} (moles
+      of of NaCl added to 1 kg H\s{2}O) and rename output with \samp{m_Naplus},
+      \samp{m_Clminus}, \samp{m_NaCl0}, and \samp{m_HCl0} (molalities of
+      aqueous species).
+
     }
   }
 

inst/tinytest/test-mosaic.R

@@ -88,8 +88,8 @@ basis("H2S", -2)
 # Set a low logfO2 to get into H2S - HS- fields
 basis("O2", -40)
 # Calculate solution composition for 1 mol/kg NaCl
-NaCl <- NaCl(T = T, P = P, m_tot = 1)
-basis("Cl-", log10(NaCl$m_Cl))
+NaCl <- NaCl(m_NaCl = 1, T = T, P = P)
+basis("Cl-", log10(NaCl$m_Clminus))
 # Calculate affinity with changing basis species
 bases <- c("H2S", "HS-", "HSO4-", "SO4-2")
 m <- mosaic(bases, pH = c(2, 10), T = 250, P = 500, IS = NaCl$IS)