[cig-commits] [commit] inversion, master, validate_MT_params: Finished ideal activity formulation (7929428)

[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 7929428af3021bf6318fa451f8eec4d272161dec
Author: Bob Myhill <myhill.bob at gmail.com>
Date:   Wed Sep 3 23:50:00 2014 +0200

    Finished ideal activity formulation


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

7929428af3021bf6318fa451f8eec4d272161dec
 burnman/solutionmodel.py | 16 ++++++++++------
 1 file changed, 10 insertions(+), 6 deletions(-)

diff --git a/burnman/solutionmodel.py b/burnman/solutionmodel.py
index 3544a36..93022b5 100644
--- a/burnman/solutionmodel.py
+++ b/burnman/solutionmodel.py
@@ -89,9 +89,8 @@ class IdealSolution ( SolutionModel ):
 
         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? 
+                if endmember_occupancy[occ] != 0.: 
                     self.endmember_configurational_entropies[idx] = \
                         self.endmember_configurational_entropies[idx] - \
                         R*self.site_multiplicities[occ]*endmember_occupancy[occ]*np.log(endmember_occupancy[occ])
@@ -108,13 +107,19 @@ class IdealSolution ( SolutionModel ):
             if occupancy > 1e-10:
                 conf_entropy=conf_entropy-R*occupancy*self.site_multiplicities[idx]*np.log(occupancy)
 
+        # Alternative entropy construction
+        #activities=self.ideal_activities(molar_fractions)
+        #conf_entropy=0.
+        #for i, activity in enumerate(activities):
+        #    if activity > 1e-10:
+        #        conf_entropy=conf_entropy - R*molar_fractions[i]*np.log(activity) + molar_fractions[i]*self.endmember_configurational_entropies[i]
+
         return conf_entropy
 
     def ideal_gibbs_excess( self, temperature, 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))
 
@@ -122,11 +127,10 @@ 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],site_occupancies[occ]*self.site_multiplicities[occ])
+                if self.endmember_occupancies[e][occ] > 1e-10:
+                    activities[e]=activities[e]*np.power(site_occupancies[occ],self.endmember_occupancies[e][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