Application of Density Functional Theory in Soil Science

Soil is the basis for life and soil science is regarded as the final frontier; however, as compared to chemistry, physics, biology, and other disciplines, soil science undergoes an obviously slower development and remains almost stagnant in the past few decades, mainly due to two reasons: (1) wrong and outdated perceptions for a large portion of soil researchers; (2) complexity of soil systems that are difficult to characterize by current experimental techniques. Computer simulations have unique advantages to handle complex systems while currently, its role during soil researches is far from being recognized. In this chapter, several examples are given with respect to application of density functional theory (DFT) calculations to soil science, focusing on the adsorption of uranyl ion and SO2 onto mineral surfaces and reaction mechanisms to form acid rain. In this way, insightful clues at the atomic level are provided for the adsorption, interaction, and reactions regarding soil systems. We believe that computer simulations including DFT are the right key to unravel the complicated processes occurring in soils. More efforts of computer simulations are anticipated for soil science with aim to decipher the experimental results and probe the uncharted principles that may result in a revolutionary in the near future.


Introduction
According to Natural Resources Conservation Service (NRCS), soil is defined as a natural body comprised of solids (mainly minerals and organic matters), liquids, and gases that occurs at the intermediate surface of the Earth, occupies space and is characterized by one or both of the following properties: horizons and layers, which are distinguishable from the initial materials as a result of addition, loss, transfer, and transformation of energy and matter or the ability to support rooted plants in natural circumstances. Soil constitutes the basis for life and bridges the biosphere, atmosphere, hydrosphere, and geosphere. Despite these facts, apparently less attention has been given to soil science than to other disciplines such as physics, chemistry, and biology. As said by Gardner (the past president of Soil Science Society of America) [1], "not a few people mistakenly perceive that everything worth knowing about soils has already been understood, and all we need to do is merely to apply that knowledge properly. Even knowledgeable scientists assume that principles and theories developed from other disciplines can be applied to the researches of soil science in a straightforward way, requiring little imagination or creativity. In their opinions, soil science is just one of expressions for the applied physics, chemistry, or biology." The situation of soil science research is alarming. The core concepts of current soil science textbooks remain almost unchanged as compared to those of half a century ago. Obsolete or even incorrect standpoints are a commonplace [2]. Fortunately, a few researchers have recognized such a crisis. On the other hand, because of the complexity of systems and co-function of multiple factors, it seems challenging for experimental techniques to in-situ characterize "real" soil properties and processes; in addition, the experimental results from one lab may not be reproducible by others, since soil samples of different areas or even different batches may vary significantly. Computer simulations have unique advantages within this context: (1) probing the various influencing factors one by one; e.g., six factors (identity of heteroatoms, crystallographically distinct T sites, structural alterations, quantity of negative charges, distance from charge centers to metal ions and source of negative charges) were identified to affect the adsorption of metal ions at clay surfaces, and their respective contributions were estimated by density functional theory (DFT) calculations. The quantity of negative charges is the foremost factor that controls the adsorption processes, while other factors in certain circumstances can also play a critical role [3]. The adsorption strengths and numbers of all metal ions increase in a direct proportion to the intensities of electric fields [4]; (2) providing useful and detailed information at the femtosecond scale such as how ions from aqueous solutions diffuse toward to clay surfaces [5]; (3) understanding the adsorption, interaction and reaction processes at the atomic level such as how metal ions interact with surface-O atoms and respond to the increase of electric fields. Based on Hirshfeld, Mulliken, and NBO charge analyses, we found that polarization rather than electrostatic interactions are more likely to result in the pronounced cation-specific effects at clay surfaces [6]; (4) unraveling the exact reaction mechanism by comparing the structural and (especially) activation barriers of competing paths. This can be considered as an extension of (3). There are a plethora of competing reactions occurring in soils; e.g., with DFT calculations, it was clarified that Mn 4+ rather than Mn 3+ sites are more reactive for the oxidation of As 3+ and the oxidation processes are significantly blocked by As 5+ complexes [7]. A more convictive example is the mechanistic study of Brönsted acid-catalyzed conversion of biomass sugars [8]. More than 120 reaction paths were explored, and the low reactivity and selectivity of glucose conversions were clearly addressed: unlike fructose that prefers to dehydrate at the anomeric O 2 H group and initiates a sequence of facile reaction steps toward 5-hydroxymethyl-2-furfural (HMF), the less reactive sites in glucose (O 2 H and O 3 H) produces levulic acid not involving fructose and HMF intermediates, while the most reactive O 1 H site leads to humin precursors or reversion products [8]. The relationship between system size and computational accuracy for representative theoretical levels is shown in Figure 1. For huge systems (100,000 or even more atoms), classical methods (Monte Carlo and Molecular Dynamics) seems to be the good choice, although the computational accuracy is relatively low; In contrast, ab initio quantum mechanical methods are restricted to relatively small systems (up to several hundred atoms or even fewer atoms) while achieve high accuracy. The computational accuracy of quantum mechanics/molecular mechanics (QM/MM) [9] and semi-empirical methods fall in-between, and QM/MM methods have recently become increasingly popular due to the satisfying computational accuracy (active sites handled by QM methods) and the easy extension to large systems (all atoms except active sites disposed by classical methods). 2013 Noble Prize in Chemistry was awarded to Karplus, Levitt, and Warshel for the development of "multiscale methods for complex systems." Let us look back to ab initio quantum mechanical methods. Strictly speaking, a large portion of DFT methods fall outside this scope; e.g., B3LYP is probably the most popular DFT method and its exchange-correlation functional is written as [10,11] where a, b, and c are empirical parameters.
Computer simulations, as discussed above, are particularly useful to handle the complex soil systems, and two applications were illustrated in this chapter. Section 2 summarized DFT Figure 1. The theoretical landscape regarding to the balance between system size and computational accuracy. calculated results of uranyl adsorption onto mineral surfaces and Section 3 elaborately discussed DFT calculated results of SO 2 (one of the main precursors for acid rain) adsorption and formation mechanisms of acid rain. In the end, concluding remarks were provided and some suggestions for DFT applications to the future soil researches were presented as well. A more critical role of computer simulations should have played in soil science, and this chapter aims to arouse the attention of general soil researchers regarding to the applications of computer simulations. More efforts in this regard are greatly beneficial to decipher the experimental results and probe the uncharted principles that will result in a revolutionary for soil science in the near future.

Uranyl adsorption onto mineral surfaces
Release of radionuclides into the environment seriously threatens the ecosystem and human health, and adsorption of radionuclides onto mineral surfaces significantly affects their migration and transport into the environment [13]. Accordingly, knowledge about the interaction between radionuclides and minerals is essential for the long-term risk assessment of radioactive waste repositories. Uranium is usually present in the uranyl (UO 2
Perron et al. [16] examined UO 2 2+ adsorption onto rutile(110) surface, the most stable face of natural rutile reported by Jones and Hockey [26,27]. UO 2 2+ forms a bidentate inner sphere complex with three H 2 O molecules to fill its first hydration shell, see In addition to the [UO 2 (H 2 O) n ] 2+ adsorption complexes, other uranyl species containing anionic ligands such as OH À and CO 3 2À also play a critical role in the environmental circumstances. Pan and co-workers [17] investigated the effects of different ligands (H 2 O, OH À , CO 3 2À ) on the adsorption on rutile(110) surface. Note that only the most stable adsorption mode (bb) has been taken into account therein. The uranyl ion (UO 2

2+
) interacts more strongly with anionic ligands  ) that strengthen the adsorption processes.
The adsorption of uranyl species on tetrahedral SiO 4 surfaces is thermodynamically unfavorable, with the adsorption energies of inner-and outer-sphere complexes being, respectively, 239 and 206 kJ/mol [20]; in contrast, the adsorption at octahedral AlO 6 surfaces is preferred due to the presence of upright (perpendicular to the surface) and lying (parallel to the surface) OH groups. The most stable configuration corresponds to an inner-sphere monodentate complex, where UO 2 2+ bonds directly to the lying surface-OH group and has an adsorption energy of À155 kJ/mol. Martorell et al. [21] continued to study the adsorption of uranyl ion on the bare and solvated octahedral AlO 6 surfaces of kaolinite, where UO 2 2+ forms two direct bonds with the surface-O atoms of deprotonated OH groups. Two different adsorption structures are assigned according to the coordination of surface-O atoms: two surface-O atoms connected to one and two Al atoms are designated to be AlOO (short-bridge site) and AlO-AlO (long-bridge site), respectively. The formation of adsorption complex at the short-bridge site (AlOO) requires an energy of 195 kJ/mol and is obviously less than at the long-bridge site (AlO-AlO, 261 kJ/mol). Thus, the uranyl ion prefers to adsorption at the short-bridge site; furthermore, similar trends for uranyl adsorption remain when including solvent effects by adding a monolayer of water molecules.
Gibbsite is a primary mineral form of aluminum hydroxide, and adsorption of the uranyl ion onto its (001) surface [18,24] and edge surface [25] have been investigated, similar to the situation of kaolinite discussed above.

Acid rain
In 1852, Smith demonstrated the relationship between acid rain and atmospheric pollution in Manchester, and after 20 years (i.e., 1872), he coined the term "acid rain." Now acid rain has become a popular term and one of the world's biggest environmental concerns, especially in North America, Europe, and China [28]. Acid rain refers to any form of precipitations with acidic components that fall to the ground from the atmosphere, including rain, snow, flog, hail, or even dust that is acidic. It results mainly from SO 2 and NO x emissions to the atmosphere and the further transport by wind and air currents, during when SO 2 and NO x react with water, oxygen, and other substances leading to the formation of sulfuric (H 2 SO 4 ) and nitric (HNO 3 ) acids. The pH of acid rain is approximately 5.6 and since 1940s, researchers began to recognize its strong impacts on the ecosystem and human health so that soils, freshwaters, forests, and buildings will be damaged. With regard to soils, acid rain inhibits the decomposition of organic matter [29], fixation of nitrogen [30], elution of calcium (Ca 2+ ), magnesium (Mg 2+ ), potassium (K + ), and other nutrients [31]. As a result, soil fertility and microbial activity show an obvious reduction [32]. Geochemical modeling indicated that Ca 2+ leaching in marble due to acid rain neutralization approximates 0.158 mmol/L, in contrast to 10.5 mmol/L by dry deposition, and the corresponding Cu 2+ losses in bronze are ca. 0.21 and 47.3 mmol/L, respectively [33].
As aforementioned, SO 2 emissions are one of the principal causes of acid rain and also represent a primary source of atmospheric aerosols, which can lead to respiratory diseases, premature deaths, and even climate changes by affecting the properties of clouds and the balance of solar radiation. Therefore, it is of great significance to convert SO 2 to other less contaminated compounds, and a number of measures to control SO 2 emissions have been proposed. During 1982-1999, SO 2 emissions have reduced by approximately 65% in Europe and 40% in the United States, and SO 2 emissions in China decline in the late 1990s while again increase after then. DFT calculations provide useful information about the adsorption of SO 2 onto mineral surfaces as well as reaction mechanisms that seem difficult to capture by current experimental techniques [34][35][36][37][38][39][40][41][42][43][44][45], which are, however, critical to understand the formation of acid rain at the molecular level and to remediate the ecosystem. Clay minerals, such as alumina (Al 2 O 3 ), iron oxides (Fe x O y ), are good candidates for the adsorption of acid components from acid rain and then convert them into less hazardous compounds. Lo et al. [35] studied the adsorption of SO 2 on clean (100), dehydrated (110), and hydrated (110) surfaces of γ-Al 2 O 3 , finding that significant adsorption differences exist for the various surfaces and the calculated adsorption energies (À13 to -85 kcal/mol) are consistent with experimental results.
The γ-Al 2 O 3 (100) surface is composed by bridging-O and five-coordinated Al atoms, and a total of five stable configurations are produced for SO 2 adsorption (Figure 4). The feeble interaction between S and surface-O atoms results in the physisorption configuration (CM3) with the S-O distance of 2.915 Å, and the corresponding binding energy is very small (À2.0 kcal/mol). The interaction between O@SO 2 and Al atoms leads to a chemisorption state named CM4, and the O-Al bond distance and adsorption energy are 2.123 Å and À23.9 kcal/mol, respectively. The other three configurations are also ascribed to chemisorption. In CM5, one O@SO 2 atom is coordinated to two Al atoms in the vicinity of an octahedral vacancy, and CM7 and CM8 can be considered to generate from CM5 conversion and recombination. As compared to CM5 and CM7, the adsorption configuration CM8 possesses a superior symmetry, and both O@SO 2 atoms participate in the formation of direct bonds with the Al atoms. The adsorption energies of SO 2 are calculated to be À39.2, À15.1, and À45.2 kcal/mol, respectively, for CM5, CM7, and CM8. In consequence, three types of sulfite (SO 3 2À ) are produced during the adsorption of SO 2 onto γ-Al 2 O 3 (100) surface. For all adsorption configurations including three with positive adsorption energies (CM1: 1.0 kcal/mol, CM2: 2.7 kcal/mol, and CM6: 21.4 kcal/mol), no direct coupling is detected between S@SO 2 and Al atoms. When γ-Al 2 O 3 (110) surface is hydrated, five stable adsorption configurations arise that are distinct from dehydrated condition: two physisorption modes (HM1 and HM2) and three chemisorption modes (HM3, HM4, and HM5), see Figure 4. HM1 and HM2 are structurally similar in that their S atoms are coordinated to a surface hydroxyl, while the coordination numbers of their Al atoms are different from each other. HM3 is produced by interaction of O@SO 2 atom with five-fold Al sites. HM4 and HM5 contain the sulfite species where the S atom is coordinated with surface-O atoms. The adsorption energies are calculated to be À20.4, À25.3, À31.1, À17.5, and À35.0 kcal/mol, respectively, for HM1, HM2, HM3, HM4, and HM5. In consequence, HM5 with formation of the sulfite species represents the lowest-energy adsorption configuration, which is the same as in dehydrated condition (CM8). Two IR peaks at 1214 and 1349 cm À1 are assigned to the sulfate species, which can be been finely interpreted by DFT calculated results.
Goethite (α-FeO(OH)), which can be found in soils and other low-temperature environments, is an iron-bearing hydroxide. Because of the considerable adsorption capacity for organic acids and anions, goethite has also been widely used in environmental remediation and protection [35]. Zubieta et al. [37] investigated the adsorption of SO 2 on partially and fully hydrated (110) surfaces of goethite and obtained eight stable products: six sulfite, one bisulfate, and one sulfate ( Figure 5). The adsorption of SO 2 onto the Cu(100), MgO and carbon surfaces was discussed as well, which may provide insightful clues for resembling processes onto mineral surfaces. It was proposed that SO 2 and H 2 O are co-adsorbed onto Cu(100) surfaces [30], through the direct coupling of Cu atom with S and O@H 2 O atoms. The Cu-O distances ascend in the order of coadsorption of SO 2 and H 2 O < adsorption of only SO 2 < < adsorption of only H 2 O. Accordingly, the interaction between SO 2 and Cu(100) surface is stronger than that of H 2 O, and the coadsorption of SO 2 conduces to the enhanced interaction of H 2 O with Cu(100) surface. At the same time, the Cu-S distance of the co-adsorption configuration is optimized at 2.385 Å and is shorter than that with only SO 2 adsorption. That is, water exhibits a promoting effect for the adsorption of SO 2 on Cu(100) surface, as corroborated by the calculated adsorption energies. Eid and collaborators [41] found that as compared to regular MgO surface, the adsorption capacity of SO 2 at MgO(Fs-center) defects is higher, and MgO(Fs-center) corresponds to an enhanced catalytic activity. With regard to pure carbon materials, SO 2 is physisorbed and van der Waals (vdW) is the driving force therein [42]. When carbon materials are modified with functional groups such as carboxyl, lactone, or/and phenolic hydroxyl, the adsorption strengths of SO 2 are enhanced pronouncedly, especially for the sites at edge surfaces. In addition, these functional groups show little effects on SO 2 adsorption, suggesting that the enhanced adsorption is mainly due to regulation of carbon surface properties. Adsorption of SO 2 is the first step for the formation of acid rain. According to our preliminary studies, the reaction mechanisms of gas phase and mineral surfaces resemble each other, and hence the gas-phase results are beneficial to understand acid rain formation at mineral surfaces. The reaction of SO 2 and H 2 O produces two isomers that have close electronic energies [43]; in addition, the two isomers have apparently lower electronic energies than H 2 SO 3 . Accordingly, the gaseous SO 2 and H 2 O mixture is likely to exist as the SO 2 ÁH 2 O complex rather than H 2 SO 3 . The activation barrier of SO 2 reacting with H 2 O to form H 2 SO 3 is so high (146.7 kJ/mol) that it becomes very difficult to produce the sulfite species (H 2 SO 3 ) in gas phase. Five years later, Stirling [44] investigated the hydrolysis of SO 2 in aqueous solutions, finding that hydrated SO 2 forms the bisulfite anion and hydronium ion after overcoming an energy barrier of about 17 kcal/mol, while the one-step formation of H 2 SO 3 has not been detected. The orientation of water molecules in the hydration shell of SO 2 implies a more facile formation of the bisulfite anion rather than H 2 SO 3 [45], in line with the results of meta-dynamics calculations [46]. When HO 2 participates in the reaction, the S atom constructs a new bond with O@H 2 O atom (S-O: 1.716 Å). The energy barrier reduces considerably and equals 56.6 kJ/mol. HO 2 exists widely in the atmosphere and participates in a variety of chemical reactions [47]. In addition to HO 2 , a number of computational studies explicitly indicated that acidic substances in the atmosphere play similar catalytic effects and reduce significantly the energy barriers for the hydrolysis of SO 2 [48][49][50]. Interaction between SO 2 and the OH radical (OH • ) has been addressed by ab initio electronic structure calculations [52], and the optimized structure of the HOSO 2 radical show differences with that predicted using a small basis set (MP4/6-31G**//HF/3-21G*) [53]: Inclusion of correlation effects results in an increase of S-O distances and OSO angle, whereas the HOS angle shows a substantial decrease. Although geometrically stable at all levels of theory, stability of the HOSO 2 radical is estimated to be 104-110 kJ/mol, suggesting that the direct dissociation to SO 3 is almost infeasible. O 3 and H 2 O 2 are two oxidizing substances in atmosphere chemistry, and their reaction mechanisms with SO 2 were studied by Jiang et al. [54].  3 . This is one-step process and requires to overcome an energy barrier of 46.5 kJ/mol. The reaction of SO 2 with H 2 O 2 proceeds via the OH-abstraction mechanism, and two OH radicals generated from the dissociation of H 2 O 2 are appended to SO 2 forming H 2 SO 4 . However, the energy barrier of this reaction is 299.5 kJ/mol, which is extremely difficult to proceed at normal conditions. Chen et al. [46] reported the reaction mechanism between SO 2 and HO 2 and showed that there exist two types of SO 2 ÁHO 2 complexes: one is to combine the terminal O@HO 2 and S atoms. In this reaction, the O-S distance decreases from 2.966 Å (the initial complex: SO 2 ÁHO 2 ) to 1.598 Å while the O-O bond of HO 2 is gradually elongated until broken. The reaction of SO 2 and H 2 O is divided into three stages, and the energy barrier of SO 3 ÁOH formation is so small (1.7 kJ/mol) that can be considered negligible. The final product of this reaction is HSO 4 . The other complex is characterized by two types of intermolecular H-bonds that form between the terminal O@HO 2 and S atoms and between the H@HO 2 and O@SO 2 atoms. It is a one-step process and the product for this reaction is HOSOÁO 2 . When a single water is added, it has a combined effect on the overall process, not only accelerating the reaction of generating HSO 4 by reducing the activation barrier of the second step from 48.8 to 44.8 kJ/mol but also inhibiting the reaction between SO 2 and HO 2 to produce HOSO due to blocking the interaction between the H@HO 2 and O@SO 2 atoms.
It can be seen from the above discussions that the SO 3 and water reaction is an integral section for the formation of acid rain. The direct reaction of SO 3 with one water to produce sulfuric acid (H 2 SO 4 ) requires to overcome a large energy barrier (28.7 kcal/mol) and seems difficult to occur at normal conditions [55][56][57]. When the second water participates, the energy barrier reduces substantially and the reaction becomes almost barrierless, which is probably due to the formation of a stable six-membered cyclic transition state, see Figure 6 [58]. The two water molecules transfer their protons in a concerted manner, and a new S-O bond is formed between the newly generated OH and SO 3 fragments. In consequence, multiple reaction pathways may co-exist as illustrated below

Concluding remarks
In contrast to the rapid development of chemistry, physics, and biology and other disciplines, soil science remains almost stagnant in the past few decades, and to best of our knowledge, no breakthroughs have been reported for rather a long time. Despite that, no one can deny the vital significance of soil to our life, and soil science has been widely acknowledged as the final frontier.
The slow progresses for soil science, we think, should be attributed to two reasons: (1) wrong and outdated perceptions. A majority of soil researchers mistakenly believe that all knowledge worth knowing about soils has already been understood and no revolutionary progresses would take place; (2) complex systems. Soils are very structurally complicated and there are multiple factors to co-function, which makes it very difficult to characterize by experimental techniques. Computer simulations have unique advantages to handle complex systems, while currently its role in soil science is far from being recognized. In this chapter, two examples are elaborately discussed with regard to application of DFT calculations in soil science: one focuses on the adsorption of uranyl onto mineral surfaces, and the other involves the adsorption of SO 2 onto mineral surfaces and reaction mechanisms to form acid rain. It can be seen from these discussions that DFT calculations are able to provide useful and detailed information about the adsorption, interaction and reactions at the atomic level that greatly promote our understanding about soil science.
With advent of high-performance computing platforms, the same DFT calculation tasks of 10 years ago can now finish within a remarkably shorter time, even if you increase model size (periodic model is also an option), consider solvent effects by adding explicit solvent molecules or/and choose more accurate methods. The methodological developments regarding to DFT calculations have also made remarkable progresses over the recent three decades, and as a result, thermodynamics and reaction barriers can now be predicted with nearly chemical accuracy (≤1 kJ/mol). In consequence, computer simulations including DFT are the right key to unravel the complicated phenomena and processes occurring within soils; e.g., with DFT calculations, the aggregation mechanisms of "real" soils and the driving force therein were unveiled at an atomic level [6].
In addition to DFT methods, there are a number of other computational methods, such as QM/ MM and Molecular Dynamics (MD). The choice of suitable computational methods is strongly recommended. We are pleased to see the birth of the ClayFF force-field [66] that was developed specially for clay minerals and the capability of the ReaxFF force-field [67] to handle reaction mechanisms. The 2013 Noble Prize in Chemistry was awarded to the development of "multiscale methods for complex systems," and now it is time to apply these methods to tackle complex soil systems.