diff --git a/burnman/classes/solution.py b/burnman/classes/solution.py index 66c44e26..8c8d131f 100644 --- a/burnman/classes/solution.py +++ b/burnman/classes/solution.py @@ -131,22 +131,20 @@ def formula(self): @material_property def site_occupancies(self): """ - :returns: The fractional occupancies of species on each site. - :rtype: list of OrderedDicts + :returns: The fractional occupancies of species on each site + as a flat vector. + :rtype: numpy.ndarray """ f = self.molar_fractions - occs = np.einsum("ij, i", self.solution_model.endmember_occupancies, f) - site_occs = [] - k = 0 - for site in self.solution_model.sites: - site_occs.append(OrderedDict()) - for species in site: - site_occs[-1][species] = occs[k] - k += 1 - - return site_occs + occs = f.dot(self.solution_model.endmember_occupancies) + return occs - def site_formula(self, precision=2): + def site_formula( + self, + precision=2, + print_absent_species=True, + site_occupancies=None, + ): """ Returns the molar chemical formula of the solution with site occupancies. For example, [Mg0.4Fe0.6]2SiO4. @@ -154,18 +152,41 @@ def site_formula(self, precision=2): :param precision: Precision with which to print the site occupancies :type precision: int + :param print_absent_species: If True, prints site occupancies that are absent (i.e., zero). + :type print_absent_species: bool + + :param site_occupancies: Optional site occupancies to use instead of + the current ones. + :type site_occupancies: numpy.ndarray + :returns: Molar chemical formula of the solution with site occupancies :rtype: str """ - split_empty = self.solution_model.empty_formula.split("[") - formula = split_empty[0] - for i, site_occs in enumerate(self.site_occupancies): + if site_occupancies is not None: + occs = site_occupancies + else: + occs = self.site_occupancies + + i = 0 + # Split the empty formula into site sections. + # This allows us to reconstruct the formula including + # occupancies and anything outside the sites. + split_formula = self.solution_model.empty_formula.split("[") + formula = split_formula[0] + for j, site in enumerate(self.solution_model.sites): formula += "[" - for species, occ in site_occs.items(): - if np.abs(occ) < 1.0e-12: - occ = np.abs(occ) + for species in site: + # avoid printing -0.00 + occ = occs[i] if np.abs(occs[i]) >= 1.0e-12 else 0.0 + + # skip absent species if requested + if np.abs(occ) < 1.0e-12 and not print_absent_species: + i += 1 + continue formula += f"{species}{occ:0.{precision}f}" - formula += split_empty[i + 1] + i += 1 + formula += split_formula[j + 1] + return formula @material_property diff --git a/burnman/tools/thermobarometry.py b/burnman/tools/thermobarometry.py index b38cd8f7..6e534d15 100644 --- a/burnman/tools/thermobarometry.py +++ b/burnman/tools/thermobarometry.py @@ -107,6 +107,12 @@ def estimate_conditions( for a given mineral assemblage. Algorithm based on Powell and Holland (1994). :param assemblage: The mineral assemblage for which to perform the inversion. + Each solution phase in the assemblage must have its composition set along with + its compositional covariance matrix (called `compositional_covariances`). + If there are compositional degrees of freedom, they can be added by setting + the attribute `free_compositional_vectors` on the relevant solution phases, where each + row of the array corresponds to a free compositional vector, and the columns correspond + to the amounts of endmembers of that phase in each vector. :type assemblage: Assemblage :param dataset_covariances: The covariance data from the thermodynamic dataset. @@ -146,13 +152,21 @@ def estimate_conditions( and the fit quality (fit, given by the square root of the reduced chi-squared value). :rtype: OptimizeResult """ - # TODO: Implement compositional unknown parameters like redox state - # to join P and T in the optimization + # Count the number of free compositional vectors across all phases + n_free_vectors = 0 + for phase in assemblage.phases: + if hasattr(phase, "free_compositional_vectors"): + n_free_vectors += phase.free_compositional_vectors.shape[0] + phase.baseline_composition = phase.molar_fractions.copy() + + # Check if P and/or T are fixed P_fixed = np.isclose(pressure_bounds[0], pressure_bounds[1]) T_fixed = np.isclose(temperature_bounds[0], temperature_bounds[1]) - if P_fixed and T_fixed: - raise Exception("Both pressure and temperature cannot be fixed!") + if P_fixed and T_fixed and n_free_vectors == 0: + raise Exception( + "Both pressure and temperature cannot be fixed if there are no free compositional vectors!" + ) # Build the reaction matrix R, possibly reducing the number of endmembers # by excluding components with small molar fractions @@ -189,7 +203,18 @@ def estimate_conditions( raise Exception("No independent reactions found for the given assemblage!") def calculate_cov_a(assemblage, theta, dataset_covariances): - assemblage.set_state(*theta) + if n_free_vectors > 0: + vi = 2 + for phase in assemblage.phases: + if hasattr(phase, "free_compositional_vectors"): + n_free_vectors_in_phase = phase.free_compositional_vectors.shape[0] + v = np.array(theta[vi : vi + n_free_vectors_in_phase]) + dc = phase.free_compositional_vectors.T.dot(v) + new_composition = phase.baseline_composition + dc + phase.set_composition(new_composition) + vi += n_free_vectors_in_phase + + assemblage.set_state(*theta[:2]) cov1 = _build_endmember_covariance_matrix(assemblage, dataset_covariances) cov2 = _build_solution_covariance_matrix(assemblage) cov_mu = cov1 + cov2 @@ -207,18 +232,36 @@ def chisqr(args): """ # These are the chemical potential covariances for all endmembers cov_a = calculate_cov_a( - assemblage, [args[0] * P_scaling, args[1]], dataset_covariances + assemblage, [args[0] * P_scaling, *args[1:]], dataset_covariances ) a = R.dot(assemblage.endmember_partial_gibbs) - chisqr = a.T.dot(np.linalg.inv(cov_a)).dot(a) + chisqr = a.T.dot(np.linalg.pinv(cov_a)).dot(a) return chisqr # Minimize the chi-squared function - x0 = [guessed_conditions[0] / P_scaling, guessed_conditions[1]] + x0 = [ + guessed_conditions[0] / P_scaling, + guessed_conditions[1], + *[0.0] * n_free_vectors, + ] bounds = [ (pressure_bounds[0] / P_scaling, pressure_bounds[1] / P_scaling), (temperature_bounds[0], temperature_bounds[1]), + *[(None, None)] * n_free_vectors, ] + + if n_free_vectors > 0: + # Nelder-Mead is more robust for problems where the feasible region + # may be complex due to compositional degrees of freedom + res = minimize( + chisqr, + x0, + method="Nelder-Mead", + bounds=bounds, + options={"maxiter": max_it}, + ) + x0 = res.x + res = minimize( chisqr, x0, diff --git a/burnman/utils/chemistry.py b/burnman/utils/chemistry.py index 579511a7..078f6551 100644 --- a/burnman/utils/chemistry.py +++ b/burnman/utils/chemistry.py @@ -561,11 +561,20 @@ def ordered_compositional_array(formulae, elements): return formula_array -def formula_to_string(formula): +def formula_to_string(formula, use_fractions=True, significant_figures=2): """ :param formula: Chemical formula :type formula: dict or Counter + :param use_fractions: Whether or not to use fractions in the output string. + If False, numbers are rounded to the number of significant figures + given by the significant_figures parameter. + :type use_fractions: bool + + :param significant_figures: Number of significant figures to use if + use_fractions=False. + :type significant_figures: int + :returns: A formula string, with element order as given in the list IUPAC_element_order. If one or more keys in the dictionary are not one of the elements @@ -574,21 +583,40 @@ def formula_to_string(formula): """ formula_string = "" - for e in IUPAC_element_order: - if e in formula and np.abs(formula[e]) > 1.0e-12: - if np.abs(formula[e] - 1.0) < 1.0e-12: - formula_string += e - else: - formula_string += e + str(nsimplify(formula[e])) - - for e in formula: - if e not in IUPAC_element_order: + if use_fractions: + for e in IUPAC_element_order: if e in formula and np.abs(formula[e]) > 1.0e-12: if np.abs(formula[e] - 1.0) < 1.0e-12: formula_string += e else: formula_string += e + str(nsimplify(formula[e])) + for e in formula: + if e not in IUPAC_element_order: + if e in formula and np.abs(formula[e]) > 1.0e-12: + if np.abs(formula[e] - 1.0) < 1.0e-12: + formula_string += e + else: + formula_string += e + str(nsimplify(formula[e])) + + else: + for e in IUPAC_element_order: + if e in formula and np.abs(formula[e]) > 1.0e-12: + if np.abs(formula[e] - 1.0) < 1.0e-12: + formula_string += e + else: + formula_string += e + str(round(formula[e], significant_figures)) + + for e in formula: + if e not in IUPAC_element_order: + if e in formula and np.abs(formula[e]) > 1.0e-12: + if np.abs(formula[e] - 1.0) < 1.0e-12: + formula_string += e + else: + formula_string += e + str( + round(formula[e], significant_figures) + ) + return formula_string diff --git a/examples/example_optimal_thermobarometry.py b/examples/example_optimal_thermobarometry.py index 449d6e55..7a22e490 100644 --- a/examples/example_optimal_thermobarometry.py +++ b/examples/example_optimal_thermobarometry.py @@ -37,6 +37,7 @@ from burnman.minerals import mp50MnNCKFMASHTO, HP_2011_ds62 from burnman.optimize.composition_fitting import fit_composition_to_solution from burnman.tools.thermobarometry import estimate_conditions +from burnman.utils.chemistry import formula_to_string if __name__ == "__main__": @@ -70,7 +71,9 @@ mu.compositional_covariances = pcov # b) Biotite (requires parameter constraining order state of - # Mg and Fe on the first two sites) + # Mg and Fe on the first two sites). + # Here, we will allow the Mg-Fe ordering to vary to improve the fit + # by adding a free compositional vector. fitted_species = ["Mn", "Fe", "Mg", "Al", "Si", "Ti", "Mg_M3"] species_amounts = np.array([0.01, 1.50, 1.00, 1.65, 2.65, 0.20, 0.10]) species_covariances = np.diag( @@ -82,6 +85,7 @@ ) bi.set_composition(popt) bi.compositional_covariances = pcov + bi.free_compositional_vectors = Composite([bi]).reaction_basis # c) Garnet fitted_species = ["Mn", "Fe", "Mg", "Ca", "Al", "Si"] @@ -94,6 +98,8 @@ g.compositional_covariances = pcov # d) Ilmenite (requires parameter constraining order state of Fe and Ti) + # Here, we will allow the Fe-Ti ordering to vary to improve the fit + # by adding a free compositional vector. fitted_species = ["Mn", "Fe", "Ti", "Mg", "Fe2+_A"] species_amounts = np.array([0.05, 1.0, 0.90, 0.05, 0.4]) species_covariances = np.diag(np.array([0.01, 0.01, 0.01, 0.01, 0.2]) ** 2) @@ -103,6 +109,7 @@ ) ilmm.set_composition(popt) ilmm.compositional_covariances = pcov + ilmm.free_compositional_vectors = Composite([ilmm]).reaction_basis # e) Staurolite fitted_species = ["Mn", "Fe", "Mg", "Al", "Si", "Ti"] @@ -128,20 +135,74 @@ # 4) Print the estimated conditions along with their uncertainties # and correlations. + print("Results of thermobarometric inversion:") print( - f"Estimated Pressure: {assemblage.pressure/1.e9:.2f} " + f" Estimated Pressure: {assemblage.pressure/1.e9:.2f} " f"+/- {np.sqrt(res.xcov[0, 0])/1.e9:.2f} GPa" ) print( - f"Estimated Temperature: {assemblage.temperature:.0f} " + f" Estimated Temperature: {assemblage.temperature:.0f} " f"+/- {np.sqrt(res.xcov[1, 1]):.0f} K" ) - print(f"Correlation between P and T: {res.xcorr:.2f}") - print(f"Number of Reactions: {res.n_reactions}") - print(f"Number of Parameters: {res.n_params}") - print(f"Degrees of Freedom: {res.degrees_of_freedom}") - print(f"Reduced Chi-squared: {res.reduced_chisqr:.2f}") - print(f"Fit (sqrt reduced chi-squared): {res.fit:.2f}") + print(f" Correlation between P and T: {res.xcorr:.2f}") + print(f" Number of Reactions: {res.n_reactions}") + print(f" Number of Parameters: {res.n_params}") + print(f" Degrees of Freedom: {res.degrees_of_freedom}") + print(f" Reduced Chi-squared: {res.reduced_chisqr:.2f}") + print(f" Fit (sqrt reduced chi-squared): {res.fit:.2f}") + + # The names of the sites in these models are a bit cumbersome, + # so we will do some string replacements to make them more + # readable in the output. + replace_strings = [ + ["Fethree", "Fe3+"], + ["monetwo", ""], + ["mtwoa", ""], + ["mtwob", ""], + ["mthree", ""], + ["x", ""], + ["y", ""], + ["tone", ""], + ["t", ""], + ["Ka", "K"], + ["Naa", "Na"], + ["Caa", "Ca"], + ["Fea", "Fe"], + ["Mga", "Mg"], + ["Mna", "Mn"], + ["Tia", "Ti"], + ["Fe3+a", "Fe3+"], + ["Oh", "(OH)"], + ["b", ""], + ["v", ""], + ] + + print("\nMineral compositions and site occupancies:") + for phase in assemblage.phases: + np.set_printoptions(precision=2, formatter={"all": lambda x: f"{x:0.2f}"}) + print(f"{phase.name}:") + + if hasattr(phase, "free_compositional_vectors"): + free_site_occupancy_vectors = phase.free_compositional_vectors.dot( + phase.solution_model.endmember_occupancies + ) + for v in free_site_occupancy_vectors: + site_formula = phase.site_formula( + precision=2, site_occupancies=v, print_absent_species=False + ) + for old, new in replace_strings: + site_formula = site_formula.replace(old, new) + print(f" free site occupancy vector: {site_formula}") + + print( + f" composition: {formula_to_string(phase.formula, use_fractions=False)}" + ) + if hasattr(phase, "molar_fractions"): + site_formula = phase.site_formula(2) + for old, new in replace_strings: + site_formula = site_formula.replace(old, new) + print(f" site occupancies: {site_formula}") + print("") # Uncomment to see the weighted reaction affinities # print("Weighted reaction affinities:") diff --git a/misc/ref/example_optimal_thermobarometry.py.out b/misc/ref/example_optimal_thermobarometry.py.out index e94427c3..cb704478 100644 --- a/misc/ref/example_optimal_thermobarometry.py.out +++ b/misc/ref/example_optimal_thermobarometry.py.out @@ -1,8 +1,31 @@ -Estimated Pressure: 0.43 +/- 0.04 GPa -Estimated Temperature: 887 +/- 11 K -Correlation between P and T: 0.11 -Number of Reactions: 11 -Number of Parameters: 2 -Degrees of Freedom: 9 -Reduced Chi-squared: 0.20 -Fit (sqrt reduced chi-squared): 0.45 +Results of thermobarometric inversion: + Estimated Pressure: 0.43 +/- 0.04 GPa + Estimated Temperature: 888 +/- 11 K + Correlation between P and T: 0.11 + Number of Reactions: 11 + Number of Parameters: 4 + Degrees of Freedom: 7 + Reduced Chi-squared: 0.25 + Fit (sqrt reduced chi-squared): 0.50 + +Mineral compositions and site occupancies: +mu: + composition: K0.56Na0.41Ca0.02Mg0.01Fe0.01Al2.99Si2.99O12.0H2.0 + site occupancies: [K0.56Na0.41Ca0.02][Al0.98Mg0.01Fe0.01][Al1.00Fe3+0.00][Si0.50Al0.50]2 +bi: + free site occupancy vector: [Mg-0.67Fe0.67][Mg0.33Fe-0.33]2[]2[]2 + composition: KMg1.0Ti0.2Mn0.01Fe1.5Al1.65Si2.65O12.0H1.6 + site occupancies: [Mg0.12Fe0.32Al0.29Ti0.20Fe3+0.06Mn0.00][Mg0.44Fe0.56Mn0.00]2[Si0.32Al0.68]2[(OH)0.80O0.20]2 +g: + composition: Ca0.05Mg0.4Mn0.25Fe2.3Al2.0Si3.0O12.0 + site occupancies: [Mg0.13Fe0.77Mn0.08Ca0.02]3[Al1.00Fe3+0.00]2 +ilmm: + free site occupancy vector: [Fe-0.50Ti0.50][Ti-0.50Fe0.50] + composition: Mg0.05Ti0.9Mn0.05FeO3.0 + site occupancies: [Fe0.37Ti0.43Fe3+0.10Mg0.05Mn0.05][Ti0.47Fe0.43Fe3+0.10] +st: + composition: Mg0.72Ti0.12Mn0.05Fe3.3Al17.77Si7.5O48.0H4.0 + site occupancies: [Mg0.18Fe0.81Mn0.01]4[Al0.89Fe3+0.03Ti0.06V0.02]2 +q: + composition: SiO2.0 + diff --git a/tests/test_solution.py b/tests/test_solution.py index 1bb804fc..4824435c 100644 --- a/tests/test_solution.py +++ b/tests/test_solution.py @@ -521,6 +521,36 @@ def test_ol_ss_formulae(self): self.assertEqual(ol_ss.site_formula(), "[Mg0.40Fe0.60]2SiO4") + def test_ol_ss_formulae_high_precision(self): + P = 1.0e5 + T = 1000.0 + ol_ss = olivine_ss() + f = 0.43673678 + ol_ss.set_composition([f, 1 - f]) + ol_ss.set_state(P, T) + self.assertEqual(ol_ss.site_formula(precision=4), "[Mg0.4367Fe0.5633]2SiO4") + + def test_ol_ss_formulae_absent_species(self): + P = 1.0e5 + T = 1000.0 + ol_ss = olivine_ss() + ol_ss.set_composition([0.0, 1.0]) + ol_ss.set_state(P, T) + self.assertEqual(ol_ss.site_formula(), "[Mg0.00Fe1.00]2SiO4") + self.assertEqual( + ol_ss.site_formula(print_absent_species=False), "[Fe1.00]2SiO4" + ) + + def test_ol_ss_formulae_custom_site_occupancies(self): + P = 1.0e5 + T = 1000.0 + ol_ss = olivine_ss() + ol_ss.set_composition([0.4, 0.6]) + ol_ss.set_state(P, T) + self.assertEqual( + ol_ss.site_formula(site_occupancies=[1.0, 2.0]), "[Mg1.00Fe2.00]2SiO4" + ) + def test_1_gibbs(self): fo, fo_ss = self.setup_1min_ss() with warnings.catch_warnings(record=True) as w: diff --git a/tests/test_tools.py b/tests/test_tools.py index 53441dd3..2fafe0f3 100644 --- a/tests/test_tools.py +++ b/tests/test_tools.py @@ -466,6 +466,25 @@ def test_estimate_conditions_fail_fixed_P_and_T(self): temperature_bounds=[500.0, 500.0], ) + def test_estimate_conditions_free_compositional_vector(self): + assemblage = setup_assemblage() + theta = np.array([0.5e9, 500.0]) + dataset_covariances = HP_2011_ds62.cov() + + # set free compositional vectors for each phase + # corresponding to isochemical reactions + for phase in assemblage.phases: + r = Composite([phase]).reaction_basis + if r.shape[0] > 0: + phase.free_compositional_vectors = Composite([phase]).reaction_basis + + res = estimate_conditions( + assemblage, dataset_covariances, guessed_conditions=theta + ) + self.assertTrue(res.success) + self.assertEqual(res.n_reactions, 18) + self.assertEqual(res.n_params, 4) + if __name__ == "__main__": unittest.main() diff --git a/tests/test_utils.py b/tests/test_utils.py index aaae46a1..158f49a3 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -7,6 +7,7 @@ from sympy import Matrix, Rational import scipy.integrate as integrate +from burnman.utils.chemistry import formula_to_string from burnman.utils.misc import extract_lines_between_markers from burnman.utils.misc import run_cli_program_with_input from burnman.utils.math import smooth_array @@ -300,6 +301,26 @@ def test_visual_center_area_zero(self): center = visual_center_of_polygon(tiny) self.assertTrue(np.allclose(center, [0.005, 0.05], atol=1e-4)) + def test_formula_to_string(self): + formula = {"Si": 1, "O": 4, "Mg": 1.5, "Fe": 0.5} + result = formula_to_string(formula, use_fractions=False) + self.assertEqual(result, "Mg1.5Fe0.5SiO4") + + result_frac = formula_to_string(formula) + self.assertEqual(result_frac, "Mg3/2Fe1/2SiO4") + + formula_zero = {"Si": 1.0, "O": 4.0, "Mg": 0.0, "Fe": 2.0} + result_zero = formula_to_string( + formula_zero, use_fractions=False, significant_figures=1 + ) + self.assertEqual(result_zero, "Fe2.0SiO4.0") + + formula_sf = {"Si": 2, "O": 4, "Mg": 1.123456, "Fe": 0.5} + result_sf = formula_to_string( + formula_sf, use_fractions=False, significant_figures=2 + ) + self.assertEqual(result_sf, "Mg1.12Fe0.5Si2O4") + if __name__ == "__main__": unittest.main()