Source code for pymatgen.analysis.interface_reactions

# coding: utf-8
# Copyright (c) Pymatgen Development Team.
# Distributed under the terms of the MIT License.

This module provides class to generate and analyze interfacial reactions.

import warnings
import numpy as np
import matplotlib.pylab as plt

from pymatgen import Composition
from pymatgen.analysis.phase_diagram import GrandPotentialPhaseDiagram
from pymatgen.analysis.reaction_calculator import Reaction

__author__ = "Yihan Xiao"
__copyright__ = "Copyright 2011, The Materials Project"
__version__ = "1.0"
__maintainer__ = "Yihan Xiao"
__email__ = ""
__status__ = "Production"
__date__ = "Aug 15 2017"

[docs]class InterfacialReactivity: """ An object encompassing all relevant data for interface reactions. """ EV_TO_KJ_PER_MOL = 96.4853 def __init__(self, c1, c2, pd, norm=True, include_no_mixing_energy=False, pd_non_grand=None, use_hull_energy=False): """ Args: c1 (Composition): Composition object for reactant 1. c2 (Composition): Composition object for reactant 2. pd (PhaseDiagram): PhaseDiagram object or GrandPotentialPhaseDiagram object built from all elements in composition c1 and c2. norm (bool): Whether or not the total number of atoms in composition of reactant will be normalized to 1. include_no_mixing_energy (bool): No_mixing_energy for a reactant is the opposite number of its energy above grand potential convex hull. In cases where reactions involve elements reservoir, this param determines whether no_mixing_energy of reactants will be included in the final reaction energy calculation. By definition, if pd is not a GrandPotentialPhaseDiagram object, this param is False. pd_non_grand (PhaseDiagram): PhaseDiagram object but not GrandPotentialPhaseDiagram object built from elements in c1 and c2. use_hull_energy (bool): Whether or not use the convex hull energy for a given composition for reaction energy calculation. If false, the energy of ground state structure will be used instead. Note that in case when ground state can not be found for a composition, convex hull energy will be used associated with a warning message. """ self.grand = isinstance(pd, GrandPotentialPhaseDiagram) # if include_no_mixing_energy is True, pd should be a # GrandPotentialPhaseDiagram object and pd_non_grand should be given. if include_no_mixing_energy and not self.grand: raise ValueError('Please provide grand phase diagram to compute' ' no_mixing_energy!') if include_no_mixing_energy and not pd_non_grand: raise ValueError('Please provide non-grand phase diagram to ' 'compute no_mixing_energy!') if self.grand and use_hull_energy and not pd_non_grand: raise ValueError('Please provide non-grand phase diagram if' ' you want to use convex hull energy.') # Keeps copy of original compositions. self.c1_original = c1 self.c2_original = c2 # Two sets of composition attributes for two processing conditions: # normalization with and without exluding element(s) from reservoir. self.c1 = c1 self.c2 = c2 self.comp1 = c1 self.comp2 = c2 self.norm = norm self.pd = pd self.pd_non_grand = pd_non_grand self.use_hull_energy = use_hull_energy # Factor is the compositional ratio between composition self.c1 and # processed composition self.comp1. E.g., factor for # Composition('SiO2') and Composition('O') is 2.0. This factor will # be used to convert mixing ratio in self.comp1 - self.comp2 # tie line to that in self.c1 - self.c2 tie line. self.factor1 = 1 self.factor2 = 1 if self.grand: # Excludes element(s) from reservoir. self.comp1 = Composition({k: v for k, v in c1.items() if k not in pd.chempots}) self.comp2 = Composition({k: v for k, v in c2.items() if k not in pd.chempots}) # Calculate the factors in case where self.grand = True and # self.norm = True. factor1 = self.comp1.num_atoms / c1.num_atoms factor2 = self.comp2.num_atoms / c2.num_atoms if self.norm: self.c1 = c1.fractional_composition self.c2 = c2.fractional_composition self.comp1 = self.comp1.fractional_composition self.comp2 = self.comp2.fractional_composition if self.grand: # Only when self.grand = True and self.norm = True # will self.factor be updated. self.factor1 = factor1 self.factor2 = factor2 # Computes energies for reactants in different scenarios. if not self.grand: if self.use_hull_energy: self.e1 = self.pd.get_hull_energy(self.comp1) self.e2 = self.pd.get_hull_energy(self.comp2) else: # Use entry energy as reactant energy if no reservoir # is present. self.e1 = InterfacialReactivity._get_entry_energy( self.pd, self.comp1) self.e2 = InterfacialReactivity._get_entry_energy( self.pd, self.comp2) else: if include_no_mixing_energy: # Computing grand potentials needs compositions containing # element(s) from reservoir, so self.c1 and self.c2 are used. self.e1 = self._get_grand_potential(self.c1) self.e2 = self._get_grand_potential(self.c2) else: self.e1 = self.pd.get_hull_energy(self.comp1) self.e2 = self.pd.get_hull_energy(self.comp2) @staticmethod def _get_entry_energy(pd, composition): """ Finds the lowest entry energy for entries matching the composition. Entries with non-negative formation energies are excluded. If no entry is found, use the convex hull energy for the composition. Args: pd (PhaseDiagram): PhaseDiagram object. composition (Composition): Composition object that the target entry should match. Returns: The lowest entry energy among entries matching the composition. """ candidate = [i.energy_per_atom for i in pd.qhull_entries if i.composition.fractional_composition == composition.fractional_composition] if not candidate: warnings.warn("The reactant " + composition.reduced_formula + " has no matching entry with negative formation" " energy, instead convex hull energy for this" " composition will be used for reaction energy " "calculation. ") return pd.get_hull_energy(composition) else: min_entry_energy = min(candidate) return min_entry_energy * composition.num_atoms def _get_grand_potential(self, composition): """ Computes the grand potential Phi at a given composition and chemical potential(s). Args: composition (Composition): Composition object. Returns: Grand potential at a given composition at chemical potential(s). """ if self.use_hull_energy: grand_potential = self.pd_non_grand.get_hull_energy(composition) else: grand_potential = InterfacialReactivity._get_entry_energy( self.pd_non_grand, composition) grand_potential -= sum([composition[e] * mu for e, mu in self.pd.chempots.items()]) if self.norm: # Normalizes energy to the composition excluding element(s) # from reservoir. grand_potential /= sum([composition[el] for el in composition if el not in self.pd.chempots]) return grand_potential def _get_energy(self, x): """ Computes reaction energy in eV/atom at mixing ratio x : (1-x) for self.comp1 : self.comp2. Args: x (float): Mixing ratio x of reactants, a float between 0 and 1. Returns: Reaction energy. """ return self.pd.get_hull_energy(self.comp1 * x + self.comp2 * (1 - x)) - self.e1 * x - self.e2 * (1 - x) def _get_reaction(self, x): """ Generates balanced reaction at mixing ratio x : (1-x) for self.comp1 : self.comp2. Args: x (float): Mixing ratio x of reactants, a float between 0 and 1. Returns: Reaction object. """ mix_comp = self.comp1 * x + self.comp2 * (1 - x) decomp = self.pd.get_decomposition(mix_comp) # Uses original composition for reactants. if np.isclose(x, 0): reactant = [self.c2_original] elif np.isclose(x, 1): reactant = [self.c1_original] else: reactant = list(set([self.c1_original, self.c2_original])) if self.grand: reactant += [Composition(e.symbol) for e, v in self.pd.chempots.items()] product = [Composition( for k, v in decomp.items()] reaction = Reaction(reactant, product) x_original = self._get_original_composition_ratio(reaction) if np.isclose(x_original, 1): reaction.normalize_to(self.c1_original, x_original) else: reaction.normalize_to(self.c2_original, 1 - x_original) return reaction def _get_elmt_amt_in_rxt(self, rxt): """ Computes total number of atoms in a reaction formula for elements not in external reservoir. This method is used in the calculation of reaction energy per mol of reaction formula. Args: rxt (Reaction): a reaction. Returns: Total number of atoms for non_reservoir elements. """ return sum([rxt.get_el_amount(e) for e in self.pd.elements])
[docs] def get_products(self): """ List of formulas of potential products. E.g., ['Li','O2','Mn']. """ products = set() for _, _, _, react, _ in self.get_kinks(): products = products.union(set([k.reduced_formula for k in react.products])) return list(products)
@staticmethod def _convert(x, factor1, factor2): """ Converts mixing ratio x in comp1 - comp2 tie line to that in c1 - c2 tie line. Args: x (float): Mixing ratio x in comp1 - comp2 tie line, a float between 0 and 1. factor1 (float): Compositional ratio between composition c1 and processed composition comp1. E.g., factor for Composition('SiO2') and Composition('O') is 2.0. factor2 (float): Compositional ratio between composition c2 and processed composition comp2. Returns: Mixing ratio in c1 - c2 tie line, a float between 0 and 1. """ return x * factor2 / ((1 - x) * factor1 + x * factor2) @staticmethod def _reverse_convert(x, factor1, factor2): """ Converts mixing ratio x in c1 - c2 tie line to that in comp1 - comp2 tie line. Args: x (float): Mixing ratio x in c1 - c2 tie line, a float between 0 and 1. factor1 (float): Compositional ratio between composition c1 and processed composition comp1. E.g., factor for Composition('SiO2') and Composition('O') is 2. factor2 (float): Compositional ratio between composition c2 and processed composition comp2. Returns: Mixing ratio in comp1 - comp2 tie line, a float between 0 and 1. """ return x * factor1 / ((1 - x) * factor2 + x * factor1)
[docs] def get_kinks(self): """ Finds all the kinks in mixing ratio where reaction products changes along the tie line of composition self.c1 and composition self.c2. Returns: Zip object of tuples (index, mixing ratio, reaction energy per atom in eV/atom, reaction formula, reaction energy per mol of reaction formula in kJ/mol). """ c1_coord = self.pd.pd_coords(self.comp1) c2_coord = self.pd.pd_coords(self.comp2) n1 = self.comp1.num_atoms n2 = self.comp2.num_atoms critical_comp = self.pd.get_critical_compositions(self.comp1, self.comp2) x_kink, energy_kink, react_kink, energy_per_rxt_formula = \ [], [], [], [] if all(c1_coord == c2_coord): x_kink = [0, 1] energy_kink = [self._get_energy(x) for x in x_kink] react_kink = [self._get_reaction(x) for x in x_kink] num_atoms = [(x * self.comp1.num_atoms + (1 - x) * self.comp2.num_atoms) for x in x_kink] energy_per_rxt_formula = [energy_kink[i] * self._get_elmt_amt_in_rxt( react_kink[i]) / num_atoms[i] * InterfacialReactivity.EV_TO_KJ_PER_MOL for i in range(2)] else: for i in reversed(critical_comp): # Gets mixing ratio x at kinks. c = self.pd.pd_coords(i) x = np.linalg.norm(c - c2_coord) / np.linalg.norm(c1_coord - c2_coord) # Modifies mixing ratio in case compositions self.comp1 and # self.comp2 are not normalized. x = x * n2 / (n1 + x * (n2 - n1)) n_atoms = x * self.comp1.num_atoms + (1 - x) * self.comp2.num_atoms # Converts mixing ratio in comp1 - comp2 tie line to that in # c1 - c2 tie line. x_converted = InterfacialReactivity._convert( x, self.factor1, self.factor2) x_kink.append(x_converted) # Gets reaction energy at kinks normalized_energy = self._get_energy(x) energy_kink.append(normalized_energy) # Gets balanced reaction at kinks rxt = self._get_reaction(x) react_kink.append(rxt) rxt_energy = normalized_energy * self._get_elmt_amt_in_rxt(rxt) / n_atoms energy_per_rxt_formula.append( rxt_energy * InterfacialReactivity.EV_TO_KJ_PER_MOL) index_kink = range(1, len(critical_comp) + 1) return zip(index_kink, x_kink, energy_kink, react_kink, energy_per_rxt_formula)
[docs] def get_critical_original_kink_ratio(self): """ Returns a list of molar mixing ratio for each kink between ORIGINAL (instead of processed) reactant compositions. This is the same list as mixing ratio obtained from get_kinks method if self.norm = False. Returns: A list of floats representing molar mixing ratios between the original reactant compositions for each kink. """ ratios = [] if self.c1_original == self.c2_original: return [0, 1] reaction_kink = [k[3] for k in self.get_kinks()] for rxt in reaction_kink: ratios.append(abs(self._get_original_composition_ratio(rxt))) return ratios
def _get_original_composition_ratio(self, reaction): """ Returns the molar mixing ratio between the reactants with ORIGINAL ( instead of processed) compositions for a reaction. Args: reaction (Reaction): Reaction object that contains the original reactant compositions. Returns: The molar mixing ratio between the original reactant compositions for a reaction. """ if self.c1_original == self.c2_original: return 1 c1_coeff = reaction.get_coeff(self.c1_original) \ if self.c1_original in reaction.reactants else 0 c2_coeff = reaction.get_coeff(self.c2_original) \ if self.c2_original in reaction.reactants else 0 return c1_coeff * 1.0 / (c1_coeff + c2_coeff)
[docs] def labels(self): """ Returns a dictionary containing kink information: {index: 'x= mixing_ratio energy= reaction_energy reaction_equation'}. E.g., {1: 'x= 0.0 energy = 0.0 Mn -> Mn', 2: 'x= 0.5 energy = -15.0 O2 + Mn -> MnO2', 3: 'x= 1.0 energy = 0.0 O2 -> O2'}. """ return {j: 'x= ' + str(round(x, 4)) + ' energy in eV/atom = ' + str(round(energy, 4)) + ' ' + str(reaction) for j, x, energy, reaction, _ in self.get_kinks()}
[docs] def plot(self): """ Plots reaction energy as a function of mixing ratio x in self.c1 - self.c2 tie line using pylab. Returns: Pylab object that plots reaction energy as a function of mixing ratio x. """ plt.rcParams['xtick.major.pad'] = '6' plt.rcParams['ytick.major.pad'] = '6' plt.rcParams['axes.linewidth'] = 2 npoint = 1000 xs = np.linspace(0, 1, npoint) # Converts sampling points in self.c1 - self.c2 tie line to those in # self.comp1 - self.comp2 tie line. xs_reverse_converted = InterfacialReactivity._reverse_convert( xs, self.factor1, self.factor2) energies = [self._get_energy(x) for x in xs_reverse_converted] plt.plot(xs, energies, 'k-') # Marks kinks and minimum energy point. kinks = self.get_kinks() _, x_kink, energy_kink, _, _ = zip(*kinks) plt.scatter(x_kink, energy_kink, marker='o', c='blue', s=20) plt.scatter(self.minimum()[0], self.minimum()[1], marker='*', c='red', s=300) # Labels kinks with indices. Labels are made draggable # in case of overlapping. for index, x, energy, _, _ in kinks: plt.annotate( index, xy=(x, energy), xytext=(5, 30), textcoords='offset points', ha='right', va='bottom', arrowprops=dict(arrowstyle='->', connectionstyle='arc3,rad=0')).draggable() plt.xlim([-0.05, 1.05]) if self.norm: plt.ylabel('Energy (eV/atom)') else: plt.ylabel('Energy (eV/f.u.)') plt.xlabel('$x$ in $x$ {} + $(1-x)$ {}'.format( self.c1.reduced_formula, self.c2.reduced_formula)) return plt
[docs] def minimum(self): """ Finds the minimum reaction energy E_min and corresponding mixing ratio x_min. Returns: Tuple (x_min, E_min). """ return min([(x, energy) for _, x, energy, _, _ in self.get_kinks()], key=lambda i: i[1])
[docs] def get_no_mixing_energy(self): """ Generates the opposite number of energy above grand potential convex hull for both reactants. Returns: [(reactant1, no_mixing_energy1),(reactant2,no_mixing_energy2)]. """ assert self.grand == 1, 'Please provide grand potential phase diagram for computing no_mixing_energy!' energy1 = self.pd.get_hull_energy(self.comp1) - self._get_grand_potential(self.c1) energy2 = self.pd.get_hull_energy(self.comp2) - self._get_grand_potential(self.c2) unit = 'eV/f.u.' if self.norm: unit = 'eV/atom' return [(self.c1_original.reduced_formula + ' ({0})'.format(unit), energy1), (self.c2_original.reduced_formula + ' ({0})'.format(unit), energy2)]
[docs] @staticmethod def get_chempot_correction(element, temp, pres): """ Get the normalized correction term Δμ for chemical potential of a gas phase consisting of element at given temperature and pressure, referenced to that in the standard state (T_std = 298.15 K, T_std = 1 bar). The gas phase is limited to be one of O2, N2, Cl2, F2, H2. Calculation formula can be found in the documentation of Materials Project website. Args: element (string): The string representing the element. temp (float): The temperature of the gas phase. pres (float): The pressure of the gas phase. Returns: The correction of chemical potential in eV/atom of the gas phase at given temperature and pressure. """ if element not in ["O", "N", "Cl", "F", "H"]: return 0 std_temp = 298.15 std_pres = 1E5 ideal_gas_const = 8.3144598 # Cp and S at standard state in J/(K.mol). Data from # # # # # Cp_dict = {"O": 29.376, "N": 29.124, "Cl": 33.949, "F": 31.302, "H": 28.836} S_dict = {"O": 205.147, "N": 191.609, "Cl": 223.079, "F": 202.789, "H": 130.680} Cp_std = Cp_dict[element] S_std = S_dict[element] PV_correction = ideal_gas_const * temp * np.log(pres / std_pres) TS_correction = - Cp_std * (temp * np.log(temp) - std_temp * np.log(std_temp)) \ + Cp_std * (temp - std_temp) * (1 + np.log(std_temp)) \ - S_std * (temp - std_temp) dG = PV_correction + TS_correction # Convert to eV/molecule unit. dG /= 1000 * InterfacialReactivity.EV_TO_KJ_PER_MOL # Normalize by number of atoms in the gas molecule. For elements # considered, the gas molecules are all diatomic. dG /= 2 return dG