Merge pull request geodynamics#144 from bobmyhill/activities

Added activity functions to solid solution, updated example document

burnman/solidsolution.py

@@ -10,8 +10,7 @@
 from .mineral import Mineral, material_property
 from .solutionmodel import *
 from .averaging_schemes import reuss_average_function
-
-
+from . import constants
 
 class SolidSolution(Mineral):
     """
@@ -122,17 +121,30 @@ def set_state(self, pressure, temperature):
         Mineral.set_state(self, pressure, temperature)
         for i in range(self.n_endmembers):
             self.endmembers[i][0].set_state(pressure, temperature)
+
 
+    @material_property
+    def activities(self):
+        """
+        Returns a list of endmember activities [unitless]
+        """
+        return self.solution_model.activities( self.pressure, self.temperature, self.molar_fractions)
 
+    @material_property
+    def activity_coefficients(self):
+        """
+        Returns a list of endmember activity coefficients (gamma = activity / ideal activity) [unitless]
+        """
+        return self.solution_model.activity_coefficients( self.pressure, self.temperature, self.molar_fractions)
 
+
     @material_property
     def internal_energy(self):
         """
         Returns internal energy of the mineral [J]
         Aliased with self.energy
         """
-        raise NotImplementedError("need to implement internal_energy() for solid solution!")
-        return None
+        return self.molar_helmholtz - self.pressure*self.molar_volume
 
 
     @material_property
@@ -152,8 +164,7 @@ def partial_gibbs(self):
         """
         return np.array([self.endmembers[i][0].gibbs for i in range(self.n_endmembers)]) + self.excess_partial_gibbs
 
-
-
+
     @material_property
     def excess_gibbs(self):
         """
@@ -178,8 +189,7 @@ def molar_helmholtz(self):
         Returns Helmholtz free energy of the solid solution [J]
         Aliased with self.helmholtz
         """
-        raise NotImplementedError("need to implement molar_helmholtz() for solid solution!")
-        return None
+        return self.molar_gibbs + self.temperature*self.molar_entropy
 
 
     @material_property
@@ -229,7 +239,7 @@ def excess_entropy(self):
     @material_property
     def molar_entropy(self):
         """
-        Returns enthalpy of the solid solution [J]
+        Returns entropy of the solid solution [J]
         Aliased with self.S
         """
         return sum([ self.endmembers[i][0].S * self.molar_fractions[i] for i in range(self.n_endmembers) ]) + self.excess_entropy
@@ -273,7 +283,7 @@ def adiabatic_bulk_modulus(self):
             return  self.isothermal_bulk_modulus
         else:
             return  self.isothermal_bulk_modulus*self.heat_capacity_p/self.heat_capacity_v
-
+        
 
     @material_property
     def isothermal_compressibility(self):
@@ -368,3 +378,4 @@ def heat_capacity_p(self):
         Aliased with self.C_p
         """
         return sum([ self.endmembers[i][0].heat_capacity_p * self.molar_fractions[i] for i in range(self.n_endmembers) ])
+

burnman/solutionmodel.py

@@ -229,6 +229,13 @@ def _ideal_activities (self, molar_fractions):
             activities[e]=normalisation_constant*activities[e]
         return activities
 
+    def activity_coefficients(self, pressure, temperature, molar_fractions):
+        return np.ones_like(molar_fractions)
+
+    def activities (self, pressure, temperature, molar_fractions):
+        return self._ideal_activities(molar_fractions)
+
+
 
 class AsymmetricRegularSolution (IdealSolution):
     """
@@ -314,7 +321,15 @@ def excess_enthalpy(self, pressure, temperature, molar_fractions):
         H_excess=np.dot(self.alpha.T,molar_fractions)*np.dot(phi.T,np.dot(self.Wh,phi))
         return H_excess + pressure*self.excess_volume ( pressure, temperature, molar_fractions )
 
-
+    def activity_coefficients(self, pressure, temperature, molar_fractions):
+        if temperature > 1.e-10:
+            return np.exp(self._non_ideal_excess_partial_gibbs(pressure, temperature, molar_fractions)/(constants.gas_constant*temperature))
+        else:
+            raise Exception("Activity coefficients not defined at 0 K.")
+
+    def activities(self, pressure, temperature, molar_fractions):
+        return IdealSolution.activities(self, pressure, temperature, molar_fractions) * self.activity_coefficients(pressure, temperature, molar_fractions)
+
 class SymmetricRegularSolution (AsymmetricRegularSolution):
     """
     Solution model implementing the symmetric regular solution model
@@ -401,3 +416,13 @@ def excess_entropy(self, pressure, temperature, molar_fractions):
     def excess_enthalpy(self, pressure, temperature, molar_fractions):
         H_excess=np.dot(molar_fractions, self._non_ideal_function(self.Wh, molar_fractions))
         return H_excess + pressure*self.excess_volume (pressure, temperature, molar_fractions)
+
+    def activity_coefficients(self, pressure, temperature, molar_fractions):
+        if temperature > 1.e-10:
+            return np.exp(self._non_ideal_excess_partial_gibbs(pressure, temperature, molar_fractions)/(constants.gas_constant*temperature))
+        else:
+            raise Exception("Activity coefficients not defined at 0 K.")
+
+    def activities(self, pressure, temperature, molar_fractions):
+        return IdealSolution.activities(self, pressure, temperature, molar_fractions) * self.activity_coefficients(pressure, temperature, molar_fractions)
+

examples/example_solid_solution.py

@@ -87,18 +87,38 @@ def __init__(self, molar_fractions=None):
     comp = np.linspace(0.001, 0.999, 100)
     g1_gibbs = np.empty_like(comp)
     g1_excess_gibbs = np.empty_like(comp)
+    g1_pyrope_activity = np.empty_like(comp)
+    g1_almandine_activity = np.empty_like(comp)
+    g1_pyrope_gamma = np.empty_like(comp)
+    g1_almandine_gamma = np.empty_like(comp)
     for i,c in enumerate(comp):
         molar_fractions=[1.0-c, c]
         g1.set_composition(molar_fractions)
         g1.set_state(P,T)
         g1_gibbs[i] = g1.gibbs
         g1_excess_gibbs[i] = g1.excess_gibbs
+        g1_pyrope_activity[i] = g1.activities[0]
+        g1_almandine_activity[i] = g1.activities[1]
+        g1_pyrope_gamma[i] = g1.activity_coefficients[0]
+        g1_almandine_gamma[i] = g1.activity_coefficients[1]
 
     plt.plot( [0., 1.], [py_gibbs/1000., alm_gibbs/1000.], 'b-', linewidth=1., label='Mechanical mixing')
     plt.plot( comp, g1_gibbs/1000., 'r-', linewidth=1., label='Ideal solid solution')
-    plt.title("Pyrope-almandine join")
+    plt.title("Ideal pyrope-almandine join")
     plt.ylabel("Gibbs free energy of solution (kJ/mol)")
-    plt.xlabel("Pyrope fraction")
+    plt.xlabel("Almandine fraction")
+    plt.legend(loc='lower right')
+    plt.show()
+
+    plt.plot( comp, g1_pyrope_activity, 'r-', linewidth=1., label='Pyrope activity')
+    plt.plot( comp, g1_almandine_activity, 'b-', linewidth=1., label='Almandine activity')
+
+    plt.plot( comp, g1_pyrope_gamma, 'r-.', linewidth=1., label='Pyrope activity coefficient')
+    plt.plot( comp, g1_almandine_gamma, 'b--', linewidth=1., label='Almandine activity coefficient')
+    plt.title("Ideal pyrope-almandine join")
+    plt.ylim(0.0, 1.01)
+    plt.ylabel("Activities")
+    plt.xlabel("Almandine fraction")
     plt.legend(loc='lower right')
     plt.show()
 
@@ -128,22 +148,44 @@ def __init__(self, molar_fractions=None):
 
             burnman.SolidSolution.__init__(self, molar_fractions)
 
+
     g2=mg_fe_garnet_2()
     g2_excess_gibbs = np.empty_like(comp)
+    g2_pyrope_activity = np.empty_like(comp)
+    g2_almandine_activity = np.empty_like(comp)
+    g2_pyrope_gamma = np.empty_like(comp)
+    g2_almandine_gamma = np.empty_like(comp)
     for i,c in enumerate(comp):
-        molar_fractions=[1.0-c, c]
+        molar_fractions=[1.-c, c]
         g2.set_composition(molar_fractions)
         g2.set_state(P,T)
         g2_excess_gibbs[i] = g2.excess_gibbs
+        g2_pyrope_activity[i] = g2.activities[0]
+        g2_almandine_activity[i] = g2.activities[1]
+        g2_pyrope_gamma[i] = g2.activity_coefficients[0]
+        g2_almandine_gamma[i] = g2.activity_coefficients[1]
 
     plt.plot( comp, g1_excess_gibbs, 'r-', linewidth=1., label='Ideal solution')
     plt.plot( comp, g2_excess_gibbs, 'g-', linewidth=1., label='Symmetric solution, 2.5 kJ/mol')
     plt.title("Pyrope-almandine join (model comparison)")
     plt.ylabel("Excess gibbs free energy of solution (J/mol)")
-    plt.xlabel("Pyrope fraction")
+    plt.xlabel("Almandine fraction")
     plt.legend(loc='upper left')
     plt.show()
 
+
+    plt.plot( comp, g2_pyrope_activity, 'r-', linewidth=1., label='Pyrope activity')
+    plt.plot( comp, g2_almandine_activity, 'b-', linewidth=1., label='Almandine activity')
+
+    plt.plot( comp, g2_pyrope_gamma, 'r-.', linewidth=1., label='Pyrope activity coefficient')
+    plt.plot( comp, g2_almandine_gamma, 'b--', linewidth=1., label='Almandine activity coefficient')
+    plt.title("Non-ideal pyrope-almandine join")
+    plt.ylabel("Activities")
+    plt.xlabel("Almandine fraction")
+    plt.legend(loc='lower right')
+    plt.show()
+
+
     '''
     Adding more endmembers is very straightforward.
     Interaction terms must be added in the order [[12,13,14,...],[23,24,...],[34,...],...]
@@ -187,7 +229,7 @@ def __init__(self, molar_fractions=None):
     plt.plot( comp, g3_excess_entropy, 'r-', linewidth=1., label='Excess entropy')
     plt.title("Pyrope-majorite join")
     plt.ylabel("Entropy (J/K/mol)")
-    plt.xlabel("Pyrope fraction")
+    plt.xlabel("Majorite fraction")
     plt.legend(loc='upper left')
     plt.show()
 
@@ -231,7 +273,7 @@ def __init__(self, molar_fractions=None):
     plt.plot( comp, g4_excess_gibbs_1200, 'b-', linewidth=1., label='1200 K')
     plt.title("Pyrope-grossular join (asymmetric model)")
     plt.ylabel("Excess gibbs free energy of solution (J/mol)")
-    plt.xlabel("Pyrope fraction")
+    plt.xlabel("Grossular fraction")
     plt.legend(loc='lower left')
     plt.show()
 
@@ -260,17 +302,19 @@ def __init__(self, molar_fractions=None):
             burnman.SolidSolution.__init__(self, molar_fractions)
 
     g5=mg_fe_ca_garnet_Ganguly()
-    g5_excess_enthalpy = np.empty_like(comp)
-
+    g5_excess_gibbs = np.empty_like(comp)
     for i,c in enumerate(comp):
-        molar_fractions=[1.0-c, 0., c, 0.]
+        molar_fractions=[1.-c, 0., c, 0.]
         g5.set_composition(molar_fractions)
         g5.set_state(1.e5,298.15)
-        g5_excess_enthalpy[i] = g5.excess_enthalpy
+        g5_excess_gibbs[i] = g5.excess_gibbs
+
+    plt.plot( comp, g5_excess_gibbs/3., 'r-', linewidth=1., label='Py-Gr excess enthalpy (J/cation-mole)')
 
-    plt.plot( comp, g5_excess_enthalpy/3., 'r-', linewidth=1., label='Py-Gr excess enthalpy (J/cation-mole)')
     plt.title("Asymmetric py-gr join (Ganguly et al., 1996; Figure 5)")
     plt.ylabel("Excess enthalpy of solution (J/cation-mol)")
-    plt.xlabel("Pyrope fraction")
+    plt.xlabel("XCa")
     plt.legend(loc='lower left')
     plt.show()
+
+

tests/test_solidsolution.py

@@ -70,6 +70,15 @@ def __init__(self, molar_fractions=None):
 
         burnman.SolidSolution.__init__(self, molar_fractions)
 
+# Ideal solid solution
+class olivine_ideal_ss(burnman.SolidSolution):
+    def __init__(self, molar_fractions=None):
+        self.name='Fo-Fo solid solution'
+        self.type='ideal'
+        self.endmembers = [[forsterite(), '[Mg]2SiO4'], [fayalite(), '[Fe]2SiO4']]
+
+        burnman.SolidSolution.__init__(self, molar_fractions)
+
 # Olivine solid solution
 class olivine_ss(burnman.SolidSolution):
     def __init__(self, molar_fractions=None):
@@ -86,8 +95,8 @@ def __init__(self, molar_fractions=None):
         # Name
         self.name='orthopyroxene'
         self.type='symmetric'
-        self.endmembers = [[forsterite(), 'Mg[Mg]Si2O6'], [forsterite(), '[Mg1/2Al1/2]2AlSiO6']]
-        self.enthalpy_interaction=[[10.0e3]]
+        self.endmembers = [[forsterite(), '[Mg][Mg]Si2O6'], [forsterite(), '[Mg1/2Al1/2][Mg1/2Al1/2]AlSiO6']]
+        self.enthalpy_interaction=[[burnman.constants.gas_constant*1.0e3]]
 
         burnman.SolidSolution.__init__(self, molar_fractions)
 
@@ -217,6 +226,24 @@ def test_subregular(self):
 
         self.assertArraysAlmostEqual(ss0.excess_partial_gibbs, ss1.excess_partial_gibbs)
 
+    def test_activities_ideal(self):
+        ol = olivine_ideal_ss()
+        ol.set_composition( np.array([0.5, 0.5]) )
+        ol.set_state(1.e5,1000.)
+        self.assertArraysAlmostEqual(ol.activities, [0.25, 0.25])
 
+
+    def test_activity_coefficients_ideal(self):
+        ol = olivine_ideal_ss()
+        ol.set_composition( np.array([0.5, 0.5]) )
+        ol.set_state(1.e5,1000.)
+        self.assertArraysAlmostEqual(ol.activity_coefficients, [1., 1.])
+
+    def test_activity_coefficients_non_ideal(self):
+        opx = orthopyroxene()
+        opx.set_composition( np.array([0.0, 1.0]) )
+        opx.set_state(1.e5,1000.)
+        self.assertArraysAlmostEqual(opx.activity_coefficients, [np.exp(1.), 1.])
+
 if __name__ == '__main__':
     unittest.main()