Soils are known to be “skin of the earth” and carry the material basis for human beings. As one of the basic disciplines of agriculture, soil science deals with the natural resources of the earth surface relating to soil formation, classification, morphology as well as associated physical, chemical, biological and fertility properties. Soils are a unique array of different substances naming clay minerals (e.g., montmorillonite and kaolinite), organic matters (e.g., humus), microorganisms (e.g., bacteria) as well as liquids (mainly moisture) and gases that are intercalated among clay mineral layers. Clay minerals are definitely the main ingredient for soils, and even within a single sample, they vary continuously in crystalline degree, chemical composition, charge distribution and other properties [1], and such complicated and variable chemical compositions endow soils with magic “physiological fertility” [2]. As indicated, soils have the ability of self-regulation, which is a natural result of the composite complexes of clay minerals interacted with soil inorganics, organics, enzymes and other constituents. It is evident that unravelling the physiological function of soil systems has to rely strongly on chemistry discipline. What frustrates us is that up to date, the direct and quantitative measurements of soil systems remains a challenge for currently advanced experimental techniques.

Computational chemistry is the application of various theories incorporated into efficient computer programs to resolve chemical problems, such as total energy, electronic charge distribution, dipoles, spectroscopic modeling as well as reaction reactivity, kinetics and mechanism. Researchers have gradually realized that chemistry is no longer a purely experimental discipline since 1998 when Kohn and Pople were awarded with the Nobel Prize of Chemistry for their outstanding contributions in the development and application of computational methods. This is not fortuitous but because of the momentous status of computational chemistry, and last year (2013) the Nobel Prize of Chemistry again patronized computational chemists. A number of computational methods have successively emerged and each has its respective merits: density functional and ab initio quantum chemical methods are precise but rather slow and therefore are restricted to relatively small systems, while molecular mechanics and molecular dynamics are prompt and effective to treat very large systems but the results thus obtained are comparatively coarse. Recently, the hybrid quantum mechanics/molecular mechanics (QM/MM) approach has been developed that combines the merits of QM (accuracy) and MM (superability in speed and model size) [3]. In QM/MM, a small portion of systems that is of central interest is treated by QM while the rest by MM, thus realizing high accuracy and efficiency at the same time. In contrast to the experimental counterpart, computational chemistry is capable of “easily” selecting one target (e.g., a specific montmorillonite surface from the various clay soils) and investigate the numerous influencing factors one by one; in addition, it can in situ “observe” the structural change and dynamics on the femtosecond time scale and at an atomic level. Accordingly, computational chemistry should have made substantial achievements in soil science, which is discouraged by the current situations, however. Even very few people are engaged in the computational studies of soils in contrast to other fields.

As commented by Prof. Borden (Editor) [4], manuscripts that focus on the methodology development are rarely Accepted for publication in Journal of the American Chemical Society (JACS), but this does not belittle its significance at all. The methodology development is a prerequisite for performing molecular dynamics, density functional or other calculation tasks. The DLVO theory, which is proposed by the Soviet scientists Boris Derjaguin and Lev Landau as well as independently by the Dutch scientists Evert Verwey and Theo Overbeek, has long been regarded as one cornerstone in colloid sciences and helps greatly to interpret the experimental results such as colloidal stability, heterocoagulation and ion adsorption. It combines the van der Walls (vdW) attraction and electrostatic (double-layer) repulsion that are respectively described by the Lifshitz force and Poisson-Boltzmann equation. In regard to the DLVO theory, the interaction potential energy (W_T) between two colloids or surfaces should be W_T = W_vdW + W_EDL, where W_vdW is not affected by presence of salts while W_EDL is sensitive to salt concentrations [5]. With increase of salt concentrations, W_EDL decreases accordingly and results in the more attractive total interaction energy (W_T).

The DLVO theory had ever achieved tremendous success but in the past two decades, it receives more criticism than praise due to the treatment of finite ion sizes and self-atmosphere effects (ion correlation) in an inconsistent manner. As matters of fact, many other forces (non- DLVO forces) that may be significant are found to be absent in the DLVO expression, such as steric repulsion, depletion, oscillatory structural undulation, hydrophobic, hydration and ionic correlation. Some researchers consider that the deviation of the DLVO theory is due to the more pronounced collective Coulomb interactions among charged species for concentrated systems, while others insist that the DLVO theory governed by the ordinary Poisson-Boltzmann equation is not accurate even for dilute systems [6]. A number of corrections have been made aiming to improve the DLVO theory. Even with five adjustable parameters, the DLVO theory remains contradictory with the experimental results at relatively short ranges. According to the report of Boström et al. [7], the inconsistency of the DLVO theory with experiments should be caused by the improper account of specific ion effects, which are ubiquitous in colloidal, biological and soil systems. Very different results from the traditional DLVO theory will be obtained by inclusion of dispersion forces accounting for the ion specificities. The Poisson-Boltzmann equation corrected with steric, polarization and valence effects describes well the selective adsorption of different ions on the charged surfaces, showing potential for demonstrating specific ion effects in soil, colloidal and biological systems [8]. Very recently, the aggregation kinetics of colloidal particles in the various alkali ion solutions has been studied, and the ion specificities implicated therein are quantitatively characterized by use of activation energies [9], which are further ascribed to polarization effect. To our surprise, the classical theory even predicts an opposite trend for the experimentally observed ion specificities and hence, the theoretical basis for specific ion effects requires urgent corrections.

Clay minerals such as kaolinite and montmorillonite are now within the power of periodic ab initio calculations [10]. In our group, QM/MM methods are being used to model the adsorption behaviors of metal ions onto the surfaces of clay minerals. As is known to us, all processes relating to clay minerals should proceed in fully watersolvated environments, which are readily available from classical molecular dynamics (MD) simulations with the use of supercells. The qualities of classical MD are significantly dependent on the choice of force fields. Note that the force field parameters are derived from experiments or high-level QM calculations. Aluminum (Al) in solidstate materials is usually four-coordinated to O atoms, which can be finely modeled by a variety of force fields that are familiar to us; nonetheless, Al in clay minerals adopts an octahedral coordination. Take kaolinite Al₂Si₂O₅(OH)₄ for example. Kaolinite is a 1:1 layer clay mineral composed of a repeated layer of gibbsite Al octahedral sheet and siloxane silicon tetrahedral sheet, and these two sheets are linked by bridging-O atoms. Accordingly, the proper modeling of octahedrally coordinated Al renders to be the bottleneck for extending classical MD simulations to clay minerals. To the best of our knowledge, CLAYFF [11] is the only force field that is developed especially for clay minerals. Classical MD simulations using the CLAYFF force field have been widely used to model clay mineral systems. Igor et al. [12] studied and compared the adsorption behaviors of different ions on the basal surfaces of kaolinite. The cations and anions are respectively absorbed on the siloxane and gibbsite surfaces, and the Cs+ and Na+ rather than Cd2+ and Pb2+ ions are driven by Cl- anions forming the inner-sphere adsorption model. No inner-sphere adsorption of Cd²⁺ or Pb²⁺ ion with kaolinite has been detected, and instead relatively strong outersphere complexes are formed for these two transition metal cations. The composite materials of montmorillonite and intercalated poly (ethylene glycol) (PEG) with or without CO₂ have been simulated by Krishnan et al. [13]. Albeit fruitful results have been obtained thus far, the CLAYFF force field is far from being perfect: Only a limited number of elements have been included, and more importantly, it seems rather robust with only bonded interactions being within water molecules and hydroxyl groups, as claimed by the developers of CLAYFF force field [11]. Accordingly, the task of top priority for computational soil chemists is to develop an appropriate force field capable of correctly modeling clay mineral associated systems, and only in this way, can MD simulations be used to interpret experimental results and make new predictions that deserve our trust.

The authors thank the National Natural Science Foundation (D010503), and Fundamental Research Funds for the Central Colleges (SWU113049 and 2362014xk01) for support of this research.