We have employed atomistic simulation techniques to investigate the thermodynamics of mixing in the solid solutions of calcite (CaCO3) and rhodochrosite (MnCO3). Our calculations show that the fully disordered solid solution has positive enthalpies of mixing for the entire range of compositions, which confirm recent experiments. The consideration of a small degree of ordering in the simulations leads to mixing enthalpies in quantitative agreement with experimental measurements. We argue that earlier measurements of negative mixing enthalpies for the Mn-rich solid solution were probably due to relatively high degrees of ordering in the samples. Our calculations show that the lowest energy configuration for each composition is always the one that maximizes the homogeneity of the cations within (0001) layers but maximizes the heterogeneity across layers. In particular, for Mn/Ca = 1, the most stable configuration corresponds to the ordered structure of kutnahorite, where layers of Ca and Mn ions alternate along the c axis, similar to the Ca/Mg ordering in dolomite. Our simulations predict that kutnahorite becomes less stable than the fully disordered 50:50 solid solution at similar to 850 K, and this disordering temperature decreases to a value in better agreement with experiment (695 K), if a transition to a partially ordered structure is considered. Our results thus suggest that the "disordered" (Mn,Ca)CO3 solid solutions, which are known to be favored by kinetic factors, are actually not fully disordered, but contain a higher abundance of lower-energy cation arrangements than that expected from a completely random distribution.