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')