[cig-commits] [commit] add_gibbs_energy: Disorganized work on getting dGdN stuff working (eb34fce)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Dec 11 17:13:02 PST 2014
Repository : https://github.com/geodynamics/burnman
On branch : add_gibbs_energy
Link : https://github.com/geodynamics/burnman/compare/0000000000000000000000000000000000000000...2148b324d3e8aa7b527f831eb397590942563008
>---------------------------------------------------------------
commit eb34fcee312687f6eb7fc227ffd03b2cf67a8442
Author: Ian Rose <ian.r.rose at gmail.com>
Date: Wed Oct 22 11:44:25 2014 -0700
Disorganized work on getting dGdN stuff working
Conflicts:
burnman/solidsolution.py
>---------------------------------------------------------------
eb34fcee312687f6eb7fc227ffd03b2cf67a8442
burnman/solidsolution.py | 10 ++---
burnman/solutionmodel.py | 89 +++++++++++++++++++++++++--------------------
solidsolution_benchmarks.py | 4 +-
3 files changed, 57 insertions(+), 46 deletions(-)
diff --git a/burnman/solidsolution.py b/burnman/solidsolution.py
index e03676c..3c488bc 100644
--- a/burnman/solidsolution.py
+++ b/burnman/solidsolution.py
@@ -52,11 +52,14 @@ class SolidSolution(Mineral):
for i in range(self.n_endmembers):
self.base_material[i][0].set_state(pressure, temperature)
+ self.excess_partial_gibbs = self.solution_model.excess_partial_gibbs_free_energies( pressure, temperature, self.molar_fraction)
self.excess_gibbs = self.solution_model.excess_gibbs_free_energy( pressure, temperature, self.molar_fraction)
+ self.partial_gibbs = np.array([self.base_material[i][0].gibbs for i in range(self.n_endmembers)]) + self.excess_partial_gibbs
+ self.gibbs= sum([ self.base_material[i][0].gibbs * self.molar_fraction[i] for i in range(self.n_endmembers) ]) + self.excess_gibbs
+
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.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
'''
@@ -74,7 +77,4 @@ class SolidSolution(Mineral):
return sum([ self.base_material[i][0].calcgibbs(pressure, temperature) * molar_fractions[i] for i in range(self.n_endmembers) ]) + self.solution_model.excess_gibbs_free_energy( pressure, temperature, molar_fractions)
def calcpartialgibbsexcesses(self, pressure, temperature, molar_fractions):
- Hint, Sint, Vint = self.solution_model.non_ideal_interactions(molar_fractions)
- partialgibbsexcesses=np.empty(len(molar_fractions))
- partialgibbsexcesses = np.array([0.+R*temperature*self.solution_model.ln_ideal_activities(molar_fractions)[i] + Hint[i] - temperature*Sint[i] + pressure*Vint[i] for i in range(self.n_endmembers)])
- return partialgibbsexcesses
+ return self.solution_model.excess_partial_gibbs_free_energies( pressure, temperature, molar_fractions)
diff --git a/burnman/solutionmodel.py b/burnman/solutionmodel.py
index ad97188..a812b49 100644
--- a/burnman/solutionmodel.py
+++ b/burnman/solutionmodel.py
@@ -25,8 +25,9 @@ class SolutionModel:
def excess_gibbs_free_energy( self, pressure, temperature, molar_fractions):
"""
Given a list of molar fractions of different phases,
- compute the excess Gibbs free energy of the solution,
- relative to an ideal model.
+ compute the excess Gibbs free energy of the solution.
+ The base class implementation assumes that the excess gibbs
+ free energy is zero.
Parameters
----------
@@ -44,14 +45,38 @@ class SolutionModel:
G_excess : float
The excess Gibbs free energy
"""
- return 0.0
+ return np.dot(np.array(molar_fractions), self.excess_partial_gibbs_free_energies( pressure, temperature, molar_fractions))
+
+ def excess_partial_gibbs_free_energies( self, pressure, temperature, molar_fractions):
+ """
+ Given a list of molar fractions of different phases,
+ compute the excess Gibbs free energy for each endmember of the solution.
+ The base class implementation assumes that the excess gibbs
+ free energy is zero.
+
+ Parameters
+ ----------
+ pressure : float
+ Pressure at which to evaluate the solution model. [Pa]
+
+ temperature : float
+ Temperature at which to evaluate the solution. [K]
+
+ molar_fractions : list of floats
+ List of molar fractions of the different endmembers in solution
+
+ Returns
+ -------
+ partial_G_excess : numpy array
+ The excess Gibbs free energy of each endmember
+ """
+ return np.empty_like( np.array(molar_fractions) )
def excess_volume( self, pressure, temperature, molar_fractions):
"""
Given a list of molar fractions of different phases,
compute the excess Gibbs free energy of the solution,
- relative to an ideal model. The base class implementation
- assumes that the excess volume is zero.
+ The base class implementation assumes that the excess volume is zero.
Parameters
----------
@@ -87,6 +112,12 @@ class IdealSolution ( SolutionModel ):
self.solution_formulae, self.n_sites, self.sites, self.n_occupancies, self.endmember_occupancies, self.site_multiplicities = \
ProcessSolidSolutionChemistry(self.formulas)
+ self._calculate_endmember_configurational_entropies()
+
+ def excess_partial_gibbs_free_energies( self, pressure, temperature, molar_fractions ):
+ return self._ideal_excess_partial_gibbs( temperature, molar_fractions )
+
+ def _calculate_endmember_configurational_entropies( self ):
self.endmember_configurational_entropies=np.zeros(shape=(self.n_endmembers))
for idx, endmember_occupancy in enumerate(self.endmember_occupancies):
for occ in range(self.n_occupancies):
@@ -95,34 +126,23 @@ class IdealSolution ( SolutionModel ):
self.endmember_configurational_entropies[idx] - \
R*self.site_multiplicities[occ]*endmember_occupancy[occ]*np.log(endmember_occupancy[occ])
+ def _endmember_configurational_entropy_contribution(self, molar_fractions):
+ return np.dot(molar_fractions, self.endmember_configurational_entropies)
-
- def excess_gibbs_free_energy( self, pressure, temperature, molar_fractions ):
- return self.ideal_gibbs_excess( temperature, molar_fractions )
-
- def configurational_entropy (self, molar_fractions):
+ def _configurational_entropy (self, molar_fractions):
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)
- # 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 endmember_configurational_entropy_contribution(self, molar_fractions):
- return np.dot(molar_fractions, self.endmember_configurational_entropies)
- def ideal_gibbs_excess( self, temperature, molar_fractions ):
- return 0.0-temperature*(self.configurational_entropy(molar_fractions) - self.endmember_configurational_entropy_contribution(molar_fractions))
+ def _ideal_excess_partial_gibbs( self, temperature, molar_fractions ):
+ return R * temperature * self._log_ideal_activities(molar_fractions)
- def ln_ideal_activities ( self, molar_fractions ):
+ def _log_ideal_activities ( self, molar_fractions ):
site_occupancies=np.dot(molar_fractions, self.endmember_occupancies)
lna=np.empty(shape=(self.n_endmembers))
@@ -137,7 +157,7 @@ class IdealSolution ( SolutionModel ):
return lna
- def ideal_activities ( self, molar_fractions ):
+ def _ideal_activities ( self, molar_fractions ):
site_occupancies=np.dot(molar_fractions, self.endmember_occupancies)
activities=np.empty(shape=(self.n_endmembers))
@@ -186,10 +206,8 @@ class AsymmetricRegularSolution ( IdealSolution ):
#initialize ideal solution model
IdealSolution.__init__(self, endmembers )
- def configurational_entropy( self, molar_fractions ):
- return IdealSolution.configurational_entropy( self, molar_fractions )
- def non_ideal_interactions( self, molar_fractions ):
+ def _non_ideal_interactions( self, molar_fractions ):
# -sum(sum(qi.qj.Wij*)
# equation (2) of Holland and Powell 2003
phi=np.array([self.alpha[i]*molar_fractions[i] for i in range(self.n_endmembers)])
@@ -210,24 +228,17 @@ class AsymmetricRegularSolution ( IdealSolution ):
return Hint, Sint, Vint
+ def _non_ideal_excess_partial_gibbs( self, pressure, temperature, molar_fractions) :
- def excess_gibbs_free_energy( self, pressure, temperature, molar_fractions ):
+ Hint, Sint, Vint = self._non_ideal_interactions( molar_fractions )
+ return Hint - temperature*Sint + pressure*Vint
- ideal_gibbs = IdealSolution.ideal_gibbs_excess( self, temperature, molar_fractions )
-
-
- phi=np.array([self.alpha[i]*molar_fractions[i] for i in range(self.n_endmembers)])
- phi=np.divide(phi, np.sum(phi))
-
- H_excess=np.dot(self.alpha.T,molar_fractions)*np.dot(phi.T,np.dot(self.Wh,phi))
- S_excess=np.dot(self.alpha.T,molar_fractions)*np.dot(phi.T,np.dot(self.Ws,phi))
- V_excess=np.dot(self.alpha.T,molar_fractions)*np.dot(phi.T,np.dot(self.Wv,phi))
-
- non_ideal_gibbs = H_excess - temperature*S_excess + pressure*V_excess
+ def excess_partial_gibbs_free_energies( self, pressure, temperature, molar_fractions ):
+ ideal_gibbs = IdealSolution._ideal_excess_partial_gibbs (self, temperature, molar_fractions )
+ non_ideal_gibbs = self._non_ideal_excess_partial_gibbs( pressure, temperature, molar_fractions)
return ideal_gibbs + non_ideal_gibbs
-
def excess_volume ( self, pressure, temperature, molar_fractions ):
phi=np.array([self.alpha[i]*molar_fractions[i] for i in range(self.n_endmembers)])
diff --git a/solidsolution_benchmarks.py b/solidsolution_benchmarks.py
index 92fef5d..b1a16e1 100644
--- a/solidsolution_benchmarks.py
+++ b/solidsolution_benchmarks.py
@@ -55,7 +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 )
+ sp_entropies[i] = sp.solution_model._configurational_entropy( molar_fractions )
sp_entropies_NK1967[i] = -8.3145*(c*np.log(c) + (1.-c)*np.log(1.-c) + c*np.log(c/2.) + (2.-c)*np.log(1.-c/2.)) # eq. 7 in Navrotsky and Kleppa, 1967.
#fig1 = mpimg.imread('configurational_entropy.png') # Uncomment these two lines if you want to overlay the plot on a screengrab from SLB2011
@@ -177,7 +177,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 )
+ opx_entropies[idx][i] = model.solution_model._configurational_entropy( molar_fractions )
More information about the CIG-COMMITS
mailing list