We have used a combination of classical and quantum-mechanical atomistic calculations, together with Monte Carlo simulations, to study order-disorder phenomena in the spinel mineral qandilite, Mg2TiO4. Using an interatomic potential model akin to those previously used for 2–3 spinels yielded a general increase in energy E as a function of inversion parameter x, and thus incorrectly predicted a normal-spinel ground state, whereas the E(x) behavior as modeled by density-functional theory exhibited a maximum at an intermediate degree of inversion and correctly predicted an inverse-spinel ground state. We therefore used the quantum-mechanical simulations to derive pair interaction parameters (for nearest-neighbor tetrahedral-tetrahedral, octahedral-octahedral, and tetrahedral-octahedral interactions) and chemical potential to use in Monte Carlo simulations of order-disorder in qandilite. The simulated cation distributions compared favorably with those obtained experimentally, although the long-range ordering transition to the tetragonal P4122 phase was not observed when using only nearest-neighbor interactions. However, this transition could be observed following the addition of two extra parameters to the model. The simulations were used to calculate the effect of short- and long-range cation order on the configurational entropy of qandilite as a function of temperature. The calculated entropy of the high-temperature cubic phase was in very good agreement with the experimental value recently determined, supporting the suggestion that there is considerable short-range order in qandilite.