diff --git a/burnman/calibrants/Decker_1971.py b/burnman/calibrants/Decker_1971.py index 4eb9bc03..60a73018 100644 --- a/burnman/calibrants/Decker_1971.py +++ b/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, diff --git a/burnman/calibrants/Fei_2007.py b/burnman/calibrants/Fei_2007.py index ac8e529d..56e77bbb 100644 --- a/burnman/calibrants/Fei_2007.py +++ b/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) diff --git a/burnman/calibrants/Holmes_1989.py b/burnman/calibrants/Holmes_1989.py new file mode 100644 index 00000000..b0efb85f --- /dev/null +++ b/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) diff --git a/burnman/calibrants/Matsui_2012.py b/burnman/calibrants/Matsui_2012.py index c4437183..1b712ac5 100644 --- a/burnman/calibrants/Matsui_2012.py +++ b/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"] diff --git a/burnman/calibrants/Tsuchiya_2003.py b/burnman/calibrants/Tsuchiya_2003.py new file mode 100644 index 00000000..23d2ad18 --- /dev/null +++ b/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}) diff --git a/burnman/calibrants/__init__.py b/burnman/calibrants/__init__.py index 65104b82..3ceca19f 100644 --- a/burnman/calibrants/__init__.py +++ b/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 diff --git a/burnman/classes/anisotropicsolution.py b/burnman/classes/anisotropicsolution.py index 28681cfc..b6fc15f0 100644 --- a/burnman/classes/anisotropicsolution.py +++ b/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. diff --git a/burnman/utils/misc.py b/burnman/utils/misc.py index cfb36d6b..6b7a7baa 100644 --- a/burnman/utils/misc.py +++ b/burnman/utils/misc.py @@ -70,7 +70,7 @@ 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 @@ -78,8 +78,10 @@ def pretty_print_values(popt, pcov, params): 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))) @@ -87,7 +89,7 @@ def pretty_print_values(popt, pcov, params): 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 diff --git a/misc/benchmarks/calibrant_benchmarks.py b/misc/benchmarks/calibrant_benchmarks.py new file mode 100644 index 00000000..7ecb40d6 --- /dev/null +++ b/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() diff --git a/misc/benchmarks/figures/Fei_2007_Au.png b/misc/benchmarks/figures/Fei_2007_Au.png new file mode 100644 index 00000000..5bad1efe Binary files /dev/null and b/misc/benchmarks/figures/Fei_2007_Au.png differ diff --git a/misc/benchmarks/figures/Fei_2007_Pt.png b/misc/benchmarks/figures/Fei_2007_Pt.png new file mode 100644 index 00000000..c61e3a2c Binary files /dev/null and b/misc/benchmarks/figures/Fei_2007_Pt.png differ diff --git a/misc/benchmarks/figures/Holmes_1989_Pt.png b/misc/benchmarks/figures/Holmes_1989_Pt.png new file mode 100644 index 00000000..59150805 Binary files /dev/null and b/misc/benchmarks/figures/Holmes_1989_Pt.png differ diff --git a/misc/ref/calibrant_benchmarks.py.out b/misc/ref/calibrant_benchmarks.py.out new file mode 100644 index 00000000..7cffaede --- /dev/null +++ b/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... diff --git a/misc/ref/example_calibrants.py.out b/misc/ref/example_calibrants.py.out index 06eebe13..d0f96700 100644 --- a/misc/ref/example_calibrants.py.out +++ b/misc/ref/example_calibrants.py.out @@ -3,15 +3,15 @@ Volume: 19.0287 +/- 0.1903 m^3/mol Temperature: 1073.15 +/- 10.0 K Pressure estimated using the Decker equation of state: -21.6417 GPa +21.7355 GPa Pressure with propagated uncertainty estimated using the Decker equation of state: -21.6417 +/- 1.0010 GPa +21.7355 +/- 1.0254 GPa PVT correlation matrix: -[[ 1. -0.99955588 0.02979992] - [-0.99955588 1. 0. ] - [ 0.02979992 0. 1. ]] +[[ 1. -0.99957671 0.02909288] + [-0.99957671 1. 0. ] + [ 0.02909288 0. 1. ]] The uncertainties and correlation matrix only take into account the uncertainties in the measured volume and temperature, not the uncertainties in the calibrated equation of state. @@ -19,8 +19,8 @@ not the uncertainties in the calibrated equation of state. Consistency checks: Volume: 19.0287 +/- 0.1903 m^3/mol Temperature: 1073.15 +/- 10.0000 K -V-T correlation: -0.000000 +V-T correlation: 0.000000 Result of conversion from Decker calibration to Decker calibration: -21.6417 +/- 1.0010 GPa +21.7355 +/- 1.0254 GPa (This should be the same as the pressures above)