Nothing
# solubility.R: vectorized solubility calculations without uniroot
# 20181031 first version jmd
# 20181106 work on the output from affinity(); no "equilibrate()" needed!
# 20190117 add find.IS and test for dissociation reaction
# 20210319 use vector 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:
#source("equilibrate.R")
#source("util.misc.R")
#source("species.R")
#source("util.args.R")
#source("util.character.R")
#source("mosaic.R")
#source("basis.R")
# Function to calculate solubilities of multiple minerals 20210303
# 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)
solubility <- function(iaq, ..., in.terms.of = NULL, dissociate = FALSE, find.IS = FALSE) {
# If iaq is the output of affinity(), use old calling style 20210318
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()
affargs <- ddd <- list(...)
is.mosaic <- FALSE
if(identical(names(ddd)[1], "bases")) {
is.mosaic <- TRUE
# For getting 'sout' from affinity(), drop arguments specific for mosaic()
affargs <- ddd[-1]
affargs <- affargs[!names(affargs) %in% c("bases", "stable", "blend", "loga_aq")]
}
# Save current thermodynamic system settings
thermo <- get("thermo", CHNOSZ)
# Use current basis species as a template for the solubility calculations
ispecies <- basis()$ispecies
logact <- basis()$logact
# The current formed species are the minerals to be dissolved
mineral <- species()
if(is.null(mineral)) stop("please load minerals or gases with species()")
if(!find.IS) {
# Get subcrt() output for all aqueous species and minerals 20210322
# Add aqueous species here - the minerals are already present
lapply(iaq, species, add = TRUE)
# Also add basis species for mosaic()!
if(is.mosaic) lapply(unlist(ddd$bases), species, add = TRUE)
sout <- suppressMessages(do.call(affinity, c(affargs, return.sout = TRUE)))
}
# Make a list to store the calculated solubilities for each mineral
slist <- list()
# Loop over minerals
for(i in seq_along(mineral$ispecies)) {
# Print message
message(paste("solubility: calculating for", mineral$name[i]))
# Define basis species with the mineral first (so it will be dissolved)
ispecies[1] <- mineral$ispecies[i]
logact[1] <- mineral$logact[i]
# Use numeric values first and put in buffer names second (needed for demo/gold.R)
loga.numeric <- suppressWarnings(as.numeric(logact))
basis(ispecies, loga.numeric)
is.na <- is.na(loga.numeric)
if(any(is.na)) basis(rownames(basis())[is.na], logact[is.na])
# Add aqueous species (no need to define activities here - they will be calculated by solubility_calc)
species(iaq)
if(find.IS) {
if(is.mosaic) aout <- suppressMessages(mosaic(...)) else aout <- suppressMessages(affinity(...))
} else {
if(is.mosaic) aout <- suppressMessages(mosaic(..., sout = sout)) else aout <- suppressMessages(affinity(..., sout = sout))
}
# Calculate solubility of this mineral
scalc <- solubility_calc(aout, in.terms.of = in.terms.of, dissociate = dissociate, find.IS = find.IS)
# Store the solubilities in the list
slist[[i]] <- scalc$loga.balance
}
# Restore the original thermodynamic system settings
assign("thermo", thermo, CHNOSZ)
if(length(mineral$ispecies) == 1) {
# For one mineral, return the results of the solubility calculation
scalc
} else {
# For multiple minerals, the overall solubility is the *minimum* among all the minerals
smin <- do.call(pmin, slist)
# Put this into the last-computed 'solubility' object
scalc$loga.balance <- smin
scalc$loga.equil <- slist
scalc$species <- mineral
# Change the function name stored in the object so diagram() plots loga.balance automatically
scalc$fun <- "solubilities"
# Return the object
scalc
}
}
# The "nuts and bolts" of solubility calculations
# Moved from solubility() to solubility_calc() 20210318
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
## Is aout the output from mosaic() instead of affinity()?
aout.save <- aout
thisfun <- aout$fun
if(thisfun == "mosaic") aout <- aout$A.species
# Get starting ionic strength (probably zero, but could be anything set by user)
IS <- aout$IS
if(is.null(IS)) IS <- 0
## To output loga.balance we need the balancing coefficients
bout <- balance(aout)
n.balance <- bout$n.balance
balance <- bout$balance
# Get logarithm of total activity of the conserved basis species
logabfun <- function(loga.equil, n.balance) {
# Exponentiate, multiply by n.balance, sum, logarithm
a.equil <- mapply("^", 10, loga.equil, SIMPLIFY = FALSE)
a.balance <- mapply("*", a.equil, n.balance, SIMPLIFY = FALSE)
a.balance <- Reduce("+", a.balance)
log10(a.balance)
}
# For find.IS = TRUE, iterate to converge on ionic strength
niter <- 1
while(TRUE) {
## The values in aout can be calculated for other than
## unit activities of species, so we have to take away the activites
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
## (can be used to reproduce Fig. 4 of Manning et al., 2013)
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)
# Recalculate things for a dissociation reaction (like CaCO3 = Ca+2 + CO3+2)
if(isTRUE(dissociate)) {
ndissoc <- 2
# The number of dissociated products is the exponent in the activity product
loga.dissoc <- loga.balance / ndissoc
# The contribution to affinity
Adissoc <- lapply(n.balance, "*", loga.dissoc)
# Adjust the affinity and get new equilibrium activities
aout$values <- mapply("-", aout$values, Adissoc, SIMPLIFY = FALSE)
loga.equil <- lapply(1:length(aout$values), Astar)
loga.balance <- logabfun(loga.equil, n.balance)
# Check that the new loga.balance == loga.dissoc
# TODO: does this work for non-1:1 species?
stopifnot(all.equal(loga.balance, loga.dissoc))
}
# Get the old IS
IS.old <- rep(IS, length.out = length(aout$values[[1]]))
# Calculate the ionic strength for unpaired ions: IS = sum(molality * Z^2) / 2
sum.mZ2 <- 0
ion.names <- character()
# is there an ion in the basis species?
if(dissociate) {
basis.ion <- rownames(aout$basis)[2]
Z.basis.ion <- makeup(basis.ion)["Z"]
if(!is.na(Z.basis.ion)) {
sum.mZ2 <- sum.mZ2 + 10^loga.balance * Z.basis.ion^2
ion.names <- c(ion.names, basis.ion)
}
}
# Add ions present in the species of interest
# Instead of using aout$species$name, use info() to get formulas 20190309
species.formulas <- suppressMessages(info(aout$species$ispecies, check.it = FALSE)$formula)
for(i in 1:length(loga.equil)) {
species.ion <- species.formulas[i]
Z.species.ion <- makeup(species.ion)["Z"]
if(!is.na(Z.species.ion)) {
sum.mZ2 <- sum.mZ2 + 10^loga.equil[[i]] * Z.species.ion^2
ion.names <- c(ion.names, species.ion)
}
}
IS <- sum.mZ2 / 2
# Report ions used in the ionic strength calculation
if(find.IS & niter == 1) message("solubility: ionic strength calculated for ", paste(ion.names, collapse = " "))
# Report the current ionic strength
if(find.IS) message("solubility: (iteration ", niter, ") ionic strength range is ", paste(round(range(IS), 4), collapse = " "))
# Stop iterating if we reached the tolerance (or find.IS = FALSE)
if(!find.IS | all(IS - IS.old < 1e-4)) break
# On the first iteration, expand argument values for affinity() or mosaic()
# (e.g. convert pH = c(0, 14) to pH = seq(0, 14, 256) so that it has the same length as the IS values)
# We don't do this if aout$vals is NA 20190731
if(niter == 1 & !all(is.na(aout$vals))) {
if(thisfun == "affinity") for(i in 1:length(aout$vals)) {
aout.save$args[[i]] <- aout$vals[[i]]
}
else if(thisfun == "mosaic") {
for(i in 1:length(aout$vals)) {
argname <- names(aout$args)[i]
aout.save$args[[argname]] <- aout$vals[[i]]
}
}
}
# Recalculate the affinity using the new IS
aout <- suppressMessages(do.call(thisfun, list(aout.save, IS = IS)))
if(thisfun == "mosaic") aout <- aout$A.species
niter <- niter + 1
}
# Do we want the solubility expressed in terms of
# something other than the first basis species? 20190123
if(!is.null(in.terms.of)) {
# write the reaction between the basis species and the new species
sbasis <- species.basis(in.terms.of)
# divide the activity of the conserved basis species by the coefficient in the formation reaction
ibalance <- which.balance(aout$species)
coeff <- sbasis[, ibalance][1]
loga.balance <- loga.balance - log10(coeff)
}
# Make the output (we don't deal with normalized formulas yet, so for now m.balance == n.balance)
# Indicate the function used to make this output 20190525
aout$fun <- "solubility"
# Add names to loga.equil 20190731
names(loga.equil) <- aout$species$name
c(aout, list(balance = bout$balance, m.balance = bout$n.balance, n.balance = bout$n.balance, in.terms.of = in.terms.of,
loga.balance = loga.balance, Astar = loga.equil, loga.equil = loga.equil))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.