How to calculate the defect formation energy by pymatgen

Dear all,

I am trying to calculate the defect formation energy of Al doped LLZO (Li24V16O40 + Al3+(Li)-3Li → Li21AlV16O40). The defect formation energy is defined below in Chem. Mater. 2015, 27, 4040−4047. Such a definition seems to inverse a normal one: E_{defect}(normal) = -E_{defect}.
1701758877714

Besides, I am not so sure if my method to obtain chemical potential and calculate the defect energy is correct or not:

from pymatgen.core.composition import Composition
from pymatgen.core.periodic_table import Element
from pymatgen.ext.matproj import MPRester
from pymatgen.analysis.phase_diagram import PDPlotter, PhaseDiagram, PDEntry

with MPRester("API") as m:
	entries = m.get_entries_in_chemsys(['Al', 'V', 'Li', 'O'])

pd = PhaseDiagram(entries)
plotter = PDPlotter(pd, show_unstable=0.2, backend="matplotlib")
plotter.show()

pd.get_all_chempots(comp=Composition('Li21AlV16O40'))
# Obtain the one combination of chemical potential, 
pd.get_all_chempots(comp=Composition('Li24V16O40'))
# Obtain the three different combinations of chemical potentials (from different facets), and choose the chemical potentials that are the same as the values for Li21AlV16O40 (same facets)

E_pure = -528.82263 # pure structure (Li24V16O40) energy eV
E_doped = -524.61966 # doped structure (Li21AlV16O40) energy eV 

a = E_pure - E_doped + (-3)*-4.163013836249989 + -8.284360923749983 # eV/atom
print(f"The defect energy is {a} eV/atom") 

The result shows that “The defect energy is 0.0017105849999339284 eV/atom”.
Besides, I also tried to calculate the Energy above the hull for Li24V16O40 and Li21AlV16O40. Please see the below code:

pure = PDEntry(composition=Composition('Li24V16O40'), energy=-528.82263)
defect = PDEntry(composition=Composition('Li21AlV16O40'), energy=-524.61966)

# note that e_above_hull: eV/atom
decomp_pure, e_above_hull_pure = pd.get_decomp_and_e_above_hull(pure)
decomp_defect, e_above_hull_defect = pd.get_decomp_and_e_above_hull(defect)

# The energy difference obtained by e_above_hull
e_above_hull_diff = e_above_hull_pure*80 - e_above_hull_defect*78
print(f"The defect energy from e_above hull is {e_above_hull_diff} eV/atom")

The results show that “The defect energy from e_above hull is 0.0017105849999339284 eV/atom”, which is the same as the defect energy calculated by chemical potential. So I was wondering if these two methods to obtain defect energy are both right.

Best regards,
Germanicus

The obtained chemical potentials are shown below:

'LiVO2-Li3VO4-V2O3-LiAlO2': {Element Li: -4.163013836249989,
  Element Al: -8.284360923749983,
  Element V: -11.26628340624998,
  Element O: -7.717288936250014}

And the number of atoms in Li24V16O40 and Li21AlV16O40 is 80 and 78, respectively.