Open Access Article
This Open Access Article is licensed under a Creative Commons Attribution-Non Commercial 3.0 Unported Licence

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

A. Gijón and E. R. Hernández *
Instituto de Ciencia de Materiales de Madrid (ICMM-CSIC), Campus de Cantoblanco, 28049 Madrid, Spain. E-mail: Eduardo.Hernandez@csic.es

Received 5th March 2022 , Accepted 19th May 2022

First published on 26th May 2022


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,1 the hydrated electron, eaq, 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, eaq 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 eaq. At the theoretical level, modelling of eaq 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 eaq 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.6 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 eaq, 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 eaq 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–eaq 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.23 to perform simulations of (H2O)5 and (D2O)5, 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 (qO = −0.82e = −2qH, 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,13 or the SPC/Fw26), 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–Borgis27,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:

 
image file: d2cp01088g-t1.tif(1)
where r(n)i = |rr(n)i|, and r(n)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–Parrinello29 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 (H2O)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[thin space (1/6-em)]100[thin space (1/6-em)]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.


image file: d2cp01088g-f1.tif
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 (H2O)N clusters.

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

 
ΔEn = E[(H2O)n] − nE(H2O),(2)
where E(H2O) 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 (H2O)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.


image file: d2cp01088g-f2.tif
Fig. 2 Averages of the intra-molecular O–H distance (a) and H–O–H angle (b), and cluster radius (c) of neutral (H2O)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−9 m2 s−1, a value comparable to that obtained from simulations of super-cooled water at the same temperature.30

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, gOO(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 gOO 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.


image file: d2cp01088g-f3.tif
Fig. 3 Short (a) and long (b) range regions of the radial pair distribution function of the oxygen–oxygen (gOO) distance for quantum (solid line) and classical (dashed line) simulations of the neutral (H2O)237 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, (H2O)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 (H2O)100 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 (H2O)573 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 (H2O)100 the electron is always observed to remain inside the cluster, but moving off-center towards the surface, as can be appreciated in panel (d).


image file: d2cp01088g-f4.tif
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, image file: d2cp01088g-t2.tif, with [r with combining macron] 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:

 
image file: d2cp01088g-t3.tif(3)
where Rc is the gyration radius of the cluster and Re is that for the electron [see eqn (A9)]. The quantity de + Re 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 de + Re, 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 de + Re 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.


image file: d2cp01088g-f5.tif
Fig. 5 (a) Comparison of the excess electron reach, de + Re and cluster radius, Rc, as a function of the cluster size for different temperatures. The solid and dashed lines represent Rc 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 rw–e = [r with combining macron]erw at short distances. This is the case for both internal and external electronic states.


image file: d2cp01088g-f6.tif
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 (H2O)100, 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].


image file: d2cp01088g-f7.tif
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, ε0, 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.37 This was indeed observed in the experimental results of Coe et al.,34 whose linear fit is shown in Fig. 8 as a dashed black line. Our average values for ε0 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.38 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.,38 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.


image file: d2cp01088g-f8.tif
Fig. 8 Classical (a) and quantum (b) average of the simulated ground-state energy of the excess electron for (H2O)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 = ε1ε0, are shown together with the experimental fit obtained from the maxima of the optical adsorption spectra by Ayotte et al.41 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


image file: d2cp01088g-f9.tif
Fig. 9 Classical (a) and quantum (b) average of the simulated energy gap of the excess electron for (H2O)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,44 who carefully analyzed the experimental VDE data previously assigned to internal states,35 concluding that the transition occurs in the size range n = 25–35. Furthermore, this assessment is in full agreement with other low temperature experiments36,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:
 
image file: d2cp01088g-t4.tif(A1)
Here, QP is the partition function of the P-bead classical isomorph approximating that of its quantum counterpart; β = (kBT)−1, 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α1, rα2,…,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 derived48–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 VSPC/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.52 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 ε52 using a Trotter-like factorisation:

 
image file: d2cp01088g-t5.tif(A3)
where ck, ajk, and bjk 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 transformations53 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α) = Ee(rα) + VSPC/F(rα),(A4)
where Ee(rα) is the electron ground state energy for bead α, which depends parametrically on rα, and VSPC/F(rα) is the SPC/F model potential for water, as in the neutral clusters.

Centroid approximation

Strictly speaking, the ground-state energy, Ee(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:
 
image file: d2cp01088g-t6.tif(A5)
 
image file: d2cp01088g-t7.tif(A6)
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.


image file: d2cp01088g-f10.tif
Fig. 10 (left) Schematic representation of the beads, rα, centroid position, [r with combining macron], and gyration radius, rg, 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:

 
Ee(rα) = Ee([r with combining macron]) + ∇Ee([r with combining macron])·(rα[r with combining macron]) + [scr O, script letter O](|rα[r with combining macron]|2).(A7)
The centroid approximation implies taking Ee(rα) ≈ Ee([r with combining macron]), which will be justified when the second and subsequent terms on the rhs of eqn (A7) are negligible compared to Ee([r with combining macron]). 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 by
 
image file: d2cp01088g-t8.tif(A8)
as can be seen using the chain rule, where [F with combining macron]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–Parrinello29 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−3, including all water molecules within spheres of varying radii.54 For each calculation, an equilibration run of 1 × 105 steps (100 ps) was executed before the production simulation of 5 × 105 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, image file: d2cp01088g-t9.tif, with k = 10−8 a.u. The relaxation was performed at each temperature during 104 steps, and after that the confining term was switched off. Typical calculations consisted of an equilibration run of 5 × 104 steps (50 ps), previous to the production run with length between 105 and 2 × 105 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[thin space (1/6-em)]800 K value, which implies using 128 replicas at T = 100 K, 64 at T = 200 K, and 48 at T = 300 K.
image file: d2cp01088g-f11.tif
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−4 hartree ≈ 2.7 × 10−3 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−2 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 Np = 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 Np, 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 × 104 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

 
image file: d2cp01088g-t10.tif(A9)
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 (∝T2), 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%).


image file: d2cp01088g-f12.tif
Fig. 12 Relative error made in the averages when using the centroid approximation to study the (H2O)100 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, ε0, and its spatial width, Re.

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 kBT 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 CrossRef.
  2. R. M. Young and D. M. Neumark, Chem. Rev., 2012, 112, 5553 CrossRef CAS PubMed.
  3. L. Turi and P. J. Rossky, Chem. Rev., 2012, 112, 5641 CrossRef CAS PubMed.
  4. J. M. Herbert and M. P. Coons, Annu. Rev. Phys. Chem., 2017, 68, 447 CrossRef CAS PubMed.
  5. J. M. Herbert, Phys. Chem. Chem. Phys., 2019, 21, 20538 RSC.
  6. R. E. Larsen, W. J. Glover and B. J. Schwartz, Science, 2010, 329, 65 CrossRef CAS PubMed.
  7. L. Turi and A. Madarász, Science, 2011, 331, 1387 CrossRef PubMed.
  8. L. D. Jacobson and J. M. Herbert, Science, 2011, 331, 1387 CrossRef PubMed.
  9. R. E. Larsen, W. J. Glover and B. J. Schwartz, Science, 2011, 331, 1387 CrossRef.
  10. J. R. Casey, A. Kahros and B. J. Schwartz, J. Phys. Chem. B, 2013, 117, 14173 CrossRef CAS PubMed.
  11. R. A. Kuharski and P. J. Rossky, J. Chem. Phys., 1985, 82, 5164 CrossRef CAS.
  12. F. Paesani, W. Zhang, D. A. Case, T. E. Cheetham and G. A. Voth, J. Chem. Phys., 2006, 125, 184507 CrossRef PubMed.
  13. S. Habershon, T. E. Markland and D. E. Manolopoulos, J. Chem. Phys., 2009, 131, 024501 CrossRef PubMed.
  14. F. Paesani, S. Yoo, H. J. Bakker and S. S. Xantheas, J. Phys. Chem. Lett., 2010, 1, 2316 CrossRef CAS.
  15. B. S. González, E. G. Noya, C. Vega and L. M. Sesé, J. Phys. Chem. B, 2010, 114, 2484 CrossRef PubMed.
  16. C. Vega and J. L.-F. Abascal, Phys. Chem. Chem. Phys., 2011, 13, 19663 RSC.
  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 CrossRef CAS PubMed.
  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 CrossRef CAS PubMed.
  19. D. Thirumalai, A. Wallqvist and B. J. Berne, J. Stat. Phys., 1986, 43, 973 CrossRef.
  20. A. Wallqvist, D. Thirumalai and B. J. Berne, J. Chem. Phys., 1987, 86, 6404 CrossRef CAS.
  21. R. N. Barnett, U. Landman, C. L. Cleveland and J. Jortner, Phys. Rev. Lett., 1987, 59, 811 CrossRef CAS PubMed.
  22. R. N. Barnett, U. Landman, C. L. Cleveland and J. Jortner, J. Chem. Phys., 1988, 88, 65 Search PubMed.
  23. T. Takayanagi, T. Yoshikawa, H. Motegi and M. Shiga, Chem. Phys. Lett., 2009, 482, 195 CrossRef CAS.
  24. K. Toukan and A. Rahman, Phys. Rev. B, 1985, 31, 2643 CrossRef CAS PubMed.
  25. J. Lobaugh and G. A. Voth, J. Chem. Phys., 1997, 106, 2400 CrossRef CAS.
  26. Y. Wu, H. L. Tepper and G. A. Voth, J. Chem. Phys., 2006, 124, 024503 CrossRef PubMed.
  27. L. Turi, M.-P. Gaigeot, N. Levy and D. Borgis, J. Chem. Phys., 2001, 114, 7805 CrossRef CAS.
  28. L. Turi and D. Borgis, J. Chem. Phys., 2002, 117, 6186 CrossRef CAS.
  29. G. Bussi and M. Parrinello, Phys. Rev. E: Stat., Nonlin., Soft Matter Phys., 2007, 75, 056707 CrossRef PubMed.
  30. D. K. Park, B. W. Yoon and S. H. Lee, Bull. Korean Chem. Soc., 2015, 36, 492 CAS.
  31. A. Wallqvist and B. J. Berne, Chem. Phys. Lett., 1985, 117, 214 CrossRef CAS.
  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 Search PubMed.
  33. D. M. Bartels, J. Chem. Phys., 2001, 115, 4404 CrossRef CAS.
  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 CrossRef CAS.
  35. J. R.-R. Verlet, A. E. Bragg, A. Kammrath, O. Cheshnovsky and D. M. Neumark, Science, 2005, 307, 93 LP CrossRef PubMed.
  36. L. Ma, K. Majer, F. Chirot and B. von Issendorff, J. Chem. Phys., 2009, 131, 144303 CrossRef PubMed.
  37. G. Makov and A. Nitzan, J. Phys. Chem., 1994, 98, 3459 CrossRef CAS.
  38. L. Turi, W.-S. Sheu and P. J. Rossky, Science, 2005, 309, 914 LP CrossRef PubMed.
  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 CrossRef CAS PubMed.
  40. Y. Tang, H. Shen, K. Sekiguchi, N. Kurahashi, Y.-I. Mizuno, T. Suzuki and T. Suzuki, Phys. Chem. Chem. Phys., 2010, 12, 3653 RSC.
  41. P. Ayotte and M. A. Johnson, J. Chem. Phys., 1997, 106, 811 CrossRef CAS.
  42. F.-Y. Jou and G. R. Freeman, J. Phys. Chem., 1977, 81, 909 CrossRef CAS.
  43. F.-Y. Jou and G. R. Freeman, J. Phys. Chem., 1979, 83, 2383 CrossRef CAS.
  44. D. M. Neumark, Mol. Phys., 2008, 106, 2183 CrossRef CAS.
  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 CrossRef CAS PubMed.
  46. T. Sommerfeld and K. D. Jordan, J. Am. Chem. Soc., 2006, 128, 5828 CrossRef CAS PubMed.
  47. T. Frigato, J. VandeVondele, B. Schmidt, P. Schütte and C. Jungwirth, J. Phys. Chem. A, 2008, 112, 6125 CrossRef CAS PubMed.
  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 Search PubMed.
  49. M. E. Tuckerman, Statistical Mechanics: Theory and Molecular Simulation, Oxford University Press, 2010 Search PubMed.
  50. M. F. Herman, E. J. Bruskin and B. J. Berne, J. Chem. Phys., 1982, 76, 5150 CrossRef CAS.
  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 CrossRef CAS.
  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 CrossRef.
  55. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Oxford Scholarship Online, 2nd edn, 2017 Search PubMed.
  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 Search PubMed.
  57. M. E. Tuckerman, J. Phys. Condens. Matter, 2002, 14, R1297 CrossRef CAS.
  58. D. Marx and M. Parrinello, Z. Phys. B: Condens. Matter, 1994, 95, 143 CrossRef CAS.
  59. D. Marx and M. Parrinello, J. Chem. Phys., 1996, 104, 4077 CrossRef CAS.
  60. M. E. Tuckerman, D. Marx, M. L. Klein and M. Parrinello, J. Chem. Phys., 1996, 104, 5579 CrossRef CAS.
  61. D. Marx, M. E. Tuckerman and G. J. Martyna, Comput. Phys. Commun., 1999, 118, 166 CrossRef CAS.
  62. E. Paquet and H. L. Viktor, Adv. Chem., 2018, 2018, 9839641 Search PubMed.
  63. R. P. Feynman, Rev. Mod. Phys., 1948, 20, 367 CrossRef.
  64. R. P. Feynman and A. R. Hibbs, Quantum Mechanics and Path Integrals, McGraw-Hill, New York, 1965 Search PubMed.
  65. R. P. Feynman, Statistical Mechanics, Addison-Wesley, New York, 1972 Search PubMed.
  66. C. P. Herrero and R. Ramírez, J. Phys.: Condens. Matter, 2014, 26, 233201 CrossRef CAS PubMed.

This journal is © the Owner Societies 2022
Click here to see how this site uses Cookies. View our privacy policy here.