In [1]:

    

"""
This is slightly modified version of the notebook in the pymatgen/examples folder
"""
import numpy as np

from pymatgen.matproj.rest import MPRester
from pymatgen.core.ion import Ion
from pymatgen import Element

from pymatgen.phasediagram.pdmaker import PhaseDiagram
from pymatgen.analysis.pourbaix.entry import PourbaixEntry, IonEntry
from pymatgen.analysis.pourbaix.maker import PourbaixDiagram
from pymatgen.analysis.pourbaix.plotter import PourbaixPlotter
from pymatgen.entries.compatibility import MaterialsProjectCompatibility, AqueousCorrection

%matplotlib inline


    

In [2]:

    

def contains_entry(entry_list, entry):
    for e in entry_list:
        if e.entry_id == entry.entry_id or \
            (abs(entry.energy_per_atom 
                - e.energy_per_atom) < 1e-6 and
                    entry.composition.reduced_formula ==
             e.composition.reduced_formula):
            return True


    

In [3]:

    

#This initializes the REST adaptor. Put your own API key in.
a = MPRester("dwvz2XCFUEI9fJiR")


    

In [4]:

    

#Entries are the basic unit for thermodynamic and other analyses in pymatgen.
#This gets all entries belonging to the Fe-O-H system.
entries = a.get_entries_in_chemsys(['Fe', 'O', 'H'])


    

In [5]:

    

#exp_entries = a.get_exp_entry("SnO")
#print exp_entries


    

In [6]:

    

#Dictionary of ion:energy, where the energy is the formation energy of ions from 
#the NBS tables. (Source: NBS Thermochemical Tables; FeO4[2-]: Misawa T., Corr. Sci., 13(9), 659-676 (1973))
ion_dict = {"Fe[2+]":-0.817471, "Fe[3+]":-0.0478, "FeO2[2-]":-3.06055, "FeOH[+]":-2.8738,
            "FeOH[2+]":-2.37954, "HFeO2[-]":-3.91578, "Fe(OH)2[+]":-4.54022, "Fe2(OH)2[4+]":-4.84285,
            "FeO2[-]":-3.81653, "FeO4[2-]":-3.33946, "Fe(OH)3(aq)":-6.83418, "Fe(OH)2[+]":-4.54022}


    

In [7]:

    

elements = set()
data = []
for entry in entries:
    elements.update(entry.composition.elements)
#print elements
for entry in entries:
    comp = entry.composition
    #print comp
    row = [comp.get_atomic_fraction(el) for el in elements]
    row.append(entry.energy_per_atom)
    data.append(row)
data = np.array(data)
#print data
#print data[:, 1:]


    

In [8]:

    

# Run aqueouscorrection on the entries
aqcompat = AqueousCorrection("/home/matk/Software/pymatgen/pymatgen/entries/MPCompatibility.yaml")
entries_aqcorr = list()

for entry in entries:
    aq_corrected_entry = aqcompat.correct_entry(entry)
    if not contains_entry(entries_aqcorr, aq_corrected_entry):
        entries_aqcorr.append(aq_corrected_entry)

# Generate a phase diagram to consider only solid entries stable in water. 
#for e in entries_aqcorr:
#    print e.name
pd = PhaseDiagram(entries_aqcorr)
#for e in pd.qhull_entries:
#    print e.name
#for el in pd.el_refs:
#    print el.name
#print pd.facets
#print pd.qhull_data


    

In [9]:

    

stable_solids = pd.stable_entries
for ss in stable_solids:
    print ss.name
    print ss.composition.get_reduced_composition_and_factor()#[1]
    
stable_solids_minus_h2o = [entry for entry in stable_solids if 
                           entry.composition.reduced_formula not in ["H2", "O2", "H2O", "H2O2"]]

pbx_solid_entries = []
for entry in stable_solids_minus_h2o:
    pbx_entry = PourbaixEntry(entry)
    pbx_entry.g0_replace(pd.get_form_energy(entry))
    pbx_entry.reduced_entry()
    pbx_solid_entries.append(pbx_entry)


    

    




Fe
(Comp: Fe1, 1.0)
O2
(Comp: O2, 4.0)
Fe3O4
(Comp: Fe3 O4, 2.0)
Fe2O3
(Comp: Fe2 O3, 8.0)
H2O
(Comp: H2 O1, 12.0)
H2
(Comp: H2, 1.0)
FeO
(Comp: Fe1 O1, 2.0)

In [10]:

    

#Dictionary of reference state:experimental formation energy (from O. Kubaschewski) for reference state.
ref_state = "Fe2O3"
ref_dict = {"Fe2O3": -7.685050670886141}

# Calculate DFT reference energy for ions (See Persson et al, PRB (2012))
ref_entry = [entry for entry in stable_solids_minus_h2o if entry.composition.reduced_formula == ref_state][0]
ion_correction = pd.get_form_energy(ref_entry)/ref_entry.composition.get_reduced_composition_and_factor()[1] - ref_dict[ref_state]
 
el = Element("Fe")
pbx_ion_entries = []

# Get PourbaixEntry corresponding to each ion
for key in ion_dict:
    comp = Ion.from_formula(key)
    factor = comp.composition[el] / (ref_entry.composition[el] / ref_entry.composition.get_reduced_composition_and_factor()[1])
    #experimental energy + correction
    energy = ion_dict[key] + ion_correction * factor
    pbx_entry_ion = PourbaixEntry(IonEntry(comp, energy))
    pbx_entry_ion.name = key
    pbx_ion_entries.append(pbx_entry_ion)


    

In [11]:

    

# Generate and plot Pourbaix diagram
#Each bulk solid/ion has a free energy g of the form:
#g = g0_ref + 0.0591 * log10(conc) - nO * mu_H2O + (nH - 2nO) * pH + phi * (-nH + 2nO + q)
all_entries = pbx_solid_entries + pbx_ion_entries
pourbaix = PourbaixDiagram(all_entries)
print pourbaix._qhull_data
plotter = PourbaixPlotter(pourbaix)
plotter.plot_pourbaix(limits=[[-2, 16],[-3, 3]])


    

    




([0.0, -2.0, -1.56866370393193], [-1.0, -2.0, -1.1666927039319304], [0.0, -3.0, -0.7989927039319302], [-1.0, -3.0, -0.6724327039319302], [-2.6666666666666665, -2.6666666666666665, -0.5833912150000033], [-2.0, -2.0, -0.5625319462500022], [-3.0, -3.0, -0.5516680393750004], [-1.0, -3.0, -0.5370177039319306], [-2.0, -3.0, -0.37481270393193], [-3.0, -3.0, -0.21047270393192896], [0.0, 0.0, 0.0], [-3.0, -2.0, 0.24962729606806988], [-4.0, -3.0, 0.3488772960680695], [-4.0, -2.0, 1.1048572960680692], [-8.0, -6.0, 5.74254729606807])






    

In [11]: