A grand canonical ensemble molecular dynamics (GMD) simulation method has been adapted to examine the thermodynamics of clay-mineral hydration. In the GMD method, the number of water molecules in the system is treated as a continuous variable for which an equation of motion is established. Fluctuations in the water content at constant chemical potential are investigated using trajectories of this particle number variable. A bias potential may be used to modify the free energy contour along the particle number coordinate. This catalyzes particle fluctuations and greatly improves simulation convergence. Adaptation of the GMD method to treat hydrated clay minerals included the introduction of a local-control technique that fixes the water chemical potential in the clay interlayer region. In addition, a bias-potential feedback algorithm was implemented to improve particle fluctuation efficiency. Information pertaining to the free energy contour, generated during the course of the simulation, was used periodically to enhance the bias potential. This allowed for the utilization of a single input bias potential under a broad range of simulation conditions. The method was used to investigate swelling of a cesium-montmorillonite clay. Measured disjoining pressures showed oscillations that are indicative of crystalline-swelling phase transitions. Integration of the disjoining pressures yielded a swelling free energy profile with distinct free-energy minima for the one- and two-layer hydrates. The results may be compared qualitatively with both clay swelling and surface force apparatus experiments, and with previous simulation studies of simple fluids in slit pores.