Merge pull request geodynamics#604 from bobmyhill/update_calibrants

Update calibrants

burnman/calibrants/Decker_1971.py

@@ -33,7 +33,7 @@ def _pressure_Decker_NaCl(volume, temperature, params):
                 + 9 * params["Kprime_0"] ** 2
                 - 63 * params["Kprime_0"]
                 + 143
-            )
+            ) / 6.0
             f = 0.5 * ((volume / params["V_0"]) ** (-2 / 3) - 1)
             K_T = (
                 params["K_0"]
@@ -55,7 +55,7 @@ def _pressure_Decker_NaCl(volume, temperature, params):
         _params_Decker_NaCl = {
             "V_0": 2.7013e-05,
             "K_0": 23.7,
-            "Kprime_0": 5.04,  # 4.91 in Matsui (2012), however 5.04 is required to reproduce values in Table 4.
+            "Kprime_0": 4.91,
             "Kprime_prime_0": -0.267,
             "Debye_0": 279.0,
             "grueneisen_0": 1.59,

burnman/calibrants/Fei_2007.py

@@ -7,7 +7,7 @@
 from burnman.eos.birch_murnaghan import BirchMurnaghanBase as BM3
 from burnman.eos.mie_grueneisen_debye import MGDBase
 from burnman.classes.calibrant import Calibrant
-
+from burnman.utils.unitcell import molar_volume_from_unit_cell_volume
 
 """
 Fei_2007
@@ -33,22 +33,14 @@ def _pressure_Fei_Pt(volume, temperature, params):
             Pth0 = thermal_model._thermal_pressure(params["T_0"], volume, params)
             Pth = thermal_model._thermal_pressure(temperature, volume, params)
 
-            # Electronic pressure
-            Pel = (
-                1.1916e-15 * temperature**4.0
-                - 1.4551e-11 * temperature**3.0
-                + 1.6209e-07 * temperature**2.0
-                + 1.8269e-4 * temperature
-                - 0.069
-            ) * 1.0e09
-
             # Total pressure
-            P = P0 + Pth - Pth0 + Pel
+            P = P0 + Pth - Pth0
 
             return P
 
+        Z = 4.0
         _params_Fei_Pt = {
-            "V_0": 9.0904e-06,
+            "V_0": molar_volume_from_unit_cell_volume(60.38, Z),
             "K_0": 277.0e9,
             "Kprime_0": 5.08,
             "Debye_0": 230.0,
@@ -57,7 +49,47 @@ def _pressure_Fei_Pt(volume, temperature, params):
             "n": 1.0,
             "T_0": 300.0,
             "P_0": 0.0,
-            "Z": 4.0,
+            "Z": Z,
         }
 
         Calibrant.__init__(self, _pressure_Fei_Pt, "pressure", _params_Fei_Pt)
+
+
+class Au(Calibrant):
+    """
+    The Au pressure standard reported by
+    Fei et al. (2007; https://doi.org/10.1073/pnas.0609013104).
+    """
+
+    def __init__(self):
+        def _pressure_Fei_Au(volume, temperature, params):
+
+            # Isothermal pressure (GPa)
+            pressure_model = Vinet()
+            P0 = pressure_model.pressure(params["T_0"], volume, params)
+
+            # Thermal pressure
+            thermal_model = MGDBase()
+            Pth0 = thermal_model._thermal_pressure(params["T_0"], volume, params)
+            Pth = thermal_model._thermal_pressure(temperature, volume, params)
+
+            # Total pressure
+            P = P0 + Pth - Pth0
+
+            return P
+
+        Z = 4.0
+        _params_Fei_Au = {
+            "V_0": molar_volume_from_unit_cell_volume(67.850, Z),
+            "K_0": 167.0e9,
+            "Kprime_0": 6.0,
+            "Debye_0": 170.0,
+            "grueneisen_0": 2.97,
+            "q_0": 0.6,
+            "n": 1.0,
+            "T_0": 300.0,
+            "P_0": 0.0,
+            "Z": Z,
+        }
+
+        Calibrant.__init__(self, _pressure_Fei_Au, "pressure", _params_Fei_Au)

burnman/calibrants/Holmes_1989.py

@@ -0,0 +1,45 @@
+# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit for
+# the Earth and Planetary Sciences
+# Copyright (C) 2012 - 2024 by the BurnMan team, released under the GNU
+# GPL v2 or later.
+
+"""
+Holmes_1989
+^^^^^^^^^^^
+"""
+
+from burnman.classes.calibrant import Calibrant
+import numpy as np
+from burnman.utils.unitcell import molar_volume_from_unit_cell_volume
+
+
+class Pt(Calibrant):
+    """
+    The Pt pressure standard reported by
+    Holmes et al. (1989; https://doi.org/10.1063/1.344177).
+    """
+
+    def __init__(self):
+        def _pressure(volume, temperature, params):
+            X = np.power(volume / params["V_0"], 1.0 / 3.0)
+            P_300 = (
+                3.0
+                * params["beta_T"]
+                * (1.0 - X)
+                / (X * X)
+                * np.exp(params["eta"] * (1.0 - X))
+            )
+
+            return P_300 + params["alpha_T"] * params["beta_T"] * (temperature - 300.0)
+
+        Z = 4.0
+        _params = {
+            "V_0": molar_volume_from_unit_cell_volume(60.38, Z),
+            "beta_T": 798.31e9 / 3.0,
+            "eta": 7.2119,
+            "beta_prime_T": (7.2119 / 1.5) + 1.0,
+            "alpha_T": 2.61e-5,
+            "Z": Z,
+        }
+
+        Calibrant.__init__(self, _pressure, "pressure", _params)

burnman/calibrants/Matsui_2012.py

@@ -3,7 +3,6 @@
 # Copyright (C) 2012 - 2024 by the BurnMan team, released under the GNU
 # GPL v2 or later.
 
-from burnman.eos.vinet import Vinet
 from burnman.eos.mie_grueneisen_debye import MGDBase
 from burnman.classes.calibrant import Calibrant
 
@@ -30,7 +29,7 @@ def _pressure_Matsui_NaCl(volume, temperature, params):
                 + 9 * params["Kprime_0"] ** 2
                 - 63 * params["Kprime_0"]
                 + 143
-            )
+            ) / 6.0
             f = 0.5 * ((volume / params["V_0"]) ** (-2 / 3) - 1)
             K_T = (
                 params["K_0"]

burnman/calibrants/Tsuchiya_2003.py

@@ -0,0 +1,61 @@
+# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit for
+# the Earth and Planetary Sciences
+# Copyright (C) 2012 - 2024 by the BurnMan team, released under the GNU
+# GPL v2 or later.
+
+import numpy as np
+from burnman.classes.calibrant import Calibrant
+from scipy.interpolate import RegularGridInterpolator
+
+"""
+Tsuchiya_2003
+^^^^^^^^^^^^^
+"""
+
+
+class Au(Calibrant):
+    """
+    The Au pressure standard reported by
+    Tsuchiya (2003; https://doi.org/10.1029/2003JB002446).
+    """
+
+    def __init__(self):
+
+        grid_compressions = np.linspace(0.0, 0.34, 18)
+        grid_temperatures = np.array([300.0, 500.0, 1000.0, 1500.0, 2000.0, 2500.0])
+        grid_pressures = np.array(
+            [
+                [0.00, 1.52, 5.35, 9.19, 13.04, 16.88],
+                [3.55, 5.04, 8.78, 12.54, 16.29, 20.05],
+                [7.68, 9.13, 12.79, 16.45, 20.12, 23.79],
+                [12.42, 13.83, 17.40, 20.98, 24.56, 28.14],
+                [17.86, 19.23, 22.71, 26.20, 29.70, 33.19],
+                [24.12, 25.46, 28.85, 32.25, 35.66, 39.07],
+                [31.30, 32.60, 35.90, 39.22, 42.54, 45.86],
+                [39.52, 40.78, 43.99, 47.22, 50.45, 53.68],
+                [48.94, 50.17, 53.29, 56.43, 59.58, 62.72],
+                [59.76, 60.95, 63.98, 67.03, 70.09, 73.15],
+                [72.11, 73.26, 76.21, 79.18, 82.14, 85.11],
+                [86.36, 87.48, 90.34, 93.22, 96.10, 98.98],
+                [102.65, 103.73, 106.50, 109.29, 112.08, 114.88],
+                [121.38, 122.42, 125.10, 127.80, 130.51, 133.21],
+                [142.98, 143.99, 146.58, 149.19, 151.81, 154.43],
+                [167.77, 168.74, 171.24, 173.77, 176.30, 178.83],
+                [196.48, 197.41, 199.83, 202.26, 204.70, 207.15],
+                [229.56, 230.45, 232.78, 235.13, 237.49, 239.84],
+            ]
+        )
+
+        self.interpolate_pressure = RegularGridInterpolator(
+            (grid_compressions, grid_temperatures),
+            grid_pressures,
+            bounds_error=False,
+            fill_value=None,
+            method="cubic",
+        )
+
+        def _pressure(volume, temperature, params):
+            compression = 1.0 - volume / params["V_0"]
+            return self.interpolate_pressure([compression, temperature])[0] * 1.0e9
+
+        Calibrant.__init__(self, _pressure, "pressure", {"V_0": 10.207e-06})

burnman/calibrants/__init__.py

@@ -20,6 +20,7 @@
   - :mod:`~burnman.calibrants.Dubrovinsky_1998`
   - :mod:`~burnman.calibrants.Fei_2007`
   - :mod:`~burnman.calibrants.Fei_2016`
+  - :mod:`~burnman.calibrants.Holmes_1989`
   - :mod:`~burnman.calibrants.Huang_2016`
   - :mod:`~burnman.calibrants.LeGodec_2000`
   - :mod:`~burnman.calibrants.Litasov_2013`
@@ -35,6 +36,7 @@
   - :mod:`~burnman.calibrants.Speziale_2001`
   - :mod:`~burnman.calibrants.Tange_2009`
   - :mod:`~burnman.calibrants.Tateno_2019`
+  - :mod:`~burnman.calibrants.Tsuchiya_2003`
   - :mod:`~burnman.calibrants.Walker_2002`
   - :mod:`~burnman.calibrants.Zeng_2010`
   - :mod:`~burnman.calibrants.Zha_2004`
@@ -56,6 +58,7 @@
 from . import Dubrovinsky_1998
 from . import Fei_2007
 from . import Fei_2016
+from . import Holmes_1989
 from . import Huang_2016
 from . import LeGodec_2000
 from . import Litasov_2013
@@ -71,6 +74,7 @@
 from . import Speziale_2001
 from . import Tange_2009
 from . import Tateno_2019
+from . import Tsuchiya_2003
 from . import Walker_2002
 from . import Zeng_2010
 from . import Zha_2004

burnman/classes/anisotropicsolution.py

@@ -599,7 +599,7 @@ def _d2Fdqdz(self):
     @material_property
     def _d2Fdqdq_fixed_epsT_pinv(self):
         """
-        The second derivative of the Helmholtz energy
+        The inverse of the second derivative of the Helmholtz energy
         with respect to the structure parameters at constant
         strain and temperature. Often referred to as the
         susceptibility matrix.

burnman/utils/misc.py

@@ -70,24 +70,26 @@ def flatten(arr):
     )
 
 
-def pretty_print_values(popt, pcov, params):
+def pretty_print_values(popt, pcov, params, extra_decimal_places=0):
     """
     Takes a numpy array of parameters, the corresponding covariance matrix
     and a set of parameter names and prints the parameters and
     principal 1-s.d.uncertainties (np.sqrt(pcov[i][i]))
     in a nice text based format.
     """
     for i, p in enumerate(params):
-        p_rnd = round_to_n(popt[i], np.sqrt(pcov[i][i]), 1)
-        c_rnd = round_to_n(np.sqrt(pcov[i][i]), np.sqrt(pcov[i][i]), 1)
+        p_rnd = round_to_n(popt[i], np.sqrt(pcov[i][i]), 1 + extra_decimal_places)
+        c_rnd = round_to_n(
+            np.sqrt(pcov[i][i]), np.sqrt(pcov[i][i]), 1 + extra_decimal_places
+        )
 
         if p_rnd != 0.0:
             p_expnt = np.floor(np.log10(np.abs(p_rnd)))
         else:
             p_expnt = 0.0
 
         scale = np.power(10.0, p_expnt)
-        nd = p_expnt - np.floor(np.log10(np.abs(c_rnd)))
+        nd = p_expnt - np.floor(np.log10(np.abs(c_rnd))) + extra_decimal_places
         print(
             "{0:s}: ({1:{4}{5}f} +/- {2:{4}{5}f}) x {3:.0e}".format(
                 p, p_rnd / scale, c_rnd / scale, scale, 0, (nd) / 10.0

misc/benchmarks/calibrant_benchmarks.py

@@ -0,0 +1,67 @@
+from __future__ import absolute_import
+
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.image as mpimg
+from burnman import calibrants
+from burnman.tools.unitcell import molar_volume_from_unit_cell_volume
+
+
+def make_VPT_figure(calibrant, temperatures, figure, figure_extent, plot_extent):
+    name = str(calibrant).split(" ")[0][1:]
+    print(f"Checking {name}...")
+
+    fig = plt.figure(figsize=(6, 4))
+    fig.suptitle(f"{name}")
+
+    ax = [fig.add_subplot(1, 1, 1)]
+
+    fig1 = mpimg.imread(figure)
+    ax[0].imshow(fig1, extent=figure_extent, aspect="auto")
+
+    pressures = np.linspace(plot_extent[0] * 1.0e9, plot_extent[1] * 1.0e9, 101)
+    volumes = np.empty_like(pressures)
+    for T in temperatures:
+        for i, P in enumerate(pressures):
+            volumes[i] = calibrant.volume(P, T) / molar_volume_from_unit_cell_volume(
+                1.0, calibrant.params["Z"]
+            )
+
+        plt.plot(pressures / 1.0e9, volumes)
+
+    ax[0].set_xlim(plot_extent[0], plot_extent[1])
+    ax[0].set_ylim(plot_extent[2], plot_extent[3])
+    plt.show()
+
+
+def check_figures():
+    make_VPT_figure(
+        calibrants.Fei_2007.Au(),
+        [300.0, 1473.0, 2173.0],
+        "figures/Fei_2007_Au.png",
+        [0, 139.5, 50, 68],
+        [0, 139.5, 50, 68],
+    )
+
+    make_VPT_figure(
+        calibrants.Fei_2007.Pt(),
+        [300.0, 1473.0, 1873.0],
+        "figures/Fei_2007_Pt.png",
+        [-2, 125, 46.98, 61.02],
+        [0, 125, 47, 61],
+    )
+
+    V_0 = calibrants.Holmes_1989.Pt().params[
+        "V_0"
+    ] / molar_volume_from_unit_cell_volume(1.0, 4)
+    make_VPT_figure(
+        calibrants.Holmes_1989.Pt(),
+        [300.0],
+        "figures/Holmes_1989_Pt.png",
+        [0, 600, 0.6 * V_0, 1.1 * V_0],
+        [0, 600, 0.6 * V_0, 1.1 * V_0],
+    )
+
+
+if __name__ == "__main__":
+    check_figures()

misc/ref/calibrant_benchmarks.py.out

@@ -0,0 +1,3 @@
+Checking burnman.calibrants.Fei_2007.Au...
+Checking burnman.calibrants.Fei_2007.Pt...
+Checking burnman.calibrants.Holmes_1989.Pt...