QM/MM molecular dynamics investigation of the binding of organic phosphates to the 100 diaspore surface
The fate of phosphorus (P) in the eco-system is strongly affected by the interaction of phosphates with soil components and especially reactive soil mineral surfaces. As a consequence, P immobilization occurs which eventually leads to P inefficiency and thus unavailability to plants with strong implications on the global food system. A molecular level understanding of the mechanisms of the P binding to soil mineral surfaces could be a key for the development of novel strategies for more efficient P application. Much experimental work has been done to understand P binding to several reactive and abundant minerals especially goethite (α-FeOOH). Complementary, atomistic modeling of the P-mineral molecular systems using molecular dynamics (MD) simulations is emerging as a new tool, which provides more detailed information regarding the mechanisms, nature, and strength of these binding processes. The present study characterize the binding of the most abundant organic phosphates in forest soils, inositol hexaphosphate (IHP) and glycerolphosphate (GP), to the 100 diaspore (α-AlOOH) surface plane. Here, different molecular models have been introduced to simulate typical situations for the P-binding at the diaspore/water interface. For all models, quantum mechanics/molecular mechanics (QM/MM) based MD simulations have been performed to explore the diaspore-IHP/GP-water interactions. The results provide evidence for the formation of monodentate (M) and bidentate (B) motifs for GP and M and as well as two monodentate (2M) motifs for IHP with the surface. The calculated interaction energies suggest that GP and IHP prefer to form the B and 2M motif, respectively. Moreover, IHP exhibited stronger binding than GP with diaspore and water. Further, the role of water in controlling binding strengths via promoting of specific binding motifs, formation of H-bonds, adsorption and dissociation at the surface, as well as proton transfer processes is demonstrated. Finally, the P-binding at the 100 diaspore surface plane is weaker than that at the 010 plane, studied previously (Ganta et al., 2019), highlighting the influential role of the coordination number of Al atoms at the top surface of diaspore.
INTRODUCTION<!>MOLECULAR MODELING APPROACH<!>RESULTS AND DISCUSSION<!>GP M motif<!>GP B motif<!>IHP M motifs<!>IHP 2M motif<!>Effect of surface saturation<!>SUMMARY AND CONCLUSIONS
<p>Phosphorus (P) is essential for plant growth and plays an important role in photosynthesis, energy storage, cell growth and many other plant processes. It has been pointed out that P scarcity could arise in near future (Cordell and Neset, 2014) and to cope with this situation there is a need to understand the P cycle in forest and agro-ecosystems (Bol et al., 2016(Bol et al., , 2018;;Missong et al., 2018). In general, phosphates bind to soil components like soil organic matter (Ahmed et al., 2018a;Gros et al., 2017Gros et al., , 2019) ) and to soil minerals like Fe/Al(oxyhydr)oxides (Hens and Merckx, 2001;Kruse et al., 2015;Jiang et al., 2015;Ahmed et al., 2019) and amorphous Fe/Al hydroxide mixtures (Gypser et al., 2018). The P bound to soil minerals forms colloidal P complexes and consequently becomes unavailable to plants causing P inefficiency (Bol et al., 2016;Holzmann et al., 2015). These colloidal P complexes disperse during heavy rains and accumulate in specific regions resulting in P leaching which further reduces P availability to plants (Boy et al., 2008). Molecular level understanding of P adsorption onto these soil minerals could support efforts to improve P availability to plants (Bol et al., 2016).</p><p>Goethite (α-FeOOH) is one of the most common and abundant soil minerals that interacts strongly with phosphates (Parfitt and Atkinson, 1976;Torrent et al., 1992;Chitrakar et al., 2006;Kubicki et al., 2012;Ahmed et al., 2019). It is a highly reactive soil mineral containing ferric ions (Fe +3 ) with common surface planes as 010, 100, 110, 021 (according to pbnm space group) (Cornell and Schwertmann, 2003). The surface iron atoms are coordinated by 3, 4, 5 or more atoms depending on the surface plane as well as the pH of the environment. Consequently, goethite exhibits different levels of saturation according to the interaction with its environment. The same holds true for most minerals, i.e. minerals exhibit a net positive or negative surface charge based on surface (un)saturation and pH (Cornell and Schwertmann, 2003). Hence, the type of surface plane, its termination and saturation are important factors that influence the adsorption of phosphates onto soil minerals. For instance, Ahmed et al. (2018b) studied glyphosate adsorption at goethite surface with three different degrees of (un)saturation (Fe surface atoms coordinated by 3, 4, and 5 O −2 /OH − groups) and showed the effect of the surface's (un)saturation on phosphate binding stability.</p><p>In addition to surface saturation and pH, the Fe and Al ratio in amorphous Fe/Al hydroxides is also vital for understanding the phosphates' interaction with soil minerals. Gypser et al. (2018) showed the influence of the Fe:Al ratio in amorphous Fe/Al hydroxide mixtures on phosphate adsorption/desorption rates. The omnipresent Al in weathering environment results in most of Fe oxides in soils being substituted by Al and goethite is no exception (Cornell and Schwertmann, 2003). Diaspore (α-AlOOH) is isomorphous with goethite with Al +3 oxidation state exhibiting a higher surface energy compared to goethite (Haibo and S., 2011). Since amorphous Fe/Al hydroxide mixtures exist in soils, analyzing phosphate binding to diaspore provides additional insight into the P interaction with these amorphous mixtures.</p><p>Orthophosphates (Newman and Tate, 1980), inositolhexaphosphate (IHP) (Turner et al., 2002;Doolette et al., 2009;Gerke, 2015) and glycerolphosphate (GP) (Pant et al., 1999;Vincent et al., 2013;Missong et al., 2016) are the most abundant phosphates in soils. Orthophosphate interaction with goethite has been studied extensively (Parfitt and Atkinson, 1976;Torrent et al., 1992;Chitrakar et al., 2006;Ahmed et al., 2019). Yan et al. (2014) studied sorption of different phosphates involving GP and IHP on aluminium (oxyhydr)oxides. They found that the maximum adsorbed phosphate normalized to the mass of adsorbent increases with decreasing crystallinity of the minerals: α − Al 2 O 3 < boehmite < amorphous Al(OH) 3 . Moreover, they concluded that the phosphate adsorption, interfacial reactions, and phosphate fate in the environment are strongly affected by molecular structure and size of phosphates, and degree of crystallinity and crystal structure of mineral surfaces. Li et al. (2017) suggested that GP adsorbs onto the goethite surface through its phosphate group forming inner-sphere complexes. IHP has six phosphate groups, and in general it exhibits strong binding with P-fixing minerals compared to other phosphates with fewer phosphate groups. Anderson and Arlidge (1962) suggested that the total number of phosphate groups in a compound determines the stability of its interaction with minerals. Ognalaga et al. (1994) showed that IHP forms inner-sphere complexes with goethite through its phosphate groups and suggested that up to four phosphate groups could be involved in binding with mineral surface; the remaining non-interacting phosphate groups could alter the electrochemical properties of the surroundings. However, Guan et al. (2006) revealed that only three phosphate groups were bound to aluminium hydroxide while the other three groups remained free. This was based on adsorption experiments of IHP on amorphous aluminium hydroxide, FTIR characterization, and quantum chemical calculations.</p><p>Quantum chemical calculations become increasingly important when it comes to develop a mechanistic understanding of chemical processes in general.</p><p>Although not yet widely used, computational chemistry approaches to environmentally relevant questions are recognized as tools to complement experimental investigations. Interestingly, as early as 1973Tossell et al. (1973) ) studied electronic structure and bonding in iron oxide minerals with molecular simulations and validated this approach with experimental studies. Kwon and Kubicki (2004) used molecular simulations to resolve controversies in experimental studies related to phosphate surface complexes on iron hydroxides. In another study, Kubicki et al. (2012) demonstrated that phosphate interaction with goethite involves a variety of surface complexes in multiple configurations, which explained the difficulties one faces when interpreting, e.g., IR spectra. Moreover, Ahmed et al. (2018b) explored the possible binding mechanisms for glyphosate (GLP) with three goethite surface planes (010, 001, and 100) in the presence of water via ab initio molecular dynamics simulations. The results showed the prominence of water in controlling the GLP-goethite-water interactions. Further (Ahmed et al., 2019) investigated the molecular level mechanism of phosphate binding at the goethite-water interface referring to the possible phosphate binding motifs formed at the modeled goethite surface planes. Moreover, the theoretical assignment of IR spectra in the latter study introduced a benchmark for characterizing experimental IR data for a distribution of adsorbed phosphate species. This incomplete list already indicates the huge potential of computational chemistry as an emerging powerful tool for detailed investigations of complex geochemical reactions and especially reaction mechanisms of P-species in soil (for a more complete overview, see also (Kubicki, 2016)). Which computational methods are available today? Thinking of it in terms of a hierarchy (Ozboyaci et al., 2016;Kubicki, 2016), quantum mechanics (QM) methods are at the top because they provide, at least in principle, an unbiased ab initio description in terms of the molecular Schrödinger equation. In practice, such a treatment is not feasible except for the simplest cases. Density Functional Theory (DFT) has emerged as a low-cost alternative, which despite of the involved approximations often provides reliable results even for systems with hundreds of atoms (Kubicki, 2016). At the bottom of the hierarchy there are molecular mechanics (MM) methods, which are based on parametrized empirical functions (force fields) describing bonded and non-bonded interactions. This allows to treat millions of atoms, but unless special purpose force fields are used chemical reactions, i.e. bond making and breaking, cannot be simulated (González, M.A., 2011). Here, the hybrid QM/MM approach comes into play, combining the accuracy of QM methods with the efficiency of MM methods. In this approach, the reactive region of the molecular system is treated at the QM level while the remaining part of the system enters at the MM level (Senn and Thiel, 2009). By this, the bond changes, proton transfer events and hydrogen bonds (HBs), e.g., at water-mineral interfaces can be properly described. Besides providing energies in dependence on nuclear positions, forces on the nuclei can be calculated as well. This is prerequisite for molecular dynamics simulations which gives statistical information including thermally accessible configurations (Marx and Hutter, 2009).</p><p>In an earlier hybrid QM/MM study of IHP and GP binding to the 010 diaspore surface plane we have demonstrated a strong interaction of IHP/GP with the diaspore surface (Ganta et al., 2019). Here, IHP forms multiple intramolecular HBs with three of its phosphate groups bound to the surface, while GP is bound through its single phosphate group only. Overall, it has been found that proton transfers from phosphate to water or surface have a stabilizing effect, most likely due to the interaction of the HBs dipole with surface charges. Moreover, in case of IHP intramolecular HBs can be formed, which lead to a steric constraint that could weaken the binding to the surface. Since the interaction of IHP and GP with diaspore is not yet fully explored, in the present work we extend our previous study in two directions, i.e. we consider a chemically different surface plane and incorporate the effect of saturation of the diaspore surface on phosphates adsorption. In our previous work (Ganta et al., 2019), IHP/GP and water showed strong and spontaneous interactions with an unsaturated diaspore surface (010 in pnma) where the surface Al atoms are coordinated by four oxygens (O −2 /OH − groups). Here, a more saturated diaspore surface (100 in pnma) is selected where the surface Al atoms are coordinated by five oxygens i.e. O −2 /OH − groups. The main objective of current work is to characterize the binding mechanism of IHP and GP at this diaspore-water interface and also to understand the effect of (un)saturation of the diaspore surface on this binding mechanism.</p><!><p>In general, the surface charge of a certain mineral can be determined as a function of pH via its point of zero charge (PZC) (Tan, 2011). For pH > PZC, the mineral surface is saturated with negative surface charges which attract cations. In contrast, for pH < PZC, the mineral surface is unsaturated with positive surface charges which attract anions. The phosphates in general exhibit overall negative charge and hence can be adsorbed at the partially unsaturated 100 diaspore surface (according to the Pnma space group) (Tan, 2011). The diaspore unit cell has four AlO(OH) units i.e. total of 16 atoms with lattice constants a = 9.4253, b = 2.8452, c = 4.4007 Å, see Fig. 1c. The 100 surface plane model is generated by repetition of the diaspore unit cell as 1a × 8b × 5c along x, y, z axes respectively, see Fig. 1d. In total, the used diaspore slab consists of 640 atoms (160 Al,160 H,and 320 O atoms). Observe that the surface Al atoms are coordinated by five oxygen atoms, see Supplementary Material (Suppl. Mat.) Fig. S3b. The phosphates IHP and GP (see Fig. 1a-b) are modeled to have their phosphate group(s) interacting via inner-sphere complexes with surface Al atoms of diaspore, see Fig. 1f-h and Suppl. Mat. Fig. S2a-c. The diaspore-IHP/GP complexes are then solvated using the solvate plugin from the VMD package (Humphrey et al., 1996) with a water layer of about 18 Å perpendicular to the surface along the x axis and with a density of ≈ 1g cm −3 , see Fig. 1e. Since we are interested in IHP/GP interaction at diaspore-water interface, the QM part of the system (see dashed box in Fig. 1e) includes the top layer of diaspore (160 atoms), IHP (54 atoms)/GP (19 atoms) and a few water molecules (≈ 53 molecules depending on the setup) surrounding IHP/GP within layer of ≈ 10 Å perpendicular to the diaspore surface. The enclosing QM box has a size of 22 × 8b × 5c Å i.e. The initial motifs of diaspore-GP complexes include the monodentate motif M (1Al+1O), see Fig. 1f and bidentate motif B (2Al+2O, here both oxygens are from same phosphate group), see Fig. 1g. Note that in contrast to this setup, in the 010 case the two oxygens bind to the same Al atom. In addition to these two, the diaspore-IHP complexes include the four monodentate motif 4M (4Al+4O) as experimental studies suggest that IHP forms multiple bonds with the goethite surface (Ognalaga et al., 1994;Guan et al., 2006), see Fig. 1h and Suppl. Mat. Fig. S2c. Note that in principle even more initial conditions/motifs could be sampled. But considering the size of the modeled systems here and the used computationally expensive QM/MM level of theory the initial configurations for the MD simulations are limited to the most common and experimentally observed binding motifs. More technical details about QM and MM methods and their mutual interaction adopted here is given in Suppl. Mat.. The QM/MM based MD simulations are performed for 25 ps with 0.5 fs time step and with an average temperature of 300 K maintained using canonical sampling through the velocity rescaling thermostat (CSVR) (Bussi et al., 2007).</p><p>Here, for each molecular model the first 10 ps of the MD trajectory is assigned for equilibration. The last 15 ps of the trajectory is considered as the production trajectory which is used for analysis of interactions at diaspore-water interface with IHP/GP.</p><p>To analyze the interaction energies of the complexes, snapshots are taken at every 100 fs of the production trajectory and interaction energies between diaspore and IHP/GP (E diaspore−IHP/GP ), IHP/GP and water (E IHP/GP−water ), and diaspore and water (E diaspore−water ) are calculated. For example, the interaction energy E int between diaspore and GP for a certain diaspore-GP-water snapshot is calculated as follows:</p><p>here E diaspore−GP , E GP , and E diaspore denote electronic energies of the diaspore-GP complex, GP and diaspore surface, respectively. Likewise, the interaction energies for each pair of diaspore, IHP/GP and water are calculated at every 100 fs during the corresponding production trajectory. The interaction energies with water are divided by the total number of water molecules involved in the simulation box for better comparison. More details regarding the calculation of interaction energies are given in (Ganta et al., 2019).</p><p>The HBs strength between IHP/GP and water as well as for the intramolecular HBs of one IHP motif are analyzed using geometrical correlations of distances between atoms in HB as discussed in (Strassner, 2006;Limbach et al., 2009;Yan and Kühn, 2010;Zentel and Kühn, 2017). The quantities q 1 and q 2 in Figs. 4 and 6 below are defined as the deviation of the hydrogen from HB center assuming a linear HB (q 1 ) and the total HB length (q 2 ), see Fig. 2a. Geometrically q 1 and q 2 are denoted as q 1 = 1 2 (r 1 − r 2 ) Å and q 2 = (r 1 + r 2 ) Å where r 1 , r 2 denote the distance between donor oxygen and hydrogen (O D -H) and distance between hydrogen and acceptor oxygen (H• • • O A ), respectively, see Fig. 2a. The same holds true for intramolecular HBs between phosphate groups for the IHP case, see Fig. 2b.</p><p>Figure 2. Definition of coordinates q 1 and q 2 in a HB between water's oxygen (O A ) and IHP's oxygen (O D ) from a phosphate group (a) (Limbach et al., 2009), intramolecular HB in IHP (b).</p><p>A HB will be called strong if q 1 ≈ 0 and q 2 is in the range 2.2-2.5 Å. Similarly, moderate and weak HBs have q 2 distances ranging from 2.5-3.2 Å and 3.2-4 Å respectively. Also if q 1 < 0 then hydrogen stays with the donor oxygen and if q 1 > 0 then hydrogen transferred to the acceptor oxygen, i.e. a proton transfer event is observed. In the following a HB analysis is performed for the QM part of the system only as the emphasis of this study is on the interface region where the binding/adsorption process takes place.</p><!><p>3.1 Diaspore-GP-water interactions</p><!><p>For the GP-M initial condition a stable monodentate (1Al+1O) motif is observed between GP and the diaspore surface over the course of the production trajectory with the Al-O1 average bond length of 2.3 Å, see Fig. 3a-d. The average geometry of the PO 4 moiety here has root mean square deviation (RMSD) value of 0.17 Å with respect to free tetrahedral PO 3− 4 , see Suppl. Mat. Fig. S2d. Proton transfer events from GP to the diaspore surface are not observed, instead three proton transfer events are found from O3, O5 and O2 oxygens of GP to water, see Fig. 3a-c respectively. On the average, eight HBs are observed between GP and water in the production trajectory. Here, the GPs oxygen atoms act as HB donors (Ox D ) as well as acceptors (Ox A ). Exemplary analysis of six of the above eight HBs shows that four (O2 A , O3 A , O5 A , O5 D ) are strong to moderately strong HBs and two (O1 A , O6 A ) moderately strong to weak HBs, see Fig. 4a-c.</p><p>Regarding the diaspore-water interaction, an average of 17 water molecules (out of 40 surface Al atoms) formed M binding motifs (Al-O H 2 O ) with the diaspore surface and have average bond length of 1.9-2.3 Å. Also moderately strong HBs are observed between water and diaspore, see Suppl. Mat. Fig. S1a. This scheme of diaspore-water interactions is also observed for the other diaspore-IHP/GP-water models studied below. For the average diaspore-water interaction energy per water molecule one obtains about -3 kcal/mol for all considered models.</p><p>The time averaged interaction energy per surface bond between diaspore and GP is around -23 kcal/mol, see Table 1. The average GP-water interaction energy per water molecule is -2.6 kcal/mol.</p><p>. HB correlation q 2 vs. q 1 of HBs formed between GP and water in GP's M motif (a)-(c) and B motif (d)-(f), respectively. The average q 1 and q 2 coordinate pairs are shown as points in white square boxes. A strong to moderately strong HB is denoted as * while a moderately strong HB and a moderately strong to weak HB are denoted as • and , respectively.</p><p>Table 1. The per bond time averaged interaction energies (calculated with Eq. 1) of diaspore-IHP/GP complexes and some selected bond lengths and distances. The RMSD difference between average phosphate geometry and isolated tetrahedral PO 3− 4 (see Suppl. Mat. Fig. S2d) is also provided. Here, P o denotes phosphate and motif denotes average motif observed during production trajectory.</p><!><p>The B motif (2Al+2O i.e, Al1-O1 and Al2-O3) is observed over the course of the production trajectory with Al1-O1 and Al2-O3 covalent bond length ranging from 2-2.7 and 1.9-2.7 Å with an average value of 2.4 and 2.3 Å, respectively, see Fig. 3e-h. Most notably the Al2-O3 and Al1-O1 bonds are elongated and compressed in an alternating see-saw fashion as seen in Figs. 3f and 3g. The B motif's average geometry of the PO 4 moiety has a RMSD value of 0.17 Å with respect to the free tetrahedral PO 3− 4 .</p><p>Proton transfer events are observed from O3 and O2 oxygens to water, see Figs. 3e and 3f respectively. The B motif features on average a total of seven HBs between GP and water.</p><p>According to Fig. 4d-f, four (O1 A , O2 A , O2 D , O5 D ) strong to moderately strong HBs and two (O1 A , O6 D ) moderately strong to weak HBs are formed between GP and water.</p><p>The average interaction energy per surface bond between diaspore and GP for the B motif is around -15 kcal/mol. The per surface bond interaction energy here is smaller than the M motif due to GP's see-saw type of motion over surface. Nevertheless, the total interaction energy observed here is larger than for the M motif and hence B motif is more likely to form.</p><p>Note that due to the formation of a strong to moderately strong HB with water, the O3 oxygen in M cannot easily transform into the B motif, i.e. the barrier is too high to be sampled in the present trajectory (see Fig. 4b). The average GP-water interaction energy per water molecule is around -2.3 kcal/mol. The smaller value as compared with M could be due to the additional proton transfer event observed in that case.</p><!><p>Here, two initial configurations (M, B motifs) of the diaspore-IHP-water model resulted in two different M final motifs. In the first case, M(1), the initial configuration was an M motif, wherein O11 oxygen is aligned to form a M motif with a surface Al atom, see Fig. S2a. A stable M motif is observed through out the production trajectory with average Al-O11 bond length of 2.4 Å, see Fig. 5a-d. The series of events that are observed during the formation of the M(1) motif are: a proton transfer from O12 to water (see Fig. 5a), followed by intramolecular proton transfer from O62 to O12 and formation of O13-H-O61 intramolecular HB (see Fig. 5b). After a few femtoseconds, the O12-H-O62 HB is formed and a proton transfer is observed from O32 to water (see Fig. 5c) to reach the final M motif in Fig. 5d. The events in Figs. 5a-c occured within 2 ps of simulation trajectory. Overall, a total of three protons transfer events are observed from IHP to water from O12, O53 and O32, see Figs. 5a-c respectively. On an average IHP has formed 19 HBs with water over the course of the production trajectory. Analyzing for illustration 12 out of these 19 HBs, IHP has formed nine (O12 A , O21 A , O22 A , O23 A , O31 A , O43 A , O51 A , O53 A , O63 D ) strong to moderately strong HBs, two (O33 D , O41 A ) moderately strong HBs and one (O31 A ) moderately strong to weak HB with water, see Fig. 6a-e. Interestingly, IHP also forms multiple intramolecular HBs, for instance, O13-H-O61, O41-H-O52 and O12-H-O62, see Fig. 5a-c. Analyzing the strength of the intramolecular HBs between P1 and P6 phosphate groups, one finds two (O13-H-O61, O12-H-O62) moderately strong HBs, see Fig. 6f. The average geometry of PO 4 moiety has RMSD value of 0.16 Å with respect to the free tetrahedral PO 3− 4 . The time averaged interaction energy per surface bond between diaspore and IHP is around -33 kcal/mol, see Table 1 and between IHP and water is around -5.5 kcal/mol per water molecule. Notice that IHP exhibits a larger interaction energy with water as well with the diaspore surface compared to GP.</p><p>For the M(2) motif, the initial configuration had the O11 and O13 oxygens aligned such as to form a B motif with adjacent surface Al atoms Å, see Suppl. Mat. Fig. S2b. The Al2-O13 covalent bond is dissociated during the trajectory and the B motif is transformed into the M(2) motif, see Fig. 5e-f. Over the course of the production trajectory the Al1-O11 bond length ranges from 2.4-3.2 Å with an average of 2.7 Å. In more detail, the series of events that unfold in this case are as follows: From the O23 and O33 oxygens, two protons are transferred to water and the</p><p>. Diaspore-IHP-water snapshots displaying proton transfer events and IHP dynamics along trajectories. The phosphate groups are color-coded and the oxygen atom label color denotes to the phosphate group it belongs. In M(1) motif, proton transfer from O12 to water (a), proton transfer from O53 to water and intramolecular proton transfer from O62 to O12 (b), moderately strong O12-H-O62 HB and proton sharing between O12 to O62 oxygens and proton transfer from O32 to water (d), M(1) motif at 25 ps. In M(2) motif, dissociation of Al2-O13 covalent bond due to O13-H-O61 HB and proton transfers from O23 and O33 to water (e), proton transfer from O53 to water and formation of O53-H-O63 HB (f), dissociation of Al1-O11 bond due to formation of O11-H-O21 HB (g), M(2) motif at 25ps (h). In 2M motif, dissociation of Al4-O51 bond and proton transfer from O13 to water (i), proton transfer from O43 to diaspore and from O32 to water, also proton sharing and HB between O22 and O32 (j), and between O33 and O42 followed by dissociation of Al3-O41 bond (k), proton transfer from O42 to water and 2M motif at 25 ps (i).</p><p>Al2-O13 covalent bond is dissociated, see Fig. 5e. After a few femtoseconds, a proton transfer is observed from O53 to water followed by an intramolecular HB between O53 and O63 oxygens, see Fig. 5f. Also formation of the intramolecular O13-H-O61 HB is observed. With progressing simulation time an intramolecular proton transfer event is observed from O12 to O21, followed by formation of O11-H-O21 HB. Afterwards, the Al1-O11 covalent bond weakens at around 6 ps and its bond length ranges from 2.4-3.2 Å further on, see Figs. 5g and S1b. In addition, a proton transfer is observed from O62 to water, see Fig. 5g. The snapshot at 25 ps shows that multiple inter-and intramolecular HBs are observed for IHP, see Fig. 5h. Their characterization in term of HB geometries leads to a similar distribution of HB strengths as for M(1). Overall the average geometry of PO 4 moiety has RMSD value of 0.19 Å with respect to the free tetrahedral PO 3− 4 .</p><p>2.2 2.5</p><p>. HB correlation q 2 vs. q 1 of HBs formed between IHP and water (a)-(e) and between P1 and P6 phosphate groups (f) along trajectory of IHP's M(1) motif. The average q 1 and q 2 coordinate pairs are shown as points in white square boxes. A strong to moderately strong HB is denoted as * while a moderately strong HB and a moderately strong to weak HB are denoted as • and , respectively.</p><p>The average interaction energy per surface bond between the diaspore surface and IHP in this case is around -18 kcal/mol. The interaction energy observed here is smaller than for the M(1) case as the Al1-O11 bond length is longer due to formation of O13-H-O61 and O12-H-O21 intramolecular HBs. The observed interaction energy between IHP and water here is -5.7 kcal/mol per water molecule which is slightly higher than for the M(1) case, probably due to additional proton transfer from IHP to water.</p><!><p>Here, IHP is initially aligned parallel to surface with the non-protonated oxygens of the four phosphate groups forming a 4M motif, i.e. Al1-O11, Al2-O21, Al3-O41, Al4-O51 covalent bonds with the surface, see Suppl. Mat. Fig. S2c. However, only a stable 2M motif is observed along the production trajectory with average bond lengths of Al1-O11, Al2-O21 as 2.13 and 2.22 Å, respectively. The events observed during simulation that led to the formation of 2M motif are as follows: within a few picoseconds, the Al-O51 bond is dissociated transforming 4M into a 3M motif. A proton transfer is observed from O13 oxygen to water and from O43 oxygen to diaspore, see Figs. 5i and 5j. Further, an intramolecular HB is observed between O22 and O32 oxygens and two proton transfer events from O32 and O52 to water, see Fig. 5j. The Al3-O41 covalent bond is disassociated due to intramolecular HB formed between O33 and O42 oxygens, see Fig. 5k, followed by a proton transfer event from O42 to water, see Fig. 5l. Totally four proton transfer events are observed from IHP to water and an average of 19 HBs are formed between IHP and water. The inter-and intramolecular HBs have a similar distribution of strengths as for M(1) and M(2). The average geometry of the PO 4 moieties deviate from that of free tetrahedral PO 3− 4 with RMSD values of 0.18 and 0.17 Å.</p><p>The interaction energy between IHP and diaspore in the 2M motif is -109 kcal/mol per bond, see Table 1. The interaction energy is larger here compared to M(1) and M(2) motifs due to additional covalent bond and a proton transfer from IHP to the diaspore surface. Hence, the 2M motif is more likely to form compared to both the M(1) and M(2) motifs. The average interaction energy per water molecule with IHP is -5.9 kcal/mol which is slightly larger than the IHP-water interaction energy observed in M(1) motif case due to additional proton transfer from IHP to water here.</p><!><p>In the following we will compare the present results with those of our previous work for the 010 surface plane (Ganta et al., 2019). The 010 diaspore surface plane is relatively unsaturated, i.e. the surface Al atoms are coordinated by only four oxygen atoms, see Suppl. Mat. Fig. S3a. In contrast, for the present more saturated 100 plane, surface Al atoms are coordinated by five oxygens, see Suppl. Mat. Fig. S3b. Also the 010 diaspore surface plane exhibits higher electrostatic potential compared to the 100 surface plane as shown in Suppl. Mat. Fig. S6.</p><p>For the 010 plane the largest total interaction energy was observed for the B and 3M motif in case of GP and IHP, respectively, see Table . 1. In case of 100 plane, the B and 2M motifs dominate the total interaction energies. Comparing the two B motifs for GP one finds that the binding to 010 being 10 times stronger than to 100 surface plane. The reason for the weaker interaction energy in case of 100 is due to see-saw type of motion of GP yielding a weakening/strengthening of Al1-O1 and Al2-O3 bonds which is not observed for the 010 plane. This can be attributed to the fact that in case of 010 plane the two oxygens are coordinated to the same Al, whereas for 100 plane the coordination is with two neighbouring Al atoms, whose distances is such as to require unfavourably large O P -P-O P angles for strong binding. Further stabilization of the GPs B motif in the 010 plane case comes from two additional proton transfers observed from GP to the diaspore surface.</p><p>Regarding the total interaction energies, the dominant binding motifs for the diaspore-IHP are also different for 010 (3M) and 100 (2M). In case of 100 the total interaction energy is about two times smaller than for 010. Comparing the two motifs we note in particular that the interaction with the two surfaces is different. This is nicely illustrated by the fact that no stable 4M motif could be observed for the 100 case. One reason for the transformation of the 4M motif to 2M motif is that the Al-O P bonds (regions R1 and R2, see Suppl. Mat. Fig. S5a-c) at the 100 diaspore surface are inclined due to O P and surface hydroxyl oxygen repulsion. Hence the movement of oxygens in the Al-O P bonds (regions R1 and R2) is restricted to the space between consecutive surface hydroxyl oxygens or to move away from surface, see Suppl. Mat. Fig. S5a-c. Consequently, upon equilibration the oxygens in the Al-O P bonds could dissociate from diaspore as they are confined between consecutive surface hydroxyl oxygens, see Suppl. Mat. Fig. S5a-c. In contrast for the 010 diaspore-IHP case, the oxygens in the Al-O P bonds (regions R3 and R4, see Suppl. Mat. Fig. S5e-g) are not restricted and they are free to move. Hence a stable 3M motif is observed over the course of production trajectory, see Fig. S5e-g. Looking at it from a geometric point of view, Al-Al distances on the 010 surface are about 4.4 to 5.4 Å which is much larger than the 2.4 Å for the 100 surface where the 2M motif forms, see Suppl. Mat. Figs. S5d and h. Given the typical distances between the phosphate groups in IHP, bonding to 100 surface plane yields a higher strain and thus it becomes weaker as compared with 010.</p><p>In case of the 100 surface the diaspore-water interaction energy is 3.4 times smaller compared to the 010 case. In fact less than half of surface Al atoms formed M motifs with water compared to the 010 diaspore surface. Also the radial distribution function of diaspore surface oxygens with water hydrogens in Suppl. Mat. Fig. S4 shows higher water accumulation near the 010 diaspore surface compared to the 100 diaspore surface which suggests stronger interaction for the 010 diaspore surface with water.</p><!><p>In our previous study (Ganta et al., 2019), a strong and spontaneous binding of IHP and GP with the 010 diaspore surface has been described, which provided the motivation for studying the effect of surface saturation on these interactions. Therefore, the more saturated 100 diaspore surface has been investigated here using periodic boundary QM/MM based MD simulations. The analysis of the MD trajectories showed the importance of inter-and intramolecular HBs in the formation of final motifs and also shed light on effects that lead to disassociation and association of P-O-Al bonds in the diaspore-IHP/GP-water complexes.</p><p>In case of the diaspore-GP-water complexes, the B motif's interaction energy per bond is 1.5 times smaller than the M motif. But considering the total average interaction energy, GP is more likely to form a B motif with the 100 diaspore surface.</p><p>Regarding the diaspore-IHP-water complexes, the interaction energy per bond follows the order 2M > M(1) > M(2). Here, the M(2) motif's interaction energy is 1.8 times smaller than the M(1) motif due to longer Al-O p bond length, i.e, due to movement of IHP away from the diaspore surface. Thus the 2M motif will be also dominating considering the total interaction energy. This is due to the additional covalent bond and a proton transfer to the diaspore surface. Hence IHP is likely to form a 2M motif with the 100 diaspore surface.</p><p>Regarding the water interaction with 100 diaspore and IHP/GP, it can be concluded that the average IHP-water interaction energy is about 2.3 times larger than the GP-water one due to IHP's higher water accessible surface area. Both IHP and GP show proton transfer events to water and formation of strong to moderately strong HBs with water. The diaspore-water interaction energy is only 1.1 times that of the GP-water case, but 2.8 times smaller than IHP-water one. Thus, water has a stronger interaction with IHP than to the 100 diaspore surface.</p><p>Of course, studying a particular perfect surface plane can at best give qualitative trends if compared to real surfaces of mineral particles in soil. Studying two abundant surface planes, however, it is possible to pinpoint important factors which influence the behaviour of P-compounds at the mineral/water interface. The present investigation focused on the effects of surface saturation and thus electrostatic potential on IHP/GP adsorption. The soil minerals exhibit a positive charge with an unsaturated surface and active sites for pH < PZC. For a pH far below PZC a higher positive charge, i.e, more unsaturated surface is observed (Tan, 2011;Cornell and Schwertmann, 2003). The IHP and GP adsorption onto goethite (with PZC around 9-10) is decreased with increasing pH and reached near zero values for pH around 10 (Celi et al., 1999;Li et al., 2017). This shows that the mineral surface saturation varies with pH and thus the ability of a mineral to adsorb phosphates. Higher surface saturation leads to more negative charges on the 100 surface as compared to the 010 case. Phosphates in water are partially deprotonated and thus have an effective negative charge. Thus the phosphate groups will be attracted stronger to the 010 diaspore</p>