diff --git a/.gitignore b/.gitignore index 9644749fc..bff73a585 100644 --- a/.gitignore +++ b/.gitignore @@ -120,3 +120,4 @@ ENV/ _autosummary uv.lock +JANAF_*_data.json diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index d7983abcf..a2ffbc9d0 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -39,4 +39,4 @@ repos: - id: codespell stages: [pre-commit, commit-msg] exclude_types: [json, bib, svg] - args: [--ignore-words-list, "mater,fwe,te"] + args: [--ignore-words-list, "mater,fwe,te,alo"] diff --git a/mp_api/client/core/_oxygen_evolution.py b/mp_api/client/core/_oxygen_evolution.py new file mode 100644 index 000000000..57121de67 --- /dev/null +++ b/mp_api/client/core/_oxygen_evolution.py @@ -0,0 +1,348 @@ +"""Compute the oxygen evolution of a phase.""" +from __future__ import annotations + +import json +from collections.abc import Sequence +from io import StringIO +from pathlib import Path +from warnings import warn + +import numpy as np +import pandas as pd +import requests +from pymatgen.analysis.phase_diagram import PhaseDiagram +from pymatgen.core import Composition, Element +from scipy.constants import Avogadro, Boltzmann, atm, elementary_charge +from scipy.interpolate import make_splrep, splev + +DEFAULT_CACHE_FILE = Path(__file__).absolute().parent / "JANAF_O2_data.json" + + +class OxygenEvolution: + """Compute the oxygen evolution using experimental data and atomistic phase diagrams. + + Parameters + ----------- + cache_file : Path or None + If a Path, the path to the file at which NIST JANAF data is cached / retrieved. + If None, no caching is performed. + """ + + def __init__( + self, + cache_file: Path | None = DEFAULT_CACHE_FILE, + ): + self._spline_pars = None + self.cache_file = cache_file + + def get_chempot_temp_data( + self, + nist_url: str = "https://janaf.nist.gov/tables/O-029.txt", + ref_p: float = 0.21 + * atm + * 1e-6, # O2 partial pressure at ambient conditions, in MPa + meas_p: float = 0.1, # The reference pressure reported in JANAF, 0.1 MPa + ) -> tuple[np.ndarray, np.ndarray]: + """Get the approximate relationship between the O2 chemical potential and temperature. + + The user is advised not to change the default parameters below unless they + know what they control. + + Parameters + ----------- + nist_url : str + The link to the txt / TSV data on NIST's JANAF tables for gaseous O2 + ref_p : float + The reference pressure for O2 at ambient conditions. + The default value is 0.21 atmospheres, which we convert to MPa + meas_p : float + The O2 partial pressure at which the experimental data was collected. + Per NIST, this is 0.1 MPa + + Returns: + ----------- + tuple of + Temperature in K (np.ndarray) + O2 chemical potential in eV/atom (np.ndarray) + """ + if self.cache_file and self.cache_file.exists(): + data = json.loads(self.cache_file.read_text()) + return tuple(np.array(data[k]) for k in ("temperature", "mu-mu_0K")) + + response = requests.get(nist_url) + + with StringIO(initial_value=response.content.decode()) as stio: + janaf_data = pd.read_csv(stio, sep="\t", header=1) + + kB_eV = Boltzmann / elementary_charge + + temp_K = janaf_data["T(K)"].to_numpy() + entropy_dimless = janaf_data["S"].to_numpy() / (Boltzmann * Avogadro) + + # Note that the ideal gas contribution to the chemical potential is + # (Gibbs-Helmholtz relation) + # G(p, T)/N = G(p_0, T_0)/N + k_B T ln(p/p_0) + # where G(p_0, T_0)/N is a reference chemical potential for the system + # We use the zero temperature energy per atom to approximate this later + mu_mu0 = kB_eV * temp_K * (1.0 + np.log(meas_p / ref_p) - entropy_dimless) + if self.cache_file: + self.cache_file.write_text( + json.dumps( + {"mu-mu_0K": mu_mu0.tolist(), "temperature": temp_K.tolist()}, + ) + ) + + return temp_K, mu_mu0 + + def mu_to_temp_spline_params(self): + if self._spline_pars is None: + temp = np.array(NIST_JANAF_O2_MU_T["temperature"]) + mu_m0 = np.array(NIST_JANAF_O2_MU_T["mu-mu_0K"]) + # spline needs the x vars to be in increasing order + mu_sort_idx = np.argsort(mu_m0) + self._spline_pars = make_splrep(mu_m0[mu_sort_idx], temp[mu_sort_idx]) + return self._spline_pars + + def mu_to_temp_spline( + self, + mu: Sequence[float], + ): + """Interpolate the O2 chemical potential vs. T relationship. + + Parameters + ----------- + mu : np.ndarray + The O2 chemical potentials to sample at. + """ + mu_arr = np.array(mu) + if np.any( + (mu_arr < min(NIST_JANAF_O2_MU_T["mu-mu_0K"])) + | (mu_arr > max(NIST_JANAF_O2_MU_T["mu-mu_0K"])) + ): + warn( + "Some of the input chemical potential values are " + "outside the fitting range - extrapolation will be inaccurate.", + stacklevel=2, + ) + return splev(mu_arr, self.mu_to_temp_spline_params()) + + @staticmethod + def stairstep(x, y, z): + x = np.array(x) + y = np.array(y) + z = np.array(z) + isort = np.argsort(x) + + new_x = np.zeros(2 * x.shape[0] - 1) + new_y = np.zeros(2 * y.shape[0] - 1) + new_z = np.zeros(2 * z.shape[0] - 1, dtype=object) + + new_x[::2] = x[isort] + new_x[1::2] = x[isort][1:] + + new_y[::2] = y[isort] + new_y[1::2] = y[isort][:-1] + + new_z[::2] = z[isort] + new_z[1::2] = z[isort][:-1] + + return new_x, new_y, new_z + + def get_oxygen_evolution_from_phase_diagram( + self, + phase_diagram: PhaseDiagram, + compositions: list[Composition] | set[Composition], + ) -> dict[str, np.ndarray]: + by_dict = { + composition.formula: [ + { + k: profile[v] + for k, v in { + "mu": "chempot", + "reaction": "reaction", + "evolution": "evolution", + }.items() + } + for profile in phase_diagram.get_element_profile("O", composition) + ] + for composition in set(compositions) + } + + for formula, data in by_dict.items(): + for idx, entry in enumerate(data): + # Normalize all reactions to have integer coefficients + scale = entry["reaction"].normalized_repr_and_factor()[1] + by_dict[formula][idx]["reaction"]._coeffs = [ + c * scale for c in entry["reaction"]._coeffs + ] + + oxy_evo_data = { + formula: { + k: [data[idx][k] for idx in range(len(data))] + for k in ( + "mu", + "reaction", + "evolution", + ) + } + for formula, data in by_dict.items() + } + + mu_0K = phase_diagram.el_refs[Element.O].energy_per_atom + + for formula, data in oxy_evo_data.items(): + ( + oxy_evo_data[formula]["mu"], + oxy_evo_data[formula]["evolution"], + oxy_evo_data[formula]["reaction"], + ) = self.stairstep(*(data[k] for k in ("mu", "evolution", "reaction"))) + oxy_evo_data[formula][ + "evolution" + ] *= -0.5 # This is the normalization convention we adopt for MP + oxy_evo_data[formula]["temperature"] = self.mu_to_temp_spline( + data["mu"] - mu_0K + ) + return oxy_evo_data + + +# This data is generated by running: +# `OxygenEvolution().get_chempot_temp_data()` +# The contents of `DEFAULT_CACHE_FILE` will contain this data. +# Because the last update of the JANAF tables was in 1998 +# (writing in 2025), we see no need to dynamically retrieve this data +NIST_JANAF_O2_MU_T = { + "mu-mu_0K": [ + 0.0, + -0.15766752295964667, + -0.35716109254410505, + -0.4633062592086276, + -0.5684747265209704, + -0.5725679617587883, + -0.6845046823072302, + -0.7988365855732343, + -0.9153600136580551, + -1.0338982557640533, + -1.2764843857253183, + -1.5256155519745132, + -1.7804988878829335, + -2.0404658980577506, + -2.304982822611666, + -2.5736050343764276, + -2.8459749660488933, + -3.1217816895393757, + -3.4007764623761245, + -3.6827416348963578, + -3.9674927231001944, + -4.254892918628172, + -4.544805413120827, + -4.837112053904078, + -5.13170297971957, + -5.428488021421295, + -5.7273687184475115, + -6.028299468011739, + -6.3311932102488635, + -6.635957703158943, + -6.9425981288768055, + -7.2509807883238855, + -7.561092207949628, + -7.872879529978791, + -8.186271240950743, + -8.501236248056516, + -8.817722729947826, + -9.135721358781987, + -9.4551295282894, + -9.775953457031855, + -10.098137177953216, + -10.421667217502918, + -10.746499009321454, + -11.072592132757157, + -11.399910312866226, + -11.728421420412731, + -12.058141001801161, + -12.389003762132672, + -12.720990009294923, + -13.05408419688343, + -13.388322599841986, + -13.723603648327966, + -14.059921123779567, + -14.397323738264184, + -14.735764852568353, + -15.075250685253877, + -15.415683812185982, + -15.757182386038748, + -16.099653128385278, + -16.44311262205701, + -16.787639635503336, + -17.13308492324121, + -17.47953139942788, + -17.826947971254377, + -18.175365731529666, + ], + "temperature": [ + 0.0, + 100.0, + 200.0, + 250.0, + 298.15, + 300.0, + 350.0, + 400.0, + 450.0, + 500.0, + 600.0, + 700.0, + 800.0, + 900.0, + 1000.0, + 1100.0, + 1200.0, + 1300.0, + 1400.0, + 1500.0, + 1600.0, + 1700.0, + 1800.0, + 1900.0, + 2000.0, + 2100.0, + 2200.0, + 2300.0, + 2400.0, + 2500.0, + 2600.0, + 2700.0, + 2800.0, + 2900.0, + 3000.0, + 3100.0, + 3200.0, + 3300.0, + 3400.0, + 3500.0, + 3600.0, + 3700.0, + 3800.0, + 3900.0, + 4000.0, + 4100.0, + 4200.0, + 4300.0, + 4400.0, + 4500.0, + 4600.0, + 4700.0, + 4800.0, + 4900.0, + 5000.0, + 5100.0, + 5200.0, + 5300.0, + 5400.0, + 5500.0, + 5600.0, + 5700.0, + 5800.0, + 5900.0, + 6000.0, + ], +} diff --git a/mp_api/client/core/settings.py b/mp_api/client/core/settings.py index 674552044..200b67785 100644 --- a/mp_api/client/core/settings.py +++ b/mp_api/client/core/settings.py @@ -3,7 +3,7 @@ from typing import List from pydantic import Field -from pydantic_settings import BaseSettings +from pydantic_settings import BaseSettings, SettingsConfigDict from pymatgen.core import _load_pmg_settings from mp_api.client import __file__ as root_dir @@ -87,5 +87,4 @@ class MAPIClientSettings(BaseSettings): _MAX_LIST_LENGTH, description="Maximum length of query parameter list" ) - class Config: - env_prefix = "MPRESTER_" + model_config = SettingsConfigDict(env_prefix="MPRESTER_") diff --git a/mp_api/client/mprester.py b/mp_api/client/mprester.py index 50fd53261..c3fffcbb3 100644 --- a/mp_api/client/mprester.py +++ b/mp_api/client/mprester.py @@ -22,6 +22,7 @@ from requests import Session, get from mp_api.client.core import BaseRester, MPRestError +from mp_api.client.core._oxygen_evolution import OxygenEvolution from mp_api.client.core.settings import MAPIClientSettings from mp_api.client.core.utils import _compare_emmet_ver, load_json, validate_ids from mp_api.client.routes import GeneralStoreRester, MessagesRester, UserSettingsRester @@ -1741,41 +1742,7 @@ def get_oxygen_evolution( inserted_chemsys, thermo_type=ThermoType(thermo_type), ) - - by_dict = { - composition.formula: [ - { - k: profile[v] - for k, v in { - "mu": "chempot", - "reaction": "reaction", - "evolution": "evolution", - }.items() - } - for profile in phase_diagram.get_element_profile("O", composition) - ] - for composition in unique_composition - } - - target_comp = Composition({"O": 1}) - for formula, data in by_dict.items(): - for idx, entry in enumerate(data): - # Normalize all reactions to have integer coefficients - scale = entry["reaction"].normalized_repr_and_factor()[1] - by_dict[formula][idx]["reaction"]._coeffs = [ - c * scale for c in entry["reaction"]._coeffs - ] - - by_dict[formula][idx]["O2_produced"] = ( - entry["reaction"].get_coeff(target_comp) - if target_comp in entry["reaction"].products - else 0 - ) - - return { - formula: { - k: [data[idx][k] for idx in range(len(data))] - for k in ("mu", "reaction", "evolution", "O2_produced") - } - for formula, data in by_dict.items() - } + return OxygenEvolution().get_oxygen_evolution_from_phase_diagram( + phase_diagram, + unique_composition, + ) diff --git a/tests/core/test_oxygen_evolution.py b/tests/core/test_oxygen_evolution.py new file mode 100644 index 000000000..a83247bf8 --- /dev/null +++ b/tests/core/test_oxygen_evolution.py @@ -0,0 +1,68 @@ +"""Test the oxygen evolution analysis features in the client.""" + +import pytest +import numpy as np +from scipy.interpolate import splev + +from pymatgen.analysis.phase_diagram import PhaseDiagram +from pymatgen.core import Composition +from pymatgen.entries.computed_entries import ComputedEntry + +from mp_api.client.core._oxygen_evolution import NIST_JANAF_O2_MU_T, OxygenEvolution + + +def test_interp(): + oxyevo = OxygenEvolution() + + # Spline interpolation is exact for all fit points: + assert np.allclose( + oxyevo.mu_to_temp_spline(np.array(NIST_JANAF_O2_MU_T["mu-mu_0K"])), + NIST_JANAF_O2_MU_T["temperature"], + ) + + # mu(T) is a generally decreasing function, at least for fit range + deriv = splev( + np.linspace( + min(NIST_JANAF_O2_MU_T["mu-mu_0K"]), + max(NIST_JANAF_O2_MU_T["mu-mu_0K"]), + 1000, + ), + oxyevo.mu_to_temp_spline_params(), + der=1, + ) + assert len(deriv[deriv > 0]) == 0 + + # ensure user is warned when data is outside fit range + for badval in (1.0, -50.0): + with pytest.warns(UserWarning, match="outside the fitting range"): + oxyevo.mu_to_temp_spline(badval) + + +def test_oxy_evo(): + """Very fake Al-O chemical system for testing purposes.""" + entries = [ + ComputedEntry({"Al": 1}, -2.0), + ComputedEntry({"O": 1}, -5.0), + ComputedEntry({"Al": 1, "O": 1}, -25.0), + ComputedEntry({"Al": 2, "O": 3}, -55.0), + ComputedEntry({"Al": 3, "O": 2}, -5.0), + ] + phase_diag = PhaseDiagram(entries) + + ref_comp = Composition({"Al": 4, "O": 1}) + evo = OxygenEvolution().get_oxygen_evolution_from_phase_diagram( + phase_diag, [ref_comp] + ) + assert ref_comp.formula in evo + assert np.allclose(evo[ref_comp.formula]["evolution"], [0.5, 0.5, -1.5]) + assert np.allclose(evo[ref_comp.formula]["temperature"], [5949.70500942, 0.0, 0.0]) + assert all( + str(evo[ref_comp.formula]["reaction"][idx]) == rxn + for idx, rxn in enumerate( + [ + "2 Al4O -> 8 Al + O2", + "2 Al4O -> 8 Al + O2", + "2 Al4O + 3 O2 -> 8 AlO", + ] + ) + )