"""This module implements added features to the ChemicalPotentialDiagram class from
pymatgen.
"""
from __future__ import annotations
from functools import cached_property
from typing import TYPE_CHECKING
import numpy as np
from pymatgen.analysis.chempot_diagram import ChemicalPotentialDiagram as ChempotDiagram
from pymatgen.analysis.phase_diagram import PhaseDiagram
from scipy.spatial import HalfspaceIntersection, KDTree
from rxn_network.entries.entry_set import GibbsEntrySet
if TYPE_CHECKING:
from pymatgen.analysis.phase_diagram import PDEntry
from pymatgen.core.periodic_table import Element
[docs]
class ChemicalPotentialDiagram(ChempotDiagram):
"""This class is an extension of the ChemicalPotentialDiagram class from pymatgen.
Several features have been added to the original class for the purpose of
efficiently calculating the shortest distance between two chemical potential
domains.
For more information on this specific implementation of the algorithm, please
cite/reference the paper below:
Todd, P. K., McDermott, M. J., Rom, C. L., Corrao, A. A., Denney, J. J.,
Dwaraknath, S. S., Khalifah, P. G., Persson, K. A., & Neilson, J. R. (2021).
Selectivity in Yttrium Manganese Oxide Synthesis via Local Chemical Potentials
in Hyperdimensional Phase Space. Journal of the American Chemical Society,
143(37), 15185-15194. https://doi.org/10.1021/jacs.1c06229
"""
def __init__( # pylint: disable=super-init-not-called
self,
entries: list[PDEntry],
limits: dict[Element, float] | None = None,
default_min_limit: float | None = -100.0,
):
"""Initialize a ChemicalPotentialDiagram object.
Args:
entries: List of PDEntry-like objects containing a composition and
energy. Must contain elemental references and be suitable for typical
phase diagram construction. Entries must be within a chemical system
with 2 or more elements.
limits: Bounds of elemental chemical potentials (min, max), which are
used to construct the border hyperplanes used in the
HalfspaceIntersection algorithm; these constrain the space over which
the domains are calculated and also determine the size of the plotted
diagram. Any elemental limits not specified are covered in the
default_min_limit argument.
default_min_limit (float): Default minimum chemical potential limit for
unspecified elements within the "limits" argument. This results in
default limits of (-100, 0).
"""
self.entries = sorted(entries, key=lambda e: e.composition.reduced_formula)
self._entry_set = GibbsEntrySet(self.entries)
self.limits = limits
self.default_min_limit = default_min_limit
self.elements = sorted({els for e in self.entries for els in e.composition.elements})
self.dim = len(self.elements)
self._min_entries, self._el_refs = self._get_min_entries_and_el_refs(self.entries)
self._entry_dict = {e.composition.reduced_formula: e for e in self._min_entries}
self._border_hyperplanes = self._get_border_hyperplanes()
(
self._hyperplanes,
self._hyperplane_entries,
) = self._get_hyperplanes_and_entries()
if self.dim < 2:
raise ValueError("ChemicalPotentialDiagram currently requires phase diagrams with 2 or more elements!")
if len(self.el_refs) != self.dim:
missing = set(self.elements).difference(self.el_refs.keys())
raise ValueError(f"There are no entries for the terminal elements: {missing}")
self._hs_int = self._get_halfspace_intersection()
num_hyperplanes = len(self._hyperplanes)
num_border_hyperplanes = len(self._border_hyperplanes)
self._border_hyperplane_indices = list(range(num_hyperplanes, num_hyperplanes + num_border_hyperplanes))
self._metastable_domains: dict[str, list] = {} # for caching
[docs]
def shortest_domain_distance(self, f1: str, f2: str, offset: float = 0.0) -> float:
"""Returns the chemical potential distance between two phase domains. Also works
for metastable phases (see metastable_domains property).
Args:
f1: chemical formula of phase 1
f2: chemical formula of phase 2
offset: an optional offset (eV/atom) to add to the calculated distance. See
get_offset() method.
Returns:
Shortest distance between domain boundaries in the full (hyper)dimensional
space, calculated using KDTree.
"""
pts1 = self.domains[f1] if f1 in self.domains else self._get_metastable_domain(f1)
pts2 = self.domains[f2] if f2 in self.domains else self._get_metastable_domain(f2)
tree = KDTree(pts1)
return min(tree.query(pts2)[0]) + offset
[docs]
def get_offset(self, entry: PDEntry) -> float:
"""For a given entry, returns the distance between its hyperplane and the surface
of the chemical potential diagram. This allows one to represent the energy above
hull in chemical potential space. Returns zero for stable entries.
Args:
entry: A stable or metastable entry within the chemical potential diagram
Returns:
Offset in chemical potential distance (eV/atom)
"""
if entry in self._min_entries and entry.composition.reduced_formula in self.domains:
offset = 0.0
else:
e_above_hull = self._entry_set.get_e_above_hull(entry)
hyperplane = self._get_hyperplane(entry)
offset = self._get_distance_between_parallel_hyperplanes(hyperplane[:-1], e_above_hull)
return offset
@cached_property
def domains(self) -> dict[str, np.ndarray]:
"""Mapping of formulas to array of domain boundary points. Cached for quicker
calculations.
"""
return self._get_domains()
@property
def metastable_domains(self) -> dict[str, np.ndarray]:
"""Gets a dictionary of the chemical potential domains for metastable chemical
formulas. This corresponds to the domains of the relevant phases if they were
just barely thermodynamically stable (on the hull).
"""
return {
e.composition.reduced_formula: self._get_metastable_domain(e.composition.reduced_formula)
for e in self._min_entries
if e.composition.reduced_formula not in self.domains
}
@property
def hs_int(self) -> HalfspaceIntersection:
"""Returns the scipy HalfSpaceIntersection object used to calculate all domains."""
return self._hs_int
def _get_halfspace_intersection(self):
hs_hyperplanes = np.vstack([self._hyperplanes, self._border_hyperplanes])
interior_point = np.min(self.lims, axis=1) + 1e-1
return HalfspaceIntersection(hs_hyperplanes, interior_point)
def _get_domains(self) -> dict[str, np.ndarray]:
"""Returns a dictionary of chemical potential domains as {formula:
np.ndarray}.
"""
domains: dict[str, list] = {entry.composition.reduced_formula: [] for entry in self._hyperplane_entries}
entries = self._hyperplane_entries
for intersection, facet in zip(self.hs_int.intersections, self.hs_int.dual_facets):
for v in facet:
if v not in self._border_hyperplane_indices:
this_entry = entries[v]
formula = this_entry.composition.reduced_formula
domains[formula].append(intersection)
return {k: np.array(v) for k, v in domains.items() if v}
def _get_hyperplanes_and_entries(self) -> tuple[np.ndarray, list[PDEntry]]:
"""Returns both the array of hyperplanes, as well as a list of the minimum
entries.
"""
data = np.array([self._get_hyperplane(e) for e in self._min_entries])
vec = [self.el_refs[el].energy_per_atom for el in self.elements] + [1]
form_e = -np.dot(data, vec)
inds = np.where(form_e < -PhaseDiagram.formation_energy_tol)[0].tolist()
inds.extend([self._min_entries.index(el) for el in self.el_refs.values()])
hyperplanes = data[inds]
hyperplane_entries = [self._min_entries[i] for i in inds]
return hyperplanes, hyperplane_entries
def _get_hyperplane(self, entry):
return np.array([entry.composition.get_atomic_fraction(el) for el in self.elements] + [-entry.energy_per_atom])
def _get_metastable_domain(self, formula, tol=1e-5):
"""Returns the metastable domain for a given formula. Tol is passed to
GibbsEntrySet.get_stabilized_entry and will affect the size of the domain.
"""
if formula in self._metastable_domains:
return self._metastable_domains[formula]
orig_entry = self._entry_set.get_min_entry_by_formula(formula)
new_entry = self._entry_set.get_stabilized_entry(orig_entry, tol=tol)
self._entry_set.add(new_entry)
cpd = ChemicalPotentialDiagram(self._entry_set, default_min_limit=-500)
try:
metastable_domain = cpd.domains[formula]
except KeyError as exc:
# sometimes if the entry is exactly on the hull it fails, so set force=True
# and make bigger tolerance
self._entry_set.remove(new_entry)
new_entry = self._entry_set.get_stabilized_entry(orig_entry, tol=1e-2, force=True)
self._entry_set.add(new_entry)
cpd = ChemicalPotentialDiagram(self._entry_set, default_min_limit=-500)
try:
metastable_domain = cpd.domains[formula]
except KeyError:
raise ValueError(
f"Failed even after attempted fix. Metastable domain for {formula} can not be created!"
) from exc
self._metastable_domains[formula] = metastable_domain
self._entry_set.remove(new_entry)
return metastable_domain
@staticmethod
def _get_distance_between_parallel_hyperplanes(a, delta_b):
"""Returns the distance between two parallel hyperplanes."""
return np.abs(delta_b) / np.linalg.norm(a)