R/util.misc.R

Defines functions which.balance unitize GHS_Tr Ttr dPdTtr

Documented in dPdTtr GHS_Tr Ttr unitize

# CHNOSZ/util.misc.R
# some utility functions for the CHNOSZ package
# speciate/thermo.R 20051021 jmd

dPdTtr <- function(ispecies) {
  # calculate dP/dT for a phase transition
  # (argument is index of the lower-T phase)
  thermo <- get("thermo")
  t1 <- subcrt(ispecies, P=0, T=thermo$obigt$z.T[ispecies], convert=FALSE, check.Ttr=FALSE)
  t2 <- subcrt(ispecies+1, P=0, T=thermo$obigt$z.T[ispecies], convert=FALSE, check.Ttr=FALSE)
  # if these aren't the same mineral all we can say is zero
  # actually, should be infinity ... the volume change is zero
  if(as.character(t1$species$name) != as.character(t2$species$name)) return(Inf)
  # we really hope the G's are the same ...
  #if(abs(t2$out[[1]]$G - t2$out[[1]]$G) > 0.1) warning('dP.dT: inconsistent values of G for different phases of ',ispecies,call.=FALSE)
  dP.dT <- convert( ( t2$out[[1]]$S - t1$out[[1]]$S ) / ( t2$out[[1]]$V - t1$out[[1]]$V ), 'cm3bar' )
  return(dP.dT)
}

Ttr <- function(ispecies,P=1,dPdT=NULL) {
  # calculate a phase transition temperature
  # taking account of P (from dP/dT)
  T <- get("thermo")$obigt$z.T[ispecies]
  if(is.null(dPdT)) dPdT <- dPdTtr(ispecies)
  return(T + P / dPdT)
}

GHS_Tr <- function(ispecies, Htr) {
  # calculate G, H, and S at Tr for cr2, cr3, ... phases 20170301
  # Htr: enthalpy(ies) of transition
  # ispecies: the species index for cr1
  thisinfo <- info(ispecies)
  name <- thisinfo$name
  # start from Tr (T=298.15 K)
  Tprev <- 298.15
  # the GHS at T
  Gf <- thisinfo$G
  Hf <- thisinfo$H
  S <- thisinfo$S
  # where to store the calculated GHS at Tr
  Gf_Tr <- Hf_Tr <- S_Tr <- numeric()
  for(i in 1:(length(Htr)+1)) {
    # check that we have the correct one of cr1, cr2, ...
    thiscr <- paste0("cr", i)
    if(thisinfo$state!=thiscr | thisinfo$name!=name) stop(paste("species", thisis, "is not", name, thiscr))
    # if we're above cr1, calculate the equivalent GHS at Tr
    if(i > 1) {
      # set the starting GHS to 0 (in case they're NA - we only need the increments over temperature)
      thisinfo$G <- thisinfo$H <- thisinfo$S <- 0
      # the HS increments from 298.15 to Ttr
      HSinc <- cgl(c("H", "S"), T=c(298.15, Ttr), ghs=thisinfo, eos=thisinfo)
      Hf_Tr <- c(Hf_Tr, Hf - diff(HSinc[[1]]$H))
      S_Tr <- c(S_Tr, S - diff(HSinc[[1]]$S))
      # plug in the calculated S_Tr to calculate the G increment correctly
      thisinfo$S <- tail(S_Tr, 1)
      Ginc <- cgl("G", T=c(298.15, Ttr), ghs=thisinfo, eos=thisinfo)
      Gf_Tr <- c(Gf_Tr, Gf - diff(Ginc[[1]]$G))
    }
    # the temperature of the next transition
    Ttr <- thisinfo$T
    # the GHS increments from Tprev to Ttr
    GHCinc <- cgl(c("G", "H", "S"), T=c(Tprev, Ttr), ghs=thisinfo, eos=thisinfo)
    # the GHS + transition at Tr
    Gf <- Gf + diff(GHCinc[[1]]$G)
    Hf <- Hf + diff(GHCinc[[1]]$H) + Htr[i]
    S <- S + diff(GHCinc[[1]]$S) + Htr[i] / Ttr
    # prepare next phase
    thisis <- ispecies + i
    thisinfo <- info(thisis)
    Tprev <- Ttr
  }
  list(Gf_Tr=Gf_Tr, Hf_Tr=Hf_Tr, S_Tr=S_Tr)
}

unitize <- function(logact=NULL,length=NULL,logact.tot=0) {
  # scale the logarithms of activities given in loga
  # so that the logarithm of total activity of residues
  # is zero (i.e. total activity of residues is one),
  # or some other value set in loga.tot.
  # length indicates the number of residues in each species.
  # if loga is NULL, take the logarithms of activities from
  # the current species definition. if any of those species
  # are proteins, get their lengths using protein.length.
  thermo <- get("thermo")
  if(is.null(logact)) {
    if(is.null(thermo$species)) stop("loga is NULL and no species are defined")
    ts <- thermo$species
    logact <- ts$logact
    length <- rep(1,length(logact))
    ip <- grep("_",ts$name)
    if(length(ip) > 0) length[ip] <- protein.length(ts$name[ip])
  }
  # the lengths of the species
  if(is.null(length)) length <- 1
  length <- rep(length,length.out=length(logact)) 
  # remove the logarithms
  act <- 10^logact
  # the total activity
  act.tot <- sum(act*length)
  # the target activity
  act.to.get <- 10^logact.tot
  # the factor to apply
  act.fact <- act.to.get/act.tot
  # apply the factor
  act <- act * act.fact
  # take the logarithms
  log10(act)
  # done!
}

### unexported functions ###

# return, in order, which column(s) of species all have non-zero values.
which.balance <- function(species) {
  # find the first basis species that
  # is present in all species of interest
  # ... it can be used to balance the system
  nbasis <- function(species) return(ncol(species)-4)
  ib <- NULL
  nb <- 1
  nbs <- nbasis(species)
  for(i in 1:nbs) {
    coeff <- species[,i]
    if(length(coeff)==length(coeff[coeff!=0])) {
      ib <- c(ib,i)
      nb <- nb + 1
    } else ib <- c(ib,NA)
  }
  return(ib[!is.na(ib)])
}

Try the CHNOSZ package in your browser

Any scripts or data that you put into this service are public.

CHNOSZ documentation built on July 27, 2017, 3:01 a.m.