Merge pull request #1558 from awvio/master

Uncertainty quantification in energy corrections

docs/genindex.html

@@ -634,7 +634,7 @@ <h2 id="A">A</h2>
 </li>
       <li><a href="pymatgen.analysis.chemenv.coordination_environments.chemenv_strategies.html#pymatgen.analysis.chemenv.coordination_environments.chemenv_strategies.NormalizedAngleDistanceNbSetWeight.angninvndist">angninvndist() (NormalizedAngleDistanceNbSetWeight method)</a>
 </li>
-      <li><a href="pymatgen.entries.compatibility.html#pymatgen.entries.compatibility.AnionCorrection">AnionCorrection (class in pymatgen.entries.compatibility)</a>
+      <li><a href="pymatgen.entries.compatibility.html#pymatgen.entries.compatibility.CompositionCorrection">CompositionCorrection (class in pymatgen.entries.compatibility)</a>
 </li>
       <li><a href="pymatgen.analysis.wulff.html#pymatgen.analysis.wulff.WulffShape.anisotropy">anisotropy() (WulffShape property)</a>
 </li>

docs/pymatgen.entries.compatibility.html

@@ -206,6 +206,10 @@
 <span id="pymatgen-entries-compatibility-module"></span><h1>pymatgen.entries.compatibility module<a class="headerlink" href="#module-pymatgen.entries.compatibility" title="Permalink to this headline">¶</a></h1>
 <p>This module implements Compatibility corrections for mixing runs of different
 functionals.</p>
+<dl class="class">
+<dt id="pymatgen.entries.compatibility.CompositionCorrection">
+<em class="property">class </em><code class="sig-name descname">CompositionCorrection</code><span class="sig-paren">(</span><em class="sig-param">*args</em>, <em class="sig-param">**kwargs</em><span class="sig-paren">)</span><a class="reference internal" href="_modules/pymatgen/entries/compatibility.html#CompositionCorrection"><span class="viewcode-link">[source]</span></a><a class="headerlink" href="#pymatgen.entries.compatibility.CompositionCorrection" title="Permalink to this definition">¶</a></dt>
+<dd><p>Bases: <a class="reference internal" href="#pymatgen.entries.compatibility.CompositionCorrection" title="pymatgen.entries.compatibility.CompositionCorrection"><code class="xref py py-class docutils literal notranslate"><span class="pre">pymatgen.entries.compatibility.CompositionCorrection</span></code></a></p>
 <dl class="py class">
 <dt id="pymatgen.entries.compatibility.AnionCorrection">
 <em class="property">class </em><code class="sig-name descname">AnionCorrection</code><span class="sig-paren">(</span><em class="sig-param"><span class="o">*</span><span class="n">args</span></em>, <em class="sig-param"><span class="o">**</span><span class="n">kwargs</span></em><span class="sig-paren">)</span><a class="reference external" href="https://github.com/materialsproject/pymatgen/blob/v2020.7.18/pymatgen/entries/compatibility.py#L200-L287"><span class="viewcode-link">[source]</span></a><a class="headerlink" href="#pymatgen.entries.compatibility.AnionCorrection" title="Permalink to this definition">¶</a></dt>

pymatgen/analysis/reaction_calculator.py

@@ -14,6 +14,7 @@
 from monty.fractions import gcd_float
 from monty.json import MSONable
 from monty.json import MontyDecoder
+from uncertainties import ufloat
 
 from pymatgen.core.composition import Composition
 from pymatgen.entries.computed_entries import ComputedEntry
@@ -111,8 +112,10 @@ def normalize_to_element(self, element, factor=1):
         """
         all_comp = self._all_comp
         coeffs = self._coeffs
-        current_el_amount = (sum([all_comp[i][element] * abs(coeffs[i])
-                                  for i in range(len(all_comp))]) / 2)
+        current_el_amount = (
+            sum([all_comp[i][element] * abs(coeffs[i]) for i in range(len(all_comp))])
+            / 2
+        )
         scale_factor = factor / current_el_amount
         self._coeffs = [c * scale_factor for c in coeffs]
 
@@ -126,8 +129,15 @@ def get_el_amount(self, element):
         Returns:
             Amount of that element in the reaction.
         """
-        return (sum([self._all_comp[i][element] * abs(self._coeffs[i])
-                    for i in range(len(self._all_comp))]) / 2)
+        return (
+            sum(
+                [
+                    self._all_comp[i][element] * abs(self._coeffs[i])
+                    for i in range(len(self._all_comp))
+                ]
+            )
+            / 2
+        )
 
     @property
     def elements(self):
@@ -155,16 +165,18 @@ def reactants(self):
         """
         List of reactants
         """
-        return [self._all_comp[i] for i in range(len(self._all_comp))
-                if self._coeffs[i] < 0]
+        return [
+            self._all_comp[i] for i in range(len(self._all_comp)) if self._coeffs[i] < 0
+        ]
 
     @property
     def products(self):
         """
         List of products
         """
-        return [self._all_comp[i] for i in range(len(self._all_comp))
-                if self._coeffs[i] > 0]
+        return [
+            self._all_comp[i] for i in range(len(self._all_comp)) if self._coeffs[i] > 0
+        ]
 
     def get_coeff(self, comp):
         """
@@ -240,8 +252,9 @@ def as_entry(self, energies):
         Returns a ComputedEntry representation of the reaction.
         :return:
         """
-        relevant_comp = [comp * abs(coeff) for coeff, comp
-                         in zip(self._coeffs, self._all_comp)]
+        relevant_comp = [
+            comp * abs(coeff) for coeff, comp in zip(self._coeffs, self._all_comp)
+        ]
         comp = sum(relevant_comp, Composition())
         entry = ComputedEntry(0.5 * comp, self.calculate_energy(energies))
         entry.name = self.__str__()
@@ -334,7 +347,9 @@ def __init__(self, reactants, products):
         diff = self._num_comp - rank
         num_constraints = diff if diff >= 2 else 1
 
-        self._lowest_num_errors = np.inf  # an error = a component changing sides or disappearing
+        self._lowest_num_errors = (
+            np.inf
+        )  # an error = a component changing sides or disappearing
 
         self._coeffs = self._balance_coeffs(comp_matrix, num_constraints)
         self._els = all_elems
@@ -373,11 +388,16 @@ def _balance_coeffs(self, comp_matrix, max_num_constraints):
 
             coeffs = np.matmul(np.linalg.pinv(comp_and_constraints), b)
 
-            if np.allclose(np.matmul(comp_matrix, coeffs), np.zeros((self._num_elems, 1))):
+            if np.allclose(
+                np.matmul(comp_matrix, coeffs), np.zeros((self._num_elems, 1))
+            ):
                 balanced = True
-                expected_signs = np.array([-1] * len(self._input_reactants) +
-                                          [+1] * len(self._input_products))
-                num_errors = np.sum(np.multiply(expected_signs, coeffs.T) < self.TOLERANCE)
+                expected_signs = np.array(
+                    [-1] * len(self._input_reactants) + [+1] * len(self._input_products)
+                )
+                num_errors = np.sum(
+                    np.multiply(expected_signs, coeffs.T) < self.TOLERANCE
+                )
 
                 if num_errors == 0:
                     self._lowest_num_errors = 0
@@ -458,9 +478,13 @@ def __init__(self, reactant_entries, product_entries):
         self._reactant_entries = reactant_entries
         self._product_entries = product_entries
         self._all_entries = reactant_entries + product_entries
-        reactant_comp = [e.composition.get_reduced_composition_and_factor()[0] for e in reactant_entries]
-        product_comp = [e.composition.get_reduced_composition_and_factor()[0] for e in product_entries]
-        super().__init__(reactant_comp, product_comp)
+        reactant_comp = [e.composition.get_reduced_composition_and_factor()[0]
+                         for e in reactant_entries]
+
+        product_comp = [e.composition.get_reduced_composition_and_factor()[0]
+                        for e in product_entries]
+
+        super().__init__(list(reactant_comp), list(product_comp))
 
     @property
     def all_entries(self):
@@ -487,10 +511,31 @@ def calculated_reaction_energy(self):
         for entry in self._reactant_entries + self._product_entries:
             (comp, factor) = entry.composition.get_reduced_composition_and_factor()
             calc_energies[comp] = min(
-                calc_energies.get(comp, float("inf")), entry.energy / factor)
+                calc_energies.get(comp, float("inf")), entry.energy / factor
+            )
 
         return self.calculate_energy(calc_energies)
 
+    @property
+    def calculated_reaction_energy_uncertainty(self):
+        """
+        Calculates the uncertainty in the reaction energy based on the uncertainty in the
+        energies of the products and reactants
+        """
+
+        calc_energies = {}
+
+        for entry in self._reactant_entries + self._product_entries:
+            (comp, factor) = entry.composition.get_reduced_composition_and_factor()
+            energy_ufloat = ufloat(
+                entry.energy, entry.correction_uncertainty
+            )
+            calc_energies[comp] = min(
+                calc_energies.get(comp, float("inf")), energy_ufloat / factor
+            )
+
+        return self.calculate_energy(calc_energies).std_dev
+
     def as_dict(self):
         """
         Returns:

pymatgen/analysis/structure_analyzer.py

@@ -16,7 +16,6 @@
 
 import numpy as np
 from scipy.spatial import Voronoi
-
 from monty.dev import deprecated
 
 from pymatgen import Element, Specie, Composition
@@ -596,13 +595,13 @@ def oxide_type(structure, relative_cutoff=1.1, return_nbonds=False):
 
 def sulfide_type(structure):
     """
-    Determines if a structure is a sulfide/polysulfide
+    Determines if a structure is a sulfide/polysulfide/sulfate
 
     Args:
         structure (Structure): Input structure.
 
     Returns:
-        (str) sulfide/polysulfide/sulfate
+        (str) sulfide/polysulfide or None if structure is a sulfate.
     """
     structure = structure.copy()
     structure.remove_oxidation_states()