The free energies for the allotropes were also evaluated using quasiharmonic approximation (QHA) method implemented in Phonopy program[1]. They were calculated from 0 to 1000 K at ambient pressure as follows,
F=U(V_0)+\frac{1}{2}\sum\limits_{\boldsymbol{q}\nu}\hbar\omega_{\boldsymbol{q}\nu}+k_BT\sum\limits_{\boldsymbol{q}\nu}\ln[1-\exp(-\hbar\omega_{\boldsymbol{q}\nu}/k_BT)] where U(V_0) is the ground state energy at the relaxed volume V_0, T is temperature, and \omega_{\boldsymbol{q}\nu} is the phonon frequency for the band \nu at wave vector \boldsymbol{q}.
F_e(V,T)=E_e(V,T)-TS_e(V,T)
S_e(V,T)=-k_B\int n(\varepsilon,V)[fln f+(1-f)ln(1-f)]d\varepsilon
E_e(V,T)=\int n(\varepsilon,V)f\varepsilon d\varepsilon-\int^{\varepsilon_F}n(\varepsilon,V)\varepsilon d\varepsilon
, where n({\varepsilon,V}) is DOS, \varepsilon_F is Fermi level, and k_B is boltzmann constant (8.6173324E-5 eV/K), f=1/(e^{(\varepsilon-\varepsilon _F)/k_BT}+1) is the Fermi–Dirac distribution function.
ISMEAR = -1
SIGMA = 0.025852 # 300K: 8.6173324E-5*300=0.025852, to use 10K as the referenceplot a curve between free energy change relative to the reference (TOTEN in OUTCAR) vs. temperature
但是我发现使用上面公式的计算结果和VASP对不上，能对的上的公式应该是：
F_e(V,T)=E_e(V,T)-TS_e(V,T)
S_e(V,T)=-k_B\int n(\varepsilon,V)[fln f+(1-f)ln(1-f)]d\varepsilon
E_e(V,T)=\int n(\varepsilon,V)f\varepsilon d\varepsilon-\int^{0}n(\varepsilon,V)\varepsilon d\varepsilon
, where \varepsilon is the energy relative to Fermi Level E_F, k_B is boltzmann constant (8.6173324E-5 eV/K), n({\varepsilon,V}) is DOS, f=1/(e^{\varepsilon/k_BT}+1) is the Fermi–Dirac distribution function.
使用这个公式的结果与VASP ISMEAR=-1的E_{el}, S_{el} 和F_{el} 严格相等。