Colloidal interactions between clay nanoparticles have been studied extensively because of their strong influence on the hydrology and mechanics of many soils and sedimentary media. The predominant theory used to describe these interactions is the Derjaguin-Landau-Verwey-Overbeek (DLVO) model, a framework widely applied in colloidal and interfacial science that accurately predicts the interactions between charged surfaces across water films at distances greater than similar to 3 nm (i.e., ten water monolayers). Unfortunately, the DLVO model is inaccurate at the shorter interparticle distances that predominate in most subsurface environments. For example, it inherently cannot predict the existence of equilibrium states wherein clay particles adopt interparticle distances equal to the thickness of one, two, or three water monolayers. Molecular dynamics (MD) simulations have the potential to provide detailed information on the free energy of interaction between clay nanoparticles; however, they have only been used to examine clay swelling and aggregation at interparticle distances below 1 nm. We present the first MD simulation predictions of the free energy of interaction of smectite clay nanoparticles in the entire range of interparticle distances from the large interparticle distances where the DLVO model is accurate (>3 nm) to the short-range swelling states where non-DLVO interactions predominate (<1 nm). Our simulations examine a range of salinities (0.0 to 1.0 M NaCl) and counterion types (Na, K, Ca) and establish a detailed picture of the breakdown of the DLVO model. In particular, they confirm previous theoretical suggestions of the existence of a strong non-DLVO attraction with a range of similar to 3 nm arising from specific ion-clay Coulomb interactions in the electrical double layer. (C) 2020 The Authors. Published by Elsevier Inc.