Quantum simulations of neutral water clusters and singly-charged water cluster anions

Abstract

We report a computational study of the structural and energetic properties of water clusters and singly-charged water cluster anions containing from 20 to 573 water molecules. We have used both a classical and a quantum description of the molecular degrees of freedom. Water intra and inter-molecular interactions have been modelled through the SPC/F model, while the water-excess electron interaction has been described via the well-known Turi–Borgis potential. We find that in general the quantum effects of the water degrees of freedom are small, but they do influence the cluster-size at which the excess electron stabilises inside the cluster, which occurs at smaller cluster sizes when quantum effects are taken into consideration.

---

I. Introduction

Sixty years after its first unequivocal spectroscopic detection,¹ the hydrated electron, e_aq⁻, continues to fascinate physicists and chemists alike. Since then a great deal of research effort has aimed at gaining a better understanding of its nature and properties, as testified by the large number of review articles devoted to this topic.^2–5 The reasons for this are several: firstly, e_aq⁻ is the simplest reducing agent in solution, and as such it can be expected to play an important role in chemical reactivity in aqueous media. It is a species known to intervene in electron transfer, electrochemistry, radiation effects in water, etc. Its presence has been detected not only in bulk water but also at water–gas interfaces, in thin water layers on metallic substrates, and in water clusters (water cluster anions). One could regard a solvated electron as one of the simplest, if not indeed the simplest, realization of a quantum system interacting with a (quasi)-classical bath, a model system that has been of interest to theoreticians since the early days of quantum mechanics. But perhaps the most obvious reason for the ongoing interest in the hydrated electron is the fact that, in spite of the great efforts, both experimental and theoretical, devoted to understanding this system, there are still several open questions concerning its microscopic nature and properties. At the experimental level, it has proved difficult to correlate spectroscopic signals with a unique and consistent microscopic picture of e_aq⁻. At the theoretical level, modelling of e_aq⁻ is challenging both because of the need to account for very different degrees of freedom (electron and water molecules), as well as accurately describing their mutual interaction. In spite of these challenges, a widely accepted paradigm of e_aq⁻ has emerged over the years, according to which the solvated electron resides in a roughly spherical cavity of excluded volume, being coordinated by approximately four water molecules in the nearest-neighbor shell (see e.g. ref. 5 and references therein). The cavity structural model has been challenged by a non-cavity one, put forward by Larsen et al.⁶ According to this model the solvated electron induces an enhanced water-density region of ∼1 nm diameter over which the electron is spread. Both cavity and non-cavity interpretations have been criticised and defended.^7–10 It should be noted that the one essential difference between the simulations that support either model lies in the way in which the water-excess electron interaction is modelled, a fact that testifies to the difficulty of the problem.

Experiments and simulations of singly-charged water cluster anions have been frequently employed as a means to gain a better understanding of the properties of e_aq⁻, and at the same time have raised interesting questions of their own, such as where does the excess electron reside, inside the cluster or at its surface, how does this depend on cluster size and temperature, or what is the stability range of the cluster. We will be addressing these questions below.

One of the most common strategies employed to simulate e_aq⁻ both in bulk and cluster geometries involves a combination of classical molecular dynamics (MD) for the water degrees of freedom with a quantum mechanical treatment of the electron, which is generally known as quantum-classical MD. While this may be expected to be a good approximation, it is nevertheless desirable to quantify the effects of treating the water molecules classically, given the well-known importance of nuclear quantum effects in water and aqueous systems,^11–18 and especially in view of the fact that employing different models for the water–e_aq⁻ interaction can lead to structurally different motifs (cavity and non-cavity models). It is also known that the excess electron can occupy either an internal or an external state, depending on such factors as cluster size and temperature; the transition from one type of state to the other may be (and indeed we will show below that it is) affected by the inclusion or otherwise of quantum effects in the molecular dynamics. Fully quantum simulations of water cluster anions have been reported before, such as in the pioneering work of Berne and coworkers,^19,20 and Barnett et al.,^21,22 works which treated both water molecules and excess electron within the Path Integral (PI) formalism. Although feasible, this approach is technically inconvenient, given the disparity of the masses involved: the adequate PI description of the excess electron requires many more beads than the oxygen and hydrogen atoms. An alternative approach is to solve the Schrödinger equation for the electron moving in the field resulting from its interaction with the water molecules, and use the PI formalism only for the molecular degrees of freedom. This methodology was, to our knowledge, first used by Takanayagi et al.²³ to perform simulations of (H₂O)₅⁻ and (D₂O)₅⁻, and in this sense our work is a continuation of theirs. Here we report a series of simulations of both neutral and singly-charged water clusters of different sizes and at different temperatures. In order to gauge the relevance of quantum effects in these systems, we perform both classical and quantum simulations of the neutral clusters, and QCMD as well as fully quantum simulations for the water cluster anions.

The paper is structured as follows: in Section II we provide a succinct description of our computational procedure; full details as well as convergence tests are provided in Appendix A. Our results are presented in Section III, focusing first on neutral water clusters, and then on the singly-charged water cluster anions. Finally, we draw our conclusions in Section IV.

II. Computational methodology

A detailed description of our computational methodology, including validation and convergence tests is deferred to Appendix A. Here we will confine ourselves to providing a birds-eye view of our simulation procedure.

One of our aims has been to perform both classical and quantum simulations of the equilibrium properties of water clusters and singly-charged water cluster anions, i.e. water clusters incorporating an excess electron. In the latter class of system the excess electron is always treated quantum mechanically by numerical solution of the Schrödinger equation for an electron moving in the field generated by its interaction with the water molecules. The molecules were simulated either classically or quantum mechanically within the Feynman Path Integral (PI) formalism. As is well known, the PI formalism approximately maps the partition function of a quantum system onto that of a classical isomorph, namely a ring polymer with P monomers, in which each monomer is a classical counterpart of the original system affected by a scaled potential, V/P, and coupled to its two nearest neighbors on the ring via harmonic springs. In the limit P → ∞ one can sample the Boltzmann statistics of the quantum system through its classical isomorph, while in the opposite limit, P → 1 one falls back onto the classical description. Convergence with respect to the value of P used in the simulations of the quantum systems has been monitored and is reported in Appendix A.

The water intra- and inter-molecular interactions are described by means of the SPC/F model.^24,25 Harmonic springs account for the O–H bond and H–O–H bond angle dynamics; the intermolecular interactions consist of a Coulomb contribution between partial charges located at the atoms (q_O = −0.82e = −2q_H, where e is the quantum of charge), and a Lennard-Jones potential acting between oxygen atoms. While more recent and accurate models exist to describe this system (e.g. the q-TIP4P/F,¹³ or the SPC/Fw²⁶), we have nevertheless opted for the SPC/F model, in part because there is a large body of classical simulation work with this model against which we can compare our results, and also because the influence of quantum effects on the ionic degrees of freedom should not depend critically on the model used.

In cluster anions the excess electron interaction with the oxygen and hydrogen atoms of water molecules is described via the Turi–Borgis^27,28 local pseudopotential, which is parametrized to reproduce the ground-state energy and the optical absorption spectrum of the hydrated electron with classical MD simulations of bulk water at 298 K. In the limit of long distances this model reproduces the Coulomb interaction between electron and atomic partial charges, while at short distances it mimics the interaction of the electron with corresponding neutral atoms, thus avoiding the Coulomb divergence. The model also includes a term accounting for the induced polarization of oxygen atoms brought about by the presence of the excess electron. In total, the potential energy accounting for the electron-water interaction has the following expression:

where r⁽ⁿ⁾_i = |r − r⁽ⁿ⁾_i|, and r⁽ⁿ⁾_i is the position of atom i in molecule n; the sum extends over all molecules in the cluster. The explicit form of the terms appearing in eqn (1) can be found in ref. 27 and 28.

We used standard Molecular Dynamics (MD) to sample the thermal properties of either the classical systems or their quantum counterparts, using a time step of 1 fs and coupling a Bussi–Parrinello²⁹ thermostat to each degree of freedom in order to ensure appropriate sampling of the canonical ensemble.

III. Results and discussion

In this section we report the results of the classical and quantum molecular dynamics simulations of neutral and negatively charged water clusters with n from 20 to 573 and temperatures from 50 to 400 K. First, in Section IIIA, we consider the case of neutral clusters, focusing in particular on the observed differences in energetic and structural equilibrium properties between the classical and quantum treatment. Then, in Section IIIB we report the results obtained for clusters charged with an excess electron.

A. Neutral water clusters

We have computed the energetic and structural properties of neutral water clusters (H₂O)_n with n in the range 20 to 237, both by means of classical and PI MD simulations at temperatures between 50 and 400 K.

Let us first consider the dependence of the energetic properties on temperature for some representative cluster sizes, as shown in Fig. 1. Panel (a) displays the kinetic energy per molecule for n = 1 and n = 237. It can be seen that in the classical treatment the kinetic energy is independent of cluster size, displaying a linear temperature behavior, as expected. In contrast, the quantum results do show a small but noticeable size dependence, with the results for intermediate cluster sizes (not shown) falling between the curves for n = 1 and 237. It can also be seen that, while the classical results extrapolate to zero at T = 0 K, this is not so for the quantum results, which tend to a size-dependent constant value, namely the zero-point energy (ZPE). In panel (b) of Fig. 1 we show the potential energy, again as obtained from both classical and quantum simulations, for sizes n = 47 100 237. The difference in average potential energy between the quantum and classical results for a given size and temperature corresponds closely to the ZPE, implying that the quantum clusters are predicted to be somewhat less energetically stable than their classical counterparts, at least at low temperatures.

Fig. 1 Averages of the kinetic energy (a), potential energy (b) and binding energy (c) per water molecule as a function of temperature, for quantum (solid line and filled symbols) and classical (dashed line and empty symbols) simulations of neutral (H₂O)_N clusters.

The formation or binding energy of a neutral cluster is defined as

ΔE_n = E[(H₂O)_n] − nE(H₂O), (2)

where E(H₂O) is the average total energy of a single water molecule. Binding energies per molecule are shown in Fig. 1c. As can be noted, the classical binding energies are always lower (higher in absolute value) than the quantum ones for the same temperature and cluster size, testifying to the lower stability of the latter. Although the classical and quantum curves approach one another at higher temperatures, the differences are appreciable in all the temperature range considered here. We note that the quantum binding energy decreases (increases in absolute value) at temperatures between 50 and 250 K, even if the quantum potential energy [see Fig. 1b] is monotonically increasing at all temperatures. In contrast, the classical binding energy remains essentially constant up to 150 K, increasing linearly above 200 K.

Next we consider some structural properties of (H₂O)_n clusters. Fig. 2a and b display the calculated temperature dependence of intra-molecular O–H distances and H–O–H bond angles. As can be seen, both the values as well as their temperature behavior are different in the classical and quantum results. Although the differences are not large, they are nevertheless appreciable; it is particularly noticeable that the temperature dependence of the O–H distance displays opposite trends, decreasing for the classical results, while increasing in the quantum simulations. The bond-angle temperature behavior is increasing in both cases, although the average values are almost two degrees larger in the quantum case at low temperatures. The results of both types of simulations can be seen to converge to the same values in the limit of high temperatures, as expected, but differences are still appreciable at the highest temperature considered here, namely 400 K.

Fig. 2 Averages of the intra-molecular O–H distance (a) and H–O–H angle (b), and cluster radius (c) of neutral (H₂O)_N clusters, as a function of the temperature, for quantum (solid line and filled symbols) and classical (dashed line and empty symbols) cases.

In order to quantify the size of clusters at different temperatures we calculate their gyration radius, shown in Fig. 2c for the particular case n = 237. At low temperatures the effective gyration radius of the quantum cluster is larger than the classical one, which is consistent with the higher degree of delocalization and lower stability in the quantum counterpart. Interestingly, however this trend is reversed at temperatures near and above 250 K. Although not shown in Fig. 2c, this effect is also noticeable for other cluster sizes at slightly different temperatures. We attribute this change in the classical results to a solid–liquid phase transition. Indeed, this seems to be corroborated by the fact that the slope of the mean squared displacements of the water molecules in the classical system increase by nearly an order of magnitude upon raising the temperature from 200 to 250 K, giving a diffusion coefficient of ∼4.3 × 10⁻⁹ m² s⁻¹, a value comparable to that obtained from simulations of super-cooled water at the same temperature.³⁰

Quantum and thermal effects are also observable in the pair radial distribution functions (RDFs). In Fig. 3 we report the short- and long-range regions of the oxygen-oxygen RDF, g_OO(r), of the n = 237 cluster for the classical and quantum cases at temperatures ranging from 100 to 400 K. The first peak of the g_OO function at 100 and 200 K [see panel (a)] at r ≈ 2.75 Å and the second at r ≈ 4.5 Å correspond to the first and second neighbor oxygen–oxygen distances in ice, respectively. As temperature is increased the intensity of the peaks reduces and they become broader; also the minima located between the first- and second-nearest neighbour peaks become gradually filled. Likewise, higher temperatures result in more extended tails at long distances, as seen in panel (b). These effects are present in both the quantum and classical simulations; however, if we compare the quantum and classical results in more detail, we can see that at short distances the quantum analysis results in broader peaks of lower intensity, an effect that is particularly noticeable at lower temperatures. This broadening is to some extent equivalent to an increased temperature effect. At higher temperatures, the distinction is more apparent in the long-distance tails of the RDF distributions (panel b), where the classical clusters are seen to be slightly more extended, which is consistent with a greater effective radius.

Fig. 3 Short (a) and long (b) range regions of the radial pair distribution function of the oxygen–oxygen (g_OO) distance for quantum (solid line) and classical (dashed line) simulations of the neutral (H₂O)₂₃₇ cluster.

Thus, taking into account the previous analysis of the RDFs and that of the cluster effective radius seen in Fig. 2c, we can conclude that the classical description leads to more highly structured and compact clusters at temperatures below 200 K, a general conclusion that is in agreement with previous work,^11,31 whereas the reverse is true for temperatures above 200 K.

B. Water cluster anions

Now we turn to consider the properties of singly-charged water cluster anions, (H₂O)_n⁻; specifically, we report a series of classical and quantum simulations with n = 20, 32, 47, 76, 100, 139, 237, 573 at temperatures of 100, 200 and 300 K.

In Fig. 4 we show some representative configurations obtained from the quantum simulations of cluster anions. Panels (a) and (b) show two configurations of the (H₂O)₁₀₀⁻ anion at 100 and 300 K, respectively. As can be seen, the electron delocalises over a volume that is comparable to that of the entire cluster. At low temperature the excess electron is accommodated inside the cluster, but at 300 K it is attached to its surface. At 200 K (not shown) the electron remains inside the cluster, but locates itself closer to the surface. Larger cluster sizes seem to stabilize the electron inside, as shown in panels (c) and (d), where we display instantaneous configurations of (H₂O)₅₇₃⁻ at 100 and 300 K. It should be noted that at the higher temperature these structures are probably unstable; during the course of our simulations we see that clusters of size n ≤ 100 occasionally evaporate molecules form the surface at 300 K; the same probably happens to larger clusters over longer periods of time, but all cluster sizes remain intact over the course of our simulations for temperatures of 200 K and lower. For cluster anions larger than (H₂O)₁₀₀⁻ the electron is always observed to remain inside the cluster, but moving off-center towards the surface, as can be appreciated in panel (d).

Fig. 4 Representative configurations of water cluster ions obtained from quantum simulations. Water molecules are shown at their centroid positions; the electron charge density is represented as concentric spheres, the outer one containing ∼90% of the charge density, the intermediate one ∼50%, and the inner one ∼2%. (a) 100 molecules, T = 100 K; (b) 100 molecules, T = 300 K; (c) 573 molecules, T = 100 K; (d) 573 molecules, T = 300 K.

Throughout each simulation, we compute the average distance between the electron mean position and the cluster's center of mass, , with r̄ being the mean position of the electron density. As mentioned in the introduction section, and seen in Fig. 4, the excess electron can be found either in an interior state (within the cluster) or attached to the surface, depending on factors such as cluster size and temperature. In order to classify the electron as being in one or the other state, it is helpful to adopt the following criteria:

where R_c is the gyration radius of the cluster and R_e is that for the electron [see eqn (A9)]. The quantity d_e + R_e can thus be interpreted as the outer reach of the electron, that is, the maximum distance from the cluster's center of mass where the electron density has an appreciable value.

Concerning differences between the fully quantum and quantum-classical treatment of the cluster anions, we can say that in general we obtain qualitatively similar trends for the excess electron properties, but some quantitative differences can be identified. First, we classify the electronic states as interior or surface states, following the criteria of eqn (3). In panels (a) and (c) of Fig. 5, we display the sum d_e + R_e, which, as argued above gives an indication of the outer reach of the excess electron along the cluster radius, as a function of the cluster size. The crossing of d_e + R_e and the cluster radius indicates a transition from an exterior to an interior state. For classical simulations of Fig. 5c, small clusters (n ≤ 47) can only accommodate an exterior state at all temperatures, while bigger clusters (n ≥ 190) can host the excess electron in an interior configuration. The transition from surface to interior states occurs at different sizes depending on the temperature. In general, when increasing the temperature, the transition takes place at larger sizes. At 100, 200 and 300 K, the transitions appear for n in the ranges 47–76, 76–100 and 139–190, respectively. Outside the transition region of intermediate sizes of n = 47–190, the differences between temperatures are very small. For the quantum simulations, for which results are plotted in Fig. 5a, the most remarkable difference with the classical results is that the transition from surface to interior states occurs at smaller cluster sizes. At 100, 200 and 300 K, the discontinuity appears for n in the ranges 20–32, 32–47 and 100–139. We attribute this to the fact that the quantum simulations allow atoms to explore classically forbidden regions, thus increasing the possibility of stabilization of the excess electron into the more strongly bound interior state. We should point out that the results plotted in Fig. 5 should be taken as qualitative trends, because the criteria of eqn (3) rely on the assumption that the cluster geometry is spherical, which may not be the case for small clusters and/or high temperatures. Nevertheless the observed trends appear to be robust.

Fig. 5 (a) Comparison of the excess electron reach, d_e + R_e and cluster radius, R_c, as a function of the cluster size for different temperatures. The solid and dashed lines represent R_c at 100 and 300 K of temperature, respectively. (b) Gyration radius of the excess electron (in color) along with the simulated value for bulk water at 298 K (dashed black line) of ref. 28. Panels (c) and (d) show the same quantities as (a) and (b) for the classical case.

The excess electron gyration radius, which measures its dispersion around its mean position, is also characteristic of each type of configuration. In panels (b) and (d) of Fig. 5 we see how the gyration radius evolves towards the bulk value of 2.42 Å obtained by Turi and Borgis for water at 298 K in ref. 28, with the same model used here. This value agrees well with the experimental moment analysis of the absorption spectrum,^32,33 which results in a radius of 2.5 Å in bulk water. As can be observed, the gyration radii for exterior configurations approach the bulk value from above as the cluster grows, whereas interior states have a nearly constant value, regardless the temperature or the cluster size.

In order to gain insight into the mechanisms resulting in the different configurations of the excess electron, in Fig. 6 we plot the projection of water molecule dipoles onto the cluster radial direction vs. distance from the cluster center of mass for a neutral cluster (panel a), and onto the axis separating water molecule and electron center of mass for the cluster anion of the same size (panel b), at two different temperatures. In the SPC/F water model the molecular dipole lies close to the H–O–H angle bisector, pointing in the direction of the hydrogen atoms. Firstly, let us consider the orientation of water molecule dipoles in the neutral n = 100 cluster, in Fig. 6a. At T = 100 K the dipole projection onto the cluster radial direction displays an oscillating character consistent with the interpretation of the cluster being solid. Indeed, at 300 K the oscillating character has been considerably washed out, although still appreciable to some extent, indicating that the cluster is liquid. It can be seen that in both cases the outer layer of the cluster has molecules preferentially oriented such that the dipole, i.e. the hydrogen atoms, point away from the cluster. In the case of the cluster anions, Fig. 6b, the nearest water molecules to the electron are oriented such that their dipoles point towards the electron center of mass, as indicated by the positive values of the projection in the direction r_w–e = r̄_e − r_w at short distances. This is the case for both internal and external electronic states.

Fig. 6 Projection of the dipolar moment of water molecules for the n = 100 neutral cluster (a) and in the presence of the excess electron (b). For the neutral cluster, the projection is on the radial direction, that is, the line between the center of mass of the cluster and each specific water molecule. For the negatively charged cluster, the dipolar moment is projected on the direction linking the mean position of the excess electron with each water molecule.

The orientation of the molecular dipoles in the outer layer of the neutral clusters, see Fig. 6a, generates a surface dipole that can serve as a trap for an excess electron. This is shown by the radial electronic potential profiles of panel (a) and (b) of Fig. 7, obtained from classical and quantum simulations, respectively. Both exhibit a potential energy well close to the cluster surface. In the presence of the excess electron (panels c and d), there is a distortion of the molecular network that can generate a stable region inside the cluster, at shorter distances from the cluster center; this is akin to the well-known polaron effect in extended systems. For the classical simulation of (H₂O)₁₀₀⁻, see Fig. 7c, an interior state is only stable at temperatures of 100 and 200 K, but becomes metastable at 300 K. In the quantum case (panel d) the situation is similar, though the internal state remains slightly more stable as the temperature is increased than in the classical case. In contrast, the external potential trap seems to be more robust against temperature increase: the surface potential minimum becomes shallower as the temperature is increased, but is less affected and remains attractive at all temperatures considered here. Comparison between the results obtained from classical and quantum simulations [Fig. 7c and d] indicate that in the latter case the interior state is somewhat more robust against temperature increase, which is consistent with our earlier observations [see Fig. 5].

Fig. 7 Radial electron potential from the center of the water cluster for the classical neutral (a), quantum neutral (b), classical charged (c) and quantum charged (d) clusters with n = 100 water molecules.

Let us now focus on the energetic properties of the excess electron. Its average ground state energy, ε₀, can be compared with the negative of the experimental vertical detachment energy (VDE) as obtained via photoelectron spectroscopy.^34–36 The VDE is expected to have a linear dependence with the inverse cluster radius (or equivalently, n^−1/3) for both surface or interior states.³⁷ This was indeed observed in the experimental results of Coe et al.,³⁴ whose linear fit is shown in Fig. 8 as a dashed black line. Our average values for ε₀ for surface states appear above the experimental line, while those of internal states are found below it, in agreement with the simulations of Turi et al.³⁸ This is the case for both the classical (a) and quantum (b) simulations, the only difference being that the latter evidence a transition to the internal state at smaller cluster sizes, consistent with our earlier discussion [see Fig. 5]. Extrapolation of our calculated average ground state energies for the interior states to the bulk limit results in an energy close to −4 eV, somewhat lower than the value of −3.1 eV found by Turi et al.,³⁸ which is closer to the experimental value of −3.3 eV.^39,40 This discrepancy could at least in part arise from the fact that the model we are using was originally designed to describe the excess electron in bulk water, and has not been readjusted for the case of cluster geometries. Out of the classical and quantum results, the latter are somewhat closer to the experimental value, but still below it.

Fig. 8 Classical (a) and quantum (b) average of the simulated ground-state energy of the excess electron for (H₂O)_n⁻ clusters at different temperatures. The experimental fit of Coe from ref. 34 is shown with dashed black lines, while the linear fits from our interior-state data are represented with dotted black lines.

In Fig. 9 the calculated energy gaps between first excited and ground state, ε_gap = ε₁ − ε₀, are shown together with the experimental fit obtained from the maxima of the optical adsorption spectra by Ayotte et al.⁴¹ It is evident from the figure that only the external-configuration states agree with the experimental fit, which indicates that in the experiments of Ayotte et al., conducted in cluster sizes in the range 6 to 50, only surface states were observed. It is noteworthy that the extrapolations of the linear fit to the calculated band-gap for both interior and exterior states to infinite size fall very close to the experimentally observed value of 1.72 eV at 298 K.^42,43

Fig. 9 Classical (a) and quantum (b) average of the simulated energy gap of the excess electron for (H₂O)_n⁻ clusters at different temperatures. The experimental fit of Ayotte from ref. 41 is shown with dashed black lines, while the linear fits from our interior-state data are represented with dotted black lines.

The most remarkable difference between the classical and quantum simulations is the appearance of a transition from surface to internal states at smaller sizes in the latter case. This trend is observed at all temperatures considered here. Specifically, at 100 K, no interior configurations of the excess electron are observed for clusters smaller than n = 76, with the transition occurring somewhere between 47 and 76. In the quantum simulations, however, the transition appears to occur in the range n = 20–32. This observation is in good agreement with the conclusions of Neumark,⁴⁴ who carefully analyzed the experimental VDE data previously assigned to internal states,³⁵ concluding that the transition occurs in the size range n = 25–35. Furthermore, this assessment is in full agreement with other low temperature experiments^36,45 and consistent with theoretical works,^46,47 which observe internal states for clusters with n < 50.

IV. Conclusions

To summarize, we have reported results of simulations of neutral water clusters and cluster anions containing an excess electron, with cluster sizes in the range 20 to 573 molecules. We have treated the water molecules both classically and quantum mechanically via the path integral formalism. Our results indicate that quantum effects in the molecular degrees of freedom generally have a small but noticeable influence on the cluster properties, and in particular they affect the cluster size at which the transition from external to internal state of the excess electron occurs, allowing it to happen at smaller cluster sizes when they are included. Therefore, this work reinforces the importance of including nuclear quantum effects to obtain a precise picture when modelling water clusters in the presence of an excess electron.

Author contributions

A. G. developed the PIMD code and carried out the simulations. E. R. H. conceived and supervised the project, and implemented the solution of the Schrödinger equation. Both authors jointly wrote the manuscript.

Conflicts of interest

The authors declare that they have no conflict of interests.

Appendix A: computational details

Simulation procedure

One of the central aims of this work is to account for the quantum effects of the atomic dynamics in the structural and energetic properties of water clusters, and especially those of water cluster anions, where we anticipate that such effects could be most appreciable, particularly on the excess electron properties. To do so, we employ the Feynman Path Integral (PI) approach.^57–66 As is well known, the PI formalism establishes a quantum-classical isomorphism, through which the partition function of a quantum system can be approximately mapped onto that of a classical system consisting of P replicas (also referred to as beads) of the original system, but subject to a modified potential, as follows:Here, Q_P is the partition function of the P-bead classical isomorph approximating that of its quantum counterpart; β = (k_BT)⁻¹, and V(r^α) is the potential energy. Latin sub-indices label atoms, while Greek super-indices label beads; if only a super-index is given, the variable is meant to represent the entire configuration for that bead, i.e.r^α ≡ (r^α₁, r^α₂,…,r^α_n). The isomorphism is exact in the limit P → ∞. For a finite number of replicas, thermodynamic estimators for the internal energy E, its kinetic and potential contributions, and a series of other properties can be easily derived^48–50 from the partition function eqn (A1). Obtaining average equilibrium properties of the quantum system then reduces to sampling the classical isomorph in the canonical ensemble, either with Monte Carlo or with molecular dynamics techniques.

In the case of the neutral clusters V(r) in eqn (A1) reduces to V_SPC/F(r), i.e. the SPC/F potential energy expression. However, for the cluster anions the interaction of water molecules with the excess electron must be incorporated into the model. Next we describe how we have done this in our simulations, more details can be found in ref. 51.

Schrödinger equation for the excess electron

In principle it would be possible to address the dynamics of water molecules and excess electron on an equal footing within the PI formalism. However, the light mass of the electron would require a large number of beads in the PI ring polymer, much larger than would be required in the absence of the electron. While it is possible to effectively use different numbers of beads for different degrees of freedom,^19,20 we have found it more convenient to invoke the Born–Oppenheimer approximation and treat separately the dynamics of molecules and excess electron. Thus we use the PI method to account for the quantum effects in the dynamics of the water molecules subject to their mutual interaction and that of the electron density charge, while the latter is obtained for each cluster configuration by direct numerical solution of the electron's Schrödinger equation.

To obtain the electron charge distribution we proceed as follows: first, the interaction potential of the electron with the water molecules, eqn (1), is discretized on a regular cubic grid covering a volume large enough to contain the water cluster, with its center fixed at the cluster center of mass. In this work we have used a grid of size L = 80 Å with a grid of 64 × 64 × 64 grid points, which provides sufficient resolution, as we discuss below in the convergence tests section. The lowest 4 electron eigen-pairs are then obtained iteratively using the imaginary time propagation method.⁵² In this method, a set of n trial wave functions ψ_i is driven towards the n lowest lying eigenstates of Hamiltonian H by acting on them with the imaginary time evolution operator:

ψ_i(ε) = e^−εHψ_i(0). (A2)

Here ε is the imaginary time-step. Notice that the complex time evolution operator is not unitary, and hence does not preserve the orthonormality of the trial set; it is therefore necessary to orthonormalise the set following each time-step. After a sufficiently large number of evolution steps (followed by orthonormalisation), this process converges to the lowest n eigenstates of H. Once these eigenstates have been obtained, the excess electron charge density is obtained.

The action of operator e^−εH on the trial set cannot be calculated exactly, but it can be approximated up to any desired order in ε⁵² using a Trotter-like factorisation:

where c_k, a_jk, and b_jk are numerical coefficients, and H = T + V, where T and V are the kinetic and potential energy operators, respectively. Operators e^αT are diagonal in momentum space, while operators e^βV are diagonal in real space (provided V is local, as is the case in our model), so the action of eqn (A3) on the trial set involves a succession of Fast Fourier transformations⁵³ between momentum and real space (see ref. 52 for details). In practice we find that a 4th-order approximation to e^−εH provides a good compromise between computational effort and speed of convergence.

In the case of the cluster anions V(r) in eqn (A1) takes the form:

V(r^α) = E_e(r^α) + V_SPC/F(r^α), (A4)

where E_e(r^α) is the electron ground state energy for bead α, which depends parametrically on r^α, and V_SPC/F(r^α) is the SPC/F model potential for water, as in the neutral clusters.

Centroid approximation

Strictly speaking, the ground-state energy, E_e(r^α), must be evaluated separately for each replica α, and then accumulated with weight 1/P. However, as argued below, it is possible to take a further approximation in our scheme, which we refer to as the centroid approximation. Each quantum atomic nucleus will be delocalized with certain mean position and spatial extension, characterized by its centroid (instantaneous center of the chain) and its gyration radius:The centroid approximation amounts to solving the electron Schrödinger equation for the ring polymer centroid instead of doing it for each individual bead.

A schematic representation of a quantum particle with 8 beads around its centroid position can be seen in left panel of Fig. 10.

Fig. 10 (left) Schematic representation of the beads, r^α, centroid position, r̄, and gyration radius, r_g, of a classical ring polymer corresponding to a quantum particle with 8 replicas. (right) Typical geometry of a 3 water molecules cluster, simulated with 40 beads.

Like the energy, the forces on the molecular atoms in the cluster anion simulations have two contributions: one resulting from the water–water interactions, modelled via the SPC/F model, and the second from their interaction with the excess electron charge density, which can be easily evaluated from the Hellman–Feynman theorem. Let us consider the electron energy at bead α and write it in terms of that evaluated for the centroid; we would have:

E_e(r^α) = E_e(r̄) + ∇E_e(r̄)·(r^α − r̄) + [scr O, script letter O](|r^α − r̄|²). (A7)

The centroid approximation implies taking E_e(r^α) ≈ E_e(r̄), which will be justified when the second and subsequent terms on the rhs of eqn (A7) are negligible compared to E_e(r̄). This will happen if the Hellman–Feynman forces are small and/or the gyration radius is small. Typically, the distances between bead coordinates of the same quantum ion are much smaller than the internuclear distances, as can be appreciated in the right panel of Fig. 10, and so we can expect the approximation to hold valid. Within this approximation, the electronic contribution to force acting on the ith atom in replica α due to the electron cloud is given byas can be seen using the chain rule, where F̄_e,i is the electronic force acting on the centroid of atom i.

From a practical point of view, the centroid approximation reduces considerably the computational cost needed to obtain the electronic energy and forces, as this is done only for the centroid configuration and not for every bead on the ring polymer. Nevertheless, the approximation needs to be justified by comparing its predictions against those of the rigorous procedure, which will be done later.

Molecular dynamics sampling

In order to sample the partition function eqn (A1) and obtain thermal averages derived from it we employ molecular dynamics simulations at constant temperature. To ensure sampling of the canonical ensemble we combine the ring polymer dynamics with a Bussi–Parrinello²⁹ thermostat attached to each degree of freedom. We used a time step of 1 fs and a friction parameter γ = 0.001 a.u., which provided effective sampling.

For classical simulations of the neutral clusters, initial geometries were obtained from large molecular dynamics simulations of bulk water at 300 K and density 1 g cm⁻³, including all water molecules within spheres of varying radii.⁵⁴ For each calculation, an equilibration run of 1 × 10⁵ steps (100 ps) was executed before the production simulation of 5 × 10⁵ steps. Quantum PIMD simulations were carried out with 128 replicas for all temperatures between 50 and 400 K, using classical equilibrated configurations as initial conditions. As will be justified later, the election of P = 128 provides converged properties and it is computationally affordable for neutral clusters.

For water cluster anions, the simulations were initiated from interior states, obtained from neutral configurations with a relaxation simulation in the presence of a smooth confining potential to keep the excess electron in the vicinity of the cluster center, , with k = 10⁻⁸ a.u. The relaxation was performed at each temperature during 10⁴ steps, and after that the confining term was switched off. Typical calculations consisted of an equilibration run of 5 × 10⁴ steps (50 ps), previous to the production run with length between 10⁵ and 2 × 10⁵ steps, which correspond to 100 and 200 ps of time. This procedure, combined with the length of our simulations, allows the average total energy to converge to its equilibrium value. For the quantum simulations of the negatively charged water clusters, we used 128, 64 and 48 replicas for temperatures of 100, 200 and 300 K, respectively, trying to keep constant the product PT as close as possible. In addition, we employed the centroid approximation, which provides accurate results for both the water cluster and the excess electron properties and its computational needs are up to two orders of magnitude lower than the exact calculation at 100 K. The election of those number of replicas and the use of the centroid approximation is justified below in the next section.

In all cases, the associated uncertainties of equilibrium averages were computed with block averaging, as the plateau of the standard error of the mean among the blocks, when increasing the block size.^55,56

Convergence tests

Here we analyze the convergence of the total energy with the number of replicas for a cluster with N = 100 water molecules at different temperatures. As anticipated, the number of replicas required to converge increases with decreasing temperature, as can be observed in the top panel of Fig. 11. While P = 32 is sufficient at T = 400 K, a similar level of convergence at T = 50 K requires 256 replicas. The differences between using 128 and 256 replicas are approximately 5% or less for temperatures equal to or higher than 100 K. Consequently, we performed all subsequent simulations at a fixed PT = 12 800 K value, which implies using 128 replicas at T = 100 K, 64 at T = 200 K, and 48 at T = 300 K.

Fig. 11 Total energy per molecule as a function of the number of replicas, for several temperatures.

For cluster anions we also require a suitable grid with appropriate size and grid spacing on which to discretize the electron wave functions and density. For this purpose we analysed the convergence of the ground and first excited state energies for an electron in a surface state in the n = 100 cluster; interior states are more localized and show better convergence with the grid parameters. The results, obtained with the complex time evolution using an energy tolerance of 10⁻⁴ hartree ≈ 2.7 × 10⁻³ eV, are shown in Table 1. As can be seen there, a cubic grid extending over 80 Å and with 64 points along each dimension is capable of providing converged eigenenergies with an error smaller than 10⁻² eV. A grid twice as large in volume but with equal grid point density provides essentially indistinguishable results. Therefore, for all subsequent simulations of cluster anions we employed L = 80 Å and N_p = 64 for all sizes and temperatures.

Table 1 Ground and first excited state energies of the excess electron as a function of the grid size L and the number of grid points N_p, in each Cartesian direction. All the calculations correspond to the n = 100 neutral cluster. Energies are in eV and lengths in Å

Surface state, N =100, T = 300 K
	N p = 32		N p = 64		N p = 128
	ε 0	ε 1	ε 0	ε 1	ε 0	ε 1
L = 40	−1.2596	−0.4268	−1.2393	−0.4009	−1.2285	−0.3872
L = 80	−1.3454	−0.5995	−1.3560	−0.5944	−1.3560	−0.5944
L = 160	−1.5857	−0.8469	−1.3511	−0.6025	−1.3561	−0.5961

---

As discussed above, the centroid approximation allows us to reduce the computational cost of the simulations, but its validity must be checked by comparing its predictions against those of the rigorous procedure. To do so, we have performed 20 ps-long simulations (2 × 10⁴ time steps) for the n = 100 cluster anion at different temperatures with 128 replicas, starting from previously equilibrated configurations. In Fig. 12 we plot the relative differences in electron energy and spread, defined as

where ρ is the electron density, between the simulations with and without centroid approximation. The relative difference decreases from 2.5 to 0.5% between temperatures from 50 to 400 K, due to the fact that the kinetic energy spring constant increases with temperature (∝T²), resulting in more localized beads around the centroid. For structural and energetic properties of the entire cluster anion the errors are even smaller (<0.5%).

Fig. 12 Relative error made in the averages when using the centroid approximation to study the (H₂O)₁₀₀⁻ cluster, with respect to the exact quantum simulation. Panel (a) shows the error in the intramolecular HOH angle and the OH distance, and the estimated cluster radius. Panel (b) displays the error in the total energy of the system and the gyration radius of oxygen and hydrogen atoms. Lastly, panel (c) presents the differences in the excess electron ground-state energy, ε₀, and its spatial width, R_e.

As can be seen in Table 1, the energy gap between ground and first excited state of the excess electron in cluster anions remains significantly larger than k_BT even at the largest temperature considered for these systems (300 K); thus it is safe to assume that the electron remains in its ground state throughout.

Acknowledgements

We thank J. M. Soler, M. V. Fernández-Serra, R. Ramírez and C. P. Herrero for helpful discussions. The authors acknowledge the allocation of computer time provided by the CESGA supercomputing centre and the Spanish Supercomputing Network (RES) through project FI-2020-2-0009. This work has been funded by MCIN/AEI/10.13039/501100011033 through project PGC2018-096955-B-C44.

References

1. E. J. Hart and J. W. Boag, J. Am. Chem. Soc., 1962, 81, 4090.
2. R. M. Young and D. M. Neumark, Chem. Rev., 2012, 112, 5553.
3. L. Turi and P. J. Rossky, Chem. Rev., 2012, 112, 5641.
4. J. M. Herbert and M. P. Coons, Annu. Rev. Phys. Chem., 2017, 68, 447.
5. J. M. Herbert, Phys. Chem. Chem. Phys., 2019, 21, 20538.
6. R. E. Larsen, W. J. Glover and B. J. Schwartz, Science, 2010, 329, 65.
7. L. Turi and A. Madarász, Science, 2011, 331, 1387.
8. L. D. Jacobson and J. M. Herbert, Science, 2011, 331, 1387.
9. R. E. Larsen, W. J. Glover and B. J. Schwartz, Science, 2011, 331, 1387.
10. J. R. Casey, A. Kahros and B. J. Schwartz, J. Phys. Chem. B, 2013, 117, 14173.
11. R. A. Kuharski and P. J. Rossky, J. Chem. Phys., 1985, 82, 5164.
12. F. Paesani, W. Zhang, D. A. Case, T. E. Cheetham and G. A. Voth, J. Chem. Phys., 2006, 125, 184507.
13. S. Habershon, T. E. Markland and D. E. Manolopoulos, J. Chem. Phys., 2009, 131, 024501.
14. F. Paesani, S. Yoo, H. J. Bakker and S. S. Xantheas, J. Phys. Chem. Lett., 2010, 1, 2316.
15. B. S. González, E. G. Noya, C. Vega and L. M. Sesé, J. Phys. Chem. B, 2010, 114, 2484.
16. C. Vega and J. L.-F. Abascal, Phys. Chem. Chem. Phys., 2011, 13, 19663.
17. J. O. Richardson, C. Pérez, S. Lobsiger, A. A. Reid, B. Temelso, G. C. Shields, Z. Kisiel, D. J. Wales, B. H. Pate and S. C. Althorpe, Science, 2016, 351, 1310.
18. M. Ceriotti, W. Fang, P. G. Kusalik, R. H. McKenzie, A. Michaelides, M. A. Morales and T. E. Markland, Chem. Rev., 2016, 116, 7529.
19. D. Thirumalai, A. Wallqvist and B. J. Berne, J. Stat. Phys., 1986, 43, 973.
20. A. Wallqvist, D. Thirumalai and B. J. Berne, J. Chem. Phys., 1987, 86, 6404.
21. R. N. Barnett, U. Landman, C. L. Cleveland and J. Jortner, Phys. Rev. Lett., 1987, 59, 811.
22. R. N. Barnett, U. Landman, C. L. Cleveland and J. Jortner, J. Chem. Phys., 1988, 88, 65.
23. T. Takayanagi, T. Yoshikawa, H. Motegi and M. Shiga, Chem. Phys. Lett., 2009, 482, 195.
24. K. Toukan and A. Rahman, Phys. Rev. B, 1985, 31, 2643.
25. J. Lobaugh and G. A. Voth, J. Chem. Phys., 1997, 106, 2400.
26. Y. Wu, H. L. Tepper and G. A. Voth, J. Chem. Phys., 2006, 124, 024503.
27. L. Turi, M.-P. Gaigeot, N. Levy and D. Borgis, J. Chem. Phys., 2001, 114, 7805.
28. L. Turi and D. Borgis, J. Chem. Phys., 2002, 117, 6186.
29. G. Bussi and M. Parrinello, Phys. Rev. E: Stat., Nonlin., Soft Matter Phys., 2007, 75, 056707.
30. D. K. Park, B. W. Yoon and S. H. Lee, Bull. Korean Chem. Soc., 2015, 36, 492.
31. A. Wallqvist and B. J. Berne, Chem. Phys. Lett., 1985, 117, 214.
32. D. J. Lavrich, P. J. Campagnola and M. A. Johnson, in Linking the Gaseous and Condensed Phases of Matter, ed. L. G. Christophorou, E. Illenberger and W. F. Schmidt, SpringerNew York NY, 1994, pp. 183–202.
33. D. M. Bartels, J. Chem. Phys., 2001, 115, 4404.
34. J. V. Coe, G. H. Lee, J. G. Eaton, S. T. Arnold, H. W. Sarkas, K. H. Bowen, C. Ludewigt, H. Haberland and D. R. Worsnop, J. Chem. Phys., 1990, 92, 3980.
35. J. R.-R. Verlet, A. E. Bragg, A. Kammrath, O. Cheshnovsky and D. M. Neumark, Science, 2005, 307, 93 LP.
36. L. Ma, K. Majer, F. Chirot and B. von Issendorff, J. Chem. Phys., 2009, 131, 144303.
37. G. Makov and A. Nitzan, J. Phys. Chem., 1994, 98, 3459.
38. L. Turi, W.-S. Sheu and P. J. Rossky, Science, 2005, 309, 914 LP.
39. K. R. Siefermann, Y. Liu, E. Lugovoy, O. Link, M. E. Faubel, U. Buck, B. Winter and B. Abel, Nat. Chem., 2010, 2, 274.
40. Y. Tang, H. Shen, K. Sekiguchi, N. Kurahashi, Y.-I. Mizuno, T. Suzuki and T. Suzuki, Phys. Chem. Chem. Phys., 2010, 12, 3653.
41. P. Ayotte and M. A. Johnson, J. Chem. Phys., 1997, 106, 811.
42. F.-Y. Jou and G. R. Freeman, J. Phys. Chem., 1977, 81, 909.
43. F.-Y. Jou and G. R. Freeman, J. Phys. Chem., 1979, 83, 2383.
44. D. M. Neumark, Mol. Phys., 2008, 106, 2183.
45. A. Herburger, E. Barwa, M. Oncák, J. Heller, C. van der Linde, D. M. Neumark and M. K. Beyer, J. Am. Chem. Soc., 2019, 141, 18000.
46. T. Sommerfeld and K. D. Jordan, J. Am. Chem. Soc., 2006, 128, 5828.
47. T. Frigato, J. VandeVondele, B. Schmidt, P. Schütte and C. Jungwirth, J. Phys. Chem. A, 2008, 112, 6125.
48. M. J. Gillan, in Computer Modelling of Fluids Polymers and Solids, ed. C. R. A. Catlow, S. C. Parker and M. P. Allen, Springer Netherlands, Dordrecht, 1990, pp. 155–188.
49. M. E. Tuckerman, Statistical Mechanics: Theory and Molecular Simulation, Oxford University Press, 2010.
50. M. F. Herman, E. J. Bruskin and B. J. Berne, J. Chem. Phys., 1982, 76, 5150.
51. A. Gijón, Classical and quantum molecular dynamics simulations of condensed aqueous systems, PhD thesis, Digital CSIC, 2021.
52. S. A. Chin, S. Janecek and E. Krotscheck, Chem. Phys. Lett., 2009, 470, 342.
53. M. Frigo and S. G. Johnson, The Fastest Fourier Transform in the West, Tech. Rep. MIT-LCS-TR-728, Massachusetts Institute of Technology, 1997.
54. E. Rudberg, E. H. Rubensson, P. Salek and A. Kruchinina, SoftwareX, 2018, 7, 107.
55. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Oxford Scholarship Online, 2nd edn, 2017.
56. A. Grossfield, P. N. Patrone, D. R. Roe, A. J. Schultz, D. Siderius and D. M. Zuckerman, Living J. Comput. Mol. Sci., 2018, 1, 5067.
57. M. E. Tuckerman, J. Phys. Condens. Matter, 2002, 14, R1297.
58. D. Marx and M. Parrinello, Z. Phys. B: Condens. Matter, 1994, 95, 143.
59. D. Marx and M. Parrinello, J. Chem. Phys., 1996, 104, 4077.
60. M. E. Tuckerman, D. Marx, M. L. Klein and M. Parrinello, J. Chem. Phys., 1996, 104, 5579.
61. D. Marx, M. E. Tuckerman and G. J. Martyna, Comput. Phys. Commun., 1999, 118, 166.
62. E. Paquet and H. L. Viktor, Adv. Chem., 2018, 2018, 9839641.
63. R. P. Feynman, Rev. Mod. Phys., 1948, 20, 367.
64. R. P. Feynman and A. R. Hibbs, Quantum Mechanics and Path Integrals, McGraw-Hill, New York, 1965.
65. R. P. Feynman, Statistical Mechanics, Addison-Wesley, New York, 1972.
66. C. P. Herrero and R. Ramírez, J. Phys.: Condens. Matter, 2014, 26, 233201.

---