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
First published on 26th May 2022
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.
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.
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:
![]() | (1) |
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.
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 = 47100
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.
The formation or binding energy of a neutral cluster is defined as
ΔEn = E[(H2O)n] − nE(H2O), | (2) |
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.
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.
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.
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).
Throughout each simulation, we compute the average distance between the electron mean position and the cluster's center of mass, , with
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:
![]() | (3) |
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.
![]() | ||
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 = e − rw at short distances. This is the case for both internal and external electronic states.
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].
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.
![]() | ||
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
![]() | ||
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.
![]() | (A1) |
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.
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) |
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:
![]() | (A3) |
In the case of the cluster anions V(r) in eqn (A1) takes the form:
V(rα) = Ee(rα) + VSPC/F(rα), | (A4) |
![]() | (A5) |
![]() | (A6) |
A schematic representation of a quantum particle with 8 beads around its centroid position can be seen in left panel of Fig. 10.
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(![]() ![]() ![]() ![]() ![]() | (A7) |
![]() | (A8) |
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.
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, , 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
![]() | ||
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.
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
![]() | (A9) |
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.
This journal is © the Owner Societies 2022 |