Source code for rxn_network.entries.gibbs

"""A computed entry object for estimating the Gibbs free energy of formation. Note that
this is similar to the implementation within pymatgen, but has been refactored here to
add extra functionality.
"""

from __future__ import annotations

from copy import deepcopy
from itertools import combinations
from typing import TYPE_CHECKING

import numpy as np
from monty.json import MontyDecoder
from pymatgen.analysis.phase_diagram import GrandPotPDEntry
from pymatgen.entries.computed_entries import ComputedEntry, ConstantEnergyAdjustment
from scipy.interpolate import interp1d

from rxn_network.core import Composition
from rxn_network.data import G_ELEMS

if TYPE_CHECKING:
    from pymatgen.core.periodic_table import Element
    from pymatgen.core.structure import Structure
    from pymatgen.entries.computed_entries import EnergyAdjustment


[docs] class GibbsComputedEntry(ComputedEntry): """An extension to ComputedEntry which estimates the Gibbs free energy of formation of solids using energy adjustments from the machine-learned SISSO descriptor from Bartel et al. (2018). Note that this is similar to the implementation within pymatgen, but has been refactored here to add extra functionality. WARNING: This descriptor only applies to solids. See entries.nist.NISTReferenceEntry for common gases (e.g. CO2). If you use this entry class in your work, please consider citing the following paper: Bartel, C. J., Millican, S. L., Deml, A. M., Rumptz, J. R., Tumas, W., Weimer, A. W., … Holder, A. M. (2018). Physical descriptor for the Gibbs energy of inorganic crystalline solids and temperature-dependent materials chemistry. Nature Communications, 9(1), 4168. https://doi.org/10.1038/s41467-018-06682-4. """ def __init__( self, composition: Composition, formation_energy_per_atom: float, volume_per_atom: float, temperature: float, energy_adjustments: list[EnergyAdjustment] | None = None, parameters: dict | None = None, data: dict | None = None, entry_id: object | None = None, ): """A new computed entry object is returned with a supplied energy correction representing the difference between the formation enthalpy at T=0K and the Gibbs formation energy at the specified temperature. Args: composition: The composition object (pymatgen) formation_energy_per_atom: Calculated formation enthalpy, dH, at T = 298 K, normalized to the total number of atoms in the composition. NOTE: since this is a _formation_ energy, it must be calculated using a phase diagram construction. volume_per_atom: The total volume of the associated structure divided by its number of atoms. temperature: Temperature [K] by which to acquire dGf(T); must be selected from a range of [300, 2000] K. If temperature is not selected from one of [300, 400, 500, ... 2000 K], then free energies will be interpolated. energy_adjustments: Optional list of energy adjustments. parameters: Optional list of calculation parameters. data: Optional dictionary containing entry data. entry_id: An identifying string for the entry, such as the entry's mpid (e.g., "mp-25025"). While optional, this is recommended to set as several downstream classes depend on its use. If an entry_id is not provided, a combination of the composition and volume will be used (e.g., "Li2O_8.4266"). """ composition = Composition(composition) self._composition = composition self.volume_per_atom = volume_per_atom num_atoms = composition.num_atoms if temperature < 300 or temperature > 2000: raise ValueError("Temperature must be selected from range: [300, 2000] K.") if energy_adjustments is not None: energy_adjustments = [ adjustment for adjustment in energy_adjustments if adjustment.name != "Gibbs SISSO Correction" ] else: energy_adjustments = [] energy_adjustments.append( ConstantEnergyAdjustment( self.gibbs_adjustment(temperature), uncertainty=0.05 * num_atoms, # descriptor has ~50 meV/atom MAD name="Gibbs SISSO Correction", description=f"Gibbs correction: dGf({temperature} K) - dHf (298 K)", ) ) formation_energy = num_atoms * formation_energy_per_atom if entry_id is None: # should set an entry_id for downstream processing entry_id = f"{composition.reduced_formula}_{volume_per_atom:.4f}" super().__init__( composition=composition, energy=formation_energy, energy_adjustments=energy_adjustments, parameters=parameters, data=data, entry_id=entry_id, ) self._composition = composition self.formation_energy_per_atom = formation_energy_per_atom self.temperature = temperature
[docs] def get_new_temperature(self, new_temperature: float) -> GibbsComputedEntry: """Return a copy of the GibbsComputedEntry at the new specified temperature. Args: new_temperature: The new temperature to use [K] Returns: A copy of the GibbsComputedEntry, initialized at the new specified temperature. """ new_entry_dict = self.as_dict() new_entry_dict["temperature"] = new_temperature return self.from_dict(new_entry_dict)
[docs] def gibbs_adjustment(self, temperature: float) -> float: """Returns the difference between the predicted Gibbs formation energy and the formation enthalpy at 298 K, i.e., dGf(T) - dHf(298 K). Calculated using SISSO descriptor from Bartel et al. (2018) and elemental chemical potentials (FactSage). Units: eV (not normalized) Args: temperature: The absolute temperature [K]. Returns: The correction to Gibbs free energy of formation (eV) from DFT energy. """ if self._composition.is_element: return 0 num_atoms = self._composition.num_atoms reduced_mass = self._reduced_mass(self._composition) return num_atoms * self._g_delta_sisso(self.volume_per_atom, reduced_mass, temperature) - self._sum_g_i( self._composition, temperature )
[docs] def to_grand_entry(self, chempots: dict[Element, float]) -> GrandPotPDEntry: """Convert a GibbsComputedEntry to a GrandComputedEntry. Args: chempots: A dictionary of {element: chempot} pairs. Returns: A GrandComputedEntry. """ return GrandPotPDEntry(self, chempots)
[docs] def copy(self) -> GibbsComputedEntry: """Returns a deepcopy of the GibbsComputedEntry object.""" return deepcopy(self)
@staticmethod def _g_delta_sisso(volume_per_atom: float, reduced_mass: float, temp: float) -> float: """G^delta as predicted by SISSO-learned descriptor from Eq. (4) in Bartel et al. (2018). Args: volume_per_atom: volume per atom [Å^3/atom] reduced_mass: reduced mass as calculated with pair-wise sum formula [amu] temp: Temperature [K] Returns: float: G^delta """ return ( (-2.48e-4 * np.log(volume_per_atom) - 8.94e-5 * reduced_mass / volume_per_atom) * temp + 0.181 * np.log(temp) - 0.882 ) @staticmethod def _sum_g_i(composition: Composition, temperature: float) -> float: """Sum of the stoichiometrically weighted chemical potentials [eV] of the elements at specified temperature, as acquired from data/elements.json. """ elems = composition.get_el_amt_dict() if temperature % 100 > 0: sum_g_i = 0 for elem, amt in elems.items(): g_interp = interp1d( [float(t) for t in G_ELEMS], [g_dict[elem] for g_dict in G_ELEMS.values()], ) sum_g_i += amt * g_interp(temperature) else: sum_g_i = sum(amt * G_ELEMS[str(temperature)][elem] for elem, amt in elems.items()) return sum_g_i @staticmethod def _reduced_mass(composition: Composition) -> float: """Reduced mass [amu] as calculated via Eq. 6 in Bartel et al. (2018), to be used in SISSO descriptor equation. """ reduced_comp = composition.reduced_composition num_elems = len(reduced_comp.elements) elem_dict = reduced_comp.get_el_amt_dict() denominator = (num_elems - 1) * reduced_comp.num_atoms all_pairs = combinations(elem_dict.items(), 2) mass_sum = 0.0 for pair in all_pairs: m_i = Composition(pair[0][0]).weight m_j = Composition(pair[1][0]).weight alpha_i = pair[0][1] alpha_j = pair[1][1] mass_sum += (alpha_i + alpha_j) * (m_i * m_j) / (m_i + m_j) return (1 / denominator) * mass_sum
[docs] @classmethod def from_structure( cls, structure: Structure, formation_energy_per_atom: float, temperature: float, **kwargs, ) -> GibbsComputedEntry: """Constructor method for building a GibbsComputedEntry from a structure, formation enthalpy, and temperature. Args: structure: Structure object (pymatgen) formation_energy_per_atom: Formation enthalpy at T = 298 K associated with structure temperature: Desired temperature [K] for acquiring dGf(T) **kwargs: Optional kwargs to be passed to GibbsComputedEntry.__init__, such as entry_id, energy_adjustments, parameters, and data. Returns: A new GibbsComputedEntry object """ composition = Composition(structure.composition).element_composition volume_per_atom = structure.volume / structure.num_sites return cls( composition=composition, formation_energy_per_atom=formation_energy_per_atom, volume_per_atom=volume_per_atom, temperature=temperature, **kwargs, )
@property def is_experimental(self) -> bool: """Returns True if self.data contains {"theoretical": False}. If theoretical is not specified but there is greater than 1 icsd_id provided, assumes that the presence of an icsd_id means the entry is experimental. """ if "theoretical" in self.data: return not self.data["theoretical"] if "icsd_ids" in self.data: return len(self.data["icsd_ids"]) >= 1 return False @property def unique_id(self) -> str: """Returns a unique ID for the entry based on its entry-id and temperature. This is useful because the same entry-id can be used for multiple entries at different temperatures. """ return f"{self.entry_id}_{self.temperature}"
[docs] def as_dict(self) -> dict: """Returns an MSONable dict.""" data = super().as_dict() data["volume_per_atom"] = self.volume_per_atom data["formation_energy_per_atom"] = self.formation_energy_per_atom data["temperature"] = self.temperature return data
[docs] @classmethod def from_dict(cls, d: dict) -> GibbsComputedEntry: """Returns a GibbsComputedEntry object from MSONable dictionary.""" dec = MontyDecoder() return cls( composition=d["composition"], formation_energy_per_atom=d["formation_energy_per_atom"], volume_per_atom=d["volume_per_atom"], temperature=d["temperature"], energy_adjustments=dec.process_decoded(d["energy_adjustments"]), parameters=d["parameters"], data=d["data"], entry_id=d["entry_id"], )
def __repr__(self): output = [ ( f"GibbsComputedEntry | {self.entry_id} | {self.composition.formula} " f"({self.composition.reduced_formula})" ), f"Gibbs Energy ({self.temperature} K) = {self.energy:.4f}", ] return "\n".join(output) def __eq__(self, other): if type(other) is not type(self): return False if not np.isclose(self.temperature, other.temperature): return False if not np.isclose(self.energy, other.energy): return False if getattr(self, "entry_id", None) and getattr(other, "entry_id", None): return self.entry_id == other.entry_id if self.composition != other.composition: return False return True def __hash__(self): return hash( ( self.composition, self.energy, self.entry_id, self.temperature, ) )