R/basis.R

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
# CHNOSZ/basis.R
# set up the basis species of a thermodynamic system

basis <- function(species=NULL, state=NULL, logact=NULL, delete=FALSE) {
  thermo <- get("thermo")
  ## delete the basis species if requested
  oldbasis <- thermo$basis
  if(delete) {
    thermo$basis <- NULL
    assign("thermo", thermo, "CHNOSZ")
  }
  ## return the basis definition if requested
  if(is.null(species)) return(oldbasis)
  ## from now on we need something to work with
  if(length(species)==0) stop("species argument is empty")
  # is the species one of the preset keywords?
  if(species[1] %in% preset.basis()) return(preset.basis(species[1]))
  # the species names/formulas have to be unique
  if(!length(unique(species))==length(species)) stop("species names are not unique")
  ## processing 'state' and 'logact' arguments
  # they should be same length as species
  if(!is.null(state)) state <- rep(state, length.out=length(species))
  if(!is.null(logact)) logact <- rep(logact, length.out=length(species))
  # results should be identical for
  # basis(c('H2O','CO2','H2'), rep('aq',3), c(0,-3,-3))
  # basis(c('H2O','CO2','H2'), c(0,-3,-3), rep('aq',3))
  # first of all, do we have a third argument?
  if(!is.null(logact)) {
    # does the 3rd argument look like states?
    if(is.character(logact[1])) {
      # swap the arguments into their correct places
      tmp <- logact
      logact <- state
      state <- tmp
    }
  } else {
    # if the second argument is numeric, treat it like logacts
    if(is.numeric(state[1])) {
      logact <- state
      state <- NULL
    }
  }
  ## processing 'species' argument
  # pH transformation
  if("pH" %in% species) {
    logact[species=="pH"] <- -logact[species=="pH"]
    if(!is.null(logact)) species[species=="pH"] <- "H+"
  }
  # Eh and pe transformations
  if("pe" %in% species) {
    logact[species=="pe"] <- -logact[species=="pe"]
    if(!is.null(logact)) species[species=="pe"] <- "e-"
  }
  if("Eh" %in% species) {
    # 20090209 should be careful with this conversion as it's only for 25 deg C
    # to be sure, just don't call species("Eh")
    if(!is.null(logact)) logact[species=="Eh"] <- -convert(logact[species=="Eh"],"pe")
    species[species=="Eh"] <- "e-"
  }
  ## if all species are in the existing basis definition, 
  ## *and* at least one of state or logact is not NULL
  ## modify the states and/or logacts of the existing basis species
  if(all(species %in% rownames(oldbasis)) | all(species %in% oldbasis$ispecies)) 
    if(!is.null(state) | !is.null(logact))
      return(mod.basis(species, state, logact))
  ## we're on to making a new basis definition
  # use default logacts if they aren't present
  if(is.null(logact)) logact <- rep(0, length(species))
  # if species argument is numeric, it's species indices
  if(is.numeric(species[1])) {
    ispecies <- species
    ina <- ispecies > nrow(thermo$obigt)
  } else {
    # get species indices using states from the argument, or default states
    if(!is.null(state)) ispecies <- suppressMessages(info(species, state, check.it=FALSE))
    else ispecies <- suppressMessages(info(species, check.it=FALSE))
    # check if we got all the species
    ina <- is.na(ispecies)
    # info() returns a list if any of the species had multiple approximate matches
    # we don't accept any of those
    if(is.list(ispecies)) ina <- ina | sapply(ispecies,length) > 1
  }
  if(any(ina)) stop(paste("species not available:", paste(species[ina], "(", state[ina], ")", sep="", collapse=" ")))
  # load new basis species
  return(put.basis(ispecies, logact))
}

### unexported functions ###

# to add the basis to thermo$obigt
put.basis <- function(ispecies, logact = rep(NA, length(ispecies))) {
  thermo <- get("thermo")
  state <- thermo$obigt$state[ispecies]
  # make the basis matrix, revised 20120114
  # get the elemental makeup of each species,
  # counting zero for any element that only appears in other species in the set
  comp <- makeup(ispecies, count.zero=TRUE)
  # turn the list into a matrix
  comp <- sapply(comp, c)
  # transpose to get put basis species on the rows
  comp <- t(comp)
  # note, makeup(count.zero=TRUE) above gave elements (colnames) sorted alphabetically
  # rownames identify the species
  rownames(comp) <- as.character(thermo$obigt$formula[ispecies])
  # FIXME: the electron doesn't look like a chemical formula
  # this is needed for affinity() to understand a 'pe' or 'Eh' variable
  if("Z0-1" %in% rownames(comp)) rownames(comp)[rownames(comp)=="Z0-1"] <- "e-"
  # now check it for validity of basis species
  # the first test: matrix is square
  if( nrow(comp) > ncol(comp) ) stop("overdetermined system; square stoichiometric matrix needed")
  if( nrow(comp) < ncol(comp) ) stop("underdetermined system; square stoichiometric matrix needed")
  # the second test: matrix is invertible
  if(class(try(solve(comp), silent=TRUE))=='try-error') 
    stop("singular stoichiometric matrix")
  # store the basis definition in thermo$basis, including
  # both numeric and character data, so we need to use a data frame
  comp <- cbind(as.data.frame(comp), ispecies, logact, state, stringsAsFactors=FALSE)
  # ready to assign to the global thermo object
  thermo$basis <- comp
  assign("thermo", thermo, "CHNOSZ")
  # remove the species since there's no guarantee the
  # new basis includes all their elements
  species(delete=TRUE)
  return(thermo$basis)
}

# modify the states or logact values in the existing basis definition
mod.basis <- function(species, state=NULL, logact=NULL) {
  thermo <- get("thermo")
  # the basis must be defined
  if(is.null(thermo$basis)) stop("basis is not defined")
  # loop over each species to modify
  for(i in 1:length(species)) {
    # which basis species are we looking at?
    if(is.numeric(species)) {
      ib <- match(species[i], thermo$basis$ispecies)
      if(is.na(ib)) stop(paste(species[i],"is not a species index of one of the basis species"))
    } else {
      ib <- match(species[i], rownames(thermo$basis))
      if(is.na(ib)) stop(paste(species[i],"is not a formula of one of the basis species"))
    }
    # first modify the state
    if(!is.null(state)) {
      if(state[i] %in% thermo$buffer$name) {
        # this is the name of a buffer
        ibuff <- which(as.character(thermo$buffers$name)==state[i])
        # check that each species in the buffer is compositionally
        # compatible with the basis definition
        for(k in 1:length(ibuff)) {
          ispecies <- suppressMessages(info(as.character(thermo$buffers$species)[ibuff[k]],
            as.character(thermo$buffers$state)[ibuff[k]], check.it=FALSE))
          bufmakeup <- makeup(ispecies)
          inbasis <- rownames(bufmakeup) %in% colnames(basis()) 
          if(FALSE %in% inbasis) {
            stop(paste("the elements '",c2s(rownames(bufmakeup)[!inbasis]),
              "' of species '",thermo$buffers$species[ibuff[k]],"' in buffer '",state[i],
              "' are not in the basis\n",sep=""))
          }
        }
        thermo$basis$logact[ib] <- state[i]
      } else {
        # look for a species with this name in the requested state
        ispecies <- suppressMessages(info(thermo$obigt$name[thermo$basis$ispecies[ib]], state[i], check.it=FALSE))
        if(is.na(ispecies) | is.list(ispecies)) 
          stop(paste("state or buffer '", state[i], "' not found for ", thermo$obigt$name[thermo$basis$ispecies[ib]], "\n", sep=""))
        thermo$basis$ispecies[ib] <- ispecies
        thermo$basis$state[ib] <- state[i]
      }
    } 
    # then modify the logact
    if(!is.null(logact)) thermo$basis$logact[ib] <- as.numeric(logact[i])
    # assign the result to the CHNOSZ environment
    assign("thermo", thermo, "CHNOSZ")
  }
  return(thermo$basis)
} 

# to load a preset basis definition by keyword
preset.basis <- function(key=NULL) {
  # the available keywords
  basis.key <- c("CHNOS", "CHNOS+", "CHNOSe", "CHNOPS+", "MgCHNOPS+", "FeCHNOS", "FeCHNOS+", "QEC")
  # just list the keywords if none is specified
  if(is.null(key)) return(basis.key)
  # delete any previous basis definition
  basis(delete=TRUE)
  # match the keyword to the available ones
  ibase <- match(key, basis.key)
  if(is.na(ibase)) stop(paste(key, "is not a keyword for preset basis species"))
  if(ibase==1) species <- c("CO2", "H2O", "NH3", "H2S", "oxygen")
  else if(ibase==2) species <- c("CO2", "H2O", "NH3", "H2S", "oxygen", "H+")
  else if(ibase==3) species <- c("CO2", "H2O", "NH3", "H2S", "e-", "H+")
  else if(ibase==4) species <- c("CO2", "H2O", "NH3", "H3PO4", "H2S", "e-", "H+")
  else if(ibase==5) species <- c("Mg+2", "CO2", "H2O", "NH3", "H3PO4", "H2S", "e-", "H+")
  else if(ibase==6) species <- c("Fe2O3", "CO2", "H2O", "NH3", "H2S", "oxygen")
  else if(ibase==7) species <- c("Fe2O3", "CO2", "H2O", "NH3", "H2S", "oxygen", "H+")
  else if(ibase==8) species <- c("glutamine", "glutamic acid", "cysteine", "H2O", "oxygen")
  # get the preset logact
  logact <- preset.logact(species)
  # load the species and return the result
  return(basis(species, logact))
}

# logarithms of activities for preset basis definitions
preset.logact <- function(species) {
  bases <- c("H2O", "CO2", "NH3", "H2S", "oxygen", "H+", "e-", "Fe2O3",
             "glutamine", "glutamic acid", "cysteine")
  # values for QEC amino acids from Dick, 2017 (http://doi.org/10.1101/097667)
  logact <- c(0, -3, -4, -7, -80, -7, -7, 0,
              -3.2, -4.5, -3.6)
  ibase <- match(species, bases)
  logact <- logact[ibase]
  # any unmatched species gets a logarithm of activity of -3
  logact[is.na(logact)] <- -3
  return(logact)
}

Questions? Problems? Suggestions? or email at ian@mutexlabs.com.

Please suggest features or report bugs with the GitHub issue tracker.

All documentation is copyright its authors; we didn't write any of that.