next up previous
Next: util.affinity Up: CHNOSZ examples Previous: species

[3] affinity

## Don't show: data(thermo)
thermo$obigt: 1809 aqueous, 3368 total species
## End(Don't show)## set up a system and calculate ## chemical affinities of formation reactions basis(c("SiO2","MgO","H2O","O2"),c(-5,-5,0,999))
H Mg O Si ispecies logact state SiO2 0 0 2 1 72 -5 aq MgO 0 1 1 0 2002 -5 cr H2O 2 0 1 0 1 0 liq O2 0 0 2 0 67 999 aq
species(c("quartz","enstatite","forsterite"))
SiO2 MgO H2O O2 ispecies logact state name 1 1 0 0 0 2014 0 cr1 quartz 2 1 1 0 0 1906 0 cr1 enstatite 3 1 2 0 0 1929 0 cr forsterite
# chemical affinities (A/2.303RT) at 25 deg C and 1 bar affinity()
energy.args: temperature is 25 C energy.args: pressure is Psat subcrt: 7 species at 298.15 K and 1 bar (wet) $sout $sout$SiO2 logK G 1 146.0076 -199190 $sout$periclase logK G 1 99.75194 -136086 $sout$water logK G 1 41.55247 -56687.71 $sout$O2 logK G 1 -2.898308 3954 $sout$quartz logK G 1 150.0069 -204646 $sout$enstatite logK G 1 255.768 -348930 $sout$forsterite logK G 1 360.3197 -491564 $property [1] "A" $basis H Mg O Si ispecies logact state SiO2 0 0 2 1 72 -5 aq MgO 0 1 1 0 2002 -5 cr H2O 2 0 1 0 1 0 liq O2 0 0 2 0 67 999 aq $species SiO2 MgO H2O O2 ispecies logact state name 1 1 0 0 0 2014 0 cr1 quartz 2 1 1 0 0 1906 0 cr1 enstatite 3 1 2 0 0 1929 0 cr forsterite $T [1] 298.15 $P [1] "Psat" $vars character(0) $vals $vals[[1]] [1] NA $values $values$`2014` [1] -1.000716 $values$`1906` [1] 0.008472902 $values$`1929` [1] -0.1917995
# at higher temperature and pressure affinity(T=500,P=2000)
energy.args: temperature is 500 C energy.args: pressure is 2000 bar subcrt: 7 species at 773.15 K and 2000 bar (wet) $sout $sout$SiO2 logK G 1 58.67394 -207570.8 $sout$periclase logK G 1 39.92007 -141225.2 $sout$water logK G 1 19.25955 -68134.5 $sout$O2 logK G 1 5.245429 -18556.76 $sout$quartz logK G 1 59.83361 -211673.4 $sout$enstatite logK G 1 102.0606 -361059.8 $sout$forsterite logK G 1 143.8758 -508989.4 $property [1] "A" $basis H Mg O Si ispecies logact state SiO2 0 0 2 1 72 -5 aq MgO 0 1 1 0 2002 -5 cr H2O 2 0 1 0 1 0 liq O2 0 0 2 0 67 999 aq $species SiO2 MgO H2O O2 ispecies logact state name 1 1 0 0 0 2014 0 cr1 quartz 2 1 1 0 0 1906 0 cr1 enstatite 3 1 2 0 0 1929 0 cr forsterite $T [1] 773.15 $P [1] 2000 $vars character(0) $vals $vals[[1]] [1] NA $values $values$`2014` [1] -3.840332 $values$`1906` [1] -6.533419 $values$`1929` [1] -9.638279
# ten different temperatures at one pressure affinity(T=c(500,1000,10),P=2000)
energy.args: pressure is 2000 bar energy.args: variable 1 is T at 10 values from 773.15 to 1273.15 K subcrt: 7 species at 10 values of T and P (wet) subcrt: some points above transition temperature for quartz cr1 (using NA for G) subcrt: some points above transition temperature for enstatite cr1 (using NA for G) $sout $sout$SiO2 T P logK G 1 773.1500 2000 58.67394 -207570.8 2 828.7056 2000 55.13264 -209057.8 3 884.2611 2000 52.05246 -210610.0 4 939.8167 2000 49.35128 -212226.1 5 995.3722 2000 46.96669 -213910.8 6 1050.9278 2000 44.85023 -215672.5 7 1106.4833 2000 42.96315 -217519.5 8 1162.0389 2000 41.27363 -219457.6 9 1217.5944 2000 39.75507 -221489.2 10 1273.1500 2000 38.38505 -223614.0 $sout$periclase T P logK G 1 773.1500 2000 39.92007 -141225.2 2 828.7056 2000 37.48789 -142150.5 3 884.2611 2000 35.37183 -143118.3 4 939.8167 2000 33.51531 -144126.4 5 995.3722 2000 31.87442 -145172.8 6 1050.9278 2000 30.41462 -146255.6 7 1106.4833 2000 29.10829 -147373.3 8 1162.0389 2000 27.93316 -148524.5 9 1217.5944 2000 26.85723 -149707.7 10 1273.1500 2000 25.90691 -150922.0 $sout$water T P logK G 1 773.1500 2000 19.25955 -68134.50 2 828.7056 2000 18.46115 -70002.94 3 884.2611 2000 17.78034 -71941.22 4 939.8167 2000 17.19508 -73944.30 5 995.3722 2000 16.68824 -76006.96 6 1050.9278 2000 16.24630 -78124.02 7 1106.4833 2000 15.85852 -80290.63 8 1162.0389 2000 15.51633 -82502.46 9 1217.5944 2000 15.21281 -84755.78 10 1273.1500 2000 14.94236 -87047.43 $sout$O2 T P logK G 1 773.1500 2000 5.245429 -18556.76 2 828.7056 2000 5.958308 -22593.34 3 884.2611 2000 6.642488 -26876.25 4 939.8167 2000 7.298209 -31384.61 5 995.3722 2000 7.922007 -36080.95 6 1050.9278 2000 8.509278 -40918.80 7 1106.4833 2000 9.056779 -45853.86 8 1162.0389 2000 9.563717 -50851.60 9 1217.5944 2000 10.031490 -55888.87 10 1273.1500 2000 10.462907 -60952.17 $sout$quartz T P logK G 1 773.1500 2000 59.83361 -211673.4 2 828.7056 2000 56.17083 -212994.5 3 884.2611 2000 52.98451 -214381.2 4 939.8167 2000 NA NA 5 995.3722 2000 NA NA 6 1050.9278 2000 NA NA 7 1106.4833 2000 NA NA 8 1162.0389 2000 NA NA 9 1217.5944 2000 NA NA 10 1273.1500 2000 NA NA $sout$enstatite T P logK G 1 773.1500 2000 102.06059 -361059.8 2 828.7056 2000 95.80160 -363270.6 3 884.2611 2000 90.35447 -365584.2 4 939.8167 2000 NA NA 5 995.3722 2000 NA NA 6 1050.9278 2000 NA NA 7 1106.4833 2000 NA NA 8 1162.0389 2000 NA NA 9 1217.5944 2000 NA NA 10 1273.1500 2000 NA NA $sout$forsterite T P logK G 1 773.1500 2000 143.87580 -508989.4 2 828.7056 2000 135.06988 -512172.2 3 884.2611 2000 127.40729 -515504.0 4 939.8167 2000 120.68351 -518977.3 5 995.3722 2000 114.73986 -522585.2 6 1050.9278 2000 109.45135 -526321.7 7 1106.4833 2000 104.71823 -530181.4 8 1162.0389 2000 100.45991 -534159.2 9 1217.5944 2000 96.61057 -538250.7 10 1273.1500 2000 93.11598 -542451.7 $property [1] "A" $basis H Mg O Si ispecies logact state SiO2 0 0 2 1 72 -5 aq MgO 0 1 1 0 2002 -5 cr H2O 2 0 1 0 1 0 liq O2 0 0 2 0 67 999 aq $species SiO2 MgO H2O O2 ispecies logact state name 1 1 0 0 0 2014 0 cr1 quartz 2 1 1 0 0 1906 0 cr1 enstatite 3 1 2 0 0 1929 0 cr forsterite $T numeric(0) $P [1] 2000 $vars [1] "T" $vals $vals[[1]] [1] 500.0000 555.5556 611.1111 666.6667 722.2222 777.7778 833.3333 888.8889 944.4444 1000.0000 $values $values$`2014` [1] -3.840332 -3.961812 -4.067943 NA NA NA NA NA NA NA $values$`1906` [1] -6.533419 -6.818924 -7.069816 NA NA NA NA NA NA NA $values$`1929` [1] -9.638279 -10.038539 -10.388830 -10.698385 -10.975681 -11.228104 -11.461496 -11.680047 -11.886569 [10] -12.082887
# at 25 temperatures and pressures affinity(T=c(500,1000,5),P=c(1000,5000,5))
energy.args: variable 1 is T at 5 values from 773.15 to 1273.15 K energy.args: variable 2 is P at 5 values from 1000 to 5000 bar subcrt: 7 species at 25 values of T and P (wet) subcrt: some points above transition temperature for quartz cr1 (using NA for G) subcrt: some points above transition temperature for enstatite cr1 (using NA for G) $sout $sout$SiO2 T P logK G 1 773.15 1000 58.61638 -207367.2 2 773.15 2000 58.67394 -207570.8 3 773.15 3000 58.64421 -207465.6 4 773.15 4000 58.59738 -207300.0 5 773.15 5000 58.54380 -207110.4 6 898.15 1000 51.07367 -209895.5 7 898.15 2000 51.34436 -211008.0 8 898.15 3000 51.36294 -29034.3 9 898.15 4000 51.34149 -210996.2 10 898.15 5000 51.30634 -210851.7 11 1023.15 1000 45.47259 -212885.6 12 1023.15 2000 45.87755 -214781.5 13 1023.15 3000 45.94656 -215104.6 14 1023.15 4000 45.94863 -215114.2 15 1023.15 5000 45.92870 -215020.9 16 1148.15 1000 41.24339 -216675.7 17 1148.15 2000 41.67903 -218964.4 18 1148.15 3000 41.78293 -219510.2 19 1148.15 4000 41.80251 -219613.1 20 1148.15 5000 41.79311 -219563.7 21 1273.15 1000 37.96343 -221157.8 22 1273.15 2000 38.38505 -223614.0 23 1273.15 3000 38.50388 -224306.2 24 1273.15 4000 38.53209 -224470.6 25 1273.15 5000 38.52669 -224439.1 $sout$periclase T P logK G 1 773.15 1000 39.99606 -141494.1 2 773.15 2000 39.92007 -141225.2 3 773.15 3000 39.84408 -140956.4 4 773.15 4000 39.76809 -140687.6 5 773.15 5000 39.69210 -140418.7 6 898.15 1000 34.95068 -143635.5 7 898.15 2000 34.88527 -143366.7 8 898.15 3000 34.81985 -143097.8 9 898.15 4000 34.75444 -142829.0 10 898.15 5000 34.68902 -142560.2 11 1023.15 1000 31.18117 -145978.6 12 1023.15 2000 31.12375 -145709.7 13 1023.15 3000 31.06633 -145440.9 14 1023.15 4000 31.00891 -145172.1 15 1023.15 5000 30.95148 -144903.2 16 1148.15 1000 28.26687 -148502.4 17 1148.15 2000 28.21570 -148233.6 18 1148.15 3000 28.16453 -147964.8 19 1148.15 4000 28.11335 -147695.9 20 1148.15 5000 28.06218 -147427.1 21 1273.15 1000 25.95306 -151190.8 22 1273.15 2000 25.90691 -150922.0 23 1273.15 3000 25.86076 -150653.1 24 1273.15 4000 25.81461 -150384.3 25 1273.15 5000 25.76847 -150115.5 $sout$water T P logK G 1 773.15 1000 19.45590 -68829.13 2 773.15 2000 19.25955 -68134.50 3 773.15 3000 19.09379 -67548.09 4 773.15 4000 18.94215 -67011.63 5 773.15 5000 18.79959 -66507.32 6 898.15 1000 17.85178 -73364.78 7 898.15 2000 17.62581 -72436.12 8 898.15 3000 17.45754 -71744.59 9 898.15 4000 17.31002 -71138.35 10 898.15 5000 17.17446 -70581.21 11 1023.15 1000 16.71777 -78266.33 12 1023.15 2000 16.45989 -77059.00 13 1023.15 3000 16.28488 -76239.70 14 1023.15 4000 16.13803 -75552.21 15 1023.15 5000 16.00624 -74935.18 16 1148.15 1000 15.87914 -83422.44 17 1148.15 2000 15.59801 -81945.48 18 1148.15 3000 15.41570 -80987.70 19 1148.15 4000 15.26800 -80211.76 20 1148.15 5000 15.13826 -79530.14 21 1273.15 1000 15.23897 -88775.36 22 1273.15 2000 14.94236 -87047.43 23 1273.15 3000 14.75404 -85950.35 24 1273.15 4000 14.60513 -85082.91 25 1273.15 5000 14.47661 -84334.20 $sout$O2 T P logK G 1 773.15 1000 5.847896 -20688.10 2 773.15 2000 5.245429 -18556.76 3 773.15 3000 4.922519 -17414.40 4 773.15 4000 4.660520 -16487.52 5 773.15 5000 4.425326 -15655.48 6 898.15 1000 8.000166 -32877.98 7 898.15 2000 6.809175 -27983.41 8 898.15 3000 6.399786 -26300.96 9 898.15 4000 6.120452 -25152.99 10 898.15 5000 5.888545 -24199.93 11 1023.15 1000 9.775211 -45763.87 12 1023.15 2000 8.220464 -38485.12 13 1023.15 3000 7.702246 -36059.02 14 1023.15 4000 7.395508 -34622.99 15 1023.15 5000 7.160975 -33525.00 16 1148.15 1000 11.053381 -58069.90 17 1148.15 2000 9.440739 -49597.74 18 1148.15 3000 8.850810 -46498.49 19 1148.15 4000 8.524758 -44785.56 20 1148.15 5000 8.291816 -43561.77 21 1273.15 1000 12.004199 -69931.04 22 1273.15 2000 10.462907 -60952.17 23 1273.15 3000 9.855069 -57411.19 24 1273.15 4000 9.529847 -55516.59 25 1273.15 5000 9.311444 -54244.27 $sout$quartz T P logK G 1 773.15 1000 59.98689 -212215.6 2 773.15 2000 59.83361 -211673.4 3 773.15 3000 59.68033 -211131.1 4 773.15 4000 59.52705 -210588.8 5 773.15 5000 59.37377 -210046.6 6 898.15 1000 NA NA 7 898.15 2000 52.25196 -214737.9 8 898.15 3000 52.12001 -214195.6 9 898.15 4000 51.98806 -213653.4 10 898.15 5000 51.85612 -213111.1 11 1023.15 1000 NA NA 12 1023.15 2000 NA NA 13 1023.15 3000 NA NA 14 1023.15 4000 NA NA 15 1023.15 5000 NA NA 16 1148.15 1000 NA NA 17 1148.15 2000 NA NA 18 1148.15 3000 NA NA 19 1148.15 4000 NA NA 20 1148.15 5000 NA NA 21 1273.15 1000 NA NA 22 1273.15 2000 NA NA 23 1273.15 3000 NA NA 24 1273.15 4000 NA NA 25 1273.15 5000 NA NA $sout$enstatite T P logK G 1 773.15 1000 102.27189 -361807.3 2 773.15 2000 102.06059 -361059.8 3 773.15 3000 101.84929 -360312.3 4 773.15 4000 101.63799 -359564.7 5 773.15 5000 101.42669 -358817.2 6 898.15 1000 89.28363 -366925.6 7 898.15 2000 89.10174 -366178.0 8 898.15 3000 88.91984 -365430.5 9 898.15 4000 88.73795 -364683.0 10 898.15 5000 88.55606 -363935.5 11 1023.15 1000 NA NA 12 1023.15 2000 NA NA 13 1023.15 3000 NA NA 14 1023.15 4000 NA NA 15 1023.15 5000 NA NA 16 1148.15 1000 NA NA 17 1148.15 2000 NA NA 18 1148.15 3000 NA NA 19 1148.15 4000 NA NA 20 1148.15 5000 NA NA 21 1273.15 1000 NA NA 22 1273.15 2000 NA NA 23 1273.15 3000 NA NA 24 1273.15 4000 NA NA 25 1273.15 5000 NA NA $sout$forsterite T P logK G 1 773.15 1000 144.17164 -510036.0 2 773.15 2000 143.87580 -508989.4 3 773.15 3000 143.57995 -507942.8 4 773.15 4000 143.28411 -506896.2 5 773.15 5000 142.98827 -505849.6 6 898.15 1000 125.89987 -516026.0 7 898.15 2000 125.64520 -516359.3 8 898.15 3000 125.39053 -515312.7 9 898.15 4000 125.13586 -514266.1 10 898.15 5000 124.88119 -513219.5 11 1023.15 1000 112.24401 -525484.3 12 1023.15 2000 112.02046 -524437.7 13 1023.15 3000 111.79690 -523391.1 14 1023.15 4000 111.57335 -522344.5 15 1023.15 5000 111.34979 -521297.9 16 1148.15 1000 101.68301 -534200.5 17 1148.15 2000 101.48379 -533153.9 18 1148.15 3000 101.28458 -532107.3 19 1148.15 4000 101.08536 -531060.7 20 1148.15 5000 100.88614 -530014.1 21 1273.15 1000 93.29563 -543498.3 22 1273.15 2000 93.11598 -542451.7 23 1273.15 3000 92.93632 -541405.1 24 1273.15 4000 92.75666 -540358.5 25 1273.15 5000 92.57700 -539311.8 $property [1] "A" $basis H Mg O Si ispecies logact state SiO2 0 0 2 1 72 -5 aq MgO 0 1 1 0 2002 -5 cr H2O 2 0 1 0 1 0 liq O2 0 0 2 0 67 999 aq $species SiO2 MgO H2O O2 ispecies logact state name 1 1 0 0 0 2014 0 cr1 quartz 2 1 1 0 0 1906 0 cr1 enstatite 3 1 2 0 0 1929 0 cr forsterite $T numeric(0) $P numeric(0) $vars [1] "T" "P" $vals $vals[[1]] [1] 500 625 750 875 1000 $vals[[2]] [1] 1000 2000 3000 4000 5000 $values $values$`2014` [,1] [,2] [,3] [,4] [,5] [1,] -3.629498 -3.840332 -3.963881 -4.070328 -4.170036 [2,] NA -4.092402 -4.242933 -4.353429 -4.450220 [3,] NA NA NA NA NA [4,] NA NA NA NA NA [5,] NA NA NA NA NA $values$`1906` [,1] [,2] [,3] [,4] [,5] [1,] -6.340555 -6.533419 -6.638997 -6.727473 -6.809210 [2,] -6.602723 -7.127890 -7.262952 -7.357977 -7.439299 [3,] NA NA NA NA NA [4,] NA NA NA NA NA [5,] NA NA NA NA NA $values$`1929` [,1] [,2] [,3] [,4] [,5] [1,] -9.436862 -9.638279 -9.75241 -9.849439 -9.939729 [2,] -10.075164 -10.469693 -10.61212 -10.714505 -10.803190 [3,] -10.590926 -11.104597 -11.28232 -11.393097 -11.481877 [4,] -11.094114 -11.626627 -11.82602 -11.943864 -12.031333 [5,] -11.573909 -12.082887 -12.28908 -12.404656 -12.486619
# as a function of logarithm of activity of MgO affinity(MgO=c(-10,-5,10))
energy.args: temperature is 25 C energy.args: pressure is Psat energy.args: variable 1 is log_a(MgO) at 10 values from -10 to -5 subcrt: 7 species at 298.15 K and 1 bar (wet) $sout $sout$SiO2 logK G 1 146.0076 -199190 $sout$periclase logK G 1 99.75194 -136086 $sout$water logK G 1 41.55247 -56687.71 $sout$O2 logK G 1 -2.898308 3954 $sout$quartz logK G 1 150.0069 -204646 $sout$enstatite logK G 1 255.768 -348930 $sout$forsterite logK G 1 360.3197 -491564 $property [1] "A" $basis H Mg O Si ispecies logact state SiO2 0 0 2 1 72 -5 aq MgO 0 1 1 0 2002 -5 cr H2O 2 0 1 0 1 0 liq O2 0 0 2 0 67 999 aq $species SiO2 MgO H2O O2 ispecies logact state name 1 1 0 0 0 2014 0 cr1 quartz 2 1 1 0 0 1906 0 cr1 enstatite 3 1 2 0 0 1929 0 cr forsterite $T [1] 298.15 $P [1] "Psat" $vars [1] "MgO" $vals $vals[[1]] [1] -10.000000 -9.444444 -8.888889 -8.333333 -7.777778 -7.222222 -6.666667 -6.111111 -5.555556 [10] -5.000000 $values $values$`2014` [1] -1.000716 -1.000716 -1.000716 -1.000716 -1.000716 -1.000716 -1.000716 -1.000716 -1.000716 -1.000716 $values$`1906` [1] -4.991527098 -4.435971543 -3.880415987 -3.324860432 -2.769304876 -2.213749321 -1.658193765 -1.102638210 [9] -0.547082654 0.008472902 $values$`1929` [1] -10.1917995 -9.0806884 -7.9695773 -6.8584662 -5.7472861 -4.6362439 -3.5251328 -2.4140217 [9] -1.3029106 -0.1917995
## equilibrium constants of formation reactions affinity(property="logK")
energy.args: temperature is 25 C energy.args: pressure is Psat subcrt: 7 species at 298.15 K and 1 bar (wet) $sout $sout$SiO2 logK G 1 146.0076 -199190 $sout$periclase logK G 1 99.75194 -136086 $sout$water logK G 1 41.55247 -56687.71 $sout$O2 logK G 1 -2.898308 3954 $sout$quartz logK G 1 150.0069 -204646 $sout$enstatite logK G 1 255.768 -348930 $sout$forsterite logK G 1 360.3197 -491564 $property [1] "logK" $basis H Mg O Si ispecies logact state SiO2 0 0 2 1 72 -5 aq MgO 0 1 1 0 2002 -5 cr H2O 2 0 1 0 1 0 liq O2 0 0 2 0 67 999 aq $species SiO2 MgO H2O O2 ispecies logact state name 1 1 0 0 0 2014 0 cr1 quartz 2 1 1 0 0 1906 0 cr1 enstatite 3 1 2 0 0 1929 0 cr forsterite $T [1] 298.15 $P [1] "Psat" $vars character(0) $vals $vals[[1]] [1] NA $values $values$`2014` [1] 3.999284 $values$`1906` [1] 10.00847 $values$`1929` [1] 14.8082
# Standard molal Gibbs energies of species, # user units (default: cal/mol) affinity(property="G.species")
energy.args: temperature is 25 C energy.args: pressure is Psat subcrt: 7 species at 298.15 K and 1 bar (wet) $sout $sout$SiO2 G 1 -199190 $sout$periclase G 1 -136086 $sout$water G 1 -56687.71 $sout$O2 G 1 3954 $sout$quartz G 1 -204646 $sout$enstatite G 1 -348930 $sout$forsterite G 1 -491564 $property [1] "G.species" $basis H Mg O Si ispecies logact state SiO2 0 0 2 1 72 -5 aq MgO 0 1 1 0 2002 -5 cr H2O 2 0 1 0 1 0 liq O2 0 0 2 0 67 999 aq $species SiO2 MgO H2O O2 ispecies logact state name 1 1 0 0 0 2014 0 cr1 quartz 2 1 1 0 0 1906 0 cr1 enstatite 3 1 2 0 0 1929 0 cr forsterite $T [1] 298.15 $P [1] "Psat" $vars character(0) $vals $vals[[1]] [1] NA $values $values$`2014` [1] -204646 $values$`1906` [1] -348930 $values$`1929` [1] -491564
# Standard molal Gibbs energies of reactions affinity(property="G")
energy.args: temperature is 25 C energy.args: pressure is Psat subcrt: 7 species at 298.15 K and 1 bar (wet) $sout $sout$SiO2 G 1 -199190 $sout$periclase G 1 -136086 $sout$water G 1 -56687.71 $sout$O2 G 1 3954 $sout$quartz G 1 -204646 $sout$enstatite G 1 -348930 $sout$forsterite G 1 -491564 $property [1] "G" $basis H Mg O Si ispecies logact state SiO2 0 0 2 1 72 -5 aq MgO 0 1 1 0 2002 -5 cr H2O 2 0 1 0 1 0 liq O2 0 0 2 0 67 999 aq $species SiO2 MgO H2O O2 ispecies logact state name 1 1 0 0 0 2014 0 cr1 quartz 2 1 1 0 0 1906 0 cr1 enstatite 3 1 2 0 0 1929 0 cr forsterite $T [1] 298.15 $P [1] "Psat" $vars character(0) $vals $vals[[1]] [1] NA $values $values$`2014` [1] -5456 $values$`1906` [1] -13654 $values$`1929` [1] -20202
## affinities of metabolic reactions ## after Amend and Shock, 2001, Fig. 7 # use aq state for all basis species (including O2) basis(c("CO2", "H2", "NH3", "O2", "H2S", "H+"), "aq")
C H N O S Z ispecies logact state CO2 1 0 0 2 0 0 69 0 aq H2 0 2 0 0 0 0 66 0 aq NH3 0 3 1 0 0 0 68 0 aq O2 0 0 0 2 0 0 67 0 aq H2S 0 2 0 0 1 0 70 0 aq H+ 0 1 0 0 0 1 3 0 aq
# we're going to make H2O species("H2O")
CO2 H2 NH3 O2 H2S H+ ispecies logact state name 1 0 1 0 0.5 0 0 1 0 liq water
# a function to create the plots doplot <- function(T) { res <- 20 # calculate affinity/2.303RT as a function of loga(H2) and loga(O2) a <- affinity(H2=c(-10, 0, res), O2=c(-10, 0, res), T=T) T.K <- convert(T, "K") # temperature in Kelvin acal <- convert(a$values[[1]], "G", T.K) # affinity (cal/mol) akJ <- convert(acal, "J")/1000 # affinity (kJ/mol) # now contour the values xyvals <- seq(-10, 0, length.out=res) contour(x=xyvals, y=xyvals, z=t(akJ), levels=seq(-150, -250, -20), labcex=1, xlab=axis.label("H2"), ylab=axis.label("O2")) # show the temperature legend("topleft", bg="white", cex=1, legend=describe.property("T", T, digits=0, ret.val=TRUE) ) } # plot layout with space for title at top layout(matrix(c(1, 1, 2, 3, 4, 5), ncol=2, byrow=TRUE), heights=c(1, 4, 4)) par(mar=c(0, 0, 0, 0)) plot.new() # we use subcrt() to generate a reaction for titling the plot rxnexpr <- describe.reaction(subcrt("H2O", 1)$reaction, states="all")
info.character: found H2O(liq), also available in gas subcrt: 1 species at 15 values of T and P (wet) subcrt: reaction is not balanced; it is missing this composition: H O -2 -1 subcrt: adding missing composition from basis definition and restarting... subcrt: 3 species at 15 values of T and P (wet)
# also in the title is the property with its units E.units("J")
[1] "J"
Gexpr <- axis.label("DGr", prefix="k")[[2]] text(0.5, 0.6, substitute(paste(G~~"for"~~r), list(G=Gexpr, r=rxnexpr)), cex=2) text(0.5, 0.2, "after Amend and Shock, 2001 Figure 7", cex=2) # now make the plots par(mar=c(3, 3, 0.5, 0.5), cex=1.3, mgp=c(2, 1, 0)) sapply(c(25, 55, 100, 150), doplot)
energy.args: temperature is 25 C energy.args: pressure is Psat energy.args: variable 1 is log_a(H2) at 20 values from -10 to 0 energy.args: variable 2 is log_a(O2) at 20 values from -10 to 0 subcrt: 7 species at 298.15 K and 1 bar (wet) energy.args: temperature is 55 C energy.args: pressure is Psat energy.args: variable 1 is log_a(H2) at 20 values from -10 to 0 energy.args: variable 2 is log_a(O2) at 20 values from -10 to 0 subcrt: 7 species at 328.15 K and 1 bar (wet) energy.args: temperature is 100 C energy.args: pressure is Psat energy.args: variable 1 is log_a(H2) at 20 values from -10 to 0 energy.args: variable 2 is log_a(O2) at 20 values from -10 to 0 subcrt: 7 species at 373.15 K and 1.01 bar (wet) energy.args: temperature is 150 C energy.args: pressure is Psat energy.args: variable 1 is log_a(H2) at 20 values from -10 to 0 energy.args: variable 2 is log_a(O2) at 20 values from -10 to 0 subcrt: 7 species at 423.15 K and 4.76 bar (wet) [,1] [,2] [,3] [,4] rect List,4 List,4 List,4 List,4 text List,2 List,2 List,2 List,2

Image affinity1

 

# affinity() can handle the three dimensions simultaneously print(affinity(H2=c(-10, 0, 3), O2=c(-10, 0, 3), T=c(25, 150, 4))$values)
energy.args: pressure is Psat energy.args: variable 1 is log_a(H2) at 3 values from -10 to 0 energy.args: variable 2 is log_a(O2) at 3 values from -10 to 0 energy.args: variable 3 is T at 4 values from 298.15 to 423.15 K subcrt: 7 species at 4 values of T and P (wet) $`1` , , 1 [,1] [,2] [,3] [1,] 31.10664 33.60664 36.10664 [2,] 36.10664 38.60664 41.10664 [3,] 41.10664 43.60664 46.10664 , , 2 [,1] [,2] [,3] [1,] 25.10194 27.60194 30.10194 [2,] 30.10194 32.60194 35.10194 [3,] 35.10194 37.60194 40.10194 , , 3 [,1] [,2] [,3] [1,] 20.29182 22.79182 25.29182 [2,] 25.29182 27.79182 30.29182 [3,] 30.29182 32.79182 35.29182 , , 4 [,1] [,2] [,3] [1,] 16.34186 18.84186 21.34186 [2,] 21.34186 23.84186 26.34186 [3,] 26.34186 28.84186 31.34186
# this is so the plots in the next examples show up OK E.units("cal")
[1] "cal"
layout(matrix(1)) par(mar=c(5.1, 4.1, 4.1, 2.1)) ## amino acid synthesis at low and high temperatures ## after Amend and Shock, 1998 # select the basis species and species of interest # and set their activities, first for the 18 degree C case basis(c("H2O", "CO2", "NH4+", "H2", "H+", "H2S"), log10(c(1, 1e-4, 5e-8, 2e-9, 5e-9, 1e-15)))
C H N O S Z ispecies logact state H2O 0 2 0 1 0 0 1 0.00000 liq CO2 1 0 0 2 0 0 69 -4.00000 aq NH4+ 0 4 1 0 0 1 18 -7.30103 aq H2 0 2 0 0 0 0 66 -8.69897 aq H+ 0 1 0 0 0 1 3 -8.30103 aq H2S 0 2 0 0 1 0 70 -15.00000 aq
species(sort(aminoacids("Z")), log10(c(3.9, 0.7, 1.1, 3.3, 0.5, 3.8, 1.0, 5.8, 1.2, 0.7, 0.8, 1.0, 2.8, 0.5, 0.5, 4.6, 5.8, 0.6, 0.9, 2.8)/1e9))
H2O CO2 NH4+ H2 H+ H2S ispecies logact state name 1 -4 3 1 6 -1 0 1504 -8.408935 aq alanine 2 -10 6 4 11 -3 0 1506 -9.154902 aq argininium 3 -5 4 2 6 -2 0 1508 -8.958607 aq asparagine 4 -4 4 1 6 -2 0 1510 -8.481486 aq aspartate 5 -4 3 1 5 -2 1 1512 -9.301030 aq cysteinate 6 -6 5 1 9 -2 0 1515 -8.420216 aq glutamate 7 -7 5 2 9 -2 0 1513 -9.000000 aq glutamine 8 -2 2 1 3 -1 0 1516 -8.236572 aq glycine 9 -10 6 3 10 -2 0 1519 -8.920819 aq histidinium 10 -10 6 1 15 -1 0 1520 -9.154902 aq isoleucine 11 -10 6 1 15 -1 0 1521 -9.096910 aq leucine 12 -10 6 2 14 -1 0 1523 -9.000000 aq lysinium 13 -8 5 1 11 -1 1 1525 -8.552842 aq methionine 14 -16 9 1 20 -1 0 1526 -9.301030 aq phenylalanine 15 -8 5 1 11 -1 0 1527 -9.301030 aq proline 16 -3 3 1 5 -1 0 1528 -8.337242 aq serine 17 -5 4 1 8 -1 0 1529 -8.236572 aq threonine 18 -20 11 2 23 -2 0 1530 -9.221849 aq tryptophan 19 -15 9 1 19 -2 0 1532 -9.045757 aq tyrosinate 20 -8 5 1 12 -1 0 1533 -8.552842 aq valine
T <- 18 TK <- convert(T, "K") # calculate A/2.303RT (dimensionless), convert to G of reaction (cal/mol) a <- affinity(T=T)
energy.args: temperature is 18 C energy.args: pressure is Psat subcrt: 26 species at 291.15 K and 1 bar (wet)
G.18.cal <- convert(unlist(a$values), "G", T=TK) # covvert to kJ/mol G.18.kJ <- convert(G.18.cal, "J")/1000 # the 100 degree C case basis(c("H2O", "CO2", "NH4+", "H2", "H+", "H2S"), log10(c(1, 2.2e-3, 2.9e-6, 3.4e-4, 1.9e-6, 1.6e-3)))
C H N O S Z ispecies logact state H2O 0 2 0 1 0 0 1 0.000000 liq CO2 1 0 0 2 0 0 69 -2.657577 aq NH4+ 0 4 1 0 0 1 18 -5.537602 aq H2 0 2 0 0 0 0 66 -3.468521 aq H+ 0 1 0 0 0 1 3 -5.721246 aq H2S 0 2 0 0 1 0 70 -2.795880 aq
species(1:20, log10(c(2.8e-9, 5.0e-10, 7.9e-10, 2.4e-9, 3.6e-10, 2.7e-9, 7.2e-10, 4.2e-9, 8.6e-10, 5.0e-10, 5.7e-10, 7.2e-10, 2.0e-9, 3.6e-10,3.6e-10, 3.3e-9, 4.2e-9, 4.3e-10, 6.5e-10, 2.0e-9)))
H2O CO2 NH4+ H2 H+ H2S ispecies logact state name 1 -4 3 1 6 -1 0 1504 -8.552842 aq alanine 2 -10 6 4 11 -3 0 1506 -9.301030 aq argininium 3 -5 4 2 6 -2 0 1508 -9.102373 aq asparagine 4 -4 4 1 6 -2 0 1510 -8.619789 aq aspartate 5 -4 3 1 5 -2 1 1512 -9.443697 aq cysteinate 6 -6 5 1 9 -2 0 1515 -8.568636 aq glutamate 7 -7 5 2 9 -2 0 1513 -9.142668 aq glutamine 8 -2 2 1 3 -1 0 1516 -8.376751 aq glycine 9 -10 6 3 10 -2 0 1519 -9.065502 aq histidinium 10 -10 6 1 15 -1 0 1520 -9.301030 aq isoleucine 11 -10 6 1 15 -1 0 1521 -9.244125 aq leucine 12 -10 6 2 14 -1 0 1523 -9.142668 aq lysinium 13 -8 5 1 11 -1 1 1525 -8.698970 aq methionine 14 -16 9 1 20 -1 0 1526 -9.443697 aq phenylalanine 15 -8 5 1 11 -1 0 1527 -9.443697 aq proline 16 -3 3 1 5 -1 0 1528 -8.481486 aq serine 17 -5 4 1 8 -1 0 1529 -8.376751 aq threonine 18 -20 11 2 23 -2 0 1530 -9.366532 aq tryptophan 19 -15 9 1 19 -2 0 1532 -9.187087 aq tyrosinate 20 -8 5 1 12 -1 0 1533 -8.698970 aq valine
T <- 100 TK <- convert(T, "K") a <- affinity(T=T)
energy.args: temperature is 100 C energy.args: pressure is Psat subcrt: 26 species at 373.15 K and 1.01 bar (wet)
G.100.cal <- convert(unlist(a$values), "G", T=TK) G.100.kJ <- convert(G.100.cal, "J")/1000 # the average oxidation states of carbon Z.C <- ZC(thermo$obigt$formula[thermo$species$ispecies]) # put everything together a la Table 3 in the paper print(out <- data.frame(G.18=G.18.kJ, G.100=G.100.kJ, Z.C=Z.C))
G.18 G.100 Z.C 1504 113.62972 -8.862465 0.0000000 1506 409.40186 202.950494 0.3333333 1508 201.53021 87.204836 1.0000000 1510 146.70540 36.990300 1.0000000 1512 225.85075 74.224858 0.6666667 1515 172.08452 4.507084 0.4000000 1513 223.31409 49.159517 0.4000000 1516 80.46687 16.572202 1.0000000 1519 362.22122 164.727944 0.6666667 1520 213.88013 -88.795949 -1.0000000 1521 204.97626 -97.981989 -1.0000000 1523 258.50429 -21.063315 -0.6666667 1525 293.76048 11.800832 -0.4000000 1526 303.57593 -104.416187 -0.4444444 1527 192.79146 -33.119766 -0.4000000 1528 173.65949 72.592217 0.6666667 1529 216.47082 58.124509 0.0000000 1530 431.08196 -26.742386 -0.1818182 1532 339.43679 -33.563139 -0.2222222 1533 177.95803 -64.020557 -0.8000000
# make a plot; set units to get correct label E.units("J")
[1] "J"
plot(out$Z.C, out$G.18, pch=20, xlim=c(-1.1, 1.1), ylim=c(-200, 500), xlab=axis.label("ZC"), ylab=axis.label("DGr")) points(out$Z.C, out$G.100, col="red", pch=20) legend("topleft", pch=c(20, 20), col=c("black", "red"), legend=describe.property(c("T", "T"), c(18, 100))) title(main="Amino acid synthesis, after Amend and Shock, 1998")

Image affinity2

 

# 9 amino acids have negative delta Gr under hydrothermal conditions # (cf. AS98 with 11; we are using more recent thermodynamic data) stopifnot(sum(out$G.100 < 0)==9) # reset units and species to run next examples E.units("cal")
[1] "cal"
species(delete=TRUE)
NULL
## calculations along a transect: methanogenesis and biosynthetic ## reactions in hydrothermal systems, after Shock and Canovas, 2010 # this file has their mixing path results for Rainbow hydrothermal field file <- system.file("extdata/cpetc/SC10_Rainbow.csv", package="CHNOSZ") rb <- read.csv(file, check.names=FALSE) # write all synthesis reactions in terms of these basis species # it's okay not to set the activities of the basis species now # because they'll be changing along with temperature basis(c("CO2", "H2", "NH4+", "H2O", "H2S", "H+"))
C H N O S Z ispecies logact state CO2 1 0 0 2 0 0 69 0 aq H2 0 2 0 0 0 0 66 0 aq NH4+ 0 4 1 0 0 1 18 0 aq H2O 0 2 0 1 0 0 1 0 liq H2S 0 2 0 0 1 0 70 0 aq H+ 0 1 0 0 0 1 3 0 aq
# now a selection of the species from SC10, with activities equal to 1e-6 species(c("CH4", "formaldehyde", "ethylene", "glycolic acid", "n-nonanoic acid", "leucine", "aspartic acid", "tryptophan", "deoxyribose", "adenine", "cytosine"), -6)
CO2 H2 NH4+ H2O H2S H+ ispecies logact state name 1 1 4 0 -2 0 0 80 -6 aq methane 2 1 2 0 -1 0 0 255 -6 aq formaldehyde 3 2 6 0 -4 0 0 88 -6 aq ethylene 4 2 3 0 -1 0 0 573 -6 aq glycolic acid 5 9 25 0 -16 0 0 522 -6 aq n-nonanoic acid 6 6 15 1 -10 0 -1 1521 -6 aq leucine 7 4 6 1 -4 0 -1 1509 -6 aq aspartic acid 8 11 23 2 -20 0 -2 1530 -6 aq tryptophan 9 5 11 0 -6 0 0 1584 -6 aq deoxyribose 10 5 5 5 -10 0 -5 1578 -6 aq adenine 11 4 5 3 -7 0 -3 1580 -6 aq cytosine
# the exception is methane; unlike SC10 we use a constant activity 1e-3 # (accounting for variable activities of the species of interest here # is possible but would require longer code ....) species("CH4", -3)
CO2 H2 NH4+ H2O H2S H+ ispecies logact state name 1 1 4 0 -2 0 0 80 -6 aq methane 2 1 2 0 -1 0 0 255 -6 aq formaldehyde 3 2 6 0 -4 0 0 88 -6 aq ethylene 4 2 3 0 -1 0 0 573 -6 aq glycolic acid 5 9 25 0 -16 0 0 522 -6 aq n-nonanoic acid 6 6 15 1 -10 0 -1 1521 -6 aq leucine 7 4 6 1 -4 0 -1 1509 -6 aq aspartic acid 8 11 23 2 -20 0 -2 1530 -6 aq tryptophan 9 5 11 0 -6 0 0 1584 -6 aq deoxyribose 10 5 5 5 -10 0 -5 1578 -6 aq adenine 11 4 5 3 -7 0 -3 1580 -6 aq cytosine
# synchronized change of temperature and five basis activities a <- affinity(T=rb$T, CO2=rb$CO2, H2=rb$H2, `NH4+`=rb$`NH4+`, H2S=rb$H2S, pH=rb$pH)
energy.args: pressure is Psat energy.args: variable 1 is T at 18 values from 278.596 to 622.8 K energy.args: variable 2 is log_a(CO2) at 18 values from -3.6506 to -1.8288 energy.args: variable 3 is log_a(H2) at 18 values from -47.4549 to -1.8288 energy.args: variable 4 is log_a(NH4+) at 18 values from -6.2165 to -4.5641 energy.args: variable 5 is log_a(H2S) at 18 values from -5.2891 to -2.9213 energy.args: variable 6 is pH at 18 values from 3.0312 to 7.1784 subcrt: 17 species at 18 values of T and P (wet)
# the tricky part: affinity() uses dimensionless values (A/2.303RT) # but we want to show the values in cal/mol a$values <- lapply(a$values, function(val) { -convert(val, "G", T=convert(a$vals[[1]], "K")) }) # if we didn't have balance=1 here the values would be # divided by the number of moles of CO2 in the reactions ... diagram(a, balance=1, ylim=c(-100000, 100000), ylab=axis.label("A"), col=topo.colors(4), lwd=2)
balance: coefficients are unity diagram: plotting A/2.303RT from affinity(), divided by balancing coefficients
# add a zero-affinity line and a title abline(h=0, lty=2, lwd=2) title(main="Affinities of organic synthesis, after Shock and Canovas, 2010")

Image affinity3

 


next up previous
Next: util.affinity Up: CHNOSZ examples Previous: species