[CHNOSZ-commits] r213 - in pkg/CHNOSZ: . R inst man/macros tests/testthat
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Sep 25 10:11:41 CEST 2017
Author: jedick
Date: 2017-09-25 10:11:41 +0200 (Mon, 25 Sep 2017)
New Revision: 213
Added:
pkg/CHNOSZ/R/DEW.R
pkg/CHNOSZ/tests/testthat/test-DEW.R
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/NAMESPACE
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/macros/macros.Rd
Log:
add equations for Deep Earth Water (DEW) model
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2017-09-25 06:49:53 UTC (rev 212)
+++ pkg/CHNOSZ/DESCRIPTION 2017-09-25 08:11:41 UTC (rev 213)
@@ -1,6 +1,6 @@
Date: 2017-09-25
Package: CHNOSZ
-Version: 1.1.0-11
+Version: 1.1.0-12
Title: Thermodynamic Calculations for Geobiochemistry
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
Modified: pkg/CHNOSZ/NAMESPACE
===================================================================
--- pkg/CHNOSZ/NAMESPACE 2017-09-25 06:49:53 UTC (rev 212)
+++ pkg/CHNOSZ/NAMESPACE 2017-09-25 08:11:41 UTC (rev 213)
@@ -57,7 +57,7 @@
"water.SUPCRT92", "water.props", "eqdata", "plot_findit",
"nonideal", "anim.TCA", "uniprot.aa", "run.guess",
# added 20170301 or later
- "GHS_Tr"
+ "GHS_Tr", "calculateDensity", "calculateGibbsOfWater", "calculateEpsilon", "calculateQ"
)
# Load shared objects
Added: pkg/CHNOSZ/R/DEW.R
===================================================================
--- pkg/CHNOSZ/R/DEW.R (rev 0)
+++ pkg/CHNOSZ/R/DEW.R 2017-09-25 08:11:41 UTC (rev 213)
@@ -0,0 +1,202 @@
+# CHNOSZ/DEW.R
+# R functions for the Deep Earth Water model
+# 20170924 jmd
+
+# Most code was translated from VBA macros in the DEW spreadsheet (DEW_May19_2017_11.0.2 .xlsm)
+# Comments starting with ' were transferred from the DEW spreadsheet
+# In the original, functions return zero for invalid input; here, NA is used instead
+
+# Equations were selected using default options in the DEW spreadsheet:
+# Density of water equation (1) - Zhang & Duan (2005)
+# Dielectric Constant Equation (4) - Sverjensky et al. (2014)
+# Water Free Energy Equation (2) - Integral of Volume
+
+# 'Returns the density of water at the input pressure and temperature, in units of g/cm^3.
+# 'pressure - The pressure to calculate the density of water at, in bars
+# 'temperature - The temperature to calculate the density of water at, in Celsius
+# 'error - The density returned will calculate a pressure which differs from the input pressure by the value of "error" or less.
+# the default value of error is taken from equations in DEW spreadsheet (table "Calculations")
+calculateDensity <- function(pressure, temperature, error = 0.01) {
+ myfunction <- function(pressure, temperature) {
+ minGuess <- 1E-05
+ guess <- 1E-05
+ equation <- 1 # 'The maxGuess is dependent on the value of "equation"
+ maxGuess <- 7.5 * equation - 5
+ calcP <- 0
+ # 'Loop through and find the density
+ calculateDensity <- NA
+ for(i in 1:50) {
+ # 'Calculates the pressure using the specified equation
+ calcP <- calculatePressure(guess, temperature)
+ # 'If the calculated pressure is not equal to input pressure, this determines a new
+ # 'guess for the density based on current guess and how the calculated pressure
+ # 'relates to the input pressure. In effect, this a form of a bisection method.
+ if(abs(calcP - pressure) > error) {
+ if(calcP > pressure) {
+ maxGuess <- guess
+ guess <- (guess + minGuess) / 2
+ } else if(calcP < pressure) {
+ minGuess <- guess
+ guess <- (guess + maxGuess) / 2
+ }
+ } else {
+ calculateDensity <- guess
+ break
+ }
+ }
+ calculateDensity
+ }
+ # make input pressure and temperature the same length
+ if(length(pressure) < length(temperature)) pressure <- rep(pressure, length.out=length(temperature))
+ if(length(temperature) < length(pressure)) temperature <- rep(temperature, length.out=length(pressure))
+ # use a loop to process vectorized input
+ sapply(1:length(pressure), function(i) myfunction(pressure[i], temperature[i]))
+}
+
+# 'Returns the Gibbs Free Energy of water in units of cal/mol.
+# 'pressure - The pressure to calculate the Gibbs Free Energy at, in bars
+# 'temperature - The temperature to calculate the Gibbs Free Energy at, in Celsius
+calculateGibbsOfWater <- function(pressure, temperature) {
+
+ # 'Equation created by Brandon Harrison. This models data for the Gibbs energy at 1 kb as a function of temperature,
+ # 'then defines the gibbs free energy as the integral over the volume as a function of temperature.
+
+ myfunction <- function(pressure, temperature) {
+ # 'Gibbs Free Energy of water at 1 kb. This equation is a polynomial fit to data as a function of temperature.
+ # 'It is valid in the range of 100 to 1000 C.
+ GAtOneKb <- 2.6880734E-09 * temperature^4 + 6.3163061E-07 * temperature^3 -
+ 0.019372355 * temperature^2 - 16.945093 * temperature - 55769.287
+
+ if(pressure < 1000) { # 'Simply return zero, this method only works at P >= 1000 bars
+ integral <- NA
+ } else if(pressure == 1000) { # 'Return the value calculated above from the polynomial fit
+ integral <- 0
+ } else if(pressure > 1000) { # 'Integrate from 1 kb to P over the volume
+ integral <- 0
+ # 'Integral is sum of rectangles with this width. This function in effect limits the spacing
+ # 'to 20 bars so that very small pressures do not have unreasonably small widths. Otherwise the width
+ # 'is chosen such that there are always 500 steps in the numerical integration. This ensures that for very
+ # 'high pressures, there are not a huge number of steps calculated which is very computationally taxing.
+ spacing <- ifelse((pressure - 1000) / 500 < 20, 20, (pressure - 1000) / 500)
+ for(i in seq(1000, pressure, by = spacing)) {
+ # 'This integral determines the density only down to an error of 100 bars
+ # 'rather than the standard of 0.01. This is done to save computational
+ # 'time. Tests indicate this reduces the computation by about a half while
+ # 'introducing little error from the standard of 0.01.
+ integral <- integral + (18.01528 / calculateDensity(i, temperature, 100) / 41.84) * spacing
+ }
+ }
+ GAtOneKb + integral
+ }
+ # make input pressure and temperature the same length
+ if(length(pressure) < length(temperature)) pressure <- rep(pressure, length.out=length(temperature))
+ if(length(temperature) < length(pressure)) temperature <- rep(temperature, length.out=length(pressure))
+ # use a loop to process vectorized input
+ sapply(1:length(pressure), function(i) myfunction(pressure[i], temperature[i]))
+}
+
+# 'Returns the Dielectric constant of water at the given density and temperature.
+# 'density - The density of water to use in calculating epsilon, in g/cm^3
+# 'temperature - The temperature to calculate epsilon with, in Celsius
+calculateEpsilon <- function(density, temperature) {
+ # 'Power Function - Created by Dimitri Sverjensky and Brandon Harrison
+ # 'Relevant parameters
+ a1 <- -0.00157637700752506
+ a2 <- 0.0681028783422197
+ a3 <- 0.754875480393944
+ b1 <- -8.01665106535394E-05
+ b2 <- -0.0687161761831994
+ b3 <- 4.74797272182151
+
+ A <- a1 * temperature + a2 * sqrt(temperature) + a3
+ B <- b1 * temperature + b2 * sqrt(temperature) + b3
+
+ exp(B) * density ^ A
+}
+
+# 'Outputs the value of Q in units of bar^-1
+# 'pressure - The pressure to calculate Q at, in bars
+# 'temperature - The temperature to calculate Q at, in Celsius
+# 'density - The density at the input pressure and temperature, input simply to save time, in g/cm^3
+calculateQ <- function(density, temperature) {
+ eps <- calculateEpsilon(density, temperature)
+ depsdrho <- calculate_depsdrho(density, temperature)
+ drhodP <- calculate_drhodP(density, temperature)
+
+ depsdrho * drhodP / eps ^2
+}
+
+### unexported functions ###
+
+# 'Returns the pressure of water corresponding to the input density and temperature, in units of bars.
+# 'density - The density to use in finding a pressure, in g/cm^3
+# 'temperature - The temperature to use in finding a pressure, in Celsius
+calculatePressure <- function(density, temperature) {
+ m <- 18.01528 # 'Molar mass of water molecule in units of g/mol
+ ZD05_R <- 83.144 # 'Gas Constant in units of cm^3 bar/mol/K
+ ZD05_Vc <- 55.9480373 # 'Critical volume in units of cm^3/mol
+ ZD05_Tc <- 647.25 # 'Critical temperature in units of Kelvin
+
+ TK <- temperature + 273.15 # 'Temperature must be converted to Kelvin
+ Vr <- m / density / ZD05_Vc
+ Tr <- TK / ZD05_Tc
+
+ B <- 0.349824207 - 2.91046273 / (Tr * Tr) + 2.00914688 / (Tr * Tr * Tr)
+ C <- 0.112819964 + 0.748997714 / (Tr * Tr) - 0.87320704 / (Tr * Tr * Tr)
+ D <- 0.0170609505 - 0.0146355822 / (Tr * Tr) + 0.0579768283 / (Tr * Tr * Tr)
+ E <- -0.000841246372 + 0.00495186474 / (Tr * Tr) - 0.00916248538 / (Tr * Tr * Tr)
+ f <- -0.100358152 / Tr
+ g <- -0.00182674744 * Tr
+
+ delta <- 1 + B / Vr + C / (Vr * Vr) + D / Vr^4 + E / Vr^5 + (f / (Vr * Vr) + g / Vr^4) * exp(-0.0105999998 / (Vr * Vr))
+
+ ZD05_R * TK * density * delta / m
+}
+
+# 'Calculates the partial derivative of density with respect to pressure, i.e. (d(rho)/dP)_T
+# 'density - The density of water, in g/cm^3
+# 'temperature - The temperature of water, in Celsius
+calculate_drhodP <- function(density, temperature) {
+ m <- 18.01528 # 'Molar mass of water molecule in units of g/mol
+ ZD05_R <- 83.144 # 'Gas Constant in units of cm^3 bar/mol/K
+ ZD05_Vc <- 55.9480373 # 'Critical volume in units of cm^3/mol
+ ZD05_Tc <- 647.25 # 'Critical temperature in units of Kelvin
+
+ TK <- temperature + 273.15 # 'temperature must be converted to Kelvin
+ Tr <- TK / ZD05_Tc
+ cc <- ZD05_Vc / m # 'This term appears frequently in the equation and is defined here for convenience
+ Vr <- m / (density * ZD05_Vc)
+
+ B <- 0.349824207 - 2.91046273 / (Tr * Tr) + 2.00914688 / (Tr * Tr * Tr)
+ C <- 0.112819964 + 0.748997714 / (Tr * Tr) - 0.87320704 / (Tr * Tr * Tr)
+ D <- 0.0170609505 - 0.0146355822 / (Tr * Tr) + 0.0579768283 / (Tr * Tr * Tr)
+ E <- -0.000841246372 + 0.00495186474 / (Tr * Tr) - 0.00916248538 / (Tr * Tr * Tr)
+ f <- -0.100358152 / Tr
+ g <- 0.0105999998 * Tr
+
+ delta <- 1 + B / Vr + C / (Vr^2) + D / Vr^4 + E / Vr^5 + (f / (Vr^2) + g / Vr^4) * exp(-0.0105999998 / Vr^2)
+
+ kappa <- B * cc + 2 * C * (cc^2) * density + 4 * D * cc^4 * density^3 + 5 * E * cc^5 * density^4 +
+ (2 * f * (cc^2) * density + 4 * g * cc^4 * density^3 - (f / (Vr^2) + g / Vr^4) * (2 * 0.0105999998 * (cc^2) * density)) * exp(-0.0105999998 / (Vr^2))
+
+ m / (ZD05_R * TK * (delta + density * kappa))
+}
+
+# 'Returns the partial derivative of the dielectric constant with respect to density in units of cm^3/g.
+# 'density - The density of water to calculate with, in g/cm^3
+# 'temperature - The temperature to calculate with, in Celsius
+calculate_depsdrho <- function(density, temperature) {
+ # 'Power Function - Created by Dimitri Sverjensky and Brandon Harrison
+ # 'Relevant parameters
+ a1 <- -0.00157637700752506
+ a2 <- 0.0681028783422197
+ a3 <- 0.754875480393944
+ b1 <- -8.01665106535394E-05
+ b2 <- -0.0687161761831994
+ b3 <- 4.74797272182151
+
+ A <- a1 * temperature + a2 * sqrt(temperature) + a3
+ B <- b1 * temperature + b2 * sqrt(temperature) + b3
+
+ A * exp(B) * density ^ (A - 1)
+}
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2017-09-25 06:49:53 UTC (rev 212)
+++ pkg/CHNOSZ/inst/NEWS 2017-09-25 08:11:41 UTC (rev 213)
@@ -1,5 +1,5 @@
-CHANGES IN CHNOSZ 1.1.0-11 (2017-09-25)
---------------------------------------
+CHANGES IN CHNOSZ 1.1.0-12 (2017-09-25)
+---------------------------------------
- water.lines() now works for diagrams of Eh, pe, logfO2, logaO2,
logfH2, or logaH2 vs pH, T, or P. It is possible to have T or P on
@@ -33,6 +33,10 @@
to calculate (at P, T) so that the calculations aren't duplicated by
subcrt().
+- Add equations from the Deep Earth Water (DEW) model (Sverjensky
+ et al., 2014) in the exported functions calculateDensity(),
+ calculateGibbsofWater(), calculateEpsilon(), and calculateQ().
+
CHANGES IN CHNOSZ 1.1.0 (2017-05-04)
------------------------------------
Modified: pkg/CHNOSZ/man/macros/macros.Rd
===================================================================
--- pkg/CHNOSZ/man/macros/macros.Rd 2017-09-25 06:49:53 UTC (rev 212)
+++ pkg/CHNOSZ/man/macros/macros.Rd 2017-09-25 08:11:41 UTC (rev 213)
@@ -23,3 +23,8 @@
% use \nH2O{̄} to call this macro (the html code can't be defined in the macro,
% which interprets '#' followed by a number as a placeholder for an argument)
\newcommand{\nH2O}{\ifelse{latex}{\eqn{\bar{n}_{\mathrm{H_{2}O}}}}{\ifelse{html}{\out{<i>n</i>#1<sub>H<sub>2</sub>O</sub>}}{nH2O}}}
+
+% other symbols
+\newcommand{\degC}{\ifelse{latex}{\eqn{^{\circ}}C}{\ifelse{html}{\out{°}C}{°C}}}
+\newcommand{\le}{\ifelse{latex}{\eqn{\le}}{\ifelse{html}{\out{≤}}{≤}}}
+\newcommand{\ge}{\ifelse{latex}{\eqn{\ge}}{\ifelse{html}{\out{≥}}{≥}}}
Added: pkg/CHNOSZ/tests/testthat/test-DEW.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-DEW.R (rev 0)
+++ pkg/CHNOSZ/tests/testthat/test-DEW.R 2017-09-25 08:11:41 UTC (rev 213)
@@ -0,0 +1,43 @@
+context("DEW")
+
+test_that("density of water is calculated correctly", {
+ pressure <- c(5000, 5000, 10000, 10000, 20000, 20000, 50000, 50000)
+ temperature <- c(100, 1000, 100, 1000, 100, 1000, 100, 1000)
+ # density from R functions
+ RDensity <- calculateDensity(pressure, temperature)
+ # density from DEW spreadsheet
+ DEWDensity <- c(1.108200, 0.597623, 1.196591, 0.798331, 1.321050, 1.000735, 1.578116, 1.287663)
+ expect_equal(RDensity, DEWDensity, tolerance=1e-6)
+})
+
+test_that("Gibbs energy of water is calculated correctly", {
+ pressure <- c(5000, 5000, 10000, 10000, 20000, 20000, 50000, 50000)
+ temperature <- c(100, 1000, 100, 1000, 100, 1000, 100, 1000)
+ # Gibbs energies from R functions
+ RGibbs <- calculateGibbsOfWater(pressure, temperature)
+ # Gibbs energies from DEW spreadsheet
+ DEWGibbs <- c(-56019.85419280258, -84262.028821198, -54155.004480575895, -81210.38766217149,
+ -50735.122222685815, -76433.07602205424, -41823.26077175943, -65187.48113532527)
+ expect_equal(RGibbs, DEWGibbs)
+})
+
+test_that("dielectric constant of water is calculated correctly", {
+ pressure <- c(5000, 5000, 10000, 10000, 20000, 20000, 50000, 50000)
+ temperature <- c(100, 1000, 100, 1000, 100, 1000, 100, 1000)
+ # epsilon from R functions
+ Repsilon <- calculateEpsilon(calculateDensity(pressure, temperature), temperature)
+ # epsilon from DEW spreadsheet
+ DEWepsilon <- c(65.63571, 6.10465, 72.40050, 8.97800, 82.16244, 12.13131, 103.12897, 16.97266)
+ expect_equal(Repsilon, DEWepsilon, tolerance=1e-7)
+})
+
+test_that("Born coefficient Q is calculated correctly", {
+ pressure <- c(5000, 5000, 10000, 10000, 20000, 20000, 50000, 50000)
+ temperature <- c(100, 1000, 100, 1000, 100, 1000, 100, 1000)
+ # Q from R functions
+ RQ <- calculateQ(calculateDensity(pressure, temperature), temperature)
+ # Q from DEW spreadsheet
+ DEWQ <- c(0.32319817, 14.50286092, 0.19453478, 3.12650897,
+ 0.10918151, 0.87729257, 0.05068788, 0.20640645) / 1e6
+ expect_equal(RQ, DEWQ)
+})
More information about the CHNOSZ-commits
mailing list