compound that due to its hygroscopic nature has to be discussed in context of water content, and which can polymerize at drying conditions.

Chapter 3 gives a comprehensive insight into how more than 1000 geometries were explored by density functional theory (DFT) calculations. Clearly, the real surface species cannot be represented by a single reference adsorption state. However, it will be shown in this chapter, how a relatively small number of structures can be chosen unambiguously.

The resulting geometries were selected by their free enthalpies. The second challenge is to verify those model structures with experimental results.

In chapter 4, chemical shifts were calculated for these selected geometries and compared with 31 P-MAS NMR experiments, that were done in collaboration with researchers at Centre de RMN à Très Hauts Champs (Université de Lyon) and IFP Energies nouvelles. This work is published in ref. [18]. It is also illustrated here, how well the most stable systems agree with experimental NMR spectra. Overall, excluding all but the most favorable systems is a straightforward method of information reduction and allows to extract conclusions by direct comparison.

It was also used later in chapter 5 in order to choose the adsorption mode for biased molecular dynamics calculations. The questions in this chapter are, how does phosphate behave under impregnation conditions? And which adsorption states are accessible by kinetics? Biased molecular dynamics (metadynamics) were used in order to test the adsorption and desorption mechanism as well as several more complex phosphorylation reactions between phosphates, water and the most important surface sites. It is not feasible to evaluate all geometries produced with the initial exploration, thus in chapters 3 to 5, trends are only identified within a hand full of structures. However, different evaluation methods permit the use of hundreds or all of the test structures and extract data in histograms and averaged structure parameters (bond lengths, angles, etc). Such methods, today often referred to as data mining and machine learning, allow to search for

List of Tables

Corral Valero for all the support and advice. Without your dedicated involvement this work would have never been possible. Thank you, Carine, for your tireless encouragement and your engaged focus on everyone's strengths. I am grateful to encounter your mindset. Thank you, Pascal, for sharing your experience and incredibly detailed knowledge with me. You knew every detail to enrich a discussion and I always enjoyed your enthusiasm. Thank you, Manuel, for calm guidance and your great patience with me on a personal level. You helped me a lot to communicate technical details in an understandable manner and allowed me to substantially improve my discussions. I honestly humbles me having had the chance to work with you. Having three supervisors can certainly be complicated, so I consider myself very lucky to have caught this great team.

To the jury of my defence: I want to thank Anne Milet and Dominique Costa, who agreed to evaluate this thesis, as well as Paul Fleurat-Lessart and Christel Gervais for taking their time to participate as jury members. As expected for an independent jury, we do not know each other allowing you to examine my work without any conflict of interest. Because of this, I appreciate the time you are going to spend for this and hope I can sustain the initial interest you had in my work.

I also am greatly indebted to Dorothea Wisser, who carried out the complex NMR experiments which were crucial in this work, and to Tao Jiang who introduced me to machine learning methods and explored various different options for NMR predictions training on the basis of DFT calculations. For all the productive discussions in the past three years I am grateful to Jérôme Rey, Daniel Walden, Thomas Pigeon, Teddy Roy, Stephan Steinmann, David Loffreda and certainly many other people.

Not least I want to thank my friends and colleagues I made here, including my reliable friends Amit Sahu and Sharmin Sharna who gave me great company and personal support over the last three years, and especially my colleagues Sarah Blanck and Laureline Treps who helped me on every request with official paperwork in the French language. I must also thank Malika, my wife and companion, who supported and encouraged me and who carefully read this work for additional corrections.

Finally, I also have to thank the ROAD4CAT industrial chair, part of the project IDEXLYON, as well as ENS Lyon for funding my position and this work.

Abstract

Phosphate adsorption on oxides is of significant interest in heterogeneous catalysis. The CoMoS/γ-Al 2 O 3 hydrodesulfuration catalyst which is used at the scale of petroleum refinement, uses phosphate additives for preparation. Its preparation follows the stages of wet incipient impregnation, drying and activation by thermal treatment. Understanding the chemistry involved in these processes requires models of the oxide-water and oxide-air interfaces at different temperatures, which are also of interest in different fields. Oxide surfaces like γ-Al 2 O 3 are by themselves very complex, which renders the characterization of surface phosphate species even more intricate.

Yet by today, improved NMR spectroscopic methods complemented with proven computational models allowed us now to mold an atomistic model of phosphate speciation. The present work starts with the creation of a systematic database by exploring more than 1000 geometries at different adsorption modes, adsorption sites, degrees of coverage and an optimized hydrogen bond environment for two polyphosphate species. For the most stable results magnetic shielding was calculated at DFT level. Good agreement was achieved with [START_REF] Luis | Hydrolysis rates of inorganic polyphosphates in aqueous solution as well as in soils and effects on P availability[END_REF] P NMR experiments at drying conditions. Liquid conditions were simulated using metadynamics and showed a different set of most stable species. The explored reaction energies and kinetic barriers consistently show that elevated temperatures are needed to form the experimentally observed species. A 31 P NMR prediction model based on machine learning of chemical shifts from the local structure was developed to access chemical shifts in dynamic simulations.

Resumé

L'adsorption des phosphates sur les oxydes présente un intérêt significatif pour les catalyseurs hétérogènes. Le catalyseur d'hydrodésulfuration CoMoS/γ-Al 2 O 3 qui est utilisé à l'échelle du raffinage du pétrole, utilise des additifs phosphatés pour sa préparation. Sa préparation suit les étapes d'imprégnation humide, de séchage et d'activation par traitement thermique. La compréhension de la chimie impliquée dans ces processus nécessite des modèles des interfaces oxyde-eau et oxyde-air à différentes températures, qui présentent également un intérêt dans différents domaines. Les surfaces d'oxyde comme la γ-Al 2 O 3 sont en elles-mêmes très complexes, ce qui rend la caractérisation des espèces de phosphate de surface encore plus compliquée.

Pourtant, aujourd'hui, les méthodes améliorées de spectroscopie RMN complétées par des modèles de calcul éprouvés nous ont permis d'élaborer un modèle atomique de spéciation des phosphates. Le présent travail commence par la création d'une base de données systématique en explorant plus de 1000 géométries à différents modes d'adsorption, sites d'adsorption, degrés de couverture et un environnement de liaison hydrogène optimisé pour deux espèces de polyphosphates. Pour obtenir les résultats les plus stables, le blindage magnétique a été calculé au niveau DFT. Un bon accord a été obtenu avec les expériences de [START_REF] Luis | Hydrolysis rates of inorganic polyphosphates in aqueous solution as well as in soils and effects on P availability[END_REF] P RMN dans des conditions de séchage. Les conditions liquides ont été simulées en utilisant la métadynamique et ont montré un ensemble différent d'espèces les plus stables. Les énergies de réaction et les barrières cinétiques explorées montrent systématiquement que des températures élevées sont nécessaires pour former les espèces observées expérimentalement. Un modèle de prédiction [START_REF] Luis | Hydrolysis rates of inorganic polyphosphates in aqueous solution as well as in soils and effects on P availability[END_REF] P RMN basé sur l'apprentissage automatique des déplacements chimiques à partir de la structure locale a été développé pour accéder aux déplacements chimiques dans les simulations dynamiques.

Introduction

Catalysts consist of active sites at which the reactions take place that would not occur otherwise due to considerable energy barriers. For a catalyst to be economically viable, it has to be ensured that the amount of active sites is maximized. In heterogeneous catalysts, this is typically done by dispersing the active phase on a support with a high specific surface area. The preparation of these catalysts is delicate and demands well-controlled processing steps. It starts from a suitable support material which is impregnated with a solution containing, for example, metal precursors and additives. The impregnated support is then dried to remove the solvent and activated by thermal treatment, whereby catalytically active sites are formed.

The advantage of heterogeneous catalysts is that the support can be a cheap mineral or oxide like alumina or silicates which carries small amounts of a highly dispersed active phase. However, it also means that several components are involved in the course of a reaction -at least the support, active phase and reactant. On top of this, several additives are commonly used during preparation. At the mesoscale, material transport within the support influences the preparation as well as the performance of the final catalyst. Much of the knowledge about the preparation of heterogeneous catalysts is based on empirical knowledge and it is difficult to rationalize the effect of each component without accurate models. For example, several organic and inorganic compounds, e. g. phosphoric acid, are known to improve the final catalyst, but it is yet unclear what the effect of those additives is during preparation and how they can impact the catalytic Introduction performance. Moreover, organic additives are typically removed during the drying or activation steps of preparation, while phosphates would remain as part of the catalyst. This raises more questions about the actual role of phosphates in every preparation step and the final catalyst.

An important example for such a catalyst is the γ-Al 2 O 3 -supported cobalt molybdenum sulfide catalyst (CoMoS/γ-Al 2 O 3 ) that is used in petroleum refinement for hydrodesulfurization (HDS) and hydrodenitrogenation (HDN) treatments. [START_REF] Jaap | Envisaging the Physicochemical Processes during the Preparation of Supported Catalysts: Raman Microscopy on the Impregnation of Mo onto Al 2 O 3 Extrudates[END_REF][START_REF] Toulhoat | Catalysis by Transition Metal Sulphides: From Molecular Theory to Industrial Application[END_REF][START_REF] Van Cleve | Enhanced Hydrothermal Stability of γ-Al 2 O 3 Catalyst Supports with Alkyl Phosphonate Coatings[END_REF] At this scale, it is critical to optimize the efficiency of the catalyst as well as the use of materials which also includes ingredients for the catalyst. More aspects such as aging determine the overall life span of the catalyst. An optimum phosphate concentration at the oxide support must be targeted to avoid the formation of undesirable and inactive phases during activation. [START_REF] Toulhoat | Catalysis by Transition Metal Sulphides: From Molecular Theory to Industrial Application[END_REF] Precise control over the phosphorus content and its dispersion on the catalytic support will not only impact the catalyst performance but can also minimize the consumption of phosphates.

Van Cleve et al. have demonstrated that phosphonate deposition increases the hydrothermal stability of γ-Al 2 O 3 , hindering the phase transition to boehmite and thus allowing the use of this support in aqueous phase reforming processes. [START_REF] Van Cleve | Enhanced Hydrothermal Stability of γ-Al 2 O 3 Catalyst Supports with Alkyl Phosphonate Coatings[END_REF] Phosphates may be added to the impregnation solution in order to modify the interplay between various metal ions in the solution and their behavior inside the alumina pores. [START_REF] Jaap | Envisaging the Physicochemical Processes during the Preparation of Supported Catalysts: Raman Microscopy on the Impregnation of Mo onto Al 2 O 3 Extrudates[END_REF] Here, the authors argue that phosphorus influences the interaction of the metallic precursors with the support and impacts the dispersion, size and morphology of the active phase. [START_REF] Iwamoto | Influence of Phosphorus on the Properties of Alumina-Based Hydrotreating Catalysts[END_REF] More recently, van Haandel et al. reported that the activity of hydrodesulfurization CoMoS catalysts was improved by adding phosphate since it reduced the formation of the undesired Co 9 S 8 phase. [START_REF] Van Haandel | The effect of organic additives and phosphoric acid on sulfidation and activity of (Co)Mo/Al 2 O 3 hydrodesulfurization catalysts[END_REF] This is again showing that phosphorus impacts the synthesis of the active phase.

Van Veen et al. suggest that an undesirable aluminium phosphate phase appears at high phosphate loading, when only using phosphates and molybdates for the catalyst preparation. [START_REF] Van Veen | Chemistry of phosphomolybdate adsorption on alumina surfaces. 1. The molybdate/alumina system[END_REF][START_REF] Van Veen | Chemistry of phosphomolybdate adsorption on alumina surfaces. 2. The molybdate/phosphated alumina and phosphomolybdate/alumina systems[END_REF] However, those are added together with other additives and metal precursors during catalyst preparation, leading to a much more complex multicomponent mixture. Combined spectroscopic studies (Raman, IR, UV/vis and Raman/NMR) show that the phosphate-γ-Al 2 O 3 interaction is most critical. [START_REF] Jaap | Spatially Resolved Raman and UV-visible-NIR Spectroscopy on the Preparation of Supported Catalyst Bodies: Controlling the Formation of H 2 PMo 11 CoO 5 -40 Inside Al 2 O 3 Pellets During Impregnation[END_REF][START_REF] Catita | Application of Magnetic Resonance Imaging and Raman Imaging to study the impact of phosphorus in impregnation of hydrotreatment catalysts[END_REF] Phosphoric acid is usually co-impregnated with other salts during catalyst preparation and Bergwerff et al. show that the speciation of these co-impregnation solutions can be extremely complex. [START_REF] Jaap | Envisaging the Physicochemical Processes during the Preparation of Supported Catalysts: Raman Microscopy on the Impregnation of Mo onto Al 2 O 3 Extrudates[END_REF][START_REF] Jaap | Spatially Resolved Raman and UV-visible-NIR Spectroscopy on the Preparation of Supported Catalyst Bodies: Controlling the Formation of H 2 PMo 11 CoO 5 -40 Inside Al 2 O 3 Pellets During Impregnation[END_REF] In some cases, the chemical equilibrium can shift towards the formation of precursor species of aluminium phosphates at the surface of alumina. [START_REF] Griboval | Alumina supported HDS catalysts prepared by impregnation with new heteropolycompounds. Comparison with catalysts prepared by conventional Co-Mo-P coimpregnation[END_REF][START_REF] Griboval | Hydrotreatment Catalysts Prepared with Heteropolycompound: Characterisation of the Oxidic Precursors[END_REF][START_REF] Decanio | FT-IR and Solid-State NMR Investigation of Phosphorus Promoted Hydrotreating Catalyst Precursors[END_REF] Moreover, it was shown that phosphate interacts more strongly with the alumina surface than molybdenum species by various characterization techniques such as spatially resolved Raman, UV/vis and NIR spectroscopy as well as by Magnetic Resonance Imaging combined to Streamline Raman Imaging. [START_REF] Jaap | Spatially Resolved Raman and UV-visible-NIR Spectroscopy on the Preparation of Supported Catalyst Bodies: Controlling the Formation of H 2 PMo 11 CoO 5 -40 Inside Al 2 O 3 Pellets During Impregnation[END_REF][START_REF] Catita | Application of Magnetic Resonance Imaging and Raman Imaging to study the impact of phosphorus in impregnation of hydrotreatment catalysts[END_REF] According to the results of in situ EXAFS experiments by Nicosia et al., [START_REF] Nicosia | The effect of phosphate and glycol on the sulfidation mechanism of CoMo/Al 2 O 3 hydrotreating catalysts: an in situ QEXAFS study[END_REF] phosphates increase the saturation of active edges in CoMoS/γ-Al 2 O 3 catalysts. The authors rationalize their findings with support deactivation by phosphates. Several characterization studies have already been devoted to the intricate connection between the speciation of phosphate on the surface of the support as function of concentration, and its beneficial or deleterious effects in catalysis. X-ray absorption spectroscopy revealed that the presence of phosphate modifies the local environment of Mo at the metal-support interface. [START_REF] Nicosia | The effect of phosphate and glycol on the sulfidation mechanism of CoMo/Al 2 O 3 hydrotreating catalysts: an in situ QEXAFS study[END_REF] However, no direct insight into the phosphate local environment could be deduced from those studies. Although IR-based experiments on dry samples and samples in solution demonstrate an interaction between phosphate anions and γ-Al 2 O 3 carriers, it remains difficult to unambiguously assign the adsorption mode of phosphate anions. [START_REF] Zheng | Sorption of phosphate onto mesoporous γ-alumina studied with in-situ ATR-FTIR spectroscopy[END_REF][START_REF] Del | Surface reactivity of α-Al 2 O 3 and mechanisms of phosphate sorption: In situ ATR-FTIR spectroscopy and ζ potential studies[END_REF][START_REF] Chao Ren | Molecular-scale investigation of fluoride sorption mechanism by nanosized hydroxyapatite using 19 F solid-state NMR spectroscopy[END_REF] At the H 2 O/corundum interface, it was proposed that phosphate binds to corundum in a bidentate adsorption mode (two phosphate oxygen bind to the surface). [START_REF] Li | Molecular level investigations of phosphate sorption on corundum (α-Al 2 O 3 ) by 31 P solid state NMR, ATR-FTIR and quantum chemical calculation[END_REF] These studies rely on the comparison against the experimental spectra of known reference samples or against DFT studies, but this assignment remains difficult to achieve in the case of adsorbed phosphate anions, since the vibrational frequencies cannot be deconvoluted unambiguously.

Computational models of the atomic scale of the materials can be of great help to disentangle the findings of these studies. However, in order to do this, these models must be developed in close connection with experiments. It is in fact straightforward to construct a few model geometries and find patterns within them that match some of the experimental observations. The problem with this approach is that it might be biased towards the experimental observation, and this is ultimately not the same as an open-ended systematic search. In contrast, it is much more complicated to demonstrate that the model is plausible from a theoretical and experimental perspective. Even more important, this approach paves the way for predictions of experimental results and suggestions for new experiments. In the future, it might allow for the generation of new catalyst designs by computational modelling. For this purpose, a catalyst model needs to include all relevant components, and it must be shown that it is likely to form at synthesis conditions, but -most importantly -it needs to be consistent with all the experimental observations. At the same time, the model design needs to be feasible for testing with the available computational resources. For the ambitious goal -to fully understand the heterogeneous catalyst -one needs to know all interactions between the individual components. This work intends to elaborate on one of these interactions, that between phosphates and γ-Al 2 O 3 during impregnation and at drying conditions.

The first challenge in this endeavor is to construct a model that accurately reflects the experimental conditions. In order to avoid missing any potentially relevant structures and to keep the research open-ended, the initial search was focussed on finding a computational protocol that is feasible with the available resources, minimizes bias and is testable. For phosphate adsorption on an oxide surface, the number of adsorption sites largely depends on the complexity of the alumina surface. Moreover, phosphate itself is a complex complex patterns in the data. The difference between these two methods is quite subtle as they are both used to give a rationale for certain trends and to obtain predictions. This aims at a third challenge, which is to identify general trends in the produced data. In chapter 6 we tested several chemical descriptors and machine learning methods to construct NMR chemical shift prediction schemes. This prediction scheme is built using the geometries from chapter 3 by identifying the complex structure-chemical shift relationships. Chemical shift estimates are of comparable accuracy as DFT but orders of magnitude faster, reducing the calculation time from several hours to fractions of a second. This would allow [START_REF] Luis | Hydrolysis rates of inorganic polyphosphates in aqueous solution as well as in soils and effects on P availability[END_REF] P NMR predictions for larger adsorption systems and trajectories from molecular dynamics. Beyond the mere construction of a faster prediction system based on geometries, the presented machine learning methods can also give insight into the leading components that influence the chemical shift.

In its broader context, this work is based on ongoing scientific efforts that aim at an improved understanding of alumina surfaces and phosphate chemistry. Much work has been done in the area of solid state NMR spectroscopy, while experimental method development allows for more and more complex experiments. For computational modelling, a variety of bulk and surface models of transition alumina are available nowadays, some of them well established. [START_REF] Wolverton | Phase stability and structure of spinel-based transition aluminas[END_REF][START_REF] Gutiérrez | Theoretical structure determination of γ-Al 2 O 3[END_REF][START_REF] Ushakov | Structure of low-temperatureγ-andη-Al 2 O 3[END_REF][START_REF] Krokidis | Theoretical Study of the Dehydration Process of Boehmite to γ-Alumina[END_REF][START_REF] Paglia | Determination of the structure of gamma-alumina from interatomic potential and first-principles calculations: The requirement of significant numbers of nonspinel positions to achieve an accurate structural model[END_REF][START_REF] Digne | Hydroxyl Groups on γ-Alumina Surfaces: A DFT Study[END_REF] While several adsorption studies of phosphates and phosphonates on oxide surfaces exist, the mechanistic insights are still limited.

Mechanistic studies of phosphates are much better covered in molecular biology and one of the ideas tested in this thesis is if their methods can be transferred to mineral surfaces. The following chapter gives an overview of what has been done in the area of phosphate chemistry and applications in general as well as adsorption in particular.

That includes the currently available methods for sensitive [START_REF] Luis | Hydrolysis rates of inorganic polyphosphates in aqueous solution as well as in soils and effects on P availability[END_REF] P NMR experiments and the current knowledge about reaction kinetics of several reactions involving phosphates.

State of the Art

This chapter will first present the role of phosphates and transition alumina in different research areas as well as applications in several industrial sectors. The following sections will give an overview of many fields where the different phosphate interactions are important, followed by a review of the phosphate chemistry relevant in this work and the currently available spectroscopic methods. In particular, recent improvements in surface enhanced NMR techniques permitted complex [START_REF] Luis | Hydrolysis rates of inorganic polyphosphates in aqueous solution as well as in soils and effects on P availability[END_REF] P quantum correlation experiments and the generation of two-dimensional 27 Al -31 P NMR spectra. [START_REF] Hühn | Structural Characterization of Phosphate Species Adsorbed on γ-Alumina by Combining DNP Surface Enhanced NMR Spectroscopy and DFT Calculations[END_REF] The last sections review the current state of the available γ-Al 2 O 3 bulk and surface models.

Importance of phosphate-mineral interactions

Phosphate interactions with very different environments play a critical role in natural and human-made materials but our understanding of these systems is still limited. These materials include catalyst supports and minerals, which are important in the chemical industry and agriculture, as well as macromolecules in biological systems. Most phosphorus-containing chemicals are derived from phosphates, in which phosphorus binds to four oxygen atoms. The simplest and most common phosphates are based on the orthophosphate anion PO 3 - 4 and its protonated derivative HPO 2 - 4 , H 2 PO - 4 and H 3 PO 4 .

Orthophosphates are common in minerals and esters, typically with all protons substituted by cations (Ca 2+ in apatite, Al 3+ in berlinite) or organic moieties. Phosphorus is rather 1. State of the Art abundant (0.1 wt% of the Earth's crust, mostly in sediments of phosphate rock), [START_REF] William | CRC Handbook of Chemistry and Physics[END_REF][START_REF] Canfield | The Phosphorus Cycle[END_REF] and is involved in countless applications and natural processes. Its use in fertilizers, agrochemicals, detergents, corrosion inhibitors, food additives and pharmaceuticals gives it a critical role in modern civilization. Due to its range and scale of usage, phosphorus is listed as a critical raw material by the European Commission. [27] The synthesis of phosphorus chemicals is still wasteful and has a high potential for optimization towards more energy-efficient processes and minimization of by-products. [START_REF] Michael | Let's Make White Phosphorus Obsolete[END_REF] Moreover, phosphorus-based agrochemicals and detergents pose a yet incalculable risk to surface waters. [START_REF] Sandy | Oxygen isotope signature of UV degradation of glyphosate and phosphonoacetate: Tracing sources and cycling of phosphonates[END_REF][START_REF] Saidani | Toxicity assessment of organophosphorus in Ruditapes decussatus via physiological, chemical and biochemical determination: A case study with the compounds γ-oximo-and γ-amino-phosphonates and phosphine oxides[END_REF] In particular, the assessment of the environmental impact depends on our understanding of adsorption and degradation processes on soil minerals. Agriculture is certainly the most ubiquitous sector where the interaction of soluble phosphates with aluminosilica in soil has an essential contribution to crop yields and fertilizer consumption.

Torres et al. studied the bioavailability of phosphate, that is the availability of monomers (orthophosphate) in aqueous solution. [START_REF] Luis | Hydrolysis rates of inorganic polyphosphates in aqueous solution as well as in soils and effects on P availability[END_REF] They found different rates for polyphosphate hydrolysis and adsorption depending on the soil and the involved polyphosphates. These bioavailable phosphates are free inorganic phosphates in groundwater and water bodies like lakes, rivers and oceans. Eutrophication, due to fertilizers from mined phosphates that accumulate in those water bodies, also becomes an emerging concern. [START_REF] Michael | Let's Make White Phosphorus Obsolete[END_REF] Its symptoms are in some cases treated with alumina compounds. Phosphate management of lakes is already done for more than 50 years by precipitation with Al(OH) 3 . [START_REF] Agstam-Norlin | A 25-year retrospective analysis of factors influencing success of aluminum treatment for lake restoration[END_REF] Its success, however, depends a lot on the presence of other ions and aging characteristics which demands a well tuned, holistic strategy. [START_REF] De Vicente | Factors affecting phosphate adsorption to aluminum in lake water: Implications for lake restoration[END_REF][START_REF] Inmaculada De Vicente | Phosphate Adsorption by Fresh and Aged Aluminum Hydroxide. Consequences for Lake Restoration[END_REF][START_REF] Huser | Changes to water quality and sediment phosphorus forms in a shallow, eutrophic lake after removal of common carp (Cyprinus carpio)[END_REF] In this context it would certainly be of help to understand the adsorption kinetics of phosphates on Al(OH) 3 . Another recently published study emphasized the importance of free phosphates in order to produce accurate Earth system models for climate change prediction. [START_REF] Martiny | Biogeochemical controls of surface ocean phosphate[END_REF] Phosphate compounds are also one of many standard compounds needed in large scale industrial processes. For heterogeneous catalysts based on alumina substrate,

Importance of phosphate-mineral interactions

phosphoric acid is generally added during the impregnation step. It is reported that phosphoric acid increases the final catalytic performance in hydrodenitrogenation (HDN) as well as hydrodesulfurization (HDS) treatments. [START_REF] Iwamoto | Influence of Phosphorus on the Properties of Alumina-Based Hydrotreating Catalysts[END_REF] This behaviour is rationalized as a surface effect with the commonly used γ-Al 2 O 3 substrate as well as inducing changes to metalpolyanions in the impregnation solution. There is a general agreement in the literature that phosphates and phosphonates strongly interact with γ-Al 2 O 3 surfaces but it was also reported that phosphate is incorporated in dissolved metal precursors. [START_REF] Jaap | Envisaging the Physicochemical Processes during the Preparation of Supported Catalysts: Raman Microscopy on the Impregnation of Mo onto Al 2 O 3 Extrudates[END_REF][START_REF] Van Cleve | Enhanced Hydrothermal Stability of γ-Al 2 O 3 Catalyst Supports with Alkyl Phosphonate Coatings[END_REF][START_REF] Jaap | Spatially Resolved Raman and UV-visible-NIR Spectroscopy on the Preparation of Supported Catalyst Bodies: Controlling the Formation of H 2 PMo 11 CoO 5 -40 Inside Al 2 O 3 Pellets During Impregnation[END_REF] The presence of phosphorus in biomolecules like ATP, biopolymers like DNA or biominerals such as tooth enamel, underlines its existential biological role. [START_REF] Diskowski | Phosphorus". In: Ullmann's Encyclopedia of Industrial Chemistry[END_REF][START_REF] Schrödter | Phosphoric Acid and Phosphates[END_REF] Nucleic acids like RNA and DNA are phosphate condensation products, as well as adenosine polyphosphate like ATP and ADP for which there still remain many open questions regarding their reactivity, formation and contributions to biochemical processes. [START_REF] Westheimer | Why nature chose phosphates[END_REF][START_REF] Shina | Why nature really chose phosphate[END_REF] A more fundamental question that may also involve mineral surfaces and phosphates is, which chemical processes may have lead to the first proto-biochemical structures on

Earth. [START_REF] Georgelin | Formation of Activated Biomolecules by Condensation on Mineral Surfaces -A Comparison of Peptide Bond Formation and Phosphate Condensation[END_REF][START_REF] Kee | Phosphate Activation via Reduced Oxidation State Phosphorus (P). Mild Routes to Condensed-P Energy Currency Molecules[END_REF] According to Georgelin et al. a key piece for those early proto-biochemical phosphate condensation reactions could be reactions on alternating wet and dry silicate surfaces as produced at ocean shores by tides. [START_REF] Georgelin | Formation of Activated Biomolecules by Condensation on Mineral Surfaces -A Comparison of Peptide Bond Formation and Phosphate Condensation[END_REF] Kee et al. propose a reaction with a phosphite precursor in presence of Ca 2+ ions. [START_REF] Kee | Phosphate Activation via Reduced Oxidation State Phosphorus (P). Mild Routes to Condensed-P Energy Currency Molecules[END_REF] All these examples show research in very different contexts but the problems that were faced in them are similar throughout the areas of interest. Phosphates can strongly interact with mineral surfaces or even be a part of a mineral, while they simultaneously interact with several other components in the system, in solution or at dry surfaces.

One integral part of understanding any of those systems is to understand the role of phosphates as adsorbate species. Based on accurate phosphate adsorption models, one can add further components depending on the system of interest. This strategy should enable one to eventually build comprehensive models and understand these complicated systems in industry or nature in much more detail. (in this work it will always be referred to as metaphosphate). [START_REF] Thilo | Die kondensierten Phosphate[END_REF] Ring-shaped phosphate polymers are called polymetaphosphates based on the fact that the molar weights of their ions, [PO Overall, they are less prone to hydrolysis which occurs preferably at the chain ends. [START_REF] Luis | Hydrolysis rates of inorganic polyphosphates in aqueous solution as well as in soils and effects on P availability[END_REF][START_REF] Crowther | The Hydrolysis of the Condensed Phosphates: III. Sodium Tetrametaphosphate and Sodium Tetraphosphate[END_REF] Moreover, long polyphosphate chains slowly decompose at 60°C, thereby also forming trimetaphosphate. Thilo theorizes that therefore the polyphosphate chains are folded or in a helical structure. [START_REF] Thilo | Die kondensierten Phosphate[END_REF] Even higher degrees of condensation exist in form of cross-linked polyphosphates where ring or chain phosphates are phosphorylated in a third direction.

Importance of phosphate-mineral interactions

Those are characterized by a local chemical environment where, instead of one or two P -O -P bridging bonds, three bridging bonds start from the same phosphorus atom.

Consequently, these bonds are higher in energy and significantly less stable. They may occur in polyphosphate glasses at high temperatures and in absence of water. One small molecular representative of these phosphates is the phosphorus oxide which forms P 4 O 10 molecules, and the isotetrametaphosphate (fig. 1.1e) that is formed during hydrolysis of the former. [START_REF] Thilo | Die kondensierten Phosphate[END_REF] There is considerable experimental research done in the field of phosphate condensation, also including several accurate measurements of thermodynamic and kinetic data that were carried out over the last decades. [START_REF] Luis | Hydrolysis rates of inorganic polyphosphates in aqueous solution as well as in soils and effects on P availability[END_REF][START_REF] Crowther | The Hydrolysis of the Condensed Phosphates: III. Sodium Tetrametaphosphate and Sodium Tetraphosphate[END_REF][START_REF] Pedley | The Hydrolysis of the Condensed Phosphates: I. Sodium Pyrophosphate and Sodium Triphosphate[END_REF][START_REF] George | Squiggle-H2O". An enquiry into the importance of solvation effects in phosphate ester and anhydride reactions[END_REF][START_REF] Stockbridge | Enhancement of the Rate of Pyrophosphate Hydrolysis by Nonenzymatic Catalysts and by Inorganic Pyrophosphatase*[END_REF] 1.1.

Computational studies of phosphate condensation

From a computational perspective, phosphate condensation studies can involve a wide variety of condensation products. Early attempts (25 years ago) to computationally reproduce the reaction energy for pyrophosphate hydrolysis only showed qualitative agreement with experiments. [START_REF] Michael | Quantum chemical studies of pyrophosphate hydrolysis[END_REF] While accurate quantum mechanical methods were ) metadynamics approach. [START_REF] Tripathi | The GTPase hGBP1 converts GTP to GMP in two steps via proton shuttle mechanisms[END_REF] Compared to the early models, their study includes the statistical dynamic properties and resolves free energy differences as well as barriers.

Similar studies by Prasad et al. and Glaves et al. address the same hydrolysis reaction in

different biochemical environments and discuss the mechanistic details, like concerted or stepwise reactions. [START_REF] Glaves | Mechanistic Insights into the Hydrolysis of a Nucleoside Triphosphate Model in Neutral and Acidic Solution[END_REF][START_REF] Prasad | Addressing Open Questions about Phosphate Hydrolysis Pathways by Careful Free Energy Mapping[END_REF] Adsorption studies with inorganic material interfaces are, in comparison to biochemical and organic chemical studies, still at their infancy, also because

State of the Art

complex modelling techniques such as QM/MM are not always available due to the lack of parameters for heavy elements in force field parts. On alumina or oxidized aluminium surfaces, phosphate or phosphonate adsorption is commonly modelled as part of a range of compounds; for example, in context of lubricant or corrosion inhibition studies. [START_REF] Luschtinetz | Adsorption of phosphonic and ethylphosphonic acid on aluminum oxide surfaces[END_REF][START_REF] Blanck | Adhesion of lubricant on aluminium through adsorption of additive head-groups on γalumina: A DFT study[END_REF][START_REF] Milošev | Editors' Choice-The Effect of Anchor Group and Alkyl Backbone Chain on Performance of Organic Compounds as Corrosion Inhibitors for Aluminum Investigated Using an Integrative Experimental-Modeling Approach[END_REF] Hence, phosphate adsorption is embedded in context with several different adsorbing groups, but mechanistic details of adsorption were not in the scope of these studies.

To my best knowledge, no studies exist that simultaneously include oxide adsorption and condensation phenomena of phosphates. Phosphate condensation reactions are in general phosphorylation reactions, but these reactions could not only occur at the alumina surface but also with the surface itself. The following section will give a comprehensive explanation of the challenges involved in phosphorylation reactions.

The phosphorylation reaction

In a simplified context, phosphates may be considered as constant structural PO 4 units which is justified in perspective of the strong P -O bonds involved. However, several elementary reactions involve P -O bond breakage and formation. Among those are polyphosphate condensation and hydrolysis as well as formation of phosphate esters, both of which are critical reactions in biochemistry. In the context of surface speciation, phosphorylation reactions have not yet been investigated. However, they might be crucial in order to explain the formation of some surface species as during some of the preparation steps, the reaction conditions are likely to allow these reactions.

Phosphorylation reactions have some noticeable similarities with nucleophilic substitution reactions known in organic chemistry. A metaphosphate group (PO - 3 ) is moved from one moiety to the other. With an attacking group Y and a leaving group X, this can be written as:

XPO 2- 3 + Y ---→ X -+ Y PO 2- 3 (1.1)

Importance of phosphate-mineral interactions

For simpler equations, X and Y are assumed to be anionic moieties with one negative charge, and the metaphosphate moiety to be deprotonated. These conditions occur in strong alkaline environments. In neutral or acidic environments, the reaction could be:

XPO 3 H -+ HY --→ HX + Y PO 3 H - (1.2)
with X and Y as well as the metaphosphate protonated. By substituting X = OH - and Y = HPO 2- 4 in eq. 1.1 or 1.2, the equations describe pyrophosphate condensation;

swapping X and Y results in pyrophosphate hydrolysis. Instead, one can also substitute Y with a surface oxygen site, which would lead to the phosphorylation of that site. This would allow to add phosphates in those positions where a substitution of oxide ions would be otherwise unlikely. The protonation states in eq. 1.1 and 1.2 are a necessary simplification to express reaction equations, since the actual protonation states are generally not known. In fact, proton transfer reactions are critical components for understanding phosphorylations. [START_REF] Shina | Why nature really chose phosphate[END_REF] Most importantly, another reaction intermediate is possible. While nucleophilic substitution reactions are well understood in terms of eliminations (S N 1) or concerted reactions (S N 2), phosphorylations also allow for initial additions. In fact, stable penta-and even hexavalent oxo-phosphate complexes are well known. [START_REF] Kumara Swamy | New Features in Pentacoordinate Phosphorus Chemistry[END_REF][START_REF] Constant | New Trends in Hexacoordinated Phosphorus Chemistry[END_REF] This leaves many more possible reaction pathways than, for example, in case of nucleophilic substitutions. The overall reaction is either limited by the P-O bond breakage/formation or by proton transfer. one needs to define up to three more axes in general, which corresponds to up to five dimensions for searching transition states and intermediates. According to ref. [51], it is also important to ask for the number of involved H 2 O molecules in the proton transport,

X PO 2 - 3 [X -PO 3 -Y ] 3 - PO - 3 Y PO 2 - 3 + Y - -X - -X - + Y - + Y - -X - (a) Alkaline conditions X PO 3 H - H 2 PO 3 XY - ± H + PO - 3 Y PO 3 H - + HY ± H + , -HX -HX + HY , ± H + + H Y , ± H + ± H + , -H X (b) Neutral/acidic conditions
which adds yet another dimension to the reaction mechanism.

The most fundamental reaction is the hydrolysis of pyrophosphoric acid (H 4 P 2 O 7 )

to form two orthophosphoric acid molecules (H 3 PO 4 ) in aqueous solution. Several experimental studies have produced thermodynamic and kinetic data for this hydrolysis reaction. The first deprotonation of orthophosphoric acid and the first two deprotonations of pyrophosphate occur rapidly, hence the hydrolysis of H 2 P 2 O 2 - 7 to 2 H 2 PO - 4 is a reasonable assumption. At 25°C, the reaction free enthalpy and free enthalpy of activation of this hydrolysis are -32 and 119 kJ mol -1 , respectively. To my best knowledge, the only mechanistic studies of this particular reaction were conducted in the 1990s, neglecting solvation and finite temperature effects. [START_REF] Michael | Quantum chemical studies of pyrophosphate hydrolysis[END_REF] Later studies based on molecular dynamics or Plot created with experimental data from Yoza et al. [START_REF] Yoza | pH-dependence of 31 P-NMR spectroscopic parameters of monofluorophosphate, phosphate, hypophosphate, phosphonate, phosphinate and their dimers and trimers[END_REF] transition state search were focussed on specific, typically biochemical problems. [START_REF] Tripathi | The GTPase hGBP1 converts GTP to GMP in two steps via proton shuttle mechanisms[END_REF][START_REF] Glaves | Mechanistic Insights into the Hydrolysis of a Nucleoside Triphosphate Model in Neutral and Acidic Solution[END_REF][START_REF] Prasad | Addressing Open Questions about Phosphate Hydrolysis Pathways by Careful Free Energy Mapping[END_REF] Thus, our understanding of the fundamental mechanism of phosphoric acid condensation in aqueous solutions is still limited. However, since thermodynamic data of the overall reaction is available, this condensation reaction is simulated as a performance benchmark in chapter 5.

Solid-state 31 P NMR of phosphates

Nuclear Magnetic Resonance (NMR) is one of the most powerful techniques to reach an atomic scale determination of phosphate interaction with alumina. It is an elementselective method which allows to draw direct conclusions about the chemical environment around the phosphorus atom. Yoza et al. produced [START_REF] Luis | Hydrolysis rates of inorganic polyphosphates in aqueous solution as well as in soils and effects on P availability[END_REF] P NMR spectra for a broad range of phosphorus species in aqueous solution at a pH range from 3 to 11. [START_REF] Yoza | pH-dependence of 31 P-NMR spectroscopic parameters of monofluorophosphate, phosphate, hypophosphate, phosphonate, phosphinate and their dimers and trimers[END_REF] NMR spectroscopy on solids can be evaluated in similar ways as solution NMR. However, compared to solution NMR, the anisotropy of the chemical environment deteriorates the signal resolution of solid state NMR. The rapid molecular movement in solution produces an average isotropic environment where momentary anisotropy is "averaged out". In contrast, the orientation of the anisotropic chemical environment is constant in a dry powder sample. As a consequence, the solid state NMR spectrum is an ensemble of random anisotropic directions instead of well resolved averages. Orientation averaging is, however, possible by rotating the probe i. e. by rapidly spinning the probe head. Using an angle of about 55°relative to the magnetic field vector, all directions in three-dimensional space are transformed into one another (c. f. fig. 1.4), provided that the spinning frequency is fast enough. This can best be understood by imaging a cube spinning on its vertex.

The angle between any edge of the cube and the surface normal, or equivalently the cube's body diagonal, describes this magic angle. from Bleam et al. [START_REF] Bleam | 31 P solid-state nuclear magnetic resonance spectroscopy of aluminum phosphate minerals[END_REF] several aluminium phosphate minerals and estimated the anisotropy of the signal. [START_REF] Bleam | 31 P solid-state nuclear magnetic resonance spectroscopy of aluminum phosphate minerals[END_REF] The principal components of the anisotropy tensor (tab. of phosphoric acid and then drying at 120°C. [START_REF] Decanio | FT-IR and Solid-State NMR Investigation of Phosphorus Promoted Hydrotreating Catalyst Precursors[END_REF] A single broad peak is observed in the { 1 H} 31 P CP NMR spectra in fig. 1.5c, which is shifted to lower ppm with higher phosphate loading. The authors attribute this peak shift to an overlap of orthophosphate, et al. [START_REF] Decanio | FT-IR and Solid-State NMR Investigation of Phosphorus Promoted Hydrotreating Catalyst Precursors[END_REF] pyrophosphate and an increasing amount of other short-chained polyphosphate species with higher loading. At higher loading between 4 and 10 wt% phosphate, a distinct peak occurs at -32 ppm in the direct polarization [START_REF] Luis | Hydrolysis rates of inorganic polyphosphates in aqueous solution as well as in soils and effects on P availability[END_REF] P spectra in fig. 1.5b. This indicates the formation of a water-free layer of aluminium phosphate, which does not contain protons and is thus cannot be observed in the { 1 H} 31 P CP NMR spectra. In the direct polarization 27 Al spectra (see fig. °C. [START_REF] Van Eck | A Solid State Double Resonance NMR Investigation of Phosphorus-Impregnated γ-Al 2 O 3[END_REF] Using 31 P{ 27 Al} REDOR and TRAPDOR experiments, the authors also suggest the formation of both polyphosphates and a layer of aluminium phosphate. Despite the insights gained by those studies, there are still fundamental questions concerning the structure of the adsorbed species. For example, DeCanio et al. demonstrate that phosphates adsorb predominantly in the form of monomeric species at phosphate loading below 1 wt%, which corresponds to a coverage of approximately 0.3 monolayers. [START_REF] Decanio | FT-IR and Solid-State NMR Investigation of Phosphorus Promoted Hydrotreating Catalyst Precursors[END_REF] For higher phosphate loading it remains difficult to quantify at which point these species form oligomers and precursors of aluminium phosphate. Moreover, the local structure at the surface is still unknown. Yet, the dentation modes of phosphates adsorbed on alumina as well as changes to the surface could give very useful insights about the role of these species in catalysis.

Using surface signal enhancement techniques like Dynamic Nuclear Polarization Surface

Enhanced NMR Spectroscopy (DNP SENS) can vastly increase the signal gain which improves the signal-to-noise ratio. DNP SENS has revolutionized solid-state NMR of surfaces as the improved gain enables complex correlation experiments, which would take too long to record with conventional methods. The idea is to add an exogenous diradical solution to the sample, thereby creating a large electronic polarization on the sample surface. This electronic polarization is then transferred to the nuclei by irradiating the sample with a high-power microwave source, which results in improved signal gain of up to a factor of 100. [START_REF] Kan-Nian | Dynamic Nuclear Polarization with Biradicals[END_REF][START_REF] Qing Zhe Ni | High Frequency Dynamic Nuclear Polarization[END_REF][START_REF] Lesage | Surface enhanced NMR spectroscopy by dynamic nuclear polarization[END_REF][START_REF] Rossini | Dynamic Nuclear Polarization Surface Enhanced NMR Spectroscopy[END_REF] In combination with CP, DNP SENS yields sufficient signal intensity to record two-dimensional 27 Al-31 P correlation spectra, as well as [START_REF] Luis | Hydrolysis rates of inorganic polyphosphates in aqueous solution as well as in soils and effects on P availability[END_REF] P double quantum and triple quantum coherence experiments, to directly detect oligomers. [START_REF] Hühn | Structural Characterization of Phosphate Species Adsorbed on γ-Alumina by Combining DNP Surface Enhanced NMR Spectroscopy and DFT Calculations[END_REF] The determination of adsorption sites and modes at the surface of the alumina support is highly challenging due to a high degree of complexity of local sites. Computational modelling can be the key in order to complement experimental techniques.

Combined NMR and DFT studies

While the experimental work is critical to narrow down the requirements to the model, at this point none of the presented publications was able to produce a coherent and complete picture for phosphate-γ-Al 2 O 3 interaction. Recently, by the combination of DNP NMR and DFT calculations, it was possible to unravel the nature of Brønsted acid sites at the surface of amorphous silica-alumina. [START_REF] Valla | Atomic Description of the Interface between Silica and Alumina in Aluminosilicates through Dynamic Nuclear Polarization Surface-Enhanced NMR Spectroscopy and First-Principles Calculations[END_REF] Lately, it was also possible to determine the location of hydroxyl groups at the edges of γ-Al 2 O 3 crystallites thanks to the combination of 1 H NMR experiments and DFT calculations of chemical shifts. [START_REF] Ana | Beyond γ-Al 2 O 3 crystallite surfaces: The hidden features of edges revealed by solid-state 1 H NMR and DFT calculations[END_REF] Those DFT calculations rely on complex γ-Al 2 O 3 models detailed in the coming sections.

The requirement for these studies is an accurate prediction of NMR chemical shifts from DFT electronic structure calculations. For periodic calculations with pseudopotentials, the most common technique is the gauge-including projector-augmented wave (GIPAW) method, [START_REF] Pickard | All-electron magnetic response with pseudopotentials: NMR chemical shifts[END_REF][START_REF] Charpentier | The PAW/GIPAW approach for computing NMR parameters: A new dimension added to NMR study of solids[END_REF] that allows accurate chemical shift estimates for all atoms in the system.

Many theoretical studies, devoted to the calculations of phosphorus NMR chemical shifts inside various bulk materials, provide an accurate assignment of phosphate sites in various minerals. [START_REF] Ashbrook | Structure and NMR assignment in calcined and as-synthesized forms of AlPO-14: a combined study by first-principles calculations and high[END_REF][START_REF] Daniel | Effects of Extraframework Species on the Structure-Based Prediction of 31 P Isotropic Chemical Shifts of Aluminophosphates[END_REF][START_REF] Pourpoint | First Principles Calculations of NMR Parameters in Biocompatible Materials Science: The Case Study of Calcium Phosphates, β-and γ-Ca(PO 3 ) 2 . Combination with MAS-J Experiments[END_REF][START_REF] Pourpoint | Calcium Phosphates and Hydroxyapatite: Solid-State NMR Experiments and First-Principles Calculations[END_REF] Wei Li et al. showed by combining REAPDOR NMR characterization and GIPAW calculations that bridging bidentate phosphates would be adsorbed at the interface of water/boehmite (γ-AlOOH). [START_REF] Li | Solid-State NMR Spectroscopic Study of Phosphate Sorption Mechanisms on Aluminum (Hydr)oxides[END_REF] Nevertheless, since the authors used a rather simplified Al2 dinuclear cluster model for describing the boehmite surface, this may provoke questioning about the transferability of the results to other types of aluminum oxides such as γ-Al 2 O 3 . More recently, F. Tielens et al. investigated by GIPAW NMR calculations the [START_REF] Li | Molecular level investigations of phosphate sorption on corundum (α-Al 2 O 3 ) by 31 P solid state NMR, ATR-FTIR and quantum chemical calculation[END_REF] O and 31 P spectroscopic features of phosphate adsorption at the H 2 O-TiO 2 (anatase) interface and identified the bidentate adsorption mode. [START_REF] Tielens | Characterization of Phosphate Species on Hydrated Anatase TiO 2 Surfaces[END_REF] While GIPAW gives quantitative predictions of NMR chemical shifts, it does not provide a rationale for its origin. Dawson and Ashbrook produced a set of empirical equations involving average P -O distances and P -O -Al angles. [START_REF] Daniel | Investigating Relationships between the Crystal Structure and 31 P Isotropic Chemical Shifts in Calcined Aluminophosphates[END_REF] Multivariate regression of a test set of DFT calculations on an aluminium phosphate model cluster was used for function fitting. The authors later extended and applied their method based on experimental data and a wider range of aluminium phosphates. [START_REF] Daniel | Effects of Extraframework Species on the Structure-Based Prediction of 31 P Isotropic Chemical Shifts of Aluminophosphates[END_REF][START_REF] Daniel | A Multinuclear NMR Study of Six Forms of AlPO-34: Structure and Motional Broadening[END_REF] It is at the moment the only known model to map structural parameters of aluminium phosphates to NMR spectroscopic results.

GIPAW calculations require calculations at DFT level, with high demands towards precision of geometry and electronic structure. Thus, these simulations can be relatively expensive, in particular if the number of relevant nuclei is small (in the simplest case of orthophosphate adsorption, one electronic structure calculation produces one chemical shift). Empirical models, like Dawson and Ashbrooks, can also be a useful tool for reducing this computational expense. An early attempt to predict 31 P NMR chemical shifts by machine learning was done by West. [START_REF] Geoffrey | Predicting phosphorus NMR shifts using neural networks[END_REF] In 1993, he tried to train neural networks using a large database of experimental chemical shifts for phosphorus compounds, but the accuracy was not satisfactory. Recently, Chaker et al. tested several machine learning methods to predict NMR iso-shifts of [START_REF] Li | Molecular level investigations of phosphate sorption on corundum (α-Al 2 O 3 ) by 31 P solid state NMR, ATR-FTIR and quantum chemical calculation[END_REF] O, 27 Al and [START_REF] Sandy | Oxygen isotope signature of UV degradation of glyphosate and phosphonoacetate: Tracing sources and cycling of phosphonates[END_REF] Si in aluminosilicate glasses. [START_REF] Chaker | NMR shifts in aluminosilicate glasses via machine learning[END_REF] Their results are instructive in how to produce an accurate machine learning model with predictable error estimation relative to GIPAW results.

Building a γ-Al

2 O 3 model 1.2.1. Characteristics of γ-Al 2 O 3
Transition alumina is a versatile group of materials with well-characterized pore and particle sizes as well as high specific surfaces. [START_REF] Euzen | Alumina". In: Handbook of Porous Solids[END_REF] Among them, γ-Al 2 O 3 is economically the most important one since it is used as a support in large scale industrial catalysis.

It has a disordered oxide structure, as characterized by X-ray diffraction (XRD) and spectroscopic methods, which is a big obstacle for atomistic modelling. Nevertheless, several surface models have been tested and established in the past decades. These models allowed for the first time to build computational models to describe surface speciation on γ-Al 2 O 3 . In this work, γ-Al 2 O 3 model surfaces were used in order to establish a robust understanding of phosphate speciation at different experimental conditions.

As one of the most common supports in heterogeneous catalysts, γ-Al 2 O 3 is a critical material in large scale industrial processes. [START_REF] Toulhoat | Catalysis by Transition Metal Sulphides: From Molecular Theory to Industrial Application[END_REF][START_REF] Euzen | Alumina". In: Handbook of Porous Solids[END_REF] The properties of γ-Al another by temperature-controlled solid phase transformations. [START_REF] Euzen | Alumina". In: Handbook of Porous Solids[END_REF][START_REF] Busca | The surface of transitional aluminas: A critical review[END_REF] Partial pressure of water also influences the reaction process; even grinding and crushing can have an effect on the reaction as it improves water release for a given temperature and it might be a different energy source for phase transformations. [START_REF] Euzen | Alumina". In: Handbook of Porous Solids[END_REF][START_REF] Threlfall | Structural and Thermodynamic Explanations of Ostwald's Rule[END_REF] γ-Al 2 O 3 is formed as the first phase of boehmite dehydration. The sequence can be illustrated like this:

boehmite 450 °C ---→ γ-Al 2 O 3 750 °C ---→ δ-Al 2 O 3 1000 °C ----→ θ-Al 2 O 3 1200 °C ----→ corundum
From γ-Al 2 O 3 , the phase transitions to δ-Al 2 O 3 , θ-Al 2 O 3 and γ-Al 2 O 3 define the stability order, with γ-Al 2 O 3 as the least stable phase in accordance with Ostwald's rule. [START_REF] Van Santen | The Ostwald step rule[END_REF] More transition aluminas (χ-, κ-and η-Al 2 O 3 ) can be obtained by direct calcination of bayerite and gibbsite, forming different phase transition pathways. [START_REF] Euzen | Alumina". In: Handbook of Porous Solids[END_REF][START_REF] Busca | The surface of transitional aluminas: A critical review[END_REF] All transition aluminas form θ-and α-Al 2 O 3 at 1000°C. γ-, δ-and θ-Al 2 O 3 all share a similar distorted fcc crystal lattice of oxygen atoms. While θ-Al 2 O 3 is fully characterized and known to crystallize in a β-Ga 2 O 3 -type lattice, γ-, and δ-Al 2 O 3 are disordered systems with irregularly distributed aluminium atoms. This was shown by various techniques including XRD, NMR, IR and Raman spectra. [START_REF] Busca | The surface of transitional aluminas: A critical review[END_REF] They nevertheless show a certain degree of crystallinity, most notably a face-centered cubic (fcc) lattice of oxygen ions. 27 Al NMR results show that in γ-Al 2 O 3 , 25 -31 % of aluminium ions are in tetrahedral sites, [START_REF] Euzen | Alumina". In: Handbook of Porous Solids[END_REF] and this proportion is increased to 35 % in δ-Al 2 O 3 . For γ-Al 2 O 3 , new evaluations show that the number of tetrahedral

1.2. Building a γ-Al 2 O 3 model
sites might be slightly underestimated, and that the amount of tetrahedral sites is more similar to δ-Al 2 O 3 . [START_REF] Prins | On the structure of γ-Al 2 O 3[END_REF] Initially, γ-Al 2 O 3 was described as a distorted spinel-like structure and the high-temperature δ-Al 2 O 3 as a less disordered, more spinel-like structure. [START_REF] Lippens | Study of phase transformations during calcination of aluminum hydroxides by selected area electron diffraction[END_REF] This early observation is consistent with the distribution of aluminium in tetrahedral and octahedral sites, but the authors did not try to rationalize the pronounced anisotropy of γ-Al 2 O 3 , which would not occur in a cubic spinel structure. After calcination, both polymorphs also contain pentacoordinated aluminium atoms. [START_REF] Busca | The surface of transitional aluminas: A critical review[END_REF][START_REF] O'dell | A 27 Al MAS NMR study of a sol-gel produced alumina: Identification of the NMR parameters of the θ-Al 2 O 3 transition alumina phase[END_REF] On a mesoscopic level, it is known that γ-Al 2 O 3 is a porous compound that typically provides surface areas of 180 to 240 m 2 g -1 . The disorder in those systems substantially compromises the efforts to describe these materials on an atomistic level. Nevertheless, in the last decades, several research groups proposed a range of models for the description of γ-Al 2 O 3 and δ-Al 2 O 3 .

The following sections provide an overview of the available γ-Al 2 O 3 bulk models, as well as on how to produce surface models from them and the current understanding about this system.

Bulk model system

Disorder is a problem for any kind of simulation as the atomic structure must be fully defined first. In other words, we need to decide where to place the atoms to create an approximate representation. For the γ-Al 2 O 3 structure we only know two bulk properties with certainty: the oxygen atoms are ordered in a distorted fcc lattice and the aluminium atoms are distributed in approximately 1/3 tetrahedral and 2/3 octahedral sites. Beyond that, structural features similar to spinels can be imposed from XRD intensities. This description already dates back to six decades ago. [START_REF] Lippens | Study of phase transformations during calcination of aluminum hydroxides by selected area electron diffraction[END_REF] Kinetic considerations and other aspects, like the topotactic transformation from its precursor boehmite, indicate that γ-Al 2 O 3 should have a structural relation to boehmite as well.

The increase of computational capacities nowadays also allows for a systematic search for structures within a large configuration space. Structures based on all these assumptions [START_REF] Wolverton | Phase stability and structure of spinel-based transition aluminas[END_REF] The defect sites could be on tetrahedral or octahedral sites but ab initio calculations showed that octahedral defects are energetically favored, and tetrahedral sites are maximized, making up 37.5 % of all aluminium sites. As a slightly different approach, hydrogen spinels (HAl 5 O 8 ) were tested as well. [START_REF] Sohlberg | Hydrogen and the Structure of the Transition Aluminas[END_REF] Wolverton and Hass compared both models and concluded that HAl 5 O 8 is thermodynamically unfavorable compared to the hydrogen-free models and boehmite but it might exist as a metastable phase. [START_REF] Wolverton | Phase stability and structure of spinel-based transition aluminas[END_REF] Based on a structural relation to boehmite, Krokidis et al. constructed a 'dehydrated boehmite' model by removing OH and H from its layered structure and manually arranging a structural collapse followed by identifying sensible tetrahedral site transitions. [START_REF] Krokidis | Theoretical Study of the Dehydration Process of Boehmite to γ-Alumina[END_REF] The resulting structure is the first nonspinel structure (see fig. 1.6). They tested structures with 0 to 50 % tetrahedral sites; among them, those with 25 % lead to the most stable configuration, followed by those with 31 % and 37 % tetrahedral sites. The structure with 25 % was considered as the best case. This unit cell has a different shape and monoclinic symmetry, but almost the same cell volume as the defect spinel models with the same number of atoms. [START_REF] Ary | Ab initio 27 Al NMR chemical shifts and quadrupolar parameters for Al 2 O 3 phases and their precursors[END_REF] In an exhaustive comparison of spinel and nonspinel models, Paglia [START_REF] Paglia | Determination of the structure of gamma-alumina from interatomic potential and first-principles calculations: The requirement of significant numbers of nonspinel positions to achieve an accurate structural model[END_REF] The authors then compared their results with more than 500k systematically generated nonspinel models and concluded that, on average, 40 % of the aluminium atoms occupy nonspinel sites. Moreover, they ruled out spinel models based on the fact that their 100k generated variations generally produced mismatched diffraction patterns. This was in qualitative agreement with Krokidis' model, however since the latter is not based on a spinel unit cell, the number of spinel vs nonspinel sites cannot be clearly partitioned. [START_REF] Wolverton | Phase stability and structure of spinel-based transition aluminas[END_REF][START_REF] Gutiérrez | Theoretical structure determination of γ-Al 2 O 3[END_REF] hydrogen spinel, [START_REF] Ushakov | Structure of low-temperatureγ-andη-Al 2 O 3[END_REF] Krokidis nonspinel [START_REF] Krokidis | Theoretical Study of the Dehydration Process of Boehmite to γ-Alumina[END_REF] and Paglia nonspinel [START_REF] Paglia | Determination of the structure of gamma-alumina from interatomic potential and first-principles calculations: The requirement of significant numbers of nonspinel positions to achieve an accurate structural model[END_REF] ) with neutron diffraction experiments and 27 Al NMR spectroscopic data. [START_REF] Ary | Ab initio 27 Al NMR chemical shifts and quadrupolar parameters for Al 2 O 3 phases and their precursors[END_REF][START_REF] Paglia | Comment on "Examination of Spinel and Nonspinel Structural Models for γ-Al 2 O 3 by DFT and Rietveld Refinement Simulations[END_REF] The model fitting to neutron diffraction data showed that Paglia's nonspinel models were in the best agreement with Rietveld refinement simulations, followed by Krokidis' nonspinel, defect spinel and hydrogen spinel in this order, but with overall comparable accuracies. 27 Al-MAS NMR iso-shifts showed similar results, but in particular the absence of a correct prediction of pentacoordinated aluminium atoms indicated that surface aluminium atoms are likely to contribute significantly to the experiments. It must be underlined, that these comparative studies neglect surface effects. [START_REF] Ary | Ab initio 27 Al NMR chemical shifts and quadrupolar parameters for Al 2 O 3 phases and their precursors[END_REF][START_REF] Paglia | Comment on "Examination of Spinel and Nonspinel Structural Models for γ-Al 2 O 3 by DFT and Rietveld Refinement Simulations[END_REF] Even though this is generally a reasonable assumption, γ-Al 2 O 3 in particular is known to have large specific surfaces, increasing the number of contributing surface atoms. For example, Łodziana et al. showed in 2004 that certain hydroxylated alumina surface models of the related θ-Al 2 O 3 exhibit a slightly "negative" surface energy, which means that the surfaces are thermodynamically more stable than the bulk structure. [START_REF] Łodziana | A negative surface energy for alumina[END_REF] This is clearly impossible for clean surface cuts, however, in principle adsorbates (here: H 2 O) can stabilize the surfaces. It must be noted here that the negative surface energy reported here concerned only one of several surface models, and at the time the computational models neglected more accurate descriptions of long range dispersion interactions which are available today. [START_REF] Grimme | A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu[END_REF] Busca summarized that none of the four computational γ-Al 2 O 3 models can be rationally excluded on the basis of thermodynamic stability, since γ-Al 2 O 3 is the first transition alumina of boehmite, and neither the γ-Al 2 O 3 nor the δ-Al 2 O 3 phase are fully characterized. A very stable model intended for γ-Al 2 O 3 might rather be a good model δ-Al 2 O 3 . [START_REF] Busca | The surface of transitional aluminas: A critical review[END_REF] These considerations are in favor of Krokidis' model, since it is the only model involving a conceptual dehydration of boehmite and thus a construction that includes a structural relation to its precursor. Krokidis' model is also considered as one of the most popular ones, [START_REF] Busca | The surface of transitional aluminas: A critical review[END_REF][START_REF] Prins | On the structure of γ-Al 2 O 3[END_REF] while, for example, hydrogen-containing γ-Al 2 O 3 models are overall rejected. [START_REF] Busca | The surface of transitional aluminas: A critical review[END_REF] Prins recently reviewed the current knowledge of the γ-Al 2 O 3 structure. [START_REF] Prins | On the structure of γ-Al 2 O 3[END_REF] It was pointed out that new 27 Al NMR experiments and the previously available [START_REF] Li | Molecular level investigations of phosphate sorption on corundum (α-Al 2 O 3 ) by 31 P solid state NMR, ATR-FTIR and quantum chemical calculation[END_REF] O

NMR data indicate that higher percentages (35 -37.5 ) of tetrahedral sites should be expected; this indicates that the existing nonspinel models have too few aluminium atoms in tetrahedral sites. On the other hand, the author discussed findings of a decreased number of tetrahedral sites in γ-Al 2 O 3 catalysts. As a consequence, the concerned nonspinel models would still be useful to produce interface models. While Prins seems to disregard Krokidis' model based on his new analysis of NMR data, he acknowledged its successful application starting from Digne's hydrated surface studies, [START_REF] Digne | Hydroxyl Groups on γ-Alumina Surfaces: A DFT Study[END_REF] followed by many computational studies in heterogeneous catalysis. [START_REF] Hühn | Structural Characterization of Phosphate Species Adsorbed on γ-Alumina by Combining DNP Surface Enhanced NMR Spectroscopy and DFT Calculations[END_REF][START_REF] Raybaud | From γ-alumina to supported platinum nanoclusters in reforming conditions: 10 years of DFT modeling and beyond[END_REF][START_REF] Larmier | Tuning the Metal-Support Interaction by Structural Recognition of Cobalt-Based Catalyst Precursors[END_REF][START_REF] Wischert | γ-Alumina: The Essential and Unexpected Role of Water for the Structure, Stability, and Reactivity of "Defect" Sites[END_REF] While the author has a valid concern regarding the underestimated fraction of tetrahedral aluminium atoms, it does not exclude the Krokidis nonspinel per se. As mentioned above, Krokidis et al.

discussed a range of structures with up to 50 % tetrahedral sites, with models composed of slightly more tetrahedral alumina as additional candidates.

Generally, Busca and Prins concluded independently that structural disorder and limits in experimental methods did not yet allow for an atomic level characterization of γ-Al 2 O 3 . [START_REF] Busca | The surface of transitional aluminas: A critical review[END_REF][START_REF] Prins | On the structure of γ-Al 2 O 3[END_REF] [START_REF] Kovarik | High temperature transition aluminas in δ-Al 2 O 3 /θ-Al 2 O 3 stability range: Review[END_REF] The authors supported their proposal by DFT models and comparison with electron microscopy data.

Their results demonstrated how involved computational modelling along with experimental investigations allows for a comprehensive characterization of complex irregular phases such as in this case δ-Al 2 O 3 -the structurally most similar polymorph to γ-Al 2 O 3 . It should also be noted, that the δ-Al 2 O 3 phases are not based on a defect spinel structure.

Surface models

While the bulk structure of γ-Al 2 O 3 is of fundamental interest, its applications as a high surface compound (i. e. as a catalyst substrate), raise even more interest in its surfaces.

It is established that γ-Al 2 O 3 crystallites are products of a topotactic dehydration of boehmite and that the most exposed facet is (1 1 0) with about 70 % (cf. fig. 1.9). [START_REF] Euzen | Alumina". In: Handbook of Porous Solids[END_REF] To construct surfaces from the γ-Al 2 O 3 bulk model (fig. 1.6), the unit cell is cleaved along one crystallographic direction, and repeated several times along the new surface normal direction in order to mimic lower bulk layers, see fig. 1.7. Its thickness and the available vacuum over the surface become parameters that determine how accurate the surface is represented. However, this process is not straightforward as it has to be decided along which surface direction the surface is cut and which atoms are exposed on the hydration. [START_REF] Digne | Hydroxyl Groups on γ-Alumina Surfaces: A DFT Study[END_REF][START_REF] Digne | Use of DFT to achieve a rational understanding of acid-basic properties of γ-alumina surfaces[END_REF] The surface hydration can be summarized as follows: H 2 O adsorption occurs at the unsaturated surface aluminium atoms, followed by water splitting at the more acidic H 2 O sites that distributes the protons to the most basic alumina oxygen. For comprehensive details, refer to chapter 3. The surface hydration is complete when the aluminium coordination sites are occupied. Typically half of the H 2 O are dissociated.

Wischert et al. studied low hydration states of the (1 1 0) facet at H 2 O loads of at most 50 % H 2 O molecules. [START_REF] Wischert | γ-Alumina: The Essential and Unexpected Role of Water for the Structure, Stability, and Reactivity of "Defect" Sites[END_REF] They discovered that at these conditions, the (1 1 0) model facet

State of the Art

is subject to reconstruction by displacement of one aluminium atom to a neighboring tetrahedral site. It was recently shown that this aluminium atom is rather labile at aqueous conditions and might be subject to dissolution. [START_REF] Réocreux | Structuration and Dynamics of Interfacial Liquid Water at Hydrated γ-Alumina Determined by ab Initio Molecular Simulations: Implications for Nanoparticle Stability[END_REF] Since then, Digne's slab models were applied in several interface studies and good agreement with experimental results was attested. [START_REF] Blanck | Adhesion of lubricant on aluminium through adsorption of additive head-groups on γalumina: A DFT study[END_REF][START_REF] Ana | Beyond γ-Al 2 O 3 crystallite surfaces: The hidden features of edges revealed by solid-state 1 H NMR and DFT calculations[END_REF][START_REF] Larmier | Tuning the Metal-Support Interaction by Structural Recognition of Cobalt-Based Catalyst Precursors[END_REF][START_REF] Lagauche | Thermodynamic Characterization of the Hydroxyl Group on the γ-Alumina Surface by the Energy Distribution Function[END_REF] Prins criticized the commonly studied (1 1 0) facet introduced by Digne et al. because it is different from the spinel-like (1 1 0) surfaces, and does not produce considerable surface reconstruction due to its stoichiometric termination. [START_REF] Prins | On the structure of γ-Al 2 O 3[END_REF] Pinto et al. observed this reconstruction on defect spinel models, however, surface hydration was neglected at the time . [START_REF] Henry | Ab initio study of gamma-Al 2 O 3 surfaces[END_REF] Prins also noted that, since γ-Al 2 O 3 is formed in a topotactic phase transition from boehmite, different surface terminations should be tested, as the cleaved model surface by Digne is not the only candidate. A recent study systematically explored these different surface terminations which should be expected by construction from boehmite dehydration. [99] While Digne's surface models are not the final answer, they remain useful to study possible surface chemistry at γ-Al 2 O 3 facets.

Surface adsorption

The goal of this study is to obtain a set of candidate structures that can account for the adsorption modes obtained after impregnation in solution and drying. For that, one first needs to know the predominant surface facets available for adsorption. About 70 % of the facets are (1 1 0) and after that (1 0 0), as established in section 1.2.3. [START_REF] Euzen | Alumina". In: Handbook of Porous Solids[END_REF] For exploration, the drying state was considered, which circumvents difficult modelling of solvent effects. The initial assumption is that after drying, phosphates only bind to aluminium atoms via oxygen sites that are accessible at the surface. The fully hydrated surface models established by Digne et al. were used throughout this work. [START_REF] Digne | Hydroxyl Groups on γ-Alumina Surfaces: A DFT Study[END_REF] In their publication, the hydration state of several surface terminations was tested at varying temperatures. According to ref. [95], the (1 0 0) and (1 1 0) facets are each hydrated at 120°C, which also shown in fig. 1.10. The more recently discovered surface reconstruction of (1 1 0) first described by Wischert et al., [START_REF] Wischert | γ-Alumina: The Essential and Unexpected Role of Water for the Structure, Stability, and Reactivity of "Defect" Sites[END_REF] will be shown as crucial part of this study.

The model adsorption of phosphates on the surface can be carried out with various strategies, but here the general idea was to substitute the surface hydroxyls in a reaction similar to

H 3 PO 4 + OH --→ H 2 PO 4 + H 2 O
where is the adsorption site of the hydroxy group. A similar computational adsorption study was done by Luschtinetz et al. involving phosphonates (H 2 RPO 3 ) and several aluminium oxides and hydroxides. [START_REF] Luschtinetz | Adsorption of phosphonic and ethylphosphonic acid on aluminum oxide surfaces[END_REF] The authors used the substitution approach with hydroxy groups on these surfaces. In particular, they explored a wide range of adsorption modes of two phosphonates on corundum, bayerite and boehmite surfaces by a systematic substitution of surface oxygen and hydroxide sites. The challenge there was similar to this adsorption study. However, adsorption on γ-Al 2 O 3 is more complex because it has

Methods

The results produced for this thesis are exclusively computational models of γ-Al 2 O 3 surfaces with the ambition to create a detailed understanding of phosphate adsorption on this material. Electronic structure calculations and construction principles of interface models are well established, however they rely on numerous approximations and assumptions that balance the accuracy of the results with computational performance.

This chapter gives a short description of the applied methods and lists all important parameters used. A comprehensive discussion of the range of available methods as well as accuracy considerations are given in appendix A.1 (p. 143).

Computational setup

The presented models have have to be assessed from three distinct point of views: (i) electronic structure, (ii) structural aspects or geometry of atom positions, (iii) system sampling. While all electronic structure calculations were done on a similar level of density functional theory (DFT), there are several important differences in the technical details.

The surface models were adjusted to answer different questions efficiently. Sampling has several meanings for the distinct challenges: for the database presented in chapter 3 and 4 it refers to a part of a global optimization to ensure that the identified geometry minimizes the energy; in contrast the sampling in chapter 5 and 6 is a central contribution in those results, where it contributes to the machine learning quality and free energy at

Methods

the solid-liquid interface. The following sections summarize the employed approximations on these levels separately.

Electronic structure of periodic DFT

All calculations were conducted using at PBE-D3 level of theory, [START_REF] Grimme | A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu[END_REF][START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF] which is a very common DFT exchange correlation potential based on the general gradient approximation that is combined with a force-field like term included to add long range dispersion interactions.

Small systems include single surface unit cells with unit cell parameters for the two γ-Al 2 O 3 facets: (1 0 0) 8.358Å×5.547Å and (1 1 0) 8.358Å×8.035Å. Each system contains around 100 atoms. The larger systems are supercells constructed with two to three times larger unit cell vectors for about 16Å×16Å sized surface unit cells. Together with a 17Å thick layer of water molecules, they contain 800-1000 atoms. The small systems allow reaching higher accuracies for of electronic structure calculations than the latter.

A high accuracy was important here, since the systems were also used to calculate NMR chemical shifts which are very sensitive the local structure. The larger systems are surface supercells build for molecular dynamics, where fast calculations are much more important.

For small systems with higher accuracy requirements, the electron structure was spanned in projector-augmented wave functions (PAW) with a kinetic energy cutoff at 500 eV and a density cutoff of 605 eV. [START_REF] Blöchl | Projector augmented-wave method[END_REF] The PAW pseudopotentials are based on a norm conserving reference density, which allowed such a smooth density cutoff. [START_REF] Kresse | From ultrasoft pseudopotentials to the projector augmented-wave method[END_REF] For efficient calculations, VASP (v. 5.3.5) was used. [START_REF] Kresse | Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set[END_REF] Other parameters that can influence the electronic structure are smearing (set to 10 meV Gaussian-type smearing which is negligible at the end of a geometry optimization but stabilizes certain stating configurations) and two switches determining overall precision (PREC = accurate) and for The large systems were constructed as supercells with 16×16Å wide surface area, on which a 17Å thick layer of water was equilibrated using molecular dynamics. The Brillouin for these supercells is decreased such that calculations based only on the Γ-point have equivalent accuracy as the 2×3 and 2×2 grids in the small systems. The calculations were carried out with the QUICKSTEP module which is part of the CP2K software package. [START_REF] Kühne | CP2K: An electronic structure and molecular dynamics software package -Quickstep: Efficient and accurate electronic structure calculations[END_REF] QUICKSTEP uses a Gaussian and plane wave (GPW) method to represent Kohn-Sham orbitals and electron density. The difference to PAW calculations is explained in more details in appendix A. 1.3 (p. 151). Norm-conserving GTH pseudopotentials were used together with molopt-DZVP basis sets. [START_REF] Kühne | CP2K: An electronic structure and molecular dynamics software package -Quickstep: Efficient and accurate electronic structure calculations[END_REF][START_REF] Goedecker | Separable dual-space Gaussian pseudopotentials[END_REF] The auxiliary density cutoff is set to 400 Ry (5442 eV). which is a bit higher than the software default of 280 Ry. The value has to be for higher accurate density transformation because the local basis used for orbitals in CP2K is significantly harder; see appendix A.1.3 (p. 151) for a detailed explanation.

DFT chemical shift calculations

NMR chemical shift calculations were carried out using the gauge-including projectoraugmented plane wave method (GIPAW) by Yates, Pickard and Mauri. [START_REF] Pickard | All-electron magnetic response with pseudopotentials: NMR chemical shifts[END_REF][START_REF] Yates | Calculation of NMR chemical shifts for extended systems using ultrasoft pseudopotentials[END_REF] The chemical shift tensor is defined by the ratio of induced magnetic field and external field:

σ(r N ) = B ind (r N ) B ext (2.1)
After careful convergence tests, the Monkhorst-Pack k-point grids were increased to 3×5×1 and 3×3×1 for the surface models. The calculation involves numerical k-space 41

Methods

derivatives, for which step sizes of DQ = 0.003 were used and symmetry operations were discarded where necessary (LNMR_SYM_RED=.TRUE). Furthermore, quadrupolar coupling constants (C q ) were estimated using quadrupole moments of 146.6 and -25.58 for 27 Al and 17 O. To ensure accurate electronic convergence SCF convergence was continued until consecutive values differed less than 10 -8 eV (EDIFF).

Phase space sampling

Theoretical models for phosphoric acid adsorption on γ-Al 2 O 3 model facets require as minimum a notion on how a phosphate group can adsorb to this interfaces. This notion can be translated into several adsorption candidates for which different sites are identified on the surface models. Depending on the number of candidates and our knowledge about the system, there are many different sampling methods for atomistic structures.

Starting from manually testing individual adsorption mode with hand-made geometries is common, and typically followed by exploring all remaining adsorption modes if feasible, or at least several of them of the exploration space is too large. As alternative, one can test automated sampling methods. Popular sampling methods include ab-initio molecular dynamics (AIMD), Monte Carlo sampling or genetic algorithms, [START_REF] Deaven | Molecular Geometry Optimization with a Genetic Algorithm[END_REF][START_REF] Andriy | New developments in evolutionary structure prediction algorithm USPEX[END_REF] which are more or less suitable for the actual challenge. AIMD is often used to sample gas phase or liquid phases, whereas Monte Carlo sampling and genetic algorithms can be used for global optimizations with little to no initial knowledge about the system. When including structural knowledge into the sampling, the situation is much more complicated, since the choice on which parts to sample (i. e. adsorption sites and reactions) and which to restrict (i. e. the surface model) have to be explicitly implemented. AIMD allows accurate sampling of water at the γ-Al 2 O 3 model interfaces, [START_REF] Réocreux | Structuration and Dynamics of Interfacial Liquid Water at Hydrated γ-Alumina Determined by ab Initio Molecular Simulations: Implications for Nanoparticle Stability[END_REF] which also allows to test adsorption models, however at high computational cost.

In this work, the adsorption sites were sampled by manual structure generation based on a detailed analysis of surface patterns and a resulting combinatoric sample space.

In-depth considerations are given in sections 3.2 and 3.3 (p. 58). For the solid-liquid interface, some of the most stable adsorption sites were simulated by means of AIMD.

In this context each sampled structure as well their AIMD simulations represent one spot of the phase space. Even though AIMD in principle allows for dynamic change between adsorption sites, the time scales of these simulations are in a range of 10 ps to 50 ps which is just enough to sample liquid water but not to observe scarce transitions.

In order to connect these spots, or in other words, to study the speciation mechanism at the impregnation conditions, biased AIMD calculations were carried out. The biasing approach applied here is metadynamics or metaD for short. This section also addresses the free energy and free enthalpy estimates, since in contrast to the electronic energy, which is only defined for a single geometry, it is directly related to the phase space by the degrees of freedom available to the system.

Geometry optimization

More than 1000 geometry optimizations were carried out within the structural exploration phase of this work. To focus on structure generation rather than correcting unsuccessful optimizations, the optimization algorithm had to be more reliable than any of the standard implementations in VASP. The fast inertial relaxation engine (FIRE) [START_REF] Bitzek | Structural Relaxation Made Simple[END_REF] was implemented in VASP by Henkelman et al. for their nudged-elastic-band (NEB) [START_REF] Henkelman | Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points[END_REF] based transition state search methods, but it can also be used separately. It is the only method that efficiently optimizes internal rotamers like free Al -OH and P -OH without user intervention, since it is only follows the force vector and adjusts its step size during the optimization. The quasi-Newton method in VASP (IBRION=1) [START_REF] Pulay | Convergence acceleration of iterative sequences. the case of scf iteration[END_REF] does not succeed for reasons explained in ref. [109]. The success of the optimizing algorithm also depends on the precision with which the force were determined. This depends the

Methods

electronic optimization, which in this case done until 10 -6 eV (EDIFF), while the forces were optimized until the largest force on all atoms was below 10 meVÅ -1 .

Partial Hessian matrix and harmonic free energies

Frequency calculations were carried out on the best candidate geometry for each adsorption site. Since the data was acquired to compute vibrational component of the free energy computations, the parameters were optimized as well. For the description of the relevant vibrational surface modes, partial Hessian matrices were calculated by finite differences.

The electronic optimization was done until 10 -8 eV, displacements of 0.015Å were used for all atoms. To limit the number of atom displacements to the relevant area only surface atoms and adsorbates were included. Without further processing, this partial Hessian matrices contain frustrated translations due to the neglected subsurface atoms. These are difficult to identify by hand, but they can be removed by projection, where they were treated as free translations. The remaining partial Hessian matrix were diagonalized and the resulting vibrational modes were checked for imaginary modes. If imaginary modes occurred, the structures were reoptimized with tighter force criteria until the imaginary mode vanished.

The vibrational frequencies were furthermore used to calculate thermodynamic correction terms for the Gibbs free enthalpy G. The overall free enthalpy is defined as

G = E el + U + pV -T S (2.2)
where contributions to the inner energy U and entropy S of free molecules are separated in translational, rotational and vibrational components:

U = U trans + U rot + U vib (2.3) S = S trans + S rot + S vib (2.4)
Each component can be determined by constructing their partition function. [START_REF] Mcquarrie | Statistical Mechanics[END_REF] They whole expression for the free enthalpy is thus:

G = E el + U trans + U rot + U vib + pV -T S trans + S rot + S vib (2.5)
It is necessary to calculate all terms for gas-phase reference molecules like H 2 O and H 3 PO 4 . For the surfaces, translational and rotational components vanish, as well as pV since surfaces do not have a partial pressure, leaving only vibrational contributions:

G surf = E el + U vib -T S vib (2.6)
The lowest vibrational frequencies have the most significant contribution to the overall thermodynamic correction terms, which is the reason why it was important to remove the frustrated translations discussed above. In connection with the phase space, the calculation and diagonalization of a Hessian matrix for harmonic vibrational mode is nothing else than a simple analytical alternative to sampling.

Molecular dynamics

The explicit solvation of the surface and surface species requires to include dynamic, statistical effects of the solvent. For this, AIMD simulations were carried out at DFT level as described in section 2.1.1. Room temperature calculations where simulated with a slightly elevated simulation temperature of 330 K using the CSVR thermostat (100 fs time constant) in order to compensate for systematic errors known about PBE-D3. [START_REF] Bussi | Canonical sampling through velocity rescaling[END_REF] In order to increase simulation throughput, the hydrogen atomic weight was increased to 3 g mol -1 which allowed to carry out the simulations with a step size of 1 fs. These parameters were in all AIMD. Unbiased AIMD calculations were carried out for at least 10 ps to equilibrate the system before starting metadynamics simulations. Four parameters, d 0 , r 0 , n and m define its behavior.

s(r) = 1- r-d 0 r 0 n 1- r-d 0 r 0 m for r > d 0

Bias AIMD along collective variables

Biased AIMD simulations are a critical tool to study rare events. AIMD simulations are carried out for time scales in the range of picoseconds, where the probability of most chemical reactions is very low. Biases are used to nudge the system over reaction barriers.

There are many possible methods to do so, but regardless of the actual method, one first needs to define a direction within phase space in which a bias should be applied. It is important to note that really any differentiable system quantity can be used to add biases to the system. When defined, a bias exerts an additional force at the atoms in the system. The simplest variable to bias for chemical bond formation and cleavage is the distance between the atoms of interest. For large systems, the difficult question is, which distance to consider. For example in case of P -O bond breakage in PO 3 - 4 , one can define any of the four bonds, whereas for bond formation with H 2 O, oxygen atoms in several surrounding H 2 O molecules could be chosen. This complexity cannot be added by means of a distance, which is only defined for a single pair of atoms, but coordination is a much more suitable concept. In general, the coordination c is a collective variable (CV) which counts the number of bonds within the

same group: c = N i=1 N j>i s(r ij ) or different groups: c = N i=1 M j=1 s(r ij ) (2.7)
A bond can be defined as an atom pair with a distance r ij below a user-defined cutoff distance. For a continuously differentiable definition, a smooth transition between bonded and non-bonded pairs, a switching function has to be established. A common switching function for the coordination CV is:

s(r) =        1- r-d 0 r 0 n 1- r-d 0 r 0 m , if r > d 0 1 , otherwise (2.8)
The function is plotted in fig. can be used as a cutoff value: s(d 0 +r 0 ) = n m = 0.5. The definition of parameters is an involved process which required months of testing, before identifying a set of parameters that performs well. The most flexible set of parameters identified only depended on r 0 ; the other parameters were fixed at d 0 = 0, n = 6, m = 12. For the coordination of P -O, r 0 is set to 2.2Å, for Al -O to 2.4Å, and for O -H to 1.4Å. This parameters are quite specific to the targeted element pairs to be split or formed but with local effects, like the coordination of octahedral and tetrahedral aluminium sites, or differences between µ 1 -OH and µ 3 -O, these parameters are not necessarily ideal in all situations.

These functions are implemented in the PLUMED code, [START_REF] Bonomi | PLUMED: A portable plugin for free-energy calculations with molecular dynamics[END_REF][START_REF] Gareth | PLUMED 2: New feathers for an old bird[END_REF][START_REF] Bonomi | Promoting transparency and reproducibility in enhanced molecular simulations[END_REF] which is compatible with CP2K. The functionality is available with the COORDINATION keyword. It should also be noted, that it is used as a general method to count the number of bonds within 2. Methods groups of atoms, which different from coordination numbers. Those describe the number of other atoms to the same atomic center, however the coordination number of a single atom can be described for one center, if this one is the only element in the group. For brevity, this section just covered the most basic CVs used the most within chapter 5.

switching coordination function parametrized in this work to mathematically describe coordinations as a sum. It should be noted that any differentiable function can be used as a switching function switching function can be used and more complex quantities than the sum in eq. 2.7 are needed in some cases. Those are described in appendix A.4 (p. 162).

Metadynamics

There are several methods to add a bias to AIMD along a set of CV's. In this work, metadynamics (metaD) [START_REF] Laio | Escaping free-energy minima[END_REF] and, more specifically, well-tempered metadynamics (WT-metaD) [START_REF] Barducci | Well-Tempered Metadynamics: A Smoothly Converging and Tunable Free-Energy Method[END_REF] were carried out to study phosphate adsorption at the solid-liquid interface.

The idea in these methods is to deposit a potential V (s, t) in form of gaussian functions (in jargon hills) along positions of a CV s in order to force the system to leave this particular area in the free energy surface (FES) F (s):

V (s, t) = t τ =0 W (τ ) exp - (s -s τ ) 2 2σ 2 (2.9)
The potential is given as a sum over τ gaussian during the simulation time t, where the σ parameter sets the gaussian width. The gaussian height in standard metadynamics is a constant. It can be shown that the potential V(s,t) converges with the FES shifted by a constant. [START_REF] Dama | Well-Tempered Metadynamics Converges Asymptotically[END_REF] The process is illustrated in fig. 2.2. For WTmetaD, the gaussian height is decreased depending on the existing bias potential V (s, t) at the point of the next gaussian deposition at time t:

W (t) = W 0 exp - V (s, t) ∆T (2.10)
where ∆T is another input parameter. The bias factor allows for an overall better convergence behavior, since the growth rate of V (s, t) is limited, but even at full convergence, V (s, t) does not fully compensate the FES:

V (s, t → ∞) = - ∆T ∆T + T F (s) + C (2.11)
Instead, V (s, t) would have to be scaled by a factor γ = -(T + ∆T )/T . This factor γ is referred to as bias factor and most commonly reported in publications instead of values for ∆T . Overall, one needs to define a gaussian height W 0 and width σ, a deposition pace 1/t and a bias factor γ for a one-dimensional WTmetaD simulation. For further dimensions, i. e. to bias multi-dimensional FES projections along further CV's, one needs to define additional width parameters along the CV directions. In this work, for all

Methods

presented WTmetaD simulations, a deposition pace of 100 ps -1 was used, which is every 10 simulation steps. The bias factor was set to γ = ∆E k B T = 72.9 assuming 200 kJ mol -1 is the largest barrier needed to be overcome. This bias factor is smaller than the one used in previous work, [START_REF] Réocreux | Structuration and Dynamics of Interfacial Liquid Water at Hydrated γ-Alumina Determined by ab Initio Molecular Simulations: Implications for Nanoparticle Stability[END_REF] and therefore expected to allow slightly faster convergence. It also is the default method to determine the best bias factor in a new derivative method of WTmetaD. [START_REF] Invernizzi | Rethinking Metadynamics: From Bias Potentials to Probability Distributions[END_REF] AIMD simulations are generally very time-consuming calculations, where WTmetaD is no exception. The discussion above considers the modification of the AIMD trajectory by constructing additional energy terms which generate forces in forthcoming simulation steps. The overall bias potential can only build up in areas of this single trajectory. The computation time that the trajectory, in this context referred to as walker, requires to cover the overall area defined by the CV's is the minimum calculation time. By adding more walkers that simultaneously contribute to the bias potential, one can parallelize the whole process. In fact, the original publication of the metadynamics (metaD) method describes this parallelization technique. [START_REF] Laio | Escaping free-energy minima[END_REF] The most obvious advantage of this parallelization is reduced computation time at the cost of higher parallel workload. A more subtle advantage is that the parallel walkers, while ultimately sampling different areas of the FES, they can also limit initial oscillation effects from an alternating coverage of different areas in the FES, if the starting points are distributed accordingly. For example 100 walkers were used in ref. [120]. In case of AIMD, the limiting factor is the overall available resources, hence for this work, the number of walkers was always limited to 4.

Structural exploration

The aim of the work presented in this chapter is a detailed exploration of all relevant adsorption patterns of H 3 PO 4 on hydrated γ-Al 2 O 3 . This exploration also sets the foundation for all derived work, that will be discussed in chapters 4, 5 and 6. Compared to other adsorption studies, [START_REF] Li | Molecular level investigations of phosphate sorption on corundum (α-Al 2 O 3 ) by 31 P solid state NMR, ATR-FTIR and quantum chemical calculation[END_REF][START_REF] Luschtinetz | Adsorption of phosphonic and ethylphosphonic acid on aluminum oxide surfaces[END_REF] in which all surface oxygen sites were µ 2 -OH sites, the hydrogen networks of γ-Al 2 O 3 surface models are more complex. They involve partially protonated µ 3 -O sites and µ 1 -OH as well as µ 1 -H 2 O sites with rotational degrees of freedom. This requires an extended search for the protonation state of every generated system, from reference to adsorption candidate. Overall, three reference structures were used as a starting point for the adsorption states of phosphates on γ-Al 2 O 3 : two Digne's original surface models, [START_REF] Digne | Use of DFT to achieve a rational understanding of acid-basic properties of γ-alumina surfaces[END_REF] here referred to as (1 0 0) and n(1 1 0), and Wischert's reconstruction, R(1 1 0). [START_REF] Wischert | γ-Alumina: The Essential and Unexpected Role of Water for the Structure, Stability, and Reactivity of "Defect" Sites[END_REF] These reference structures were reoptimized at the applied level of theory, including an initial, thorough exploration of the proton network. The most stable configurations are illustrated in fig. 3.1b. The sites relevant for chemisorption in these structures are oxygen atoms highlighted in red for sites that originate from surface hydration and blue for the topmost layer of alumina oxygen atoms. This distinction is arbitrary but a quite useful first hierarchical order in which phosphates may adsorb by substitution at these sites: first substitution at hydration sites, then eventually sites of the alumina network. The structures in fig. 3.1b reflect the chemical bonding situation accurately and even vaguely map the hydrogen bond network of each surface, however, but they are too complex to systematically address the exploration of different adsorption patterns. Fig. 3.1a is an abstract representation of the top view of the three surfaces that includes the coordination polyhedrons of the aluminium atoms and labels for the highlighted oxygen atoms; lowercase labels for hydration sites, uppercase labels for alumina network sites. Moreover, this surface representation includes symmetry aspects of Digne's (1 0 0) and n(1 1 0), see fig. 3.1i and ii. Except for the hydrogen network, these symmetries overall reduce th number of possible hydration states, which considerably reduces the number of adsorption modes. Thus also the labels were adjusted. Symmetryequivalent sites were provided with the same letters and one marked with a prime. For example on fig. With the involved challenges in finding the relevant adsorption modes, the next three sections explain in detail, how all adsorption structures were sampled. First, some general considerations are laid out, including an explanation on how structural exploration is related to phosphate speciation at drying conditions and how the adsorption energy can be interpreted at the solid-liquid interface. This is followed by a detailed description and a notation for all the adsorption sites. After that, the most stable adsorption sites are analyzed and construction principles for them are identified.

Drying conditions and structural exploration

For the purpose of covering all the targeted adsorption patterns, one needs to choose the relevant environmental conditions. This choice is also strictly necessary in order to quantitatively compare different adsorption candidates by their stability. The preparation steps of a heterogeneous catalyst starting from the calcined γ-Al 2 O 3 support are wetincipient impregnation, drying and activation. Of those, we are mostly interested in impregnation (solid-liquid interface, 25°C, 1 atm) and drying conditions (solid-gas interface, 120°C and 2 kPa for H 2 O -air moisture levels of 60 % at ambient conditions). (1 0 0) [START_REF] Digne | Hydroxyl Groups on γ-Alumina Surfaces: A DFT Study[END_REF] (leftmost structures), n(1 1 0) [START_REF] Digne | Hydroxyl Groups on γ-Alumina Surfaces: A DFT Study[END_REF] (centered structures) and R(1 1 0) [START_REF] Wischert | γ-Alumina: The Essential and Unexpected Role of Water for the Structure, Stability, and Reactivity of "Defect" Sites[END_REF] (rightmost structures). Oxygen atoms are highlighted in all structures to distinguish the sites of preadsorbed water (red) and alumina oxygen (blue). The structures in (a) provide the necessary level of abstraction from (b), and in particular the labelling that is critical to follow the discussion in this work in detail.

Starting from the impregnation conditions is impractical for structural exploration, since it would require to simulate liquid water at the surface. Sampling of liquids using molecular dynamics and comparing stabilities is computationally demanding and complicated. In contrast, at drying conditions the surface model only involves adsorbed water with comparably well-defined positions. The interaction with the gas phase over the surface can be neglected, and the translational and rotational degrees of freedom in the gas phase itself can be described with analytical terms within ideal gas and rigid rotor approximations (cf. eq 2.5, p. 45 and ref. [112]). Thus, drying conditions can be simulated by creating static models of surface adsorption, and comparing their energies. Only H 2 O (g) has to be considered as a gas-phase reference state, as will be shown below. The adsorption enthalpy of H 2 O was well described by Digne et al. using only (nowadays) inexpensive static calculations. [START_REF] Digne | Hydroxyl Groups on γ-Alumina Surfaces: A DFT Study[END_REF] Since these are isobaric conditions, i. e. the partial pressure of H 2 O vapor is considered constant, the relevant energetic quantity is the free enthalpy.

After impregnation with H 3 PO 4 , the stability of the adsorbed phosphate species at drying conditions depends on a competition for sites between H 2 O and H 3 PO 4 adsorbates.

It seems trivial to describe this by means of adsorption enthalpies, but it brings out a subtle problem: while during impregnation, the adsorption of phosphate species can be described in chemical balance with phosphates in solution, no such reference exists in the dried state, since the solution had evaporated and phosphates or phosphoric acid are not volatile at drying conditions. Hence, for a thermodynamic description of phosphates at the dried interface, adsorption free enthalpies raise questions because no consistent phosphate reservoir such as H 3 PO 4 (aq) exists. Consequently, it is not possible to define an unambiguous adsorption free enthalpy, the only relevant quantities being free enthalpy differences between model adsorption candidates. As a pragmatic choice (i. e. simple reproducibility, knowability from the beginning of the study), the theoretical adsorption free enthalpy is estimated with respect to H 3 PO 4 (g) using the same parameters as for H 2 O (120°C and 2 kPa). An equally valid choice would be any related phosphate-containing reference, for example the most stable adsorption state.

Since the surfaces are fully hydrated, there are two cases of H 3 PO 4 adsorption: physisorption or chemisorption after substitution of chemisorbed water, H 2 O. Moreover, a single H 3 PO 4 can substitute more than one H 2 O molecule, thereby forming higher dentation modes which will be explained in detail in section 3.3. The adsorption reaction for any such reaction of a single H 3 PO 4 can be written in the following simplified way:

H 3 PO 4 + x H 2 O --→ H 3 PO 4 + x H 2 O ↑ (3.1)
where the " " indicates surface adsorbed species and the only variable x indicates the amount of substituted water. Physisorption is described with x = 0. Even without inspecting the reference surfaces in fig.

∆G ads = G[ H 3 PO 4 ] + x G[H 2 O(g)] -G[ ] -G[H 3 PO 4 (g)] (3.3)
The phosphate reservoir G[H 3 PO 4 (g)] is always constant but the surface reference, G[ ] depends on the facet. This is important when comparing two adsorption free enthalpies

[n] ----[R] + H 2 O , ∆G = +29 kJ mol -1 (3.4)
The reconstructed structural reference R(1 1 0) is with a free enthalpy of this reaction of 29 kJ mol -1 less stable than the unreconstructed n(1 1 0). To my best knowledge this parameter is reported by us for the first time, since Wischert discovered the reconstruction in a study of high-temperature models with incomplete hydration, [START_REF] Wischert | γ-Alumina: The Essential and Unexpected Role of Water for the Structure, Stability, and Reactivity of "Defect" Sites[END_REF][START_REF] Wischert | Dinitrogen: a selective probe for tri-coordinate Al "defect" sites on alumina[END_REF][START_REF] Wischert | Visibility of Al Surface Sites of γ-Alumina: A Combined Computational and Experimental Point of View[END_REF] and other studies involving Wischert's R(1 1 0) did not compare it with Digne's n(1 1 0). [START_REF] Réocreux | Structuration and Dynamics of Interfacial Liquid Water at Hydrated γ-Alumina Determined by ab Initio Molecular Simulations: Implications for Nanoparticle Stability[END_REF] By combining the reaction eq. 3.1 and eq. 3.4, one can calculate free enthalpies for adsorption on both structural references with respect to the same surface reference:

∆G [R] ads = G[ H 3 PO 4 ] + (x+1) G[H 2 O(g)] -G[n] -G[H 3 PO 4 (g)] (3.5)
where G[n] is the unreconstructed reference. The reconstruction reaction in eq. 3.4 adds 1 H 2 O to the equation. Overall, these considerations allow for an extensive exploration of the most exposed (1 1 0) facet even including the possibility of reconstruction while avoiding limits in the comparability.

For adsorption of more than one phosphate, as for orthophosphate coadsorption or pyrophosphate adsorption, eq. 3.2 needs to be applied for both phosphate atoms, however, this way we introduce a second phosphate from the arbitrarily chosen reference of H 3 PO 4 in gas phase. For the sake of a good global comparison, average adsorption enthalpies will be given, as the effect of the phosphate reference remains constant this way. For the adsorption of two H 3 PO 4 , leading to a substitution of x and x H 2 O molecules, the average adsorption free enthalpy is defined as:

∆G 2 P ads = 1 2 G[2 H 3 PO 4 ] + (x+x ) G[H 2 O(g)] -G[ ] -G[H 3 PO 4 (g)] (3
∆G P 2 ads = 1 2 G [ H 4 P 2 O 7 ] + (x+x +1) G[H 2 O(g)] -G[ ] -G[H 3 PO 4 (g)] (3.7)
where P 2 is a short form for pyrophosphates. Finally, for coadsorption and pyrophosphate adsorption on R(1 1 0), the free number of H 2 O molecules is again adjusted:

∆G 2 P ads = 1 2 G[2 H 3 PO 4 ] + (x+x +1) G[H 2 O(g)] -G[ ] -G[H 3 PO 4 (g)] (3.8) ∆G P 2 ads = 1 2 G [ H 4 P 2 O 7 ] + (x+x +2) G[H 2 O(g)] -G[ ] -G[H 3 PO 4 (g)] (3.9)
Eq. 3.3 and 3.5-3.9 describe the material balance based on x and x . Structural considerations are described by y and z in eq. 3.2, i. e. the substituted surface groups and proton transfer from orthophosphoric acid to the surface. The latter is a general problem for any adsorption structure that starts from a reference structure with an involved hydrogen bond network. There is in fact no reason to assume that the adsorbate can be constructed without significantly impacting the hydrogen bonds. In principle, it has to be recreated for every adsorption. The next section will discuss how this can be done based on a study of the reference states in fig. 3.1b.

Hydrogen bond network

In order to determine accurate phosphate adsorption free enthalpies, the correct hydrogen bond network needs to be identified. The similarity between µ 2 and µ 3 is surprising, since µ 2 should be more basic than µ 3 .

It seems that, while the local environment of the site matters, it is also important that the protons are regularly dispersed. With H 3 PO 4 adsorption the number of protons in the system changes depending on the number of substituted H 2 O. The capability of the sites to yield or remove protons depends on their proton affinities. Most information for this lies in the reference surfaces, but after the first adsorption tests of phosphate it was found that this adsorbate typically dissociates and transfers one or two protons to the surface (more for higher dentation). All these considerations lead to proton affinities that can be ranked in the following way:

µ 1 -O > > µ 2 -O > P-O ≈ µ 3 -O ≈ µ Al O 1 -OH mobile > µ Al T 1 -OH > > µ 2 -OH, µ 3 -OH, µ 4 -O (3.10)
Here, µ Al O 1 and µ Al T 1 are oxygen atoms at octahedral and tetrahedral aluminium atoms.

For most adsorption candidates, this allowed to limit the range of protonation states to only a 3-4 test calculations which were calculated by DFT. Still, in some cases the number of combinations was more than 30, thus the number of calculations was limited to about 15 cases for which the protons were well dispersed over the surface and layers.

On average, six hydrogen bond networks were tested per considered site, which lead to around 1000 geometry optimizations for the whole phosphate adsorption database. The following section describes which sites were actually considered.

Surfaces and adsorption modes

To understand all tested phosphate surface species in detail, one has to combine all possible adsorption modes with the given surface models. The adsorption modes of chemisorbed molecules or ions can be understood by the number of bonds they form species that only binds to the surface via a single atom is commonly described in a notation like µ x , where x refers to the number of surface atoms it binds to. For more complex adsorption cases, like for phosphates binding via several oxygen atoms to the surface, there is no generally accepted notation available, but the bond patterns are described as dentations. [START_REF] Luschtinetz | Adsorption of phosphonic and ethylphosphonic acid on aluminum oxide surfaces[END_REF][START_REF] Li | Solid-State NMR Spectroscopic Study of Phosphate Sorption Mechanisms on Aluminum (Hydr)oxides[END_REF] In this work, a short notation is introduced based on the Greek letter ν. For one orthophosphate (PO 4 ) unit, x in ν x refers to the number of phosphate oxygen atoms binding to the surface. Notably, the variable x in the reaction eq. 3.2 and x in ν x are equivalent, because by the definition of the reference surfaces, the phosphate adsorption always involves substitution of OH or H 2 O, releasing an equivalent amount of water. Fig. 3.2 shows all the considered adsorption modes from physisorption (ν 0 ) to tridentates (ν 3 ). Except for physisorption, all adsorption modes are constructed by substitution of preadsorbed water sites (cf. red colored surface sites in fig. 3.1), which are mostly µ 1 -OH and µ 1 -H 2 O. For the higher dentation modes (ν 2 and ν 3 ), the actual adsorption requires to identify matching adsorption site patterns. In case of ν 2 , a minimum requirement is two neighboring adsorption sites with a matching distance, whereas ν 3 needs an arrangement of three adsorption sites close to a regular triangle.

Based on the findings of Luschtinetz et al., [START_REF] Luschtinetz | Adsorption of phosphonic and ethylphosphonic acid on aluminum oxide surfaces[END_REF] the distance criterion was held on relatively analyzing the free enthalpies of all candidates on (1 0 0) and n(1 1 0), the following criteria were derived:

(1 0 0) ν a 1 , ν b 1 , ν c 1 ν ab 2 , ν aa 2 * , ν bc 2 , ν cb 2 * ν abc 3 , ν acb 3 * , ν aa B 3 * , ν aa C 3 * , ν abA 3 n(1 1 0) ν a 1 , ν b 1 , ν d 1 , ν e 1 ν ab 2 , ν bb 2 , ν de 2 , ν ee 2 ν abA 3 , ν abD 3 , ν bb B 3 , ν bb E 3 , ν deA 3 , ν deD 3 , ν ee B 3 , ν ee E 3 R(1 1 0) ν a 1 , ν b 1 , ν c 1 , ν d 1 , ν e 1 ν ab 2 , ν ac 2 , ν bc 2 , ν de 2 , ν abA 3 , ν abD 3 , ν acC 3 , ν acF 3 , ν bcB 3 , ν dC 2 , ν dF 2 , ν eB 2 , ν eE 2 ν bcE 3 , ν deD 3 , ν deA 3 , ν cBC 3 , ν cEF 3 * sites
1. Two phosphate oxygen atoms never bind to the same aluminium atom.

2. Substitution of µ 1 -O is generally preferred over µ 2 -O and µ 3 -O.

3. Free phosphate oxygen atoms are either protonated or stabilized by hydrogen bonds.

The first criterion is consistent with one of the findings in ref. [52]. It simplified the structural exploration on R(1 1 0), where several of the candidates in tab. 3.1 were discarded. However, it should be noted that some of the unfavorable adsorption structures lead to a critical further insight. For example, n(1 1 0)ν ee B 3 leads to the reconstruction of R(1 1 0)ν eB 3 , which is one of the most stable adsorption modes discovered; it in fact showed that it was necessary to include models based on Wischert's reconstruction.

(1 0 0)ν abA 3 , while clearly violating the first criterion, is only 10 kJ mol -1 less stable than

(1 0 0)ν abc 3 . After all, while there are some exceptions, many adsorption candidates can be excluded only based on the first criterion. As shown in tab. 3.1, only eight ν 3 candidates remain out of originally 23. Overall it reduces the number of relevant sites to 30 from originally 51 candidates. The second criterion excludes monodentates on alumina oxygen sites, which in most cases were not considered due to steric limitations, but it would also exclude two water adsorption sites, ν a 1 on n(1 1 0) and R(1 1 0). As a consequence of the third criterion, the phosphate protonation state for most cases was HPO 2 - 4 , as at least one oxygen protrudes from the surface and hence cannot be stabilized in a hydrogen bond.

Based on the results for adsorption of a single orthophosphate on the three γ-Al 2 O 3 model facets, two additional cases were studied: coadsorption of two orthophosphates and pyrophosphate adsorption. This was critical since it is well established by [START_REF] Luis | Hydrolysis rates of inorganic polyphosphates in aqueous solution as well as in soils and effects on P availability[END_REF] P NMR that the phosphate species change with coverage. [START_REF] Decanio | FT-IR and Solid-State NMR Investigation of Phosphorus Promoted Hydrotreating Catalyst Precursors[END_REF][START_REF] Van Eck | A Solid State Double Resonance NMR Investigation of Phosphorus-Impregnated γ-Al 2 O 3[END_REF] While there is no tight control over coverage for the computational model, the actual modelled coverage can be calculated based on the surface unit cell. For one phosphate on the (1 0 0) facet, the coverage is 2.2 P nm -2 , on (1 1 0) it is 1.5 P nm -2 . These values double when adding a second phosphate to the surface model (4.3 and 3.0 P nm -2 ). The ν x notation for coadsorption is concatenated in case of coadsorption like in ν x +ν x . The new variable x also leads to water desorption as given in reaction eq. 3.2, depending on the adsorption mode ν x of the additional phosphate. For polyphosphates like pyrophosphate, the ν strict limitation (two phosphate oxygen atoms never bind to the same aluminium atom).

While -conservatively estimated -up to 203 structures could have been included (i. e. by application of this first criterion), the other three criteria also allow for a more selective approach, within which all plausible sites and site combinations are included in the data set.

The most stable adsorption modes

As pointed out in sections 3.2 and 3.3 with considerations for creating hydrogen bond networks and with site selection criteria, there are several general rules which can be used to identify a manageable number of candidate structures. Nevertheless, the most stable results comprise exceptional structural features that will be discussed here. It should also be noted that the reported absolute free enthalpy values heavily depend on the environmental conditions, whereas enthalpies are less dependent on temperature. Tab. Only the free enthalpy of the most favorable system at a given degree of water abstraction was included.

water is released as no substitution occurred but on R(1 1 0)ν 3 , three H 2 O are released for the ν 3 adsorption and one for surface reconstruction (cf. eq. 3.5). The enthalpy difference between n(1 1 0)ν 0 and R(1 1 0)ν 3 is -188 kJ mol -1 but the free enthalpy 112 kJ mol -1 , which corresponds to the difference of the entropy contribution T •∆S. In this example it is T •∆S = 300 kJ mol -1 for a vapor entropy gain after releasing four H 2 O. When

(1 0 0)ν 0 is compared with R(1 1 0)ν 3 , the result is similar, with T •∆S = 315 kJ mol -1 .

While identifying the most stable structures, at drying conditions they are ν 2 and ν 3 on R(1 1 0) -indicating that entropy gain by water release is the most decisive term in the free enthalpy results. ∆H and ∆S also slightly depend on the temperature, but the most favorable adsorption sites still remain the same. This balance of enthalpy ∆H and entropy ∆S is linear dependent on the temperature due to T •∆S, which means that the order by free enthalpy of the tabulated structures in tab. 3.2 is sensitive to the temperature. For example at ambient conditions, i. e. the free enthalpy of R(1 1 0)ν 3 is 118 kJ mol -1 and 152 kJ mol -1 for (1 0 0)ν 0 (cf. tab. A. Energy comparisons between the facets (1 0 0) and (1 1 0) are not ideal for two reasons. In the following sections the coverage on both facets will be considered as equivalent.

The coverage effect is qualitatively separated in high and low model coverage domains.

Orthophosphate coadsorption and pyrophosphates effectively represent the same coverage, thus the coadsorption is the preferred reference to quantify the stability of pyrophosphates relative to orthophosphates.

Single orthophosphates

Since the order of most stable structures can vary with environmental conditions, the focus on most stable few structures can become arbitrary. In fig. 3.4 three adsorption species from each reference surface were selected based on their exceptional structural features. They are also within the most stable candidates, but at different temperatures.

Physisorption is the most stable adsorption mode on (1 0 0) at ambient conditions (see fig. It is a remarkable fact that the most stable adsorption mode is found on the reconstructed surface reference, since the hydrated structural reference of n(1 1 0) is 29 kJ mol -1 more stable than R(1 1 0) before phosphate adsorption. This stabilization is clearly triggered by the adsorption, but the question is, how exactly. Tab. 3.3 lists all stable bidentate and tridentate adsorption modes ordered by their location on the surface on n(1 1 0) and R(1 1 0). In most of these cases, i. e. ν bc 2 , ν de 2 , ν bcE surface model that can occur at low water coverage due to high temperature, [START_REF] Wischert | γ-Alumina: The Essential and Unexpected Role of Water for the Structure, Stability, and Reactivity of "Defect" Sites[END_REF] or in solution, provided that it can be stabilized with a suitable adsorbate at alumina oxygen site B. In general, it underlines the complexity of these adsorption phenomena, since adsorbate-induced surface reorganization like this one is still very difficult to predict.

Coadsorption and pyrophosphate adsorption

Coadsorption and pyrophosphate enthalpies are determined as described before with eq. 3.6-3.9, where for better comparison with single orthophosphate adsorption, average phosphate adsorption enthalpies are reported. The two best coadsorption and pyrophosphate cases are tabulated in tab. 3.4. Compared to the adsorption of a single orthophosphate on the reference surfaces, coadsorption and pyrophosphate adsorption reach a limit in terms of water abstraction. For single orthophosphate adsorptions, the free enthalpy was dominated by the entropy term, with R(1 1 0)ν 3 (4 H 2 O) being the most stable case at 120°C. In comparison, the coadsorption case with most water desorption is R(1 1 0)ν 3 +ν 3 (7 H 2 O), which is with ∆G = -113 kJ mol -1 less stable than the most favorable coadsorption cases on n(1 1 0) and R(1 1 0). Removal of 4 or 5 H 2 O seems to be the limit for phosphate adsorption at these conditions. A similar trend is observed on the (1 0 0) facet. Regarding adsorption site combinations, it is clear that the second adsorbing phosphate is limited to the remaining sites. One would generally expect that the most favorable sites are occupied first followed by the second most favorable site. However, the sites in tab. A similar overall trend can be observed for pyrophosphates. The available range of sites are arranged in such a way that they can bind a pyrophosphate on four sites, while still staying with the selection criteria discussed in section 3.3.

Summary and perspectives

The result of this structural exploration is a database of overall more than 1000 DFT calculations based on 159 adsorption site candidates. To ensure that the most stable cases were identified, the exploration was done on multiple levels of complexity. With a defined set of structural references, candidate sites had to be defined. This first involved labelling of the relevant surface oxygen sites given in fig. 3.1a, followed by defining a set of relevant adsorption modes based on oxygen site combinations. These adsorption modes were recombined to obtain candidate structures for orthophosphate coadsorption and pyrophosphates. All structures identified in this process required a test for several relevant hydrogen bond networks by applying combinatorics described in section 3.2. Several general trends and features were identified just by comparing their thermodynamic properties. One of the important findings is that the presence of phosphates leads to surface reorganization, here observed by the stabilization of Wischert's reconstruction. It was proposed in several experimental studies, [START_REF] Decanio | FT-IR and Solid-State NMR Investigation of Phosphorus Promoted Hydrotreating Catalyst Precursors[END_REF][START_REF] Van Eck | A Solid State Double Resonance NMR Investigation of Phosphorus-Impregnated γ-Al 2 O 3[END_REF] but this work is the first successful attempt to reconstruct these effects in computational atomistic models.

The database collected as described in this chapter was the basis for further studies.

As an attempt to verify the result database, GIPAW chemical shift simulations were conducted for the most stable structures. The results and performance in comparison to the experiment are discussed in chapter 4. With the ambition to better understand the structure-NMR relationship, the chemical shift simulations of this database were used to train a machine learning model. The outcome of this is given in chapter 6.

Finally, in order to cover phosphate species at impregnation conditions, the best database structures were used as target geometries in state-of-the-art metadynamics simulations (chapter 5). In particular, high dentation patterns were probed in order to see whether they can form by means of a phosphorylation. Moreover, based on thermodynamic estimates, the condensation of pyrophosphates was excluded within the tested models.

This has important consequences, as experimental evidence demonstrates their presence.

Pyrophosphate and other polyphosphates either come from the solution, they condensate faster than the most stable adsorption modes can form or they are a consequence of inhomogeneous distribution in the material. These three alternatives are discussed in detail in section 4.3 (p. 4.3), after a detailed analysis of the prediction performance of these models for 31 P NMR spectroscopy.

Combining NMR and DFT

Experimental and theoretical work have always been complementary. The discussed [START_REF] Luis | Hydrolysis rates of inorganic polyphosphates in aqueous solution as well as in soils and effects on P availability[END_REF] P and 27 Al NMR experiments involve -even without computational models -a considerable theoretical foundation in order to infer chemically relevant information from the raw spectral data. Many regular single component systems can be well characterized with different experimental methods, allowing accurate experimentally derived models. In contrast, disordered multicomponent mixtures are usually too complex to permit a complete atomic scale description that is exclusively based on experimental insight. In the past, these knowledge gaps had been commonly filled with more or less speculative hypotheses in lack of better alternatives. Computational modelling is different from the theory needed to interpret experiments, since it is carried out separately. As such, computational results can be evaluated independently, but since they do not involve any measurements of a real chemical system, their relevance is sometimes disputed.

Combining NMR cross-polarization experimental studies with computational modelling can therefore improve the discussion at both sides.

In the structural exploration (chapter 3), we presented computational results, identified several trends within the most stable adsorption modes and sites and rationalized them by the balance between enthalpy and H 2 O vapor entropy. As described there, the exploration was comprehensive and it is therefore unlikely that the most representative structures were overlooked -but that only includes the three model surfaces tested. From a computational perspective, it is very difficult to proof that the explored phosphate adsorption geometries are relevant models to describe the real systems, thus the structures have to be validated against the experimental results. On the other hand, the unambiguous discussion of the experimental findings requires a lot of results to be reconciled, where drawn figures of adsorption patterns are often speculative. In computational models, however, these presumptions can be evaluated and compared quantitatively; essentially by composing structures in silico, and benchmarking it by computation of spectroscopic data from the model geometry. For this work, it is particularly beneficial to compare calculated chemical shifts of different adsorption modes and sites and compare them with the experiment. This is also considerably different to some studies with a strong focus on experimental results and only few computational models that are produced to illustrate the conclusions.

The goal in this chapter is to determine the extent to which the model geometries represent the real phosphate speciation on γ-Al 2 O 3 . It requires to identify which computational models reproduce the experimental 31 P NMR signal band of dried H 3 PO 4 -impregnated γ-Al 2 O 3 . This also requires a computational NMR reference, because -equivalently to NMR measurements -a reference substance needs to be computed in order to compare the computational findings with literature data. Since simulations of 85 % phosphoric acid, the 31 P NMR standard reference, are not feasible, aluminium phosphate minerals with known crystal structures and NMR iso-shifts were used instead. While a single one could be used as a reference, three compounds -berlinite, variscite and wavellitewere used in order to minimize systematic errors. The reference is explained in detail in appendix A.3. In order to compare the most stable candidates discussed in section 3.4 and 3.4.2, correlation graphs were produced where the chemical shifts from GIPAW calculations are plotted against the phosphate adsorption free enthalpy.

Experimental results

As discussed in section 1.1.5 (p. 23), DeCanio et al. showed with [START_REF] Luis | Hydrolysis rates of inorganic polyphosphates in aqueous solution as well as in soils and effects on P availability[END_REF] P NMR, that two aluminium phosphate phases can form when preparing samples with 4 to 10 wt% H 3 PO 4

Experimental results

into γ-Al 2 O 3 , which corresponds to coverages of approximately 3.6 to 9 P nm -2 . They rationalized the initial signal shift from 1 to 4 wt% (≈0.9-3.6 P nm -2 ) with polyphosphate formation. [START_REF] Decanio | FT-IR and Solid-State NMR Investigation of Phosphorus Promoted Hydrotreating Catalyst Precursors[END_REF] Moreover, DeCanio et al. observed a sudden drop in γ-Al 2 O 3 pore volume between phosphate loadings of 6 and 10 wt% (≈5.4-9 P nm -2 ), which is one reason why our experimental collaborators focussed only on loadings up to 4.1 P nm -2 . Some experimental work on similar systems also attempted to resolve the adsorption mode of phosphates. Li et al. studied the adsorption of phosphates on boehmite using 31 P REAPDOR CP MAS NMR experiments. They identified bidentation (ν 2 ) as a more favorable adsorption mode than monodentation (ν 1 ) of orthophosphates on boehmite on dried and wet samples. [START_REF] Li | Solid-State NMR Spectroscopic Study of Phosphate Sorption Mechanisms on Aluminum (Hydr)oxides[END_REF][START_REF] Li | Surface Speciation of Phosphate on Boehmite (γ-AlOOH) Determined from NMR Spectroscopy[END_REF] While this is an important finding and a comparable behaviour is expected on γ-Al 2 O 3 , the authors did not include other adsorption modes like tridentates (ν 3 ) or physisorption (ν 0 ). Overall, these studies set a range of species and adsorption modes to be identified, and also a "confined" range of requirements that the computational models need to satisfy. The focus here is on phosphate loadings at which DeCanio et al. proposed that condensed phosphates would influence the 31 P NMR signal the most. Lesage on γ-Al 2 O 3 samples provided by IFPEN with the ROAD4CAT project. [START_REF] Hühn | Structural Characterization of Phosphate Species Adsorbed on γ-Alumina by Combining DNP Surface Enhanced NMR Spectroscopy and DFT Calculations[END_REF] does not concern all phosphate centers. The comparison between through-bond and through-space [START_REF] Luis | Hydrolysis rates of inorganic polyphosphates in aqueous solution as well as in soils and effects on P availability[END_REF] by experiments. However, currently there is no alternative to computational methods and simulations in order produce a coherent model which can be benchmarked against the experimental observations. Moreover, future experiments on even more complex systems could rely on such models as foundation.

Computational chemical shifts

In order to completely characterize phosphate adsorption at the γ-Al 2 O 3 surface, one should take into account the following aspects:

• Formation of aluminium phosphate phases [START_REF] Decanio | FT-IR and Solid-State NMR Investigation of Phosphorus Promoted Hydrotreating Catalyst Precursors[END_REF][START_REF] Van Eck | A Solid State Double Resonance NMR Investigation of Phosphorus-Impregnated γ-Al 2 O 3[END_REF] • Polyphosphate formation [START_REF] Decanio | FT-IR and Solid-State NMR Investigation of Phosphorus Promoted Hydrotreating Catalyst Precursors[END_REF] • Coadsorption or coverage effect (interaction with other phosphates in vicinity)

• Adsorption mode effect (i. e. dentation, for boehmite discussed in ref. [123])

• Adsorption site effect

• Stability (discussed for calcination in ref. [12])

Several of these have been studied before, but it was not yet possible to include these aspects in a coherent model. The coming chapters demonstrate how DFT calculations can enhance our understanding by comparing these results with model chemical shifts. The GIPAW method [START_REF] Pickard | All-electron magnetic response with pseudopotentials: NMR chemical shifts[END_REF][START_REF] Yates | Calculation of NMR chemical shifts for extended systems using ultrasoft pseudopotentials[END_REF] employed in this work explicitly calculates the shielding tensor from the electronic structure for every nucleus of the system. Nevertheless, in the following sections, only the chemical shift of 31 P nuclei will be considered.

Chemical shift of orthophosphate

For a single adsorbed orthophosphate species, the coverage of each model facet is determined by the surface area of the unit cell, that is 0.46 nm 2 for the (1 0 0) facet and 0.66 nm 2 for the (1 1 0) facet. This corresponds to local coverages of 2.2 and 1.5 P nm -2 .

These local coverage values cannot be related to experimental phosphate loadings, as the phosphates are not expected to be equally distributed over the different facets. When comparing the results with experimental data, one also needs to bear in mind that theoretical models are built at the pH corresponding to the point of zero charge. For the computational models, it just means that the calculated systems are uncharged, but for the experimental part it is required to ensure the absence of ions except for the phosphates during sample preparation. GIPAW calculations were conducted for the most stable sites for each adsorption mode -ν 0 , ν 1 , ν 2 and ν 3 -and structural reference, (1 0 0), n(1 1 0), R(1 1 0). The stability of the modelled phosphate species was already discussed in detail in chapter 3. In this chapter, only the most stable systems are considered for each adsorption mode. In fig. 4.2, the simulated 31 P chemical shifts are plotted against the phosphate adsorption free enthalpies defined in section 3.1 (p. 52).

The adsorption on the (1 0 0) facet is less stable than on (1 1 0), thus more phosphates are expected to bind to the latter. The signals of phosphates on (1 0 0) are between 12 and 4 ppm, which is not within the experimental signal range that starts from 0 ppm. Thus, the phosphate does not seem to adsorb on the (1 0 0) at low coverage. For the (1 1 0) facet, the most stable adsorption modes describe the upfield half of the experimental

Orthophosphate coadsorption

As explained in section 3.1 (p. 52), the coadsorption free enthalpy is the average free enthalpy of adsorption for both phosphates. This is the only value that can be directly compared with single orthophosphate adsorption. Nevertheless, each coadsorption system produces two results for the GIPAW 31 P chemical shift simulations. Fig. 4.3a presents the around -5 ppm albeit all of them are coadsorbed with another predicted signal that is outside the experiment range. Even though it should in general be expected to find phosphates on the (1 0 0) facet, as soon as the most stable adsorption sites on (1 1 0) facet are occupied, these monodentation (ν 1 ) signals could as well be characteristic to any of the monodentates. Monodentates on all facets produce calculated 31 P chemical shifts between 6 and -6 ppm. This also concludes our findings for the influence of the adsorption mode. The predominant species, bidentate (ν 2 ) and tridentate (ν 3 ), can contribute to an upfield shift of the 31 P NMR signal, whereas lower dentation modes (ν 0 , ν 1 ) should have a minor influence to the NMR signal. Another aspect is the decreased free enthalpy difference between adsorption based on the n(1 1 0) reference structure and R(1 1 0). The chemical shifts on n(1 1 0) are clearly at more negative values than those on R(1 1 0). However these differences are not easy to disentangle with the general view on the best adsorption mode combinations, but require a closer look at individual adsorption sites. This is given in fig. 4.4 and will be discussed in the next section.

Coverage and site effect

The previous sections validated the models by comparison with the chemical shift calculations and furthermore analyzed the effect of the predominant ν 2 and ν 3 adsorption modes. In this section, the chemical shifts of the six most stable adsorption candidates at varying coverages will be compared. The same sites are connected by lines. At the highest coverage, the site of the coadsorbing phosphate cannot be chosen unambiguously, so several cases were plotted.

The highest coverages (coadsorption) depend on the choice of the other coadsorbing phosphate. Since there is no simple criterion to pick this other phosphate, all adsorption mode combinations were plotted. Fig. 4.4a shows the effect of coverage while keeping the site constant and fig. 4.4b provides structures of the tested sites. Note that this plot has a focus on sites which means that all connected points not only share the same structural reference and the adsorption mode, but they are bound to the model surface at the same position and orientation. They are different in where and how their surrounding rearranges depending on the presence of neighboring phosphates. It seems that there is a small trend of decreasing chemical shifts at higher coverage, which could be called "direct" coverage effect. The only exception is the monodentate (1 0 0)ν a 1 , which shows a strong deshielding for the highest coverage. As this observation is unique to monodentates, which are much less stable based on the thermodynamic results (cf. fig.

Adsorbed pyrophosphates

Pyrophosphates are the only condensed phosphates discussed in this work, due to the excessively large phase space that would need to be explored otherwise for additional polyphosphates. They are evaluated in a similar way as the coadsorption cases, where also each system produces two chemical shifts and one free enthalpy.

The two most stable pyrophosphates, R(1 1 0)ν bcB,e orthophosphate species on (1 0 0). The thermodynamic prediction on the model surfaces disfavor the presence of pyrophosphate, which is a weakness of the model at this stage.

This is likely to improve by exploring higher coverages (4.5-6.5 P nm -2 ), i. e. a comparison of three orthophosphate coadsorption ν x +ν x +ν x vs. pyrophosphate/orthophosphate coadsorption ν xx +ν x vs. triphosphates ν xx x . Since the most stable adsorption sites are generally occupied with two phosphates, and physisorption is in most cases much less stable, pyrophosphates might be more favorable at this stage. Overall the pyrophosphate adsorption cannot yet be conclusively explained, but different rationales for the presence of pyrophosphates are presented in the discussion in section 4.3.

Terminal pyrophosphates or physisorbed orthophosphates

Summary and perspectives

The computed chemical shifts of the most stable model adsorption candidates agree well with the experiments. In particular the most stable adsorption modes -bidentates (ν 2 ) and tridentates (ν 3 ) on (1 1 0) -consistently cover experimental range of [START_REF] Luis | Hydrolysis rates of inorganic polyphosphates in aqueous solution as well as in soils and effects on P availability[END_REF] The thermodynamic prediction for phosphate adsorption does not permit the presence of pyrophosphates, because coadsorbed orthophosphates are overall more stable. This result is not in alignment with the experimental findings, as pyrophosphate are clearly observed. There are three rationales for this. One possibility is that pyrophosphates come from the solution and adsorbed before being hydrolyzed. This is probable if they were part of the impregnation solution due to their kinetic stability. [START_REF] Pedley | The Hydrolysis of the Condensed Phosphates: I. Sodium Pyrophosphate and Sodium Triphosphate[END_REF][START_REF] Stockbridge | Enhancement of the Rate of Pyrophosphate Hydrolysis by Nonenzymatic Catalysts and by Inorganic Pyrophosphatase*[END_REF] It is not expected that they form in aqueous solution, due to the positive reaction enthalpy; [START_REF] Stockbridge | Enhancement of the Rate of Pyrophosphate Hydrolysis by Nonenzymatic Catalysts and by Inorganic Pyrophosphatase*[END_REF] significant phosphate condensation only occurs at high phosphate concentrations with trace amounts of water. [START_REF] Thilo | Zur Chemie der kondensierten Phosphate und Arsenate. XVII. Der Verlauf und die Produkte der Entwässerung der Monophosphorsäure H 3 PO 4[END_REF] This possibility can be tested by reevaluation of the synthesis.

Another rational could be found in a more complex relationship between local coverages and phosphate loading. If phosphate chemisorption is rapid and the adsorbates at the solid-liquid interface are too stable, then the adsorption process is exclusively under diffusion control. In turn, phosphates are immobile and remain at the sites they bound to first would lead to an inhomogeneous surface distribution of phosphates, between outermost areas and pores. The same trend is expected if the most stable adsorption sites are concentrated at relatively small facets, which in turn produce high local phosphate coverages there. Since this initial rapid chemisorption would only occur on the reactive surface sites, this effect should be more stronger in case of low phosphate loadings since the remaining free phosphates would adsorb elsewhere after the most stable sites are occupied. Testing this hypothesis would involve additional structure/free enthalpy explorations at higher coverage and likely needs to be extended by different surface models. Atomic scale and mesoscale models for phosphate diffusion in the pores might be needed as well. As last alternative, pyrophosphates may also form at the surface at lower coverage, if they are kinetically more stabilized than coadsorption. However, this also requires that the reaction barrier of the most stable orthophosphate adsorption modes is too high for these cases to occur, which would exclude some of the most stable model adsorption structure based on kinetics.

In general the GIPAW calculations for [START_REF] Luis | Hydrolysis rates of inorganic polyphosphates in aqueous solution as well as in soils and effects on P availability[END_REF] P NMR demonstrate that different structural aspects as well as facet distribution have to be considered simultaneously. The computational models provide a considerably expanded foundation to discuss these effects.

However, we are still at the beginning of understanding all aspects of phosphate surface speciation. Reaction kinetics need to be studied in order improve the comparison between orthophosphates and pyrophosphates. This is much more complicated and computationally demanding, but as a first step in this direction, the kinetics of one adsorption reaction H 2 O/γ-Al 2 O 3 interface and several phosphorylation reactions (cf. section 1.1.3, p. 18) will be studied in chapter 5.

Towards alumina impregnation

After exploration of phosphate adsorption on dried γ-Al 2 O 3 , one can define a set of target geometries and ask how they formed through wet incipient impregnation. The most stable adsorption modes in dry samples, bidentates (ν 2 ) and tridentates (ν 3 ), may not be the most stable modes during impregnation. [START_REF] Roy | Phosphate Adsorption on γ-Alumina: A Surface Complex Model Based on Surface Character-ization and Zeta Potential Measurements[END_REF] However, experiments showed that once impregnated, the phosphates cannot be removed easily without specific leaching treatment. This chapter provides insight into stabilities of the explored adsorption modes in chapter 3 and demonstrates how well-tempered metadynamics (WTmetaD) can be used to explore reaction paths producing some of the described dentations. The most favorable adsorption sites under drying conditions were in focus. In particular, the presence of adsorbed pyrophosphates is evident and they likely have a key role to rationalize the coverage-dependent 31 P NMR shift. Here, we want to take a closer look at how the initial adsorption could take place and how pyrophosphate can be involved in some reaction steps. The goal is in particular to find out whether the higher dentation states, i. e.

bidentates and tridentates that involve alumina oxygen atoms, are kinetically possible and if they are still stable compared to the other adsorption modes. This chapter also focussed on the (1 1 0) facet, since chapter 4 showed that the (1 0 0) facet is less active for H 3 PO 4 adsorption.

The adsorption of orthophosphates can in general be thought of a reaction sequence from physisorption to the highest possible dentation mode:

H 3 PO 4 (aq) --→ ν 0 -H 2 O
--------

+ H 2 O ν 1 -H 2 O -------- + H 2 O ν 2 -H 2 O
--------

+ H 2 O ν 3 (5.1)
Starting from phosphoric acid in solution, the adsorption species are formed in the sequence: physisorption (ν 0 ), monodentate (ν 1 ), bidentate (ν 2 ), and tridentate (ν 3 ). The first reaction of physisorption of H 3 PO 4 (aq) from bulk water was not tested here. On All WTmetaD calculations were done at a simulation temperature of 330 K with a bias factor γ of 72.9 and gaussian depositions with a height of 3 kJ mol -1 every 10 fs if not noted otherwise. Furthermore parallelization using four walkers, with two starting at the initial state and two at the final state of the reaction.

Phosphate condensation in bulk water

The condensation and hydrolysis of pyrophosphate is the most fundamental phosphorylation reaction (cf. section 1.1.3, p. 18) because it only involves the two smallest polyphosphates, orthophosphate and pyrophosphate. At 120°C, 10 % of water-free orthophosphoric acid is converted to pyrophosphate, and trace amounts of longer polyphosphate chains can be measured. [START_REF] Thilo | Zur Chemie der kondensierten Phosphate und Arsenate. XVII. Der Verlauf und die Produkte der Entwässerung der Monophosphorsäure H 3 PO 4[END_REF] However, pyrophosphate is not thermodynamically stable in aqueous solutions, thus they do not form such solutions, but existing phosphates hydrolyze at very slow rates (more than 21 h half-life at 65°C, [START_REF] Pedley | The Hydrolysis of the Condensed Phosphates: I. Sodium Pyrophosphate and Sodium Triphosphate[END_REF] or several days at room temperature). [START_REF] Luis | Hydrolysis rates of inorganic polyphosphates in aqueous solution as well as in soils and effects on P availability[END_REF] In particular the cyclic polymetaphosphates can remain in solution for weeks, before significant decomposition. [START_REF] Luis | Hydrolysis rates of inorganic polyphosphates in aqueous solution as well as in soils and effects on P availability[END_REF] In this section, the pyrophosphate hydrolysis reaction is discussed in aqueous solution. It is used as a testbed of the method since the barrier and reaction energy are known by experiments for the hydrolysis. [START_REF] Stockbridge | Enhancement of the Rate of Pyrophosphate Hydrolysis by Nonenzymatic Catalysts and by Inorganic Pyrophosphatase*[END_REF] Then in the next sections, these findings are compared to phosphate speciation reactions via phosphorylation on

γ-Al 2 O 3 .
Condensation and hydrolysis in phosphorylation reactions are quite complex as they involve two P -O bond formations and may also include several proton shifts, involving the solvation shell. [START_REF] Shina | Why nature really chose phosphate[END_REF][START_REF] Glaves | Mechanistic Insights into the Hydrolysis of a Nucleoside Triphosphate Model in Neutral and Acidic Solution[END_REF][START_REF] Prasad | Addressing Open Questions about Phosphate Hydrolysis Pathways by Careful Free Energy Mapping[END_REF] Several theoretical studies were done using AIMD and metadynamics, but typically within the context of biochemistry, involving nucleoside polyphosphates in presence of magnesium and within enzymes. [START_REF] Tripathi | The GTPase hGBP1 converts GTP to GMP in two steps via proton shuttle mechanisms[END_REF][START_REF] Prasad | Addressing Open Questions about Phosphate Hydrolysis Pathways by Careful Free Energy Mapping[END_REF] To our best knowledge, the reaction of the pure compounds (the phosphates in water) was not yet simulated by AIMD/WTmetaD. For brevity, the following discussion only refers to the condensation reaction and neglects the hydrolysis since it is equivalent. The system consists of 124 H 2 O and two H 3 PO 4 which reflects a solution of 1 M phosphoric acid with pH=1-2. It is generally complicated to identify good collective variables (CV). Since the system is symmetric, starting from two orthophosphates, one would ideally find CV that respect this symmetry. However, the symmetry is broken in the course of the reaction, in which only

(0, 8) H 2 PO - 4 (1, 8) [H 2 O 3 PO -PO 4 H 2 ] 2 - (0, 7) PO - 3
(1, 7) one phosphate provides the bridging oxygen atom. For a two-dimensional well-tempered metadynamics simulation, that was discussed in section 1. 1.3 (p. 18), the following two collective variables were identified: (i) a variable which counts the number of phosphorusbridging oxygen atoms keeps track of the pyrophosphate bond formation, CN(P -O bridge ), and (ii) the number of phosphate-binding oxygen atoms, CN(P -O all ), counts the total number of oxygen atom binding to phosphorus. The CN(P -O bridge ) should vary from 0 (orthophosphate) to 1 (pyrophosphate) and CN(P -O all ) from 8 (2 H 3 PO 4 ) to 7 (H 4 P 2 O 7 ).

H 2 P 2 O 2 - 7 + H 2 PO - 4 + H 2 PO - 4 + H 2 P O - 4 -O H - , -H + -OH -, -H + -OH -, -H +
Both variables are explained in detail in appendix A. 4, p. 162. This means that in the basis of these CV, one can find the orthophosphates at (0, 8) and the pyrophosphate at (1, 7). The scheme of a phosphorylation reaction (cf. sec- These values were chosen after careful analysis of (unbiased) AIMD simulations such that those are not affected by the walls.

As the FES in fig. 5.2c shows, the reaction is activated by initial elimination of water from one phosphate leading to a short-lived metaphosphate ion PO - 3 :

2 H 2 PO - 4 ----PO - 3 + H 2 PO - 4 + H 2 O ----H 2 P 2 O 2- 7 + H 2 O (5.2)
With a free energy of 95 kJ mol -1 , the PO - 3 (positioned at (0.1, 6.7) in fig. 5.2c) is far more likely to react back to orthophosphate (F ‡ = 9 kJ mol -1 ) than to pyrophosphate (F ‡ = 33 kJ mol -1 ). In fact it cannot be excluded that the small barrier of F ‡ = 9 kJ mol -1 vanishes at longer simulation time. The activation by initial addition of water to orthophosphate was also covered at (1.0, 7.7) but with a higher free energy of 119 kJ mol -1 . Overall, H 2 PO - 4 is 27 kJ mol -1 more stable than the H 2 P 2 O 2 - 7 , which is consistent with experimental findings of 32 kJ mol -1 . [START_REF] George | Squiggle-H2O". An enquiry into the importance of solvation effects in phosphate ester and anhydride reactions[END_REF] The total barrier for condensation is F ‡ = 127 kJ mol -1 , and 100 kJ mol -1 for the hydrolysis. The hydrolysis barrier is 19 kJ mol -1 lower than the result of G ‡ = 119 kJ mol -1 reported by Stockbridge and

Wolfenden. [START_REF] Stockbridge | Enhancement of the Rate of Pyrophosphate Hydrolysis by Nonenzymatic Catalysts and by Inorganic Pyrophosphatase*[END_REF] This is still an excellent agreement, since it is well known, that PBE, the used DFT functional, tends to underestimate barriers and ab initio WTmetaD simulation times are difficult to converge.

Besides P -O bond reformations, the reaction also consists of proton transfer reactions in order to form the activated metaphosphate ion PO - 3 . The protonation states cannot be retrieved from the FES directly, but by evaluation of the biased trajectory. Fig. 5.2d is a scatter plot where each dot refers to the protonation of the less protonated phosphate center. This plot clearly shows the most common protonation states in each area of the FES: The orthophosphoric acid is clearly deprotonated once to H 2 PO - 4 and the metaphosphate is PO - 3 ; the pyrophosphate is an interesting case, since the protonation varies between 0 and 1, which is mostly an intramolecular protonation:

[HO 3 P-O-PO 3 H] 2-----[O 3 P-O-PO 3 H 2 ] 2- (5.
3)

The right-side configuration forms after addition of H 2 PO - 4 to PO - 3 , but it is surprisingly common in the simulation. In fig. 5.2d this is visible as overlapping dark green and purple points reflecting protonations of 0 an 1 in one of the phosphate groups The overall protonation of the pyrophosphate ion can be identified in a different scatter plot (not shown) that estimates the protonation of both phosphate centers, which is in nearly all points 2 and validates eq. 5.3. This is also the case for the transition state at the highest barrier around (0.4, 6.7), which validates eq. 5. unexpected ways which are difficult to explore. [START_REF] Prasad | Addressing Open Questions about Phosphate Hydrolysis Pathways by Careful Free Energy Mapping[END_REF]

Orthophosphate adsorption in presence of water

The γ-Al 2 O 3 (1 1 0) model facet comprises dentation sites for two sites to form bidentates.

Following the labels introduced in fig. 3.1a on page 53, they can be located as neighboring oxygen atoms at sites de and bb' (or bc on R(1 1 0)). The de sites bind to tetrahedral and octahedral aluminium atoms (labelled as α and β in fig. 3.1ii and iii, p. 53), whereas bb' are sites from two symmetry-equivalent octahedral aluminium atoms (labelled with γ and γ'). The de adsorption reaction is much more complicated, due to the tetrahedral aluminium binding site α, potentially involving the adsorption of one or two additional H 2 O molecules in an intermediate octahedral complex around the aluminium atom α, similar to those shown in [96]. With these aspects, the formation of the ν de 2 bidentate would require a difficult and time-consuming simulation. The bb' site, on the other hand, is limited to a much simpler elimination-addition process as it only involves two symmetry-equivalent sites. In this chapter, we focus on the adsorption on the bb' site.

The relevant collective variables for the reaction are CN w , the coordination number of Knowing this, CN w and CN p can be combined to a single CV:

CV = CN w -CN p (5.4)
For this CV, bidentate, monodentate and physisorption occur at well distinguished values: Overall, the one-dimensional WTmetaD successfully sampled the adsorption sequence from physisorption to bidentate.

ν bb 2 -→ CV = -2 , ν b 1 -→ CV = 0 and ν 0 -→ CV = 2
Compared to these barriers, a large free energy landscape of more than 300 kJ mol -1

was carved out throughout the simulation, which is visible at the ordinate of fig. 5.3a.

Hence towards the end of it, less probable states like additional hydration barriers of aluminium atoms (CV > 3, see fig. 5.3) can be overcome. Besides parallelization, the use of multiple walkers serves the purpose of minimizing oscillations during the exploration of the free energy landscape. Nevertheless, these walkers tend to diverge and explore separate areas, given to the fact that they generate a potential that repels itself and other walkers. In this simulation, four walkers simultaneously explored the FES projection along one CV. In particular for the exploration of the more stable adsorption modes and barrier (structures I-IV) are sampled by three of the walkers, in fig. 5.3a green, purple and yellow, whereas physisorption (structure V) is only sampled by one walker (blue), which also contributes to the transition state (structure IV). At some point, the blue walker starts to determine a new local minimum at CV ≈ 4, which can occur, given that relative to structure V (ν 0 ) a potential of 250 kJ mol -1 forces this walker to proceed along the CV towards in this case irrelevant areas. In the meantime of the simulation, the remaining ones continue to improve the minima for ν 2 and ν 1 (I and III in fig. 5.3). This problem is intrinsic standard metadynamics (metaD) but also WTmetaD only in parts overcomes this problem. In particular for the overall short simulation times at DFT level, it can introduce artefacts that require excessively long simulation times to "smoothen".

In order to avoid these irrelevant areas of the free energy surface from the beginning of the simulation, additional wall biases can be introduced. These additional biases keep the walkers from diverging into areas of "irrelevant" chemistry by effectively limiting the area covered within the metadynamics simulation. For the other simulations presented in the following sections, such walls were introduced. However, in the region of interest, this does not change the overall conclusion that the adsorption process is exergonic and readily occurs at room temperature. While one can expect a similar trend for the much more complex ν de 2 bidentate formation involving two different adsorption sites, one being the tetrahedral aluminium atom α (cf. fig. 3.1ii, p. 53), the result may not be exactly the same.

According to new geochemical models based on in situ ATR-IR and zetametry, [START_REF] Roy | Phosphate Adsorption on γ-Alumina: A Surface Complex Model Based on Surface Character-ization and Zeta Potential Measurements[END_REF] a significant bonding mode of phosphoric acid at the solid-liquid interface is a monodentate binding to tetrahedral aluminium sites. This in general means that adsorption is also favorable on tetrahedral alumina, even though possibly as monodentate. Nevertheless the results in this sections indicate, that the monomer ν d 1 is in a fast equilibrium with the bidentate ν de 2 .

Phosphorylation barrier of the γ-Al 2 O 3 surface

Monodentates and most bidentates are readily accessible via substitution of µ 1 -OH hydroxy groups from preadsorbed water, as shown in section 5.2. While some of these structures were rather stable, the most stable candidates in the dry state are tridentates and a bidentate, where an alumina oxygen is substituted. However, these alumina oxygen are µ 2 -OH, µ 3 -OH or µ 3 -O, which cannot be substituted directly, since they do not easily desorb from the surface. Thus, it is important to determine how and when they can form.

Here, we test the phosphorylation reaction of those sites, which could be similar to the reaction:

ν 2 + µ 2 -O --→ ν 3 + H 2 O := (5.5)
The phosphate oxygen atoms are highlighted in red and the alumina oxygen in blue.

The H 2 O leaving the surface was formed from a P -OH moiety. This is similar to the pyrophosphate condensation, which is a phosphorylation between two orthophosphates.

It completes the reaction sequence as an exploration of reactions up to ν 3 in eq. 5.1. Since pyrophosphates are also present at the surface, this reaction could also start from an adsorbed pyrophosphate and form an orthophosphate as leaving group instead of water:

ν 20 + µ 2 -O --→ ν 3 +ν 0 := (5.6)
The same coloring of oxygen applies here with additional green labels for phosphate oxygen of the leaving orthophosphate. In this reaction, the alumina phosphorylation reaction occurs together with pyrophosphate hydrolysis, however, the only change compared to eq. 5.5 is the leaving group of H 2 PO - 4 . At least for the simulations, this leaving group cannot diffuse into bulk H 2 O, hence it is described as physisorption (ν 0 ). The presence of pyrophosphate was discussed in detail in chapter 3 and 4, with a proposed rationale for [START_REF] Luis | Hydrolysis rates of inorganic polyphosphates in aqueous solution as well as in soils and effects on P availability[END_REF] P NMR [START_REF] Decanio | FT-IR and Solid-State NMR Investigation of Phosphorus Promoted Hydrotreating Catalyst Precursors[END_REF] and direct 31 P NMR measurements [START_REF] Hühn | Structural Characterization of Phosphate Species Adsorbed on γ-Alumina by Combining DNP Surface Enhanced NMR Spectroscopy and DFT Calculations[END_REF] at the dry γ-Al The lowest reaction path is a concerted reaction regardless of the leaving group. This conclusion is not only based on the fact that the reaction path is approximately described as diagonal through the two graphs but also because it does not include intermediates in contrast to the phosphate condensation. It is clearly visible in fig. 5.5b that the phosphorylation starting from orthophosphate is endergonic (33 kJ mol -1 cf. tab. 5.1), while the reaction with pyrophosphate in fig. 5.5d is exergonic (-27 kJ mol -1 ). The free energy barriers are 146 kJ mol -1 for the orthophosphate and 115 kJ mol -1 for pyrophosphate. Overall the barrier is decreased when considering pyrophosphate but it is still slightly higher than the estimate for phosphate hydrolysis (100 kJ mol -1 ). The relevant comparison is the hydrolysis, since the pyrophosphate is split into orthophosphates (while phosphorylation, thus no H 2 O is formed here).

Regarding the elimination-activated path observed for phosphate condensation in bulk liquid, it is not expected to occur based on the high free energy region at (0.1, 1.3) in coordinates of CN a and CN w of fig. 5.5b. The free energy in the same region in fig. 5.5d seems to be lower, but one has to keep in mind that, the phosphorylation reaction with a phosphate leaving group is in competition with the analogous reaction of leaving water.

While sampling CN a , thereby cleaving the bond without immediate addition of phosphate, would lead to the elimination state around (0.1, 0.0) in coordinates of CN a and CN p .

This reacts spontaneously with water and even in another alumina oxygen site, thereby escaping the relevant FES region, i. e. the state at (0.1, 0.0) is usually not a metaphosphate (d): Free energy surface of (c). See eq. 5.5 and 5.6 on page 102 for the reaction mechanism.

structure like [(µ 1 O) 2 PO)] + (cf. fig. 5.4b). It is therefore not possible to make reliable statements about the elimination path, however, the concerted reaction seems to be favored by the system. Additional wall biases were included that artificially disfavor water adsorption in order to circumvent this problem, however these wall biases did not entirely suppress the competing reactions. In consequence, the all biases keep the walkers of the simulation "on track" and allow to sample the desired area but the elimination-activated

(0, 1) (µ B 2 -O) -H (1, 1)
O 3 P -(µ B 2 -O) -H (0, 0) (µ B 2 -O) (1, 0) O 3 P -(µ B 2 -O) + PO - 3 + PO - 3 -H + -H + Figure 5.6.:
The local coordination of µ B 2 was biased with two CV, CN a (x-axis) and CN h (yaxis), in order to describe the deprotonation and phosphorylation of this site. The coordinates (x, y) correspond to CN a and CN h . mechanism for pyrophosphate cannot be sampled reliably. It must be emphasized that this problem is intrinsic to the tested reaction for this system and independent of the added walls. This concerns all surface reactions involving pyrophosphate when using the chosen collective variables.

Bidentate ν eB

2 and tridentate ν bcB 3 on R(1 1 0)

The two most stable adsorption sites found at drying conditions are R(1 1 0)ν bcB 3 and R(1 1 0)ν eB 2 . They involve a different alumina site, the hydroxide µ B 2 -OH, which is always protonated at the water interface. This is in contrast to µ A 3 from n(1 1 0)ν deA 3 discussed in section 5.3.2; this site can be protonated as µ A 3 -OH but it is not required to be.

Unbiased AIMD simulations, which consist of 2×2 unit cells with symmetry-equivalent pairs of A/A , in total 8 sites representing the behavior of A, showed that some sites are protonated but there are always several unprotonated sites. Here, the protonation of µ A 3 was neglected, assuming that some of these sites are always deprotonated. The alumina oxygen site µ B 2 is however always protonated -it rapidly recovers this protonation in AIMD simulations, if a starting geometry involving µ B 2 -O was chosen.

This requires a deprotonation during the reaction in order to form µ That means, the deprotonation occurs independently from the phosphorylation. The reaction is possible by initial deprotonation as well as deprotonation after phosphorylation, where the latter is preferred, since the lowest free energy path leads through point ( it has to be considered that this WTmetaD simulation was much shorter than the others. As a consequence of the third CV, CN h , the overall convergence time of these simulations increases by an order of magnitude. Moreover, the carried out simulations were shorter than the simulations for n(1 1 0)ν deA 3 formation (cf. tab. 5.1). Consequently, the convergence of the presented FES differs from one another. The big differences in the barriers indicate that those are strongly affected by this, whereas the reaction free energies are consistent. This is actually expected, since the bias potential generated WTmetaD does not fully compensate for the minima in the FES ("flat target"), but converges to a "well-tempered target". It is a target which depends on the bias factor and allows sampling of the FES by flattening it up to the point where thermal overcomes the barrier (see eq. 2.11, p. 49 for details). As consequence, minima are still better sampled than transition states, because they are still visited more often throughout the simulation.

With the limited simulation time, this shows in much less accurate barriers than reaction energies.

Regarding the apparent preference of ν e 10 → ν eB 2 + ν 0 (fig. 5.8e) for an eliminationactivated reaction: this is rather likely in this case, but for the same reason as in the previous section, the minimum at (0, 0) cannot be reliably sampled due to competing reactions Therefore, this possibility was not evaluated.

Summary and perspectives

In this chapter, well-tempered metadynamics (WTmetaD) simulations were used to explore the formation of several phosphate surface species on γ-Al indicate, which of eq. 5.5 and eq. 5.6 are used. HPO 2 - 4 is considered physisorbed (ν 0 ). The hills implicitly include a simulation time, since 100 hills are added per picosecond, however, this time is the sum of four walkers producing 4 biased trajectories in parallel. dentation, only involving ν 1 adsorption sites which are accessible by water substitution, R(1 1 0)ν bc 2 , R(1 1 0)ν e 1 and n(1 1 0)ν de 2 , and the related pyrophosphate species, R(1 1 0)ν bc 20 , R(1 1 0)ν e 10 and n(1 1 0)ν de 20 . As demonstrated in section 5.2, these types of chemisorptions in form of bidentates (ν 2 ) and monodentates (ν 1 ) are the predominant species.

n(1 1 0)ν deA 3 R(1 1 0)ν bcB 3 R(1 1 0)ν eB
For benchmarking the method and as a reference for comparisons, the condensation of two orthophosphates to pyrophosphate was simulated in bulk water. The focus in all cases was also on the activation mechanism of the reaction, simultaneously covering elimination-activated, concerted and addition-activated mechanisms. Based on the results of seven simulations, the addition-activated mechanism can be excluded, whereas the elimination-activated path was only observed in bulk water. The concerted reaction is observed as the only mechanism of surface phosphorylation, however, convergence of the barriers was only satisfactory in the two simulations on n(1 1 0):

ν de 2 → ν deA 3 + H 2 O and ν de 20 → ν deA 3 + ν 0 .
In the four simulations on R(1 1 0) the barriers are much higher, but this is most likely due to inaccurate sampling of those more demanding three-dimensional WTmetaD simulations. The reaction free energy is not affected by this and shows the same trend as for ν deA 3 : the surface phosphorylation using an orthophosphate precursor is endergonic, whereas the phosphorylation using pyrophosphate is exergonic. The barrier for the pyrophosphate-based phosphorylation is with 115 kJ mol -1 similar the related pyrophosphate hydrolysis in bulk water (100 kJ mol -1 ), which means that surface phosphorylation may already occur slowly during impregnation, provided that pyrophosphate or other polyphosphates are present in solution. However, it is overall unlikely to be a driving factor of surface speciation during impregnation at room temperature, which confirms that the chemisorption described in section 5. While the results in this chapter certainly do not directly support the presence of high dentations involving the phosphorylation of alumina oxygen sites, they are not excluded either. For example in Wischert's reconstruction, [START_REF] Wischert | γ-Alumina: The Essential and Unexpected Role of Water for the Structure, Stability, and Reactivity of "Defect" Sites[END_REF] a water adsorption site is formally removed, and the R(1 1 0)µ B 2 becomes more reactive (much more basic), than n(1 1 0)µ B 3 (cf. fig. 3.1a, p. 53). More elaborate reconstructions may effectively allow the exchange of alumina oxygen sites at a lower barrier than phosphorylation, but such reactions are unknown and difficult to explore. As it seems, the surface species may change considerably during the drying step. AIMD simulations of those surfaces could be used to explore formation of such species under these conditions. Simulations of the dry surface are difficult to achieve without ambiguity, since the protonation states cannot change without solvent contact. The presented results also do not yet include the coverage considerations, in fact the model surface coverage is with 0.38 P nm -2 lower than most experiments. Further simulations of thin layers (monolayer, double layer) of different phosphoric acid content may be use to complete the atomistic description by filling the void between impregnation step and drying as well as allowing to test for surface coverage effects. This requires the use of more complex collective variables (CV) such as those described in appendix A. 4 (p. 162). Even though open questions remain, the most likely surface species were identified as bidentates and potentially monodentates which form exclusively by water substitution from µ 1 sites. These are ν bc 2 and ν de 2 but it excludes the most stable bidentate at drying conditions, ν eB 2 as it binds to an alumina oxygen site.

The formation of higher dentations involving alumina oxygen sites can be excluded during impregnation since they are endergonic at reaction conditions, except in the presence of polyphosphates in solution. As alternative to pyrophosphates during impregnation, the same types of surface phosphorylation may occur at the drying conditions, since temperatures are higher and a drop in H 2 O solvent should activate the reaction. It would be interesting to explore the drying process with a similar set of calculations in which the phosphorylation reaction is explored for orthophosphates. Different phosphate loadings could actually induce further reactions, including polyphosphate condensation. Chapter 3 indicated, which model coverages would have to be reached.

Machine learning NMR chemical shifts

During the last decade, the field of machine learning got a tremendous amount of attention in all material sciences. [START_REF] Yosipof | Data Mining and Machine Learning Tools for Combinatorial Material Science of All-Oxide Photovoltaic Cells[END_REF][START_REF] Schmidt | Recent advances and applications of machine learning in solid-state materials science[END_REF][START_REF] Kailkhura | Reliable and explainable machine-learning methods for accelerated material discovery[END_REF] The two main reasons for this growing interest are (i) the availability of computing resources for this type of research and (ii) the availability of steadily growing databases. Still, the field is very new and at this time, without generally established strategies.

There are very different problems that could be addressed by machine learning methods.

The overall idea is to build a structure-activity relationship model on complex data sets with many variables, where subtle correlations can be exploited in an automated way, without the need of human intervention to discover and quantify them. This, however, does not necessary mean that the prediction model cannot anymore be interpreted, but that interpretation follows the model construction. Interpretation of the prediction scheme is usually not straightforward and in many cases not necessary. In such cases, machine learning models are used as "black box". A good example for this would be a machine learning interpolation of the potential energy surface (PES) for geometry optimization. [START_REF] Schlexer Lamoureux | Machine Learning for Computational Heterogeneous Catalysis[END_REF] The success of this method is easily understood when considering that most geometry optimization algorithms are at best based on a parabolic approximation to the real potential hyperplane, which only applies in close proximity to the real minimum.

The machine learning approximation for the PES is used to update the structure, and the structures improve the estimate. The final result is not further transferable however, at least in case of ref. [129], because the machine learning estimates are based on cartesian coordinates, which are difficult to generalize. For that purpose, dedicated chemical descriptors exist which are explained in section 6.1.

Long before computational NMR predictions, empirical relations between [START_REF] Luis | Hydrolysis rates of inorganic polyphosphates in aqueous solution as well as in soils and effects on P availability[END_REF] P NMR and structural parameters were tested. [START_REF] Geoffrey | Predicting phosphorus NMR shifts using neural networks[END_REF][START_REF] Müller | High-resolution solid-state 27 Al and 31 P NMR: correlation between chemical shift and mean Al-O-P angle in AlPO4 polymorphs[END_REF] In case of [START_REF] Luis | Hydrolysis rates of inorganic polyphosphates in aqueous solution as well as in soils and effects on P availability[END_REF] , that both were used in previous works. [START_REF] Müller | High-resolution solid-state 27 Al and 31 P NMR: correlation between chemical shift and mean Al-O-P angle in AlPO4 polymorphs[END_REF][START_REF] Campomar | Rmn en rotation a l'angle magique du 31p. Application a l'etude de quelques phosphates[END_REF] Their work method can be applied to a range of aluminium phosphates, however, it cannot be transferred to more complex chemical environments that involve bonds other than a local structure of [P(OAl) 4 ]. Dawson and Asbrook's first study was even more specific by only studying aluminium phosphate clusters with tetrahedral aluminium atoms, which they extended for more complex compounds involving octahedral alumina. [START_REF] Daniel | Effects of Extraframework Species on the Structure-Based Prediction of 31 P Isotropic Chemical Shifts of Aluminophosphates[END_REF][START_REF] Daniel | A Multinuclear NMR Study of Six Forms of AlPO-34: Structure and Motional Broadening[END_REF] Chaker et al. recently identified a reliable method to learn chemical shifts from ab initio calculations by combining several local chemical descriptors with different machine learning methods. [START_REF] Chaker | NMR shifts in aluminosilicate glasses via machine learning[END_REF] Their work was focussed on the solid-state NMR chemical shift prediction of [START_REF] Sandy | Oxygen isotope signature of UV degradation of glyphosate and phosphonoacetate: Tracing sources and cycling of phosphonates[END_REF] Si, [START_REF] Li | Molecular level investigations of phosphate sorption on corundum (α-Al 2 O 3 ) by 31 P solid state NMR, ATR-FTIR and quantum chemical calculation[END_REF] O, [START_REF] Paglia | Determination of the structure of gamma-alumina from interatomic potential and first-principles calculations: The requirement of significant numbers of nonspinel positions to achieve an accurate structural model[END_REF]

Chemical descriptors

The value of machine learning frameworks is to accept arbitrarily complex input data and systematically process it in order to extract relevant information. Those involve regression models, i. e. machine learning algorithms like LRR, neural networks or boosted trees.

For a table with property columns and one set of data points per line, machine learning algorithms are typically applied directly, for example to train a prediction algorithm for one property based on correlation with the others. More complex, unstructured data, such as pixel information in images, require considerable adjustments in order to preprocess the information before using it as input for a machine learning algorithm.

For pixel images, those might be scanning of local image snippets on a grid which is necessary to learn features irrespective of their global position in an image. Similar considerations have to be made when processing atomistic data. There are two common types of encoding atoms in chemistry: first by defining a connectivity graph, and second by explicitly providing atomic coordinates in space. The former is most common for molecules -essentially IUPAC names of molecules are an encoding of such graphs; the latter are typically geometry representations which could be a computed or for instance crystallographic positions. Even though connectivity graphs contain much less information than the atomic coordinates, Jonas and Kuhn recently used an experimental database of 43k molecules to successfully generate accurate prediction schemes for 1 H and 13 C NMR chemical shifts using graph neural networks, with RMSE = 1.2 ppm ( 13 C). [START_REF] Jonas | Rapid prediction of NMR spectral properties with quantified uncertainty[END_REF] 119

Machine learning NMR chemical shifts

Efforts towards empiric NMR chemical shift predictions are not new, HOSE codes are well established for [START_REF] Nicosia | The effect of phosphate and glycol on the sulfidation mechanism of CoMo/Al 2 O 3 hydrotreating catalysts: an in situ QEXAFS study[END_REF] C and 1 H NMR for more than 30 years (albeit lower accuracy, for [START_REF] Nicosia | The effect of phosphate and glycol on the sulfidation mechanism of CoMo/Al 2 O 3 hydrotreating catalysts: an in situ QEXAFS study[END_REF] C RMSE = 4.2 ppm). [START_REF] Jonas | Rapid prediction of NMR spectral properties with quantified uncertainty[END_REF][START_REF] Bremser | Hose -a novel substructure code[END_REF] West explored similar graph neural network methods experimental [START_REF] Luis | Hydrolysis rates of inorganic polyphosphates in aqueous solution as well as in soils and effects on P availability[END_REF] P NMR data in 1993, but did not reach satisfactory accuracy (at best about 20 ppm). [START_REF] Geoffrey | Predicting phosphorus NMR shifts using neural networks[END_REF] The author rationalized this with errors in the included databases and weaknesses in the applied descriptors and machine learning methods at the time. This is consistent with Jonas and Kuhn, who also filtered the available database, removing about 25 % based on various criteria and added several input features (like hybridization and aromaticity) for each atom explicitly. [START_REF] Jonas | Rapid prediction of NMR spectral properties with quantified uncertainty[END_REF] It shows that the applied neural networks do not automatically derive such information, even though most of them are implicitly contained in the molecular graph. As an alternative to the connectivity graph, one can also consider atomic representations as lists of cartesian coordinates, which is physically much more accurate since these contain exact positions and with it bond distances and angles. However, even bonds are not defined in a coordinate list, this is left to be recognized by an algorithm. In some cases, machine learning algorithms can use them directly, like in interpolation schemes which produce estimates on few data points, [START_REF] Schlexer Lamoureux | Machine Learning for Computational Heterogeneous Catalysis[END_REF] but in most cases one needs to process the atom list in a way that allows for generalization over a wide range of different systems with different size. Those generalizations are chemical descriptors where the structural information in the atom list is encoded in a well defined form. A descriptor must fulfill the following requirements: [START_REF] Tanaka | Nanoinformatics[END_REF] • It must faithfully represent the properties of interest. No two different structures should generate the same descriptor.

• It should be continuous and differentiable to quantify structural similarities.

• The size and dimensionality of the descriptor needs to be constant, independent of the input size.

• The descriptor should be invariant to rotations, translations, inversion, and list permutations.

Furthermore, the descriptor can be global, capturing the whole geometry or local, atomcentered, in which case it is constructed to capture a chemical environment. It may or may not be suitable for periodic structures like crystals. A simple global descriptor is the Coulomb matrix, a matrix that contains two-body terms

Z i Z j
|r i -r j | of the Coulomb repulsion of the nuclei. The Coulomb matrix is a unique representation of the system of interest (which means that it is "faithful") and it is invariant to rotation, translation and inversion by construction, and can be made invariant to list permutations by sorting (but at the cost of differentiability). [START_REF] Faber | Crystal structure representations for machine learning models of formation energies[END_REF] Several extensions of the Coulomb matrix for periodic structures exist, [START_REF] Faber | Crystal structure representations for machine learning models of formation energies[END_REF][START_REF] Himanen | DScribe: Library of descriptors for machine learning in materials science[END_REF] as well as derived representations like the Bag of Bonds. [START_REF] Hansen | Machine Learning Predictions of Molecular Properties: Accurate Many-Body Potentials and Nonlocality in Chemical Space[END_REF] The size of a Coulomb matrix depends on the number of atoms; a constant descriptor size can be achieved by zero-padding. A different global descriptor is the many-body tensor representation (MBTR), which uses bond lengths and angles based on the chemical element. It uses one-body (g 1 , elements), two-body (g 2 , distance or inverse distances) and three-body terms (g 3 , angles in degree or cosine). These motifs are grouped by elements, element pairs and all possible element triples, and then added together on a range D.

The scalar values are then broadened with a gaussian kernel density estimate and the range is sampled on a numerical grid, similar to a histogram:

D(χ) = gn 1 σ √ 2π exp - (χ -g) 2 2σ 2
, with (6.1)

g = g 1 (Z i ), g 2 (r i , r j ), g 3 (r i , r j , r k ) (6.2)
where χ is a grid point and σ is the broadening. Altogether the ranges of g 1 , g 2 and g 3 define the descriptor. g 1 is equivalent to the system stoichiometry, g 2 and g 3 characterize element distances and angles in a way that is independent of the number of atoms in the descriptor is the smooth overlap of atomic positions (SOAP). [START_REF] Bartók | On representing chemical environments[END_REF] It encodes the local environment of a central atom as smoothened atomic positions by a local expansion in orthonormal functions based on spherical harmonics Y lm and radial basis functions g n . The coefficients c nlm of these basis functions for one atom sort are defined as inner product:

c nlm = V dV g n (r)Y lm (θ, φ)ρ(r) (6.3)
Those basis functions make the SOAP descriptor naturally a faithful descriptor of the local environment, which is differentiable, invariant to translation and permutation, and of constant size, based on a finite basis set and the number of involved elements. In order to achieve rotational invariance, the SOAP descriptor defines a partial power spectrum vector p of the angular momenta of all atoms in the vicinity. An element p Z 1 Z 2 nn l of this vector is calculated by:

p Z 1 Z 2 nn l = 8π 2 2l + 1 m c Z 1 nlm • c Z 2 n lm (6.4)
Note that in the implementation used here, real-valued (tesseral) spherical harmonics were used; [START_REF] Himanen | DScribe: Library of descriptors for machine learning in materials science[END_REF] for complex-valued spherical harmonic, the multiplication is with the complex conjugate c Z 1 * nlm . There are other local descriptors like atom-centered symmetry functions (ACSF) or angular and radial distribution functions (ARDF), [START_REF] Behler | Atom-centered symmetry functions for constructing high-dimensional neural network potentials[END_REF][START_REF] Artrith | Efficient and accurate machine-learning interpolation of atomic energies in compositions with many species[END_REF] which use a different methodology to sample the distribution of atoms in a vicinity. Both can be used to construct continuous differentiable, invariant descriptors of constant size.

The ACSF descriptor is different to the other descriptors, since its components are all 6.2. Machine learning methods user-defined allowing more flexibility, but it leaves it to the user to ensure faithfulness. [START_REF] Himanen | DScribe: Library of descriptors for machine learning in materials science[END_REF] ARDF provides a consistent framework similar to SOAP, however, the author promises a much higher efficiency for systems involving multiple chemical species. [START_REF] Artrith | Efficient and accurate machine-learning interpolation of atomic energies in compositions with many species[END_REF] Chaker et al. tested SOAP, ARDF and a simple implementation of ACSF involving only two-body terms for performance in predicting NMR chemical shifts. They demonstrated that of these descriptors, SOAP clearly leads to the most accurate prediction models.

In this work, the SOAP descriptor was compared with the LMBTR in order to test the performance of this new descriptor. In comparison to SOAP, LMBTR explicitly encodes three-body terms including all atoms within its vicinity. This should be advantageous in comparison to SOAP, since angles, such as the P -O -Al angle discussed by Dawson and Ashbrook [START_REF] Daniel | Investigating Relationships between the Crystal Structure and 31 P Isotropic Chemical Shifts in Calcined Aluminophosphates[END_REF] needed to accurately predict the 31 P NMR chemical shift.

Machine learning methods

Overall, the two descriptors SOAP and LMBTR were combined with three machine learning methods: linear ridge regression (LRR), [START_REF] Hilt | Ridge, a computer program for calculating ridge regression estimates[END_REF] the multilayer perceptron (MLP) [START_REF] Haykin | Neural Networks and Learning Machines[END_REF] and extreme gradient boosting (XGBoost). [START_REF] Chen | XGBoost: A Scalable Tree Boosting System[END_REF] LRR is a simple technique in which a regularization term is added to a least square fit:

LRR = min       i y i - k (m k x ik -b) linear least squares + λ k m 2 k regularization       (6.5)
where x ik are k independent variables for data points i for regression and y i values are to be predicted. The linear least square fit part needs to be overdetermined (i > > k) in order to quantify uncertainties. On the other hand, least square fit is not possible if the system is underdetermined (i < k). The regularization term is effectively a bias, leading to a systematic deviation from the least square fit results, but it provides solutions to 123 the minimization in eq. 6.5 even if the system is underdetermined. While this seems ill-advised for the low-dimensional case or in case the independent variables are truly independent, [START_REF] Campbell | A Critical Analysis of Ridge Regression[END_REF] ridge regression can be very useful to identify the most important components in underdetermined data sets, provided that the variables are strongly correlated, i. e. linear dependent. [START_REF] Gruber | Improving Efficiency by Shrinkage: The James-Stein and Ridge Regression Estimators[END_REF] Like indicated, LRR works extraordinary well for linear regressions in which the input is heavily correlated, but it fails where non-linear components are crucial. The parameter λ is a so-called hyperparameter, since it is fixed by the user, as opposed to the parameters m k and b which are optimized in the regression.

MLP and XGBoost are more advanced techniques; the former is an artificial neural network and the latter based on boosted decision tree algorithms. They offer more flexibility and are suitable for very complex data structures. It means that while LRR cannot take into account nonlinear components, and hence is inherently limited, the prediction accuracy of MLP and XGBoost can be arbitrarily improved by adding more data (and adding parameters). MLP requires defining a neural network topology between input neurons (defined by the feature array) and the output neurons (here, a single one that returns a chemical shift estimate). This can be very complex, as one can choose the number of neuron layers as well as the dimensions of each layer and the connectivity between all the neurons in every layer. For simplicity, here the number of layers is limited to a single, fully connected layer. Only the number of neurons was varied.

The data set

For the presented results, the DFT GIPAW calculations of the phosphate species at single adsorption, coadsorption and pyrophosphates. An even smaller subset of those was discussed in chapter 4. This number of data points is small for machine learning purposes, hence additional GIPAW calculations were carried out for the less stable hydrogen bond networks from the previous sampling. Only orthophosphate coadsorption and pyrophosphate adsorption cases were used as they produce two data points per calculation. With the additional calculations, the number of data points was increased to 1227. All data points were used without filtering outliers. While it does not seem to disturb machine learning results, in one (faulty) system a peroxide geometry PO 3 (O 2 ) 3 - was identified, which apparently did not substantially affect the results, and was therefore kept. Fig. 6.1a provides a comparison between the original set of 361 data points of the most stable structures, and the extended set of 1227 data points. A core issue of machine learning is overfitting, which means that a machine learning algorithm can produce accurate fits of the data, but the result is not transferable to new data. The RMSE of the training set actually increases, since with more data points, the degree of overfitting is reduced. The test and validation sets (explanation is given the next section) are on the other hand significantly improved. This effect on the training accuracy is due to the fact that the same method needs to find a regression on more data points for a constant set of features, which overall reduces overfitting, but also shows that one can only reach limited accuracy with the used algorithm based on SOAP/LRR. Fig. 6.1b shows a more systematic plot of this size effect. The size effect was in this case determined by picking data points from the large data set at random, which produces clear trends. The RMSE of training and test/validation converges with additional data, but the RMSE is still larger than zero. Based on the full data set, the RMSE converges between 2.7 and 3.6 ppm.

This simply reflects the systematic error introduced by the machine learning model itself, in this case SOAP/LRR. Consequently, in order to improve the prediction, the machine learning model has to be changed, either by an improved structural descriptor or a better fitting algorithm.

Parametrization

Parametrization in machine learning is a global optimization problem in which hyperparameters of an algorithm are varied in order to maximize prediction performance. In this work, it means to minimize the RMSE of the machine learning predictions of 31 P GIPAW results. By combining the descriptors SOAP and LMBTR with the machine learning algorithms LRR, MLP and XGBoost six different methods were generated. Besides the hyperparameters, the descriptors also depend on parameters. Since the parameters depend on one another in non-trivial ways, for all six combinations test calculations were carried out. To keep track of overfitting, cross-validation and testing were done as described by Chaker et al. [START_REF] Chaker | NMR shifts in aluminosilicate glasses via machine learning[END_REF] A test set of 10 % of the database was separated at random.

For cross-validation, the remaining 90 % were split in 5 groups, 72 % for training and 18 % for validation. Training and validation were conducted for 5 times, i. e. in every case and the fraction for the validation set was swapped in the groups. Based on this procedure, descriptor parameters and machine learning hyperparameters were scanned.

Descriptor parameters

The SOAP descriptor can be manipulated by defining different basis sets of the radial (nmax) and angular (lmax) components. Furthermore, one has to define a radial cutoff which restricts vicinity and with it the number of included atoms. Lastly the atomic positions are smeared out by a gaussian kernel which can be adjusted by a parameter sigma. We tested the radial cutoff in a range from 5 to 9Å and equal values for nmax and lmax from 2 to 9. This corresponds to a feature size of 108 for lmax/nmax = 2 and 4752 for lmax/nmax = 9. Fig. This is consistent with findings by Chaker et al., [START_REF] Chaker | NMR shifts in aluminosilicate glasses via machine learning[END_REF] who also observed overall lower performance for more complex methods.

For LMBTR, a slightly different approach was used. In contrast to SOAP, in LMBTR ranges and grid resolution of such ranges have to be defined for the two-body and threebody terms. The two-body terms were defined in inverse units to improve resolution for small distances from 0.182 to 1.2Å -1 , which corresponds to a cutoff value at 5.5Å.

two-body (n 2 ) and three-body (n 3 ) ranges were tested with numeric grid densities from 20 to 100 corresponding to descriptor sizes of 800 to 4000 features. Fig. 6.3 contains the results for the tests of parameters. For n 2 < 50, or n 3 < 50 more irregularities are observed. In particular, LMBTR/MLP leads to bad predictions if n 2 and n 3 are not balanced. LRR is as before better than XGBoost, which is itself better than MLP.

Since 28 % of the database were used for test and validation, just about 880 data points are included in each training, which means that the system was in most cases underdetermined for regression, but the features were strongly correlated. The parameter scans of LMBTR seemed overall more noisy than the results for SOAP; this is misleading tested together with the SOAP parameters by performing a three-dimensional scan of the SOAP parameters and each hyperparameter. Fig. 6.4b and c provide the best RMSE results of the SOAP parameter scan for a given hyperparameter. As expected, the lowest RMSE were observed for the training set, while RMSE for validation and test set are similar, which means that in all cases a certain degree of overfitting was observed due to the finite database. The 'J' shaped curves for increasing λ indicate an initial performance improvement followed by monotonous increase of RMSE where regularization overly restricts the regression procedure. MLP converges with respect to the number of neurons while in contrast, increasing tree depth quickly leads to overfitting with XGBoost. Regarding overfitting, it also seems clear that it cannot be entirely avoided for a given data set with LRR and MLP. The difference between training RMSE and validation/testing RMSE was in both cases nearly constant 0.9 ppm. With XGBoost, the tree depth rapidly improved the training RMSE, while testing/validation RMSE only improves only up to a tree depth of 3. Increasing the tree depth therefore can improve the prediction performance, but it quickly leads to overfitting.

SOAP vs. LMBTR

The previous sections have shown that for the conducted tests, LRR outperforms the two more complex machine learning algorithms. This is consistent with the results of Chaker et al. [START_REF] Chaker | NMR shifts in aluminosilicate glasses via machine learning[END_REF] prediction scheme is based on the descriptor, SOAP, rather than the algorithm, LRR. Fig. 6.5 provides a direct comparison between both methods. As opposed to the previous cross-validation scheme, the machine learning scheme uses Each point is an LMBTR g 2 grid point (cf. eq. 6.1), lines added for guidance.

The y-axis describes the slope (m k ) in eq. 6.5. Positive coefficients describe higher NMR shielding compared to the average structures.

since their methodologies are very different. However, both provide distance and angular information. SOAP encodes it as two-dimensional power spectrum based on gaussian radial functions and spherical harmonics, hence distances and angles are grid positions in the spectrum. This is equivalent to other descriptors which include only two-body terms, including angular information but only from the central atom X with the motif A-X-B, even though implicitly this also includes other angles. LMBTR also encodes angles around neighboring atoms such as X-A-B directly. Simple machine learning algorithms like LRR cannot recognize implicit angular information, thus the SOAP descriptor is inherently restricted by these aspects. This is also in agreement with findings by Dawson and Ashbrook who parametrized distances and angles manually using a multivariate regression. [START_REF] Daniel | Investigating Relationships between the Crystal Structure and 31 P Isotropic Chemical Shifts in Calcined Aluminophosphates[END_REF] LMBTR can be used as one possible way of generalizing their findings.

As discussed earlier, both SOAP and LMBTR discard information compared to the original geometry by smoothening atomic positions, and forming smooth distributions from distances and angles. Fig.

LRR performance

It is important to note that the chemical descriptors SOAP and LMBTR are both derived from the geometry and encoded in a numerical form that allows to unambiguously capture a wide range of chemical environments. For this, both methods are excessively large, saving hundreds to thousands of values as a fingerprint of generally less than 20 atoms in the vicinity. It is therefore expected that these values, when used as independent variables for regression, are strongly correlated. As discussed, ridge regression methods perform very well in automatically identifying such correlations and scaling such correlated features simultaneously. It seems that NMR chemical shift predictions are in large parts driven by a linear correlation to structural features. Up to this point, it was more important to identify the most important structural features in order to generate an 31 P NMR prediction scheme, than to find the most suitable algorithm. Dawson et al. [START_REF] Daniel | Effects of Extraframework Species on the Structure-Based Prediction of 31 P Isotropic Chemical Shifts of Aluminophosphates[END_REF][START_REF] Daniel | Investigating Relationships between the Crystal Structure and 31 P Isotropic Chemical Shifts in Calcined Aluminophosphates[END_REF][START_REF] Daniel | A Multinuclear NMR Study of Six Forms of AlPO-34: Structure and Motional Broadening[END_REF] started from hand-selected structural features like bond distances and angles, whereas Chaker et al. [START_REF] Chaker | NMR shifts in aluminosilicate glasses via machine learning[END_REF] combined recently developed local descriptors and various more complex machine learning algorithms. Here, we also tested two algorithms despite LRR, but in agreement with Chaker et al., we realized that those do not outperform simpler regression schemes.

Nevertheless, a more elaborate machine learning algorithm should outperform LRR if the descriptor captures all the relevant structural motifs, as well as potentially other quantities such as atomic charges (e. g. by Bader or Hirshfeld) Successful examples for the use of neural networks for NMR prediction exist in the literature, however, those also heavily rely on the design of the descriptor. [START_REF] Jonas | Rapid prediction of NMR spectral properties with quantified uncertainty[END_REF]

Summary and perspectives

This machine learning project is not yet complete; however, it already provides valuable insights into the challenges for building a machine learning-based 31 P NMR chemical shift predictor. Recent work showed that both complex algorithms like neural networks and simple regression schemes like LRR likewise allow for accurate machine learning chemical shift predictions, [START_REF] Chaker | NMR shifts in aluminosilicate glasses via machine learning[END_REF][START_REF] Jonas | Rapid prediction of NMR spectral properties with quantified uncertainty[END_REF] At the current state many perspectives for future work are possible. One important task is to further ensure the consistence of the database. Optimizing all hyperparameters of the tested descriptor/algorithm combinations is necessary in order to complete this work. Further adjustments to the LMBTR descriptor like exchanging the representation of the angular part (g 3 ), using the angle directly instead of the cosine, has now proven to be more accurate. The search could also include squared components of the g 2 component of the descriptor, in alignment with the work of Dawson and Ashbrook. [START_REF] Daniel | Investigating Relationships between the Crystal Structure and 31 P Isotropic Chemical Shifts in Calcined Aluminophosphates[END_REF] In order to improve the system description, atomic charges might improve accuracy. They can be derived from DFT results as Bader or Hirshfeld charges. Testing more complex algorithms could also improve the prediction performance. Elastic net regression might improve the learning model, since the method is more suited for removing irrelevant components. Since the XGBoost worked best with the shallowest trees, so-called stumps (tree depth = 1), a simpler boosted tree method, AdaBoost, could be an alternative.

Based on the available GIPAW data, a similar machine learning method could also be applied to 27 Al NMR. This would provide insights into NMR shifts of aluminium surface species. Overall, using machine learning for [START_REF] Luis | Hydrolysis rates of inorganic polyphosphates in aqueous solution as well as in soils and effects on P availability[END_REF] facet, phosphate species stabilize the reconstruction of the surface, which was earlier identified at high-temperature models in which undersaturation of water adsorbates is accommodated by shifting a surface aluminium atom from an octahedral to a tetrahedral site, thereby reducing the degree of undersaturation. The stabilization by phosphates is clearly different, since the discussed systems are always coordinatively saturated. Hence it could also be expected that the reconstruction is (in the absence of chemisorbed phosphates) less stable at drying conditions. We also identified an alumina oxygen site, which is destabilized in the reconstruction and which is much more stable after phosphorylation. This initial exploration in chapter 3 was the basis for all following chapters, which was also a motivation for the phase space search that was carried out.

While exclusively theoretical models can improve our understanding of the more complex real systems, its foundation should be regularly justified on experimental grounds.

In chapter 4, magnetic shielding tensors of the explored model surface species were calculated using the GIPAW method. The magnetic shielding was used in order to estimate experimental 31 P chemical shifts. Overall, the model structures could reproduce most of the experimental observations. 31 P chemical shifts of phosphate species are well-known to produce a broad signal band that is broadened and shifted to lower values when increasing phosphate loading. The comparison between the coverage-dependent chemical shifts allowed us to elaborate on changes that actually occur at the atomic scale.

We can distinguish surface between the interaction effects from (i) the adsorption site, (ii) the adsorption mode and (iii) direct coverage effects (influence of phosphate increase in the vicinity of the adsorbed phosphate). The strongest effect identified was the site effect (i), which means that the relative weight of the site distribution changes with the phosphate loading. The free enthalpy estimates from DFT confirm this in comparison between single orthophosphate adsorption versus coadsorption of two orthophosphates on the (1 1 0) facet. After this the adsorption mode (ii) influences the chemical shift, in particular bidentates (ν 2 ) show lower chemical shift values than tridentates (ν 3 ). For the direct coverage effect (iii), additional chemical shift calculations with larger unit cells (lower coverage) were carried out for the selected adsorption sites. The interaction of orthophosphate with the surface can rationalize 31 P chemical shifts down to -15 ppm, but the experimental chemical shift bands can reach values down to around -30 ppm. [START_REF] Decanio | FT-IR and Solid-State NMR Investigation of Phosphorus Promoted Hydrotreating Catalyst Precursors[END_REF][START_REF] Hühn | Structural Characterization of Phosphate Species Adsorbed on γ-Alumina by Combining DNP Surface Enhanced NMR Spectroscopy and DFT Calculations[END_REF][START_REF] Van Eck | A Solid State Double Resonance NMR Investigation of Phosphorus-Impregnated γ-Al 2 O 3[END_REF] Two explanations for this trend exist: polyphosphate formation and aluminium phosphate [START_REF] Van Eck | A Solid State Double Resonance NMR Investigation of Phosphorus-Impregnated γ-Al 2 O 3[END_REF] and DeCanio et al. [START_REF] Decanio | FT-IR and Solid-State NMR Investigation of Phosphorus Promoted Hydrotreating Catalyst Precursors[END_REF] also depends on the overall phosphate loading. Our experimental data and our computational models show that up to a phosphate loading of 4.5 wt% or 4.1 P nm -2 , the adsorption site distribution and An important open question from chapters 3 and 4 is how the most stable adsorption structure can actually form at the involved experimental conditions. The initial impregnation is done in aqueous solution, in which phosphate from bulk water reaches the surface through the solid-liquid interface. Do the most stable structures at drying conditions already form during impregnation or does the solid-liquid interface produce different adsorption modes? This was studied in chapter 5. As an initial benchmark, the reaction mechanism of orthophosphate condensation to pyrophosphate was tested by well-tempered metadynamics. Even though the reaction was not yet simulated by this approach, its reaction free energy and the barriers are known experimentally and were successfully reproduced with this method. While the initial adsorption could also be assumed to happen, based on experimental observation, it was tested by a well-tempered metadynamics exploration. It was confirmed for octahedral aluminium sites, that water substitution by phosphates readily occurs for the explored case, leading to strongly adsorbed bidentates (ν 2 ). The phosphorylation of alumina oxygen sites in order to form tridentates was explored by well-tempered metadynamics for the three most stable adsorption sites identified at drying conditions on the (1 1 0) facet. It also included one case of a bidentate constructed from an alumina oxygen. Two configurations for the reactions were used for each site: preadsorbed orthophosphates ν 2 ----ν 3 +H 2 O and pyrophosphates ν 20 ----ν 3 +ν 0 , which only differ in the leaving group, i. e. water and another orthophosphate. The bidentate (ν 2 ) starts from ν 1 and ν 10 . The trend in all cases is that the reaction free energy is positive for all cases involving orthophosphates, and negative for the reaction starting from pyrophosphate. This means that, based on reaction energies, none of the reaction products are stable, as even after the pyrophosphate hydrolysis, the phosphorylated alumina oxygen can be hydrolyzed as well, i. e.

ν 20 --→ -ν 0 ν 3 -→ ν 2 .
Only two of the six surface phosphorylation reaction barriers were converged; compared to pyrophosphate hydrolysis in water, with a calculated value of 100 kJ mol -1 , they showed a higher phosphorylation barrier (115 kJ mol -1 starting from n(1 1 0)ν 20 , 146 kJ mol -1 starting from n(1 1 0)ν2). Consequently, since pyrophosphate hydrolysis is a slow reaction at room temperature, taking overall several days, higher dentations cannot even form temporarily during impregnation based on such barriers.

These barriers also indicate that, once formed, the tridentates should be stable enough to resist hydrolysis for some time. It is also consistent with the results of Van Cleve et al. who demonstrated how phosphonate adsorption on γ-Al 2 O 3 improves hydrothermal stability. [START_REF] Van Cleve | Enhanced Hydrothermal Stability of γ-Al 2 O 3 Catalyst Supports with Alkyl Phosphonate Coatings[END_REF] Overall, dentations involving alumina oxygens do not occur at the solid-liquid interface, but they are more likely under drying conditions. Once formed, these bonds are quite stable. At this point, the mechanistic studies under drying conditions are a very interesting topic for future studies.

As a result, we clearly see that higher dentations involving alumina oxygen are limited to drying conditions. This fact can also hold as an a posteriori justification for using 120°C as reference of the dry model surfaces (chapter 3) for comparison with the experimental involving many more data points per GIPAW calculation or cluster models. [START_REF] Daniel | Effects of Extraframework Species on the Structure-Based Prediction of 31 P Isotropic Chemical Shifts of Aluminophosphates[END_REF][START_REF] Daniel | Investigating Relationships between the Crystal Structure and 31 P Isotropic Chemical Shifts in Calcined Aluminophosphates[END_REF][START_REF] Daniel | A Multinuclear NMR Study of Six Forms of AlPO-34: Structure and Motional Broadening[END_REF][START_REF] Chaker | NMR shifts in aluminosilicate glasses via machine learning[END_REF] The recent local chemical descriptor LMBTR (local many-body tensor representation) [START_REF] Huo | Unified Representation of Molecules and Crystals for Machine Learning[END_REF] was applied for the regression, leading to a significant improvement over one of the most popular local descriptors (smooth overlap of atomic positions, SOAP). [START_REF] Bartók | On representing chemical environments[END_REF] This opens the possibility for future work: a study of the chemical shift including dynamic properties of the system, potentially quantifying the proportions at which different adsorption dynamics, for which GIPAW calculations would be prohibitively expensive. A similar study could also be conducted for 27 that is applied on a wave function |Ψ in order to determine the total energy, E, of the system: [START_REF] Schrödinger | An Undulatory Theory of the Mechanics of Atoms and Molecules[END_REF] Ĥ |Ψ(r N , r e ) = E |Ψ(r N , r e ) (A.1)

For non-stationary states with its time-dependent form,

ih d dt |Ψ(r e , r N , t) = Ĥ |Ψ(r e , r N , t) (A.2)
The wave function |Ψ is a descriptor that -in accordance with Heisenberg's uncertainty principle -represents the many-body system of particles by probability amplitudes of their positions r e (electrons) and r N (nuclei). By taking advantage of the Dirac notation, [START_REF] Szabo | Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory[END_REF] we can use a short form to multiply with the conjugate transpose of the wave function (implying that it is represented by some matrix), in order to separate observables and wave functions. For example, eq. A.1 becomes

Ψ(r N , r e )| Ĥ |Ψ(r N , r e ) = E (A.3)
Analogously to classical mechanics, the Hamilton operator for charged particles can be written as a sum of all operators that determine the kinetic energy of translation T and Coulomb potential energy V operators for every particle and particle pair: [START_REF] Born | Zur Quantentheorie der Molekeln[END_REF][148] results by even current standards, [152] while it was published in 1934 where calculations were limited to the simplest atomic or diatomic cases. Only in 1951 the Hartree-Fock method was supplied with an improved mathematical framework by Roothaan as one of the first steps that allowed it to be run on computers. [153] The Møller-Plesset perturbation method is the earliest post-Hartree-Fock method and was long seen as "gold-standard" for ab-initio models of organic molecules, however its approximation is only accurate for large energy gaps between occupied and virtual orbitals, and, as opposed to the Hartree-Fock method, it is not variational. As a consequence, perturbation methods are difficult to improve other than by including higher order correction terms. It is common to interpret the additional terms as excitations with first order terms as single excitations, second order terms as double excitations, etc, however, rather than real excitation phenomena, they are terms formed by excited state Slater determinants from the ground state optimized Hartree-Fock orbitals. Since the single excitations do not overlap with the ground state, the first perturbative correction is achieved by inclusion of double excitations, within the second order Møller-Plesset perturbation method (MP2).

Ĥ = i Ti + i Vi = - h2 2 i ∆ i m i + 1 4π 0 i j<i q i q j |r i -r j | (A.
Higher order corrections like triple excitations in MP3 or beyond may improve the accuracy of the result at computational expense and numerical stability. The fact that virtual excited states are included also leads to high demands for the used basis functions, needing a more "complete basis" set in order to be accurate.

In modern jargon, Hartree-Fock and in particular post-Hartree-Fock methods termed as wave function methods. They include methods that can be exact in principle, like configuration interaction and are rather computationally demanding, limiting its applicability to a smaller scale. The reason is the electron-electron repulsion term which rapidly increases computational cost with the number of electrons which leads to an asymptotic computation time of O(n 4 ). Many alternatives were proposed, which are related with wave function methods in the sense that they solve the electronic Schrödinger equation while replacing some of the terms in eq. A.5 by numerically simpler ones. The most common is based on density functional theory (DFT) in which the explicit calculation of the largest sum Vee can be substituted. While otherwise iteratively minimizing the total energy. Density functional theory (DFT) is presented in the following section in more detail. Empirical or semiempirical methods can be characterized by replacing the iterative procedure by empirical parameters, for example to build a wave function and access orbital energies and total energies by a single diagonalization, or by replacing the whole electronic wave function in eq. A.6 by classic spring potentials and force fields that mimic forces of the electron structure on the nuclei.

A.1.2. Density functional theory

DFT arisen with the study of the inhomogeneous electron gas of Hohenberg and Kohn, who proved that the density distribution ρ 0 can completely describe the physics of a quantum mechanical ground state system, like the electronic wave function |ψ e (r e ) . [START_REF] Hohenberg | Inhomogeneous Electron Gas[END_REF] The Hohenberg and Kohn theorem was reformulated and generalized over the years, [START_REF] Garrigue | Hohenberg-Kohn theorems for interactions, spin and temperature[END_REF] but it is just a proof of existence and does not give any guidance on how to replace a wave function by a density. It was also demonstrated to break down on excited states; [START_REF] Gaudoin | Lack of Hohenberg-Kohn Theorem for Excited States[END_REF] while this does not disqualify DFT as a method, it shows its limits. It nevertheless proved to be superior in terms of computational efficiency and comparable accuracy.

The general idea is that the ground state density of a system is, like the total energy, an observable of the system. It can be determined from the wave function by summation over the orbital densities: Coulomb interactions, which can be expressed as integrals over the density: [START_REF] Kohn | Self-Consistent Equations Including Exchange and Correlation Effects[END_REF] E el [ρ 0 ] = T e [ρ 0 ] + V Ne [ρ 0 ] + V ee [ρ 0 ]

= T e [ρ 0 ] + dr V N (r)ρ 0 (r) + e 2 2 drdr ρ 0 (r)ρ 0 (r ) |r -r | (A.9)

V N is the static Coulomb field of the nuclei, which is commonly generalized as external field. The integration over r and r is over the density and pair density and is usually implemented as numerical integration. With these integrals it is possible to calculate the potential energy components of the electronic structure exactly, and only based on the electron density. However, there is currently no exact density functional known for the kinetic energy T e [ρ 0 ]. Thomas and Fermi independently found an analytical approximation for the kinetic energy in 1927: [START_REF] Thomas | The calculation of atomic fields[END_REF][START_REF] Elliott | The Thomas-Fermi theory of atoms, molecules and solids[END_REF] T TF e [ρ 0 ] = 3h 2 40m e 3 π 2 3

[ρ 0 (r)] While this model was also extended to account for exchange energy correctly, [START_REF] Weizsäcker | Zur Theorie der Kernmassen[END_REF] the approximation is not even accurate enough to predict covalent bonds in molecules. [START_REF] Elliott | The Thomas-Fermi theory of atoms, molecules and solids[END_REF] This problem can in principle be overcome by finding a better approximation for T e [ρ 0 ].

Recently Snyder et al. demonstrated that it is possible to find an improved kinetic energy term by machine learning, [START_REF] Snyder | Finding Density Functionals with Machine Learning[END_REF] and several other groups explored different possible machine learning approaches. [START_REF] Yao | Kinetic Energy of Hydrocarbons as a Function of Electron Density and Convolutional Neural Networks[END_REF][START_REF] Meyer | Machine Learning Approaches toward Orbital-free Density Functional Theory: Simultaneous Training on the Kinetic Energy Density Functional and Its Functional Derivative[END_REF] While these may improve in the future, and allow systematic improvement of a "pure" kinetic energy functional, the most established method involves the same mean-field kinetic energy terms as in the Hartree-Fock method. Kohn and Sham combined the DFT terms in eq. A.9 with a wave function build from single Slater determinant of non-interacting electrons in a fictitious potential. The Kohn-Sham kinetic energy can be expressed by the orbitals |ϕ i (r) of this wave function:

T KS e [ρ 0 ] = i ϕ i (r)| - h2 ∆ 2m |ϕ i (r) (A.11)
The Kohn-Sham kinetic energy is much more accurate than any analytical approximation, but introduces the subtle errors, inherited from the mean field approximation. It is generally collected in an unknown term, the exchange-correlation potential V xc [ρ 0 ], which is empirically approximated based on the density. It is the only approximative term of the Kohn-Sham energy:

E KS el [ρ 0 ] = T e [ρ 0 ] + V Ne [ρ 0 ] + V ee [ρ 0 ] + V xc [ρ 0 ] (A.12)
The non-interacting electrons in the wave function part are produced by SCF minimization.

Their orbital energy levels ε i are given by:

ε i = ϕ i (r)| - h2 ∆ 2m + V s (r) |ϕ i (r) (A.13)
The fictitious potential V s (r) depends on the density itself: As alternative, so-called hybrid DFT methods substitute varying amounts of their approximate exchange by exact exchange from Hartree-Fock terms. Perdew described this process of including increasingly complex terms into a successively more accurate method with the figure of speech of Jacob's ladder. [START_REF] Perdew | Jacob's ladder of density functional approximations for the exchange-correlation energy[END_REF] In more recent years the focus was on improving the description of dispersion interaction, which is included in wave function methods like MP2 but at high computational cost.

Grimme et al. created several method agnostic correction potentials which, irrespective of the functional, add an empirical term based on the atomic positions. [START_REF] Grimme | A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu[END_REF] These schemes aim to determine long range correlation effects by parametrized force fields independent of the density. The method was also adapted for charge-dependent polarizability by estimating atom charges from the density. [START_REF] Steinmann | Comprehensive Benchmarking of a Density-Dependent Dispersion Correction[END_REF] On the other hand Klimeš et al. built whole new functionals from scratch. [START_REF] Klimeš | Chemical accuracy for the van der Waals density functional[END_REF] Overall, DFT methods are still evolving, with continuous development of more accurate functionals and efficient implementations.

A.1.3. Basis functions and pseudopotentials

The electronic structure is described by a set of mathematical functions that do not only have to be well suited to describe its topology but they also need to ensure numerical stability for ground state optimization. At the same time it must be possible to evaluate their accuracy and transferability, while in particular small basis sets should be reasonably accurate.

In computational chemistry, which emerged from wave function methods, in the past almost exclusively Gaussian orbitals were used. Gaussian orbitals are localized at the atomic positions and are parametrized to describe atomic orbitals and molecular orbitals alike. Angular components are included by multiplying with spherical harmonics. They are commonly "contracted" in fixed linear combinations to approximate more physically sound Slater-type orbitals without increasing the number of available basis functions.

Since Gaussian orbitals are not orthogonal, but the determined molecular orbitals must be orthogonal, one has to ensure, that linear dependencies are avoided. This is unproblematic as long as the Gaussian basis set does not include too diffuse basis functions.

In periodic boundary conditions the most "natural" functions are plane waves (PW).

They are likely the most popular basis set in computational solid state physics for decades.

Plane waves are similar to a particle in a box, where a quantum mechanical particle is described by a superposition of a infinite discrete set of trigonometric functions. In contrast to Gaussian basis sets, PW basis sets are defined by the unit cell dimensions and independent of the number of atoms in the system. To actually allow computations with PW basis sets, it has to be truncated to a finite size, usually by choosing a maximum frequency. Plane waves are ideal to describe smooth continuous densities, however, local density spikes like atom nuclei would require a very large number of plane waves. For the interfacial systems described here, it is also important to note that even vacuum volume in surface models contributes to the basis set size.

To effectively save computation time without affecting the accuracy of the chemically most relevant valence states (which are relatively smooth), the core electrons are substituted by pseudopotentials. Pseudopotentials are in essence radial functions that describe the same electrostatic field as the nucleus at long distances, while the core region is optimized in such a way that the remaining electron orbitals as well as the overall density remains the same outside of the core region given by the pseudopotential. Thus, effectively the number of explicitly calculated electrons is reduced to valence electrons.

Calculations with Gaussian basis sets can be done without pseudopotentials, but in particular heavy elements are not accurately represented without relativistic models of the core electrons. To avoid more costly relativistic methods, the pseudopotentials are constructed including relativistic effects. For PW methods, pseudopotentials are crucial in order to smoothen the electron density in close proximity to the nuclei. This in fact changes the density description in close vicinity to the nucleus without physical justification. A simple method for this is the use of norm-conserving pseudopotentials. [START_REF] Hamann | Norm-Conserving Pseudopotentials[END_REF] While the shape of the density close to the nuclei is changed, the number of valence electrons is maintained. Ideally, one would be able to recover the all-electron density after carrying out pseudopotential calculations. This can be in principle be achieved within Blöchl's projector-augmented wave (PAW) formalism. [START_REF] Blöchl | Projector augmented-wave method[END_REF] The method uses a much smaller PW basis to built Kohn-Sham orbitals than for the electron density. With the method, all-electron densities can be produced from pseudo-orbitals. Nevertheless the most common application of PAW is as theoretical generalization of Vanderbilt's ultrasoft pseudopotentials (USPP), [START_REF] Vanderbilt | Soft self-consistent pseudopotentials in a generalized eigenvalue formalism[END_REF] as for example implemented in VASP. [START_REF] Kresse | From ultrasoft pseudopotentials to the projector augmented-wave method[END_REF] The motivation of USPP is, rather than constructing an all-electron density, to soften the requirement for the orbital plane waves beyond limitations due to norm conservation by introducing local components to the plane wave orbitals.

In Kohn-Sham DFT, the kinetic energy is calculated by the construction of noninteracting Kohn-Sham orbitals in an effective field that is depends on the electron density. The other terms of the total energy are directly determined as functionals of the electron density. For the choice of basis functions and pseudopotentials, this combination of orbitals and electron density leads to inconsistent numeric requirements.

Calculations based on the density improve with a larger, more flexible basis set, but the computational cost for orbital energy calculations in plane wave basis functions increases rapidly, ultimately leading to the bottleneck of this method. As described above, within Gaussian basis sets this problem is addressed by keeping the number of available basis functions constant, and improving their quality in contractions that are more suitable to describe orbitals than single Gaussian orbital functions. For the PAW method, a small PW basis set is used for description of the orbitals, and a larger auxiliary basis for the density.

The Gaussian and plane wave (GPW) method does the same separation between Kohn-Sham orbitals and density calculations, where the former is calculated with Gaussian orbitals and the latter -after transformation of the Gaussian basis -in an auxiliary plane wave basis. The Gaussian basis sets are adapted for the use of the same pseudopotential as the auxiliary basis. [START_REF] Kühne | CP2K: An electronic structure and molecular dynamics software package -Quickstep: Efficient and accurate electronic structure calculations[END_REF] Compared to PW methods, GPW is similar to the use of normconserving pseudopotentials, both basis sets rely on the same pseudopotentials. However it is much more efficient and scalable since orbitals and densities can be much better evaluated in those basis sets, while additional computational load by transformation is usually negligible. as 1/2 occupancies. These were also fully relaxed at the DFT level. The other two geometry tests were full optimizations at DFT level, one with the standard functional of this work, PBE-D3, and another test of TPSS-D3(BJ), a more accurate but also more costly method. Tab. A.4 summarizes the results for the GIPAW calculations. A limitation of the VASP code is the lack of pseudopotentials for any functional that is not PBE. This problem might have unpredictable effects when using a different functional to optimize the pseudized orbitals, since we then mix the description of core-electrons (described by PBE) and valence electrons (described by TPSS). Grimme's dispersion corrections were used in all tests. [START_REF] Grimme | A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu[END_REF] Hence for the test with TPSS both cases are presented, one that uses TPSS-D3(BJ) only for geometry optimization, and another that is based on the electronic structure of TPSS (last two rows in tab. A.4).

Overall, the different structures lead to significant iso-shift variations, up to 13.3 ppm for wavellite and 11.1 on average. While the clear choice is to use the same functional that

1. 1 . 31 P

 131 NMR iso-shift and anisotropy . . . . . . . . . . . . . . . . . . . . . . . 3.1. Tested adsorption sites on reference surface models . . . . . . . . . . . . . 3.2. Lowest adsorption free enthalpies at drying condition. . . . . . . . . . . . 3.3. Site comparison between n(1 1 0) and R(1 1 0) . . . . . . . . . . . . . . . . 3.4. Best two pyrophosphate and coadsorption cases . . . . . . . . . . . . . . . 5.1. WTmetaD simulation results of γ-Al 2 O 3 surface phosphorylation . . . . . A.1. Free enthalpy of orthophosphate adsorption and 31 P NMR at drying conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.2. Free enthalpy of orthophosphate coadsorption and 31 P NMR at drying conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.3. Free enthalpy of surface pyrophosphate formation and 31 P NMR at drying conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.4. Computational 31 P NMR reference . . . . . . . . . . . . . . . . . . . . . . vii First, I want to thank my thesis supervisors Carine Michel, Pascal Raybaud and Manuel

Figure 1 . 1 .:

 11 Figure 1.1.: Condensed phosphoric acid molecules: a) orthophosphoric acid, H 3 PO 4 ; b)

Fig. 1 .

 1 Fig. 1.2 gives an overview of the range of possible reaction types. The two reaction schemes summarize the range of reactions and intermediates for P -O bond formation and cleavage as a starting point to systematically analyze the whole phosphorylation reaction as a function of two bonds. Ref. [40] and [51] use a third axis to extend fig. 1.2 for proton transfer. However, by including protonation states of PO - 3 , X-and Y-separately,

Figure 1 . 2 .:

 12 Figure 1.2.: Phosphorylation reaction between sites X and Y, where vertical and horizontal directions refer to bond cleavage and formation with each site, respectively. Either concerted reaction or two types of intermediates can be expected: the addition complex with pentacoordinated phosphorus (top right) or elimination with a metaphosphate intermediate (bottom left). In strong alkaline environments (a), the reaction is only driven by formation and cleavage of phosphorus bonds. In neutral or acidic conditions (b), the reaction also depends on proton transfer reactions that have to be considered in every reaction path.

Figure 1 . 3 .:

 13 Figure 1.3.: Polyphosphate 31 P NMR chemical shift in aqueous solution depending on pH with respect to 85 % phosphoric acid. a: orthophosphoric acid, H 3 PO 4 (black line); b: pyrophosphoric acid, H 4 P 2 O 7 (red line); c: triphosphoric acid, H 5 P 3 O 10 (two blue lines).

Fig. 1 . 3 shows

 13 their experimental results on the three smallest polyphosphates: orthophosphate, pyrophosphate and triphosphate. The signal of orthophosphate is at 0 ppm at high acidity and for increasing pH values the signal is shifted to slightly more positive values up to 3 -5 ppm. This is related to an increased shielding of the deprotonated species. The pyrophosphate and triphosphate signals are clearly shifted to negative values, and triphosphate produces two signals. Pyrophosphate consists of two groups connected by oxygen, which produces a signal at -10 ppm at pH=3. Triphosphate H 5 P 3 O 10 is the first polyphosphate with one chain group connecting its two terminal groups: H 2 O 3 P -O -(HPO 2 ) -O -PO 3 H 2 . This leads to another signal for the central phosphate group at about -22 ppm (pH=3). The shift is related to a deshielding when substituting hydroxy groups in orthophosphate by one or two more electronegative phosphate groups. As for orthophosphate, signals by terminal and chain phosphate groups are gradually shifted to positive ppm values at higher pH as a consequence of deprotonation.

Fig. 1 .

 1 4 provides another graphical explanation of the angle and spinning operation. The spinning method is called magic angle spinning (MAS) and, only by including this spinning operation, well resolved solid state spectra are accessible. Bleam et al. conducted MAS 31 P NMR experiments on

Figure 1 . 4 .:

 14 Figure 1.4.: Probe spinning geometry in MAS NMR spectroscopy. The probe spins around a defined axis with respect to the magnetic field B 0 .

  1.1) vary over dramatic ranges of 50 to 100 ppm while the isotropic shift range of the presented minerals is only up to 20 ppm. The error estimates of isotropic chemical shifts (±0.1 ppm) are much lower than anisotropic components. 1.1.5. Surface-sensitive 27 Al and 31 P NMR For signal enhancement, cross-polarization (CP) experiments are performed, where an abundant nucleus (usually 1 H) is polarized, and which is brought into spin-thermal contact via dipolar coupling. DeCanio et al. recorded 27 Al and 31 P NMR spectra, with and without CP, of phosphate on γ-Al 2 O 3 after impregnation with an aqueous solution

Figure 1 . 5 .:

 15 Figure 1.5.: Direct polarization 27 Al-and 31 P-MAS NMR and cross polarization NMR of dried phosphate species on γ-Al 2 O 3 . Plot created with experimental data from DeCanio

Figure 1 . 6 .:

 16 Figure 1.6.: Unit cell content of the γ-Al 2 O 3 bulk model by Krokidis et al.[START_REF] Krokidis | Theoretical Study of the Dehydration Process of Boehmite to γ-Alumina[END_REF] . Oxygen and aluminium atoms are colored in red and blue. Cuts along the horizontal lead to (1 1 0) model facets (illustrated as plane and normal vector), whereas cuts along the vertical produce (1 0 0) model facets. Note that the Miller indices are with respect to the fcc lattice.

Figure 1 . 7 .:

 17 Figure 1.7.: Cutting a surface slab from a bulk unit cell (a). Several copies of the bulk unit cell produce a surface slab (b) with a thickness that depends on the bulk geometry.

Figure 1 . 8 .:

 18 Figure 1.8.: Cutting a surface slab along arbitrary surface directions in (a): (0 1 0) with small dashed line, (1 0 1) with dotted line, (1 1 1) with rounded dashed line. After choosing the Miller indices, the slab can be cut along different layers (b), which lead to different surface terminations.

Figure 1 . 9 .:

 19 Figure 1.9.: Shape of a γ-Al 2 O 3 platelet after topotactic transformation from boehmite.The (1 1 0) and (1 0 0) facets are most exposed. Miller indices are given with respect to the oxygen fcc lattice (often referred to as "spinel" symmetry).

Figure 1 . 10 .:

 110 Figure 1.10.: Surface free enthalpy of γ-Al 2 O 3 (1 0 0) (a) and γ-Al 2 O 3 (1 1 0) (b) for different degrees of hydration with coverage values θ given OH nm -2 . The inset graphs show the differential free enthalpy of H 2 O adsorption as a function of θ. This is the free enthalpy of each adsorbing H 2 O molecule in the sequence of full hydration. Figures are from Digne et al. in ref.[95].

2. 1 .

 1 Computational setup optimized projector evaluation in real space (LREAL = automatic). The Brillouin zone was sampled by Γ-centered Monkhorst-Pack grids of 2×3×1 and 2×2×1 on the (1 0 0) and (1 1 0) model facets, reflecting the different surface cell vectors and the aperiodic surface normal vector. These parameters were thoroughly benchmarked and lead to well converged structures and energies.

Figure 2 . 1 .:

 21 Figure 2.1.: Rational switching function to define a coordination collective variable.

2 . 1 .

 21 The cutoff distance in this function is given d 0 + r 0 , where r 0 stretches the overall function and d 0 is a radial shift. The exponents n and m change the stiffness and particularly influence the decay at long distances. For convenience, they are typically chosen such that m = 2n. Only this ensures that the sum d 0 + r 0

( a )Figure 2 . 2 .:

 a22 Figure 2.2.: Illustration of how a collective variable (CV) is biased by metadynamics. The free energy surface (FES) is based on an arbitrarily chosen polynomial (125 -45x + 35x 2 + 3x 3 -7x 4 + 0.1x 5 + 0.3x 6 ), whereas the biasing potential is a sum of gaussians fixed on 25 grid points along the graph.

  3.1ii the site b and b These labels are important throughout the whole thesis, because they are used to refer to individual adsorption sites where necessary.

( a )

 a Abstract structural representation involving local symmetry patterns and global symmetry features. Al atoms are in the center of their coordination octahedra (blue) and tetrahedra (green) and labelled with Greek letters when necessary. Oxygen sites from adsorbed water are labelled with red lowercase while alumina oxygen sites are labelled with blue upper case letters. Mirror planes are shown as thick black lines. Top view (iv, v, vi) and side view (vii, viii, ix) structures including the protonation states of all individual sites. Subsurface atoms are given in grey to complete the local coordination of the surface oxygen atoms.

Figure 3 . 1 .:

 31 Figure 3.1.: Three model reference surfaces of γ-Al 2 O 3 (1 0 0)[START_REF] Digne | Hydroxyl Groups on γ-Alumina Surfaces: A DFT Study[END_REF] (leftmost structures), n(1 1 0)[START_REF] Digne | Hydroxyl Groups on γ-Alumina Surfaces: A DFT Study[END_REF]

  from different model facets, because systematic errors in both reference calculations can be different, limiting the accuracy of this comparison. Consequently, free enthalpy comparisons between the (1 0 0) and (1 1 0) facet have to be discussed with more caution. It also adds an important aspect to the comparison of the unreconstructed and reconstructed structural references of (1 1 0), since for this facet it is desirable to use the same G[ ] in order to avoid this problem altogether. Wischert's reconstruction was identified in a study about low water coverage at high temperatures, which lead to the transition of the β-aluminium atom (see fig. 3.1 ii and iii) from an octahedral to a tetrahedral site. This structural change in fact reduces the number of available adsorption sites by one, hence the number of preadsorbed H 2 O on R(1 1 0) is reduced compared to n(1 1 0) as well. The reaction from one (1 1 0) structural reference [n] to another [R] is given by:

3 Figure 3 . 2 .:

 332 Figure 3.2.: Adsorption modes of orthophosphates on γ-Al 2 O 3 . Physisorption (ν 0 ) and several dentation patterns are possible: monodentates (ν 1 ), bidentates (ν 2 ) and tridentates (ν 3 ). All phosphate oxygen atoms are coordinated with at least one more atom other than phosphorus. The dentation sites are highlighted in red for substituted water adsorption sites, and blue for alumina oxygen sites.

22 . 2 = 168 .

 222168 x notation was extended by additional digits, such as ν xx . The superscript site information is comma-separated, e. g. ν bc,eB Even after elimination of several candidates by the first criterion for the adsorption of single orthophosphates, the number of adsorption combinations is overwhelming. The full scope of coadsorption includes all combinations of x and x in ν x +ν x . With values from 0 to 3, this gives 10 combinations from lowest to highest overall adsorptions, i. e. ν 00 to ν 33 . As shown in tab. 3.1, four adsorption modes lead to at least 30 adsorption candidates, after selecting a subset of structures from 51 initial chemisorption candidates. One also needs to add physisorption here, which adds one candidate for each structural reference, and thus 33 single adsorption cases for recombination to coadsorption and pyrophosphate adsorption. These 33 sites are ten on (1 0 0), ten on n(1 1 0) and 13 on R(1 1 0). The upper limit of combinations for a coadsorption is first estimated by calculating the number of inequivalent site combinations with n=10,10,13 n(n-1) This higher than the actual number of combinations since, for example on R(11 0), ν b 1 , ν bc 2 and ν bcB 3 all involve site b and thus cannot coexist. Nevertheless, by counting only possible combinations, 107 new coadsorption structures remain. This can be slightly reduced by excluding the less stable site µ a 2 -OH site on n(1 1 0) and R(1 1 0) leaving 86 candidates. The same considerations apply for pyrophosphates ν xx , with combinations for x and x . Several more combinations can be excluded due to the limited flexibility of the H 4 P 2 O 7 molecule, but still up to 73 structures would require testing. In the end this leaves us with 180 additional site combinations for coadsorption and pyrophosphates, which would each involve additional sampling for the hydrogen bond network (cf. section 3.2). We chose to restrict the number of combinations even further, only including those systems systematically, that involve one of the best single orthophosphate adsorption modes. These are combinations of 12 best adsorption geometries with one another, or if this is not possible, involving second best structures of each adsorption mode. Additional calculations were done for the more interesting cases. Overall 118 site combinations were tested. The tested sites are listed in tab. A.2 and A.3 (appendix, p. 156). Based on this exploration a fourth selection criterion was identified:4. If possible, coadsorbed phosphates are well distributed over the surface. That means coadsorbed phosphates maximize their distance to one another and avoid neighboring sites if possible. Identifying all important adsorption site candidates on (1 0 0), n(1 1 0) and R(1 1 0) started from single orthophosphate adsorbates, which were evaluated based on their free enthalpies. The best adsorption modes were used to study coadsorption and pyrophosphate, which gave additional insight into coverage effects. Overall 159 of these sites and site combinations were optimized. This way, selection criteria for new adsorption sites were developed, of which only the first is a

( 1 0 3 - 122 3 - 150 ( 3 -

 1312231503 0) ν 0 -126 (-232) ν a 1 -121 (-148) ν ab 2 -156 (-106) ν abc ν 0 -56 (-147) ν e 1 -105 (-121) ν de 2 -132 (-75) ν deA ν 0 -66 (-79) ν b 1 -120 (-62) ν eB 2 -160 (-30) ν bcB 168 (+41)

3 . 2 summarizesFigure 3 . 3 .:

 3233 Figure 3.3.: Gibbs free enthalpy (a) and enthalpy (b) of H 3 PO 4 adsorption per mole of phosphate at 120°C as a function of the number of released water molecules from the reference surfaces. For R(1 1 0), the common energy reference of the n(1 1 0) surface was used, hence adsorption occurs after removal of one H 2 O. Only the lowest enthalpy system of each category is included. Left panels: adsorption on the (1 0 0) surface. Right panels: adsorption on the (1 1 0) facets, n(1 1 0) with filled symbols, R(1 1 0) with empty symbols.

First, those energies

  are based on calculations of two different reference surfaces, neither of which is free of model errors. Second, both systems describe different coverages due to different areas of their surface unit cell. These two aspects are actually very complex to include into a quantitative model, since the coverages on both facets would differ significantly even for slightly different adsorption free enthalpies. It would require to calculate accurate adsorption energies for a range of coverages which would allow us to generate estimates for two local coverages based on the phosphate loading in the sample.

  fig. 3.4a). It is adsorbed as HPO 2 - 4 ion, in which the three free P -O oxygen atoms are stabilized by a hydrogen bond pocket consisting of six protons, i. e. two hydrogen bonds per phosphate oxygen. In contrast, the much less stable physisorption cases on n(1 1 0) and R(1 1 0) (cf. tab. 3.2) are H 2 PO - 4 ions, stabilized by hydrogen bonds from additional three protons (for two phosphate oxygen, i. e. 1.5 hydrogen bonds each). The high dentation cases n(1 1 0)ν deA 3

Figure 3 . 4 . 1 and µ Al O 1 are µ 1

 34111 Figure 3.4.: Ball-stick models of the most stable cases on n(1 1 0) and R(1 1 0) at 120°C and of (1 0 0) at 25°C. These species are some of the overall best orthophosphate adsorption cases, depending on environmental conditions. µ Al T 1 and µ Al O 1 are µ 1 sites from tetrahedral or octahedral oxygen atoms. Hydrogen atoms are white, oxygen atoms red, aluminium atoms blue and phosphorus atoms yellow. Hydrogen bonds are illustrated with dots.

3 and ν deA 3 ,and ν bcB 3 involveTable 3 . 3 .:

 3333 the trend is the same as for the reference, with n(1 1 0) being consistently 16 to 35 kJ mol -1 more stable than R(1 1 0). Only ν bcB 3 is stabilized by 46 kJ mol -1 on R(1 1 0). A second adsorption mode, R(1 1 0)ν eB 2 , is considerably more stable than any of the adsorption modes on n(1 1 0), despite the reconstruction. Both sites ν eB 2 the alumina oxygen site B which transforms from a µ 3 -O site on n(1 1 0) to a µ 2 -OH site on R(1 1 0). For illustration see fig. 3.1 and note that the reconstruction breaks the symmetry in n(1 1 0); thus b Comparison of bidentate on tridentate free enthalpies, ∆G in kJ mol -1 , of the same sites on n(1 1 0) and R(1 1 0) at drying conditions, 120°C and 2 kPa H 2 O partial pressure. Labels are given for R(1 1 0), on n(1 1 0) c is equivalent to b . It shows, that Wischert's reconstruction is a phenomenon of Digne's (1 1 0)

+ν a 1 and n(1 1 0)ν bc 2 +ν de 2 and ν bcB 3 .

 123 3.3 and tab. 3.4 are not consistent in this regard. In fact, only (1 0 0)ν ab 2 specifically fit this expectation. On (1 0 0), the adsorption of two bidentates is only possible by one pattern ν aa 2 +ν bc 2 that excludes the most stable single bidentate ν ab 2 , since the remaining sites a and c are not direct neighbors. The free adsorption energy of the monodentate R(1 1 0)ν d 1 is ∆G = -86 kJ mol -1 , which makes it very unfavorable on R(1 1 0); ν c 1 (-120 kJ mol -1 ) and ν e 1 (-98 kJ mol -1 ) are better alternatives. Nevertheless, ν d 1 is the most stable candidate as coadsorbate with ν eB 2 and ν bcB 3 , where in each case better sites would be still available. This is likely related to the hydrogen bonds formed with those sites to stabilize ν eB 2 The ν d 1 adsorption isfor construction of the most favorable site -the best choice, as it does not destabilize the other site. The structures comprise different hydrogen bond patterns but the differences are subtle.

Figure 3 . 5 .:

 35 Figure 3.5.: The most stable pyrophosphate species, R(1 1 0)ν bcB,e 31

Fig. 4 .

 4 Fig. 4.1 shows NMR spectra acquired by D. Wisser and A. Lesage on γ-Al 2 O 3 samples provided by IFPEN based on Dynamic Nuclear Polarization Surface Enhanced NMR Spectroscopy (DNP SENS) techniques. [18] The 31 P CP NMR signal in fig. 4.1a confirms the trend observed by DeCanio et al.: an almost featureless signal ranging from 0 to -20 ppm at low coverage. With increasing coverage, the signal broadens and its maximum is shifted to lower values. At the highest experimentally tested coverage of 4.1 P nm -2 , the signal band covers an area with an upfield shift to -30 ppm. For the first time they detected the presence of pyrophosphates and other short-chained polyphosphates by using double and triple quantum correlation spectra (fig. 4.1b and 4.1c). The 31 Al -27 Al-INEPT NMR experiments in fig. 4.1d showed clearly that the phosphates are chemisorbed to the γ-Al 2 O 3 surfaces by P -O -Al bonds. However, it is surprising that this apparently

Figure 4 . 1 .

 41 Figure 4.1.: Different DNP enhanced NMR techniques recorded at 9.4 T (400 MHz 1 H resonance, 263 GHz microwave frequency). The 31 P NMR signal (a) is shifted with the overall phosphate coverage. The results of the 1D INADEQUATE experiments (b) and (c) provide the evidence for the presence of condensed phosphates. Vertical line follows the signal maximum at a coverage of 2.8 P nm -2 . The 2D 31 P -27 Al INEPT correlation experiments (d) lead to new conclusions regarding the nature of P -O -Al bonds with the surface (discussed in main text). NMR spectra obtained at CRMN by D. Wisser and A.

Figure 4 . 2 .

 42 Figure 4.2.: (a) Plot of the relationship between chemical shift and free enthalpy for single adsorbed phosphates (1.5-2.2 P nm -2 ). (b) Structures of the four most stable adsorption sites. Oxygen sites from H 2 O adsorption are labelled in red and alumina oxygen sites in blue.

31 P

 31 NMR spectrum. The best candidate is the tridentate R(1 1 0)ν 3 with a chemical shift of -2 ppm, which is at the most positive end of the experimental range. The other tridentate, n(1 1 0)ν 3 , leads to a chemical shift of -11 ppm, just close to the band maximum of the lowest coverage in fig. 4.1a. The bidentates R(1 1 0)ν 2 and n(1 1 0)ν 2 are both deshielded compared to their tridentate counterparts. Hence, if the amount of bidentate increases, the overall signal would be shifted to lower values. The most stable monodentates (ν 1 ) and physisorption cases (ν 0 ) are 20 kJ mol -1 less stable than the most favorable site on each reference structure, and in most cases not in the chemical shift experimental range. The fact that the most stable adsorption structures are in agreement with the experiment is a clear validation of our computational method, and allows us to analyze further trends. The next section takes a more elaborate view on 31 P chemical shift for coadsorption.

Figure 4 . 3 .

 43 Figure 4.3.: (a) Plot between chemical shift and free enthalpy of orthophosphate coadsorption (3.0-4.3 P nm -2 ). Colors and symbols refer to the reference surfaces (cf. fig 3.1) and adsorption modes, respectively. The chemical shifts of two 31 P nucleus in the same system are connected by a horizontal line. (b) Chemical environment of the three most stable adsorption sites combinations.

Figure 4 . 4 .

 44 Figure 4.4.: (a) Coverage effect on the 31 P chemical shift for six selected sites, i. e. the two most stable adsorption sites of each reference structure: (1 0 0)ν a 1 , (1 0 0)ν ab 2 , n(1 1 0)ν de 2 , n(1 1 0)ν deA 3 , R(1 1 0)ν eB 2 and R(1 1 0)ν bcB 2 (for labels cf. fig. 3.1a, p. 53). The phosphate coverages are given as ranges, because the two model facets have different surface areas.

4 . 2 )

 42 , this aspect was not further examined.The five considered bidentates (ν 2 ) and tridentates (ν 3 ) lead to fairly separated chemical shifts all over the range. The chemical shift of each site does not change systematically with increased coverage but might be different within a range of 5-10 ppm depending on their surrounding. The chemical shift difference between the adsorption sites is in comparison much more significant. For all bidentates (ν 2 ) and tridentates (ν 3 ), the signals are in distinct areas, even at high coverage where overlapping ranges are overall in the same order. Combined with the stabilization of adsorption on n(1 1 0) sites, this would introduce a significant change of chemical shift to negative values for increased coverage, if these sites become more populated. At low coverage, the predominant adsorption species are found at reconstructed R(1 1 0) surface sites. With increasing the coverage, contributions of adsorption sites identified on the n(1 1 0) structural references become more and more important.

31 and R(1 1 0)ν bc,eB 22 (

 3122 fig. 3.1a, p. 53), share the same oxygen surface sites and have almost the same free enthalpy. Their calculated signals are in a range from -16 to -27 ppm (cf. fig. 4.5). The third most stable pyrophosphate is on the (1 0 0) facet. The chemical shift of both groups is at about -10 ppm, which is in agreement with the experimental range. Of all calculated chemical shifts on the (1 0 0) facet, pyrophosphates are the only results that align well with the experiment. Hence it is possible that pyrophosphates are more present than

Figure 4 . 5 .

 45 Figure 4.5.: (a) Correlation between chemical shift and stability of pyrophosphates (1.5-2.2 P nm -2 ). The chemical shift of two 31 P nuclei of the pyrophosphate are connected by a horizontal line. (b) Chemical environment of the three most stable adsorption sites.

31 P 31 P

 3131 P NMR chemical shifts from 0 to -15. The computed pyrophosphate model systems complement the calculated 31 P chemical shifts for orthophosphate by covering a 31 P NMR signal range from -10 to -28 ppm. Several different aspects were tested that change with coverage. The preferred adsorption mode shifts from tridentates (ν 3 ) to bidentates (ν 2 ) at increased coverage, given the change in thermodynamic stability and simultaneously in agreement with the downfield shift of the 31 P NMR signal. Fig. 4.1d clearly shows the site dependence of chemical shifts, whereas the effect of different phosphates and protonation networks does not change the 31 P chemical shift in a clear direction. Pyrophosphates and longer polyphosphate chains were for the first time directly detected by our collaborators' INADEQUATE double and triple quantum correlation experiments. The calculated chemical shifts of orthophosphate on the (1 0 0) facet are in most cases at considerably more positive values, indicating that the facet is almost negligible compared to the (1 1 0)facet. Pyrophosphates models on (1 0 0) are in a much better agreement, which is expected due to the downfield chemical shift of condensed phosphates. Similar to the cases on the (1 0 0) facet, the chemical shifts of physisorbed orthophosphates do not agree well with the experiments. While the difference between through-space and through-bond[START_REF] Luis | Hydrolysis rates of inorganic polyphosphates in aqueous solution as well as in soils and effects on P availability[END_REF] P -27 Al INEPT spectra, is evidence for nonbonded phosphate sites, rather than physisorbed phosphates, those are more likely dangling groups of ν x0 bonded pyrophosphates.

  the other hand the reaction to form monodentates (ν 1 ) involves the substitution of one adsorbed H 2 O molecule as discussed in section 3.1 (p. 52). This can occur a second time, in which a free P -O oxygen in the monodentate substitutes another preadsorbed H 2 O.The preadsorbed water sites are either µ 1 -H 2 O or µ 1 -OH, the latter of which has to be protonated before desorption. Unlike in the drying state, this desorbed water molecules become part of the solution, hence the entropy effect is expected to be significantly lower. The H 2 O abstraction mechanism to form tridentates (ν 3 ) on the (1 1 0) facet is considerable more complex, since the formation of another bond with the γ-Al 2 O 3 surface now involves a triangular pattern of surface oxygen sites. On the (1 1 0) facet this always involves at least one alumina oxygen atom.The next section demonstrates the possibility of predicting barriers with WTmetaD simulations by computationally reproducing the experimentally determined orthophosphate condensation barrier. Then the focus is on the first chemisorption steps for the formation of mono-and bidentates (ν 1 , ν 2 ) from physisorbed orthophosphate (ν 0 ). At last the question whether tridentates (ν 3 ) can form at the solid liquid interface is addressed.

Figure 5 . 1 .:

 51 Figure 5.1.: Schematic of the orthophosphate condensation of phosphoric acid in aqueous solution. H 2 PO - 4 is more stable than H 3 PO 4 . The coordinates of the different structures of each path in the basis of CN(P -O all ) and CN(P -O bridge ) (see main text) on x and y-axis are given as 2D vectors (x, y).

Figure 5 . 2 . 3 (

 523 Figure 5.2.: 2D WTmetaD simulation of two H 3 PO 4 (2 P) condensation to H 4 P 2 O 7 (P 2 ) in a box with 124 H 2 O molecules. (a) Optimized paths of H 2 O elimination to form PO - 3 (metaphosphate, mP) followed by H 2 PO - 4 addition. (b) H 2 PO - 4 addition to PO - 3 . (c) The free energy surface is defined by the number of bridging oxygen CN(P -O bridge ) and the overall number of phosphorus-coordinating oxygen CN(P -O all ) The well in the top left corner of the graph at (0.1, 7.6) represents the situation with two H 3 PO 4 with a total of 8 oxygen binding to phosphorus but none bridging. In the bottom right, a well appears around (1.0, 6.8) which shows a total of 7 phosphorus binding oxygen and 1 bridging as in H 2 P 2 O 2 - 7 . (d) The number of H + in the least protonated phosphate, 0 (purple), 1 (dark green), 2 (bright green), 3 (yellow).

  the two relevant aluminium atoms (with γ-label in 3.1ii on page 53) with water oxygen (Al γ -O b w ), and CN p , number of phosphate oxygen atoms binding to the same site b (Al γ -O b p ). Since the H 3 PO 4 adsorption requires the initial desorption of H 2 O from b, there is only a single possible pathway to form the monodentate ν b 1 . This monodentate can, again, only form the bidentate ν bb 2 after desorption of the second H 2 O from site b . Both, the reactions ν 0 -→ ν b 1 and ν b 1 -→ ν bb 2 decrease CN w by 1 and increase CN p by 1.

Fig. 5 .

 5 Fig. 5.3 shows the results of a one-dimensional well-tempered metadynamics simulation. At the end of the simulation, the bidentate (ν bb 2 , structure I in fig. 5.3) is clearly the most stable state followed by monodentate (ν b 1 , structure III) and physisorption (ν 0 , structure V). The free energies of ν b 1 and ν bb 2 formation are, -58 and -88 kJ mol -1 , respectively, relative to the physisorbed state ν 0 . The first and second adsorption barrier

Figure 5 . 3 .

 53 Figure 5.3.: (a) 1D WTmetaD simulation of phosphate adsorption to γ-Al 2 O 3 by H 2 O substitution at the (1 1 0)bb' sites. The CV is the difference of coordination numbers of Al -O with phosphate oxygen and water oxygen: CN w -CN p (a) Consecutive FES in steps of 1000 gaussian hills each (hair lines), the last FES is highlighted as thick line. (b) Projection of trajectories on CN w +CN p . This reveals the progression of the simulation along CN w and CN p (equivalent to diagonals, north-east: CN w , north-west: CN p ). Four different colors indicate each walker (section 2.2.5), p. 48). (c) Structures at minima and transition states. Both aluminium atoms are octahedral sites. Gaussian hills: σ = 0.04, W 0 = 3 kJ mol -1 ; deposition rate: 100 ps -1 .

5 . 3 . 2 . Tridentate ν deA 3 on n(1 1 0 )

 53230 fig.3.1iii (p. 53). The formation of this species involves the alumina oxygen site µ B 2 -OH, that needs to be phosphorylated in order to form this adsorption structure. Similarly

(a) ν de 2 3 + H 2 O 1 Figure 5 . 5 .

 232155 Figure 5.5.: Two 2D WTmetaD simulations for formation of ν deA 3 on n(1 1 0). Bias along coordination number of P -O alumina (CN a ), P -O water (CN w ) and P -O phosphate (CN p ). Isolines in free energy surface in kJ mol -1 . (a): Minimum energy path of the phosphorylation reaction of alumina oxygen "A" (cf. fig. 3.1ii, p. 53), starting from the bidentate ν de 2 with H 2 O as leaving group. (b): Free energy surface of (a). (c): Minimum energy path of the same phosphorylation with a ν 20 pyrophosphate, such that the leaving group is orthophosphate.

B 2 - 3 + H 2 OFigure 5 . 7 .

 23257 Figure 5.7.: Two 3D WTmetaD simulations as candidate reactions to form ν bcB 3 on R(1 1 0). Bias along coordination number of P -O alumina (CN a ), P -O water (CN w ), P -O phosphate (CN p ) and µ B 2 -O -H (CN h ). Iso-lines in free energy surface in kJ mol -1 . (a): Minimum energy path of the phosphorylation of alumina oxygen "B" (cf. fig. 3.1ii, p. 53), starting from the bidentate ν bc 2 with H 2 O as leaving group. (b): Free energy surface of (a). (c): Free energy surface of proton transfer in (a). (d): Minimum energy path of the same phosphorylation with a ν 20 pyrophosphate, such that the leaving group is orthophosphate. Path incomplete in order to include correct transition state (software limitation). (e): Free energy surface of (d). (f): Free energy surface of proton transfer in (d). See eq. 5.5 and 5.6 on page 102 for the reaction mechanism.

and ν bcB 3 fromFigure 5 . 8 .

 358 Figure 5.8.: Two 3D WTmetaD simulations as candidate reactions to form ν eB 2 on R(1 1 0). Bias along coordination number of P -O alumina (CN a ), P -O water (CN w ), P -O phosphate (CN p ) and µ B 2 -O -H (CN h ). Iso-lines in free energy surface in kJ mol -1 . (a): Minimum energy path of the phosphorylation of alumina oxygen "B" (cf. fig. 3.1ii, p. 53), starting from the bidentate ν e 2 with H 2 O as leaving group. (b): Free energy surface of (a). (c): Free energy surface of proton transfer in (a). (d): Minimum energy path of the same phosphorylation with a ν 20 pyrophosphate, such that the leaving group is orthophosphate. (e): Free energy surface of (d). (f): Free energy surface of proton transfer in (d).

deA 3 (

 3 2 O 3 (1 1 0) that were identified as the most stable model adsorption structures at drying conditions: R(1 cf. fig. 3.1, p. 53). The idea was a comprehensive study of surface phosphorylation reactions targeting the aforementioned most stable of orthophosphate species. The precursor of these reactions are orthophosphates at the next lower Table 5.1.: WTmetaD results of the six conducted explorations for γ-Al 2 O 3 surface phosphorylation: Reaction free energies and free energy barriers by the end of the simulations and the number of included CV and simulation size in hills are given. The leaving groups

  2 indeed leads to the most relevant species during impregnation: bidentates and monodentates. Even though the three-dimensional WTmetaD simulations on R(1 1 0) (fig. 5.7 and fig. fig:reB) are not accurately converged, they show trends regarding the site deprotonation of µ B 2 -OH: the deprotonation typically occurs after phosphorylation, which means the deprotonation does not occur readily by itself.

  P NMR, typically only single parameters like average bond distance or bond angles were correlated with an NMR chemical shift. Dawson and Asbrook unravelled the combined effect P -O bond distances and P -O -Al angles in aluminium phosphates by screening calculations at DFT level and GIPAW chemical shift simulation and multivariate regression. Their results demonstrated that the NMR prediction depends on multiple structural features, of which P -O bonds and P -O -Al angles are the leading contributions. The authors proposed several regression functions with mean absolute errors (MAE) of 0.49 ppm -much more accurate than simpler correlations of 31 P NMR with either P -O bonds (quadratic fit: MAE = 3.9 ppm) or P -O -Al angles (linear fit: MAE = 7.1 ppm)

  Na and 27 Al nuclei in aluminosilica glasses based on DFT geometry data and GIPAW calculations. The authors combined three atom centered descriptors with seven common machine learning algorithms and found that the SOAP descriptor (smooth overlap of atomic positions) in combination with linear ridge regression (LRR) produces accurate results for their purpose ( 29 Si: RMSE = 1.3; 17 O: RMSE = 1.5; 23 Na: RMSE = 1.5; 27 Al: RMSE = 1.5) This chapter applies the methodology of Chaker et al. to the calculated 31 P chemical shifts results from the database generated in chapter 3. As starting point prediction schemes based on SOAP/LRR were used to generate a prediction model for 31 P chemical shifts. This work was only possible in collaboration with Tao Jiang, Research Engineer at ENSL, who carried out the machine learning parametrizations and helped during evaluation.

6 .

 6 system. It is straightforward to derive a local variant of MBTR (LMBTR) by defining a 121 Machine learning NMR chemical shifts central atom and using a subset of ranges: g 1 is removed, since stoichiometry is difficult to establish as a local system property, and only g 2 and g 3 involving the central atom are included. The central atom is typically treated separately, independent from the element. Since the NMR chemical shift is a local property of the chemical shielding by the local electron density around a nucleus, local descriptors are a natural choice. A popular local

  the γ-Al 2 O 3 model surfaces were used to train different machine learning algorithms to predict the 31 P chemical shift based on the geometry. The original set of GIPAW calculations only included the most favorable cases of the sampled hydrogen bond networks described in section 3.2 (p. 58), which were 361 data points in total including (a) RMSE comparison of smaller data set which only includes the most favorable adsorption sites (361 data points) and the full data set (1227 data points). (b) Randomized fractions of large data set. RMSE of training, validation and test set.

Figure 6 . 1 .

 61 Figure 6.1.: RMSE change with respect to data set size, in case of the use of on SOAP/LRR. (a) Scaled histograms of RMSE (x-axis) for small an large data set; y-axis in arbitrary units. (b) RMSE (y-axis) dependent on the data set size (x-axis).

( a ) 1 Figure 6 . 2 .: 1 Figure 6 . 3 .:

 a162163 Figure 6.2.: Test of the SOAP descriptor against three different machine learning algorithms for (a) LRR, (b) MLP, (c) XG-Boost; rcut and lmax/nmax see the text. Comparison of training data (left), validation (middle) and test data (right). Color coded range describes the RMSE, ranges defined individual for each graph. Validation and test RMSE were higher for MLP and XGBoost.

( a )Figure 6 . 5 .:

 a65 Figure 6.5.: Performance of machine learning result with respect to GIPAW training data. LMBTR/LRR produces more accurate predictions, than SOAP/LRR. Value for λ has been optimized separately.

Figure 6 . 6 .

 66 Figure 6.6.: LMBTR/LRR coefficients between the central phosphorus atom and distances for all elements in the system.

  P chemical shifts allows for a different pathway of evaluating the chemical shifts at drying conditions. By evaluating the chemical shift in whole AIMD trajectories, it might be possible to quantify the contributions of each adsorption mode individually based on the experimental band instead of the model free enthalpy estimates.

  polyphosphate condensation are crucial to rationalize the change of chemical shift for different phosphate loadings.

7 )

 7 According to the Hohenberg and Kohn theorem, in principle, this relationship can be reversed, i. e. |ψ e (r e ) = |ψ e [ρ 0 ] , where the brackets indicate that the wave function is a functional of the ground state density. If the wave function is fully determined by its density, the Schrödinger equation can be expressed as density functional as well:E el = ψ[ρ 0 ]| Ĥe |ψ e [ρ 0 ] = E el [ρ 0 ] (A.8)This implies that for most cases should be possible to replace a many body wave function, a complex 3N-dimensional construction comprised of densities for every particle N, it contains by one single 3-dimensional density distribution. Consequently calculations based on the density are much less time consuming, and, regarding numerical stability for convergence, much more robust. Most of the terms shown in eq. A.4 and eq A.5 are

V

  s (r) = V N (r) + e 2 ρ(r ) |r -r | dr + V xc (r) (A.14) Starting from an initial guess, the density and orbitals are iteratively optimized, and new densities produced by eq. A.7. Kohn-Sham DFT is a compromise between wave function method and pure DFT. While the wave function based kinetic energy introduces a single reference determinant, it allows for an accurate estimate of the kinetic energy; hence the correction added in V xc (r) is smaller and would therefore introduce smaller inaccuracies. The mix between both sides sides also leads to a clear computational advantage, since the Vee part could be substituted by a much smaller integral V ee [ρ 0 ]. Its asymptotic computation times scales with O(n 3 ) for a system of size n. Several approximations are available for the exchange-correlation potential which can be ordered by a dependence on the local density, its first derivatives, and further corrective measures: the local density approximation (LDA) only includes terms that scale with the density, while general gradient approximations (GGAs) include density changes and meta-GGAs add kinetic energy corrections to V xc [ρ].

Table 1 .1.: 31

 1 P NMR isotropic chemical shifts (δ iso ) and principal components of chemical shift tensors (δ ii , i = 1, 2, 3) and the anisotropy ∆δ for several aluminium phosphate minerals. ∆δ = δ 33 -δ 11 , δ 11 < δ 22 < δ 33 .Table and caption (rephrased)

Table 3 .

 3 1.: 51 adsorption site candidates on the three surface models. The site labels in superscript are given in fig. 3.1i, ii and iii. The full structure exploration included all possible adsorption combinations based on neighboring sites. The struck through structures are excluded later on, since two phosphate oxygen atoms bind to the same aluminium atom.

	Facet	Monodentate	Bidentate	Tridentate

  What followed based on these considerations was the systematic listing of all possible sites. Even after excluding symmetry-equivalent cases, this left at least 51 chemisorbed candidate sites which are tabulated in tab. 3.1. Physisorption (ν 0 ) was handled separately by positional screening of the surfaces, leading up to 16 test structures (including hydrogen bond optimization, see section 3.2) for each physisorption result. Without further knowledge, it is difficult to reduce this number. It should be noted that for each site, sampling the hydrogen bond network are necessary, leading to hundreds of calculations, as discussed in section 3.2. Moreover, for the purpose of processing the results and extending this data base with relevant candidates, further selection criteria are needed. A full list of free enthalpies is given in tab. A.1 (appendix, p. 154). After

	2
	and ν abc 3 , respectively. Due to symmetry, the overall number of adsorption sites can be
	slightly reduced, which if reflected in the labels, e. g. b and b on n(1 1 0) are equivalent.

are not in direct view of fig.

3

.1i due to periodic boundary conditions strictly, ignoring "distant" neighbors with distances larger than 3Å. Since these patterns, in particular the triangular patterns are not easily attained, alumina oxygen sites (cf. blue colored surface sites in fig.

3

.1) were included in the construction of several higher dentation modes. A known reaction mechanism to form such adsorption species is the phosphorylation, as discussed in section 1.

1.3 (p. 18)

. Now, the systematic exploration of all possible adsorption is a major challenge, because despite four distinct adsorption modes, the three structural references provide a considerable number of adsorption sites. In order to constrain the number of explored structures to the most relevant cases, selection criteria were identified and applied for next explored structures. Overall the three model surfaces are comprised of 15 water adsorption sites and 18 alumina oxygen sites. In order to keep track of every tested adsorption site, the ν x notation is extended by the site with the labels given in fig. 3.1a. For example, on (1 0 0) the monodentate, bidentate and tridentate adsorption modes are labelled ν a 1 , ν ab

Table 3 .2.:

 3 Most favorable adsorption sites and free enthalpies ∆G and enthalpies ∆H in kJ mol -1 at drying conditions, 120°C and 2 kPa H 2 O partial pressure.

	Facet	Physisorption	Monodentate	Bidentate	Tridentate
		∆G (∆H)	∆G (∆H)	∆G (∆H)	∆G (∆H)

  1, appendix, p. 154). It is also not the case for higher coverages, such as the calculated coadsorption and pyrophosphate adsorption cases. For a general overview in fig. 3.3, ∆G and ∆H were plotted against the number of released H 2 O. With eq. 3.3 and 3.5-3.9, the number of released H 2 O

varies from 0 to 7 relative to the water-covered reference structures. Coadsorption is less favorable than single orthophosphates, which means that phosphate adsorbates compete for surface area, and there is no cooperative adsorption found. This is consistent with the observation, that structures with maximum phosphate distance minimize the free enthalpy. Pyrophosphates are generally less stable even though there are some cases in which pyrophosphates are slightly better stabilized. At drying conditions, up to three or, respectively, four H 2 O are substituted on the (1 0 0) and (1 1 0) facets. Substitution of further H 2 O is endergonic at 120°C, but becomes exergonic at higher temperatures (for more data, cf. tab. A.2, appendix, p. 156).

Table 3 .4.:

 3 Free enthalpies of the most favorable adsorption mode for coadsorption and pyrophosphates in comparison to the nextmost stable case with different overall dentation x+x in ν x +ν x and ν xx , sorted by dentation. ∆G in kJ mol -1 at drying conditions, 120°C and 2 kPa H 2 O partial pressure. ν e,b11 is a pyrophosphate species that binds with its P -O -P bond positioned over alumina oxygen sites A and B, leaving boundaries of unit cell depictions in fig.3.1ab.

	Facet	Orthophosphate coadsorption	Pyrophosphate adsorption
		sites	∆G	sites	∆G	sites	∆G	sites	∆G
	(1 0 0)	ν ab 2 +ν a 1 -119 ν aa 2 +ν bc 2	-116	ν a,a 11	-99	ν bc,a 21	-106
	n(1 1 0)	ν bc 2 +ν de 2	-127 ν bcE 3 +ν de 2	-128	ν e,b * 11	-100	ν de,b 21	-97
	R(1 1 0)	ν eB 2 +ν d 1 -132 ν bcB 3 +ν d 1 -121	ν eB,b 21	-105	ν bcB,e 31	-117

*

  Al -31 P-INEPT experiments (left side of the 2D graphs in fig. 4.1d) reveals that the signal bands deviate significantly, which indicates that the down-field part of the overall signal is represented by physisorbed phosphates or dangling groups of terminally binding polyphosphates. A more subtle observation is that also the 27 Al signals are slightly shifted to positive values for increasing phosphate loading. As previous work and these results show, many relevant aspects of phosphate adsorption are accessible

  The 27 Al -[START_REF] Luis | Hydrolysis rates of inorganic polyphosphates in aqueous solution as well as in soils and effects on P availability[END_REF] P-INEPT NMR experiments (see fig. 4.1d) showed that a significant part of phosphates does not directly bind to the alumina surface. The through-bond INEPT signal projected for 31 P nucleus, is shifted by around -6 ppm relative to the through-space INEPT signal, which shows that different phosphate species are involved in both cases.

	The best candidates in for those nonbonded cases are
	• physisorbed phosphates, ν 0 and
	• terminal pyrophosphates, ν x0 .
	The calculated 31 P chemical shifts for physisorption are in most cases at significantly
	more positive values than the experimental range. The few exception that match the
	experiment are significantly less stable than the former. In contrast both computed 31 P

NMR chemical shifts for any of the ν x0 are in agreement with experiment. While it might still be possible that physisorbed phosphates are misrepresented due to their more dynamic behaviour (without covalent anchors they are more mobile), from the current data, pyrophosphates provide a good rationale for these experimental findings.

  6.2 contains the results for the tests of SOAP parameters. From the ranges of the RMSE, it becomes clear that the LRR produces the most accurate predictions. Training and validation graphs of MLP and XGBoost contain notable noise, which indicates that both methods need more data points in order to work reliably.

  who found that LRR works better than any other tested algorithm for

	chemical shift estimations. Their data also indicate that similar algorithms like gaussian
	kernel ridge regression (GKRR) and elastic net regression (ENR) produce similar results
	with SOAP, while different local descriptors with LRR, like Behler-Parinello symmetry
	functions (BPSF) and the angular and radial distribution function (ARDF), were clearly
	inferior to SOAP. However, our results show that the prediction accuracy is improved by
	replacing SOAP/LRR with LMBTR/LRR. These trends imply that the limitation in the

  90 % of all data for training and 10 % for testing, thus test and training RMSE are more similar than in section 6.3. Due to overfitting, the training RMSE is still lower than the validation/test RMSE. As shown in fig. 6.1b, the training RMSE and validation/test RMSE asymptotically converge with the data set size. While it is certainly possible to extrapolate the RMSE trend, at this stage we can assume that the training RMSE as optimistic and test RMSE as conservative estimate of the accuracy of the same method for a much larger data set.The two values provide a lower and upper boundary to the real RMSE. Since the test RMSE value of LMBTR/LRR is lower than the training RMSE of SOAP/LRR, we can conclude that LMBTR/LRR is a better prediction scheme than SOAP/LRR. This trend is consistent for all the tested machine learning algorithms, which demonstrates that the LMBTR descriptor is more suitable for predicting[START_REF] Luis | Hydrolysis rates of inorganic polyphosphates in aqueous solution as well as in soils and effects on P availability[END_REF] P NMR chemical shifts. For chemical shift estimates, the descriptor which provides the most relevant information should perform best. It is not trivial to compare both descriptors based on their features,

  6.6 is a plot of the radial part of the numerical grid. The P -O distance around 1.5Å reflects P -O bonds. The smallest P -O bonds are correlated with deshielding which is due to few significantly distorted phosphate adsorbates with free P -O bonds. They are actually too unstable to be relevant here, and should thus be excluded from the data set. At slightly higher values, the P -O bonds lead to slightly positive contributions from P -O -X. The large values between 1.7 and 2.0Å are again Al IV and 3.5Å for P -Al VI , the angles approach 180°which leads to deshielding of phosphorus due to bond competition in the O-p σ orbital of the phosphate (electron pulling of Lewis acidic moieties is not effective at small angles). At further distances the effects are overall more subtle except for the P -O distance, which leads to a deshielding from 3.1 to 4.5Å.

	exceptionally distorted geometries. For the smallest P -H, P -Al and P -P distances,
	another positive correlation is observed, which is followed by an even stronger negative
	correlation at a slightly larger distance in case of P -Al and P -P. This is related to the
	P -O -(P, Al) angles that are directly connected with P -(P, Al) distances. The lowest
	P -Al and P -P distances at 2.5 and 2.7Å minimize P -O -Al and P -O -P angles, which
	leads to high amounts of electron density around phosphorus. At larger distances, which
	are 3.1Å for P -P, 3.0Å for P -

  but this largely depends on the descriptors and data set size. The comprehensive tests byChaker et al. indicate that the choice of a good descriptor, in their case SOAP, leads to the biggest improvement in NMR prediction accuracy. The simplest machine learning method, LRR, produced the best result. The success of LRR is overall surprising, but it might mean that better descriptors have to be identified before more sophisticated machine learning methods can improve the prediction model. In this study, it was demonstrated that with LMBTR, a clear improvement is observed for all machine learning methods tested. One rationale for this improvement is that LMBTR describes angles and distances independently.

  formation. Our experimental collaborators recently produced evidence for the presence of polyphosphates at dry surfaces. Pyrophosphate adsorbates can actually fill the gap from -15 to -30 ppm based on the presented DFT models; however, based on the free enthalpy estimates, they might not be stable enough. Moreover, the presence of terminally bound pyrophosphates can also rationalize discrepancies in two-dimensional

31

P -27 Al NMR INEPT experiments, which show two different distributions of phosphate species for through-bond and through-space correlation. The alternative rationale for this discrepancy are physisorbed (ν 0 ) orthophosphates, although those are much less stable than the chemisorbed species on (1 1 0). Aluminium phosphate formation, which was discussed earlier by Van Eck et al.

  MAS NMR results. Those experiments were done with samples dried at 120°C but they were kept at room temperature in contact with air. Hence, it was an open question which temperature is the better reference. In fact, the most stable adsorption structure identified was (1 0 0)ν 0 at room temperature, but this structure is on a minority facet and does not reproduce the 31 P chemical shift (chapter 4). Only by assuming the thermodynamic stability ordering at the drying conditions themselves, free enthalpy estimates and experimental results are consistent. The barriers identified in chapter 5 rationalize why this is the case: at room temperature, surface hydration of P -O -Al bonds is kinetically inhibited.

31 P

The last chapter, chapter 6, demonstrates that it is possible to generate efficient prediction schemes of

[START_REF] Luis | Hydrolysis rates of inorganic polyphosphates in aqueous solution as well as in soils and effects on P availability[END_REF]

P chemical shifts based on atomic positions based on linear ridge regression, a simple machine learning regression algorithm. This study was only possible due to the large structural exploration carried out in chapter 3 (p. 51), which allowed to generate 1227 data points for 31 P GIPAW chemical shifts as the basis for this work. To the best of our knowledge, this work is the first application of machine learning on 31 P NMR of phosphate surface species; nevertheless, several previous studies dealt with bulk materials

  Al NMR in order to computationally reproduce the two-dimensional INEPT correlation NMR experiments. Overall, the discussed structures are the first realistic adsorption models available for the γ-Al 2 O 3 /H 3 PO 4 /H 2 O system at the time, and as such they should be used as the basis for future DFT studies of heterogeneous catalyst models involving γ-Al 2 O 3 , phosphate additives and other active compounds. Their plausibility was demonstrated by comparison with 31 P NMR experimental results (chapter 4, p. 75), with a perspective on further quantification of the spectra by generating a machine learning NMR prediction scheme and using it to test NMR chemical shifts on trajectories of molecular dynamics. The metadynamics simulations presented in chapter 5 (p. 91) allowed to disentangle which adsorption structures occurred at which of the two tested conditions, dried or wet surface. In the context of future research, this work vastly improves our understanding of the effects of phosphates at γ-Al 2 O 3 surfaces which paves the way to construct new, more holistic models of heterogeneous catalysts.

4 )

 4 Since the Coulomb terms are sums over particle pairs, they rapidly grow with the number of particles. If we are furthermore interested in dynamic information of, i. e. structural changes over time, we would have to solve the time-dependent Schrödinger equation. As referred to as Hartree-Fock method, does not directly produce solutions to a ground state wave functions but is capable of iteratively improving a trial wave function by optimizing the one-electron states or orbitals in the average electric field of the other orbitals. Hence it is also called self-consistent field (SCF) method. The most critical approximation here stems from the fact that each electron orbital is expressed in the mean-field of the other electrons. In consequence, electron correlation phenomena other than covalent bonds, like van-der-Waals interaction cannot be reproduced. To alleviate this problem, many improvements were proposed, for example based on perturbation theory, where the Hartree-Fock solution which is exact for single electrons is improved by adding perturbation terms to approximate the neglected Coulomb correlation effects. Most noteworthy is Møller-Plesset perturbation method since it can efficiently produce accurate

Table A .

 A 2.: Average free enthalpy of adsorption of two orthophosphates in kJ mol -1 and 31 P NMR shift in ppm with respect to to 85 % H 3 PO 4 for the (at 393 K) most stable system of every tested single adsorption site, labelled according to fig.3.1a (p. 53).

					ΔG ads		31 P NMR
	facet sites			(25°C) (120°C) (300°C)	1 st	2 nd
	(1 0 0) ν 0 +ν 0	-98	-74	-27	2.1	1.3
	(1 0 0) ν a 1 +ν 0	-93	-77	-47	-3.2	5.5
	(1 0 0) ν b 1 +ν 0	-105	-91	-64	-2.0	4.1
	(1 0 0) ν c 1 +ν 0	-84	-69	-41	-8.8	9.0
	(1 0 0) ν a 1 +ν a 1	-82	-75	-62	5.5	-6.2
	(1 0 0) ν a 1 +ν b 1	-102	-97	-86	11.6	-5.3
	(1 0 0) ν a 1 +ν c 1	-90	-85	-75	-11.1	4.9
	(1 0 0) ν ab 2 +ν 0	-107	-101	-88	1.6	-6.0
	(1 0 0) ν bc 2 +ν 0	-119	-114	-102	1.5	12.1
	(1 0 0) ν ab 2 +ν a 1	-116	-119	-126	12.7	-6.0
	(1 0 0) ν ab 2 +ν c 1	-91	-95	-101	8.8	0.4
	(1 0 0) ν bc 2 +ν a 1	-105	-110	-118	-2.6	-0.6
	(1 0 0) ν aa 2 +ν bc 2	-102	-116	-141	1.7	8.3
	(1 0 0) ν aa 2 +ν cb 2	-66	-80	-106	8.8	12.4
	(1 0 0) ν abc 3 +ν 0	-77	-80	-85	9.6	-4.7
	(1 0 0) ν abc 3 +ν a 1	-91	-105	-132	-0.8	-4.9
	(1 0 0) ν aaC 3		+ν bc 2	-84	-106	-148	9.4	6.4
	(1 0 0) ν aaC 3		+ν cb 2	-76	-98	-139	1.0	7.7
	n(1 1 0) ν 0 +ν 0	-101	-77	-33	4.7	0.1
	n(1 1 0) ν b 1 +ν 0	-112	-99	-72	10.6	8.1
	n(1 1 0) ν e 1 +ν 0	-89	-75	-48	-10.7	0.3
	n(1 1 0) ν b 1 +ν e 1	-111	-106	-96	-8.3	3.4
	n(1 1 0) ν e 1 +ν e 1	-83	-78	-69	12.3	9.6
	n(1 1 0) ν de 2 +ν 0	-102	-97	-88	-15.0	13.4
	n(1 1 0) ν bb 2 +ν e 1	-89	-92	-99	-13.7	9.3
	n(1 1 0) ν de 2 +ν b 1	-101	-106	-115	-14.2	-0.2
	n(1 1 0) ν bb 2 +ν de 2	-114	-127	-153	-13.2 -10.9
	n(1 1 0) ν deA 3	+ν 0	-89	-93	-101	-9.4	13.3
	n(1 1 0) ν bb E 3	+ν d 1	-76	-90	-116	-6.4	0.0
				continued on next page

Table A .

 A 3.: Free enthalpy of surface pyrophosphate formation in kJ mol -1 and 31 P NMR shift in ppm with respect to to 85 % H 3 PO 4 for the (at 393 K) most stable system of every tested single adsorption site, labelled according to fig.3.1a (p. 53). Free enthalpy reported with respect to a single orthophosphate in gas phase, in order to allow comparison with orthophosphate adsorption.

			ΔG ads		31 P NMR
	facet site	(25°C) (120°C) (300°C)	1 st	2 nd
	(1 0 0) ν 00	-49	-35	-8	-7.9	1.7
	(1 0 0) ν a 10	-96	-89	-76	-19.2 -23.8
	(1 0 0) ν a,a 11

Table A .

 A

		3: continued from previous page
			ΔG ads		31 P NMR
	facet site	(25°C) (120°C) (300°C)	1 st	2 nd
	n(1 1 0) ν de,ab 22	-51	-74	-117	-25.9 -35.9
	n(1 1 0) ν ee ,bc 22	-35	-59	-103	-32.6	4.7
	n(1 1 0) ν bb E 30	-20	-34	-59	-16.1	-5.0
	n(1 1 0) ν deA 30	-43	-56	-83	-27.0	-5.6
	n(1 1 0) ν bb E,e 31	-61	-84	-125	-21.6 -13.8
	n(1 1 0) ν deA,b 31	-57	-79	-121	-25.2	-8.6
	n(1 1 0) ν bb B,ee 32	-35	-67	-127	-18.3	8.3
	n(1 1 0) ν bb E,ee 32	4	-28	-89	-14.1	12.0
	n(1 1 0) ν deA,ab 32	-16	-47	-107	-21.7 -13.7
	R(1 1 0) ν 00	-84	-79	-69	-6.3 -13.6
	R(1 1 0) ν c 10

Table A .

 A 4.: Tested 31 P NMR references at different levels of theory. The experimental MAS NMR iso-shifts are given with respect to 85 % H 3 PO 4 , while the computational results are given as calculated. The zero-shifts are the difference between experimental and computational value, i. e. the extrapolated values for 85 % H 3 PO 4 . A: GIPAW calculation at PBE-D3 level on the experimental X-ray structure; B: Full optimization and GIPAW calculation at PBE-D3 level; C: Full optimization and GIPAW calculation at TPSS-D3(BJ) level; D: Full optimization at TPSS-D3(BJ) level, GIPAW calculation at PBE-D3 level.

		31 P iso-shift, δ		zero-shift to 85 % H 3 PO 4 , δ
		Berlinite Variscite Wavellite	Berlinite Variscite Wavellite	average
	exp.	-25.3	-19.2	-11.2
	A	-331.4 -324.7	-318.6	-306.1 -305.5	-307.4	-306.4
	B	-322.7 -313.7	-305.3	-297.4 -294.5	-294.1	-295.3
	C	-325.1 -318.6	-314.1	-299.8 -299.4	-302.9	-300.7
	D	-324.3 -318.7	-310.4	-299.0 -299.5	-299.2	-299.2

ions with only one P -O oxygen atom in hydrogen bonds.

Acknowledgements

solid-liquid interface, this vapor entropy is substituted for a H 2 O liquid entropy that contributes the overall FES. At drying conditions the free enthalpy difference between ν b 1 and ν bb 2 is 32 kJ mol -1 (cf. tab. A.1, p. 154), which is almost equal to the to the solid-liquid interface (30 kJ mol -1 ). While it is possible that ν bb 2 is also the most stable adsorption mode at the solid-liquid interface it is unexpected to observe it with the same free energy as for drying conditions. One likely reason here is that the WTmetaD simulation is not ideally converged. Nevertheless, the same trend, i. e. stability order is clearly visible for the last 10 ps of the simulation, hence at least the stability is expected to be correct.

In order to verify the combined CV in eq. 5.4, the positions of the minima are studied in more detail. The minima are not exactly at 2, 0 and -2 because CN w and CN p are defined by smooth functions and surrounding water and phosphate oxygen atoms still produce subtle contributions at distances around 4Å, which would definitely not be considered coordination. Note that these smooth functions are necessary in order to generate the forces that pull the to-be-bonded oxygen atoms to its site. It is possible to achieve numeric values close to 2, 0 and -2 by using stiff functions for the coordination, but this is at the expense of a less effective metadynamics sampling, since the forces vanish too quickly around the cutoff distance (cf. section 2. structures can be explored by adapting the reaction in eq. 5.5 for the specific site. In

For phosphorylation of µ 2 -OH it was thus also necessary to add a CV for deprotonation of this site. Furthermore, the starting structures can be specified: ν de 2 for ν deA 3 , ν bc 2 for ν bcB

3

and ν e 1 for ν eB 2 .

The collective variables and bias

Three distinct sites were chosen to carry out WTmetaD simulations of the surface phosphorylation, using two different reactions (eq. 5.5 and eq. 5.6). The principle is the same in all cases: the collective variables (CV) for the surface phosphorylations in fig. and ν 3 structures (also illustrated in eq. 5.5 and fig. 5.4a). For the formation of ν eB 2 from ν e 1 the value varies from 3 to 2. In the case of the pyrophosphate reaction (eq. 5.6) CN w is replaced by the P -O coordination CN p with the four oxygen atoms from the leaving phosphate group; in all cases CN p varies between 0 and 1, since pyrophosphate to (a) eq. 5.5, with the formation of H 2 O as leaving group, and (b) eq. 5.6, with instead the formation of a phosphate ion at the surface. Concerning proton transfers and charges, the figures are simplified by only involving µ x O sites with x = 1 for water substitution sites (cf. section 5.2) and x = 2 for alumina oxygen, each assumed to carry one negative formal charge. The released HPO 2 - 4 can be protonated to H 2 PO - 4 but it is assumed that it stays close to the interface. The desorption into bulk H 2 O is neglected. is characterized by a single P -O -P bridge, which is hydrolyzed during the reaction.

Regardless of the different integer values of the different coordinations, their range is always ΔCN = 1 and equivalent WTmetaD simulations can be carried out to explore the FES of those phosphorylations. As before, the minima of the resulting FES are not exactly at the described integer values, due to the smooth switching functions employed to define a P -O coordination. The bias is applied by depositing gaussians in the visited FES regions, thereby "flooding" minima and automatically exploring the FES area of interest. Here, a gaussian height of 3 kJ mol -1 was used, and for CN a , CN w and CN p gaussian widths of 0.04. In section 5. in case of LRR and XGBoost, since the RMSE range for both methods is very narrow.

In contrast, the MLP neural network was too small for the much larger feature input, leading to substantial fluctuations.

Hyperparameters

For LRR, the regularization parameter λ was scanned, MLP was limited to a single hidden layer, for which the number of all-to-all connected neurons was varied, and for XGBoost, both the tree depth and number of trees were scanned. The scans of all three hyperparameters give different insight into the machine learning performance.

The effect of different values for λ were tested at fixed SOAP parameters (rcut = 5.9, lmax/nmax = 8) (fig. 6.4a), whereas the hyperparameters for MLP and XGBoost were

Conclusions and Perspectives

The main goal of this work was to generate a better understanding of how phosphates A. Appendix

A.1. Electronic structure methods

Today, computational chemists have access to a wide variety of methods in order to produce atomistic models for their systems of interest. These methods can largely vary in computational expense and general accuracy as well as in robustness and practicality.

To start from the fundamentals, the most basic approximations are introduced in this section with a short description of the underlying physics, followed by an outline of a few notable methods. Since quantum mechanical particles are represented by wave functions which are continuous probability amplitudes, the underlying set of basis functions is a crucial part for the accuracy of the systems. Both aspects are explained in more detail below.

A.1.1. The total energy

A system's energy is a conserved quantity of any physical system from large interacting bodies down to molecules and crystals. The latter are described in atomistic models of atom nuclei and electrons. An accurate determination of the behavior of either of these particles requires quantum mechanics. The total energy of a quantum mechanical system is given by the Schrödinger's equation which defines the Hamilton operator Ĥ first approximation, we separate nuclear and electronic terms:

The terms each refer to the sum of operators of involved electrons e and nuclei N, i. e.

VNN refers to the electrostatic repulsion of all nuclei. Nuclear movement is much smaller than electronic movement given their mass ratio of m N /m e > 10 3 , as considered in the Born-Oppenheimer approximation. [START_REF] Born | Zur Quantentheorie der Molekeln[END_REF] Nuclei are then given as fixed point charges, and

TN is zero by that definition. VNN is included to the electron part as constant term. Even though this is not strictly necessary, it simplifies the construction of energy derivatives like forces and Hessian matrices. The full Schrödinger equation can be approximated by its (stationary) electronic part and an independent term T N for nuclear motion.

In prospect of molecular dynamics, this equation provides the total energy of a dynamic system, where the nuclear positions change over time, r N (t), while the electronic ground state is stationary, that means the electrons instantly adapt to the moving coordinates.

In general T N is the expectation value of a nuclear wave function that is only defined by nuclear movement in the potential given by the electronic wave function and VNN . In order to solve the electronic part of the Schrödinger equation, a mathematical expression for the wave function is needed.

The first approximative solution was given by Hartree in 1927 a year after publication of the Schrödinger equation itself. [150][148] He expressed the many-electron wave function as product of single-electron wave functions, the Hartree product. Since it violates the principle of antisymmetry of proton exchange in the wave function, [START_REF] Slater | The Theory of Complex Spectra[END_REF] it was extended in the following few years by constructing a determinant of Hartree products, the Slater determinant, which is a construction that satisfies this principle. This method, later

A.2. Database tables

A.3. Computational 31 P NMR reference

Experimental NMR data is by convention tabulated with respect to a well defined reference compound. These compounds are usually liquids as NMR spectroscopy is particularly suitable for molecules in solution. In case of 31 P NMR, 85 % H 3 PO 4 is used as reference which defines the origin of the 31 P NMR scale by convention. GIPAW calculations calculate the chemical shift from the electronic structure by estimating the shielding tensor of the all-electron wave function. Its reference for the chemical shift is hence the "naked" nucleus, P 31+ , which is experimentally impractical. Thus, in order to compare experiment and theory, a similar reference needs to be established.

Computational modelling is much easier for static structures, like isolated molecules (gas phase) or crystalline solids. Here aluminium phosphate minerals were considered, since they should have a rather similar bonding as the discussed surface adsorptions. Three minerals were considered: berlinite, variscite and wavellite. All of them are characterized by X-ray crystallography and solid state MAS NMR. [START_REF] Bleam | 31 P solid-state nuclear magnetic resonance spectroscopy of aluminum phosphate minerals[END_REF][START_REF] Goiffon | Etude comparée à diverses températures (173, 293 et 373°K) des structures de type quartz α des phases M III X V O 4 (M III = Al, Ga et X V = P, As)[END_REF][START_REF] Kniep | Variscite[END_REF][START_REF] Araki | The crystal structure of wavellite[END_REF] Berlinite is the only waterfree aluminium phosphate with a AlPO 4 stoichiometry and which has the same structure as β-quartz (both aluminium and phosphorus are in tetrahedral coordination). Variscite (AlPO 4 • H 2 O) is structurally similar to berlinite in that it consists of a three-dimensional mesh of alternating P and Al units; the Al units are in octahedral coordination due to hydration. Wavellite is an aluminium phosphate hydroxide, Al 3 (PO 4 ) 2 (OH) 3 • 5 H 2 O, with a reduced relative phosphate content. Hence, wavellite might be structurally most similar to the discussed surface adsorption states. The phosphate sites in all three minerals are equivalent, producing only a single signal.

It is well known that NMR results are very sensitive to the local geometry. Three different geometries were tested starting from the experimental X-ray crystal structure.

Since X-ray diffraction cannot resolve hydrogen positions, they were optimized at the DFT level. In case of Wavellite, free water molecules are in cavities, which were reported is also used throughout the work (PBE-D3) the structural effects provide information about the overall accuracy that can be expected. The deviations within the zero-shift should be negligible if any of the three minerals provided an accurate reference system.

For PBE-D3 the deviations from berlinite to variscite and wavellite is notable with 3 ppm.

PBE results based on the TPSS-D3(BJ) geometry seems to be most accurate, however, when using the electronic structure of TPSS, the zero-shift prediction of wavellite changes by 3.7 ppm compared to the PBE electronic structure on the TPSS geometry. 7. These numbers are accessible in PLUMED [START_REF] Bonomi | PLUMED: A portable plugin for free-energy calculations with molecular dynamics[END_REF][START_REF] Gareth | PLUMED 2: New feathers for an old bird[END_REF][START_REF] Bonomi | Promoting transparency and reproducibility in enhanced molecular simulations[END_REF] with the keyword COORDINATION-NUMBER, introduced in ref. [172]. Instead of summing the whole coordination within one or between two groups of atoms, as described in section 2.2.4, the coordination number cn i of each atom i in one group is calculated and accessible for further operation:

Any scalar quantities can be derived from these values, like e. g. the mean value or the quantity of cn i smaller than a particular value, or within a range. In order to reach a smooth transition from 0 to 1 bond, a different switching function (based on a gaussian curve) was used:

The smoothness is set with x 0 as gaussian width. This switching function vanishes rapidly with the difference x-d 0 , whereas a switching function such as in eq. 2.8 produces a "tail end" which slowly converges with 0 (see also fig. 2.1, p. 46). Such a tail end is undesirable as we want to use the result from a switching function in cn i and this does not permit convergence with 0 or 1 as a result, since if cn i = 1, then s(cn i ) > 0.

Here, the quantity of cn i of one or more was used in order to count the oxygen atoms that bind to phosphorus, irrespective of the number of P -O bonds. The parameters