Nothing
# CHNOSZ/util.formula.R
# Functions to compute some properties of chemical formulas
i2A <- function(formula) {
## Assemble the stoichiometric matrix (A)
## for the given formulas 20120108 jmd
if(is.matrix(formula)) {
# dDo nothing if the argument is already a matrix
A <- formula
} else if(is.numeric(formula) & !is.null(names(formula))) {
# Turn a named numeric object into a formula matrix
A <- t(as.matrix(formula))
} else {
# Get the elemental makeup of each formula, counting
# zero for elements that appear only in other formulas
msz <- makeup(formula, count.zero = TRUE)
# Convert formulas into a stoichiometric matrix with elements on the columns
A <- t(sapply(msz, c))
# Add names from character argument
# or from thermo()$OBIGT for numeric argument
if(is.numeric(formula[1])) rownames(A) <- get("thermo", CHNOSZ)$OBIGT$name[formula]
else rownames(A) <- formula
}
return(A)
}
as.chemical.formula <- function(makeup, drop.zero = TRUE) {
# Make a formula character string from the output of makeup()
# or from a stoichiometric matrix (output of i2A() or protein.formula())
# First define a function to work with a single makeup object
cffun <- function(makeup) {
# First strip zeroes if needed
if(drop.zero) makeup <- makeup[makeup != 0]
# Z always goes at end
makeup <- c(makeup[names(makeup) != "Z"], makeup[names(makeup) == "Z"])
# The elements and coefficients
elements <- names(makeup)
coefficients <- as.character(makeup)
# Any 1's get zapped
coefficients[makeup == 1] <- ""
# Any Z's get zapped (if they're followed by a negative number)
# or turned into a plus sign (to indicate a positive charge)
elements[elements == "Z" & makeup < 0] <- ""
elements[elements == "Z" & makeup >= 0] <- "+"
# Put the elements and coefficients together
formula <- paste(elements, coefficients, sep = "", collapse = "")
# If the formula is uncharged, and the last element has a negative
# coefficient, add an explicit +0 at the end
if(!"Z" %in% names(makeup) & tail(makeup,1) < 0)
formula <- paste(formula, "+0", sep = "")
return(formula)
}
# Call cffun() once for a single makeup, or loop for a matrix
if(is.matrix(makeup)) out <- sapply(1:nrow(makeup), function(i) {
mkp <- makeup[i, ]
return(cffun(mkp))
}) else out <- cffun(makeup)
return(out)
}
mass <- function(formula) {
# Calculate the mass of elements in chemical formulas
thermo <- get("thermo", CHNOSZ)
formula <- i2A(get.formula(formula))
ielem <- match(colnames(formula), thermo$element$element)
if(any(is.na(ielem))) stop(paste("element(s)",
colnames(formula)[is.na(ielem)], "not available in thermo()$element"))
mass <- as.numeric(formula %*% thermo$element$mass[ielem])
return(mass)
}
entropy <- function(formula) {
# Calculate the standard molal entropy at Tref of elements in chemical formulas
thermo <- get("thermo", CHNOSZ)
formula <- i2A(get.formula(formula))
ielem <- match(colnames(formula), thermo$element$element)
if(any(is.na(ielem))) warning(paste("element(s)",
paste(colnames(formula)[is.na(ielem)], collapse = " "), "not available in thermo()$element"))
# Entropy per atom
Sn <- thermo$element$s[ielem] / thermo$element$n[ielem]
# If there are any NA values of entropy, put NA in the matrix, then set the value to zero
# this allows mixed finite and NA values to be calculated 20190802
ina <- is.na(Sn)
if(any(ina)) {
for(i in which(ina)) {
hasNA <- formula[, i] != 0
formula[hasNA, i] <- NA
}
Sn[ina] <- 0
}
entropy <- as.numeric( formula %*% Sn )
# Convert to Joules 20220325
convert(entropy, "J")
}
GHS <- function(formula, G = NA, H = NA, S = NA, T = 298.15, E_units = "J") {
# For all NA in G, H and S, do nothing
# For no NA in G, H and S, do nothing
# For one NA in G, H and S, calculate its value from the other two:
# G - standard molal Gibbs energy of formation from the elements
# H - standard molal enthalpy of formation from the elements
# S - standard molal entropy
# Argument checking
if(!all(diff(sapply(list(formula, G, H, S), length)) == 0))
stop("formula, G, H and S arguments are not same length")
# Calculate Se (entropy of elements)
Se <- entropy(formula)
if(E_units == "cal") Se <- convert(Se, "cal")
# Calculate one of G, H, or S if the other two are given
GHS <- lapply(seq_along(G), function(i) {
G <- G[i]
H <- H[i]
S <- S[i]
Se <- Se[i]
if(is.na(G)) G <- H - T * (S - Se)
else if(is.na(H)) H <- G + T * (S - Se)
else if(is.na(S)) S <- (H - G) / T + Se
list(G, H, S)
})
# Turn the list into a matrix and add labels
GHS <- t(sapply(GHS, c))
colnames(GHS) <- c("G", "H", "S")
rownames(GHS) <- formula
GHS
}
ZC <- function(formula) {
# Calculate average oxidation state of carbon
# from chemical formulas of species
# If we haven't been supplied with a stoichiometric matrix, first get the formulas
formula <- i2A(get.formula(formula))
# Is there carbon there?
iC <- match("C", colnames(formula))
if(is.na(iC)) stop("carbon not found in the stoichiometric matrix")
# The nominal charges of elements other than carbon
# FIXME: add more elements, warn about missing ones
knownelement <- c("H", "N", "O", "S", "Z")
charge <- c(-1, 3, 2, 2, 1)
# Where are these elements in the formulas?
iknown <- match(knownelement, colnames(formula))
# Any unknown elements in formula get dropped with a warning
iunk <- !colnames(formula) %in% c(knownelement, "C")
if(any(iunk)) warning(paste("element(s)", paste(colnames(formula)[iunk], collapse = " "),
"not in", paste(knownelement, collapse = " "), "so not included in this calculation"))
# Contribution to charge only from known elements that are in the formula
formulacharges <- t(formula[, iknown[!is.na(iknown)]]) * charge[!is.na(iknown)]
# Sum of the charges; the arrangement depends on the number of formulas
if(nrow(formula) == 1) formulacharge <- rowSums(formulacharges)
else formulacharge <- colSums(formulacharges)
# Numbers of carbons
nC <- formula[, iC]
# Average oxidation state
ZC <- as.numeric(formulacharge/nC)
return(ZC)
}
### Unexported functions ###
# Accept a numeric or character argument; the character argument can be mixed
# (i.e. include quoted numbers). as.numeric is tested on every value; numeric values
# are then interpreted as species indices in the thermodynamic database (rownumbers of thermo()$OBIGT),
# and the chemical formulas for those species are returned.
# Values that can not be converted to numeric are returned as-is.
get.formula <- function(formula) {
# Return the argument if it's a matrix or named numeric
if(is.matrix(formula)) return(formula)
if(is.numeric(formula) & !is.null(names(formula))) return(formula)
# Return the argument as matrix if it's a data frame
if(is.data.frame(formula)) return(as.matrix(formula))
# Return the values in the argument, or chemical formula(s) for values that are species indices
# For numeric values, get the formulas from those rownumbers of thermo()$OBIGT
i <- suppressWarnings(as.integer(formula))
# We can't have more than the number of rows in thermo()$OBIGT
thermo <- get("thermo", CHNOSZ)
iover <- i > nrow(thermo$OBIGT)
iover[is.na(iover)] <- FALSE
if(any(iover)) stop(paste("species number(s)", paste(i[iover], collapse = " "),
"not available in thermo()$OBIGT"))
# We let negative numbers pass as formulas
i[i < 0] <- NA
# Replace any species indices with formulas from thermo()$OBIGT
formula[!is.na(i)] <- thermo$OBIGT$formula[i[!is.na(i)]]
return(formula)
}
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.