Next: info
Up: CHNOSZ examples
Previous: thermo
exmpls> ## Don't show:
exmpls> data(thermo)
thermo: loaded 1997 aqueous, 3089 total species to thermo$obigt
thermo: loaded 5264 proteins to thermo$ECO
thermo: loaded 6717 proteins to thermo$SGD
thermo: loaded 4155 localizations and 3570 abundances to thermo$yeastgfp
exmpls> ## End Don't show
exmpls> ## string to character
exmpls> s2c("hello world")
[1] "h" "e" "l" "l" "o" " " "w" "o" "r" "l" "d"
exmpls> s2c("hello world",sep=" ",keep.sep=FALSE)
[1] "hello" "world"
exmpls> s2c("3.141592",sep=c(".","9"))
[1] "3" ".1415" "92"
exmpls> s2c("3.141592",sep=c(".","9"),move.sep=TRUE)
[1] "3." "14159" "2"
exmpls> # character to string
exmpls> c2s(aminoacids())
[1] "A C D E F G H I K L M N P Q R S T V W Y"
exmpls> c2s(aminoacids(),sep=".")
[1] "A.C.D.E.F.G.H.I.K.L.M.N.P.Q.R.S.T.V.W.Y"
exmpls> ## Spearman's rho
exmpls> spearman(c(1,2,3),c(2,3,4)) # 1
[1] 1
exmpls> spearman(c(1,2,3),c(4,3,2)) # -1
[1] -1
exmpls> ## argument processing
exmpls> eos.args("hkf",c("g","H","S","cP","V","kT","e"))
$props
[1] "G" "H" "S" "Cp" "V" "kT" "E"
$prop
[1] "g" "h" "s" "cp" "v" "kt" "e"
$Prop
[1] "G" "H" "S" "Cp" "V" "kT" "E"
exmpls> ## produces an error because "Q" is not allowed in water.SUPCRT92
exmpls> ## Not run: eos.args("hkf",c("G","H","S","Cp","V","kT","E","Q"))
exmpls> thermo$opt$water <- "IAPWS" # needed for p and n in next line
exmpls> eos.args("water",c("p","u","cv","psat","rho","n","q","x","y","epsilon"))
$props
[1] "G" "H" "S" "Cp" "V" "kT" "E" "A" "U"
[10] "Cv" "Psat" "rho" "Q" "X" "Y" "epsilon" "w" "P"
[19] "N" "UBorn" "de.dT" "de.dP"
$prop
[1] "p" "u" "cv" "psat" "rho" "n" "q" "x" "y"
[10] "epsilon"
$Prop
[1] "P" "U" "Cv" "Psat" "rho" "N" "Q" "X" "Y"
[10] "epsilon"
exmpls> TP.args(c(273.15,373.15))
$T
[1] 273.16 373.15
$P
NULL
exmpls> TP.args(c(273.15,373.15),"Psat")
$T
[1] 273.16 373.15
$P
[1] 1.00000 1.01418
exmpls> TP.args(c(273.15,373.15),c(100,100,200,200))
$T
[1] 273.16 373.15 273.16 373.15
$P
[1] 100 100 200 200
exmpls> state.args(c("AQ","GAS"))
[1] "AQ" "GAS"
exmpls> state.args(c("a","l","liq"))
[1] "aq" "liq" "liq"
exmpls> ## converting among Gibbs, enthalpy, entropy
exmpls> GHS("H") # entropy of H (element)
[1] 15.61663
exmpls> # calculate enthalpy of formation of arsenopyrite
exmpls> GHS("FeAsS",DG=-33843,S=68.5)
[1] -20149.05
exmpls> # return the value of DG calculated from DH and S
exmpls> # cf. -56687.71 from subcrt("water")
exmpls> GHS("H2O",DH=-68316.76,S=16.7123)
[1] -56677.81
exmpls> ## count selected elements in a formula
exmpls> t <- makeup("H2O")
exmpls> expand.formula(c("H","O"),t)
[1] 2 1
exmpls> expand.formula(c("C","H","S"),t)
[1] 0 2 0
exmpls> ## count amino acids in a sequence
exmpls> aminoacids("GGSGG")
A C D E F G H I K L M N P Q R S T V W Y
1 0 0 0 0 0 4 0 0 0 0 0 0 0 0 0 1 0 0 0 0
exmpls> aminoacids("WhatAmIMadeOf?")
A C D E F G H I K L M N P Q R S T V W Y
1 3 0 1 1 1 0 1 1 0 0 2 0 0 0 0 0 1 0 1 0
exmpls> ## count nucleobases in a sequence
exmpls> nucleicacids("ACCGGGTTT")
A C G T
1 1 2 3 3
exmpls> # the DNA complement of that sequence
exmpls> nucleicacids("ACCGGGTTT",comp="DNA")
T G C A
1 1 2 3 3
exmpls> # the RNA complement of the DNA complement
exmpls> n <- nucleicacids("ACCGGGTTT",comp="DNA",comp2="RNA")
exmpls> # the formula of the RNA complement
exmpls> nucleicacids(n,type="RNA")
info: uracil (C4H4N2O2) available in aq, cr.
info: guanine (C5H5N5O) available in aq, cr.
info: cytosine (C4H5N3O) available in aq, cr.
info: adenine (C5H5N5) available in aq, cr.
info: 1709 refers to uracil, C4H4N2O2 aq (LH06a, 3.Feb.03).
info: 1706 refers to guanine, C5H5N5O aq (LH06a, 3.Feb.03).
info: 1707 refers to cytosine, C4H5N3O aq (LH06a, 3.Feb.03).
info: 1705 refers to adenine, C5H5N5 aq (LH06a, 3.Feb.03).
[1] "C41H44N36O7"
exmpls> ## calculate protein length
exmpls> protein.length("LYSC_CHICK")
[1] 129
exmpls> # another way to do it
exmpls> basis("CHNOS")
C H N O S ispecies logact state
CO2 1 0 0 2 0 69 -3 aq
H2O 0 2 0 1 0 1 0 liq
NH3 0 3 1 0 0 68 -4 aq
H2S 0 2 0 0 1 70 -7 aq
O2 0 0 0 2 0 2852 -80 gas
exmpls> species("LYSC_CHICK")
protein: found LYSC_CHICK (C613H959N193O185S10, 129 residues)
exmpls> protein.length(species()$ispecies)
[1] 129
exmpls> # another way to do it
exmpls> ip <- protein("LYSC","CHICK")
exmpls> protein.length(-ip)
[1] 129
exmpls> ## heat capacity as a function of temperature
exmpls> ## (Makhatadze & Privalov, 1990) units: J mol-1
exmpls> MP90.cp(c(5,25,50,75,100,125),"LYSC_CHICK")
[1] 27536.9 29714.4 31570.3 32239.3 32967.6 33114.8
exmpls> ## properties of phase transitions
exmpls> t <- info("enstatite")
info: enstatite (MgSiO3) available in cr1, cr2, cr3.
info: 2094 refers to enstatite, MgSiO3 cr1 (HDN+78, 5.May.78).
exmpls> # (dP/dT) of transitions
exmpls> dPdTtr(t) # first transition
[1] 384.5063
exmpls> dPdTtr(t+1) # second transition
[1] 11.90108
exmpls> # temperature of transitions (Ttr) as a function of P
exmpls> Ttr(t,P=c(1,10,100,1000))
[1] 903.0026 903.0260 903.2601 905.6007
exmpls> Ttr(t,P=c(1,10,100,1000))
[1] 903.0026 903.0260 903.2601 905.6007
exmpls> ## nominal carbon oxidation states
exmpls> ZC("CHNOSZ")
[1] 7
exmpls> t <- info(info("LYSC_CHICK"))
info: 3090 refers to LYSC_CHICK, C613H959N193O185S10 aq (BBA+03).
exmpls> ZC(t$formula)
[1] 0.01631321
exmpls> ## the basis stoichiometry of a made-up species
exmpls> # warns because P isn't in our basis
exmpls> basis("CHNOS")
C H N O S ispecies logact state
CO2 1 0 0 2 0 69 -3 aq
H2O 0 2 0 1 0 1 0 liq
NH3 0 3 1 0 0 68 -4 aq
H2S 0 2 0 0 1 70 -7 aq
O2 0 0 0 2 0 2852 -80 gas
exmpls> basis.comp("SPONCH")
basis.comp: missing element P from species SPONCH.
CO2 H2O NH3 H2S O2
missing 1 -2 1 1 0.5
exmpls> ## describing the basis species
exmpls> basis("CHNOSe")
C H N O S Z ispecies logact state
CO2 1 0 0 2 0 0 69 -3 aq
H2O 0 2 0 1 0 0 1 0 liq
NH3 0 3 1 0 0 0 68 -4 aq
H2S 0 2 0 0 1 0 70 -7 aq
e- 0 0 0 0 0 -1 2 -7 aq
H+ 0 1 0 0 0 1 3 -7 aq
exmpls> describe(thermo$basis)
[1] "25 C, 1 bar, aCO2(-3), aH2O(0), aNH3(-4), aH2S(-7), ae-(-7), aH+(-7)"
exmpls> describe(thermo$basis,T=NULL,P=NULL)
[1] "aCO2(-3), aH2O(0), aNH3(-4), aH2S(-7), ae-(-7), aH+(-7)"
exmpls> ## mass and entropy of compounds of elements
exmpls> element("CH4")
mass entropy
1 16.04276 63.83843
exmpls> element(c("CH4","H2O"),"mass")
mass
1 16.04276
2 18.01528
exmpls> element("Z") # charge
mass entropy
1 0 -15.61663
exmpls> # same mass, opposite entropy as charge
exmpls> element("Z-1") # i.e., electron
mass entropy
1 0 15.61663
exmpls> ## scale logarithms of activity
exmpls> ## suppose we have two proteins whose lengths
exmpls> ## are 100 and 200; what are the logarithms
exmpls> ## of activity of the proteins that are equal to
exmpls> ## each other and that give a total activity of
exmpls> ## residues equal to unity?
exmpls> logact <- c(-3,-3) # could be any two equal numbers
exmpls> length <- c(100,200)
exmpls> logact.tot <- 0
exmpls> loga <- unitize(logact,length,logact.tot)
exmpls> # the proteins have equal activity
exmpls> stopifnot(identical(loga[1],loga[2]))
exmpls> # the sum of activity of the residues is unity
exmpls> stopifnot(isTRUE(all.equal(sum(10^loga * length),1)))
exmpls> ## now, what if the activity of protein 2 is ten
exmpls> ## times that of protein 1?
exmpls> logact <- c(-3,-2)
exmpls> loga <- unitize(logact,length,logact.tot)
exmpls> # the proteins have unequal activity
exmpls> stopifnot(isTRUE(all.equal(loga[2]-loga[1],1)))
exmpls> # but the activities of residues still add up to one
exmpls> stopifnot(isTRUE(all.equal(sum(10^loga * length),1)))
exmpls> ## modify/add species
exmpls> info(t <- info("alanine","cr"))
info: 2315 refers to alanine, C3H7NO2 cr (HOK+98, 30.Aug.06).
name abbrv formula state source1 source2 date G H S Cp V a
2315 alanine <NA> C3H7NO2 cr HOK+98 <NA> 30.Aug.06 -88440 -134500 30.88 29.21 107.556 5.772
b c d e f lambda T
2315 0.07862 0 0 0 0 0 587
exmpls> mod.obigt("alanine",state="cr",G=0,H=0,S=0)
mod.obigt: updating alanine cr (2315).
exmpls> # now the values of G, H, and S are inconsistent
exmpls> # with the elemental composition of alanine
exmpls> info(t)
info: G (from H and S) of alanine cr differs by 55266 from tabulated value.
name abbrv formula state source1 source2 date G H S Cp V a b c d e
2315 alanine <NA> C3H7NO2 cr USER <NA> 28.Nov.09 0 0 0 29.21 107.556 5.772 0.07862 0 0 0
f lambda T
2315 0 0 587
exmpls> # add a species
exmpls> mod.obigt("myname",formula="CHNOSZ",G=0,H=0)
mod.obigt: adding myname aq (3091).
exmpls> info(t <- info("myname"))
info: 3091 refers to myname, CHNOSZ aq (USER, 28.Nov.09).
info: S of myname aq is NA; set to 56.45.
name abbrv formula state source1 source2 date G H S Cp V a1 a2 a3 a4 c1 c2
3091 myname <NA> CHNOSZ aq USER <NA> 28.Nov.09 0 0 56.44706 NA NA NA NA NA NA NA NA
omega Z
3091 NA NA
exmpls> # values of G, H, S as a function of T
exmpls> # (without any equations of state parameters)
exmpls> subcrt(t)
subcrt: 1 species at 15 values of T and P (wet)
$species
name formula state ispecies
3091 myname CHNOSZ aq 3091
$out
$out$myname
T P rho logK G H S V Cp
1 0.01 1.000000 0.9998431 -1.1285822 1410.612 0 56.44706 NA NA
2 25.00 1.000000 0.9970470 0.0000000 0.000 0 56.44706 NA NA
3 50.00 1.000000 0.9880345 0.9543769 -1411.177 0 56.44706 NA NA
4 75.00 1.000000 0.9748423 1.7716897 -2822.353 0 56.44706 NA NA
5 100.00 1.014180 0.9583491 2.4794872 -4233.530 0 56.44706 NA NA
6 125.00 2.322344 0.9390238 3.0983990 -5644.706 0 56.44706 NA NA
7 150.00 4.761587 0.9170077 3.6441792 -7055.883 0 56.44706 NA NA
8 175.00 8.926011 0.8922826 4.1290669 -8467.059 0 56.44706 NA NA
9 200.00 15.549392 0.8646581 4.5627142 -9878.236 0 56.44706 NA NA
10 225.00 25.497499 0.8337473 4.9528357 -11289.412 0 56.44706 NA NA
11 250.00 39.762043 0.7988943 5.3056714 -12700.589 0 56.44706 NA NA
12 275.00 59.463969 0.7590050 5.6263229 -14111.765 0 56.44706 NA NA
13 300.00 85.878675 0.7121356 5.9190016 -15522.942 0 56.44706 NA NA
14 325.00 120.510200 0.6543280 6.1872150 -16934.118 0 56.44706 NA NA
15 350.00 165.293399 0.5747058 6.4339077 -18345.295 0 56.44706 NA NA
Next: info
Up: CHNOSZ examples
Previous: thermo