Bottom-up modeling of clay behavior from the molecular scale requires a detailed understanding of the free energy between pairs of clay platelets. We investigate the potential of mean force (PMF) for hydrated clays in face-to-face interactions with free energy perturbation (FEP) methods through molecular dynamics simulations using simple overlap sampling (SOS). We show that PMF results for open systems with one finite in-plane dimension are affected by migration of counterions from within the interlayer space compared with fully confined closed system conditions. We compare PMFs for two common 2:1 clay sheet minerals Illite (IMt-1) and Na-smectite. The PMFs for the open illite systems exhibit a strong attractive energy well at a basal layer separation, d = 11 Å and interlayer water content, wIL = ∼0.4% while the attractive minimum for the closed system occurs at d = 12 Å, wIL = 3.5%. In contrast, net repulsion occurs between pairs of Na-smectite platelets for both open and closed systems (for d < 15-16 Å). The free energy is closely related to the distribution of counterions; while K+ ions are bound closely to the surfaces of the illite platelets, Na+ ions are more spatially disperse. This PMF results contradict prior findings for Na-smectite and prompted further comparisons with other published results. We find that most of the published results do not represent accurately the free energy for face-face interactions between pairs of clay platelets that are effectively infinite (with width/thickness O[104]). The PMF results presented in this paper form a reliable basis for mesoscale, coarse-grained modeling of illite and smectite particle assemblies. We show that the Gay-Berne potential provides a reasonable first-order model for upscaling, while the solvation potential proposed by Masoumi enables a more accurate representation of the computed PMFs.