MULTI-SCALE MODELLING OF WATER AND HYDROXIDE IN SOLIDS AND SOLUTIONS ●

This report discusses some of the most pressing challenges that need to be overcome for computational condensed-matter chemistry to become fully accepted, at par with experiments. The prospects are rather bright. By means of a few examples, all connected to the bound water molecule and hydroxide ion, and their mysteries, the unique capabilities of theoretical calculations will be demonstrated. They provide new insights and details, and can even surpass experiments in accuracy.


INTRODUCTION
Chemical industries worldwide make use of catalytic surfaces to produce enormous amounts of thousands of different chemicals.The very same factories use catalytic surfaces to mitigate the consequences of this production through pollution control processes.Many of the key processes in environmental and atmospheric chemistry, electrochemistry and materials chemistry are even governed by ion-water interactions and water/solid interfaces.As water is omnipresent, it affects a range of important chemical processes at functional surfaces and interfaces, with beneficial or damaging consequences.To improve and develop the materials themselves and the accompanying molecular processes, new knowledge and insights will be needed.These developments are unfortunately hampered by the fact that the atomic-level mechanisms that govern the key functionalities of materials and their interfaces are usually unknown.
Here Computational Materials Chemistry can be of immense help as it provides results of unmatched detail, as well as the needed atomic-level understandingif the computational models and methods are accurate and realistic enough.Also, on the experimental side, new powerful characterization techniques for surface and interface systems have emerged in the last decades and new infrastructures are under development in Europe and in particular in Sweden (Max-IV and ESS).The role of modelling will now become even more important as the dimensions of experimental and computational targets approach each other.
A main scientific challenge in focus of my own research concerns the exploration of the links between microscopic structural features and the properties and (re)activity of solid surfaces and nanoparticles.We also devote much effort to the design of modelling strategies that can make "calculations meet reality", or reach the slightly less ambitious goal of making "calculations meet experiments".Water has a prominent position in this effort of ours and this molecule insists on much attention because of its importance, intricacies and ubiquity, and because it is such a delicate and difficult molecule to model.
In these endeavors we are grateful for the many important contributions made by Professor Bojan Šoptrajanov and his research team regarding molecular vibrations and structure-property relationships in crystalline hydrates and other hydrogen-bonded and aqueous systems.These findings are of great scientific value not only to myself and my colleagues in the former Hydrogen-bond project led by professor Ivar Olovsson at Uppsala University, but also to the international scientific community at large.Below I will give a brief overview of the research I and my group are involved in, the challenges that we are up against, and some promising results.I will focus on our studies involving water and vibrational calculations.

COMPUTATIONAL CHALLENGES
Computational materials chemistry is a challenging topic because: • Chemistry deals with interactions, reactions, and bonding between species, and not just with their individual properties.Complicated electronic structures and their changes often need to be described.
• The thermodynamic conditions, in particular temperature effects, as well as dynamical processes, also often need to be taken into account.
• The structures of real materials systems are generally very complex, containing imperfections and defects, interfaces or multiple phases, and attached molecules.Large-scale computational models are needed to accommodate such diversity.
Given this complexity, the development of methods and strategies to build relevant models is clearly central.There exists no single model or work-flow capable of treating simultaneously the full length-, time-and energy-ranges, from the electronic level to macroscopic devices.Instead a set of models and work-flows along a multi-scale ladder, illustrated in Figure 1a, need to be interconnected in clever ways.At the core of this daunting task lies the issue of obtaining "sufficient accuracy" for the problem at hand.Together with my co-workers, my own piece in this puzzle is the development of multi-scale methods and models to bring modelling closer to the complexity of realistic applications.With combinations of high-and low-level electronic structure and force-field models, we explore structure-activity relations for inorganic crystals, surfaces and nano-particles, and in aqueous media.

RESULTS -METHOD AND MODEL DEVELOPMENT
Some methodological areas in need of close attention and developments that we contributed to are: • Multiscale modelling techniques: on the one hand protocols to move "sequentially" and seamless-ly along the multiscale ladder (see Figure 1a) and on the other "concurrent" approaches such as the embedding, or QM/MM, models (see Figure 1b).
• High-quality alternatives to standard quantum-mechanical methods, e.g.advanced force-fields or approximate quantum-mechanical (QM) methods.
• Improved models to describe the shortand long-range interactions with a system's surroundings, such as solvent effects of molecules in solution or at solid-liquid interfaces (see Figure 1c).
• Strategies and software that allow us to link simulated data to experimental results (spectra, images, etc.); see Figure 1c.
Our multi-scale modelling efforts focus on the "chemical scales" within the multi-scale ladder in Figure 1a.In the lower left-hand corner (QM), we use various flavors of density functional theory (DFT) and wavefunction-based methods.One step up the ladder, we currently use the self-consistent charge density functional based tight-binding approximation (SCC-DFTB), which is an approximate DFT method, about two orders of magnitude faster than standard DFT calculations.The accuracy and transferability of SCC-DFTB models crucially depend on a set of parameters, which have to be carefully optimized, an area to which we devote considerable effort.
At the next step up on the ladder, we dispose of the explicit electronic information and use/develop parametrized force-fields (FFs), often of a sophisticated form that allows for bond breaking and formation (e.g. of the ReaxFF type).Forcefields allow to reach experimentally relevant system sizes and time scales.As indicated in Figure 1, we use data obtained at the QM-level to parametrize models higher up in the ladder.A great deal of research and insight is required to determine which degrees of freedom are safe to sacrifice when approximations need to be made.Some recent publications here are our SCC-DFTB work on H 2 O/ZnO (10-10) [1] and CeO 2 [2], and our Re-axFF work on CeO 2 [3].
Data generation is one side of the coin, data analysis for property calculations and materials characterization is the other.Both aspects require adequate models and workflows, and both present theoretical and computational challenges.Our simulated IRRAS spectra of the CO/TiO 2 (110) system compared to experiment is a good example [4].Another one is our generation of simulated STM images of a defective CeO 2 (111) surface, where we mimicked experimental images for three cases of surface defectsoxygen vacancies, F substituents and H substituentsand challenged the interpretation of experimental STM images prevailing in the literature [5].

RESULTS -SCIENTIFIC APPLICATIONS
Three examples from our current research are presented below.They are connected to experimental investigations, either as collaborations, or using data from the literature.

Water in crystalline hydrates
As already commented on above, water is an important and treacherous molecule: seemingly so simple, but in practice so difficult to model satisfactorily.We have found that theoretical calculations for ionic hydrate crystals present a unique opportunity to gather novel information about the water molecule and its structure-property relations [6].Ample high-quality structural data are available from diffraction experiments for a large number of crystalline hydrates where water coexists with metal ions, sulphates, nitrates, halides etc., and we, the modelers, can then focus on modelling the properties themselves.Our ongoing study of structureproperty relations for crystalline hydrates can be seen as a "Materials informatics" project for the water molecule.
We recently proposed that highly hydrated crystals, represent an "economical" advantage, allowing us to collect information about many water molecules "in one go"; for example the nine structurally different water molecules in the Al(NO 3 ) 3 • 9H 2 O crystal. Figure 2 displays a property which is difficult to extract from experiment, namely the in situ polarisation of water inside a crystal.The contour diagrams display two of the nine unique water molecules and reveal that the nature of each of them is strongly altered by the surroundings; the reference in the figure is the gas-phase molecule, which would thus have no contours at all.The crystal-induced polarisation in Figure 2 is especially dramatic for the molecule labeled W5, which is bound to an Al 3+ ion.
This result ties back to the multi-scale picture in Figure 1c, which displayed the computational scheme that we used to generate the infrared vibrational spectrum for the first-shell water molecules around the Al 3+ ion in an aqueous solution.Figure 2c demonstrates that water molecules binding to Al 3+ in the solid state become highly polarized and we may infer that this should happen in aqueous solution as well, thereby offering a likely explanation for the large frequency downshift observed in the vibrational spectrum in Figure 1c.3a.The result is not impressive.In Figure 3b the results from the X-ray diffraction data-set have been removed, leaving only the three neutron diffraction studies.The result is somehow improved, although no correlation is discernible.We also note that the intra-molecular water OH distances in the crystals are sometimes shorter than the gas-phase value, a result which is obviously unphysical.The reason behind these systematic experimental distance "errors" is known and was elaborated on in Ref. 6.The resulting 'r(OH) vs. R(O•••O)' correlation curve from the quantum-mechanical crystal calculations is shown in Figure 3c.As should be, the scatter of the points is reasonable and all points lie above the free-water value.Clearly, the calculations perform better than experiment in locating the H atom.
We are of course concerned with the accuracy of the methods and models we use.In a current (unpublished [7]) study we investigate the performance of a range of dispersion-corrected DFT methods concerning the structure of the four crystalline hydrates mentioned above.Without going into details here, we note from the dart board in Figure 4 that, while the traditional LDA method is not competitive, many of the methods reproduce the experimentally determined hydrogen-bond distances very well.
Next we will turn to the OH vibrational frequencies.We routinely calculate anharmonic intra-molecular OH vibrational frequencies as the OH vibrational anharmonicity is large and very frequency-dependent.Moreover, we routinely calculate uncoupled OH vibrational frequencies to make use of the fact that an uncoupled OH frequency (i.e.corresponding to a one-legged OH vibration) only pertains to one OH oscillator; its gas-to-crystal frequency shift is a descriptor of only that oscillator's environment.The uncoupled OH vibration is thus more informative than the coupled ones, as far as intermolecular interactions go.For the very same reason experimentalists often try, whenever possible, to perform isotope-isolated experiments.Figure 5 shows our calculated OH frequencies for the four crystalline hydrates mentioned above against the electric field strength component along the vibrating OH bond, probed at the equilibrium H position.The figure displays quite a good correlation and a rather narrow distribution around the least-squaresfitted line.A comparison of our curve with that of Auer and Skinner [8] for instantaneous structures in liquid water reveals that our correlation is tighter, in spite of the fact that it it covers larger frequency and field ranges.

Water on surfaces
We are also interested in water molecules on surfaces: we have studied intact and dissociated water molecules and the transfer between them on metal oxide surfaces.Figure 6 is a snapshot from a Molecular dynamics simulation of a thick water film on wurtzite ZnO (10-10) [9]; it illustrates the dynamic nature of the system and shows that many of the water molecules are dissociated also on the common ZnO  surface (without extra steps; stepped surfaces were also investigated in Ref. 9).The Molecular Dynamics simulations in Figure 6 were based on a ReaxFF force-field that was derived from electronic QM calculations in a consistent manner, as indicated by the red arrows in the right-hand part of the multi-scale scheme in Figure 1a (in a collaborative project between A. van Duin (Penn State U) and the Uppsala group).We have also studied the ZnO and water-ZnO systems in the approximate (tightbinding) DFT approach in the SCC-DFTB formulation and developed a parameter set for these systems, again based on QM calculations in a consistent setup.[1] Such an endeavor can be illustrated by the light-green arrows in the right-hand part of Figure 1a.This work was developed as a collaboration with the group of Th.Frauenheim (Bremen Center for Computational Materials Science).
Using a range of dispersion-corrected DFT functionals and water coverages, we recently studied water on the archetypical NaCl(001) and MgO(001) surfaces [10], and found that, for both NaCl(001) and MgO(001), the dispersion-flavored functionals stabilize the water-surface interface by 20%-40% compared to the PBE-results.A monolayer water coverage on MgO(001) leads to a mixed overlayer with both intact and dissociated water molecules in an intricate hydrogen-bonded scheme.For NaCl(001), the water molecules remain intact for all water coverages, which allows a more elaborate analysis of "the nature of a surface water molecule" compared to the gas-phase.Thus to probe the strength of the perturbations from the surface and the rest of the water layer on an adsorbed water molecule, we calculated water dipole moments and found an increase of up to 85% for water at the MgO(001) surface and 70% at the NaCl(001) surface, compared to the gas-phase dipole moment.Likewise, for intact water molecules, it was meaningful to divide the total adsorption energy into water-surface and surface-surface contributions: they were found to be approximately of the same magnitude but the dispersion correction affected the water-surface interactions more than the water-water interactions.

Water in aqueous solution
It is particularly challenging to calculate vibrational OH spectra in a liquid because the liquid configurations have to be properly sampled or the spectrum will be scewed and misleading.We have continued to develop this approach ever since our first computational study of the (anharmonic) infrared OH spectrum for liquid water in 1991 [11].Our approach contains several stages.A Molecular Dynamics or a Monte Carlo simulation is performed, a suitable number of snapshots are collected for analysis, and for each of these a large number of electronic QM calculations are performed to generate potential energy curves for the OH stretching modes for each of the first-shell water molecules around the ion in solution.In our approach, this step is in fact in itself a QM/MM calculation, where our preferred QM method for water is B3LYP.This step is followed by a quantum-mechanical calculation for each potential energy curve, but here it is a nuclear QM calculation, which generates the vibrational energy levels.The energy difference between the ground and first excited vibrational energy levels is collected into a frequency spectrum.The scheme is illustrated in Figure 1c.In practice we sometimes do not include the first step (a QMbased force-field generation) but instead go straight to the MD stage using an available force-field or an ab initio-MD approach.
Figure 7 shows results from two collaborative projects with Lj.Pejov (Sts.Cyril and Metho-dius University).Figure 7a displays a snapshot from the force-field-based MD simulation of Li + (aq) [12]; all water molecules except those closest to the cation have been removed for visual clarity.There resulting gas-to-solution frequency shift is given in Table 1, together with similar results for a divalent and a trivalent ion.The agreement with experiment is overall good, and this is actually also true for the absolute frequencies (not shown here).The water OH-frequency of bound water molecules are always downshifted with respect to the gasphase and this effect is very large in the cases of Mg 2+ (aq) and Al 3+ (aq); in both cases the OH band origin even lies much below that of liquid water.
The last line in Table 1 does not refer to water but to the OH vibrations of the hydroxide ion itself, when immersed in water.A snapshot from the ab initio MD simulation is given in Figure 7b and illustrates the strong hydrogen-bonds donated to the hydroxide's O atom, and the much loser and cage-like water structure around the hydroxide's H atom, as discussed in Ref. [13].The eperimentally measured gas-to-solution frequency shift is seen to be well reproduced by our calculations.The shift is quite modest, much smaller even than that of the water molecules around the monovalent Li + ion, in spite the strong hydrogen bonds donated to the hydroxide ion.How can this be?The explanation was given in Ref. [13] and confirms that the hydroxide ion in water conforms with the behavior that we had earlier found in the gas-phase [14]: when the ion is exposed to a small electric field (from an ex-ternal source, e.g. its surrounding molecules) the frequency is up-shifted (opposite to what happens for the water molecule), reaches a maximum as the field is increased and then decreases and, for very large electric fields, gives a large downshift compared to the gas-phase OH -frequency.The effect is illustrated in Figure 8 for the ideal cases of a water molecule and a hydroxide ion (separately) exposed to a uniform electric field [14].We conclude that the study in [13] both demonstrated that the calculations were able to reproduce the experimentelly determined gas-to-liquid frequency shift (and actually here also the absolute values), and provided an explanation for the observedpositive and modest frequency shift observed for the hydroxide ion, which is contrary to our experience from the behavior of water molecules.
In summary, water molecules in aqueous solutions, in crystals, and on surfaces are always found to be downshifted (red-shifted) in OH frequency by their surroundings.Examples given in this mini-review, demonstrate, however, that while the hydroxide ion can be red-shifted when bound in a strongly polarizing environment, it is in fact often found to be up-shifted (blue-shifted).We found this parabola-like behavior for the hydroxide ion in crystalline hydroxides as well [15].The explanation underlying this qualitatively different behavior of the water molecule and the hydroxide ion is related to the relative signs and magnitudes of the permanent and induced dipole moment derivatives along the OH-stretching coordinate: the sign relations are different for water and OH - [14].The simulation is described in Ref. [12], where references to the force-field generation is given.(b) Snapshot from an ab initio-MD simulation for OH -(aq).The simulation is described in Ref. [13]. a References to the experimental work: see references within [12] and [13], respectively.b One of the two experimental values that have been presented in the literature; see [12].This is the energetically favored direction.

CONCLUSION, IMPACT AND OUTLOOK
The examples above demonstrate that Computational Materials Chemistry delivers unique quantitative and qualitative information, and knowledge that is unattainable from experiments.In my opinion, the most valuable and long-lasting results from modelling are likely to be the scientific insights that can be formulated in terms of rules-ofthumb and structure-property relations: they will save time and money for many categories of endusers of modelling results, from high-school science students to industrial stakeholders.
Considerable methodological challenges remain, however.In fact, one of the narrowest bottlenecks, and one which currently limits the wider use of modelling in e.g.industry, is the lack of appropriate models to treat large-scale systems of realistic complexity, and their time evolution.This was recently demonstrated in a survey [16] conducted by the European Materials Modelling Council (EMMC, https://emmc.info)where interaction models of high accuracy and computational efficiency were requested as well as techniques to per-form multi-scale modeling of materials.In fact, the Royal Swedish Academy of Sciences awarded the Nobel Prize in Chemistry for 2013 to Martin Karplus, Michael Levitt and Arieh Warshel "for the development of multiscale models for complex chemical systems" but it was the development of QM/MM models for biological molecules that was in focus in this prize.In the field of materials modelling, however, the use of QM/MM models and other multiscale approaches (cf. Figure 1) is much less mature.Significant development efforts are needed.
In this brief review we have discussed diffraction experiments, vibrational spectroscopy data and, of course, theoretical calculations.The importance of combining different techniques was early on realized by Professor Bojan Šotrajanov and his research team at the Sts.Cyril and Methodius University.Many of my colleagues and I at the Structural chemistry program at Uppsala University are grateful for many years of inspiring collaboration and friendship.Some of the first contacts were shaped at the fourth Horizons in Hydrogen Bond Research meeting at Sånga-Säby in Sweden in 1980.This has been followed by research visits in both directions, including those of Professors Gligor (Glišo) Jovanovski and Ljupčo Pejov, who both spent a research year in Uppsala.The scientific discussions between Uppsala and Skopje are still very much alive.

Figure 1 .
Figure 1.Overview of the multi-scale modelling efforts in our group.(a) The multi-scale ladder to the left indicates systems with electronic and atomistic resolution, while more coarse-grained (mesoscopic) models and continuum models, which lie in the gap between atomistic models and the Real World, are not included in the figure.The figure to the right lists the methods most frequently used in our group.The arrows pointing upwards along the ladder indicate that fine-grained model are used to parametrize more coarse-grained models.(b) Our QM/MM approach for solid ionic surfaces.(c) An illustration of our multistage strategy to calculate vibrational spectra in liquid solutions: QM elec => FF-based MD => QM elec /MM potential energy curve => QM nuclear calculations of the vibrational energies.The sample system in this figure is Al 3+ (aq) with 10000 water molecules interacting through an ab initio-generated force-field with effective three-body terms [12].

Figure 2 .
Figure 2. Polarization of water molecules inside a crystal of the water-rich Al(NO 3 ) 3 • 9H 2 O crystal, from quantum-mechanical calculations [6].The figures display both the difference electron density (total electron density of the crystal minus the sum of the electron densities of the isolated building blocks) and the total dipole moment of the water molecules.Red areas in the difference electron density maps indicate "more electrons" than in an isolated molecule, blue means "fewer electrons".

Figure 3 .
Figure 3. r(OH) and R(O•••O) correlations from experiments and from our calculations for four highly hydrated crystals, Na 2 CO 3 •10H 2 O, MgSO 4 •7H 2 O, MgSO 4 •11H 2 O, and Al(NO 3 ) 3 •9H 2 O. (a) Correlations based on X-ray and neutron diffraction experiments (references are given in Ref. [6]).(b) The same figure as in (a) except that the X-ray diffraction data have been excluded.(c) Correlations based on the optimized structures from our DFT calculations described in Ref. [6].

Figure 4 .
Figure 4. Assessment of the performance of various DFT functionals with respect to reproducing the experimental R(O•••O) hydrogen-bond distances for four crystalline hydrates.A point at the center of the dart-board would mean that the functional produces all distances in total agreement with experiments.

Figure 5 .
Figure 5. ν(OH) vs. "external" electric field.Each electric field value in the plot is the electric field calculated at one of the 74 water water H atoms in our sample of 37 water molecules in crystalline hydrates listed in the caption of Figure 3.For each water OH bond, the electric field generated by the whole crystalline surroundings outside the probed molecule itself was calculated at the equilibrium H position, and the component of this electric field along the OH bond was used in the plot.The frequency plotted is the corresponding anharmonic, uncoupled OH frequency.

Figure 6 .
Figure 6.Snapshot from the ReaxFF-based MD simulation of a water-covered wurtzite ZnO(1 0 -1 0) surface, reported in Ref. [9].The picture was drawn from the trajectory of the simulation.

Figure 7 .
Figure 7. (a)Snapshot from an MD simulation for Li + (aq) based on an ab initio-generated many-body force-field.The simulation is described in Ref.[12], where references to the force-field generation is given.(b) Snapshot from an ab initio-MD simulation for OH -(aq).The simulation is described in Ref.[13].

Figure 8 .
Figure8.Ab initio-calculated anharmonic OH frequency vs. electric field strength for the uncoupled OH stretching vibration of an isolated HDO molecule exposed to a uniform electric field directed along the vibrating OH bond, and the corresponding result for an isolated OH -ion in a uniform electric field.The calculations are described in Ref.[14].The positive field direction is from O, towards H, i.e. the field is directed as if there is a positive charge far away on the O side and a negative charge far away on the H side.This is the energetically favored direction.

Table 1 .
Calculated and experimental gas-to-solution OH-frequency shifts for water molecules in the first hydration shell of a series of metal cations (lines 1-3) and for the hydroxide ion surrounded by water molecules in an aqueous solution (line 4)