diff --git a/burnman/classes/solutionmodel.py b/burnman/classes/solutionmodel.py index fe6ef84f8..d7736696e 100644 --- a/burnman/classes/solutionmodel.py +++ b/burnman/classes/solutionmodel.py @@ -108,7 +108,7 @@ def _non_ideal_interactions_subreg(p, n_endmembers, Wijk): return Wint -def logish(x, eps=1.0e-5): +def logish(x, eps=1.0e-7): """ 2nd order series expansion of log(x) about eps: log(eps) - sum_k=1^infty (f_eps)^k / k @@ -434,7 +434,13 @@ def __init__(self, endmembers): def _calculate_endmember_configurational_entropies(self): S_conf = -( constants.gas_constant - * (self.endmember_noccupancies * logish(self.endmember_occupancies)).sum(-1) + * ( + self.endmember_noccupancies + * ( + logish(self.endmember_noccupancies) + - logish(self.site_multiplicities) + ) + ).sum(-1) ) self.endmember_configurational_entropies = S_conf 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 new file mode 100644 index 000000000..bec3e914a --- /dev/null +++ b/burnman/data/input_raw_endmember_datasets/HeFESTo_2024_to_burnman.py @@ -0,0 +1,607 @@ +# 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. + + +# This is a script that converts a version of the +# HeFESTo data format into the standard burnman +# format (printed to stdout) +# This script is used to generate the 2024 dataset, +# taking parameters from https://github.com/stixrude/HeFESTo_parameters_010123 +# via git clone https://github.com/stixrude/HeFESTo_parameters_010123.git +# The SLB_2024 BurnMan file is then created via the commands +# python HeFESTo_2024_to_burnman.py > ../../minerals/SLB_2024.py; black ../../minerals/SLB_2024.py +# and from . import SLB_2024 added to the minerals __init__.py + + +import os +import pathlib +import pandas as pd +import numpy as np +from burnman.utils.chemistry import dictionarize_formula +from burnman.constants import gas_constant +from itertools import groupby +import re + + +def all_equal(iterable): + g = groupby(iterable) + return next(g, True) and not next(g, False) + + +def rfloat(x, m=1.0): + return round(float(x) * m, ndigits=10) + + +# Hard-coded into HeFESTo +""" +Ferropericlase: +Crystallography of wustites shows that the locations of the tetrahedral site +and the octahedral vacancies are highly correlated, with the tetrahedral +site being surrounded by a tetrahedral array of octahedral vacancies +(Gavarri et al. 1979 ).Because of this highly structured arrangement, +we ignore the contribution to the ideal entropy of mixing from the +tetrahedral site and the vacant octahedral sites. + +...the formula of mag : [Fe]Fe2O4, with the square brackets indicating +that the first site does not participate in mixing. Thus, the ideal entropy +of mixing of an equimolar solution of wu and mag is 2Rln2 per 4 oxygens, +where R is the gas constant. + +Spinel: +Pure magnetite is known to have an inverse arrangement of cations at +ambient conditions, rather than normal, as we have assumed. + +Bridgmanite: +hepv and hlpv where only the cation on the +smaller B site undergoes the high spin to low spin transition... +""" + +hardcoded_site_occupancies = { + "wu": "[Fe]2[Fe]2", + "wuls": "[Fels]2[Fels]2", + "mag": "[Fef]0[Fef]2", # first site not counted + "smag": "[Fe][Fef][Fef]", + "hepv": "[Fef][Fef]", + "hlpv": "[Fef][Fefls]", + "fapv": "[Fef][Al]", + "hppv": "[Fef][Fef]", + "hmag": "[Fe][Fef][Fef]", + "acm": "[Na][Fef][Si]2", + "andr": "[Ca]3[Fef][Fef]", + "hem": "[Fef][Fef]", + "lppv": "[Fefls][Fefls]", +} + +magnetic_structural_parameters = { + "fea": {"p": 0.4}, + "mag": {"p": 0.4}, + "smag": {"p": 0.4}, + "hmag": {"p": 0.4}, +} + +hefesto_path = "HeFESTo_parameters_010123" +mbrdir = pathlib.Path(hefesto_path) +phasedir = pathlib.Path(f"{hefesto_path}/phase") + +param_dict = {} +with open("hefesto_parameter_names.txt") as file: + lines = [line.rstrip() for line in file] + for line in lines: + param_dict[" ".join(line.split()[1:])] = line.split()[0] + +ignored_mbrs = [] + +solution_aliases = { + "c2c": "c2c_pyroxene", + "cf": "calcium_ferrite_structured_phase", + "cpx": "clinopyroxene", + "gt": "garnet", + "il": "ilmenite", + "mw": "ferropericlase", + "nal": "new_aluminous_phase", + "ol": "olivine", + "opx": "orthopyroxene", + "plg": "plagioclase", + "ppv": "post_perovskite", + "pv": "bridgmanite", + "ri": "ringwoodite", + "sp": "mg_fe_aluminous_spinel", + "wa": "wadsleyite", +} + +# solutions +sol_params = {} +for f in phasedir.iterdir(): + # Do not process liquid!! + if os.stat(f).st_size > 0 and str(f).split("/")[-1] != "liq": + name = str(f).split("/")[-1] + d = pd.read_csv(f, sep="\\s+") + + endmembers = list(d.columns) + + n_mbrs = len(endmembers) + + data = d.values.tolist() + + n_data = len(data) + v0 = n_data - n_mbrs + + # Preprocess interaction matrices + # to get rid of ignored endmembers + use_indices = np.array( + [i for i, mbr in enumerate(endmembers) if mbr not in ignored_mbrs] + ) + use_endmembers = np.array( + [mbr for mbr in endmembers if mbr not in ignored_mbrs] + ) + n_use_mbrs = len(use_indices) + + energy_matrix = np.array( + np.array(data)[np.ix_(use_indices, use_indices)], dtype=float + ) + volume_matrix = np.array( + np.array(data)[np.ix_(v0 + use_indices, use_indices)], dtype=float + ) + + if n_data == 2 * n_mbrs + 1: + energy_interactions = [] + energy_interactions_2 = [] + volume_interactions = [] + for i in range(n_use_mbrs - 1): + j = i + 1 + energy_interactions.append( + list(np.around(energy_matrix[i][j:] * 1000.0, decimals=6)) + ) + energy_interactions_2.append( + list(np.around(energy_matrix[i][:j] * 1000.0, decimals=6)) + ) + volume_interactions.append( + list(np.around(volume_matrix[i][j:] * 1.0e-6, decimals=20)) + ) + + max_abs_e = np.max( + np.abs([item for sublist in energy_interactions for item in sublist]) + ) + max_abs_v = np.max( + np.abs([item for sublist in volume_interactions for item in sublist]) + ) + max_abs_lower_triangular_e2 = np.max( + np.abs([item for sublist in energy_interactions_2 for item in sublist]) + ) + assert max_abs_lower_triangular_e2 < 1.0e-3 + + # print(endmembers, n_mbrs) + # print(energy_interactions) + if max_abs_v > 1.0e-20: + # print(volume_interactions) + pass + + alphas = [] + mbr_sites = [] + for mbr in use_endmembers: + with open(f"{hefesto_path}/{mbr}") as file: + lines = [line.rstrip() for line in file] + i_laar = [i for i, line in enumerate(lines) if "Laar" in line][0] + alphas.append(float(lines[i_laar].split()[0])) + + formula_raw = lines[0].split()[0] + + # Catch for MgTs mixing + if mbr == "mgts": + formula_raw = "Mg_1Al_1Si_2O_6" + + # Catch for plag mixing + if mbr == "ab": + formula_raw = "Na_1Al_2Si_2O_8" + + # Catch for mt mixing + if mbr == "mag": + formula_raw = "(Vac_1Fe_1)Fe_2O_4" + + starts = [i for i, char in enumerate(formula_raw) if char.isupper()] + mixstarts = [i for i, char in enumerate(formula_raw) if char == "("] + mixends = [i for i, char in enumerate(formula_raw) if char == ")"] + + for i in range(len(mixstarts)): + ms = mixstarts[i] + me = mixends[i] + starts = [i for i in starts if i < ms or i > me] + + starts.extend(mixstarts) + starts.sort() + sites = [ + formula_raw[i:j] for i, j in zip(starts, starts[1:] + [None]) + ] + mbr_sites.append(sites) + + n_sites = [len(sites) for sites in mbr_sites] + + assert all_equal(n_sites) + + active_sites = [] + for i in range(n_sites[0]): + occupancies = [sites[i] for sites in mbr_sites] + if not all_equal(occupancies): + active_sites.append(i) + + site_formula_strings = [] + totals = [] + for mbr_site in mbr_sites: + totals.append([]) + site_formula_strings.append("") + for i in active_sites: + # Get number of atoms on sites + n_atoms = [int(s) for s in re.findall(r"\d+", mbr_site[i])] + + species_starts = [ + i for i, char in enumerate(mbr_site[i]) if char.isupper() + ] + species_ends = [ + i for i, char in enumerate(mbr_site[i]) if char == "_" + ] + total_atoms = np.sum(n_atoms) + totals[-1].append(total_atoms) + + site_string = "[" + if len(n_atoms) > 1: + species = "" + for j in range(len(n_atoms)): + species += mbr_site[i][species_starts[j] : species_ends[j]] + species += f"{n_atoms[j]}/{total_atoms}" + site_string += species + else: + site_string += mbr_site[i][species_starts[0] : species_ends[0]] + + site_string += "]" + if total_atoms != 1: + site_string += f"{total_atoms}" + + site_formula_strings[-1] += site_string + + for i in range(len(totals[0])): + t = [total[i] for total in totals] + assert all_equal(t) + + solution_is_symmetric = np.max(np.abs(np.array(alphas) - 1.0)) < 1.0e-5 + + sol_params[name] = { + "name": name, + "n_mbrs": n_use_mbrs, + "mbr_names": use_endmembers, + "mbr_site_formulae": site_formula_strings, + "alphas": alphas, + } + + if max_abs_e < 1.0: + sol_params[name]["solution_type"] = "IdealSolution" + elif solution_is_symmetric: + sol_params[name]["solution_type"] = "SymmetricRegularSolution" + else: + sol_params[name]["solution_type"] = "AsymmetricRegularSolution" + + if max_abs_e > 1.0: + sol_params[name]["energy_interaction"] = energy_interactions + if max_abs_v > 1.0e-20: + sol_params[name]["volume_interaction"] = volume_interactions + + else: + pass + +# endmembers +mbr_params = {} +mbr_mods = {} +mbr_names = {} +for f in mbrdir.iterdir(): + if os.stat(f).st_size > 0: + name = str(f).split("/")[-1] + + process = False + try: + with open(f) as file: + lines = [line.rstrip().split() for line in file] + if lines[-1][-1] == "prime": + process = True + if lines[0][-1] == "Liquid": + process = False + except IsADirectoryError: + pass + + if process: + form = lines[0][0].replace("_", "").replace("(", "").replace(")", "") + formula = dictionarize_formula(form) + full_name = "_".join(lines[0][1:]) + # formula = formula_to_string(formula) + # print(formula) + # print(name) + idict = {} + for line in lines[1:]: + key = param_dict[" ".join(line[1:])] + idict[key] = line[0] + + if "T_crit" not in idict: + idict["T_crit"] = "0.00000" + + assert idict["debye"] == "1.00000" + + if idict["zpp"] != "0.00000": + process = False + + if process: + mbr_params[name] = { + "name": full_name, + "formula": formula, + "equation_of_state": "slb3", + "F_0": rfloat(idict["F_0"], 1.0e3), + "V_0": rfloat(idict["V_0"], 1.0e-6), + "K_0": rfloat(idict["K_0"], 1.0e9), + "Kprime_0": rfloat(idict["Kprime_0"]), + "Debye_0": rfloat(idict["Theta_0"]), + "grueneisen_0": rfloat(idict["gamma_0"]), + "q_0": rfloat(idict["q_0"]), + "G_0": rfloat(idict["G_0"], 1.0e9), + "Gprime_0": rfloat(idict["Gprime_0"]), + "eta_s_0": rfloat(idict["eta_s_0"]), + "n": rfloat(idict["n"]), + "Z": rfloat(idict["Z"]), + "molar_mass": rfloat(idict["molar_mass"], 1.0e-3), + } + + # electronic terms + if np.abs(rfloat(idict["bel_0"])) > 1.0e-10: + mbr_params[name]["equation_of_state"] = "slb3-conductive" + mbr_params[name]["bel_0"] = rfloat(idict["bel_0"]) + mbr_params[name]["gel"] = rfloat(idict["gel"]) + + # The positivity constraint for T_crit below + # is hardcoded into HeFESTo. + if idict["T_crit"] != "0.00000" and float(idict["T_crit"]) > 0.0: + if name in magnetic_structural_parameters: + # The magnetic model in SLB2024 is identical + # to the chs model already coded in BurnMan + # but with the zero Gibbs energy and entropy + # assigned to the fully ordered form, + # not the disordered form,. + 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) + fT0 = (1.0 / A) * float(idict["T_crit"]) * 79.0 / (140.0 * p) + E0 = Scrit * fT0 + + mbr_mods[name] = [ + [ + "magnetic_chs", + { + "structural_parameter": p, + "curie_temperature": [float(idict["T_crit"]), 0.0], + "magnetic_moment": [magnetic_moment, 0.0], + }, + ], + [ + "linear", + { + "delta_E": E0, + "delta_S": Scrit, + "delta_V": 0.0, + }, + ], + ] + else: + mbr_mods[name] = [ + [ + "landau_slb_2022", + { + "Tc_0": float(idict["T_crit"]), + "S_D": float(idict["S_crit"]), + "V_D": float(idict["V_crit"]) / 1.0e6, + }, + ] + ] + mbr_names[name] = ( + full_name.lower() + .replace("-", "_") + .replace(" ", "_") + .replace("'", "") + .replace('\\"', "") + ) + +# WRITE FILE + +print( + "# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit\n" + "# for the Earth and Planetary Sciences\n" + "# Copyright (C) 2012 - 2024 by the BurnMan team,\n" + "# released under the GNU GPL v2 or later.\n" + "\n\n" + '"""\n' + "SLB_2024\n" + "Minerals from Stixrude & Lithgow-Bertelloni 2024 and references therein\n" + f"File autogenerated using HeFESTO_2024_to_burnman.py and {hefesto_path}\n" + '"""\n\n' + "from __future__ import absolute_import\n" + "\n\n" + "from ..eos.slb import SLB3\n" + "from ..eos import debye\n" + "from ..eos import birch_murnaghan as bm\n" + "from ..classes.mineral import Mineral\n" + "from ..classes.solution import Solution, RelaxedSolution\n" + "from ..classes.solutionmodel import (\n" + " IdealSolution,\n" + " SymmetricRegularSolution,\n" + " AsymmetricRegularSolution,\n" + ")\n\n" +) + + +st_class = ( + "class SLB3Stishovite(SLB3):\n" + " def shear_modulus(self, pressure, temperature, volume, params):\n" + ' """\n' + " Returns the shear modulus. :math:`[Pa]`\n" + " SLB2024 first calculate the shear modulus at T_0 using\n" + " the third order Birch-Murnaghan equation of state,\n" + " then apply the softening relevant to the temperature\n" + " and pressure of interest, and then apply the usual thermal term.\n" + ' """\n' + "\n" + " G_bare = bm.shear_modulus_third_order(volume, params)\n" + " Pc = 51.6e9 + 11.1 * (temperature - 300.0) * 1.0e6 # Eq B4\n" + " PcQ = (\n" + " Pc + 50.7e9\n" + " ) # From the HeFESTo code in stishtran.f. PcQ is T-independent in paper.\n" + " Cs0 = 128.0e9 * PcQ / Pc\n" + " Delta = 0.0 # 100.e9 in the HeFESTo code in stishtran.f, 0. in paper.\n" + " PminusPc = pressure - Pc\n" + " if PminusPc < 0.0:\n" + " Cs = Cs0 * PminusPc / (pressure - PcQ) # Eq B2\n" + " else:\n" + " Ast = (Cs0 - Delta) * (PcQ - Pc)\n" + " Cs = Cs0 - Ast / (PcQ - pressure + 3.0 * PminusPc) # Eq B3\n" + "\n" + " # VRH-like average of the bare shear modulus with softened (C11-C12)/2.\n" + " G = 0.5 * (13.0 / 15.0 * G_bare + 2.0 / 15.0 * Cs) + 0.5 / (\n" + " 13.0 / 15.0 / G_bare + 2.0 / 15.0 / Cs\n" + " ) # Eq B5\n" + "\n" + " debye_T = self._debye_temperature(params['V_0'] / volume, params)\n" + " eta_s = self._isotropic_eta_s(params['V_0'] / volume, params)\n" + "\n" + " E_th = debye.thermal_energy(temperature, debye_T, params['n'])\n" + " E_th_ref = debye.thermal_energy(params['T_0'], debye_T, params['n'])\n" + "\n" + " return G - eta_s * (E_th - E_th_ref) / volume\n" +) + +fper_relaxed_class = ( + "class ferropericlase_relaxed(RelaxedSolution):\n" + " def __init__(self):\n" + ' """RelaxedSolution model for ferropericlase (mw).\n' + " Endmembers (and site species distributions) are given in the order:\n" + " - pe ([Mg]2[Mg]2)\n" + " - wu and wuls ([Fe,Fels]2[Fe,Fels]2)\n" + " - anao ([Na]2[Al]2)\n" + " - mag ([Fe]0[Fef]2)\n" + " The entropy from the first site in magnetite is not counted.\n" + " The equilibrium spin state is calculated automatically.\n" + ' """\n' + " solution = ferropericlase()\n" + " vrel = [[0.0, 0.5, -0.5, 0.0, 0.0]]\n" + " vunrel = [\n" + " [1.0, 0.0, 0.0, 0.0, 0.0],\n" + " [0.0, 0.5, 0.5, 0.0, 0.0],\n" + " [0.0, 0.0, 0.0, 1.0, 0.0],\n" + " [0.0, 0.0, 0.0, 0.0, 1.0],\n" + " ]\n" + " 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,Fefls])\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": + print(st_class) + mbr_prm = mbr_prm.__repr__().replace("'slb3'", "SLB3Stishovite()") + mbr_prm = mbr_prm.replace("slb3", "fault") + print( + f"\nclass {key}(Mineral):\n" + " def __init__(self):\n" + f" self.params = {mbr_prm}\n" + ) + if key in mbr_mods: + print( + f" self.property_modifiers = {mbr_mods[key]}", + ) + print("") + print(" Mineral.__init__(self)") + print("") + +print('"""\n' "SOLUTIONS\n" '"""\n') + +for key, prm in sorted(sol_params.items()): + for i in range(prm["n_mbrs"]): + if prm["mbr_names"][i] in hardcoded_site_occupancies: + prm["mbr_site_formulae"][i] = hardcoded_site_occupancies[ + prm["mbr_names"][i] + ] + + docstring = '"""' + docstring += f'{prm["solution_type"]} model for {solution_aliases[key]} ({key}).\n' + docstring += "Endmembers (and site species distributions) are given in the order:\n" + for i in range(prm["n_mbrs"]): + docstring += f'- {prm["mbr_names"][i]} ({prm["mbr_site_formulae"][i]})\n' + if key == "mw": + docstring += "The entropy from the first site in magnetite is not counted.\n" + docstring += '"""' + + print( + f"\nclass {solution_aliases[key]}(Solution):\n" + " def __init__(self, molar_fractions=None):\n" + f" {docstring}\n" + f' self.name = "{solution_aliases[key]}"\n' + f" self.solution_model = {prm['solution_type']}(\n" + " endmembers=[" + ) + for i in range(prm["n_mbrs"]): + print( + f' [{prm["mbr_names"][i]}(), "{prm["mbr_site_formulae"][i]}"],' + ) + print(" ],") + + if prm["solution_type"] == "AsymmetricRegularSolution": + print(f" alphas={prm['alphas']},") + + if prm["solution_type"] != "IdealSolution": + print(f' energy_interaction={prm["energy_interaction"]},') + if "volume_interaction" in prm: + print(f' volume_interaction={prm["volume_interaction"]},') + + print(" )\n") + print(" Solution.__init__(self, molar_fractions=molar_fractions)\n") + + if key == "mw": + print(fper_relaxed_class) + if key == "pv": + print(bdg_relaxed_class) + + +print('\n"""\n' "ENDMEMBER ALIASES\n" '"""\n') + +for key, value in sorted(mbr_names.items()): + print(f"{value.replace('(', '').replace(')', '')} = {key}") + +print("") + + +print('\n"""\n' "SOLUTION ALIASES\n" '"""\n') + +for key, value in sorted(solution_aliases.items()): + if key != "sp": + print(f"{key} = {value}") diff --git a/burnman/data/input_raw_endmember_datasets/hefesto_parameter_names.txt b/burnman/data/input_raw_endmember_datasets/hefesto_parameter_names.txt index 7a0a49328..265707ed3 100644 --- a/burnman/data/input_raw_endmember_datasets/hefesto_parameter_names.txt +++ b/burnman/data/input_raw_endmember_datasets/hefesto_parameter_names.txt @@ -25,13 +25,13 @@ unknown Upper Limit of Optic Continuum unknown Lower Limit of Optic Continuum gamma_0 gamma_0 q_0 q_0 -beta beta (electronic contribution) -gammael_0 gammael_0 (electronic contribution) -q_2A_2 q_2A_2 (anharmonic contribution) +bel_0 beta (electronic contribution) +gel gammael_0 (electronic contribution) +q_2A_2 q_2A_2 (anharmonic contribution) unknown High Temperature Approximation (1=yes) unknown Birch-Murnaghan (0) or Vinet (1) -debye Einstein (0) or Debye (1) -zpp Zero Point Pressure (1=yes) +debye Einstein (0) or Debye (1) +zpp Zero Point Pressure (1=yes) G_0 Ambient Shear Modulus Gprime_0 Pressure Derivative eta_s_0 Temperature Derivative diff --git a/burnman/eos/__init__.py b/burnman/eos/__init__.py index 8b2e16115..3fe86584d 100644 --- a/burnman/eos/__init__.py +++ b/burnman/eos/__init__.py @@ -17,7 +17,7 @@ from .birch_murnaghan import BM2, BM3 from .birch_murnaghan_4th import BM4 from .mie_grueneisen_debye import MGD2, MGD3 -from .slb import SLB2, SLB3 +from .slb import SLB2, SLB3, SLB3Conductive from .modified_tait import MT from .hp import HP_TMT from .hp import HP_TMTL diff --git a/burnman/eos/bukowinski_electronic.py b/burnman/eos/bukowinski_electronic.py new file mode 100644 index 000000000..825fa63d4 --- /dev/null +++ b/burnman/eos/bukowinski_electronic.py @@ -0,0 +1,63 @@ +# 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. + +from __future__ import absolute_import +import numpy as np + +""" +Functions for the Bukowinski (1977) model for the electronic component +of the Helmholtz energy as used by Stixrude an Lithgow-Bertelloni (2024). +""" + + +def helmholtz(temperature, volume, T_0, V_0, bel_0, gel): + return ( + -0.5 + * bel_0 + * np.power(volume / V_0, gel) + * (temperature * temperature - T_0 * T_0) + ) + + +def pressure(temperature, volume, T_0, V_0, bel_0, gel): + """ + P = -dF/dV + """ + return ( + 0.5 + * gel + * bel_0 + * np.power(volume / V_0, gel) + * (temperature * temperature - T_0 * T_0) + / volume + ) + + +def KToverV(temperature, volume, T_0, V_0, bel_0, gel): + """ + KT = -V dP/dV + """ + return -(gel - 1.0) * pressure(temperature, volume, T_0, V_0, bel_0, gel) / volume + + +def entropy(temperature, volume, V_0, bel_0, gel): + """ + S = -dF/dT + """ + return bel_0 * np.power(volume / V_0, gel) * temperature + + +def CVoverT(volume, V_0, bel_0, gel): + """ + CV = T dS/dT + """ + return bel_0 * np.power(volume / V_0, gel) + + +def aKT(temperature, volume, V_0, bel_0, gel): + """ + aKT = dP/dT + """ + return gel * bel_0 * np.power(volume / V_0, gel) * temperature / volume diff --git a/burnman/eos/helper.py b/burnman/eos/helper.py index 3a9178dda..235f1f668 100644 --- a/burnman/eos/helper.py +++ b/burnman/eos/helper.py @@ -58,6 +58,8 @@ def create(method): return mgd.MGD3() elif method == "slb3": return slb.SLB3() + elif method == "slb3-conductive": + return slb.SLB3Conductive() elif method == "murnaghan": return murnaghan.Murnaghan() elif method == "bm2": diff --git a/burnman/eos/slb.py b/burnman/eos/slb.py index 78e04c0fc..094384199 100644 --- a/burnman/eos/slb.py +++ b/burnman/eos/slb.py @@ -1,6 +1,6 @@ # This file is part of BurnMan - a thermoelastic and thermodynamic toolkit # for the Earth and Planetary Sciences -# Copyright (C) 2012 - 2017 by the BurnMan team, released under the GNU +# Copyright (C) 2012 - 2024 by the BurnMan team, released under the GNU # GPL v2 or later. from __future__ import absolute_import @@ -23,11 +23,12 @@ def jit(fn): from . import birch_murnaghan as bm from . import debye from . import equation_of_state as eos +from . import bukowinski_electronic as el from ..utils.math import bracket @jit(nopython=True) -def _grueneisen_parameter_fast(V_0, volume, gruen_0, q_0): +def _grueneisen_parameter_slb(V_0, volume, gruen_0, q_0): """global function with plain parameters so jit will work""" x = V_0 / volume f = 1.0 / 2.0 * (pow(x, 2.0 / 3.0) - 1.0) @@ -39,7 +40,19 @@ def _grueneisen_parameter_fast(V_0, volume, gruen_0, q_0): @jit(nopython=True) def _delta_pressure( - x, pressure, temperature, V_0, T_0, Debye_0, n, a1_ii, a2_iikk, b_iikk, b_iikkmm + x, + pressure, + temperature, + V_0, + T_0, + Debye_0, + n, + a1_ii, + a2_iikk, + b_iikk, + b_iikkmm, + bel_0, + gel, ): f = 0.5 * (pow(V_0 / x, 2.0 / 3.0) - 1.0) nu_o_nu0_sq = 1.0 + a1_ii * f + 1.0 / 2.0 * a2_iikk * f * f @@ -53,11 +66,21 @@ def _delta_pressure( nu_o_nu0_sq = 1.0 + a1_ii * f + (1.0 / 2.0) * a2_iikk * f * f # EQ 41 gr = 1.0 / 6.0 / nu_o_nu0_sq * (2.0 * f + 1.0) * (a1_ii + a2_iikk * f) + Pel = ( + 0.5 + * gel + * bel_0 + * np.power(x / V_0, gel) + * (temperature * temperature - T_0 * T_0) + / x + ) + return ( (1.0 / 3.0) * (pow(1.0 + 2.0 * f, 5.0 / 2.0)) * ((b_iikk * f) + (0.5 * b_iikkmm * f * f)) + gr * (E_th - E_th_ref) / x + + Pel - pressure ) # EQ 21 @@ -149,14 +172,6 @@ def _isotropic_eta_s(self, x, params): return eta_s - # calculate isotropic thermal pressure, see - # Matas et. al. (2007) eq B4 - def _thermal_pressure(self, T, V, params): - Debye_T = self._debye_temperature(params["V_0"] / V, params) - gr = self.grueneisen_parameter(0.0, T, V, params) # P not important - P_th = gr * debye.thermal_energy(T, Debye_T, params["n"]) / V - return P_th - def volume(self, pressure, temperature, params): """ Returns molar volume. :math:`[m^3]` @@ -177,6 +192,10 @@ def volume(self, pressure, temperature, params): b_iikk = 9.0 * params["K_0"] # EQ 28 b_iikkmm = 27.0 * params["K_0"] * (params["Kprime_0"] - 4.0) # EQ 29z + bel_0, gel = 0.0, 1.0 + if self.conductive: + bel_0, gel = params["bel_0"], params["gel"] + # Finding the volume at a given pressure requires a # root-finding scheme. Here we use brentq to find the root. @@ -193,6 +212,8 @@ def volume(self, pressure, temperature, params): a2_iikk, b_iikk, b_iikkmm, + bel_0, + gel, ) try: @@ -241,9 +262,10 @@ def pressure(self, temperature, volume, params): [Pa] """ debye_T = self._debye_temperature(params["V_0"] / volume, params) - gr = self.grueneisen_parameter( - 0.0, temperature, volume, params - ) # does not depend on pressure + gr = _grueneisen_parameter_slb( + params["V_0"], volume, params["grueneisen_0"], params["q_0"] + ) + # does not depend on pressure # thermal energy at temperature T E_th = debye.thermal_energy(temperature, debye_T, params["n"]) # thermal energy at reference temperature @@ -257,24 +279,26 @@ def pressure(self, temperature, volume, params): ) + gr * ( E_th - E_th_ref ) / volume # EQ 21 - + if self.conductive: + P += el.pressure( + temperature, + volume, + params["T_0"], + params["V_0"], + params["bel_0"], + params["gel"], + ) return P - def grueneisen_parameter(self, pressure, temperature, volume, params): - """ - Returns grueneisen parameter :math:`[unitless]` - """ - return _grueneisen_parameter_fast( - params["V_0"], volume, params["grueneisen_0"], params["q_0"] - ) - def isothermal_bulk_modulus_reuss(self, pressure, temperature, volume, params): """ Returns isothermal bulk modulus :math:`[Pa]` """ T_0 = params["T_0"] debye_T = self._debye_temperature(params["V_0"] / volume, params) - gr = self.grueneisen_parameter(pressure, temperature, volume, params) + gr = _grueneisen_parameter_slb( + params["V_0"], volume, params["grueneisen_0"], params["q_0"] + ) # thermal energy at temperature T E_th = debye.thermal_energy(temperature, debye_T, params["n"]) @@ -294,18 +318,17 @@ def isothermal_bulk_modulus_reuss(self, pressure, temperature, volume, params): - (pow(gr, 2.0) / volume) * (C_v * temperature - C_v_ref * T_0) ) + if self.conductive: + K += volume * el.KToverV( + temperature, + volume, + params["T_0"], + params["V_0"], + params["bel_0"], + params["gel"], + ) return K - def isentropic_bulk_modulus_reuss(self, pressure, temperature, volume, params): - """ - Returns adiabatic bulk modulus. :math:`[Pa]` - """ - K_T = self.isothermal_bulk_modulus_reuss(pressure, temperature, volume, params) - alpha = self.thermal_expansivity(pressure, temperature, volume, params) - gr = self.grueneisen_parameter(pressure, temperature, volume, params) - K_S = K_T * (1.0 + gr * alpha * temperature) - return K_S - def shear_modulus(self, pressure, temperature, volume, params): """ Returns shear modulus. :math:`[Pa]` @@ -335,47 +358,32 @@ def molar_heat_capacity_v(self, pressure, temperature, volume, params): Returns heat capacity at constant volume. :math:`[J/K/mol]` """ debye_T = self._debye_temperature(params["V_0"] / volume, params) - return debye.molar_heat_capacity_v(temperature, debye_T, params["n"]) + C_v = debye.molar_heat_capacity_v(temperature, debye_T, params["n"]) - def molar_heat_capacity_p(self, pressure, temperature, volume, params): - """ - Returns heat capacity at constant pressure. :math:`[J/K/mol]` - """ - alpha = self.thermal_expansivity(pressure, temperature, volume, params) - gr = self.grueneisen_parameter(pressure, temperature, volume, params) - C_v = self.molar_heat_capacity_v(pressure, temperature, volume, params) - C_p = C_v * (1.0 + gr * alpha * temperature) - return C_p + if self.conductive: + C_v += temperature * el.CVoverT( + volume, params["V_0"], params["bel_0"], params["gel"] + ) + return C_v def thermal_expansivity(self, pressure, temperature, volume, params): """ Returns thermal expansivity. :math:`[1/K]` """ - C_v = self.molar_heat_capacity_v(pressure, temperature, volume, params) - gr = self.grueneisen_parameter(pressure, temperature, volume, params) - K = self.isothermal_bulk_modulus_reuss(pressure, temperature, volume, params) - alpha = gr * C_v / K / volume - return alpha - - def gibbs_free_energy(self, pressure, temperature, volume, params): - """ - Returns the Gibbs free energy at the pressure and temperature - of the mineral [J/mol] - """ - G = ( - self.helmholtz_free_energy(pressure, temperature, volume, params) - + pressure * volume + debye_T = self._debye_temperature(params["V_0"] / volume, params) + C_v = debye.molar_heat_capacity_v(temperature, debye_T, params["n"]) + gr_slb = _grueneisen_parameter_slb( + params["V_0"], volume, params["grueneisen_0"], params["q_0"] ) - return G + K = self.isothermal_bulk_modulus_reuss(pressure, temperature, volume, params) + alpha = gr_slb * C_v / volume / K - def molar_internal_energy(self, pressure, temperature, volume, params): - """ - Returns the internal energy at the pressure and temperature - of the mineral [J/mol] - """ - return self.helmholtz_free_energy( - pressure, temperature, volume, params - ) + temperature * self.entropy(pressure, temperature, volume, params) + if self.conductive: + aKTel = el.aKT( + temperature, volume, params["V_0"], params["bel_0"], params["gel"] + ) + alpha += aKTel / K + return alpha def entropy(self, pressure, temperature, volume, params): """ @@ -384,19 +392,12 @@ def entropy(self, pressure, temperature, volume, params): """ Debye_T = self._debye_temperature(params["V_0"] / volume, params) S = debye.entropy(temperature, Debye_T, params["n"]) - return S - def enthalpy(self, pressure, temperature, volume, params): - """ - Returns the enthalpy at the pressure and temperature - of the mineral [J/mol] - """ - - return ( - self.helmholtz_free_energy(pressure, temperature, volume, params) - + temperature * self.entropy(pressure, temperature, volume, params) - + pressure * volume - ) + if self.conductive: + S += el.entropy( + temperature, volume, params["V_0"], params["bel_0"], params["gel"] + ) + return S def helmholtz_free_energy(self, pressure, temperature, volume, params): """ @@ -421,8 +422,89 @@ def helmholtz_free_energy(self, pressure, temperature, volume, params): + F_quasiharmonic ) + if self.conductive: + F += el.helmholtz( + temperature, + volume, + params["T_0"], + params["V_0"], + params["bel_0"], + params["gel"], + ) return F + # Derived properties from here + # These functions can use pressure as an argument, + # as pressure should have been determined self-consistently + # by the point at which these functions are called. + def grueneisen_parameter(self, pressure, temperature, volume, params): + """ + Returns grueneisen parameter :math:`[unitless]` + """ + if self.conductive: + temperature = max(temperature, 1.0e-6) + K_T = self.isothermal_bulk_modulus_reuss( + pressure, temperature, volume, params + ) + alpha = self.thermal_expansivity(pressure, temperature, volume, params) + C_v = self.molar_heat_capacity_v(pressure, temperature, volume, params) + return alpha * K_T * volume / C_v + else: + return _grueneisen_parameter_slb( + params["V_0"], volume, params["grueneisen_0"], params["q_0"] + ) + + def isentropic_bulk_modulus_reuss(self, pressure, temperature, volume, params): + """ + Returns adiabatic bulk modulus. :math:`[Pa]` + """ + K_T = self.isothermal_bulk_modulus_reuss(pressure, temperature, volume, params) + alpha = self.thermal_expansivity(pressure, temperature, volume, params) + gr = self.grueneisen_parameter(pressure, temperature, volume, params) + return K_T * (1.0 + gr * alpha * temperature) + + def molar_heat_capacity_p(self, pressure, temperature, volume, params): + """ + Returns heat capacity at constant pressure. :math:`[J/K/mol]` + """ + alpha = self.thermal_expansivity(pressure, temperature, volume, params) + gr = self.grueneisen_parameter(pressure, temperature, volume, params) + C_v = self.molar_heat_capacity_v(pressure, temperature, volume, params) + C_p = C_v * (1.0 + gr * alpha * temperature) + return C_p + + def molar_internal_energy(self, pressure, temperature, volume, params): + """ + Returns the internal energy at the pressure and temperature + of the mineral [J/mol] + """ + return self.helmholtz_free_energy( + pressure, temperature, volume, params + ) + temperature * self.entropy(pressure, temperature, volume, params) + + def gibbs_free_energy(self, pressure, temperature, volume, params): + """ + Returns the Gibbs free energy at the pressure and temperature + of the mineral [J/mol] + """ + G = ( + self.helmholtz_free_energy(pressure, temperature, volume, params) + + pressure * volume + ) + return G + + def enthalpy(self, pressure, temperature, volume, params): + """ + Returns the enthalpy at the pressure and temperature + of the mineral [J/mol] + """ + + return ( + self.helmholtz_free_energy(pressure, temperature, volume, params) + + temperature * self.entropy(pressure, temperature, volume, params) + + pressure * volume + ) + def validate_parameters(self, params): """ Check for existence and validity of the parameters @@ -445,6 +527,9 @@ def validate_parameters(self, params): # Now check all the required keys for the # thermal part of the EoS are in the dictionary expected_keys = ["molar_mass", "n", "Debye_0", "grueneisen_0", "q_0", "eta_s_0"] + if self.conductive: + expected_keys.extend(["bel_0", "gel"]) + for k in expected_keys: if k not in params: raise KeyError("params object missing parameter : " + k) @@ -470,11 +555,12 @@ class SLB3(SLBBase): """ SLB equation of state with third order finite strain expansion for the shear modulus (this should be preferred, as it is more thermodynamically - consistent.) + consistent). """ def __init__(self): self.order = 3 + self.conductive = False class SLB2(SLBBase): @@ -487,3 +573,16 @@ class SLB2(SLBBase): def __init__(self): self.order = 2 + self.conductive = False + + +class SLB3Conductive(SLBBase): + """ + SLB equation of state with third order finite strain expansion for the + shear modulus and a contribution to the Helmholtz free energy + that arises from the thermal excitation of electrons (Bukowinski, 1977). + """ + + def __init__(self): + self.order = 3 + self.conductive = True diff --git a/burnman/minerals/SLB_2024.py b/burnman/minerals/SLB_2024.py new file mode 100644 index 000000000..68b67d9af --- /dev/null +++ b/burnman/minerals/SLB_2024.py @@ -0,0 +1,2598 @@ +# 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. + + +""" +SLB_2024 +Minerals from Stixrude & Lithgow-Bertelloni 2024 and references therein +File autogenerated using HeFESTO_2024_to_burnman.py and HeFESTo_parameters_010123 +""" + +from __future__ import absolute_import + + +from ..eos.slb import SLB3 +from ..eos import debye +from ..eos import birch_murnaghan as bm +from ..classes.mineral import Mineral +from ..classes.solution import Solution, RelaxedSolution +from ..classes.solutionmodel import ( + IdealSolution, + SymmetricRegularSolution, + AsymmetricRegularSolution, +) + + +""" +ENDMEMBERS +""" + + +class ab(Mineral): + def __init__(self): + self.params = { + "name": "Albite", + "formula": {"Na": 1.0, "Al": 1.0, "Si": 3.0, "O": 8.0}, + "equation_of_state": "slb3", + "F_0": -3717909.79, + "V_0": 0.000100452, + "K_0": 59752590000.0, + "Kprime_0": 2.77846, + "Debye_0": 719.0831, + "grueneisen_0": 0.57877, + "q_0": 1.0, + "G_0": 36000000000.0, + "Gprime_0": 1.38571, + "eta_s_0": 1.02954, + "n": 13.0, + "Z": 4.0, + "molar_mass": 0.262223, + } + + Mineral.__init__(self) + + +class acm(Mineral): + def __init__(self): + self.params = { + "name": "Acmite", + "formula": {"Na": 1.0, "Fe": 1.0, "Si": 2.0, "O": 6.0}, + "equation_of_state": "slb3", + "F_0": -2419124.39, + "V_0": 6.4606e-05, + "K_0": 116100000000.0, + "Kprime_0": 4.4, + "Debye_0": 726.94116, + "grueneisen_0": 0.77467, + "q_0": 0.60142, + "G_0": 72454000000.0, + "Gprime_0": 1.71384, + "eta_s_0": 1.7391, + "n": 10.0, + "Z": 1.0, + "molar_mass": 0.231, + } + + self.property_modifiers = [ + ["landau_slb_2022", {"Tc_0": 4.0, "S_D": 14.89723, "V_D": 0.0}] + ] + + Mineral.__init__(self) + + +class al(Mineral): + def __init__(self): + self.params = { + "name": "Almandine", + "formula": {"Fe": 3.0, "Al": 2.0, "Si": 3.0, "O": 12.0}, + "equation_of_state": "slb3", + "F_0": -4930861.659999999, + "V_0": 0.00011543, + "K_0": 173896370000.0, + "Kprime_0": 4.91341, + "Debye_0": 741.38227, + "grueneisen_0": 1.06493, + "q_0": 1.42169, + "G_0": 96000000000.0, + "Gprime_0": 1.40927, + "eta_s_0": 2.09289, + "n": 20.0, + "Z": 4.0, + "molar_mass": 0.49776, + } + + self.property_modifiers = [ + ["landau_slb_2022", {"Tc_0": 7.5, "S_D": 40.14405, "V_D": 0.0}] + ] + + Mineral.__init__(self) + + +class alpv(Mineral): + def __init__(self): + self.params = { + "name": "Al-Perovskite", + "formula": {"Al": 2.0, "O": 3.0}, + "equation_of_state": "slb3", + "F_0": -1517031.0, + "V_0": 2.4944e-05, + "K_0": 242400000000.0, + "Kprime_0": 4.1, + "Debye_0": 856.18212, + "grueneisen_0": 1.54466, + "q_0": 0.83352, + "G_0": 169199620000.0, + "Gprime_0": 1.55703, + "eta_s_0": 2.27978, + "n": 5.0, + "Z": 4.0, + "molar_mass": 0.101961, + } + + Mineral.__init__(self) + + +class an(Mineral): + def __init__(self): + self.params = { + "name": "Anorthite", + "formula": {"Ca": 1.0, "Al": 2.0, "Si": 2.0, "O": 8.0}, + "equation_of_state": "slb3", + "F_0": -4012381.1300000004, + "V_0": 0.00010061, + "K_0": 84093390000.0, + "Kprime_0": 6.73397, + "Debye_0": 754.46887, + "grueneisen_0": 0.38512, + "q_0": 1.0, + "G_0": 39900000000.0, + "Gprime_0": 1.09129, + "eta_s_0": 1.63405, + "n": 13.0, + "Z": 2.0, + "molar_mass": 0.278211, + } + + Mineral.__init__(self) + + +class anao(Mineral): + def __init__(self): + self.params = { + "name": "alpha-NaO2_phase", + "formula": {"Na": 2.0, "Al": 2.0, "O": 4.0}, + "equation_of_state": "slb3", + "F_0": -2124562.13, + "V_0": 4.542e-05, + "K_0": 161143930000.0, + "Kprime_0": 3.90838, + "Debye_0": 744.51451, + "grueneisen_0": 1.45033, + "q_0": 1.5487, + "G_0": 108465250000.0, + "Gprime_0": 2.14668, + "eta_s_0": 0.77305, + "n": 8.0, + "Z": 1.0, + "molar_mass": 0.16394023, + } + + Mineral.__init__(self) + + +class andr(Mineral): + def __init__(self): + self.params = { + "name": "Grossular", + "formula": {"Ca": 3.0, "Fe": 2.0, "Si": 3.0, "O": 12.0}, + "equation_of_state": "slb3", + "F_0": -5414950.41, + "V_0": 0.00013199, + "K_0": 153580040000.0, + "Kprime_0": 4.78447, + "Debye_0": 750.98472, + "grueneisen_0": 1.04336, + "q_0": 1.42169, + "G_0": 89700000000.0, + "Gprime_0": 0.97013, + "eta_s_0": 2.84503, + "n": 20.0, + "Z": 4.0, + "molar_mass": 0.50817, + } + + self.property_modifiers = [ + ["landau_slb_2022", {"Tc_0": 11.7, "S_D": 29.79445, "V_D": 0.0}] + ] + + Mineral.__init__(self) + + +class apbo(Mineral): + def __init__(self): + self.params = { + "name": "alpha-PbO_2-SiO_2", + "formula": {"Si": 1.0, "O": 2.0}, + "equation_of_state": "slb3", + "F_0": -792997.39, + "V_0": 1.367e-05, + "K_0": 327162360000.0, + "Kprime_0": 4.0166, + "Debye_0": 1132.97205, + "grueneisen_0": 1.55723, + "q_0": 2.21141, + "G_0": 227412210000.0, + "Gprime_0": 1.77076, + "eta_s_0": 4.56617, + "n": 3.0, + "Z": 4.0, + "molar_mass": 0.060085, + } + + Mineral.__init__(self) + + +class appv(Mineral): + def __init__(self): + self.params = { + "name": "Al-Post-Perovskite", + "formula": {"Al": 2.0, "O": 3.0}, + "equation_of_state": "slb3", + "F_0": -1434367.18, + "V_0": 2.45e-05, + "K_0": 247740000000.0, + "Kprime_0": 4.0, + "Debye_0": 752.02929, + "grueneisen_0": 1.86524, + "q_0": 1.76454, + "G_0": 91897150000.0, + "Gprime_0": 1.81823, + "eta_s_0": 2.70624, + "n": 5.0, + "Z": 4.0, + "molar_mass": 0.101961, + } + + Mineral.__init__(self) + + +class capv(Mineral): + def __init__(self): + self.params = { + "name": "Ca-Perovskite", + "formula": {"Ca": 1.0, "Si": 1.0, "O": 3.0}, + "equation_of_state": "slb3", + "F_0": -1461201.28, + "V_0": 2.745e-05, + "K_0": 236000000000.0, + "Kprime_0": 3.9, + "Debye_0": 800.29043, + "grueneisen_0": 1.88997, + "q_0": 0.89608, + "G_0": 126000000000.0, + "Gprime_0": 1.6, + "eta_s_0": 1.35422, + "n": 5.0, + "Z": 1.0, + "molar_mass": 0.116164, + } + + Mineral.__init__(self) + + +class cats(Mineral): + def __init__(self): + self.params = { + "name": "Lime_Tschermak'", + "formula": {"Ca": 1.0, "Al": 2.0, "Si": 1.0, "O": 6.0}, + "equation_of_state": "slb3", + "F_0": -3119890.8000000003, + "V_0": 6.3574e-05, + "K_0": 113759900000.0, + "Kprime_0": 4.8061, + "Debye_0": 804.36068, + "grueneisen_0": 0.82288, + "q_0": 0.60142, + "G_0": 75337410000.0, + "Gprime_0": 1.71384, + "eta_s_0": 1.77238, + "n": 10.0, + "Z": 1.0, + "molar_mass": 0.218123, + } + + Mineral.__init__(self) + + +class cen(Mineral): + def __init__(self): + self.params = { + "name": "Clinoenstatite", + "formula": {"Mg": 2.0, "Si": 2.0, "O": 6.0}, + "equation_of_state": "slb3", + "F_0": -2904584.98, + "V_0": 6.25e-05, + "K_0": 113759900000.0, + "Kprime_0": 4.8061, + "Debye_0": 805.59286, + "grueneisen_0": 1.00921, + "q_0": 0.60142, + "G_0": 76810310000.0, + "Gprime_0": 1.71384, + "eta_s_0": 1.42387, + "n": 10.0, + "Z": 1.0, + "molar_mass": 0.2007774, + } + + Mineral.__init__(self) + + +class co(Mineral): + def __init__(self): + self.params = { + "name": "Corundum", + "formula": {"Al": 2.0, "O": 3.0}, + "equation_of_state": "slb3", + "F_0": -1582443.1400000001, + "V_0": 2.5577e-05, + "K_0": 252585720000.0, + "Kprime_0": 3.88671, + "Debye_0": 932.21586, + "grueneisen_0": 1.3081, + "q_0": 1.71245, + "G_0": 163200000000.0, + "Gprime_0": 1.64704, + "eta_s_0": 2.84761, + "n": 5.0, + "Z": 2.0, + "molar_mass": 0.101961, + } + + Mineral.__init__(self) + + +class coes(Mineral): + def __init__(self): + self.params = { + "name": "Coesite", + "formula": {"Si": 1.0, "O": 2.0}, + "equation_of_state": "slb3", + "F_0": -853834.92, + "V_0": 2.0657e-05, + "K_0": 103537990000.0, + "Kprime_0": 2.9007, + "Debye_0": 875.22323, + "grueneisen_0": 0.29043, + "q_0": 1.0, + "G_0": 61600000000.0, + "Gprime_0": 0.49686, + "eta_s_0": 2.75631, + "n": 3.0, + "Z": 16.0, + "molar_mass": 0.060085, + } + + Mineral.__init__(self) + + +class cppv(Mineral): + def __init__(self): + self.params = { + "name": "Cr-Post-Perovskite", + "formula": {"Cr": 2.0, "O": 3.0}, + "equation_of_state": "slb3", + "F_0": -1110023.23, + "V_0": 2.6949e-05, + "K_0": 247740000000.0, + "Kprime_0": 4.0, + "Debye_0": 755.01863, + "grueneisen_0": 1.64015, + "q_0": 1.76454, + "G_0": 187873500000.0, + "Gprime_0": 1.98845, + "eta_s_0": 3.16235, + "n": 5.0, + "Z": 2.0, + "molar_mass": 0.15199, + } + + self.property_modifiers = [ + ["landau_slb_2022", {"Tc_0": 5.0, "S_D": 23.05213, "V_D": 0.0}] + ] + + Mineral.__init__(self) + + +class crcf(Mineral): + def __init__(self): + self.params = { + "name": "Cr_Ca-Ferrite", + "formula": {"Mg": 1.0, "Cr": 2.0, "O": 4.0}, + "equation_of_state": "slb3", + "F_0": -1736682.29, + "V_0": 3.9421e-05, + "K_0": 185400000000.0, + "Kprime_0": 4.0, + "Debye_0": 684.9543, + "grueneisen_0": 1.56656, + "q_0": 1.0, + "G_0": 141453800000.0, + "Gprime_0": 1.93591, + "eta_s_0": 1.93943, + "n": 7.0, + "Z": 8.0, + "molar_mass": 0.19229, + } + + self.property_modifiers = [ + ["landau_slb_2022", {"Tc_0": 5.0, "S_D": 23.05213, "V_D": 0.0}] + ] + + Mineral.__init__(self) + + +class crpv(Mineral): + def __init__(self): + self.params = { + "name": "Cr-Perovskite", + "formula": {"Cr": 2.0, "O": 3.0}, + "equation_of_state": "slb3", + "F_0": -1186685.95, + "V_0": 2.8189e-05, + "K_0": 250565350000.0, + "Kprime_0": 4.13438, + "Debye_0": 758.1187, + "grueneisen_0": 1.54466, + "q_0": 0.83352, + "G_0": 152053560000.0, + "Gprime_0": 1.73254, + "eta_s_0": 2.83537, + "n": 5.0, + "Z": 2.0, + "molar_mass": 0.15199, + } + + self.property_modifiers = [ + ["landau_slb_2022", {"Tc_0": 5.0, "S_D": 23.05213, "V_D": 0.0}] + ] + + Mineral.__init__(self) + + +class di(Mineral): + def __init__(self): + self.params = { + "name": "Diopside", + "formula": {"Ca": 1.0, "Mg": 1.0, "Si": 2.0, "O": 6.0}, + "equation_of_state": "slb3", + "F_0": -3029614.78, + "V_0": 6.6039e-05, + "K_0": 113759900000.0, + "Kprime_0": 4.8061, + "Debye_0": 782.57306, + "grueneisen_0": 1.00921, + "q_0": 0.60142, + "G_0": 72700000000.0, + "Gprime_0": 1.71384, + "eta_s_0": 1.06175, + "n": 10.0, + "Z": 1.0, + "molar_mass": 0.2165504, + } + + Mineral.__init__(self) + + +class en(Mineral): + def __init__(self): + self.params = { + "name": "Enstatite", + "formula": {"Mg": 2.0, "Si": 2.0, "O": 6.0}, + "equation_of_state": "slb3", + "F_0": -2912202.5, + "V_0": 6.2676e-05, + "K_0": 107076810000.0, + "Kprime_0": 7.02751, + "Debye_0": 812.21227, + "grueneisen_0": 0.78477, + "q_0": 3.43847, + "G_0": 76800000000.0, + "Gprime_0": 1.54596, + "eta_s_0": 2.5045, + "n": 10.0, + "Z": 4.0, + "molar_mass": 0.2007774, + } + + Mineral.__init__(self) + + +class esk(Mineral): + def __init__(self): + self.params = { + "name": "Eskolaite", + "formula": {"Cr": 2.0, "O": 3.0}, + "equation_of_state": "slb3", + "F_0": -1206914.71, + "V_0": 2.8904e-05, + "K_0": 233340360000.0, + "Kprime_0": 4.01705, + "Debye_0": 766.73627, + "grueneisen_0": 1.15191, + "q_0": 2.22481, + "G_0": 123200000000.0, + "Gprime_0": 1.81492, + "eta_s_0": 2.55521, + "n": 5.0, + "Z": 2.0, + "molar_mass": 0.15199, + } + + self.property_modifiers = [ + ["landau_slb_2022", {"Tc_0": 305.5, "S_D": 23.05213, "V_D": 0.0}] + ] + + Mineral.__init__(self) + + +class fa(Mineral): + def __init__(self): + self.params = { + "name": "Fayalite", + "formula": {"Fe": 2.0, "Si": 1.0, "O": 4.0}, + "equation_of_state": "slb3", + "F_0": -1371708.3399999999, + "V_0": 4.629e-05, + "K_0": 136485580000.0, + "Kprime_0": 4.88157, + "Debye_0": 618.96116, + "grueneisen_0": 1.08388, + "q_0": 2.88055, + "G_0": 51220000000.0, + "Gprime_0": 0.85893, + "eta_s_0": 1.65937, + "n": 7.0, + "Z": 4.0, + "molar_mass": 0.203777, + } + + self.property_modifiers = [ + ["landau_slb_2022", {"Tc_0": 65.0, "S_D": 26.7627, "V_D": 0.0}] + ] + + Mineral.__init__(self) + + +class fapv(Mineral): + def __init__(self): + self.params = { + "name": "FeAlO3-Perovskite_HS", + "formula": {"Fe": 1.0, "Al": 1.0, "O": 3.0}, + "equation_of_state": "slb3", + "F_0": -1131674.78, + "V_0": 2.709e-05, + "K_0": 223325500000.0, + "Kprime_0": 4.1, + "Debye_0": 755.62079, + "grueneisen_0": 1.54466, + "q_0": 0.83352, + "G_0": 159881770000.0, + "Gprime_0": 1.73254, + "eta_s_0": 2.8548, + "n": 5.0, + "Z": 4.0, + "molar_mass": 0.1308249, + } + + self.property_modifiers = [ + ["landau_slb_2022", {"Tc_0": 5.0, "S_D": 14.89723, "V_D": 0.0}] + ] + + Mineral.__init__(self) + + +class fea(Mineral): + def __init__(self): + self.params = { + "name": "alpha_(bcc)_Iron", + "formula": {"Fe": 1.0}, + "equation_of_state": "slb3-conductive", + "F_0": 11.82, + "V_0": 7.093e-06, + "K_0": 163421140000.0, + "Kprime_0": 6.01732, + "Debye_0": 398.00486, + "grueneisen_0": 1.66729, + "q_0": 0.90658, + "G_0": 81470000000.0, + "Gprime_0": 1.9223, + "eta_s_0": 5.96207, + "n": 1.0, + "Z": 2.0, + "molar_mass": 0.05584515, + "bel_0": 0.00388, + "gel": 1.4796, + } + + self.property_modifiers = [ + [ + "magnetic_chs", + { + "structural_parameter": 0.4, + "curie_temperature": [1043.01, 0.0], + "magnetic_moment": [2.1199287186834, 0.0], + }, + ], + [ + "linear", + {"delta_E": 8932.739871162863, "delta_S": 9.46028, "delta_V": 0.0}, + ], + ] + + Mineral.__init__(self) + + +class fec2(Mineral): + def __init__(self): + self.params = { + "name": "HP-Clinoferrosilite", + "formula": {"Fe": 2.0, "Si": 2.0, "O": 6.0}, + "equation_of_state": "slb3", + "F_0": -2219723.5300000003, + "V_0": 6.38541e-05, + "K_0": 116025560000.0, + "Kprime_0": 6.23685, + "Debye_0": 691.84626, + "grueneisen_0": 1.12478, + "q_0": 0.20409, + "G_0": 79289860000.0, + "Gprime_0": 1.84119, + "eta_s_0": 1.14272, + "n": 10.0, + "Z": 1.0, + "molar_mass": 0.2638614, + } + + self.property_modifiers = [ + ["landau_slb_2022", {"Tc_0": 5.0, "S_D": 26.7627, "V_D": 0.0}] + ] + + Mineral.__init__(self) + + +class fecf(Mineral): + def __init__(self): + self.params = { + "name": "Fe-Ca-Ferrite", + "formula": {"Fe": 1.0, "Al": 2.0, "O": 4.0}, + "equation_of_state": "slb3", + "F_0": -1768312.12, + "V_0": 3.7216e-05, + "K_0": 213000000000.0, + "Kprime_0": 4.1, + "Debye_0": 715.88779, + "grueneisen_0": 1.56656, + "q_0": 1.0, + "G_0": 164115350000.0, + "Gprime_0": 1.93591, + "eta_s_0": 2.46745, + "n": 7.0, + "Z": 4.0, + "molar_mass": 0.173806, + } + + self.property_modifiers = [ + ["landau_slb_2022", {"Tc_0": 5.0, "S_D": 13.38135, "V_D": 0.0}] + ] + + Mineral.__init__(self) + + +class fee(Mineral): + def __init__(self): + self.params = { + "name": "epsilon_(hcp)_Iron", + "formula": {"Fe": 1.0}, + "equation_of_state": "slb3-conductive", + "F_0": 3428.48, + "V_0": 6.7651e-06, + "K_0": 165000000000.0, + "Kprime_0": 4.97, + "Debye_0": 379.74281, + "grueneisen_0": 2.13092, + "q_0": 1.02253, + "G_0": 91000000000.0, + "Gprime_0": 2.24804, + "eta_s_0": 5.20603, + "n": 1.0, + "Z": 2.0, + "molar_mass": 0.05584515, + "bel_0": 0.00411, + "gel": 1.6927, + } + + Mineral.__init__(self) + + +class feg(Mineral): + def __init__(self): + self.params = { + "name": "gamma_(fcc)_Iron", + "formula": {"Fe": 1.0}, + "equation_of_state": "slb3-conductive", + "F_0": 4675.83, + "V_0": 6.9883e-06, + "K_0": 165000000000.0, + "Kprime_0": 5.97, + "Debye_0": 285.19315, + "grueneisen_0": 1.84924, + "q_0": 1.02176, + "G_0": 91000000000.0, + "Gprime_0": 2.33869, + "eta_s_0": 6.50014, + "n": 1.0, + "Z": 2.0, + "molar_mass": 0.05584515, + "bel_0": 0.00375, + "gel": 1.4796, + } + + Mineral.__init__(self) + + +class feil(Mineral): + def __init__(self): + self.params = { + "name": "Fe-Akimotoite", + "formula": {"Fe": 1.0, "Si": 1.0, "O": 3.0}, + "equation_of_state": "slb3", + "F_0": -1044260.0800000001, + "V_0": 2.6854e-05, + "K_0": 210692630000.0, + "Kprime_0": 5.2154, + "Debye_0": 760.91558, + "grueneisen_0": 1.19328, + "q_0": 2.22481, + "G_0": 167069220000.0, + "Gprime_0": 1.81492, + "eta_s_0": 3.6318, + "n": 5.0, + "Z": 2.0, + "molar_mass": 0.131931, + } + + self.property_modifiers = [ + ["landau_slb_2022", {"Tc_0": 5.0, "S_D": 13.38135, "V_D": 0.0}] + ] + + Mineral.__init__(self) + + +class fepv(Mineral): + def __init__(self): + self.params = { + "name": "Fe-Perovskite", + "formula": {"Fe": 1.0, "Si": 1.0, "O": 3.0}, + "equation_of_state": "slb3", + "F_0": -1002026.07, + "V_0": 2.5321e-05, + "K_0": 270582550000.00003, + "Kprime_0": 4.01, + "Debye_0": 740.39231, + "grueneisen_0": 1.54466, + "q_0": 0.83352, + "G_0": 130058789999.99998, + "Gprime_0": 1.37279, + "eta_s_0": 2.08541, + "n": 5.0, + "Z": 4.0, + "molar_mass": 0.131931, + } + + self.property_modifiers = [ + ["landau_slb_2022", {"Tc_0": 5.0, "S_D": 13.38135, "V_D": 0.0}] + ] + + Mineral.__init__(self) + + +class feri(Mineral): + def __init__(self): + self.params = { + "name": "Fe-Ringwoodite", + "formula": {"Fe": 2.0, "Si": 1.0, "O": 4.0}, + "equation_of_state": "slb3", + "F_0": -1356599.78, + "V_0": 4.186e-05, + "K_0": 213409100000.0, + "Kprime_0": 4.22035, + "Debye_0": 651.49411, + "grueneisen_0": 1.26156, + "q_0": 2.39214, + "G_0": 92000000000.0, + "Gprime_0": 1.35412, + "eta_s_0": 1.76941, + "n": 7.0, + "Z": 2.0, + "molar_mass": 0.203777, + } + + self.property_modifiers = [ + ["landau_slb_2022", {"Tc_0": 11.8, "S_D": 26.7627, "V_D": 0.0}] + ] + + Mineral.__init__(self) + + +class fewa(Mineral): + def __init__(self): + self.params = { + "name": "Fe-Wadsleyite", + "formula": {"Fe": 2.0, "Si": 1.0, "O": 4.0}, + "equation_of_state": "slb3", + "F_0": -1356211.64, + "V_0": 4.28e-05, + "K_0": 168566850000.0, + "Kprime_0": 4.12302, + "Debye_0": 636.8306, + "grueneisen_0": 1.20498, + "q_0": 2.20831, + "G_0": 72000000000.0, + "Gprime_0": 1.50973, + "eta_s_0": 0.94487, + "n": 7.0, + "Z": 8.0, + "molar_mass": 0.203777, + } + + self.property_modifiers = [ + ["landau_slb_2022", {"Tc_0": 5.0, "S_D": 26.7627, "V_D": 0.0}] + ] + + Mineral.__init__(self) + + +class fnal(Mineral): + def __init__(self): + self.params = { + "name": "Fe-NAL_Phase", + "formula": {"Na": 1.0, "Fe": 2.0, "Al": 5.0, "Si": 1.0, "O": 12.0}, + "equation_of_state": "slb3", + "F_0": -5465247.03, + "V_0": 0.000112045, + "K_0": 204006720000.0, + "Kprime_0": 4.31789, + "Debye_0": 788.03574, + "grueneisen_0": 1.43147, + "q_0": 1.0, + "G_0": 152250130000.0, + "Gprime_0": 1.74226, + "eta_s_0": 2.73801, + "n": 21.0, + "Z": 1.0, + "molar_mass": 0.48966601, + } + + self.property_modifiers = [ + ["landau_slb_2022", {"Tc_0": 5.0, "S_D": 26.7627, "V_D": 0.0}] + ] + + Mineral.__init__(self) + + +class fo(Mineral): + def __init__(self): + self.params = { + "name": "Forsterite", + "formula": {"Mg": 2.0, "Si": 1.0, "O": 4.0}, + "equation_of_state": "slb3", + "F_0": -2055345.33, + "V_0": 4.3603e-05, + "K_0": 127955500000.0, + "Kprime_0": 4.21796, + "Debye_0": 809.1977, + "grueneisen_0": 0.9928, + "q_0": 2.10671, + "G_0": 81600000000.0, + "Gprime_0": 1.46257, + "eta_s_0": 2.29968, + "n": 7.0, + "Z": 4.0, + "molar_mass": 0.140695, + } + + Mineral.__init__(self) + + +class fppv(Mineral): + def __init__(self): + self.params = { + "name": "Fe-PostPerovskite", + "formula": {"Fe": 1.0, "Si": 1.0, "O": 3.0}, + "equation_of_state": "slb3", + "F_0": -988569.5, + "V_0": 2.52963e-05, + "K_0": 247740000000.0, + "Kprime_0": 4.0, + "Debye_0": 769.31113, + "grueneisen_0": 1.64015, + "q_0": 1.76454, + "G_0": 129500000000.0, + "Gprime_0": 1.40076, + "eta_s_0": 1.8495, + "n": 5.0, + "Z": 4.0, + "molar_mass": 0.131931, + } + + self.property_modifiers = [ + ["landau_slb_2022", {"Tc_0": 5.0, "S_D": 13.38135, "V_D": 0.0}] + ] + + Mineral.__init__(self) + + +class fs(Mineral): + def __init__(self): + self.params = { + "name": "Ferrosilite", + "formula": {"Fe": 2.0, "Si": 2.0, "O": 6.0}, + "equation_of_state": "slb3", + "F_0": -2224559.86, + "V_0": 6.5941e-05, + "K_0": 100545440000.0, + "Kprime_0": 7.87556, + "Debye_0": 677.91886, + "grueneisen_0": 0.7144, + "q_0": 3.43847, + "G_0": 52000000000.0, + "Gprime_0": 1.54596, + "eta_s_0": 1.08228, + "n": 10.0, + "Z": 4.0, + "molar_mass": 0.2638614, + } + + self.property_modifiers = [ + ["landau_slb_2022", {"Tc_0": 5.0, "S_D": 26.7627, "V_D": 0.0}] + ] + + Mineral.__init__(self) + + +class gr(Mineral): + def __init__(self): + self.params = { + "name": "Grossular", + "formula": {"Ca": 3.0, "Al": 2.0, "Si": 3.0, "O": 12.0}, + "equation_of_state": "slb3", + "F_0": -6276906.68, + "V_0": 0.00012512, + "K_0": 167062260000.0, + "Kprime_0": 3.91544, + "Debye_0": 822.77062, + "grueneisen_0": 1.05402, + "q_0": 1.88886, + "G_0": 109000000000.0, + "Gprime_0": 1.16274, + "eta_s_0": 2.38415, + "n": 20.0, + "Z": 4.0, + "molar_mass": 0.450449, + } + + Mineral.__init__(self) + + +class hc(Mineral): + def __init__(self): + self.params = { + "name": "Hercynite", + "formula": {"Fe": 1.0, "Al": 2.0, "O": 4.0}, + "equation_of_state": "slb3", + "F_0": -1841813.73, + "V_0": 4.0843e-05, + "K_0": 208947270000.0, + "Kprime_0": 4.62702, + "Debye_0": 747.13664, + "grueneisen_0": 1.18794, + "q_0": 3.97087, + "G_0": 84500000000.0, + "Gprime_0": 0.62795, + "eta_s_0": 2.46339, + "n": 7.0, + "Z": 8.0, + "molar_mass": 0.17381, + } + + self.property_modifiers = [ + ["landau_slb_2022", {"Tc_0": 13.0, "S_D": 13.38134, "V_D": 0.0}] + ] + + Mineral.__init__(self) + + +class he(Mineral): + def __init__(self): + self.params = { + "name": "Hedenbergite", + "formula": {"Ca": 1.0, "Fe": 1.0, "Si": 2.0, "O": 6.0}, + "equation_of_state": "slb3", + "F_0": -2675636.34, + "V_0": 6.7867e-05, + "K_0": 119204720000.0, + "Kprime_0": 4.81927, + "Debye_0": 702.08234, + "grueneisen_0": 0.96665, + "q_0": 0.60142, + "G_0": 61000000000.0, + "Gprime_0": 1.71384, + "eta_s_0": 1.01745, + "n": 10.0, + "Z": 1.0, + "molar_mass": 0.2480924, + } + + self.property_modifiers = [ + ["landau_slb_2022", {"Tc_0": 34.5, "S_D": 13.38135, "V_D": 0.0}] + ] + + Mineral.__init__(self) + + +class hem(Mineral): + def __init__(self): + self.params = { + "name": "Hematite", + "formula": {"Fe": 2.0, "O": 3.0}, + "equation_of_state": "slb3", + "F_0": -744045.2000000001, + "V_0": 3.0287e-05, + "K_0": 204245370000.0, + "Kprime_0": 4.0997, + "Debye_0": 653.80747, + "grueneisen_0": 1.58944, + "q_0": 2.22481, + "G_0": 91000000000.0, + "Gprime_0": 1.81492, + "eta_s_0": 0.5241, + "n": 5.0, + "Z": 1.0, + "molar_mass": 0.15968852, + } + + self.property_modifiers = [ + ["landau_slb_2022", {"Tc_0": 950.01, "S_D": 29.79445, "V_D": 0.0}] + ] + + Mineral.__init__(self) + + +class hepv(Mineral): + def __init__(self): + self.params = { + "name": "Fe2O3-Perovskite_HS", + "formula": {"Fe": 2.0, "O": 3.0}, + "equation_of_state": "slb3", + "F_0": -706527.0299999999, + "V_0": 2.95768e-05, + "K_0": 204251000000.0, + "Kprime_0": 4.1, + "Debye_0": 646.79863, + "grueneisen_0": 1.54466, + "q_0": 0.83352, + "G_0": 123483250000.0, + "Gprime_0": 1.73254, + "eta_s_0": 1.88876, + "n": 5.0, + "Z": 4.0, + "molar_mass": 0.15968852, + } + + self.property_modifiers = [ + ["landau_slb_2022", {"Tc_0": 5.0, "S_D": 29.79445, "V_D": 0.0}] + ] + + Mineral.__init__(self) + + +class hlpv(Mineral): + def __init__(self): + self.params = { + "name": "Fe2O3-Perovskite_LS", + "formula": {"Fe": 2.0, "O": 3.0}, + "equation_of_state": "slb3", + "F_0": -597158.86, + "V_0": 2.75209e-05, + "K_0": 204251000000.0, + "Kprime_0": 4.1, + "Debye_0": 759.63863, + "grueneisen_0": 1.54466, + "q_0": 0.83352, + "G_0": 177577950000.0, + "Gprime_0": 1.73254, + "eta_s_0": 3.54218, + "n": 5.0, + "Z": 4.0, + "molar_mass": 0.15968852, + } + + self.property_modifiers = [ + ["landau_slb_2022", {"Tc_0": 5.0, "S_D": 20.66026, "V_D": 0.0}] + ] + + Mineral.__init__(self) + + +class hmag(Mineral): + def __init__(self): + self.params = { + "name": "High-Pressure_Magnetit", + "formula": {"Fe": 3.0, "O": 4.0}, + "equation_of_state": "slb3", + "F_0": -993123.89, + "V_0": 4.1702e-05, + "K_0": 172000000000.0, + "Kprime_0": 4.0, + "Debye_0": 542.9312, + "grueneisen_0": 1.56656, + "q_0": 0.41872, + "G_0": 120889150000.0, + "Gprime_0": 1.93591, + "eta_s_0": 1.37608, + "n": 7.0, + "Z": 1.0, + "molar_mass": 0.23153307, + } + + self.property_modifiers = [ + [ + "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) + + +class hppv(Mineral): + def __init__(self): + self.params = { + "name": "HS_Fe2O3-Post-Perovski", + "formula": {"Fe": 2.0, "O": 3.0}, + "equation_of_state": "slb3", + "F_0": -620839.7200000001, + "V_0": 2.76884e-05, + "K_0": 176500000000.0, + "Kprime_0": 4.0, + "Debye_0": 680.92363, + "grueneisen_0": 1.64015, + "q_0": 1.76454, + "G_0": 172363480000.0, + "Gprime_0": 1.98845, + "eta_s_0": 2.56327, + "n": 5.0, + "Z": 4.0, + "molar_mass": 0.15968852, + } + + self.property_modifiers = [ + ["landau_slb_2022", {"Tc_0": 5.0, "S_D": 29.79445, "V_D": 0.0}] + ] + + Mineral.__init__(self) + + +class jd(Mineral): + def __init__(self): + self.params = { + "name": "Jadeite", + "formula": {"Na": 1.0, "Al": 1.0, "Si": 2.0, "O": 6.0}, + "equation_of_state": "slb3", + "F_0": -2855160.68, + "V_0": 6.0508e-05, + "K_0": 137357290000.0, + "Kprime_0": 3.6305, + "Debye_0": 820.2985, + "grueneisen_0": 0.85734, + "q_0": 2.05453, + "G_0": 84000000000.0, + "Gprime_0": 1.71384, + "eta_s_0": 1.91023, + "n": 10.0, + "Z": 1.0, + "molar_mass": 0.2021387, + } + + Mineral.__init__(self) + + +class knor(Mineral): + def __init__(self): + self.params = { + "name": "Knorringite", + "formula": {"Mg": 3.0, "Cr": 2.0, "Si": 3.0, "O": 12.0}, + "equation_of_state": "slb3", + "F_0": -5523022.38, + "V_0": 0.000118716, + "K_0": 157000000000.0, + "Kprime_0": 4.5, + "Debye_0": 776.39637, + "grueneisen_0": 1.24672, + "q_0": 1.42169, + "G_0": 99527660000.0, + "Gprime_0": 1.35756, + "eta_s_0": 2.11433, + "n": 20.0, + "Z": 4.0, + "molar_mass": 0.45316, + } + + self.property_modifiers = [ + ["landau_slb_2022", {"Tc_0": 5.1, "S_D": 23.05213, "V_D": 0.0}] + ] + + Mineral.__init__(self) + + +class ky(Mineral): + def __init__(self): + self.params = { + "name": "Kyanite", + "formula": {"Al": 2.0, "Si": 1.0, "O": 5.0}, + "equation_of_state": "slb3", + "F_0": -2443738.48, + "V_0": 4.4227e-05, + "K_0": 160000000000.0, + "Kprime_0": 4.0, + "Debye_0": 943.19593, + "grueneisen_0": 0.92549, + "q_0": 1.0, + "G_0": 117646620000.0, + "Gprime_0": 1.69117, + "eta_s_0": 2.89863, + "n": 8.0, + "Z": 4.0, + "molar_mass": 0.1620456, + } + + Mineral.__init__(self) + + +class lppv(Mineral): + def __init__(self): + self.params = { + "name": "LS_Fe2O3-Post-Perovski", + "formula": {"Fe": 2.0, "O": 3.0}, + "equation_of_state": "slb3", + "F_0": -405645.45, + "V_0": 2.62534e-05, + "K_0": 176500000000.0, + "Kprime_0": 4.0, + "Debye_0": 713.13099, + "grueneisen_0": 1.64029, + "q_0": 1.76443, + "G_0": 220101260000.0, + "Gprime_0": 1.98843, + "eta_s_0": 3.79074, + "n": 5.0, + "Z": 4.0, + "molar_mass": 0.15968852, + } + + self.property_modifiers = [ + ["landau_slb_2022", {"Tc_0": 10.0, "S_D": 20.66026, "V_D": 0.0}] + ] + + Mineral.__init__(self) + + +class mag(Mineral): + def __init__(self): + self.params = { + "name": "Magnetite", + "formula": {"Fe": 3.0, "O": 4.0}, + "equation_of_state": "slb3", + "F_0": -1016009.99, + "V_0": 4.4528e-05, + "K_0": 183876640000.0, + "Kprime_0": 5.22819, + "Debye_0": 529.46966, + "grueneisen_0": 1.35821, + "q_0": 0.41872, + "G_0": 60300000000.0, + "Gprime_0": 0.04657, + "eta_s_0": 1.12258, + "n": 7.0, + "Z": 1.0, + "molar_mass": 0.23153307, + } + + self.property_modifiers = [ + [ + "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) + + +class mgc2(Mineral): + def __init__(self): + self.params = { + "name": "HP-Clinoenstatite", + "formula": {"Mg": 2.0, "Si": 2.0, "O": 6.0}, + "equation_of_state": "slb3", + "F_0": -2904465.49, + "V_0": 6.076e-05, + "K_0": 116025560000.0, + "Kprime_0": 6.23685, + "Debye_0": 824.44051, + "grueneisen_0": 1.12478, + "q_0": 0.20409, + "G_0": 87927170000.0, + "Gprime_0": 1.84119, + "eta_s_0": 2.14193, + "n": 10.0, + "Z": 1.0, + "molar_mass": 0.2007774, + } + + Mineral.__init__(self) + + +class mgcf(Mineral): + def __init__(self): + self.params = { + "name": "Mg-Ca-Ferrite", + "formula": {"Mg": 1.0, "Al": 2.0, "O": 4.0}, + "equation_of_state": "slb3", + "F_0": -2122885.3600000003, + "V_0": 3.6135e-05, + "K_0": 213000000000.0, + "Kprime_0": 4.1, + "Debye_0": 830.714, + "grueneisen_0": 1.56656, + "q_0": 1.0, + "G_0": 129699999999.99998, + "Gprime_0": 1.93591, + "eta_s_0": 1.30292, + "n": 7.0, + "Z": 4.0, + "molar_mass": 0.142266, + } + + Mineral.__init__(self) + + +class mgil(Mineral): + def __init__(self): + self.params = { + "name": "Mg-Akimotoite", + "formula": {"Mg": 1.0, "Si": 1.0, "O": 3.0}, + "equation_of_state": "slb3", + "F_0": -1409229.53, + "V_0": 2.6354e-05, + "K_0": 210692630000.0, + "Kprime_0": 5.2154, + "Debye_0": 928.95623, + "grueneisen_0": 1.19328, + "q_0": 2.22481, + "G_0": 132000000000.0, + "Gprime_0": 1.81492, + "eta_s_0": 3.3993, + "n": 5.0, + "Z": 2.0, + "molar_mass": 0.100389, + } + + Mineral.__init__(self) + + +class mgmj(Mineral): + def __init__(self): + self.params = { + "name": "Mg-Majorite", + "formula": {"Mg": 4.0, "Si": 4.0, "O": 12.0}, + "equation_of_state": "slb3", + "F_0": -5690008.32, + "V_0": 0.000114324, + "K_0": 165118370000.0, + "Kprime_0": 4.21183, + "Debye_0": 822.48562, + "grueneisen_0": 0.97681, + "q_0": 1.53581, + "G_0": 85000000000.0, + "Gprime_0": 1.42969, + "eta_s_0": 1.01779, + "n": 20.0, + "Z": 4.0, + "molar_mass": 0.40156, + } + + Mineral.__init__(self) + + +class mgpv(Mineral): + def __init__(self): + self.params = { + "name": "Mg-Perovskite", + "formula": {"Mg": 1.0, "Si": 1.0, "O": 3.0}, + "equation_of_state": "slb3", + "F_0": -1365338.1199999999, + "V_0": 2.4445e-05, + "K_0": 250565350000.0, + "Kprime_0": 4.13438, + "Debye_0": 892.95164, + "grueneisen_0": 1.54466, + "q_0": 0.83352, + "G_0": 172900000000.0, + "Gprime_0": 1.73254, + "eta_s_0": 1.65233, + "n": 5.0, + "Z": 4.0, + "molar_mass": 0.100389, + } + + Mineral.__init__(self) + + +class mgri(Mineral): + def __init__(self): + self.params = { + "name": "Mg-Ringwoodite", + "formula": {"Mg": 2.0, "Si": 1.0, "O": 4.0}, + "equation_of_state": "slb3", + "F_0": -2018381.1, + "V_0": 3.9493e-05, + "K_0": 184901750000.0, + "Kprime_0": 4.22035, + "Debye_0": 879.84656, + "grueneisen_0": 1.10843, + "q_0": 2.39214, + "G_0": 123000000000.0, + "Gprime_0": 1.35412, + "eta_s_0": 2.30588, + "n": 7.0, + "Z": 2.0, + "molar_mass": 0.140693, + } + + Mineral.__init__(self) + + +class mgts(Mineral): + def __init__(self): + self.params = { + "name": "Mg-Tschermak's", + "formula": {"Mg": 1.0, "Al": 2.0, "Si": 1.0, "O": 6.0}, + "equation_of_state": "slb3", + "F_0": -3000959.4699999997, + "V_0": 5.914e-05, + "K_0": 107076810000.0, + "Kprime_0": 7.02751, + "Debye_0": 788.01368, + "grueneisen_0": 0.78477, + "q_0": 3.43847, + "G_0": 93345240000.0, + "Gprime_0": 1.54596, + "eta_s_0": 2.39272, + "n": 10.0, + "Z": 4.0, + "molar_mass": 0.20235, + } + + Mineral.__init__(self) + + +class mgwa(Mineral): + def __init__(self): + self.params = { + "name": "Mg-Wadsleyite", + "formula": {"Mg": 2.0, "Si": 1.0, "O": 4.0}, + "equation_of_state": "slb3", + "F_0": -2029078.6, + "V_0": 4.0515e-05, + "K_0": 168702760000.0, + "Kprime_0": 4.12302, + "Debye_0": 849.12535, + "grueneisen_0": 1.20498, + "q_0": 2.20831, + "G_0": 112000000000.0, + "Gprime_0": 1.50973, + "eta_s_0": 2.56411, + "n": 7.0, + "Z": 8.0, + "molar_mass": 0.140695, + } + + Mineral.__init__(self) + + +class mnal(Mineral): + def __init__(self): + self.params = { + "name": "Mg-NAL_Phase", + "formula": {"Na": 1.0, "Mg": 2.0, "Al": 5.0, "Si": 1.0, "O": 12.0}, + "equation_of_state": "slb3", + "F_0": -6172514.21, + "V_0": 0.000109883, + "K_0": 204006720000.0, + "Kprime_0": 4.31789, + "Debye_0": 868.59088, + "grueneisen_0": 1.43147, + "q_0": 1.0, + "G_0": 129000000000.0, + "Gprime_0": 1.74226, + "eta_s_0": 1.9384, + "n": 21.0, + "Z": 1.0, + "molar_mass": 0.42658581, + } + + Mineral.__init__(self) + + +class mppv(Mineral): + def __init__(self): + self.params = { + "name": "Mg-PostPerovskite", + "formula": {"Mg": 1.0, "Si": 1.0, "O": 3.0}, + "equation_of_state": "slb3", + "F_0": -1332544.16, + "V_0": 2.41108e-05, + "K_0": 247740000000.0, + "Kprime_0": 4.0, + "Debye_0": 931.02549, + "grueneisen_0": 1.64015, + "q_0": 1.76454, + "G_0": 166644930000.0, + "Gprime_0": 1.98845, + "eta_s_0": 1.43815, + "n": 5.0, + "Z": 4.0, + "molar_mass": 0.100389, + } + + Mineral.__init__(self) + + +class nacf(Mineral): + def __init__(self): + self.params = { + "name": "Na-Ca-Ferrite", + "formula": {"Na": 1.0, "Al": 1.0, "Si": 1.0, "O": 4.0}, + "equation_of_state": "slb3", + "F_0": -1848182.13, + "V_0": 3.627e-05, + "K_0": 220000000000.0, + "Kprime_0": 4.1, + "Debye_0": 709.33152, + "grueneisen_0": 1.56656, + "q_0": 1.0, + "G_0": 135023869999.99998, + "Gprime_0": 1.74135, + "eta_s_0": 1.67725, + "n": 7.0, + "Z": 4.0, + "molar_mass": 0.142054, + } + + Mineral.__init__(self) + + +class namj(Mineral): + def __init__(self): + self.params = { + "name": "Na-Majorite", + "formula": {"Na": 2.0, "Mg": 1.0, "Si": 5.0, "O": 12.0}, + "equation_of_state": "slb3", + "F_0": -5309561.23, + "V_0": 0.000110842, + "K_0": 172035520000.0, + "Kprime_0": 5.20045, + "Debye_0": 845.23671, + "grueneisen_0": 1.25087, + "q_0": 0.10909, + "G_0": 114700000000.0, + "Gprime_0": 1.35756, + "eta_s_0": 2.48526, + "n": 20.0, + "Z": 4.0, + "molar_mass": 0.40270437, + } + + Mineral.__init__(self) + + +class neph(Mineral): + def __init__(self): + self.params = { + "name": "Nepheline", + "formula": {"Na": 1.0, "Al": 1.0, "Si": 1.0, "O": 4.0}, + "equation_of_state": "slb3", + "F_0": -1996488.39, + "V_0": 5.38684e-05, + "K_0": 53055500000.0, + "Kprime_0": 4.0, + "Debye_0": 743.57985, + "grueneisen_0": 0.6969, + "q_0": 1.0, + "G_0": 30700000000.0, + "Gprime_0": 1.33087, + "eta_s_0": 0.6241, + "n": 7.0, + "Z": 2.0, + "molar_mass": 0.14205431, + } + + self.property_modifiers = [ + [ + "landau_slb_2022", + {"Tc_0": 467.0, "S_D": 10.0, "V_D": 8.000000000000001e-07}, + ] + ] + + Mineral.__init__(self) + + +class nnal(Mineral): + def __init__(self): + self.params = { + "name": "Na-NAL_Phase", + "formula": {"Na": 3.0, "Al": 3.0, "Si": 3.0, "O": 12.0}, + "equation_of_state": "slb3", + "F_0": -5570853.449999999, + "V_0": 0.000109401, + "K_0": 204006720000.0, + "Kprime_0": 4.31789, + "Debye_0": 846.08425, + "grueneisen_0": 1.43147, + "q_0": 1.0, + "G_0": 144164270000.0, + "Gprime_0": 1.74226, + "eta_s_0": 2.40665, + "n": 21.0, + "Z": 1.0, + "molar_mass": 0.42616294, + } + + Mineral.__init__(self) + + +class odi(Mineral): + def __init__(self): + self.params = { + "name": "Ortho-Diopside", + "formula": {"Ca": 1.0, "Mg": 1.0, "Si": 2.0, "O": 6.0}, + "equation_of_state": "slb3", + "F_0": -3015701.23, + "V_0": 6.8054e-05, + "K_0": 107076810000.0, + "Kprime_0": 7.02751, + "Debye_0": 744.48915, + "grueneisen_0": 0.78477, + "q_0": 3.43847, + "G_0": 58247330000.0, + "Gprime_0": 1.54596, + "eta_s_0": 1.35202, + "n": 10.0, + "Z": 1.0, + "molar_mass": 0.2165504, + } + + Mineral.__init__(self) + + +class pe(Mineral): + def __init__(self): + self.params = { + "name": "Periclase", + "formula": {"Mg": 4.0, "O": 4.0}, + "equation_of_state": "slb3", + "F_0": -2278119.6300000004, + "V_0": 4.4976e-05, + "K_0": 161143930000.0, + "Kprime_0": 3.90838, + "Debye_0": 770.90151, + "grueneisen_0": 1.45033, + "q_0": 1.5487, + "G_0": 130900000000.0, + "Gprime_0": 2.14668, + "eta_s_0": 2.56123, + "n": 8.0, + "Z": 1.0, + "molar_mass": 0.16121782, + } + + Mineral.__init__(self) + + +class picr(Mineral): + def __init__(self): + self.params = { + "name": "Pircochromite", + "formula": {"Mg": 1.0, "Cr": 2.0, "O": 4.0}, + "equation_of_state": "slb3", + "F_0": -1808028.0, + "V_0": 4.3564e-05, + "K_0": 184400000000.0, + "Kprime_0": 5.7, + "Debye_0": 750.72523, + "grueneisen_0": 0.99168, + "q_0": 3.97087, + "G_0": 90290490000.0, + "Gprime_0": 0.62795, + "eta_s_0": 3.07014, + "n": 7.0, + "Z": 8.0, + "molar_mass": 0.19229, + } + + self.property_modifiers = [ + ["landau_slb_2022", {"Tc_0": 12.55, "S_D": 23.05213, "V_D": 0.0}] + ] + + Mineral.__init__(self) + + +class pwo(Mineral): + def __init__(self): + self.params = { + "name": "Pseudo-Wollastonite", + "formula": {"Ca": 1.0, "Si": 1.0, "O": 3.0}, + "equation_of_state": "slb3", + "F_0": -1547095.09, + "V_0": 4.0272e-05, + "K_0": 86000000000.0, + "Kprime_0": 3.8, + "Debye_0": 703.00475, + "grueneisen_0": 0.95232, + "q_0": 1.0, + "G_0": 28991610000.0, + "Gprime_0": 0.77536, + "eta_s_0": 0.76187, + "n": 5.0, + "Z": 1.0, + "molar_mass": 0.116164, + } + + Mineral.__init__(self) + + +class py(Mineral): + def __init__(self): + self.params = { + "name": "Pyrope", + "formula": {"Mg": 3.0, "Al": 2.0, "Si": 3.0, "O": 12.0}, + "equation_of_state": "slb3", + "F_0": -5932095.930000001, + "V_0": 0.00011308, + "K_0": 170239640000.0, + "Kprime_0": 4.11067, + "Debye_0": 823.23783, + "grueneisen_0": 1.01422, + "q_0": 1.42169, + "G_0": 93700000000.0, + "Gprime_0": 1.35756, + "eta_s_0": 0.98186, + "n": 20.0, + "Z": 4.0, + "molar_mass": 0.40313, + } + + Mineral.__init__(self) + + +class qtz(Mineral): + def __init__(self): + self.params = { + "name": "Quartz", + "formula": {"Si": 1.0, "O": 2.0}, + "equation_of_state": "slb3", + "F_0": -858043.33, + "V_0": 2.24211e-05, + "K_0": 61425960000.0, + "Kprime_0": 19.78014, + "Debye_0": 883.46813, + "grueneisen_0": -0.03957, + "q_0": 1.0, + "G_0": 44857750000.0, + "Gprime_0": -0.04277, + "eta_s_0": 2.40464, + "n": 3.0, + "Z": 16.0, + "molar_mass": 0.060085, + } + + self.property_modifiers = [ + [ + "landau_slb_2022", + {"Tc_0": 847.0, "S_D": 5.76, "V_D": 1.3593599999999998e-06}, + ] + ] + + Mineral.__init__(self) + + +class smag(Mineral): + def __init__(self): + self.params = { + "name": "Magnetite", + "formula": {"Fe": 3.0, "O": 4.0}, + "equation_of_state": "slb3", + "F_0": -1016009.99, + "V_0": 4.4528e-05, + "K_0": 183876640000.0, + "Kprime_0": 5.22819, + "Debye_0": 529.46966, + "grueneisen_0": 1.35821, + "q_0": 0.41872, + "G_0": 60300000000.0, + "Gprime_0": 0.04657, + "eta_s_0": 1.12258, + "n": 7.0, + "Z": 1.0, + "molar_mass": 0.23153307, + } + + self.property_modifiers = [ + [ + "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) + + +class sp(Mineral): + def __init__(self): + self.params = { + "name": "Spinel", + "formula": {"Mg": 1.0, "Al": 2.0, "O": 4.0}, + "equation_of_state": "slb3", + "F_0": -2171789.32, + "V_0": 3.9762e-05, + "K_0": 195103490000.0, + "Kprime_0": 4.62702, + "Debye_0": 801.86054, + "grueneisen_0": 0.97405, + "q_0": 3.97087, + "G_0": 109000000000.0, + "Gprime_0": 0.62795, + "eta_s_0": 2.4035, + "n": 7.0, + "Z": 8.0, + "molar_mass": 0.142278, + } + + Mineral.__init__(self) + + +class SLB3Stishovite(SLB3): + def shear_modulus(self, pressure, temperature, volume, params): + """ + Returns the shear modulus. :math:`[Pa]` + SLB2024 first calculate the shear modulus at T_0 using + the third order Birch-Murnaghan equation of state, + then apply the softening relevant to the temperature + and pressure of interest, and then apply the usual thermal term. + """ + + G_bare = bm.shear_modulus_third_order(volume, params) + Pc = 51.6e9 + 11.1 * (temperature - 300.0) * 1.0e6 # Eq B4 + PcQ = ( + Pc + 50.7e9 + ) # From the HeFESTo code in stishtran.f. PcQ is T-independent in paper. + Cs0 = 128.0e9 * PcQ / Pc + Delta = 0.0 # 100.e9 in the HeFESTo code in stishtran.f, 0. in paper. + PminusPc = pressure - Pc + if PminusPc < 0.0: + Cs = Cs0 * PminusPc / (pressure - PcQ) # Eq B2 + else: + Ast = (Cs0 - Delta) * (PcQ - Pc) + Cs = Cs0 - Ast / (PcQ - pressure + 3.0 * PminusPc) # Eq B3 + + # VRH-like average of the bare shear modulus with softened (C11-C12)/2. + G = 0.5 * (13.0 / 15.0 * G_bare + 2.0 / 15.0 * Cs) + 0.5 / ( + 13.0 / 15.0 / G_bare + 2.0 / 15.0 / Cs + ) # Eq B5 + + debye_T = self._debye_temperature(params["V_0"] / volume, params) + eta_s = self._isotropic_eta_s(params["V_0"] / volume, params) + + E_th = debye.thermal_energy(temperature, debye_T, params["n"]) + E_th_ref = debye.thermal_energy(params["T_0"], debye_T, params["n"]) + + return G - eta_s * (E_th - E_th_ref) / volume + + +class st(Mineral): + def __init__(self): + self.params = { + "name": "Stishovite", + "formula": {"Si": 1.0, "O": 2.0}, + "equation_of_state": SLB3Stishovite(), + "F_0": -817124.2, + "V_0": 1.4017e-05, + "K_0": 305839630000.0, + "Kprime_0": 4.02918, + "Debye_0": 1096.06023, + "grueneisen_0": 1.55723, + "q_0": 2.21141, + "G_0": 250293920000.0, + "Gprime_0": 2.2962, + "eta_s_0": 5.39736, + "n": 3.0, + "Z": 2.0, + "molar_mass": 0.060085, + } + + Mineral.__init__(self) + + +class wo(Mineral): + def __init__(self): + self.params = { + "name": "Wollastonite", + "formula": {"Ca": 1.0, "Si": 1.0, "O": 3.0}, + "equation_of_state": "slb3", + "F_0": -1548541.46, + "V_0": 3.9901e-05, + "K_0": 93900000000.0, + "Kprime_0": 4.0, + "Debye_0": 713.58788, + "grueneisen_0": 1.05734, + "q_0": 1.0, + "G_0": 50300000000.0, + "Gprime_0": 1.23206, + "eta_s_0": 1.2711, + "n": 5.0, + "Z": 1.0, + "molar_mass": 0.116164, + } + + Mineral.__init__(self) + + +class wu(Mineral): + def __init__(self): + self.params = { + "name": 'W\\"ustite', + "formula": {"Fe": 4.0, "O": 4.0}, + "equation_of_state": "slb3", + "F_0": -962050.35, + "V_0": 4.9024e-05, + "K_0": 160700000000.0, + "Kprime_0": 4.0, + "Debye_0": 427.00102, + "grueneisen_0": 1.45033, + "q_0": 1.5487, + "G_0": 59000000000.0, + "Gprime_0": 1.44764, + "eta_s_0": 0.0596, + "n": 8.0, + "Z": 1.0, + "molar_mass": 0.28737822, + } + + self.property_modifiers = [ + ["landau_slb_2022", {"Tc_0": 191.0, "S_D": 53.5254, "V_D": 0.0}] + ] + + Mineral.__init__(self) + + +class wuls(Mineral): + def __init__(self): + self.params = { + "name": 'W\\"ustite_Low_Spin', + "formula": {"Fe": 4.0, "O": 4.0}, + "equation_of_state": "slb3", + "F_0": -609335.18, + "V_0": 4.33997e-05, + "K_0": 199700000000.0, + "Kprime_0": 4.0, + "Debye_0": 492.99392, + "grueneisen_0": 1.45033, + "q_0": 1.5487, + "G_0": 59000000000.0, + "Gprime_0": 1.44073, + "eta_s_0": -0.14773, + "n": 8.0, + "Z": 1.0, + "molar_mass": 0.28737822, + } + + Mineral.__init__(self) + + +""" +SOLUTIONS +""" + + +class c2c_pyroxene(Solution): + def __init__(self, molar_fractions=None): + """IdealSolution model for c2c_pyroxene (c2c). + Endmembers (and site species distributions) are given in the order: + - mgc2 ([Mg]2) + - fec2 ([Fe]2) + """ + self.name = "c2c_pyroxene" + self.solution_model = IdealSolution( + endmembers=[ + [mgc2(), "[Mg]2"], + [fec2(), "[Fe]2"], + ], + ) + + Solution.__init__(self, molar_fractions=molar_fractions) + + +class calcium_ferrite_structured_phase(Solution): + def __init__(self, molar_fractions=None): + """AsymmetricRegularSolution model for calcium_ferrite_structured_phase (cf). + Endmembers (and site species distributions) are given in the order: + - mgcf ([Mg][Al][Al]) + - fecf ([Fe][Al][Al]) + - nacf ([Na][Al][Si]) + - hmag ([Fe][Fef][Fef]) + - crcf ([Mg][Cr][Cr]) + """ + self.name = "calcium_ferrite_structured_phase" + self.solution_model = AsymmetricRegularSolution( + endmembers=[ + [mgcf(), "[Mg][Al][Al]"], + [fecf(), "[Fe][Al][Al]"], + [nacf(), "[Na][Al][Si]"], + [hmag(), "[Fe][Fef][Fef]"], + [crcf(), "[Mg][Cr][Cr]"], + ], + alphas=[1.0, 1.0, 3.97647, 1.0, 1.0], + energy_interaction=[ + [0.0, 67194.52, 0.0, 0.0], + [67194.52, 0.0, 0.0], + [0.0, 67194.52], + [0.0], + ], + ) + + Solution.__init__(self, molar_fractions=molar_fractions) + + +class clinopyroxene(Solution): + def __init__(self, molar_fractions=None): + """AsymmetricRegularSolution model for clinopyroxene (cpx). + Endmembers (and site species distributions) are given in the order: + - di ([Ca][Mg][Si]2) + - he ([Ca][Fe][Si]2) + - cen ([Mg][Mg][Si]2) + - cats ([Ca][Al][Si1/2Al1/2]2) + - jd ([Na][Al][Si]2) + - acm ([Na][Fef][Si]2) + """ + self.name = "clinopyroxene" + self.solution_model = AsymmetricRegularSolution( + endmembers=[ + [di(), "[Ca][Mg][Si]2"], + [he(), "[Ca][Fe][Si]2"], + [cen(), "[Mg][Mg][Si]2"], + [cats(), "[Ca][Al][Si1/2Al1/2]2"], + [jd(), "[Na][Al][Si]2"], + [acm(), "[Na][Fef][Si]2"], + ], + alphas=[1.0, 1.0, 1.0, 3.5, 1.0, 1.0], + energy_interaction=[ + [0.0, 24740.0, 26000.0, 24300.0, 24300.0], + [24740.0, 26000.0, 24300.0, 24300.0], + [63718.12, 46037.09, 46037.09], + [10000.0, 10000.0], + [0.0], + ], + ) + + Solution.__init__(self, molar_fractions=molar_fractions) + + +class garnet(Solution): + def __init__(self, molar_fractions=None): + """SymmetricRegularSolution model for garnet (gt). + Endmembers (and site species distributions) are given in the order: + - py ([Mg]3[Al][Al]) + - al ([Fe]3[Al][Al]) + - gr ([Ca]3[Al][Al]) + - mgmj ([Mg]3[Mg][Si]) + - namj ([Na2/3Mg1/3]3[Si][Si]) + - andr ([Ca]3[Fef][Fef]) + - knor ([Mg]3[Cr][Cr]) + """ + self.name = "garnet" + self.solution_model = SymmetricRegularSolution( + endmembers=[ + [py(), "[Mg]3[Al][Al]"], + [al(), "[Fe]3[Al][Al]"], + [gr(), "[Ca]3[Al][Al]"], + [mgmj(), "[Mg]3[Mg][Si]"], + [namj(), "[Na2/3Mg1/3]3[Si][Si]"], + [andr(), "[Ca]3[Fef][Fef]"], + [knor(), "[Mg]3[Cr][Cr]"], + ], + energy_interaction=[ + [0.0, 19088.57, 23164.32, 23164.32, 53000.0, 0.0], + [0.0, 23164.32, 23164.32, 44000.0, 10000.0], + [64682.45, 64682.45, 0.0, -20000.0], + [71172.45, 64682.45, 25879.13], + [64682.45, 25879.13], + [75000.0], + ], + volume_interaction=[ + [0.0, 1.03e-06, 0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0], + [1.03e-06, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0], + [0.0, 0.0], + [0.0], + ], + ) + + Solution.__init__(self, molar_fractions=molar_fractions) + + +class ilmenite(Solution): + def __init__(self, molar_fractions=None): + """AsymmetricRegularSolution model for ilmenite (il). + Endmembers (and site species distributions) are given in the order: + - mgil ([Mg][Si]) + - feil ([Fe][Si]) + - co ([Al][Al]) + - hem ([Fef][Fef]) + - esk ([Cr][Cr]) + """ + self.name = "ilmenite" + self.solution_model = AsymmetricRegularSolution( + endmembers=[ + [mgil(), "[Mg][Si]"], + [feil(), "[Fe][Si]"], + [co(), "[Al][Al]"], + [hem(), "[Fef][Fef]"], + [esk(), "[Cr][Cr]"], + ], + alphas=[1.0, 1.0, 1.0, 1.3, 1.0], + energy_interaction=[ + [0.0, 62487.4, 65000.0, 93184.02], + [62487.4, 65000.0, 93184.02], + [65000.0, 40000.0], + [40000.0], + ], + ) + + Solution.__init__(self, molar_fractions=molar_fractions) + + +class ferropericlase(Solution): + def __init__(self, molar_fractions=None): + """AsymmetricRegularSolution model for ferropericlase (mw). + Endmembers (and site species distributions) are given in the order: + - pe ([Mg]2[Mg]2) + - wu ([Fe]2[Fe]2) + - wuls ([Fels]2[Fels]2) + - anao ([Na]2[Al]2) + - mag ([Fef]0[Fef]2) + The entropy from the first site in magnetite is not counted. + """ + self.name = "ferropericlase" + self.solution_model = AsymmetricRegularSolution( + endmembers=[ + [pe(), "[Mg]2[Mg]2"], + [wu(), "[Fe]2[Fe]2"], + [wuls(), "[Fels]2[Fels]2"], + [anao(), "[Na]2[Al]2"], + [mag(), "[Fef]0[Fef]2"], + ], + alphas=[1.0, 1.0, 1.0, 1.0, 0.08293], + energy_interaction=[ + [44000.0, -87120.47, 120000.0, 302745.2], + [-60219.09, 120000.0, 106151.49], + [120000.0, 106151.49], + [120000.0], + ], + volume_interaction=[ + [4.4e-07, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0], + [0.0, 0.0], + [0.0], + ], + ) + + Solution.__init__(self, molar_fractions=molar_fractions) + + +class ferropericlase_relaxed(RelaxedSolution): + def __init__(self): + """RelaxedSolution model for ferropericlase (mw). + Endmembers (and site species distributions) are given in the order: + - pe ([Mg]2[Mg]2) + - wu and wuls ([Fe,Fels]2[Fe,Fels]2) + - anao ([Na]2[Al]2) + - mag ([Fe]0[Fef]2) + The entropy from the first site in magnetite is not counted. + The equilibrium spin state is calculated automatically. + """ + solution = ferropericlase() + vrel = [[0.0, 0.5, -0.5, 0.0, 0.0]] + vunrel = [ + [1.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, 1.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 1.0], + ] + RelaxedSolution.__init__(self, solution, vrel, vunrel) + + +class new_aluminous_phase(Solution): + def __init__(self, molar_fractions=None): + """SymmetricRegularSolution model for new_aluminous_phase (nal). + Endmembers (and site species distributions) are given in the order: + - mnal ([Mg]2[Al5/6Si1/6]6) + - fnal ([Fe]2[Al5/6Si1/6]6) + - nnal ([Na]2[Al3/6Si3/6]6) + """ + self.name = "new_aluminous_phase" + self.solution_model = SymmetricRegularSolution( + endmembers=[ + [mnal(), "[Mg]2[Al5/6Si1/6]6"], + [fnal(), "[Fe]2[Al5/6Si1/6]6"], + [nnal(), "[Na]2[Al3/6Si3/6]6"], + ], + energy_interaction=[[0.0, -62082.92], [-62082.92]], + ) + + Solution.__init__(self, molar_fractions=molar_fractions) + + +class olivine(Solution): + def __init__(self, molar_fractions=None): + """SymmetricRegularSolution model for olivine (ol). + Endmembers (and site species distributions) are given in the order: + - fo ([Mg]2) + - fa ([Fe]2) + """ + self.name = "olivine" + self.solution_model = SymmetricRegularSolution( + endmembers=[ + [fo(), "[Mg]2"], + [fa(), "[Fe]2"], + ], + energy_interaction=[[2474.06]], + ) + + Solution.__init__(self, molar_fractions=molar_fractions) + + +class orthopyroxene(Solution): + def __init__(self, molar_fractions=None): + """SymmetricRegularSolution model for orthopyroxene (opx). + Endmembers (and site species distributions) are given in the order: + - en ([Mg][Mg]) + - fs ([Fe][Fe]) + - mgts ([Mg][Al]) + - odi ([Ca][Mg]) + """ + self.name = "orthopyroxene" + self.solution_model = SymmetricRegularSolution( + endmembers=[ + [en(), "[Mg][Mg]"], + [fs(), "[Fe][Fe]"], + [mgts(), "[Mg][Al]"], + [odi(), "[Ca][Mg]"], + ], + energy_interaction=[[0.0, 0.0, 32213.03], [0.0, 32213.03], [46640.38]], + ) + + Solution.__init__(self, molar_fractions=molar_fractions) + + +class plagioclase(Solution): + def __init__(self, molar_fractions=None): + """SymmetricRegularSolution model for plagioclase (plg). + Endmembers (and site species distributions) are given in the order: + - an ([Ca]) + - ab ([Na]) + """ + self.name = "plagioclase" + self.solution_model = SymmetricRegularSolution( + endmembers=[ + [an(), "[Ca]"], + [ab(), "[Na]"], + ], + energy_interaction=[[13000.0]], + ) + + Solution.__init__(self, molar_fractions=molar_fractions) + + +class post_perovskite(Solution): + def __init__(self, molar_fractions=None): + """SymmetricRegularSolution model for post_perovskite (ppv). + Endmembers (and site species distributions) are given in the order: + - mppv ([Mg][Si]) + - fppv ([Fe][Si]) + - appv ([Al][Al]) + - hppv ([Fef][Fef]) + - cppv ([Cr][Cr]) + """ + self.name = "post_perovskite" + self.solution_model = SymmetricRegularSolution( + endmembers=[ + [mppv(), "[Mg][Si]"], + [fppv(), "[Fe][Si]"], + [appv(), "[Al][Al]"], + [hppv(), "[Fef][Fef]"], + [cppv(), "[Cr][Cr]"], + ], + energy_interaction=[ + [-39831.77, 31662.95, 93867.53, 93184.02], + [31662.95, 93867.53, 93184.02], + [65000.0, 40000.0], + [40000.0], + ], + ) + + Solution.__init__(self, molar_fractions=molar_fractions) + + +class bridgmanite(Solution): + def __init__(self, molar_fractions=None): + """SymmetricRegularSolution model for bridgmanite (pv). + Endmembers (and site species distributions) are given in the order: + - mgpv ([Mg][Si]) + - fepv ([Fe][Si]) + - alpv ([Al][Al]) + - hepv ([Fef][Fef]) + - hlpv ([Fef][Fefls]) + - fapv ([Fef][Al]) + - crpv ([Cr][Cr]) + """ + self.name = "bridgmanite" + self.solution_model = SymmetricRegularSolution( + endmembers=[ + [mgpv(), "[Mg][Si]"], + [fepv(), "[Fe][Si]"], + [alpv(), "[Al][Al]"], + [hepv(), "[Fef][Fef]"], + [hlpv(), "[Fef][Fefls]"], + [fapv(), "[Fef][Al]"], + [crpv(), "[Cr][Cr]"], + ], + energy_interaction=[ + [-12467.39, 31662.95, 93867.53, 49859.96, 0.0, 93184.02], + [0.0, 93867.53, 49859.96, 0.0, 93184.02], + [65000.0, 65000.0, 65000.0, 40000.0], + [-5875.1, 65000.0, 40000.0], + [65000.0, 40000.0], + [40000.0], + ], + ) + + 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,Fefls]) + - 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). + Endmembers (and site species distributions) are given in the order: + - mgri ([Mg]2) + - feri ([Fe]2) + """ + self.name = "ringwoodite" + self.solution_model = SymmetricRegularSolution( + endmembers=[ + [mgri(), "[Mg]2"], + [feri(), "[Fe]2"], + ], + energy_interaction=[[6229.27]], + ) + + Solution.__init__(self, molar_fractions=molar_fractions) + + +class mg_fe_aluminous_spinel(Solution): + def __init__(self, molar_fractions=None): + """SymmetricRegularSolution model for mg_fe_aluminous_spinel (sp). + Endmembers (and site species distributions) are given in the order: + - sp ([Mg][Al][Al]) + - hc ([Fe][Al][Al]) + - smag ([Fe][Fef][Fef]) + - picr ([Mg][Cr][Cr]) + """ + self.name = "mg_fe_aluminous_spinel" + self.solution_model = SymmetricRegularSolution( + endmembers=[ + [sp(), "[Mg][Al][Al]"], + [hc(), "[Fe][Al][Al]"], + [smag(), "[Fe][Fef][Fef]"], + [picr(), "[Mg][Cr][Cr]"], + ], + energy_interaction=[ + [-1310.9, 63000.0, 20993.94], + [55000.0, 21000.0], + [42000.0], + ], + ) + + Solution.__init__(self, molar_fractions=molar_fractions) + + +class wadsleyite(Solution): + def __init__(self, molar_fractions=None): + """SymmetricRegularSolution model for wadsleyite (wa). + Endmembers (and site species distributions) are given in the order: + - mgwa ([Mg]2) + - fewa ([Fe]2) + """ + self.name = "wadsleyite" + self.solution_model = SymmetricRegularSolution( + endmembers=[ + [mgwa(), "[Mg]2"], + [fewa(), "[Fe]2"], + ], + energy_interaction=[[12868.75]], + ) + + Solution.__init__(self, molar_fractions=molar_fractions) + + +""" +ENDMEMBER ALIASES +""" + +albite = ab +acmite = acm +almandine = al +al_perovskite = alpv +anorthite = an +alpha_nao2_phase = anao +grossular = andr +alpha_pbo_2_sio_2 = apbo +al_post_perovskite = appv +ca_perovskite = capv +lime_tschermak = cats +clinoenstatite = cen +corundum = co +coesite = coes +cr_post_perovskite = cppv +cr_ca_ferrite = crcf +cr_perovskite = crpv +diopside = di +enstatite = en +eskolaite = esk +fayalite = fa +fealo3_perovskite_hs = fapv +alpha_bcc_iron = fea +hp_clinoferrosilite = fec2 +fe_ca_ferrite = fecf +epsilon_hcp_iron = fee +gamma_fcc_iron = feg +fe_akimotoite = feil +fe_perovskite = fepv +fe_ringwoodite = feri +fe_wadsleyite = fewa +fe_nal_phase = fnal +forsterite = fo +fe_postperovskite = fppv +ferrosilite = fs +grossular = gr +hercynite = hc +hedenbergite = he +hematite = hem +fe2o3_perovskite_hs = hepv +fe2o3_perovskite_ls = hlpv +high_pressure_magnetit = hmag +hs_fe2o3_post_perovski = hppv +jadeite = jd +knorringite = knor +kyanite = ky +ls_fe2o3_post_perovski = lppv +magnetite = mag +hp_clinoenstatite = mgc2 +mg_ca_ferrite = mgcf +mg_akimotoite = mgil +mg_majorite = mgmj +mg_perovskite = mgpv +mg_ringwoodite = mgri +mg_tschermaks = mgts +mg_wadsleyite = mgwa +mg_nal_phase = mnal +mg_postperovskite = mppv +na_ca_ferrite = nacf +na_majorite = namj +nepheline = neph +na_nal_phase = nnal +ortho_diopside = odi +periclase = pe +pircochromite = picr +pseudo_wollastonite = pwo +pyrope = py +quartz = qtz +magnetite = smag +spinel = sp +stishovite = st +wollastonite = wo +wustite = wu +wustite_low_spin = wuls + + +""" +SOLUTION ALIASES +""" + +c2c = c2c_pyroxene +cf = calcium_ferrite_structured_phase +cpx = clinopyroxene +gt = garnet +il = ilmenite +mw = ferropericlase +nal = new_aluminous_phase +ol = olivine +opx = orthopyroxene +plg = plagioclase +ppv = post_perovskite +pv = bridgmanite +ri = ringwoodite +wa = wadsleyite diff --git a/burnman/minerals/__init__.py b/burnman/minerals/__init__.py index 0e3caa660..85c16cc48 100644 --- a/burnman/minerals/__init__.py +++ b/burnman/minerals/__init__.py @@ -11,6 +11,7 @@ - :mod:`~burnman.minerals.SLB_2011_ZSB_2013` - :mod:`~burnman.minerals.SLB_2011` - :mod:`~burnman.minerals.SLB_2022` + - :mod:`~burnman.minerals.SLB_2024` - :mod:`~burnman.minerals.DKS_2013_liquids` - :mod:`~burnman.minerals.DKS_2013_solids` - :mod:`~burnman.minerals.RS_2014_liquids` @@ -27,6 +28,7 @@ from __future__ import absolute_import # Stixrude and Lithgow-Bertelloni +from . import SLB_2024 from . import SLB_2022 from . import SLB_2011 from . import SLB_2011_ZSB_2013 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 000000000..69a56c111 Binary files /dev/null and b/misc/benchmarks/figures/SLB_2024_Fe2O3_V.png differ diff --git a/misc/benchmarks/figures/SLB_2024_Fe_Cp_V.png b/misc/benchmarks/figures/SLB_2024_Fe_Cp_V.png new file mode 100644 index 000000000..eebdc6ab6 Binary files /dev/null and b/misc/benchmarks/figures/SLB_2024_Fe_Cp_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 000000000..e6261c158 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 000000000..8b0716d15 Binary files /dev/null and b/misc/benchmarks/figures/SLB_2024_Fe_O_fO2_T.png differ diff --git a/misc/benchmarks/figures/SLB_2024_Fe_O_phase_diagram.png b/misc/benchmarks/figures/SLB_2024_Fe_O_phase_diagram.png new file mode 100644 index 000000000..e63059a71 Binary files /dev/null and b/misc/benchmarks/figures/SLB_2024_Fe_O_phase_diagram.png differ diff --git a/misc/benchmarks/figures/SLB_2024_Fe_phase_diagram.png b/misc/benchmarks/figures/SLB_2024_Fe_phase_diagram.png new file mode 100644 index 000000000..ff7667e93 Binary files /dev/null and b/misc/benchmarks/figures/SLB_2024_Fe_phase_diagram.png differ diff --git a/misc/benchmarks/figures/SLB_2024_fper_K_T.png b/misc/benchmarks/figures/SLB_2024_fper_K_T.png new file mode 100644 index 000000000..04e6c5c32 Binary files /dev/null and b/misc/benchmarks/figures/SLB_2024_fper_K_T.png differ diff --git a/misc/benchmarks/figures/SLB_2024_fper_V.png b/misc/benchmarks/figures/SLB_2024_fper_V.png new file mode 100644 index 000000000..3bbf5153f Binary files /dev/null and b/misc/benchmarks/figures/SLB_2024_fper_V.png differ diff --git a/misc/benchmarks/figures/SLB_2024_fper_fO2.png b/misc/benchmarks/figures/SLB_2024_fper_fO2.png new file mode 100644 index 000000000..7137ce985 Binary files /dev/null and b/misc/benchmarks/figures/SLB_2024_fper_fO2.png differ diff --git a/misc/benchmarks/figures/SLB_2024_fper_spin_20_50_80.png b/misc/benchmarks/figures/SLB_2024_fper_spin_20_50_80.png new file mode 100644 index 000000000..4c8d9e315 Binary files /dev/null and b/misc/benchmarks/figures/SLB_2024_fper_spin_20_50_80.png differ diff --git a/misc/benchmarks/figures/SLB_2024_fper_spin_percents.png b/misc/benchmarks/figures/SLB_2024_fper_spin_percents.png new file mode 100644 index 000000000..5fffe700d Binary files /dev/null and b/misc/benchmarks/figures/SLB_2024_fper_spin_percents.png differ diff --git a/misc/benchmarks/figures/SLB_2024_stv_VS.png b/misc/benchmarks/figures/SLB_2024_stv_VS.png new file mode 100644 index 000000000..c2ffce005 Binary files /dev/null and b/misc/benchmarks/figures/SLB_2024_stv_VS.png differ diff --git a/misc/benchmarks/slb_2024_Fe_endmembers.dat b/misc/benchmarks/slb_2024_Fe_endmembers.dat new file mode 100644 index 000000000..c268b7cb8 --- /dev/null +++ b/misc/benchmarks/slb_2024_Fe_endmembers.dat @@ -0,0 +1,13 @@ + # id spec Pi Ti Enthalpy Entropy Cp Gibbs Volume Smix Sconf CV + 1 wu 0.00 833.00 -787.6607 460.0267 207.0322 -1170.8630 50.0110 6.8091 -0.0000 197.0969 414.6426 + 2 mag 0.00 833.00 -863.5328 345.8858 288.7717 -1151.6557 45.2057 18.1363 -0.0000 282.8666 518.6847 + 3 wu 0.00 833.00 -787.6607 460.0267 207.0322 -1170.8630 50.0110 11.5263 -0.0000 197.0969 414.6426 + 4 mag 0.00 833.00 -863.5328 345.8858 288.7717 -1151.6557 45.2057 11.5263 -0.0000 282.8666 518.6847 + 5 hem 0.00 833.00 -646.7135 221.0119 164.2131 -830.8164 30.7971 0.0000 -0.0000 158.3723 636.3544 + 6 hepv 0.00 833.00 -618.1526 232.5549 126.4491 -811.8709 30.0608 5.7632 -0.0000 121.2145 630.6739 + 7 hlpv 0.00 833.00 -515.1353 204.1707 125.4698 -685.2095 27.9927 5.7632 -0.0000 119.9336 739.8142 + 8 hppv 0.00 833.00 -533.6867 227.6978 128.9032 -723.3590 28.2999 0.0000 -0.0000 120.9283 656.4890 + 9 fea 0.00 833.00 25.0948 59.0890 39.5161 -24.1263 7.2463 0.0000 -0.0000 37.7495 383.9275 + 10 feg 0.00 833.00 30.4729 64.4276 30.2506 -23.1952 7.1594 0.0000 -0.0000 28.0479 272.5673 + 11 fee 0.00 833.00 27.4641 58.1046 31.3715 -20.9371 6.9617 0.0000 -0.0000 28.3096 356.9375 + 12 hmag 0.00 833.00 -840.9691 343.3957 291.6978 -1127.0177 42.5359 0.0000 -0.0000 282.7694 526.2834 diff --git a/misc/benchmarks/slb_2024_benchmarks.py b/misc/benchmarks/slb_2024_benchmarks.py new file mode 100644 index 000000000..ae6fd7fc1 --- /dev/null +++ b/misc/benchmarks/slb_2024_benchmarks.py @@ -0,0 +1,519 @@ +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 +from burnman.eos.property_modifiers import calculate_property_modifications +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.image as mpimg + + +def check_iron_bearing_mineral_property_values(): + print("Checking iron bearing mineral properties...") + m_names = np.genfromtxt("slb_2024_Fe_endmembers.dat", dtype=str)[:, 1] + minerals = [eval("SLB_2024." + m)() for m in m_names] + d = np.genfromtxt("slb_2024_Fe_endmembers.dat", dtype=float)[:, 2:] + for i, m in enumerate(minerals): + P, T, _, entropy, Cp, gibbs, volume, _, _, _, _ = d[i] + m.set_state(P * 1.0e4, T) + delta_gibbs = gibbs * 1.0e3 - m.gibbs + if np.abs(delta_gibbs) > 1.0: + print(m_names[i]) + print("delta Gibbs:", delta_gibbs) + print("delta S:", entropy - m.S) + print("delta Cp:", Cp - m.C_p) + print("delta V:", volume * 1.0e-6 - m.V) + else: + rel_error = max( + [ + (entropy - m.S) / entropy, + (Cp - m.C_p) / Cp, + (volume - m.V * 1.0e6) / volume, + ] + ) + print( + f"{m_names[i]} in agreement with HeFESTo (max rel error {rel_error*100.:.1g} %)" + ) + + +def check_bcc_iron_consistency(): + print("") + fe_bcc = SLB_2024.alpha_bcc_iron() + assert check_eos_consistency(fe_bcc, verbose=True) + + +def check_fper_entropy(): + print("\nChecking fper configurational entropy...") + fper = SLB_2024.ferropericlase() + print("Endmember entropies (should all be zero):") + print(np.abs(fper.solution_model.endmember_configurational_entropies)) + + fper.set_composition([0.0, 0.5, 0.0, 0.0, 0.5]) + Sxs = -2 * gas_constant * np.log(0.5) + print( + f"Equimolar wu-mag: {fper.excess_entropy:.5f} J/K/mol (should be {Sxs:.5f} J/K/mol)" + ) + + +def check_fig_1_fper_relaxed(): + print("\nChecking Figure 1...") + fper = SLB_2024.ferropericlase_relaxed() + c = molar_volume_from_unit_cell_volume(1.0, SLB_2024.periclase().params["Z"]) + + pressures = np.linspace(1.0e5, 140.0e9, 141) + temperatures = pressures * 0.0 + 300.0 + + fig = plt.figure(figsize=(12, 8)) + fig.suptitle("Figure 1 (Stixrude and Lithgow-Bertelloni, 2024)") + ax = [fig.add_subplot(2, 2, i) for i in range(1, 5)] + + fig1 = mpimg.imread("figures/SLB_2024_fper_V.png") + fig2 = mpimg.imread("figures/SLB_2024_fper_K_T.png") + fig3 = mpimg.imread("figures/SLB_2024_fper_spin_20_50_80.png") + fig4 = mpimg.imread("figures/SLB_2024_fper_spin_percents.png") + + ax[0].imshow(fig1, extent=[0.0, 140.0, 50.0, 80.0], aspect="auto") + ax[1].imshow(fig2, extent=[0.0, 140.0, 100.0, 700.0], aspect="auto") + ax[2].imshow(fig3, extent=[0.0, 1.0, 0.0, 100.0], aspect="auto") + ax[3].imshow(fig4, extent=[0.0, 200.0, 0.0, 4000.0], aspect="auto") + + for x_fe in [0.1, 0.3, 0.5, 1.0]: + fper.set_composition([1.0 - x_fe, x_fe, 0.0, 0.0]) + Vs, K_Ss = fper.evaluate(["V", "K_S"], pressures, temperatures) + + ax[0].plot(pressures / 1.0e9, Vs / c, linestyle=":", linewidth=3.0) + ax[1].plot(pressures / 1.0e9, K_Ss / 1.0e9, linestyle=":", linewidth=3.0) + + ax[0].set_ylim(50.0, 80.0) + ax[1].set_ylim(100.0, 700.0) + + fper.set_composition([0.5, 0.5, 0.0, 0.0]) + fper.set_state(70.0e9, 300.0) + print("Mg2Fe2O4 ferropericlase at 70 GPa, 300 K:") + fstr = ", ".join([f"{f:.2f}" for f in fper.molar_fractions]) + print(f"Molar proportions: {fstr}") + print(f"Volume: {fper.V*1.e6:.2f} cm^3/mol") + + fper = SLB_2024.ferropericlase() + assemblage = simplify_composite_with_composition( + Composite([fper]), {"Mg": 0.5, "Fe": 0.5, "O": 1.0} + ) + assemblage.set_state(50.0e9, 300.0) + fper = assemblage.phases[0] + x_MgOs = np.linspace(1.0e-5, 0.9999, 101) + pressures = np.empty((3, 101)) + for i, fLS in enumerate([0.2, 0.5, 0.8]): + for j, 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) / 2.0, (1.0 - x_MgO) / 2.0]) + equality_constraints = [ + ["T", 300.0], + [ + "phase_composition", + ( + fper, + ( + ["Fe_A", "Fels_A", "Fe_B", "Fels_B"], + [0.0, 1.0, 0.0, 1.0], + [1.0, 1.0, 1.0, 1.0], + fLS, + ), + ), + ], + ] + equilibrate(composition, assemblage, equality_constraints) + pressures[i, j] = assemblage.pressure + ax[2].plot( + 1.0 - x_MgOs, + (pressures[2] - pressures[0]) / 1.0e9, + linestyle=":", + linewidth=3.0, + ) + ax[2].plot(1.0 - x_MgOs, (pressures[1]) / 1.0e9, linestyle=":", linewidth=3.0) + + composition = {"Mg": 0.75, "Fe": 0.25, "O": 1.0} + fractions = np.linspace(0.1, 0.9, 9) + temperatures = np.linspace(1.0, 4000.0, 21) + equality_constraints = [ + ["T", temperatures], + [ + "phase_composition", + ( + fper, + ( + ["Fe_A", "Fels_A", "Fe_B", "Fels_B"], + [0.0, 1.0, 0.0, 1.0], + [1.0, 1.0, 1.0, 1.0], + fractions, + ), + ), + ], + ] + solss = equilibrate(composition, assemblage, equality_constraints)[0] + pressures = np.array([[sol.assemblage.pressure for sol in sols] for sols in solss]) + for i, f in enumerate(fractions): + ax[3].plot(pressures[:, i] / 1.0e9, temperatures, linestyle=":", linewidth=3.0) + ax[3].set_xlim(0.0, 200.0) + ax[3].set_ylim(0.0, 4000.0) + + plt.show() + + +def check_fig_3_fcc_ferric_fper(): + print("\nChecking Figure 3...") + + fig = plt.figure(figsize=(10, 5)) + fig.suptitle("Figure 3 (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_O_phase_diagram.png") + fig2 = mpimg.imread("figures/SLB_2024_fper_fO2.png") + + ax[0].imshow(fig1, extent=[-1, -0.75, 700.0, 1800.0], aspect="auto") + ax[1].imshow(fig2, extent=[0.0, 1.0, 0.0, 0.12], aspect="auto") + + # Subplot a + bcc = SLB_2024.alpha_bcc_iron() + fcc = SLB_2024.gamma_fcc_iron() + fper = SLB_2024.ferropericlase() + mag = SLB_2024.smag() + + assemblage = simplify_composite_with_composition( + Composite([fper]), {"Fe": 1.0, "O": 1.1} + ) + fper = assemblage.phases[0] + + composition = {"Fe": 1.0, "O": 1.0} + assemblage = Composite([bcc, fper, mag], [0.5, 0.49, 0.01]) + fper.set_composition([0.05, 0.0, 0.95]) + equality_constraints = [["P", 1.0e5], ["phase_fraction", (mag, 0.0)]] + sol = equilibrate(composition, assemblage, equality_constraints, tol=1.0e-5) + T_a_wu_mag = sol[0].assemblage.temperature + print(f"BCC-fper-mag triple point at 1 bar: {T_a_wu_mag:.2f} K") + + assemblage = Composite([bcc, fcc], [0.5, 0.5]) + equality_constraints = [["P", 1.0e5], ["phase_fraction", (bcc, 0.0)]] + sol = equilibrate(composition, assemblage, equality_constraints, tol=1.0e-5) + T_a_g = sol[0].assemblage.temperature + + assemblage = Composite([bcc, fper], [0.5, 0.5]) + temperatures = np.linspace(T_a_wu_mag, T_a_g, 31) + equality_constraints = [["P", 1.0e5], ["T", temperatures]] + sol = equilibrate(composition, assemblage, equality_constraints, tol=1.0e-5) + xs = np.empty_like(temperatures) + for i, s in enumerate(sol[0]): + xs[i] = -s.assemblage.phases[1].formula["Fe"] / 4.0 + ax[0].plot(xs, temperatures, linestyle=":", linewidth=3.0) + + assemblage = Composite([fcc, fper], [0.5, 0.5]) + temperatures = np.linspace(T_a_g, 1800.0, 31) + equality_constraints = [["P", 1.0e5], ["T", temperatures]] + sol = equilibrate(composition, assemblage, equality_constraints, tol=1.0e-5) + xs = np.empty_like(temperatures) + for i, s in enumerate(sol[0]): + xs[i] = -s.assemblage.phases[1].formula["Fe"] / 4.0 + ax[0].plot(xs, temperatures, linestyle=":", linewidth=3.0) + + assemblage = Composite([mag, fper], [0.5, 0.5]) + temperatures = np.linspace(T_a_wu_mag, 1800.0, 31) + equality_constraints = [["P", 1.0e5], ["T", temperatures]] + sol = equilibrate(composition, assemblage, equality_constraints, tol=1.0e-5) + xs = np.empty_like(temperatures) + for i, s in enumerate(sol[0]): + xs[i] = -s.assemblage.phases[1].formula["Fe"] / 4.0 + ax[0].plot(xs, temperatures, linestyle=":", linewidth=3.0) + + # Subplot b + fcc = SLB_2024.gamma_fcc_iron() + fper = SLB_2024.ferropericlase() + + x_MgOs = np.linspace(0.01, 0.99, 33) + x_Mgs = np.empty_like(x_MgOs) * np.nan + x_Fe3oversumFes = np.empty_like(x_MgOs) * np.nan + composition = {"Mg": 0.01, "Fe": 2.0, "O": 1.0} + assemblage = Composite([fper, fcc], [0.5, 0.5]) + x_MgO = 0.01 + fper.set_composition([x_MgO, 0.95 * (1.0 - x_MgO), 0.0, 0.0, 0.05 * (1.0 - x_MgO)]) + assemblage = simplify_composite_with_composition(assemblage, composition) + + 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": 2.0, "O": 1.0} + try: + sol = equilibrate( + composition, assemblage, equality_constraints, tol=1.0e-5 + ) + if sol[0].success: + f = assemblage.phases[0].formula + x_Mgs[i] = f["Mg"] / 4.0 + x_Fe3oversumFes[i] = 2.0 * ((f["O"] - f["Mg"]) / f["Fe"] - 1.0) + except AssertionError: + pass + ax[1].plot(x_Mgs, x_Fe3oversumFes, linestyle=":", linewidth=3.0) + plt.show() + + +def check_fig_6a_iron_Cp_V(): + print("\nChecking Figure 6a...") + fe_bcc = SLB_2024.alpha_bcc_iron() + fe_fcc = SLB_2024.gamma_fcc_iron() + + temperatures = np.linspace(1.0, 2000.0, 101) + 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") + ax[0].imshow(fig1, extent=[0.0, 2000.0, -15.0, 70.0], aspect="auto") + ax[1].imshow(fig1, extent=[0.0, 2000.0, 6.9, 8.4], aspect="auto") + + for m in [fe_bcc, fe_fcc]: + Cp, V = m.evaluate(["C_p", "V"], pressures, temperatures) + ax[0].plot(temperatures, Cp, linestyle=":", linewidth=3.0) + ax[1].plot(temperatures, V * 1.0e6, linestyle=":", linewidth=3.0) + + fe_bcc.set_state(1.0e5, 1000.0) + print( + f"BCC iron Cp, V at 1 bar, 1000 K: {fe_bcc.C_p:.2f} J/K/mol, {fe_bcc.V*1.e6:.2f} cm^3/mol" + ) + ax[0].set_ylim(0.0, 70.0) + ax[1].set_ylim(6.9, 7.8) + plt.show() + + +def check_fig_6c_iron_phase_diagram(): + print("\nChecking Figure 6c...") + 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=(8, 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") + + # 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) + equilibrate( + {"Fe": 1.0}, + irons, + equality_constraints=[ + ["phase_fraction", [fe_bcc, 0.0]], + ["phase_fraction", [fe_hcp, 0.0]], + ], + ) + P_invariant = irons.pressure + T_invariant = irons.temperature + + print( + f"The bcc-fcc-hcp invariant is at {P_invariant/1.e9:.2f} GPa, {T_invariant:.1f} K" + ) + + for Tmin, Tmax, m1, m2 in [ + [300.0, T_invariant, fe_bcc, fe_hcp], + [T_invariant, 2000.0, fe_bcc, fe_fcc], + [T_invariant, 4800.0, fe_fcc, fe_hcp], + ]: + + temperatures = np.linspace(Tmin, Tmax, 101) + irons = Composite([m1, m2], [1.0 / 2.0, 1.0 / 2.0]) + sols = equilibrate( + {"Fe": 1.0}, + irons, + equality_constraints=[["phase_fraction", [m1, 0.0]], ["T", temperatures]], + )[0] + Ps, Ts = np.array( + [ + [sol.assemblage.pressure, sol.assemblage.temperature] + for sol in sols + if sol.success + ] + ).T + ax[0].plot(Ps / 1.0e9, Ts, linestyle=":", linewidth=3.0) + + plt.show() + + +def check_fig_7_fO2(): + print("\nChecking Figure 7...") + 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=(8, 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, tol=1.0e-5) + 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=(8, 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) * np.nan + T = 1500.0 + for i, P in enumerate(pressures): + assemblage.set_state(P, T) + try: + sol = equilibrate( + assemblage.formula, assemblage, [["P", P], ["T", T]], tol=1.0e-7 + )[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.0) + except AssertionError: + pass + ax[0].plot(pressures / 1.0e9, logfO2, linestyle=":", linewidth=3.0) + + ax[0].set_ylim(-14, 22) + plt.show() + + +def check_fig_a1_fe2o3(): + print("\nChecking Figure A1...") + 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=(8, 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) + + bdg.set_state(60.0e9, 300.0) + print(f"(Mg0.5Fe0.5)(Fe0.5Si0.5)O3 V at 60 GPa, 300 K: {bdg.V*1.e6:.3f} cm^3/mol") + plt.show() + + +def check_fig_a4_stishovite(): + print("\nChecking Figure A4...") + stv = SLB_2024.stishovite() + + fig = plt.figure(figsize=(8, 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) + + stv.set_state(1.0e10, 2000.0) + print(f"Stv Vs at 10 GPa, 2000 K: {stv.shear_wave_velocity/1.e3:.3f} km/s") + plt.show() + + +if __name__ == "__main__": + check_iron_bearing_mineral_property_values() + check_bcc_iron_consistency() + check_fper_entropy() + 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/performance.py b/misc/performance.py index 4f09b083d..910f1391a 100644 --- a/misc/performance.py +++ b/misc/performance.py @@ -46,9 +46,10 @@ def test_grueneisen(): for aa in range(0, 500): a = burnman.eos.SLB3() m = minerals.SLB_2011.fayalite() + m.set_state(1.0e9, 1.0e3) x = 0 for i in range(0, 10000): - x += a.grueneisen_parameter(1e9, 1e4, 1.1, m.params) + x += m.grueneisen_parameter return x test_grueneisen() diff --git a/misc/ref/slb_2024_benchmarks.py.out b/misc/ref/slb_2024_benchmarks.py.out new file mode 100644 index 000000000..4a846b151 --- /dev/null +++ b/misc/ref/slb_2024_benchmarks.py.out @@ -0,0 +1,67 @@ +Checking iron bearing mineral properties... +wu in agreement with HeFESTo (max rel error 0.0001 %) +mag in agreement with HeFESTo (max rel error 0.0001 %) +wu in agreement with HeFESTo (max rel error 0.0001 %) +mag in agreement with HeFESTo (max rel error 0.0001 %) +hem in agreement with HeFESTo (max rel error 0.0001 %) +hepv in agreement with HeFESTo (max rel error 0.0001 %) +hlpv in agreement with HeFESTo (max rel error 0.0001 %) +hppv in agreement with HeFESTo (max rel error 0.0001 %) +fea in agreement with HeFESTo (max rel error 0.0002 %) +feg in agreement with HeFESTo (max rel error 0.0006 %) +fee in agreement with HeFESTo (max rel error 0.0003 %) +hmag in agreement with HeFESTo (max rel error 0.0001 %) + +Checking EoS consistency for 'burnman.minerals.SLB_2024.fea' +Expressions within tolerance of 0.000100 +G = F + PV : True +G = H - TS : True +G = E - TS + PV : True +S = -dG/dT : True +alpha = 1/V dV/dT : True +C_p = T dS/dT : True +V = dG/dP : True +K_T = -V dP/dV : True +C_v = Cp - alpha^2*K_T*V*T : True +K_S = K_T*Cp/Cv : True +gr = alpha*K_T*V/Cv : True +Vphi = np.sqrt(K_S/rho) : True +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' + +Checking fper configurational entropy... +Endmember entropies (should all be zero): +[0. 0. 0. 0. 0.] +Equimolar wu-mag: 11.52629 J/K/mol (should be 11.52629 J/K/mol) + +Checking Figure 1... +Mg2Fe2O4 ferropericlase at 70 GPa, 300 K: +Molar proportions: 0.50, 0.24, 0.26, 0.00, 0.00 +Volume: 35.49 cm^3/mol + +Checking Figure 3... +BCC-fper-mag triple point at 1 bar: 825.98 K + +Checking Figure 6a... +BCC iron Cp, V at 1 bar, 1000 K: 52.29 J/K/mol, 7.30 cm^3/mol + +Checking Figure 6c... +The bcc-fcc-hcp invariant is at 11.27 GPa, 825.7 K + +Checking Figure 7... +fO2 values at 1500 K: +['Fayalite', 'gamma_(fcc)_Iron', 'Quartz'] +log10fO2 = -12.076 +['Fayalite', 'alpha_(bcc)_Iron', 'Quartz'] +log10fO2 = -12.077 +['Fayalite', 'Magnetite', 'Quartz'] +log10fO2 = -8.233 +['Magnetite', 'Hematite', 'Hematite'] +log10fO2 = -2.635 + +Checking Figure A1... +(Mg0.5Fe0.5)(Fe0.5Si0.5)O3 V at 60 GPa, 300 K: 21.634 cm^3/mol + +Checking Figure A4... +Stv Vs at 10 GPa, 2000 K: 6.421 km/s diff --git a/misc/ref/table_mineral_library.py.out b/misc/ref/table_mineral_library.py.out index 7d560cbf3..c0341cd0d 100644 --- a/misc/ref/table_mineral_library.py.out +++ b/misc/ref/table_mineral_library.py.out @@ -277,7 +277,7 @@ HHPH_2013.maj 0.00011457 160000000000.0 4.56 -2.8e-11 HHPH_2013.mak 2.635e-05 211000000000.0 4.55 -2.2e-11 0.1003887 5.0 [147.8, 0.002015, -2395000.0, -801.8] HHPH_2013.manal 0.00011166 184000000000.0 4.0 -2.2e-11 0.42679680000000003 21.0 [600.0, 0.018756, -8989200.0, -2665.2] HHPH_2013.mcor 2.635e-05 211000000000.0 4.55 -2.2e-11 0.1003887 5.0 [147.8, 0.002015, -2395000.0, -801.8] -HHPH_2013.mgts 6.05e-05 102800000000.0 8.55 -8.3e-11 0.20234990000000003 10.0 [371.4, -0.004082, -398400.0, -3547.1] +HHPH_2013.mgts 6.05e-05 102800000000.0 8.55 -8.3e-11 0.2023499 10.0 [371.4, -0.004082, -398400.0, -3547.1] HHPH_2013.mpv 2.445e-05 251000000000.0 4.14 -1.6e-11 0.1003887 5.0 [149.3, 0.002918, -2983000.0, -799.1] HHPH_2013.mrw 3.949e-05 178100000000.0 4.35 -2.4e-11 0.14069310000000002 7.0 [213.3, 0.00269, -1410400.0, -1495.9] HHPH_2013.mscf 3.649e-05 185000000000.0 4.0 -1.7e-11 0.14069310000000002 7.0 [213.3, 0.00269, -1410400.0, -1495.9] @@ -824,7 +824,7 @@ SLB_2011.mg_majorite 0.000114324 165118300000.0 4 SLB_2011.mg_perovskite 2.4445e-05 250526400000.0 4.14 172900000000.0 1.69037 0.1003887 5.0 905.9412 1.56508 1.10945 2.56536 SLB_2011.mg_post_perovskite 2.4419e-05 231200000000.0 4.0 150167000000.0 1.97874 0.1003887 5.0 855.8173 1.89155 1.09081 1.16704 SLB_2011.mg_ringwoodite 3.9493e-05 184900900000.0 4.22035 123000000000.0 1.35412 0.14069310000000002 7.0 877.7094 1.10791 2.3914 2.30461 -SLB_2011.mg_tschermaks 5.914e-05 107076800000.0 7.02751 95950860000.0 1.54596 0.20234990000000003 10.0 783.8404 0.78479 3.43846 2.49099 +SLB_2011.mg_tschermaks 5.914e-05 107076800000.0 7.02751 95950860000.0 1.54596 0.2023499 10.0 783.8404 0.78479 3.43846 2.49099 SLB_2011.mg_wadsleyite 4.0515e-05 168694800000.0 4.3229 112000000000.0 1.44424 0.14069310000000002 7.0 843.4973 1.2061 2.0188 2.63683 SLB_2011.mgc2 (SLB_2011.hp_clinoenstatite) 6.076e-05 116025400000.0 6.23685 87927170000.0 1.84119 0.2007774 10.0 824.4439 1.12473 0.20401 2.14181 SLB_2011.mgcf (SLB_2011.mg_ca_ferrite) 3.6177e-05 210666300000.0 4.0528 129826000000.0 1.75878 0.1422656 7.0 838.6291 1.31156 1.0 2.1073 @@ -832,7 +832,7 @@ SLB_2011.mgil (SLB_2011.mg_akimotoite) 2.6354e-05 210706000000.0 5 SLB_2011.mgmj (SLB_2011.mg_majorite) 0.000114324 165118300000.0 4.21183 84999990000.0 1.42969 0.4015548 20.0 822.458 0.97682 1.53581 1.0178 SLB_2011.mgpv (SLB_2011.mg_perovskite) 2.4445e-05 250526400000.0 4.14 172900000000.0 1.69037 0.1003887 5.0 905.9412 1.56508 1.10945 2.56536 SLB_2011.mgri (SLB_2011.mg_ringwoodite) 3.9493e-05 184900900000.0 4.22035 123000000000.0 1.35412 0.14069310000000002 7.0 877.7094 1.10791 2.3914 2.30461 -SLB_2011.mgts (SLB_2011.mg_tschermaks) 5.914e-05 107076800000.0 7.02751 95950860000.0 1.54596 0.20234990000000003 10.0 783.8404 0.78479 3.43846 2.49099 +SLB_2011.mgts (SLB_2011.mg_tschermaks) 5.914e-05 107076800000.0 7.02751 95950860000.0 1.54596 0.2023499 10.0 783.8404 0.78479 3.43846 2.49099 SLB_2011.mgwa (SLB_2011.mg_wadsleyite) 4.0515e-05 168694800000.0 4.3229 112000000000.0 1.44424 0.14069310000000002 7.0 843.4973 1.2061 2.0188 2.63683 SLB_2011.mppv (SLB_2011.mg_post_perovskite) 2.4419e-05 231200000000.0 4.0 150167000000.0 1.97874 0.1003887 5.0 855.8173 1.89155 1.09081 1.16704 SLB_2011.na_ca_ferrite 3.627e-05 161338500000.0 4.32479 122004900000.0 2.07687 0.1420544 7.0 812.4769 0.69428 1.0 2.79016 @@ -988,6 +988,145 @@ SLB_2022.wu 4.9024e-05 160700000000.0 SLB_2022.wuls 4.33997e-05 199700000000.0 4.0 59000000000.0 1.44073 0.28737822 8.0 524.57881 1.45033 1.5487 -0.13801 SLB_2022.wustite (SLB_2022.wu) 4.9024e-05 160700000000.0 4.0 59000000000.0 1.44764 0.28737822 8.0 454.1752 1.45033 1.5487 0.06776 SLB_2022.wustite_low_spin (SLB_2022.wuls) 4.33997e-05 199700000000.0 4.0 59000000000.0 1.44073 0.28737822 8.0 524.57881 1.45033 1.5487 -0.13801 +Name (slb3 equation of state) V_0 K_0 Kprime_0 G_0 Gprime_0 molar_mass n Debye_0 grueneisen_0 q_0 eta_s_0 +SLB_2024.ab 0.000100452 59752590000.0 2.77846 36000000000.0 1.38571 0.262223 13.0 719.0831 0.57877 1.0 1.02954 +SLB_2024.acm 6.4606e-05 116100000000.0 4.4 72454000000.0 1.71384 0.231 10.0 726.94116 0.77467 0.60142 1.7391 +SLB_2024.acmite (SLB_2024.acm) 6.4606e-05 116100000000.0 4.4 72454000000.0 1.71384 0.231 10.0 726.94116 0.77467 0.60142 1.7391 +SLB_2024.al 0.00011543 173896370000.0 4.91341 96000000000.0 1.40927 0.49776 20.0 741.38227 1.06493 1.42169 2.09289 +SLB_2024.al_perovskite (SLB_2024.alpv) 2.4944e-05 242400000000.0 4.1 169199620000.0 1.55703 0.101961 5.0 856.18212 1.54466 0.83352 2.27978 +SLB_2024.al_post_perovskite (SLB_2024.appv) 2.45e-05 247740000000.0 4.0 91897150000.0 1.81823 0.101961 5.0 752.02929 1.86524 1.76454 2.70624 +SLB_2024.albite (SLB_2024.ab) 0.000100452 59752590000.0 2.77846 36000000000.0 1.38571 0.262223 13.0 719.0831 0.57877 1.0 1.02954 +SLB_2024.almandine (SLB_2024.al) 0.00011543 173896370000.0 4.91341 96000000000.0 1.40927 0.49776 20.0 741.38227 1.06493 1.42169 2.09289 +SLB_2024.alpha_nao2_phase (SLB_2024.anao) 4.542e-05 161143930000.0 3.90838 108465250000.0 2.14668 0.16394023 8.0 744.51451 1.45033 1.5487 0.77305 +SLB_2024.alpha_pbo_2_sio_2 (SLB_2024.apbo) 1.367e-05 327162360000.0 4.0166 227412210000.0 1.77076 0.060085 3.0 1132.97205 1.55723 2.21141 4.56617 +SLB_2024.alpv 2.4944e-05 242400000000.0 4.1 169199620000.0 1.55703 0.101961 5.0 856.18212 1.54466 0.83352 2.27978 +SLB_2024.an 0.00010061 84093390000.0 6.73397 39900000000.0 1.09129 0.278211 13.0 754.46887 0.38512 1.0 1.63405 +SLB_2024.anao 4.542e-05 161143930000.0 3.90838 108465250000.0 2.14668 0.16394023 8.0 744.51451 1.45033 1.5487 0.77305 +SLB_2024.andr 0.00013199 153580040000.0 4.78447 89700000000.0 0.97013 0.50817 20.0 750.98472 1.04336 1.42169 2.84503 +SLB_2024.anorthite (SLB_2024.an) 0.00010061 84093390000.0 6.73397 39900000000.0 1.09129 0.278211 13.0 754.46887 0.38512 1.0 1.63405 +SLB_2024.apbo 1.367e-05 327162360000.0 4.0166 227412210000.0 1.77076 0.060085 3.0 1132.97205 1.55723 2.21141 4.56617 +SLB_2024.appv 2.45e-05 247740000000.0 4.0 91897150000.0 1.81823 0.101961 5.0 752.02929 1.86524 1.76454 2.70624 +SLB_2024.ca_perovskite (SLB_2024.capv) 2.745e-05 236000000000.0 3.9 126000000000.0 1.6 0.116164 5.0 800.29043 1.88997 0.89608 1.35422 +SLB_2024.capv 2.745e-05 236000000000.0 3.9 126000000000.0 1.6 0.116164 5.0 800.29043 1.88997 0.89608 1.35422 +SLB_2024.cats 6.3574e-05 113759900000.0 4.8061 75337410000.0 1.71384 0.218123 10.0 804.36068 0.82288 0.60142 1.77238 +SLB_2024.cen 6.25e-05 113759900000.0 4.8061 76810310000.0 1.71384 0.2007774 10.0 805.59286 1.00921 0.60142 1.42387 +SLB_2024.clinoenstatite (SLB_2024.cen) 6.25e-05 113759900000.0 4.8061 76810310000.0 1.71384 0.2007774 10.0 805.59286 1.00921 0.60142 1.42387 +SLB_2024.co 2.5577e-05 252585720000.0 3.88671 163200000000.0 1.64704 0.101961 5.0 932.21586 1.3081 1.71245 2.84761 +SLB_2024.coes 2.0657e-05 103537990000.0 2.9007 61600000000.0 0.49686 0.060085 3.0 875.22323 0.29043 1.0 2.75631 +SLB_2024.coesite (SLB_2024.coes) 2.0657e-05 103537990000.0 2.9007 61600000000.0 0.49686 0.060085 3.0 875.22323 0.29043 1.0 2.75631 +SLB_2024.corundum (SLB_2024.co) 2.5577e-05 252585720000.0 3.88671 163200000000.0 1.64704 0.101961 5.0 932.21586 1.3081 1.71245 2.84761 +SLB_2024.cppv 2.6949e-05 247740000000.0 4.0 187873500000.0 1.98845 0.15199 5.0 755.01863 1.64015 1.76454 3.16235 +SLB_2024.cr_ca_ferrite (SLB_2024.crcf) 3.9421e-05 185400000000.0 4.0 141453800000.0 1.93591 0.19229 7.0 684.9543 1.56656 1.0 1.93943 +SLB_2024.cr_perovskite (SLB_2024.crpv) 2.8189e-05 250565350000.0 4.13438 152053560000.0 1.73254 0.15199 5.0 758.1187 1.54466 0.83352 2.83537 +SLB_2024.cr_post_perovskite (SLB_2024.cppv) 2.6949e-05 247740000000.0 4.0 187873500000.0 1.98845 0.15199 5.0 755.01863 1.64015 1.76454 3.16235 +SLB_2024.crcf 3.9421e-05 185400000000.0 4.0 141453800000.0 1.93591 0.19229 7.0 684.9543 1.56656 1.0 1.93943 +SLB_2024.crpv 2.8189e-05 250565350000.0 4.13438 152053560000.0 1.73254 0.15199 5.0 758.1187 1.54466 0.83352 2.83537 +SLB_2024.di 6.6039e-05 113759900000.0 4.8061 72700000000.0 1.71384 0.2165504 10.0 782.57306 1.00921 0.60142 1.06175 +SLB_2024.diopside (SLB_2024.di) 6.6039e-05 113759900000.0 4.8061 72700000000.0 1.71384 0.2165504 10.0 782.57306 1.00921 0.60142 1.06175 +SLB_2024.en 6.2676e-05 107076810000.0 7.02751 76800000000.0 1.54596 0.2007774 10.0 812.21227 0.78477 3.43847 2.5045 +SLB_2024.enstatite (SLB_2024.en) 6.2676e-05 107076810000.0 7.02751 76800000000.0 1.54596 0.2007774 10.0 812.21227 0.78477 3.43847 2.5045 +SLB_2024.esk 2.8904e-05 233340360000.0 4.01705 123200000000.0 1.81492 0.15199 5.0 766.73627 1.15191 2.22481 2.55521 +SLB_2024.eskolaite (SLB_2024.esk) 2.8904e-05 233340360000.0 4.01705 123200000000.0 1.81492 0.15199 5.0 766.73627 1.15191 2.22481 2.55521 +SLB_2024.fa 4.629e-05 136485580000.0 4.88157 51220000000.0 0.85893 0.203777 7.0 618.96116 1.08388 2.88055 1.65937 +SLB_2024.fapv 2.709e-05 223325500000.0 4.1 159881770000.0 1.73254 0.1308249 5.0 755.62079 1.54466 0.83352 2.8548 +SLB_2024.fayalite (SLB_2024.fa) 4.629e-05 136485580000.0 4.88157 51220000000.0 0.85893 0.203777 7.0 618.96116 1.08388 2.88055 1.65937 +SLB_2024.fe2o3_perovskite_hs (SLB_2024.hepv) 2.95768e-05 204251000000.0 4.1 123483250000.0 1.73254 0.15968852 5.0 646.79863 1.54466 0.83352 1.88876 +SLB_2024.fe2o3_perovskite_ls (SLB_2024.hlpv) 2.75209e-05 204251000000.0 4.1 177577950000.0 1.73254 0.15968852 5.0 759.63863 1.54466 0.83352 3.54218 +SLB_2024.fe_akimotoite (SLB_2024.feil) 2.6854e-05 210692630000.0 5.2154 167069220000.0 1.81492 0.131931 5.0 760.91558 1.19328 2.22481 3.6318 +SLB_2024.fe_ca_ferrite (SLB_2024.fecf) 3.7216e-05 213000000000.0 4.1 164115350000.0 1.93591 0.173806 7.0 715.88779 1.56656 1.0 2.46745 +SLB_2024.fe_nal_phase (SLB_2024.fnal) 0.000112045 204006720000.0 4.31789 152250130000.0 1.74226 0.48966601 21.0 788.03574 1.43147 1.0 2.73801 +SLB_2024.fe_perovskite (SLB_2024.fepv) 2.5321e-05 270582550000.00003 4.01 130058789999.99998 1.37279 0.131931 5.0 740.39231 1.54466 0.83352 2.08541 +SLB_2024.fe_postperovskite (SLB_2024.fppv) 2.52963e-05 247740000000.0 4.0 129500000000.0 1.40076 0.131931 5.0 769.31113 1.64015 1.76454 1.8495 +SLB_2024.fe_ringwoodite (SLB_2024.feri) 4.186e-05 213409100000.0 4.22035 92000000000.0 1.35412 0.203777 7.0 651.49411 1.26156 2.39214 1.76941 +SLB_2024.fe_wadsleyite (SLB_2024.fewa) 4.28e-05 168566850000.0 4.12302 72000000000.0 1.50973 0.203777 7.0 636.8306 1.20498 2.20831 0.94487 +SLB_2024.fealo3_perovskite_hs (SLB_2024.fapv) 2.709e-05 223325500000.0 4.1 159881770000.0 1.73254 0.1308249 5.0 755.62079 1.54466 0.83352 2.8548 +SLB_2024.fec2 6.38541e-05 116025560000.0 6.23685 79289860000.0 1.84119 0.2638614 10.0 691.84626 1.12478 0.20409 1.14272 +SLB_2024.fecf 3.7216e-05 213000000000.0 4.1 164115350000.0 1.93591 0.173806 7.0 715.88779 1.56656 1.0 2.46745 +SLB_2024.feil 2.6854e-05 210692630000.0 5.2154 167069220000.0 1.81492 0.131931 5.0 760.91558 1.19328 2.22481 3.6318 +SLB_2024.fepv 2.5321e-05 270582550000.00003 4.01 130058789999.99998 1.37279 0.131931 5.0 740.39231 1.54466 0.83352 2.08541 +SLB_2024.feri 4.186e-05 213409100000.0 4.22035 92000000000.0 1.35412 0.203777 7.0 651.49411 1.26156 2.39214 1.76941 +SLB_2024.ferrosilite (SLB_2024.fs) 6.5941e-05 100545440000.0 7.87556 52000000000.0 1.54596 0.2638614 10.0 677.91886 0.7144 3.43847 1.08228 +SLB_2024.fewa 4.28e-05 168566850000.0 4.12302 72000000000.0 1.50973 0.203777 7.0 636.8306 1.20498 2.20831 0.94487 +SLB_2024.fnal 0.000112045 204006720000.0 4.31789 152250130000.0 1.74226 0.48966601 21.0 788.03574 1.43147 1.0 2.73801 +SLB_2024.fo 4.3603e-05 127955500000.0 4.21796 81600000000.0 1.46257 0.140695 7.0 809.1977 0.9928 2.10671 2.29968 +SLB_2024.forsterite (SLB_2024.fo) 4.3603e-05 127955500000.0 4.21796 81600000000.0 1.46257 0.140695 7.0 809.1977 0.9928 2.10671 2.29968 +SLB_2024.fppv 2.52963e-05 247740000000.0 4.0 129500000000.0 1.40076 0.131931 5.0 769.31113 1.64015 1.76454 1.8495 +SLB_2024.fs 6.5941e-05 100545440000.0 7.87556 52000000000.0 1.54596 0.2638614 10.0 677.91886 0.7144 3.43847 1.08228 +SLB_2024.gr 0.00012512 167062260000.0 3.91544 109000000000.0 1.16274 0.450449 20.0 822.77062 1.05402 1.88886 2.38415 +SLB_2024.grossular (SLB_2024.gr) 0.00012512 167062260000.0 3.91544 109000000000.0 1.16274 0.450449 20.0 822.77062 1.05402 1.88886 2.38415 +SLB_2024.hc 4.0843e-05 208947270000.0 4.62702 84500000000.0 0.62795 0.17381 7.0 747.13664 1.18794 3.97087 2.46339 +SLB_2024.he 6.7867e-05 119204720000.0 4.81927 61000000000.0 1.71384 0.2480924 10.0 702.08234 0.96665 0.60142 1.01745 +SLB_2024.hedenbergite (SLB_2024.he) 6.7867e-05 119204720000.0 4.81927 61000000000.0 1.71384 0.2480924 10.0 702.08234 0.96665 0.60142 1.01745 +SLB_2024.hem 3.0287e-05 204245370000.0 4.0997 91000000000.0 1.81492 0.15968852 5.0 653.80747 1.58944 2.22481 0.5241 +SLB_2024.hematite (SLB_2024.hem) 3.0287e-05 204245370000.0 4.0997 91000000000.0 1.81492 0.15968852 5.0 653.80747 1.58944 2.22481 0.5241 +SLB_2024.hepv 2.95768e-05 204251000000.0 4.1 123483250000.0 1.73254 0.15968852 5.0 646.79863 1.54466 0.83352 1.88876 +SLB_2024.hercynite (SLB_2024.hc) 4.0843e-05 208947270000.0 4.62702 84500000000.0 0.62795 0.17381 7.0 747.13664 1.18794 3.97087 2.46339 +SLB_2024.high_pressure_magnetit (SLB_2024.hmag) 4.1702e-05 172000000000.0 4.0 120889150000.0 1.93591 0.23153307 7.0 542.9312 1.56656 0.41872 1.37608 +SLB_2024.hlpv 2.75209e-05 204251000000.0 4.1 177577950000.0 1.73254 0.15968852 5.0 759.63863 1.54466 0.83352 3.54218 +SLB_2024.hmag 4.1702e-05 172000000000.0 4.0 120889150000.0 1.93591 0.23153307 7.0 542.9312 1.56656 0.41872 1.37608 +SLB_2024.hp_clinoenstatite (SLB_2024.mgc2) 6.076e-05 116025560000.0 6.23685 87927170000.0 1.84119 0.2007774 10.0 824.44051 1.12478 0.20409 2.14193 +SLB_2024.hp_clinoferrosilite (SLB_2024.fec2) 6.38541e-05 116025560000.0 6.23685 79289860000.0 1.84119 0.2638614 10.0 691.84626 1.12478 0.20409 1.14272 +SLB_2024.hppv 2.76884e-05 176500000000.0 4.0 172363480000.0 1.98845 0.15968852 5.0 680.92363 1.64015 1.76454 2.56327 +SLB_2024.hs_fe2o3_post_perovski (SLB_2024.hppv) 2.76884e-05 176500000000.0 4.0 172363480000.0 1.98845 0.15968852 5.0 680.92363 1.64015 1.76454 2.56327 +SLB_2024.jadeite (SLB_2024.jd) 6.0508e-05 137357290000.0 3.6305 84000000000.0 1.71384 0.2021387 10.0 820.2985 0.85734 2.05453 1.91023 +SLB_2024.jd 6.0508e-05 137357290000.0 3.6305 84000000000.0 1.71384 0.2021387 10.0 820.2985 0.85734 2.05453 1.91023 +SLB_2024.knor 0.000118716 157000000000.0 4.5 99527660000.0 1.35756 0.45316 20.0 776.39637 1.24672 1.42169 2.11433 +SLB_2024.knorringite (SLB_2024.knor) 0.000118716 157000000000.0 4.5 99527660000.0 1.35756 0.45316 20.0 776.39637 1.24672 1.42169 2.11433 +SLB_2024.ky 4.4227e-05 160000000000.0 4.0 117646620000.0 1.69117 0.1620456 8.0 943.19593 0.92549 1.0 2.89863 +SLB_2024.kyanite (SLB_2024.ky) 4.4227e-05 160000000000.0 4.0 117646620000.0 1.69117 0.1620456 8.0 943.19593 0.92549 1.0 2.89863 +SLB_2024.lime_tschermak (SLB_2024.cats) 6.3574e-05 113759900000.0 4.8061 75337410000.0 1.71384 0.218123 10.0 804.36068 0.82288 0.60142 1.77238 +SLB_2024.lppv 2.62534e-05 176500000000.0 4.0 220101260000.0 1.98843 0.15968852 5.0 713.13099 1.64029 1.76443 3.79074 +SLB_2024.ls_fe2o3_post_perovski (SLB_2024.lppv) 2.62534e-05 176500000000.0 4.0 220101260000.0 1.98843 0.15968852 5.0 713.13099 1.64029 1.76443 3.79074 +SLB_2024.mag 4.4528e-05 183876640000.0 5.22819 60300000000.0 0.04657 0.23153307 7.0 529.46966 1.35821 0.41872 1.12258 +SLB_2024.magnetite (SLB_2024.smag) 4.4528e-05 183876640000.0 5.22819 60300000000.0 0.04657 0.23153307 7.0 529.46966 1.35821 0.41872 1.12258 +SLB_2024.mg_akimotoite (SLB_2024.mgil) 2.6354e-05 210692630000.0 5.2154 132000000000.0 1.81492 0.100389 5.0 928.95623 1.19328 2.22481 3.3993 +SLB_2024.mg_ca_ferrite (SLB_2024.mgcf) 3.6135e-05 213000000000.0 4.1 129699999999.99998 1.93591 0.142266 7.0 830.714 1.56656 1.0 1.30292 +SLB_2024.mg_majorite (SLB_2024.mgmj) 0.000114324 165118370000.0 4.21183 85000000000.0 1.42969 0.40156 20.0 822.48562 0.97681 1.53581 1.01779 +SLB_2024.mg_nal_phase (SLB_2024.mnal) 0.000109883 204006720000.0 4.31789 129000000000.0 1.74226 0.42658581 21.0 868.59088 1.43147 1.0 1.9384 +SLB_2024.mg_perovskite (SLB_2024.mgpv) 2.4445e-05 250565350000.0 4.13438 172900000000.0 1.73254 0.100389 5.0 892.95164 1.54466 0.83352 1.65233 +SLB_2024.mg_postperovskite (SLB_2024.mppv) 2.41108e-05 247740000000.0 4.0 166644930000.0 1.98845 0.100389 5.0 931.02549 1.64015 1.76454 1.43815 +SLB_2024.mg_ringwoodite (SLB_2024.mgri) 3.9493e-05 184901750000.0 4.22035 123000000000.0 1.35412 0.140693 7.0 879.84656 1.10843 2.39214 2.30588 +SLB_2024.mg_tschermaks (SLB_2024.mgts) 5.914e-05 107076810000.0 7.02751 93345240000.0 1.54596 0.20235 10.0 788.01368 0.78477 3.43847 2.39272 +SLB_2024.mg_wadsleyite (SLB_2024.mgwa) 4.0515e-05 168702760000.0 4.12302 112000000000.0 1.50973 0.140695 7.0 849.12535 1.20498 2.20831 2.56411 +SLB_2024.mgc2 6.076e-05 116025560000.0 6.23685 87927170000.0 1.84119 0.2007774 10.0 824.44051 1.12478 0.20409 2.14193 +SLB_2024.mgcf 3.6135e-05 213000000000.0 4.1 129699999999.99998 1.93591 0.142266 7.0 830.714 1.56656 1.0 1.30292 +SLB_2024.mgil 2.6354e-05 210692630000.0 5.2154 132000000000.0 1.81492 0.100389 5.0 928.95623 1.19328 2.22481 3.3993 +SLB_2024.mgmj 0.000114324 165118370000.0 4.21183 85000000000.0 1.42969 0.40156 20.0 822.48562 0.97681 1.53581 1.01779 +SLB_2024.mgpv 2.4445e-05 250565350000.0 4.13438 172900000000.0 1.73254 0.100389 5.0 892.95164 1.54466 0.83352 1.65233 +SLB_2024.mgri 3.9493e-05 184901750000.0 4.22035 123000000000.0 1.35412 0.140693 7.0 879.84656 1.10843 2.39214 2.30588 +SLB_2024.mgts 5.914e-05 107076810000.0 7.02751 93345240000.0 1.54596 0.20235 10.0 788.01368 0.78477 3.43847 2.39272 +SLB_2024.mgwa 4.0515e-05 168702760000.0 4.12302 112000000000.0 1.50973 0.140695 7.0 849.12535 1.20498 2.20831 2.56411 +SLB_2024.mnal 0.000109883 204006720000.0 4.31789 129000000000.0 1.74226 0.42658581 21.0 868.59088 1.43147 1.0 1.9384 +SLB_2024.mppv 2.41108e-05 247740000000.0 4.0 166644930000.0 1.98845 0.100389 5.0 931.02549 1.64015 1.76454 1.43815 +SLB_2024.na_ca_ferrite (SLB_2024.nacf) 3.627e-05 220000000000.0 4.1 135023869999.99998 1.74135 0.142054 7.0 709.33152 1.56656 1.0 1.67725 +SLB_2024.na_majorite (SLB_2024.namj) 0.000110842 172035520000.0 5.20045 114700000000.0 1.35756 0.40270437 20.0 845.23671 1.25087 0.10909 2.48526 +SLB_2024.na_nal_phase (SLB_2024.nnal) 0.000109401 204006720000.0 4.31789 144164270000.0 1.74226 0.42616294 21.0 846.08425 1.43147 1.0 2.40665 +SLB_2024.nacf 3.627e-05 220000000000.0 4.1 135023869999.99998 1.74135 0.142054 7.0 709.33152 1.56656 1.0 1.67725 +SLB_2024.namj 0.000110842 172035520000.0 5.20045 114700000000.0 1.35756 0.40270437 20.0 845.23671 1.25087 0.10909 2.48526 +SLB_2024.neph 5.38684e-05 53055500000.0 4.0 30700000000.0 1.33087 0.14205431 7.0 743.57985 0.6969 1.0 0.6241 +SLB_2024.nepheline (SLB_2024.neph) 5.38684e-05 53055500000.0 4.0 30700000000.0 1.33087 0.14205431 7.0 743.57985 0.6969 1.0 0.6241 +SLB_2024.nnal 0.000109401 204006720000.0 4.31789 144164270000.0 1.74226 0.42616294 21.0 846.08425 1.43147 1.0 2.40665 +SLB_2024.odi 6.8054e-05 107076810000.0 7.02751 58247330000.0 1.54596 0.2165504 10.0 744.48915 0.78477 3.43847 1.35202 +SLB_2024.ortho_diopside (SLB_2024.odi) 6.8054e-05 107076810000.0 7.02751 58247330000.0 1.54596 0.2165504 10.0 744.48915 0.78477 3.43847 1.35202 +SLB_2024.pe 4.4976e-05 161143930000.0 3.90838 130900000000.0 2.14668 0.16121782 8.0 770.90151 1.45033 1.5487 2.56123 +SLB_2024.periclase (SLB_2024.pe) 4.4976e-05 161143930000.0 3.90838 130900000000.0 2.14668 0.16121782 8.0 770.90151 1.45033 1.5487 2.56123 +SLB_2024.picr 4.3564e-05 184400000000.0 5.7 90290490000.0 0.62795 0.19229 7.0 750.72523 0.99168 3.97087 3.07014 +SLB_2024.pircochromite (SLB_2024.picr) 4.3564e-05 184400000000.0 5.7 90290490000.0 0.62795 0.19229 7.0 750.72523 0.99168 3.97087 3.07014 +SLB_2024.pseudo_wollastonite (SLB_2024.pwo) 4.0272e-05 86000000000.0 3.8 28991610000.0 0.77536 0.116164 5.0 703.00475 0.95232 1.0 0.76187 +SLB_2024.pwo 4.0272e-05 86000000000.0 3.8 28991610000.0 0.77536 0.116164 5.0 703.00475 0.95232 1.0 0.76187 +SLB_2024.py 0.00011308 170239640000.0 4.11067 93700000000.0 1.35756 0.40313 20.0 823.23783 1.01422 1.42169 0.98186 +SLB_2024.pyrope (SLB_2024.py) 0.00011308 170239640000.0 4.11067 93700000000.0 1.35756 0.40313 20.0 823.23783 1.01422 1.42169 0.98186 +SLB_2024.qtz 2.24211e-05 61425960000.0 19.78014 44857750000.0 -0.04277 0.060085 3.0 883.46813 -0.03957 1.0 2.40464 +SLB_2024.quartz (SLB_2024.qtz) 2.24211e-05 61425960000.0 19.78014 44857750000.0 -0.04277 0.060085 3.0 883.46813 -0.03957 1.0 2.40464 +SLB_2024.smag 4.4528e-05 183876640000.0 5.22819 60300000000.0 0.04657 0.23153307 7.0 529.46966 1.35821 0.41872 1.12258 +SLB_2024.sp 3.9762e-05 195103490000.0 4.62702 109000000000.0 0.62795 0.142278 7.0 801.86054 0.97405 3.97087 2.4035 +SLB_2024.spinel (SLB_2024.sp) 3.9762e-05 195103490000.0 4.62702 109000000000.0 0.62795 0.142278 7.0 801.86054 0.97405 3.97087 2.4035 +SLB_2024.wo 3.9901e-05 93900000000.0 4.0 50300000000.0 1.23206 0.116164 5.0 713.58788 1.05734 1.0 1.2711 +SLB_2024.wollastonite (SLB_2024.wo) 3.9901e-05 93900000000.0 4.0 50300000000.0 1.23206 0.116164 5.0 713.58788 1.05734 1.0 1.2711 +SLB_2024.wu 4.9024e-05 160700000000.0 4.0 59000000000.0 1.44764 0.28737822 8.0 427.00102 1.45033 1.5487 0.0596 +SLB_2024.wuls 4.33997e-05 199700000000.0 4.0 59000000000.0 1.44073 0.28737822 8.0 492.99392 1.45033 1.5487 -0.14773 +SLB_2024.wustite (SLB_2024.wu) 4.9024e-05 160700000000.0 4.0 59000000000.0 1.44764 0.28737822 8.0 427.00102 1.45033 1.5487 0.0596 +SLB_2024.wustite_low_spin (SLB_2024.wuls) 4.33997e-05 199700000000.0 4.0 59000000000.0 1.44073 0.28737822 8.0 492.99392 1.45033 1.5487 -0.14773 Name (hp_tmt equation of state) V_0 K_0 Kprime_0 Kdprime_0 molar_mass n Cp Sundman_1991.bcc_iron 7.09e-06 164000000000.0 5.16 -3.1e-11 0.055845 1.0 [21.09, 0.0101455, -221508.0, 47.1947] Sundman_1991.fcc_iron 6.93863394593e-06 153900000000.0 5.2 -3.37e-11 0.055845 1.0 [22.24, 0.0088656, -221517.0, 47.1998] diff --git a/tests/test_minerals.py b/tests/test_minerals.py index 9a1face29..7a1d4d413 100644 --- a/tests/test_minerals.py +++ b/tests/test_minerals.py @@ -27,6 +27,7 @@ def test_instantiate_minerals(self): for m in member_list if inspect.isclass(m) and issubclass(m, burnman.Mineral) + and not issubclass(m, burnman.RelaxedSolution) and m is not burnman.Mineral and m is not burnman.SolidSolution and m is not burnman.CombinedMineral diff --git a/tests/test_solidsolution.py b/tests/test_solidsolution.py index 223214f95..9476236bd 100644 --- a/tests/test_solidsolution.py +++ b/tests/test_solidsolution.py @@ -810,6 +810,36 @@ def test_subregular_model_ternary_hessian_multicomponent_change(self): self.assertArraysAlmostEqual(H0.dot(df), dGdx2 - dGdx1) + def test_slb_ferropericlase_partial_gibbs_multicomponent_change(self): + ss = burnman.minerals.SLB_2024.ferropericlase() + f0 = [0.15, 0.25, 0.35, 0.05, 0.2] + ss.set_composition(f0) + ss.set_state(1.0e5, 300.0) + mu0 = ss.partial_gibbs + + df = np.array([2.0e-6, -1.5e-6, -0.5e-6, 0.25e-6, -0.25e-6]) + ss.set_composition(f0 - df / 2.0) + G1 = ss.gibbs + ss.set_composition(f0 + df / 2.0) + G2 = ss.gibbs + + self.assertAlmostEqual(mu0.dot(df), G2 - G1) + + def test_slb_ferropericlase_hessian_multicomponent_change(self): + ss = burnman.minerals.SLB_2024.ferropericlase() + f0 = [0.15, 0.25, 0.35, 0.05, 0.2] + ss.set_composition(f0) + ss.set_state(1.0e5, 300.0) + H0 = ss.gibbs_hessian + + df = np.array([2.0e-6, -1.5e-6, -0.5e-6, 0.25e-6, -0.25e-6]) + ss.set_composition(f0 - df / 2.0) + dGdx1 = ss.partial_gibbs + ss.set_composition(f0 + df / 2.0) + dGdx2 = ss.partial_gibbs + + self.assertArraysAlmostEqual(H0.dot(df), dGdx2 - dGdx1) + def test_polynomial_model_partial_entropy_multicomponent_change(self): ss = two_site_ss_polynomial_ternary() f0 = np.array([0.25, 0.35, 0.4])