Pourbaix binary fix (materialsproject#941)

* update maker to sieve initial entries based on phase stability

* edited unittests, added test plot

* cleanup test_maker

* mod analyzer test to make selection of entries a bit more robust

* cleanup analyzer

* fixes entry issue with normalization

* update the example notebook

* refactor PhaseDiagram filtering to make it optional for multicomponent diagrams

* minor tweak to example

* add proper warnings for multielement ehull finding

* change ehull count back and add binary analyzer ghull test

* few cosmetic things and added entry finder for multi

* added binary unittest for plotter

pymatgen/analysis/pourbaix/analyzer.py

@@ -11,6 +11,7 @@
 from scipy.spatial import HalfspaceIntersection, ConvexHull
 from pymatgen.analysis.pourbaix.entry import MultiEntry
 from six.moves import zip
+import warnings
 
 """
 Class for analyzing Pourbaix Diagrams. Similar to PDAnalyzer
@@ -238,13 +239,13 @@ def get_all_decomp_and_e_above_hull(self, single_entry):
         decomp_entries, hull_energies, decomp_species, entries = [], [], [], []
 
         # for all entries where the material is the only solid
-        if not isinstance(self._pd.all_entries[0], MultiEntry):
+        if not self._pd._multielement:
            possible_entries = [e for e in self._pd.all_entries
-                             if single_entry == e]
+                               if single_entry == e]
         else:
            possible_entries = [e for e in self._pd.all_entries
-                         if e.phases.count("Solid") == 1
-                         and single_entry in e.entrylist]
+                               if e.phases.count("Solid") == 1
+                               and single_entry in e.entrylist]
 
         for possible_entry in possible_entries:
             # Find the decomposition details if the material
@@ -341,16 +342,36 @@ def get_gibbs_free_energy(self, pH, V):
         stable_entry = [k for k, v in data.items() if v == gibbs_energy]
         return (gibbs_energy, stable_entry)
 
+    def _min_multientry_from_single_entry(self, single_entry):
+        """
+        Gives lowest energy multi-entry from single entry
+
+        Args:
+            single_entry (PourbaixEntry): pourbaix entry to find valid
+                multientries from
+        """
+        de, ehulls, entries = self.get_all_decomp_and_e_above_hull(single_entry)
+        if not ehulls:
+            raise ValueError("No entries where {} is the only solid".format(
+                             single_entry.name))
+        return entries[np.argmin(ehulls)]
+
     def get_entry_stability(self, entry, pH, V):
         """
         Get the energy difference between an entry and the
         most stable decomposition product (i.e. the pourbaix-stable
         entry) at a given pH and voltage.
 
         Args:
-            entry (PourbaixEntry): Pourbaix entry corresponding to the
-                stability to be calculated
+            entry (PourbaixEntry): Pourbaix entry or MultiEntry
+                corresponding to the stability to be calculated
             pH (float): pH at which to calculate stability of entry
             V (float): voltage at which to calculate stability of entry
         """
+        if self._pd._multielement and not isinstance(entry, MultiEntry):
+            _, _, entries = self.get_all_decomp_and_e_above_hull(entry)
+            warnings.warn("{} is not a multi-entry, calculating stability of "
+                          "representative {} multientry")
+            gs = [self.g(e, pH, V) for e in entries]
+            return min(gs) - self.get_gibbs_free_energy(pH, V)[0]
         return self.g(entry, pH, V) - self.get_gibbs_free_energy(pH, V)[0]

pymatgen/analysis/pourbaix/entry.py

@@ -81,6 +81,13 @@ def __init__(self, entry, correction=0.0, entry_id=None):
     def energy(self):
         return self.uncorrected_energy + self.correction
 
+    @property
+    def energy_per_atom(self):
+        # This is slightly different because it's normalized by formula unit
+        # TODO: remember to fix this if this module ever gets refactored
+        #       to be more consistent with the other PD infrastructure
+        return self.energy / self.composition.reduced_composition.num_atoms
+
     @property
     def g0(self):
         """

pymatgen/analysis/pourbaix/maker.py

@@ -12,7 +12,9 @@
 from pymatgen.analysis.pourbaix.entry import MultiEntry, ion_or_solid_comp_object
 from pymatgen.core.periodic_table import Element
 from pymatgen.core.composition import Composition
+from pymatgen.entries.computed_entries import ComputedEntry
 from pymatgen.analysis.reaction_calculator import Reaction, ReactionError
+from pymatgen.analysis.phase_diagram import PhaseDiagram
 
 """
 Module containing analysis classes which compute a pourbaix diagram given a
@@ -35,7 +37,7 @@
 
 PREFAC = 0.0591
 MU_H2O = -2.4583
-
+elements_HO = {Element('H'), Element('O')}
 # TODO: There's a lot of functionality here that diverges
 #   based on whether or not the pbx diagram is multielement
 #   or not.  Could be a more elegant way to
@@ -44,71 +46,75 @@
 class PourbaixDiagram(object):
     """
     Class to create a Pourbaix diagram from entries
+
     Args:
-        entries: Entries list containing both Solids and Ions
-        comp_dict: Dictionary of compositions
+        entries [Entry]: Entries list containing both Solids and Ions
+        comp_dict {str: float}: Dictionary of compositions, defaults to
+            equal parts of each elements
+        conc_dict {str: float}: Dictionary of ion concentrations, defaults
+            to 1e-6 for each element
+        filter_multielement (bool): applying this filter to a multi-
+            element pourbaix diagram makes generates it a bit more
+            efficiently by filtering the entries used to generate
+            the hull.  This breaks some of the functionality of
+            the analyzer, though, so use with caution.
     """
-    def __init__(self, entries, comp_dict=None):
-        self._solid_entries = list()
-        self._ion_entries = list()
-        for entry in entries:
-            if entry.phase_type == "Solid":
-                self._solid_entries.append(entry)
-            elif entry.phase_type == "Ion":
-                self._ion_entries.append(entry)
-            else:
-                raise Exception("Incorrect Phase type - needs to be \
-                Pourbaix entry of phase type Ion/Solid")
-        if len(self._ion_entries) == 0:
-            raise Exception("No ion phase. Equilibrium between ion/solid "
-                            "is required to make a Pourbaix Diagram")
-        self._unprocessed_entries = self._solid_entries + self._ion_entries
+    def __init__(self, entries, comp_dict=None, conc_dict=None,
+                 filter_multielement=False):
+        # Get non-OH elements
+        pbx_elts = set(itertools.chain.from_iterable(
+            [entry.composition.elements for entry in entries]))
+        pbx_elts = list(pbx_elts - elements_HO)
+
+        # Set default conc/comp dicts
+        if not comp_dict:
+            comp_dict = {elt.symbol : 1. / len(pbx_elts) for elt in pbx_elts}
+        if not conc_dict:
+            conc_dict = {elt.symbol : 1e-6 for elt in pbx_elts}
+
         self._elt_comp = comp_dict
+        self.pourbaix_elements = pbx_elts
 
-        if comp_dict and len(comp_dict) > 1:
-            self._multielement = True
-            pbx_elements = set()
-            for comp in comp_dict.keys():
-                for el in [el for el in
-                           ion_or_solid_comp_object(comp).elements
-                           if el not in ["H", "O"]]:
-                    pbx_elements.add(el.symbol)
-            self.pourbaix_elements = pbx_elements
-            w = [comp_dict[key] for key in comp_dict]
-            A = []
-            for comp in comp_dict:
-                comp_obj = ion_or_solid_comp_object(comp)
-                Ai = []
-                for elt in self.pourbaix_elements:
-                    Ai.append(comp_obj[elt])
-                A.append(Ai)
-            A = np.array(A).T.astype(float)
-            w = np.array(w)
-            A /= np.dot([a.sum() for a in A], w)
-            x = np.linalg.solve(A, w)
-            self._elt_comp = dict(zip(self.pourbaix_elements, x))
+        solid_entries = [entry for entry in entries
+                         if entry.phase_type == "Solid"]
+        ion_entries = [entry for entry in entries
+                       if entry.phase_type == "Ion"]
+
+        for entry in ion_entries:
+            ion_elts = list(set(entry.composition.elements) - elements_HO)
+            if len(ion_elts) != 1:
+                raise ValueError("Elemental concentration not compatible "
+                                 "with multi-element ions")
+            entry.conc = conc_dict[ion_elts[0].symbol]
+
+        if not len(solid_entries + ion_entries) == len(entries):
+            raise ValueError("All supplied entries must have a phase type of "
+                             "either \"Solid\" or \"Ion\"")
 
+        self._unprocessed_entries = entries
+
+        if len(comp_dict) > 1:
+            self._multielement = True
+            if filter_multielement:
+                # Add two high-energy H/O entries that ensure the hull
+                # includes all stable solids.
+                entries_HO = [ComputedEntry('H', 10000), ComputedEntry('O', 10000)]
+                solid_pd = PhaseDiagram(solid_entries + entries_HO)
+                solid_entries = list(set(solid_pd.stable_entries) - set(entries_HO))
+            self._processed_entries = self._generate_multielement_entries(
+                    solid_entries + ion_entries)
         else:
             self._multielement = False
-            pb_els = set(itertools.chain.from_iterable(
-                            [e.composition for e in entries]))
-            pb_els -= {Element('H'), Element('O')}
-            self.pourbaix_elements = [e.name for e in list(pb_els)]
-            # TODO: document the physical meaning
-            # of ternary oxide of pbx diagrams in both modes
-            self._elt_comp = {self.pourbaix_elements[0]: 1.0}
+
+            self._processed_entries = solid_entries + ion_entries
         self._make_pourbaix_diagram()
 
     def _create_conv_hull_data(self):
         """
         Make data conducive to convex hull generator.
         """
-        if self._multielement:
-            self._all_entries = self._process_multielement_entries()
-        else:
-            self._all_entries = self._unprocessed_entries
         entries_to_process = list()
-        for entry in self._all_entries:
+        for entry in self._processed_entries:
             entry.scale(entry.normalization_factor)
             entry.correction += (- MU_H2O * entry.nH2O + entry.conc_term)
             entries_to_process.append(entry)
@@ -129,23 +135,25 @@ def _process_conv_hull_data(self, entries_to_process):
         [data, self._qhull_entries] = list(zip(*temp))
         return data
 
-    def _process_multielement_entries(self):
+    def _generate_multielement_entries(self, entries):
         """
         Create entries for multi-element Pourbaix construction.
 
         This works by finding all possible linear combinations
         of entries that can result in the specified composition
         from the initialized comp_dict.
+
+        Args:
+            entries ([PourbaixEntries]): list of pourbaix entries
+                to process into MultiEntries
         """
         N = len(self._elt_comp)  # No. of elements
-        entries = self._unprocessed_entries
-        dummy_prod = Composition(self._elt_comp) 
         total_comp = Composition(self._elt_comp)
 
         # generate all possible combinations of compounds that have all elts
         entry_combos = [itertools.combinations(entries, j+1) for j in range(N)]
         entry_combos = itertools.chain.from_iterable(entry_combos)
-        entry_combos = filter(lambda x: dummy_prod < MultiEntry(x).total_composition, 
+        entry_combos = filter(lambda x: total_comp < MultiEntry(x).total_composition, 
                               entry_combos)
 
         # Generate and filter entries
@@ -311,9 +319,9 @@ def unstable_entries(self):
     @property
     def all_entries(self):
         """
-        Return all entries
+        Return all entries used to generate the pourbaix diagram
         """
-        return self._all_entries
+        return self._processed_entries
 
     @property
     def vertices(self):

pymatgen/analysis/pourbaix/plotter.py

@@ -16,6 +16,7 @@
 import six
 from six.moves import map
 from six.moves import zip
+import warnings
 
 __author__ = "Sai Jayaraman"
 __copyright__ = "Copyright 2011, The Materials Project"
@@ -211,9 +212,21 @@ def plot_entry_stability(self, entry, resolution=100, e_hull_max=1,
 
         # Get stability map
         phs, vs = np.mgrid[ph_min:ph_max:resolution*1j, 
-                         v_min:v_max:resolution*1j]
-        stability = [self._analyzer.get_entry_stability(entry, ph, v)
-                     for ph, v in zip(phs.ravel(), vs.ravel())]
+                           v_min:v_max:resolution*1j]
+        if self._pd._multielement and not isinstance(entry, MultiEntry):
+            _, _, m_entries = self._analyzer.\
+                    get_all_decomp_and_e_above_hull(entry)
+            warnings.warn("{} is not a multi-entry, calculating stability of "
+                          "representative multientries".format(entry.name))
+            stability = [self._analyzer.get_entry_stability(m_entries[0], ph, v)
+                         for ph, v in zip(phs.ravel(), vs.ravel())]
+            for m_entry in m_entries[1:]:
+                new_stab = [self._analyzer.get_entry_stability(m_entry, ph, v)
+                            for ph, v in zip(phs.ravel(), vs.ravel())]
+                stability = np.min(np.array([stability, new_stab]), axis=0)
+        else:
+            stability = [self._analyzer.get_entry_stability(entry, ph, v)
+                         for ph, v in zip(phs.ravel(), vs.ravel())]
         stability = np.array(stability).reshape(phs.shape)
 
         # Plot stability map
@@ -390,7 +403,6 @@ def get_pourbaix_plot(self, limits=[[-2, 14], [-3, 3]],
             plt:
                 matplotlib plot object
         """
-#        plt = pretty_plot(24, 14.4)
         plt = pretty_plot(16)
         (stable, unstable) = self.pourbaix_plot_data(limits)
         if limits:

pymatgen/analysis/pourbaix/tests/test_analyzer.py

@@ -14,6 +14,7 @@
 
 test_dir = os.path.join(os.path.dirname(__file__), '..', '..', '..', '..', 'test_files')
 
+# TODO: refactor with more sensible binary/ternary test data
 class TestPourbaixAnalyzer(unittest.TestCase):
 
     def setUp(self):
@@ -36,7 +37,8 @@ def test_get_facet_chempots(self):
             self.assertEqual(len(range_map_dict[entry]), self.num_simplices[entry])
 
     def test_get_decomp(self):
-        for entry in [entry for entry in self.pd.all_entries if entry not in self.pd.stable_entries]:
+        for entry in [entry for entry in self.pd.all_entries
+                      if entry not in self.pd.stable_entries]:
             decomp_entries = self.analyzer.get_decomposition(entry)
             for entr in decomp_entries:
                 self.assertEqual(decomp_entries[entr], self.decomp_test[entry.name][entr.name])
@@ -49,11 +51,15 @@ def test_binary(self):
                                     comp_dict = {"Ag": 0.5, "Te": 0.5})
         analyzer_binary = PourbaixAnalyzer(pd_binary)
 
-        te_entry = pd_binary._unprocessed_entries[4]
+        te_entry = [e for e in pd_binary._unprocessed_entries
+                    if e.composition.formula == "Te3"][0]
         de, hull_e, entries = analyzer_binary.get_all_decomp_and_e_above_hull(te_entry)
         self.assertEqual(len(de), 10)
         self.assertEqual(len(hull_e), 10)
-        self.assertAlmostEqual(hull_e[0], 5.4419792326439893)
+        # Find a specific multientry to test
+        tuples = zip(de, hull_e, entries)
+        test_tuple = [t for t in tuples if t[2].name=='Te(s) + Ag[2+]'][0]
+        self.assertAlmostEqual(test_tuple[1], 5.1396968548627315)
 
     def test_ternary(self):
         # Ternary
@@ -66,20 +72,20 @@ def test_ternary(self):
         de, hull_e, entries = analyzer_ternary.get_all_decomp_and_e_above_hull(te_entry)
         self.assertEqual(len(de), 116)
         self.assertEqual(len(hull_e), 116)
-        self.assertAlmostEqual(hull_e[0], 29.2520325229)
-
-    def test_v_fe(self):
-        v_fe_entries = self.multi_data['v_fe']
-        pd = PourbaixDiagram(v_fe_entries, comp_dict = {"Fe": 0.5, "V": 0.5})
-        analyzer_vfe = PourbaixAnalyzer(pd)
-        entries = sorted(pd._all_entries, key=lambda x: x.energy)
-        self.assertAlmostEqual(entries[1].weights[1], 0.6666666666)
-        self.assertAlmostEqual(entries[1].energy, -110.77582628499995)
-        self.assertAlmostEqual(entries[100].energy, -20.685496454465955)
+        tuples = zip(de, hull_e, entries)
+        test_tuple = [t for t in tuples if t[2].name=='N2(s) + TeO4[2-] + Ag[2+]'][0]
+        self.assertAlmostEqual(test_tuple[1], 50.337069095866745)
 
     def test_get_entry_stability(self):
         stab = self.analyzer.get_entry_stability(self.pd.all_entries[0], pH=0, V=1)
         self.assertAlmostEqual(stab, 3.88159999)
 
+        # binary
+        pd_binary = PourbaixDiagram(self.multi_data['binary'],
+                                    comp_dict = {"Ag": 0.5, "Te": 0.5})
+        analyzer_binary = PourbaixAnalyzer(pd_binary)
+        ghull = analyzer_binary.get_entry_stability(self.multi_data['binary'][16], 
+                                                    0, -5)
+
 if __name__ == '__main__':
     unittest.main()

pymatgen/analysis/pourbaix/tests/test_maker.py

@@ -19,13 +19,14 @@ def setUp(self):
                                                     "test_entries.csv"))
         self._entries = entries
         self._pd = PourbaixDiagram(entries)
-        self.list_of_stable_entries = ["ZnO(s)", "Zn[2+]", "ZnO2(s)", "ZnHO2[-]", "ZnO2[2-]", "Zn(s)"]
+        self.list_of_stable_entries = ["ZnO(s)", "ZnO2(s)", "Zn[2+]", "ZnHO2[-]", "ZnO2[2-]", "Zn(s)"]
 
 
     def test_pourbaix_diagram(self):
         self.assertEqual(len(self._pd.facets), 6, "Incorrect number of facets")
-        for entry in self._pd.stable_entries:
-            self.assertIn(entry.name, self.list_of_stable_entries, "List of stable entries does not match")
+        self.assertEqual(set([e.name for e in self._pd.stable_entries]), 
+                         set(self.list_of_stable_entries), "List of stable entries does not match")
+
 
 if __name__ == '__main__':
     unittest.main()