††thanks: email<EMAIL_ADDRESS>
# Quantum simulations of neutral water clusters and singly-charged water
cluster anions

A. Gijón  Instituto de Ciencia de Materiales de Madrid (ICMM–CSIC), Campus de
Cantoblanco, 28049 Madrid, Spain  E. R. Hernández  Instituto de Ciencia de
Materiales de Madrid (ICMM–CSIC), Campus de Cantoblanco, 28049 Madrid, Spain

###### 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 Hart and Boag
(1962), 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 Young and Neumark
(2012); Turi and Rossky (2012); Herbert and Coons (2017); Herbert (2019). 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.
Herbert, 2019 and references therein). The cavity structural model has been
challenged by a non-cavity one, put forward by Larsen et al. Larsen  _et al._
(2010). According to this model the solvated electron induces an enhanced
water-density region of $\sim$1 nm diameter over which the electron is spread.
Both cavity and non-cavity interpretations have been criticised and defended
Turi and Madarász (2011); Jacobson and Herbert (2011); Larsen  _et al._
(2011); Casey  _et al._ (2013). 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 Kuharski and Rossky (1985); Paesani  _et al._ (2006);
Habershon  _et al._ (2009); González  _et al._ (2010); Vega and Abascal
(2011); Richardson  _et al._ (2016); Ceriotti  _et al._ (2016), 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 Thirumalai  _et al._ (1986); Wallqvist  _et al._
(1987), and Barnett et al. Barnett  _et al._ (1987, 1988), 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. Takayanagi
_et al._ (2009) to perform simulations of $(\mbox{H}_{2}\mbox{O})_{5}^{-}$ and
$(\mbox{D}_{2}\mbox{O})_{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 Sec. 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
Sec. III, focusing first on neutral water clusters, and then on the singly-
charged water cluster anions. Finally, we draw our conclusions in Sec. 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\rightarrow\infty$ one can sample the Boltzmann
statistics of the quantum system through its classical isomorph, while in the
opposite limit, $P\rightarrow 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 Toukan and Rahman (1985); Lobaugh and Voth (1997). 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_{\text{O}}=-0.82e=-2q_{\text{H}}$, where $e$
is the quantum of charge), and a Lennard-Jones potential acting between oxygen
atoms.

In cluster anions the excess electron interaction with the oxygen and hydrogen
atoms of water molecules is described via the Turi-Borgis Turi  _et al._
(2001); Turi and Borgis (2002) 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:

$\begin{split}V_{e-\text{W}}({\bf
r}){}=\sum_{n=1}^{N_{\text{W}}}\Big{[}V_{e-\text{O}}\big{(}r^{(n)}_{\text{O}}\big{)}+V_{\text{pol}}\big{(}r^{(n)}_{\text{O}}\big{)}+\\\
{}V_{e-\text{H}}\big{(}r^{(n)}_{\text{H}_{1}}\big{)}+V_{e-\text{H}}\big{(}r^{(n)}_{\text{H}_{2}}\big{)}\Big{]},\end{split}$
(1)

where $r^{(n)}_{i}=|{\bf r}-{\bf r}^{(n)}_{i}|$, and ${\bf 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 Eq. (1) can be found
in refs. Turi  _et al._ , 2001; Turi and Borgis, 2002.

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 Bussi and Parrinello (2007) 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 III.1, 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 III.2 we report the results obtained
for clusters charged with an excess electron.

###  III.1 Neutral water clusters

We have computed the energetic and structural properties of neutral water
clusters $(\text{H}_{2}\text{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.

Figure 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 $(\text{H}_{2}\text{O})_{N}$ clusters.

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=0K$, 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.

Figure 2: Averages of the intra-molecular O-H distance (a) and H-O-H angle
(b), and cluster radius (c) of neutral $(\text{H}_{2}\text{O})_{N}$ clusters,
as a function of the temperature, for quantum (solid line and filled symbols)
and classical (dashed line and empty symbols) cases.

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

$\Delta E_{n}=E[(\text{H}_{2}\text{O})_{n}]-nE(\text{H}_{2}\text{O})\,,$ (2)

where $E(\text{H}_{2}\text{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 $(\text{H}_{2}\text{O})_{n}$
clusters. Fig. (2-a) 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. (2-c) 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. 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 $\sim
4.3\times 10^{-9}\mbox{m}^{2}/\mbox{s}$, a value comparable to that obtained
from simulations of super-cooled water at the same temperature Park  _et al._
(2015).

Figure 3: Radial pair distribution function of the oxygen-oxygen
($g_{\text{OO}}$) distance for quantum (solid line) and classical (dashed
line) simulations of the neutral $(\text{H}_{2}\text{O})_{237}$ cluster.

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_{\text{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_{\text{OO}}$ function at 100 and 200 K [see panel
(a)] at $r\approx 2.75$ Å and the second at $r\approx 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. (2-c), 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 Kuharski and Rossky (1985); Wallqvist and Berne (1985), whereas
the reverse is true for temperatures above 200 K.

###  III.2 Water cluster anions

Now we turn to consider the properties of singly-charged water cluster anions,
$(\text{H}_{2}\text{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.

Figure 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 $\sim$ 90 % of the charge density, the intermediate one $\sim$
50 %, and the inner one $\sim$ 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.

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 $(\mbox{H}_{2}\mbox{O})^{-}_{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 $(\mbox{H}_{2}\mbox{O})^{-}_{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\leq 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
$(\mbox{H}_{2}\mbox{O})^{-}_{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,
$d_{e}=\sqrt{\langle|\bar{\bm{r}}|^{2}\rangle}$, with $\bar{\bm{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:

$\displaystyle\begin{split}d_{e}+R_{e}<R_{c}&\implies\text{Interior
state}\,,\\\ d_{e}+R_{e}\gtrsim R_{c}&\implies\text{Surface
state}\,,\end{split}$ (3)

where $R_{\text{c}}$ is the gyration radius of the cluster and $R_{e}$ is that
for the electron [see Eq. (12)]. 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.

Figure 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. Turi and Borgis, 2002. Panels (c) and (d) show the same
quantities as (a) and (b) for the classical case.

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 Eq. (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. (5-c), small clusters
($n\leq 47$) can only accommodate an exterior state at all temperatures, while
bigger clusters ($n\geq 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. (5-a), 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 Eq. (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.

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
Turi and Borgis, 2002, with the same model used here. This value agrees well
with the experimental moment analysis of the absorption spectrum Lavrich  _et
al._ (1994); Bartels (2001), 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.

Figure 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.

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. (6-a). 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
$\bm{r}_{\text{w-}e}=\bar{\bm{r}}_{e}-\bm{r}_{\text{w}}$ at short distances.
This is the case for both internal and external electronic states.

Figure 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.

The orientation of the molecular dipoles in the outer layer of the neutral
clusters, see Fig. (6-a), 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 $(\text{H}_{2}\text{O})_{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)].

Figure 8: Classical (a) and quantum (b) average of the simulated ground-state
energy of the excess electron for $(\text{H}_{2}\text{O})_{n}^{-}$ clusters at
different temperatures. The experimental fit of Coe from Ref. Coe  _et al._ ,
1990 is shown with dashed black lines, while the linear fits from our
interior-state data are represented with dotted black lines.

Let us now focus on the energetic properties of the excess electron. Its
average ground state energy, $\varepsilon_{0}$, can be compared with the
negative of the experimental vertical detachment energy (VDE) as obtained via
photoelectron spectroscopy Coe  _et al._ (1990); Verlet  _et al._ (2005); Ma
_et al._ (2009). 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 Makov and Nitzan (1994). This was indeed observed in the
experimental results of Coe et al. Coe  _et al._ (1990), whose linear fit is
shown in Fig. (8) as a dashed black line. Our average values for
$\varepsilon_{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. Turi  _et al._ (2005). 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. Turi  _et al._ (2005), which is closer to the
experimental value of -3.3 eV Siefermann  _et al._ (2010); Tang  _et al._
(2010). 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.

Figure 9: Classical (a) and quantum (b) average of the simulated energy gap of
the excess electron for $(\text{H}_{2}\text{O})_{n}^{-}$ clusters at different
temperatures. The experimental fit of Ayotte from Ref. Ayotte and Johnson,
1997 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,
$\varepsilon_{\text{gap}}=\varepsilon_{1}-\varepsilon_{0}$, are shown together
with the experimental fit obtained from the maxima of the optical adsorption
spectra by Ayotte et al. Ayotte and Johnson (1997). 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 Jou and R. Freeman
(1977, 1979).

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 Neumark
(2008), who carefully analyzed the experimental VDE data previously assigned
to internal states Verlet  _et al._ (2005), concluding that the transition
occurs in the size range $n=25-35$. Furthermore, this assessment is in full
agreement with other low temperature experiments Ma  _et al._ (2009);
Herburger  _et al._ (2019) and consistent with theoretical works Sommerfeld
and D. Jordan (2006); Frigato  _et al._ (2008), 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 cluster in the presence of an excess
electron, apart from employing a suitable water-electron interaction.

##  V Acknowledgments

We thank J. M. Soler, M. V. Fernández-Serra, R. Ramírez and C. P. Herrero for
helpful discussions. This work has been supported by the Spanish State Agency
for scientific research through project PGC2018-096955-B-C44.

##  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. 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:

$Q_{P}=\left[\prod_{i=1}^{N}\left(\frac{m_{i}P}{2\pi\beta\hbar^{2}}\right)^{3P/2}\int
d\bm{r}_{i}^{1}\cdots
d\bm{r}_{i}^{P}\right]\times\exp{\left\\{-\beta\sum_{\alpha=1}^{P}\left[\sum_{i=1}^{N}\frac{m_{i}P}{2\beta^{2}\hbar^{2}}(\bm{r}_{i}^{\alpha}-\bm{r}_{i}^{\alpha+1})^{2}+\frac{V(\bm{r}^{\alpha})}{P}\right]\right\\}}_{\bm{r}_{i}^{P+1}=\bm{r}_{i}^{1}}.$
(4)

Here, $Q_{P}$ is the partition function of the P-bead classical isomorph
approximating that of its quantum counterpart; $\beta=(k_{B}T)^{-1}$, and
$V(\bm{r}^{\alpha})$ 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.
$\bm{r}^{\alpha}\equiv(\bm{r}_{1}^{\alpha},\bm{r}_{2}^{\alpha}\ldots\bm{r}_{n}^{\alpha})$.
The isomorphism is exact in the limit $P\rightarrow\infty$. 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 Gillan (1990); Tuckerman (2010); Herman  _et al._ (1982) from
the partition function Eq. (4). 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(\bm{r})$ in Eq. (4) reduces to
$V_{\text{SPC/F}}(\bm{r})$, i.e. the SPC/F potential energy expression.
However, for the 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.

### 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 Thirumalai  _et al._ (1986); Wallqvist  _et
al._ (1987), 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, Eq. (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\times
64\times 64$ grid points, which provides sufficient resolution, as we discuss
below in the convergence tests subsection. The lowest 4 electron eigen-pairs
are then obtained iteratively using the imaginary time propagation method Chin
_et al._ (2009). In this method, a set of n trial wave functions $\psi_{i}$ is
driven towards the n lowest lying eigenstates of Hamiltonian $H$ by acting on
them with the imaginary time evolution operator:

$\psi_{i}(\epsilon)=e^{-\epsilon H}\psi_{i}(0).$ (5)

Here $\epsilon$ 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^{-\epsilon H}$ on the trial set cannot be calculated
exactly, but it can be approximated up to any desired order in $\epsilon$ Chin
_et al._ (2009) using a Trotter-like factorisation:

$e^{-\epsilon H}\approx\sum_{k}c_{k}\prod_{j}e^{-a_{jk}T}e^{-b_{jk}V},$ (6)

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^{\alpha T}$ are diagonal in momentum space, while
operators $e^{\beta V}$ are diagonal in real space (provided $V$ is local, as
is the case in our model), so the action of Eq. (6) on the trial set involves
a succession of Fast Fourier transformations Frigo and Johnson (1997) between
momentum and real space (see Ref. Chin  _et al._ , 2009 for details). In
practice we find that a $4^{th}$-order approximation to $e^{-\epsilon H}$
provides a good compromise between computational effort and speed of
convergence.

In the case of the cluster anions $V(\bm{r})$ in Eq. (4) takes the form:

$V(\bm{r}^{\alpha})=E_{e}(\bm{r}^{\alpha})+V_{\text{SPC/F}}(\bm{r}^{\alpha}),$
(7)

where $E_{e}(\bm{r}^{\alpha})$ is the electron ground state energy for bead
$\alpha$, which depends parametrically on $\bm{r}^{\alpha}$, and
$V_{\text{SPC/F}}(\bm{r}^{\alpha})$ is the SPC/F model potential for water, as
in the neutral clusters.

### Centroid approximation

Strictly speaking, the ground-state energy, $E_{e}(\bm{r}^{\alpha})$, must be
evaluated separately for each replica $\alpha$, 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:

$\displaystyle\bar{\bm{r}}_{i}$
$\displaystyle=\frac{1}{P}\sum_{\alpha=1}^{P}\bm{r}_{i}^{\alpha}\,,$ (8)
$\displaystyle r_{g,i}^{2}$
$\displaystyle=\frac{1}{P}\sum_{\alpha=1}^{P}(\bm{r}_{i}^{\alpha}-\bar{\bm{r}}_{i})^{2}\,.$
(9)

The centroid approximation amounts to solving the electron Schrödinger
equation for the ring polymer centroid instead of doing it for each individual
bead.

Figure 10: (Left) Schematic representation of the beads, $\bm{r}^{\alpha}$,
centroid position, $\bar{\bm{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.

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 $\alpha$ and write it in terms of that evaluated for the
centroid; we would have:

$E_{e}(\bm{r}^{\alpha})=E_{e}(\bar{\bm{r}})+\nabla
E_{e}(\bar{\bm{r}})\cdot(\bm{r}^{\alpha}-\bar{\bm{r}})+\mathcal{O}(|\bm{r}^{\alpha}-\bar{\bm{r}}|^{2})\,.$
(10)

The centroid approximation implies taking $E_{e}(\bm{r}^{\alpha})\approx
E_{e}(\bar{\bm{r}})$, which will be justified when the second and subsequent
terms on the rhs of Eq. (10) are negligible compared to $E_{e}(\bar{\bm{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 $i$th atom in replica $\alpha$ due to the
electron cloud is given by

$\bm{F}_{e,i}^{\alpha}=\frac{1}{P}\bar{\bm{F}}_{e,i}\,,$ (11)

as can be seen using the chain rule, where $\bar{\bm{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 Eq. (4) 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 Bussi and Parrinello (2007)
thermostat attached to each degree of freedom. We used a time step of 1 fs and
a friction parameter $\gamma=0.001$ au, 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 $\text{g/cm}^{3}$, including all water molecules within spheres of
varying radii Rudberg  _et al._ (2018). For each calculation, an equilibration
run of $1\times 10^{5}$ steps (100 ps) was executed before the production
simulation of $5\times 10^{5}$ 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,
$V_{\text{conf}}=\frac{1}{2}k(x^{8}+y^{8}+z^{8})$, with $k=10^{-8}$ au. The
relaxation was performed at each temperature during $10^{4}$ steps, and after
that the confining term was switched off. Typical calculations consisted of an
equilibration run of $5\times 10^{4}$ steps (50 ps), previous to the
production run with length between $10^{5}$ and $2\times 10^{5}$ steps, which
correspond to 100 and 200 ps of time. 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 subsection.

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 Allen and Tildesley
(2017); Grossfield  _et al._ (2018).

### Convergence tests

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

Surface state, $N$=100, $T=300\,$K

| $N_{p}=32$ | $N_{p}=64$ | $N_{p}=128$
---|---|---|---
| $\varepsilon_{0}$ | $\varepsilon_{1}$ | $\varepsilon_{0}$ | $\varepsilon_{1}$ | $\varepsilon_{0}$ | $\varepsilon_{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

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 Å.

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=12800$ K value, which implies using 128 replicas at $T=100$ K, 64 at
$T=200$ K, and 48 at $T=300$ K.

Figure 12: Relative error made in the averages when using the centroid
approximation to study the $(\text{H}_{2}\text{O})_{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, $\varepsilon_{0}$, and
its spatial width, $R_{e}$.

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}\,\text{Hartree}\approx 2.7\times 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 $N_{p}=64$ for all sizes and
temperatures.

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\times 10^{4}$ 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

$R_{e}=\left[\int|\bm{r}-\bar{\bm{r}}|^{2}\rho(\bm{r})d\bm{r}\right]^{1/2}\,,$
(12)

where $\rho$ 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 ($\propto T^{2}$), 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$ %).

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_{\text{B}}T$ 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.

## References

  * Hart and Boag (1962) E. J. Hart and J. W. Boag, J. Am. Chem. Soc. 81, 4090 (1962).
  * Young and Neumark (2012) R. M. Young and D. M. Neumark, Chemical Reviews 112, 5553 (2012).
  * Turi and Rossky (2012) L. Turi and P. J. Rossky, Chemical Reviews 112, 5641 (2012).
  * Herbert and Coons (2017) J. M. Herbert and M. P. Coons, Annual Reviews of Physical Chemistry 68, 447 (2017).
  * Herbert (2019) J. M. Herbert, Physical Chemistry and Chemical Physics 21, 20538 (2019).
  * Larsen  _et al._ (2010) R. E. Larsen, W. J. Glover, and B. J. Schwartz, Science 329, 65 (2010).
  * Turi and Madarász (2011) L. Turi and A. Madarász, Science 331, 1387 (2011).
  * Jacobson and Herbert (2011) L. D. Jacobson and J. M. Herbert, Science 331, 1387 (2011).
  * Larsen  _et al._ (2011) R. E. Larsen, W. J. Glover, and B. J. Schwartz, Science 331, 1387 (2011).
  * Casey  _et al._ (2013) J. R. Casey, A. Kahros, and B. J. Schwartz, Journal of Physical Chemistry B 117, 14173 (2013).
  * Kuharski and Rossky (1985) R. A. Kuharski and P. J. Rossky, The Journal of Chemical Physics 82, 5164 (1985).
  * Paesani  _et al._ (2006) F. Paesani, W. Zhang, D. A. Case, T. E. Cheatham, and G. A. Voth, The Journal of Chemical Physics 125, 184507 (2006), https://doi.org/10.1063/1.2386157 .
  * Habershon  _et al._ (2009) S. Habershon, T. E. Markland, and D. E. Manolopoulos, The Journal of Chemical Physics 131, 024501 (2009), https://doi.org/10.1063/1.3167790 .
  * González  _et al._ (2010) B. S. González, E. G. Noya, C. Vega, and L. M. Sesé, The Journal of Physical Chemistry B 114, 2484 (2010), pMID: 20121175, https://doi.org/10.1021/jp910770y .
  * Vega and Abascal (2011) C. Vega and J. L. F. Abascal, Phys. Chem. Chem. Phys. 13, 19663 (2011).
  * Richardson  _et al._ (2016) 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 351, 1310 (2016).
  * Ceriotti  _et al._ (2016) M. Ceriotti, W. Fang, P. G. Kusalik, R. H. McKenzie, A. Michaelides, M. A. Morales, and T. E. Markland, Chemical Reviews 116, 7529 (2016).
  * Thirumalai  _et al._ (1986) D. Thirumalai, A. Wallqvist, and B. J. Berne, The Journal of Statistical Physics 43, 973 (1986).
  * Wallqvist  _et al._ (1987) A. Wallqvist, D. Thirumalai, and B. J. Berne, The Journal of Chemical Physics 86, 6404 (1987).
  * Barnett  _et al._ (1987) R. N. Barnett, U. Landman, C. L. Cleveland, and J. Jortner, Physical Review Letters 59, 811 (1987).
  * Barnett  _et al._ (1988) R. N. Barnett, U. Landman, C. L. Cleveland, and J. Jortner, Journal of Chemical Physics 88, 65 (1988).
  * Takayanagi  _et al._ (2009) T. Takayanagi, T. Yoshikawa, H. Motegi, and M. Shiga, Chemical Physics Letters 482, 195 (2009).
  * Toukan and Rahman (1985) K. Toukan and A. Rahman, Physical Review B (1985), 10.1103/PhysRevB.31.2643.
  * Lobaugh and Voth (1997) J. Lobaugh and G. A. Voth, Journal of Chemical Physics (1997), 10.1063/1.473151.
  * Turi  _et al._ (2001) L. Turi, M.-P. Gaigeot, . N. Levy, and D. Borgis, The Journal of Chemical Physics 114, 7805 (2001).
  * Turi and Borgis (2002) L. Turi and D. Borgis, The Journal of Chemical Physics 117, 6186 (2002).
  * Bussi and Parrinello (2007) G. Bussi and M. Parrinello, Physical Review E - Statistical, Nonlinear, and Soft Matter Physics 75 (2007), 10.1103/PhysRevE.75.056707.
  * Park  _et al._ (2015) D. K. Park, B. W. Yoon, and S. H. Lee, Bulletin of the Korean Chemical Society 36, 492 (2015).
  * Wallqvist and Berne (1985) A. Wallqvist and B. J. Berne, Chemical Physics Letters 117, 214 (1985).
  * Lavrich  _et al._ (1994) D. J. Lavrich, P. J. Campagnola, , and M. A. Johnson (Plenum, New York, 1994).
  * Bartels (2001) D. M. Bartels, The Journal of Chemical Physics 115, 4404 (2001).
  * Coe  _et al._ (1990) 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, The Journal of Chemical Physics 92, 3980 (1990).
  * Verlet  _et al._ (2005) J. R. R. Verlet, A. E. Bragg, A. Kammrath, O. Cheshnovsky, and D. M. Neumark, Science 307, 93 LP (2005).
  * Ma  _et al._ (2009) L. Ma, K. Majer, F. Chirot, and B. von Issendorff, The Journal of Chemical Physics 131, 144303 (2009).
  * Makov and Nitzan (1994) G. Makov and A. Nitzan, The Journal of Physical Chemistry 98, 3459 (1994).
  * Turi  _et al._ (2005) L. Turi, W.-S. Sheu, and P. J. Rossky, Science 309, 914 LP (2005).
  * Siefermann  _et al._ (2010) K. R. Siefermann, Y. Liu, E. Lugovoy, O. Link, M. e. Faubel, U. Buck, B. Winter, and B. Abel, Nature chemistry 2, 274 (2010).
  * Tang  _et al._ (2010) Y. Tang, H. Shen, K. Sekiguchi, N. Kurahashi, Y.-I. Mizuno, Tomoya an d Suzuki, and T. Suzuki, Physical Chemistry Chemical Physics 12, 3653 (2010).
  * Ayotte and Johnson (1997) P. Ayotte and M. A. Johnson, The Journal of Chemical Physics 106, 811 (1997).
  * Jou and R. Freeman (1977) F.-Y. Jou and G. R. Freeman, The Journal of Physical Chemistry 81, 909 (1977).
  * Jou and R. Freeman (1979) F.-Y. Jou and G. R. Freeman, The Journal of Physical Chemistry 83, 2383 (1979).
  * Neumark (2008) D. M. Neumark, Molecular Physics 106, 2183 (2008).
  * Herburger  _et al._ (2019) A. Herburger, E. Barwa, M. Ončák, J. Heller, C. van dẽr Linde, D. M. Neumark, and M. K. Beyer, Journal of the American Chemical Society 141, 18000 (2019).
  * Sommerfeld and D. Jordan (2006) T. Sommerfeld and K. D. Jordan, Journal of the American Chemical Society 128, 5828 (2006).
  * Frigato  _et al._ (2008) T. Frigato, J. VandeVondele, B. Schmidt, and P. Schütte, Christof a nd Jungwirth, The Journal of Physical Chemistry A 112, 6125 (2008).
  * Gillan (1990) M. J. Gillan, in  _Computer Modelling of Fluids Polymers and Solids_, edited by C. R. A. Catlow, S. C. Parker, and M. P. Allen (Springer Netherlands, Dordrecht, 1990) pp. 155–188.
  * Tuckerman (2010) M. E. Tuckerman,  _Statistical Mechanics: Theory and Molecular Simulation_ (Oxford University Press, 2010).
  * Herman  _et al._ (1982) M. F. Herman, E. J. Bruskin, and B. J. Berne, The Journal of Chemical Physics 76, 5150 (1982).
  * Chin  _et al._ (2009) S. A. Chin, S. Janecek, and E. Krotscheck, Chemical Physics Letters 470, 342 (2009).
  * Frigo and Johnson (1997) M. Frigo and S. G. Johnson,  _The Fastest Fourier Transform in the West_ , Tech. Rep. MIT-LCS-TR-728 (Massachusetts Institute of Technology, 1997).
  * Rudberg  _et al._ (2018) E. Rudberg, E. H. Rubensson, P. Sałek, and A. Kruchinina, SoftwareX 7, 107 (2018).
  * Allen and Tildesley (2017) M. P. Allen and D. J. Tildesley,  _Computer Simulation of Liquids: Second Edition_ (Oxford Scholarship Online, 2017).
  * Grossfield  _et al._ (2018) A. Grossfield, P. N. Patrone, D. R. Roe, A. J. Schultz, D. Siderius, and D. M. Zuckerman, Living Journal of Computational Molecular Science 1, 5067 (2018).
  * Tuckerman (2002) M. E. Tuckerman, Journal of Physics Condensed Matter 14, R1297 (2002).
  * Marx and Parrinello (1994) D. Marx and M. Parrinello, Zeitschrift für Physik B Condensed Matter 95, 143 (1994).
  * Marx and Parrinello (1996) D. Marx and M. Parrinello, The Journal of Chemical Physics 104, 4077 (1996).
  * Tuckerman  _et al._ (1996) M. E. Tuckerman, D. Marx, M. L. Klein, and M. Parrinello, The Journal of Chemical Physics 104, 5579 (1996).
  * Marx  _et al._ (1999) D. Marx, M. E. Tuckerman, and G. J. Martyna, Computer Physics Communications 118, 166 (1999).
  * Paquet and Viktor (2018) E. Paquet and H. L. Viktor, Advances in Chemistry 2018 (2018), https://doi.org/10.1155/2018/9839641.
  * Feynman (1948) R. P. Feynman, Reviews of Modern Physics 20, 367 (1948).
  * Richard P Feynman; Albert R Hibbs (1965) Richard P Feynman; Albert R Hibbs,  _Quantum Mechanics and Path Integrals_ (New York : McGraw-Hill, ©1965, 1965).
  * Feynman (1972) R. P. Feynman,  _Statistical Mechanics_ (Addison-Wesley, New York, 1972).
  * Herrero and Ramírez (2014) C. P. Herrero and R. Ramírez, Journal of Physics: Condensed Matter 26, 233201 (2014).