R/util.misc.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
# 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)])
}

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.