[cig-commits] [commit] add_thermodynamic_potentials: Added endmember Gibbs with Landau, and a new test script. (cdb8c44)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue Dec 9 09:53:14 PST 2014
Repository : https://github.com/geodynamics/burnman
On branch : add_thermodynamic_potentials
Link : https://github.com/geodynamics/burnman/compare/2e5646d93cedbbf8eae54cc37cffc14e0aa85180...d5ddad03ff9f30f5a4efaddb4e3ec585ea1a7c51
>---------------------------------------------------------------
commit cdb8c44f9d1eafd69f748f7f8505f1a20fef8606
Author: Bob Myhill <myhill.bob at gmail.com>
Date: Mon Aug 11 19:23:10 2014 +0200
Added endmember Gibbs with Landau, and a new test script.
>---------------------------------------------------------------
cdb8c44f9d1eafd69f748f7f8505f1a20fef8606
burnman/mineral.py | 8 +-
burnman/minerals/HP_2011.py | 73 ++++++++++--
burnman/minerals/SLB_2011.py | 11 +-
.../minerals/{SLB_2011.py => SLB_2011_working.py} | 0
burnman/minerals/__init__.py | 1 +
burnman/modified_tait.py | 95 +++++++++++-----
burnman/processchemistry.py | 39 +++++++
burnman/testing.py | 126 +++++++++++++++++++++
8 files changed, 310 insertions(+), 43 deletions(-)
diff --git a/burnman/mineral.py b/burnman/mineral.py
index d6c4d87..8d2a8ed 100644
--- a/burnman/mineral.py
+++ b/burnman/mineral.py
@@ -122,12 +122,18 @@ class Mineral(Material):
self.C_v = self.method.heat_capacity_v(self.pressure, self.temperature, self.V, self.params)
self.C_p = self.method.heat_capacity_p(self.pressure, self.temperature, self.V, self.params)
self.alpha = self.method.thermal_expansivity(self.pressure, self.temperature, self.V, self.params)
+ self.gibbs = self.method.gibbs(self.pressure, self.temperature, self.params)
if (self.params.has_key('G_0') and self.params.has_key('Gprime_0')):
self.G = self.method.shear_modulus(self.pressure, self.temperature, self.V, self.params)
else:
self.G = float('nan') #nan if there is no G, this should propagate through calculations to the end
- warnings.warn(('Warning: G_0 and or Gprime_0 are undefined for ' + self.to_string()))
+ if self.params['equation_of_state'] != 'mtait':
+ warnings.warn(('Warning: G_0 and or Gprime_0 are undefined for ' + self.to_string()))
+
+ # The following gibbs function avoids having to calculate a bunch of unnecessary parameters over P-T space. This will be useful for gibbs minimisation.
+ def calcgibbs(self, pressure, temperature):
+ return self.method.gibbs(pressure, temperature, self.params)
def molar_mass(self):
"""
diff --git a/burnman/minerals/HP_2011.py b/burnman/minerals/HP_2011.py
index f4cd3af..ea35dbc 100644
--- a/burnman/minerals/HP_2011.py
+++ b/burnman/minerals/HP_2011.py
@@ -12,14 +12,14 @@ Minerals from Holland and Powell 2011 and references therein.
Note the units in Holland and Powell's Table 1 are not SI. They are:
H_0 [kJ/mol] i.e. multiply by 1e3 to get J/mol
err_H_0 [kJ/mol] i.e. multiply by 1e3 to get J/mol
-S [J/mol] i.e. S.I!
+S [J/K/mol] i.e. S.I!
V [kJ/kbar/mol] i.e. multiply by 1e-5 to get m^3/mol
Cp [kJ/K/mol] is given as a + bT + cT^-2 + dT^-0.5.
b is multiplied by 1e5 to be more readable.
Thus, conversions to SI are [1e3, 1e-2, 1e3, 1e3]
a_0 [1e5 K^-1] i.e. multiply by 1e-5 to get K^-1
K_0 [kbar] i.e. multiply by 1e8 to get Pa
-Kprime_0 [] i.e. SI!
+Kprime_0 [] i.e. S.I.!
Kdprime_0 [kbar^-1] i.e. multiply by 1e-8
The values in the database itself are NOT the same as the paper. They are strictly in the units kJ, kbar, K.
@@ -27,28 +27,85 @@ The values in the database itself are NOT the same as the paper. They are strict
There are also parameters which deal with internal order-disorder and Landau transitions. NOTE: These still need to be implemented.
"""
+from burnman.processchemistry import ProcessChemistry
from burnman.mineral import Mineral
-from burnman.solidsolution import SolidSolution
+#from burnman.solidsolution import SolidSolution
class stishovite (Mineral):
"""
Holland and Powell, 2011 and references therein
"""
def __init__(self):
+ formula='SiO2'
self.params = {
- 'formula': 'SiO2',
+ 'formula': formula,
'equation_of_state': 'mtait',
'H_0': -876.39e3,
'S_0': 24.0,
'V_0': 1.401e-5,
'Cp': [68.1,6.01e-3,-1.9782e6,-82.1],
- 'a_0': 15.8e-6
- 'K_0': 3090.e8
+ 'a_0': 15.8e-6,
+ 'K_0': 3090.e8,
'Kprime_0': 4.6,
'Kdprime_0': -0.00150e-8,
- 'molar_mass': .0601,
- 'n': 3}
+ 'n': ProcessChemistry(formula)[0],
+ 'molar_mass': ProcessChemistry(formula)[1]}
self.uncertainties = {
'err_H_0':0.49e3}
+class spinel (Mineral):
+ """
+ Holland and Powell, 2011 and references therein
+ """
+ def __init__(self):
+ formula='MgAl2O4'
+ self.params = {
+ 'formula': formula,
+ 'equation_of_state': 'mtait',
+ 'H_0': -2301.26e3,
+ 'S_0': 82.0,
+ 'V_0': 3.978e-5,
+ 'Cp': [222.9,6.127e-3,-1.686e6,-1551.0],
+ 'a_0': 1.93e-5,
+ 'K_0': 1922.e8,
+ 'Kprime_0': 4.04,
+ 'Kdprime_0': -0.00210e-8,
+ 'n': ProcessChemistry(formula)[0],
+ 'molar_mass': ProcessChemistry(formula)[1],
+ 'BW_deltaH': 8.00e3,
+ 'BW_deltaV': 0.00,
+ 'BW_W': 1.20e3,
+ 'BW_Wv': 0.00,
+ 'BW_n': 2.0,
+ 'BW_factor': 0.50}
+
+ self.uncertainties = {
+ 'err_H_0':0.84e3}
+
+class quartz (Mineral):
+ """
+ Holland and Powell, 2011 and references therein
+ """
+ def __init__(self):
+ formula='SiO2'
+ self.params = {
+ 'formula': formula,
+ 'equation_of_state': 'mtait',
+ 'H_0': -910.70e3,
+ 'S_0': 41.43,
+ 'V_0': 2.269e-5,
+ 'Cp': [92.9,-6.42e-4,-7.149e5,-716.1],
+ 'a_0': 0.0e-6,
+ 'K_0': 730.e8,
+ 'Kprime_0': 6.0,
+ 'Kdprime_0': -0.00820e-8,
+ 'n': ProcessChemistry(formula)[0],
+ 'molar_mass': ProcessChemistry(formula)[1],
+ 'landau_Tc':847.,
+ 'landau_Smax':4.95,
+ 'landau_Vmax':0.1188e-5}
+
+ self.uncertainties = {
+ 'err_H_0':0.27e3}
+
diff --git a/burnman/minerals/SLB_2011.py b/burnman/minerals/SLB_2011.py
index b4ebebf..5be815e 100644
--- a/burnman/minerals/SLB_2011.py
+++ b/burnman/minerals/SLB_2011.py
@@ -13,7 +13,6 @@ Minerals from Stixrude & Lithgow-Bertelloni 2011 and references therein
import burnman.mineral_helpers as bmb
from burnman.mineral import Mineral
-from burnman.solidsolution import SolidSolution
class stishovite (Mineral):
"""
@@ -105,15 +104,11 @@ class wuestite (Mineral):
'err_eta_s_0':1.0}
-class ferropericlase(SolidSolution):
+class ferropericlase(bmb.HelperSolidSolution):
def __init__(self, x):
- model_type='asymmetric'
- base_material = [periclase(), wuestite()]
+ base_materials = [periclase(), wuestite()]
molar_fraction = [1. - x, 0. + x] # keep the 0.0 +, otherwise it is an array sometimes
- site_occupancy = [['x(Mg)',1. - x],['x(Fe)',0. + x]]
- interaction_parameter = [[4.e3,0.0,0.0]] # Wh + T*Ws + P*Wv, J/mol
- van_laar_parameter=[1.0,1.0]
- SolidSolution.__init__(self, base_material, molar_fraction, site_occupancy, interaction_parameter, van_laar_parameter)
+ bmb.HelperSolidSolution.__init__(self, base_materials, molar_fraction)
class mg_fe_perovskite(bmb.HelperSolidSolution):
def __init__(self, fe_num):
diff --git a/burnman/minerals/SLB_2011.py b/burnman/minerals/SLB_2011_working.py
similarity index 100%
copy from burnman/minerals/SLB_2011.py
copy to burnman/minerals/SLB_2011_working.py
diff --git a/burnman/minerals/__init__.py b/burnman/minerals/__init__.py
index 4ce6e2b..cd9c388 100644
--- a/burnman/minerals/__init__.py
+++ b/burnman/minerals/__init__.py
@@ -19,5 +19,6 @@ import Matas_etal_2007
import SLB_2011
import SLB_2011_ZSB_2013
import SLB_2005
+import HP_2011
import Murakami_2013
import other
diff --git a/burnman/modified_tait.py b/burnman/modified_tait.py
index 94c21be..b65f745 100644
--- a/burnman/modified_tait.py
+++ b/burnman/modified_tait.py
@@ -10,18 +10,24 @@ import burnman.equation_of_state as eos
T_0=298.15 # Standard temperature = 25 C
def einst(S, n):
- "
+ """
Einstein temperature
base of p.346, para.1
- "
+ """
return 10636./(S/n + 6.44)
def ksi(u):
- "
+ """
Einstein function to describe behaviour of ak
EQ 11+1
- "
- return pow(u,2.)*exp(u)/pow((exp(u)-1.), 2.)
+ """
+ return pow(u,2.)*np.exp(u)/pow((np.exp(u)-1.), 2.)
+
+def tait_constants(params):
+ a=(1.+params['Kprime_0'])/(1. + params['Kprime_0'] + params['K_0']*params['Kdprime_0'])
+ b=params['Kprime_0']/params['K_0'] - params['Kdprime_0']/(1. + params['Kprime_0'])
+ c=(1. + params['Kprime_0'] + params['K_0']*params['Kdprime_0'])/(params['Kprime_0']*params['Kprime_0'] + params['Kprime_0'] - params['K_0']*params['Kdprime_0'])
+ return a, b, c
def modified_tait(x, params):
"""
@@ -31,9 +37,7 @@ def modified_tait(x, params):
EQ 2
"""
- a=(1.+params['Kprime_0'])/(1. + params['Kprime_0'] + params['K_0']*params['Kdprime_0'])
- b=(params['Kprime_0'])/(params['K_0'] - (params['Kdprime_0'])/(1. + params['Kprime_0']))
- c=(1. + params['Kprime_0'] + params['K_0']*params['Kdprime_0'])/(params['Kprime_0']*params['Kprime_0'] + params['Kprime_0'] - params['K_0']*params['Kdprime_0'])
+ a, b, c = tait_constants(params)
return (pow((x + a - 1.) / a, -1./c) - 1.)/b
@@ -60,10 +64,8 @@ class MTaitBase(eos.EquationOfState):
Returns volume [m^3] as a function of pressure [Pa] and temperature [K]
EQ 12
"""
- a=(1.+params['Kprime_0'])/(1. + params['Kprime_0'] + params['K_0']*params['Kdprime_0'])
- b=(params['Kprime_0'])/(params['K_0'] - (params['Kdprime_0'])/(1. + params['Kprime_0']))
- c=(1. + params['Kprime_0'] + params['K_0']*params['Kdprime_0'])/(params['Kprime_0']*params['Kprime_0'] + params['Kprime_0'] - params['K_0']*params['Kdprime_0'])
+ a, b, c = tait_constants(params)
Pth=self.__rel_thermal_pressure(temperature,params)
x = 1 - a*( 1. - pow(( 1. + b*(pressure-Pth)), -1.0*c))
return x*params['V_0']
@@ -73,13 +75,10 @@ class MTaitBase(eos.EquationOfState):
Returns isothermal bulk modulus [Pa] as a function of pressure [Pa],
temperature [K], and volume [m^3]. EQ 13+2
"""
- a=(1.+params['Kprime_0'])/(1. + params['Kprime_0'] + params['K_0']*params['Kdprime_0'])
- b=(params['Kprime_0'])/(params['K_0'] - (params['Kdprime_0'])/(1. + params['Kprime_0']))
- c=(1. + params['Kprime_0'] + params['K_0']*params['Kdprime_0'])/(params['Kprime_0']*params['Kprime_0'] + params['Kprime_0'] - params['K_0']*params['Kdprime_0'])
-
+ a, b, c = tait_constants(params)
Pth=self.__rel_thermal_pressure(temperature,params)
psubpth=pressure-Pth
- return params['K_0']*(1. + b(psubpth))*(a + (1.-a)*pow((1. + b(psubpth)), c))
+ return params['K_0']*(1. + b*(psubpth))*(a + (1.-a)*pow((1. + b*(psubpth)), c))
#calculate the shear modulus as a function of P, V, and T
def shear_modulus(self, pressure, temperature, volume, params):
@@ -103,9 +102,7 @@ class MTaitBase(eos.EquationOfState):
Returns thermal expansivity at the pressure, temperature, and volume [1/K]
Replace -Pth in EQ 13+1 with P-Pth for non-ambient temperature
"""
- a=(1.+params['Kprime_0'])/(1. + params['Kprime_0'] + params['K_0']*params['Kdprime_0'])
- b=(params['Kprime_0'])/(params['K_0'] - (params['Kdprime_0'])/(1. + params['Kprime_0']))
- c=(1. + params['Kprime_0'] + params['K_0']*params['Kdprime_0'])/(params['Kprime_0']*params['Kprime_0'] + params['Kprime_0'] - params['K_0']*params['Kdprime_0'])
+ a, b, c = tait_constants(params)
Pth=self.__rel_thermal_pressure(temperature,params)
psubpth=pressure-Pth
ein=einst(params['S_0'], params['n'])
@@ -123,6 +120,19 @@ class MTaitBase(eos.EquationOfState):
Cp = params['Cp'][0] + params['Cp'][1]*temperature + params['Cp'][2]*pow(temperature,-2.) + params['Cp'][3]*pow(temperature,-0.5)
return Cp
+
+ def heat_capacity_p(self, pressure, temperature, volume, params):
+ """
+ Returns heat capacity at constant pressure at the pressure, temperature, and volume [J/K/mol]
+ Not yet implemented. Returns 0.
+ """
+ #alpha = self.thermal_expansivity(pressure, temperature, volume, params)
+ #gr = self.grueneisen_parameter(pressure, temperature, volume, params)
+ #C_v = self.heat_capacity_v(pressure, temperature, volume, params)
+ #C_p = C_v*(1. + gr * alpha * temperature)
+ return 0.
+
+
def adiabatic_bulk_modulus(self,pressure,temperature,volume,params):
"""
Returns adiabatic bulk modulus [Pa] as a function of pressure [Pa],
@@ -141,20 +151,53 @@ class MTaitBase(eos.EquationOfState):
and temperature [K].
"""
# Calculate temperature and pressure integrals
- a=(1.+params['Kprime_0'])/(1. + params['Kprime_0'] + params['K_0']*params['Kdprime_0'])
- b=(params['Kprime_0'])/(params['K_0'] - (params['Kdprime_0'])/(1. + params['Kprime_0']))
- c=(1. + params['Kprime_0'] + params['K_0']*params['Kdprime_0'])/(params['Kprime_0']*params['Kprime_0'] + params['Kprime_0'] - params['K_0']*params['Kdprime_0'])
+ a, b, c = tait_constants(params)
Pth=self.__rel_thermal_pressure(temperature,params)
psubpth=pressure-Pth
- intCpdT = (params['Cp'][0]*temperature + 0.5*params['Cp'][1]*pow(temperature,2.) - params['Cp'][2]/temperature + 2.*params['Cp'][3]*sqrt(temperature)) - (params['Cp'][0]*T_0 + 0.5*params['Cp'][1]*T_0*T_0 - params['Cp'][2]/T_0 + 2.0*params['Cp'][3]*sqrt(T_0))
+ intCpdT = (params['Cp'][0]*temperature + 0.5*params['Cp'][1]*pow(temperature,2.) - params['Cp'][2]/temperature + 2.*params['Cp'][3]*np.sqrt(temperature)) - (params['Cp'][0]*T_0 + 0.5*params['Cp'][1]*T_0*T_0 - params['Cp'][2]/T_0 + 2.0*params['Cp'][3]*np.sqrt(T_0))
- intCpoverTdT = (params['Cp'][0]*log(temperature) + params['Cp'][1]*temperature - 0.5*params['Cp'][2]/pow(temperature,2.) - 2.0*params['Cp'][3]/sqrt(temperature)) - (params['Cp'][0]*log(T_0) + params['Cp'][1]*T_0 - 0.5*params['Cp'][2]/(T_0*T_0) - 2.0*params['Cp'][3]/sqrt(T_0))
+ intCpoverTdT = (params['Cp'][0]*np.log(temperature) + params['Cp'][1]*temperature - 0.5*params['Cp'][2]/pow(temperature,2.) - 2.0*params['Cp'][3]/np.sqrt(temperature)) - (params['Cp'][0]*np.log(T_0) + params['Cp'][1]*T_0 - 0.5*params['Cp'][2]/(T_0*T_0) - 2.0*params['Cp'][3]/np.sqrt(T_0))
# EQ 13
intVdP = pressure*params['V_0']*(1. - a + (a*(pow((1.-b*Pth), 1.-c) - pow((1. + b*(pressure-Pth)), 1.-c))/(b*(c-1.)*pressure)))
- return params['H_0'] + intCpdT - temperature*(params['S_0'] + intCpoverTdT) + intVdP
+ # Add order-disorder terms if required
+ if params.has_key('landau_Tc'): # For a phase transition described by Landau term
+ # From Holland and Powell, 1996, corrected using
+ # landaunote.pdf on Tim Holland's esc web page
+ Tcstar=params['landau_Tc'] + (params['landau_Vmax']/params['landau_Smax'])*pressure
+ # Q_0 is Q at T0, P0?
+ Q_0=pow((params['landau_Tc']-T_0)/params['landau_Tc'],1./4.)
+
+ # Find state of ordering
+ # Note that Q can be > 1. This is a bit puzzling, and wasn't in the original
+ # (apparently incorrect) expressions by Holland and Powell, 1996
+ if Tcstar-temperature > 0.:
+ Q=pow((Tcstar-temperature)/params['landau_Tc'],1./4.)
+ else:
+ Q=0.0
+
+ # Vt is cryptically defined in landaunote.pdf.
+ # R.M. imagines we would use the expression in brackets in
+ # EQ13 of Holland and Powell (2011) but apparently Vt==1 fits
+ # the output from thermocalc version tc337L.
+ # This is possibly a bug in thermocalc.
+ Vt=1.# EQ 13 would give: 1. - a + (a*(pow((1.-b*Pth), 1.-c) - pow((1. + b*(pressure-Pth)), 1.-c))/(b*(c-1.)*pressure))
+ Gdisord=params['landau_Tc']*params['landau_Smax']*(pow(Q_0,2) - pow(Q_0,6)/3.0) - params['landau_Smax']*(Tcstar*pow(Q,2) - params['landau_Tc']*pow(Q,6)/3.0) - temperature*(params['landau_Smax']*(pow(Q_0,2) - pow(Q,2))) + pressure*(params['landau_Vmax']*pow(Q_0,2)*Vt)
+ else:
+ if params.has_key('BW_deltaH'): # Add Bragg-Williams disordering
+ #Q_0=
+ #Q=
+ #Gdisord=params['BW_factor']*
+ Gdisord=0.0
+ print 'WARNING: This phase has B-W order-disorder which has not been implemented yet'
+ else:
+ Gdisord=0.0
+
+
+
+ return params['H_0'] + intCpdT - temperature*(params['S_0'] + intCpoverTdT) + intVdP + Gdisord
# calculate P = P(T0) + Pth
def pressure(self, temperature, volume, params):
@@ -171,7 +214,7 @@ class MTaitBase(eos.EquationOfState):
ein=einst(params['S_0'],params['n'])
u=ein/T
u_0=ein/T_0
- P_th = params['a_0']*params['k_0']*ein/ksi(T_0)*((1./(exp(u))-1.))-(1./(exp(u_0))-1.))
+ P_th = params['a_0']*params['K_0']*ein/ksi(u_0)*((1./(np.exp(u)-1.))-(1./(np.exp(u_0)-1.)))
return P_th
diff --git a/burnman/processchemistry.py b/burnman/processchemistry.py
new file mode 100644
index 0000000..41a4bbd
--- /dev/null
+++ b/burnman/processchemistry.py
@@ -0,0 +1,39 @@
+# BurnMan - a lower mantle toolkit
+# Copyright (C) 2012, 2013, Myhill, R., Heister, T., Unterborn, C., Rose, I. and Cottaar, S.
+# Released under GPL v2 or later.
+
+# This is a function which returns the number of atoms and molar mass of a compound given its unit formula as an argument.
+
+import re
+
+def ProcessChemistry(formula):
+ filename="data/input_masses/atomic_masses.dat"
+ file = open(filename, "r")
+ el_name = []
+ el_mass = []
+
+ i=0
+ for line in file:
+ data="%".join(line.split("%")[:1]).split()
+ if data != []:
+ el_name.append(data[0])
+ el_mass.append(float(data[1]))
+ file.close()
+
+# Loop over elements
+ n=0.
+ molar_mass=0.
+ for element in re.findall('[A-Z][^A-Z]*', formula):
+ list=filter(None, re.split(r'(\d+)', element))
+ # Look up number of atoms of element
+ if len(list) == 1:
+ nel=1.
+ else:
+ nel=float(list[1])
+ # Increment size of compound
+ n=n+nel
+
+ # Find atomic mass of element
+ molar_mass=molar_mass+nel*el_mass[el_name.index(list[0])]
+
+ return n, molar_mass
diff --git a/burnman/testing.py b/burnman/testing.py
new file mode 100644
index 0000000..c779f73
--- /dev/null
+++ b/burnman/testing.py
@@ -0,0 +1,126 @@
+# BurnMan - a lower mantle toolkit
+# Copyright (C) 2012, 2013, Heister, T., Unterborn, C., Rose, I. and Cottaar, S.
+# Released under GPL v2 or later.
+
+
+
+# Here we import standard python modules that are required for
+# usage of BurnMan. In particular, numpy is used for handling
+# numerical arrays and mathematical operations on them, and
+# matplotlib is used for generating plots of results of calculations
+import os, sys, numpy as np, matplotlib.pyplot as plt
+
+#hack to allow scripts to be placed in subdirectories next to burnman:
+if not os.path.exists('burnman') and os.path.exists('../burnman'):
+ sys.path.insert(1,os.path.abspath('..'))
+
+
+# Here we import the relevant modules from BurnMan. The burnman
+# module imports several of the most important functionalities of
+# the library, including the ability to make composites, and compute
+# thermoelastic properties of them. The minerals module includes
+# the mineral physical parameters for the predefined minerals in
+# BurnMan
+import burnman
+from burnman import minerals
+
+
+# The following script checks the Gibbs Free energies of *solid* Holland and Powell (2011) endmembers satisfying one of several cases
+
+print 'N.B. Small differences (<100 J/mol) between the burnman and thermocalc outputs are presumably the result of Holland and Powell missing .15 K somewhere in their code, either in the standard state temperature (298.15 K), or in the conversion from Celsius to Kelvin.'
+print ''
+
+
+# A mineral without phase transitions: stishovite
+stv=minerals.HP_2011.stishovite()
+stv.set_method("mtait")
+
+Gstv=[['P \ T',25,500,1500],[0.001, -883.54, -909.13, -1020.78],[10.0, -869.55, -895.00, -1006.26],[100.0, -745.58, -769.83, -877.86]]
+
+print 'G(stv)'
+
+Pmaxerror=0.
+Tmaxerror=0.
+maxerror=0.
+serror=0.
+for pi in range(len(Gstv)-1):
+ P=Gstv[pi+1][0]*1e8
+ for ti in range(len(Gstv[0])-1):
+ T=Gstv[0][ti+1]+273.15
+ Gburnman=stv.calcgibbs(P,T)
+ GHP=Gstv[pi+1][ti+1]*1000. # convert to J/mol
+ Gdiff=Gburnman-GHP
+# print P, T, Gburnman, GHP, Gdiff
+ serror=serror + pow(Gdiff,2)
+ if abs(Gdiff) > abs(maxerror):
+ maxerror = Gdiff
+ Tmaxerror = T
+ Pmaxerror = P
+
+rmserror = np.sqrt(serror/((len(Gstv)-1)*(len(Gstv[0])-1)))
+print 'Root mean squared error:', rmserror, 'J/mol'
+print 'Maximum error:', maxerror, 'J/mol @', P/1.e9, 'GPa and', T, 'K'
+print ''
+
+# A mineral with a phase transition described by Landau theory: quartz
+q=minerals.HP_2011.quartz()
+q.set_method("mtait")
+
+Gq=[['P \ T',25,500,1500],[0.001,-923.06,-957.10,-1088.42],[10.0,-900.62,-934.16,-1064.93],[100.0,-715.40,-746.62,-870.48]]
+
+print 'G(q)'
+
+Pmaxerror=0.
+Tmaxerror=0.
+maxerror=0.
+serror=0.
+for pi in range(len(Gq)-1):
+ P=Gq[pi+1][0]*1e8
+ for ti in range(len(Gq[0])-1):
+ T=Gq[0][ti+1]+273.15
+ Gburnman=q.calcgibbs(P,T)
+ GHP=Gq[pi+1][ti+1]*1000. # convert to J/mol
+ Gdiff=Gburnman-GHP
+# print P, T, Gburnman, GHP, Gdiff
+ serror=serror + pow(Gdiff,2)
+ if abs(Gdiff) > abs(maxerror):
+ maxerror = Gdiff
+ Tmaxerror = T
+ Pmaxerror = P
+
+rmserror = np.sqrt(serror/((len(Gq)-1)*(len(Gq[0])-1)))
+print 'Root mean squared error:', rmserror, 'J/mol'
+print 'Maximum error:', maxerror, 'J/mol @', P/1.e9, 'GPa and', T, 'K'
+print ''
+
+
+# A mineral with a phase transition described by Bragg-Williams theory: spinel
+sp=minerals.HP_2011.spinel()
+sp.set_method("mtait")
+
+Gsp=[['P \ T',25,500,1500],[0.001,-2325.64,-2401.93,-2713.86],[10.0,-2285.96,-2361.80,-2672.59],[100.0,-1937.39,-2009.65,-2311.34]]
+
+print 'G(sp)'
+
+Pmaxerror=0.
+Tmaxerror=0.
+maxerror=0.
+serror=0.
+for pi in range(len(Gsp)-1):
+ P=Gsp[pi+1][0]*1e8
+ for ti in range(len(Gsp[0])-1):
+ T=Gsp[0][ti+1]+273.15
+ Gburnman=sp.calcgibbs(P,T)
+ GHP=Gsp[pi+1][ti+1]*1000. # convert to J/mol
+ Gdiff=Gburnman-GHP
+ print P, T, Gburnman, GHP, Gdiff
+ serror=serror + pow(Gdiff,2)
+ if abs(Gdiff) > abs(maxerror):
+ maxerror = Gdiff
+ Tmaxerror = T
+ Pmaxerror = P
+
+rmserror = np.sqrt(serror/((len(Gsp)-1)*(len(Gsp[0])-1)))
+print 'Root mean squared error:', rmserror, 'J/mol'
+print 'Maximum error:', maxerror, 'J/mol @', P/1.e9, 'GPa and', T, 'K'
+print ''
More information about the CIG-COMMITS
mailing list