diff --git a/burnman/data/input_raw_endmember_datasets/HeFESTo_2024_to_burnman.py b/burnman/data/input_raw_endmember_datasets/HeFESTo_2024_to_burnman.py index e38cb433..6be7776e 100644 --- a/burnman/data/input_raw_endmember_datasets/HeFESTo_2024_to_burnman.py +++ b/burnman/data/input_raw_endmember_datasets/HeFESTo_2024_to_burnman.py @@ -74,7 +74,13 @@ def rfloat(x, m=1.0): "lppv": "[Fefls][Fefls]", } -magnetic_structural_parameters = {"fea": 0.4, "mag": 0.4, "wu": 0.6} +magnetic_structural_parameters = { + "fea": {"p": 0.4}, + "mag": {"p": 0.4}, + "smag": {"p": 0.4}, + "hmag": {"p": 0.4}, + "wu": {"p": 0.6}, +} hefesto_path = "HeFESTo_parameters_010123" mbrdir = pathlib.Path(hefesto_path) @@ -358,7 +364,7 @@ def rfloat(x, m=1.0): # but with the zero Gibbs energy and entropy # assigned to the fully ordered form, # not the disordered form,. - p = magnetic_structural_parameters[name] + p = magnetic_structural_parameters[name]["p"] Scrit = float(idict["S_crit"]) magnetic_moment = np.exp(Scrit / gas_constant) - 1.0 A = (518.0 / 1125.0) + (11692.0 / 15975.0) * ((1.0 / p) - 1.0) @@ -390,7 +396,7 @@ def rfloat(x, m=1.0): { "Tc_0": float(idict["T_crit"]), "S_D": float(idict["S_crit"]), - "V_D": float(idict["V_crit"]), + "V_D": float(idict["V_crit"]) / 1.0e6, }, ] ] @@ -471,6 +477,33 @@ def rfloat(x, m=1.0): " RelaxedSolution.__init__(self, solution, vrel, vunrel)\n" ) +bdg_relaxed_class = ( + "class bridgmanite_relaxed(RelaxedSolution):\n" + " def __init__(self):\n" + ' """RelaxedSolution model for bridgmanite (pv).\n' + " Only the spin transition is relaxed.\n" + " Endmembers (and site species distributions) are given in the order:\n" + " - mgpv ([Mg][Si])\n" + " - fepv ([Fe][Si])\n" + " - alpv ([Al][Al])\n" + " - hepv + hlpv ([Fef][Fef,Fels])\n" + " - fapv ([Fef][Al])\n" + " - crpv ([Cr][Cr])\n" + " The equilibrium spin state is calculated automatically.\n" + ' """\n' + " solution = bridgmanite()\n" + " vrel = [[0.0, 0., 0., 0.5, -0.5, 0.0, 0.0]]\n" + " vunrel = [\n" + " [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],\n" + " [0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0],\n" + " [0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0],\n" + " [0.0, 0.0, 0.0, 0.5, 0.5, 0.0, 0.0],\n" + " [0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0],\n" + " [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0],\n" + " ]\n" + " RelaxedSolution.__init__(self, solution, vrel, vunrel)\n" +) + print('"""\n' "ENDMEMBERS\n" '"""\n') for key, mbr_prm in sorted(mbr_params.items()): if key == "st": @@ -533,6 +566,8 @@ def rfloat(x, m=1.0): if key == "mw": print(fper_relaxed_class) + if key == "pv": + print(bdg_relaxed_class) print('\n"""\n' "ENDMEMBER ALIASES\n" '"""\n') diff --git a/burnman/minerals/SLB_2024.py b/burnman/minerals/SLB_2024.py index 5bc8b6c5..acaa408a 100644 --- a/burnman/minerals/SLB_2024.py +++ b/burnman/minerals/SLB_2024.py @@ -1147,7 +1147,18 @@ def __init__(self): } self.property_modifiers = [ - ["landau_slb_2022", {"Tc_0": 845.5, "S_D": 43.1758, "V_D": 0.0}] + [ + "magnetic_chs", + { + "structural_parameter": 0.4, + "curie_temperature": [845.5, 0.0], + "magnetic_moment": [178.9816942215617, 0.0], + }, + ], + [ + "linear", + {"delta_E": 33048.07971316818, "delta_S": 43.1758, "delta_V": 0.0}, + ], ] Mineral.__init__(self) @@ -1634,7 +1645,10 @@ def __init__(self): } self.property_modifiers = [ - ["landau_slb_2022", {"Tc_0": 467.0, "S_D": 10.0, "V_D": 0.8}] + [ + "landau_slb_2022", + {"Tc_0": 467.0, "S_D": 10.0, "V_D": 8.000000000000001e-07}, + ] ] Mineral.__init__(self) @@ -1810,7 +1824,10 @@ def __init__(self): } self.property_modifiers = [ - ["landau_slb_2022", {"Tc_0": 847.0, "S_D": 5.76, "V_D": 1.35936}] + [ + "landau_slb_2022", + {"Tc_0": 847.0, "S_D": 5.76, "V_D": 1.3593599999999998e-06}, + ] ] Mineral.__init__(self) @@ -1838,7 +1855,18 @@ def __init__(self): } self.property_modifiers = [ - ["landau_slb_2022", {"Tc_0": 845.5, "S_D": 43.1758, "V_D": 0.0}] + [ + "magnetic_chs", + { + "structural_parameter": 0.4, + "curie_temperature": [845.5, 0.0], + "magnetic_moment": [178.9816942215617, 0.0], + }, + ], + [ + "linear", + {"delta_E": 33048.07971316818, "delta_S": 43.1758, "delta_V": 0.0}, + ], ] Mineral.__init__(self) @@ -2372,6 +2400,32 @@ def __init__(self, molar_fractions=None): Solution.__init__(self, molar_fractions=molar_fractions) +class bridgmanite_relaxed(RelaxedSolution): + def __init__(self): + """RelaxedSolution model for bridgmanite (pv). + Only the spin transition is relaxed. + Endmembers (and site species distributions) are given in the order: + - mgpv ([Mg][Si]) + - fepv ([Fe][Si]) + - alpv ([Al][Al]) + - hepv + hlpv ([Fef][Fef,Fels]) + - fapv ([Fef][Al]) + - crpv ([Cr][Cr]) + The equilibrium spin state is calculated automatically. + """ + solution = bridgmanite() + vrel = [[0.0, 0.0, 0.0, 0.5, -0.5, 0.0, 0.0]] + vunrel = [ + [1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.5, 0.5, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0], + ] + RelaxedSolution.__init__(self, solution, vrel, vunrel) + + class ringwoodite(Solution): def __init__(self, molar_fractions=None): """SymmetricRegularSolution model for ringwoodite (ri). diff --git a/misc/benchmarks/figures/SLB_2024_Fe2O3_V.png b/misc/benchmarks/figures/SLB_2024_Fe2O3_V.png new file mode 100644 index 00000000..69a56c11 Binary files /dev/null and b/misc/benchmarks/figures/SLB_2024_Fe2O3_V.png differ diff --git a/misc/benchmarks/figures/SLB_2024_Fe_O_fO2.png b/misc/benchmarks/figures/SLB_2024_Fe_O_fO2.png new file mode 100644 index 00000000..e6261c15 Binary files /dev/null and b/misc/benchmarks/figures/SLB_2024_Fe_O_fO2.png differ diff --git a/misc/benchmarks/figures/SLB_2024_Fe_O_fO2_T.png b/misc/benchmarks/figures/SLB_2024_Fe_O_fO2_T.png new file mode 100644 index 00000000..8b0716d1 Binary files /dev/null and b/misc/benchmarks/figures/SLB_2024_Fe_O_fO2_T.png differ diff --git a/misc/benchmarks/slb_2024_benchmarks.py b/misc/benchmarks/slb_2024_benchmarks.py index 2e96c9a1..9480bbb2 100644 --- a/misc/benchmarks/slb_2024_benchmarks.py +++ b/misc/benchmarks/slb_2024_benchmarks.py @@ -1,33 +1,24 @@ from __future__ import absolute_import from burnman import Composite, equilibrate +from burnman.constants import gas_constant +from burnman.tools.polytope import simplify_composite_with_composition from burnman.tools.eos import check_eos_consistency from burnman.tools.unitcell import molar_volume_from_unit_cell_volume from burnman.minerals import SLB_2024 +from burnman.minerals.HP_2011_fluids import O2 import numpy as np import matplotlib.pyplot as plt import matplotlib.image as mpimg -def check_stishovite(): - stv = SLB_2024.stishovite() - - fig = plt.figure(figsize=(10, 5)) - ax = [fig.add_subplot(1, 1, i) for i in range(1, 2)] - - fig1 = mpimg.imread("figures/SLB_2024_stv_VS.png") - ax[0].imshow(fig1, extent=[0.0, 100.0, 0.0, 9.0], aspect="auto") - - for T in [300, 1000, 2000]: - pressures = np.linspace(1.0e5, 100.0e9, 101) - temperatures = pressures * 0.0 + T - Vs = stv.evaluate(["shear_wave_velocity"], pressures, temperatures)[0] - ax[0].plot(pressures / 1.0e9, Vs / 1.0e3, linestyle=":", linewidth=3.0) - plt.show() +def check_bcc_iron_consistency(): + fe_bcc = SLB_2024.alpha_bcc_iron() + check_eos_consistency(fe_bcc, verbose=True) -def check_fper_relaxed(): +def check_fig_1_fper_relaxed(): fper = SLB_2024.ferropericlase_relaxed() c = molar_volume_from_unit_cell_volume(1.0, SLB_2024.periclase().params["Z"]) @@ -35,6 +26,7 @@ def check_fper_relaxed(): temperatures = pressures * 0.0 + 300.0 fig = plt.figure(figsize=(10, 5)) + fig.suptitle("Figure 1 (Stixrude and Lithgow-Bertelloni, 2024)") ax = [fig.add_subplot(1, 2, i) for i in range(1, 3)] fig1 = mpimg.imread("figures/SLB_2024_fper_V.png") @@ -54,12 +46,36 @@ def check_fper_relaxed(): plt.show() -def check_bcc_iron_consistency(): - fe_bcc = SLB_2024.alpha_bcc_iron() - check_eos_consistency(fe_bcc, verbose=True) +def check_fig_3_fcc_ferric_fper(): + fcc = SLB_2024.gamma_fcc_iron() + fper = SLB_2024.ferropericlase() + + x_MgOs = np.linspace(0.001, 0.999, 101) + x_Fe3oversumFes = np.empty_like(x_MgOs) + + fig = plt.figure(figsize=(10, 5)) + fig.suptitle("Figure 3 (Stixrude and Lithgow-Bertelloni, 2024)") + ax = [fig.add_subplot(1, 1, i) for i in range(1, 2)] + + for equality_constraints in [ + [["P", 1.0e5], ["T", 1473.0]], + [["P", 21.0e9], ["T", 1673.0]], + ]: + for i, x_MgO in enumerate(x_MgOs): + composition = {"Mg": x_MgO, "Fe": 1.0 - x_MgO, "O": 1.0} + fper.set_composition([x_MgO, 1.0 - x_MgO, 0.0, 0.0, 0.0]) + assemblage = Composite([fper, fcc], [1.0, 0.0]) + assemblage = simplify_composite_with_composition(assemblage, composition) + equilibrate(composition, assemblage, equality_constraints) + + f = assemblage.phases[0].formula + x_Fe3oversumFes[i] = 2.0 * ((f["O"] - f["Mg"]) / f["Fe"] - 1.0) + + ax[0].plot(x_MgOs, x_Fe3oversumFes, linestyle=":", linewidth=3.0) + plt.show() -def check_iron_Cp_V(): +def check_fig_6a_iron_Cp_V(): fe_bcc = SLB_2024.alpha_bcc_iron() fe_fcc = SLB_2024.gamma_fcc_iron() @@ -67,6 +83,7 @@ def check_iron_Cp_V(): pressures = temperatures * 0.0 + 1.0e5 fig = plt.figure(figsize=(10, 5)) + fig.suptitle("Figure 6a (Stixrude and Lithgow-Bertelloni, 2024)") ax = [fig.add_subplot(1, 2, i) for i in range(1, 3)] fig1 = mpimg.imread("figures/SLB_2024_Fe_Cp_V.png") @@ -83,12 +100,13 @@ def check_iron_Cp_V(): plt.show() -def check_iron_phase_diagram(): +def check_fig_6c_iron_phase_diagram(): fe_bcc = SLB_2024.alpha_bcc_iron() fe_fcc = SLB_2024.gamma_fcc_iron() fe_hcp = SLB_2024.epsilon_hcp_iron() fig = plt.figure(figsize=(10, 5)) + fig.suptitle("Figure 6c (Stixrude and Lithgow-Bertelloni, 2024)") ax = [fig.add_subplot(1, 1, i) for i in range(1, 2)] fig1 = mpimg.imread("figures/SLB_2024_Fe_phase_diagram.png") ax[0].imshow(fig1, extent=[0.0, 250.0, 300.0, 5800.0], aspect="auto") @@ -96,7 +114,7 @@ def check_iron_phase_diagram(): # Calculate the invariant point irons = Composite([fe_bcc, fe_fcc, fe_hcp], [1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0]) irons.set_state(5.0e9, 600.0) - sol = equilibrate( + equilibrate( {"Fe": 1.0}, irons, equality_constraints=[ @@ -136,9 +154,159 @@ def check_iron_phase_diagram(): plt.show() +def check_fig_7_fO2(): + hem = SLB_2024.hematite() + mag = SLB_2024.smag() + hmag = SLB_2024.high_pressure_magnetit() + wu = SLB_2024.ferropericlase() + wu.set_composition([0.0, 0.97, 0.005, 0.0, 0.025]) + wu = simplify_composite_with_composition( + Composite([wu]), {"Fe": 1, "O": 1.01} + ).phases[0] + + hepv = SLB_2024.hepv() + hppv = SLB_2024.hppv() + Fe_bcc = SLB_2024.alpha_bcc_iron() + Fe_fcc = SLB_2024.gamma_fcc_iron() + Fe_hcp = SLB_2024.epsilon_hcp_iron() + q = SLB_2024.quartz() + fa = SLB_2024.fayalite() + O2_gas = O2() + + O2_gas.set_state(1.0e5, 298.15) + f = O2_gas.S * 298.15 + + fig = plt.figure(figsize=(10, 5)) + fig.suptitle("Figure 7 (Stixrude and Lithgow-Bertelloni, 2024)") + ax = [fig.add_subplot(1, 1, i) for i in range(1, 2)] + + fig1 = mpimg.imread("figures/SLB_2024_Fe_O_fO2_T.png") + ax[0].imshow(fig1, extent=[700.0, 1500.0, -32.0, 0.0], aspect="auto") + + print("fO2 values at 1500 K:") + for phases in [[fa, Fe_fcc, q], [fa, Fe_bcc, q], [fa, mag, q], [mag, hem, hem]]: + assemblage = Composite(phases, [0.2, 0.3, 0.5]) + temperatures = np.linspace(700.0, 1500.0) + logfO2 = np.empty_like(temperatures) + + for i, T in enumerate(temperatures): + assemblage.set_state(1.0e5, T) + mu_O2 = assemblage.chemical_potential([{"O": 2.0}])[0] + O2_gas.set_state(1.0e5, T) + O2_gibbs = O2_gas.gibbs + f + logfO2[i] = (mu_O2 - O2_gibbs) / (gas_constant * T) / np.log(10) + print([ph.name for ph in phases]) + print(f"log10fO2 = {logfO2[-1]:.3f}") + ax[0].plot(temperatures, logfO2, linestyle=":", linewidth=3.0) + + ax[0].set_xlim(700.0, 1500.0) + ax[0].set_ylim(-32, 0) + plt.show() + + P_bounds = [[1.0e5, 100.0e9] for i in range(9)] + + a1 = Composite([mag, hmag], [0.5, 0.5]) + a2 = Composite([hem, hepv], [0.5, 0.5]) + a3 = Composite([hepv, hppv], [0.5, 0.5]) + a4 = Composite([hmag, hppv, wu], [0.3, 0.2, 0.5]) + a5 = Composite([Fe_fcc, Fe_hcp], [0.5, 0.5]) + + for i, a in enumerate([a1, a2, a3, a4, "null", a1, "null", a5]): + if a != "null": + equality_constraints = [ + ["phase_fraction", [a.phases[0], 0.0]], + ["T", 1500.0], + ] + equilibrate(a.phases[0].formula, a, equality_constraints) + P_bounds[i][1] = a.pressure + P_bounds[i + 1][0] = a.pressure + + P_bounds[6][1] = P_bounds[3][1] + + a1 = Composite([hem, mag], [0.5, 0.5]) + a2 = Composite([hem, hmag], [0.5, 0.5]) + a3 = Composite([hepv, hmag], [0.5, 0.5]) + a4 = Composite([hppv, hmag], [0.5, 0.5]) + a5 = Composite([hppv, wu], [0.5, 0.5]) + a6 = Composite([mag, wu], [0.5, 0.5]) + a7 = Composite([hmag, wu], [0.5, 0.5]) + a8 = Composite([Fe_fcc, wu], [0.5, 0.5]) + a9 = Composite([Fe_hcp, wu], [0.5, 0.5]) + + fig = plt.figure(figsize=(10, 5)) + ax = [fig.add_subplot(1, 1, i) for i in range(1, 2)] + + fig1 = mpimg.imread("figures/SLB_2024_Fe_O_fO2.png") + ax[0].imshow(fig1, extent=[0, 100.0, -14.0, 22.0], aspect="auto") + + for i, assemblage in enumerate([a1, a2, a3, a4, a5, a6, a7, a8, a9]): + pressures = np.linspace(P_bounds[i][0], P_bounds[i][1], 11) + logfO2 = np.empty_like(pressures) + T = 1500.0 + for i, P in enumerate(pressures): + assemblage.set_state(P, T) + sol = equilibrate(assemblage.formula, assemblage, [["P", P], ["T", T]])[0] + if sol.success: + mu_O2 = assemblage.chemical_potential([{"O": 2.0}])[0] + O2_gas.set_state(1.0e5, T) + O2_gibbs = O2_gas.gibbs + f + logfO2[i] = (mu_O2 - O2_gibbs) / (gas_constant * T) / np.log(10) + ax[0].plot(pressures / 1.0e9, logfO2, linestyle=":", linewidth=3.0) + + ax[0].set_ylim(-14, 22) + plt.show() + + +def check_fig_a1_fe2o3(): + hem = SLB_2024.hematite() + hppv = SLB_2024.hppv() + bdg = SLB_2024.bridgmanite_relaxed() + bdg.set_composition([0.0, 0.0, 0.0, 1.0, 0.0, 0.0]) + + fig = plt.figure(figsize=(10, 5)) + fig.suptitle("Figure A1 (Stixrude and Lithgow-Bertelloni, 2024)") + ax = [fig.add_subplot(1, 1, i) for i in range(1, 2)] + + fig1 = mpimg.imread("figures/SLB_2024_Fe2O3_V.png") + ax[0].imshow(fig1, extent=[0.0, 120.0, 19.0, 31.0], aspect="auto") + + pressures = np.linspace(1.0e5, 120.0e9, 101) + temperatures = 300.0 + 0.0 * pressures + for phase in [hem, hppv, bdg]: + volumes = phase.evaluate(["V"], pressures, temperatures)[0] + ax[0].plot(pressures / 1.0e9, volumes * 1.0e6, linestyle=":", linewidth=3.0) + + bdg.set_composition([0.5, 0.0, 0.0, 0.5, 0.0, 0.0]) + volumes = bdg.evaluate(["V"], pressures, temperatures)[0] + ax[0].plot(pressures / 1.0e9, volumes * 1.0e6, linestyle=":", linewidth=3.0) + + plt.show() + + +def check_fig_a4_stishovite(): + stv = SLB_2024.stishovite() + + fig = plt.figure(figsize=(10, 5)) + fig.suptitle("Figure A4 (Stixrude and Lithgow-Bertelloni, 2024)") + ax = [fig.add_subplot(1, 1, i) for i in range(1, 2)] + + fig1 = mpimg.imread("figures/SLB_2024_stv_VS.png") + ax[0].imshow(fig1, extent=[0.0, 100.0, 0.0, 9.0], aspect="auto") + + for T in [300, 1000, 2000]: + pressures = np.linspace(1.0e5, 100.0e9, 101) + temperatures = pressures * 0.0 + T + Vs = stv.evaluate(["shear_wave_velocity"], pressures, temperatures)[0] + ax[0].plot(pressures / 1.0e9, Vs / 1.0e3, linestyle=":", linewidth=3.0) + plt.show() + + if __name__ == "__main__": - check_fper_relaxed() check_bcc_iron_consistency() - check_iron_phase_diagram() - check_iron_Cp_V() - check_stishovite() + check_fig_1_fper_relaxed() + check_fig_3_fcc_ferric_fper() + check_fig_6a_iron_Cp_V() + check_fig_6c_iron_phase_diagram() + check_fig_7_fO2() + check_fig_a1_fe2o3() + check_fig_a4_stishovite() diff --git a/misc/ref/slb_2024_benchmarks.py.out b/misc/ref/slb_2024_benchmarks.py.out index 0a1bf2f6..1b400a15 100644 --- a/misc/ref/slb_2024_benchmarks.py.out +++ b/misc/ref/slb_2024_benchmarks.py.out @@ -16,3 +16,12 @@ Vp = np.sqrt((K_S + 4G/3)/rho) : True Vs = np.sqrt(G_S/rho) : True All EoS consistency constraints satisfied for 'burnman.minerals.SLB_2024.fea' The bcc-fcc-hcp invariant is at 11.250 GPa, 827.042 K +fO2 values at 1500 K: +['Fayalite', 'gamma_(fcc)_Iron', 'Quartz'] +log10fO2 = -12.063 +['Fayalite', 'alpha_(bcc)_Iron', 'Quartz'] +log10fO2 = -12.063 +['Fayalite', 'Magnetite', 'Quartz'] +log10fO2 = -8.233 +['Magnetite', 'Hematite', 'Hematite'] +log10fO2 = -2.635