The Role of Quantum Effects on Structural and Electronic Fluctuations

Article pubs.acs.org/JPCB

The Role of Quantum Effects on Structural and Electronic Fluctuations in Neat and Charged Water Federico Giberti,† Ali A. Hassanali,*,‡ Michele Ceriotti,*,¶ and Michele Parrinello† †

Department of Chemistry and Applied Biosciences, ETH Zurich and Università della Svizzera Italiana, via G. Buffi 13, CH-6900 Lugano, Switzerland ‡ The Abdus Salaam International Center for Theoretical Physics, Condensed Matter Physics Section, Strada Costiera 11, I-34151 Trieste, Italy ¶ Laboratory of Computational Science and Modeling, IMX, École Polytechnique Fédérale de Lausanne, 1015 Lausanne, Switzerland S Supporting Information *

ABSTRACT: In this work, we revisit the role of nuclear quantum effects on the structural and electronic properties of the excess proton in bulk liquid water using advanced molecular dynamics techniques. The hydronium ion is known to be a weak acceptor of a hydrogen bond which gives it some hydrophobic character. Quantum effects reduce the degree of this hydrophobicity which facilitates the fluctuations of the protons along the wires compared to the classical proton. Although the Eigen and Zundel species still appear to be dominant motifs, quantum fluctuations result in rather drastic events where both transient autoionization and delocalization over extended proton wires can simultaneously occur. These wild fluctuations also result in a significant change of the electronic properties of the system such as the broadening of the electronic density of states. An analysis of the Wannier functions indicate that quantum fluctuations of neat water molecules result in transient charging with subtle similarities and differences to that of the excess proton.

1. INTRODUCTION Understanding the mechanisms by which the proton moves through liquid water is a classical problem in physical chemistry.1−3 Unlike other ions which move exclusively via hydrodynamic diffusion, the proton is capable of migrating through the hydrogen bond (HB) network in a process commonly referred to as the Grotthuss mechanism.4 Most of the insights into the molecular mechanisms associated with this process have come from a combination of both ab initio molecular dynamics (AIMD) and empirical potentials based on the multistate empirical valence bond formalism (MSEVB). The commonly accepted picture from AIMD simulations describes the proton moving from one water molecule to the next in a stepwise manner with the Eigen and Zundel cations being dominant limiting states.5−10 In the Eigen cation, the proton is localized on one water molecule forming the H3O+ ion, which is stabilized by strong hydrogen-bond interactions with the surrounding water molecules. On the other hand in the Zundel species, the proton is shared equally by two water molecules. We have recently revisited the accepted view of the Grotthuss mechanism built on previous AIMD simulations.11,12 Medium-range directional correlations of water molecules within closed rings result in the presence of proton wires, which provide conduits for proton jumps over several water molecules. Rather than undergoing uninterrupted diffusion, the proton is subject to traps in the hydrogen bond network where it resides for longer times than expected giving it a multiscale and multidynamical character. MSEVB simulations of the © 2014 American Chemical Society

excess proton in water have also elucidated the collective behavior of proton transport and shown that the movement of the excess charge involves the reorganization of many water molecules surrounding it over a broad range of time scales.13−16 In the preceding theoretical studies, situations where the proton delocalizes over more than two water molecules occur rarely, for example, during concerted hopping events. It should be emphasized that the great majority of simulations of liquid water with charge defects using both AIMD and MSEVB treat the nuclei classically. It is well-known that both the structural and dynamical properties of light particles such as protons are significantly affected by nuclear quantum effects (NQEs).17−20 Earlier ab initio path integral (PI) molecular dynamics studies (PIMD) of the excess proton in water suggested that zero-point energy effects (ZPE) significantly lower the barrier for PT and sometimes led to situations where the proton was delocalized over more than two water molecules.5 The lowering of the barrier for PT was also observed by Voth and co-workers who performed simulations using an EVB model coupled with the centroid molecular dynamics (CMD) method.21−23 Besides these studies in bulk systems, there is also some indication that in confined environments NQEs enhance the delocalization of the excess proton beyond just one hydrogen bond.24 Although these studies have given some qualitative indications of longerrange proton delocalization, the statistical significance of these Received: August 1, 2014 Revised: September 30, 2014 Published: October 6, 2014 13226

dx.doi.org/10.1021/jp507752e | J. Phys. Chem. B 2014, 118, 13226−13235

The Journal of Physical Chemistry B

Article

fluctuations has been underappreciated in most descriptions of PT in liquid water.1,3,11,16 In this work, we use state of the art methodologies that combine smart thermostatting techniques and path integral (PIGLET) approaches to include NQEs25,26 to revisit the role of quantum effects on the structural and electronic properties of the excess proton in liquid water. Similar to previous studies of the classical proton, we find that the proton can reside in one region of the network without undergoing long-range diffusion making the convergence of our three independent simulations quite challenging.11 Despite this limitation, our simulations provide valuable information on how quantum fluctuations influence the delocalization of hydrogens around the excess proton and in the surrounding hydrogen bond network in these trapped states. While NQEs do not change the average value of the examined properties, they significantly alter their fluctuations enhancing both short and long-range delocalization of hydrogen atoms along the hydrogen bonds in the vicinity of the excess proton. An important consequence of this feature is that there is a non-negligible probability for the proton to delocalize over 3−4 water molecules confirming qualitative observations made in previous path integral simulations of the excess proton.5 This feature is facilitated by the fact that the enhanced quantum fluctuations of protons being donated to the hydronium oxygen reduce its hydrophobic character making it easier to accommodate the proton across several hydrogen bonds. Quantum effects also result in important changes in the electronic properties of the system. In particular, NQEs introduce new electronic energy levels as revealed by the density of states (DOS) akin to a form of thermal broadening. An analysis of the electronic configurations as probed by the Wannier functions also shows that intrinsic quantum effects in neat water resemble those of the excess proton. Molecules involved in extreme quantum fluctuations exhibit a strongly ionic character, that would not be captured by either classical AIMD simulations or empirical potentials of water. The paper is laid out as follows: we begin by describing the systems simulated and the computational methods used in Section 2. In Section 3 we catalogue and contrast the wide diversity of different species generated from classical and quantum simulations of neat water and the excess proton. In Section 4 we provide a semiquantitative analysis on how the excess proton changes the degree of delocalization of the hydrogen bonds in the network. In Section 5 we describe the role of quantum effects on the electronic density of states and the Wannier distributions. Finally, we end with some concluding remarks on our work in Section 6.

simulations has been the PI approach,29−31 which is unfortunately computationally prohibitive for the systems sizes that we are interested in. It has recently been shown, however, that this limitation can be circumvented by combining PI using a fewer number of beads with thermostats based on the generalized Langevin equation (PIGLET).26 PIGLET simulations of neutral liquid water indicate that the radial distribution functions, mean potential and quantum kinetic energy using 6 beads match those obtained from full 32 bead path integral calculations. PIGLET thus provides a tractable way of examining the role of NQEs in larger and more complex systems. To practically include them in our AIMD, we take advantage of the recently released code i-PI, that decouples the calculation of the interatomic forces from the dynamic evolution of the nuclei.32 We performed three independent AIMD PIGLET simulations of the excess proton in liquid water. The starting configurations for these simulations were obtained from classical simulations in our previous study.11 The three simulations were run for a total of 35, 15, and 13 ps. As we will see later, the qualitative conclusions made on the role of quantum fluctuations on the delocalization do not change in each of these three trajectories. A visual inspection of our trajectories indicates that similar to the classical trajectories, the short simulation times and small system sizes make it difficult to extract information on the long-range proton diffusion. In the first trajectory (QP1), the proton resides mostly in one region of the hydrogen bond (HB) network while making significant excursions along local HB wires. In the second trajectory (QP2), the proton shuffles along a proton wire made up of four molecules but remains otherwise trapped in one region of the network. In the third trajectory (QP3), we observe long-range migration of the proton away from a trap site. However, since the process occurs on a much faster time scale than that taken to reorganize the solvent environment around this trap, the periodic boundary conditions of the box allow the proton to return to the original trap site. Despite these limitations, our simulations provide an important starting point for quantifying the type of fluctuations that the proton experiences along the HB network and the differences from the classical picture. 2.2. Graph Representation of the Hydrogen Bond Network. Inspired by recent work that looked at the HB network from a graph theoretical perspective, we consider oxygen atoms as the nodes of a network, and hydrogen bonds as the edges of the graph.11 To connect two nodes, we use a simple definition of a putative hydrogen bond, as follows. We assume that every H atom contributes to one and only one HB, which is the one for which d(O,H) + d(O′,H) is minimal. By doing so, we ignore corner-case geometries such as bifurcated HBs or highly distorted configurations in which the H can hardly be seen as forming a bond. This is, however, irrelevant for the following analysis, where we focus on well-formed HBs, which are unambiguously identified by this definition. Once the hydrogen bond has been assigned to an oxygen pair, we calculate the proton transfer coordinate (ν = d(O,H) − d(O′,H)). It has recently been shown that ν provides a good measure of well-formed hydrogen bonds, and at the same time serves to characterize quantitatively the strength of the bond which can interpolate between a broken or weak interaction, and a strong, delocalized proton along the bond.28,33 The instantaneous proton potential along νand in particular the position of the most stable configuration of the hydrogen bond

2. COMPUTATIONAL METHODS 2.1. Ab Initio MD/PIGLET Simulations. The analysis presented in this work involves 4 systems that were simulated with AIMD using the CP2K package:27 simulations of classical water (CW), quantum water (QW), classical proton (CP) and quantum protonated water (QP) in a cubic box of side 15.6404 Å consisting of 128 water molecules. The CW, QW, and CP simulations were analyzed from trajectories obtained during recent studies.11,28 In these simulations the BLYP functional was employed along with the DZVP basis set and a plane wave cutoff of 300 Ryd. The reader is referred to those papers for more details on the calculations. It is well appreciated that light atoms such as protons are greatly affected by NQEs, which alter both the structural and dynamical properties of systems such as liquid water. The standard method for the inclusion of NQEs in 13227



Article

do not represent limiting states of a particular chemical species, but rather they provide a practical way to quantify and compare the extent of the fluctuations of the hydrogen bond and the degree of delocalization of the protons in the CW, QW, CP, and QP simulations. Once the procedure is completed and the oxygen atoms have been clustered, we have immediate access to the number of clusters at each frame, as well as the number of atoms within each cluster and their charge. It is easy then to identify that a cluster composed of three atoms consists of a single water molecule. Clusters composed of 3nO atoms, where nO is the number of oxygen atoms and nO ≥ 2, are neutral clusters created from transient autolysis events, while clusters made up of 3nO + 1 atoms, where nO ≥ 1, are positively charged clusters hosting the excess proton. For CP and QP, once we have identified a single positively charged cluster in each frame, we can introduce the center of charge x+ = ∑iqixi, where xi are the atom coordinates, qi is a formal charge of +1 for hydrogen atoms and −2 for oxygen atoms, the sum runs over the atoms assigned to the cluster, and we take into account periodic boundary conditions. We will use x+ as a reference point to disentangle the interplay between (quantum) fluctuations and the perturbation to the HB network introduced by the excess proton.

for a given environmentwould be the ideal tool to characterize the degree of delocalization. Since we do not have easy access to these quantities, we considered the instantaneous value of the proton transfer coordinate to be a proxy for the level of delocalization, and we set a cutoff value to qualify a proton as shared. In practice, two oxygens are connected in the graph if the HB possesses a value of |ν| < 0.2 Å. The choice of the value of this cutoff parameter will naturally have effects on the quantitative details of the cluster analysis. In particular, focusing on the CP simulations, a cutoff larger than 0.2 Å would lead to more Zundel species, while a smaller cutoff would enhance the proportion of Eigen-like species. In the context of the current work, we have checked that our main conclusions, namely, the effect that quantum fluctuations has on causing wild delocalization events of both neutral water molecules and the excess proton, is not sensitive to the choice of this cutoff (see the Supporting Information (SI) for results of using |ν| < 0.15 Å and |ν| < 0.25 Å). 2.3. Clustering Algorithm. Quantum fluctuations of both the excess charge as well as neutral water molecules in the hydrogen bond network result in a wide variety of different structures due to the varying tendency for the protons to delocalize over several hydrogen bonds. Ceriotti et al.28 have in fact shown that in quantum simulations of neat water, large fluctuations may lead to the formation of transient proton-like species. As we will see in detail later, these types of species also spontaneously form around the excess proton. Thus, in the QP simulations, using criteria such as the local coordination number or a Wannier analysis, it is possible to identify in any simulation snapshot several candidate sites in the water network that appear both structurally and electronically as proton-like. This makes it quite challenging to distinguish between species produced by “transient autolysis” events and the true excess proton. In order to disentangle the different proton-like species generated in the CW, CP, QW, and QP systems, we performed a cluster analysis, which we will now describe. Starting from the graph constructed in the previous section, we identify the clusters applying a Depth First Search algorithm.34 Once the oxygens have been clustered, we assign the hydrogen atoms to the nearest oxygen atom following the strategy used in ref 35. We can then calculate the formal charge of each cluster assigning a nominal −2 charge to the oxygen and +1 charge to the hydrogens. By applying this criterion, we found that in the quantum simulations there is a small probability of finding more than one cluster with a net +1 charge. These false positives come from transient autolysis events in which two or more water molecules undergoing a large fluctuation create a hydronium-hydroxide like pair at close distance. These pairs can either be an Eigen-like or a Zundel-like proton with hydroxide depending on the number of water molecules that are participating in the autolysis event. We eliminate this feature by performing an iterative algorithm, as follows. The presence of two clusters in the system with net +1 charge implies that there must be one other cluster with −1 charge. For this negatively charged cluster, we increase by 10% the cutoff on the value of ν, only for the hydrogen bonds that involve one of the O atoms that resides in the cluster. We repeat the clustering and iterate until all the +1/−1 pairs are connected. This ensures that no positive or negative clusters will be found in neat water, and only one positive cluster is obtained in the simulations with the excess proton. It is important to stress that the clusters identified by this procedure

3. CLUSTER ANALYSIS Quantum fluctuations of the hydrogen bonds tend to engender large deviations of water from its molecular character. In the presence of an excess charge, these effects include both intrinsic fluctuations of the proton and that of the surrounding water network. Much to our surprise, the structural and electronic properties of the species generated from autolysis events have close resemblance to those associated with the excess proton, making it rather challenging to catalogue the plethora of species generated. In this section, we focus on understanding the role of NQEs on the delocalization of the proton in the HB network and how these effects are different for CW, CP, QW, and QP. We begin our analysis by examining the differences in cluster distributions between the CW and QW simulations. Although quantum fluctuations lead to wild excursions of the protons, configurations with ν ≈ 0 are quite rare (about 1 in 1000).28 In the following analysis we focus on quantifying the cluster properties associated with these rare configurations. The total number of clusters in the CW and QW simulations is reported in Figure 1a, in black and red, respectively. These rare configurations involve pairs of water molecules where a proton is significantly shared between them. We see that in CW, most snapshots do not contain any clusters, since the vast majority of water molecules maintain their chemical integrity. On the other hand, for the QW simulations, the distribution is peaked at around 2 and is much broader, reflecting a significant strengthening of the hydrogen bonds, which results in an increased probability of observing water molecules that deviate from their molecular character. The same distribution is reported in Figure 1b for the CP (black) and QP (red) simulations. It appears as though at least based on this metric, the presence of the excess charge has little effect on the distribution of the number of clusters. Deeper insight into the structural origins of these clusters can be obtained by examining the number of water molecules that is involved in each of these neutral clusters. These distributions are compared for CW and QW in Figure 1c. The few neutral clusters found in neat water involve two water molecules 13228



Article

Figure 2. Histogram of the number of constituents for the cluster containing the excess proton for the three quantum simulations and one classical simulation. The classical trajectory is reported in black. The three QP trajectories are reported in red, blue, and purple, respectively (QP1, QP2, and QP3).

the HB network. These results are somewhat anticipated from earlier path integral simulations (see Figure 3 of the earlier work by Marx et al.5) but their statistical significance was not quantified. As indicated earlier, our three PIGLET simulations led to different types of proton traps in the HB network. Pinpointing

Figure 1. (a) Total number of clusters in the CW (black) and QW (red) simulations. (b) Total number of clusters in the excess proton simulation with (red) and without (black) NQEs. (c) Number of constituents of the neutral clusters for CW (black) and QW (red). (d) number of constituents for the neutral cluster for the three NQEs (red, blue, and purple) and a classical simulation (black). Although the most probable cluster is composed of two water molecules undergoing an autolysis event, in the presence of NQEs there is a noticeable probability of finding even overcoordinate water and autolysis chains spanning three or four water molecules.

undergoing transient autoionization. On the other hand, in QW there is a small but non-negligible probability of finding clusters composed of three water molecules. These extended clusters results from cooperative autolysis events where the incipient proton and hydroxide form over a longer hydrogen bond wire. This is rather interesting when considered within the context of the autoionization of liquid water. By examining the timereversed process, namely, the recombination of the hydronium and hydroxide ions, we have recently found with classical AIMD simulations, that this process involves concerted proton transfer across a water wire made up of three hydrogen bonds.36 Analogously, quantum fluctuations along HB wires in neutral water are so drastic that they instantaneously create quasi-ionized species without the need for any rare electric field fluctuations. Revisiting the role of quantum effects on the mechanism of recombination would be an interesting topic for future studies. Similar trends are also observed in CP and QP which is reported in Figure 1d. Perhaps more interesting is the distribution of the number of species within the clusters involving the excess proton shown in Figure 2. In this case, the charged cluster consists of at least one water molecule and an additional H+. CP leads to peaks with one and two water molecules with an excess charge which effectively correspond to the Eigen and Zundel-like species, respectively. The QP simulations on the other hand are characterized by more complex behavior. We observe a much broader distribution in the number of water molecules in the cluster. There is, however, a significant number of 3-water molecule clusters and a small but non-negligible fraction of 4water and 5-water clusters. These correspond to situations where protons along several hydrogen bonds delocalize in a concerted manner. These species are completely absent in the classical simulations an reinforce the importance of the role of quantum effects on the delocalization of the excess proton in

Figure 3. Probability of fluctuations of hydrogen bonds in a simulation of neat water with classical nuclei. The lower panel shows the overall distribution of proton transfer coordinate ν for hydrogen bonds; negative values correspond to donated bonds, and positive values to accepted bonds. The upper panel shows the joint probability distribution. The joint probability distribution P(νA,νB) is symmetric by construction along the main diagonal, and so we only show the upper left corner. In the lower right corner, instead, we show the ratio between the joint probability and the product of the marginals, which highlights correlations between the hydrogen bonds that involve the tagged oxygen atom. 13229



Article

at ν = ± 0.75 Å, and there is negligible probability of observing a shared proton with ν = 0 (lower panel in Figure 3). The chemical nature of water molecules and the correlations between fluctuations across the HB network become more apparent when one considers the joint probability distribution P(νA,νB) of two hydrogen bonds involving the same oxygen atom O. In this construction, correlations are determined across wires made up of two hydrogen bonds, and longer range correlations are neglected. The peaks along the diagonal are about half as strong as the off-diagonal one. This is because in an ideal tetrahedral environment, each water donates two HBs and accepts two HBs, so if one hydrogen bond is donated (ν ≈ −0.75 Å), there is twice the chance that a second HB chosen at random around the same O atom is being accepted (ν ≈ 0.75 Å) than that it is the second donated HB. These kinds of correlations can be visualized more clearly by plotting the ratio between the joint probability distribution P(νA,νB) and the product of the marginals P(νA) and P(νB). This ratio quantifies how more or less likely it is to find a given configuration relative to the case in which there was no correlation between the two HBs. For instance, a ratio of one implies no correlation, a ratio of 10 would mean that the joint configuration is 10 times more likely than in the absence of correlations and a ratio of 0.1 implies that the configuration is 10 times less likely than if the HBs were uncorrelated. This quantity is shown in the lower half of the main panel in Figure 3, and demonstrates clearly that donor/acceptor or acceptor/donor configurations are more likely than donor/donor or acceptor/acceptor pairs. As discussed before, this is just a manifestation of the molecular nature of water molecules, and is reminiscent of the ice rules in crystalline phases of water. Besides these trivial correlations that occur simply because of the local tetrahedrality of water, it is also clear that within each quadrant of the figure, fluctuations around the most likely configuration are highly correlated. If an accepted hydrogen bond has a large fluctuation toward the tagged O atom (νA ≳ 0), then there is an increased probability that a donated HB will experience a large fluctuation toward zero (νB ≳ 0). On the contrary, extreme fluctuations of HBs in the same direction are negatively correlated: the probability of two donated HBs to fluctuate simultaneously toward ν = 0 is less likely than if there was no cross-correlation, and the same is the case when one considers two accepted HBs. These differences in the correlations reflect the extent of the nonlocality of the asymmetric and symmetric vibrational stretches in water. The asymmetric stretch is more susceptible to a nonlocal delocalization over at least 2 hydrogen bonds. It is interesting to note that this correlated behavior of the asymmetric stretch is reminiscent of the early events observed in concerted proton hopping of the excess proton.11 Quantum effects dramatically increase the extent of fluctuations of protons along the hydrogen bonds. Figure 4 clearly shows the large extent of fluctuations in both the distribution of a single bond P(ν) and in the joint probability distribution. Note that the absolute fraction of extreme fluctuations to ν ≈ 0 is probably overestimated by at least a factor of 2, because of the excessive propensity of generalized gradient approximation (GGA) water to undergo charge transfer. Comparison with different choices of density functionals, including hybrid functionals, shows, however, that the relative increase over a classical simulation and the qualitative features of this observation are remarkably stable.28 The same qualitative cross-correlations seen in classical simulations are

the microscopic origins of these differences is beyond the scope of the current study. However, we can make some qualitative remarks about how the differences between QP1, QP2, and QP3 affect the results in Figure 2. For example, the more localized trap (QP1, in red in the figure), possesses a higher population of Eigen structures and fewer long chains compared to QP2 and QP3. This suggests that the excess proton feels a slightly different thermodynamic potential in these three simulations, which can influence its degree of delocalization in a nontrivial way. To conclude, the preceding analysis shows that while the Eigen and Zundel-like species still appear to be the dominant structures, quantum fluctuations of the excess charge in bulk water are quite drastic and thus longer range delocalization of the proton does not appear to be an exclusive feature of protons in confined environments.24

4. FLUCTUATIONS AND CORRELATIONS IN HYDROGEN BONDS In the previous section, we have cataloged the rich variety of species resulting from quantum fluctuations of the hydrogen bonds both in water and in the presence of the excess proton. These enhanced fluctuations ultimately result from the softening of the potential of mean force for protons to shuffle along hydrogen bonds in the water network. The presence of clusters with 2, 3, and 4 water molecules in QP suggests that this delocalizing effect is quite long-range. In the ensuing discussion, we try to give a semiquantitative analysis of how quantum effects and the presence of an excess charge conspire to determine the behavior of both bound and free protons in water. The definition we have used for a hydrogen bond does not rely on knowing the chemical identity of water molecules, and can therefore be safely used also in the case of protonated water, where molecules in the vicinity of the excess charge undergo bond-changing fluctuations. We base our analysis on the distribution of the proton transfer coordinate, and we consider its value relative to a given oxygen atomso that negative values of ν correspond to donated HBs, and positive values to accepted HBs. For each oxygen atom, one can define the probability distribution P(ν) for all the HBs it is involved in (both donated and accepted), as well as the joint probability distribution P(νA,νB) for two hydrogen bonds involving the same oxygen atom O. In order to give a semiquantitative account of the effect of the presence of an excess charge, we will also comment on how these probability distributions change with the distance of an oxygen from the center of charge x+ of the cluster associated with the excess proton. 4.1. Classical and Quantum Neat Water. The QP system is characterized by both intrinsic quantum fluctuations of the excess proton, where it can delocalize over several water molecules as well as those of neutral waters undergoing transient autolysis events. To develop a deeper understanding of the complexity of this problem, it will be useful to calibrate ourselves on the nature of the hydrogen bond fluctuations and correlations that exist in CW and QW. We begin by comparing the single bond P(ν) and joint probability distributions. In neat water, sitting on an oxygen atom and looking at the hydrogen bonds that it participates in, one will find on average about four HBs: two in which the tagged oxygen is acting as a HB donor (ν < 0), and two in which it acts as a HB acceptor (ν > 0). The overall distribution is symmetric across ν = 0 by construction. In the case of a simulation with classical nuclei, P(ν) is peaked 13230



Article

hydrogen bonds. Extracting similar mechanistic details from the PIGLET simulations is not feasible, since it does not give access to dynamical properties. Hence, in the ensuing discussion we focus on the differences in structural perturbations that a classical or quantum proton would make to the hydrogen bond network as one moved further away or closer to the center of charge. Figure 5 shows that at large d the distribution of the PT

Figure 4. Probability of fluctuations of hydrogen bonds in a simulation of neat water with quantum nuclei. See Figure 3 for a detailed description of the plot.

also observed here, with simultaneous fluctuations of the same kind of bond being depleted, and simultaneous fluctuations of a donated and an accepted bond being promoted. Note that cross-correlations are being somewhat underestimated here. Since we do not assume a priori the directionality of hydrogen bonds, our analysis cannot distinguish between a donated HB that fluctuates toward ν = 0 and an accepted HB experiencing an extreme fluctuation beyond zero. Performing the same analysis assuming that the chemical identity of water molecules is preserved along the simulation shows that cross-correlations between extreme fluctuations are very strong.28 Furthermore, since they happen with considerably higher probability, they are much more relevant to the equilibrium description of liquid water than one would expect based on an analysis that neglects nuclear quantum effects. 4.2. Classical and Quantum Proton. Let us move on to consider the effect of an excess proton on the fluctuations of the HB network. We aim at verifying how P(ν) and the crosscorrelations for HBs that involve the same oxygen atom vary as a function of the distance between the O and the excess charge. We define the conditional proton-transfer distribution as P(ν,d)/P(d), where P(ν,d) is the joint probability of finding a hydrogen bond with a PT coordinate ν and involving an oxygen atom at a distance d from the center of charge, and P(d) is the distribution of distances between O atoms and the center of charge x+. As indicated earlier, we find that within the limited time scales of our current simulations, there is a strong tendency for the excess proton to be trapped in one region of the network while undergoing large fluctuations. The trapping behavior was also observed in classical simulations, and it was followed by bursts of activity involving correlated proton hops over several

Figure 5. (Upper panel) Conditional probability distribution of the proton transfer coordinate as a function of the distance between the oxygen atom and the center of excess charge. (Lower panel) Four slices of the conditional probability at the distances indicated, on a linear scale.

coordinate is symmetric about ν = 0 as in the case of neat water. At shorter distances, the presence of an excess charge breaks the symmetry, and one can recognize the signature of the different limiting structures of the hydrated proton. Configurations with very small d correspond to Eigen-like configurations, where the center of excess charge coincides with the O atom at the center of the complex. One sees a very strong asymmetry between ν > 0 and ν < 0, with a predominance of donated HBs and few, elongated accepted bonds. This feature captures the hydrophobic nature of the Eigen cation, and the strengthening of the donated HB which has a maximum probability for ν ≈ 0.5 Å. Oxygen atoms that are at 1−2 Å from the excess charge are likely to be part of a Zundel-like structure. At this distance, there is a large probability of finding a shared proton with ν = 0, donated HBs are strengthened (albeit to a lesser extent than at very short d), and there is a depletion of density of accepted HBs. At larger distances the distribution of ν becomes symmetric, but the effect of the excess charge is still visible in an increased delocalization of protons along the hydrogen bond up to d ≈ 3.5 Å. While we do not have enough statistics or large enough boxes to comment on the asymptotic behavior, it appears that by d ≈ 4 Å the probability of observing a wild fluctuation to ν = 0 has become very low. 13231



Article

≈ 0.4 Å and, in turn, relax the requirement of the “presolvation” criterion. In other words, it is possible for the proton to delocalize over several water molecules, without requiring or involving a drastic change in the coordination numbers of water molecules. Beyond d ≈ 4 Å, the potential of mean force along the HBs is still somewhat stronger than in neat water, but has become essentially symmetric. In summary, if we compare Figures 3, 4, 5, and 6, it is apparent that the effects of quantum nuclei and of the presence of an excess proton in the hydrogen bond network are, to a first approximation, additive. Both conspire to soften the potential of mean force for the motion of a H atom in between the HB partner oxygen atoms. This means that the effect of the excess proton is noticeable up to well beyond 3.5 Å, even though the strong “chemical” correlations that are seen at short distances from x+ are very similar to those observed in a classical simulation. One may wonder whether this long-range increase in delocalization can be explained in terms of polarization of the HBs due to the excess charge associated with the extra proton. To probe deeper into whether the behavior of the proton is a purely electrostatic effect at play, we performed PIGLET simulations of a solvated Na ion. Figure 7 shows the potential

Joint probability distributions also vary as a function of distance from the center of excess charge. Given the heterogeneity between different trajectories and the difficulties in converging a three-dimensional histogram, we do not believe that the present calculations make it possible to comment quantitatively on these features, which are reported in the SI. At large distances from x+, a distribution very similar to the bulk is observed. As one moves closer to the excess charge, the asymmetry between ν > 0 and ν < 0 becomes more apparent, and correlations are enhanced. In the regions that we have associated with Zundel-like and Eigen-like clusters, the distribution becomes qualitatively different, and exhibits very large correlations between hydrogen bonds involving the same O atom. These strong correlations are a consequence of the fluxional character of the excess charge. Fluctuations of the HB are not simply (broad) oscillations around a stable structure, but involve the correlated motion of protons across the HB network both from autolysis events as seen in the CW and QW analysis and the intrinsic fluctuations of the excess charge. Figure 6 elucidates the role of quantum effects on the longrange strengthening of hydrogen bonds around the center of

Figure 7. (Upper panel) Conditional probability distribution of the proton transfer coordinate as a function of the distance between the oxygen atom and the Na+ ion in a simulation of a solvated sodium cation including NQEs. (Lower panel) Two slices of the conditional probability at the distances indicated, on a linear scale.

Figure 6. (Upper panel) Conditional probability distribution of the proton transfer coordinate as a function of the distance between the oxygen atom and the center of excess charge, for a simulation including NQEs. (Lower panel) Four slices of the conditional probability at the distances indicated, on a linear scale.

of mean force computed as a function of the distance between the reference oxygen and Na+. Other than the fact that oxygen atoms in the first solvation shell preferentially donate hydrogen bonds (the lone pairs tend to point toward the ion) there is barely any effect on the potential of mean force. This result suggests that the effect of the excess proton involves effects that are more complicated than that given by classical electrostatics. In particular, it is very likely, and we will see evidence of this later through an examination of the electronic structure, that these effects involve a mixture of electrostatics, polarization, and charge transfer, which simultaneously act together to give the proton its peculiar behavior. The long-range HB strengthening and the extended delocalized states that are

excess charge. Qualitatively, we observe similar features as in the classical proton. However, the combination of intrinsic fluctuations of the excess proton generating clusters with three to four water molecules hosting the proton, as well as quantum fluctuations of neutral waters, results at all distances in significantly more delocalized hydrogen bonds compared the classical proton in Figure 5. In particular, quantum effects reduce the hydrophobicity of both Eigen and Zundel like configurations, as probed by the extent of the donated HB at ν 13232



Article

to lowest unoccupied molecular orbital (HOMO−LUMO) gap of liquid water, changes in the band gap with the same method are more meaningful. In particular, the overlay of the QP and CP DOS at the Fermi level ∼0 suggests that NQEs change the band gap of water. This effect originates from changes in the population of both the occupied and nonoccupied levels. These results suggest that there should be an isotope effect on the optical properties of water whose origin would be related to NQEs. In fact, X-ray Raman spectroscopy experiments on liquid water have found differences in the spectra of hydrogenated and deuterated water, which was interpreted as originating from differences in asymmetries of the hydrogenbonding configurations between the two systems due to NQEs.38 Further studies by Fuchs and co-workers also found strong isotope effects in both the X-ray absorption and emission spectra of liquid water.39 Since our calculations are all performed on the ground state, a one-to-one comparison with the experiments is not possible at this point. However, it is clear that NQEs will have a sizable effect on both the absorption and emission spectra of liquid water and are likely to have an impact in excited state properties of systems characterized by hydrogen bonded networks. 5.2. Wannier Center and Spread Distributions. Earlier in the manuscript, we examined the clusters that existed in the CW, QW, CP, and QP simulations and showed that quantum fluctuations of the proton resulted in both short and long-range delocalization of the protons along hydrogen bonds. We will now examine the local electronic character of these structures by using the distribution of the Wannier functions associated with them. The Wannier functions provide a physically intuitive representation of the Kohn−Sham orbitals from density functional theory. For example, a single water has four Wannier centers associated with it: two involved in covalent bonds with the hydrogen atoms, and two forming lone pairs.40,41 Using the results from the clustering algorithm, we separately computed the distribution of the Wannier center and spreads (WC and WS) of clusters involved in transient autolysis events and those involved in the delocalization of the true excess charge. In the ensuing results, we focus on the QP simulations. The top panel of Figure 9a shows the density plot of WS− WC averaged over all the water molecules in the system. The distribution is characterized by two peaks in WC−WS space: one at ∼ (0.48,0.7) for the bonding centers and another at ∼ (0.3,0.75) for the lone-pairs. These results are consistent with those found by Silvestrelli et al. for classical neat water,40,41 even though quantum fluctuations result in considerably broader distributions. The second panel of Figure 9b shows the WC−WS density plot obtained by exclusively examining the clusters consisting of two water molecules exhibiting autolysis character, where the proton transfer coordinate along the hydrogen bond is ∼0. Ceriotti et al. recently examined the Wannier distributions associated with the transient autolysis events by focusing solely on the WCtheir analysis shows extensive polarization of the WC caused by quantum effects. Here, we disentangle the contributions from both the excess charge and those from quantum fluctuations in neat water in both WC and WS space. Several interesting features are found in this distribution. We observe four distinct regions in WC− WS space. It is rather striking that the electronic structure as probed by the WC−WS space of the autolysis clusters exhibit some important differences compared to the most of the other water molecules in the system. First, we observe that for the bonding centers an additional shoulder, where the WC is

seen in our simulations can be understood in terms of multiple defects shuffling around the hydrogen bond network. This involves a complex coupling between the fluxonial nature of the solvated proton and the quantum fluctuations of all the hydrogen bonds in the network.

5. NQES ON ELECTRONIC PROPERTIES The aforementioned results show that quantum effects result in a large perturbation to the fluctuations of the hydrogen bonds in the network. This is a manifestation of the long-range correlations of the HB’s away from the center of excess charge, as well as the formation of a rich diversity of water species characterized by nonmolecular or rather significant ionic character. We will now discuss the electronic perturbations resulting from these effects. 5.1. Electronic Density of States. We begin by examining the role of NQEs on the electronic DOS. An approximate DOS is easily obtained from solving the Kohn−Sham equations and obtaining the distribution of eigenvalues εi from which a trajectory averaged DOS is defined by f (ε) = ⟨∑ δ(ε − εi)⟩ i

(1)

Comparison of this quantity evaluated for the QP and CP simulations probes to what extent NQEs, which essentially perturb the configurations of the ionic nuclei, affect the electronic Hamiltonian. Figure 8 below compares the DOS

Figure 8. Electronic DOS for CP (black) and QP (red) simulations.

obtained for the QP and CP systems. The CP DOS is similar to previously reported CPMD simulations of neat and protonated liquid water.37 In particular for the valence band, one obtains four peaks corresponding to the 2a1, 1b2, 3a1, and 1b1 bands from more negative to positive values on the x-axis. While the peak positions for the CP and QP systems are quite similar, in the latter each of the 4 bands is broader. There is thus a significant population of electronic states that are completely forbidden in CP, but when NQEs are included the quantum fluctuations adiabatically introduce new types of occupied levels. The changes in all four bands of these occupied levels reflect changes in both the lone-pair and bonding electrons and thus have direct effects on the properties of the hydrogen bonds described earlier. The widening of each of these peaks is very similar to effects of thermal broadening on the line-spectra. In this case, the results of Figure 8 show that broadening is introduced by quantum rather than thermal fluctuations. While it is well-known that GGA density functional theory (DFT) underestimates the highest occupied molecular orbital 13233



Article

present and earlier results for the classical proton strongly suggest the presence of a slower collective reorganization of the network, which is required to release the proton from the trapped states and is not being adequately sampled within the time scales of the ab initio molecular dynamics. Quantum effects do not appear to change the average properties of the system, but instead result in significantly enhanced fluctuations. These fluctuations in neat water result in rare but non-negligible events where transiently formed hydronium-hydroxide-like pairs appear. At the same time, quantum fluctuations of the excess proton result in wild delocalization events where the charge can spread over 2−5 water molecules. The hydronium ion is a weak acceptor of a hydrogen bond and is hence thought to be hydrophobic. In this context, the transfer or delocalization of the proton over several water molecules has been thought to be a high barrier process since it would need the weakening of several hydrogen bonds being donated to regions where the proton delocalizes. We instead find here that quantum effects reduce the hydrophobicity of the excess proton and at the same time weaken the “pre-solvation” requirement. This phenomenon is strongly coupled to the collective delocalization of the protons along hydrogen bonds in the neighborhood of the excess proton, and results in considerably longer-range delocalization than it is observed in classical simulations. Besides the perturbations made to structural properties, quantum effects also alter the electronic properties of the system. In particular, NQEs appear to have a similar effect on the electronic density of states as thermal broadening does, resulting in the introduction of new previously inaccessible electronic levels and a change in the band gap of water. Thus, quantum effects have quite a drastic impact on both the bonding and antibonding properties in a molecular orbital picture and must be taken into account to develop a deeper knowledge of the coupling of nuclear motions, electronic degrees of freedom and hence chemistry. Understanding the role of all these effects in tuning the behavior of water at hydrophobic and biological interfaces will be an interesting area to explore in the near future.

Figure 9. Wannier center (WC) and Wannier spread (WS) distributions for QP simulations. Top panel shows distribution for all waters in the system, while the lower panels considers only the Wannier functions of O atoms that belong to clusters that have been identified as transient autolysis, Zundel-like structures and Eigen-like structures.

displaced inward and the WS is slightly more inflated, is formed. This feature is caused by the bonding center along the O−H bond that is elongating. On the other hand, for the lonepair we observe a shoulder where the WC is displaced outward and the WS is slightly reduced. This in turn is a perturbation induced by the water molecule that is receiving the transiently dissociating proton from the elongated O−H bond. The bottom two panels of Figure 9 show similar density plots but now for Eigen and Zundel clusters consisting of one and two water molecules plus an additional proton. The electronic structure of the true excess charge is quite different from most water molecules particularly for the Eigen species. As one might expect the presence of the three bound protons around the oxygen hosting the proton creates a stronger signal from bonding WC. What is quite interesting is that the Zundelexcess charge, which most closely resembles the structure from an autolysis event, also has some density with some overlap in both WC and WS with the shoulder regions described earlier for the autolysis clusters. Thus, what we see here is that NQEs in neutral water introduce a peculiar type of similarity between the electronic structure of the true excess charge and that of neutral waters.

■

ASSOCIATED CONTENT

S Supporting Information *

Additional information on the clustering sensibility and the variation of the joint probability distribution as a function of distance from the excess charge are reported in the Supporting Information. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: michele.ceriotti@epfl.ch.

6. CONCLUSION In this work, we have used recent developments in colored noise thermostats coupled with path integral molecular dynamics to revisit the role of nuclear quantum effects on the structural and electronic properties of the excess proton in liquid water. Despite the dramatic computational savings made possible by using a GLE to reduce the number of path integral replicas, statistical convergence of these calculations has proven to be quite a challenge. In all of our PIGLET simulations we find that the excess proton resides in one part of the network without undergoing long-range proton diffusion. Both our

Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS We would like to acknowledge generous allocation of CPU time by CSCS under the project id s466. REFERENCES

(1) Atkins, P.; de Paula, J. Physical Chemistry; Oxford University Press: Oxford, U.K., 2006.

13234



Article

(2) Eigen, M.; De Maeyer, L. Self-Dissociation and Protonic Charge Transport in Water and Ice. Proc. R. Soc. A 1958, 247, 505−533. (3) Agmon, N. The Grotthuss Mechanism. Chem. Phys. Lett. 1995, 244, 456−462. (4) de Grotthuss, C. J. T. Theory of Decomposition of Liquids by Electrical Currents. Ann. Chim. (Paris) 1806, LVIII, 54−74. (5) Marx, D.; Tuckerman, M. E.; Hutter, J.; Parrinello, M. The Nature of the Hydrated Excess Proton in Water. Nature 1999, 397, 601−604. (6) Marx, D.; Chandra, A.; Tuckerman, M. E. Aqueous Basic Solutions: Hydroxide Solvation, Structural Diffusion, and Comparison to the Hydrated Proton. Chem. Rev. 2010, 110, 2174−2216. (7) Tuckerman, M.; Laasonen, K.; Sprik, M.; Parrinello, M. Ab Initio Molecular Dynamics Simulation of the Solvation and Transport of H3O+ and OH-Ions in Water. J. Phys. Chem. 1995, 99, 5749−5752. (8) Tuckerman, M.; Laasonen, K.; Sprik, M.; Parrinello, M. Ab Initio Molecular Dynamics Simulation of the Solvation and Transport of Hydronium and Hydroxyl Ions in Water. J. Phys. Chem. 1995, 103, 150−161. (9) Hynes, J. T. Physical Chemistry: The Protean Proton in Water. Nat. 1999, 397, 565−567. (10) Berkelbach, T. C.; Lee, H.-S.; Tuckerman, M. E. Concerted Hydrogen-Bond Dynamics in the Transport Mechanism of the Hydrated Proton: A First-Principles Molecular Dynamics Study. Phys. Rev. Lett. 2009, 103, 238302. (11) Hassanali, A.; Giberti, F.; Cuny, J.; Kühne, T. D.; Parrinello, M. Proton Transfer Through the Water Gossamer. Proc. Natl. Acad. Sci. U.S.A. 2013, 110, 13723−13728. (12) Hassanali, A. A.; Giberti, F.; Sosso, G. C.; Parrinello, M. The Role of the Umbrella Inversion Mode in Proton Diffusion. Chem. Phys. Lett. 2014, 599, 133−138. (13) Lapid, H.; Agmon, N.; Petersen, M. K.; Voth, G. A. A BondOrder Analysis of the Mechanism for Hydrated Proton Mobility in Liquid Water. J. Chem. Phys. 2005, 122, 014506. (14) Wu, Y.; Chen, H.; Wang, F.; Paesani, F.; Voth, G. A. An improved Multistate Empirical Valence Bond Model for Aqueous Proton Solvation and Transport. J. Phys. Chem. B 2008, 112, 467−482. (15) Markovitch, O.; Chen, H.; Izvekov, S.; Paesani, F.; Voth, G. A.; Agmon, N. Special Pair Dance and Partner Selection: Elementary Steps in Proton Transport in Liquid Water. J. Phys. Chem. B 2008, 112, 9456−9466. (16) Knight, C.; Voth, G. A. The Curious Case of the Hydrated Proton. Acc. Chem. Res. 2011, 45, 101−109. (17) Morrone, J. A.; Car, R. Nuclear Quantum Effects in Water. Phys. Rev. Lett. 2008, 101, 017801. (18) Li, X.-Z.; Walker, B.; Michaelides, A. Quantum Nature of the Hydrogen Bond. Proc. Natl. Acad. Sci. U.S.A. 2011, 108, 6369−6373. (19) Giuliani, A.; Bruni, F.; Ricci, M.; Adams, M. Isotope Quantum Effects on the Water Proton Mean Kinetic Energy. Phys. Rev. Lett. 2011, 106, 255502. (20) Romanelli, G.; Ceriotti, M.; Manolopoulos, D. E.; Pantalei, C.; Senesi, R.; Andreani, C. Direct Measurement of Competing Quantum Effects on the Kinetic Energy of Heavy Water upon Melting. j. Phys. Chem. Lett. 2013, 4, 3251−3256. (21) Lobaugh, J.; Voth, G. A. The Quantum Dynamics of an Excess Proton in Water. J. Chem. Phys. 1996, 104, 2056−2069. (22) Pavese, M.; Chawla, S.; Lu, D.; Lobaugh, J.; Voth, G. A. Quantum Effects and the Excess Proton in Water. J. Chem. Phys. 1997, 107, 7428−7432. (23) Pavese, M.; Voth, G. A. Quantum and Classical Simulations of an Excess Proton in Water. Ber. Bunsenges. Phys. Chem. 1998, 102, 527−532. (24) Chen, J.; Li, X.-Z.; Zhang, Q.; Michaelides, A.; Wang, E. Nature of Proton Transport in a Water-Filled Carbon Nanotube and in Liquid Water. Phys. Chem. Chem. Phys. 2013, 15, 6344−6349. (25) Ceriotti, M.; Manolopoulos, D. E.; Parrinello, M. Accelerating the Convergence of Path Integral Dynamics with a Generalized Langevin Equation. J. Chem. Phys. 2011, 134, 084104.

(26) Ceriotti, M.; Manolopoulos, D. E. Efficient First-Principles Calculation of the Quantum Kinetic Energy and Momentum Distribution of Nuclei. Phys. Rev. Lett. 2012, 109, 100604. (27) VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. Quickstep: Fast and Accurate Density Functional Calculations Using a Mixed Gaussian and Plane Waves Approach. Comput. Phys. Commun. 2005, 167, 103−128. (28) Ceriotti, M.; Cuny, J.; Parrinello, M.; Manolopoulos, D. E. Nuclear Quantum Effects and Hydrogen Bond Fluctuations in Water. Proc. Natl. Acad. Sci. U.S.A. 2013, 110, 15591−15596. (29) Feynman, R. P.; Hibbs, A. R. Quantum Mechanics and Path Integrals; McGraw-Hill: New York, 1965. (30) Chandler, D.; Wolynes, P. G. Exploiting the Isomorphism Between Quantum Theory and Classical Statistical Mechanics of Polyatomic Fluids. J. Chem. Phys. 1981, 74, 4078−4095. (31) Parrinello, M.; Rahman, A. Study of an F Center in Molten KCl. J. Chem. Phys. 1984, 80, 860−867. (32) Ceriotti, M.; More, J.; Manolopoulos, D. E. i-PI: A Python Interface for Ab Initio Path Integral Molecular Dynamics Simulations. Comput. Phys. Commun. 2014, 185, 1019−1026. (33) McKenzie, R. H.; Bekker, C.; Athokpam, B.; Ramesh, S. G. Effect of Quantum Nuclear Motion on Hydrogen Bonding. J. Chem. Phys. 2014, 140, 174508. (34) Donald, E. K. The Art of Computer Programming. Sorting and searching 1999, 3, 426−458. (35) Schwegler, E.; Galli, G.; Gygi, F.; Hood, R. Q. Dissociation of Water under Pressure. Phys. Rev. Lett. 2001, 87, 265501. (36) Hassanali, A.; Prakash, M. K.; Eshet, H.; Parrinello, M. On the Recombination of Hydronium and Hydroxide Ions in Water. Proc. Natl. Acad. Sci. U.S.A. 2011, 108, 20410−20415. (37) Laasonen, K.; Sprik, M.; Parrinello, M.; Car, R. “Ab initio” Liquid Water. J. Chem. Phys. 1993, 99, 9080−9089. (38) Bergmann, U.; Nordlund, D.; Wernet, P.; Odelius, M.; Pettersson, L. G.; Nilsson, A. Isotope Effects in Liquid Water Probed by X-ray Raman Spectroscopy. Phys. Rev. B 2007, 76, 024202. (39) Fuchs, O.; Zharnikov, M.; Weinhardt, L.; Blum, M.; Weigand, M.; Zubavichus, Y.; Bär, M.; Maier, F.; Denlinger, J. D.; Heske.; et al. Isotope and Temperature Effects in Liquid Water Probed by X-ray Absorption and Resonant X-ray Emission Spectroscopy. Phys. Rev. Lett. 2008, 100, 027801. (40) Silvestrelli, P. L.; Parrinello, M. Water Molecule Dipole in the Gas and in the Liquid Phase. Phys. Rev. Lett. 1999, 82, 3308. (41) Silvestrelli, P. L.; Parrinello, M. Structural, Electronic, and Bonding Properties of Liquid Water from First Principles. J. Chem. Phys. 1999, 111, 3572−3580.

13235


The Role of Quantum Effects on Structural and Electronic Fluctuations

Recommend Documents