[cig-commits] [commit] inversion, master, validate_MT_params: Simplified configurational entropy calculation (ce35392)

[cig_noreply at geodynamics.org]

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

On branches: inversion,master,validate_MT_params
Link       : https://github.com/geodynamics/burnman/compare/80c2a295c42dfdb38f83f6c1334bf7d8f97a8463...409647ff05dfad6a686198cac1481bd46b5e2e62

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

commit ce353928cda56245c2f6f736b709378759760a60
Author: Bob Myhill <myhill.bob at gmail.com>
Date:   Wed Sep 3 23:01:07 2014 +0200

    Simplified configurational entropy calculation


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

ce353928cda56245c2f6f736b709378759760a60
 burnman/solutionmodel.py    | 22 ++++++++++------------
 solidsolution_benchmarks.py |  6 ++----
 2 files changed, 12 insertions(+), 16 deletions(-)

diff --git a/burnman/solutionmodel.py b/burnman/solutionmodel.py
index 4c67757..3544a36 100644
--- a/burnman/solutionmodel.py
+++ b/burnman/solutionmodel.py
@@ -102,21 +102,19 @@ class IdealSolution ( SolutionModel ):
         return self.ideal_gibbs_excess( temperature, molar_fractions )
 
     def configurational_entropy (self, molar_fractions):
-        activities = self.ideal_activities( molar_fractions )
-        conf_entropy = 0.0
+        site_occupancies=np.dot(molar_fractions, self.endmember_occupancies)
+        conf_entropy=0
+        for idx, occupancy in enumerate(site_occupancies):
+            if occupancy > 1e-10:
+                conf_entropy=conf_entropy-R*occupancy*self.site_multiplicities[idx]*np.log(occupancy)
 
-        tol = 1.e-10
-        for i in range(self.n_endmembers):
-            if molar_fractions[i] > tol:
-                conf_entropy = conf_entropy -  \
-                               molar_fractions[i] * R * np.log(activities[i])
         return conf_entropy
 
     def ideal_gibbs_excess( self, temperature, molar_fractions ): 
-        return -temperature*self.configurational_entropy(molar_fractions)
+        return 0.0-temperature*self.configurational_entropy(molar_fractions)
 
     def ideal_activities ( self, molar_fractions ):
-
+        # XXX Need to check this against configurational entropy!!!
         site_occupancies=np.dot(molar_fractions, self.endmember_occupancies)
         activities=np.empty(shape=(self.n_endmembers))
 
@@ -124,9 +122,9 @@ class IdealSolution ( SolutionModel ):
             activities[e]=1.0
             normalisation_constant=1.0
             for occ in range(self.n_occupancies):
-                if self.endmember_occupancies[e][occ] != 0: #integer?  
-                    activities[e]=activities[e]*np.power(site_occupancies[occ],self.site_multiplicities[occ])
-                    normalisation_constant=normalisation_constant/np.power(self.endmember_occupancies[e][occ],self.site_multiplicities[occ])
+                if self.endmember_occupancies[e][occ] != 0.: #integer?  
+                    activities[e]=activities[e]*np.power(site_occupancies[occ],site_occupancies[occ]*self.site_multiplicities[occ])
+                    normalisation_constant=normalisation_constant/np.power(self.endmember_occupancies[e][occ],self.endmember_occupancies[e][occ]*self.site_multiplicities[occ])
             activities[e]=normalisation_constant*activities[e]
 
         return activities
diff --git a/solidsolution_benchmarks.py b/solidsolution_benchmarks.py
index ea8bb85..eb57592 100644
--- a/solidsolution_benchmarks.py
+++ b/solidsolution_benchmarks.py
@@ -55,8 +55,7 @@ for i,c in enumerate(comp):
         molar_fractions=[1.0-c, c]
         sp.set_composition( np.array(molar_fractions) )
         sp.set_state( 1e5, 298.15 )
-        sp_entropies[i] = sp.solution_model.configurational_entropy( molar_fractions ) + \
-            (sum(molar_fractions[mbr]*sp.solution_model.endmember_configurational_entropies[mbr] for mbr in range(sp.solution_model.n_endmembers)))
+        sp_entropies[i] = sp.solution_model.configurational_entropy( molar_fractions )
         x=c/1.5
         if i > 0:
             sp_S[i] = -8.3145*(x*np.log(x) + (1.-x)*np.log(1.-x) + x*np.log(x/2.) + (2.-x)*np.log(1.-x/2.))
@@ -182,8 +181,7 @@ for idx, model in enumerate(opx_models):
         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_entropies[idx][i] = model.solution_model.configurational_entropy( molar_fractions )