Skip to content

Commit

Permalink
Improve formatting of ΔG°f and ΔH°f
Browse files Browse the repository at this point in the history
git-svn-id: svn://scm.r-forge.r-project.org/svnroot/chnosz/pkg/CHNOSZ@797 edb9625f-4e0d-4859-8d74-9fd3b1da38cb
  • Loading branch information
jedick committed Jul 6, 2023
1 parent 04aaee8 commit 7b76817
Show file tree
Hide file tree
Showing 6 changed files with 18 additions and 18 deletions.
4 changes: 2 additions & 2 deletions DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
Date: 2023-07-01
Date: 2023-07-05
Package: CHNOSZ
Version: 2.0.0-16
Version: 2.0.0-17
Title: Thermodynamic Calculations and Diagrams for Geochemistry
Authors@R: c(
person("Jeffrey", "Dick", , "[email protected]", role = c("aut", "cre"),
Expand Down
4 changes: 2 additions & 2 deletions R/info.R
Original file line number Diff line number Diff line change
Expand Up @@ -236,7 +236,7 @@ info.numeric <- function(ispecies, check.it = TRUE) {
calcCp <- check.EOS(this, this$model, "Cp", return.difference = FALSE)
# Fill in NA heat capacity
if(!is.na(calcCp) & is.na(this$Cp)) {
message("info.numeric: Cp of ", this$name, "(", this$state, ") is NA; set by EOS parameters to ", round(calcCp, 2), " ", this$E_units, " K-1 mol-1")
message("info.numeric: Cp\u00B0 of ", this$name, "(", this$state, ") is NA; set by EOS parameters to ", round(calcCp, 2), " ", this$E_units, " K-1 mol-1")
this$Cp <- as.numeric(calcCp)
} else if(isBerman) {
# Calculate Cp for Berman minerals 20220208
Expand All @@ -248,7 +248,7 @@ info.numeric <- function(ispecies, check.it = TRUE) {
calcV <- check.EOS(this, this$model, "V", return.difference = FALSE)
# Fill in NA volume
if(!is.na(calcV) & is.na(this$V)) {
message("info.numeric: V of ", this$name, "(", this$state, ") is NA; set by EOS parameters to ", round(calcV, 2), " cm3 mol-1")
message("info.numeric: V\u00B0 of ", this$name, "(", this$state, ") is NA; set by EOS parameters to ", round(calcV, 2), " cm3 mol-1")
this$V <- as.numeric(calcV)
}
}
Expand Down
4 changes: 2 additions & 2 deletions R/util.data.R
Original file line number Diff line number Diff line change
Expand Up @@ -209,7 +209,7 @@ check.EOS <- function(eos, model, prop, return.difference = TRUE) {
if(!is.na(calcval)) {
if(!is.na(refval)) {
if(abs(diff) > tol) {
message(paste("check.EOS: calculated ", prop, " of ", eos$name, "(", eos$state,
message(paste("check.EOS: calculated ", prop, "\u00B0 of ", eos$name, "(", eos$state,
") differs by ", round(diff,2), " ", units, " from database value", sep = ""))
return(calcval)
}
Expand Down Expand Up @@ -243,7 +243,7 @@ check.GHS <- function(ghs, return.difference = TRUE) {
if(!is.na(refval)) {
diff <- calcval - refval
if(abs(diff) > thermo$opt$G.tol) {
message(paste("check.GHS: calculated G of ", ghs$name, "(", ghs$state,
message(paste("check.GHS: calculated \u0394G\u00B0f of ", ghs$name, "(", ghs$state,
") differs by ", round(diff), " ", ghs$E_units, " mol-1 from database value", sep = ""))
return(calcval)
}
Expand Down
4 changes: 2 additions & 2 deletions inst/tinytest/test-info.R
Original file line number Diff line number Diff line change
Expand Up @@ -14,9 +14,9 @@ expect_equal(CHNOSZ:::info.character("H2O", "aq"), CHNOSZ:::info.character("H2O"
info <- "info.numeric() produces expected errors and messages"
expect_error(CHNOSZ:::info.numeric(9999), "species index 9999 not found in thermo\\(\\)\\$OBIGT", info = info)
iargon <- info("argon", "gas")
expect_message(CHNOSZ:::info.numeric(iargon), "Cp of argon\\(gas\\) is NA; set by EOS parameters to 4.97", info = info)
expect_message(CHNOSZ:::info.numeric(iargon), "Cp° of argon\\(gas\\) is NA; set by EOS parameters to 4.97", info = info)
iMgSO4 <- info("MgSO4")
expect_message(CHNOSZ:::info.numeric(iMgSO4), "V of MgSO4\\(aq\\) is NA; set by EOS parameters to 1.34", info = info)
expect_message(CHNOSZ:::info.numeric(iMgSO4), "V° of MgSO4\\(aq\\) is NA; set by EOS parameters to 1.34", info = info)

info <- "info.approx() produces expected messages"
expect_message(CHNOSZ:::info.approx("lactic"), "is similar to lactic acid", info = info)
Expand Down
4 changes: 2 additions & 2 deletions inst/tinytest/test-util.data.R
Original file line number Diff line number Diff line change
Expand Up @@ -9,9 +9,9 @@ expect_equal(d1, d2, info = info)

info <- "check.GHS() and check.EOS() (via info()) produce messages"
i1 <- info("S2O3-2")
expect_message(info(i1), "calculated G of S2O3-2\\(aq\\) differs by 939 cal mol-1 from database value", info = info)
expect_message(info(i1), "calculated ΔG°f of S2O3-2\\(aq\\) differs by 939 cal mol-1 from database value", info = info)
i2 <- info("Cu+2")
expect_message(info(i2), "calculated Cp of Cu\\+2\\(aq\\) differs by 3.62 cal K-1 mol-1 from database value", info = info)
expect_message(info(i2), "calculated Cp° of Cu\\+2\\(aq\\) differs by 3.62 cal K-1 mol-1 from database value", info = info)

info <- "check.GHS() and check.EOS() respond to thermo()$opt$*.tol"
i1 <- info("SO4-2")
Expand Down
16 changes: 8 additions & 8 deletions vignettes/FAQ.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -250,7 +250,7 @@ mtitle(c("Inorganic and organic species", "C[total] = 3 molal"))

Many minerals undergo phase transitions upon heating.
The stable phases at different temperatures are called polymorphs.
Each polymorph for a given mineral should have its own entry in OBIGT, including values of the standard thermodynamic properties (Δ*G*°, Δ*H*°, and *S*°) at 25 °C.
Each polymorph for a given mineral should have its own entry in OBIGT, including values of the standard thermodynamic properties (Δ*G*°~*f*~, Δ*H*°~*f*~, and *S*°) at 25 °C.
The 25 °C (or 298.15 K) values for high-temperature polymorphs are often not listed in thermodynamic tables, but they can be calculated.
This thermodynamic cycle shows how: we calculate the changes of a thermodyamic property (pictured here as `DS1` and `DS2`) between 298.15 K and the transition temperature (`Ttr`) for two polymorphs, then combine those with the property of the polymorphic transition (`DStr`) to obtain the difference of the property between the polymorphs at 298.15 K (`DS298`).

Expand Down Expand Up @@ -312,7 +312,7 @@ DStr <- TDStr / Ttr
```

Start new database entries that include basic information, volume, and heat capacity coefficients for each polymorph.
@PMW87 don't list `Cp` of the high-temperature polymorph extrapolated to 298.15 K, so leave it out.
@PMW87 don't list *C~p~*° of the high-temperature polymorph extrapolated to 298.15 K, so leave it out.
If the properties were in Joules, we would omit `E_units = "cal"` or change it to `E_units = "J"`.
```{r pyrrhotite_Cp, results = "hide", message = FALSE}
mod.OBIGT("pyrrhotite_new", formula = formula, state = "cr", ref1 = ref1,
Expand All @@ -324,14 +324,14 @@ mod.OBIGT("pyrrhotite_new", formula = formula, state = "cr2", ref1 = ref1,
```

For the time being, we set `G`, `H`, and `S` (i.e., the properties at 25 °C) to zero in order to easily calculate the temperature integrals of the properties from 298.15 K to `Ttr`.
Values of zero are placeholders that don't satisfy Δ*G*° = Δ*H*° − *T*Δ*S*° for the formation properties from the elements for either polymorph, as indicated by the following messages.
Values of zero are placeholders that don't satisfy Δ*G*°~*f*~ = Δ*H*°~*f*~*T*Δ*S*°~*f*~ for either polymorph (the subscript *f* represents formation from the elements), as indicated by the following messages.
We will check again for consistency of the thermodynamic parameters at the end of the example.
```{r pyrrhotite_info, results = "hide", collapse = TRUE}
info(info("pyrrhotite_new", "cr"))
info(info("pyrrhotite_new", "cr2"))
```

In order to calculate the temperature integral of Δ*G*°, we need not only the heat capacity coefficients but also actual values of *S*°.
In order to calculate the temperature integral of Δ*G*°~*f*~, we need not only the heat capacity coefficients but also actual values of *S*°.
Therefore, we start by calculating the entropy changes of each polymorph from 298.15 to `r Ttr` K (`DS1` and `DS2`) and combining those with the entropy of transition to obtain the difference of entropy between the polymorphs at 298.15 K.
For polymorph 1 (with `state = "cr"`) it's advisable to include `use.polymorphs = FALSE` to prevent `subcrt()` from trying to identify the most stable polymorph at the indicated temperature.
```{r pyrrhotite_S, message = FALSE}
Expand All @@ -353,7 +353,7 @@ It's a good idea to check that the entropy of transition is calculated correctly
stopifnot(all.equal(D2$S - D1$S, DStr))
```

Now we're ready to add up the contributions to Δ*G*° and Δ*H*° from the left, top, and right sides of the cycle.
Now we're ready to add up the contributions to Δ*G*°~*f*~ and Δ*H*°~*f*~ from the left, top, and right sides of the cycle.
This gives us the differences between the polymorphs at 298.15 K, which we use to make the final changes to the database.
```{r pyrrhotite_DG298_DH298, results = "hide", message = FALSE}
DG298 <- D1$G + DGtr - D2$G
Expand All @@ -365,7 +365,7 @@ mod.OBIGT("pyrrhotite_new", state = "cr2", G = G1 + DG298, H = H1 + DH298)
It's a good idea to check that the values of `G`, `H`, and `S` at 25 °C for a given polymorph are consistent with each other.
Here we use `check.GHS()` to calculate the difference between the value given for `G` and the value calculated from `H` and `S`.
The difference of less than 1 `r E.units()`/mol can probably be attributed to small differences in the entropies of the elements used by @PMW87 and in CHNOSZ.
We still get a message that the database value of `Cp` at 25 °C for the high-temperature polymorph is NA; this is OK because the (extrapolated) value can be calculated from the heat capacity coefficients.
We still get a message that the database value of *C~p~*° at 25 °C for the high-temperature polymorph is NA; this is OK because the (extrapolated) value can be calculated from the heat capacity coefficients.
```{r pyrrhotite_info2, collapse = TRUE}
cr_parameters <- info(info("pyrrhotite_new", "cr"))
stopifnot(abs(check.GHS(cr_parameters)) < 1)
Expand All @@ -381,15 +381,15 @@ cr2_parameters

Finally, let's look at the thermodynamic properties of the newly added pyrrhotite as a function of temperature around `Ttr`.
Here, we use the feature of `subcrt()` that identifies the stable polymorph at each temperature.
Note that Δ*G*° is a continuous function -- visual confirmation that the parameters yield zero for `DGtr` -- but Δ*H*°, *S*°, and *Cp*° are discontinuous at the transition temperature.
Note that Δ*G*°~*f*~ is a continuous function -- visual confirmation that the parameters yield zero for `DGtr` -- but Δ*H*°~*f*~, *S*°, and *C~p~*° are discontinuous at the transition temperature.

<button id="B-pyrrhotite_T" onclick="ToggleDiv('pyrrhotite_T')">Show code</button>
<div id="D-pyrrhotite_T" style="display: none">
```{r pyrrhotite_T, eval = FALSE}
T <- seq(550, 650, 1)
sout <- subcrt("pyrrhotite_new", T = T, P = 1)$out[[1]]
par(mfrow = c(1, 4), mar = c(4, 4.2, 1, 1))
labels <- c(G = "DG0", H = "DH0", S = "S0", Cp = "Cp0")
labels <- c(G = "DG0f", H = "DH0f", S = "S0", Cp = "Cp0")
for(property in c("G", "H", "S", "Cp")) {
plot(sout$T, sout[, property], col = sout$polymorph,
xlab = axis.label("T"), ylab = axis.label(labels[property])
Expand Down

0 comments on commit 7b76817

Please sign in to comment.