Reliable prediction of the behaviour of radionuclides and their transport and retention in clayey formations at nuclear waste repositories requires detailed molecular scale understanding of these complex multicomponent systems. As the first step in our study of the effects of organic molecules on the adsorption and transport of radionuclides in hydrated clay systems we have investigated the effects of the ordering in charge distributions on the swelling behavior of montmorillonite (a smectite clay). Montmorillonite layered structure consists of aluminum-oxygen octahedral sheet sandwiched between two opposing silicon-oxygen tetrahedral sheets giving rise to a 2:1 clay mineral. Isomorphic substitutions in the tetrahedral and octahedral sheets are responsible of the negative layer charge of montmorillonite clay minerals having the chemical composition (Si8-xXx)(Al4-yYy)O20(OH)4 where X = Al3+, Y = Mg2+, Fe2+...[1]. The montmorillonite models for our study are based on a pyrophillite unit cell structure (5.16Å×8.966Å×9.347Å) obtained from the crystallographic data of Lee et al. [2]. The 4×4×2 simulation supercells were built and substitutions were made in the pyrophillite structure in order to approximate as close as possible the chemical composition of Wyoming montmorillonite [1] M24(Si248Al8)(Al112Mg16)O640(OH)128, where M is either Cs+, or K+. We have explored three different models of the substitution distributions. In the first model, the substitutions were uniformly and orderly distributed within the tetrahedral and octahedral sheets. In the second model, the substituted positions were kept ordered in the octahedral sheets but made disordered in the tetrahedral one. In the third model, the substituted positions of the octahedral sites were additionally made disordered. In order to study the swelling behavior of these montmorillonites, NPT-ensemble molecular dynamics (MD) simulations were run at T = 298 K and P = 1 bar for each of the three different substitution models and with 23 different hydration states ranging from 0 to 700 mgwater/gclay (from 0 to 42 H2O molecules per one monovalent cation). All MD runs were performed for a total of 2 ns using the CLAYFF force field [3]. After the system reached equilibrium, the last 1ns of each MD trajectory was used to compute the clay basal spacing and the swelling thermodynamic properties: hydration energy, immersion energy, isosteric heat of adsorption. These calculations indicate that in addition to the commonly observed 1-layer and 2-layer hydrates, stable hydration states corresponding to 3-layer and 4-layer hydrates can also be distinguished. These stable states (minima of hydration and immersion energies) were then selected to run further 500 ps NVT-ensemble MD simulations at the same temperature. The equilibrium parts of these NVT-simulated trajectories were then used to calculate radial distribution functions and atomic density profiles of the interlayer species in hydrated montmorillonites. References [1] Tsipursky, S.I., Drits, V.A. Clay Minerals, 19, 177-193 (1984). [2] Lee, J.H. and Guggenheim, S. American Mineralogist, 66, 350-357 (1981). [3] Cygan, R.T., Liang, J.J., Kalinichev, A.G. Journal of Physical Chemistry B, 108, 1255-1266 (2004).