Exploring the Non-Covalent Bonding in Water Clusters

QTAIM and source function analysis were used to explore the non-covalent bonding in twelve different water clusters (H2O)n obtained by considering n = 2–7 and various geometrical arrangements. A total of seventy-seven O−H⋯O hydrogen bonds (HBs) were identified in the systems under consideration, and the examination of the electron density at the bond critical point (BCP) of these HBs revealed the existence of a great diversity of O−H⋯O interactions. Furthermore, the analysis of quantities, such as |V(r)|/G(r) and H(r), allowed a further description of the nature of analogous O−H⋯O interactions within each cluster. In the case of 2-D cyclic clusters, the HBs are nearly equivalent between them. However, significant differences among the O−H⋯O interactions were observed in 3-D clusters. The assessment of the source function (SF) confirmed these findings. Finally, the ability of SF to decompose the electron density (ρ) into atomic contributions allowed the evaluation of the localized or delocalized character of these contributions to ρ at the BCP associated to the different HBs, revealing that weak O−H⋯O interactions have a significant spread of the atomic contributions, whereas strong interactions have more localized atomic contributions. These observations suggest that the nature of the O−H⋯O hydrogen bond in water clusters is determined by the inductive effects originated by the different spatial arrangements of the water molecules in the studied clusters.


Introduction
One of the most important principles in chemistry and material science states that the physicochemical properties of the matter do not depend solely on the type of atoms composing it, but they are also the result of the main interactions between its constituents at a given geometry [1]. Within this context, it can be claimed that, as the interatomic bonds determine the properties of the molecules [2,3], the interaction between molecules (and other building blocks) defines the main characteristics of the substances [4][5][6][7][8]. The latter is particularly evident in condensed phases where it has been deduced that the intermolecular interactions obey an interplay between dispersive and electrostatic forces [9]. Therefore, different types of intermolecular interactions can be readily defined in terms of the dispersive or electrostatic degree of the involved forces. The latter, also referred to as non-covalent interactions, are significantly important since they play a fundamental role in several scenarios, such as: the conformation of molecular crystals [10,11], the pairing of amino acids and their protein folding [12,13], the supramolecular assembly [14], the DNA structural stabilization [15], and others.
Hydrogen bonds, referred to as HBs hereinafter, are a particular case of non-covalent interactions that have gained a great deal of attention due to their peculiar nature and their Int. J. Mol. Sci. 2023, 24, 5271 2 of 16 ubiquitous presence in many chemical and biochemical systems. On the one hand, the definition of the HB as a classical bond is still at the center of intense scientific debate [16], but on the other hand, HBs with strengths comparable to weak covalent bonds have been reported [17]. This contradictory behavior has been explained by considering that, besides the dispersive and electrostatic components governing HBs, these interactions are highly dependent on the inductive (i.e., non-local) character of the chemical environment in which they are formed [9]. It is important to point out that the latter behavior explains the electron delocalization existing between the HB donor and acceptor water molecules, which in turn results in the partial shared-closed shell character of HB interactions [18]. Due to the aforementioned, HB networks are known to possess a cooperative (i.e., non-additivity) behavior. In this regard, it must be indicated that a thorough comprehension of the latter is an essential requirement for the development of reliable models for solvents. However, its full description through experimental as well as theoretical means remains a very challenging task [19,20].
Water is a solvent, considered as essential to biological and chemical processes [15], where HBs have certainly a key role. Thus, exploration of the structure and binding properties of small water clusters is of significant importance since it provides insights on the properties and behavior associated to the different condensed phases of water. Liquid water, the most common form on earth, has several unique properties, such as a maximum density occurring at 4 • C, a relatively high vapor pressure, the expansion of its volume upon freezing, and its exceptionally high surface tension [21]. It is also well-known that long-range order does not exist in liquid water [21,22]. Thus, the understanding of its local structure, i.e., the HB networks in water clusters, is inextricably tied to properly describing its distinctive properties. In addition to liquid phase, multiple ice phases, which dwell in the complicated phase diagram of this system, are also derived from distinct local environments, resulting in a wide range of densities, lattice energies, and other thermodynamic properties [21][22][23][24].
As commented before, the study of the HB networks present in water clusters is essential to obtain a proper description of water as a solvent. In this vein, small water clusters, (H 2 O) n with n < 10, have been intensively studied at different levels of calculation [25][26][27]. These reports have suggested that, on a microscopic scale, the density heterogeneity in liquid water is directly related to its local structure [28][29][30][31]. Related to the latter, Chaplin proposed a water model capable of switching from lower to higher density forms without breaking HBs [30]. For this dynamic behavior to be achieved, the model was constructed from various HB arrangements formed by a combination of solely hexamer and pentamer substructures [28,30]. Thus, a major finding of these works is that a precise description of the non-covalent bonding in small water clusters is accurate enough to investigate the structural changes as well as the dynamical behavior of water in its liquid state.
Before closing the present section, it is important to point out that many of the interesting properties of water have their origin in the non-additivity nature of the HBs (see above). In this context, Koehler et al., conducted ab initio calculations on the ground state of the linear water dimer with Cs symmetry and the cyclic water tetramer with S4 symmetry. This work demonstrated that an energy gain of 29% results from the cooperativity effects as estimated by using parameters based on two-body non-neighbor interaction energy as well as three-and four-body contributions [32]. Later, Suhai demonstrated that HB networks in ice are also the result of highly cooperative behavior. Moreover, this study also found that the cohesive energy of ice is the consequence of a complex interplay between repulsive and attractive terms, which exhibits different long-term interactions [33]. Likewise, other researchers investigated the significance of all higher-order interaction energy components, particularly the three-body factor among the non-additive terms for water clusters ranging from the trimer to the hexamer [34,35]. From the latter, it was suggested that non-additivity is more significant in HB networks belonging to small water clusters, where donor-acceptor arrangements involve the largest number of water molecules. At this point, it is important to comment that studies, such as those mentioned above as well as works cited in Refs. [36][37][38], are focused on structural or energetic aspects of the water clusters, but a detailed description of the nature of the non-covalent bonding on water clusters is still missing.
In this work, the HB network of small water clusters, ranging from the dimer to the heptamer, with different geometrical arrangements will be described by means of the quantum theory of atoms in molecules, QTAIM, and further analysis performed by employing the source function, SF. From there, the nature of the different HBs, O−H· · · O, found in the studied clusters will also be highlighted using quantities related to the Laplacian of the electron density. Regarding the adopted methodology, it must be indicated that it has been demonstrated to produce an accurate description of the intermolecular interactions in ion-solvent systems [39].

Geometric and Energetic Features of Water Clusters
The clusters adopted in the present study correspond to the models previously employed byŘezáč et al. [40], which have been recognized to provide an adequate representation of the water system. The atoms of each cluster involved in the HB interaction were denoted as O d -H· · · O a , referring to donor and acceptor oxygen atoms, as depicted in Figure 1. This notation will be upheld along the discussion. To understand the differences between the HBs formed in each cluster, the geometrical aspects of all the clusters are to be described. Moreover, for the sake of clearness, the models are divided into two subgroups: (i) 2-D clusters comprising the cyclic clusters (H 2 O) n , with n = 3-6 and (ii) 3-D clusters comprising the cage-like clusters (H 2 O) n , this being exclusively for n = 6-7.
water clusters ranging from the trimer to the hexamer [34,35]. From the latter, it was suggested that non-additivity is more significant in HB networks belonging to small water clusters, where donor-acceptor arrangements involve the largest number of water molecules. At this point, it is important to comment that studies, such as those mentioned above as well as works cited in Refs. [36][37][38], are focused on structural or energetic aspects of the water clusters, but a detailed description of the nature of the non-covalent bonding on water clusters is still missing.
In this work, the HB network of small water clusters, ranging from the dimer to the heptamer, with different geometrical arrangements will be described by means of the quantum theory of atoms in molecules, QTAIM, and further analysis performed by employing the source function, SF. From there, the nature of the different HBs, O−H···O, found in the studied clusters will also be highlighted using quantities related to the Laplacian of the electron density. Regarding the adopted methodology, it must be indicated that it has been demonstrated to produce an accurate description of the intermolecular interactions in ion-solvent systems [39].

Geometric and Energetic Features of Water Clusters
The clusters adopted in the present study correspond to the models previously employed by Řezáč et al. [40], which have been recognized to provide an adequate representation of the water system. The atoms of each cluster involved in the HB interaction were denoted as Od-H···Oa, referring to donor and acceptor oxygen atoms, as depicted in Figure  1. This notation will be upheld along the discussion. To understand the differences between the HBs formed in each cluster, the geometrical aspects of all the clusters are to be described. Moreover, for the sake of clearness, the models are divided into two subgroups: (i) 2-D clusters comprising the cyclic clusters (H2O)n, with n = 3-6 and (ii) 3-D clusters comprising the cage-like clusters (H2O)n, this being exclusively for n = 6-7. Analysis of the equilibrium geometry of the models showed in Figure 2 (see Section 3) shows that the 2-D clusters are characterized by the shortening of the average Od···Oa distance as the cluster size increases. As reported in Table S1 (see Supplementary Materials), the Od···Oa interatomic distance changes from 2.78 Å to 2.70 Å for (H2O)2 and (H2O)6, respectively. Conversely, the donor covalent Od-H bond distance increases with the ring size, going from 0.87 Å to 0.98 Å for (H2O)2 and (H2O)6, respectively. Previous reports [34,35,41,42] have attributed this behavior to the non-additive, many-body nature of the systems. However, this trend is no longer observed in 3-D clusters. In those cases, the molecular arrangement allows the formation of a larger number of HBs with respect to the number of water molecules of the system. For instance, in the 2-D cluster 6a, the water molecules donate only one hydrogen atom, whereas in the 3-D cluster 6d three water molecules donate their two hydrogen atoms, increasing by three the number of total HBs (see Table 1). It is worth of mentioning that the additional HBs adopt a linear arrangement, which, as reported by Gosh et al. [43] for the water dimer, represents the most favorable geometry. However, the additionally formed HBs have, on average, longer distances than the ones observed in cyclic systems. For instance, in the case of cluster 6b, the HB Od···Oa Analysis of the equilibrium geometry of the models showed in Figure 2 (see Section 3) shows that the 2-D clusters are characterized by the shortening of the average O d · · · O a distance as the cluster size increases. As reported in  6 , respectively. Previous reports [34,35,41,42] have attributed this behavior to the non-additive, many-body nature of the systems. However, this trend is no longer observed in 3-D clusters. In those cases, the molecular arrangement allows the formation of a larger number of HBs with respect to the number of water molecules of the system. For instance, in the 2-D cluster 6a, the water molecules donate only one hydrogen atom, whereas in the 3-D cluster 6d three water molecules donate their two hydrogen atoms, increasing by three the number of total HBs (see Table 1). It is worth of mentioning that the additional HBs adopt a linear arrangement, which, as reported by Gosh et al. [43] for the water dimer, represents the most favorable geometry. However, the additionally formed HBs have, on average, longer distances than the ones observed in cyclic systems. For instance, in the case of cluster 6b, the HB O d · · · O a distance is larger by 0.04 Å than in 6a. Despite this, three HBs are shorter than the average O d · · · O a distance of 6a. Moreover, in 6b, the donor covalent O d -H bond distance is, on average, larger by 0.06 Å than 6a. Similar observations are made on cage-like clusters 6c and 6d, where the O d · · · O a distances are, on average, 0.08 Å and 0.10 Å larger than those observed in 6a, whereas the donor covalent O d -H bond distances are larger by 0.12 Å and 0.18 Å. Finally, clusters 7a-c, corresponding to 3-D cage-like structures (Figure 2), are characterized by O d · · · O a distances larger than those in 6a but shorter than those in 6d. These observations point out that clusters 6b-6d and 7a-7c have a greater diversity of O d -H· · · O a HBs. It is important to mention that the O· · · O distances observed in clusters (H 2 O) 6 and (H 2 O) 7 are between the O d · · · O a distances reported for the liquid water [44] and hexagonal Ice I h [45]. distance is larger by 0.04 Å than in 6a. Despite this, three HBs are shorter than the average Od···Oa distance of 6a. Moreover, in 6b, the donor covalent Od-H bond distance is, on average, larger by 0.06 Å than 6a. Similar observations are made on cage-like clusters 6c and 6d, where the Od···Oa distances are, on average, 0.08 Å and 0.10 Å larger than those observed in 6a, whereas the donor covalent Od-H bond distances are larger by 0.12 Å and 0.18 Å. Finally, clusters 7a-c, corresponding to 3-D cage-like structures (Figure 2), are characterized by Od···Oa distances larger than those in 6a but shorter than those in 6d. These observations point out that clusters 6b-6d and 7a-7c have a greater diversity of Od-H···Oa HBs. It is important to mention that the O···O distances observed in clusters (H2O)6 and (H2O)7 are between the Od···Oa distances reported for the liquid water [44] and hexagonal Ice Ih [45]. Here 2a shows the dimer geometry, 3a, 4a, 4b, 5a and 6a shows the geometries of the 2D ring clusters. Finally, 6b, 6c, 6d, 7a, 7b and 7c depict the geometries for the 3D clusters. Here 2a shows the dimer geometry, 3a, 4a, 4b, 5a and 6a shows the geometries of the 2D ring clusters. Finally, 6b, 6c, 6d, 7a, 7b and 7c depict the geometries for the 3D clusters.
Regarding the energetic aspects, Table 1 (1) From Table 1, it can be noticed that our counterpoised corrected binding energies computed at the M062X/aug-cc-pVTZ level are comparable with the values obtained at CCSD(T)/CBS/CBS level reported in Ref. [46], a fact that validates the method employed in the present study. It is also observed in Table 1 that the BE per water molecule increases with respect to the number of molecules of the system, as expected. Although the latter definition of BE is customarily used to describe the energetic features of molecular aggregates, in the case of water clusters, it is also interesting to consider the BE per HB, because it provides an average of the stabilization energy gained per HB formed in each system. As reported in Table S2 (see Supplementary Materials), the BE per HB systematically increases with the system size for the case of 2-D cyclic clusters. However, when the more complex 3-D structures are considered, this value decreases as the number of HBs increases. The latter is particularly clear in the hexamers, where, given a fixed number of molecules (i.e., n = 6), the HBs go from 6 to 9. The data summarized in Table S2 shows that the BE per HB is 7.90 in the 2-D cyclic cluster 6a, whereas this quantity decreases to 5.71 in the 3-D cluster 6d. The last result allows us to conclude that the more complex the geometrical arrangement of the cluster, the lower the average contribution of each HB to the total energy of the system, suggesting that, while in 2-D structures the HBs formed are equivalent, in the 3-D systems a larger diversity of HBs can be expected. Despite the significant differences between the BE per water molecule and its counterpart with respect to HBs formed, the plot depicted in Figure 3 shows that BE[(H 2 O) n ]/n increases with the total number of HB formed. Thus, these results indicate that the BE[(H 2 O) n ] depends on two main factors: (i) the cluster size and (ii) the number of formed HBs, a behavior that has been also attributed to the non-additive nature of the HB network ( Figure 3) [34,35,41,42]. The data summarized in Table 1 show that, in the case of the 2-D tetramers (H 2 O) 4 , a difference of 0.26 kcal mol −1 is observed between 4a and 4b, this being attributed to the relative orientation of free hydrogen atoms: in 4a, the free hydrogen atoms of neighboring molecules are oriented at the same side of the main plane (Figure 2c), whereas the hydrogen atoms show an alternate orientation in 4b (Figure 2d). On the other hand, clusters 6a, 6b, and 6c have slightly lower binding energies than 6d (0.26 kcal mol −1 , 0.40 kcal mol −1 and 0.03 kcal mol −1 , respectively). Finally, the 3-D clusters 7a, 7b, and 7c, can be considered as quasi-degenerated structures. In the next sections, the non-covalent bonding in all the considered water clusters will be analyzed. whereas the hydrogen atoms show an alternate orientation in 4b ( Figure 2d). On the other hand, clusters 6a, 6b, and 6c have slightly lower binding energies than 6d (0.26 kcal mol −1 , 0.40 kcal mol −1 and 0.03 kcal mol −1 , respectively). Finally, the 3-D clusters 7a, 7b, and 7c, can be considered as quasi-degenerated structures. In the next sections, the non-covalent bonding in all the considered water clusters will be analyzed.

QTAIM
As the second stage of this work, the HB interactions in each cluster were analyzed using the quantum theory of atoms in molecules (QTAIM) [47]. This theory represents a powerful tool for describing different types of chemical interactions, including non-covalent interactions. Some examples of this affirmation can be found in Refs. [48][49][50][51]. Within the QTAIM, a zero-flux surface defines each atom's boundary in the electron density gradient vector field. This partition allows the separation of regions, Ω, identified as atoms in molecules. Moreover, the topological analysis of the electron density allows the identification of a set of critical points that occurs when the gradient vanishes (i.e., ∇ = 0). This set of points can be classified using another quantity derived from the electron density, that is the Laplacian of the density ∇ . This quantity can be written as the sum of contributions along the three principal axes of maximum variation, namely the eigenvalues ( ) of the Hessian matrix: ∇ = + + . Furthermore, the algebraic sum of the signs of , usually referred to as the signature (σ), together with the number of the non-zero curvatures, also known as rank (ω), allows the classification of the electron density critical points. Four different types of critical points can be found: (i) local maxima (3, −3), associated to nuclear critical points; (ii) saddle points (3, −1), associated to bond critical points; (iii) saddle points (3, +1), associated to ring critical points; and (iv) local minima (3, +3), associated to cage critical points. In this work, the HBs will be characterized by evaluating different quantities derived from the electron density at the bond critical points (BCP). The importance of the latter ones can be understood if two important aspects are considered: (i) at the BCP, ( ) has the smallest value of the electron density along the bond path, and (ii) it corresponds to the maximum of density at the interatomic surface between the two atoms. Therefore, these two important properties of ( ) will be used to characterize the HBs in the considered water clusters. In all the clusters, two different BCPs can be found, one corresponding to the covalent Od-H polar bonds and the second related to the non-covalent HB interaction Oa···H, following the notation depicted in Figure 1. As we mention before, the uniqueness of the

QTAIM
As the second stage of this work, the HB interactions in each cluster were analyzed using the quantum theory of atoms in molecules (QTAIM) [47]. This theory represents a powerful tool for describing different types of chemical interactions, including noncovalent interactions. Some examples of this affirmation can be found in Refs. [48][49][50][51]. Within the QTAIM, a zero-flux surface defines each atom's boundary in the electron density gradient vector field. This partition allows the separation of regions, Ω, identified as atoms in molecules. Moreover, the topological analysis of the electron density allows the identification of a set of critical points that occurs when the gradient vanishes (i.e., ∇ρ = 0). This set of points can be classified using another quantity derived from the electron density, that is the Laplacian of the density ∇ 2 ρ. This quantity can be written as the sum of contributions along the three principal axes of maximum variation, namely the eigenvalues (λ i ) of the Hessian matrix: ∇ 2 ρ = λ 1 + λ 2 + λ 3 . Furthermore, the algebraic sum of the signs of λ i , usually referred to as the signature (σ), together with the number of the non-zero curvatures, also known as rank (ω), allows the classification of the electron density critical points. Four different types of critical points can be found: (i) local maxima (3, −3), associated to nuclear critical points; (ii) saddle points (3, −1), associated to bond critical points; (iii) saddle points (3, +1), associated to ring critical points; and (iv) local minima (3, +3), associated to cage critical points. In this work, the HBs will be characterized by evaluating different quantities derived from the electron density at the bond critical points (BCP). The importance of the latter ones can be understood if two important aspects are considered: (i) at the BCP, ρ(r BCP ) has the smallest value of the electron density along the bond path, and (ii) it corresponds to the maximum of density at the interatomic surface between the two atoms. Therefore, these two important properties of ρ(r BCP ) will be used to characterize the HBs in the considered water clusters.
In all the clusters, two different BCPs can be found, one corresponding to the covalent O d -H polar bonds and the second related to the non-covalent HB interaction O a · · · H, following the notation depicted in Figure 1. As we mention before, the uniqueness of the BCP allows the characterization of the whole bond in terms of different properties of ρ(r BCP ). Figure 4a depicts the ρ(r BCP ) dependence with the two oxygen-hydrogen interatomic distances; namely: O d -H and O a · · · H. From Figure 4, two sets of points are clearly identified, one at shorter distances, associated with the O d -H bonds, and the other one, at larger distances, corresponding to the O a · · · H HBs. Moreover, these ρ(r BCP ) vs. d OH points fit an exponential function [ρ(r BCP ) = 5.380e −2.810d OH , r 2 = 0.9996]. This result shows that both types of bond O d -H and O a · · · H fit an empirical model previously reported for other systems [52,53]. From Figure 4b, an inverse relationship between the ρ(r BCP ) value computed at the O d -H BCP and its O a · · · H counterpart is observed. This last result suggests that a charge transfer occurs from the covalent O d -H bond to the O a · · · H HB. It is important to comment that this has also been observed in other systems characterized by N−H· · · N and F−H· · · F HB interactions [18,54]. BCP allows the characterization of the whole bond in terms of different properties of ( ). Figure 4a depicts the ( ) dependence with the two oxygen-hydrogen interatomic distances; namely: Od-H and Oa···H. From Figure 4, two sets of points are clearly identified, one at shorter distances, associated with the Od-H bonds, and the other one, at larger distances, corresponding to the Oa···H HBs. Moreover, these ( ) vs.
points fit an exponential function [ ( ) = 5.380 . , = 0.9996]. This result shows that both types of bond Od-H and Oa···H fit an empirical model previously reported for other systems [52,53]. From Figure 4b, an inverse relationship between the ( ) value computed at the Od-H BCP and its Oa···H counterpart is observed. This last result suggests that a charge transfer occurs from the covalent Od-H bond to the Oa···H HB. It is important to comment that this has also been observed in other systems characterized by N−H···N and F−H···F HB interactions [18,54].
(a) (b) The previously discussed results show that a wide variety of HBs are present in the studied models, being characterized by different Oa···H distances and ( ) values computed at the BCP. Therefore, to obtain a more complete description of the HB network, other topological properties are needed. One quantity that can be employed for this purpose is the Laplacian of the density, ∇ , because it not only provides a classification of the BCPs, but its sign can be used to distinguish between regions of concentration or depletion of the electron density with respect to the surroundings. In fact, for pure sharedshell interactions, ∇ < 0, whereas, for pure closed-shell interactions, ∇ > 0. However, in the transit region from closed-shell to shared-shell, ∇ presents a particular behavior, first increasing until reach a maximum and then decreasing [55]. The latter makes ∇ alone not fully adequate to study HB interactions that indeed belong to the aforementioned transit region. Thus, a more convenient quantity to characterize the Oa···H HBs would be the total electron energy density, ( ), defined as follows: where ( ) and ( ) are the local kinetic and potential density energies. It is worth of mentioning that ( ) is connected to the Laplacian through the local virial theorem: Because ( ) is always negative and ( ) is always positive, a local concentration of the electron density implies that the local potential energy density exceeds twice the kinetic energy density, whereas a local depletion of the electron energy corresponds to the The previously discussed results show that a wide variety of HBs are present in the studied models, being characterized by different O a · · · H distances and ρ(r BCP ) values computed at the BCP. Therefore, to obtain a more complete description of the HB network, other topological properties are needed. One quantity that can be employed for this purpose is the Laplacian of the density, ∇ 2 ρ, because it not only provides a classification of the BCPs, but its sign can be used to distinguish between regions of concentration or depletion of the electron density with respect to the surroundings. In fact, for pure shared-shell interactions, ∇ 2 ρ < 0, whereas, for pure closed-shell interactions, ∇ 2 ρ > 0. However, in the transit region from closed-shell to shared-shell, ∇ 2 ρ presents a particular behavior, first increasing until reach a maximum and then decreasing [55]. The latter makes ∇ 2 ρ alone not fully adequate to study HB interactions that indeed belong to the aforementioned transit region. Thus, a more convenient quantity to characterize the O a · · · H HBs would be the total electron energy density, H(r), defined as follows: where G(r) and V(r) are the local kinetic and potential density energies. It is worth of mentioning that H(r) is connected to the Laplacian through the local virial theorem: Because V(r) is always negative and G(r) is always positive, a local concentration of the electron density implies that the local potential energy density exceeds twice the kinetic energy density, whereas a local depletion of the electron energy corresponds to the opposite situation. The use of the local virial theorem in the H(r) definition (Equation (2)), results in the following expression: where the second term measures the concentration or depletion of the charge in the interatomic region. From Equation (4), it can be deduced that the more dominant the character of V(r) the more negative are the values of H(r), which results in a further stabilization of the molecular system. In this context, analysis of the local properties of V(r BCP ) and G(r BCP ) are adequate for the description of chemical bonding and intermolecular interactions [56,57]. From all above, it can be deduced that H(r BCP ) < 0, and ∇ 2 ρ(r BCP ) > 0 are conditions associated to the transit region from closed-shell to share-shell interactions; whereas, H(r BCP ) > 0, and ∇ 2 ρ(r BCP ) > 0 are conditions associated to pure closed-shell interactions. On the other hand, the value of the |V(r BCP )|/G(r BCP ) quotient can be divided in three types: (i) |V(r BCP )|/G(r BCP ) > 2, associated to shared-shell interactions; (ii) 1 < |V(r BCP )|/G(r BCP ) < 2, associated to the transit region from closed-shell to share-shell interaction; and (iii) |V(r BCP )|/G(r BCP ) < 1 associated to pure closed-shell interactions. It is important to point out that the described method has been successfully employed to classify HBs in biuret-water systems [58]. As reported in Table S1 (see Supplementary Materials), HBs in water dimer (2a) and trimer (3a) are characterized by H(r BCP ) > 0 (0.0007 a.u. and 0.0008 a.u., respectively) and |V(r BCP )|/G(r BCP ) < 1 (0.966 and 0.988, respectively). On the contrary, in the 4a, 4b, 5a, and 6a clusters, H(r BCP ) < 0 and 1 < |V(r BCP )|/G(r BCP ) < 2. In these clusters, a change in the value of the ratio |V(r BCP )|/G(r BCP ) suggests a strong reorganization of ρ(r BCP ) associated to the shortening of the O a · · · H distance. Nevertheless, in three dimensional hexamers, 6b, 6c, and 6d, H(r BCP ) shows both positive and negative values, being the extreme cases those of the 6d cluster where H(r BCP ) varies form −0.0161 a.u. to 0.0028 a.u. In agreement, the |V(r BCP )|/G(r BCP ) ratio takes values between 0.833 and 1.387. These observations indicate that the HBs in 3-D clusters have different nature as compared to its 2-D cyclic counterparts, and they are highly dependent on the electronic and geometric features of the systems. As the O a · · · H distance shortens and the O a · · · H-O d angle becomes more linear (see Table S1), H(r BCP ) takes more negative values, whereas the |V(r BCP )|/G(r BCP ) ratio increases, suggesting that the interaction becomes stronger as observed by the increase of ρ(r BCP ) (see Figure 5). Moreover, despite the ∇ 2 ρ(r BCP ) > 0 for all the O a · · · H interactions, it can be noticed that the larger the magnitude of |V(r BCP )| with respect to G(r BCP ), the lower the ∇ 2 ρ(r BCP ) value. opposite situation. The use of the local virial theorem in the ( ) definition (Equation (2)), results in the following expression: where the second term measures the concentration or depletion of the charge in the interatomic region. From Equation (4), it can be deduced that the more dominant the character of ( ) the more negative are the values of ( ), which results in a further stabilization of the molecular system. In this context, analysis of the local properties of ( ) and ( ) are adequate for the description of chemical bonding and intermolecular interactions [56,57]. From all above, it can be deduced that ( ) < 0, and ∇ ( ) > 0 are conditions associated to the transit region from closed-shell to share-shell interactions; whereas, ( ) > 0, and ∇ ( ) > 0 are conditions associated to pure closed-shell interactions. On the other hand, the value of the | ( )|/ ( ) quotient can be divided in three types: (i) | ( )|/ ( ) > 2, associated to shared-shell interactions; (ii) 1 < | ( )|/ ( ) < 2, associated to the transit region from closed-shell to share-shell interaction; and (iii) | ( )|/ ( ) < 1 associated to pure closed-shell interactions. It is important to point out that the described method has been successfully employed to classify HBs in biuret-water systems [58]. As reported in Table S1 ( In agreement, the | ( )|/ ( ) ratio takes values between 0.833 and 1.387. These observations indicate that the HBs in 3-D clusters have different nature as compared to its 2-D cyclic counterparts, and they are highly dependent on the electronic and geometric features of the systems. As the Oa···H distance shortens and the Oa···H-Od angle becomes more linear (see Table S1), ( ) takes more negative values, whereas the | ( )|/ ( ) ratio increases, suggesting that the interaction becomes stronger as observed by the increase of ( ) (see Figure 5). Moreover, despite the ∇ ( ) > 0 for all the Oa···H interactions, it can be noticed that the larger the magnitude of | ( )| with respect to ( ), the lower the ∇ ( ) value.  Additional insights on the nature of the HBs in the different water clusters can be gained by employing the delocalization index, DI, defined as follows: where ρ(r) and Γ(r 1 , r 2 ) correspond to the one and pair densities, respectively. As its name states, the DI quantity is related to the number of electrons delocalized between atoms A and B. This quantity can be interpreted as the fraction of electrons that are shared in the classical Lewis model, and it has been successfully used as a descriptor to analyze HBs [59,60]. In this context, the DI is associated with the shared-shell character of an intermolecular bond [60,61]. DI(O a ,H) computed for the different HBs in the studied water clusters ranges from 0.041 to 0.121. In the particular case of 2-D clusters, the DI(O a ,H) values obtained for the HBs belonging to a given cluster are very close (see Table S3). However, as the number of water molecules grows, the average DI value also increases going from 0.  Figure 6 shows the plots of the bond descriptors H(r BCP ) and |V(r BCP )|/G(r BCP ) versus the DI(O a ,H) for the intermolecular O a · · · H interactions. From Figure 6, a clear relationship between the DI(O a ,H) and the two descriptors is observed, where the shared-shell character of the interaction grows as the DI(O a ,H) values increases. The latter shows that different types of O a · · · H interactions are present in water clusters as suggested before. It is worth noting that this behavior also supports the idea that HBs in water clusters do not obey a pairwise additive rule. Indeed, at this stage, it is possible to establish that the results presented above indicate that each individual water molecule affects the entire HB network. Moreover, the influence of each water molecule can vary significantly in 3-D clusters, suggesting that an accurate enough molecular potential must consider the inductive (i.e., non-local) effect of each molecule of the system over the whole HBs network. To further explore the latter conjectures, the non-local effects of the water clusters on their HB network will be analyzed by means of the source function in the following section.
Additional insights on the nature of the HBs in the different water clusters can be gained by employing the delocalization index, DI, defined as follows: where ( ) and Γ( , ) correspond to the one and pair densities, respectively. As its name states, the DI quantity is related to the number of electrons delocalized between atoms A and B. This quantity can be interpreted as the fraction of electrons that are shared in the classical Lewis model, and it has been successfully used as a descriptor to analyze HBs [59,60]. In this context, the DI is associated with the shared-shell character of an intermolecular bond [60,61]. DI(Oa,H) computed for the different HBs in the studied water clusters ranges from 0.041 to 0.121. In the particular case of 2-D clusters, the DI(Oa,H) values obtained for the HBs belonging to a given cluster are very close (see Table S3). However, as the number of water molecules grows, the average DI value also increases going from 0.066 for the dimer 2a to 0.094 for the hexamer 6a. On the other hand, as the Oa···H interatomic distance takes different values in 3-D clusters, the DI(Oa,H) also spans a wide range of values. For instance, in the cluster 6d, the strongest Oa···H has a DI(Oa,H) value of 0.115, being 2.5 times the DI(Oa,H) observed for the weakest Oa···H. Moreover, clusters, 7a, 7b, and 7c possess the highest DI(Oa,H) values, 0.199, 0.117, and 0.121, respectively, but these clusters also have the lowest DI(Oa,H) values, being 0.044, 0.051, and 0.041 for 7a, 7b, and 7c, respectively. Figure 6 shows the plots of the bond descriptors ( ) and | ( )|/ ( ) versus the DI(Oa,H) for the intermolecular Oa···H interactions. From Figure 6, a clear relationship between the DI(Oa,H) and the two descriptors is observed, where the shared-shell character of the interaction grows as the DI(Oa,H) values increases. The latter shows that different types of Oa···H interactions are present in water clusters as suggested before. It is worth noting that this behavior also supports the idea that HBs in water clusters do not obey a pairwise additive rule. Indeed, at this stage, it is possible to establish that the results presented above indicate that each individual water molecule affects the entire HB network. Moreover, the influence of each water molecule can vary significantly in 3-D clusters, suggesting that an accurate enough molecular potential must consider the inductive (i.e., non-local) effect of each molecule of the system over the whole HBs network. To further explore the latter conjectures, the non-local effects of the water clusters on their HB network will be analyzed by means of the source function in the following section.

Source Function Analysis
The non-local effects that stabilize the O a · · · H interactions present in the considered water clusters were analyzed by adopting the source function (SF) introduced by Bader and Gatti in 1998 [62]. This function can be applied to decompose the electron density at r as the sum of contributions of different local sources LS r, r , operating at all the other points of the space [63,64]. Within this context, the electron density can be written as: where the local source is given by LS r, r = − 4π r − r −1 · ∇ 2 ρ r . As observed in its definition, LS represents how the Laplacian of the electron density at the position r affects the electron density at a different position r [25]. Accordingly, ρ(r) can be written as a sum of atomic source contributions: where S, referred to as source function (SF), is defined as: Equation (7) implies that ρ(r) is given by the contribution, S(r, Ω), of each atomic basin of the complete system, indicating that SF provides a measure of the relative contribution of an atom or group of atoms to the electron density at r [47]. The physical meaning of this contributions can be understood employing the local expression of the virial theorem (Equation (3)) in the LS r, r definition to obtain the following: where G(r ) is the positively defined kinetic energy density, and V(r ) is the electronic potential energy. According to Equation (9), molecular regions, where the electron density is concentrated (V(r ) > 2G(r )), are sources for the electron density at r. On the other hand, regions where electron density is depleted (V(r ) < 2G(r )), acts as sinks for ρ(r). It is worth to mention that the capability of the electron density at r to be a source (or a sink) at another point r is related to the magnitude of its charge concentration or depletion (magnitude of ∇ 2 ρ(r )), weighted by the inverse of the distance between these two points [64]. Here, for the sake of simplicity, the source function values will be expressed as the percentage contribution to the electron density as: It is important to point out that S%(r, Ω) represents the percentage share of electron density from Ω to ρ(r). Thus, in principle, S%(r, Ω) reflects the delocalized or localized character of a chemical interaction. For two bonded atoms S%(r, Ω) will be large if the interaction is localized and will be small for highly delocalized interactions.
The SF contributions to the O a · · · H electronic density at the bond critical point, ρ(r BCP ), were examined to investigate the nature of the HB in water clusters.  Table S4). In 2-D clusters, S%(H) increases with the number of water molecules, suggesting that the positive contributions of ∇ 2 ρ r within the hydrogen basin decreases in comparison to negative contributions to O a · · · H ρ(r BCP ). In other words, in the H basin, the charge concentration regions become more important than the charge depletion regions; thus, hydrogen atoms became a more significant source as the number of water molecules increases in the system. The average source contribution at ρ(r BCP ) for the HB interactions in 2-D clusters are shown in Table 2. From this Table, it can be concluded that the larger the |V(r BCP )|/G(r BCP ) ratio, the larger S%(H) and S%(O a ), but the lower S%(O d ). The latter results also indicate that as the HB increase its shared-shell character, the ρ(r BCP ) is mainly due to the O d and O a atomic contributions; because, as was mentioned before, the hydrogen atoms act as sinks. Furthermore, the increase in S%(O d , H, O a ) values can be attributed to a localization of the atomic density charge contributions.  In other words, in the H basin, the charge concentration regions become more important than the charge depletion regions; thus, hydrogen atoms became a more significant source as the number of water molecules increases in the system. The average source contribution at ( ) for the HB interactions in 2-D clusters are shown in Table 2. From this Table,  which, in turn, can be associated to an increase in the interaction's shared-shell nature.   The SF was employed also to unveil the character of the HBs in 3-D clusters. This analysis was focused on the 3-D cluster 6d as a representative case, because it contains HBs that belong to the pure closed-shell and the transit regions, as shown by the evaluation of |V(r BCP )|/G(r BCP ), H(r BCP ), and DI(O,H) (see tables in Supplementary Material). In this case, the SF was employed to examine the localized/delocalized character of the atomic contributions to the ρ(r BCP ) of the O a · · · H interaction. These atomic contributions are calculated by integrating the local source, LS(r, r ), on the atomic basin (see Equations (7)- (9). It is important to mention that, since the topography of the electron density at any location in space affects the electron density at the O a · · · H BCP, the atomic SF contribution can be used as an alternative to explore electron delocalization [64][65][66][67]. Figure 8 displays a pictorial description of the atomic source contributions to the electron density at the BCP of different bonds and interactions. For the O d -H covalent bond, depicted in Figure 8a, the contributions of the atoms that form the bond account for 96.8% of the electron density at the BCP, whereas 1.5% of the electron density refers to the other H atom in the molecule and 1.7% corresponds to the remaining atoms in the cluster. In a strong HB (see Figure 8b), a similar behavior is exhibited, where the atomic SF contribution of the interacting atoms, S%(O d , H, O a ), is 79.7% and 20.3% is spread along the rest of the water molecules. On the other hand, Figure 8c shows the atomic SF contributions for a weak HB (i.e., |V(r BCP )|/G(r BCP ) < 1), where the SF triad contribution S%(O d , H, O a ) is only 48.1%, and the remaining 51.9% contribution is associated to the remaining atoms in the cluster. The latter observations mean that as the shared-shell character of the HB increases, the SF contributions become more localized and vice versa. The SF was employed also to unveil the character of the HBs in 3-D clusters. This analysis was focused on the 3-D cluster 6d as a representative case, because it contains HBs that belong to the pure closed-shell and the transit regions, as shown by the evaluation of | ( )|/ ( ), ( ), and DI(O,H) (see tables in Supplementary Material).
In this case, the SF was employed to examine the localized/delocalized character of the atomic contributions to the ( ) of the Oa···H interaction. These atomic contributions are calculated by integrating the local source, ( , ), on the atomic basin (see Equations (7)- (9). It is important to mention that, since the topography of the electron density at any location in space affects the electron density at the Oa···H BCP, the atomic SF contribution can be used as an alternative to explore electron delocalization [64][65][66][67].

Models and Methods
The present study considered twelve water clusters ranging from water dimer (H2O)2 to heptamers (H2O)7. The initial geometries were obtained from the Refs. [40,46] and further optimized at the M062X/aug-cc-pVTZ level of theory. In this regard, it is important to mention that a recently published benchmark study, focused on the performance of different DFT functionals in the description of hydrogen-bonded systems [68], has determined that the M062X is among the methods providing the lowest RMSE as compared with results obtained by means of highly correlated methods, i.e., CCSD(T). Moreover, this functional has been successfully tested in non-covalent benchmark databases, S66 [69,70]. At this level of calculation, no significant geometric differences were found after the optimization. The error due to the basis set superposition (BSSE) was corrected during the minimization of the energy with respect to the geometry by employing the Boys-Bernardi counterpoise method as implemented in the GAUSSIAN16 suite of programs [71]. Moreover, subsequent vibrational analysis was performed to check the nature of the Figure 8. Pictorial representation of the SF contributions to ρ(r BCP ) at (a) a O d -H covalent bond, (b) strong O a · · · H HB (i.e, 1 < |V(r BCP )|/G(r BCP ) < 2 ) and (c) weak O a · · · H HB (i.e., |V(r BCP )|/G(r BCP ) < 1 ). Blue color represents positive contributions (sources), whereas yellow color represents negative SF contributions (sinks).

Models and Methods
The present study considered twelve water clusters ranging from water dimer (H 2 O) 2 to heptamers (H 2 O) 7 . The initial geometries were obtained from the Refs. [40,46] and further optimized at the M062X/aug-cc-pVTZ level of theory. In this regard, it is important to mention that a recently published benchmark study, focused on the performance of different DFT functionals in the description of hydrogen-bonded systems [68], has determined that the M062X is among the methods providing the lowest RMSE as compared with results obtained by means of highly correlated methods, i.e., CCSD(T). Moreover, this functional has been successfully tested in non-covalent benchmark databases, S66 [69,70]. At this level of calculation, no significant geometric differences were found after the optimization. The error due to the basis set superposition (BSSE) was corrected during the minimization of the energy with respect to the geometry by employing the Boys-Bernardi counterpoise method as implemented in the GAUSSIAN16 suite of programs [71]. Moreover, subsequent vibrational analysis was performed to check the nature of the stationary points. The QTAIM and source function analysis were performed using the AIMALL package [72].

Conclusions
Non-covalent bonding interactions in twelve different water clusters were investigated in this work. An analysis of the electron density at the BCP between the O and H atoms revealed an inverse relationship between the ρ(r BCP ) value at O d -H BCP and its analogous O a · · · H, suggesting a charge transfer from the covalent O d -H bond to the O a · · · H HB. Furthermore, examination of the ∇ 2 ρ(r BCP ) values evidenced important differences between analogous O a · · · H interactions within each cluster. Whereas no noticeable changes were observed in 2D clusters, substantial changes (closed-shell vs. shared-shell character) were detected between different O a · · · H interactions in 3D clusters.
SF analysis was also employed to characterize the O a · · · H interactions. Unlike 2D clusters, in 3D clusters the S%(H) and S%(O a ) values varied greatly, which can be explained by the different geometric arrangements that adopt the water molecules. SF analyses also permitted to evaluate the localized/non-localized nature of the O a · · · H interactions via ρ(r BCP ) decomposition into atomic contributions. These findings revealed that each individual water molecule affects the entire HB network. Therefore, an accurate description of water cluster interactions beyond the dimer is crucial to obtain accurate potentials for the condensate water simulation.