[CHNOSZ-commits] r754 - in pkg/CHNOSZ: . R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Nov 25 10:44:10 CET 2022
Author: jedick
Date: 2022-11-25 10:44:09 +0100 (Fri, 25 Nov 2022)
New Revision: 754
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/NaCl.R
pkg/CHNOSZ/inst/NEWS.Rd
pkg/CHNOSZ/man/extdata.Rd
Log:
NaCl() now works with 'm_tot = 0'
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2022-11-18 00:40:29 UTC (rev 753)
+++ pkg/CHNOSZ/DESCRIPTION 2022-11-25 09:44:09 UTC (rev 754)
@@ -1,6 +1,6 @@
-Date: 2022-11-18
+Date: 2022-11-22
Package: CHNOSZ
-Version: 1.9.9-45
+Version: 1.9.9-46
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-11-18 00:40:29 UTC (rev 753)
+++ pkg/CHNOSZ/R/NaCl.R 2022-11-25 09:44:09 UTC (rev 754)
@@ -1,61 +1,67 @@
# NaCl.R
-# calculate ion activities and ionic strength
-# given a total molality of NaCl
+# Calculate ion activities and ionic strength
+# for a total molality of NaCl,
# taking account of ion association: Na+ + Cl- = NaCl(aq)
# 20181102 jmd first version
# 20181105 use activity coefficient of Na+
# 20181106 use activity coefficient of NaCl
+# 20221122 make it work with m_tot = 0
-NaCl <- function(T=seq(100, 500, 100), P=1000, m_tot=2, ...) {
- # define a function for the reaction quotient
+NaCl <- function(T = seq(100, 500, 100), P = 1000, m_tot = 2, ...) {
+ # Define a function for the reaction quotient
logQ <- function(m_Cl, gam_NaCl, gam_Na, gam_Cl) {
- # starting with Q = a_NaCl / (a_Na+ * a_Cl-),
+ # Starting with Q = a_NaCl / (a_Na+ * a_Cl-),
# substitute m_tot = m_NaCl + m_Cl and m_Cl = m_Na
# to write:
log10( (m_tot - m_Cl) * gam_NaCl / (m_Cl ^ 2 * gam_Na * gam_Cl) )
}
- # define a function for affinity = log(K / Q)
+ # Define a function for affinity = log(K / Q)
A <- function(m_Cl, gam_NaCl, gam_Na, gam_Cl, logK) logK - logQ(m_Cl, gam_NaCl, gam_Na, gam_Cl)
- # calculate equilibrium constant at all temperatures (standard conditions: IS = 0)
+ # Calculate equilibrium constant at all temperatures (standard conditions: IS = 0)
logK <- subcrt(c("Na+", "Cl-", "NaCl"), c(-1, -1, 1), T = T, P = P, ...)$out$logK
- # calculate Debye-Huckel parameters at all temperatures
+ # Calculate Debye-Huckel parameters at all temperatures
wout <- water(c("A_DH", "B_DH"), T = convert(T, "K"), P = P)
- # initialize output variables
+ # Initialize output variables
N <- length(T)
ISout <- a_Cl <- numeric(N)
- # initial guess for m_Cl and ionic strength assuming complete dissociation of NaCl
+ # Initial guess for m_Cl and ionic strength assuming complete dissociation of NaCl
IS <- m_Cl <- rep(m_tot, N)
- # the corresponding total molality of dissolved species (NaCl + Cl- + Na+)
+ # The corresponding total molality of dissolved species (NaCl + Cl- + Na+)
m_star <- (m_tot - m_Cl) + 2*m_Cl
- # the species indices for Na+, Cl-, and NaCl(aq)
+ # The species indices for Na+, Cl-, and NaCl(aq)
ispecies <- info(c("Na+", "Cl-", "NaCl"))
- # we start by doing calculations for all temperatures
- doit <- !logical(N)
- while(any(doit)) {
- # calculate activity coefficient at given ionic strength
- speciesprops <- rep(list(data.frame(G=numeric(N))), length(ispecies))
- gammas <- suppressMessages(nonideal(ispecies, speciesprops, IS=IS, T=convert(T, "K"), P=P, A_DH=wout$A_DH, B_DH=wout$B_DH, m_star=m_star))
- gam_Na <- 10^gammas[[1]]$loggam
- gam_Cl <- 10^gammas[[2]]$loggam
- gam_NaCl <- 10^gammas[[3]]$loggam
- # in case Setchenow calculations are turned off 20191209
- if(length(gam_NaCl) == 0) gam_NaCl <- rep(1, length(T))
- # solve for m_Cl
- for(i in which(doit)) m_Cl[i] <- uniroot(A, c(0, m_tot), gam_NaCl=gam_NaCl[i], gam_Na=gam_Na[i], gam_Cl=gam_Cl[i], logK=logK[i])$root
- # calculate new total molality
- m_star <- (m_tot - m_Cl) + 2*m_Cl
- # calculate new ionic strength and deviation
- ISnew <- m_Cl
- dIS <- ISnew - IS
- # set new ionic strength
- IS <- ISnew
- # keep going until the deviation in ionic strength at any temperature is less than 0.01
- doit <- abs(dIS) > 0.01
+ # Data frame needed for nonideal()
+ speciesprops <- rep(list(data.frame(G = numeric(N))), length(ispecies))
+ # Only try to find root for non-zero NaCl concentration 20221122
+ if(m_tot > 0) {
+ # Start by doing calculations for all temperatures
+ doit <- !logical(N)
+ while(any(doit)) {
+ # Calculate activity coefficient at given ionic strength
+ gammas <- suppressMessages(nonideal(ispecies, speciesprops, IS = IS, T = convert(T, "K"), P = P, A_DH = wout$A_DH, B_DH = wout$B_DH, m_star = m_star))
+ gam_Na <- 10^gammas[[1]]$loggam
+ gam_Cl <- 10^gammas[[2]]$loggam
+ gam_NaCl <- 10^gammas[[3]]$loggam
+ # In case Setchenow calculations are turned off 20191209
+ if(length(gam_NaCl) == 0) gam_NaCl <- rep(1, length(T))
+ # Solve for m_Cl
+ for(i in which(doit)) m_Cl[i] <- uniroot(A, c(0, m_tot), gam_NaCl = gam_NaCl[i], gam_Na = gam_Na[i], gam_Cl = gam_Cl[i], logK = logK[i])$root
+ # Calculate new total molality
+ m_star <- (m_tot - m_Cl) + 2*m_Cl
+ # Calculate new ionic strength and deviation
+ ISnew <- m_Cl
+ dIS <- ISnew - IS
+ # Set new ionic strength
+ IS <- ISnew
+ # Keep going until the deviation in ionic strength at any temperature is less than 0.01
+ doit <- abs(dIS) > 0.01
+ }
}
- # assemble final molality of Cl- and gammas
- gammas <- suppressMessages(nonideal(ispecies, speciesprops, IS=IS, T=convert(T, "K"), P=P, A_DH=wout$A_DH, B_DH=wout$B_DH))
+ # Assemble final molality of Cl- and gammas
+ gammas <- suppressMessages(nonideal(ispecies, speciesprops, IS = IS, T = convert(T, "K"), P = P, A_DH = wout$A_DH, B_DH = wout$B_DH))
gam_Na <- 10^gammas[[1]]$loggam
gam_Cl <- 10^gammas[[2]]$loggam
- # return the calculated values
- list(IS=IS, m_Cl=m_Cl, gam_Na=gam_Na, gam_Cl=gam_Cl, gam_NaCl=gam_NaCl)
+ gam_NaCl <- 10^gammas[[3]]$loggam
+ # Return the calculated values
+ list(IS = IS, m_Cl = m_Cl, gam_Na = gam_Na, gam_Cl = gam_Cl, gam_NaCl = gam_NaCl)
}
Modified: pkg/CHNOSZ/inst/NEWS.Rd
===================================================================
--- pkg/CHNOSZ/inst/NEWS.Rd 2022-11-18 00:40:29 UTC (rev 753)
+++ pkg/CHNOSZ/inst/NEWS.Rd 2022-11-25 09:44:09 UTC (rev 754)
@@ -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-45 (2022-11-18)}{
+\section{Changes in CHNOSZ version 1.9.9-46 (2022-11-22)}{
\subsection{MAJOR USER-VISIBLE CHANGES}{
\itemize{
@@ -99,6 +99,9 @@
aqueous basis species. This obviates a workaround that was previously
used in the Mosaic Stacking 2 section of \strong{multi-metal.Rmd}.
+ \item Add \code{add_alpha()} to add alpha value to a color specified in
+ any notation accepted by \code{col2rgb()}.
+
}
}
@@ -195,9 +198,11 @@
\item \code{axis.label()} produces labels with units delimited by
parentheses instead of a comma.
- \item \code{seq2aa} now has the \strong{sequence} argument first and a
+ \item \code{seq2aa()} now has the \strong{sequence} argument first and a
default of NA for \strong{protein} (the protein name).
+ \item \code{NaCl()} now works with \code{mtot = 0}.
+
}
}
Modified: pkg/CHNOSZ/man/extdata.Rd
===================================================================
--- pkg/CHNOSZ/man/extdata.Rd 2022-11-18 00:40:29 UTC (rev 753)
+++ pkg/CHNOSZ/man/extdata.Rd 2022-11-25 09:44:09 UTC (rev 754)
@@ -17,7 +17,7 @@
\item Other files with names like \code{xxx_yyyy.csv} contain thermodynamic data from other sources; xxx in the filename corresponds to the reference in \code{\link{thermo}$OBIGT} and yyyy gives the year of publication.
\code{\link{Berman}} uses these data for the calculation of thermodynamic properties at specified \P and \T, which are then available for use in \code{\link{subcrt}}.
If there are any duplicated mineral names in the files, only the most recent data are used, as determined by the year in the file name.
- Following conventions used in other data files, the names of sanidine and microcline were changed to K-feldspar,high and K-feldspar,low.
+ Following conventions used SUPCRT92 (see Helgeson et al., 1978), the names of sanidine and microcline were changed to K-feldspar,high and K-feldspar,low (by using the same names in all data files, loading the optional SUPCRT92 data file updates these minerals rather than makes new ones).
\item \code{sympy.R} is an R script that uses \CRANpkg{rSymPy} to symbolically integrate Bermans's equations for heat capacity and volume to write experessions for enthalpy, entropy and Gibbs energy.
\item The \code{testing} directory contains data files based on Berman and Aranovich (1996). These are used to demonstrate the addition of data from a user-supplied file (see \code{\link{Berman}}).
}
@@ -83,6 +83,8 @@
Gaucher, E. A., Thomson, J. M., Burgan, M. F. and Benner, S. A (2003) Inferring the palaeoenvironment of ancient bacteria on the basis of resurrected proteins. \emph{Nature} \bold{425}(6955), 285--288. \doi{10.1038/nature01977}
+Helgeson, H. C., Delany, J. M., Nesbitt, H. W. and Bird, D. K. (1978) Summary and critique of the thermodynamic properties of rock-forming minerals. \emph{Am. J. Sci.} \bold{278-A}, 1--229. \url{https://www.worldcat.org/oclc/13594862}
+
Hnědkovský, L., Wood, R. H. and Majer, V. (1996) Volumes of aqueous solutions of \CH4, \CO2, \H2S, and \NH3 at temperatures from 298.15 K to 705 K and pressures to 35 MPa. \emph{J. Chem. Thermodyn.} \bold{28}, 125--142. \doi{10.1006/jcht.1996.0011}
Hnědkovský, L. and Wood, R. H. (1997) Apparent molar heat capacities of aqueous solutions of \CH4, \CO2, \H2S, and \NH3 at temperatures from 304 K to 704 K at a pressure of 28 MPa. \emph{J. Chem. Thermodyn.} \bold{29}, 731--747. \doi{10.1006/jcht.1997.0192}
More information about the CHNOSZ-commits
mailing list