Integrating Density Functional Theory Modeling with Experimental Data to Understand and Predict Sorption Reactions: Exchange of Salicylate for Phosphate on Goethite

: Density functional theory (DFT) calculations are a quantum mechanical approach that can be used to model chemical reactions on an atomistic scale. DFT provides predictions on structures, thermodynamics, spectroscopic parameters and kinetics that can be compared against experimentally determined data. This paper is a primer on the basics of utilizing DFT for applications in mineral-water interfaces. In our case-study, we use DFT to model the surface complexes of phosphate and salicylate adsorbed onto the (101) and (210) surfaces of α -FeOOH (goethite), as an example of combining DFT and experiment. These three components are important in the phosphorus-organic matter interactions in soils, and by comparing the energies of the two surface complexes, the exchange energy of salicylate for phosphate onto goethite can be estimated. The structures of the surface complexes are predicted and the resulting vibrational frequencies calculated based on these structures are compared to previous observations. Upon veriﬁcation of reasonable surface complex models, the potential energy of exchanging salicylate for phosphate is calculated and shown to be signiﬁcantly exothermic. This model result is consistent with observations of plant exudates, such as salicylate freeing adsorbed phosphate in soils under P-limited conditions.


Introduction
Scientists studying the Earth's processes are often interested in the macroscopic chemical behavior of a field site or region. For example, the problem of eutrophication of water bodies can be caused by nutrient runoff from agricultural fields or concentrated animal feeding areas [1]. The transport of N and P from soils into rivers, lakes and bays can be monitored via sample collection and models (e.g., SWAT [2]; and pedon-scale behavior is ultimately influenced by pore scale (millimeters to microns), nanoscale, and molecular-level processes. Consequently, molecular-level mineral-water interface chemistry must be understood in a manner that can transfer information to larger scales in order to apply it at the field scale.
Equilibrium and reaction rate constants are two critical parameters that affect the transport and fate of elements and compounds in soils and aquifers and have molecular-level origins. In addition, knowing aqueous and surface speciation is imperative to connect the model parameters to the chemical species present in a given environment. Otherwise, the results are empirical fitting models that will be less useful for prediction and designing intelligent management practices.
An excellent example of this molecular-level approach is Hanna et al. [3]. The authors used batch and flow-through adsorption experiments and attenuated total-reflectance Fourier-transform infrared spectroscopy (ATR FTIR) to observe phthalic acid (C 8 O 4 H 6 ) adsorption onto goethite (α-FeOOH)-coated sand. Phthalic acid is a naturally occurring compound and representative of functional groups present in fulvic and humic acids. Goethite-coatings are common on mineral surfaces due to the insoluble nature of Fe 3+ in most aerobic soil and natural water systems, and these coatings are active adsorbents, especially for carboxylic acids and oxyanions [4]. Adsorption was measured as a function of time to ensure equilibration and as a function of pH. ATR FTIR revealed the relative concentrations of covalent and H-bonded species of phthalate with the goethite. Surface complexation models were constructed and the results were applied to predict breakthrough curves in flow-through experiments. The authors concluded that the concentrations of non-specifically adsorbed species are a function of pH. Such non-specific (i.e., H-bonded or outer-sphere complexes) have more labile kinetics than covalently bonded surface complexes, so knowledge of the extent of both non-specific and covalent bonding is necessary in order to predict macroscopic adsorption/desorption behavior. The authors point out that these results indicate that reactive transport models need to be revised to include this complex chemistry.
The concentration of P in soil solution, the dominant factor controlling the bioavailability of soil P [5,6], is determined in many soils by adsorption and desorption processes of P onto metal (oxyhydr)oxides. These minerals, with their preponderance of surface hydroxyl groups, have high affinity sorption sites for both P and OM [7,8]. Thus, the bioavailability of P is likely to be sensitive to both the level of OM in the soil and the surface properties of the metal (oxyhydr)oxides (i.e., surface coverage with respect to both P and OM). The presence of organic acids, such as root exudates, in soils has been shown to increase the bioavailability of soil P [9,10]. Hue [11] reported that three low molecular weight organic acids decreased P sorption by acidic soils and increased the yields of a lettuce (Lactuca sativa L.) bioassay crop. Bolan et al. [12] reported that addition of organic acids typically found in soil, leaf litter, and poultry manure decreased P sorption in soils. These findings suggest that root exudate-derived organic acids may increase the bioavailability of soil P fertilizer.

Macroscopic Data-Adsorption Isotherms, pH Edges, Calorimetry, Adsorption/Desorption Kinetics
Macroscopic data such as extent of adsorption as a function of solution adsorbate concentration or pH, calorimetry of the adsorption reaction, and adsorption/desorption rates are the primary constraints on the molecular-level chemistry that occur in soils and aquifers. The aqueous/surface partition constant, the enthalpy of adsorption, and the adsorption/desorption rates determine the observed field-scale behavior of elements and compounds in soils. However, these macroscopic techniques neither determine nor require knowledge of the molecular species involved. For example, aqueous silica can be represented thermodynamically as SiO 2(aq) or Si(OH) 4 , but the chemistry of these two as written would be different. Another instance is that the toxicity of free metal ions such may be orders of magnitude higher than the same metal ions in organic complexes [13][14][15]. Kinetics and biological impacts depend heavily on the molecular speciation, so informing the larger scale models with the appropriate species is desirable.

Molecular Data-IR/Raman, EXAFS/XANES
Molecular speciation can be determined with a variety of methods, and generally it is a best practice to employ more than one method when practical. Each type of spectroscopy has strengths and weaknesses, so a combination of techniques can provide for a more complete picture of chemical bonding. Vibrational spectroscopies such as IR and Raman are sensitive to the bonds in a compound that undergo changes in their dipole moments and polarizabilities, respectively [16]. Hence, IR and Raman detect changes in bond strengths and the types of bonds present. Nuclear magnetic resonance spectroscopy (NMR) is sensitive to the electronic environment around nuclei with an odd number of spins (e.g., 13 C, 27 Al, 29 Si) [17]. Thus, the coordination environment around selected nuclei can be probed. Extended X-ray adsorption fine structure and X-ray adsorption near-edge spectroscopies (EXAFS and XANES) can be used to determine interatomic distances and electronic orbital energy Soil Syst. 2020, 4, 27 3 of 14 levels, respectively, that are important in soil chemistry [18]. Bonding environments and oxidation states are derived from this type of data.
Each technique may not capture all the important aspects of the molecular structure. IR and Raman do not detect all vibrational modes and many may be weak or result in frequencies that overlap with other vibrational modes, making interpretation problematic. NMR works with many nuclei but not all, and the presence of Fe in the sample can disrupt NMR spectra [19,20]. EXAFS is insensitive to lighter elements, so H positions cannot be identified. Changes in the electronic orbital energies may be due to various effects, so the interpretation of XANES spectra can be ambiguous. To complement these spectroscopies, molecular simulations are commonly employed [21].

Classical Molecular Mechanics Simulations
Many studies of soil mineral surface chemistry employ classical molecular mechanics (see [22] and references therein). By classical, we mean that Newtonian mechanics are applied to the atoms and molecules in the system. Each atom is represented by a series of parameters, collectively termed a "force field" or "interatomic potential". Most parameterizations include an atomic charge (formal or otherwise), atomic repulsion parameters, so that cations and anions do not attempt to undergo fusion, van der Waals terms, and bonding parameters, such as harmonic force constants for bond stretches and angle bends [23].
The advantage of classical molecular mechanics is that model systems of millions of atoms can be simulated and the timeframe for these the simulations may be on the millisecond scale. For complex biogeochemical systems, these are highly desirable characteristics. The limitations of the classical approach are mainly that the accuracy depends on the parameterization and that most parameterizations focus on reproducing structures and energies near minima. This results in higher uncertainties when describing the tails of the Boltzmann distribution of configurations for a system and commonly means that bond-breaking and bond-making are not allowed during the simulation. Reactive force fields such as ReaxFF [24] exist and create a promising avenue for simulating complex chemistry when the accuracy of the parameterization is high [25].

DFT
Another computational approach is to use quantum mechanics to calculate properties of the system of interest. The main limitations are that models of 1000 atoms stretch the computational capacity of today's supercomputers and that molecular dynamics simulations are limited to hundreds of picoseconds, rather than the milliseconds possible with classical methods. The extra computational demand to calculate the system's electron density distribution as opposed to using simple analytical expressions for bonding decreases the spatial and temporal domains one may explore.
However, even within these constraints, the quantum mechanical approach, most often using density functional theory (DFT), can provide insights into soil chemistry. A key is to utilize the experimental and spectroscopic data to constrain the possible structures one needs to examine and define questions that DFT can shed light upon. Because the DFT approach provides electronic and molecular structures, thermodynamic, kinetic and spectroscopic properties can all be predicted from the calculations. The DFT method is not dependent upon fitting parameters to match the experimental data, so when the DFT results are consistent with experimental observables, one can be reasonably certain that the resulting DFT molecular model is reliable. This also gives the researcher the ability to estimate the error in the computational results, so that one can distinguish among possibilities, or recognize that two results are equivalent within error. However, we note that the DFT models are calculating ∆E exch values, not ∆G exch . The calculation of ∆G values is possible (e.g., [26]) but requires long duration MD simulations that sample configuration space to an extent not currently practical with DFT on systems of this size. Entropy can be significant, because ∆S for desolvation reactions of ions in solution are 70 J/mol-K [27]. This contribution to ∆G would be (70 J/mol-K)(298K) ≈ 21 kJ/mol Soil Syst. 2020, 4, 27 4 of 14 at 298K. We rely on the hypothesis that the ∆S exch is relatively small compared to ∆E, because this is an exchange reaction. The larger ∆E exch is, the more likely this is to be true. This paper will illustrate a case-study of combining DFT with experimental data to better understand an important soil chemistry reaction-the exchange of organic acids for phosphate on goethite.

Adsorption and Exchange Experiments
A batch study was conducted to determine the salicylate adsorption isotherm onto α-FeOOH (goethite). Thirty mL of solutions containing 0, 40, 80, 120, 160, 320, 480, and 640 mg L −1 salicylate adjusted to pH 5.0 with a background matrix of 0.02 mol L −1 KCl were added to 0.50 g of goethite in 50-mL centrifuge tubes. The tubes were placed on an orbital shaker for 4 h at 4 • C, to minimize microbial degradation of the metabolite and vacuum filtered through a 0.4 µm polycarbonate filter. The total dissolved organic carbon was determined using Shimadzu TOC-Vcsh Analyzer and the quantity of salicylate adsorbed was determined by subtracting the salicylate concentration in the filtrate from the initial salicylate concentration. The salicylate adsorption was fit to the Langmuir model: where Q is the quantity adsorbed (mg kg −1 ), C is the equilibrium concentration (mg L −1 ), Q max is maximum monolayer adsorption, and K is the binding strength constant. The non-linear regression fitting was conducted with MATLAB to obtain the best-fit values of Q max and K.

DFT
For a review of basic quantum mechanics and density functional theory, see books such as [28] or a review chapters such as [29]. This section will describe some important practical aspects about using DFT to model mineral surface-mediated reactions. This author primarily utilizes Materials Studio (Biovia Inc.) to build and visualize models, the Vienna Ab-initio Simulation Package (VASP; [30][31][32][33][34] for periodic calculations, and Gaussian 16 [35] for vibrational frequencies and NMR chemical shieldings on clusters. However, there are numerous other codes, such as Crystal [36], Siesta [37], Quantum Espresso [38], and CASTEP [39], that can be used depending on the system under study and user preferences. Desktop computers and laptops are capable of running DFT calculations, but computational and memory limits mean that larger scale models are impractical. Generally, one would prefer to run DFT calculations on a Linux cluster using the order of tens to thousands of processors, to model systems of hundreds to thousands of atoms.

Basis Sets
A primary factor controlling the accuracy of the DFT results is the choice of basis set or energy cutoff used in the calculation. The basis of the DFT approach is to mimic the actual electron density distribution of the system, so the larger the set of functions used to model the electrons, the more accurate the calculation can be. Unfortunately, increasing the basis set size or energy cutoff used also increases computational demand. Thus, one must be judicious in finding the most efficient method that can produce the accuracy necessary for the problem at hand. When there are large differences between observables, then less model accuracy is acceptable. For example, if a wavenumber shift in a vibrational frequency is over 100 cm −1 , then an error or ±10 cm −1 is not a major concern. If subtle changes such as a few wavenumbers are being modeled, it may be difficult to find any method that can discern such differences. This author recommends running tests with increasing basis set size or energy cutoff on small test systems similar to the system of interest. When results converge, the most computationally efficient method can be selected for modeling larger, more complex systems of interest.

Exchange-Correlation Functionals
The "functional" components of DFT are the electron exchange and correlation functionals developed to describe these quantum phenomena that go beyond the electrostatic interactions of electrons with nuclei and with other electrons. New functionals are continually being developed to increase accuracy, as judged by the reproduction of experimental data. Many are developed with dispersion corrections that account for van der Waals forces between atoms that are not accounted for in DFT. The investigator may need to run method testing calculations with a matrix of basis sets/energy cutoffs and exchange-correlation functionals before diving into the question of interest; however, it would be wise to check the literature for such studies first, because this is a time-consuming project. An excellent example is the study of Demichelis et al. [40], who examined exchange-correlation functional methods for accuracy in reproducing silicate structures.

Model Construction
In parallel with finding a practical and accurate computational methodology, care must be given to the construction of the model system to account for all the factors possible that may affect the results. Final states of molecular modeling can be heavily influenced by the initial configuration constructed (i.e., metastable states may be predicted that are similar in structure to the starting condition). Furthermore, the farther the initial state is from the true structure, the longer the calculations will take to achieve the final result. Consequently, spending time to build a model or models based on as much experimental data as possible is worthwhile. Generating a mineral surface may appear to be a trivial exercise as a plane along a Miller index can be used to cleave the known bulk crystal structure. Where this plane is placed within the crystal structure and any possible reconstruction should not be left to automation within a program, however. An example of where this issue was critical can be found in Lo et al. [41], where construction of the hematite (α-Fe 2 O 3 ) (1-102) surface required study in its own right. Similar issues can exist with complex mineral structures such as goethite [42].
Solute and adsorbate initial structures may be guided by experimental data on coordination numbers and bond lengths in the aqueous state (see [43] for example). Although surface complexes may deviate significantly from aqueous coordination environments [44], starting from a known solvation structure is a reasonable choice. Many early calculations related to mineral surface chemistry avoided the presence of water altogether, mainly for computational efficiency. Today, many researchers use polarized continuum models for solvation (e.g., IEFPCM- [45] and references therein; COSMO- [46][47][48]) to save computer time. Continuum or "implicit solvation" models are an improvement over gas-phase calculations and can provide accurate results; however, implicit solvation does not account for explicit H-bonding effects. When H-bonding is strong, it can have significant effects on modeled structures, energies, vibrational frequencies and NMR chemical shifts [21].
Adding H 2 O molecules to the model reaction box can be problematic. One generally determines the volume remaining in the system after the mineral slab, solutes and adsorbates have been accounted for, then adds H 2 O molecules to create a density desired (typically ≈ 1 g/cm 3 ). As mentioned above, the initial configuration will influence the final result, so creating an initial guess for the water structure is a step that creates uncertainty. This author creates the slab-adsorbate-solute model based on experimental information, freezes these atoms, adds H 2 O to the appropriate density, and then runs classical energy minimizations (EM) and molecular dynamics (MD) simulations to create a reasonably H-bonded network for a starting configuration. Most force fields do not accurately reflect water structure for complex systems with surfaces and solutes, but the configurations can be realistic enough that the DFT calculations will optimize to a configuration close to reality. All these assumptions and approximations must be tested in each case against experimental data on the system of interest to ensure agreement within the needed accuracy. When large discrepancies exist, the procedure must be revisited.
Periodic models can be created to realistically represent mineral surfaces and solution interfaces, but often it is useful to create molecular clusters extracted from these periodic models. Selecting the atoms of interest around the surface complex (including H 2 O molecules in the solvation sphere), one can create a cluster to subject to vibrational and NMR calculations [49].

Energy Minimization and Molecular Dynamics
After an initial model is created, the atomic structure is subjected to energy minimization (EM) and/or molecular dynamics (MD) calculations. MD has commonly been used with classical force fields, so many non-modelers equate MD with classical methods. This is not the case, however. Both EM and MD can be applied to classical and quantum calculations. The main difference is that the quantum calculations are typically smaller models and run for shorter durations, as mentioned above. EM is used to determine the most stable structure of a system. Due to the presence of local energy minima for complex systems, finding the global energy minimum is problematic. A practical approach is to model the system in a variety of possible configurations hypothesized from experimental data, and then see which model has the lowest predicted energy and best matches experimental observables.
Even if one determines a global minimum using EM, real systems at finite temperature sample configurations around minima and metastable minima if they exist without crossing prohibitive activation energy barriers. In these cases, MD is useful for exploring configuration space and sampling structural variations around the minimum. MD is run at a finite temperature so the atoms take on a Boltzmann distribution of energies [50] and can overcome potential energy barriers that are small compared to the kinetic energy available. MD is run with small time steps (e.g., 0.1 to 1 × 10 −15 s), so the duration of these simulations in DFT-MD may be on the order of tens of picoseconds (1 × 10 −11 s). The modeler must recognize that this is a severe constraint on sampling the system, so comparison of the final results to experimentally observable quantities is encouraged wherever possible.

Frequency Analysis
MD simulations can be used to calculate vibrational densities of states (VDOS) through the velocity autocorrelation function [51]. The VDOS is not directly comparable to IR and Raman spectra, however, because it does not include the IR and Raman intensity of each vibrational mode. Therefore, this author prefers to use energy-minimized structures and subject them to harmonic frequency analysis [52]. Anharmonic frequency calculations are possible, and in some cases necessary; however, often the faster harmonic frequency approximation and correction for DFT method and anharmonicity is sufficient to interpret spectra. Gaussian 16 [35] has the capability of calculating the IR and Raman intensities of each vibrational frequency, as well as the vibrational mode associated with each frequency (see [53] for more details). The comparison of modeled with observed frequencies is straightforward, so discerning which model best matches experiment is a relatively simple exercise. Vibrational frequencies can also be calculated for periodic systems, but IR and Raman intensities are not always produced. In addition, the numerical procedures for calculating vibrational frequencies in large periodic systems may make this approach impractical. In this paper, we will show how periodic models are used to generate molecular cluster models that are then subjected to frequency analyses.

Model and Computational Details
In this study, the (210) and (101) surfaces of goethite were cleaved from the experimental crystal structure using the "Cleave Surface" tool of Materials Studio 2016 (Biovia Inc., San Diego, CA, USA) The unsatisfied valences of the surface Fe atoms created by this cleavage were satisfied by adding H 2 O molecules, such that Fe-OH 2 bonds completed the octahedral coordination, except for one surface Fe atom, where the valence was satisfied by an O atom of the phosphate or salicylate ligand. The goethite slab was three Fe layers' thick, to create a symmetric representation of the surface. A vacuum space of 2 nm was added to the slab to provide space for the phosphate and salicylate groups, as well as solvating H 2 O molecules. The resulting periodic simulation cells were 9.24 × 11.63 × 30.26 Å for After construction of the model, the atoms in the goethite, phosphate and salicylate were fixed and the H 2 O molecules were subjected to energy minimization and 100 ps MD simulation at 300K, using the universal force field [54], and energy re-minimization, to obtain an initial structure for energy minimization using VASP [30]. The initial structure of the goethite (210)-phosphate surface complex was based on the structure produced in [49]. The salicylate structure was taken from previous calculations on salicylate by Trout and Kubicki [55]. These initial structures were converted to VASP input (POSCAR) format via a perl script written by A.V. Bandura (St. Petersburg State University). Energy minimizations in VASP were then performed using the PBE functional, a 500 eV energy cutoff and the Grimme D2 dispersion correction [56]. A single k-point at the point in the Brillouin zone was used with the Monkhorst-Pack scheme, and energy convergence criteria of 1.0 × 10 −7 eV and −0.02 eV/Å were used for the electron density and geometry convergence cutoffs, respectively.
The molecular cluster models were extracted from the energy-minimized periodic models by selecting the surface complex atoms of interest (i.e., the surface Fe, phosphate or salicylate); the O atoms bonded to these atoms and any H 2 O molecules within 2 Å of these atoms to account for H-bonding. This was done for both the surface and aqueous species. These molecular clusters were then subjected to energy minimization and frequency analyses in Gaussian 16 [35]. The B3LYP/6-31G(d,p) exchange-correlation functional and basis set were used [57,58]. The frequencies were scaled by 0.96 to account for basis set and anharmonicity effects to compare to observed frequencies [49].

Adsorption Data
The adsorption isotherm relating the quantity of salicylate adsorbed, as a function of the equilibrium concentration in solution, is shown in Figure 1. The isotherm fits the L-type classification where the slope decreases with increasing salicylate concentration in the equilibrium solution [59]. The decreasing high initial slope characteristic of the isotherm reflects the high affinity of salicylate for the goethite surface, coupled with decreasing adsorptive sites with increasing adsorption [60]. The salicylate adsorption maximum on the goethite was calculated using a one-site Langmuir model was calculated to be 8641 mg kg −1 goethite. The affinity of organic acids, such as salicylate, for goethite can be due to "non-specific" electrostatic interactions or through "specific" interactions involving the formation of covalent bonding between the organic molecule and Fe atom of the mineral. Computational chemistry and FT-IR spectroscopy can be used to help establish the nature of the bonding mechanism present in the system under investigation [61].

Energetics
Adsorption to minerals will be controlled by the type of surface groups present which, for a given mineral, is controlled by the habit of the particle (i.e., the crystal faces present). Some faces are likely to be more reactive towards certain compounds than others. This was demonstrated clearly by Villalobos and coworkers [62,63] for the adsorption of various ions onto goethite. Based on this work and a previous DFT modeling study of phosphate adsorption onto goethite [49], two of the more reactive faces are the (210) and (101) surfaces. Consequently, rather than modeling all possible goethite surfaces, this study reports on phosphate and salicylate adsorption on goethite (210) and (101) only.
Previous studies have used cluster models to estimate adsorption enthalpies (e.g., [64], and others have used periodic models of a single adsorbate in various configurations (e.g., [49]). Periodic models should be more realistic in terms of modeling the adsorption reaction as it actually occurs. Compared to clusters, the periodic model surface is a better representation of real surfaces, and the solvation of compounds in solution is more complete. However, periodic DFT calculations are not without their limitations. A major issue is the calculation of the potential energy differences without calculating Gibbs free energies of adsorption (ΔGads). To circumvent this issue, this study modeled systems, with the salicylate and phosphate included simultaneously. Thus, the relative energy differences reflect the surface bonding and solvation, and the ΔS terms should be minimized, because they will tend to cancel each other out when the two configurations are compared.
The four models for the (210) and (101) surfaces with phosphate and salicylate adsorbed and in solution respectively are shown in Figure 2. The (210) model consists of 24 FeOOH, HPO4 2− , C7H5O3 − (salicylate) and 68 H2O and 3 H + . Either the HPO4 2− or the salicylate are bonded to the surface in a monodentate configuration with the other anion in solution. The (101) model consists of 24 FeOOH, HPO4 2− , C7H5O3 − (salicylate) and 57 H2O and 3 H3O + , with the anions as bidentate bridging innersphere complexes. The choice for phosphate surface complex bonding was guided by the lowest energy configurations found in Kubicki et al. [49]. The small number of H2O molecules means that the anion concentrations are high and the pH low compared to a real system, but the current computational facilities make increasing the number of atoms in the system substantially impractical. Similar simulations on larger scale systems using an accurate, reactive force field could address this issue.

Energetics
Adsorption to minerals will be controlled by the type of surface groups present which, for a given mineral, is controlled by the habit of the particle (i.e., the crystal faces present). Some faces are likely to be more reactive towards certain compounds than others. This was demonstrated clearly by Villalobos and coworkers [62,63] for the adsorption of various ions onto goethite. Based on this work and a previous DFT modeling study of phosphate adsorption onto goethite [49], two of the more reactive faces are the (210) and (101) surfaces. Consequently, rather than modeling all possible goethite surfaces, this study reports on phosphate and salicylate adsorption on goethite (210) and (101) only.
Previous studies have used cluster models to estimate adsorption enthalpies (e.g., [64], and others have used periodic models of a single adsorbate in various configurations (e.g., [49]). Periodic models should be more realistic in terms of modeling the adsorption reaction as it actually occurs. Compared to clusters, the periodic model surface is a better representation of real surfaces, and the solvation of compounds in solution is more complete. However, periodic DFT calculations are not without their limitations. A major issue is the calculation of the potential energy differences without calculating Gibbs free energies of adsorption (∆G ads ). To circumvent this issue, this study modeled systems, with the salicylate and phosphate included simultaneously. Thus, the relative energy differences reflect the surface bonding and solvation, and the ∆S terms should be minimized, because they will tend to cancel each other out when the two configurations are compared.
The four models for the (210) and (101) surfaces with phosphate and salicylate adsorbed and in solution respectively are shown in Figure 2. The (210) model consists of 24 FeOOH, HPO 4 2− ,  [49]. The small number of H 2 O molecules means that the anion concentrations are high and the pH low compared to a real system, but the current computational facilities make increasing the number of atoms in the system substantially impractical. Similar simulations on larger scale systems using an accurate, reactive force field could address this issue.
ΔEads values should still dominate the ΔGads, because the ΔSads is expected to be significantly smaller as they largely cancel out when the salicylate-phosphate exchange occurs. Figure 3 illustrates how the energy changes for each configuration from the initial to the final step when the system is at a potential energy minimum. The curves show the convergence for each model and the discernible differences in energies between the phosphate and salicylate surface complexes.   Nevertheless, the DFT results are consistent with observation of salicylate replacing phosphate on goethite. For the (210) monodentate configuration, a ∆E ads = −93 kJ/mol is predicted and for the (101) bidentate bridging model, ∆E ads = −45 kJ/mol. As mentioned above, entropic factors could affect the equilibrium for salicylate-phosphate exchange on goethite, but these relatively large negative ∆E ads values should still dominate the ∆G ads , because the ∆S ads is expected to be significantly smaller as they largely cancel out when the salicylate-phosphate exchange occurs. Figure 3 illustrates how the energy changes for each configuration from the initial to the final step when the system is at a potential energy minimum. The curves show the convergence for each model and the discernible differences in energies between the phosphate and salicylate surface complexes. Nevertheless, the DFT results are consistent with observation of salicylate replacing phosphate on goethite. For the (210) monodentate configuration, a ΔEads = −93 kJ/mol is predicted and for the (101) bidentate bridging model, ΔEads = −45 kJ/mol. As mentioned above, entropic factors could affect the equilibrium for salicylate-phosphate exchange on goethite, but these relatively large negative ΔEads values should still dominate the ΔGads, because the ΔSads is expected to be significantly smaller as they largely cancel out when the salicylate-phosphate exchange occurs. Figure 3 illustrates how the energy changes for each configuration from the initial to the final step when the system is at a potential energy minimum. The curves show the convergence for each model and the discernible differences in energies between the phosphate and salicylate surface complexes.

Vibrational Frequencies
In order to verify that the model surface complexes represent the configurations present on goethite, molecular clusters representing the short-range bonding of these surface complexes were created and subjected to energy minimization and frequency analyses using Gaussian 16 [35] (Frisch et al., 2016). The extracted molecular cluster method results in similar structural parameters between the periodic and molecular models studied in Paul et al. [65]. In this case, P-Fe distances are 3.25 Å in the bidentate bridging periodic model and 3.13 Å in the molecular model consistent with observation (3.2 to 3.3 Å; [66]). For the monodentate model of the (210) surface complex, the P-Fe distance is 3.15 and 3.14 Å in the periodic and molecular cluster, respectively, however, which is significantly shorter than the interpreted EXAFS result of 3.6 [66]. These authors do state "The shortest P-Fe distance of 2.83-2.87 Å was indicative of a bidentate mononuclear bonding configuration." So, it may be that monodentate configurations can have shorter P-Fe distances than commonly thought [67]. This method has worked well in the past to produce good correlations of calculated and observed frequencies, but the comparisons can be complicated due to the possibility that numerous surface complexes form giving rise to various vibrational modes [49]. Nonetheless, model complexes that give rise to vibrational frequencies similar to those observed in IR spectra are likely to be present at some concentration in the samples.
After correcting the calculated harmonic frequencies for anharmonicity and basis set effects using the 0.96 factor obtain from NIST (https://cccbdb.nist.gov/vibscale.asp), the correlations of modeled and observed vibrational frequencies in the range 800 to 1800 cm −1 are plotted in Figure 4. Correlation statistics are good for both the phosphate-and salicylate-goethite clusters. These results do not preclude other complexes (e.g., outer-sphere) existing on other surfaces of goethite, but they do suggest that the modeled complexes have a finite concentration on the surface.

Vibrational Frequencies
In order to verify that the model surface complexes represent the configurations present on goethite, molecular clusters representing the short-range bonding of these surface complexes were created and subjected to energy minimization and frequency analyses using Gaussian 16 [35] (Frisch et al., 2016). The extracted molecular cluster method results in similar structural parameters between the periodic and molecular models studied in Paul et al. [65]. In this case, P-Fe distances are 3.25 Å in the bidentate bridging periodic model and 3.13 Å in the molecular model consistent with observation (3.2 to 3.3 Å; [66]). For the monodentate model of the (210) surface complex, the P-Fe distance is 3.15 and 3.14 Å in the periodic and molecular cluster, respectively, however, which is significantly shorter than the interpreted EXAFS result of 3.6 [66]. These authors do state "The shortest P-Fe distance of 2.83-2.87 Å was indicative of a bidentate mononuclear bonding configuration." So, it may be that monodentate configurations can have shorter P-Fe distances than commonly thought [67]. This method has worked well in the past to produce good correlations of calculated and observed frequencies, but the comparisons can be complicated due to the possibility that numerous surface complexes form giving rise to various vibrational modes [49]. Nonetheless, model complexes that give rise to vibrational frequencies similar to those observed in IR spectra are likely to be present at some concentration in the samples.
After correcting the calculated harmonic frequencies for anharmonicity and basis set effects using the 0.96 factor obtain from NIST (https://cccbdb.nist.gov/vibscale.asp), the correlations of modeled and observed vibrational frequencies in the range 800 to 1800 cm −1 are plotted in Figure 4. Correlation statistics are good for both the phosphate-and salicylate-goethite clusters. These results do not preclude other complexes (e.g., outer-sphere) existing on other surfaces of goethite, but they do suggest that the modeled complexes have a finite concentration on the surface.

Implications and Conclusions
DFT serves to complement the experimental data in this case-study. By providing an independent test on proposed species and reactions, DFT eliminates ambiguity in analyzing complex datasets regarding surface adsorption reactions. Because the DFT method results in predictions that can be directly compared to experimentally observed spectroscopic and thermodynamic data, the accuracy of the model results can be estimated. Limitations in the system size and length of molecular dynamics simulations using DFT hinder application of DFT to complex soil systems. However, deconstruction into simpler components in combined experimental/computational studies can provide insights into the molecular-level processes that are critical in soils.
The adsorption of phosphate by Fe and Al oxy(hydr)oxide minerals affects P solubility in soils [7,8]. These minerals also adsorb organic molecules, suggesting that competition between phosphate and organic acids affects soil P availability. Plants have evolved the ability to manipulate the chemical environment surrounding the root to enhance its ability to obtain the essential elements needed for growth, by expending up to 30% of its assimilated C in the release of root exudates [68]. It is difficult to experimentally evaluate individually the large suite of compounds identified as root exudates, because the compounds may occur at low concentrations and may have short half-lives, with microbes using them as a C source [69]. Our results show that DFT can be used to calculate the free energy of adsorption of root exudates known to compete with phosphate for adsorption to goethite. Computation methods provide the opportunity to predict the potential of plant metabolites to increase the bioavailability of soil P.