Vacancy formation energy and chemical potential TLDR: What exactly is the difference between vacancy formation energy and chemical potential? Does the vacancy formation energy include, e.g., the energetic cost of bond-breaking, or is this contribution substracted from it via the chemical potential?
I have a conceptual problem on the notion of vacancy formation energy and/or chemical potential. I will try to expose my (mis)understanding of the problem, feel free to point out my mistakes.
My question will be focused on a defect with neutral charge state in a solid crystal for simplification. 
Chemical potential:
the chemical potential of a species A can be defined by the differential relation:
$$dU = TdS - PdV + \sum_A{dN_A \mu_A}$$
with $N_A$ the number of A particles. Let's simplify in a particular context by assuming that:


*

*We are at $0K$ (ground state calculation)

*The enthalpic term PdV is negligible for a solid OR we are working at $0$ pressure

*We are looking at the formation of a single vacancy (only one species A in the sum, and $\Delta N_A = -1$ between initial and final state)

*constant chemical potential between initial and final state


We then simply have (after integration between perfect crystal and 1-vacancy crystal):
$$-\mu_A = \Delta U $$ 
Vacancy formation energy from DFT calculations
In many articles (for example: https://aip.scitation.org/doi/abs/10.1063/1.1682673), vacancy formation energy is calculated from DFT results in the following way:
$$E_f  = E^{tot}_{vacant} - E^{tot}_{bulk} - \sum_A n_A \mu_A + qE_{Fermi}$$
where $E^{tot}_{vacant}$ is the total energy of a crystal supercell with 1 vacancy and $E^{tot}_{bulk}$ for a  perfect crystal.
The last term is for charged defects, and can be seen as yet another chemical potential term, but for additional/missing electrons (the fermi energy is the chemical potential of an electron at $0K$). 
Again, by assuming neutral defect ($q=0$) and single vacancy ($n_A=-1$), we can rewrite this as:
$$E_f  = \Delta E_{DFT} + \mu_A $$
Now, My problem should appear pretty obvious. The formation energy calculated in this way should always be $0$ according to the thermodynamics relation written above... (provided $\Delta E_{DFT}$ has been calculated precisely enough to coincide with the actual $\Delta U$)
In short: I thought the formation energy of a vacancy was telling you how much energy is needed (in terms of bond-breaking, etc) for a species to escape the crystal potential and create a vacancy. However, this is also my understanding of the chemical potential ! The variation of energy according to particle number...
Specifically, the fact that the chemical potential depends on the compound that a species is in, reflects that it encompasses the contribution of interaction (e.g. bonding) with the neighbouring particles.
Therefore, I am well confused as to why the chemical potential is being substracted, when it seems that it encompasses all the information we are looking for. 
Thanks for reading!
 A: The chemical potentials used in the equation for the defect formation energy are the energies at which atoms are taken from / given to reservoirs to form the defect in the solid. The electronic chemical potential similarly is the energy at which electrons are exchanged to form charged defects (i.e. the Fermi level, which is mid-band gap for insulators and depends on doping for semi-conductors).
Explaining the equation for the formation energy $E_{\rm f}$ in words: it is the total energy of the solid with defect (also accounting for the energy associated with distortion of the host lattice) minus the total energy of the solid without defects minus the total energy of the defect back in its reservoir (including the equivalent term for electrons if the defect is charged). One can also consider it as the calculation of a reaction energy: product is the solid with defect, reactants are the defect free solid and the defect still in its reservoir. If dealing e.g. with a neutral phosphorus defect, the chemical potential for it used in the formation energy is the total energy of bulk phosphorus per phosphorus atom; for a neutral nitrogen defect it would be half the total energy of N$_2$ in gas phase (bulk phosphorus and N$_2$ gas are the thermodynamic particle reservoirs for the corresponding defects in the host solid).
The chemical potential used here is the energy of the defect in the reservoir not in the host solid. If the chemical potential referred to the energy in the host solid, indeed the definition of the defect formation energy would be nonsensical and always be zero.
To relate this to the specific case mentioned in the question: creation of a vacancy in a crystal. The chemical potential being subtracted from the vacancy formation energy is accounting for the fact that once you have removed an atom from the crystal (leaving behind the vacancy), this atom has to go somewhere at a certain energy. Depending on what species (and conditions such as temperature and pressure or applied potentials in case of charged defects one is dealing with), this atom will typically bind with other (removed) atoms outside the crystal it was taken from. Metal atoms for example will form a bulk metal phase. The energy the removed atom has in this other phase, is what the subtracted chemical potentials in the vacancy formation energy formula account for.
A reference (picked out of many possible ones since it is not behind a pay wall): Jund et al. arXiv:1309.7246 "Lattice stability and formation energies of intrinsic defects in Mg2Si and Mg2Ge via first principles simulations", where, e.g., as chemical potential reference for intrinsic Si defects in Mg2Si $\mu({\rm Si}) = E({\rm Si})$ is used ($E({\rm Si})$ is the energy of bulk Si per Si atom; above linked PDF, bottom of page 4).
