Nothing
# CHNOSZ/retrieve.R
# Retrieve species with given elements
# 20190214 initial version
# 20190224 define a chemical system using multiple arguments [defunct; use list() for chemical system]
# 20190304 update the stoichiometric matrix instead of doing a full recalculation when the database changes
# 20190305 use c() for combination of elements, list() for chemical system,
# and add 'ligands' argument to retrieve element-bearing species
retrieve <- function(elements = NULL, ligands = NULL, state = NULL, T = NULL, P = "Psat", add.charge = TRUE, hide.groups = TRUE) {
## Empty argument handling
if(is.null(elements)) return(integer())
## Stoichiometric matrix
# What are the formulas of species in the current database?
formula <- thermo()$OBIGT$formula
# Get a previously calculated stoichiometric matrix
stoich <- thermo()$stoich
# If it doesn't match the current database, update it
if(!identical(rownames(stoich), formula)) {
message("retrieve: updating stoichiometric matrix")
# First put rows in the right order for the current database
istoich <- match(formula, rownames(stoich))
stoich <- stoich[na.omit(istoich), ]
# Deal with any missing formulas
if(any(is.na(istoich))) {
# Convert the stoichiometric matrix to data frame with a 'formula' column
oldstoich <- cbind(formula = rownames(stoich), as.data.frame(stoich))
# Get the stoichiometry of the additional formulas
# and suppress warning messages about missing elements
addstoich <- suppressWarnings(i2A(formula[is.na(istoich)]))
addstoich <- cbind(formula = rownames(addstoich), as.data.frame(addstoich))
# Merge the old and added stoichiometries, and assign 0 for missing elements in either one
newstoich <- merge(oldstoich, addstoich, all = TRUE)
newstoich[is.na(newstoich)] <- 0
# Convert the data frame to matrix with rownames corresponding to formulas
stoich <- as.matrix(newstoich[, 2:ncol(newstoich)])
rownames(stoich) <- newstoich$formula
# Put rows in the right order for the current database
istoich <- match(formula, rownames(stoich))
stoich <- stoich[istoich, ]
}
# Store the stoichiometric matrix for later calculations
thermo("stoich" = stoich)
}
## Generate error for missing element(s)
allelements <- c(unlist(elements), unlist(ligands))
not.present <- ! allelements %in% c(colnames(stoich), "all")
if(any(not.present)) {
if(sum(not.present) == 1) stop('"', allelements[not.present], '" is not an element that is present in any species in the database')
else stop('"', paste(allelements[not.present], collapse = '", "'), '" are not elements that are present in any species in the database')
}
## Handle 'ligands' argument
if(!is.null(ligands)) {
# If 'ligands' is cr, liq, gas, or aq, use that as the state
if(any(ligands %in% c("cr", "liq", "gas", "aq")) & length(ligands) == 1) {
state <- ligands
ispecies <- retrieve(elements, add.charge = add.charge)
} else {
# Include the element in the system defined by the ligands list
ligands <- c(elements, as.list(ligands))
# Call retrieve() for each argument and take the intersection
r1 <- retrieve(elements, add.charge = add.charge)
r2 <- retrieve(ligands, add.charge = add.charge)
ispecies <- intersect(r1, r2)
}
} else {
## Species identification
ispecies <- list()
# Automatically add charge to a system 20190225
if(add.charge & is.list(elements) & !"Z" %in% elements) elements <- c(elements, "Z")
# Proceed element-by-element
for(i in seq_along(elements)) {
element <- unlist(elements[i])
if(identical(element, "all")) {
ispecies[[i]] <- 1:nrow(thermo()$OBIGT)
} else {
# Identify the species that have the element
has.element <- rowSums(stoich[, element, drop = FALSE] != 0) == 1
ispecies[[i]] <- which(has.element)
}
}
# Now we have a list of ispecies (one vector for each element)
# What we do next depends on whether the argument is a list() or c()
if(is.list(elements)) {
# For a chemical system, all species are included that do not contain any other elements
ispecies <- unique(unlist(ispecies))
ielements <- colnames(thermo()$stoich) %in% elements
if(any(!ielements)) {
otherstoich <- thermo()$stoich[, !ielements]
iother <- rowSums(otherstoich[ispecies, ] != 0) > 0
ispecies <- ispecies[!iother]
}
} else {
# Get species that have all the elements; the species must be present in each vector
# Reduce() hint from https://stackoverflow.com/questions/27520310/union-of-intersecting-vectors-in-a-list-in-r
ispecies <- Reduce(intersect, ispecies)
}
}
# Exclude groups
if(hide.groups) {
igroup <- grepl("^\\[.*\\]$", thermo()$OBIGT$name[ispecies])
ispecies <- ispecies[!igroup]
}
# Filter on state
if(!is.null(state)) {
istate <- thermo()$OBIGT$state[ispecies] %in% state
ispecies <- ispecies[istate]
}
# Require non-NA Delta G0 at specific temperature 20200825
if(!is.null(T)) {
G <- sapply(suppressMessages(subcrt(ispecies, T = T, P = P))$out, "[[", "G")
ispecies <- ispecies[!is.na(G)]
}
# Assign names; use e- instead of (Z-1)
names(ispecies) <- thermo()$OBIGT$formula[ispecies]
names(ispecies)[names(ispecies) == "(Z-1)"] <- "e-"
# If there's nothing, don't give it a name
if(length(ispecies) == 0) ispecies <- integer()
ispecies
}
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.