Many-Body Contributions in Water Nanoclusters

Many-body interactions in water are known to be important but difficult to treat in atomistic models and often are included only as a correction. Polarizable models treat them explicitly, with long-range many-body potentials, within their classical approximation. However, their calculation is computationally expensive. Here, we evaluate how relevant the contributions to the many-body interaction associated with different coordination shells are. We calculate the global energy minimum, and the corresponding configuration, for nanoclusters of up to 20 water molecules. We find that including the first coordination shell, i.e., the five-body term of the central molecule, is enough to approximate within 5% the global energy minimum and its structure. We show that this result is valid for three different polarizable models, the Dang–Chang, the MB-pol, and the Kozack–Jordan potentials. This result suggests a strategy to develop many-body potentials for water that are reliable and, at the same time, computationally efficient.


Introduction
Water is an object of intense research for its unusual properties and central role in many areas of science and technology. 1 In the last 50 years, computer simulations have contributed to understanding some of these peculiar phenomena.On one hand, ab-initio calculations have been used to predict structural [2][3][4][5] and dynamical properties 6,7,[7][8][9] of water from quantum calculations.However, this technique requires a high computational cost and generally treats small systems of the order of ≈ 100 molecules.[12][13][14][15] However, in classical Molecular Dynamics or Monte Carlo simulations, the choice of the force field is critical.Many of the classical force fields for water are based on pairwise dispersion-repulsion and electrostatic interactions, and the many-body contributions are neglected.One of the most popular potential models is the TIP4P 16 and its family of TIP4P-like models. 17,18In these models, each water molecule is considered rigid, with four sites, including oxygen, two hydrogens, and one, often called the M site, located along the bisector of the oxygen-hydrogen vectors.The water-water interaction is given by Lennard-Jones and Coulomb pairwise potentials.The parameters of the TIP4P potential are chosen to replicate the structural properties of bulk water at standard temperature and pressure.
However, this nonpolarizable potential cannot reproduce, for example, high-density properties where the many-body interactions play a fundamental role. 19Thus, the polarizable models are built to overcome these deficiencies and are based on explicitly incorporating a non-additive term.
In this work, we aim to elucidate how relevant are the many-body effects on the energetic and structural properties of water nano-clusters.To achieve this goal, we first employ the rigid-body polarizable Dang-Chang (DC) potential. 20This model, defined in the Methods section, is characterized by two terms.One is associated with the pairwise additive and the other with the non-additive polarization term.
To describe the importance of the many-body effects in water, we introduce a cut-off radius in the polarization term and, employing the Basin-Hopping global optimization technique, 21 we identify the putative global energy minimum for several selected water clusters with nm-size.We expect to obtain global minimum structures similar to the known TIP4P configurations for a minimal value of the cut-off radius, whereas, for a larger cut-off, the DC minimum structures.We ask if we can find an intermediate cut-off value that could reproduce the DC results within a reasonable approximation.

Energy dependence
Based on previous work, 22 we focus on water nano-clusters with 6 to 20 molecules (Fig. 1).
We minimize the energy of each cluster, as described in the Methods section, by applying a cut-off exclusively to the non-additive polarization term of the DC potential.
Figure 1: The lowest-energy configuration for a cluster with 6 (a,f), 8 (b, g), 10 (c,h), 16 (d, i), and 20 (e, j) water molecules calculated for the model with the shortest cut-off 1 Å, as an approximation (see text) of the TIP4P-like model (a -e), and with the largest cut-off 20 Å, as in the DC limit (f -j).For N = 16, our method recovers the same water molecule orientations in the two limits.Surprisingly, we observe that the cut-off also influences the Lennard-Jones and the Coulomb energy contributions as an indirect effect due to the structural changes of the global-energy minima induced by the polarization.
We find that, for all cluster sizes, the resulting binding energy for the minimum-energy configurations is non-monotonic as a function of the cut-off radius r ( Fig. 1 and Fig. 2 in the Supporting Information).
When r is larger than d max , the largest O-O distance in the cluster, all the energy contributions must converge to a constant value.For the octamer, for example, this is true at r > 5 Å because it is d max 4.85 Å.Thus, we recover the energy calculated for the DC model 22 and the corresponding configurations (Fig. 1 f-j).
When r is shorter than the water first-coordination shell distance, r 3 Å, we expect that the polarization energy vanishes and the other terms are constant.Under these circumstances, although the DC parameters for the isotropic potentials are slightly different from those of the TIP4P, one could expect that the DC model approximates well the TIP4P-like's minimum energy configurations (Fig. 1 a-e).We verify it is so for all the cases we considered except for N = 20.However, there is no apparent change between the results for the two extreme cut-off radii for 16 molecules (Fig. 1 d-i).We will discuss the surprising result for 16 molecules in a separate section.
Therefore, by increasing the cut-off radius, we tune the global-minimum energy from the unpolarizable to the polarizable model.In particular, for all the cluster sizes, we observe a switching behavior in the binding energy when we cross the r 3 Å threshold as a consequence of the sudden change in the number of water molecules interacting via the polarization potential.To illustrate this point, we calculate the average number N i of interacting water molecules as a function of r in each cluster (Fig. 3).The average is over all the molecules of the same clusters.Above r 3 Å, the number N i jumps from 1 to 4 or 5 for the three smaller and the two larger clusters, respectively (in bulk water, the number of molecules within the first shell would be 5 in a tetrahedral configuration and 6 in a local high-density configuration with an interstitial molecule).
By further increasing r, N i has step-like increases at each coordination-shell distance for the different clusters.However, the steps smooth out for larger clusters due to the broadening of the O-O distances distribution over which we calculate N i .As a consequence of the variation of N i , the calculated binding energy has a non-monotonic behavior with r, e.g., with a minimum at r = 4.0 Å, for the octamer (Fig. 2) or maxima for the other clusters (Fig. 1 in the Supporting Information), due to partial contributions of the coordination shells, as we discuss next.

Minimum-energy configurations
The visual comparison of our minimum-energy configurations at the two extreme cut-off radii, r 3 Å and r > d max , with the configurations of the global-energy minima found in the literature for TIP4P 23 and DC, 22 confirms that, by tuning r, we move between the space of minima of the two limiting models.For example, for the hexamer (6 molecules), we modulate between the cage structure of the TIP4P 24 for r < 3Å (Fig. 1a) to the trigonal prism of the DC 22 for r > d max = 4.16 Å (Fig. 1d).
To make this comparison more quantitative, we calculate how close the configuration at a given cut-off r (the cut-off cluster) is to the DC reference cluster with the long-range many-body contributions by computing the Root-mean-square deviation (RMSD) between the two configurations where R i and r i are the O-O distances within the reference and the cut-off cluster, respectively.The index i labels the different distances between molecules in the cluster, and N is the size of the cluster.We evaluate the RMSD over all the possible permutations of distances within the clusters and consider the minimum as our estimate (Fig. 4).
As expected, we recover the minimum-energy configuration of the full DC case for large cut-off values.Interestingly, we find that the dependence of the RMSD is not monotonic with the cut-off, displaying several minima and maxima.In particular, all the clusters have a minimum RMSD < 0.05 Å at r 3 Å, i.e., the distance of the first coordination shell.Thus, the minimum energy configuration becomes very similar to the DC limit when we include the first coordination shell (see Fig. 3 in the Supporting Information to observe the minimum energy configurations for r ∼ 3 Å).
At cut-off radii that do not correspond to the distance of a coordination shell, the RMSD increases, i.e., the agreement between the configuration of minimum energy and the DC reference case is reduced.This happens although the total energy of the cluster reaches a global minimum, as for the octamer at r = 4.0 Å (Fig. 2).Indeed, this minimum is a consequence of the many-body interaction acting on a number of molecules N i that is intermediate between two consecutive coordination shells (Fig. 4 b).
Next, we study the total energy deviation from the DC reference results as a function of the cut-off radius r.We find that, for all cluster sizes, the energy deviations at the first coordination shell ( 3 Å) are 5% than the reference DC energy (Fig. 2 in Supporting Information).The drop in the energy deviation at the first coordination-shell distance is evident when we represent it as a function of N i , observing that for our clusters, the first coordination shell includes from 3 to 4.6 molecules (Fig. 5).
Furthermore, the energy deviation increases for larger cut-off radii that are intermediate between two consecutive coordination shells, and it does not improve significantly when we include more shells.Overall, we conclude that both the RMSD and the energy deviation of the cut-off cluster drop below 5% compared with the DC reference cluster when the cutoff coincides with the first coordination shell.To prove that these results are not model dependent, we perform the same analysis using the MB-pol potential, [25][26][27] which has been shown to correctly predict the properties of water across a wide range of thermodynamic conditions, 28,29 and the Kozack-Jordan (KJ) potential. 30In the Supporting Information, we show that we find the same behavior as a function of the cut-off radius, proving the generality of our result.

The nano-cluster with 16 molecules
The RMSD is an accurate observable to differentiate similar clusters.For example, for 8 and 10 molecules, the clusters with the cut-off at the first coordination shell are similar to the reference DC clusters.However, they have different ordering in the hydrogen bond directions.Consequently, the RMSD of the corresponding clusters is finite (Fig. 4 b, and c).
The accuracy of RMSD allows us to understand also the case with 16 molecules, displaying minimum-energy configurations with no polarization contribution (Fig. 1 d) and full polarization contribution (Fig. 1 i) that are indistinguishable by the naked eye.We find this surprising result also for 12 molecules (not shown) that likewise minimize their energy by clustering as fused cubes.Nevertheless, although the configurations of fused cubes look the same in both short and large cut-off limits, they are not.Indeed, the RMSD for 16 molecules at the first coordination shell is > 1% larger than the reference DC cluster (Fig. 4 d).A more refined analysis (not presented here) reveals that these differences are due to minimal variations in the O-O distances of the short-cut-off cluster that account for an energy deviation > 6% larger than the reference DC cluster (Fig. 2 in the Supporting Information).

Conclusions
We analyze the effect of introducing a cut-off r in three many-body, long-range, polarizable potentials: the DC model, which adds polarizability to a TIP4P-like potential, the MB-pol, and the KJ potential.First, we consider DC-water nano-clusters of up to 20 molecules and calculate their minimum-energy configuration.We evaluate the root-mean-square deviation of the structure and the deviation of the energy of the minima compared to the reference unrestricted simulations with the full DC polarizable potential. 22We repeat the calculations for the MB-pol and the KJ potentials.To minimize the computational cost, we 1) replicate the entire procedure only for two representative cluster sizes (8 and 20) of the KJ water, and 2) use the minimum-energy configurations of DC-water as a starting point for the MB-pol analysis.
Surprisingly, we find that the cut-off, although applied only to the polarization term of the potential and not the other interactions, induces variations in the energies of each term of the potential.This is a consequence of the change in the global minimum structure dominated by the many-body interaction.
For the DC model, we find that the deviations are not monotonic with the cut-off radius r.When r does not correspond to the average distance of a coordination shell, the nanoclusters reach artificial minima with energy below the full DC case but with a significant structural deviation from the correct minimum.When the r corresponds to the first shell, we find an agreement within 5% for the RMSD of the configuration and its total energy relative to the reference values.For the larger 1-nm clusters considered here, the first shell comprises five molecules.The same behavior is found for the MB-pol and KJ models.Furthermore, the same number of first coordinated molecules is also characteristic of the bulk water in a local tetrahedral structure.Therefore, our results show that approximating the long-range, many-body (polarization) interaction with a short-range, five-body interaction is better than using fewer-body terms, what may lead to improvements over recent potential energy functions that account for explicit short-range interactions only up to three-body contribution. 318][39][40] These approaches find support in our present results.It is likely that contribution beyond the first coordination shell would be necessary to determine the correct minimum binding energy and other properties of water.However, we have shown that, for three different polarizable water models at the cost of (less than) 5% error in interaction energies and the right global minimum structure, it is enough to consider just the contribution of the first coordination shell.This approximation would vastly reduce the computational cost of large-scale simulations of hydrated systems, including an effective approximation of the long-range many-body interactions.

Methods
Pairwise-additive Lennard-Jones and Coulomb terms plus a many-body polarization contribution describe the rigid-body polarization DC potential. 20The Lennard-Jones term is applied between oxygen atoms, whereas the Coulomb interaction is on partial charges on the hydrogen and M sites.The polarization term is characterized by the isotropic molecular polarizability on the M site and the induced dipole moments due to the electric field produced by fixed charges in the system.
The putative global energy minima of water clusters were located using the Basin-Hopping method. 212][43] To treat the rigid-body orientational degrees of freedom, we employed the angle-axis scheme. 44The advantage of this coordinate system is that it does not suffer from the "gimbal lock" problem, which can occur with the use of Euler angles. 45 perform four independent trajectories of 1×10 5 basin-hopping steps for each cluster, starting with random geometries and a constant optimization temperature of k B T 3 kcal/mol.We attempt blocks of 100 translational and 200 angular moves with an acceptance ratio of 20%.
To evaluate how the global energy minimum changes with the cut-off radius, first, we calculate the global minimum for 20 water molecules with a cut-off radius of 20.0 Å.We check that this distance is enough to account for all many-body energy contributions for the cluster sizes considered here.Under this condition, we recover the global minimum obtained with the Dang-Chang Model. 22This preliminary analysis allows us to calculate the maximum O-O distance in each cluster, corresponding to the minimum cut-off radius needed to include the total contribution of the many-body potential.
Finally, for each cut-off radius r, we find the minimum energy configurations with the Basin-Hopping method and calculate the average, over all the cluster molecules, of the number N i of molecules interacting.Since the global energy structure can change with the cut-off radius, the relative distance between the molecules varies, resulting in a nonmonotonic N i function, highlighting intriguing features of the global energy minimum configurations.

Figure 2 :
Figure 2: Contributions to binding energy for the minimum-energy configuration as a function of the cut-off radius r for a water cluster of 8 molecules (octamer).The Lennard-Jones contribution (blue circles) is always repulsive, while Coulomb (orange squares) and polarization (green triangles) terms are attractive.The Coulomb contribution is three times larger than the other two terms, dominating the total energy (red crosses).The three contributions show a clear correlation and a non-monotonic dependence on the cut-off radius, with a significant drop at r = 3 Å (gray dashed line).The minimum-energy configurations for the shortest and the largest cut-off are shown in Fig.1 b and g, respectively.

Figure 3 :
Figure 3: Average number N i of molecules interacting via the many-body potential for the global minimum configuration as a function of the cut-off radius r for the different water cluster sizes.

Figure 4 :
Figure 4: Root mean square deviation (RMSD) between a cut-off cluster and the reference DC cluster as a function of the cut-off radius r for (a) 6, (b) 8, (c) 10, (d) 16 and (e) 20 molecules.Vertical dashed lines mark the first value of the average number of interacting molecules N i computed where the coordination shells are complete: gray for the first and red for the second coordination shells.

Figure 5 :
Figure 5: Parametric plot of the energy deviation, compared to the reference DC value, as a function of coordination shell number for clusters with 6 (blue circles), 8 (orange squares), 10 (green pluses), 16 (red triangles) and 20 (purple crosses) molecules.The deviations are within 5% when the first coordination shell is completed, i.e., when N i 4 and r 3 (see Fig. 2 in Supporting Information).