High-Density “Windowpane” Coordination Patterns of Water Clusters and Their NBO/NRT Characterization

Cluster mixture models for liquid water at higher pressures suggest the need for water clusters of higher coordination and density than those commonly based on tetrahedral H-bonding motifs. We show here how proton-ordered water clusters of increased coordination and density can assemble from a starting cyclic tetramer or twisted bicyclic (Möbius-like) heptamer to form extended Aufbau sequences of stable two-, three-, and four-coordinate “windowpane” motifs. Such windowpane clusters exhibit sharply reduced (~90°) bond angles that differ appreciably from the tetrahedral angles of idealized crystalline ice Ih. Computed free energy and natural resonance theory (NRT) bond orders provide quantitative descriptors for the relative stabilities of clusters and strengths of individual coordinative linkages. The unity and consistency of NRT description is demonstrated to extend from familiar supra-integer bonds of the molecular regime to the near-zero bond orders of the weakest linkages in the present H-bond clusters. Our results serve to confirm that H-bonding exemplifies resonance–covalent (fractional) bonding in the sub-integer range and to further discount the dichotomous conceptions of “electrostatics” for intermolecular bonding vs. “covalency” for intramolecular bonding that still pervade much of freshman-level pedagogy and force-field methodology.


Introduction
The earliest applications of ab initio natural bond orbital (NBO) analysis [1][2][3][4] consistently revealed a "donor-acceptor" (resonance-covalency-type "charge transfer") picture of hydrogen bonding that was sharply at odds with then-prevalent "electrostatic" conceptions of intermolecular interactions [5,6]. Although the IUPAC Gold Book definition of H-bonding was subsequently revised to acknowledge the importance of covalency in Hbonding [7], superficial "dipole-dipole" rationalizations of H bonding continue to survive in many freshman-level expositions [8]. Arguments against the charge-transfer picture or in support of classical-type long-range, multipole, or "electrostatically driven" conceptions of H-bonding continue to appear [9,10] (vs. replies in [11][12][13]) in the research literature, and similar simplifying approximations persist in the empirical force fields of popular molecular dynamics (MD) simulation methods [14] that are commonly adopted to describe H-bonding in condensed phases.
The daunting task of describing macroscopic phases of liquid water or other H-bonded fluids may seem to demand the drastic long-range approximations of intermolecular ("noncovalent") interactions as compared to the exchange-type ("covalent") interactions of the short-range molecular regime. However, a more practical and accurate approach to describing intermolecular H-bonding is achieved by adopting supramolecular clusters [15] {C n } as the conceptual "building blocks" of the macroscopic liquid-phase description, based on the known continuity of high-density liquid and low-density gaseous phases around the fluid critical point [16]. More specifically, quantum cluster equilibrium (QCE) theory [17][18][19] provides a practical numerical implementation of such "cluster mixture" [20][21][22][23][24][25] modeling of macroscopic phase properties, based on accurate values of electronic and vibrational properties of H-bonded {C n } clusters that can be obtained at any chosen ab initio or density functional theory (DFT) level. The key input for QCE-based thermodynamic modeling of an aqueous phase is the data set of supramolecular clusters whose self-consistent (T,P)dependent equilibrium populations are determined from the computed partition functions for each cluster by the standard methods of quantum statistical thermodynamics [26].
Among the many H-bonded fluids of practical interest, water itself presents the most studied yet still most perplexing phase behavior of the terrestrial regime [27]. Even the microscopic structure and properties of "ordinary" liquid water under near-ambient conditions remain matters of controversy [28]. Further mysteries surround the phase behavior of water at higher temperatures and pressures, where both theory [29][30][31][32] and experiments [33][34][35][36] have suggested the existence of an alternative high-density phase of liquid water that could lead to a liquid-liquid critical point and an exotic new domain of thermodynamic behavior near 220 K and 1-2 kbar.
The primary goal of present work was to computationally search for a new class of water clusters {W n } based on the quadrilateral ("windowpane") coordination motif of the cyclic tetramer ( Figure 1) that might contribute to equilibrium QCE populations in the neighborhood of the proposed high-density phase. In each case, we restricted attention to clusters that maintain maximal Grotthuss-type proton ordering for the powerful effects of cooperative stabilization [37,38], as exemplified by the clockwise ordering of in-ring OH bonds in the view of Figure 1. The near −90 • coordination angles of the windowpane class correspond to reduced next-neighbor distances and increased mass/volume ratios compared to the characteristic tetrahedral angles and chair-hexagon coordination motifs of ice-I-like clusters. The search for cooperatively stabilized windowpane clusters is organized in Aufbau fashion toward increasing numbers of fully four-coordinate sites that more adequately sample the intermolecular interactions expected to dominate in the phase behavior of the low-temperature and high-pressure regime. The resulting windowpane clusters can serve as computational input for subsequent QCE studies to examine their possible role in the equilibrium cluster distributions of the water-phase diagram. A secondary goal of this study was to characterize each computed cluster in deeper conceptual terms that can clarify distinctive features of the underlying H-bond interactions. Such characterization should include aspects of overall cluster stability, strengths of individual coordinative linkages, shifts in atomic charge distribution, and other orbital-level features of free vs. coordinated water molecules. For these purposes, we employed NBO analysis [39,40] to obtain localized descriptors of molecular and intermolecular bonding features. Of particular interest are natural resonance theory (NRT) bond orders [41], which are expected to exhibit useful correlations with bond lengths [42,43], bond energies [44,45], bond stretching frequencies [46][47][48], NMR 1 J and 1h J spin-coupling constants [49], and other experimentally measurable properties.
Although the present study of novel water clusters was primarily directed toward equilibrium thermodynamic properties, it is important to note that such studies can also yield information on the kinetics and mechanisms of water cluster reactions. This is particularly true when, as in the present case, each cluster of the class is created in a sequential Aufbau manner from a previous member, e.g., by successive dimer additions of the form where W k = (H 2 O) k is a k-mer of a chosen coordination pattern. Analogous to elementary A + B C chemical reactions, one can compute the transition state (W k ···W 2 ) ‡ and other features of the intrinsic reaction coordinate [50] (IRC) for each such cluster reaction. Similarly, for other cluster species satisfying the simultaneous QCE equilibrium conditions, standard quantum chemical methods can be employed to determine transition-state features and associated absolute rate constants along the associated IRC [51]. However, such deeper mechanistic aspects of cluster formation were not addressed in the present work.

Computational Methods
For direct comparisons with many previous chemical applications in the NBO/NRT literature [52,53], we employed the familiar B3LYP/6-311++G** level of hybrid density functional theory for all geometry optimizations and energy evaluations of the present work. As shown elsewhere [54,55], realistic treatment of thermodynamic properties requires balanced treatment of energetic (primarily electronic) and entropic (primarily vibrational) contributions to free energy. All species were fully optimized and checked for vibrational stability with standard options of the Gaussian-16 program [56]. NBO/NRT analyses were completed with the NBO7 program [57,58] in interactive G16/NBO7 configuration. Structural and orbital graphics were obtained with the NBOPro7@Jmol utility program [59]. For NRT analyses of larger clusters, keyword selections for enlarged dynamic memory and the number of resonance structures were required to obtain fully converged bond orders. Ready-to-run input files containing optimized cartesian coordinates and keyword input for each cluster are included in the Supporting Information (SI). As shown particularly in ref. [54], many DFT variants and additional "corrections" (for dispersion, counterpoise, etc.) give qualitatively similar results for individual cluster structures and relative energies, even if some choices prove "best" for a particular thermodynamic comparison. The provided SI files allow re-optimization of cluster structures for alternative method/basis levels of choice.

Sequential Aufbau of 2-, 3-, 4-Coordinate Windowpane Water Clusters
The properties of each water cluster W k of an envisioned class are dictated by its specific H-bond coordination pattern. As primary descriptors of this pattern, we expect that each water molecule may generally be involved in two-, three-, or four-coordinate H-bonding to other molecules of the cluster (with singly coordinated "dangling" molecules excluded in leading clusters of the equilibrium thermodynamic distribution). For label-ing purposes, the coordination pattern of each cluster may be usefully described by the number of quadruply (q), triply (t), or doubly (d) coordinated sites, appended as presuperscripts (viz., q,t,d W n ) to the cluster symbol. In this notation, the cyclic water tetramer of Figure 1 is labeled 0,0,4 W 4 , with each monomer doubly coordinated in chain-like linkages to the substrate.
The structural logic for sequential Aufbau construction of windowpane clusters is straightforward. Starting from an existing cluster of this class, such as the cyclic water tetramer of Figure 1, one can choose any edge-type coordination (such as that between O(1) and O(10) in Figure 1) as a "base" for a new windowpane by attaching a water dimer in parallel fashion with two new H-bonds, as shown in the left panel of Figure 2. For maximum stabilization in forming this new H-bond attachment (e.g., from emanating H (12) at O(10)), the Grotthuss-type proton ordering should be continued around the edges of the newly formed windowpane that joins to O(1). The net result of this particular attachment is that sites O(10) and O(1) become tri-coordinate (t → t + 2), while other sites remain di-coordinate, leading to an overall 0,0,4 W 4 → 0,2,4 W 6 change in labeling. Some of these clusters, such as 0,0,4 W 4 itself or the cubane-like 0,8,0 W 8 described below, are featured in many previous cluster investigations, but the emphasis here is on hierarchical families of clusters that can be associated with a well-defined mechanistic Aufbau sequence of dimer additions, particularly leading to higher four-coordinate (q-type) motifs.
By alternating the sign of folding angles between panes, such additions can be continued indefinitely in "ladder-like" procession, as shown in successive panels of Figure 2. Each panel of Figure 2 includes (in parentheses) the per-monomer energy and standard-state Gibbs free energy change with respect to free water molecules, which serve to exhibit the important cooperative (nonadditive) effects of Grotthuss-ordered coordination patterns. The first four panels ( 0,0,4 W 4 , 0,2,4 W 6 , 0,4,4 W 8 , 0,6,4 W 10 ) show the addition of successive rungs to the ladder pattern, up to the four-pane member. The ensuing 1,4,4 W 9 (row 3, left) is the alternative "2 × 2" four-pane cluster, which adopts a buckled saddle-shape deformation from planarity with a central four-coordinate monomer. From the starting two-pane ladder ( 0,2,4 W 6 ) at the upper right, one can also attempt to add another rung that curls backward (E-like) rather than forward (Z-like), but this optimizes to the cubane-like 0,8,0 W 8 cluster (row 3, right). The cubane motif becomes an evident building block for extensions to two-cube ( 4,8,0 W 12 ), three-cube ( 8,8,0 W 16 ), or longer rod-like clusters, as illustrated in the final row of the figure.  Gibbs free energy (∆G (0) ) from free water molecules in each panel.
An alternative Aufbau starting point is provided by the twisted two-pane ( 1,0,6 W 7 ) cluster shown in the upper-left panel of Figure 3. This cluster features "Möbius-like" coordination with a continuous Grotthuss-ordered chain passing twice through the unique four-coordinate central monomer to form a closed loop. Remaining panels of Figure 3 show selected clusters that are obtained by successive Grotthuss-ordered dimer additions to 1,0,6 W 7 , aimed at increasing q numbers of saturated four-coordinate sites. The resulting structures all incorporate the higher density coordination angles of the windowpane motif, but they exhibit irregular overall shapes that appear suitable as possible contributions to bulk liquid or amorphous solid phases. As seen in Figures 2 and 3, the 8,8,0 W 16 cluster ( Figure 2, lower right) achieves the largest number of three-and four-coordinate sites (q = t = 8) and the deepest per-monomer energy (−10.62 kcal/mol) in the depicted sequences. However, whether some or all of these clusters contribute significantly to known roots of the QCE equations, or whether (like the buckyball-type clathrate clusters previously studied [60]) they can serve as leading contributors to entirely new roots (phases) of the QCE phase diagram remains to be investigated.  It is evident that each Aufbau cluster shown in Figures 2 and 3 may have alternative isomeric rearrangements of the proton network without altering the q/t/d descriptors of O···(H)···O coordination linkages. Such alternative q,t,d W n (alt) isomers may have higher point group symmetry, different proton orderings (e.g., Grotthuss cycles around individual panes rather than overall periphery), and higher or lower energy than the Aufbau-derived clusters described above. Figure 4 displays two such alternative high-symmetry forms of the 0,2k,4 W 2k+4 (sym) sequence (k = 1, 2), with respective C s (k = 1), C i (k = 2) symmetry. The C s -symmetric 0,2,4 W 6 (Cs) structure (Figure 4, left) is slightly higher in energy than 0,2,4 W 6 of Figure 2, but C i -symmetric 0,4,4 W 8 (Ci) (Figure 4, right) is slightly lower in energy than its low-symmetry counterpart in Figure 2. The inherent chirality of the coordination pattern about each O atom of higher-coordinated water clusters of Figures 2 and 3 indicates that reduced symmetry (net chirality) is a high-probability feature of equilibrium water cluster distributions in any phase involving their participation. Note that although H-bonds are considered weak noncovalent attractions, the cumulative energy release from larger cluster formation (viz., ∆E ≈ 170 kcal/mol for the 8,8,0 W 16 cluster) can readily exceed that necessary to dissociate a strong covalent bond, as in the ion pair clusters involved in self-dissociation (pH) of liquid water [54,55]. The per-monomer free energies of formation shown in Figures 2 and 3 remain slightly positive under standard-state conditions, but the windowpane clusters are expected to gain increased stability relative to the ice-like clusters of the near-ambient regime as pressure increases. Full thermochemical and vibrational spectroscopic values for each cluster are included with the optimized coordinates in SI.

Natural Atomic Charge and Bond Order Characterizations
Among the many descriptors provided by NBO analysis, the natural atomic charges {Q A } and interatomic bond orders {b AB } are most intimately associated with traditional empirical concepts of chemical bonding theory. Long-held perceptions of dichotomy between intra-vs. intermolecular forces (viz., "covalency" for chemical bond formation (b AB = 1, 2, 3,...) vs. "electrostatics" for H-bond formation (b H···O ≈ 0.1-0.2)) have long impeded true progress in the supramolecular domain. Demonstrations of how quantal Q A , b AB descriptors extend seamlessly across the supposed divide can therefore serve to refute the obsolete dipole-dipole conceptions of H-bonding (and other so-called "noncovalent" interactions) that still pervade freshman-level pedagogy and classical force-field methodology. In the present section, we wish to test the usefulness of NBO/NRT-based Q A , b AB descriptors when applied to the large data base of windowpane water clusters as described above.

General Features of Donor-Acceptor Interactions in Water Clusters
In every H-bond of every water cluster, NBO analysis reveals the characteristic n O →σ* OH donor-acceptor ("charge transfer") interaction that transfers a slight electronic charge (Q CT ) from the oxygen lone pair (n O ) of the Lewis base (LB) site into the valence antibond (σ* OH ) of the proximal Lewis acid (LA) site. Figure 5 depicts the n O -σ* OH interaction for one of the H-bonds of W 4c , showing the strongly overlapping forms of pre-orthogonal PNBOs deep inside van der Waals contact. The insets show details of the interaction that are routinely provided in NBO output, including (in kcal/mol; upper right) the second-order perturbative estimate of n O -σ* OH donor-acceptor attraction (∆E CT (2) ), the corresponding steric opposition of n O -σ OH donor-donor repulsion (∆E steric ), and the net binding energy (∆E net ). The known high transferability of NBOs [61] then assures that the individual n O , σ* OH orbitals are quite similar to those in water monomer and dimer as well as other windowpane clusters. However, one can also recognize the slight misalignments of ring strain that lower PNBO overlaps throughout the windowpane series and lead to the nuances in charge distribution, structure, and bond strength discussed below. Alternatively, the effects of n O(4) →σ* O(1)H(2) interaction can be quantified by deleting this single specific matrix element from the DFT calculation (with standard $DEL keylist options [62]) and recalculating the energy and reoptimized geometry as though it were absent in nature. As shown in Figure 6, this single deletion "breaks" the O(4)···H(2)−O(1) hydrogen bond (and initial S 4 symmetry) to give an open-chain structure with R O(1)···O(4) separation increased by~0.5 Å to near-van der Waals contact distance. The monomers at each chain terminus also reorient to near coplanarity (contrary to the~120 • dihedral twisting of the two remaining monomers), thereby allowing partial re-gain of n  All such NBO-based energetic and $DEL deletion descriptors can be obtained for other windowpane clusters of Figures 2 and 3. In the following, we focus instead on subtleties of the charge distributions and H-bond strengths that relate to the interesting cooperative effects of the highly ordered proton patterns ("water wires") formed by the H-bond networks.

Natural Atomic Charge Distributions
In principle, the simple water dimer (W 2 ) might be seen as the fundamental conceptual building block for studies of electronic charge distribution and stability in clusters of higher complexity. However, Figure 7 exhibits the detailed comparisons of H atom (italic) and O atom (plain text) natural charges in W 2 vs. cubane-like 0,8,0 W 8 to show the surprising contrasts between these species. In the two panels of Figure   How can the conflicting charge patterns of Figure 7 be rationalized? At the termini of each H-bond are two water monomers that can be identified as the LB (formal e-pair "donor") and LA (formal σ* OH "acceptor" vacancy). In the simple water dimer, the n O → σ* OH donor-acceptor interaction necessarily results in net charge transfer (ca. 0.017e) from LB to LA (Figure 7, left), resulting in the LB δ+ ···LA δ− charge pattern. However, in more complex water clusters, the surroundings of any chosen H-bond may be seen as a network of "water wires" that allow charge to redistribute as necessary to optimize overall cluster stability. Specifically, the multiple network connections allow electronic charge to be redistributed to achieve near neutrality at qor d-coordinated sites, whose equal numbers of donor and acceptor interactions can be tuned to avoid capacitive build-up. However, at t-coordinated sites, which necessarily have an imbalance of donor (t d , LB) vs. acceptor (t a , LA) connections, it becomes advantageous to confer excess anionic charge on t d sites (increasing LB strength) and cationic charge on t a sites (increasing LA strength), thus leading to the commonly observed LB δ− ···LA δ+ charge pattern.
To illustrate these propensities of cluster charge distribution, Figure 8 displays selected Q O (plain text) and Q H (italic) atomic charges of the 8,0,8 W 16 cluster for two q-type sites (centered at O(1), O(13)) and one d-type site (at O(46)), showing the significantly reduced net monomer charges compared to those of the water dimer. More complete listings of monomer charge values and coordination type at each O atom for all clusters of this study are included in SI, as illustrated for the 8,0,8 W 16 cluster in Table 1. The subtle variations in molecular charge indicate the extreme "feedback" sensitivity to every detail of the surrounding H-bond network, showing that overall network topology has taken precedence over characteristics of the water dimer (single H-bond) "building block" of which the network is composed.

Natural Bond Order Correlations
The distended shapes of windowpane clusters provide clear evidence of the severe effects of "ring strain" in altering the network O−H···O bonds from the idealized geometries of isolated H-bonds in binary complexes. Nevertheless, one expects that network H-bonds should continue to exhibit the robust correlations with NBO/NRT measures of bond order and charge transfer that were previously demonstrated for free binary H-bonded species [63]. We now turn to examining the supramolecular extension of such correlations for the classical bond order-bond length (BOBL) relationships that have long been fruitfully employed in the integer (single-, double-, triple-, etc., bond) range of covalent bonding in molecules [42,43].
A simple example of such BOBL correlations is illustrated in Figure 9 for the 1,4,4 W 9 windowpane cluster of Figure 2. For each O···H−O linkage, the total b O···O bond order is obtained as the sum of b O···H (major) and "long-bond" [64] b OˆO (minor) contributions, with sub-integer values ranging from 0.02 to 0.18 in this simple cluster. As shown in the right panel, the BOBL correlation is of excellent quality (Pearson correlation coefficient χ ≈ −0.97), and the least-squares regression line (shown in the inset) allows close prediction of R O···O distances to near the 0.01Å level(!), despite the fact that NBO/NRT descriptors receive no input from real-space molecular geometry or spatial distribution of electron density. Thus, the resonance-covalency concepts underlying NRT bond order evaluations appear to extend seamlessly into this sub-integer range of weak H-bonding in clusters, practically as well as the familiar supra-integer range of strong covalent bonding and resonance in molecules.  Table 2. In this case, the b O···O -R O···O correlation is found to be weaker, but still of reasonably high quality (χ ≈ −0.91), reflecting the heterogeneities of higher-coordination motifs. It is also of interest to examine the global BOBL correlations for all windowpane clusters of the present work, covering ca. 250 individual b O···O -R O···O H-bonded pairs in a broad variety of coupled coordination motifs. Figure 10 displays the BOBL scatter plot, least-squares regression line, and Pearson correlation coefficient for this entire data set of hydrogen bonds, showing the strong correlation (χ = −0.90) that persists in spite of increasingly heterogeneous cluster topologies. The degraded accuracy of the linear least-squares regression fit in Figure 10 (compared, e.g., to that in Figure 9) can be primarily attributed to the upward deviations from linearity that are evident near b ij → 0. However, it is important to recognize that these deviations are required on physical grounds, because intermonomer separation should asymptotically diverge (R ij → ∞) as bond order vanishes (b ij → 0). Indeed, only the higher-order connectivity of the H-bond network prevents such asymptotic dissociation when any single H-bond is severed, so the proper appearance of such nonlinearity in the b ij → 0 limit serves to further confirm the resonance-covalent nature of H-bonding even in this range of interaction strengths near the limit of chemical interest.

Conclusions
In the present work, we have employed standard density functional methods to computationally characterize a broad variety of unusual "windowpane" clusters that may play a role in the high-density fluid phase(s) of water. Despite their diverse topological forms and unusual angular features, we have demonstrated that these clusters are fully compliant with water's known facility in forming doubly (d-type), triply (t-type), and quadruply (q-type) coordinative linkages to other water molecules, leading to multiply connected ("water-wired") networks of increasing energetic stability when proper Grotthuss-type proton ordering is maintained. The Aufbau construction approach also suggests the mechanistic sequence by which such Grotthuss-ordered clusters can readily form from successive aggregation with water dimers.
We have also employed natural bond orbital (NBO) and natural resonance theory (NRT) analysis tools to demonstrate the consistency and accuracy with which H-bonding in these clusters conforms to the general conceptual picture of resonance-covalency ("charge transfer") as the authentic origin of intermolecular O−H···O attractions. The charge flows and adaptive bond order and structural shifts in these clusters are shown to obey familiar bond order-bond length (BOBL) correlations with high accuracy (|χ| > 0.9). Moreover, the BOBL correlations also exhibit the expected deviations from linearity in the asymptotic limit of vanishing bond order where R O···O distance becomes divergent. Although connections can be shown between NBO and Bader-type descriptors [65], we believe that the NRT bond orders of the present work provide broader predictive utility and more nuanced inclusion of resonance effects than the topological descriptors as employed in previous studies of water clustering (e.g., [66]).
The reader is reminded that "correlation is not causation." Nevertheless, the continuity of robust BOBL correlations that stretch across the broad extremes of supramolecular (subinteger) vs. molecular (multi-integer) bond orders strongly implies their shared origin in unified "covalency" concepts, contrary to the dichotomous viewpoint that still dominates freshman-level teaching of chemical principles and many facets of force-field methodology.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/molecules27134218/s1. The Supporting Information (SI) file contains optimized geometrical coordinates, NBO/NRT keyword input, and other computational details in ready-to-run Gaussian input files for all equilibrium water clusters described in the paper. The file also contains tables of computed natural atomic charges and natural bond orders for all water clusters of the work. Acknowledgments: Thanks are due to Eric Glendening (Indiana State University) for assistance in wrangling consistent NRT bond orders for the many challenging H-bonding interactions in this study.