next up previous
Next: [3] diagram (1 variable) Up: CHNOSZ examples Previous: util.affinity

equil.boltzmann

## Don't show: data(thermo)
thermo$obigt: 1809 aqueous, 3368 total species
## End(Don't show) ## equilibrium in a simple system: ## ionization of carbonic acid basis("CHNOS+")
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 O2 0 0 0 2 0 0 3095 -80 gas H+ 0 1 0 0 0 1 3 -7 aq
species(c("CO2", "HCO3-", "CO3-2"))
CO2 H2O NH3 H2S O2 H+ ispecies logact state name 1 1 0 0 0 0 0 69 -3 aq CO2 2 1 1 0 0 0 -1 13 -3 aq HCO3- 3 1 1 0 0 0 -2 14 -3 aq CO3-2
# set unit activity of the species (0 = log10(1)) species(1:3, 0)
CO2 H2O NH3 H2S O2 H+ ispecies logact state name 1 1 0 0 0 0 0 69 0 aq CO2 2 1 1 0 0 0 -1 13 0 aq HCO3- 3 1 1 0 0 0 -2 14 0 aq CO3-2
# calculate Astar (for unit activity) res <- 100 Astar <- affinity(pH=c(0, 14, res))$values
energy.args: temperature is 25 C energy.args: pressure is Psat energy.args: variable 1 is pH at 100 values from 0 to 14 subcrt: 9 species at 298.15 K and 1 bar (wet)
# the logarithms of activity for a total activity # of the balanced quantity (C or CO2) equal to 0.001 loga.boltz <- equil.boltzmann(Astar, c(1, 1, 1), 0.001) # calculated another way loga.react <- equil.reaction(Astar, c(1, 1, 1), 0.001) # probably close enough for most purposes stopifnot(all.equal(loga.boltz, loga.react)) # the first ionization constant (pKa) ipKa <- which.min(abs(loga.boltz[[1]] - loga.boltz[[2]])) pKa.equil <- seq(0, 14, length.out=res)[ipKa] # calculate logK directly logK <- subcrt(c("CO2","H2O","HCO3-","H+"), c(-1, -1, 1, 1), T=25)$out$logK
info.character: found CO2(aq), also available in gas info.character: found H2O(liq), also available in gas subcrt: 4 species at 298.15 K and 1 bar (wet)
# we could decrease tolerance here by increasing res stopifnot(all.equal(pKa.equil, -logK, tolerance=1e-2))


next up previous
Next: [3] diagram (1 variable) Up: CHNOSZ examples Previous: util.affinity