[cig-commits] [commit] add_thermodynamic_potentials: Implemented endmember configurational entropies for disordered phases (a4dc8be)

[cig_noreply at geodynamics.org]

Repository : https://github.com/geodynamics/burnman

On branch  : add_thermodynamic_potentials
Link       : https://github.com/geodynamics/burnman/compare/2e5646d93cedbbf8eae54cc37cffc14e0aa85180...d5ddad03ff9f30f5a4efaddb4e3ec585ea1a7c51

>---------------------------------------------------------------

commit a4dc8be8baffc837113f4a5a61f8eeab5db03274
Author: Bob Myhill <myhill.bob at gmail.com>
Date:   Wed Sep 3 20:05:37 2014 +0200

    Implemented endmember configurational entropies for disordered phases


>---------------------------------------------------------------

a4dc8be8baffc837113f4a5a61f8eeab5db03274
 burnman/modified_tait.py    |  2 +-
 burnman/solutionmodel.py    | 15 ++++++++++--
 solidsolution_benchmarks.py | 57 ++++++++++++++++++++++++++++++++-------------
 3 files changed, 55 insertions(+), 19 deletions(-)

diff --git a/burnman/modified_tait.py b/burnman/modified_tait.py
index 36e1635..0049726 100644
--- a/burnman/modified_tait.py
+++ b/burnman/modified_tait.py
@@ -13,7 +13,7 @@ import burnman.einstein as einstein
 
 T_0=298.15 # Standard temperature = 25 C
 P_0=1.e5 # Standard pressure = 1.e5 Pa
-R=8.3145 # J/K/mol
+R=8.31446 # J/K/mol
 
 # see Holland and Powell, 2011
 def einstein_temperature(S, n):
diff --git a/burnman/solutionmodel.py b/burnman/solutionmodel.py
index 7f6ada7..4c67757 100644
--- a/burnman/solutionmodel.py
+++ b/burnman/solutionmodel.py
@@ -7,7 +7,7 @@ import warnings
 import burnman
 from burnman.processchemistry import *
 
-R = 8.3145 # J/K/mol
+R = 8.31446 # J/K/mol
 
 class SolutionModel:
     """
@@ -87,6 +87,17 @@ class IdealSolution ( SolutionModel ):
         self.solution_formulae, self.n_sites, self.sites, self.n_occupancies, self.endmember_occupancies, self.site_multiplicities = \
             ProcessSolidSolutionChemistry(self.formulas)
 
+        self.endmember_configurational_entropies=np.zeros(shape=(self.n_endmembers))
+        for idx, endmember_occupancy in enumerate(self.endmember_occupancies):
+            S_conf=0.
+            for occ in range(self.n_occupancies):
+                if endmember_occupancy[occ] != 0: #integer? 
+                    self.endmember_configurational_entropies[idx] = \
+                        self.endmember_configurational_entropies[idx] - \
+                        R*self.site_multiplicities[occ]*endmember_occupancy[occ]*np.log(endmember_occupancy[occ])
+
+
+
     def excess_gibbs_free_energy( self, pressure, temperature, molar_fractions ):
         return self.ideal_gibbs_excess( temperature, molar_fractions )
 
@@ -101,7 +112,7 @@ class IdealSolution ( SolutionModel ):
                                molar_fractions[i] * R * np.log(activities[i])
         return conf_entropy
 
-    def  ideal_gibbs_excess( self, temperature, molar_fractions ): 
+    def ideal_gibbs_excess( self, temperature, molar_fractions ): 
         return -temperature*self.configurational_entropy(molar_fractions)
 
     def ideal_activities ( self, molar_fractions ):
diff --git a/solidsolution_benchmarks.py b/solidsolution_benchmarks.py
index 4b27be6..21c8917 100644
--- a/solidsolution_benchmarks.py
+++ b/solidsolution_benchmarks.py
@@ -100,30 +100,55 @@ class orthopyroxene_blue(burnman.SolidSolution):
         burnman.SolidSolution.__init__(self, base_material, \
                           burnman.solutionmodel.SymmetricRegularSolution(base_material, enthalpy_interaction) )
 
-opx_red = orthopyroxene_red()
-opx_red.set_method('slb3')
-opx_blue = orthopyroxene_blue()
-opx_blue.set_method('slb3')
+class orthopyroxene_long_dashed(burnman.SolidSolution):
+    def __init__(self):
+        # Name
+        self.name='orthopyroxene'
+
+        # Endmembers (cpx is symmetric)
+        base_material = [[enstatite(), 'Mg[Mg]Si2O6'],[mg_tschermaks_molecule(), '[Mg1/2Al1/2]2AlSiO6'] ]
+
+        # Interaction parameters
+        enthalpy_interaction=[[0.0]]
+
+        burnman.SolidSolution.__init__(self, base_material, \
+                          burnman.solutionmodel.SymmetricRegularSolution(base_material, enthalpy_interaction) )
+
+class orthopyroxene_short_dashed(burnman.SolidSolution):
+    def __init__(self):
+        # Name
+        self.name='orthopyroxene'
+
+        # Endmembers (cpx is symmetric)
+        base_material = [[enstatite(), 'Mg[Mg][Si]2O6'],[mg_tschermaks_molecule(), 'Mg[Al][Al1/2Si1/2]2O6'] ]
+
+        # Interaction parameters
+        enthalpy_interaction=[[0.0]]
+
+        burnman.SolidSolution.__init__(self, base_material, \
+                          burnman.solutionmodel.SymmetricRegularSolution(base_material, enthalpy_interaction) )
 
 comp = np.linspace(0, 1.0, 100)
-entropy_red = np.empty_like(comp)
-entropy_blue = np.empty_like(comp)
+opx_models=[orthopyroxene_red(), orthopyroxene_blue(), orthopyroxene_long_dashed(), orthopyroxene_short_dashed()]
+opx_entropies = [ np.empty_like(comp) for model in opx_models ]
+for idx, model in enumerate(opx_models):
+    model.set_method('slb3')
 
-for i,c in enumerate(comp):
-   opx_red.set_composition( np.array([1.0-c, c]) )
-   opx_red.set_state( 0.0, 0.0 )
-   entropy_red[i] = opx_red.solution_model.configurational_entropy( [1.0-c, c] )
+    for i,c in enumerate(comp):
+        molar_fractions=[1.0-c, c]
+        model.set_composition( np.array(molar_fractions) )
+        model.set_state( 0.0, 0.0 )
+        opx_entropies[idx][i] = model.solution_model.configurational_entropy( molar_fractions ) + \
+            (sum(molar_fractions[mbr]*model.solution_model.endmember_configurational_entropies[mbr] for mbr in range(model.solution_model.n_endmembers)))
 
-   opx_blue.set_composition( np.array([1.0-c, c]) )
-   opx_blue.set_state( 0.0, 0.0 )
-   entropy_blue[i] = opx_blue.solution_model.configurational_entropy( [1.0-c, c] )
 
 
 #fig1 = mpimg.imread('configurational_entropy.png')  # Uncomment these two lines if you want to overlay the plot on a screengrab from SLB2011
 #plt.imshow(fig1, extent=[0.0, 1.0,0.,17.0], aspect='auto')
-
-plt.plot( comp, entropy_red, 'r--', linewidth=3.)
-plt.plot( comp, entropy_blue, 'b--', linewidth=3.)
+plt.plot( comp, opx_entropies[0], 'r--', linewidth=3.)
+plt.plot( comp, opx_entropies[1], 'b--', linewidth=3.)
+plt.plot( comp, opx_entropies[2], 'g--', linewidth=3.)
+plt.plot( comp, opx_entropies[3], 'g-.', linewidth=3.)
 plt.xlim(0.0,1.0)
 plt.ylim(0.,17.0)
 plt.ylabel("Configurational enthalpy of solution")