By using atomistic simulation and configurational statistics techniques, the thermodynamics of mixing for calcite-smithsonite solid solutions have been investigated. By employing a 2 × 2 × 1 supercell, the configuration with the lowest energy for the solid solution with a certain composition was determined. The incorporated Zn2+ tends to occur at the sites neighboring to another substituted Zn2+ within the (0001) layer, but the substituted layers are preferentially segregated by calcite layers, and vice versa. The supercells with compositions around the two end-members stand positive enthalpies at any temperatures, whereas those supercells with composition of about Ca0.5Zn0.5CO3 prominently exhibit negative values in various temperatures of reality (e.g., <1000 K). The free energies are prominently negative at high temperatures (>1500 K) for the whole range of compositions, only those around both end-members have positive values at some low temperatures (<1200 K). In the derived phase relations of this solid solution system, the potential incorporation content of ZnCO3 into calcite is only 0–2.5% mole fraction (i.e., Zn content of 0–1.6 wt%) in most geochemistry equilibrium processes, and vice versa.