diff --git a/burnman/data/input_raw_endmember_datasets/hpx_eos_to_burnman.py b/burnman/data/input_raw_endmember_datasets/hpx_eos_to_burnman.py new file mode 100644 index 00000000..d4aa46e2 --- /dev/null +++ b/burnman/data/input_raw_endmember_datasets/hpx_eos_to_burnman.py @@ -0,0 +1,418 @@ +# 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 HPx-eos data format into the standard burnman +# format (printed to stdout). For example, the dataset file +# tc-thermoinput-igneous-2022-01-23/tc-ig50NCKFMASHTOCr.txt +# can be processed via the commands +# python hpx_eos_to_burnman.py > ../../minerals/ig50NCKFMASHTOCr.py; black ../../minerals/ig50NCKFMASHTOCr.py +# the text "from . import ig50NCKFMASHTOCr must then be added to the minerals __init__.py + +import numpy as np +import logging +import sys +from sympy import Symbol, prod, sympify +from sympy.parsing.sympy_parser import parse_expr +from burnman.constants import gas_constant +from burnman.minerals import HGP_2018_ds633 +from datetime import date + +# Change the following lines to process a different dataset +dataset = HGP_2018_ds633 +mbr_dataset = "HGP_2018_ds633" +solution_file = "tc-thermoinput-igneous-2022-01-23/tc-ig50NCKFMASHTOCr.txt" +solution_file = "tc-thermoinput-igneous-2022-01-23/tc-ig50NCKFMASTOCr.txt" +first_model = "pl4tr" # melt and fl are (currently unreadable) Temkin models + +# logging.basicConfig(stream=sys.stderr, level=logging.DEBUG) +logging.basicConfig(stream=sys.stderr, level=logging.ERROR) + + +current_year = str(date.today())[:4] + + +out_ss = "" +out_make = "" + + +def reverse_polish(lines): + # get symbols + syms = [] + for line in lines: + n_terms = int(line[0]) + i = 2 + for term in range(n_terms): + n_add = int(line[i]) + syms.extend(line[i + 2 : i + 2 * n_add + 2][::2]) + i += 2 * n_add + 2 + unique_syms = list(np.unique(sorted(syms))) + sympy_syms = {u: Symbol(u) for u in unique_syms} + + # get expressions + expr = 0.0 + for line in lines: + n_terms = int(line[0]) + i = 2 + terms = [] + for term in range(n_terms): + n_add = int(line[i]) + en = parse_expr(line[i - 1]) # this is the constant + for j in range(n_add): + en += parse_expr(line[i + 1 + 2 * j]) * sympy_syms[line[i + 2 + 2 * j]] + terms.append(en) + # syms.extend(l[i+2:i+2*n_add + 2][::2]) + # print("a", n_add) + i += 2 * n_add + 2 + product = prod(terms) + expr += product + return expr + + +def ordering_modifier(property_modifier, ordered): + mtype, params = property_modifier + if mtype == "bragg_williams": + if ordered: + return np.array([0.0, 0.0, 0.0]) + else: + n = params["n"] + if params["factor"] > 0.0: + f = [params["factor"], params["factor"]] + else: + f = [1.0, -params["factor"]] + + S_disord = ( + -gas_constant + * ( + f[0] * (np.log(1.0 / (n + 1.0)) + n * np.log(n / (n + 1.0))) + + f[1] + * (n * np.log(1.0 / (n + 1.0)) + n * n * np.log(n / (n + 1.0))) + ) + / (n + 1.0) + ) + return np.array([params["deltaH"], S_disord, params["deltaV"]]) + else: + if params["T_0"] < params["Tc_0"]: + Q_0 = np.power((params["Tc_0"] - params["T_0"]) / params["Tc_0"], 0.25) + else: + Q_0 = 0.0 + + if ordered: + E_ord = ( + params["Tc_0"] * params["S_D"] * (Q_0 * Q_0 - np.power(Q_0, 6.0) / 3.0) + - ( + params["S_D"] * params["Tc_0"] * 2.0 / 3.0 + - params["V_D"] * params["P_0"] + ) + - params["P_0"] * params["V_D"] * Q_0 * Q_0 + ) + S_ord = params["S_D"] * (Q_0 * Q_0 - 1.0) + V_ord = params["V_D"] * (Q_0 * Q_0 - 1.0) + return np.array([E_ord, S_ord, V_ord]) + + else: + E_disord = ( + params["Tc_0"] * params["S_D"] * (Q_0 * Q_0 - np.power(Q_0, 6.0) / 3.0) + - params["P_0"] * params["V_D"] * Q_0 * Q_0 + ) + S_disord = params["S_D"] * Q_0 * Q_0 + V_disord = params["V_D"] * Q_0 * Q_0 + return np.array([E_disord, S_disord, V_disord]) + + +increment_i = True +i = 0 +data = [] +with open(solution_file, "r") as file: + # Read each line in the file one by one + process = True + for line in file: + # Process the line (for example, print it) + ln = line.strip().split() + if len(ln) > 0: + if ln[0] == "verbatim": + process = not process + if ( + ln[0][0] != "%" + and ln[0][0] != "*" + and ln[0] != "header" + and ln[0] != "verbatim" + and process + ): + data.append(ln) + if ln[0] == first_model: + increment_i = False + if increment_i: + i += 1 + +n_data = len(data) + +noods = [] +while True: + if i > n_data - 3: + break + name = data[i][0] + n_mbrs = int(data[i][1]) + logging.debug(f"Solution: {name}") + logging.debug(f"Number of endmembers: {n_mbrs}") + + vars = data[i + 1 : i + n_mbrs] + logging.debug(vars) + i += n_mbrs + + mbr_proportions = [] + mbr_names = [] + for j in range(n_mbrs): + n_lines = int(data[i][1]) + mbr_names.append(data[i][0].split("(")[1][:-1]) + mbr_proportions.append(data[i : i + n_lines]) + i += n_lines + logging.debug(mbr_proportions) + + formulation = data[i][0] + logging.debug(f"formulation: {formulation}") + i += 1 + n_ints = int((n_mbrs * (n_mbrs - 1.0)) / 2.0) + interactions = data[i : i + n_ints] + logging.debug(interactions) + + nonideal_entropies = False + nonideal_volumes = False + + We = [[0.0 for j in mbr_names[:-i]] for i in range(1, n_mbrs)] + Ws = [[0.0 for j in mbr_names[:-i]] for i in range(1, n_mbrs)] + Wv = [[0.0 for j in mbr_names[:-i]] for i in range(1, n_mbrs)] + for ints in interactions: + m = ints[0].split("(")[1].replace(")", "").split(",") + m0 = mbr_names.index(m[0]) + m1 = mbr_names.index(m[1]) + We[m0][m1 - m0 - 1] = float(ints[1]) * 1.0e3 + Ws[m0][m1 - m0 - 1] = -float(ints[2]) + Wv[m0][m1 - m0 - 1] = float(ints[3]) * 1.0e-5 + + if np.abs(Ws[m0][m1 - m0 - 1]) > 1.0e-10: + nonideal_entropies = True + if np.abs(Wv[m0][m1 - m0 - 1]) > 1.0e-10: + nonideal_volumes = True + + i += n_ints + + if formulation == "asf": + text_alphas = data[i : i + n_mbrs] + logging.debug(text_alphas) + i += n_mbrs + + Ae = [0.0 for i in range(n_mbrs)] + As = [0.0 for i in range(n_mbrs)] + Av = [0.0 for i in range(n_mbrs)] + for alpha in text_alphas: + try: + m = alpha[0].split("(")[1].replace(")", "").split(",")[0] + except IndexError: + m = alpha[0] + m0 = mbr_names.index(m) + Ae[m0] = float(alpha[1]) + As[m0] = -float(alpha[2]) + Av[m0] = float(alpha[3]) * 1.0e6 + + n_site_species = int(data[i][0]) + + logging.debug(f"Number of site species: {n_site_species}") + i += 1 + + site_fractions = [] + for j in range(n_site_species): + n_lines = int(data[i][1]) + site_fractions.append(data[i : i + n_lines]) + i += n_lines + + summed = 0.0 + sites = [[]] + for site_species_and_expression in site_fractions: + site_species = site_species_and_expression[0][0] + expression = site_species_and_expression + expression[0] = expression[0][2:] + expression = reverse_polish(expression) + summed += expression + summed = sympify(summed) + sites[-1].append(site_species) + if summed == 1: + summed = 0 + sites.append([]) + sites = sites[:-1] + site_species = {} + for n_site, site in enumerate(sites): + for species in site: + site_species[species] = n_site + n_sites = len(sites) + + endmember_site_fractions_and_checks = [] + for j in range(n_mbrs): + n_lines = 1 + while True: + if ( + data[i + n_lines][0] == "make" + or data[i + n_lines][0] == "delG(tran)" + or data[i + n_lines][0] == "check" + or data[i + n_lines][0] == "delG(od)" + or data[i + n_lines][0] == "delG(make)" + or data[i + n_lines][0] == "delG(mod)" + ): + n_lines += 1 + else: + break + endmember_site_fractions_and_checks.append(data[i : i + n_lines]) + i += n_lines + + mbr_initializations = [] + for f in endmember_site_fractions_and_checks: + mbr = f[0][0] + if len(f) == 1: + mbr = f"{mbr_dataset}.{mbr}()" + elif len(f) == 2: + mbr = f"{mbr_dataset}.{mbr}()" + assert f[1][0] == "check" + else: + make_idx = [g[0] for g in f].index("make") + make = f[make_idx][1:] + delG = [np.array(g[1:4], dtype=float) for g in f if "delG" in g[0]] + delG = np.array([0.0, 0.0, 0.0]) if len(delG) == 0 else delG[0] + delG *= [1.0e3, -1.0e3, 1.0e-5] + n_make = int(make[0]) + el = 1 + combined_mbrs = "[" + combined_amounts = "[" + for j in range(n_make): + if make[el] == "ordered": + el += 1 + m_string = make[el] + m_amount = float(parse_expr(make[el + 1])) + mineral = getattr(dataset, m_string)() + property_modifier = mineral.property_modifiers[0] + delG += ordering_modifier(property_modifier, True) * m_amount + combined_mbrs += f"{make[el]}_nood, " + noods.append(make[el]) + elif make[el] == "disordered": + el += 1 + m_string = make[el] + m_amount = float(parse_expr(make[el + 1])) + mineral = getattr(dataset, m_string)() + property_modifier = mineral.property_modifiers[0] + delG += ordering_modifier(property_modifier, False) * m_amount + combined_mbrs += f"{make[el]}_nood, " + noods.append(make[el]) + elif make[el] == "equilibrium": + el += 1 + combined_mbrs += f"{mbr_dataset}.{make[el]}(), " + else: + combined_mbrs += f"{mbr_dataset}.{make[el]}(), " + combined_amounts += f"{float(parse_expr(make[el+1]))}, " + el += 2 + combined_mbrs = combined_mbrs[:-2] + combined_amounts = combined_amounts[:-2] + combined_mbrs += "]" + combined_amounts += "]" + out_make += f'{mbr} = CombinedMineral({combined_mbrs}, {combined_amounts}, {list(delG)}, "{mbr}")\n' + n_occs = int(f[0][2]) + site_occ = [[] for ln in range(n_sites)] + site_occ_list = f[0][3:] + ss = [[] for idx in range(n_sites)] + site_amounts = [[] for idx in range(n_sites)] + for j in range(n_occs): + site_sp = site_occ_list[2 * j] + i_site = site_species[site_sp] + # rename for insertion into string + if site_sp[0] == "x": + site_sp = site_sp[1:] + site_sp = "".join(filter(str.isalpha, site_sp)).title() + + ss[i_site].append(site_sp) + site_amounts[i_site].append(parse_expr(site_occ_list[2 * j + 1])) + total_site_amounts = [sum(amounts) for amounts in site_amounts] + site_fractions = [ + [amount / total_site_amounts[j] for amount in amounts] + for j, amounts in enumerate(site_amounts) + ] + + site_formula = "" + for m, site in enumerate(ss): + site_formula += "[" + for n, species in enumerate(site): + site_formula += species + if site_fractions[m][n] != 1: + site_formula += str(site_fractions[m][n]) + site_formula += "]" + if str(total_site_amounts[m]) != "1": + site_formula += str(total_site_amounts[m]) + + mbr_initializations.append(f'[{mbr}, "{site_formula}"]') + + out_ss += f"class {name}(Solution):\n" + out_ss += " def __init__(self, molar_fractions=None):\n" + out_ss += f' self.name = "{name}"\n' + if formulation == "asf": + out_ss += " self.solution_model = AsymmetricRegularSolution(\n" + else: + out_ss += " self.solution_model = SymmetricRegularSolution(\n" + out_ss += " endmembers=[\n" + for mbr_initialization in mbr_initializations: + out_ss += f" {mbr_initialization},\n" + out_ss += " ],\n" + if formulation == "asf": + out_ss += f" alphas={Ae},\n" + out_ss += f" energy_interaction={We},\n" + if nonideal_entropies: + out_ss += f" entropy_interaction={Ws},\n" + if nonideal_volumes: + out_ss += f" volume_interaction={Wv},\n" + + out_ss += " )\n" + out_ss += " Solution.__init__(self, molar_fractions=molar_fractions)\n\n" + + logging.debug(site_fractions) + logging.debug(endmember_site_fractions_and_checks) + logging.debug("") + + +out_noods = "" +noods = np.unique(sorted(noods)) +for nood in noods: + mineral = getattr(dataset, nood)() + params = mineral.params + out_noods += f"{nood}_nood = Mineral({params})\n\n" + +# Preamble +solution_dataset = solution_file.split("-")[-1].split(".")[0] +underline = "^" * len(solution_dataset) +preamble = ( + "# This file is part of BurnMan - a thermoelastic\n" + "# and thermodynamic toolkit for the Earth and Planetary Sciences\n" + f"# Copyright (C) 2012 - {current_year} by the BurnMan team, released under the GNU\n" + "# GPL v2 or later.\n" + "\n" + '"""\n' + f"{solution_dataset}\n" + f"{underline}\n" + "\n" + "HPx-eos solutions using endmembers from\n" + f"dataset {mbr_dataset}.\n" + "The values in this document are all in S.I. units,\n" + "unlike those in the original THERMOCALC file.\n" + "This file is autogenerated using process_HPX_eos.py\n" + '"""\n\n' + f"from numpy import array, nan\n" + f"from . import {mbr_dataset}\n" + "from ..classes.mineral import Mineral\n" + "from ..classes.solution import Solution\n" + "from ..classes.solutionmodel import SymmetricRegularSolution\n" + "from ..classes.solutionmodel import AsymmetricRegularSolution\n" + "from ..classes.combinedmineral import CombinedMineral\n\n" +) + + +print(preamble) +print(out_noods) +print(out_make) +print(out_ss) diff --git a/burnman/minerals/__init__.py b/burnman/minerals/__init__.py index 85c16cc4..767c920c 100644 --- a/burnman/minerals/__init__.py +++ b/burnman/minerals/__init__.py @@ -52,6 +52,8 @@ from . import HHPH_2013 from . import JH_2015 from . import HGP_2018_ds633 +from . import ig50NCKFMASHTOCr +from . import ig50NCKFMASTOCr # Kurnosov et al. 2017 from . import KMFBZ_2017 diff --git a/burnman/minerals/ig50NCKFMASHTOCr.py b/burnman/minerals/ig50NCKFMASHTOCr.py new file mode 100644 index 00000000..2721545c --- /dev/null +++ b/burnman/minerals/ig50NCKFMASHTOCr.py @@ -0,0 +1,748 @@ +# 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. + +""" +ig50NCKFMASHTOCr +^^^^^^^^^^^^^^^^ + +HPx-eos solutions using endmembers from +dataset HGP_2018_ds633. +The values in this document are all in S.I. units, +unlike those in the original THERMOCALC file. +This file is autogenerated using process_HPX_eos.py +""" + +from numpy import array, nan +from . import HGP_2018_ds633 +from ..classes.mineral import Mineral +from ..classes.solution import Solution +from ..classes.solutionmodel import SymmetricRegularSolution +from ..classes.solutionmodel import AsymmetricRegularSolution +from ..classes.combinedmineral import CombinedMineral + + +cats_nood = Mineral( + { + "name": "cats", + "formula": {"Al": 2.0, "Ca": 1.0, "O": 6.0, "Si": 1.0}, + "equation_of_state": "hp_tmt", + "H_0": -3310050.0, + "S_0": 135.0, + "V_0": 6.356e-05, + "Cp": [347.6, -0.006974, -1781600.0, -2757.5], + "a_0": 2.08e-05, + "K_0": 119200000000.0, + "Kprime_0": 5.19, + "Kdprime_0": -4.4e-11, + "n": 10.0, + "molar_mass": 0.2181229, + "T_0": 298.15, + "E_0": 0.0, + "P_0": 100000.0, + "G_0": nan, + "Gprime_0": nan, + "T_einstein": 533.4002006018054, + } +) + +hem_nood = Mineral( + { + "name": "hem", + "formula": {"Fe": 2.0, "O": 3.0}, + "equation_of_state": "hp_tmt", + "H_0": -825420.0, + "S_0": 87.4, + "V_0": 3.027e-05, + "Cp": [163.9, 0.0, -2257200.0, -657.6], + "a_0": 2.79e-05, + "K_0": 223000000000.0, + "Kprime_0": 4.04, + "Kdprime_0": -1.8e-11, + "n": 5.0, + "molar_mass": 0.1596882, + "T_0": 298.15, + "E_0": 0.0, + "P_0": 100000.0, + "G_0": nan, + "Gprime_0": nan, + "T_einstein": 444.6488294314381, + } +) + +herc_nood = Mineral( + { + "name": "herc", + "formula": {"Al": 2.0, "Fe": 1.0, "O": 4.0}, + "equation_of_state": "hp_tmt", + "H_0": -1949470.0, + "S_0": 113.9, + "V_0": 4.075e-05, + "Cp": [184.9, 0.01417, -3674800.0, -404.0], + "a_0": 2.06e-05, + "K_0": 192200000000.0, + "Kprime_0": 4.04, + "Kdprime_0": -2.1e-11, + "n": 7.0, + "molar_mass": 0.1738056, + "T_0": 298.15, + "E_0": 0.0, + "P_0": 100000.0, + "G_0": nan, + "Gprime_0": nan, + "T_einstein": 468.310479305573, + } +) + +ilm_nood = Mineral( + { + "name": "ilm", + "formula": {"Fe": 1.0, "O": 3.0, "Ti": 1.0}, + "equation_of_state": "hp_tmt", + "H_0": -1232320.0, + "S_0": 107.5, + "V_0": 3.169e-05, + "Cp": [138.9, 0.005081, -1288800.0, -463.7], + "a_0": 2.4e-05, + "K_0": 170000000000.0, + "Kprime_0": 8.3, + "Kdprime_0": -4.9e-11, + "n": 5.0, + "molar_mass": 0.1517102, + "T_0": 298.15, + "E_0": 0.0, + "P_0": 100000.0, + "G_0": nan, + "Gprime_0": nan, + "T_einstein": 380.67287043664993, + } +) + +sp_nood = Mineral( + { + "name": "sp", + "formula": {"Al": 2.0, "Mg": 1.0, "O": 4.0}, + "equation_of_state": "hp_tmt", + "H_0": -2300180.0, + "S_0": 80.63, + "V_0": 3.978e-05, + "Cp": [200.5, 0.006252, -2996400.0, -888.4], + "a_0": 1.93e-05, + "K_0": 192200000000.0, + "Kprime_0": 4.04, + "Kdprime_0": -2.1e-11, + "n": 7.0, + "molar_mass": 0.1422656, + "T_0": 298.15, + "E_0": 0.0, + "P_0": 100000.0, + "G_0": nan, + "Gprime_0": nan, + "T_einstein": 592.252008591202, + } +) + + +abhI = CombinedMineral([HGP_2018_ds633.abh()], [1.0], [570.0, 4.12, 0.0], "abhI") +anC = CombinedMineral([HGP_2018_ds633.an()], [1.0], [7030.0, 4.66, 0.0], "anC") +cfm = CombinedMineral( + [HGP_2018_ds633.fa(), HGP_2018_ds633.fo()], [0.5, 0.5], [0.0, -0.0, 0.0], "cfm" +) +anC = CombinedMineral([HGP_2018_ds633.an()], [1.0], [7030.0, 4.66, 0.0], "anC") +mam = CombinedMineral([HGP_2018_ds633.ma()], [1.0], [6500.0, -0.0, 0.0], "mam") +fmu = CombinedMineral( + [HGP_2018_ds633.andr(), HGP_2018_ds633.gr(), HGP_2018_ds633.mu()], + [0.5, -0.5, 1.0], + [25000.0, -0.0, 0.0], + "fmu", +) +annm = CombinedMineral([HGP_2018_ds633.ann()], [1.0], [-6000.0, -0.0, 0.0], "annm") +obi = CombinedMineral( + [HGP_2018_ds633.ann(), HGP_2018_ds633.phl()], + [0.3333333333333333, 0.6666666666666666], + [-6000.0, -0.0, 0.0], + "obi", +) +tbi = CombinedMineral( + [HGP_2018_ds633.br(), HGP_2018_ds633.phl(), HGP_2018_ds633.ru()], + [-1.0, 1.0, 1.0], + [55000.0, -0.0, 0.0], + "tbi", +) +fbi = CombinedMineral( + [HGP_2018_ds633.andr(), HGP_2018_ds633.east(), HGP_2018_ds633.gr()], + [0.5, 1.0, -0.5], + [-3000.0, -0.0, 0.0], + "fbi", +) +knom = CombinedMineral([HGP_2018_ds633.knor()], [1.0], [18200.0, -0.0, 0.0], "knom") +tig = CombinedMineral( + [ + HGP_2018_ds633.py(), + HGP_2018_ds633.per(), + HGP_2018_ds633.ru(), + HGP_2018_ds633.cor(), + ], + [1.0, 0.5, 0.5, -0.5], + [46700.0, 17.3, 0.0], + "tig", +) +fm = CombinedMineral( + [HGP_2018_ds633.en(), HGP_2018_ds633.fs()], [0.5, 0.5], [-6600.0, -0.0, 0.0], "fm" +) +odi = CombinedMineral( + [HGP_2018_ds633.di()], [1.0], [2800.0, -0.0, 5.0000000000000004e-08], "odi" +) +cren = CombinedMineral( + [HGP_2018_ds633.mgts(), HGP_2018_ds633.kos(), HGP_2018_ds633.jd()], + [1.0, 1.0, -1.0], + [-25900.0, -15.5, 5.000000000000001e-07], + "cren", +) +obuf = CombinedMineral( + [ + HGP_2018_ds633.mgts(), + HGP_2018_ds633.per(), + HGP_2018_ds633.ru(), + HGP_2018_ds633.cor(), + ], + [1.0, 0.5, 0.5, -0.5], + [-5000.0, 5.1000000000000005, -6.1e-08], + "obuf", +) +mess = CombinedMineral( + [HGP_2018_ds633.mgts(), HGP_2018_ds633.acm(), HGP_2018_ds633.jd()], + [1.0, 1.0, -1.0], + [4800.0, -0.0, -8.900000000000001e-07], + "mess", +) +ojd = CombinedMineral([HGP_2018_ds633.jd()], [1.0], [18800.0, -0.0, 0.0], "ojd") +cfs = CombinedMineral( + [HGP_2018_ds633.fs()], [1.0], [2100.0, 2.0, 4.5000000000000003e-07], "cfs" +) +crdi = CombinedMineral( + [cats_nood, HGP_2018_ds633.kos(), HGP_2018_ds633.jd()], + [1.0, 1.0, -1.0], + [-1100.0, 2.881573160768881, 1e-07], + "crdi", +) +cess = CombinedMineral( + [cats_nood, HGP_2018_ds633.acm(), HGP_2018_ds633.jd()], + [1.0, 1.0, -1.0], + [350.0, 2.881573160768881, 1e-07], + "cess", +) +cbuf = CombinedMineral( + [cats_nood, HGP_2018_ds633.per(), HGP_2018_ds633.ru(), HGP_2018_ds633.cor()], + [1.0, 0.5, 0.5, -0.5], + [-12400.0, 4.0815731607688805, 4.999999999999999e-08], + "cbuf", +) +cen = CombinedMineral( + [HGP_2018_ds633.en()], [1.0], [3500.0, 2.0, 4.800000000000001e-07], "cen" +) +cfm = CombinedMineral( + [HGP_2018_ds633.en(), HGP_2018_ds633.fs()], + [0.5, 0.5], + [-1600.0, 2.0, 4.6500000000000005e-07], + "cfm", +) +kjd = CombinedMineral( + [HGP_2018_ds633.jd(), HGP_2018_ds633.abh(), HGP_2018_ds633.san()], + [1.0, -1.0, 1.0], + [11700.0, -0.0, 6e-06], + "kjd", +) +nsp = CombinedMineral([sp_nood], [1.0], [0.0, 0.0, 0.0], "nsp") +isp = CombinedMineral([sp_nood], [1.0], [23600.0, 5.76303, 0.0], "isp") +nhc = CombinedMineral([herc_nood], [1.0], [0.0, 0.0, 0.0], "nhc") +ihc = CombinedMineral([herc_nood], [1.0], [23600.0, 5.76303, 0.0], "ihc") +nmt = CombinedMineral([HGP_2018_ds633.mt()], [1.0], [0.0, -5.76303, 0.0], "nmt") +imt = CombinedMineral([HGP_2018_ds633.mt()], [1.0], [300.0, -0.0, 0.0], "imt") +pcr = CombinedMineral([HGP_2018_ds633.picr()], [1.0], [0.0, -0.0, 0.0], "pcr") +qndm = CombinedMineral([HGP_2018_ds633.qnd()], [1.0], [-30000.0, -0.0, 0.0], "qndm") +tsm = CombinedMineral([HGP_2018_ds633.ts()], [1.0], [10000.0, -0.0, 0.0], "tsm") +prgm = CombinedMineral([HGP_2018_ds633.parg()], [1.0], [-10000.0, -0.0, 0.0], "prgm") +glm = CombinedMineral([HGP_2018_ds633.gl()], [1.0], [-3000.0, -0.0, 0.0], "glm") +grnm = CombinedMineral([HGP_2018_ds633.grun()], [1.0], [-3000.0, -0.0, 0.0], "grnm") +a = CombinedMineral( + [HGP_2018_ds633.cumm(), HGP_2018_ds633.grun()], + [0.42857142857142855, 0.5714285714285714], + [-11200.0, -0.0, 0.0], + "a", +) +b = CombinedMineral( + [HGP_2018_ds633.cumm(), HGP_2018_ds633.grun()], + [0.2857142857142857, 0.7142857142857143], + [-13800.0, -0.0, 0.0], + "b", +) +mrb = CombinedMineral( + [HGP_2018_ds633.gl(), HGP_2018_ds633.gr(), HGP_2018_ds633.andr()], + [1.0, -1.0, 1.0], + [0.0, -0.0, 0.0], + "mrb", +) +kprg = CombinedMineral( + [HGP_2018_ds633.mu(), HGP_2018_ds633.pa(), HGP_2018_ds633.parg()], + [1.0, -1.0, 1.0], + [-7060.0, -20.0, 0.0], + "kprg", +) +tts = CombinedMineral( + [HGP_2018_ds633.dsp(), HGP_2018_ds633.ru(), HGP_2018_ds633.ts()], + [-2.0, 2.0, 1.0], + [95000.0, -0.0, 0.0], + "tts", +) +oilm = CombinedMineral( + [ilm_nood], + [1.0], + [1444.0572257227777, 1.5923196732102785, 1.836386612201713e-07], + "oilm", +) +dilm = CombinedMineral( + [ilm_nood], + [1.0], + [17044.35722572278, 13.118319673210278, 1.836386612201713e-07], + "dilm", +) +dhem = CombinedMineral( + [hem_nood], [1.0], [9522.770803030953, 12.937668369038724, 0.0], "dhem" +) + + +class pl4tr(Solution): + def __init__(self, molar_fractions=None): + self.name = "pl4tr" + self.solution_model = AsymmetricRegularSolution( + endmembers=[ + [HGP_2018_ds633.ab(), "[Naa][Altb1/4Sitb3/4]"], + [HGP_2018_ds633.an(), "[Caa][Altb1/2Sitb1/2]"], + [HGP_2018_ds633.san(), "[Ka][Altb1/4Sitb3/4]"], + ], + alphas=[0.674, 0.55, 1.0], + energy_interaction=[[14600.0, 24100.0], [48500.0]], + entropy_interaction=[[0.00935, 0.00957], [-0.0]], + volume_interaction=[ + [-4.0000000000000003e-07, 3.3800000000000007e-06], + [-1.3e-06], + ], + ) + Solution.__init__(self, molar_fractions=molar_fractions) + + +class k4tr(Solution): + def __init__(self, molar_fractions=None): + self.name = "k4tr" + self.solution_model = AsymmetricRegularSolution( + endmembers=[ + [HGP_2018_ds633.ab(), "[Naa][Altb1/4Sitb3/4]"], + [HGP_2018_ds633.an(), "[Caa][Altb1/2Sitb1/2]"], + [HGP_2018_ds633.san(), "[Ka][Altb1/4Sitb3/4]"], + ], + alphas=[0.674, 0.55, 1.0], + energy_interaction=[[14600.0, 24100.0], [48500.0]], + entropy_interaction=[[0.00935, 0.00957], [-0.0]], + volume_interaction=[ + [-4.0000000000000003e-07, 3.3800000000000007e-06], + [-1.3e-06], + ], + ) + Solution.__init__(self, molar_fractions=molar_fractions) + + +class pli(Solution): + def __init__(self, molar_fractions=None): + self.name = "pli" + self.solution_model = AsymmetricRegularSolution( + endmembers=[ + [abhI, "[Na]"], + [HGP_2018_ds633.an(), "[Ca]"], + [HGP_2018_ds633.san(), "[K]"], + ], + alphas=[0.643, 1.0, 1.0], + energy_interaction=[[15000.0, 25100.0], [40000.0]], + entropy_interaction=[[-0.0, 0.0108], [-0.0]], + volume_interaction=[[0.0, 3.3800000000000007e-06], [0.0]], + ) + Solution.__init__(self, molar_fractions=molar_fractions) + + +class plc(Solution): + def __init__(self, molar_fractions=None): + self.name = "plc" + self.solution_model = AsymmetricRegularSolution( + endmembers=[ + [HGP_2018_ds633.abh(), "[Na]"], + [anC, "[Ca]"], + [HGP_2018_ds633.san(), "[K]"], + ], + alphas=[0.643, 1.0, 1.0], + energy_interaction=[[3100.0, 25100.0], [40000.0]], + entropy_interaction=[[-0.0, 0.0108], [-0.0]], + volume_interaction=[[0.0, 3.3800000000000007e-06], [0.0]], + ) + Solution.__init__(self, molar_fractions=molar_fractions) + + +class ol(Solution): + def __init__(self, molar_fractions=None): + self.name = "ol" + self.solution_model = SymmetricRegularSolution( + endmembers=[ + [HGP_2018_ds633.mont(), "[Mgm][Cam]"], + [HGP_2018_ds633.fa(), "[Fem][Fem]"], + [HGP_2018_ds633.fo(), "[Mgm][Mgm]"], + [cfm, "[Mgm][Fem]"], + ], + energy_interaction=[ + [24000.0, 38000.0, 24000.0], + [9000.0, 4500.0], + [4500.0], + ], + ) + Solution.__init__(self, molar_fractions=molar_fractions) + + +class ksp(Solution): + def __init__(self, molar_fractions=None): + self.name = "ksp" + self.solution_model = AsymmetricRegularSolution( + endmembers=[ + [HGP_2018_ds633.san(), "[K]"], + [HGP_2018_ds633.abh(), "[Na]"], + [anC, "[Ca]"], + ], + alphas=[1.0, 0.643, 1.0], + energy_interaction=[[25100.0, 40000.0], [3100.0]], + entropy_interaction=[[0.0108, -0.0], [-0.0]], + volume_interaction=[[3.3800000000000007e-06, 0.0], [0.0]], + ) + Solution.__init__(self, molar_fractions=molar_fractions) + + +class mu(Solution): + def __init__(self, molar_fractions=None): + self.name = "mu" + self.solution_model = AsymmetricRegularSolution( + endmembers=[ + [HGP_2018_ds633.mu(), "[Ka][Alma][Almb][Sit1/2Alt1/2]2"], + [HGP_2018_ds633.cel(), "[Ka][Mgma][Almb][Sit]2"], + [HGP_2018_ds633.fcel(), "[Ka][Fema][Almb][Sit]2"], + [HGP_2018_ds633.pa(), "[Naa][Alma][Almb][Sit1/2Alt1/2]2"], + [mam, "[Caa][Alma][Almb][Alt]2"], + [fmu, "[Ka][Alma][Femb][Sit1/2Alt1/2]2"], + ], + alphas=[0.63, 0.63, 0.63, 0.37, 0.63, 0.63], + energy_interaction=[ + [0.0, 0.0, 10120.0, 35000.0, 0.0], + [0.0, 45000.0, 50000.0, 0.0], + [45000.0, 50000.0, 0.0], + [15000.0, 30000.0], + [35000.0], + ], + entropy_interaction=[ + [-0.0, -0.0, -0.0034, -0.0, -0.0], + [-0.0, -0.0, -0.0, -0.0], + [-0.0, -0.0, -0.0], + [-0.0, -0.0], + [-0.0], + ], + volume_interaction=[ + [2.0000000000000003e-06, 2.0000000000000003e-06, 3.53e-06, 0.0, 0.0], + [0.0, 2.5e-06, 0.0, 0.0], + [2.5e-06, 0.0, 0.0], + [0.0, 0.0], + [0.0], + ], + ) + Solution.__init__(self, molar_fractions=molar_fractions) + + +class bi(Solution): + def __init__(self, molar_fractions=None): + self.name = "bi" + self.solution_model = SymmetricRegularSolution( + endmembers=[ + [HGP_2018_ds633.phl(), "[Mgm][Mgm]2[Sit1/2Alt1/2]2[Ohv]2"], + [annm, "[Fem][Fem]2[Sit1/2Alt1/2]2[Ohv]2"], + [obi, "[Fem][Mgm]2[Sit1/2Alt1/2]2[Ohv]2"], + [HGP_2018_ds633.east(), "[Alm][Mgm]2[Alt]2[Ohv]2"], + [tbi, "[Tim][Mgm]2[Sit1/2Alt1/2]2[Ov]2"], + [fbi, "[Fem][Mgm]2[Alt]2[Ohv]2"], + ], + energy_interaction=[ + [12000.0, 4000.0, 10000.0, 30000.0, 8000.0], + [8000.0, 5000.0, 32000.0, 13600.0], + [7000.0, 24000.0, 5600.0], + [40000.0, 1000.0], + [40000.0], + ], + ) + Solution.__init__(self, molar_fractions=molar_fractions) + + +class g(Solution): + def __init__(self, molar_fractions=None): + self.name = "g" + self.solution_model = AsymmetricRegularSolution( + endmembers=[ + [HGP_2018_ds633.py(), "[Mgm]3[Alm]2"], + [HGP_2018_ds633.alm(), "[Fem]3[Alm]2"], + [HGP_2018_ds633.gr(), "[Cam]3[Alm]2"], + [HGP_2018_ds633.andr(), "[Cam]3[Fem]2"], + [knom, "[Mgm]3[Crm]2"], + [tig, "[Mgm]3[Alm1/2Mgm1/4Tim1/4]2"], + ], + alphas=[1.0, 1.0, 2.5, 2.5, 1.0, 1.0], + energy_interaction=[ + [4000.0, 45400.0, 107000.0, 2000.0, 0.0], + [17000.0, 65000.0, 6000.0, 0.0], + [2000.0, 1000.0, 0.0], + [63000.0, 0.0], + [0.0], + ], + entropy_interaction=[ + [-0.0, 0.01, 0.01, -0.0, -0.0], + [0.01, 0.01, -0.0, -0.0], + [-0.0, 0.01, -0.0], + [0.01, -0.0], + [-0.0], + ], + volume_interaction=[ + [1.0000000000000002e-06, 4.0000000000000003e-07, -3.6e-07, 0.0, 0.0], + [1.0000000000000002e-06, 3.9e-07, 1.0000000000000001e-07, 0.0], + [0.0, 1.8000000000000001e-06, 0.0], + [1.0000000000000002e-06, 0.0], + [0.0], + ], + ) + Solution.__init__(self, molar_fractions=molar_fractions) + + +class ep(Solution): + def __init__(self, molar_fractions=None): + self.name = "ep" + self.solution_model = SymmetricRegularSolution( + endmembers=[ + [HGP_2018_ds633.cz(), "[Alm][Alm]"], + [HGP_2018_ds633.ep(), "[Alm][Fem]"], + [HGP_2018_ds633.fep(), "[Fem][Fem]"], + ], + energy_interaction=[[1000.0, 3000.0], [1000.0]], + ) + Solution.__init__(self, molar_fractions=molar_fractions) + + +class cd(Solution): + def __init__(self, molar_fractions=None): + self.name = "cd" + self.solution_model = SymmetricRegularSolution( + endmembers=[ + [HGP_2018_ds633.crd(), "[Mgx]2[Vh]"], + [HGP_2018_ds633.fcrd(), "[Fex]2[Vh]"], + [HGP_2018_ds633.hcrd(), "[Mgx]2[Hoh]"], + ], + energy_interaction=[[6000.0, 0.0], [0.0]], + ) + Solution.__init__(self, molar_fractions=molar_fractions) + + +class opx(Solution): + def __init__(self, molar_fractions=None): + self.name = "opx" + self.solution_model = AsymmetricRegularSolution( + endmembers=[ + [HGP_2018_ds633.en(), "[Mgm][Mgm][Sit]1/2"], + [HGP_2018_ds633.fs(), "[Fem][Fem][Sit]1/2"], + [fm, "[Mgm][Fem][Sit]1/2"], + [odi, "[Mgm][Cam][Sit]1/2"], + [HGP_2018_ds633.mgts(), "[Alm][Mgm][Sit1/2Alt1/2]1/2"], + [cren, "[Crm][Mgm][Sit1/2Alt1/2]1/2"], + [obuf, "[Mgm1/2Tim1/2][Mgm][Sit1/2Alt1/2]1/2"], + [mess, "[Fem][Mgm][Sit1/2Alt1/2]1/2"], + [ojd, "[Alm][Nam][Sit]1/2"], + ], + alphas=[1.0, 1.0, 1.0, 1.2, 1.0, 1.0, 1.0, 1.0, 1.2], + energy_interaction=[ + [7000.0, 4000.0, 29400.0, 12500.0, 8000.0, 6000.0, 8000.0, 35000.0], + [4000.0, 21500.0, 11000.0, 10000.0, 7000.0, 10000.0, 35000.0], + [18000.0, 15000.0, 12000.0, 8000.0, 12000.0, 35000.0], + [75500.0, 20000.0, 40000.0, 20000.0, 35000.0], + [2000.0, 10000.0, 2000.0, 7000.0], + [6000.0, 2000.0, -11000.0], + [6000.0, 20000.0], + [-11000.0], + ], + volume_interaction=[ + [0.0, 0.0, 0.0, -4.0000000000000003e-07, 0.0, 0.0, 0.0, 0.0], + [0.0, 8.000000000000001e-07, -1.5e-06, 0.0, 0.0, 0.0, 0.0], + [8.000000000000001e-07, -1.5e-06, 0.0, 0.0, 0.0, 0.0], + [-8.400000000000001e-06, 0.0, 0.0, 0.0, 0.0], + [0.0, 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 cpx(Solution): + def __init__(self, molar_fractions=None): + self.name = "cpx" + self.solution_model = AsymmetricRegularSolution( + endmembers=[ + [HGP_2018_ds633.di(), "[Mgm][Cam][Sit]1/2"], + [cfs, "[Fem][Fem][Sit]1/2"], + [HGP_2018_ds633.cats(), "[Alm][Cam][Sit1/2Alt1/2]1/2"], + [crdi, "[Crm][Cam][Sit1/2Alt1/2]1/2"], + [cess, "[Fem][Cam][Sit1/2Alt1/2]1/2"], + [cbuf, "[Mgm1/2Tim1/2][Cam][Sit1/2Alt1/2]1/2"], + [HGP_2018_ds633.jd(), "[Alm][Nam][Sit]1/2"], + [cen, "[Mgm][Mgm][Sit]1/2"], + [cfm, "[Mgm][Fem][Sit]1/2"], + [kjd, "[Alm][Km][Sit]1/2"], + ], + alphas=[1.2, 1.0, 1.9, 1.9, 1.9, 1.9, 1.2, 1.0, 1.0, 1.2], + energy_interaction=[ + [ + 25800.0, + 13000.0, + 8000.0, + 8000.0, + 8000.0, + 26000.0, + 29800.0, + 20600.0, + 26000.0, + ], + [25000.0, 38300.0, 43300.0, 24000.0, 24000.0, 2300.0, 3500.0, 24000.0], + [2000.0, 2000.0, 6000.0, 6000.0, 45200.0, 27000.0, 6000.0], + [2000.0, 6000.0, 3000.0, 52300.0, 40300.0, 3000.0], + [6000.0, 3000.0, 57300.0, 45300.0, 3000.0], + [16000.0, 24000.0, 22000.0, 16000.0], + [40000.0, 40000.0, 28000.0], + [4000.0, 40000.0], + [40000.0], + ], + volume_interaction=[ + [0.0, -6.000000000000001e-07, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [-1.0000000000000002e-06, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, -3.5e-06, -1.0000000000000002e-06, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0], + [0.0, 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 spn(Solution): + def __init__(self, molar_fractions=None): + self.name = "spn" + self.solution_model = SymmetricRegularSolution( + endmembers=[ + [nsp, "[Mgt][Alm]"], + [isp, "[Alt][Mgm1/2Alm1/2]"], + [nhc, "[Fet][Alm]"], + [ihc, "[Alt][Fem1/2Alm1/2]"], + [nmt, "[Fet][Fem]"], + [imt, "[Fet][Fem1/2Fem1/2]"], + [pcr, "[Mgt][Crm]"], + [qndm, "[Mgt][Mgm1/2Tim1/2]"], + ], + energy_interaction=[ + [-8200.0, 3500.0, -13000.0, 43200.0, 49100.0, -5000.0, 22500.0], + [4400.0, -6000.0, 36800.0, 20000.0, 14000.0, 21500.0], + [-8200.0, 18100.0, 49000.0, -19000.0, 35100.0], + [-4000.0, 7600.0, -11000.0, 9000.0], + [18100.0, 11900.0, 62200.0], + [-6400.0, 24300.0], + [60000.0], + ], + ) + Solution.__init__(self, molar_fractions=molar_fractions) + + +class hb(Solution): + def __init__(self, molar_fractions=None): + self.name = "hb" + self.solution_model = AsymmetricRegularSolution( + endmembers=[ + [HGP_2018_ds633.tr(), "[Va][Mgm]3[Mgm]2[Cam]2[Sit][Ohv]2"], + [tsm, "[Va][Mgm]3[Alm]2[Cam]2[Sit1/2Alt1/2][Ohv]2"], + [prgm, "[Naa][Mgm]3[Mgm1/2Alm1/2]2[Cam]2[Sit1/2Alt1/2][Ohv]2"], + [glm, "[Va][Mgm]3[Alm]2[Nam]2[Sit][Ohv]2"], + [HGP_2018_ds633.cumm(), "[Va][Mgm]3[Mgm]2[Mgm]2[Sit][Ohv]2"], + [grnm, "[Va][Fem]3[Fem]2[Fem]2[Sit][Ohv]2"], + [a, "[Va][Mgm]3[Fem]2[Fem]2[Sit][Ohv]2"], + [b, "[Va][Fem]3[Mgm]2[Fem]2[Sit][Ohv]2"], + [mrb, "[Va][Mgm]3[Fem]2[Nam]2[Sit][Ohv]2"], + [kprg, "[Ka][Mgm]3[Mgm1/2Alm1/2]2[Cam]2[Sit1/2Alt1/2][Ohv]2"], + [tts, "[Va][Mgm]3[Tim]2[Cam]2[Sit1/2Alt1/2][Ov]2"], + ], + alphas=[1.0, 1.5, 1.7, 0.8, 1.0, 1.0, 1.0, 1.0, 0.8, 1.7, 1.5], + energy_interaction=[ + [ + 20000.0, + 25000.0, + 65000.0, + 45000.0, + 75000.0, + 57000.0, + 63000.0, + 52000.0, + 30000.0, + 85000.0, + ], + [ + -40000.0, + 25000.0, + 70000.0, + 80000.0, + 70000.0, + 72500.0, + 20000.0, + -40000.0, + 35000.0, + ], + [ + 50000.0, + 90000.0, + 106700.0, + 94800.0, + 94800.0, + 40000.0, + 8000.0, + 15000.0, + ], + [100000.0, 113500.0, 100000.0, 111200.0, 0.0, 54000.0, 75000.0], + [33000.0, 18000.0, 23000.0, 80000.0, 87000.0, 100000.0], + [12000.0, 8000.0, 91000.0, 96000.0, 65000.0], + [20000.0, 80000.0, 94000.0, 95000.0], + [90000.0, 94000.0, 95000.0], + [50000.0, 50000.0], + [35000.0], + ], + ) + Solution.__init__(self, molar_fractions=molar_fractions) + + +class ilm(Solution): + def __init__(self, molar_fractions=None): + self.name = "ilm" + self.solution_model = SymmetricRegularSolution( + endmembers=[ + [oilm, "[Fea][Tib]"], + [dilm, "[Fea1/2Tia1/2][Feb1/2Tib1/2]"], + [dhem, "[Fea][Feb]"], + ], + energy_interaction=[[15600.0, 26600.0], [11000.0]], + ) + Solution.__init__(self, molar_fractions=molar_fractions) diff --git a/burnman/minerals/ig50NCKFMASTOCr.py b/burnman/minerals/ig50NCKFMASTOCr.py new file mode 100644 index 00000000..e30e2d6b --- /dev/null +++ b/burnman/minerals/ig50NCKFMASTOCr.py @@ -0,0 +1,535 @@ +# 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. + +""" +ig50NCKFMASTOCr +^^^^^^^^^^^^^^^ + +HPx-eos solutions using endmembers from +dataset HGP_2018_ds633. +The values in this document are all in S.I. units, +unlike those in the original THERMOCALC file. +This file is autogenerated using process_HPX_eos.py +""" + +from numpy import array, nan +from . import HGP_2018_ds633 +from ..classes.mineral import Mineral +from ..classes.solution import Solution +from ..classes.solutionmodel import SymmetricRegularSolution +from ..classes.solutionmodel import AsymmetricRegularSolution +from ..classes.combinedmineral import CombinedMineral + + +cats_nood = Mineral( + { + "name": "cats", + "formula": {"Al": 2.0, "Ca": 1.0, "O": 6.0, "Si": 1.0}, + "equation_of_state": "hp_tmt", + "H_0": -3310050.0, + "S_0": 135.0, + "V_0": 6.356e-05, + "Cp": [347.6, -0.006974, -1781600.0, -2757.5], + "a_0": 2.08e-05, + "K_0": 119200000000.0, + "Kprime_0": 5.19, + "Kdprime_0": -4.4e-11, + "n": 10.0, + "molar_mass": 0.2181229, + "T_0": 298.15, + "E_0": 0.0, + "P_0": 100000.0, + "G_0": nan, + "Gprime_0": nan, + "T_einstein": 533.4002006018054, + } +) + +hem_nood = Mineral( + { + "name": "hem", + "formula": {"Fe": 2.0, "O": 3.0}, + "equation_of_state": "hp_tmt", + "H_0": -825420.0, + "S_0": 87.4, + "V_0": 3.027e-05, + "Cp": [163.9, 0.0, -2257200.0, -657.6], + "a_0": 2.79e-05, + "K_0": 223000000000.0, + "Kprime_0": 4.04, + "Kdprime_0": -1.8e-11, + "n": 5.0, + "molar_mass": 0.1596882, + "T_0": 298.15, + "E_0": 0.0, + "P_0": 100000.0, + "G_0": nan, + "Gprime_0": nan, + "T_einstein": 444.6488294314381, + } +) + +herc_nood = Mineral( + { + "name": "herc", + "formula": {"Al": 2.0, "Fe": 1.0, "O": 4.0}, + "equation_of_state": "hp_tmt", + "H_0": -1949470.0, + "S_0": 113.9, + "V_0": 4.075e-05, + "Cp": [184.9, 0.01417, -3674800.0, -404.0], + "a_0": 2.06e-05, + "K_0": 192200000000.0, + "Kprime_0": 4.04, + "Kdprime_0": -2.1e-11, + "n": 7.0, + "molar_mass": 0.1738056, + "T_0": 298.15, + "E_0": 0.0, + "P_0": 100000.0, + "G_0": nan, + "Gprime_0": nan, + "T_einstein": 468.310479305573, + } +) + +ilm_nood = Mineral( + { + "name": "ilm", + "formula": {"Fe": 1.0, "O": 3.0, "Ti": 1.0}, + "equation_of_state": "hp_tmt", + "H_0": -1232320.0, + "S_0": 107.5, + "V_0": 3.169e-05, + "Cp": [138.9, 0.005081, -1288800.0, -463.7], + "a_0": 2.4e-05, + "K_0": 170000000000.0, + "Kprime_0": 8.3, + "Kdprime_0": -4.9e-11, + "n": 5.0, + "molar_mass": 0.1517102, + "T_0": 298.15, + "E_0": 0.0, + "P_0": 100000.0, + "G_0": nan, + "Gprime_0": nan, + "T_einstein": 380.67287043664993, + } +) + +sp_nood = Mineral( + { + "name": "sp", + "formula": {"Al": 2.0, "Mg": 1.0, "O": 4.0}, + "equation_of_state": "hp_tmt", + "H_0": -2300180.0, + "S_0": 80.63, + "V_0": 3.978e-05, + "Cp": [200.5, 0.006252, -2996400.0, -888.4], + "a_0": 1.93e-05, + "K_0": 192200000000.0, + "Kprime_0": 4.04, + "Kdprime_0": -2.1e-11, + "n": 7.0, + "molar_mass": 0.1422656, + "T_0": 298.15, + "E_0": 0.0, + "P_0": 100000.0, + "G_0": nan, + "Gprime_0": nan, + "T_einstein": 592.252008591202, + } +) + + +abhI = CombinedMineral([HGP_2018_ds633.abh()], [1.0], [570.0, 4.12, 0.0], "abhI") +anC = CombinedMineral([HGP_2018_ds633.an()], [1.0], [7030.0, 4.66, 0.0], "anC") +cfm = CombinedMineral( + [HGP_2018_ds633.fa(), HGP_2018_ds633.fo()], [0.5, 0.5], [0.0, -0.0, 0.0], "cfm" +) +anC = CombinedMineral([HGP_2018_ds633.an()], [1.0], [7030.0, 4.66, 0.0], "anC") +knom = CombinedMineral([HGP_2018_ds633.knor()], [1.0], [18200.0, -0.0, 0.0], "knom") +tig = CombinedMineral( + [ + HGP_2018_ds633.py(), + HGP_2018_ds633.per(), + HGP_2018_ds633.ru(), + HGP_2018_ds633.cor(), + ], + [1.0, 0.5, 0.5, -0.5], + [46700.0, 17.3, 0.0], + "tig", +) +fm = CombinedMineral( + [HGP_2018_ds633.en(), HGP_2018_ds633.fs()], [0.5, 0.5], [-6600.0, -0.0, 0.0], "fm" +) +odi = CombinedMineral( + [HGP_2018_ds633.di()], [1.0], [2800.0, -0.0, 5.0000000000000004e-08], "odi" +) +cren = CombinedMineral( + [HGP_2018_ds633.mgts(), HGP_2018_ds633.kos(), HGP_2018_ds633.jd()], + [1.0, 1.0, -1.0], + [-25900.0, -15.5, 5.000000000000001e-07], + "cren", +) +obuf = CombinedMineral( + [ + HGP_2018_ds633.mgts(), + HGP_2018_ds633.per(), + HGP_2018_ds633.ru(), + HGP_2018_ds633.cor(), + ], + [1.0, 0.5, 0.5, -0.5], + [-5000.0, 5.1000000000000005, -6.1e-08], + "obuf", +) +mess = CombinedMineral( + [HGP_2018_ds633.mgts(), HGP_2018_ds633.acm(), HGP_2018_ds633.jd()], + [1.0, 1.0, -1.0], + [4800.0, -0.0, -8.900000000000001e-07], + "mess", +) +ojd = CombinedMineral([HGP_2018_ds633.jd()], [1.0], [18800.0, -0.0, 0.0], "ojd") +cfs = CombinedMineral( + [HGP_2018_ds633.fs()], [1.0], [2100.0, 2.0, 4.5000000000000003e-07], "cfs" +) +crdi = CombinedMineral( + [cats_nood, HGP_2018_ds633.kos(), HGP_2018_ds633.jd()], + [1.0, 1.0, -1.0], + [-1100.0, 2.881573160768881, 1e-07], + "crdi", +) +cess = CombinedMineral( + [cats_nood, HGP_2018_ds633.acm(), HGP_2018_ds633.jd()], + [1.0, 1.0, -1.0], + [350.0, 2.881573160768881, 1e-07], + "cess", +) +cbuf = CombinedMineral( + [cats_nood, HGP_2018_ds633.per(), HGP_2018_ds633.ru(), HGP_2018_ds633.cor()], + [1.0, 0.5, 0.5, -0.5], + [-12400.0, 4.0815731607688805, 4.999999999999999e-08], + "cbuf", +) +cen = CombinedMineral( + [HGP_2018_ds633.en()], [1.0], [3500.0, 2.0, 4.800000000000001e-07], "cen" +) +cfm = CombinedMineral( + [HGP_2018_ds633.en(), HGP_2018_ds633.fs()], + [0.5, 0.5], + [-1600.0, 2.0, 4.6500000000000005e-07], + "cfm", +) +kjd = CombinedMineral( + [HGP_2018_ds633.jd(), HGP_2018_ds633.abh(), HGP_2018_ds633.san()], + [1.0, -1.0, 1.0], + [11700.0, -0.0, 6e-06], + "kjd", +) +nsp = CombinedMineral([sp_nood], [1.0], [0.0, 0.0, 0.0], "nsp") +isp = CombinedMineral([sp_nood], [1.0], [23600.0, 5.76303, 0.0], "isp") +nhc = CombinedMineral([herc_nood], [1.0], [0.0, 0.0, 0.0], "nhc") +ihc = CombinedMineral([herc_nood], [1.0], [23600.0, 5.76303, 0.0], "ihc") +nmt = CombinedMineral([HGP_2018_ds633.mt()], [1.0], [0.0, -5.76303, 0.0], "nmt") +imt = CombinedMineral([HGP_2018_ds633.mt()], [1.0], [300.0, -0.0, 0.0], "imt") +pcr = CombinedMineral([HGP_2018_ds633.picr()], [1.0], [0.0, -0.0, 0.0], "pcr") +qndm = CombinedMineral([HGP_2018_ds633.qnd()], [1.0], [-30000.0, -0.0, 0.0], "qndm") +oilm = CombinedMineral( + [ilm_nood], + [1.0], + [1444.0572257227777, 1.5923196732102785, 1.836386612201713e-07], + "oilm", +) +dilm = CombinedMineral( + [ilm_nood], + [1.0], + [17044.35722572278, 13.118319673210278, 1.836386612201713e-07], + "dilm", +) +dhem = CombinedMineral( + [hem_nood], [1.0], [9522.770803030953, 12.937668369038724, 0.0], "dhem" +) + + +class pl4tr(Solution): + def __init__(self, molar_fractions=None): + self.name = "pl4tr" + self.solution_model = AsymmetricRegularSolution( + endmembers=[ + [HGP_2018_ds633.ab(), "[Naa][Altb1/4Sitb3/4]"], + [HGP_2018_ds633.an(), "[Caa][Altb1/2Sitb1/2]"], + [HGP_2018_ds633.san(), "[Ka][Altb1/4Sitb3/4]"], + ], + alphas=[0.674, 0.55, 1.0], + energy_interaction=[[14600.0, 24100.0], [48500.0]], + entropy_interaction=[[0.00935, 0.00957], [-0.0]], + volume_interaction=[ + [-4.0000000000000003e-07, 3.3800000000000007e-06], + [-1.3e-06], + ], + ) + Solution.__init__(self, molar_fractions=molar_fractions) + + +class k4tr(Solution): + def __init__(self, molar_fractions=None): + self.name = "k4tr" + self.solution_model = AsymmetricRegularSolution( + endmembers=[ + [HGP_2018_ds633.ab(), "[Naa][Altb1/4Sitb3/4]"], + [HGP_2018_ds633.an(), "[Caa][Altb1/2Sitb1/2]"], + [HGP_2018_ds633.san(), "[Ka][Altb1/4Sitb3/4]"], + ], + alphas=[0.674, 0.55, 1.0], + energy_interaction=[[14600.0, 24100.0], [48500.0]], + entropy_interaction=[[0.00935, 0.00957], [-0.0]], + volume_interaction=[ + [-4.0000000000000003e-07, 3.3800000000000007e-06], + [-1.3e-06], + ], + ) + Solution.__init__(self, molar_fractions=molar_fractions) + + +class pli(Solution): + def __init__(self, molar_fractions=None): + self.name = "pli" + self.solution_model = AsymmetricRegularSolution( + endmembers=[ + [abhI, "[Na]"], + [HGP_2018_ds633.an(), "[Ca]"], + [HGP_2018_ds633.san(), "[K]"], + ], + alphas=[0.643, 1.0, 1.0], + energy_interaction=[[15000.0, 25100.0], [40000.0]], + entropy_interaction=[[-0.0, 0.0108], [-0.0]], + volume_interaction=[[0.0, 3.3800000000000007e-06], [0.0]], + ) + Solution.__init__(self, molar_fractions=molar_fractions) + + +class plc(Solution): + def __init__(self, molar_fractions=None): + self.name = "plc" + self.solution_model = AsymmetricRegularSolution( + endmembers=[ + [HGP_2018_ds633.abh(), "[Na]"], + [anC, "[Ca]"], + [HGP_2018_ds633.san(), "[K]"], + ], + alphas=[0.643, 1.0, 1.0], + energy_interaction=[[3100.0, 25100.0], [40000.0]], + entropy_interaction=[[-0.0, 0.0108], [-0.0]], + volume_interaction=[[0.0, 3.3800000000000007e-06], [0.0]], + ) + Solution.__init__(self, molar_fractions=molar_fractions) + + +class ol(Solution): + def __init__(self, molar_fractions=None): + self.name = "ol" + self.solution_model = SymmetricRegularSolution( + endmembers=[ + [HGP_2018_ds633.mont(), "[Mgm][Cam]"], + [HGP_2018_ds633.fa(), "[Fem][Fem]"], + [HGP_2018_ds633.fo(), "[Mgm][Mgm]"], + [cfm, "[Mgm][Fem]"], + ], + energy_interaction=[ + [24000.0, 38000.0, 24000.0], + [9000.0, 4500.0], + [4500.0], + ], + ) + Solution.__init__(self, molar_fractions=molar_fractions) + + +class ksp(Solution): + def __init__(self, molar_fractions=None): + self.name = "ksp" + self.solution_model = AsymmetricRegularSolution( + endmembers=[ + [HGP_2018_ds633.san(), "[K]"], + [HGP_2018_ds633.abh(), "[Na]"], + [anC, "[Ca]"], + ], + alphas=[1.0, 0.643, 1.0], + energy_interaction=[[25100.0, 40000.0], [3100.0]], + entropy_interaction=[[0.0108, -0.0], [-0.0]], + volume_interaction=[[3.3800000000000007e-06, 0.0], [0.0]], + ) + Solution.__init__(self, molar_fractions=molar_fractions) + + +class g(Solution): + def __init__(self, molar_fractions=None): + self.name = "g" + self.solution_model = AsymmetricRegularSolution( + endmembers=[ + [HGP_2018_ds633.py(), "[Mgm]3[Alm]2"], + [HGP_2018_ds633.alm(), "[Fem]3[Alm]2"], + [HGP_2018_ds633.gr(), "[Cam]3[Alm]2"], + [HGP_2018_ds633.andr(), "[Cam]3[Fem]2"], + [knom, "[Mgm]3[Crm]2"], + [tig, "[Mgm]3[Alm1/2Mgm1/4Tim1/4]2"], + ], + alphas=[1.0, 1.0, 2.5, 2.5, 1.0, 1.0], + energy_interaction=[ + [4000.0, 45400.0, 107000.0, 2000.0, 0.0], + [17000.0, 65000.0, 6000.0, 0.0], + [2000.0, 1000.0, 0.0], + [63000.0, 0.0], + [0.0], + ], + entropy_interaction=[ + [-0.0, 0.01, 0.01, -0.0, -0.0], + [0.01, 0.01, -0.0, -0.0], + [-0.0, 0.01, -0.0], + [0.01, -0.0], + [-0.0], + ], + volume_interaction=[ + [1.0000000000000002e-06, 4.0000000000000003e-07, -3.6e-07, 0.0, 0.0], + [1.0000000000000002e-06, 3.9e-07, 1.0000000000000001e-07, 0.0], + [0.0, 1.8000000000000001e-06, 0.0], + [1.0000000000000002e-06, 0.0], + [0.0], + ], + ) + Solution.__init__(self, molar_fractions=molar_fractions) + + +class opx(Solution): + def __init__(self, molar_fractions=None): + self.name = "opx" + self.solution_model = AsymmetricRegularSolution( + endmembers=[ + [HGP_2018_ds633.en(), "[Mgm][Mgm][Sit]1/2"], + [HGP_2018_ds633.fs(), "[Fem][Fem][Sit]1/2"], + [fm, "[Mgm][Fem][Sit]1/2"], + [odi, "[Mgm][Cam][Sit]1/2"], + [HGP_2018_ds633.mgts(), "[Alm][Mgm][Sit1/2Alt1/2]1/2"], + [cren, "[Crm][Mgm][Sit1/2Alt1/2]1/2"], + [obuf, "[Mgm1/2Tim1/2][Mgm][Sit1/2Alt1/2]1/2"], + [mess, "[Fem][Mgm][Sit1/2Alt1/2]1/2"], + [ojd, "[Alm][Nam][Sit]1/2"], + ], + alphas=[1.0, 1.0, 1.0, 1.2, 1.0, 1.0, 1.0, 1.0, 1.2], + energy_interaction=[ + [7000.0, 4000.0, 29400.0, 12500.0, 8000.0, 6000.0, 8000.0, 35000.0], + [4000.0, 21500.0, 11000.0, 10000.0, 7000.0, 10000.0, 35000.0], + [18000.0, 15000.0, 12000.0, 8000.0, 12000.0, 35000.0], + [75500.0, 20000.0, 40000.0, 20000.0, 35000.0], + [2000.0, 10000.0, 2000.0, 7000.0], + [6000.0, 2000.0, -11000.0], + [6000.0, 20000.0], + [-11000.0], + ], + volume_interaction=[ + [0.0, 0.0, 0.0, -4.0000000000000003e-07, 0.0, 0.0, 0.0, 0.0], + [0.0, 8.000000000000001e-07, -1.5e-06, 0.0, 0.0, 0.0, 0.0], + [8.000000000000001e-07, -1.5e-06, 0.0, 0.0, 0.0, 0.0], + [-8.400000000000001e-06, 0.0, 0.0, 0.0, 0.0], + [0.0, 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 cpx(Solution): + def __init__(self, molar_fractions=None): + self.name = "cpx" + self.solution_model = AsymmetricRegularSolution( + endmembers=[ + [HGP_2018_ds633.di(), "[Mgm][Cam][Sit]1/2"], + [cfs, "[Fem][Fem][Sit]1/2"], + [HGP_2018_ds633.cats(), "[Alm][Cam][Sit1/2Alt1/2]1/2"], + [crdi, "[Crm][Cam][Sit1/2Alt1/2]1/2"], + [cess, "[Fem][Cam][Sit1/2Alt1/2]1/2"], + [cbuf, "[Mgm1/2Tim1/2][Cam][Sit1/2Alt1/2]1/2"], + [HGP_2018_ds633.jd(), "[Alm][Nam][Sit]1/2"], + [cen, "[Mgm][Mgm][Sit]1/2"], + [cfm, "[Mgm][Fem][Sit]1/2"], + [kjd, "[Alm][Km][Sit]1/2"], + ], + alphas=[1.2, 1.0, 1.9, 1.9, 1.9, 1.9, 1.2, 1.0, 1.0, 1.2], + energy_interaction=[ + [ + 25800.0, + 13000.0, + 8000.0, + 8000.0, + 8000.0, + 26000.0, + 29800.0, + 20600.0, + 26000.0, + ], + [25000.0, 38300.0, 43300.0, 24000.0, 24000.0, 2300.0, 3500.0, 24000.0], + [2000.0, 2000.0, 6000.0, 6000.0, 45200.0, 27000.0, 6000.0], + [2000.0, 6000.0, 3000.0, 52300.0, 40300.0, 3000.0], + [6000.0, 3000.0, 57300.0, 45300.0, 3000.0], + [16000.0, 24000.0, 22000.0, 16000.0], + [40000.0, 40000.0, 28000.0], + [4000.0, 40000.0], + [40000.0], + ], + volume_interaction=[ + [0.0, -6.000000000000001e-07, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [-1.0000000000000002e-06, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, -3.5e-06, -1.0000000000000002e-06, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0], + [0.0, 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 spn(Solution): + def __init__(self, molar_fractions=None): + self.name = "spn" + self.solution_model = SymmetricRegularSolution( + endmembers=[ + [nsp, "[Mgt][Alm]"], + [isp, "[Alt][Mgm1/2Alm1/2]"], + [nhc, "[Fet][Alm]"], + [ihc, "[Alt][Fem1/2Alm1/2]"], + [nmt, "[Fet][Fem]"], + [imt, "[Fet][Fem1/2Fem1/2]"], + [pcr, "[Mgt][Crm]"], + [qndm, "[Mgt][Mgm1/2Tim1/2]"], + ], + energy_interaction=[ + [-8200.0, 3500.0, -13000.0, 43200.0, 49100.0, -5000.0, 22500.0], + [4400.0, -6000.0, 36800.0, 20000.0, 14000.0, 21500.0], + [-8200.0, 18100.0, 49000.0, -19000.0, 35100.0], + [-4000.0, 7600.0, -11000.0, 9000.0], + [18100.0, 11900.0, 62200.0], + [-6400.0, 24300.0], + [60000.0], + ], + ) + Solution.__init__(self, molar_fractions=molar_fractions) + + +class ilm(Solution): + def __init__(self, molar_fractions=None): + self.name = "ilm" + self.solution_model = SymmetricRegularSolution( + endmembers=[ + [oilm, "[Fea][Tib]"], + [dilm, "[Fea1/2Tia1/2][Feb1/2Tib1/2]"], + [dhem, "[Fea][Feb]"], + ], + energy_interaction=[[15600.0, 26600.0], [11000.0]], + ) + Solution.__init__(self, molar_fractions=molar_fractions) diff --git a/tests/test_hp_solutions.py b/tests/test_hp_solutions.py new file mode 100644 index 00000000..a6274208 --- /dev/null +++ b/tests/test_hp_solutions.py @@ -0,0 +1,49 @@ +from __future__ import absolute_import +from __future__ import print_function +import unittest +from util import BurnManTest +from burnman import Mineral, CombinedMineral, Solution +from burnman.minerals import ig50NCKFMASHTOCr, ig50NCKFMASTOCr +import inspect +import numpy as np +from collections import Counter + + +class hp_solutions(BurnManTest): + def test_ig50NCKFMASHTOCr(self): + for m in dir(ig50NCKFMASHTOCr): + mineral = getattr(ig50NCKFMASHTOCr, m) + if ( + inspect.isclass(mineral) + and mineral != Mineral + and mineral != CombinedMineral + and issubclass(mineral, Mineral) + ): + if issubclass(mineral, Solution) and mineral is not Solution: + mbr_names = mineral().endmember_names + n_mbrs = len(mbr_names) + molar_fractions = np.ones(n_mbrs) / n_mbrs + m = mineral() + m.set_composition(molar_fractions) + self.assertTrue(type(m.formula) is Counter) + + def test_ig50NCKFMASTOCr(self): + for m in dir(ig50NCKFMASTOCr): + mineral = getattr(ig50NCKFMASTOCr, m) + if ( + inspect.isclass(mineral) + and mineral != Mineral + and mineral != CombinedMineral + and issubclass(mineral, Mineral) + ): + if issubclass(mineral, Solution) and mineral is not Solution: + mbr_names = mineral().endmember_names + n_mbrs = len(mbr_names) + molar_fractions = np.ones(n_mbrs) / n_mbrs + m = mineral() + m.set_composition(molar_fractions) + self.assertTrue(type(m.formula) is Counter) + + +if __name__ == "__main__": + unittest.main()