[CHNOSZ-commits] r652 - in pkg/CHNOSZ: . R demo inst man tests/testthat

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Fri Mar 19 08:50:28 CET 2021


Author: jedick
Date: 2021-03-19 08:50:28 +0100 (Fri, 19 Mar 2021)
New Revision: 652

Modified:
   pkg/CHNOSZ/DESCRIPTION
   pkg/CHNOSZ/R/solubility.R
   pkg/CHNOSZ/demo/solubility.R
   pkg/CHNOSZ/demo/zinc.R
   pkg/CHNOSZ/inst/NEWS.Rd
   pkg/CHNOSZ/man/solubility.Rd
   pkg/CHNOSZ/tests/testthat/test-solubility.R
Log:
Remove automatic detection of dissociation reactions from solubility()


Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION	2021-03-19 07:15:53 UTC (rev 651)
+++ pkg/CHNOSZ/DESCRIPTION	2021-03-19 07:50:28 UTC (rev 652)
@@ -1,6 +1,6 @@
 Date: 2021-03-19
 Package: CHNOSZ
-Version: 1.4.0-21
+Version: 1.4.0-22
 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/solubility.R
===================================================================
--- pkg/CHNOSZ/R/solubility.R	2021-03-19 07:15:53 UTC (rev 651)
+++ pkg/CHNOSZ/R/solubility.R	2021-03-19 07:50:28 UTC (rev 652)
@@ -2,7 +2,6 @@
 # 20181031 first version jmd
 # 20181106 work on the output from affinity(); no "equilibrate()" needed!
 # 20190117 add find.IS and test for dissociation reaction
-# 20190730 add codeanal; test all species for dissociation reaction
 # 20210319 use list of aqueous species as main argument (with back-compatibility for affinity output) and handle multiple minerals
 
 ## if this file is interactively sourced, the following are also needed to provide unexported functions:
@@ -16,11 +15,10 @@
 # species() should be used first to load the minerals (all bearing the same metal)
 # 'iaq' lists aqueous species that can be produced by dissolution of the minerals
 # '...' contains arguments for affinity() or mosaic() (i.e. plotting variables)
-# FIXME: what to do about 'dissociation' argument?
-solubility <- function(iaq, ..., in.terms.of = NULL, dissociation = NULL, find.IS = FALSE) {
+solubility <- function(iaq, ..., in.terms.of = NULL, dissociate = FALSE, find.IS = FALSE) {
 
   # If iaq is the output of affinity(), use old method 20210318
-  if(is.list(iaq)) return(solubility_calc(aout = iaq, in.terms.of = in.terms.of, dissociation = dissociation, find.IS = find.IS))
+  if(is.list(iaq)) return(solubility_calc(aout = iaq, in.terms.of = in.terms.of, dissociate = dissociate, find.IS = find.IS))
   # Check whether to use affinity() or mosaic()
   ddd <- list(...)
   if(identical(names(ddd)[1], "bases")) is.mosaic <- TRUE else is.mosaic <- FALSE
@@ -47,7 +45,7 @@
     species(iaq)
     if(is.mosaic) a <- suppressMessages(mosaic(...)$A.species) else a <- suppressMessages(affinity(...))
     # Calculate solubility of this mineral
-    scalc <- solubility_calc(a, in.terms.of = in.terms.of, dissociation = dissociation, find.IS = find.IS)
+    scalc <- solubility_calc(a, in.terms.of = in.terms.of, dissociate = dissociate, find.IS = find.IS)
     # Store the solubilities in the list
     slist[[i]] <- scalc$loga.balance
   }
@@ -75,7 +73,7 @@
 
 # The "nuts and bolts" of solubility calculations
 # Moved from solubility() to solubility_calc() 20210318
-solubility_calc <- function(aout, dissociation = NULL, find.IS = FALSE, in.terms.of = NULL, codeanal = FALSE) {
+solubility_calc <- function(aout, dissociate = NULL, find.IS = FALSE, in.terms.of = NULL) {
   ## concept: the logarithms of activities of species at equilibrium are equal to
   ## Astar, the affinities calculated for unit activities of species
 
@@ -84,35 +82,6 @@
   thisfun <- aout$fun
   if(thisfun=="mosaic") aout <- aout$A.species
 
-  ## does the system involve a dissociation reaction?
-  if(is.null(dissociation)) {
-    # assume FALSE unless determined otherwise
-    dissociation <- FALSE
-    # we can only test for dissociation with more than one basis species 20190730
-    nbasis <- nrow(aout$basis)
-if(codeanal) print(paste("number of basis species:", nbasis))
-    if(nbasis > 1) {
-if(codeanal) print("-- testing for dissociation reactions:")
-      # if the reaction to form the first species involves the second basis species, we consider it to be a dissociation reaction
-      #if(aout$species[1, 2] != 0) {
-      # change this to test for all species 20190730
-      if(all(aout$species[2] != 0)) {
-if(codeanal) print("   reactions to form all species involve the second basis species")
-        # 20190123 (corundum calculation): if there are only H2O, H+, and e-
-        # besides the first basis species, it's not a dissociation reaction
-        nH2O <- sum(rownames(aout$basis) %in% c("H2O", "H+", "e-", "O2", "H2"))
-        if(nbasis > (nH2O + 1)) {
-          # if we got here, and the second basis species is H2O, H+, e-, O2 (or others?), we don't have enough information, so stop
-          if(rownames(aout$basis)[2] %in% c("H2O", "H+", "e-", "O2", "H2")) {
-            stop("Unsure whether the first formation reaction is a dissociation reaction.\nSet the 'dissociation' argument to TRUE or FALSE, or redefine the basis to put a product ion second.")
-          }
-          dissociation <- TRUE
-        } else if(codeanal) print("   the basis consists of only the conserved species and H2O-derived species")
-      } else if(codeanal) print("   at least one formation reaction doesn't involve the second basis species")
-    } else if(codeanal) print("-- no dissociation reactions are possible")
-    message("solubility: test for dissociation reaction returns ", dissociation)
-  }
-
   # get starting ionic strength (probably zero, but could be anything set by user)
   IS <- aout$IS
   if(is.null(IS)) IS <- 0
@@ -121,7 +90,6 @@
   bout <- balance(aout)
   n.balance <- bout$n.balance
   balance <- bout$balance
-if(codeanal) print(paste0("balancing coefficients [", balance, "]: ", paste(n.balance, collapse = " ")))
   # get logarithm of total activity of the conserved basis species
   logabfun <- function(loga.equil, n.balance) {
     # exponentiate, multiply by n.balance, sum, logarithm
@@ -139,26 +107,18 @@
     Astar <- function(i) aout$values[[i]] + aout$species$logact[i]
     loga.equil <- lapply(1:length(aout$values), Astar)
 
-    ## for a dissociation on a *per reaction* (not system) basis,
-    ## apply the divisor here and skip the if(dissociation){} part below
+    ## for a dissociation on a *per reaction* (not system) basis, apply the divisor here
     ## (can be used to reproduce Fig. 4 of Manning et al., 2013)
-    if(is.numeric(dissociation)) {
-if(codeanal) for(ii in 1:length(loga.equil)) print(paste0("loga.equil0 [", aout$species$name[ii], "]: ", round(loga.equil[[ii]], 3)))
-if(codeanal) print(paste0("applying ", dissociation, "-fold dissociation correction (no interaction)"))
-      loga.equil <- lapply(loga.equil, "/", dissociation)
+    if(is.numeric(dissociate)) {
+      loga.equil <- lapply(loga.equil, "/", dissociate)
     }
 
-
     # get the total activity of the balancing basis species
     loga.balance <- logabfun(loga.equil, n.balance)
-if(codeanal & !isTRUE(dissociation)) print(paste0("loga.balance [", balance, "]: ", round(loga.balance, 3)))
 
     # recalculate things for a dissociation reaction (like CaCO3 = Ca+2 + CO3+2)
-    if(isTRUE(dissociation)) {
+    if(isTRUE(dissociate)) {
       ndissoc <- 2
-if(codeanal) print(paste0("loga.balance0 [", balance, "]: ", round(loga.balance, 3)))
-if(codeanal) for(ii in 1:length(loga.equil)) print(paste0("loga.equil0 [", aout$species$name[ii], "]: ", round(loga.equil[[ii]], 3)))
-if(codeanal) print(paste0("applying ", ndissoc, "-fold dissociation correction (with interaction)"))
       # the number of dissociated products is the exponent in the activity product
       loga.dissoc <- loga.balance / ndissoc
       # the contribution to affinity
@@ -167,7 +127,6 @@
       aout$values <- mapply("-", aout$values, Adissoc, SIMPLIFY = FALSE)
       loga.equil <- lapply(1:length(aout$values), Astar)
       loga.balance <- logabfun(loga.equil, n.balance)
-if(codeanal) print(paste0("loga.balance [", balance, "]: ", round(loga.balance, 3)))
       # check that the new loga.balance == loga.dissoc
       # TODO: does this work for non-1:1 species?
       stopifnot(all.equal(loga.balance, loga.dissoc))
@@ -179,7 +138,7 @@
     sum.mZ2 <- 0
     ion.names <- character()
     # is there an ion in the basis species?
-    if(dissociation) {
+    if(dissociate) {
       basis.ion <- rownames(aout$basis)[2]
       Z.basis.ion <- makeup(basis.ion)["Z"]
       if(!is.na(Z.basis.ion)) {
@@ -234,7 +193,6 @@
     ibalance <- which.balance(aout$species)
     coeff <- sbasis[, ibalance][1]
     loga.balance <- loga.balance - log10(coeff)
-if(codeanal) print(paste0("loga.balance [in terms of ", in.terms.of, "]: ", round(loga.balance, 3)))
   }
 
   # make the output (we don't deal with normalized formulas yet, so for now m.balance==n.balance)

Modified: pkg/CHNOSZ/demo/solubility.R
===================================================================
--- pkg/CHNOSZ/demo/solubility.R	2021-03-19 07:15:53 UTC (rev 651)
+++ pkg/CHNOSZ/demo/solubility.R	2021-03-19 07:50:28 UTC (rev 652)
@@ -14,7 +14,7 @@
 opar <- par(no.readonly = TRUE)
 layout(matrix(1:4, nrow = 2))
 
-# set pH and T range and resolution, constant temperature and ionic strength
+# Set pH and T range and resolution, constant temperature and ionic strength
 pH <- c(0, 14)
 T <- c(0, 300)
 res <- 100
@@ -21,18 +21,18 @@
 T1 <- 25
 IS <- 0
 
-# start with CO2
+# Start with CO2
 basis(c("carbon dioxide", "H2O", "O2", "H+"))
-# ca. atmospheric PCO2
+# This is ca. atmospheric PCO2
 basis("CO2", -3.5)
 species(c("CO2", "HCO3-", "CO3-2"))
 a <- affinity(pH = c(pH, res), T = T1, IS = IS)
 s <- solubility(a)
-# first plot total activity line
+# First plot total activity line
 diagram(s, ylim = c(-10, 4), type = "loga.balance", lwd = 4, col = "green2")
-# add activities of species
+# Add activities of species
 diagram(s, ylim=c(-10, 4), add = TRUE, dy = 1)
-# add legend
+# Add legend
 lexpr <- as.expression(c("total", expr.species("CO2", state = "aq"),
   expr.species("HCO3-"), expr.species("CO3-2")))
 legend("topleft", lty = c(1, 1:3), lwd = c(4, 2, 2, 2),
@@ -40,7 +40,7 @@
 title(main = substitute("Solubility of"~what~"at"~T~degree*"C",
   list(what = expr.species("CO2"), T = T1)), line = 1.6)
 mtext("cf. Fig. 4.5 of Stumm and Morgan, 1996")
-# check the endpoints
+# Check the endpoints
 stopifnot(round(s$loga.balance[c(1, res)])==c(-5, 6))
 
 # CO2 T-pH plot
@@ -49,11 +49,12 @@
 diagram(s, type = "loga.balance")
 title(main = substitute("Solubility of"~what, list(what = expr.species("CO2"))))
 
-# now do calcite
+# Now do calcite
 basis(c("calcite", "Ca+2", "H2O", "O2", "H+"))
 species(c("CO2", "HCO3-", "CO3-2"))
 a <- affinity(pH = c(pH, res), T = T1, IS = IS)
-s <- solubility(a)
+# Change this to dissociate = 2 to reproduce straight lines in Fig. 4A of Manning et al., 2013
+s <- solubility(a, dissociate = TRUE)
 diagram(s, ylim = c(-10, 4), type = "loga.balance", lwd = 4, col = "green2")
 diagram(s, add = TRUE, dy = 1)
 legend("topright", lty = c(1, 1:3), lwd = c(4, 2, 2, 2),
@@ -61,10 +62,10 @@
 title(main = substitute("Solubility of"~what~"at"~T~degree*"C",
   list(what = "calcite", T = T1)), line = 1.6)
 mtext("cf. Fig. 4A of Manning et al., 2013")
-# check the endpoints
+# Check the endpoints
 stopifnot(round(s$loga.balance[c(1, res)])==c(4, -4))
 
-# calcite T-pH plot
+# Calcite T-pH plot
 a <- affinity(pH = c(pH, res), T = c(T, res), IS = IS)
 s <- solubility(a)
 diagram(s, type = "loga.balance")

Modified: pkg/CHNOSZ/demo/zinc.R
===================================================================
--- pkg/CHNOSZ/demo/zinc.R	2021-03-19 07:15:53 UTC (rev 651)
+++ pkg/CHNOSZ/demo/zinc.R	2021-03-19 07:50:28 UTC (rev 652)
@@ -69,8 +69,7 @@
 # Calculate *minimum* solubility among all the minerals 20201008
 # (i.e. saturation condition for the solution)
 # Use solubilities() 20210303
-# FIXME: why do we need to set dissociation = FALSE here?
-s <- solubility(iaq, bases = bases, pH = pH, O2 = O2, T = T, P = P, IS = sat$IS, in.terms.of = "Zn", dissociation = FALSE)
+s <- solubility(iaq, bases = bases, pH = pH, O2 = O2, T = T, P = P, IS = sat$IS, in.terms.of = "Zn")
 # Specify contour levels
 levels <- seq(-12, 9, 3)
 diagram(s, type = "loga.balance", levels = levels, contour.method = "flattest")

Modified: pkg/CHNOSZ/inst/NEWS.Rd
===================================================================
--- pkg/CHNOSZ/inst/NEWS.Rd	2021-03-19 07:15:53 UTC (rev 651)
+++ pkg/CHNOSZ/inst/NEWS.Rd	2021-03-19 07:50:28 UTC (rev 652)
@@ -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.4.0-21 (2021-03-19)}{
+\section{Changes in CHNOSZ version 1.4.0-22 (2021-03-19)}{
 
   \subsection{NEW FEATURES}{
     \itemize{
@@ -32,6 +32,10 @@
   \subsection{OTHER CHANGES}{
     \itemize{
 
+      \item Automatic detection of dissociation reactions was fragile and has
+      been removed from \code{solubility()}. The new default (\code{dissociate
+        = FALSE}) is to not condiser dissociation reactions.
+
       \item In the return value of \code{diagram()}, \code{predominant.values}
       previously only contained affinities (extracted from the \code{values}
       element of the \code{eout} argument). Now it has the values for the

Modified: pkg/CHNOSZ/man/solubility.Rd
===================================================================
--- pkg/CHNOSZ/man/solubility.Rd	2021-03-19 07:15:53 UTC (rev 651)
+++ pkg/CHNOSZ/man/solubility.Rd	2021-03-19 07:50:28 UTC (rev 652)
@@ -7,7 +7,7 @@
 }
 
 \usage{
-  solubility(iaq, ..., in.terms.of = NULL, dissociation = NULL, find.IS = FALSE)
+  solubility(iaq, ..., in.terms.of = NULL, dissociate = FALSE, find.IS = FALSE)
 }
 
 \arguments{
@@ -14,7 +14,7 @@
   \item{iaq}{character (names) or numeric (species indices) of aqueous species}
   \item{...}{arguments for \code{\link{affinity}} or \code{\link{mosaic}} (i.e. plotting variables)}
   \item{in.terms.of}{character, express the total solubility in terms of moles of this species}
-  \item{dissociation}{logical, does the mineral undergo a dissociation reaction?}
+  \item{dissociate}{logical, does the mineral undergo a dissociation reaction?}
   \item{find.IS}{logical, find the equilibrium ionic strength by iteration?}
 }
 
@@ -26,7 +26,6 @@
 It must be possible to obtain a valid set of basis species by substituting each of the minerals or gases in the first position of the current \code{\link{basis}} defintion, and all of the aqueous species must include that basis species in their formation reactions.
 (This essentially means that all minerals, gases and aqueous species must share a common element, which is what the reactions are balanced on.)
 
-The output of the function is structured to facilitate making useful diagrams.
 For a single mineral or gas, the output of \code{solubility} can be used by \code{\link{diagram}} with \code{type = NULL} (the default) to plot the activities of the aqueous species or with \code{type = "loga.balance"} to plot the sum of activities of aqueous species, which corresponds to the solubility of the mineral or gas.
 This value corresponds to the total extent of dissolution of the mineral or gas; \code{in.terms.of} can be used to express this value in terms of another species or element.
 For example, for dissolution of gaseous S\s{2}, \code{in.terms.of = "S"} gives the total amount of S in solution, which is twice the amount of S\s{2} dissolved.
@@ -35,7 +34,7 @@
 For multiple minerals, the function calculates the solubilities for each of the minerals separately; these are stored in the \code{loga.equil} element of the output.
 The overall \emph{minimum} solubility among all the minerals at each point is stored in \code{loga.balance}.
 This corresponds to the total activity of dissolved species in equilibrium with the most stable mineral.
-In contrast to the situation for a single mineral or gas, \code{\link{diagram}} be default plots \code{loga.balance}, and \code{type = "loga.balance"} should be used to plot the solubilities for the individual minerals or gases.
+In contrast to the situation for a single mineral or gas, \code{\link{diagram}} by default plots \code{loga.balance}; \code{type = "loga.equil"} can be used to plot the solubilities for the individual minerals or gases.
 }
 
 \section{Backward Compatibility}{
@@ -58,15 +57,9 @@
 The solution is accomplished by recalculating the affinities, essentially working backward from the assumption that the dissociation didn't occur.
 The resulting activities correspond to equilibrium considering the system-wide activity of Ca\S{+2}.
 
-The function attempts to automatically detect whether dissociation reactions are involved.
-If the formation reactions of all species includes \emph{both} the first and second basis species, the \code{dissociation} flag is set to TRUE.
-An example reaction of this type can be found in demo(solubility): CaCO3 (first basis species) = Ca+2 (second basis species) + CO3-2 (first species).
-Note that if the conserved basis species is not the first basis species, then the automatic detection of \code{dissociation} will always return FALSE.
-Therefore, a reaction corresponding to Au (fourth basis species) + ... = ... gives \code{dissociation} = FALSE (see \code{demo(gold)}).
-This algorithm for determining whether dissociation occurs is not perfect, so \code{dissociation} can be explicitly set in the arguments.
-A \emph{not recommended} alternative is to set \code{dissociation} to a numeric value corresponding to the stoichiometry of released species (i.e. 2 for a 1:1 electrolyte).
+A \emph{not recommended} alternative is to set \code{dissociate} to a numeric value corresponding to the number of dissociated species (i.e. 2 for a 1:1 electrolyte).
 This setting indicates to calculate activities on a per-reaction basis, where each reaction has its own (independent) activity of Ca\S{+2}.
-That does not give a complete equilibrium in the system, but may be required to reproduce some published diagrams.
+That does not give a complete equilibrium in the system, but may be required to reproduce some published diagrams (see comment in the calcite example of \code{demo("solubility")}).
 }
 
 \section{Ionic Strength}{
@@ -85,11 +78,10 @@
 \seealso{
 \code{\link{retrieve}} provides a way to list all of the aqueous species in the database that have the specified elements.
 
-\code{demo("solubility")} shows \T-pH diagrams related to the CO\s{2} and calcite example here.
+\code{demo("solubility")} shows solubilities of CO\s{2} and calcite calculated as a function of pH and T; note that for calcite, the \code{dissociate} argument is set to TRUE.
 \code{demo("gold")} shows solubility calculations for Au in aqueous solutions with hydroxide, chloride, and hydrosulfide complexes.
-In the latter demo, the conserved basis species is not the first one (this might have been done to prevent the incorrect automatic detection of dissociation reactions).
 
-Solubility calculations for multiple minerals are used for generating isosolubility (aka. equisolubility) lines in \code{demo("Pourbaix")} and \code{demo{"zinc"}}.
+Solubility calculations for multiple minerals are used for generating isosolubility (aka. equisolubility) lines in \code{demo("Pourbaix")} and \code{demo("zinc")}.
 The latter demo combines the calculation of solubilities with a \code{\link{mosaic}} calculation to account for the speciation of aqueous sulfur.
 
 Whereas \code{solubility} yields a stable equilibrium condition (the affinities of formation reactions of aqueous species are zero), \code{\link{equilibrate}} generates metastable equilibrium (the affinities of formation reactions are equal to each other, but not necessarily zero).
@@ -96,45 +88,6 @@
 }
 
 \examples{\dontshow{reset()}
-## Solubility of CO2 and calcite as a function of pH
-opar <- par(mfrow = c(1, 2))
-## Set pH range and resolution, constant temperature and ionic strength
-pH <- c(0, 14)
-res <- 100
-T <- 25
-IS <- 0
-## Start with CO2 (not a dissociation reaction)
-basis(c("carbon dioxide", "H2O", "O2", "H+"))
-# This is close to atmospheric PCO2
-basis("CO2", -3.5)
-species(c("CO2", "HCO3-", "CO3-2"))
-a <- affinity(pH = c(pH, res), T = T, IS = IS)
-s <- solubility(a)
-# First plot total activity line
-diagram(s, ylim = c(-10, 4), type = "loga.balance", lwd = 4, col = "green2")
-# Add activities of species
-diagram(s, ylim=c(-10, 4), add = TRUE, dy = 1)
-# Add legend
-lexpr <- as.expression(c("total", expr.species("CO2", state = "aq"),
-  expr.species("HCO3-"), expr.species("CO3-2")))
-legend("topleft", lty = c(1, 1:3), lwd = c(4, 2, 2, 2),
-  col = c("green2", rep("black", 3)), legend = lexpr)
-title(main = substitute("Solubility of"~what~"at"~T~degree*"C",
-  list(what = expr.species("CO2"), T = T)), line = 1.5)
-mtext("cf. Fig. 4.5 of Stumm and Morgan, 1996")
-## Now do calcite (a dissociation reaction)
-basis(c("calcite", "Ca+2", "H2O", "O2", "H+"))
-species(c("CO2", "HCO3-", "CO3-2"))
-a <- affinity(pH = c(pH, res), T = T, IS = IS)
-s <- solubility(a)
-diagram(s, ylim = c(-10, 4), type = "loga.balance", lwd = 4, col = "green2")
-diagram(s, add = TRUE, dy = 1)
-legend("topright", lty = c(1, 1:3), lwd = c(4, 2, 2, 2),
-  col = c("green2", rep("black", 3)), legend = lexpr)
-title(main = substitute("Solubility of"~what~"at"~T~degree*"C",
-  list(what = "calcite", T = T)))
-mtext("cf. Fig. 4A of Manning et al., 2013")
-par(opar)
 
 ## Two ways to calculate pH-dependent solubility of calcite
 ## with ionic strength determination
@@ -143,14 +96,14 @@
 species(c("CO2", "HCO3-", "CO3-2"))
 # Ionic strength calculations don't converge below around pH = 3
 a <- affinity(pH = c(4, 14))
-sa0 <- solubility(a)
-saI <- solubility(a, find.IS = TRUE)
+sa0 <- solubility(a, dissociate = TRUE)
+saI <- solubility(a, dissociate = TRUE, find.IS = TRUE)
 ## Method 2: CO2 and carbonate species as basis species
 basis(c("calcite", "CO2", "H2O", "O2", "H+"))
 species(c("Ca+2"))
 m <- mosaic(c("CO2", "HCO3-", "CO3-2"), pH = c(4, 14))
-sm0 <- solubility(m)
-smI <- solubility(m, find.IS = TRUE)
+sm0 <- solubility(m, dissociate = TRUE)
+smI <- solubility(m, dissociate = TRUE, find.IS = TRUE)
 ## Plot the results
 plot(0, 0, xlab="pH", ylab="solubility, log mol", xlim = c(4, 14), ylim = c(-5, 2))
 # Method 1 with/without ionic strength
@@ -169,10 +122,4 @@
 ## more flexible; e.g. complexes with Ca+2 could be included
 }
 
-\references{
-Manning, C. E., Shock, E. L. and Sverjensky, D. A. (2013) The chemistry of carbon in aqueous fluids at crustal and upper-mantle conditions: Experimental and theoretical constraints. \emph{Rev. Mineral. Geochem.} \bold{75}, 109--148. \url{https://doi.org/10.2138/rmg.2013.75.5}
-
-Stumm, W. and Morgan, J. J. (1996) \emph{Aquatic Chemistry: Chemical Equilibria and Rates in Natural Waters}, John Wiley & Sons, New York, 1040 p. \url{https://www.worldcat.org/oclc/31754493}
-}
-
 \concept{Main workflow}

Modified: pkg/CHNOSZ/tests/testthat/test-solubility.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-solubility.R	2021-03-19 07:15:53 UTC (rev 651)
+++ pkg/CHNOSZ/tests/testthat/test-solubility.R	2021-03-19 07:50:28 UTC (rev 652)
@@ -30,7 +30,7 @@
   basis(c("calcite", "Ca+2", "H2O", "O2", "H+"))
   species(c("CO2", "HCO3-", "CO3-2"))
   a <- affinity(pH = c(pH, res), T = T, IS = IS)
-  s <- solubility(a)
+  s <- solubility(a, dissociate = TRUE)
   # here we need to also set the activity of Ca+2
   checkfun <- function(i) {
     logact <- sapply(s$loga.equil, "[", i)
@@ -43,16 +43,6 @@
   expect_equal(max(abs(unlist(checkfun(99)$values))), 0)
 })
 
-test_that("solubility() catches some error conditions", {
-  # demo(solubility) has basis(c("calcite", "Ca+2", "H2O", "O2", "H+")), but what if the user puts H2O second?
-  basis(c("calcite", "H2O", "Ca+2", "O2", "H+"))
-  species(c("CO2", "HCO3-", "CO3-2"))
-  a <- affinity()
-  expect_message(solubility(a), "test for dissociation reaction returns FALSE")
-  # TODO: write a test to get this error 20190730
-  #expect_error(solubility(a), "Unsure whether the first formation reaction is a dissociation reaction.")
-})
-
 test_that("backward compatible and new calling styles produce identical results", {
   # Test added 20210319
   # Calculate solubility of a single substance:



More information about the CHNOSZ-commits mailing list