[CHNOSZ-commits] r737 - in pkg/CHNOSZ: . R inst man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Jul 23 01:52:24 CEST 2022
Author: jedick
Date: 2022-07-23 01:52:24 +0200 (Sat, 23 Jul 2022)
New Revision: 737
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/mosaic.R
pkg/CHNOSZ/R/stack_mosaic.R
pkg/CHNOSZ/inst/NEWS.Rd
pkg/CHNOSZ/man/mosaic.Rd
pkg/CHNOSZ/man/stack_mosaic.Rd
pkg/CHNOSZ/vignettes/multi-metal.Rmd
Log:
Add loga_aq argument to mosaic() and stack_mosaic()
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2022-07-22 06:33:37 UTC (rev 736)
+++ pkg/CHNOSZ/DESCRIPTION 2022-07-22 23:52:24 UTC (rev 737)
@@ -1,4 +1,4 @@
-Date: 2022-07-22
+Date: 2022-07-23
Package: CHNOSZ
Version: 1.9.9-29
Title: Thermodynamic Calculations and Diagrams for Geochemistry
Modified: pkg/CHNOSZ/R/mosaic.R
===================================================================
--- pkg/CHNOSZ/R/mosaic.R 2022-07-22 06:33:37 UTC (rev 736)
+++ pkg/CHNOSZ/R/mosaic.R 2022-07-22 23:52:24 UTC (rev 737)
@@ -5,15 +5,15 @@
# and improve speed by pre-calculating subcrt values (sout)
# 20190505 bug fix: adjust affinities of species formation reactions for mole fractions of basis species
-## if this file is interactively sourced, the following are also needed to provide unexported functions:
+## If this file is interactively sourced, the following are also needed to provide unexported functions:
#source("basis.R")
#source("util.character.R")
#source("util.args.R")
-# function to calculate affinities with mosaic of basis species
-mosaic <- function(bases, blend = TRUE, stable = list(), ...) {
+# Function to calculate affinities with mosaic of basis species
+mosaic <- function(bases, blend = TRUE, stable = list(), loga_aq = NULL, ...) {
- # argument recall 20190120
+ # Argument recall 20190120
# if the first argument is the result from a previous mosaic() calculation,
# just update the remaining arguments
if(is.list(bases)) {
@@ -31,7 +31,7 @@
}
}
- # backward compatibility 20190131:
+ # Backward compatibility 20190131:
# bases can be a vector instead of a list
if(!is.list(bases)) {
bases <- list(bases)
@@ -38,21 +38,21 @@
otherargs <- list(...)
allargs <- c(list(bases = bases, blend = blend, stable = stable), otherargs)
out <- do.call(mosaic, allargs)
- # replace A.bases (affinity calculations for all groups of basis species) with backwards-compatbile A.bases
+ # Replace A.bases (affinity calculations for all groups of basis species) with backwards-compatbile A.bases
out$A.bases <- out$A.bases[[1]]
return(out)
}
- # save starting basis and species definition
+ # Save starting basis and species definition
basis0 <- get("thermo", CHNOSZ)$basis
species0 <- get("thermo", CHNOSZ)$species
- # get species indices of requested basis species
+ # Get species indices of requested basis species
ispecies <- lapply(bases, info)
if(any(is.na(unlist(ispecies)))) stop("one or more of the requested basis species is unavailable")
- # identify starting basis species
+ # Identify starting basis species
ispecies0 <- sapply(ispecies, "[", 1)
ibasis0 <- match(ispecies0, basis0$ispecies)
- # quit if starting basis species are not present
+ # Quit if starting basis species are not present
ina <- is.na(ibasis0)
if(any(ina)) {
names0 <- unlist(lapply(bases, "[", 1))
@@ -66,7 +66,7 @@
sout <- ddd$sout
} else {
ddd_has_sout <- FALSE
- # run subcrt() calculations for all basis species and formed species 20190131
+ # Run subcrt() calculations for all basis species and formed species 20190131
# this avoids repeating the calculations in different calls to affinity()
# add all the basis species here - the formed species are already present
lapply(bases, species, add = TRUE)
@@ -73,7 +73,7 @@
sout <- affinity(..., return.sout = TRUE)
}
- # calculate affinities of the basis species themselves
+ # Calculate affinities of the basis species themselves
A.bases <- list()
for(i in 1:length(bases)) {
message("mosaic: calculating affinities of basis species group ", i, ": ", paste(bases[[i]], collapse=" "))
@@ -80,67 +80,78 @@
mysp <- species(bases[[i]])
# 20191111 include only aq species in total activity
iaq <- mysp$state == "aq"
- # use as.numeric in case a buffer is active 20201014
+ # Use as.numeric in case a buffer is active 20201014
if(any(iaq)) species(which(iaq), as.numeric(basis0$logact[ibasis0[i]]))
if(ddd_has_sout) A.bases[[i]] <- suppressMessages(affinity(...))
else A.bases[[i]] <- suppressMessages(affinity(..., sout = sout))
}
- # get all combinations of basis species
+ # Get all combinations of basis species
newbases <- as.matrix(expand.grid(ispecies))
allbases <- matrix(basis0$ispecies, nrow = 1)[rep(1, nrow(newbases)), , drop = FALSE]
allbases[, ibasis0] <- newbases
- # calculate affinities of species for all combinations of basis species
+ # Calculate affinities of species for all combinations of basis species
aff.species <- list()
message("mosaic: calculating affinities of species for all ", nrow(allbases), " combinations of the basis species")
- # run backwards so that we put the starting basis species back at the end
+ # Run backwards so that we put the starting basis species back at the end
for(i in nrow(allbases):1) {
- # use logact = 0 for solids 20191111
+ # Get default loga from starting basis species
thislogact <- basis0$logact
+ # Use logact = 0 for solids 20191111
states <- sout$species$state[match(allbases[i, ], sout$species$ispecies)]
icr <- grepl("cr", states)
thislogact[icr] <- 0
+ # Use loga_aq for log(activity) of mosaiced aqueous basis species 20220722
+ if(!is.null(loga_aq)) {
+ iaq <- grep("aq", states)
+ # Loop over sets of mosaiced basis species
+ for(j in 1:length(ibasis0)) {
+ if(is.na(loga_aq[j])) next
+ iaqbasis0 <- intersect(iaq, ibasis0[j])
+ thislogact[iaqbasis0] <- loga_aq[j]
+ }
+ }
put.basis(allbases[i, ], thislogact)
- # we have to define the species using the current basis
+ # We have to define the species using the current basis
species(species0$ispecies, species0$logact)
if(ddd_has_sout) aff.species[[i]] <- suppressMessages(affinity(...))
else aff.species[[i]] <- suppressMessages(affinity(..., sout = sout))
}
- # calculate equilibrium mole fractions for each group of basis species
+ # Calculate equilibrium mole fractions for each group of basis species
group.fraction <- list()
blend <- rep(blend, length(A.bases))
E.bases <- list()
for(i in 1:length(A.bases)) {
if(blend[i] & is.null(stable[i][[1]])) {
- # this isn't needed (and doesn't work) if all the affinities are NA 20180925
+ # This isn't needed (and doesn't work) if all the affinities are NA 20180925
if(any(!sapply(A.bases[[1]]$values, is.na))) {
# 20190504: when equilibrating the changing basis species, use a total activity equal to the activity from the basis definition
# 20191111 use equilibrate(loga.balance = ) instead of setting activities in species definition
e <- equilibrate(A.bases[[i]], loga.balance = as.numeric(basis0$logact[ibasis0[i]]))
- # exponentiate to get activities then divide by total activity
+ # Exponentiate to get activities then divide by total activity
a.equil <- lapply(e$loga.equil, function(x) 10^x)
a.tot <- Reduce("+", a.equil)
group.fraction[[i]] <- lapply(a.equil, function(x) x / a.tot)
- # include the equilibrium activities in the output of this function 20190504
+ # Include the equilibrium activities in the output of this function 20190504
E.bases[[i]] <- e
} else {
group.fraction[[i]] <- A.bases[[i]]$values
}
} else {
- # for blend = FALSE, we just look at whether
+ # For blend = FALSE, we just look at whether
# a basis species predominates within its group
if(is.null(stable[i][[1]])) {
d <- diagram(A.bases[[i]], plot.it = FALSE, limit.water = FALSE)
predom <- d$predominant
} else {
- # get the stable species from the argument 20200715
+ # Get the stable species from the argument 20200715
predom <- stable[i][[1]]
}
group.fraction[[i]] <- list()
for(j in 1:length(bases[[i]])) {
- # if a basis species predominates, it has a mole fraction of 1, or 0 otherwise
+ # If a basis species predominates, it has a mole fraction of 1, or 0 otherwise
yesno <- predom
yesno[yesno != j] <- 0
yesno[yesno == j] <- 1
@@ -149,46 +160,46 @@
}
}
- # make an indexing matrix for all combinations of basis species
+ # Make an indexing matrix for all combinations of basis species
ind.mat <- list()
for(i in 1:length(ispecies)) ind.mat[[i]] <- 1:length(ispecies[[i]])
ind.mat <- as.matrix(expand.grid(ind.mat))
- # loop over combinations of basis species
+ # Loop over combinations of basis species
for(icomb in 1:nrow(ind.mat)) {
- # loop over groups of changing basis species
+ # Loop over groups of changing basis species
for(igroup in 1:ncol(ind.mat)) {
- # get mole fractions for this particular basis species
+ # Get mole fractions for this particular basis species
basisx <- group.fraction[[igroup]][[ind.mat[icomb, igroup]]]
- # loop over species
+ # Loop over species
for(jspecies in 1:length(aff.species[[icomb]]$values)) {
- # get coefficient of this basis species in the formation reaction for this species
+ # Get coefficient of this basis species in the formation reaction for this species
nbasis <- aff.species[[icomb]]$species[jspecies, ibasis0[igroup]]
- # adjust affinity of species for mole fractions (i.e. lower activity) of basis species 20190505
+ # Adjust affinity of species for mole fractions (i.e. lower activity) of basis species 20190505
aff.adjust <- nbasis * log10(basisx)
- # avoid infinite values (from log10(0))
+ # Avoid infinite values (from log10(0))
isfin <- is.finite(aff.adjust)
aff.species[[icomb]]$values[[jspecies]][isfin] <- aff.species[[icomb]]$values[[jspecies]][isfin] + aff.adjust[isfin]
}
- # multiply fractions of basis species from each group to get overall fraction
+ # Multiply fractions of basis species from each group to get overall fraction
if(igroup==1) groupx <- basisx
else groupx <- groupx * basisx
}
- # multiply affinities by the mole fractions of basis species
+ # Multiply affinities by the mole fractions of basis species
aff.species[[icomb]]$values <- lapply(aff.species[[icomb]]$values, function(values) values * groupx)
}
- # get total affinities for the species
+ # Get total affinities for the species
A.species <- aff.species[[1]]
for(i in 1:length(A.species$values)) {
- # extract the affinity contributions from each basis species
+ # Extract the affinity contributions from each basis species
A.values <- lapply(lapply(aff.species, "[[", "values"), "[[", i)
- # sum them to get total affinities for this species
+ # Sum them to get total affinities for this species
A.species$values[[i]] <- Reduce("+", A.values)
}
- # for argument recall, include all arguments in output 20190120
+ # For argument recall, include all arguments in output 20190120
allargs <- c(list(bases = bases, blend = blend), list(...))
- # return the affinities for the species and basis species
+ # Return the affinities for the species and basis species
return(list(fun = "mosaic", args = allargs, A.species = A.species, A.bases = A.bases, E.bases = E.bases))
}
Modified: pkg/CHNOSZ/R/stack_mosaic.R
===================================================================
--- pkg/CHNOSZ/R/stack_mosaic.R 2022-07-22 06:33:37 UTC (rev 736)
+++ pkg/CHNOSZ/R/stack_mosaic.R 2022-07-22 23:52:24 UTC (rev 737)
@@ -4,27 +4,37 @@
# Adapted from vignettes/multi-metal.Rmd
# col: Colors for species1, species2, and species12
# (default of NA for col[3] means to plot species12 boundaries with color for species2)
-# ...: Arguments for mosaic() (incl. affinity() arguments)
+# ...: Arguments for mosaic() (including affinity() arguments)
stack_mosaic <- function(bases, species1, species2, species12, names = NULL, col = list(4, 3, 6), col.names = list(4, 3, 6),
fill = NULL, dx = list(0, 0, 0), dy = list(0, 0, 0), srt = list(0, 0, 0), lwd = list(1, 1, 1), lty = list(1, 1, 1),
- plot.it = TRUE, ...) {
+ loga_aq = NULL, plot.it = TRUE, ...) {
# Default is to use semi-transparent fill for bimetallic species
if(is.null(fill)) fill <- list(NA, NA, add_alpha(col.names[3], "50"))
# Load species1 (first metal-bearing species)
- species(species1)
+ isp1 <- species(species1)
+ # Set activities 20220722
+ if(!is.null(loga_aq)) {
+ iaq1 <- which(isp1$state == "aq")
+ species(iaq1, loga_aq)
+ }
# Calculate affinity of species1 while speciating bases (e.g. aqueous S species)
- mosaic1 <- mosaic(bases, ...)
+ mosaic1 <- mosaic(bases, loga_aq = c(NA, loga_aq), ...)
# Show predominance fields
diagram1 <- diagram(mosaic1$A.species, names = names[[1]], col = col[[1]], col.names = col.names[[1]], fill = fill[[1]],
dx = dx[[1]], dy = dy[[1]], srt = srt[[1]], lwd = lwd[[1]], lty = lty[[1]], plot.it = plot.it)
# Load species12 (bimetallic species) and species2 (second metal-bearing species)
- species(c(species12, species2))
+ isp2 <- species(c(species12, species2))
+ # Set activities 20220722
+ if(!is.null(loga_aq)) {
+ iaq2 <- which(isp2$state == "aq")
+ species(iaq2, loga_aq)
+ }
# Speciate bases again (NULL)
# Take the predominant members of species1 (diagram1$predominant)
- mosaic2 <- mosaic(list(bases, species1), stable = list(NULL, diagram1$predominant), ...)
+ mosaic2 <- mosaic(list(bases, species1), stable = list(NULL, diagram1$predominant), loga_aq = c(NA, loga_aq), ...)
# Set colors
col <- c(rep_len(col[[3]], length(species12)), rep_len(col[[2]], length(species2)))
col.names <- c(rep_len(col.names[[3]], length(species12)), rep_len(col.names[[2]], length(species2)))
Modified: pkg/CHNOSZ/inst/NEWS.Rd
===================================================================
--- pkg/CHNOSZ/inst/NEWS.Rd 2022-07-22 06:33:37 UTC (rev 736)
+++ pkg/CHNOSZ/inst/NEWS.Rd 2022-07-22 23:52:24 UTC (rev 737)
@@ -9,8 +9,10 @@
% subscript and superscript
\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}}}
+% 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-27 (2022-07-15)}{
+\section{Changes in CHNOSZ version 1.9.9-29 (2022-07-23)}{
\subsection{MAJOR CHANGE}{
\itemize{
@@ -55,7 +57,7 @@
(\samp{G}, \samp{S}, and \samp{Cp} at 25 \degC) to formation constants of
aqueous species as a function of temperature.
- \item Add vignette \strong{customizing.Rmd} with description of database
+ \item Add vignette \viglink{customizing} with description of database
format, data-entry conventions, and examples of customizing the
thermodynamic database using \code{add.OBIGT()}, \code{mod.OBIGT()}, and
\code{logB_to_OBIGT()}.
@@ -68,10 +70,15 @@
\item Add \code{stack_mosaic()} to create stacked mosaic diagrams, where
the species formed in the first layer become the basis species for the
- species formed in the second layer. This function implements the main
+ species formed in the second layer. This function serves as a wrapper for the
steps for making mosaic diagrams for bimetallic systems described in
- \strong{multi-metal.Rmd}.
+ \viglink{multi-metal}.
+ \item Add \code{loga_aq} argument to \code{mosiac()} (also present in new
+ function \code{stack_mosaic()}) to control the activity of mosaiced
+ aqueous basis species. This obviates a workaround that was previously
+ used in the Mosaic Stacking 2 section of \strong{multi-metal.Rmd}.
+
}
}
Modified: pkg/CHNOSZ/man/mosaic.Rd
===================================================================
--- pkg/CHNOSZ/man/mosaic.Rd 2022-07-22 06:33:37 UTC (rev 736)
+++ pkg/CHNOSZ/man/mosaic.Rd 2022-07-22 23:52:24 UTC (rev 737)
@@ -7,7 +7,7 @@
}
\usage{
- mosaic(bases, blend = TRUE, stable = list(), ...)
+ mosaic(bases, blend = TRUE, stable = list(), loga_aq = NULL, ...)
}
\arguments{
@@ -14,12 +14,13 @@
\item{bases}{character, basis species to be changed in the calculation, or list, each element of which defines an independent group of changing basis species}
\item{blend}{logical, use relative abundances of basis species?}
\item{stable}{list, previously determined stable species}
+ \item{loga_aq}{numeric, activities of aqueous species (overrides current values in \code{\link{basis}})}
\item{...}{additional arguments to be passed to \code{\link{affinity}}}
}
\details{
-\code{mosaic} can be used to calculate the affinities of formation of species when the relative abundances of the basis species listed in \code{bases} change over the range of conditions, due to e.g. ionization, complexation or redox reactions.
+\code{mosaic} calculates the affinities of formation of species when the relative abundances of the basis species listed in \code{bases} change over the range of conditions, due to e.g. ionization, complexation or redox reactions.
This is a way to \dQuote{speciate the basis species}.
For example, the speciation of sulfur (\samp{SO4-2}, \samp{HSO4-}, \samp{HS-} and \samp{H2S}) as a function of Eh and pH affects the formation affinities, and therefore relative stabilities of iron oxide and sulfide minerals.
Chemical activity diagrams constructed by assembling sub-diagrams corresponding to the predominant (i.e. most stable) basis species can described as \dQuote{mosaic diagrams}.
@@ -39,11 +40,16 @@
This is appropriate when minerals, rather than aqueous species, are used as the changing basis species.
Note, however, that \code{mosaic} is not internally recursive: the stabilities of one group of basis species (e.g. minerals) are not affected by changes in another group (e.g. aqueous species).
-This limitation can be overcome by supplying the calculated stabilities (from the \code{predominant} element of the output of \code{\link{diagram}}) as an element of the list in the \code{stable} argument whose position corresponds to the appropriate group of basis species.
-By using the stabilities of minerals in one calculation to identify the basis species in a subsequent calculation, a series of linked \code{mosaic} diagrams with increasing complexity can be made.
-This is useful for diagrams for minerals with multiple metallic elements.
+(\dQuote{Stacked mosaic diagrams}) are useful for making diagrams for multi-metal systems.
+By using the stable minerals in one calculation as the new basis species in a subsequent calculation, a series of stacked \code{mosaic} diagrams with increasing complexity can be made.
+Specifically, this is done by supplying previously calculated stabilities (from the \code{predominant} element of the output of \code{\link{diagram}}) as an element of the list in the \code{stable} argument whose position corresponds to the appropriate group of basis species.
Note that a value in any position of the \code{stable} list forces \code{blend = FALSE} for the corresponding group of basis species, so there is no need to explicity change the \code{blend} argument.
+The activities of mosaiced basis species in each group are taken from the current \code{\link{basis}} definition.
+Generally it makes sense to set the activity of minerals to 1 (logact = 0) and the activity of aqueous species to some smaller value.
+For mosaic stacking calculations where the mosaiced basis species include both minerals and aqueous species, the \code{loga_aq} argument specifies the activity of aqueous species to be used \emph{in each group}.
+That is, there should be one value of \code{loga_aq} for each group of basis species; use NA to indicate that the activity comes from the current \code{\link{basis}} definition.
+See the Mosaic Stacking 2 section of the vignette \viglink{multi-metal} for an example.
}
\value{
Modified: pkg/CHNOSZ/man/stack_mosaic.Rd
===================================================================
--- pkg/CHNOSZ/man/stack_mosaic.Rd 2022-07-22 06:33:37 UTC (rev 736)
+++ pkg/CHNOSZ/man/stack_mosaic.Rd 2022-07-22 23:52:24 UTC (rev 737)
@@ -12,7 +12,8 @@
stack_mosaic(bases, species1, species2, species12, names = NULL,
col = list(4, 3, 6), col.names = list(4, 3, 6), fill = NULL,
dx = list(0, 0, 0), dy = list(0, 0, 0), srt = list(0, 0, 0),
- lwd = list(1, 1, 1), lty = list(1, 1, 1), plot.it = TRUE, ...)
+ lwd = list(1, 1, 1), lty = list(1, 1, 1), loga_aq = NULL,
+ plot.it = TRUE, ...)
add_alpha(col, alpha)
}
@@ -30,6 +31,7 @@
\item{srt}{label rotation}
\item{lwd}{line width}
\item{lty}{line type}
+ \item{loga_aq}{numeric, activity of aqueous species}
\item{plot.it}{make plots?}
\item{...}{arguments for \code{\link{mosaic}} and \code{\link{affinity}}}
\item{alpha}{character, hexadecimal value of color transparency (alpha)}
@@ -38,8 +40,13 @@
\details{
\code{stack_mosaic} creates a stacked mosaic diagram following steps that are described in detail in the vignette \viglink{multi-metal}.
Briefly, the first layer of the diagram is made by speciating the species in \code{bases} across the diagram to form the first set of species in \code{species1}.
-Then, both \code{bases} \emph{and} \code{species1} (the stable species at each point on the diagram) are used to form the second set of species, consisting of both \code{species2} \emph{and} \code{species12}.
+Then, both \code{bases} \emph{and} \code{species1} (the stable species at each point on the diagram) are used to form the second set of species, including those in both \code{species2} \emph{and} \code{species12}.
+Note that \code{basis} has aqueous S species in the examples provided, and \code{species1} consists of minerals and/or aqueous species with a single metal (e.g. Fe).
+\code{species2} has minerals and/or aqueous species with a second metal (e.g. Cu), and \code{species12} has bimetallic minerals.
+For \dQuote{mixed} diagrams (where \code{species1} or \code{species2} has both minerals and aqueous species), use \code{loga_aq} to set the logarithms of activities of aqueous species.
+Here, only a single value of \code{loga_aq} is needed, unlike in \code{\link{mosaic}}, where a value for each set of basis species is required.
+
The plot parameters \code{col}, \code{col.names}, \code{fill}, \code{dx}, \code{dy}, \code{srt}, \code{lwd}, and \code{lty} should be length-3 lists (not vectors).
The values of elements 1--3 of the list are recycled to the number of species in \code{species1}, \code{species2}, and \code{species12}, respectively.
@@ -51,7 +58,7 @@
\section{Warning}{
The bimetallic species in \code{species12} are shown as part of the second layer, although their formation is sensitive to the presence of stable species in the first layer.
It follows that changing the order of layers (i.e., swapping \code{species1} and \code{species2}) can affect the depiction of mineral assemblages that have \code{species12}.
-Only one of the alternatives is thermodynamically correct, but currently there is no check to determine which one it is.
+It is likely that only one of the alternatives is thermodynamically correct, but currently there is no check to determine which one it is.
}
\value{
Modified: pkg/CHNOSZ/vignettes/multi-metal.Rmd
===================================================================
--- pkg/CHNOSZ/vignettes/multi-metal.Rmd 2022-07-22 06:33:37 UTC (rev 736)
+++ pkg/CHNOSZ/vignettes/multi-metal.Rmd 2022-07-22 23:52:24 UTC (rev 737)
@@ -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 <- TRUE
+hires <- FALSE
res1.lo <- 150
res1.hi <- 256
res1 <- ifelse(hires, res1.hi, res1.lo)
@@ -568,18 +568,6 @@
# Add aqueous species 20210220
species(iCu.aq, logm_aq, add = TRUE)
-# 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_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_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 Fe species as basis species (to speed up calculation) 20210224
@@ -587,7 +575,11 @@
ipredom <- sort(unique(as.numeric(predom)))
for(i in seq_along(ipredom)) predom[dFe$predominant == ipredom[i]] <- i
Fe.predom <- c(Fe.cr, Fe.aq)[ipredom]
-mFeCu <- mosaic(list(S.aq, Fe.predom), pH = pH, O2 = O2, T = T, stable = list(NULL, predom), IS = nacl$IS)
+# Use loga_aq argument to control the activity of aqueous species in mosaic calculation 20220722
+# c(NA, logm_aq) means to use:
+# basis()'s value for logact of aqueous S species
+# logm_aq for logact of aqueous Fe species
+mFeCu <- mosaic(list(S.aq, Fe.predom), pH = pH, O2 = O2, T = T, stable = list(NULL, predom), IS = nacl$IS, loga_aq = c(NA, logm_aq))
# Adjust labels
bold <- c(rep(TRUE, length(FeCu.abbrv)), rep(FALSE, length(Cu.aq)))
More information about the CHNOSZ-commits
mailing list