The order-disorder behavior of the isomorphous cation substitution of the octahedral sheet of phyllosilicates was investigated by Monte Carlo simulations based only on atomistic models in some three-species systems Al/Fe/Mg including a wide range of different octahedral compositions that can be relevant to clay compositions found in nature, especially for smectites and illites. In many cases, phase transitions do not occur, in that long-range order is not attained, but most systems exhibit short-range order at low temperature. The ordering of the octahedral cations is highly dependent on the cation composition. Variations in the tetrahedral charge (smectite vs. illite) produce slight differences in the cation distribution and the short-range and long-range order of octahedral cations do not change drastically. The average size of Fe clusters and the long-range order of Fe are not larger in illites than in smectites as previous reports concluded, but the proportion of Fe3+ cations non-clustered is higher in smectites than in illites. This behavior supports the experimental behavior of the Fe effect on the Al-NMR signal, which is lower in illites than in smectites.