[cig-commits] [commit] add_thermodynamic_potentials: Initial version of SolutionModel class, including one that does nothing, and IdealSolution class, a SymmetricVanLaar class, and an AsymmetricVanLaar class. Still work to be done (f3b6439)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue Dec 9 09:55:05 PST 2014
Repository : https://github.com/geodynamics/burnman
On branch : add_thermodynamic_potentials
Link : https://github.com/geodynamics/burnman/compare/2e5646d93cedbbf8eae54cc37cffc14e0aa85180...d5ddad03ff9f30f5a4efaddb4e3ec585ea1a7c51
>---------------------------------------------------------------
commit f3b6439d5221f4023d4a4804b98b4fe1e219f8ca
Author: ian-r-rose <ian.r.rose at gmail.com>
Date: Tue Sep 2 11:57:35 2014 -0700
Initial version of SolutionModel class, including one that does nothing, and IdealSolution class, a SymmetricVanLaar class, and an AsymmetricVanLaar class. Still work to be done
>---------------------------------------------------------------
f3b6439d5221f4023d4a4804b98b4fe1e219f8ca
burnman/minerals/HP_2011_ds62.py | 10 +--
burnman/minerals/SLB_2011_w_solid_solutions.py | 79 +++++++++-----------
burnman/solidsolution.py | 99 +++-----------------------
burnman/test_solidsolution.py | 22 +++---
4 files changed, 62 insertions(+), 148 deletions(-)
diff --git a/burnman/minerals/HP_2011_ds62.py b/burnman/minerals/HP_2011_ds62.py
index 436768d..20a6354 100644
--- a/burnman/minerals/HP_2011_ds62.py
+++ b/burnman/minerals/HP_2011_ds62.py
@@ -18,6 +18,7 @@ N.B. VERY IMPORTANT: The excess entropy term in the regular solution model has t
from burnman.mineral import Mineral
from burnman.solidsolution import SolidSolution
from burnman.processchemistry import read_masses, dictionarize_formula, formula_mass
+from burnman.solutionmodel import *
atomic_masses=read_masses('data/input_masses/atomic_masses.dat')
@@ -28,15 +29,14 @@ class garnet(SolidSolution):
self.name='garnet'
# Endmembers
- base_material = [[py(), '[Mg]3[Al]2Si3O12', 1.0],[alm(), '[Fe]3[Al]2Si3O12', 1.0],[gr(), '[Ca]3[Al]2Si3O12', 2.7],[maj(), '[Mg]3[Mg1/2Si1/2]2Si3O12', 1.0]]
-
- # Interaction parameters
+ endmembers = [[py(), '[Mg]3[Al]2Si3O12'], [alm(), '[Fe]3[Al]2Si3O12'], [gr(), '[Ca]3[Al]2Si3O12'], [maj(), '[Mg]3[Mg1/2Si1/2]2Si3O12']]
+ alphas = [1.0, 1.0, 2.7, 1.0]
excess_enthalpy=[[2.5e3, 29.1e3, 15e3],[10e3,18e3],[48e3]]
excess_entropy=[[0., 0., 0.],[0., 0.],[0.]]
excess_volume=[[0., 0.164e-5, 0.],[0., 0.],[0.]]
- interaction_parameter=[excess_enthalpy,excess_entropy,excess_volume]
- SolidSolution.__init__(self, base_material, interaction_parameter)
+ sm = AsymmetricVanLaar( endmembers, alphas, excess_enthalpy, excess_volume, excess_entropy)
+ SolidSolution.__init__(self, endmembers, sm)
class fo (Mineral):
diff --git a/burnman/minerals/SLB_2011_w_solid_solutions.py b/burnman/minerals/SLB_2011_w_solid_solutions.py
index 9e68e73..48aeb46 100644
--- a/burnman/minerals/SLB_2011_w_solid_solutions.py
+++ b/burnman/minerals/SLB_2011_w_solid_solutions.py
@@ -14,6 +14,7 @@ Solid solutions from inv251010 of HeFESTo
from burnman.mineral import Mineral
from burnman.solidsolution import SolidSolution
+from burnman.solutionmodel import *
from burnman.processchemistry import read_masses, dictionarize_formula, formula_mass
atomic_masses=read_masses('data/input_masses/atomic_masses.dat')
@@ -30,10 +31,7 @@ class c2c_pyroxene(SolidSolution):
# Endmembers (C2/c is symmetric)
base_material = [[hp_clinoenstatite(), '[Mg]2Si2O6'],[hp_clinoferrosilite(), '[Fe]2Si2O6']]
- # Interaction parameters (C2c is ideal)
- interaction_parameter=[]
-
- SolidSolution.__init__(self, base_material, interaction_parameter)
+ SolidSolution.__init__(self, base_material, IdealSolution(base_material))
@@ -45,10 +43,7 @@ class ca_ferrite_structured_phase(SolidSolution):
# Endmembers (CF is symmetric)
base_material = [[mg_calcium_ferrite(), '[Mg]Al[Al]O4'],[fe_calcium_ferrite(), '[Fe]Al[Al]O4'],[na_calcium_ferrite(), '[Na]Al[Si]O4']]
- # Interaction parameters (CF is ideal)
- interaction_parameter=[]
-
- SolidSolution.__init__(self, base_material, interaction_parameter)
+ SolidSolution.__init__(self, base_material, IdealSolution(base_material))
class ca_rich_perovskite(SolidSolution):
def __init__(self):
@@ -58,10 +53,7 @@ class ca_rich_perovskite(SolidSolution):
# Endmembers (cpv is symmetric)
base_material = [[ca_perovskite(), '[Ca][Ca]Si2O6'],[jd_perovskite(), '[Na][Al]Si2O6']]
- # Interaction parameters (cpv is ideal)
- interaction_parameter=[]
-
- SolidSolution.__init__(self, base_material, interaction_parameter)
+ SolidSolution.__init__(self, base_material, IdealSolution(base_materia))
class clinopyroxene(SolidSolution):
def __init__(self):
@@ -72,9 +64,9 @@ class clinopyroxene(SolidSolution):
base_material = [[diopside(), '[Ca][Mg][Si]2O6'],[hedenbergite(), '[Ca][Fe][Si]2O6'],[clinoenstatite(), '[Mg][Mg][Si]2O6'],[ca_tschermaks_molecule(), '[Ca][Al][Si1/2Al1/2]2O6'],[jadeite(), '[Na][Al][Si]2O6']]
# Interaction parameters
- interaction_parameter=[[0., 24.74e3, 26.e3, 24.3e3],[24.74e3, 0., 0.e3], [60.53136e3, 0.0], [10.e3]]
+ enthalpy_interaction=[[0., 24.74e3, 26.e3, 24.3e3],[24.74e3, 0., 0.e3], [60.53136e3, 0.0], [10.e3]]
- SolidSolution.__init__(self, base_material, interaction_parameter)
+ SolidSolution.__init__(self, base_material, SymmetricVanLaar(base_material, enthalpy_interaction) )
class garnet(SolidSolution):
def __init__(self):
@@ -84,9 +76,9 @@ class garnet(SolidSolution):
# Endmembers (garnet is symmetric)
base_material = [[pyrope(), '[Mg]3[Al][Al]Si3O12'],[almandine(), '[Fe]3[Al][Al]Si3O12'],[grossular(), '[Ca]3[Al][Al]Si3O12'],[mg_majorite(), '[Mg]3[Mg][Si]Si3O12'],[jd_majorite(), '[Na2/3Al1/3]3[Al][Si]Si3O12']]
# Interaction parameters
- interaction_parameter=[[0.0, 30.e3, 21.20278e3, 0.0],[0.0,0.0,0.0],[57.77596e3, 0.0],[0.0]]
+ enthalply_interaction=[[0.0, 30.e3, 21.20278e3, 0.0],[0.0,0.0,0.0],[57.77596e3, 0.0],[0.0]]
- SolidSolution.__init__(self, base_material, interaction_parameter)
+ SolidSolution.__init__(self, base_material, SymmetricVanLaar(base_material, enthalpy_interaction) )
class akimotoite(SolidSolution):
@@ -97,9 +89,9 @@ class akimotoite(SolidSolution):
# Endmembers (ilmenite/akimotoite is symmetric)
base_material = [[mg_akimotoite(), '[Mg][Si]O3'],[fe_akimotoite(), '[Fe][Si]O3'],[corundum(), '[Al][Al]O3']]
# Interaction parameters
- interaction_parameter=[[0.0, 66.e3],[66.e3]]
+ enthalpy_interaction=[[0.0, 66.e3],[66.e3]]
- SolidSolution.__init__(self, base_material, interaction_parameter)
+ SolidSolution.__init__(self, base_material, SymmetricVanLaar(base_material, enthalpy_interaction) )
class ferropericlase(SolidSolution):
def __init__(self):
@@ -109,9 +101,9 @@ class ferropericlase(SolidSolution):
# Endmembers (ferropericlase is symmetric)
base_material = [[periclase(), '[Mg]O'],[wuestite(), '[Fe]O']]
# Interaction parameters
- interaction_parameter=[[13.e3]]
+ enthalpy_interaction=[[13.e3]]
- SolidSolution.__init__(self, base_material, interaction_parameter)
+ SolidSolution.__init__(self, base_material, SymmetricVanLaar(base_material, enthalpy_interaction) )
class mg_fe_olivine(SolidSolution):
def __init__(self):
@@ -121,9 +113,9 @@ class mg_fe_olivine(SolidSolution):
# Endmembers (olivine is symmetric)
base_material = [[forsterite(), '[Mg]2SiO4'],[fayalite(), '[Fe]2SiO4']]
# Interaction parameters
- interaction_parameter=[[7.81322e3]]
+ enthalpy_interaction=[[7.81322e3]]
- SolidSolution.__init__(self, base_material, interaction_parameter)
+ SolidSolution.__init__(self, base_material, SymmetricVanLaar(base_material, enthalpy_interaction) )
class orthopyroxene(SolidSolution):
def __init__(self):
@@ -134,9 +126,9 @@ class orthopyroxene(SolidSolution):
base_material = [[enstatite(), '[Mg][Mg][Si]SiO6'],[ferrosilite(), '[Fe][Fe][Si]SiO6'],[mg_tschermaks_molecule(), '[Mg][Al][Al]SiO6'],[orthodiopside(), '[Ca][Mg][Si]SiO6']]
# Interaction parameters
- interaction_parameter=[[0.0, 0.0, 32.11352e3],[0.0, 0.0],[48.35316e3]]
+ enthalpy_interaction=[[0.0, 0.0, 32.11352e3],[0.0, 0.0],[48.35316e3]]
- SolidSolution.__init__(self, base_material, interaction_parameter)
+ SolidSolution.__init__(self, base_material, SymmetricVanLaar(base_material, enthalpy_interaction))
class plagioclase(SolidSolution):
def __init__(self):
@@ -146,9 +138,9 @@ class plagioclase(SolidSolution):
# Endmembers (plagioclase is symmetric)
base_material = [[anorthite(), '[Ca][Al]2Si2O8'],[albite(), '[Na][Al1/2Si1/2]2Si2O8']]
# Interaction parameters
- interaction_parameter=[[26.0e3]]
+ enthalpy_interaction=[[26.0e3]]
- SolidSolution.__init__(self, base_material, interaction_parameter)
+ SolidSolution.__init__(self, base_material, SymmetricVanLaar(base_material, enthalpy_interaction))
class post_perovskite(SolidSolution):
def __init__(self):
@@ -159,9 +151,9 @@ class post_perovskite(SolidSolution):
base_material = [[mg_post_perovskite(), '[Mg][Si]O3'],[fe_post_perovskite(), '[Fe][Si]O3'],[al_post_perovskite(), '[Al][Al]O3']]
# Interaction parameters
- interaction_parameter=[[0.0, 60.0e3],[0.0]]
+ enthalpy_interaction=[[0.0, 60.0e3],[0.0]]
- SolidSolution.__init__(self, base_material, interaction_parameter)
+ SolidSolution.__init__(self, base_material, SymmetricVanLaar(base_material, enthalpy_interaction))
class mg_fe_perovskite(SolidSolution):
def __init__(self):
@@ -172,9 +164,9 @@ class mg_fe_perovskite(SolidSolution):
base_material = [[mg_perovskite(), '[Mg][Si]O3'],[fe_perovskite(), '[Fe][Si]O3'],[al_perovskite(), '[Al][Al]O3']]
# Interaction parameters
- interaction_parameter=[[0.0, 116.0e3],[0.0]]
+ enthalpy_interaction=[[0.0, 116.0e3],[0.0]]
- SolidSolution.__init__(self, base_material, interaction_parameter)
+ SolidSolution.__init__(self, base_material, SymmetricVanLaar(base_material, enthalpy_interaction))
class mg_fe_ringwoodite(SolidSolution):
def __init__(self):
@@ -185,9 +177,9 @@ class mg_fe_ringwoodite(SolidSolution):
base_material = [[mg_ringwoodite(), '[Mg]2SiO4'],[fe_ringwoodite(), '[Fe]2SiO4']]
# Interaction parameters
- interaction_parameter=[[9.34084e3]]
+ enthalpy_interaction=[[9.34084e3]]
- SolidSolution.__init__(self, base_material, interaction_parameter)
+ SolidSolution.__init__(self, base_material, SymmetricVanLaar(base_material, enthalpy_interaction))
class mg_fe_aluminous_spinel(SolidSolution):
def __init__(self):
@@ -198,9 +190,9 @@ class mg_fe_aluminous_spinel(SolidSolution):
base_material = [[spinel(), '[Mg3/4Al1/4]4[Al7/8Mg1/8]8O16'],[hercynite(), '[Fe3/4Al1/4]4[Al7/8Fe1/8]8O16']]
# Interaction parameters
- interaction_parameter=[[5.87646e3]]
+ enthalpy_interaction=[[5.87646e3]]
- SolidSolution.__init__(self, base_material, interaction_parameter)
+ SolidSolution.__init__(self, base_material, SymmetricVanLaar(base_material, enthalpy_interaction))
class mg_fe_wadsleyite(SolidSolution):
def __init__(self):
@@ -211,10 +203,9 @@ class mg_fe_wadsleyite(SolidSolution):
base_material = [[mg_wadsleyite(), '[Mg]2SiO4'],[fe_wadsleyite(), '[Fe]2SiO4']]
# Interaction parameters
- interaction_parameter=[[16.74718e3]]
-
- SolidSolution.__init__(self, base_material, interaction_parameter)
+ enthalpy_interaction=[[16.74718e3]]
+ SolidSolution.__init__(self, base_material, SymmetricVanLaar(base_material, enthalpy_interaction) )
'''
ENDMEMBERS
@@ -240,7 +231,7 @@ class anorthite (Mineral):
'Gprime_0': 1.1,
'eta_s_0': 1.6}
- self.uncertainties = {
+ self.uncertainties = {
'err_F_0': 4.e3,
'err_V_0': 0.,
'err_K_0': 5.e9,
@@ -274,7 +265,7 @@ class albite (Mineral):
'Gprime_0': 1.4,
'eta_s_0': 1.0}
- self.uncertainties = {
+ self.uncertainties = {
'err_F_0': 5.e3,
'err_V_0': 0.0,
'err_K_0': 5.e9,
@@ -308,7 +299,7 @@ class spinel (Mineral):
'Gprime_0': 0.4,
'eta_s_0': 2.7}
- self.uncertainties = {
+ self.uncertainties = {
'err_F_0': 32.e3,
'err_V_0': 0.0,
'err_K_0': 1.e9,
@@ -342,7 +333,7 @@ class hercynite (Mineral):
'Gprime_0': 0.4,
'eta_s_0': 2.8}
- self.uncertainties = {
+ self.uncertainties = {
'err_F_0': 35.e3,
'err_V_0': 0.0,
'err_K_0': 2.e9,
@@ -376,7 +367,7 @@ class forsterite (Mineral):
'Gprime_0': 1.5,
'eta_s_0': 2.3}
- self.uncertainties = {
+ self.uncertainties = {
'err_F_0': 2.e3,
'err_V_0': 0.0,
'err_K_0': 2.e9,
@@ -410,7 +401,7 @@ class fayalite (Mineral):
'Gprime_0': 1.5,
'eta_s_0': 1.0}
- self.uncertainties = {
+ self.uncertainties = {
'err_F_0': 1.e3,
'err_V_0': 0.0,
'err_K_0': 2.e9,
@@ -423,7 +414,6 @@ class fayalite (Mineral):
'err_eta_s_0': 0.6
}
-'''
class mg_wadsleyite (Mineral):
def __init__(self):
formula=''
@@ -1818,7 +1808,6 @@ class nepheline (Mineral):
'err_eta_s_0':
}
-'''
# Mineral aliases
# Feldspars
diff --git a/burnman/solidsolution.py b/burnman/solidsolution.py
index 8828394..f24757d 100644
--- a/burnman/solidsolution.py
+++ b/burnman/solidsolution.py
@@ -5,6 +5,7 @@
import numpy as np
from burnman.mineral import Mineral
from burnman.processchemistry import ProcessSolidSolutionChemistry
+from burnman.solutionmodel import SolutionModel
import warnings
R = 8.3145 # J/K/mol
@@ -27,110 +28,32 @@ class SolidSolution(Mineral):
"""
# init sets up matrices to speed up calculations for when P, T, X is defined.
- def __init__(self, base_material, interaction_parameter):
+ def __init__(self, base_material, solution_model = SolutionModel() ):
# Initialise the solid solution inputs
self.base_material = base_material
- self.interaction_parameter = interaction_parameter
+ self.solution_model = solution_model
# 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]) == self.n_endmembers-1)
- for j in range(self.n_endmembers-1):
- assert(len(interaction_parameter[i][j]) == (self.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]
-
- # Process solid solution chemistry
- self.solution_formulae, self.n_sites, self.sites, self.n_occupancies, self.endmember_occupancies, self.site_multiplicities = ProcessSolidSolutionChemistry([base_material[i][1] for i in range(self.n_endmembers)])
-
- # Check the chemical composition of each endmember is consistent with the dataset
- for i in range(self.n_endmembers):
- endmember_formula=self.base_material[i][0].params['formula']
- solution_formula=self.solution_formulae[i]
- if endmember_formula != solution_formula:
- print 'Formula of endmember', self.base_material[i][0].params['name'], 'does not agree with formula in the', self.name, 'SolidSolution model'
- print endmember_formula, self.solution_formulae[i]
- exit()
-
- # Create array of van Laar parameters
- self.alpha=np.array([base_material[i][2] for i in range(self.n_endmembers)])
-
- # Create 2D arrays of interaction parameters
- self.Wh=np.zeros(shape=(self.n_endmembers,self.n_endmembers))
- self.Ws=np.zeros(shape=(self.n_endmembers,self.n_endmembers))
- self.Wv=np.zeros(shape=(self.n_endmembers,self.n_endmembers))
- for i in range(self.n_endmembers):
- for j in range(i+1, self.n_endmembers):
- self.Wh[i][j]=2.*self.excess_enthalpy[i][j-i-1]/(self.alpha[i]+self.alpha[j])
- self.Ws[i][j]=2.*self.excess_entropy[i][j-i-1]/(self.alpha[i]+self.alpha[j])
- self.Wv[i][j]=2.*self.excess_volume[i][j-i-1]/(self.alpha[i]+self.alpha[j])
-
- # END INITIALISATION
-
-
- def set_composition(self, molar_fraction):
- # Check that the composition declared is consistent with the solution model
+ 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)
+ self.molar_fraction = molar_fraction
- # Make molar fraction an attribute of the solid solution
- self.molar_fraction=molar_fraction
-
- # Ideal activities
- self.site_occupancies=np.dot(self.molar_fraction, self.endmember_occupancies)
- self.ideal_activity=np.empty(shape=(self.n_endmembers))
- for endmember in range(self.n_endmembers):
- self.ideal_activity[endmember]=1.0
- normalisation_constant=1.0
- for element in range(self.n_occupancies):
- if self.endmember_occupancies[endmember][element] != 0:
- self.ideal_activity[endmember]=self.ideal_activity[endmember]*pow(self.site_occupancies[element],self.site_multiplicities[element])
- normalisation_constant=normalisation_constant/pow(self.endmember_occupancies[endmember][element],self.site_multiplicities[element])
- self.ideal_activity[endmember]=normalisation_constant*self.ideal_activity[endmember]
-
- # Nonideal contributions
- phi=np.array([self.alpha[i]*molar_fraction[i] for i in range(self.n_endmembers)])
- phi=np.divide(phi, np.sum(phi))
- self.H_excess=np.dot(self.alpha.T,molar_fraction)*np.dot(phi.T,np.dot(self.Wh,phi))
- self.S_excess=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))
+ def set_state(self, pressure, temperature):
- 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].set_method(self.base_material[i][0].params['equation_of_state'])
+ 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)
-
- # Ideal contribution (configurational entropy)
- # Don't include endmembers with negligable molar fractions as the log term blows up
- tol=1e-10
- self.gibbs_excess_ideal=0
- for i in range(self.n_endmembers):
- if molar_fraction[i] > tol:
- self.gibbs_excess_ideal=self.gibbs_excess_ideal + molar_fraction[i]*R*temperature*np.log(self.ideal_activity[i])
-
- # Non-ideal contribution
- self.gibbs_excess_nonideal=self.H_excess - temperature*self.S_excess + pressure*self.V_excess
-
- # Total excess Gibbs
- self.gibbs_excess=self.gibbs_excess_ideal + self.gibbs_excess_nonideal
+ self.excess_gibbs = self.solution_model.excess_gibbs_free_energy( pressure, temperature, self.molar_fraction)
+ self.excess_volume = self.solution_model.excess_volume( pressure, temperature, self.molar_fraction)
- self.V= sum([ self.base_material[i][0].V * self.molar_fraction[i] for i in range(self.n_endmembers) ]) + self.V_excess
- self.gibbs= sum([ self.base_material[i][0].gibbs * self.molar_fraction[i] for i in range(self.n_endmembers) ]) + self.gibbs_excess
+ self.V= sum([ self.base_material[i][0].V * self.molar_fraction[i] for i in range(self.n_endmembers) ]) + self.excess_volume
+ self.gibbs= sum([ self.base_material[i][0].gibbs * self.molar_fraction[i] for i in range(self.n_endmembers) ]) + self.excess_gibbs
'''
diff --git a/burnman/test_solidsolution.py b/burnman/test_solidsolution.py
index 6159d4a..45aa289 100644
--- a/burnman/test_solidsolution.py
+++ b/burnman/test_solidsolution.py
@@ -13,20 +13,19 @@ from burnman import minerals
garnet=minerals.HP_2011_ds62.garnet()
composition=np.array([ 0.5, 0.2, 0.1, 0.2 ])
+garnet.set_method('mtait')
garnet.set_composition(composition)
+garnet.set_state(0.0, 300.0)
print 'Molar fraction'
print garnet.molar_fraction
print ''
print 'Site occupancies'
-print garnet.sites
-print garnet.site_occupancies
-print ''
-print 'Ideal activities'
-print garnet.ideal_activity
+print garnet.solution_model.sites
+#print garnet.solution_model.site_occupancies
print ''
print 'Volume excess'
-print garnet.V_excess, 'm^3/mol'
+print garnet.excess_volume, 'm^3/mol'
print ''
# Excess volumes for the pyrope-grossular join
@@ -37,12 +36,14 @@ for i in range(n+1):
pyrope_proportion[i]=float(i)/n
composition=([ pyrope_proportion[i], 0.0, 1.-pyrope_proportion[i], 0.0 ])
garnet.set_composition(composition)
- garnet_excess_volume[i]=garnet.V_excess
+ garnet.set_state(0.0, 300.)
+ garnet_excess_volume[i]=garnet.excess_volume
pressure=1.e9 # Pa
temperature=573.15 # K
composition=[0.9, 0.0, 0.1, 0.0]
-garnet.set_state(pressure, temperature, composition)
+garnet.set_composition(composition)
+garnet.set_state(pressure, temperature)
# Excess gibbs for the pyrope-grossular join
n=100
@@ -51,8 +52,9 @@ garnet_excess_gibbs= np.empty(shape=(n+1))
for i in range(n+1):
pyrope_proportion[i]=float(i)/n
composition=([ pyrope_proportion[i], 0.0, 1.-pyrope_proportion[i], 0.0 ])
- garnet.set_state(pressure, temperature, composition)
- garnet_excess_gibbs[i]=garnet.gibbs_excess
+ garnet.set_composition(composition)
+ garnet.set_state(pressure, temperature)
+ garnet_excess_gibbs[i]=garnet.excess_gibbs
import matplotlib.pyplot as plt
More information about the CIG-COMMITS
mailing list