Merge branch 'synthesis' of github.com:materialsproject/pymatgen into…

… synthesis

pymatgen/analysis/structure_prediction/dopant_predictor.py

@@ -0,0 +1,186 @@
+import warnings
+import numpy as np
+
+from pymatgen.analysis.structure_prediction.substitution_probability import \
+    SubstitutionPredictor
+
+from pymatgen.core.periodic_table import Specie, Element
+
+
+def get_dopants_from_substitution_probabilities(structure, num_dopants=5,
+                                                threshold=0.001,
+                                                match_oxi_sign=False):
+    """
+    Get dopant suggestions based on substitution probabilities.
+
+    Args:
+        structure (Structure): A pymatgen structure decorated with
+            oxidation states.
+        num_dopants (int): The number of suggestions to return for
+            n- and p-type dopants.
+        threshold (float): Probability threshold for substitutions.
+        match_oxi_sign (bool): Whether to force the dopant and original species
+            to have the same sign of oxidation state. E.g. If the original site
+            is in a negative charge state, then only negative dopants will be
+            returned.
+
+    Returns:
+        (dict): Dopant suggestions, given as a dictionary with keys "n_type" and
+        "p_type". The suggestions for each doping type are given as a list of
+        dictionaries, each with they keys:
+
+        - "probability": The probability of substitution.
+        - "dopant_species": The dopant species.
+        - "original_species": The substituted species.
+    """
+    els_have_oxi_states = [hasattr(s, "oxi_state") for s in structure.species]
+
+    if not all(els_have_oxi_states):
+        raise ValueError("All sites in structure must have oxidation states to "
+                         "predict dopants.")
+
+    sp = SubstitutionPredictor(threshold=threshold)
+
+    subs = [sp.list_prediction([s]) for s in set(structure.species)]
+    subs = [{'probability': pred['probability'],
+             'dopant_species': list(pred['substitutions'].keys())[0],
+             'original_species': list(pred['substitutions'].values())[0]}
+            for species_preds in subs for pred in species_preds]
+    subs.sort(key=lambda x: x['probability'], reverse=True)
+
+    return _get_dopants(subs, num_dopants, match_oxi_sign)
+
+
+def get_dopants_from_shannon_radii(bonded_structure, num_dopants=5,
+                                   match_oxi_sign=False):
+    """
+    Get dopant suggestions based on Shannon radii differences.
+
+    Args:
+        bonded_structure (StructureGraph): A pymatgen structure graph
+            decorated with oxidation states. For example, generated using the
+            CrystalNN.get_bonded_structure() method.
+        num_dopants (int): The nummber of suggestions to return for
+            n- and p-type dopants.
+        match_oxi_sign (bool): Whether to force the dopant and original species
+            to have the same sign of oxidation state. E.g. If the original site
+            is in a negative charge state, then only negative dopants will be
+            returned.
+
+    Returns:
+        (dict): Dopant suggestions, given as a dictionary with keys "n_type" and
+        "p_type". The suggestions for each doping type are given as a list of
+        dictionaries, each with they keys:
+
+        - "radii_diff": The difference between the Shannon radii of the species.
+        - "dopant_spcies": The dopant species.
+        - "original_species": The substituted species.
+    """
+    # get a list of all Specie for all elements in all their common oxid states
+    all_species = [Specie(el, oxi) for el in Element
+                   for oxi in el.common_oxidation_states]
+
+    # get a series of tuples with (coordination number, specie)
+    cn_and_species = set((bonded_structure.get_coordination_of_site(i),
+                          bonded_structure.structure[i].specie)
+                         for i in range(bonded_structure.structure.num_sites))
+
+    cn_to_radii_map = {}
+    possible_dopants = []
+
+    for cn, species in cn_and_species:
+        cn_roman = _int_to_roman(cn)
+
+        try:
+            species_radius = species.get_shannon_radius(cn_roman)
+        except KeyError:
+            warnings.warn("Shannon radius not found for {} with coordination "
+                          "number {}.\nSkipping...".format(species, cn))
+            continue
+
+        if cn not in cn_to_radii_map:
+            cn_to_radii_map[cn] = _shannon_radii_from_cn(
+                all_species, cn_roman, radius_to_compare=species_radius)
+
+        shannon_radii = cn_to_radii_map[cn]
+
+        possible_dopants += [{'radii_diff': p['radii_diff'],
+                              'dopant_species': p['species'],
+                              'original_species': species}
+                             for p in shannon_radii]
+
+    possible_dopants.sort(key=lambda x: abs(x['radii_diff']))
+
+    return _get_dopants(possible_dopants, num_dopants, match_oxi_sign)
+
+
+def _get_dopants(substitutions, num_dopants, match_oxi_sign):
+    """
+    Utility method to get n- and p-type dopants from a list of substitutions.
+    """
+    n_type = [pred for pred in substitutions
+              if pred['dopant_species'].oxi_state >
+              pred['original_species'].oxi_state
+              and (not match_oxi_sign or
+                   np.sign(pred['dopant_species'].oxi_state) ==
+                   np.sign(pred['original_species'].oxi_state))]
+    p_type = [pred for pred in substitutions
+              if pred['dopant_species'].oxi_state <
+              pred['original_species'].oxi_state
+              and (not match_oxi_sign or
+                   np.sign(pred['dopant_species'].oxi_state) ==
+                   np.sign(pred['original_species'].oxi_state))]
+
+    return {'n_type': n_type[:num_dopants], 'p_type': p_type[:num_dopants]}
+
+
+def _shannon_radii_from_cn(species_list, cn_roman, radius_to_compare=0):
+    """
+    Utility func to get Shannon radii for a particular coordination number.
+
+    As the Shannon radii depends on charge state and coordination number,
+    species without an entry for a particular coordination number will
+    be skipped.
+
+    Args:
+        species_list (list): A list of Species to get the Shannon radii for.
+        cn_roman (str): The coordination number as a roman numeral. See
+            Specie.get_shannon_radius for more details.
+        radius_to_compare (float, optional): If set, the data will be returned
+            with a "radii_diff" key, containing the difference between the
+            shannon radii and this radius.
+
+    Returns:
+        (list of dict): The Shannon radii for all Species in species. Formatted
+        as a list of dictionaries, with the keys:
+
+        - "species": The species with charge state.
+        - "radius": The Shannon radius for the species.
+        - "radius_diff": The difference between the Shannon radius and the
+            radius_to_compare optional argument.
+    """
+    shannon_radii = []
+
+    for s in species_list:
+        try:
+            radius = s.get_shannon_radius(cn_roman)
+            shannon_radii.append({
+                'species': s, 'radius': radius,
+                'radii_diff': radius - radius_to_compare})
+        except KeyError:
+            pass
+
+    return shannon_radii
+
+
+def _int_to_roman(number):
+    """Utility method to convert an int (less than 20) to a roman numeral."""
+    roman_conv = [(10, "X"), (9, "IX"), (5, "V"), (4, "IV"), (1, "I")]
+
+    result = []
+    for (arabic, roman) in roman_conv:
+        (factor, number) = divmod(number, arabic)
+        result.append(roman * factor)
+        if number == 0:
+            break
+    return "".join(result)

pymatgen/analysis/structure_prediction/tests/test_dopant_predictor.py

@@ -0,0 +1,88 @@
+# coding: utf-8
+# Copyright (c) Pymatgen Development Team.
+# Distributed under the terms of the MIT License.
+
+from __future__ import unicode_literals
+
+import unittest
+
+from pymatgen.core.periodic_table import Specie
+from pymatgen.core.structure import Structure
+from pymatgen.analysis.local_env import CrystalNN
+
+from pymatgen.analysis.structure_prediction.dopant_predictor import \
+    get_dopants_from_shannon_radii, get_dopants_from_substitution_probabilities
+
+
+class DopantPredictionTest(unittest.TestCase):
+
+    def setUp(self):
+        self.tin_dioxide = Structure(
+            [3.24, 0, 0, 0, 4.83, 0, 0, 0, 4.84],
+            ['O', 'O', 'O', 'O', 'Sn', 'Sn'],
+            [[0.5, 0.19, 0.80], [0.5, 0.80, 0.19], [0, 0.30, 0.30],
+             [0, 0.69, 0.69], [0.5, 0.50, 0.50], [0, 0, 0]]
+        )
+        self.tin_dioxide.add_oxidation_state_by_element(
+            {'Sn': 4, 'O': -2})
+
+    def test_dopants_from_substitution_probabilities(self):
+        dopants = get_dopants_from_substitution_probabilities(
+            self.tin_dioxide, num_dopants=5)
+
+        self.assertTrue("n_type" in dopants)
+        self.assertTrue("p_type" in dopants)
+
+        self.assertTrue(len(dopants['n_type']) <= 5)
+
+        self.assertTrue(len(dopants['p_type']) <= 5)
+
+        self.assertAlmostEqual(dopants['n_type'][0]['probability'],
+                               0.06692682583342474)
+        self.assertEqual(dopants['n_type'][0]['dopant_species'],
+                         Specie('F', -1))
+        self.assertEqual(dopants['n_type'][0]['original_species'],
+                         Specie('O', -2))
+
+        self.assertAlmostEqual(dopants['p_type'][0]['probability'],
+                               0.023398867249112935)
+        self.assertEqual(dopants['p_type'][0]['dopant_species'],
+                         Specie('Co', 2))
+        self.assertEqual(dopants['p_type'][0]['original_species'],
+                         Specie('Sn', 4))
+
+        # test oxidation sign matching
+        dopants = get_dopants_from_substitution_probabilities(
+            self.tin_dioxide, num_dopants=15, match_oxi_sign=False)
+        self.assertEqual(dopants['n_type'][14]['dopant_species'],
+                         Specie('Li', 1))
+        self.assertEqual(dopants['n_type'][14]['original_species'],
+                         Specie('O', -2))
+
+        dopants = get_dopants_from_substitution_probabilities(
+            self.tin_dioxide, num_dopants=15, match_oxi_sign=True)
+        self.assertNotEqual(dopants['n_type'][14]['dopant_species'],
+                            Specie('Li', 1))
+
+    def test_dopants_from_shannon_radii(self):
+        bonded_structure = CrystalNN().get_bonded_structure(self.tin_dioxide)
+
+        dopants = get_dopants_from_shannon_radii(bonded_structure,
+                                                 num_dopants=5)
+
+        self.assertTrue("n_type" in dopants)
+        self.assertTrue("p_type" in dopants)
+
+        self.assertTrue(len(dopants['n_type']) <= 5)
+        self.assertTrue(len(dopants['p_type']) <= 5)
+
+        self.assertAlmostEqual(dopants['n_type'][0]['radii_diff'], 0.04)
+        self.assertEqual(dopants['n_type'][0]['dopant_species'], Specie('U', 6))
+        self.assertEqual(dopants['n_type'][0]['original_species'],
+                         Specie('Sn', 4))
+
+        self.assertAlmostEqual(dopants['p_type'][0]['radii_diff'], 0.)
+        self.assertEqual(dopants['p_type'][0]['dopant_species'],
+                         Specie('Ni', 2))
+        self.assertEqual(dopants['p_type'][0]['original_species'],
+                         Specie('Sn', 4))

pymatgen/analysis/synthesis.py

@@ -9,6 +9,7 @@
 from pymatgen import MPRester
 from pymatgen.analysis.phase_diagram import PhaseDiagram, PDPlotter, PDEntry
 from pymatgen.analysis.reaction_calculator import ComputedReaction
+from pymatgen.entries.computed_entries import ComputedEntry
 import numpy as np
 from anytree import Node, RenderTree, PreOrderIter
 
@@ -40,7 +41,10 @@ def __init__(self, entries, target, max_nelements=2):
         self.pd = PhaseDiagram(entries)
         self.stable_entries = self.pd.stable_entries
         self.max_nelements = max_nelements
-        self.target = PDEntry(target, 0)
+        if not isinstance(target, ComputedEntry):
+            self.target = PDEntry(target, 0)
+        else:
+            self.target = target
 
         def _get_tree(parent, to_remove):
             # Recursive algo to get all reactions
@@ -94,7 +98,7 @@ def get_pathways(self):
                 pathways.append(path)
 
         for p in pathways:
-            ref = p[0].rxn.calculated_reaction_energy
+            ref = p[-1].rxn.calculated_reaction_energy
             for k in p:
                 print("%s, %.3f" % (k.name, k.rxn.calculated_reaction_energy - ref))
 
@@ -114,6 +118,19 @@ def print(self):
             print(output)
 
 
+def plot_pathways(pathways):
+    from pymatgen.util.plotting import pretty_plot
+    colors = ["r", "g", "b", "c", "m", "y", "k"]
+    plt = pretty_plot(12, 8)
+    for i, p in enumerate(pathways):
+        for j, k in enumerate(p):
+            energy = k.rxn.calculated_reaction_energy
+            plt.plot([j, j+1], [-energy, -energy], color=colors[i % len(colors)], linestyle='solid')
+            plt.text(j, -energy, k.name)
+        if i > 6:
+            break
+    plt.show()
+
 from pymatgen.util.testing import PymatgenTest
 
 
@@ -140,7 +157,14 @@ def test_func(self):
 
         for rxn in a.get_unique_reactions():
             print("%s (avg_nelements = %.2f)" % (rxn.name, rxn.avg_nelements))
-        rxn_tree.get_pathways()
+
+        lfpo = MPRester().get_entries_in_chemsys(["Li", "Fe", "P", "O"])
+        lifepo4 = [e for e in lfpo if e.composition.reduced_formula == "LiFePO4"]
+        lifepo4 = min(lifepo4, key=lambda e: e.energy_per_atom)
+        rxn_tree = PDSynthesisTree(lfpo, lifepo4)
+        pathways = rxn_tree.get_pathways()
+        plot_pathways(pathways)
+
 
 
 if __name__ == "__main__":

pymatgen/core/structure.py

@@ -366,9 +366,12 @@ def __init__(self, lattice, species, coords, charge=None,
             coords (Nx3 array): list of fractional/cartesian coordinates of
                 each species.
             charge (int): overall charge of the structure. Defaults to behavior
-                in SiteCollection where total charge is the sum of the oxidation states
+                in SiteCollection where total charge is the sum of the oxidation
+                states.
             validate_proximity (bool): Whether to check if there are sites
                 that are less than 0.01 Ang apart. Defaults to False.
+            to_unit_cell (bool): Whether to map all sites into the unit cell,
+                i.e., fractional coords between 0 and 1. Defaults to False.
             coords_are_cartesian (bool): Set to True if you are providing
                 coordinates in cartesian coordinates. Defaults to False.
             site_properties (dict): Properties associated with the sites as a
@@ -2358,9 +2361,9 @@ class Structure(IStructure, collections.MutableSequence):
     """
     __hash__ = None
 
-    def __init__(self, lattice, species, coords, charge=None, validate_proximity=False,
-                 to_unit_cell=False, coords_are_cartesian=False,
-                 site_properties=None):
+    def __init__(self, lattice, species, coords, charge=None,
+                 validate_proximity=False, to_unit_cell=False,
+                 coords_are_cartesian=False, site_properties=None):
         """
         Create a periodic structure.
 
@@ -2379,9 +2382,15 @@ def __init__(self, lattice, species, coords, charge=None, validate_proximity=Fal
                 ii. List of dict of elements/species and occupancies, e.g.,
                     [{"Fe" : 0.5, "Mn":0.5}, ...]. This allows the setup of
                     disordered structures.
-            fractional_coords: list of fractional coordinates of each species.
+            coords (Nx3 array): list of fractional/cartesian coordinates of
+                each species.
+            charge (int): overall charge of the structure. Defaults to behavior
+                in SiteCollection where total charge is the sum of the oxidation
+                states.
             validate_proximity (bool): Whether to check if there are sites
                 that are less than 0.01 Ang apart. Defaults to False.
+            to_unit_cell (bool): Whether to map all sites into the unit cell,
+                i.e., fractional coords between 0 and 1. Defaults to False.
             coords_are_cartesian (bool): Set to True if you are providing
                 coordinates in cartesian coordinates. Defaults to False.
             site_properties (dict): Properties associated with the sites as a