subcrt {CHNOSZ} | R Documentation |
Calculate the thermodynamic properties of one or more species or a reaction between species as a function of temperature and pressure.
subcrt(species, coeff = 1, state = NULL, property = c('logK','G','H','S','V','Cp'), T = seq(273.15,623.15,25), P = 'Psat', grid = NULL, convert = TRUE, do.phases = TRUE, logact = NULL, action.unbalanced = 'warn', IS = 0) read.supcrt(file) write.supcrt(file = 'supcrt.dat', obigt = thermo$obigt)
species |
character, name or formula of species, or numeric, rownumber of species in thermo$obigt . |
coeff |
numeric, reaction coefficients on species. |
state |
character, state(s) of species. |
property |
character, property(s) to calculate. |
T |
numeric, temperature(s) of the calculation. |
P |
numeric, pressure(s) of the calculation, or character, Psat. |
grid |
character, type of P xT grid to produce (NULL, the default, means no gridding). |
do.phases |
logical, replace Gibbs energies of unstable phases of minerals and of other species beyond their upper temperature limits of stability or extrapolation with 999999? |
logact |
numeric, logarithms of activities of species in reaction. |
file |
character, name of SUPCRT92 data file. |
obigt |
dataframe, thermodynamic data. |
convert |
logical, are input and output units of T and P those of the user (TRUE ), or are they Kelvin and bar (FALSE )? |
action.unbalanced |
character warn or NULL, what action to take if unbalanced reaction is provided. |
IS |
numeric, ionic strength(s) at which to calculate apparent molal properties, mol kg^-1. |
subcrt
calculates the standard molal thermodynamic properties of species and reactions as a function of temperature and pressure. For each of the species
(a formula or name), optionally identified in a given state
, the standard molal thermodynamic properties and equations-of-state parameters are retrieved via info
(except for H2O liquid). The standard molal properties of the species are computed using equations-of-state functions for aqueous species (hkf
), crystalline, gas, and liquid species (cgl
) and liquid or supercritical H2O (water
).
T
and P
denote the temperature and pressure conditions for the calculations and should generally be of the same length, unless P
is Psat (the default; see water
). Argument grid
if present can be one of T
or P
to flag the computation of a T
xP
or P
xT
grid. The property
s to be calculated can be one or more of those shown below:
rho | Density of water | g cm^-3 |
logK | Logarithm of equilibrium constant | dimensionless |
G | Gibbs energy | (cal | J) mol^-1 |
H | Enthalpy | (cal | J) mol^-1 |
S | Entropy | (cal | J) K^-1 mol^-1 |
V | Volume | cm^3 mol^-1 |
Cp | Heat capacity | (cal | J) K^-1 mol^-1 |
E | Exapansibility | cm^3 K^-1 |
kT | Isothermal compressibility | cm^3 bar^-1 |
(Note that E
and kT
can only be calculated for aqueous species and only if the option (thermo$opt$water
) for calculations of properties using water
is set to IAPWS
. On the other hand, if the water
option is SUPCRT (the default), E
and kT
can be calculated for water but not for aqueous species. (This is not an inherent limitation in either formulation, but it is just a matter of implementation.))
Depending on the global units definition (see nuts
) the values of G
, H
, S
and Cp
are returned in calories or Joule, but only if convert
is TRUE
. Likewise, the input values of T
and P
should be in units specified through nuts
, but setting convert
to FALSE
forces subcrt
to treat them as Kelvin and bar, respectively.
A chemical reaction is defined if coeff
is given. In this mode the standard molal properties of species are summed according to the stoichiometric coeff
icients, where negative values denote reactants. Reactions that do not conserve elements are permitted; subcrt
prints the missing composition needed to balance the reaction and produces a warning but computes anyway. Alternatively, if the basis
species of a system were previously defined, and if the species being considered are within the compositional range of the basis species, an unbalanced reaction given in the arguments to subcrt
will be balanced automatically, without altering the coefficients on the species identified in the arguments (unless perhaps they correspond to basis species), and without a warning. However, if a reaction is unbalanced and action.unbalanced
is set to NULL, no warnings are generated and no attempt is made to balance the reaction.
Minerals with phase transitions (denoted by having states cr1, cr2 etc.) can be defined generically by cr in the state
argument. subcrt
in this case simultaneously calculates the requested properties of all the species representing the phases of each such mineral (phase species) and, using the values of the transition temperatures calculated from those at P = 1 bar given in the thermodynamic database together with functions of the entropies and volumes of transitions (see dPdTtr
), determines the stable phase of the mineral at any grid point and substitutes the properties of this phase at that point for further calculations. If individual phase species of minerals are specified (by cr1, cr2 etc. in state
), the Gibbs energies of these minerals are assigned values of 999999 at conditions where another of the phases is stable (only if do.phases
is TRUE
). (NOTE: if you set do.phases
to FALSE
while investigating the properties of phases of minerals identified generically (by cr
) the function will identify the stable phase on the basis not of the transition temperatures but of the calculated Gibbs energies of the phase species. This is generally not advised. Computations of the phase transition properties of minerals are aided by the concept of phase species, but the computed standard molal properties of a phase species lose their physical meaning if computed beyond the transition temperatures of the corresponding phase.)
The chemical affinities of reactions are calculated if logact
is provided. This argument indicates the logarithms of activities (fugacities for gases) of species in the reaction; if there are fewer values of logact
than number of species those values are repeated as necessary. If the reaction was unbalanced to start, the logarithms of activities of any basis species added to the reaction are taken from the global basis definition. Columns appended to the output are logQ
for the activity product of the reaction, and A
for the chemical affinity (the latter in units corresponding to those of Gibbs energy). Note that affinity
provides related functionality but is geared toward the properties of formation reactions of species from the basis species and can be performed in more dimensions. Calculations of chemical affinity in subcrt
can be performed for any reaction of interest; however, they are currently limited to constant values of the logarithms of activities of species in the reactions, and hence of logQ
, across the computational range.
If IS
is set to a single value other than zero, function nonideal
is used to calculate the apparent properties (G
, H
, S
and Cp
) of charged aqueous species at the given ionic strength. To perform calculations at a single P
and T
and for multiple values of ionic strength, supply these values in IS
. Calculations can also be performed on a P
-IS
, T
-IS
or P,T
-IS
grid. Values of logK
of reactions calculated for IS
not equal to zero are consistent with the apparent Gibbs energies of the charged aqueous species.
subcrt
is modeled after the SUPCRT92 package (Johnson et al., 1992). Certain features of SUPCRT92 are not available here, for example, calculations as a function of density of H2O instead of pressure, or calculations of temperatures of univariant curves (i.e. for which logK
is zero), or calculations of the molar volumes of quartz and coesite as a function of temperature (taken here to be constant). The informative messages produced by SUPCRT92
when temperature or pressure limits of the equations of state are exceeded generally are not reproduced here. However, NA
s may be produced in the output of subcrt
if the requisite thermodynamic or electrostatic properties of water can not be calculated at given conditions.
Comparison of the output of subcrt
with that of SUPCRT92 indicates the two give similar values of properties for neutral aqueous species up to 1000 degrees C and 5000 bar and for charged aqueous species over the temperature range 0 to 300 degrees C. At higher temperatures, there are significant divergences for charged species (e.g., less than a 2 percent difference in the value of G for K+(aq) at 1000 deg C and 5000 bar). Further testing is necessary in CHNOSZ in connection with the g- and f-function equations (Shock et al., 1992) for the partial derivatives of the omega parameter which have only recently been incorporated into the hkf
function. (Note in particular the NaCl dissociation example under water
, the results of which are noticeably different from the calculations of Shock et al., 1992.)
read.supcrt
and write.supcrt
are used to read and write thermodynamic data file
s in the format of SUPCRT92
. read.supcrt
parses the thermodynamic data contained in a SUPCRT92-format file into a dataframe that is compatible with the format used in CHNOSZ (see thermo$obigt
). write.supcrt
creates a SUPCRT92 file
using the dataframe given in the obigt
argument (if missing the entire species database is processed). An edited version of the slop98.dat data file (Shock et al., 1998) that is compatible with read.supcrt
is available at http://chnosz.net/download/ (small formatting inconsistencies in the original version cause glitches with this function).
For subcrt
, a list of length two or three. If the properties of a reaction were calculated, the first element of the list (named reaction) contains a dataframe with the reaction parameters; the second element, named out, is a dataframe containing the calculated properties. Otherwise, the properties of species (not reactions) are returned: the first element, named species, contains a dataframe with the species identification; the second element, named out, is itself a list, each element of which is a dataframe of properties for a given species. If minerals with phase transitions are present, a third element (a dataframe) in the list indicates for all such minerals the stable phase at each grid point.
read.supcrt
returns a dataframe of properties and parameters of species.
Alberty, R. A., 2003. Thermodynamics of Biochemical Reactions. John Wiley & Sons, Hoboken, New Jersey, 397 p.
Johnson, J. W., Oelkers, E. H. and Helgeson, H. C., 1992. SUPCRT92: A software package for calculating the standard molal thermodynamic properties of minerals, gases, aqueous species, and reactions from 1 to 5000 bar and 0 to 1000degrees C. Comp. Geosci., 18, 899-947.
LaRowe, D. E. and Helgeson, H. C., 2007. Quantifying the energetics of metabolic reactions in diverse biogeochemical systems: electron flow and ATP synthesis. Geobiology, 5, 153-168.
Schulte, M. D. and Shock, E. L., 1995. Thermodynamics of Strecker synthesis in hydrothermal systems. Orig. Life Evol. Biosph., 25, 161-173.
Shock, E. L., Oelkers, E. H., Johnson, J. W., Sverjensky, D. A. and Helgeson, H. C., 1992. Calculation of the thermodynamic properties of aqueous species at high pressures and temperatures: Effective electrostatic radii, dissociation constants and standard partial molal properties to 1000 degrees C and 5 kbar. J. Chem. Soc. Faraday Trans., 88, 803-826.
Shock, E. L. et al., 1998. slop98.dat (computer data file). http://geopig.asu.edu/supcrt92_data/slop98.dat, accessed on 2005-11-05.
info
for finding species data, makeup
for parsing formulas to check reaction balance, water
, hkf
and cgl
for equations of state calculations.
## properties of species subcrt('water') # calculating at different temperatures subcrt('water',T=seq(0,100,10)) # calculating at even increments subcrt('water',T=seq(500,1000,length.out=10), P=seq(5000,10000,length.out=10)) # calculating on a temperature-pressure grid subcrt('water',T=c(500,1000),P=c(5000,10000),grid='P') # to calculate only selected properties subcrt('water',property=c('G','E')) # the properties of multiple species subcrt(c('glucose','ethanol','CO2')) ## input/output in different units nuts(c('K','J')) subcrt('water') subcrt('water',T=298.15) nuts(c('cal','C')) ## properties of reactions subcrt(c('H2O','H+','k-feldspar','kaolinite','K+','SiO2'), c(-1,-2,-2,1,2,4)) subcrt(c('glucose','ethanol','CO2'),c(-1,2,2)) # the states are specified subcrt(c('glucose','ethanol','CO2'),c(-1,2,2),c('aq','aq','gas')) # this warns about an unbalanced reaction subcrt(c('glucose','ethanol'),c(-1,1)) ## auto balancing reactions basis(c('CO2','H2O','NH3','H2S','O2')) subcrt(c('glucose','ethanol'),'aq',c(-1,3)) # properties of a species and a formation # reaction for that species subcrt('C2H5OH') # species basis('CHNOS') subcrt('C2H5OH',1) # reaction ## properties of mineral phases subcrt(c('pyrrhotite','pyrrhotite'),state=c('cr1','cr2')) # calculation of phase transitions of mineral subcrt('pyrrhotite') # phase transitions in a reaction subcrt(c('pyrite','pyrrhotite','H2O','H2S','O2'),c(-1,1,-1,1,0.5)) ## heat capacity of Fe(cr) # compare calculated values of heat capacity with # values from Robie and Hemingway, 1995, (from which # the parameters in the database were derived) nuts(c('J','K')) # set the units # we set pressure here otherwise subcrt goes for # Psat (saturation vapor pressure of H2O above # 100 degrees C) which can not be calculated above # the critical point of H2O (~647 K) t <- subcrt('Fe',T=seq(300,1800,50),P=1) plot(t$out[[1]]$T,t$out[[1]]$Cp,type='l', xlab=axis.label('T'),ylab=axis.label('Cp')) # add points from RH95 RH95 <- thermo$expt$RH95 points(RH95[,1],RH95[,2]) title(main=paste('Heat capacity of Fe(cr)\n', '(points - Robie and Hemingway, 1995)')) # reset the units nuts(c('C','cal')) ## these produce NAs and warnings # above the limits for calculating Psat subcrt('alanine',T=seq(0,5000,by=1000)) # above the T, P limits for the H2O equations of state subcrt('alanine',T=c(2250,2251),P=c(30000,30001),grid='T') ## Skarn example after Johnson et al., 1992 P <- seq(500,5000,500) # this is like the temperature specification used # in the example by Johnson et al., 1992 # T <- seq(0,1000,100) # we use this one to avoid warnings at 0 deg C, 5000 bar T <- c(2,seq(100,1000,100)) subcrt(c('andradite','carbon dioxide','H2S','Cu+','quartz','calcite', 'chalcopyrite','H+','H2O'),coeff=c(-1,-3,-4,-2,3,3,2,2,3), state=c('cr','g','aq','aq','cr','cr','cr','aq','liq'), P=P,T=T,grid='P') ## Standard Gibbs energy of reactions with HCN and ## formaldehyde after Schulte and Shock, 1995 Fig. 2 rxn1 <- subcrt(c('formaldehyde','HCN','H2O','glycolic acid','NH3'), c(-1,-1,-2,1,1),P=300) rxn2 <- subcrt(c('formaldehyde','HCN','H2O','glycine'), c(-1,-1,-1,1),P=300) plot(x=rxn1$out$T,rxn1$out$G/1000,type='l',ylim=c(-40,-10), xlab=axis.label('T'),ylab=axis.label('DG0r','k')) lines(rxn1$out$T,rxn2$out$G/1000) # write the reactions on the plot text(150,-14,describe(rxn1$reaction, use.name=c(TRUE,FALSE,FALSE,TRUE,FALSE))) text(200,-35,describe(rxn2$reaction, use.name=c(TRUE,FALSE,FALSE,TRUE))) title(main=paste('Standard Gibbs energy of reactions\n', 'after Schulte and Shock, 1995')) ## Calculation of chemical affinities # after LaRowe and Helgeson, 2007 # Fig. 3 (a): reduction of nicotinamide adenine # dinucleotide (NAD) coupled to oxidation of glucose # list the available NAD species info('NAD ') T <- seq(0,120,10) # oxidation of glucose (C6H12O6) basis(c('glucose','H2O','NH3','CO2','H+'),c(-3,0,999,-3,-7)) t <- subcrt(c('NAD(ox)-','NAD(red)-2'),c(-12,12),logact=c(0,0),T=T) # LH07's diagrams are shown per mole of electron (24 e- per 12 NAD) A <- t$out$A/24/1000 plot(x=T,y=A,xlim=range(T),ylim=c(1.4,5.4), xlab=axis.label('T'),ylab=axis.label('A',opt='k'),type='l') text('NAD(ox)-/NAD(red)-2 = 1',x=median(T),y=median(A)) # different activity ratio t <- subcrt(c('NAD(ox)-','NAD(red)-2'),c(-12,12),logact=c(1,0),T=T) A <- t$out$A/24/1000 lines(x=T,y=A) text('NAD(ox)-/NAD(red)-2 = 10',x=median(T),y=median(A)) # different activity ratio t <- subcrt(c('NAD(ox)-','NAD(red)-2'),c(-12,12),logact=c(0,1),T=T) A <- t$out$A/24/1000 lines(x=T,y=t$out$A/24/1000) text('NAD(ox)-/NAD(red)-2 = 0.1',x=median(T),y=median(A)) # this command prints the reaction on the plot text(40,4.5,c2s(s2c(describe(t$reaction, use.name=c(TRUE,TRUE,FALSE,FALSE,FALSE,FALSE)), sep='=',move.sep=TRUE),sep='\n')) # label the plot title(main=paste('Chemical affinity of NAD reduction', 'after LaRowe and Helgeson, 2007',sep='\n'), sub=describe(thermo$basis,T=NULL)) ### ### non-ideality calculations -- activity coefficients of ### aqueous species as a function of charge, temperature, ### and ionic strength -- after Alberty, 2003 ## ## p. 16 Table 1.3 apparent pKa of acetic acid with ## changing ionic strength subcrt(c('acetic acid','acetate','H+'),c(-1,1,1), IS=c(0,0.1,0.25),T=25,property='logK') # note that these *apparent* values of G and logK approach # their *standard* counterparts as IS goes to zero. ## ## p. 95: basis and elemental stoichiometries of species ## (a digression here from the nonideality calculations) # note coefficient of O2 and NH3 will be zero for these species basis(c('ATP-4','H+','H2O','HPO4-2','O2','NH3')) # cf Eq. 5.1-33: (basis composition) species(c('ATP-4','H+','H2O','HPO4-2','ADP-3','HATP-3','HADP-2','H2PO4-')) lb <- length(basis()) # cf Eq. 5.1-32: (elemental composition) as.matrix(species()[,1:lb]) ## ## p. 273-275: activity coefficient (gamma) ## as a function of ionic strength and temperature ## (as of 20080304, these do look quantitatively different ## from the plots in Alberty's book.) iplotfun <- function(T,col,add=TRUE) { IS <- seq(0,0.25,0.0025) s <- subcrt(c('H2PO4-','HADP-2','HATP-3','ATP-4'),IS=IS,grid='IS',T=T) if(!add) thermo.plot.new(xlim=range(IS),ylim=c(0,1), xlab=axis.label('IS'),ylab='gamma') for(i in 1:4) lines(IS,10^s$out[[i]]$loggam,col=col) } iplotfun(0,'blue',add=FALSE) iplotfun(25,'black') iplotfun(40,'red') title(main=paste('activity coefficient (gamma) of -1,-2,-3,-4', 'charged species at 0, 25, 40 deg C, after Alberty, 2003',sep='\n'))