Although multiple types of adsorption sites have long been observed in montmorillonite, a consistent explanation about the chemical structure of these adsorption sites has not yet been established. Identifying the cation interlayer adsorption sites based on the octahedral cation distribution on montmorillonite was investigated in this study by using a Density Functional Theory (DFT) simulation. A clay structural model (H[Al6MgFe]Si16O40(OH)8) with a similar composition to Wyoming SWy-1 montmorillonite was built, where two octahedral Al were respectively substituted by Fe and Mg, and H+ was the charge compensating cation. This model had twenty-one different possible configurations as a function of the distribution of octahedral Al, Fe, and Mg cations. The DFT simulations of 15 of these different configurations showed no preference for the formation of any configuration with a specific octahedral Fe-Mg distance. However, the H+ adsorption energy was separated into three distinct groups based on the number of octahedral jumps from Fe to Mg atoms. The H+ adsorption energy significantly decreased with increasing number of octahedral jumps from Fe to Mg. Assuming an even probability of occurrence of 21 octahedral structures in montmorillonite, the percentages of these three groups are 43, 43, and 14%, respectively, which are very close to the three major sites on montmorillonite from published cation adsorption data. These DFT simulations offer an entirely new explanation for the location and chemical structure of the three major adsorption sites on montmorillonite, namely, all three sites are in the interlayer, and their adsorption strengths are a function of the number of octahedral jumps from Fe to Mg atoms.