[CHNOSZ-commits] r41 - in pkg/CHNOSZ: . R data inst inst/tests man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Feb 6 16:17:21 CET 2013
Author: jedick
Date: 2013-02-06 16:17:21 +0100 (Wed, 06 Feb 2013)
New Revision: 41
Added:
pkg/CHNOSZ/R/IAPWS95.R
pkg/CHNOSZ/inst/tests/test-IAPWS95.R
pkg/CHNOSZ/inst/tests/test-util.seq.R
pkg/CHNOSZ/man/IAPWS95.Rd
Removed:
pkg/CHNOSZ/inst/tests/test-water.IAPWS95.R
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/iprotein.R
pkg/CHNOSZ/R/subcrt.R
pkg/CHNOSZ/R/util.args.R
pkg/CHNOSZ/R/util.data.R
pkg/CHNOSZ/R/util.expression.R
pkg/CHNOSZ/R/util.fasta.R
pkg/CHNOSZ/R/util.seq.R
pkg/CHNOSZ/R/water.R
pkg/CHNOSZ/data/xxx.R
pkg/CHNOSZ/data/yyy.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/inst/tests/test-iprotein.R
pkg/CHNOSZ/inst/tests/test-subcrt.R
pkg/CHNOSZ/inst/tests/test-util.R
pkg/CHNOSZ/man/eos.Rd
pkg/CHNOSZ/man/iprotein.Rd
pkg/CHNOSZ/man/protein.Rd
pkg/CHNOSZ/man/protein.info.Rd
pkg/CHNOSZ/man/revisit.Rd
pkg/CHNOSZ/man/subcrt.Rd
pkg/CHNOSZ/man/util.fasta.Rd
pkg/CHNOSZ/man/util.seq.Rd
pkg/CHNOSZ/man/water.Rd
pkg/CHNOSZ/src/H2O92D.f
Log:
allow subzero (to -20 degrees C) calculations in H2O92D.f
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2013-02-02 06:30:41 UTC (rev 40)
+++ pkg/CHNOSZ/DESCRIPTION 2013-02-06 15:17:21 UTC (rev 41)
@@ -1,6 +1,6 @@
-Date: 2013-02-02
+Date: 2013-02-06
Package: CHNOSZ
-Version: 0.9-9.3
+Version: 0.9-9.4
Title: Chemical Thermodynamics and Activity Diagrams
Author: Jeffrey M. Dick
Maintainer: Jeffrey M. Dick <jmdick at asu.edu>
Added: pkg/CHNOSZ/R/IAPWS95.R
===================================================================
--- pkg/CHNOSZ/R/IAPWS95.R (rev 0)
+++ pkg/CHNOSZ/R/IAPWS95.R 2013-02-06 15:17:21 UTC (rev 41)
@@ -0,0 +1,299 @@
+# functions for properties of water using
+# the IAPWS-95 formulation (Wagner and Pruss, 2002)
+
+IAPWS95.idealgas <- function(p, delta, tau) {
+ ## the ideal gas part in the IAPWS-95 formulation
+ # from Table 6.1 of Wagner and Pruss, 2002
+ n <- c( -8.32044648201, 6.6832105268, 3.00632, 0.012436,
+ 0.97315, 1.27950, 0.96956, 0.24873 )
+ gamma <- c( NA, NA, NA, 1.28728967,
+ 3.53734222, 7.74073708, 9.24437796, 27.5075105 )
+ # Equation 6.5
+ phi <- function() log(delta) + n[1] + n[2]*tau + n[3]*log(tau) +
+ sum( n[4:8] * log(1-exp(-gamma[4:8]*tau)) )
+ # derivatives from Table 6.4
+ phi.delta <- function() 1/delta+0+0+0+0
+ phi.delta.delta <- function() -1/delta^2+0+0+0+0
+ phi.tau <- function() 0+0+n[2]+n[3]/tau+sum(n[4:8]*gamma[4:8]*((1-exp(-gamma[4:8]*tau))^-1-1))
+ phi.tau.tau <- function() 0+0+0-n[3]/tau^2-sum(n[4:8]*gamma[4:8]^2 *
+ exp(-gamma[4:8]*tau)*(1-exp(-gamma[4:8]*tau))^-2)
+ phi.delta.tau <- function() 0+0+0+0+0
+ return(get(p)())
+}
+
+IAPWS95.residual <- function(p, delta, tau) {
+ ## the residual part in the IAPWS-95 formulation
+ # from Table 6.2 of Wagner and Pruss, 2002
+ c <- c(rep(NA,7),rep(1,15),rep(2,20),rep(3,4),4,rep(6,4),rep(NA,5))
+ d <- c(1,1,1,2,2,3,4,1,1,1,2,2,3,4,
+ 4,5,7,9,10,11,13,15,1,2,2,2,3,4,
+ 4,4,5,6,6,7,9,9,9,9,9,10,10,12,
+ 3,4,4,5,14,3,6,6,6,3,3,3,NA,NA)
+ t <- c(-0.5,0.875,1,0.5,0.75,0.375,1,4,6,12,1,5,4,2,
+ 13,9,3,4,11,4,13,1,7,1,9,10,10,3,
+ 7,10,10,6,10,10,1,2,3,4,8,6,9,8,
+ 16,22,23,23,10,50,44,46,50,0,1,4,NA,NA)
+ n <- c( 0.12533547935523E-1, 0.78957634722828E1 ,-0.87803203303561E1 ,
+ 0.31802509345418 ,-0.26145533859358 ,-0.78199751687981E-2,
+ 0.88089493102134E-2,-0.66856572307965 , 0.20433810950965 ,
+ -0.66212605039687E-4,-0.19232721156002 ,-0.25709043003438 ,
+ 0.16074868486251 ,-0.40092828925807E-1, 0.39343422603254E-6,
+ -0.75941377088144E-5, 0.56250979351888E-3,-0.15608652257135E-4,
+ 0.11537996422951E-8, 0.36582165144204E-6,-0.13251180074668E-11,
+ -0.62639586912454E-9,-0.10793600908932 , 0.17611491008752E-1,
+ 0.22132295167546 ,-0.40247669763528 , 0.58083399985759 ,
+ 0.49969146990806E-2,-0.31358700712549E-1,-0.74315929710341 ,
+ 0.47807329915480 , 0.20527940895948E-1,-0.13636435110343 ,
+ 0.14180634400617E-1, 0.83326504880713E-2,-0.29052336009585E-1,
+ 0.38615085574206E-1,-0.20393486513704E-1,-0.16554050063734E-2,
+ 0.19955571979541E-2, 0.15870308324157E-3,-0.16388568342530E-4,
+ 0.43613615723811E-1, 0.34994005463765E-1,-0.76788197844621E-1,
+ 0.22446277332006E-1,-0.62689710414685E-4,-0.55711118565645E-9,
+ -0.19905718354408 , 0.31777497330738 ,-0.11841182425981 ,
+ -0.31306260323435E2 , 0.31546140237781E2 ,-0.25213154341695E4 ,
+ -0.14874640856724 , 0.31806110878444)
+ alpha <- c(rep(NA,51),20,20,20,NA,NA)
+ beta <- c(rep(NA,51),150,150,250,0.3,0.3)
+ gamma <- c(rep(NA,51),1.21,1.21,1.25,NA,NA)
+ epsilon <- c(rep(NA,51),1,1,1,NA,NA)
+ a <- c(rep(NA,54),3.5,3.5)
+ b <- c(rep(NA,54),0.85,0.95)
+ B <- c(rep(NA,54),0.2,0.2)
+ C <- c(rep(NA,54),28,32)
+ D <- c(rep(NA,54),700,800)
+ A <- c(rep(NA,54),0.32,0.32)
+ # from Table 6.5
+ i1 <- 1:7
+ i2 <- 8:51
+ i3 <- 52:54
+ i4 <- 55:56
+ # deriviatives of distance function
+ Delta <- function(i) { Theta(i)^2 + B[i] * ((delta-1)^2)^a[i] }
+ Theta <- function(i) { (1-tau) + A[i] * ((delta-1)^2)^(1/(2*beta[i])) }
+ Psi <- function(i) { exp ( -C[i]*(delta-1)^2 - D[i]*(tau-1)^2 ) }
+ dDelta.bi.ddelta <- function(i) { b[i]*Delta(i)^(b[i]-1)*dDelta.ddelta(i) }
+ d2Delta.bi.ddelta2 <- function(i) { b[i]*( Delta(i)^(b[i]-1) * d2Delta.ddelta2(i) +
+ (b[i]-1)*Delta(i)^(b[i]-2)*dDelta.ddelta(i)^2 ) }
+ dDelta.bi.dtau <- function(i) { -2*Theta(i)*b[i]*Delta(i)^(b[i]-1) }
+ d2Delta.bi.dtau2 <- function(i) { 2*b[i]*Delta(i)^(b[i]-1) + 4*Theta(i)^2*b[i]*(b[i]-1)*Delta(i)^(b[i]-2) }
+ d2Delta.bi.ddelta.dtau <- function(i) { -A[i]*b[i]*2/beta[i]*Delta(i)^(b[i]-1)*(delta-1) *
+ ((delta-1)^2)^(1/(2*beta[i])-1) - 2*Theta(i)*b[i]*(b[i]-1)*Delta(i)^(b[i]-2)*dDelta.ddelta(i) }
+ dDelta.ddelta <- function(i) { (delta-1) * ( A[i]*Theta(i)*2/beta[i]*((delta-1)^2)^(1/(2*beta[i])-1) +
+ 2*B[i]*a[i]*((delta-1)^2)^(a[i]-1) ) }
+ d2Delta.ddelta2 <- function(i) { 1/(delta-1)*dDelta.ddelta(i) + (delta-1)^2 * (
+ 4*B[i]*a[i]*(a[i]-1)*((delta-1)^2)^(a[i]-2) + 2*A[i]^2*(1/beta[i])^2 *
+ (((delta-1)^2)^(1/(2*B[i])-1))^2 + A[i]*Theta(i)*4/beta[i]*(1/(2*B[i])-1) *
+ ((delta-1)^2)^(1/(2*beta[i])-2) ) }
+ # derivatives of exponential function
+ dPsi.ddelta <- function(i) { -2*C[i]*(delta-1)*Psi(i) }
+ d2Psi.ddelta2 <- function(i) { ( 2*C[i]*(delta-1)^2 - 1 ) * 2*C[i]*Psi(i) }
+ dPsi.dtau <- function(i) { -2*D[i]*(tau-1)*Psi(i) }
+ d2Psi.dtau2 <- function(i) { (2*D[i]*(tau-1)^2 - 1) * 2*D[i]*Psi(i) }
+ d2Psi.ddelta.dtau <- function(i) { 4*C[i]*D[i]*(delta-1)*(tau-1)*Psi(i) }
+ # dimensionless Helmholtz free energy and derivatives
+ phi <- function() {
+ sum(n[i1]*delta^d[i1]*tau^t[i1]) +
+ sum(n[i2]*delta^d[i2]*tau^t[i2]*exp(-delta^c[i2])) +
+ sum(n[i3]*delta^d[i3]*tau^t[i3] *
+ exp( -alpha[i3]*(delta-epsilon[i3])^2 - beta[i3]*(tau-gamma[i3])^2 ) ) +
+ sum(n[i4]*Delta(i4)^b[i4]*delta*Psi(i4))
+ }
+ phi.delta <- function() {
+ sum(n[i1]*d[i1]*delta^(d[i1]-1)*tau^t[i1]) +
+ sum(n[i2]*exp(-delta^c[i2])*(delta^(d[i2]-1)*tau^t[i2]*(d[i2]-c[i2]*delta^c[i2]))) +
+ sum(n[i3]*delta^d[i3]*tau^t[i3] *
+ exp( -alpha[i3]*(delta-epsilon[i3])^2 - beta[i3]*(tau-gamma[i3])^2 ) *
+ (d[i3]/delta - 2 * alpha[i3]*(delta-epsilon[i3])) ) +
+ sum(n[i4] * ( Delta(i4)^b[i4] * (Psi(i4)+delta*dPsi.ddelta(i4)) + dDelta.bi.ddelta(i4)*delta*Psi(i4) ) )
+ }
+ phi.delta.delta <- function() {
+ sum(n[i1]*d[i1]*(d[i1]-1)*delta^(d[i1]-2)*tau^t[i1]) +
+ sum(n[i2]*exp(-delta^c[i2])*(delta^(d[i2]-2)*tau^t[i2]*((d[i2]-c[i2]*delta^c[i2]) *
+ (d[i2]-1-c[i2]*delta^c[i2])-c[i2]^2*delta^c[i2]))) +
+ sum(n[i3]*tau^t[i3]*exp(-alpha[i3]*(delta-epsilon[i3])^2 - beta[i3]*(tau-gamma[i3])^2) * (
+ -2*alpha[i3]*delta^d[i3]+4*alpha[i3]^2*delta^d[i3]*(delta-epsilon[i3])^2 -
+ 4*d[i3]*alpha[i3]*delta^(d[i3]-1)*(delta-epsilon[i3])+d[i3]*(d[i3]-1)*delta^(d[i3]-2) ) ) +
+ sum(n[i4]*( Delta(i4)^b[i4]*(2*dPsi.ddelta(i4)+delta*d2Psi.ddelta2(i4)) +
+ 2*dDelta.bi.ddelta(i4)*(Psi(i4)+delta*dPsi.ddelta(i4)) + d2Delta.bi.ddelta2(i4)*delta*Psi(i4) ) )
+ }
+ phi.tau <- function() {
+ sum(n[i1]*t[i1]*delta^d[i1]*tau^(t[i1]-1)) +
+ sum(n[i2]*t[i2]*delta^d[i2]*tau^(t[i2]-1)*exp(-delta^c[i2])) +
+ sum(n[i3]*delta^d[i3]*tau^t[i3]*exp(-alpha[i3]*(delta-epsilon[i3])^2-beta[i3]*(tau-gamma[i3])^2) *
+ (t[i3]/tau-2*beta[i3]*(tau-gamma[i3]))) +
+ sum(n[i4]*delta*(dDelta.bi.dtau(i4)*Psi(i4)+Delta(i4)^b[i4]*dPsi.dtau(i4)))
+ }
+ phi.tau.tau <- function() {
+ sum(n[i1]*t[i1]*(t[i1]-1)*delta^d[i1]*tau^(t[i1]-2)) +
+ sum(n[i2]*t[i2]*(t[i2]-1)*delta^d[i2]*tau^(t[i2]-2)*exp(-delta^c[i2])) +
+ sum(n[i3]*delta^d[i3]*tau^t[i3]*exp(-alpha[i3]*(delta-epsilon[i3])^2-beta[i3]*(tau-gamma[i3])^2) *
+ (((t[i3]/tau)-2*beta[i3]*(tau-gamma[i3]))^2-t[i3]/tau^2-2*beta[i3])) +
+ sum(n[i4]*delta*(d2Delta.bi.dtau2(i4)*Psi(i4)+2*dDelta.bi.dtau(i4)*dPsi.dtau(i4) +
+ Delta(i4)^b[i4]*d2Psi.dtau2(i4)))
+ }
+ phi.delta.tau <- function() {
+ sum(n[i1]*d[i1]*t[i1]*delta^(d[i1]-1)*tau^(t[i1]-1)) +
+ sum(n[i2]*t[i2]*delta^(d[i2]-1)*tau^(t[i2]-1)*(d[i2]-c[i2]*delta^c[i2])*exp(-delta^c[i2])) +
+ sum(n[i3]*delta^d[i3]*tau^t[i3]*exp(-alpha[i3]*(delta-epsilon[i3])^2-beta[i3]*(tau-gamma[i3])^2) *
+ ((d[i3]/delta)-2*alpha[i3]*(delta-epsilon[i3]))*(t[i3]/tau-2*beta[i3]*(tau-gamma[i3])) ) +
+ sum(n[i4]*(Delta(i4)^b[i4]*(dPsi.dtau(i4)+delta*d2Psi.ddelta.dtau(i4)) +
+ delta*dDelta.bi.ddelta(i4)*dPsi.dtau(i4)+dDelta.bi.dtau(i4) * (Psi(i4)+delta*dPsi.ddelta(i4)) +
+ d2Delta.bi.ddelta.dtau(i4)*delta*Psi(i4) ))
+ }
+ return(get(p)())
+}
+
+IAPWS95 <- function(property,T=298.15,rho=1000) {
+ ## the IAPWS-95 formulation for ordinary water substance
+ ## Wagner and Pruss, 2002
+ property <- tolower(property)
+ # triple point
+ T.triple <- 273.16 # K
+ P.triple <- 611.657 # Pa
+ rho.triple.liquid <- 999.793
+ rho.triple.vapor <- 0.00485458
+ # normal boiling point
+ T.boiling <- 373.124
+ P.boiling <- 0.101325
+ rho.boiling.liquid <- 958.367
+ rho.boiling.vapor <- 0.597657
+ # critical point constants
+ T.critical <- 647.096 # K
+ rho.critical <- 322 # kg m-3
+ # specific and molar gas constants
+ R <- 0.46151805 # kJ kg-1 K-1
+ # R.M <- 8.314472 # J mol-1 K-1
+ # molar mass
+ M <- 18.015268 # g mol-1
+ ## define functions idealgas and residual, supplying arguments delta and tau
+ idealgas <- function(p) IAPWS95.idealgas(p, delta, tau)
+ residual <- function(p) IAPWS95.residual(p, delta, tau)
+ ## relation of thermodynamic properties to Helmholtz free energy
+ a <- function() {
+ x <- idealgas('phi')+residual('phi')
+ return(x*R*T)
+ }
+ # Table 6.3
+ p <- function() {
+ x <- 1 + delta*residual('phi.delta')
+ return(x*rho*R*T/1000) # for MPa
+ }
+ s <- function() {
+ x <- tau * (idealgas('phi.tau')+residual('phi.tau'))-idealgas('phi')-residual('phi')
+ return(x*R)
+ }
+ u <- function() {
+ x <- tau * (idealgas('phi.tau')+residual('phi.tau'))
+ return(x*R*T)
+ }
+ h <- function() {
+ x <- 1 + tau * (idealgas('phi.tau')+residual('phi.tau')) + delta*residual('phi.delta')
+ return(x*R*T)
+ }
+ g <- function() {
+ x <- 1 + idealgas('phi') + residual('phi') + delta*residual('phi.delta')
+ return(x*R*T)
+ }
+ cv <- function() {
+ x <- -tau^2*(idealgas('phi.tau.tau')+residual('phi.tau.tau'))
+ return(x*R)
+ }
+ cp <- function() {
+ x <- -tau^2*(idealgas('phi.tau.tau')+residual('phi.tau.tau')) +
+ (1+delta*residual('phi.delta')-delta*tau*residual('phi.delta.tau'))^2 /
+ (1+2*delta*residual('phi.delta')+delta^2*residual('phi.delta.delta'))
+ return(x*R)
+ }
+# 20090420 speed of sound calculation is incomplete
+# (delta.liquid and drhos.dT not visible)
+# cs <- function() {
+# x <- -tau^2*(idealgas('phi.tau.tau')+residual('phi.tau.tau')) +
+# (1+delta*residual('phi.delta')-delta*tau*residual('phi.delta.tau'))^2 /
+# (1+2*delta*residual('phi.delta')+delta^2*residual('phi.delta.delta')) *
+# ((1+delta.liquid*residual('phi.delta')-delta.liquid*tau*residual('phi.tau.tau'))-rho.critical/(R*delta.liquid)*drhos.dT)
+# return(x*R)
+# }
+ w <- function() {
+ x <- 1 + 2*delta*residual('phi.delta') + delta^2*residual('phi.delta.delta') -
+ (1+delta*residual('phi.delta')-delta*tau*residual('phi.delta.tau'))^2 /
+ tau^2*(idealgas('phi.tau.tau')+residual('phi.tau.tau'))
+ return(sqrt(x*R*T))
+ }
+ mu <- function() {
+ x <- -(delta*residual('phi.delta')+delta^2*residual('phi.delta.delta')+delta*tau*residual('phi.delta.tau')) /
+ ( ( 1+delta*residual('phi.delta')-delta*tau*residual('phi.delta.tau')^2 ) - tau^2 *
+ (idealgas('phi.tau.tau')+residual('phi.tau.tau'))*(1+2*delta*residual('phi.delta')+delta^2*residual('phi.delta.delta')) )
+ return(x/(R*rho))
+ }
+ ## run the calculations
+ ww <- NULL
+ my.T <- T
+ my.rho <- rho
+ for(j in 1:length(property)) {
+ t <- numeric()
+ for(i in 1:length(my.T)) {
+ T <- my.T[i]
+ rho <- my.rho[i]
+ # Equation 6.4
+ delta <- rho / rho.critical
+ tau <- T.critical / T
+ t <- c(t,get(property[j])())
+ }
+ t <- data.frame(t)
+ if(j==1) ww <- t else ww <- cbind(ww,t)
+ }
+ colnames(ww) <- property
+ return(ww)
+}
+
+WP02.auxiliary <- function(property='rho.liquid',T=298.15) {
+ # auxiliary equations for liquid-vapor phase boundary
+ # from Wagner and Pruss, 2002
+ # critical point
+ T.critical <- 647.096 # K
+ P.critical <- 22.064 # MPa
+ rho.critical <- 322 # kg m-3
+
+ if(property %in% c("P.sigma","dP.sigma.dT")) {
+ # vapor pressure
+ V <- 1 - T / T.critical # theta (dimensionless)
+ a1 <- -7.85951783
+ a2 <- 1.84408259
+ a3 <- -11.7866497
+ a4 <- 22.6807411
+ a5 <- -15.9618719
+ a6 <- 1.80122502
+ ln.P.sigma.P.critical <- T.critical / T *
+ ( a1*V + a2*V^1.5 + a3*V^3 + a4*V^3.5 + a5*V^4 + a6*V^7.5 )
+ P.sigma <- P.critical * exp(ln.P.sigma.P.critical)
+ if(property=="dP.sigma.dT") out <- - P.sigma / T * ( ln.P.sigma.P.critical +
+ a1 + 1.5*a2*V^0.5 + 3*a3*V^2 + 3.5*a4*V^2.5 + 4*a5*V^3 + 7.5*a6*V^6.5 )
+ else out <- P.sigma
+ } else if(property=="rho.liquid") {
+ # saturated liquid density
+ V <- 1 - T / T.critical
+ b1 <- 1.99274064
+ b2 <- 1.09965342
+ b3 <- -0.510839303
+ b4 <- -1.75493479
+ b5 <- -45.5170352
+ b6 <- -6.74694450E5
+ rho.liquid <- rho.critical * (
+ 1 + b1*V^(1/3) + b2*V^(2/3) + b3*V^(5/3) + b4*V^(16/3) + b5*V^(43/3) + b6*V^(110/3) )
+ out <- rho.liquid
+ } else if(property=="rho.vapor") {
+ # saturated vapor density
+ V <- 1 - T / T.critical
+ c1 <- -2.03150240
+ c2 <- -2.68302940
+ c3 <- -5.38626492
+ c4 <- -17.2991605
+ c5 <- -44.7586581
+ c6 <- -63.9201063
+ rho.vapor <- rho.critical * exp (
+ c1*V^(2/6) + c2*V^(4/6) + c3*V^(8/6) + c4*V^(18/6) + c5*V^(37/6) + c6*V^(71/6) )
+ out <- rho.vapor
+ } else stop(paste('i can not calculate',property))
+ return(out)
+}
+
Modified: pkg/CHNOSZ/R/iprotein.R
===================================================================
--- pkg/CHNOSZ/R/iprotein.R 2013-02-02 06:30:41 UTC (rev 40)
+++ pkg/CHNOSZ/R/iprotein.R 2013-02-06 15:17:21 UTC (rev 41)
@@ -6,7 +6,6 @@
# ip2aa - select amino acid counts (data frame) from thermo$protein
# aa2eos - perform group additivity calculations
# seq2aa - calculate amino acid counts from a sequence
-# dl.aa - get amino acid counts from SWISS-PROT
# aasum - combine amino acid counts (sum, average, or weighted sum by abundance)
# read.aa - read amino acid counts from a file
# add.protein - add amino acid counts to thermo$protein (returns iprotein)
@@ -116,48 +115,6 @@
return(aa)
}
-dl.aa <- function(protein) {
- # download protein sequence information from SWISS-PROT
- iprotein <- numeric()
- # construct the initial URL
- proteinURL <- paste("http://www.uniprot.org/uniprot/", protein, sep="")
- msgout("dl.aa: trying ", proteinURL, " ...")
- # try loading the URL, hiding any warnings
- oldopt <- options(warn=-1)
- URLstuff <- try(readLines(proteinURL),TRUE)
- options(oldopt)
- if(class(URLstuff)=="try-error") {
- msgout(" failed\n")
- return(NA)
- }
- # 20091102: look for a link to a fasta file
- linkline <- URLstuff[[grep("/uniprot/.*fasta", URLstuff)[1]]]
- # extract accession number from the link
- linkhead <- strsplit(linkline, ".fasta", fixed=TRUE)[[1]][1]
- accession.number <- tail(strsplit(linkhead, "/uniprot/", fixed=TRUE)[[1]], 1)
- msgout(" accession ", accession.number, " ...\n")
- # now download the fasta file
- fastaURL <- paste("http://www.uniprot.org/uniprot/", accession.number, ".fasta", sep="")
- URLstuff <- readLines(fastaURL)
- # show the name of the protein to the user
- header <- URLstuff[[1]]
- header2 <- strsplit(header, paste(protein, ""))[[1]][2]
- header3 <- strsplit(header2, " OS=")[[1]]
- protein.name <- header3[1]
- header4 <- strsplit(header3[2], " GN=")[[1]][1]
- header5 <- strsplit(header4[1], " PE=")[[1]]
- organism.name <- header5[1]
- msgout("dl.aa: ", protein.name, " from ", organism.name)
- # get rid of the header before counting amino acid letters
- URLstuff[[1]] <- ""
- aa <- count.aa(c2s(URLstuff, sep=""))
- msgout(" (length ", sum(aa[1,]), ")\n", sep="")
- colnames(aa) <- colnames(thermo$protein)[6:25]
- po <- strsplit(protein, "_")[[1]]
- out <- data.frame(protein=po[1], organism=po[2], ref=NA, abbrv=NA, chains=1, aa)
- return(out)
-}
-
aasum <- function(aa, abundance=1, average=FALSE, protein=NULL, organism=NULL) {
# returns the sum of the amino acid counts in aa,
# multiplied by the abundances of the proteins
@@ -191,21 +148,14 @@
}
read.aa <- function(file="protein.csv") {
- # if its a fasta file, read the sequences
- if(is.fasta(file)) {
- aa <- read.fasta(file)
- msgout("read.aa: first line in FASTA file is\n")
- msgout(readLines(file, n=1), "\n")
- } else {
- # 20090428 added colClasses here
- aa <- read.csv(file,colClasses=c(rep("character",4),rep("numeric",21)))
- }
+ # 20090428 added colClasses here
+ aa <- read.csv(file,colClasses=c(rep("character",4),rep("numeric",21)))
if(!identical(colnames(aa), colnames(thermo$protein)))
- stop("format of", file, "is incompatible with thermo$protein")
+ stop(paste("format of", file, "is incompatible with thermo$protein"))
return(aa)
}
-add.protein <- function(aa, print.existing=FALSE) {
+add.protein <- function(aa) {
# add a properly constructed data frame of
# amino acid counts to thermo$protein
if(!identical(colnames(aa), colnames(thermo$protein)))
@@ -215,15 +165,13 @@
ip <- suppressMessages(iprotein(po))
ipdup <- !is.na(ip)
# now we're ready to go
- thermo$protein <<- rbind(thermo$protein, aa[!ipdup, ])
+ if(!all(ipdup)) thermo$protein <<- rbind(thermo$protein, aa[!ipdup, ])
+ if(any(ipdup)) thermo$protein[ip[ipdup], ] <<- aa[ipdup, ]
rownames(thermo$protein) <<- NULL
# return the new rownumbers
ip <- iprotein(po)
# make some noise
- msgout("add.protein: added ", nrow(aa)-sum(ipdup), " of ", nrow(aa), " proteins\n")
- if(!all(is.na(ipdup)) & print.existing) {
- potext <- paste(aa$protein[ipdup], aa$organism[ipdup], sep="_", collapse=" ")
- msgout("add.protein: skipped existing ", potext, "\n")
- }
+ if(!all(ipdup)) msgout("add.protein: added ", nrow(aa)-sum(ipdup), " new protein(s) to thermo$protein\n")
+ if(any(ipdup)) msgout("add.protein: replaced ", sum(ipdup), " existing protein(s) in thermo$protein\n")
return(ip)
}
Modified: pkg/CHNOSZ/R/subcrt.R
===================================================================
--- pkg/CHNOSZ/R/subcrt.R 2013-02-02 06:30:41 UTC (rev 40)
+++ pkg/CHNOSZ/R/subcrt.R 2013-02-06 15:17:21 UTC (rev 41)
@@ -180,6 +180,7 @@
}
if(length(P)==1) {
if(can.be.numeric(P)) P.text <- paste(round(as.numeric(P),2),'bar')
+ else P.text <- "P"
} else P.text <- 'P'
#} else P.text <- paste(length(P),'values of P')
if(identical(P[[1]],'Psat')) P.text <- P
Modified: pkg/CHNOSZ/R/util.args.R
===================================================================
--- pkg/CHNOSZ/R/util.args.R 2013-02-02 06:30:41 UTC (rev 40)
+++ pkg/CHNOSZ/R/util.args.R 2013-02-06 15:17:21 UTC (rev 41)
@@ -6,7 +6,6 @@
props <- c('G','H','S','Cp','V','kT','E')
if(eos=='water') {
# things we also get with water
- #props <- c(colnames(thermo$water)[4:length(colnames(thermo$water))])
props <- c(props,'A','U','Cv','Psat','rho','Q','X','Y','epsilon','w')
# they keep on coming: things we also get with SUPCRT92
if(length(agrep(tolower(thermo$opt$water),'supcrt9',max.distance=0.3))>0)
Modified: pkg/CHNOSZ/R/util.data.R
===================================================================
--- pkg/CHNOSZ/R/util.data.R 2013-02-02 06:30:41 UTC (rev 40)
+++ pkg/CHNOSZ/R/util.data.R 2013-02-06 15:17:21 UTC (rev 41)
@@ -318,7 +318,7 @@
# we use the value of X consistent with SUPCRT
X <- -3.055586E-7
refval <- eos$Cp
- calcval <- eos$c1 + eos$c2/(298.15-228)^2 + eos$omega*298.15*X
+ calcval <- eos$c1 + eos$c2/(298.15-thermo$opt$Theta)^2 + eos$omega*298.15*X
tol <- thermo$opt$Cp.tol
units <- "cal K-1 mol-1"
} else if(prop=="V") {
@@ -328,7 +328,7 @@
Q <- 0.00002775729
refval <- eos$V
calcval <- 41.84*eos$a1 + 41.84*eos$a2/2601 +
- (41.84*eos$a3 + 41.84*eos$a4/2601) / (298.15-228) - Q * eos$omega
+ (41.84*eos$a3 + 41.84*eos$a4/2601) / (298.15-thermo$opt$Theta) - Q * eos$omega
tol <- thermo$opt$V.tol
units <- "cm3 mol-1"
}
Modified: pkg/CHNOSZ/R/util.expression.R
===================================================================
--- pkg/CHNOSZ/R/util.expression.R 2013-02-02 06:30:41 UTC (rev 40)
+++ pkg/CHNOSZ/R/util.expression.R 2013-02-06 15:17:21 UTC (rev 41)
@@ -36,7 +36,7 @@
# write a designation of physical state
# use the state given in log if it's a gas or neutral aqueous species
if(log %in% c("g", "gas")) state <- "g"
- if(!"Z" %in% names(elements)) state <- log
+ else if(!"Z" %in% names(elements)) state <- log
if(state != "") {
# subscript it if we're not in a log expression
if(log != "") expr <- substitute(a*group('(',italic(b),')'),list(a=expr, b=state))
@@ -54,8 +54,6 @@
expr <- substitute(a==b, list(a=expr, b=value))
}
}
- # turn the label into an expression, and we're done!
- expr <- as.expression(expr)
return(expr)
}
@@ -66,15 +64,15 @@
propchar <- s2c(property)
expr <- ""
# some special cases
- if(property=="logK") return(expression(log~italic(K)))
+ if(property=="logK") return(quote(log~italic(K)))
# grepl here b/c diagram() uses "loga.equil" and "loga.basis"
- if(grepl("loga", property)) return(expression(log~italic(a)))
- if(property=="alpha") return(expression(alpha))
- if(property=="Eh") return(expression(Eh))
- if(property=="pH") return(expression(pH))
- if(property=="pe") return(expression(pe))
- if(property=="IS") return(expression(IS))
- if(property=="ZC") return(expression(bar(italic(Z))[C]))
+ if(grepl("loga", property)) return(quote(log~italic(a)))
+ if(property=="alpha") return(quote(alpha))
+ if(property=="Eh") return("Eh")
+ if(property=="pH") return("pH")
+ if(property=="pe") return("pe")
+ if(property=="IS") return("IS")
+ if(property=="ZC") return(quote(bar(italic(Z))[C]))
# process each character in the property abbreviation
prevchar <- character()
for(i in 1:length(propchar)) {
@@ -97,7 +95,7 @@
# put it together
expr <- substitute(a*b, list(a=expr, b=thisexpr))
}
- return(as.expression(expr))
+ return(expr)
}
expr.units <- function(property, prefix="", per="mol") {
@@ -135,7 +133,7 @@
if(!any(sapply(c("P", "T", "Eh", "IS"), function(x) grepl(x, property))))
expr <- substitute(a~b^-1, list(a=expr, b=per))
}
- return(as.expression(expr))
+ return(expr)
}
axis.label <- function(label, units=NULL, basis=thermo$basis, prefix="") {
@@ -154,14 +152,14 @@
} else {
# the label is for a chemical property or condition
# make the label by putting a comma between the property and the units
- property <- expr.property(label)[[1]]
- if(is.null(units)) units <- expr.units(label, prefix=prefix)[[1]]
+ property <- expr.property(label)
+ if(is.null(units)) units <- expr.units(label, prefix=prefix)
# no comma needed if there are no units
if(units=="") desc <- substitute(a, list(a=property))
else desc <- substitute(list(a, b), list(a=property, b=units))
}
# done!
- return(as.expression(desc))
+ return(desc)
}
describe.basis <- function(basis=thermo$basis, ibasis=1:nrow(basis), digits=1, oneline=FALSE) {
@@ -194,7 +192,7 @@
for(i in 1:length(property)) {
propexpr <- c(propexpr, expr.property(property[i]))
thisvalue <- format(round(value[i], digits), nsmall=digits)
- thisunits <- expr.units(property[i])[[1]]
+ thisunits <- expr.units(property[i])
thisvalexpr <- substitute(a~b, list(a=thisvalue, b=thisunits))
valexpr <- c(valexpr, as.expression(thisvalexpr))
}
@@ -223,8 +221,8 @@
if(i %in% iname) species <- reaction$name[i]
else {
# should the chemical formula have a state?
- if(identical(states,"all")) species <- expr.species(reaction$formula[i], state=reaction$state[i])[[1]]
- else species <- expr.species(reaction$formula[i])[[1]]
+ if(identical(states,"all")) species <- expr.species(reaction$formula[i], state=reaction$state[i])
+ else species <- expr.species(reaction$formula[i])
}
# get the absolute value of the reaction coefficient
abscoeff <- abs(reaction$coeff[i])
Modified: pkg/CHNOSZ/R/util.fasta.R
===================================================================
--- pkg/CHNOSZ/R/util.fasta.R 2013-02-02 06:30:41 UTC (rev 40)
+++ pkg/CHNOSZ/R/util.fasta.R 2013-02-06 15:17:21 UTC (rev 41)
@@ -1,13 +1,6 @@
# CHNOSZ/util.fasta.R
# read and manipulate FASTA sequence files
-is.fasta <- function(file) {
- # check if the file is in FASTA format
- # read two lines in case the first one is blank
- l <- readLines(file,n=2)
- if(length(grep("^>",l)) == 0) return(FALSE) else return(TRUE)
-}
-
grep.file <- function(file,pattern="",y=NULL,ignore.case=TRUE,startswith=">",lines=NULL,grep="grep") {
# return the line numbers of the file that contain
# the search term x and optionally don't contain y
@@ -51,7 +44,8 @@
return(as.numeric(out))
}
-read.fasta <- function(file,i=NULL,ret="count",lines=NULL,ihead=NULL,pnff=FALSE) {
+read.fasta <- function(file, i=NULL, ret="count", lines=NULL, ihead=NULL,
+ pnff=FALSE, start=NULL, stop=NULL) {
# read sequences from a fasta file
if(file != "") msgout("read.fasta: reading ",basename(file),"\n")
# all of them or only those indicated by i
@@ -83,13 +77,14 @@
if(is.null(ihead)) ihead <- which(substr(lines,1,1)==">")
linefun <- function(i1,i2) lines[i1:i2]
}
+ # identify the lines that begin and end each sequence
if(is.null(i)) {
i <- ihead
- start <- i + 1
+ begin <- i + 1
end <- i - 1
end <- c(end[-1], nlines)
} else {
- start <- i + 1
+ begin <- i + 1
iend <- match(i,ihead)
# we have to be careful about the last record
iend[iend==ihead[length(ihead)]] <- NA
@@ -99,11 +94,12 @@
# just return the lines from the file
if(ret=="fas") {
iline <- numeric()
- for(i in 1:length(start)) iline <- c(iline,(start[i]-1):end[i])
+ for(i in 1:length(begin)) iline <- c(iline,(begin[i]-1):end[i])
return(lines[iline])
}
- seqfun <- function(i) paste(linefun(start[i],end[i]),collapse="")
- sequences <- palply(1:length(i), seqfun)
+ # get each sequences from the begin to end lines
+ seqfun <- function(i) paste(linefun(begin[i],end[i]),collapse="")
+ sequences <- lapply(1:length(i), seqfun)
# process the header line for each entry
# (strip the ">" and go to the first space or underscore)
nomfun <- function(befund) {
@@ -140,7 +136,7 @@
organism <- bnf
}
if(ret=="count") {
- aa <- count.aa(sequences)
+ aa <- count.aa(sequences, start, stop)
colnames(aa) <- aminoacids(3)
ref <- abbrv <- NA
chains <- 1
@@ -150,33 +146,79 @@
} else return(sequences)
}
-splitline <- function(line,length) {
- # to split a line into multiple lines with a specified length
- out <- character()
- count <- 0
- n <- nchar(line)
- while(count < n) {
- split <- substr(line,count+1,count+length)
- out <- c(out,split)
- count <- count + length
+uniprot.aa <- function(protein, start=NULL, stop=NULL) {
+ # download protein sequence information from UniProt
+ iprotein <- numeric()
+ # construct the initial URL
+ proteinURL <- paste("http://www.uniprot.org/uniprot/", protein, sep="")
+ msgout("uniprot.aa: trying ", proteinURL, " ...")
+ # try loading the URL, hiding any warnings
+ oldopt <- options(warn=-1)
+ URLstuff <- try(readLines(proteinURL),TRUE)
+ options(oldopt)
+ if(class(URLstuff)=="try-error") {
+ msgout(" failed\n")
+ return(NA)
}
- return(out)
+ # 20091102: look for a link to a fasta file
+ linkline <- URLstuff[[grep("/uniprot/.*fasta", URLstuff)[1]]]
+ # extract accession number from the link
+ linkhead <- strsplit(linkline, ".fasta", fixed=TRUE)[[1]][1]
+ accession.number <- tail(strsplit(linkhead, "/uniprot/", fixed=TRUE)[[1]], 1)
+ msgout(" accession ", accession.number, " ...\n")
+ # now download the fasta file
+ fastaURL <- paste("http://www.uniprot.org/uniprot/", accession.number, ".fasta", sep="")
+ URLstuff <- readLines(fastaURL)
+ # show the name of the protein to the user
+ header <- URLstuff[[1]]
+ header2 <- strsplit(header, paste(protein, ""))[[1]][2]
+ header3 <- strsplit(header2, " OS=")[[1]]
+ protein.name <- header3[1]
+ header4 <- strsplit(header3[2], " GN=")[[1]][1]
+ header5 <- strsplit(header4[1], " PE=")[[1]]
+ organism.name <- header5[1]
+ msgout("uniprot.aa: ", protein.name, " from ", organism.name)
+ # 20130206 use read.fasta with lines, start, stop arguments
+ aa <- read.fasta(file="", lines=URLstuff, start=start, stop=stop)
+ msgout(" (length ", sum(aa[1, 6:25]), ")\n", sep="")
+ po <- strsplit(protein, "_")[[1]]
+ aa$protein <- po[1]
+ aa$organism <- po[2]
+ return(aa)
}
-trimfas <- function(file,start,stop) {
- # to extract certain positions from an (aligned) fasta file
- lines <- readLines(file)
- fas <- read.fasta(file="",lines=lines,ret="seq")
- # the length of lines to use
- ll <- nchar(lines[2])
- ihead <- grep("^>",lines)
- head <- lines[ihead]
- out <- character()
- for(i in 1:length(head)) {
- extract <- substr(fas[i],start,stop)
- out <- c(out,head[i],splitline(extract,ll),"")
+count.aa <- function(seq, start=NULL, stop=NULL) {
+ # count amino acids in one or more sequences
+ # sequences are given as elements of the list seq
+ aa <- aminoacids(1)
+ # to count the letters in each sequence
+ countfun <- function(seq, start, stop) {
+ count <- numeric(20)
+ # get a substring if one or both of start or stop are given
+ # if only one of start or stop is given, get a default value for the other
+ if(!is.null(start)) {
+ if(is.null(stop)) stop <- nchar(seq)
+ seq <- substr(seq, start, stop)
+ } else if(!is.null(stop)) {
+ seq <- substr(seq, 1, stop)
+ }
+ # the actual counting
+ naa <- table(strsplit(toupper(seq), "")[[1]])
+ # put them in the same order as in aa (i.e. thermo$protein)
+ iaa <- match(names(naa), aa)
+ # in case any letters don't match some amino acid
+ ina <- is.na(iaa)
+ count[iaa[!ina]] <- naa[!ina]
+ if(any(ina)) msgout("count.aa: unrecognized amino acid code(s): ",
+ paste(names(naa)[ina], collapse=" "), "\n")
+ return(count)
}
- #write.table(out,paste(file,"trim",sep="."),row.names=FALSE,col.names=FALSE,quote=FALSE)
- return(out)
+ # count amino acids in each sequence
+ a <- palply(seq, countfun, start, stop)
+ a <- t(as.data.frame(a, optional=TRUE))
+ # clean up row/column names
+ colnames(a) <- aa
+ rownames(a) <- 1:nrow(a)
+ return(a)
}
Modified: pkg/CHNOSZ/R/util.seq.R
===================================================================
--- pkg/CHNOSZ/R/util.seq.R 2013-02-02 06:30:41 UTC (rev 40)
+++ pkg/CHNOSZ/R/util.seq.R 2013-02-06 15:17:21 UTC (rev 41)
@@ -31,32 +31,6 @@
else if(nchar=="Z") return(aacharged[iaa])
}
-count.aa <- function(seq) {
- # count amino acids in one or more sequences
- # sequences are given as elements of the list seq
- aa <- aminoacids(1)
- countfun <- function(seq) {
- # count the letters in each sequence, putting them in the same order as aa
- count <- numeric(20)
- # the actual counting
- naa <- table(strsplit(toupper(seq), "")[[1]])
- iaa <- match(names(naa), aa)
- # in case any letters don't match some amino acid
- ina <- is.na(iaa)
- count[iaa[!ina]] <- naa[!ina]
- if(any(ina)) msgout("count.aa: unrecognized amino acid code(s): ",
- paste(names(naa)[ina], collapse=" "), "\n")
- return(count)
- }
- # count amino acids in each sequence
- a <- palply(seq, countfun)
- a <- t(as.data.frame(a, optional=TRUE))
- # clean up row/column names
- colnames(a) <- aa
- rownames(a) <- 1:nrow(a)
- return(a)
-}
-
nucleicacids <- function(seq=NULL,type="DNA",comp=NULL,comp2=NULL) {
# count bases or compute the formula, e.g.
# n <- nucleicacids(list("AGCT","TTTT")) # a dataframe of counts
Modified: pkg/CHNOSZ/R/water.R
===================================================================
--- pkg/CHNOSZ/R/water.R 2013-02-02 06:30:41 UTC (rev 40)
+++ pkg/CHNOSZ/R/water.R 2013-02-06 15:17:21 UTC (rev 41)
@@ -51,302 +51,6 @@
return(t)
}
-idealgas.IAPWS95 <- function(p, delta, tau) {
- ## the ideal gas part in the IAPWS-95 formulation
- # from Table 6.1 of Wagner and Pruss, 2002
- n <- c( -8.32044648201, 6.6832105268, 3.00632, 0.012436,
- 0.97315, 1.27950, 0.96956, 0.24873 )
- gamma <- c( NA, NA, NA, 1.28728967,
- 3.53734222, 7.74073708, 9.24437796, 27.5075105 )
- # Equation 6.5
- phi <- function() log(delta) + n[1] + n[2]*tau + n[3]*log(tau) +
- sum( n[4:8] * log(1-exp(-gamma[4:8]*tau)) )
- # derivatives from Table 6.4
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 41
More information about the CHNOSZ-commits
mailing list