Thank you for answer for my questions! That is really useful. I find a problem in phb result. When we consider the chemical potential of ground states. We simply just need to calculate (E(gs=0)-E(gs=1))/(x(gs=0)-x(gs=1)) to obtain chemical potential that stablize the two phases, right? I compare my result to phb result. However there is small differences. Is this a bug? Or I made mistake in obtaining the chemical potential? Thanks for your help!
Yes, at 0K this seems correct. Note that phb (and all other atat codes) use the energies after subtracting the reference energies of the pure end members. This would be done automatrically if you use the energies in fit.out (from maps), but not if you take the VASP energies.
Does the reference energies matter in phase boundary? I try to set ref-energy to be 0 and do phase boundary to compare with those already been subtracted the reference energies, and I find those two groups of phase boundary are quiet similar. This may indicates ref-energy does not inflence much. However I am thinking about one problem. The nature of phase boundary is that, for phase 0 and 1, thermaldynamic potential phi_0 = phi_1. phi_0 = Helmholz free energy - mu*x. At 0K, Helmholz free energy = Internal energy from VASP. So in principle, mu = (E1 - E0)/(x1 - x0) which should be difference of internal energy divide by difference of composition x. There is the problem, if we consider ref_energy, the energy input of phb code is actually enthalpy of formation H. Then the mu is calculated by difference of enthalpy of formation divided by difference of composition. In this case, reference energy should matter a lot in determing the chemcial potential and because difference of energy is not equal to difference of enthalpy of formation. I wonder whether this will effect the phase boundary. What is your opinion about this?