Thermo-chemical properties and T–X phase relations diagram of the (Mg,Fe)O solid solution are modelled using mixing Helmholtz energy, ΔF(T,x)mixing, calculated by quantum mechanical and semi-empirical techniques. The sub-solidus MgO–FeO binary has been explored as a function of composition, with iron either in high-spin (HS) or low-spin (LS) configuration. Only the HS model provides physically sound results at room pressure, yielding a correct trend of cell edge versus composition, whereas LS’s issues are at variance with observations. Mixing Helmholtz energy has been parametrized by the following relationship: ΔF(T,x)mixing = x × y × [U0(T) + U1(T) × (x – y) + U2(T) × (x − y)2]−T × S(x,y)config, where y = 1−x and Uj(T) are polynomials in T of the second order. ΔF(T,x)mixing exhibits a quasi-symmetric behaviour and allows one to build the T–X phase relations diagram over the MgO–FeO join. The HS model including vibrational contribution to the Helmholtz energy predicts a solid solution’s critical temperature of some 950 K, remarkably larger than olivine’s and Mg–Fe garnet’s. All this points to a more difficult Mg–Fe mixing in periclase-like structure than olivine and garnet, which, in turn, provide more structure degrees of freedom for atomic relaxation. From ΔF(T,x)mixing, we have then derived ΔH(T,x)excess and ΔS(T,x)excess. The former, characterized by a quasi-regular behaviour, has been parametrized through W × x × (1−x), obtaining WH,Mg–Fe of 17.7(5) kJ/mol. ΔS(T,x)excess, in turn, increases as a function of temperature, showing absolute figures confined within 0.1 J/mol/K. Mixing Gibbs energy, calculated combining the present issues with earlier theoretical determinations of the magnesio-wüstite’s elastic properties, has shown that the HS configuration is stable and promote Mg–Fe solid solution up to ≈15 GPa