Diversity at the Water-Metal Interface: Metal, Water Thickness, and Confinement Effects.

The structure and properties of water films in contact with metal surfaces are crucial to understand the chemical and electrochemical processes involved in energy-related technologies. The nature of thin water films on Pd, Pt, and Ru has been investigated by first-principles molecular dynamics to assess how the chemistry at the water-metal surface is responsible for the diversity in the behavior of the water layers closer to the metal. The characteristics of liquid water: the radial distribution functions, coordination, and fragment speciation appear only for unconfined water layers of a minimum of 1.4 nm thick. In addition, the water layer is denser in the region closest to the metal for Pd and Pt, where seven- and five-membered ring motifs appear. These patterns are identical to those identified by scanning tunneling microscopy for isolated water bilayers. On Ru densification at the interface is not observed, water dissociates, and protons and hydroxyl groups are locked at the surface. Therefore, the acid-base properties in the area close to the metal are not perturbed, in agreement with experiments, and the bulk water resembles an electric double layer. Confinement affects water making it closer to ice for both structural and dynamic properties, thus being responsible for the higher viscosity experimentally found at the nanoscale. All these contributions modify the solvation of reactants and products at the water-metal interface and will affect the catalytic and electrocatalytic properties of the surface.

Introduction

Water surface interactions are ubiquitous in chemistry, physics, materials, and planetary sciences and for the development of life itself. A better knowledge of the water–metal structure will enable us to tackle the energy challenges that lie ahead of us,[1] as this interface is key to control the electrochemical[2] and chemical[3] processes that will generate the future energy vectors by transforming inactive molecules into fuels and back to energy with the lowest possible losses. The description of the water–metal interface relies mainly on the classical models describing electric double layer that goes back to the studies by Helmholtz. Only recently have molecular studies reached enough accuracy to present an increasingly detailed analysis of low-coverage structures formed at the interface between water and the metal by scanning tunneling microscopy (STM) and density functional theory (DFT).[1,4] Still, the nature, structure, and dynamics of water on metal surfaces and the emergence of liquid-like behavior are far from being understood.

Water is among the most studied and yet most complex compounds as both solid and liquid states present multiple phases. Three main features differentiate water phases: the density, the atom–pair radial distribution function, and the oxygen coordination shell.[5] Rearrangements in the water coordination shell and the topological ring structure have been taken as suitable descriptors for the liquid to ice transition.[6] The ice density at 273 K is 0.934 g cm–3, the X-ray diffraction gives peaks in oxygen–oxygen radial distribution function at 2.8, 3.5, and 4.5 Å,[7] and each oxygen is characterized by an ideal tetrahedrality. Both Ih and XI forms of ice, the common ones at atmospheric pressure, share these features. In contrast, the structure of liquid water still faces important challenges.[6,8−10] The most popular interpretation of the macroscopic state of liquid water is that it is constituted by a complex mixture[11] of two microscopic domains in variable proportions.[12,13] At room temperature the minor domain, low density liquid (LDL), corresponds to an enthalpy-favored tetrahedral arrangement where water molecules are four-fold coordinated. The LDL arrangement is fairly similar to hexagonal Ih both regarding the atomic coordination and the presence of two different distances: the direct chemical O–H bond at 1.0 Å and the hydrogen bond at 1.8 Å. This structure is very rigid resulting in a small entropic contribution. In turn, the major domain, high density liquid (HDL) is characterized by a lower number of hydrogen bonds with weakened strength.[5] The thermodynamic penalties of lesser H-bonds are counterbalanced by the higher entropy associated with a larger fraction of free librational modes. Therefore, in the HDL microstructure water molecules present a fraction of trihedral, instead of tetrahedral, arrangements. The experimental detection of density fluctuations in the bulk liquid is limited to around 1 nm.[13] The macroscopic features of liquid water at around room temperature are (i) the presence of a peak in O–O radial distribution function (RDF) around 2.8–2.9 and a plateau at 4.5 Å, corresponding to first and second coordination spheres in LDL; (ii) and only a single peak between 2.8 and 2.9 for HDL.[7,14,15] LDL has been claimed to be responsible for ice nucleation due to slow dynamics[6] and strong bond directionality.[16] The tetrahedrality (q) of oxygen in water is between the limiting values of the ideal LDL 1.00 and that of the isolated water value 0.75.

The water–metal interface, at room temperature, is expected to show intermediate features between those of the ordered bilayer or ice structures[1] and the labile bulk water.[7,15,17] STM experiments and complementary DFT-based simulations have provided accurate descriptions for the first wetting layer,[1] but extensions to multiple layers exist.[18,19] Notice that several of these studies are devoted to the study of ice and its nucleation,[15] and many of the experiments were performed at low temperatures. These studies have shown different degrees of wettability on metals[1,18] and the ability of some of them, like Ru, to split water into hydroxides and protons.[20] This feature adds more structural complexity to the interface. Theoretical simulations hold the key to understand the nature of the boundary between water and metals. For extended water–metal systems early molecular dynamics (MD) simulations have highlighted the need for dispersion interactions in the simulations,[21,22] and when those are included then it is possible to explain subtle details such as the structural changes for water adsorbed on gold associated with the decrease in the fraction of flat and planar molecules at the interface.[23] In addition, classical simulations, have identified the hydrophobic character of water on Pt,[24] and the role of extended time and size dynamics in water exchange close to Pt.[25−27] Unfortunately, such long time and length scales are not accessible by ab initio simulations, and classical MD cannot address the properties on reactive surfaces like Ru. Thus, only with ab initio MD it is possible to understand how the chemistry at the surface induces different interfacial water–metal structures.

In the present work, we have employed first-principles Born–Oppenheimer molecular dynamics (BOMD) to unravel the intimate structure of the interface between water and different metals typically employed as electrodes or catalyst (Pd, Pt, and Ru) driven by the water reactivity at the interface in their most common surface and identify the differences between confined and free-standing water layers of different thickness, ranging from 4 to 16 water layers (6–24 Å width). The ultimate goal is to (i) define the minimum amount of water on a metal surface for which liquid-like features appear, (ii) illustrate the differences between metals, (iii) highlight how the interfaces will affect the catalytic properties of the metals, and (iv) identify the differences induced by confinement of the water structures.

Results and Discussion

BOMD simulations for adsorption of water on metals were carried out for at least 11 ps at 300 K with 1 fs as the time step. The model systems investigated present a metal surface where different water slabs following the ice structure with an H-down configuration are added. The structures were labeled as MN(c) where M states for the metal surface, N is the number of water bilayers and (c) indicates that the water slab is confined. For Ru(0001) dissociation of 50% of the water molecules adsorbed on the surface was taken as the initial configuration.[20] A schematic representation of the models can be found in Figure . Movies for the trajectories (Pd4, Ru4, and Pt4) are available as Supporting Information.

Figure 1

Models employed for the adsorption of N-water bilayers on closed packed M = Pd, Pt, and Ru metal surfaces. Longitudinal views of systems for MN systems: Pd and Pt: two (four) free bilayers Pd2, Pt2 in a (Pd4, Pt4 in b), four (eight) confined bilayers Pd4c, Pt4c in c (and Pd8c, Pt8c in d). Equivalent systems are shown for Ru(0001): Ru2 in e, Ru4 in f, and Ru4c in g. Axial views for ice-like adsorbed Pd/Pt in h, and Ru in i. Metal, oxygen, and hydrogen atoms are represented by white, red, and black spheres; bulk water is shown as sticks, respectively.

We start by analyzing the most general features, like density and radial distribution function, and then move to more atomistically detailed terms as the local coordination of the oxygen atoms described by the tetrahedrality (q) and speciation of fragments for the reactive Ru(0001) surface. Notice that these parameters taken one by one are not conclusive, but instead the full toolbox presented here is complementary and describes the landscape in a detailed manner, as already discussed in refs (28 and 29). The final aim is to assess when these properties start to correspond to those of the liquid water, described in terms of both HDL and LDL contributions. In all simulations, no evaporation is observed; i.e., no molecule stays freely in the vacuum region. The water density (see Supporting Information for the definition) in our models is 0.95–1.06 g cm–3 for unconfined systems and 0.88–1.07 g cm–3 for the confined ones. This demonstrates that there is not a significant effect from the use of PBE-D or the lack of continuity of water in the MN models. However, we have observed that some water molecules might displace from the average z-position in the vacuum configurations. This agrees with the displacements observed by Michaelides and van de Vondele for the outermost water molecules in an ice slab.[30]

As for the analysis, we start from the most general patterns, like the radial distribution functions and then evaluate more microscopic features. The radial distribution function for oxygen, gO–O(r), averaged for the last 2 ps of simulation is presented together with available experiments of bulk ice, 122 K, and water, 298 K, at 0.1 MPa,[7,31] in Figure . All the water molecules were considered, and although the frontier (with the surface or the vacuum) might present a different behavior, we consider that this is meaningful when trying to establish the perturbations induced by the interfaces. Complementary gO–H and gH–H are presented and compared to a slab of ice following that in ref (30) in Figure S.I. 9–10. The results for Pd and Pt parallel each other provided that the same N-bilayer is compared. For both the 2-bilayer does not present the characteristic long-distance peak due to the small size of the model. In turn, M4 shows a first peak around 2.80 Å and a wide region at about 4.5 Å with a higher radial distribution value but very broad. The lack of clear signals at high distances hints to only one defined coordination sphere, typical of HDL structures. On the other side, the confined systems M4c and M8c show a second broader but clear peak at 4.8 Å, as well as weaker coordination spheres at 5.8 and 6.6 Å. The case of Ru(0001) is peculiar, because it displays a wide peak at 2.8 Å, starting at rather short distances. This shape is an indicator of a plethora of configurations where oxygen atoms accommodate around each other in different ways as a consequence of water dissociation on Ru. Again the Ru4 model is closer to the gO–O liquid water pattern than the confined ones.

Figure 2

Radial oxygen–oxygen distribution function, gO–O(r), for the water molecules in the multilayers adsorbed on Pd (orange), Pt (dark yellow), and Ru (wine) metal (M) surfaces, for the confined M4c, M8c and unconfined M2, M4 unconfined water slabs. Experimental data for liquid (water 298 K and 0.1 MPa solid black)[31] and ice (122 K dashed black)[7] are shown for comparison (gO–H and gH–H can be found in Figure S.I. 9–10). The comparison between the computed and the experimental liquid results can be found in the Supporting Information Figure S.I. 11.

From the differences between the experimentally retrieved values[31] and the results we obtained with the different models it is possible to identify that a minimum water thickness of four bilayers, 1.4 nm, is needed to retrieve the complex interplay between high and low density liquid features characteristic of liquid water. This agrees with the correlation regime of 1 nm reported in the experiments.[12] Instead confined water slabs even if much thicker show a low density structure LDL domain after 11 ps. As a consequence in the following we will mainly center on the analysis of the M4 models, retrieving the results from the other models when required only and reporting in the Supporting Information the rest of the results.

The second step is to characterize the coordination sphere of each water molecule in the slab along the BOMD trajectories. To inspect the HDL–LDL water domains we have mapped the H-related tetrahedrality (q) for each oxygen with respect to its four nearest neighboring hydrogen atoms, following the equation by Chau and Hardwick,[32]where r is the positions of the oxygen atoms in the j-molecule; see Supporting Information, Section 1.2 for further details and tests. With this definition q = 1 corresponds to a fully tetrahedral oxygen environment (LDL-like), while lower values, close to 0.84 and lower indicate less ordered structures (see right column in Figure ). The time evolution of the tetrahedrality index for each molecule in M4 models along the simulation time is shown in Figure a–c. The molecules were assigned to a layer according to their initial (t = 0 ps) z-positions. Darker colors indicate higher coordination. The accumulated values for all molecules during the entire simulation are presented in Figure d–f. The layer decomposition can be found in the Supporting Information, Figure S.I. 15. In these panels the results from the equivalent four bilayer confined slab are presented for comparison (black line). Additional values for all the other systems are reported in Figure S.I. 12−13.

Figure 3

Left panel: tetrahedrality evolution for each of the water molecules in the slab as a function of time: (a) Pd4, (b) Pt4, and (c) Ru4 models. Darker colors stand for tetrahedral coordination. Central panel: cumulative tetrahedrality index (all simulation time and all molecules) shown as normalized histogram for the unconfined (color) referenced to the corresponding confined (black). Right panel, local water configuration environments: tetrahedral (top), trihedral (middle), and dicoordinated (bottom). Complementary figures can be found in the Supporting Information Figure S.I. 12–14.

We start the analysis by the water–vacuum termination (fourth layer) and then go deep into the water layers down to the interface (top block in each panel). Because of the presence of the vacuum roughly 1:3 on the top bilayer has a trihedral configuration. The central second and third layers present subtle changes indicative of a fast and more labile environment. After the dynamics, a large fraction of the water molecules still preserve the full tetrahedrality, thus pointing toward an LDL-like behavior. However, the accumulated index in Figure d–f indicates that a significant amount have reduced their tetrahedrality to ca. 0.84. Indeed, the fraction of trihedral oxygens remained 14% in the bulk layers. Regarding the Pd- or Pt-water interfacial bilayer the light-colored areas forming long stripes of low tetrahedrality along the simulation time correspond to flat water molecules where the oxygen atoms were originally coordinated to the surface that reorient to be more in contact with the bulk water molecules. This agrees with the recent work that indicates that 49% of the water molecules on the interface with Au(111) lie flat.[23]

In comparison, the accumulated q index for the confined system, marked by the black line in Figure d–f, shows that confinement induces a much higher amount of tetrahedral molecules. The scarcity of trihedral oxygens in confined configurations is confirmed with the 8c models, Figure S.I. 12. The few trihedral molecules that appear in the confined mode revert to their tetrahedral environments in less than 0.5 ps. Therefore, all these four structures, Pd4c, Pd8c, Pt4c, Pt8c, keep their high-ordered ice-like structures during the whole trajectory. This discards confined systems for the theoretical modeling of a bulk water–metal (electrode-like) interface since they show induced rigidity and absence of HDL-like features. In addition, previous studies have assigned larger water viscosity[33,34] to an increase of the average tetrahedrality index. Our results would indicate that confined systems are more rigid and viscous in agreement with experimental results on confined water layers that have reported a more viscous liquid at the nanoscale[35] and proposed to exhibit improved electrocatalytic properties from static theoretical methods[36] (only ligand effects without dynamics were taken into account).

For Ru, the bulk layers are slightly more compact than their Pd and Pt counterparts in agreement with the observation that the first bilayer on Ru is hydrophobic.[18] From our simulations, the main reason is that water dissociation at the interface perturbs the layers on top of it. This is very clearly seen in Figure a–c; the interface layer is mainly formed by alternate dark/light colors for Pd and Pt, while the structure at the Ru interface is more homogeneous and the coordination pattern darker, which are indicative of values closer to 1. For instance, the second bilayer shows an average coordination shell smaller than the ideal liquid value. Indeed, around 35% of the water molecules present a non-tetrahedral arrangement, q values closer to HDL. In turn, the third bilayer is slightly more rigid showing higher q and LDL-like features. This structure of alternated regions with different characteristics is reminiscent of the electrical double layer.

To further understand the nature of the interface we have inspected the interfacial bilayer that is in contact with the surface by analyzing the topology of the water rings. Once the initial hexagonal symmetry is thermally relaxed the structures at the interface might display unusual water motifs. These rings have been identified through the visual inspection of the last 2 ps of the dynamics on Pd4, Pt4, and Ru4, and are exemplified by the last configuration analyzed in Figure . Pd and Pt effectively increase the water interfacial coverage to 0.74 and 0.70 ML; this is 20 (two more molecules) and 19 (one more), respectively. This densification occurs in the very early stages of the simulation, during equilibration, and once it occurs the z-displacement of these molecules is within 2 Å; see Supporting Information Figure S.I. 14. Note that although the lattice employed in the present work, and that reported in experiments, , are not commensurate, our simulations are able to retrieve both (i) the higher water coverage at the interface and (ii) the multiple motifs present. For Pd and Pt the interface water layer breaks the hexagonal lattice yielding a tessellation that includes five-, six-, and seven-membered rings. Some kind of smaller, severely distorted, arrangements appear as well. Similar structures, in particular the 6–5–7 linked patterns, have been experimentally identified when increasing coverage from ice,[37] or in water bilayers[38] on Pt, and they are theoretically predicted for Pd.[39] This is in contrast with previous first-principles dynamics simulations due to the geometry constraints (very small lateral size) employed.[22]

Figure 4

Longitudinal views of water layer closest to the metal surface in the BOMD simulations for the Pd4 (a), Pt4 (b), and Ru4 (c). The axial projection shown in the top indicates the local motifs in which the first layer water molecules coordinate to each other as tessellations (Voronoi-like) following the color code: 4,5,6,7-member cycles in gray; dark yellow, orange, wine, and larger in purple darker colors. Metal, oxygen, and hydrogen atoms are represented by white spheres and gray sticks in the longitudinal views. Hydrogen atoms are hidden in the axial projection.

The situation is completely different for Ru(0001) where the hexagonal lattice is less disrupted. In Figure c the surface hydroxyls attract one of the water molecules belonging to the adlayer bending the hexagonal pattern. The ultimate consequence is that large pseudo-rings with missing units appear. For instance in Figure c an open water hexagon (one monomer missing) and an open nine-membered cycle are identified. This patch is surrounded by mostly regular hexagonal cycles thanks to the higher binding strength of the half-dissociated water layer.

In summary, the local structure at the water−metal interface, i.e., within 5 Å, is denser and contains unusual coordination patterns for metals like Pt and Pd. In contrast, for Ru the structure at the interface keeps the hexagonal arrangement without showing any density fluctuations in the vicinity of the metal. These terms will affect the solvation sphere and how the solvent surrounding reactants and products will need to rearrange to reach the surface active site. According to Marcus theories these solvent redistributions are crucial to explain the rates of electron transfer processes.

Compared to recent results in the literature we have found some differences. For instance, the coordination patterns for water molecules at the interfaces do not retrieve five- or seven-membered rings, and mainly on-top adsorption is observed for all water molecules on Pt(111) with classical MD.[24] A similar effect together with a large fraction of flat water molecules at the interface was reported by Cao et al. when illustrating proton diffusion by reactive force fields.[25] Both effects can be caused by a too strong directionality and the strength of the water–metal bond interaction in the force fields. Extrapolation of the classical MD behavior of Pt(111) to other surfaces was suggested,[24] but our results for Ru show that for reactive surfaces the chemistry at the interface dominates. Our results are also in line with recent first-principles MD simulations,[22] performed for Pd(111) that illustrated the need for dispersion and large unit cells with a large number of water molecules. Moreover, we reinforce the need for large lateral cells to retrieve the interfacial water patterns.

A final aspect that requires attention is the study of the fragments that appear in the simulations for the water interface on Ru(0001). In static simulations and for STM images at low coverages it is normally found that at the first bilayer the number of dissociated water molecules is close to 50%.[20] Recent simulations have described the water/Ru(0001) phase diagram showing a wide range of configurations that include hydroxyl/water and hydride domains.[40] However, in larger water thicknesses the question regarding the number of ions and their mobility once formed at the interface has not been addressed. Results are reported in Figure for the systems Ru4 and Ru4c. Ru2 shows similar patterns to Ru4, and thus it is only presented in Figure S.I. 8 together with more details on the computational settings.

Figure 5

Most common ionic species resulting from the interaction of water with Ru(0001) surfaces found during the 10 ps of the MD. Free standing Ru4 (bright colors) and confined Ru4c (pale colors) are shown for comparison. Isolated protons H+ are reported in dark yellow, whereas anions are in shown in wine (OH–) and orange (H3O2–). H3O+ species are seldom observed during the simulation. The local configurations for each of these ions are shown on the right.

The species identified during the simulations are isolated adsorbed protons, hydronia in the liquid, adsorbed hydroxides, and H3O2–. Others, like the water clusters H4O2 and H6O3, have been observed, but their lifetime is very short (below 4 fs behave as transient species) and represent less than 0.01% of the total simulation time. Therefore, they have not been further considered. All the isolated H atoms (dark yellow lines in Figure ) are adsorbed on the metal surface; they do not sink into the metal or are released to the surrounding bulk of water, and the total number was constant during the run. This agrees with the tritration experiments showing that the solvent pH is preserved in the close contact of Ru electrodes.[41] Thus, they are uncorrelated to H3O+ species which only have been observed in the bulk and for very short periods of 3 fs (that actually do not enable us to fully identify them as nontransient species). In turn, the hydroxides left on the surface interact with the closest water molecule, forming H3O2– species via a shared H that shuttles to and from the two oxygen atoms. Thus, these structures oscillate from a pure H3O2– with one H exactly in the middle between the two oxygens, to polarized H2O and OH– fragments. This interconversion is reflected in Figure by the trend of the H3O2– species (orange) that almost perfectly mirror the OH– lines (wine). The H3O2– species is quite common, counting around 350 times over 1 ps.

The same fragments have also been identified in the confined Ru4c model. The presence of two surfaces in contact with the water slab increases the number of ionic species as three more water molecules were split on the upper surface. Whereas unconfined systems confirm the stability of the half dissociated layer configuration (0.50 dissociated water/Ru), in confined ones this ratio is sensibly decreased (0.33 dissociated water counting both surfaces/Ru). However, Ru4c is less dynamic than Ru4: the frequency of water splitting and reforming is much smaller than for Ru4 (13.3 ps–1 for the protons in Ru4 versus 7.9 ps–1 for Ru4c). This is due to the more ice-like structure of confined systems that let a minor number of rearrangements occur. Therefore, confined systems present differences regarding the number and dynamics of the ions at the interface and thus can be explored in an alternative way to benefit from these results of impeded mobility and transport again adding an extra contribution to the impeded mobility responsible for high viscosity.[35]

In summary, the species formed at the surface stay there and are not transferred to the bulk liquid. Thus, reactants need to reach the surface to interact with acid or basic species anchored there.

Conclusions

Simulations hold the key for the understanding of the water structure at the interface with metals, and the complex behavior including different time and length scale phenomena requires different computational approaches. The reactivity of water molecules on the surface has a large impact on the structure of the region in close contact with the metal. This phenomenon can only be observed through ab initio MD simulations. With this computational tool and analyzing the radial distribution function, tetrahedrality, interfacial configuration, and relevant fragments we have shown that the water–metal interface is different for reactive and more noble metals. At the interface between water and Pd and Pt structural patterns of five- and seven-membered rings are identified resulting in the densification of the water interface layer close to the surface. The patterns that were identified in surface science studies[38] for the first wetting layer are thus also present under water-thick conditions. Dissociation only occurs on Ru as a consequence of the electronic structure of this reactive surface; we have found the fraction of dissociated molecules and again this is in good agreement with the values reported from dissociation of the bilayer.[20] Water dissociation on Ru does not imply the acidification of the surrounding interfacial water layer as protons and hydroxyls are bound to the surface, in agreement with experimental observations.[41] These differences will severely affect the transport of active species from the bulk solvent to the surface where chemical and electrochemical reactions take place. Finally, we have observed that in order to retrieve the main characteristics of liquid water a minimum of four water bilayers (1.4 nm) interleaved by vacuum are required. The confinement induced by sandwiched configurations reduces the mobility of the layer, and this has important implications in the viscosity.[35] The present results identify new key contributions that are fundamental to reach an atomistic understanding and control of chemical, electrochemical, and photoelectrochemical processes.

Methods

The metal surfaces were modeled with slabs of four (Pd and Pt(111)) and five metal layers for Ru(0001) and a supercell. A multilayered proton-ordered form of hexagonal ice, ice XI,[42] was adsorbed on top. The different models systems were labeled in terms of the metal, M, number of water bilayers, N, and confinement, c, MN(c) shown in Figure . The single water bilayer corresponding to 18 water molecules (0.67 ML) is found to be more stable in H-down configurations. However, when growing multiple layers the H-down configurations for Pd and Pt turn in 500 ps of simulation to H-up ones. This agrees with the experimental observation that water layers reorient to form the water-tetrahedral structure. Thus, the simplest model considered in the present work, is the N = 2 double bilayer (coverage 1.33 ML), 36 water molecules in H-up configurations for Pd and Pt and the H-down for Ru, with a vacuum larger than 10 Å, Figure a,e. By replicating the double bilayer, larger coverages of 2.67 ML (N = 4) can be retrieved (72 water molecules) and were prepared both free (10 Å vacuum) and confined, Figure b,f and 1c,g. For Pt and Pd a thicker water slab with 144 water molecules (5.33 ML, (8)) in the confined mode was studied, Figure d.

DFT calculations were carried out with the Vienna Ab initio Simulation Package VASP.[43,44] Core electrons were described using the projector-augmented-wave PAW formalism.[45] The plane-wave set contained components with energies up to 450 eV. The Perdew, Burke, and Ernzerhof (PBE) functional was used.[46] Dispersion energies were accounted for by the semiempirical DFT-D2 potential with modified coefficients for the surface.[47,48] With this setup the density of water is slightly overestimated. For the unconfined systems we have investigated whether the presence of the vacuum boundary induces a tetragonal distortion to relax the density, and we have found that this is not the case; see Supporting Information. The threshold for electronic convergence was set to 10–6 eV. All initial configurations were relaxed until the forces acting on water atoms were lower than 25 meV/Å. These relaxed structures were taken as input for 11 ps BOMD in the NVT ensemble at 300 K controlled by a Nosé–Hoover thermostat.[49,50] The computational setting of BOMD was simplified to a gamma-only k-point sampling and a cutoff energy of 400 eV. The initial 1 ps of the run was taken as equilibration, leaving 10 ps of productive run with 1 ps as time step. Test with shorter time steps demonstrate no difference in the parameters investigated Supporting Information Section 2.2. For systems comprising four unconfined water bilayers a larger run was performed. In that case, the equilibration time was 4 ps, and productive dynamics were extended up to 20 ps. Extensive tests done for the structural properties: radial distribution function plotted similarly to ref (51) and tetrahedrality show that the features obtained with the present simulations are maintained if improving the k-point sampling; see Supporting Information Section 2.1.