Included the 4th order Birch murnaghan EOS. Added 2 Epsilon Fe minerals

in minerals.other and a liquid Fe EOS to use the 4th order BM EOS

burnman/eos/birch_murnaghan.py

@@ -21,7 +21,22 @@ def bulk_modulus(volume, params):
     K = pow(1. + 2.*f, 5./2.)* (params['K_0'] + (3. * params['K_0'] * params['Kprime_0'] -  \
            5*params['K_0'] ) * f + 27./2. * (params['K_0']*params['Kprime_0'] - 4.* params['K_0'])*f*f)
     return K
+def bulk_modulus_fourth(volume, params):
+    """
+    compute the bulk modulus as per the third order
+    birch-murnaghan equation of state.  Returns bulk
+    modulus in the same units as the reference bulk
+    modulus.  Pressure must be in :math:`[Pa]`.
+    """
+    B1 = (params['K_0']*params['Kprime_prime_0'])+((params['Kprime_0']-4.)*(params['Kprime_0']-5.))+(59./9.)
+    B2 = (3.*params['K_0']*params['Kprime_prime_0'])+((params['Kprime_0']-4.)*((3.*params['Kprime_0'])-13.))+(129./9.)
+    B3 = (3.*params['K_0']*params['Kprime_prime_0'])+((params['Kprime_0']-4.)*((3.*params['Kprime_0'])-11.))+(105./9.)
+    B4 = (params['K_0']*params['Kprime_prime_0'])+((params['Kprime_0']-4.)*(params['Kprime_0']-3.))+(35./9.)
 
+    x = params['V_0']/volume
+
+    K = (9./16.)*params['K_0']*( (-5.*B1/3.)*pow(x,-5./3.)+(7.*B2/3.)*pow(x,-7./3.)-3*B3*pow(x,-3.)+(11.*B4/3.)*pow(x,-11./3.))
+    return K
 
 def birch_murnaghan(x, params):
     """
@@ -44,6 +59,21 @@ def volume(pressure, params):
     V = opt.brentq(func, 0.5*params['V_0'], 1.5*params['V_0'])
     return V
 
+def volume_fourth_order(pressure,params):
+
+    func = lambda x: birch_murnaghan_fourth(x/params['V_0'], params) - pressure
+    #print "down",birch_murnaghan_fourth(.1*params['V_0']/params['V_0'], params) - pressure
+    #print "up",birch_murnaghan_fourth(200.*params['V_0']/params['V_0'], params) - pressure
+    V = opt.brentq(func, .1*params['V_0'], 1.5*params['V_0'])
+    return V
+
+def birch_murnaghan_fourth(x, params):
+    B1 = (params['K_0']*params['Kprime_prime_0'])+((params['Kprime_0']-4.)*(params['Kprime_0']-5.))+(59./9.)
+    B2 = (3.*params['K_0']*params['Kprime_prime_0'])+((params['Kprime_0']-4.)*((3.*params['Kprime_0'])-13.))+(129./9.)
+    B3 = (3.*params['K_0']*params['Kprime_prime_0'])+((params['Kprime_0']-4.)*((3.*params['Kprime_0'])-11.))+(105./9.)
+    B4 = (params['K_0']*params['Kprime_prime_0'])+((params['Kprime_0']-4.)*(params['Kprime_0']-3.))+(35./9.)
+
+    return (9.*params['K_0']/16.)*((-B1*pow(x, -5./3.)+(B2*pow(x, -7./3.))) - (B3*pow(x, -3.))+(B4*pow(x, -11./3.)))
 def shear_modulus_second_order(volume, params):
     """
     Get the birch murnaghan shear modulus at a reference temperature, for a
@@ -78,18 +108,26 @@ def volume(self,pressure, temperature, params):
         """
         Returns volume :math:`[m^3]` as a function of pressure :math:`[Pa]`.
         """
-        return volume(pressure,params)
+        if(self.order == 4):
+          return volume_fourth_order(pressure,params)
+        else:
+            return volume(pressure,params)
 
     def pressure(self, temperature, volume, params):
-        return birch_murnaghan(params['V_0']/volume, params)
+        if(self.order == 4):
+            return birch_murnaghan_fourth(volume/params['V_0'], params)
+        else:
+            return birch_murnaghan(params['V_0']/volume, params)
 
     def isothermal_bulk_modulus(self,pressure,temperature, volume, params):
         """
         Returns isothermal bulk modulus :math:`K_T` :math:`[Pa]` as a function of pressure :math:`[Pa]`,
         temperature :math:`[K]` and volume :math:`[m^3]`. 
         """
-        return bulk_modulus(volume, params)
-
+        if(self.order == 4):
+            return bulk_modulus_fourth(volume,params)
+        else:
+            return bulk_modulus(volume,params)
     def adiabatic_bulk_modulus(self,pressure, temperature, volume, params):
         """
         Returns adiabatic bulk modulus :math:`K_s` of the mineral. :math:`[Pa]`.
@@ -104,7 +142,8 @@ def shear_modulus(self,pressure, temperature, volume, params):
           return shear_modulus_second_order(volume,params)
         elif(self.order == 3):
           return shear_modulus_third_order(volume,params)
-
+        elif(self.order == 4):
+          return shear_modulus_third_order(volume,params)
     def heat_capacity_v(self,pressure, temperature, volume, params):
         """
         Since this equation of state does not contain temperature effects, simply return a very large number. :math:`[J/K/mol]`
@@ -183,3 +222,11 @@ class BM2(BirchMurnaghanBase):
     """
     def __init__(self):
         self.order=2
+
+class BM4(BirchMurnaghanBase):
+    """
+    Fourth order Birch Murnaghan isothermal equation of state.
+    This currently does not include shear modulus
+    """
+    def __init__(self):
+        self.order=4

burnman/eos/helper.py

@@ -22,8 +22,8 @@ def create(method):
     if isinstance(method, str):
         if method == "slb2":
             return slb.SLB2()
-        elif method == "V":
-            return v.Vinet()
+        elif method == "Vinet":
+            return vinet.Vinet()
         elif method == "mgd2":
             return mgd.MGD2()
         elif method == "mgd3":

burnman/eos/vinet.py

@@ -6,6 +6,7 @@
 import equation_of_state as eos
 import warnings
 from math import exp
+
 def bulk_modulus(volume, params):
     """
     compute the bulk modulus as per the third order
@@ -20,7 +21,6 @@ def bulk_modulus(volume, params):
     K = (params['K_0']*pow(x,-2./3.))*(1+((eta*pow(x,1./3.)+1.)*(1.-pow(x,1./3.))))*exp(eta*(1.-pow(x,1./3.)))
     return K
 
-
 def vinet(x, params):
     """
     equation for the third order Vinet equation of state, returns
@@ -31,13 +31,6 @@ def vinet(x, params):
     return 3.*params['K_0'] * (pow(x, -2./3.))*(1.-(pow(x,1./3.))) \
     * exp(eta*(1-pow(x,1/3.)))
 
-def density(pressure, params):
-    """
-    Get the Vinet density at a reference temperature for a given
-    pressure :math:`[Pa]`. Returns density in :math:`[kg/m^3]'
-    """
-    return params['molar_mass']/volume(pressure,params)
-
 def volume(pressure, params):
     """
     Get the Vinet volume at a reference temperature for a given
@@ -48,30 +41,6 @@ def volume(pressure, params):
     V = opt.brentq(func, 0.1*params['V_0'], 1.5*params['V_0'])
     return V
 
-def shear_modulus_second_order(volume, params):
-    """
-    Get the Vinet shear modulus at a reference temperature, for a
-    given volume.  Returns shear modulus in :math:`[Pa]` (the same units as in
-    params['G_0']).  This uses a second order finite strain expansion
-    """
-
-    x = params['V_0']/volume
-    G=params['G_0'] * pow(x,5./3.)*(1.-0.5*(pow(x,2./3.)-1.)*(5.-3.*params['Gprime_0']*params['K_0']/params['G_0']))
-    return G
-
-def shear_modulus_third_order(volume, params):
-    """
-    Get the Vinet shear modulus at a reference temperature, for a
-    given volume.  Returns shear modulus in :math:`[Pa]` (the same units as in
-    params['G_0']).  This uses a third order finite strain expansion
-    """
-
-    x = params['V_0']/volume
-    f = 0.5*(pow(x, 2./3.) - 1.0)
-    G = pow((1. + 2*f), 5./2.)*(params['G_0']+(3.*params['K_0']*params['Gprime_0'] - 5.*params['G_0'])*f + (6.*params['K_0']*params['Gprime_0']-24.*params['K_0']-14.*params['G_0']+9./2. * params['K_0']*params['Kprime_0'])*f*f)
-    return G
-
-
 class Vinet(eos.EquationOfState):
     """
     Base class for the isothermal Vinet equation of state.  This is third order in strain, and
@@ -134,17 +103,15 @@ def validate_parameters(self, params):
         """
         Check for existence and validity of the parameters
         """
-
-        #if G and Gprime are not included this is presumably deliberate,
-        #as we can model density and bulk modulus just fine without them,
-        #so just add them to the dictionary as nans
+
+        #G is not included in the Vinet EOS so we shall set them to NaN's
         if 'G_0' not in params:
             params['G_0'] = float('nan')
         if 'Gprime_0' not in params:
             params['Gprime_0'] = float('nan')
 
         #check that all the required keys are in the dictionary
-        expected_keys = ['V_0', 'K_0', 'Kprime_0', 'G_0', 'Gprime_0']
+        expected_keys = ['V_0', 'K_0', 'Kprime_0']
         for k in expected_keys:
             if k not in params:
                 raise KeyError('params object missing parameter : ' + k)
@@ -158,7 +125,3 @@ def validate_parameters(self, params):
             warnings.warn( 'Unusual value for K_0' , stacklevel=2)
         if params['Kprime_0'] < -5. or params['Kprime_0'] > 10.:
             warnings.warn( 'Unusual value for Kprime_0', stacklevel=2 )
-        if params['G_0'] < 0.0 or params['G_0'] > 1.e13:
-            warnings.warn( 'Unusual value for G_0', stacklevel=2 )
-        if params['Gprime_0'] < -5. or params['Gprime_0'] > 10.:
-            warnings.warn( 'Unusual value for Gprime_0', stacklevel=2 )

burnman/minerals/other.py

@@ -132,7 +132,44 @@ def __init__(self):
                         'q_0': 1.2  }
         Mineral.__init__(self)
 
+class Liquid_Fe_Anderson(Mineral):
+    """
+    Anderson & Ahrens, 1994 JGR
+    """
+    def __init__(self):
+        self.params = {
+            'equation_of_state': 'bm4',
+            'V_0': 7.95626e-6,
+            'K_0': 109.7e9,
+            'Kprime_0': 4.66,
+            'Kprime_prime_0': -0.043e-9,
+            'molar_mass': 0.055845,
+        }
+        Mineral.__init__(self)
+class Fe_Dewaele(Mineral):
+    """
+    Dewaele et al., 2006, Physical Review Letters
+    """
+    def __init__(self):
+        self.params = {
+                'equation_of_state': 'Vinet',
+                'V_0': 6.75e-6,
+                'K_0': 163.4e9,
+                'Kprime_0': 5.38,
+                'molar_mass': 0.055845,
+                'n': 1}
+        Mineral.__init__(self)
 
-
-
-
+class Fe_Anderson(Mineral):
+    """
+    Speziale et al. 2007, Mg#=83
+    """
+    def __init__(self):
+        self.params = {
+                'equation_of_state':'Vinet',
+                'V_0': (6.72e-6),
+                'K_0': 156.2e9,
+                'Kprime_0': 6.08,
+                'molar_mass': 0.055845,
+                'n': 1,} #number of atoms per formula unit
+        Mineral.__init__(self)