Next: utilities
Up: CHNOSZ examples
Previous: chnosz
thermo> ## Don't show:
thermo> 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
thermo> ## End Don't show
thermo> ## exploring thermo$obigt
thermo> # what physical states there are
thermo> unique(thermo$obigt$state)
[1] "liq" "aq" "cr1" "cr2" "cr3" "cr" "cr4" "cr5" "cr6" "cr7" "cr8" "cr9" "gas"
thermo> # formulas of ten random species
thermo> n <- nrow(thermo$obigt)
thermo> thermo$obigt$formula[runif(10)*n]
[1] "Ni(CH3CH2CO2)2" "H2CO2" "Fe+3"
[4] "Na(Ca2Mg5)(AlSi7)O22(F)2" "C9H20S2" "CdC2O4"
[7] "K(AlSi3)O8" "C21H28N7O14P2-" "Eu+3"
[10] "H2N2O2"
thermo> ## cross-checking sources
thermo> # the reference sources
thermo> ref.source <- thermo$source$source
thermo> # only take those that aren't journal abbreviations
thermo> ref.source <- ref.source[-grep('_',ref.source)]
thermo> # sources of elemental data
thermo> element.source <- thermo$element$source
thermo> # primary sources in thermodynamic database
thermo> obigt.source1 <- thermo$obigt$source1
thermo> # secondary sources; some are NA
thermo> obigt.source2 <-
thermo+ thermo$obigt$source2[!is.na(thermo$obigt$source2)]
thermo> # sources of protein compositions
thermo> protein.source <- thermo$protein$source
thermo> # sources of stress response proteins
thermo> stress.source <- as.character(thermo$stress[2,])
thermo> # if the sources are all accounted for
thermo> # these all produce character(0)
thermo> element.source[!(element.source %in% ref.source)]
character(0)
thermo> obigt.source1[!(obigt.source1 %in% ref.source)]
[1] "PKR02"
thermo> obigt.source2[!(obigt.source2 %in% ref.source)]
character(0)
thermo> protein.source[!(protein.source %in% ref.source)]
character(0)
thermo> stress.source[!(stress.source %in% ref.source)]
character(0)
thermo> # determine if all the reference sources are cited
thermo> my.source <- c(element.source,obigt.source1,
thermo+ obigt.source2,protein.source,stress.source)
thermo> # this should produce character(0)
thermo> ref.source[!(ref.source %in% my.source)]
character(0)
thermo> ## make a table of duplicated species
thermo> name <- thermo$obigt$name
thermo> state <- thermo$obigt$state
thermo> source <- thermo$obigt$source1
thermo> species <- paste(name,state)
thermo> dups <- species[which(duplicated(species))]
thermo> id <- numeric()
thermo> for(i in 1:length(dups)) id <- c(id,which(species %in% dups[i]))
thermo> data.frame(name=name[id],state=state[id],source=source[id])
name state source
1 NH3 aq SSW01
2 NH3 aq SHS89
3 CO2 aq SSW01
4 CO2 aq SHS89
5 H2S aq SSW01
6 H2S aq SHS89
7 SO2 aq SSW01
8 SO2 aq SHS89
9 Ga+3 aq SSW+97
10 Ga+3 aq BDP+97
11 GaOH+2 aq SSW+97
12 GaOH+2 aq BDP+97
13 NaSO4- aq SLOP98
14 NaSO4- aq PSS95
15 As(OH)3 aq PPB+08
16 As(OH)3 aq PGS+96
17 H3AsO4 aq SSW+97
18 H3AsO4 aq PPB+08
19 Al+3 aq SSW+97
20 Al+3 aq TS01
21 AlOH+2 aq SSW+97
22 AlOH+2 aq TS01
23 ribose aq LH06a
24 ribose aq AP01
25 HCl aq SLOP98
26 HCl aq TZA97
27 CaOH+ aq SSW+97
28 CaOH+ aq AZ99
29 AgCl2- aq SSH97
30 AgCl2- aq AZ01
31 Ag(HS)2- aq SSH97
32 Ag(HS)2- aq AZ01
33 AgCl aq SSH97
34 AgCl aq AZ01
35 AgOH aq SSW+97
36 AgOH aq AZ01
37 Cu+ aq SSW+97
38 Cu+ aq AZ01
39 CuCl2- aq SSH97
40 CuCl2- aq AZ01
41 CuCl aq SSH97
42 CuCl aq AZ01
43 Au+ aq SSW+97
44 Au+ aq AZ01
45 AuCl2- aq SSH97
46 AuCl2- aq AZ01
47 Au(HS)2- aq SSH97
48 Au(HS)2- aq AZ01
49 AuCl aq SSH97
50 AuCl aq AZ01
51 B(OH)3 aq SHS89
52 B(OH)3 aq AVZP06
53 B(OH)3 aq PSS95
54 AgCl3-2 aq SSH97
55 AgCl3-2 aq AVZP06
56 B(OH)3 aq SHS89
57 B(OH)3 aq AVZP06
58 B(OH)3 aq PSS95
59 B(OH)4- aq AVZP06
60 B(OH)4- aq PSS95
61 NaB(OH)4 aq AVZP06
62 NaB(OH)4 aq PSS95
63 phenol aq DSM+97
64 phenol aq CSM+07
65 o-cresol aq DSM+97
66 o-cresol aq CSM+07
67 m-cresol aq DSM+97
68 m-cresol aq CSM+07
69 p-cresol aq DSM+97
70 p-cresol aq CSM+07
71 Eu(Ac)+2 aq SK93
72 Eu(Ac)+2 aq ZKO06
73 Eu(Ac)2+ aq SK93
74 Eu(Ac)2+ aq ZKO06
75 Eu(Ac)3 aq SK93
76 Eu(Ac)3 aq ZKO06
77 Y(Ac)+2 aq SK93
78 Y(Ac)+2 aq TZS+07
79 Y(Ac)2+ aq SK93
80 Y(Ac)2+ aq TZS+07
81 Y(Ac)3 aq SK93
82 Y(Ac)3 aq TZS+07
83 arsenopyrite cr PPB+08
84 arsenopyrite cr PKR02
85 pyridine liq LH06b
86 pyridine liq HOK+98
87 2-thiabutane liq Ric01
88 2-thiabutane liq HOK+98
89 thiophene liq Ric01
90 thiophene liq HOK+98
91 2-methylthiacyclopentane liq Ric01
92 2-methylthiacyclopentane liq HOK+98
93 methane gas WEP+82
94 methane gas HOK+98
95 2-thiabutane gas Ric01
96 2-thiabutane gas HOK+98
97 thiophene gas Ric01
98 thiophene gas HOK+98
thermo> ## accessing duplicated species
thermo> # using info()
thermo> i <- info("Al+3","aq") # length = 2
info: Al+3 available in aq, aq.
info: 990 refers to Al+3 aq (SSW+97, 07.Nov.97).
info: 1588 refers to Al+3 aq (TS01, 25.Aug.06).
thermo> i1 <- info("Al+3") # length = 1
info: Al+3 available in aq, aq.
info: 990 refers to Al+3 aq (SSW+97, 07.Nov.97).
thermo> stopifnot(i1==i[1])
thermo> thermo$opt$level <- 2
thermo> i2 <- info("Al+3") # length = 1
info: Al+3 available in aq, aq.
info: 1588 refers to Al+3 aq (TS01, 25.Aug.06).
thermo> stopifnot(i2==i[2])
thermo> # using subcrt()
thermo> subcrt("Al+3","aq") # always uses the first species
subcrt: 1 species at 15 values of T and P (wet)
$species
name formula state ispecies
990 Al+3 Al+3 aq 990
$out
$out$`Al+3`
T P logK G H S V Cp
1 0.01 1.000000 94.01696 -117511.55 -125825.3 -74.15060 -42.76480 -56.66752
2 25.00 1.000000 84.74217 -115609.00 -126834.0 -77.70000 -45.24283 -30.75917
3 50.00 1.000000 76.85294 -113637.56 -127523.5 -79.92314 -47.38180 -25.98128
4 75.00 1.000000 70.06455 -111614.86 -128182.5 -81.88682 -49.75352 -27.40196
5 100.00 1.013220 64.15671 -109542.53 -128920.0 -83.93053 -52.70333 -32.04680
6 125.00 2.320144 58.96212 -107418.01 -129804.7 -86.21867 -56.54294 -39.09308
7 150.00 4.757169 54.35060 -105233.96 -130884.8 -88.83875 -61.61639 -47.51521
8 175.00 8.918049 50.22111 -102983.34 -132173.4 -91.78033 -68.19544 -56.32462
9 200.00 15.536499 46.49419 -100659.51 -133748.5 -95.17128 -76.93731 -71.25484
10 225.00 25.478603 43.10403 -98250.61 -135866.9 -99.48324 -90.80821 -96.81480
11 250.00 39.736493 39.99422 -95737.21 -138788.2 -105.11750 -116.73264 -128.52847
12 275.00 59.431251 37.11994 -93103.07 -142325.3 -111.58158 -165.90997 -137.02743
13 300.00 85.837843 34.45831 -90369.01 -144889.0 -115.95789 -251.19155 -18.89501
14 325.00 120.457572 32.02703 -87656.49 -141022.0 -109.18120 -375.23788 731.41195
15 350.00 165.211289 29.93062 -85342.55 -92253.2 -30.25258 -391.83571 11149.62853
thermo> subcrt("Al+3") # pays attention to thermo$opt$level
subcrt: 1 species at 15 values of T and P (wet)
$species
name formula state ispecies
1588 Al+3 Al+3 aq 1588
$out
$out$`Al+3`
T P rho logK G H S V Cp
1 0.01 1.000000 0.9998289 94.80932 -118501.92 -127793.95 -77.77010 -43.71978 -56.69448
2 25.00 1.000000 0.9970614 85.40261 -116510.00 -128769.02 -81.20300 -45.33933 -28.52400
3 50.00 1.000000 0.9880295 77.40481 -114453.58 -129391.90 -83.21208 -47.06624 -22.99876
4 75.00 1.000000 0.9748643 70.52702 -112351.59 -129971.56 -84.93934 -49.19221 -24.07131
5 100.00 1.013220 0.9583926 64.54525 -110205.93 -130623.04 -86.74463 -51.97450 -28.51414
6 125.00 2.320144 0.9390726 59.28925 -108013.97 -131417.70 -88.79917 -55.68851 -35.42275
7 150.00 4.757169 0.9170577 54.62652 -105768.20 -132404.76 -91.19275 -60.65963 -43.74565
8 175.00 8.918049 0.8923427 50.45423 -103461.37 -133598.29 -93.91626 -67.15044 -52.48965
9 200.00 15.536499 0.8647434 46.69146 -101086.60 -135076.63 -97.09750 -75.81165 -67.35135
10 225.00 25.478603 0.8338733 43.27127 -98631.81 -137095.88 -101.20591 -89.59236 -92.82547
11 250.00 39.736493 0.7990719 40.13632 -96077.36 -139915.51 -106.64208 -115.38094 -124.51120
12 275.00 59.431251 0.7592362 37.24101 -93406.71 -143353.00 -112.92153 -164.32908 -133.36597
13 300.00 85.837843 0.7124075 34.56165 -90640.02 -145837.96 -117.15885 -249.28186 -17.11843
14 325.00 120.457572 0.6545772 32.11490 -87896.97 -141976.95 -110.39309 -373.19024 724.02682
15 350.00 165.211289 0.5746875 30.00321 -85549.53 -93726.76 -32.30182 -393.04091 11035.24949
thermo> # .. showing the energetic differences between the
thermo> # duplicated species
thermo> subcrt(i,c(-1,1))
subcrt: 2 species at 15 values of T and P (wet)
$reaction
coeff name formula state ispecies
990 -1 Al+3 Al+3 aq 990
1588 1 Al+3 Al+3 aq 1588
$out
T P rho logK G H S V Cp
1 0.01 1.000000 0.9998289 0.79236323 -990.3728 -1968.6759 -3.619496 -0.9549867 -0.02695691
2 25.00 1.000000 0.9970614 0.66043900 -901.0000 -1935.0210 -3.503000 -0.0965015 2.23516866
3 50.00 1.000000 0.9880295 0.55186869 -816.0132 -1868.4367 -3.288943 0.3155545 2.98251778
4 75.00 1.000000 0.9748643 0.46246997 -736.7280 -1789.0641 -3.052519 0.5613137 3.33065603
5 100.00 1.013220 0.9583926 0.38853939 -663.4005 -1703.0848 -2.814104 0.7288277 3.53265968
6 125.00 2.320144 0.9390726 0.32712235 -595.9560 -1612.9825 -2.580497 0.8544348 3.67032976
7 150.00 4.757169 0.9170577 0.27592000 -534.2380 -1519.9335 -2.353997 0.9567604 3.76956808
8 175.00 8.918049 0.8923427 0.23311847 -478.0324 -1424.8535 -2.135936 1.0450023 3.83497092
9 200.00 15.536499 0.8647434 0.19727437 -427.0972 -1328.0895 -1.926219 1.1256646 3.90348204
10 225.00 25.478603 0.8338733 0.16724057 -381.2054 -1228.9522 -1.722664 1.2158497 3.98933327
11 250.00 39.736493 0.7990719 0.14209807 -340.1509 -1127.3343 -1.524576 1.3516984 4.01727180
12 275.00 59.431251 0.7592362 0.12106263 -303.6455 -1027.7426 -1.339954 1.5808968 3.66145551
13 300.00 85.837843 0.7124075 0.10333746 -271.0088 -948.9408 -1.200960 1.9096916 1.77657539
14 325.00 120.457572 0.6545772 0.08786574 -240.4844 -954.9760 -1.211887 2.0476400 -7.38513427
15 350.00 165.211289 0.5746875 0.07259233 -206.9858 -1473.5697 -2.049237 -1.2052078 -114.37904456
thermo> # using basis()
thermo> basis(c("Al+3","H+","H2"))
Al H Z ispecies logact state
Al+3 1 0 3 1588 0 aq
H+ 0 1 1 3 0 aq
H2 0 2 0 66 0 aq
thermo> thermo$opt$level <- 1
thermo> basis(c("Al+3","H+","H2"))
basis: changed basis to Al+3 H+ H2.
basis: elements unchanged. maintaining elemental chemical potentials.
Al H Z ispecies logact state
Al+3 1 0 3 990 -0.660439 aq
H+ 0 1 1 3 0.000000 aq
H2 0 2 0 66 0.000000 aq
thermo> # using species()
thermo> species("Al+3","aq") # listens to thermo$opt$level
thermo> thermo$opt$level <- 2
thermo> species("Al+3") # also listens to thermo$opt$level
thermo> species(delete=TRUE)
thermo> species(rev(i)) # can also use the species indices
Next: utilities
Up: CHNOSZ examples
Previous: chnosz