Predicting phase relations and reactions involving phyllosilicates is of critical importance when modeling geological reservoirs and making thermobarometric estimates at temperatures below 350 °C. However, it is difficult to perform experimental studies of phase relations of phyllosilicates because of very slow reaction rates. Simple methods of estimations of thermodynamic properties of minerals such as oxide summation techniques are insufficient to constrain activity models for minerals. In this study, we use a recent atomistic technique that allows constraining activity models in minerals independently from kinetics. Among the various solid solutions occurring in phyllosilicates, the magnitude and the thermodynamic significance of the pyrophyllitic substitution has strong implications on the stability of clay minerals at surface conditions. We have applied a lattice energy estimation method combined with Monte Carlo simulations to estimate the energy of mixing along the muscovite-pyrophyllite solid solution at 25 °C and 1 bar. The results suggest a strongly positive and asymmetric Gibbs free energy of mixing, concordant at first order with previous thermodynamic models issued from phase relations and with field observations. The calculated solvus implies that pyrophyllite-rich phyllosilicates are unstable, unless another phenomenon such as hydration stabilizes them. Calculated structures at low muscovite contents present large variations of interlayer occupancies due to short- and long-range ordering. The observed ordering suggests that sub-layers in illite/smectite mixed layers minerals are not independent as the alternation of K-rich and K-depleted sublayers minimizes the Gibbs free energy of the mineral.