Each formula is normalized to contain exactly one atom, after which the reaction is balanced to get coefficients. In this way the coefficients are for atomic fractions, and not formula units, so that you can use energies per atom for each compound to get e_above_hull. Let's see it on LLZO example:

Here is decomposition path from materials project: https://next-gen.materialsproject.org/materials/mp-942733?formula=Li7La3Zr2O12

(Li7La3Zr2O12)/24 -> 5⁄8 (Li₆Zr₂O₇)/15 + 1⁄16 (Li₂O)/3 + 5⁄16 (La₂O₃)/5
1/24 Li7La3Zr2O12  -> 1⁄24 Li₆Zr₂O₇ + 0.5⁄24 Li₂O + 1.5⁄24 La₂O₃
Li7La3Zr2O12  -> 1 Li₆Zr₂O₇ + 0.5 Li₂O + 1.5 La₂O₃

import numpy as np
 
Li7La3Zr2O12 = np.array([7,3,2,12])
Li6Zr2O7 = np.array([6,2,7])
Li2O = np.array([2,1])
La2O3 = np.array([2,3])
 
LLZOn = Li7La3Zr2O12/sum(Li7La3Zr2O12)
LZOn  = Li6Zr2O7/sum(Li6Zr2O7)
LiOn  = Li2O/sum(Li2O)
LaOn  = La2O3/sum(La2O3)
 
print('Li balance:', LLZOn[0] - (5/8*LZOn[0] + 1/16 * LiOn[0])  )
print('La balance:', LLZOn[1] - (5/16*LaOn[0])  )
print('Zr balance:', LLZOn[2] - (5/8*LZOn[1])   )
print('O  balance:', LLZOn[3] - (5/8*LZOn[2] + 1/16 * LiOn[1] + 5/16*LaOn[1])  )
#if you run this code, you will get all zeroes

However, if you want to calculate energy above hull it is better to use normal formulas, find coefficients and then normalize by the number of atoms in reactant.

Let's take as example mp-1020016 the reaction according to MP is

0.2381 Li7La3Zr2O12 + 0.7619 Li4Ca2Mg(SiN3)2 -> 0.1429 Li5La5(SiN3)4 + 0.4762 Li2ZrN2 + 0.9524 Li2SiN2 + 0.5714 Li2O + 0.7619 MgO + 1.524 CaO

In this case the coefficients are stoichiometric and balance the reaction. The reaction energy is calculated E(prodcuts) - E(reactants) and normalized to total number of atoms. its energy is -0.057 according to MP API

https://github.com/materialsproject/emmet/blob/main/emmet-core/emmet/core/thermo.py https://pymatgen.org/pymatgen.analysis.html#pymatgen.analysis.phase_diagram.PhaseDiagram.get_decomp_and_phase_separation_energy

get_decomp_and_phase_separation_energy(entry: PDEntry, space_limit: int = 200, stable_only: bool = False, tols: Sequence[float] = (1e-08,), maxiter: int = 1000, **kwargs: Any)→ tuple[dict[PDEntry, float], float] | tuple[None, None]

For unstable entries that are not polymorphs of stable entries (or completely novel entries) this is simply the energy above (or below) the convex hull.

For entries with the same composition as one of the stable entries in the phase diagram setting stable_only to False (Default) allows for entries not previously on the convex hull to be considered in the combination. In this case the energy returned is what is referred to as the decomposition enthalpy in:

Bartel, C., Trewartha, A., Wang, Q., Dunn, A., Jain, A., Ceder, G., A critical examination of compound stability predictions from machine-learned formation energies, npj Computational Materials 6, 97 (2020)

For stable entries setting stable_only to True returns the same energy as get_equilibrium_reaction_energy. This function is based on a constrained optimization rather than recalculation of the convex hull making it algorithmically cheaper. However, if tol is too loose there is potential for this algorithm to converge to a different solution.

Ta-doped LiGdO2 (mp-754204)

Number of atoms: 63 16 11 40 1

mp-21332 mp-755013 mp-504886 mp-50

Energies per atom: Li15Gd15Ta1O32 = -9.35497095047619 Li4Gd4O8 = -8.90854510625 #-9.34 Li5Ta1O5 = -6.553853966363636 # -6.536 Gd16O24 = -11.42174259375 # -12.0819 Ta = -11.83995368 # -11.85

Li15Gd15Ta1O32 → 0.6984 * Li4Gd4O8 + 0.13968 * Li5Ta1O5 + 0.15873 * Gd16O24 + 0.00317 * Ta