Theoretical investigations of a platinum–water interface using quantum-mechanics-molecular-mechanics based molecular dynamics simulations

R. P. Hardikara, Unmesh Mondalb, Foram M. Thakkarc, Sudip Roy§c and Prasenjit Ghosh*ad
aDepartment of Physics, Indian Institute of Science Education and Research (IISER), Pune 411008, Maharashtra, India. E-mail:;
bDepartment of Chemistry, Indian Institute of Science Education and Research (IISER), Pune 411008, Maharashtra, India
cShell Technology Centre, Bangalore, India
dCentre for Energy Science, Indian Institute of Science Education and Research (IISER), Pune 411008, Maharashtra, India

Received 24th June 2019 , Accepted 5th October 2019

First published on 7th October 2019

Pt–water interfaces have been of immense interest in the field of energy storage and conversion. Studying this interface using both experimental and theoretical tools is challenging. On the theoretical front, typically one uses classical molecular dynamics (MD) simulations to handle large system sizes or time scales while for a more accurate quantum mechanical description Born Oppenheimer MD (BOMD) is typically used. The latter is limited to smaller system sizes and time-scales. In this study using quantum-mechanics-molecular-mechanics (QMMM), we have performed atomistic MD simulations to have a microscopic understanding of the structure of the Pt–water interface using a system size that is much larger than that accessible when using BOMD simulations. In contrast to recent reports using BOMD simulations, our study reveals that the water molecules typically form two distinct layers above the Pt-surface before they form bulk like structures. Further, we also find that a significant fraction of the water molecules at the interface are pointed towards the surface thereby disrupting the H-bond network. Consistent with this observation, the layer resolved oxygen–oxygen radial distribution function for the water molecules belonging to the solvating water layer shows a high density liquid like behaviour even though the overall water behaves like a low density liquid. A charge transfer analysis reveals that this solvating water layer donates electrons to the Pt atoms in contact with it thereby resulting in the formation of an interface dipole that is pointing towards the surface. Our results suggest that, using QMMM-MD, on one hand it is possible to study more realistic models of solid–liquid interfaces that are inaccessible with BOMD, while on the other hand one also has access to information about such systems that are not obtained from conventional classical MD simulations.

1 Introduction

With the ever-increasing demand for energy, there is a need for alternative sources of clean and renewable energy. One such source of energy is hydrogen. In low-temperature fuel cells (which convert chemical energy in hydrogen to electricity) and electrolyzers (which store energy in the chemical form by splitting water to form hydrogen), the Hydrogen Oxidation Reaction (HOR) and the Hydrogen Evolution Reaction (HER) are two of the key processes. For these reactions, though Pt is still the most efficient catalyst, there are efforts to replace it with cheaper materials which exhibit efficiencies similar to Pt. Yet, to achieve this, one needs to have a microscopic understanding of the structure of the metal–water interface and the nature of the interactions between the two components, which, to the best of our knowledge, are not yet well understood.

Experimentally accessing this information is challenging and in some cases inaccessible (e.g. local structure of the interface). The typical experimental tools that are used are grazing incidence X-ray scattering,1 which gives crystallographic information and vibrational spectroscopies such as infra-red (IR) and Raman, together with non-linear second order sum frequency generation that gives vibrational information of the water molecules at the interface.2,3 Recently, there were reports on experiments performed on the Au(111)/water interface using X-ray absorption, which provides information regarding the local structure and chemical environment of water molecules at the interface.4 Further, by using ultra-low temperature (5 K) STM Guo et al. resolved the internal structure of water at the interface of NaCl(001) films supported on a Au(111) substrate.5 However, most of these methods rely on theoretical modelling of the interface to interpret their data.

On the computational front, studying these systems is equally challenging because there is not only a need for an accurate description of the potential energy surface (PES) but also for adequate treatment of the thermal fluctuations given the statistical nature of liquid water. While classical molecular dynamics (MD) simulations allow one to access large system sizes and time scales, ensuring that the thermal fluctuations are treated properly and finding an accurate description of the PES is challenging. Moreover, in a classical MD simulation, the PES on which the atoms evolve are determined by parametrization. Hence, every time the interface changes, one needs to generate a new set of parameters. Further, as mentioned in the first paragraph, the HER/HOR happens at the interface and involves bond breaking and bond formation that are not captured by conventional classical force fields. These problems can be circumvented by using BOMD simulations. However, because of the enormously large system size in this case, the time scale and system size accessible to these simulations are limited, thereby questioning the statistical accuracy of the results.

Amongst the different metal–water interfaces, Pt–water is the most studied using both computational (classical and BOMD) simulations and experiments. At room temperature, these are expected to show features between those of the ordered bilayer or ice structure6 and the labile bulk.7,8 Classical MD simulations by Willard et al. have shown the hydrophobic nature of Pt surfaces.9 Further, they have also highlighted the role of extended time and size dynamics in water exchange close to Pt. Unfortunately, such large time scales are not accessible in BOMD simulations. In contrast, in most BOMD simulations, the water–Pt interfaces have been primarily modelled in two ways: (a) a few layers of water (typically four water layers) on a Pt slab followed by vacuum on the other side10 or (b) a couple of layers of water with symmetric Pt slabs on either side.11 While in the former case the number of water layers may not be enough, the latter describes the case of confined water.

Hence to address these issues, in this work we have studied the Pt–water interface using the hybrid Quantum-Mechanics-Molecular-Mechanics (QMMM) MD simulations (QMMM-MD) where the Pt atoms and a few layers of water molecules near the interface are treated quantum mechanically and the rest of the water molecules are treated classically. We have not only used an interface whose lateral dimensions are larger than the ones reported in the literature but also a significantly thicker water layer with a vacuum on the other side. This model ensures that we are simulating an interface between bulk water and the Pt-surface. The rest of the paper is arranged as follows: Section 2 provides details of the MM and QMMM-MD simulations. The results are presented and discussed in Section 3 and finally summarized in Section 4.

2 Computational details

In this section we describe the computational protocol that has been used to set-up a Pt–H2O simulation cell for the MM calculations followed by details of the QMMM-MD simulations.

2.1 MM setup

In order to obtain the starting configuration for our QMMM-MD simulations, we have performed classical MD simulations for the interface. Water was modelled with the modified SPC/E force field.12 Typical MD calculations employ constraints to keep the water molecule rigid, resulting in accurate reproduction of macroscopic observables such as density (ρ) and radial distribution functions (g(r)). However, since we intend to avoid using any constraints at the QMMM-MD level of theory, we incorporate the effect of rigid water molecules by modifying the harmonic potentials of the form image file: c9cp03558c-t1.tif where x = b,θ for the OH bond and the H–O–H angle, respectively. We use Kb = 1600 kcal mol−1 Å−2 and Kθ = 110 kcal mol−1 rad−2, which reproduces the same structural parameters for bulk water as that of the SPC/E model. Details regarding the tests performed with these modified forcefield parameters for bulk H2O are provided in the ESI. The Pt atoms interact amongst themselves and with the water molecules at the interface through Lennard-Jones potential. The values for these parameters were adopted from MD studies by Sphor et al., where a lamina of water is constrained between two Pt surfaces.13 The resultant geometry of Pt layers at the MM level gives a Pt–Pt bond length of 2.76 Å and an average interlayer separation of 2.26 Å.

The Pt(111) surface is modelled using a slab of 4 Pt layers. Instead of the primitive hexagonal surface unit cell, we have considered the orthorhombic surface supercell containing two atoms per surface unit cell with cell dimensions of 2.77 Å × 4.80 Å. For our calculations we have used (5 × 3) and (7 × 4) unit cells containing 30 and 56 Pt atoms in each layer. In the rest of this paper we call the former model-A and the latter model-B.

Before building the H2O–Pt structure, we begin with equilibrating geometries of bulk water using our parameterized forcefield parameters. To achieve this we fill up a cubic box (sides 20.11 Å) with 268 water molecules. Separate simulations are carried out first by equilibrating the density of water at a constant temperature and pressure (NPT ensemble, at atmospheric pressure and a temperature of 300 K) followed by a canonical ensemble (NVT). The average density of water with the forcefields we consider was found to be 0.989 g cc−1 (Fig. S1(d), ESI). This value matches well with the rigid model using SPC/E parameters.14 To construct model-A, we cut an orthorhombic structure of bulk water to fit over the Pt surface. This gives 180 water molecules for model-A. For model-B, since we were also interested in a thicker water slab, we constructed a supercell of this box by doubling the cell dimension along the z-direction. This was further equilibrated. From this equilibrated supercell, orthorhombic slab structures are then cut from the bulk water to fit over the Pt surface. This results in 480 water molecules in model-B. Further, since CP2K uses periodic boundary conditions along all the directions and the system is periodic in reality only along the xy direction, we have used a vacuum of 29 and 47 Å for model-A and model-B, respectively, along the z-direction, to minimize the spurious interactions between the periodic images along the z-direction. A distance of about 3 Å is maintained at the interface of Pt–H2O. An empty space at the edges of the box is also maintained to get rid of any periodic boundary issues for the water molecules, thereby achieving proper connectivity of water molecules in the XY-dimensions. For model-B the long-range interactions are truncated at a cutoff distance set to half of the min[Lx,Ly], where Lx and Ly are the cell dimensions in the X and Y directions, respectively. Since the min[Lx,Ly] is less than the conventional cutoff radius used in bulk water, we set 9 Å as the cutoff distance for model-A. This corresponds to the minimum distance required for long-range interactions in the case of bulk water. These long-range electrostatic interactions over the periodic boundary are incorporated using the smooth particle mesh Ewald (SPME) scheme.15 Temperature is controlled using a CSVR16 thermostat with a time constant of 1 ps. Very long equilibration runs of ∼9 ns are performed at the MM level to get a well equilibrated geometry for the Pt–H2O interface structure.

2.2 QMMM-MD setup

The starting geometry for the QMMM-MD setup is chosen from a random snapshot of the equilibrated region of the MM trajectory. The total system is divided into two parts: the Pt slab and a layer of water molecules above it constitute the QM region, while the rest of the system belongs to the MM region.

To describe the interactions between the molecules in the MM region we have used the same forcefields that were formalised in Section 2.1.

For the QM region we use the QUICKSTEP17 module as implemented in the CP2K package, which is an implementation of density functional theory (DFT). The electronic wavefunctions were expanded in a double-ζ valence Gaussian basis set (DZVP-MOLOPT)18 while the charge density was expanded in an auxiliary plane wave basis, with a kinetic energy cut off of 300 Ry. The electron–electron exchange–correlation potential is approximated by the Perdew–Burke–Ernzerhof (PBE) parametrization of the generalized gradient approximation (GGA).19 It is known that the PBE functional gives rise to over structured water and that the radial distribution function (g(r)) does not match well with experiments, unless dispersion correction is included.20 In our calculations, we have included semi-emperical Grimme's van der Waals (DFTD3) corrections to account for the dispersion interactions.21 The electron–ion interactions were described using the Goedecker–Teter–Hutter (GTH) pseudopotentials.22,23

For QMMM-MD input parameters, a list of indices of atoms comprising the QM region (called the QM core) is chosen (shown with VDW representation on the right-hand side panel of Fig. 1(a)) and includes atoms from the four Pt layers and the first layer of water near the interface. Such bifurcation into different regions often results in discontinuities of forces at the QM and MM boundaries. To ensure a smooth transition from the QM to the MM region, we used the adaptive buffer force method.24–27 Accordingly, we provided radius ranges rQM and rBUF as the input. The code then selects atoms based on the distance criteria (rQM) from the edges of the water molecules in the QM core list. This region is referred to in the literature as the QM extended region (QMEXT). In our setup rQM is set to 2.0–2.5 Å from the edges of the QM core list of atoms, represented as a Licorice model in Fig. 1(a). Similarly, the distance of rBUF is calculated from the edges of the QM extended region. To understand the optimal buffer size, we first perform calculations with various buffer sizes and evaluate the difference in forces (mean and max) on the oxygen atoms in the QM core list from a full-QM calculation and QMMM-MD calculations with different buffer sizes. As shown in Fig. 1(b) a clear indication of convergence of mean and max force curve is seen beyond a buffer size of 3 Å, indicating converged forces on the atoms in the QM core list. Hence, we use a buffer size of 4 Å for further QMMM-MD calculations. An electrostatic embedding scheme using the Gaussian Expansion of Electrostatic Potential (GEEP)25,28,29 algorithm employed in the CP2K code is used for coupling the QM and MM part. The advantage of the GEEP algorithm is that it can handle both periodic (PBC) and open boundary conditions and scales linearly25 for calculation of electrostatic interactions.

image file: c9cp03558c-f1.tif
Fig. 1 (a) Schematic representation of different regions in a QM/MM model setup. (b) Convergence of the difference between mean and maximum force as a function of buffer size.

Using the above set up we have performed 10 ps of QMMM-MD simulations under the NVT ensemble at 300 K. The temperature is maintained through the use of an adaptive Langevin thermostat.30 Integrations were performed with a timestep of 1.0 fs. The system equilibrated in 2.0 ps. The rest of the 8 ps of the trajectory was used for further analysis. At this point we would like to note that the dynamics of the interfacial water molecules is slower than those in the bulk.9 Hence it is not guaranteed that accumulating 8 ps of trajectory is sufficient. However, performing simulations longer than those reported in this paper was not possible for us due to limited computational resources. In contrast, the 9 ns long MM based MD simulations are better converged. However, it cannot provide us with information which is accessible using QM simulations, for example, charge transfer. Keeping this in mind, we describe below our results and discuss their implications. Throughout the manuscript, we have also compared our QMMM-MD results with literature reports (experimental, as well as classical and ab initio based MD simulations). In Section C of the ESI, we have also provided a comparison of the results obtained from our QMMM-MD simulations with those obtained from our MM simulations, and we find that our results are in reasonably good agreement.

3 Results and discussion

In this section we present and discuss the structural properties of the Pt–water interface at room temperature as revealed from our QMMM-MD simulations using two different system sizes.

Throughout the simulation we did not find any water molecules leaving the water layer in contact with the Pt-surface to the bulk region. For model-A we have about 26 water molecules on the Pt surface resulting in a coverage of 0.87 ML [1 monolayer (ML) = 1 water molecule per surface Pt atom] while for model-B we have about 47 water molecules on the Pt surface resulting in a coverage of 0.84 ML. We note that these coverages are larger than those considered by Bellarosa et al.10 Fig. 2(a) and (b) show the side view of representative snapshots from the trajectories for model-A and -B, respectively. Careful visualization of the two snapshots shows that the water is more structured in model-A, forming four distinct layers, compared to the water in model-B. Similar structuring was also observed by Ryczko et al.11 This is further corroborated by plotting the layer resolved density of water as a function of distance from the Pt surface (Fig. 3(a) and (b)). For model-A the density of water is much larger than 1 g cc−1 (density of bulk water) until reaching about 12.5 Å away from the surface (Fig. 3(a)), while for model-B, the water density decays rapidly as one moves away from the surface and attains the bulk value when the water molecules are about 7.5 Å above the Pt surface (Fig. 3(b)). We note that our results for model-B are in reasonably good agreement with those reported by Limmer et al. where they used image file: c9cp03558c-t2.tif a Pt unit cell and a 40 Å thick layer of water.31 This suggests that model-A is not sufficiently large enough to represent the water–Pt interface correctly. Hence in the rest of this paper, we will discuss our results obtained using model-B only.

image file: c9cp03558c-f2.tif
Fig. 2 Side view of the last snapshots for model-A (a) and model-B (b) at the QMMM-MD level. Yellow, blue and red water molecules denote the first, second and third water layers above the surface while the green water molecules denote the liquid-like region.

image file: c9cp03558c-f3.tif
Fig. 3 Density of water in g cc−1 for the water layers over the Pt surface for (a) model-A and (b) model-B.

Structural parameters

One can ask two interesting questions at this point: (a) what is the distance of the water layer closest to the surface from the top-most layer of the Pt surface? and (b) how does the density of the water layer change as one moves away from the surface? To answer these we have computed the planar average of the density of water molecules as a function of the distance from the Pt surface using the Density Profile Tool plugin32 in VMD that evaluates 1D projection of density by applying the following formula:
image file: c9cp03558c-t3.tif(1)
where, Lx and Ly are the box lengths as described earlier, Δz is the bin size in the z direction, mi indicates the mass density and δb(zi) defines whether the atoms are within the spatial region of bΔzz < (b+ 1)δz, b being the slab index. The density in our case is calculated by using a bin size of 0.33 Å. The same is plotted in Fig. 3. From Fig. 3(b) we can see that the highest peak in the density profile is about 2.6 Å away from the Pt surface. We note that our results are in excellent agreement with recent ab initio molecular dynamics calculations by Sakong et al.33 where they found a large peak at about 3 Å away from the surface Pt layer in the planar average of the distribution of O atoms for the Pt(111)/water interface. The density is about 3.5 times larger that that observed in bulk water. We note that similar features are found in other studies of this interface as well as other reports on different metal–water interfaces; for example Velasco-Velez et al. found that for the case of Au(111)–water interfaces this forbidden region is about 2.5 Å away from the Au surface. However, our results are in variance with those obtained by Ryczko et al. from their BOMD simulations at 330 K.11 In addition to the largest peak in the density profile at about 3.5 Å, they also find a small peak at about 2 Å away from the surface that they have attributed to the water molecule adsorbed onto the Pt surface. From our simulations we do not find any water molecules adsorbed on Pt(111). In addition to this highest peak at 2.6 Å, we also find a second peak at about 6 Å from the surface. Following this peak, the density rapidly decays and hovers around the bulk value till one reaches the water–vapour interface. At the water–vapour interface there is a slight increase in density and then it rapidly decays to zero. We do not observe any evaporation of water molecules during the length of our simulations. Thus, from the density profile obtained from our simulations, the high-density region (HDR) extends up to about 6 Å from the surface. Further we also observe some microscopic variations of the density around the bulk water region. We believe that this variation is due to the small size of the data set used to compute the planar average of the density in our QMMM simulations because for the MM simulations, this density profile in the bulk water region becomes flat (Fig. S3 in ESI).

To have a better understanding of the local structure, we have plotted the radial distribution function (RDF) for oxygen, gO–O(r), in Fig. 4. We have shown both the total and layer resolved RDFs and compared the former with that of bulk water. The layers are defined as shown in Fig. 2. The details regarding the calculations of g(r) are given in Section B of the ESI. The total RDF for our system, the magenta line in Fig. 4, shows two peaks, a sharp one at about 2.75 Å and a broad one peaked at about 4.67 Å representing the first and second coordination spheres respectively. We note that these are characteristic features of the low density liquid (LDL) that is typically observed in bulk water7,8,34 as shown by the O–O RDF plot for bulk water (obtained from MM simulations), the black dashed line in Fig. 4. In addition to that we find another peak at around 7 Å that corresponds to a weak coordination sphere. We note that similar features are also observed from the BOMD simulations of Bellarosa et al.10 However in their case the water was confined between two Pt electrodes. Bellarosa et al. found similar features only when the thickness of the water layer exceeded 1.4 nm. These results indicate that the thickness of the water layer considered in model-B is sufficient to capture all the features.

image file: c9cp03558c-f4.tif
Fig. 4 O–O g(r) for all the water molecules obtained from the QMMM-MD trajectory. Also shown is the layer resolved g(r) for the first layer and the liquid region (LL). The total g(r) is compared with that for bulk water obtained from classical MD simulations (represented by the dashed black line). The figure in the inset shows the zoomed in region of the first peak. The dashed vertical lines indicate the peak position for each of the RDFs.

In contrast to the total RDF, the layer resolved RDFs show that the behaviour of liquid water changes as one moves away from the Pt surface. For example the O–O RDF of the water layer solvating the Pt surface (first layer) is significantly different from that of LDL. The first peak position is slightly shifted to a lower value of r compared to that observed in LDL (as shown in the inset of Fig. 4). The second peak observed at 4.67 Å in LDL, which corresponds to the second coordination shell of bulk water, collapses. Rather, we find a hump close to the first peak at 3.4 Å (Fig. 4). Further, the third peak at about 7 Å which corresponds to a weak coordination shell for LDL moves inwards to 5.69 Å and is more well defined. We note that these features in the RDF are characteristics of a high density liquid.35 Similar features are observed for water thin films (about 2.4 ML thick) on a BaF2(111) surface under ambient conditions.36 As we move away from the surface, the LDL-like behaviour of water is recovered. Our simulations suggest that bulk water on a Pt surface behaves in a heterogeneous way. Additionally, we find that for this layer the peak height is also enhanced compared to that from bulk water. These features in the g(r) indicate that the water in this solvating layer is more structured and is consistent with the enhanced density we find in this layer (Fig. 3). Though our results are consistent with the layer resolved density profile shown in Fig. 3 (we note that this agrees excellently with literature reports), we would like to point out that these are based on 10 ps simulations. Unfortunately, to the best of our knowledge, similar analyses are not reported in the literature on Pt or other metal surfaces. Hence it would be desirable to perform longer simulations, which we believe is possible using the QMMM-MD protocol that has been described in this paper.

Effect on the –OH covalent bond

We have also monitored how the covalent –OH bond lengths are affected in the presence of the Pt surface. The layer resolved probability distribution function (PDF) of the –OH bond lengths is plotted in Fig. 5. For the layer that is in contact with the surface, we find that the average bond length is about 0.99 Å suggesting that the bond lengths are slightly shortened in comparison to 1.01 Å observed in the bulk liquid region. As we move away from the surface, we find that the average bond length gradually converges to that observed for bulk liquid water.
image file: c9cp03558c-f5.tif
Fig. 5 Layer dependent distribution of –OH bond lengths.

Orientation of water molecules

Since we have got the macroscopic picture of the interface, now we proceed to perform further analysis to elucidate more on the nature of the interaction of the water molecules with the surface. A water monomer can interact with the surface either through the lone pair on the O atom or through the H atoms. Previous first principles calculations show that on Pt(111) a single water molecule is almost flat, with the geometric dipole tilted slightly away from the surface.37 However, when we have liquid water, there is a competition between the metal–water interactions and the water–water interactions. The balance of these two interactions will not only determine the structure of the water layer that is in direct contact with the surface but also of those that are away from the surface till the bulk liquid regime is reached. To understand this and also to understand how far away from the metal surface the metal layers affect these interactions, we have plotted the layer resolved probability distribution function (PDF) of the orientation of the geometric dipole moment of the water molecules with respect to the surface normal. The geometric dipole is defined as the vector that bisects the H–O–H angle and lies in the plane formed by the H, O and H atoms of each water molecule. The same for the first three layers and the liquid-like region is shown in Fig. 6. The angle θ is defined as the angle between the geometric dipole of a water molecule and the surface normal pointing towards the surface (inset of Fig. 6). Hence cos[thin space (1/6-em)]θ = 1(−1) implies that the water molecule is oriented in such a way that the H atoms point towards (away from) the surface, while cos[thin space (1/6-em)]θ = 0 implies that the molecules are parallel to the surface.
image file: c9cp03558c-f6.tif
Fig. 6 Probability distribution function of the angle between the geometric dipole vector of the water molecules in the water layer at the interface (1st layer, a), those above it (2nd layer, b) and (3rd layer, c), and in the bulk liquid-like region (LL, d) with the surface normal. The inset in (a) shows schematically how the angle is defined.

To study the convergence of the PDF, the 8 ps of the QMMM-MD trajectory was split into two parts: “first 4 ps” and “second 4 ps” (Fig. 6). We find that these PDFs agree qualitatively. From Fig. 6(a), based on the last 8 ps of our trajectory, we found that the probability distribution for the water layer just above the surface is almost bimodal with a broad peak at θ ≃ 101° and a relatively sharper peak at θ ≃ 53°. Additionally we find a shoulder at around θ ≃ 127°. This suggests that most of the water molecules are either pointing towards or away from the surface. Further, we also find that there is a sizable fraction of water molecules that lie almost parallel to the surface. The overall features are in reasonably good agreement with those reported by Ryczko et al.11 In contrast, in the second and third layers (Fig. 6(b) and (c) respectively), the PDFs shows a multimodal distribution. In the liquid like region (Fig. 6(d)) we find an almost uniform distribution. These results suggest that as we move away from the surface, the randomness in the orientation of the water molecules increases, tending towards bulk water-like behaviour. We note that our results agree reasonably well with those reported by Steinmann et al.38

H-Bond distribution

To understand how the connectivity of the H-bond network in water changes in the presence of a surface as we move from the water layer in contact with the surface towards bulk, we have calculated fraction of zero, single and double donors amongst the water molecules in the individual layers. The definition for the donor categories has been adopted from ref. 4. A H-bond is present if the O–O distance between the oxygen atoms of two water molecules lies within a radius of 3.5 Å from that of the O atom of the reference water molecule and if the O–O–H angle is within 35°. Accordingly, a double donor (DD) is a water molecule in which both the hydrogens are involved in hydrogen bonding; a single donor (SD) is one in which only one hydrogen is involved in the hydrogen bonding and zero donor (ZD), where the hydrogens are not connected with other water molecules.

To find out how similar or different the distributions of ZD, SD and DD are in the layers that are interacting with the surface compared to those that are away from the surface, we have plotted the layer resolved donor distribution statistics for the solvating water layer, one right above it and the bulk liquid like region. The same is shown in Fig. 7. We find that while the number of ZD in the different layers is similar, the number of SD and DD in each layer changes as we move away from the surface. While the water layer that is interacting with the surface has the largest number of SD, the number of SDs decrease as we move away from the surface. Accordingly, the number of DDs in the water layers also increases as we move away from the surface. The presence of large numbers of SD in the first layer indicates disruption of the H-bond network. This is in accordance with the gO–O(r) plot for this layer (Fig. 4) where we observe a collapse of the second coordination shell whose origin in the case of bulk water is due to the H-bond network.

image file: c9cp03558c-f7.tif
Fig. 7 Layer resolved zero donors (ZD), single donors (SD) and double donors (DD). Orange, blue and green bars represent the first layer, second layer and liquid like layer respectively. For definition of these regions see Fig. 2.

Charge analysis

Since these interfaces are used for electrochemical purposes, it is important to understand the electrostatics at the interface. For this purpose we have computed the charge transfer between the water–Pt(111) interface. The charge transfer within DFT is computed by computing the charge density of the total system (Pt(111) and water, ρPt+H2O), the charge density of the Pt(111) slab (ρPt) and the charge density of the water layer (ρH2O). The Pt(111) surface and the water layer should be in the same geometry as they are in the Pt(111)–water interface. The latter two are computed in the same geometry as they are in the Pt/water system. Then the charge transfer (Δρ(r)) is given by:
Δρ(r) = ρPt+H2OρPtρH2O (2)

The planar average of Δρ(r) (image file: c9cp03558c-t4.tif), where A is the area of the box in the xy plane, is then plotted as a function of distance along the direction of the interface normal. A positive value of Δρ(r) indicates accumulation of charge, while a negative value indicates depletion of charge. However, for our case, there are two types of water molecules, namely classical and quantum. Hence for ρPt+H2O and ρH2O the charge density is computed including only the water molecules in the QM region (both core and buffer) in the presence of the electrostatic potential generated by the charge distribution of the charges on the atoms comprising the water molecules in the MM region. We note that while computing the charge distribution for the QM water molecules, it is important to consider the effect of the electrostatic potential due to charges in the MM region, otherwise there will be spill-over of charges at the QM/MM interface giving rise to spurious charge transfer.

The time averaged plot of Δρ(z) (averaged over 30 snapshots) is shown in Fig. 8. We find that the solvating water layer transfers about 1.56 electrons to the surface Pt atoms. This charge transfer induces the formation of an interfacial dipole in the direction perpendicular to the surface. Similar charge transfer had been observed by Sakong et al. from their ab initio MD simulations.33 Fig. 8(c) shows Δρ(r) for a representative snapshot.

image file: c9cp03558c-f8.tif
Fig. 8 (a) The probability distribution of the z-coordinates of the Pt atoms and the O atoms of water molecules in the QM region. (b) The planar average of Δρ(r) in the direction normal to the surface. (c) Isosurface of charge transfer at an isolvalue of 0.017 e Å−3. Green lobes represent charge accumulation and violet charge depletion. The grey, red and white spheres denote Pt, O and H atoms respectively.

4 Conclusions

In summary, we have performed a detailed study of the Pt–water interface using QMMM-MD simulations. QMMM-MD allowed us to use a larger system which enabled us to obtain a more realistic model of the Pt/water interface. Our simulations show that the bulk macroscopic properties are more or less recovered at a distance of about 7.5 Å from the surface. The solvating water layer behaves as a high density liquid. At the interface, the orientation of the water molecules follows a bimodal distribution with the water molecules pointing either towards the surface or away from it, thereby resulting in a disrupted H-bond network as is evident from the large number of single donors present in this layer. However, this bimodal distribution is lifted as we move away from the surface. Moreover, we also observe a charge transfer from the solvating water layer to the Pt surface, resulting in an interface dipole in the direction normal to the surface. Our results suggest that by using QMMM-MD simulations, not only can we study more realistic models of a liquid water interface, for which, typically classical MD simulations are used, but also have access to information related to charge transfer at interfaces which can only be obtained from a quantum mechanical description of the interface.

Conflicts of interest

There are no conflicts to declare.


RPH would like to thank Shell and IISER Pune for support. PG would like to acknowledge DST-SERB, Govt. of India, Grant No.: EMR/2016/005275.

Notes and references

  1. P. Fenter, L. Cheng, S. Rihs, M. Machesky, M. Bedzyk and N. Sturchio, J. Colloid Interface Sci., 2000, 225, 154–165 CrossRef CAS PubMed.
  2. K.-I. Ataka, T. Yotsuyanagi and M. Osawa, J. Phys. Chem., 1996, 100, 10664–10672 CrossRef CAS.
  3. Y. R. Shen and V. Ostroverkhov, Chem. Rev., 2006, 106, 1140–1154 CrossRef CAS PubMed.
  4. J.-J. Velasco-Velez, C. H. Wu, T. A. Pascal, L. F. Wan, J. Guo, D. Prendergast and M. Salmeron, Science, 2014, 346, 831–834 CrossRef CAS PubMed.
  5. J. Guo, X. Meng, J. Chen, J. Peng, J. Sheng, X.-Z. Li, L. Xu, J.-R. Shi, E. Wang and Y. Jiang, Nat. Mater., 2014, 13, 184 CrossRef CAS PubMed.
  6. J. Carrasco, A. Hodgson and A. Michaelides, Nat. Mater., 2012, 11, 667–674 CrossRef CAS PubMed.
  7. A. Soper, Chem. Phys., 2000, 258, 121–137 CrossRef CAS.
  8. J. A. Sellberg, C. Huang, N. D. AU McQueen, T. A. Loh, H. Laksmono, D. Schlesinger, R. G. Sierra, D. Nordlund, C. Y. Hampton, D. Starodub, D. P. DePonte, M. Beye, C. Chen, A. V. Martin, A. Barty, K. T. Wikfeldt, T. M. Weiss, C. Caronna, J. Feldkamp, L. B. Skinner, M. M. Seibert, M. Messerschmidt, G. J. Williams, S. Boutet, L. G. M. Pettersson, M. J. Bogan and A. Nilsson, Nature, 2014, 510, 381 CrossRef CAS PubMed.
  9. A. P. Willard, D. T. Limmer, P. A. Madden and D. Chandler, J. Chem. Phys., 2013, 138, 184702 CrossRef PubMed.
  10. L. Bellarosa, R. García-Muelas, G. Revilla-López and N. López, ACS Cent. Sci., 2016, 2, 109–116 CrossRef CAS PubMed.
  11. K. Ryczko and I. Tamblyn, Phys. Rev. B, 2017, 96, 064104 CrossRef.
  12. H. J. C. Berendsen, J. R. Grigera and T. P. Straatsma, J. Phys. Chem., 1987, 91, 6269–6271 CrossRef CAS.
  13. E. Spohr and K. Heinzinger, Chem. Phys. Lett., 1986, 123, 218–221 CrossRef CAS.
  14. P. Mark and L. Nilsson, J. Phys. Chem. A, 2001, 105, 9954–9960 CrossRef CAS.
  15. U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen, J. Chem. Phys., 1995, 103, 8577–8593 CrossRef CAS.
  16. G. Bussi, D. Donadio and M. Parrinello, J. Chem. Phys., 2007, 126, 014101 CrossRef PubMed.
  17. J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing and J. Hutter, Comput. Phys. Commun., 2005, 167, 103–128 CrossRef CAS.
  18. J. VandeVondele and J. Hutter, J. Chem. Phys., 2007, 127, 114105 CrossRef PubMed.
  19. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868 CrossRef CAS PubMed.
  20. J. Schmidt, J. VandeVondele, I.-F. W. Kuo, D. Sebastiani, J. I. Siepmann, J. Hutter and C. J. Mundy, J. Phys. Chem. B, 2009, 113, 11959–11964 CrossRef CAS PubMed.
  21. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104 CrossRef PubMed.
  22. S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 1703–1710 CrossRef CAS PubMed.
  23. C. Hartwigsen, S. Goedecker and J. Hutter, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 58, 3641–3662 CrossRef CAS.
  24. L. Mones, A. Jones, A. W. Götz, T. Laino, R. C. Walker, B. Leimkuhler, G. Csányi and N. Bernstein, J. Comput. Chem., 2015, 36, 633–648 CrossRef CAS PubMed.
  25. C. Várnai, N. Bernstein, L. Mones and G. Csányi, J. Phys. Chem. B, 2013, 117, 12202–12211 CrossRef PubMed.
  26. N. Bernstein, C. Várnai, I. Solt, S. A. Winfield, M. C. Payne, I. Simon, M. Fuxreiter and G. Csányi, Phys. Chem. Chem. Phys., 2012, 14, 646–656 RSC.
  27. A. Peguiron, L. C. Ciacchi, A. D. Vita, J. R. Kermode and G. Moras, J. Chem. Phys., 2015, 142, 064116 CrossRef PubMed.
  28. T. Laino, F. Mohamed, A. Laio and M. Parrinello, J. Chem. Theory Comput., 2005, 1, 1176–1184 CrossRef CAS PubMed.
  29. T. Laino, F. Mohamed, A. Laio and M. Parrinello, J. Chem. Theory Comput., 2006, 2, 1370–1378 CrossRef CAS PubMed.
  30. A. Jones and B. Leimkuhler, J. Chem. Phys., 2011, 135, 084125 CrossRef PubMed.
  31. D. T. Limmer, A. P. Willard, P. Madden and D. Chandler, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 4200–4205 CrossRef CAS.
  32. T. Giorgino, Comput. Phys. Commun., 2014, 185, 317–322 CrossRef CAS.
  33. S. Sakong and A. Groß, J. Chem. Phys., 2018, 149, 084705 CrossRef PubMed.
  34. C. G. Venkatesh, S. A. Rice and A. H. Narten, Science, 1974, 186, 927–928 CrossRef CAS PubMed.
  35. A. K. Soper and M. A. Ricci, Phys. Rev. Lett., 2000, 84, 2881–2884 CrossRef CAS PubMed.
  36. S. Kaya, D. Schlesinger, S. Yamamoto, J. T. Newberg, H. Bluhm, H. Ogasawara, T. Kendelewicz, G. E. Brown Jr., L. G. M. Pettersson and A. Nilsson, Sci. Rep., 2013, 3, 1074 CrossRef PubMed.
  37. J. L. C. Fajín, M. N. D. S. Cordeiro and J. R. B. Gomes, J. Phys. Chem. A, 2014, 118, 5832–5840 CrossRef PubMed.
  38. S. N. Steinmann, R. Ferreira De Morais, A. W. Götz, P. Fleurat-Lessard, M. Iannuzzi, P. Sautet and C. Michel, J. Chem. Theory Comput., 2018, 14, 3238–3251 CrossRef CAS PubMed.


Electronic supplementary information (ESI) available: Details of the (a) validation of our modified SPC/E model for bulk water, (b) computation of layer resolved radial distribution function and (c) comparison of results obtained from QMMM-MD and classical MD simulations. See DOI: 10.1039/c9cp03558c
Current address: CareerTrails Academy LLP.
§ Current address: Prescience Insilico Private Limited, Bangalore, India.

This journal is © the Owner Societies 2019