[cig-commits] [commit] add_thermodynamic_potentials: First alpha of SolidSolution base class (8c5e99f)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue Dec 9 09:53:59 PST 2014
Repository : https://github.com/geodynamics/burnman
On branch : add_thermodynamic_potentials
Link : https://github.com/geodynamics/burnman/compare/2e5646d93cedbbf8eae54cc37cffc14e0aa85180...d5ddad03ff9f30f5a4efaddb4e3ec585ea1a7c51
>---------------------------------------------------------------
commit 8c5e99fa601d480389cc186e9df127ffda0588c1
Author: Bob Myhill <myhill.bob at gmail.com>
Date: Sat Aug 30 20:50:05 2014 +0200
First alpha of SolidSolution base class
>---------------------------------------------------------------
8c5e99fa601d480389cc186e9df127ffda0588c1
burnman/solidsolution.py | 156 +++++++++++++++------
burnman/test_process_solidsolution.py | 69 +++++----
burnman/testing.py | 257 ----------------------------------
3 files changed, 150 insertions(+), 332 deletions(-)
diff --git a/burnman/solidsolution.py b/burnman/solidsolution.py
index ebb75d2..44198d6 100644
--- a/burnman/solidsolution.py
+++ b/burnman/solidsolution.py
@@ -1,10 +1,12 @@
# BurnMan - a lower mantle toolkit
-# Copyright (C) 2012, 2013, Heister, T., Unterborn, C., Rose, I. and Cottaar, S.
+# Copyright (C) 2012-2014, Myhill, R., Heister, T., Unterborn, C., Rose, I. and Cottaar, S.
# Released under GPL v2 or later.
import warnings
+import re
import numpy as np
+from fractions import Fraction
kd = lambda x,y : 1 if x==y else 0
@@ -29,58 +31,130 @@ class SolidSolution(Mineral):
and P derivatives in J/K/mol and m^3/(mol molecule).
"""
- # init sets up matrices to speed up calculations when P, T, X is defined.
+ # init sets up matrices to speed up calculations for when P, T, X is defined.
def __init__(self, base_material, interaction_parameter):
+ # Initialise the solid soltion inputs
self.base_material = base_material
- self.molar_fraction = molar_fraction
- assert(len(base_material) == len(molar_fraction))
- assert(len(base_material) == len(van_laar_parameter))
+ self.interaction_parameter = interaction_parameter
+
+ # Number of endmembers in the solid solution
+ self.n_endmembers=len(base_material)
+
+ # Check that base_material and interaction parameter each have three parts
+ assert(len(map(list, zip(*base_material))) == 3)
+ assert(len(interaction_parameter) == 3)
+
+ # Check that the interaction parameters have the correct number of variables
+ for i in range(3):
+ assert(len(interaction_parameter[i]) == n_endmembers-1)
+ for j in range(n_endmembers-1):
+ assert(len(interaction_parameter[i][j]) == (n_endmembers-1)-j)
+
+ # Split interaction parameter into H, S, V terms
+ self.excess_enthalpy=interaction_parameter[0]
+ self.excess_entropy=interaction_parameter[1]
+ self.excess_volume=interaction_parameter[2]
+
+ # Number of sites
+ self.n_sites=base_material[0][1].count('[')
+
+ # Check the number of sites is the same for each endmember
+ for i in range(n_endmembers):
+ assert(base_material[i][1].count('[') == n_sites)
+
+ # Check the chemical composition of each endmember is consistent with the dataset
+ # NOT IMPLEMENTED YET
+
+ # Number of unique site occupancies (e.g.. Mg on X etc.)
+ sites=[[] for i in range(n_sites)]
+ list_occupancies=[]
+ list_multiplicity=np.empty(shape=(n_sites))
+ self.n_occupancies=0
+ for endmember in range(n_endmembers):
+ list_occupancies.append([[0]*len(sites[site]) for site in range(n_sites)])
+ s=re.split(r'\[', base_material[endmember][1])[1:]
+ for site in range(n_sites):
+ site_occupancy=re.split(r'\]', s[site])[0]
+ mult=re.split('[A-Z][^A-Z]*',re.split(r'\]', s[site])[1])[0]
+ if mult == '':
+ list_multiplicity[site]=1.0
+ else:
+ list_multiplicity[site]=mult
+ elements=re.findall('[A-Z][^A-Z]*',site_occupancy)
+ for i in range(len(elements)):
+ element_on_site=re.split('[0-9][^A-Z]*',elements[i])[0]
+ proportion_element_on_site=re.findall('[0-9][^A-Z]*',elements[i])
+ if len(proportion_element_on_site) == 0:
+ proportion_element_on_site=Fraction(1.0)
+ else:
+ proportion_element_on_site=Fraction(proportion_element_on_site[0])
+
+ if element_on_site not in sites[site]:
+ self.n_occupancies=n_occupancies+1
+ sites[site].append(element_on_site)
+ element_index=sites[site].index(element_on_site)
+ for parsed_mbr in range(len(list_occupancies)):
+ list_occupancies[parsed_mbr][site].append(0)
+ else:
+ element_index=sites[site].index(element_on_site)
+ list_occupancies[endmember][site][element_index]=proportion_element_on_site
+
+ # Site occupancies and multiplicities
+ self.site_occupancies=np.empty(shape=(n_endmembers,n_occupancies))
+ self.site_multiplicities=np.empty(shape=(n_occupancies))
+ for endmember in range(n_endmembers):
+ n_element=0
+ for site in range(n_sites):
+ for element in range(len(list_occupancies[endmember][site])):
+ self.site_occupancies[endmember][n_element]=list_occupancies[endmember][site][element]
+ self.site_multiplicities[n_element]=list_multiplicity[site]
+ n_element=n_element+1
+
+ self.alpha=np.array([base_material[i][2] for i in range(n_endmembers)])
+
+ self.Wh=np.zeros(shape=(n_endmembers,n_endmembers))
+ self.Ws=np.zeros(shape=(n_endmembers,n_endmembers))
+ self.Wv=np.zeros(shape=(n_endmembers,n_endmembers))
+ for i in range(n_endmembers):
+ for j in range(i+1, n_endmembers):
+ self.Wh[i][j]=2.*excess_enthalpy[i][j-i-1]/(alpha[i]+alpha[j])
+ self.Ws[i][j]=2.*excess_entropy[i][j-i-1]/(alpha[i]+alpha[j])
+ self.Wv[i][j]=2.*excess_volume[i][j-i-1]/(alpha[i]+alpha[j])
+
+ # END INITIALISATION
+
+
+ def set_composition(self, molar_fraction):
+ assert(len(self.base_material) == len(molar_fraction))
assert(sum(molar_fraction) > 0.9999)
assert(sum(molar_fraction) < 1.0001)
- # make sure the number of atoms in each endmember is the same
- # (when we start adding minerals with vacancies
- # this will need to be changed).
- for m in base_material:
- if(base_material[0].params.has_key('n')):
- assert(m.params['n'] == base_material[0].params['n'])
+ self.molar_fraction=molar_fraction
- # Set the state of all the endmembers
- def set_state(self, pressure, temperature):
- for mat in self.base_materials:
- mat.method = self.method
- mat.set_state(pressure, temperature)
+ phi=np.array([self.alpha[i]*molar_fraction[i] for i in range(self.n_endmembers)])
+ phi=np.divide(phi, np.sum(phi))
- itrange = range(0, len(self.base_material))
+ self.H_excess=np.dot(self.alpha.T,molar_fraction)*np.dot(phi.T,np.dot(self.Wh,phi))
+ self.S_deficit=np.dot(self.alpha.T,molar_fraction)*np.dot(phi.T,np.dot(self.Ws,phi))
+ self.V_excess=np.dot(self.alpha.T,molar_fraction)*np.dot(phi.T,np.dot(self.Wv,phi))
+ self.S_excess=0.0-S_deficit
- self.params = {}
-
- # NEEDS CHANGING FROM HERE!!!
- # Description of RTlngamma terms given by the asymmetric formalism
- # is included in Holland and Powell (2003)
-
- # the volume is calculated in the same way as Gex (p.493),
- # replacing Wij with Wvij
- Vex=-2.*sum([ sum([ self.molar_fraction[i]*self.van_laar_parameter[i]*self.molar_fraction[j]*self.van_laar_parameter[j]*self.interaction_parameter[i][2] / (self.van_laar_parameter[i] + self.van_laar_parameter[j]) for j in range(i, len(self.base_material)) ]) for i in range(0, len(self.base_material)-1) ])/sum([ self.van_laar_parameter[l]*self.molar_fraction[l] for l in itrange ])
- V= sum([ self.base_materials[i].molar_volume() * self.molar_fraction[i] for i in itrange ]) + Vex
-
- # the volume is calculated in the same way as Gex (p.493),
- # replacing Wij with Wvij
- # NEED TO CORRECT INDICES FOR INT PARAM
- # W[i][j] == W[(itrange-1)*i+(j-1)]
- Vex=-2.*sum([ sum([ self.molar_fraction[i]*self.van_laar_parameter[i]*self.molar_fraction[j]*self.van_laar_parameter[j]*self.interaction_parameter[(itrange-1)*i+(j-1)][2] / (self.van_laar_parameter[i] + self.van_laar_parameter[j]) for j in range(i, len(self.base_material)) ]) for i in range(0, len(self.base_material)-1) ])/sum([ self.van_laar_parameter[l]*self.molar_fraction[l] for l in itrange ])
- V= sum([ self.base_materials[i].molar_volume() * self.molar_fraction[i] for i in itrange ]) + Vex
-
-
- # Same for the Gibbs free energy
- # NEED TO CORRECT INDICES FOR INT PARAM
- Gex=-2.*sum([ sum([ self.molar_fraction[i]*self.van_laar_parameter[i]*self.molar_fraction[j]*self.van_laar_parameter[j]*(pressure*self.interaction_parameter[(itrange-1)*i+(j-1)][0] + temperature*self.interaction_parameter[(itrange-1)*i+(j-1)][1] + pressure*self.interaction_parameter[(itrange-1)*i+(j-1)][2]) / (self.van_laar_parameter[i] + self.van_laar_parameter[j]) for j in range(i, len(self.base_material)) ]) for i in range(0, len(self.base_material)-1) ])/sum([ self.van_laar_parameter[l]*self.molar_fraction[l] for l in itrange ])
- G= sum([ self.base_materials[i].molar_gibbs() * self.molar_fraction[i] for i in itrange ]) + Gex
+ def set_state(self, pressure, temperature, molar_fraction):
+ # Set the state of all the endmembers
+ for i in range(self.n_endmembers):
+ self.base_material[i][0].method = self.method
+ self.base_material[i][0].set_state(pressure, temperature)
+ # Find excess properties
+ self.set_composition(molar_fraction)
+ self.G_excess=self.H_excess - temperature*self.S_excess + pressure*self.V_excess
+ self.molar_volume= sum([ self.base_material[i][0].molar_volume() * self.molar_fraction[i] for i in n_endmembers ]) + self.V_excess
+ self.molar_gibbs= sum([ self.base_materials[i][0].molar_gibbs() * self.molar_fraction[i] for i in n_endmembers ]) + self.G_excess
+ '''
for prop in self.base_materials[0].params:
try:
self.params[prop] = sum([ self.base_materials[i].params[prop] * self.molar_fraction[i] for i in itrange ])
@@ -88,4 +162,6 @@ class SolidSolution(Mineral):
#if there is a type error, it is probably a string. Just go with the value of the first base_material.
self.params[prop] = self.base_materials[0].params[prop]
Mineral.set_state(self, pressure, temperature)
+ '''
+
diff --git a/burnman/test_process_solidsolution.py b/burnman/test_process_solidsolution.py
index dda3bff..92fbbde 100644
--- a/burnman/test_process_solidsolution.py
+++ b/burnman/test_process_solidsolution.py
@@ -13,12 +13,12 @@ P,T
# Things we can initialise without endmember proportions
n_sites
-n_elements
+n_occupancies
n_endmembers
alpha
Wh, Ws, Wv
-array_occupancies
-array_multiplicities
+site_occupancies
+site_multiplicities
# Endmember proportions needed
phi
@@ -49,30 +49,30 @@ interaction_parameter=[excess_enthalpy,excess_entropy,excess_volume]
# INPUT PROPORTIONS
-endmember_proportions = np.array([ 0.5, 0.2, 0.1, 0.2 ])
+molar_fraction = np.array([ 0.5, 0.2, 0.1, 0.2 ])
# "sites" is a 2D list of sites and the elements which reside on them
-# "site_occupancies" is a 3D list describing the elemental site occupancies of each endmember
-# "array_occupancies" is a 2D np.array of site_occupancies, concatenating the 2nd and 3rd dimension
+# "list_occupancies" is a 3D list describing the elemental site occupancies of each endmember
+# "site_occupancies" is a 2D np.array of list_occupancies, concatenating the 2nd and 3rd dimension
n_sites=base_material[0][1].count('[')
print 'Number of sites:', n_sites
print ''
sites=[[] for i in range(n_sites)]
-site_occupancies=[]
-site_multiplicity=np.empty(shape=(n_sites))
-n_elements=0
+list_occupancies=[]
+list_multiplicity=np.empty(shape=(n_sites))
+n_occupancies=0
for endmember in range(n_endmembers):
- site_occupancies.append([[0]*len(sites[site]) for site in range(n_sites)])
+ list_occupancies.append([[0]*len(sites[site]) for site in range(n_sites)])
s=re.split(r'\[', base_material[endmember][1])[1:]
for site in range(n_sites):
site_occupancy=re.split(r'\]', s[site])[0]
mult=re.split('[A-Z][^A-Z]*',re.split(r'\]', s[site])[1])[0]
if mult == '':
- site_multiplicity[site]=1.0
+ list_multiplicity[site]=1.0
else:
- site_multiplicity[site]=mult
+ list_multiplicity[site]=mult
elements=re.findall('[A-Z][^A-Z]*',site_occupancy)
for i in range(len(elements)):
element_on_site=re.split('[0-9][^A-Z]*',elements[i])[0]
@@ -83,30 +83,30 @@ for endmember in range(n_endmembers):
proportion_element_on_site=Fraction(proportion_element_on_site[0])
if element_on_site not in sites[site]:
- n_elements=n_elements+1
+ n_occupancies=n_occupancies+1
sites[site].append(element_on_site)
element_index=sites[site].index(element_on_site)
- for parsed_mbr in range(len(site_occupancies)):
- site_occupancies[parsed_mbr][site].append(0)
+ for parsed_mbr in range(len(list_occupancies)):
+ list_occupancies[parsed_mbr][site].append(0)
else:
element_index=sites[site].index(element_on_site)
- site_occupancies[endmember][site][element_index]=proportion_element_on_site
+ list_occupancies[endmember][site][element_index]=proportion_element_on_site
-array_occupancies=np.empty(shape=(n_endmembers,n_elements))
-array_multiplicities=np.empty(shape=(n_elements))
+site_occupancies=np.empty(shape=(n_endmembers,n_occupancies))
+site_multiplicities=np.empty(shape=(n_occupancies))
for endmember in range(n_endmembers):
n_element=0
for site in range(n_sites):
- for element in range(len(site_occupancies[endmember][site])):
- array_occupancies[endmember][n_element]=site_occupancies[endmember][site][element]
- array_multiplicities[n_element]=site_multiplicity[site]
+ for element in range(len(list_occupancies[endmember][site])):
+ site_occupancies[endmember][n_element]=list_occupancies[endmember][site][element]
+ site_multiplicities[n_element]=list_multiplicity[site]
n_element=n_element+1
# Matrix operation to return site occupancies, given proportions of endmembers
-occupancies=np.dot(endmember_proportions, array_occupancies)
+occupancies=np.dot(molar_fraction, site_occupancies)
print 'Site occupancies'
print sites
print occupancies
@@ -119,11 +119,11 @@ ideal_activity=np.empty(shape=(n_endmembers))
for endmember in range(n_endmembers):
ideal_activity[endmember]=1.0
normalisation_constant=1.0
- for element in range(n_elements):
- if array_occupancies[endmember][element] != 0:
- #print base_material[endmember][0], element, occupancies[element], array_occupancies[endmember][element]
- ideal_activity[endmember]=ideal_activity[endmember]*pow(occupancies[element],array_multiplicities[element])
- normalisation_constant=normalisation_constant/pow(array_occupancies[endmember][element],array_multiplicities[element])
+ for element in range(n_occupancies):
+ if site_occupancies[endmember][element] != 0:
+ #print base_material[endmember][0], element, occupancies[element], site_occupancies[endmember][element]
+ ideal_activity[endmember]=ideal_activity[endmember]*pow(occupancies[element],site_multiplicities[element])
+ normalisation_constant=normalisation_constant/pow(site_occupancies[endmember][element],site_multiplicities[element])
ideal_activity[endmember]=normalisation_constant*ideal_activity[endmember]
@@ -134,7 +134,7 @@ print ideal_activity
# alpha^T*p*(phi^T*W*phi)
alpha=np.array([base_material[i][2] for i in range(n_endmembers)])
-phi=np.array([alpha[i]*endmember_proportions[i] for i in range(n_endmembers)])
+phi=np.array([alpha[i]*molar_fraction[i] for i in range(n_endmembers)])
phi=np.divide(phi, np.sum(phi))
Wh=np.zeros(shape=(n_endmembers,n_endmembers))
@@ -148,25 +148,24 @@ for i in range(n_endmembers):
print ''
print 'Excess enthalpy, entropy, volume (non-configurational)'
-print np.dot(alpha.T,endmember_proportions)*np.dot(phi.T,np.dot(Wh,phi)), 'J/mol'
+print np.dot(alpha.T,molar_fraction)*np.dot(phi.T,np.dot(Wh,phi)), 'J/mol'
-print np.dot(alpha.T,endmember_proportions)*np.dot(phi.T,np.dot(Ws,phi)), 'J/K/mol'
+print np.dot(alpha.T,molar_fraction)*np.dot(phi.T,np.dot(Ws,phi)), 'J/K/mol'
-print np.dot(alpha.T,endmember_proportions)*np.dot(phi.T,np.dot(Wv,phi)), 'm^3/mol'
+print np.dot(alpha.T,molar_fraction)*np.dot(phi.T,np.dot(Wv,phi)), 'm^3/mol'
# Plot excess volumes for the pyrope-grossular join
-
ppy= np.linspace(0, 1, 101)
vex= np.empty(shape=(101))
for p in range(len(ppy)):
a=ppy[p]
- endmember_proportions = np.array([ a, 0.0, 1.-a, 0.0 ])
- phi=np.array([alpha[i]*endmember_proportions[i] for i in range(n_endmembers)])
+ molar_fraction = np.array([ a, 0.0, 1.-a, 0.0 ])
+ phi=np.array([alpha[i]*molar_fraction[i] for i in range(n_endmembers)])
phi=np.divide(phi, np.sum(phi))
- vex[p]=np.dot(alpha.T,endmember_proportions)*np.dot(phi.T,np.dot(Wv,phi))
+ vex[p]=np.dot(alpha.T,molar_fraction)*np.dot(phi.T,np.dot(Wv,phi))
import matplotlib.pyplot as plt
diff --git a/burnman/testing.py b/burnman/testing.py
deleted file mode 100644
index 7d31eb8..0000000
--- a/burnman/testing.py
+++ /dev/null
@@ -1,257 +0,0 @@
-# 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
-print 'Stishovite: a simple phase'
-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]]
-
-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 'RMS error on G:', 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 'Quartz: a phase governed by a phase transition described by Landau theory'
-
-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 'RMS error on G:', 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: iron
-iron=minerals.HP_2011.iron()
-iron.set_method("mtait")
-
-Giron=[['P \ T',25,500,1500],[0.001,-8.07,-28.25,-103.64],[10.0,-1.00,-21.05,-96.09],[100.0,60.89,41.85,-30.62]]
-
-print 'Iron: a phase governed by a phase transition described by Landau theory'
-
-Pmaxerror=0.
-Tmaxerror=0.
-maxerror=0.
-serror=0.
-for pi in range(len(Giron)-1):
- P=Giron[pi+1][0]*1e8
- for ti in range(len(Giron[0])-1):
- T=Giron[0][ti+1]+273.15
- Gburnman=iron.calcgibbs(P,T)
- GHP=Giron[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(Giron)-1)*(len(Giron[0])-1)))
-print 'RMS error on G:', rmserror, 'J/mol'
-print 'Maximum error:', maxerror, 'J/mol @', P/1.e9, 'GPa and', T, 'K'
-print ''
-
-# Another mineral with a phase transition described by Bragg-Williams theory: hercynite
-print 'Hercynite: a phase with order-disorder described by Bragg-Williams'
-herc=minerals.HP_2011.hercynite()
-herc.set_method("mtait")
-
-Gherc=[['P \ T',25, 500, 1500],[0.001,-1986.97,-2079.75,-2436.29],[10.0,-1946.33,-2038.65,-2394.06],[100.0,-1589.25,-1677.93,-2024.39]]
-
-print 'Scaling on configurational energy:', herc.params['BW_factor']
-
-Pmaxerror=0.
-Tmaxerror=0.
-maxerror=0.
-serror=0.
-for pi in range(len(Gherc)-1):
- P=Gherc[pi+1][0]*1e8
- for ti in range(len(Gherc[0])-1):
- T=Gherc[0][ti+1]+273.15
- Gburnman=herc.calcgibbs(P,T)
- GHP=Gherc[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(Gherc)-1)*(len(Gherc[0])-1)))
-print 'RMS error on G:', rmserror, 'J/mol'
-print 'Maximum error:', maxerror, 'J/mol @', P/1.e9, 'GPa and', T, 'K'
-print ''
-
-
-
-# Another mineral with a phase transition described by Bragg-Williams theory: dolomite
-print 'Dolomite: a phase with order-disorder described by Bragg-Williams'
-dol=minerals.HP_2011.dolomite()
-dol.set_method("mtait")
-
-Gdol=[['P \ T',25, 500, 1500],[0.001,-2372.73,-2496.11,-2964.31],[10.0,-2308.78,-2430.98,-2895.98],[100.0,-1759.31,-1872.97,-2315.19]]
-
-print 'Scaling on configurational energy:', dol.params['BW_factor']
-
-Pmaxerror=0.
-Tmaxerror=0.
-maxerror=0.
-serror=0.
-for pi in range(len(Gdol)-1):
- P=Gdol[pi+1][0]*1e8
- for ti in range(len(Gdol[0])-1):
- T=Gdol[0][ti+1]+273.15
- Gburnman=dol.calcgibbs(P,T)
- GHP=Gdol[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(Gdol)-1)*(len(Gdol[0])-1)))
-print 'RMS error on G:', 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
-print 'Spinel: a phase with order-disorder described by Bragg-Williams'
-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 'Scaling on configurational energy:', sp.params['BW_factor']
-
-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 'RMS error on G:', rmserror, 'J/mol'
-print 'Maximum error:', maxerror, 'J/mol @', P/1.e9, 'GPa and', T, 'K'
-print ''
-
-
-
-
-# Another mineral with a phase transition described by Bragg-Williams theory: sillimanite
-print 'Sillimanite: a phase with order-disorder described by Bragg-Williams'
-sill=minerals.HP_2011.sillimanite()
-sill.set_method("mtait")
-
-Gsill=[['P \ T',25, 500, 1500],[0.001,-2614.21,-2698.95,-3039.60],[10.0,-2564.51,-2648.92,-2988.73],[100.0,-2129.23,-2211.19,-2544.71]]
-
-print 'Scaling on configurational energy:', sill.params['BW_factor']
-
-Pmaxerror=0.
-Tmaxerror=0.
-maxerror=0.
-serror=0.
-for pi in range(len(Gsill)-1):
- P=Gsill[pi+1][0]*1e8
- for ti in range(len(Gsill[0])-1):
- T=Gsill[0][ti+1]+273.15
- Gburnman=sill.calcgibbs(P,T)
- GHP=Gsill[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(Gsill)-1)*(len(Gsill[0])-1)))
-print 'RMS error on G:', 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