Plot Phasediagram

This python script use ase to plot phasediagram

#lipai@mail.ustc.edu.cn
from ase.phasediagram import PhaseDiagram
u=[0.0,-2.92,-3.10,-3.20]

Li=0
Br=0.0
O2=0.0

Br2O3=-0.209*5

for u_Li in u:
      Li2O=-2.071*3-2*u_Li
      Li2O2=-1.65*4-2*u_Li

      LiBr=-1.573*2-u_Li

      #Li3OBr=Li2O+LiBr+0.075*5
      Li32O11Br10=11*Li2O+10*LiBr+0.059*53

      refs = [('Br',Br),
              ('O2',O2),
              ('Li',u_Li),
              ('Br2O3',Br2O3),
              ('Li2O',Li2O),
              ('Li2O2',Li2O2),
              ('LiBr',LiBr),
      #        ('Li3OBr',Li3OBr),
              ('Li32O11Br10',Li32O11Br10)]

      pd = PhaseDiagram(refs)
      fig=pd.plot(show=False)
      fig.savefig("try"+str(u_Li)+".png")
      #pd.savefig(str(u_Li)+".png")

import matplotlib.pyplot as plt
u=[0.0,-2.92,-3.10,-3.20]

for u_Li in u:
      O2=0.0
      Li2O=(-2.071*3-2*u_Li)/1
      Li2O2=(-1.65*4-2*u_Li)/2

      Br2O3=(-0.209*5)/5

      Br=0.0
      LiBr=(-1.573*2-u_Li)/1
      Li32O11Br10=(11*Li2O+10*LiBr+0.059*53)/21

      x_Br=[1,0,2/5.0,0,0,1,10/21.0]
      y=[Br,O2,Br2O3,Li2O,Li2O2,LiBr,Li32O11Br10]
      name=["Br","O2","Br2O3","Li2O","Li2O2","LiBr","Li32O11Br10"]
      plt.plot(x_Br,y,'sr')
      for i,j,k in zip(x_Br,y,name):
            #cname=
            if k=="Li2O" :
                  plt.text(i-0.09,j+0.01,k,color='g')
            else:
                  plt.text(i+0.02,j+0.01,k,color='r')
      plt.xlim([-0.1,1.1])
      plt.axis("off")
      plt.savefig(str(u_Li)+".png")
      plt.clf()

      print "x_Br     y"
      for i,j in zip(x_Br,y):
            print i,j

# coding: utf-8

# In[ ]:

from pymatgen import MPRester, Composition, Element
from pymatgen.io.vasp.outputs import Vasprun
from pymatgen.entries.compatibility import MaterialsProjectCompatibility
#from pymatgen.phasediagram.maker import PhaseDiagram
from pymatgen.analysis.phase_diagram import PhaseDiagram, PDPlotter, PDAnalyzer
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
#from pymatgen.analysis.phase_diagram.analyzer import PDAnalyzer
#from pymatgen.phasediagram.plotter import PDPlotter
#from pymatgen.util.plotting_utils import get_publication_quality_plot
import json
import re
import palettable
import matplotlib as mpl


v = Vasprun("./vasprun.xml")
#entry = v.get_computed_entry(inc_structure=True)
entry = v.get_computed_entry()

m = MPRester("g8u7XE1HNSpaI2c3")
entries = m.get_entries_in_chemsys(["Li", "O","Br"]) # pass in a list of elements in your material


#compatibility = MaterialsProjectCompatibility()
#entry = compatibility.process_entry(entry)
#allentries = compatibility.process_entries([entry] + entries)
#comp = MaterialsProjectCompatibility()

#entry = comp.process_entry(entry)
#all_entries = comp.process_entries(entries + [entry])
#allentries = entries
allentries = [entry] + entries


pd = PhaseDiagram(allentries)
#plotter = PDPlotter(pd,)
plotter = PDPlotter(pd,show_unstable=True)

plotter.write_image("Figure_2.svg", image_format="svg", label_stable=True,
                    label_unstable=True, ordering=None,
                    energy_colormap=None, process_attributes=False)

#plotter.show()
# I suppose you need the energy above hull for your material
pda = PDAnalyzer(pd)
print("Energy above hull is %.2f eV/atom" % pda.get_e_above_hull(entry))
analyzer = PDAnalyzer(pd)
ehull = analyzer.get_e_above_hull(entry)
print("The energy above hull of * is %.3f eV/atom." % ehull)


li_entries = [e for e in allentries if e.composition.reduced_formula == "Li"]
uli0 = min(li_entries, key=lambda e: e.energy_per_atom).energy_per_atom

el_profile = analyzer.get_element_profile(Element("Li"), entry.composition)
for i, d in enumerate(el_profile):
    voltage = -(d["chempot"] - uli0)
    print("Voltage: %s V" % voltage)
    print(d["reaction"])
    print("")


# plot
# Some matplotlib settings to improve the look of the plot.
mpl.rcParams['axes.linewidth']=3
mpl.rcParams['lines.markeredgewidth']=4
mpl.rcParams['lines.linewidth']=3
mpl.rcParams['lines.markersize']=15
mpl.rcParams['xtick.major.width']=3
mpl.rcParams['xtick.major.size']=8
mpl.rcParams['xtick.minor.width']=3
mpl.rcParams['xtick.minor.size']=4
mpl.rcParams['ytick.major.width']=3
mpl.rcParams['ytick.major.size']=8
mpl.rcParams['ytick.minor.width']=3
mpl.rcParams['ytick.minor.size']=4


# Plot of Li uptake per formula unit (f.u.) of Li6PS5Cl against voltage vs Li/Li+.
colors = palettable.colorbrewer.qualitative.Set1_9.mpl_colors
plt = get_publication_quality_plot(12, 8)

for i, d in enumerate(el_profile):
    v = - (d["chempot"] - uli0)
    if i != 0:
        plt.plot([x2, x2], [y1, d["evolution"] / 4.0], 'k', linewidth=3)
    x1 = v
    y1 = d["evolution"] / 4.0
    if i != len(el_profile) - 1:
        x2 = - (el_profile[i + 1]["chempot"] - uli0)
    else:
        x2 = 5.0

    if i in [0, 4, 5, 7]:
        products = [re.sub(r"(\d+)", r"$_{\1}$", p.reduced_formula)
                    for p in d["reaction"].products if p.reduced_formula != "Li"]

        plt.annotate(", ".join(products), xy=(v + 0.05, y1 + 0.3),
                     fontsize=24, color=colors[0])

        plt.plot([x1, x2], [y1, y1], color=colors[0], linewidth=5)
    else:
        plt.plot([x1, x2], [y1, y1], 'k', linewidth=3)

plt.xlim((0, 4.0))
plt.ylim((-6, 10))
plt.xlabel("Voltage vs Li/Li$^+$ (V)")
plt.ylabel("Li uptake per f.u.")
plt.tight_layout()
plt.savefig('Figure_V.png')