Source code for rxn_network.reactions.hull

"""Code for analyzing sets of reactions between two phases."""

from __future__ import annotations

from functools import cached_property, lru_cache
from itertools import combinations

import numpy as np
import plotly.express as px
from monty.json import MSONable
from plotly.graph_objs import Figure
from scipy.spatial import ConvexHull

from rxn_network.core import Composition
from rxn_network.reactions.computed import ComputedReaction


[docs] class InterfaceReactionHull(MSONable): """A class for storing and analyzing a set of reactions at an interface between two reactants. This class is more generalized than the InterfacialReactivity class and can encompass any set of reactions between two reactants, regardless of whether the reaction products are "stable" (i.e. together on the convex hull). """ def __init__( self, c1: Composition, c2: Composition, reactions: list[ComputedReaction], ): """ Args: c1: Composition of reactant 1 c2: Composition of reactant 2 reactions: List of reactions containing all enumerated reactions between the two reactants. Note that this list should not include identity reactions of the precursors. """ self.c1 = Composition(c1).reduced_composition self.c2 = Composition(c2).reduced_composition self.e1 = None self.e2 = None for rxn in reactions: for e in rxn.reactant_entries: if e.composition.reduced_composition == self.c1: self.e1 = e elif e.composition.reduced_composition == self.c2: self.e2 = e if self.e1 is not None and self.e2 is not None: break else: raise ValueError( "Provided reactions do not correspond to reactant compositons!", c1, c2, reactions, ) endpoint_reactions = [ ComputedReaction.balance([self.e1], [self.e1]), ComputedReaction.balance([self.e2], [self.e2]), ] reactions_with_endpoints = reactions + endpoint_reactions coords = np.array([(self.get_coordinate(r), r.energy_per_atom) for r in reactions]) coords = np.append(coords, [[0, 0], [1, 0]], axis=0) idx_sort = coords[:, 0].argsort() self.coords = coords[idx_sort] self.reactions = [reactions_with_endpoints[i] for i in idx_sort] self.hull = ConvexHull(self.coords) self.endpoint_reactions = endpoint_reactions
[docs] def plot(self, y_max: float = 0.2) -> Figure: """Plot the reaction hull.""" pts = self._get_scatter() lines = self._get_lines() fig = Figure(data=[*lines, pts]) fig.update_traces( hovertemplate=( "<b>%{hovertext}</b><br> <br><b>Atomic fraction</b>: %{x:.3f}<br><b>Energy</b>: %{y:.3f} (eV/atom)" ) ) fig.update_layout(yaxis_range=[min(self.coords[:, 1]) - 0.01, y_max]) fig.update_layout(xaxis_title="Mixing ratio") fig.update_layout(yaxis_title="Energy (eV/atom)") return fig
[docs] def get_energy_above_hull(self, reaction: ComputedReaction) -> float: """Get the energy of a reaction above the reaction hull.""" idx = self.reactions.index(reaction) x, y = self.coords[idx] return y - self.get_hull_energy(x)
[docs] def get_coordinate(self, reaction: ComputedReaction) -> float: """Get coordinate of reaction in reaction hull. This is expressed as the atomic mixing ratio of component 2 in the reaction. """ amt_c1 = reaction.reactant_atomic_fractions.get(self.c1, 0) amt_c2 = reaction.reactant_atomic_fractions.get(self.c2, 0) total = amt_c1 + amt_c2 # will add to 1.0 with two-component reactions try: coordinate = amt_c2 / total except ZeroDivisionError as e: raise ValueError(f"Can't compute coordinate for {reaction} with {self.c1}, {self.c2}") from e return round(coordinate, 12) # avoids numerical issues
[docs] def get_hull_energy(self, coordinate: float) -> float: """Get the energy of the reaction at a given coordinate. Args: coordinate: Coordinate of reaction in reaction hull. Returns: Energy of reaction at given coordinate. """ reactions = self.get_reactions_by_coordinate(coordinate) return sum(weight * r.energy_per_atom for r, weight in reactions.items())
[docs] def get_reactions_by_coordinate(self, coordinate: float) -> dict[ComputedReaction, float]: """Get the reaction(s) at a given coordinate.""" sorted_vertices = np.sort(self.hull_vertices) for i in range(len(sorted_vertices) - 1): v1 = sorted_vertices[i] v2 = sorted_vertices[i + 1] x1 = self.coords[v1, 0] x2 = self.coords[v2, 0] if np.isclose(coordinate, x1): return {self.reactions[v1]: 1.0} if np.isclose(coordinate, x2): return {self.reactions[v2]: 1.0} if x1 < coordinate < x2: return { self.reactions[v1]: (x2 - coordinate) / (x2 - x1), self.reactions[v2]: (coordinate - x1) / (x2 - x1), } continue raise ValueError("No reactions found!")
[docs] def get_primary_competition(self, reaction: ComputedReaction) -> float: """Calculates the primary competition, C_1, for a reaction (in eV/atom). If you use this selectivity metric in your work, please cite the following work: McDermott, M. J.; McBride, B. C.; Regier, C.; Tran, G. T.; Chen, Y.; Corrao, A. A.; Gallant, M. C.; Kamm, G. E.; Bartel, C. J.; Chapman, K. W.; Khalifah, P. G.; Ceder, G.; Neilson, J. R.; Persson, K. A. Assessing Thermodynamic Selectivity of Solid-State Reactions for the Predictive Synthesis of Inorganic Materials. arXiv August 22, 2023. https://doi.org/10.48550/arXiv.2308.11816. Args: reaction: A computed reaction. Returns: The primary competition (C1 score) in eV/atom. """ energy = reaction.energy_per_atom coord = self.get_coordinate(reaction) matching = np.where(self.coords[:, 0] == coord)[0] # remove all reactions at same coordinate idx_min = matching.min() idx_max = matching.max() competing_rxns = self.reactions[:idx_min] + self.reactions[idx_max + 1 :] # assumes ordered competing_rxns = self.reactions[:idx_min] + self.reactions[idx_max + 1 :] competing_rxn_energies = [r.energy_per_atom for r in competing_rxns] min_energy = min(competing_rxn_energies) return energy - min_energy
[docs] def get_secondary_competition( self, reaction: ComputedReaction, normalize: bool = True, include_e_hull: bool = False, recursive: bool = False, ) -> float: """Calculates the secondary competition, C_2, for a reaction (in eV/atom). If you use this selectivity metric in your work, please cite the following work: McDermott, M. J.; McBride, B. C.; Regier, C.; Tran, G. T.; Chen, Y.; Corrao, A. A.; Gallant, M. C.; Kamm, G. E.; Bartel, C. J.; Chapman, K. W.; Khalifah, P. G.; Ceder, G.; Neilson, J. R.; Persson, K. A. Assessing Thermodynamic Selectivity of Solid-State Reactions for the Predictive Synthesis of Inorganic Materials. arXiv August 22, 2023. https://doi.org/10.48550/arXiv.2308.11816. Args: reaction: A computed reaction. normalize: Whether or not to normalize the sum of secondary reaction sequence energies by the total number of sequnces. Defaults to True according to original definition. include_e_hull: Whether or not to include the energy above hull of the target reaction in the definition of C_2. According to the original paper, this defaults to False. recursive: Whether or not to perform the secondary reaction analysis via a recursive approach. This should return the same answer, but much slower. This argument exists for debugging / validation and defaults to False. Returns: The secondary competition (C2 score) in eV/atom. """ x = self.get_coordinate(reaction) if not recursive: left_energy = self.get_decomposition_energy(0, x) left_num_paths = self.count(len(self.get_coords_in_range(0, x)) - 2) right_energy = self.get_decomposition_energy(x, 1) right_num_paths = self.count(len(self.get_coords_in_range(x, 1)) - 2) else: ( left_energy, left_num_paths, ) = self.get_decomposition_energy_and_num_paths_recursive(0, x) ( right_energy, right_num_paths, ) = self.get_decomposition_energy_and_num_paths_recursive(x, 1) if left_num_paths == 0: left_num_paths = 1 if right_num_paths == 0: right_num_paths = 1 energy = left_energy * right_num_paths + right_energy * left_num_paths total = left_num_paths * right_num_paths if normalize: energy = energy / total if include_e_hull: e_above_hull = self.get_energy_above_hull(reaction) energy -= e_above_hull return -1 * energy
[docs] def get_secondary_competition_max_energy(self, reaction: ComputedReaction) -> float: """Calculates the score for a given reaction. Args: reaction: Reaction to calculate the competition score for. Returns: The competition score for the reaction """ x = self.get_coordinate(reaction) left_energy = self.get_max_decomposition_energy(0, x) right_energy = self.get_max_decomposition_energy(x, 1) return -1 * (left_energy + right_energy)
[docs] def get_secondary_competition_area(self, reaction: ComputedReaction) -> float: """Calculates the score for a given reaction. This formula is based on a methodology presented in the following paper: (TBD). Args: reaction: Reaction to calculate the competition score for. Returns: The competition score for the reaction """ x = self.get_coordinate(reaction) left_area = self.get_decomposition_area(0, x) right_area = self.get_decomposition_area(x, 1) return left_area + right_area
[docs] def get_decomposition_energy(self, x1: float, x2: float) -> float: """Calculates the energy of the reaction decomposition between two points. Args: x1: Coordinate of first point. x2: Coordinate of second point. Returns: The energy of the reaction decomposition between the two points. """ coords = self.get_coords_in_range(x1, x2) n = len(coords) - 2 energy = 0 for c in combinations(range(len(coords)), 3): i_left, i_mid, i_right = sorted(c) c_left, c_mid, c_right = coords[[i_left, i_mid, i_right]] n_left = (i_mid - i_left) - 1 n_right = (i_right - i_mid) - 1 count = self._altitude_multiplicity(n_left, n_right, n) energy += count * self._calculate_altitude(c_left, c_mid, c_right) return energy
[docs] def get_max_decomposition_energy(self, x1: float, x2: float) -> float: """Similar to get_decomposition_energy, but returns only the decomposition energy that is maximum (rather than the sum of all possible decompositions). Args: x1: Coordinate of first point. x2: Coordinate of second point. Returns: The maximum energy of the reaction decomposition between the two points. """ coords = self.get_coords_in_range(x1, x2) max_energy = 0 for c in combinations(range(len(coords)), 3): i_left, i_mid, i_right = sorted(c) c_left, c_mid, c_right = coords[[i_left, i_mid, i_right]] energy = self._calculate_altitude(c_left, c_mid, c_right) if energy < max_energy: max_energy = energy return max_energy
[docs] def get_decomposition_area(self, x1: float, x2: float) -> float: """Similar to get_decomposition_energy, but returns the area enclosed instead of an energy. **WARNING**: This may be numerically unstable. Args: x1: Coordinate of first point. x2: Coordinate of second point. Returns: The area of reaction decomposition between the two points. """ coords = self.get_coords_in_range(x1, x2) if len(coords) == 2: return 0 return ConvexHull(coords, qhull_options="QJ i").volume # this is how area is defined in scipy
[docs] def get_coords_in_range(self, x1: float, x2: float): """Get the coordinates in the range [x1, x2]. Args: x1: Start of range. x2: End of range. Returns: Array of coordinates in the range. """ x_min, x_max = sorted([x1, x2]) y_min, y_max = self.get_hull_energy(x_min), self.get_hull_energy(x_max) coords = [] if np.isclose(x_min, 0.0) and not np.isclose(y_min, 0.0): coords.append([0.0, 0.0]) coords.append([x_min, y_min]) coords.extend( [ self.coords[i] for i in self.hull_vertices if self.coords[i, 0] < x_max and self.coords[i, 0] > x_min and self.coords[i, 1] <= 0 ] ) if x_max != x_min: coords.append([x_max, y_max]) if np.isclose(x_max, 1.0) and not np.isclose(y_max, 0.0): coords.append([1.0, 0.0]) coords = np.array(coords) return coords[coords[:, 0].argsort()] # type: ignore
[docs] def count(self, num: int) -> int: """Reurns the number of decomposition pathways for the interface reaction hull based on the number of **product** vertices (i.e., total # of vertices considered - 2 reactant vertices). Precomputed for hulls up to size 15. Otherwise, calls recursive implementation of counting function. Args: num: Number of product vertices. """ counts = [ 1, 1, 2, 5, 14, 42, 132, 429, 1430, 4862, 16796, 58786, 208012, 742900, 2674440, 9694845, 35357670, 129644790, 477638700, 1767263190, 6564120420, 24466267020, 91482563640, ] return counts[num] if num <= 22 else self._count_recursive(num)[0]
[docs] @lru_cache(maxsize=128) # noqa: B019 def get_decomposition_energy_and_num_paths_recursive( self, x1: float, x2: float, use_x_min_ref: bool = True, use_x_max_ref: bool = True, ) -> tuple[float, int]: """This is a recursive implementation of the get_decomposition_energy function. It significantly slower than the non-recursive implementation but is more straightforward to understand. Both should return the same answer, however the refcursive implementation also includes "free" computation of the total number of paths. The function has been cached for speed. Args: x1: Coordinate of first point. x2: Coordinate of second point. use_x_min_ref: Useful for recursive calls. If true, uses the reactant at x=0 as the reference (sometimes there is a decomposition reaction of the reactant that is lower in energy than the reactant). use_x_max_ref: Useful for recursive calls. If true, uses the reactant at x=1.0 as the reference (sometimes there is a decomposition reaction of the reactant that is lower in energy than the reactant). Returns: Tuple of decomposition energy and the number of decomposition pathways. """ all_coords = self.get_coords_in_range(x1, x2) if not use_x_min_ref and all_coords[1, 0] == 0.0: all_coords = all_coords[1:] if not use_x_max_ref and all_coords[-1, 0] == 1.0: all_coords = all_coords[:-1] x_min, y_min = all_coords[0] x_max, y_max = all_coords[-1] coords = all_coords[1:-1, :] if len(coords) == 0: return 0, 1 if len(coords) == 1: val = self._calculate_altitude([x_min, y_min], coords[0], [x_max, y_max]) # type: ignore return val, 1 val = 0 total = 0 for c in coords: if c[0] == 0.0 and not np.isclose(c[1], 0.0): use_x_min_ref = False elif c[0] == 1.0 and not np.isclose(c[1], 0.0): use_x_max_ref = False height = self._calculate_altitude([x_min, y_min], c, [x_max, y_max]) # type: ignore ( left_decomp, left_total, ) = self.get_decomposition_energy_and_num_paths_recursive(x_min, c[0], use_x_min_ref, use_x_max_ref) ( right_decomp, right_total, ) = self.get_decomposition_energy_and_num_paths_recursive(c[0], x_max, use_x_min_ref, use_x_max_ref) val += height * (left_total * right_total) + left_decomp * right_total + right_decomp * left_total total += left_total * right_total return val, total
@lru_cache(maxsize=128) # noqa: B019 def _altitude_multiplicity(self, n_left, n_right, n): """This function is used in the non-recursive implementation of the get_decomposition_energy function. It allows for rapid computation of the number times (multiplicitly) that a particular altitude appears in the decomposition. Args: n_left: number of vertices occurring in between the leftmost and middle vertices. n_right: number of vertices occurring in between the middle and rightmost vertices. n: total number of product vertices in the full interface reaction hull (i.e., not including the 2 reactant reference vertices ) """ remainder = n - n_left - n_right - 1 if remainder < 0: return 0 return self.count(n_left) * self.count(n_right) * self.count(remainder) def _count_recursive(self, n, cache=None): """A recursive implementation of the counting function. This implementation is courtesy of @mcgalcode. """ if cache is None: cache = [] if n == 0: return 1, [1] if n == 1: return 1, [1, 1] if len(cache) >= n: return cache[n], cache total = 0 biggest_cache = [] for i in range(n): left = i right = n - i - 1 left_divs, c1 = self._count_recursive(left, biggest_cache) right_divs, c2 = self._count_recursive(right, biggest_cache) if len(c1) > len(biggest_cache): biggest_cache = c1 if len(c2) > len(biggest_cache): biggest_cache = c2 total += left_divs * right_divs return total, [*biggest_cache, total] def _get_scatter(self): marker_size = 10 pts = px.scatter( x=self.coords[:, 0], y=self.coords[:, 1], hover_name=[str(r) for r in self.reactions], labels={ "x": "Mixing Ratio", "y": (r"$\Delta G_{\mathrm{rxn}} ~ \mathrm{\left(\dfrac{\mathsf{eV}}{\mathsf{atom}}\right)}$"), }, ) pts.update_traces(marker={"size": marker_size}) return pts.data[0] def _get_lines(self): coords = self.coords[self.hull.simplices] coords = coords[(coords[:, :, 1] <= 0).all(axis=1)] coords = coords[~(coords[:, :, 1] == 0).all(axis=1)] lines = [px.line(x=c[:, 0], y=c[:, 1]) for c in coords if not (c[:, 1] == 0).all()] line_data = [] for line in lines: line.update_traces(line={"color": "black"}) line_data.append(line.data[0]) return line_data @cached_property def hull_vertices(self) -> np.ndarray: """Returns numpy array of indices of the vertices that are on the hull.""" hull_vertices = [ i for i in self.hull.vertices if self.coords[i, 1] <= 0 and np.isclose( self.coords[self.coords[:, 0] == self.coords[i, 0]][:, 1].min(), self.coords[i, 1], # make sure point is lower than others on same x ) ] return np.array(hull_vertices) @cached_property def stable_reactions(self) -> list[ComputedReaction]: """Returns the reactions that are stable (on the convex hull) of the interface reaction hull. """ return [r for i, r in enumerate(self.reactions) if i in self.hull_vertices] @cached_property def unstable_reactions(self): """Returns the reactions that are unstable (NOT on the convex hull) of the interface reaction hull. """ return [r for i, r in enumerate(self.reactions) if i not in self.hull_vertices] @staticmethod def _calculate_altitude(c_left, c_mid, c_right): """Helper geometry method: calculates the altitude of a point on a line defined by three points. Args: c_left: coordinates of point 1 c_mid: coordinates of point 2 c_right: coordinates point 3 Returns: The altitude of the point """ x1, y1 = c_left x2, y2 = c_mid x3, y3 = c_right xd = (x2 - x1) / (x3 - x1) yd = y1 + xd * (y3 - y1) return y2 - yd