Static and dynamic water structures at interfaces: A case study with focus on Pt(111)

Solid–liquid interfaces (SLIs) are ubiquitous in nature, and their study is highly relevant for understanding the (electro)chemical transformation processes that lead to natural corrosion or (electro)catalytic applications, such as in batteries, electrolyzers, or fuel cells.^1–6 While liquid properties at distances larger than ∼1 nm from the interface are already bulk-like^7,8 and approximately described with standard continuum models,^9,10 the structure and composition of the first (few) solvent layers in contact with a solid surface typically show a significant dependence on, e.g., the solid substrate and the thermodynamic conditions, e.g., applied electrode potentials in electrochemistry contexts.^11–19 In aqueous solutions, an accurate atomistic treatment of the substrate–water interface necessitates the explicit inclusion of at least the first water layer for a wide range of properties such as adsorption energies or the potential of zero charge (PZC).^15,20–23 On the other hand, the sensitive dependence of these properties on the interfacial water structure^24–27 necessitates, in principle, an appropriate sampling of these to obtain reliable thermodynamic averages.^{11,21,26,28,29} Indeed, as the accurate description of electronic degrees of freedom and the chemical reactivity of the substrate are important, ab initio molecular dynamics (AIMD) results based on density functional theory (DFT) are the only reliable benchmarks to date. The long relaxation times of interfacial water^30–32 pose, however, a significant challenge to AIMD simulations and restrict such studies to only a few selected model systems.^{14–18,21,26,33–36} As a result, many studies, e.g., on adsorption and solvation energies across different substrates and adsorbates^27,37–41 or on electrochemical reaction barriers^42–49 still rely on simplified and quasi-static interfacial water models. The water–Pt(111) interface is most likely the best studied of such systems and has been addressed by AIMD simulations^{13,14,21,29,50,51} and static low-temperature water models.^{25,50,52–54} Typical static interfacial water models are comprised of ice-like, hexagonal interfacial layers,^{20,48–50,52,54} H^down, H^up, and chain-H^down, and also involve more complex phases^52–55 such as $37×37$R25.3° and $39×39$R16.1°. The latter structure comprises five, six, and seven water rings.^52–54 Such structures are typically less prominent in room temperature AIMD studies^16,50,56 where only (predominant) local water configurations can be identified.¹⁶ This puts forward the question about the relation between theoretical results from static and dynamically sampled extended interfacial water structures and possible structural biases due to the dominant use of hexagonal water (bilayer) arrangements in applied works.^42,48,49 Furthermore, as no “standard” static water models seem available for other than the prototypical fcc-(111) surface, there is a dramatic lack of knowledge of interfacial water structures on other interfaces and substrates.

In order to fill this void, we systematically study and compare ensembles of static and dynamic water structures. We introduce first our construction protocol for obtaining a consistent set of static water structures at solid interfaces for any given surface termination, supercell size, and crystal structure. Here, we start by creating a dataset of selected, topologically different 2D water layers in vacuum. The created 2D water layers can be adsorbed on any given substrate, leveraging optimal lattice matching algorithms and consecutive geometric relaxation via DFT. Subsequently, we apply the algorithm to construct a dataset of 2D static water (2DSW) structures on low index Pt surfaces and compare the obtained structural ensemble of interfacial water with the structure of water layers from other reference works. In particular, we investigate structural similarities within quasi-2D water layers in bulk liquid water and ice, at (liquid) water–vacuum interfaces, and at Pt(111)–water interfaces using the local structure descriptor SOAP (Smooth Overlap of Atomic Positions).⁵⁷ One central aspect in this respect is the analysis of the relations between structure, adsorption energy, and work function reduction for Pt(111) with a single adsorbed water layer and as obtained from AIMD and our 2DSW ensemble. In general, the analysis reveals a bias of 2DSW structures toward three intermolecular hydrogen bonds per water and excess hydrogen bonds pointing to the substrate, which yield a maximization of the overall number of bonds within the interface in contrast to the AIMD configurations (see Secs. III B and III C). Nevertheless, we find a decent overlap between AIMD and 2DSW ensemble averages, while single selected static models exhibit large variations in their predictive quality.

Our starting point for the proposed protocol is the common understanding that the first water layer is extremely important for a variety of properties at metal–water interfaces and that H₂O–H₂O interactions dominate the total energy, at least on coinage metal low index surfaces.^8,58,59 An intuitive approach to model adsorbed water is by using water layers in vacuum as a starting guess, which directly motivates the following procedure:

• Creation of a dataset of 2D water phases in vacuum.

• Creation of the target substrate [(hkl) surface and supercell size].

• Adsorption of 2D water phases that match the substrate supercell geometry.

• Relaxation of the water adsorbate layer.

Based on DFT calculations, we initially investigated 15 different polymorphs of 2D water in vacuum consisting of 1–8 water molecules in the primitive cell as reported in the literature.^60–63 We varied lattice constants and unit cell shapes in order to get an overview of the typical 2D water orderings in vacuum. 2D water is highly flexible, and typical literature studies are in the context of confined water, e.g., in between two sheets of graphene, rationalizing size or pressure restrictions in the out-of-plane direction.⁶¹ As these are not present here, a compression mainly leads to a change in the degree of buckling and stacking of water molecules into several layers, an expansion typically to the formation of 1D water chains with varying degrees of separation. In general, energetic differences per molecule between the considered systems with different numbers of water are small, indicating that complex, large-scale water patterns are not necessary to describe the approximate energetic landscape. Furthermore, a wide variety of observed water patterns can already be achieved with only four water molecules⁶² in a rectangular primitive lattice, featuring water molecules arranged in rectangles, parallelograms, and hexagons (see Fig. 1). Note that known prototypical motives such as squares, rectangles, parallelograms, and hexagons can be understood within rectangular cells by a variation in the relative position of the molecules within the cell, as shown in Fig. 1(a). Furthermore, all observed buckling patterns (water molecules not in one plane) fall in four classes, where the water network is in a plane, exhibits horizontal or vertical chains, or forms a 2D checkerboard arrangement [see Fig. 1(b)]. In addition, the relative orientation of the water molecules can follow several patterns; the most symmetric (in-plane) representations where each water molecule can maximize the number of hydrogen bonds are drawn in Fig. 1(c).

After these observations that a wide variety of different topologies are, indeed, already possible for four molecules in rectangular cells, we decided to limit ourselves to polymorphs, which obey these restrictions, as they have several advantages: On the one hand, the relatively small water unit cells allow for a more exhaustive representation of the structural landscape due to the reduced number of degrees of freedom. On the other hand, lattice matching can be achieved by smaller increments in the (water) supercell size. As a result, it is more likely to find small substrate supercells, which is beneficial, e.g., in terms of applying the protocol to many different systems. In addition, the choice of equal lattices for all water polymorphs ensures the consistent inclusion of all the considered water topologies on any substrate.

Our final set of 2D water structures consists of eight different topologies [Fig. 1(d)] that are observed in the studied selection of 15 low-temperature polymorphs and include a large variety of position (p), buckling (b), and orientation (o) patterns. We furthermore scanned the lattice constants on a logarithmically spaced grid (99 grid points), for which we report the relative stability with respect to the optimum lattice constant in Fig. 1(e). This clarifies the small energetic differences of any of these topological variants. The first six topologies are chains in the direction of the lattice vector a such that their energy is highly sensitive to changes in that direction and largely invariant for changes in the direction of the lattice vector b. In contrast, the energy of the latter two configurations with their hexagonal geometry is sensitive to changes in both directions. A logarithmic lattice spacing is chosen as it leads to identical relative offsets of neighboring lattice vectors (here 10%), which ensures lattice matching (see below) without overlaps. Finally, we reduced this set to 20 identical lattice constant combinations for all eight polymorphs, which cover the lowest energy regions in order to construct our final dataset of 160 2D static water polymorphs.

Choosing an appropriate substrate supercell evidently depends on the quantities and systems of interest. As surface cells with maximum distances to periodic images in the lateral directions are expected to yield the least periodic boundary artifacts, we think all studies on surfaces should ideally be performed in such maximally isotropic supercells. For this purpose, we implemented an algorithm that yields (with minimal user intervention) slabs in vacuum, which are maximally isotropic in the in-plane directions for any given substrate material, (hkl) surface, and target surface area. The algorithm is outlined in the supplementary material, and a Python implementation is provided alongside this work.

With our given dataset of 2D water prototypes at the varying lattice constant and chosen substrate supercell, we use (a slightly modified version of) pymatgen^64–66 to adsorb 2D water on the substrate via leveraging its implemented lattice matching⁶⁷ functionality similarly as implemented in the pymatgen.analysis.interface module. The logarithmic lattice spacing in our dataset ensures consistency with the lattice vector mismatch constraints in pymatgen (relative length differences) and enables non-overlapping matches for neighboring lattice constants.

In addition, we include possible lateral shifts of the adsorbed water layer (scanned on a regular grid within the primitive surface unit cell and a lateral grid spacing of ∼2 Å) while also taking both possible water layer orientations into account in case the water prototype is polar in the out-of-plane direction. A self-contained Python implementation of the outlined 2DSW construction protocol is provided alongside this publication.

Using the outlined approach, we created lattice-matched 2DSW water structures on maximally isotropic substrate supercells for supercell sizes between 6 and 36 surface atoms for the low index surfaces (100), (110), and (111) of fcc-Pt. While a wide variety of possible combinations for the cell size and number of water molecules is possible, interfaces with 8 and 16 water molecules take a special role as they yield matches for all three (hkl) terminations. Here, we restrict our analysis to minimal supercell sizes with only eight water molecules. For the Pt(111) termination, which we will analyze in more detail subsequently, the algorithm automatically chooses the well-known $(12×12)R30°$ supercell and yields 80 different initial configurations for interfacial water. These are relaxed via DFT, and we compute the adsorption energetics, as well as geometric and electronic properties (eight initial structures ran into convergence problems, which were not further analyzed and discarded). The computational parameters are reported in Subsection 2 of the  Appendix, and the results for Pt(100) and Pt(110) surfaces are shown in the supplementary material.

In order to analyze the water ordering of quasi-2D water layers, we used a variety of sources, including ab initio-determined bulk ice⁶⁸ and bulk liquid water structures,⁷² liquid water slabs in vacuum,⁷² liquid water on Pt(111),^15,21 and static ice-like-layers on Pt(111) from Ref. 54 and from our 2DSW construction protocol. Table I lists the according datasets together with our chosen naming convention and the most important associated properties. Note that only the latter two sets of structures [H₂O@Pt(111) ice and H₂O@Pt 2DSW] include natively well-defined quasi-2D water layers, while all other sources also incorporate in parts or in full bulk-like water regions. However, water layers are also evidently present in 3D bulk water structures (H₂O ice, H₂O AIMD) and are, in particular, distinct and non-bulk-like for water–vacuum (H₂O/vacuum AIMD) or water–Pt interfaces [H₂O@Pt(111) AIMD(P), H₂O@Pt(111) AIMD(R)].

In order to study the structure within these layers, we create from the latter datasets first systems with water–vacuum boundaries, if needed, and subsequently extract water layers at water–vacuum interfaces. In particular, we selected snapshots from the bulk H₂O AIMD trajectory and introduced vacuum layers at random positions along the z direction while enforcing the integrity of H₂O molecules. Similarly, for H₂O ice (54 ice phases from Ref. 68), we created water slabs up to a maximum Miller index of 3 using pymatgen^64–66 again while enforcing molecular integrity. For the water–Pt AIMD trajectories [H₂O@Pt(111) AIMD(P), H₂O@Pt(111) AIMD(R)], we removed all Pt-atoms from the system, thus creating a vacuum void. Interfacial water layers were then extracted from the so-preprocessed structures by determining all interfacial water molecules only for the respective vacuum–water interfaces of interest (e.g., the original H₂O@Pt interface) via the algorithm described in Subsection 1 of the  Appendix.

The so-created structural datasets for water layers in bulk, at vacuum–water, and at Pt–water interfaces are analyzed via the use of an abstract representation for atomic environments, in particular, based on the local structure descriptor SOAP⁵⁷ as implemented in the Dscribe package⁷⁹ and as similarly performed for 3D periodic water structures.⁶⁸ Local environments are represented for each water molecule by computing a SOAP vector centered on the oxygen atom and using a cutoff radius of 4 Å, consistent with the geometrical characteristics of water. Interfacial water structures are characterized by the average over all SOAP vectors taken from each O atom in a configuration,⁷⁹ which thus also allows (global) for a structure comparison across structures that differ in the number of H₂O molecules (cf. also Subsection 3 of the  Appendix).

In order to gain an intelligible, e.g., visual, representation of the similarity between and within the different structural ensembles of extended water layers, we apply principal component analysis (PCA) to the complete dataset of average SOAP vectors. The first and second PCA components are depicted in Fig. 2, and the solid lines indicate approximate locations for each dataset as obtained via kernel density estimation.⁸⁰ Additional dataset-specific plots are provided in the supplementary material. Figures 2(a) and 2(b) and the color-coding illustrate how PCA discriminates between the average number of intermolecular hydrogen bonds (upper left to the lower right) and the average number of oxygens within the SOAP radius (lower left to the upper right). Hydrogen bonds between two water molecules are defined as in Ref. 16, where the oxygens are closer than 3.5 Å to each other and the angle between O–O–H is less than 35°. There is a broad distribution of structures in the analyzed data, ranging from rather isolated water molecules to highly coordinated ones with four hydrogen bonds to its neighbors and the bulk-ice-derived layers (green line) covering the whole range, with selected examples plotted in Fig. 2. Bulk liquid water (H₂O AIMD, red) is mainly located in the low density region at the lower left corner, while water at water–vacuum interfaces (H₂O/vacuum AIMD, violet) and [H₂O@Pt(111) AIMD(P/R), orange/brown] move upward along the diagonal, toward higher density structures. Note that the increase in density (number of nearby oxygens ⟨N_O⟩) does not dramatically alter the number of H-bonds ⟨N_O–H⋯O⟩ within the water layer, as can be rationalized from the color codings in Figs. 2(a) and 2(b) and from the reported average values provided in Table II.

For comparison, we also marked in Fig. 2 specifically the prototypical low-temperature structures⁵⁴—(hexagonal) H^up/down and $37/39$—which fall in the region of our static water dataset H₂O@Pt 2DSW, which includes natively hexagonal water layers. It is worth mentioning that the widely used H^up/down structures are observed at the extreme boundaries of the configurational space, whereby on average oxygen atoms in H^up/down structures form 0.8 more hydrogen bonds than in H₂O@Pt(111) AIMD, as shown in Table II. While there is no overlap between these and the AIMD sampled interface structures on Pt(111) (orange and brown lines), our 2DSW protocol yields at least some overlap. In particular, structures with less hydrogen bonds, which mainly consist of chain-like water arrangements, fall in the region of AIMD results (see also discussion below). These results indicate that ensemble averages based on 2DSW static water structures might approximate AIMD averages more accurately than using only selected low-temperature (e.g., hexagonal) structural models. Interestingly, Le et al. pointed out that beyond the first solvation layer, the hydrogen bonding is mainly affected by the water–water interaction rather than the water–metal interaction.¹⁶ This suggests that better results might be achieved by including more water layers in the system.

While all 2DSW water structures are provided alongside this work, we show a selection of three obtained structures in Fig. 3, namely, the lowest energy configuration after relaxation (a) and two local minima [(b) and (c)].

Interestingly, the starting structure of the least stable polymorph p3b1o3 [Fig. 3(a)] leads to the most stable relaxed structure in our set, which is characterized by hexagonal H₂O rings, in which the hydrogen network is distinguished by two chains: one constituted of only H^down oriented water molecules and the other chain comprises H₂O molecules parallel to the surface. Here, the structure is stabilized by forming the three hydrogen bonds per oxygen atom as in other ice-like water layers. Furthermore, this structure was also reported by Clabaut et al.⁵⁴ as the most stable among the hexagonal ice adlayers.

The configuration in Fig. 3(b) is also initialized from the p3b1o3 polymorph, but the start configurations differ from each other in the atomic positions of the water molecules, while Fig. 3(c) started from chain-like water structures. Both relaxed structures in Figs. 3(b) and 3(c) form chains of hexagonal H₂O rings, however, with differing water orientations. Regarding energy, both structures are only 0.03 eV/H₂O less stable relative to the lowest energy configuration [Fig. 3(a)]. Similar chain-structures have been reported for (110) metal surfaces.⁵³ 

In general, our 2DSW protocol yields a variety of hexagonal ring structures comprised of H^down, H^up oriented water molecules and water molecules parallel to the surface, similarly as shown in Fig. 1(a) at 11 ps in Ref. 51, supporting their importance for simulations in small unit cells.

For the analysis of water layers on the substrate Pt(111), we follow the previous approach and leverage the average SOAP vector as the structural descriptor, this time, however, including specifically the Pt substrate atoms in the atomic environments (see Subsection 3 of the  Appendix). Furthermore, we use a SOAP-based average distance measure $⟨D̄⟩$ to evaluate the structural similarity of a single structure with the AIMD(P) ensemble of interfacial water structures, which we took as our reference (cf. Subsection 3 of the  Appendix). The PCA map in Fig. 4 yields similar results as before and illustrates how $⟨D̄⟩$ (cf. color) is able to measure structural similarity to the reference ensemble. As expected, both AIMD simulations remain close to each other. For the H₂O@Pt(111) ice layers (diamonds), we find a clear trend in distances $D̄$⁠, namely, H^up (0.27) > H^down ≈ chain-H^down > $37$ > $39$ (0.21), where the rather large dissimilarity of the hexagonal structures is likely ascribed to the lack of disorder. Similarly as before, our 2DSW ensemble shows better agreement with the AIMD results, although some large dissimilarities $⟨D̄⟩$ are observed. Further analysis is reported in the supplementary material, where we also analyze other commonly used descriptors of (interfacial) water, such as the atomic density distribution as a function of the distance to the surface and the radial distribution functions within the water adlayers. In general, we find good qualitative agreement in the vertical density profile for the 2DSW ensemble, which is able to reproduce the bimodal oxygen peak structure at the approximately right vertical distance z_Pt–O between slab and bottom-most oxygen atoms {2.2 Å (2DSW) vs 2.0 Å [AIMD(P)]}. However, some differences are observed in the oxygen density profile, e.g., the first peak for AIMD data is wider than the peak for 2DSW, as shown in Fig. S4 of the supplementary material. This can be ascribed to the absence of the dynamic interchange of interfacial water molecules within the solvation layer and with the bulk liquid water.

To this point, we have compared structural similarity to the AIMD(P) reference trajectory; however, we have not yet analyzed how far this influences derived quantities such as average adsorption energies and work function changes, which show a strong structural sensitivity, e.g., to the distance between the metal surface and the water adlayer, as pointed out by Tripkovic et al.²⁵ Similarly, AIMD simulations are comprised of a variety of interfacial water structures due to fluctuating water orientations and ad- and desorption events, which strongly influences the work function due to interfacial charge transfer from the first water layer.^17,18,36 In this regard, Surendralal et al. reported values of absolute potentials ranging from 2.5 to roughly 7.5 eV for H₂O@Pt(111) AIMD simulations, which result in an average of 4.86 eV (see the supplementary material in Ref. 81), and it is of interest to see how far 2DSW ensemble averages relate to such AIMD-based reference results. Here, we evaluate the adsorption energy per water molecule as

where $NH2O$ is the number of water molecules at the interface and $EH2O–PtDFT$⁠, $EPtDFT$⁠, and $EH2ODFT$ are the total energy of the surface with the adsorbed water layer, the energy of a clean surface, and the energy of an isolated water molecule, respectively. The work function change ΔΦ due to the adsorption of a water film is determined by the difference in the work function between water-covered and clean slab $ΔΦ=ΦH2O–Pt−ΦPt$⁠, both of which are readily accessible as the position of the Fermi level relative to the electrostatic potential in the vacuum region. We computed the respective quantities for our 2DSW dataset but reanalyzed in the same way the AIMD(P) trajectory with all waters removed except for the first interfacial water layer. Gratifyingly, the so-computed work function reduction ⟨ΔΦ⟩ for the AIMD(P) trajectory by only the first water layer is identical to the results reported in the original study by Le et al.²¹ (−1.11 eV here vs −1.1 eV in Ref. 21). This supports once more the major role of only the first water layer.^22,23

In order to understand better the impact of different possible choices of static interfacial water structures, we studied two specific subsets of our raw 2DSW dataset, which include (i) only configurations with an energy difference ΔE ≤ 0.05 eV/H₂O relative to the lowest energy configuration and (ii) configurations with $⟨D̄⟩$ < 0.21. These two subsets are inspired by two typical assumptions in theoretical works, namely, to focus (i) on low energy structures and (ii) on structures that are close to some accurate reference (e.g., experiments or as here more accurate theoretical simulations). Table III summarizes the central properties of the considered structural datasets.

Structurally, the low energy subset (i) exhibits the highest observed number of hydrogen bonds ⟨N_O–H⋯O⟩, and the number of nearby oxygens ⟨N_O⟩ is close to 3, which likely derives from the fact that structural optimization leads to a maximized number of hydrogen bonds within the water layer and maximized number of bonds with the metal surface [cf. most stable configuration in Fig. 3(a)]. Therefore, the structures mainly consist of the hexagonal (N_O–H⋯O = 3) and some two-dimensional chain structures (N_O–H⋯O = 2.5). In contrast, the latter subset (ii), which is optimized for structural similarity with the AIMD(P) dataset, is characterized by chain structures, as shown in Fig. 3(b), where the average number of hydrogen bonds (2.23) and the average number of nearby oxygen atoms (3.41) are close to the values of the reference set (cf. Table III and Fig. 4).

For all static water ensembles, we find E_ads values that are significantly lower than the AIMD(P) average, as the latter structures were dynamically sampled with additional bonding partners in the second water layer in the original simulations, which are (artificially) removed in our analysis. Within our dataset, E_ads of the structurally most similar configuration [see Fig. 5(a)] is roughly 0.1 eV higher than the lowest energy configuration. On the other hand, our average adsorption energies fall close to the reported adsorption energies for static water models.^20,52,58,82 While the numerical differences to the AIMD averages are thus evident, a diverse set of local minima configurations as starting points for further investigations, e.g., based on short AIMD simulations or nudged elastic band calculations, might be beneficial for a better understanding of the effects of local water environments on other processes, such as solvation, diffusion, or electrochemical reactions.

In terms of the work function reduction, the better structural agreement of subset (ii) with the AIMD(P) reference dataset does not lead to better agreement in the work function reduction [⟨ΔΦ⟩ = −0.43/− 0.30 eV for 2DSW/subset (ii) vs −1.11 eV for AIMD(P)]. Indeed, the low energy subset (i) performs best with ⟨ΔΦ⟩ = −0.83 eV, which lies within the margin of errors of published AIMD reference results (−1.1 to −0.55 eV^14,21,81). Note that accurate work functions and work function reductions are relevant for applications in electrochemical contexts as they are directly related to the potential of zero charge of the electrode on an absolute potential scale.⁸³ 

In order to get a better understanding of the correlations between the structural distances $⟨D̄⟩$⁠, the adsorption energy E_ads, and the work function change ΔΦ, we plot an according analysis of the complete 2DSW dataset in Fig. 5. The beige areas in Figs. 5(a) and 5(b) correspond to our selection criteria for subsets (i) and (ii), respectively. Consistent with the described lower stability of dynamical water structures, we also observe here that the best matching structures (blue) are higher in energy [∼0.1 eV, Fig. 5(a)]. At the same time, at these E_ads values, the distance measure $⟨D̄⟩$ seems to exhibit a bimodal distribution, thus also incorporating a significant amount of structures with maximal $⟨D̄⟩$ values. The observed work function change ΔΦ exhibits a rather large spread [Fig. 5(b)], which is slightly larger at higher $⟨D̄⟩$ values. Indeed, structures with different $⟨D̄⟩$ values can still exhibit similar energies and work function change [Fig. 5(c)], indicating a high degree of complexity for the studied structure–property relationships, providing renewed evidence for the difficulty in treating SLIs, in particular water–metal surfaces with simplified, e.g., individual, static interfacial water models. As energies and work functions are not included in the PCA, the prediction can be a challenging task in some cases as here shown.⁸⁴ 

The systematic study of static interfacial water layers on Pt(111) revealed that prototypical, low energy interfacial water structures exhibit quite special orderings and lie rather at the boundaries of the observed configuration space of AIMD simulations. At variance, our 2DSW ensemble, which includes a wide variety of water topologies, does exhibit some more overlap. The main structural discrepancy is the overemphasis of a high number of in-layer H₂O–H₂O bonds for static, single layer water models, which is most prominent for the lowest energy structures (cf. Table III). These differences are, on the one hand, linked to the conceptual differences between local minimum structures and finite temperature MD trajectories and, on the other hand, to the neglect of more water layers in the present work, which can influence the configuration of interfacial water through providing donors and acceptors for hydrogen bonds.⁵¹ Some better structural agreement might be obtained already by extending/substituting the dataset of 2D water start configurations with multilayer water models. Whether an implicit solvent environment^{22,23,47,85–87} can appropriately emulate the bonding to the second water layer and whether AIMD simulations of minimal explicit–implicit hybrid models can reproduce results from fully explicit AIMD simulations remain an open question. In particular, similar results would only be expected if the implicit model would stabilize stronger non-H-bond-saturated water structures at the interfaces, e.g., H^up, which clarifies the need for such models.

Furthermore, the observed overstructuring of interfacial water might also be influenced by the choice of the substrate supercell, as the studied $12×12$ Pt(111) substrate allows for commensurate hexagonal water orderings. Our protocol thus also provides a starting point for a consistent analysis of such (artificial) periodic boundary effects, which are often discussed, but not yet well understood.

In general, while the work function change and the adsorption energy can vary dramatically for different static water structures, our 2DSW ensemble averages can provide relatively robust averages that show decent correlations with AIMD results, at least trend-wise. As many applications still necessitate the study of quasi-static water models, e.g., for evaluating reaction barriers via nudged elastic band methods, these results suggest that ensemble averages over a consistent set of static water structures might provide more robust predictions. Indeed, the application of clustering methods to our 2DSW dataset reveals a total number of 11 topologically different prototype structures (see Sec. S.5.A of the supplementary material), which indicates that such an approach would still remain computationally much more feasible than brute-force AIMD-based methods.

Finally, our approach, in particular the implementation as a Python package provided alongside this work, might serve as a good starting point for the study of low-temperature water structures (or other solvents) on other surface supercell sizes and substrate materials or as uncorrelated starting points for AIMD simulations.

See the supplementary material for (i) a discussion about the choice of supercell in the construction protocol, (ii) a description of the in-house algorithm to find water molecules at the interface, (iii) PCA maps for the water configurations by using the local SOAP descriptors, (iv) density profile and O–O radial distribution function plots for H₂O@Pt(111) systems listed in Table I, (v) an overview of the prototype structures found in H₂O@Pt(111) 2DSW, and (vi) further results for Pt(100) and Pt(110) surfaces. For access to the relevant data of this work, see the Data Availability section.

We would like to thank Professor Dr. Jun Cheng and Dr. Hendrik H. Heenen who provided the AIMD(P) and AIMD(R) trajectories, which are related to Refs. 21 and 15, respectively. Furthermore, we want to thank Ambrosio et al., Monserrat et al., and Clabaut et al. who made the AIMD trajectories and atomistic structures publicly available alongside their publications, namely, Refs. 72, 68, and 54.

N.G.H. acknowledges financial support through the EuroTech Postdoc Programme, which is co-funded by the European Commission under its framework program Horizon 2020 and Grant Agreement No. 754462, and support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy—Grant No. EXC 2089/1-390776260. This work was supported by a grant from the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) by providing computing time through the John von Neumann Institute for Computing (NIC) on the GCS Supercomputer JUWELS⁸⁸ at the Juelich Supercomputing Centre (JSC).

The authors have no conflicts to declare.

A.C.D.L and T.E. contributed equally to this work.

The data that support the findings of this study are available within the article and its supplementary material or from the original publications (cf. Acknowledgments). The newly generated H₂O@Pt 2DSW structural dataset and Python codes (analysis and determination of interface water molecules) are openly available at https://dx.doi.org/10.17617/3.7c. A Python implementation of the 2DSW protocol is available as an open-source Python package at https://github.com/computationalelectrochemistrygroup/WaterStructureCreator.

Interfacial water layers at (artificially created or naturally present) water–vacuum interfaces are selected by first determining all interfacial O atoms and then adding all chemically bound H atoms (bond cutoff = 1.2 Å). Interfacial O atoms are determined via the following algorithm that is inspired from the definition of quantum mechanical cavities in the context of implicit solvation methods:^89,90

• Define the water-filled region Ω by superimposing soft spheres localized on the atoms with given radius r₀ and that smooth decay from a value of 1 within the sphere to 0 outside it. Evaluate the so-obtained filling function Ω on a numerical grid and renormalize all values Ω > 1 to 1.

• Estimate the surface of Ω by computing the numerical derivative S = ∇Ω, and define a grid point i as an interfacial grid points i_S whenever

  $‖SiS‖≥τmaxi‖Si‖,τ<1.$

  (A1)

• Interface atoms are given by the set of closest oxygen atoms to the real space positions ${riS}$ of all interface grid points.

More details, e.g., on the numerical implementation in periodic boundary conditions and chosen numerical parameters (e.g., r₀, τ) are provided in the supplementary material.

Periodic DFT calculations are performed with the Quantum ESPRESSO package⁹¹ (PWscf) and the Perdew–Burke–Ernzerhof-generalized gradient approximation (PBE-GGA)⁷⁴ for the exchange–correlation energy functional. All atoms are represented by ultrasoft pseudopotentials from the open-source GBRV library⁹² with density and wave function cutoffs of 360 and 45 Ry, respectively. The water films, modeled with eight H₂O molecules, are placed on side of the surface, and a dipole correction as implemented in ENVIRON⁸⁵ is applied. The Pt surfaces are simulated by using a (3 × 3) supercell with three layers for Pt(100), a (3 × 2) supercell with four layers for Pt(110), and a $(12×12)R30°$ supercell with three layers for Pt(111). The two bottommost metal layers are fixed at their ideal bulk positions. Slabs are separated by ≈17 Å. Geometry optimizations are carried through until the forces on all relaxed atoms were smaller than 0.1 eV Å⁻¹. Brillouin zone integrations are performed using Γ-centered Monkhorst–Pack meshes of (4 × 4 × 1) and a Marzari–Vanderbilt smearing of 0.02 Ry.

Single point calculations of the interfacial water on Pt(111) for the AIMD simulation from the work of Le et al.²¹ were performed with the same settings, except that a (2 × 2 × 1) k-point grid was used. The symmetric setup of the original AIMD was split in two asymmetric setups with three layers of platinum each. The bottom half was rotated by 180° around the y axis to have the same orientation for both sides. No further manipulations were performed on the cell or positions of the atoms.

We selected the local structure descriptor SOAP⁵⁷ as implemented in the Dscribe package⁷⁹ and as similarly performed for 3D periodic water structures⁶⁸ for the representation of the atomic environment. The local environment $X$ around an atom is modeled by summing over Gaussians centered on each atom i in the structure within a cutoff, yielding the density $ρX(r)$⁠. After expansion of the density in a basis of orthonormal radial functions $gb|r|$ and spherical harmonics $Ylmr̂$⁠, the power spectrum $pXb1b2l$ can be expressed in terms of expansion coefficients c_blm as

The c_blm coefficients form the components of an abstract SOAP descriptor vector, which is conveniently used for measuring structural similarity throughout this work.

For the analysis of the internal structure of water layers, we evaluate the SOAP vectors centered on the oxygen atoms, including neighboring H and O atoms up to a cutoff radius of 4 Å. For the analysis of water layers on top of Pt(111), we also include the substrate Pt positions in the environment within the sphere defined by the cutoff. In both cases, crossover terms between different atomic species in the power spectrum were allowed to include the information of hydrogen and, in the later case, platinum as well. A description of the entire system was obtained by using the inner average of the SOAP vectors.⁷⁹ The hyperparameters for computing the SOAP descriptors used throughout this work are n_max = 8, l_max = 8, σ = 0.3, and r_cutoff = 4 Å.

For comparing the structural similarity of single structures A with a reference structure B, we use the average distance $D̄$⁹³ as given by the average kernel $K̄$⁹³ and defined by

where N and M are the number of atoms in structures A and B, respectively. The so-obtained distance measure $D̄$ is subsequently averaged over all structures in the reference dataset to obtain an average structural dissimilarity measure $⟨D̄⟩$ with the reference ensemble, which we call shortly average SOAP distance, here. In the present work, the reference ensemble was created by taking snapshots every 0.1 ps from the AIMD(P) trajectory²¹ that was provided by Cheng.