Python pdanalyzer.PDAnalyzer类代码示例

 def test_1d_pd(self):
     entry = PDEntry("H", 0)
     pd = PhaseDiagram([entry])
     pda = PDAnalyzer(pd)
     decomp, e = pda.get_decomp_and_e_above_hull(PDEntry("H", 1))
     self.assertAlmostEqual(e, 1)
     self.assertAlmostEqual(decomp[entry], 1.0)

 def test_get_stability(self):
     entries = self.rester.get_entries("Fe-O")
     modified_entries = []
     for entry in entries:
         # Create modified entries with energies that are 0.01eV higher
         # than the corresponding entries.
         if entry.composition.reduced_formula == "Fe2O3":
             modified_entries.append(
                 ComputedEntry(entry.composition,
                               entry.uncorrected_energy + 0.01,
                               parameters=entry.parameters,
                               entry_id="mod_{}".format(entry.entry_id)))
     rest_ehulls = self.rester.get_stability(modified_entries)
     all_entries = entries + modified_entries
     compat = MaterialsProjectCompatibility()
     all_entries = compat.process_entries(all_entries)
     pd = PhaseDiagram(all_entries)
     a = PDAnalyzer(pd)
     for e in all_entries:
         if str(e.entry_id).startswith("mod"):
             for d in rest_ehulls:
                 if d["entry_id"] == e.entry_id:
                     data = d
                     break
             self.assertAlmostEqual(a.get_e_above_hull(e),
                                    data["e_above_hull"])

def main(comp="La0.5Sr0.5MnO3", energy=-43.3610, ostart="", oend="", ostep=""):
    """Get energy above hull for a composition
        Args:
            comp <str>: Composition in string form
            energy <float>: Energy PER FORMULA UNIT of composition given
            (Leave the following arguments blank for a non-grand potential
                phase diagram.)
            ostart <float>: Starting oxygen chemical potential. 
            oend <float>: Ending oxygen chemical potential. 
            ostep <float>: Step for oxygen chemical potential
        Returns:
            Prints to screen
    """        
    #a = MPRester("<YOUR_MPREST_API_KEY_HERE>")
    a = MPRester("wfmUu5VSsDCvIrhz")
    
    mycomp=Composition(comp)
    print "Composition: ", mycomp
    myenergy=energy
    print "Energy: ", myenergy
    myPDEntry = PDEntry(mycomp, myenergy)

    elements = mycomp.elements
    ellist = map(str, elements)
    
    chemsys_entries = a.get_entries_in_chemsys(ellist)
    #For reference: other ways of getting entries
    #entries = a.mpquery(criteria={'elements':{'$in':['La','Mn'],'$all':['O']},'nelements':3})
    #entries = a.mpquery(criteria={'elements':{'$in':['La','Mn','O'],'$all':['O']}},properties=['pretty_formula'])
    #entries = a.get_entries_in_chemsys(['La', 'Mn', 'O', 'Sr'])
   
    if ostart=="": #Regular phase diagram
        entries = list(chemsys_entries)
        entries.append(myPDEntry)
        pd = PhaseDiagram(entries)
        #plotter = PDPlotter(gppd)
        #plotter.show()
        ppda = PDAnalyzer(pd)
        eabove=ppda.get_decomp_and_e_above_hull(myPDEntry)
        print "Energy above hull: ", eabove[1]
        print "Decomposition: ", eabove[0]
        return eabove
    else: #Grand potential phase diagram
        orange = np.arange(ostart, oend+ostep, ostep) #add ostep because otherwise the range ends before oend
        for o_chem_pot in orange:
            entries = list(chemsys_entries)
            myGrandPDEntry = GrandPotPDEntry(myPDEntry,{Element('O'): float(o_chem_pot)}) #need grand pot pd entry for GPPD
            entries.append(myGrandPDEntry)
            gppd = GrandPotentialPhaseDiagram(entries,{Element('O'): float(o_chem_pot)})
            gppda = PDAnalyzer(gppd)
            geabove=gppda.get_decomp_and_e_above_hull(myGrandPDEntry, True)
            print "******** Decomposition for mu_O = %s eV ********" % o_chem_pot
            print "%30s%1.4f" % ("mu_O: ",o_chem_pot)
            print "%30s%1.4f" % ("Energy above hull (eV): ",geabove[1])
            decomp=geabove[0]
            #print "Decomp: ", decomp
            print "%30s" % "Decomposition: "
            for dkey in decomp.keys():
                print "%30s:%1.4f" % (dkey.composition,decomp[dkey])
    return

class  PDAnalyzerTest(unittest.TestCase):

    def setUp(self):
        module_dir = os.path.dirname(os.path.abspath(__file__))
        (elements, entries) = PDEntryIO.from_csv(os.path.join(module_dir,
                                                              "pdentries_test.csv"))
        self.pd = PhaseDiagram(entries)
        self.analyzer = PDAnalyzer(self.pd)

    def test_get_e_above_hull(self):
        for entry in self.pd.stable_entries:
            self.assertLess(self.analyzer.get_e_above_hull(entry), 1e-11,
                            "Stable entries should have e above hull of zero!")
        for entry in self.pd.all_entries:
            if entry not in self.pd.stable_entries:
                e_ah = self.analyzer.get_e_above_hull(entry)
                self.assertGreaterEqual(e_ah, 0)
                self.assertTrue(isinstance(e_ah, Number))

    def test_get_equilibrium_reaction_energy(self):
        for entry in self.pd.stable_entries:
            self.assertLessEqual(
                self.analyzer.get_equilibrium_reaction_energy(entry), 0,
                "Stable entries should have negative equilibrium reaction energy!")

    def test_get_decomposition(self):
        for entry in self.pd.stable_entries:
            self.assertEquals(len(self.analyzer.get_decomposition(entry.composition)), 1,
                              "Stable composition should have only 1 decomposition!")
        dim = len(self.pd.elements)
        for entry in self.pd.all_entries:
            ndecomp = len(self.analyzer.get_decomposition(entry.composition))
            self.assertTrue(ndecomp > 0 and ndecomp <= dim,
                            "The number of decomposition phases can at most be equal to the number of components.")

        #Just to test decomp for a ficitious composition
        ansdict = {entry.composition.formula: amt
                   for entry, amt in
                   self.analyzer.get_decomposition(Composition("Li3Fe7O11")).items()}
        expected_ans = {"Fe2 O2": 0.0952380952380949,
                        "Li1 Fe1 O2": 0.5714285714285714,
                        "Fe6 O8": 0.33333333333333393}
        for k, v in expected_ans.items():
            self.assertAlmostEqual(ansdict[k], v)

    def test_get_transition_chempots(self):
        for el in self.pd.elements:
            self.assertLessEqual(len(self.analyzer.get_transition_chempots(el)),
                                 len(self.pd.facets))

    def test_get_element_profile(self):
        for el in self.pd.elements:
            for entry in self.pd.stable_entries:
                if not (entry.composition.is_element):
                    self.assertLessEqual(len(self.analyzer.get_element_profile(el, entry.composition)),
                                         len(self.pd.facets))

    def test_get_get_chempot_range_map(self):
        elements = [el for el in self.pd.elements if el.symbol != "Fe"]
        self.assertEqual(len(self.analyzer.get_chempot_range_map(elements)), 10)

    def test_dim1(self):
        #Ensure that dim 1 PDs can eb generated.
        for el in ["Li", "Fe", "O2"]:
            entries = [e for e in self.entries
                       if e.composition.reduced_formula == el]
            pd = PhaseDiagram(entries)
            self.assertEqual(len(pd.stable_entries), 1)

            a = PDAnalyzer(pd)
            for e in entries:
                decomp, ehull = a.get_decomp_and_e_above_hull(e)
                self.assertGreaterEqual(ehull, 0)
            plotter = PDPlotter(pd)
            lines, stable_entries, unstable_entries = plotter.pd_plot_data
            self.assertEqual(lines[0][1], [0, 0])

    def from_composition_and_pd(comp, pd, working_ion_symbol="Li"):
        """
        Convenience constructor to make a ConversionElectrode from a
        composition and a phase diagram.

        Args:
            comp:
                Starting composition for ConversionElectrode, e.g.,
                Composition("FeF3")
            pd:
                A PhaseDiagram of the relevant system (e.g., Li-Fe-F)
            working_ion_symbol:
                Element symbol of working ion. Defaults to Li.
        """
        working_ion = Element(working_ion_symbol)
        entry = None
        working_ion_entry = None
        for e in pd.stable_entries:
            if e.composition.reduced_formula == comp.reduced_formula:
                entry = e
            elif e.is_element and \
                    e.composition.reduced_formula == working_ion_symbol:
                working_ion_entry = e

        if not entry:
            raise ValueError("Not stable compound found at composition {}."
                             .format(comp))

        analyzer = PDAnalyzer(pd)

        profile = analyzer.get_element_profile(working_ion, comp)
        # Need to reverse because voltage goes form most charged to most
        # discharged.
        profile.reverse()
        if len(profile) < 2:
            return None
        working_ion_entry = working_ion_entry
        working_ion = working_ion_entry.composition.elements[0].symbol
        normalization_els = {}
        for el, amt in comp.items():
            if el != Element(working_ion):
                normalization_els[el] = amt
        vpairs = [ConversionVoltagePair.from_steps(profile[i], profile[i + 1],
                                                   normalization_els)
                  for i in range(len(profile) - 1)]
        return ConversionElectrode(vpairs, working_ion_entry, comp)

    def extract_phase_diagram_info(self,MP_phase_diagram_json_data_filename):

        computed_entries  = self._extract_MP_data(MP_phase_diagram_json_data_filename)
        processed_entries = self.compat.process_entries(computed_entries)

        pd = PhaseDiagram(processed_entries)
        self.phase_diagram_analyser = PDAnalyzer(pd)

        return

def get_decomp(o_chem_pot, mycomp, verbose=1):
    """Get decomposition from open phase diagram
        Args:
            o_chem_pot <float>: Oxygen chemical potential
            mycomp <pymatgen Composition>: Composition
            verbose <int>: 1 - verbose (default)
                           0 - silent
        Returns:
            decomposition string
    """        
    a = MPRester("<YOUR_MPREST_API_KEY_HERE>")
    elements = mycomp.elements
    ellist = map(str, elements)
    entries = a.get_entries_in_chemsys(ellist)
    #entries = a.get_entries_in_chemsys(['La', 'Mn', 'O', 'Fe'])
    pd = PhaseDiagram(entries)
    gppd = GrandPotentialPhaseDiagram(entries,{Element('O'): float(o_chem_pot)})
    print gppd
    #plotter = PDPlotter(gppd)
    #plotter.show()

    gppda = PDAnalyzer(gppd)
    #mychempots = gppda.get_composition_chempots(mycomp)
    #print "My chem pots:"
    #print mychempots
    mydecompgppd = gppda.get_decomposition(mycomp)
    #pdentry = PDEntry(mycomp, 0)
    #print "Decomp and energy:"
    #decompandenergy = gppda.get_decomp_and_e_above_hull(pdentry)
    #print decompandenergy
    #mydecomppd = pda.get_decomposition(mycomp)
    #print "Mn profile:"
    #mnprof= gppda.get_element_profile(Element('Mn'),mycomp)
    #print mnprof

    if verbose:
        for (entry,amount) in mydecompgppd.iteritems():
            print "%s: %3.3f" % (entry.name, amount)
            #mymurangegppd = gppda.getmu_range_stability_phase(Composition(entry.name),Element('O'))
            #print mymurangegppd
        #for (entry,amount) in mydecomppd.iteritems():
        #    print "%s: %3.3f" % (entry.name, amount)
        print ""
    return mydecompgppd

    def get_contour_pd_plot(self):
        """
        Plot a contour phase diagram plot, where phase triangles are colored
        according to degree of instability by interpolation. Currently only
        works for 3-component phase diagrams.

        Returns:
            A matplotlib plot object.
        """
        from scipy import interpolate
        from matplotlib import cm

        pd = self._pd
        entries = pd.qhull_entries
        data = np.array(pd.qhull_data)

        plt = self._get_2d_plot()
        analyzer = PDAnalyzer(pd)
        data[:, 0:2] = triangular_coord(data[:, 0:2]).transpose()
        for i, e in enumerate(entries):
            data[i, 2] = analyzer.get_e_above_hull(e)

        gridsize = 0.005
        xnew = np.arange(0, 1.0, gridsize)
        ynew = np.arange(0, 1, gridsize)

        f = interpolate.LinearNDInterpolator(data[:, 0:2], data[:, 2])
        znew = np.zeros((len(ynew), len(xnew)))
        for (i, xval) in enumerate(xnew):
            for (j, yval) in enumerate(ynew):
                znew[j, i] = f(xval, yval)

        plt.contourf(xnew, ynew, znew, 1000, cmap=cm.autumn_r)

        plt.colorbar()
        return plt

 def setUp(self):
     module_dir = os.path.dirname(os.path.abspath(__file__))
     (elements, entries) = PDEntryIO.from_csv(os.path.join(module_dir,
                                                           "pdentries_test.csv"))
     self.pd = PhaseDiagram(entries)
     self.analyzer = PDAnalyzer(self.pd)

class PDAnalyzerTest(unittest.TestCase):

    def setUp(self):
        module_dir = os.path.dirname(os.path.abspath(__file__))
        (elements, entries) = PDEntryIO.from_csv(os.path.join(module_dir,
                                                              "pdentries_test.csv"))
        self.pd = PhaseDiagram(entries)
        self.analyzer = PDAnalyzer(self.pd)

    def test_get_e_above_hull(self):
        for entry in self.pd.stable_entries:
            self.assertLess(self.analyzer.get_e_above_hull(entry), 1e-11,
                            "Stable entries should have e above hull of zero!")
        for entry in self.pd.all_entries:
            if entry not in self.pd.stable_entries:
                e_ah = self.analyzer.get_e_above_hull(entry)
                self.assertGreaterEqual(e_ah, 0)
                self.assertTrue(isinstance(e_ah, Number))

    def test_get_equilibrium_reaction_energy(self):
        for entry in self.pd.stable_entries:
            self.assertLessEqual(
                self.analyzer.get_equilibrium_reaction_energy(entry), 0,
                "Stable entries should have negative equilibrium reaction energy!")

    def test_get_decomposition(self):
        for entry in self.pd.stable_entries:
            self.assertEquals(len(self.analyzer.get_decomposition(entry.composition)), 1,
                              "Stable composition should have only 1 decomposition!")
        dim = len(self.pd.elements)
        for entry in self.pd.all_entries:
            ndecomp = len(self.analyzer.get_decomposition(entry.composition))
            self.assertTrue(ndecomp > 0 and ndecomp <= dim,
                            "The number of decomposition phases can at most be equal to the number of components.")

        #Just to test decomp for a ficitious composition
        ansdict = {entry.composition.formula: amt
                   for entry, amt in
                   self.analyzer.get_decomposition(Composition("Li3Fe7O11")).items()}
        expected_ans = {"Fe2 O2": 0.0952380952380949,
                        "Li1 Fe1 O2": 0.5714285714285714,
                        "Fe6 O8": 0.33333333333333393}
        for k, v in expected_ans.items():
            self.assertAlmostEqual(ansdict[k], v)

    def test_get_transition_chempots(self):
        for el in self.pd.elements:
            self.assertLessEqual(len(self.analyzer.get_transition_chempots(el)),
                                 len(self.pd.facets))

    def test_get_element_profile(self):
        for el in self.pd.elements:
            for entry in self.pd.stable_entries:
                if not (entry.composition.is_element):
                    self.assertLessEqual(len(self.analyzer.get_element_profile(el, entry.composition)),
                                         len(self.pd.facets))

    def test_get_get_chempot_range_map(self):
        elements = [el for el in self.pd.elements if el.symbol != "Fe"]
        self.assertEqual(len(self.analyzer.get_chempot_range_map(elements)), 10)

    def test_getmu_vertices_stability_phase(self):
        results = self.analyzer.getmu_vertices_stability_phase(Composition.from_formula("LiFeO2"), Element("O"))
        self.assertAlmostEqual(len(results), 6)
        test_equality = False
        for c in results:
            if abs(c[Element("O")]+7.115) < 1e-2 and abs(c[Element("Fe")]+6.596) < 1e-2 and \
                    abs(c[Element("Li")]+3.931) < 1e-2:
                test_equality = True
        self.assertTrue(test_equality,"there is an expected vertex missing in the list")


    def test_getmu_range_stability_phase(self):
        results = self.analyzer.get_chempot_range_stability_phase(
            Composition("LiFeO2"), Element("O"))
        self.assertAlmostEqual(results[Element("O")][1], -4.4501812249999997)
        self.assertAlmostEqual(results[Element("Fe")][0], -6.5961470999999996)
        self.assertAlmostEqual(results[Element("Li")][0], -3.6250022625000007)

            "Project parameters."
        sys.exit()

    syms = [el.symbol for el in entry.composition.elements]
    #This gets all entries belonging to the relevant system.
    entries = a.get_entries_in_chemsys(syms)
    entries.append(entry)

    #Process entries with Materials Project compatibility.
    entries = compat.process_entries(entries)

    print [e.composition.reduced_formula for e in entries]

    pd = PhaseDiagram(entries)

    analyzer = PDAnalyzer(pd)
    ehull = analyzer.get_e_above_hull(entry) * 1000

    print "Run contains formula {} with corrected energy {:.3f} eV.".format(
        entry.composition, entry.energy
    )
    print "Energy above convex hull = {:.1f} meV".format(ehull)
    if ehull < 1:
        print "Entry is stable."
    elif ehull < 30:
        print "Entry is metastable and could be stable at finite temperatures."
    elif ehull < 50:
        print "Entry has a low probability of being stable."
    else:
        print "Entry is very unlikely to be stable."

    def get_chempot_range_map_plot(self, elements):
        """
        Returns a plot of the chemical potential range map. Currently works
        only for 3-component PDs.

        Args:
            elements:
                Sequence of elements to be considered as independent variables.
                E.g., if you want to show the stability ranges of all Li-Co-O
                phases wrt to uLi and uO, you will supply
                [Element("Li"), Element("O")]
        Returns:
            A matplotlib plot object.
        """

        plt = get_publication_quality_plot(12, 8)
        analyzer = PDAnalyzer(self._pd)
        chempot_ranges = analyzer.get_chempot_range_map(elements)
        missing_lines = {}
        excluded_region = []
        for entry, lines in chempot_ranges.items():
            comp = entry.composition
            center_x = 0
            center_y = 0
            coords = []
            contain_zero = any([comp.get_atomic_fraction(el) == 0 for el in elements])
            is_boundary = (not contain_zero) and sum([comp.get_atomic_fraction(el) for el in elements]) == 1
            for line in lines:
                (x, y) = line.coords.transpose()
                plt.plot(x, y, "k-")

                for coord in line.coords:
                    if not in_coord_list(coords, coord):
                        coords.append(coord.tolist())
                        center_x += coord[0]
                        center_y += coord[1]
                if is_boundary:
                    excluded_region.extend(line.coords)

            if coords and contain_zero:
                missing_lines[entry] = coords
            else:
                xy = (center_x / len(coords), center_y / len(coords))
                plt.annotate(latexify(entry.name), xy, fontsize=22)

        ax = plt.gca()
        xlim = ax.get_xlim()
        ylim = ax.get_ylim()

        # Shade the forbidden chemical potential regions.
        excluded_region.append([xlim[1], ylim[1]])
        excluded_region = sorted(excluded_region, key=lambda c: c[0])
        (x, y) = np.transpose(excluded_region)
        plt.fill(x, y, "0.80")

        # The hull does not generate the missing horizontal and vertical lines.
        # The following code fixes this.
        el0 = elements[0]
        el1 = elements[1]
        for entry, coords in missing_lines.items():
            center_x = sum([c[0] for c in coords])
            center_y = sum([c[1] for c in coords])
            comp = entry.composition
            is_x = comp.get_atomic_fraction(el0) < 0.01
            is_y = comp.get_atomic_fraction(el1) < 0.01
            n = len(coords)
            if not (is_x and is_y):
                if is_x:
                    coords = sorted(coords, key=lambda c: c[1])
                    for i in [0, -1]:
                        x = [min(xlim), coords[i][0]]
                        y = [coords[i][1], coords[i][1]]
                        plt.plot(x, y, "k")
                        center_x += min(xlim)
                        center_y += coords[i][1]
                elif is_y:
                    coords = sorted(coords, key=lambda c: c[0])
                    for i in [0, -1]:
                        x = [coords[i][0], coords[i][0]]
                        y = [coords[i][1], min(ylim)]
                        plt.plot(x, y, "k")
                        center_x += coords[i][0]
                        center_y += min(ylim)
                xy = (center_x / (n + 2), center_y / (n + 2))
            else:
                center_x = sum(coord[0] for coord in coords) + xlim[0]
                center_y = sum(coord[1] for coord in coords) + ylim[0]
                xy = (center_x / (n + 1), center_y / (n + 1))

            plt.annotate(
                latexify(entry.name), xy, horizontalalignment="center", verticalalignment="center", fontsize=22
            )

        plt.xlabel("$\mu_{{{0}}} - \mu_{{{0}}}^0$ (eV)".format(el0.symbol))
        plt.ylabel("$\mu_{{{0}}} - \mu_{{{0}}}^0$ (eV)".format(el1.symbol))
        plt.tight_layout()
        return plt

    def _get_2d_plot(self, label_stable=True, label_unstable=True,
                     ordering=None, energy_colormap=None, vmin_mev=-60.0,
                     vmax_mev=60.0, show_colorbar=True,
                     process_attributes=False):
        """
        Shows the plot using pylab.  Usually I won't do imports in methods,
        but since plotting is a fairly expensive library to load and not all
        machines have matplotlib installed, I have done it this way.
        """

        plt = get_publication_quality_plot(8, 6)
        from matplotlib.font_manager import FontProperties
        if ordering is None:
            (lines, labels, unstable) = self.pd_plot_data
        else:
            (_lines, _labels, _unstable) = self.pd_plot_data
            (lines, labels, unstable) = order_phase_diagram(
                _lines, _labels, _unstable, ordering)
        if energy_colormap is None:
            if process_attributes:
                for x, y in lines:
                    plt.plot(x, y, "k-", linewidth=3, markeredgecolor="k")
                # One should think about a clever way to have "complex"
                # attributes with complex processing options but with a clear
                #  logic. At this moment, I just use the attributes to know
                # whether an entry is a new compound or an existing (from the
                #  ICSD or from the MP) one.
                for x, y in labels.keys():
                    if labels[(x, y)].attribute is None or \
                            labels[(x, y)].attribute == "existing":
                        plt.plot(x, y, "ko", linewidth=3, markeredgecolor="k",
                                 markerfacecolor="b", markersize=12)
                    else:
                        plt.plot(x, y, "k*", linewidth=3, markeredgecolor="k",
                                 markerfacecolor="g", markersize=18)
            else:
                for x, y in lines:
                    plt.plot(x, y, "ko-", linewidth=3, markeredgecolor="k",
                             markerfacecolor="b", markersize=15)
        else:
            from matplotlib.colors import Normalize, LinearSegmentedColormap
            from matplotlib.cm import ScalarMappable
            pda = PDAnalyzer(self._pd)
            for x, y in lines:
                plt.plot(x, y, "k-", linewidth=3, markeredgecolor="k")
            vmin = vmin_mev / 1000.0
            vmax = vmax_mev / 1000.0
            if energy_colormap == 'default':
                mid = - vmin / (vmax - vmin)
                cmap = LinearSegmentedColormap.from_list(
                    'my_colormap', [(0.0, '#005500'), (mid, '#55FF55'),
                                    (mid, '#FFAAAA'), (1.0, '#FF0000')])
            else:
                cmap = energy_colormap
            norm = Normalize(vmin=vmin, vmax=vmax)
            _map = ScalarMappable(norm=norm, cmap=cmap)
            _energies = [pda.get_equilibrium_reaction_energy(entry)
                         for coord, entry in labels.items()]
            energies = [en if en < 0.0 else -0.00000001 for en in _energies]
            vals_stable = _map.to_rgba(energies)
            ii = 0
            if process_attributes:
                for x, y in labels.keys():
                    if labels[(x, y)].attribute is None or \
                            labels[(x, y)].attribute == "existing":
                        plt.plot(x, y, "o", markerfacecolor=vals_stable[ii],
                                 markersize=12)
                    else:
                        plt.plot(x, y, "*", markerfacecolor=vals_stable[ii],
                                 markersize=18)
                    ii += 1
            else:
                for x, y in labels.keys():
                    plt.plot(x, y, "o", markerfacecolor=vals_stable[ii],
                             markersize=15)
                    ii += 1

        font = FontProperties()
        font.set_weight("bold")
        font.set_size(24)

        # Sets a nice layout depending on the type of PD. Also defines a
        # "center" for the PD, which then allows the annotations to be spread
        # out in a nice manner.
        if len(self._pd.elements) == 3:
            plt.axis("equal")
            plt.xlim((-0.1, 1.2))
            plt.ylim((-0.1, 1.0))
            plt.axis("off")
            center = (0.5, math.sqrt(3) / 6)
        else:
            all_coords = labels.keys()
            miny = min([c[1] for c in all_coords])
            ybuffer = max(abs(miny) * 0.1, 0.1)
            plt.xlim((-0.1, 1.1))
            plt.ylim((miny - ybuffer, ybuffer))
            center = (0.5, miny / 2)
            plt.xlabel("Fraction", fontsize=28, fontweight='bold')
            plt.ylabel("Formation energy (eV/fu)", fontsize=28,
                       fontweight='bold')
#.........这里部分代码省略.........

class AnalyseMaterialsProjectJsonDataWithComputedEntries():
    """
    Class which will wrap around boilerplate analysis of MaterialsProject-like
    json files, containing data extracted using borgs and queens.

    It will be assumed that we are providing ComputedEntries objects directly.
    """

    def __init__(self):
        # some MP analysis power tools
        self.compat  = MaterialsProjectCompatibility()

        return

    def extract_alkali_energy(self, computed_Alkali_entry ):
        processed_Alkali_entry = self.compat.process_entry(computed_Alkali_entry)
        self.E_Alkali = processed_Alkali_entry.energy
        return

    def extract_phase_diagram_info(self,MP_phase_diagram_json_data_filename):

        computed_entries  = self._extract_MP_data(MP_phase_diagram_json_data_filename)
        processed_entries = self.compat.process_entries(computed_entries)

        pd = PhaseDiagram(processed_entries)
        self.phase_diagram_analyser = PDAnalyzer(pd)

        return

    def extract_processed_entries(self,computed_entries):
        processed_entries = self.compat.process_entries(computed_entries)

        return processed_entries

    def extract_energies_above_hull(self,computed_entries,alkali):

        processed_entries = self.extract_processed_entries(computed_entries)

        list_energy_above_hull  = []
        list_alkali_content = []

        for entry in processed_entries: 
            decomposition_dict, energy_above_hull  = \
                self.phase_diagram_analyser.get_decomp_and_e_above_hull(entry, allow_negative=True)

            list_energy_above_hull.append(energy_above_hull)  
            list_alkali_content.append(entry.composition[alkali])

        list_energy_above_hull  = np.array(list_energy_above_hull)
        list_alkali_content     = np.array(list_alkali_content )

        return list_alkali_content, list_energy_above_hull  

    def extract_energies(self,computed_entries,alkali):

        processed_entries = self.extract_processed_entries(computed_entries)

        list_energy         = []
        list_alkali_content = []
        for entry in processed_entries:
            list_energy.append(entry.energy)
            list_alkali_content.append(entry.composition[alkali])

        list_energy         = np.array(list_energy)
        list_alkali_content = np.array(list_alkali_content )

        I = np.argsort(list_alkali_content )
        
        return list_alkali_content[I], list_energy[I]

    def _extract_MP_data(self,MP_data_filename):

        drone = VaspToComputedEntryDrone()
        queen = BorgQueen(drone, "dummy", 1)

        queen.load_data(MP_data_filename)
        computed_entries = queen.get_data()

        del drone
        del queen

        return computed_entries