Question about Grand Potential Interfacial Reactivity

[Kimmj]

I want to make a grand potential interfacial reactivity plot of Li3PS4.
So I imported GrandPotentialInterfacialReactivity.
Here is my python code.

from pymatgen.analysis.interface_reactions import InterfacialReactivity, GrandPotentialInterfacialReactivity
from pymatgen.core.composition import Composition
from pymatgen.analysis.phase_diagram import PhaseDiagram, GrandPotentialPhaseDiagram, PDPlotter
from pymatgen.ext.matproj import MPRester

with MPRester(“API”) as m:

entries = m.get_entries_in_chemsys(["Li","P","S"], inc_structure=True)

gpd = GrandPotentialPhaseDiagram(entries,{“Li”:2.2})
GPInfReact = GrandPotentialInterfacialReactivity(Composition(“Li3PS4”), Composition(“Li”),gpd, norm=True, include_no_mixing_energy=True, pd_non_grand=pd, use_hull_energy=False)
GPInfReact.plot()

But, error occurs at line GrandPotentialInterfacialReactivity.
I don’t know why this error occur.
I think I entered the right factors.
Please tell me what the problem is and how to solve it.

[mattmcdermott]

Hi @Kimmj,

If you are using the GrandPotentialInterfacialReactivity class, you can not use the open element (in this case, Li) as one of the reactants.

If you try typing in any other composition for reactant 2 (e.g., LiS), the code you wrote should work and you will still be able to capture any open decomposition reactions of Li3PS4 – just check the endpoint and you should see possible reactions of the form Li3PS4 + Li → … depending on the chemical potential you enter.

[Kimmj]

Thank you for your reply.

[buithitham]

Dear Mr.mattmcdermott,

May I bother you or anyone here with the Interface Reaction App in Open system? I am really confused in determining the chemical potential of the open element. I did similar for LiCoO2 and Li3PS4 reactants using the MP Interface Reaction code (version 2019.2.28). I chose Li as the open-element with input mu_Li=-2. However, when I check the mu_Li in output, its value is mu_Li=-3.9.

I think the value of mu_Li is important in calculating reaction energy (Er) for open system. For example:
0.2683 Li3PS4 + 0.7317 LiCoO2 → 0.5366 Li + 0.2439 Co3S4 + 0.09756 Li2SO4 + 0.2683 Li3PO4
Er(eV/atom)= {0.5366mu_Li + 0.2439Ef[Co3S4] + 0.09756Ef[Li2SO4] + 0.2683Ef[Li3PO4] - 0.2683Ef[Li3PS4] - 0.7317Ef[LiCoO2]}/[0.2683x8 + 0.7317x4]
Even I use mu_Li=-2 or mu_Li=-3.9, I can’t get the same Er=-0.596 (eV/atom) as shown by the Interface Reaction app. Is there anything wrong with my reaction energy calculation or determining the mu_Li?
Thank you so much in advance.

Tham

[mattmcdermott]

Hi @buithitham, I am a bit confused by your question. Could you please post the code you are using?

With regard to the calculation of reaction energy, you will have to be very careful to write out the correct equation since you are working with grand potentials. It is also easy to get tripped up by energy normalization here. You may want to start with a simpler example to get a feel for it. If you write out all the values for Ef and where you’re getting them that would help.

My post here might also help (it’s about formation energy calculations; remember for grand potential it’s always divided by number of atoms that are not open, e.g. everything but Li).

How to calculate Formation Energy / Atom on material project?

Sorry @Tanmoy_Paul for the late reply. The issue of calculating formation energies correctly comes up very often, so I wanted to clear it up for anyone who stumbles upon this in the future: When calculating formation energy by hand, you have to be extremely careful about energy normalization. When you have the energy result of a DFT calculation, there are three possible energy normalization schemes you could follow: Normalize the energy to a single atom (this is the standard eV/atom you see. …

[buithitham]

Dear Mr.mattmcdermott,

Thank for your great information,
I was use te code for calculate the electrochemical inferface energy between solid elctrolyte|cathode as bellow:

from pymatgen.ext.matproj import MPRester 

from pymatgen.analysis.interface_reactions import  GrandPotentialInterfacialReactivity
from pymatgen.analysis.phase_diagram import PhaseDiagram, GrandPotentialPhaseDiagram, PDEntry, PDPlotter

from pymatgen.core import Composition, Element

%matplotlib inline
mpr = MPRester('API')

reactant1 = 'LiCoO2'
reactant2 = 'Li3PS4'
grand = True
if grand:
   open_el = 'Li' 
relative_mu = -2

comp1 = Composition(reactant1)
comp2 = Composition(reactant2)
elements = [e.symbol for e in comp1.elements + comp2.elements]
if grand:
    elements.append(open_el)
elements = list(set(elements)) 
entries = mpr.get_entries_in_chemsys(elements)
pd = PhaseDiagram(entries)
if grand:
    mu = pd.get_transition_chempots(Element(open_el))[0]
    chempots = {open_el: relative_mu + mu}
    gpd = GrandPotentialPhaseDiagram(entries, chempots)
    interface = GrandPotentialInterfacialReactivity(
        comp1, comp2, gpd, norm=True, include_no_mixing_energy=True, pd_non_grand=pd, use_hull_energy=False)
else:
    interface = GrandPotentialInterfacialReactivity(
        comp1, comp2, pd, norm=True, include_no_mixing_energy=False, pd_non_grand=None, use_hull_energy=False)
interface.plot()
from collections import OrderedDict
import pandas as pd

critical_rxns = [
    OrderedDict([
        ("Atomic fraction", round(ratio, 3)),
        ("Reaction equation", rxn),
        ("E$_{rxt}$ per mol equation (kJ/mol)", round(rxn_energy, 1)),
        ("E$_{rxt}$ per reactant atom (eV/atom)", round(reactivity, 3)),

    ])
    for _, ratio, reactivity, rxn, rxn_energy in interface.get_kinks()]
interface_reaction_table = pd.DataFrame(critical_rxns)
interface_reaction_table

I am follow this paper (First principles study on electrochemical and chemical stability of solid electrolyte–electrode interfaces in all-solid-state Li-ion batteries - Journal of Materials Chemistry A (RSC Publishing))

With mu_Li=-2 , the code print minimum energy reaction is -641 meV/atom. However, when I apply mu_Li = -2.91 the code print minimum energy reaction is -783 meV/atom, it is different value with table 5 on paper. So I can not understand the meaning of mu_Li and the chempots in the code.
How can I apply potential of charge bateries for calculate electrochemical interface reaction?
Thank you so much for helping,

Best reagrds,
Tham

[mattmcdermott]

\Delta \mu_{Li} is what you refer to in your code as relative_mu and it should relate to the Li potential as \Phi = -\Delta\mu_{Li}. To get the actual chemical potential you can do: \mu_{Li} = \mu_{Li}^0 (ref) + \Delta \mu_{Li}. Note that \mu_{Li}^0 (ref) is just the calculated energy per atom of Li metal. It is in terms of VASP energies; if this was with standard experimental formation energies you would say it is exactly 0.

Your code seems to be calculating the chemical potential correctly, although you can access the \mu_{Li}^0 (ref) more simply by writing mu = pd.el_refs[Element("Li")].

I think you are having reproducibility issues because you are comparing your results to a paper from 2015. The MP database has changed significantly from then, and glancing at the paper, it looks like they supplemented the data with their own DFT calculations. So the chances that you will be able to get the exact same results 8 years later is pretty slim; it’s impressive that even the first value of the table matches!

[buithitham]

Dear Sir,
I used Grand potential phasediagram code to calculation fomation energy of ZnCr2O4 as bellow:

from mp_api.client import MPRester
with MPRester('API key') as mpr:
    entries =mpr.get_entries_in_chemsys("Zn-Cr-O")
    
  

ZnCr2O4_entries = [e for e in entries if e.composition.reduced_formula =="ZnCr2O4"]

ZnCr2O4_groundstate = sorted(ZnCr2O4_entries, key = lambda e: e.energy_per_atom)[0]

from pymatgen.analysis.phase_diagram import PhaseDiagram

pd= PhaseDiagram(entries)

from pymatgen.analysis.phase_diagram import PDPlotter

plotter= PDPlotter(pd)

plotter.get_plot(pd)

pd.get_form_energy_per_atom(ZnCr2O4_groundstate)

from pymatgen.entries.computed_entries import GibbsComputedStructureEntry


entries_300 = GibbsComputedStructureEntry.from_entries(entries, 300)

pd_300 = PhaseDiagram(entries_300)
PDPlotter(pd_300).get_plot()

pd_300.get_form_energy_per_atom(ZnCr2O4_groundstate)

from pymatgen.analysis.phase_diagram import GrandPotentialPhaseDiagram

from pymatgen.core.composition import Element

oxygen = Element("Zn")

mu= 0

gpd = GrandPotentialPhaseDiagram(entries_300, {oxygen : mu})

PDPlotter(gpd).get_plot()

After that, I give results of formation energy as a bellow image:

When I read this paper (DOI : 10.1088/2515-7655/acdd9c), they show the figure 7a is Formation energy, and figure 7b is voltage profile constructed from the convex hull

So How can I make the plot of voltage like the figure 7b using pymatgen code?
Thank you so much for helping !

[Rahul1]

may you please help me to understand following :
as you have mentioned μLi=μ0Li(ref)+ΔμLi ; my query is about ΔμLi. as here it’s for interface so if i am saying lets say anode and elctrolyte ,

in that case ΔμLi = μLi(anode) - μLi(electrolyte) = applied potential (ф)
is it true ?