ABSTRACT The majority of close massive binary stars with initial periods of a few days experience a contact phase, in which both stars overflow their Roche lobes simultaneously. We perform the first dedicated study of the evolution of massive contact binaries and provide a comprehensive prediction of their observed properties. We compute 2790 detailed binary models for the Large and Small Magellanic Clouds each, assuming mass transfer to be conservative. The initial parameter space for both grids span total masses from 20 to 80$\, \mathrm{M}_⊙$ , orbital periods of 0.6–2 d and mass ratios of 0.6–1.0. We find that models that remain in contact over nuclear time-scales evolve towards equal masses, echoing the mass ratios of their observed counterparts. Ultimately, the fate of our nuclear-time-scale models is to merge on the main sequence. Our predicted period–mass ratio distributions of O-type contact binaries are similar for both galaxies, and we expect 10 such systems together in both Magellanic Clouds. While we can largely reproduce the observed distribution, we overestimate the population of equal-mass contact binaries. This situation is somewhat remedied if we also account for binaries that are nearly in contact. Our theoretical distributions work particularly well for contact binaries with periods <2 d and total masses $⪅ 45\, \mathrm{M}_⊙ \,$. We expect stellar winds, non-conservative mass transfer, and envelope inflation to have played a role in the formation of the more massive and longer-period contact binaries.