Computational study of hydrogen bond interactions in water cluster-organic molecule complexes

We analyzed the interactions present in complexes that acetone, azomethane, dimethylamine, dimethyl ether, methyl acetate, and oxirane form with 39 diﬀerent (H 2 O) n clusters (n=1-10). A random generation of conﬁgurations and a subsequent screening procedure were employed to sample representative interactions. Using quantum chemical computations, we calculated the associated binding energies, ranging from -0.19 kcal/mol to -10.76 kcal/mol at the DLPNO-CCSD(T)/CBS level. It was found that this set of energies is diverse and can be elucidated in terms of various factors, including the water cluster size, the nature of the organic molecule, and type of hydrogen bond donor. We ﬁnd that the most stable complexes often arise from a combination of a strong hydrogen bond plus a secondary interaction between the organic molecule and the water cluster.


Introduction
Water can be used as a reaction medium in organic chemistry. Diels and Alder mixed furan and maleic anhydride in hot water, obtaining a cycloaddition adduct with an increased endoselectivity compared to their organic solvent counterpart. 1 However, it was not until the 1980s that a systematic and pioneering series of papers by Breslow and coworkers established the paradigm of water as useful and relevant in organic chemistry. In their research, it was demonstrated that H 2 O as a solvent increases the reaction rate of Diels-Alder reactions 700-fold, compared to reactions in organic solvents, mainly due to hydrophobic effects. [2][3][4][5] Breslow's seminal contributions set the foundations of aqueous organic chemistry (AOC) as a new field of research that is highly active today. [6][7][8] One of the key areas of activity is developing chemical procedures that reduce organic solvent usage in favour of H 2 O; 9,10 the non-toxicity, reusability, and low-cost of water perfectly align with the principles of green chemistry. 11 It is expected that further expanding research in AOC will lead to significant discoveries in chemistry.
Understanding hydrogen bonding (HB) interactions between organic molecules and water molecules is central to unravel the poorly understood mechanisms of AOC. For example, it has been proposed that a group of accelerated reactions in aqueous emulsions, called "on water" reactions, occurs due to HB stabilization of transition states relative to reactants. 12 Also, countless chemical processes relevant to biology involve HBs between water molecules and organic substrates. [13][14][15][16][17][18][19] In light of the importance of organic molecule-water cluster bonding, the present work focuses on the accurate quantum mechanical description of such interactions.
HB interactions within water clusters have been extensively studied over the past decades.
It is a well-known fact that cooperative effects are prominent in these interactions [20][21][22][23][24] and that they increase with the number of water molecules in the cluster. In addition to the cluster size, the local environment of hydrogen bonds is relevant; the nature of the HB donor and acceptor and their neighbouring molecules has been used to describe noncovalent interactions within water clusters of small size. 25,26 However, to the best of our knowledge, a comprehensive work on organic molecule interactions with water clusters of distinct sizes has not yet been explored. Some computational studies have shown that density functional theory (DFT) can reproduce experimental results in AOC reactions. 12,27,28 However, most quantum chemistry models utilize only a few water molecules to represent the complex traits of HB in these systems, completely neglecting the role of water cluster size in accurately depicting hydrogen bonding interactions with organic molecules. For these reasons, in this work, we explore the interaction strengths of water clusters of varying sizes with organic HB acceptors. In the process of our study, we developed a new benchmark data set of hydrogen bonding for organic molecule-water cluster complexes.

Computational Methods
An accurate quantum chemical description of water systems can only be achieved for small water clusters. 29,30 These limitations demand that we approximate our molecular systems with small models. Despite their reduced size, model systems can often produce insightful information for a variety of chemical problems. For our study, we wanted to understand the behaviour of binding energies between organic molecules and water clusters of increasing size. For this purpose, we retrieved the water clusters produced by Temelso and coworkers 31 that range in size from the monomer to the decamer. A summary of the structures studied in this paper is shown in Table 1, and an example of them is displayed in Figure 1. Temelso et al. obtained these structures at the RI-MP2 32,33 /aug-cc-pVDZ 34 level of theory, making them appropriate for benchmarking purposes. We also selected the following seven organic molecules as representative of species that are present in AOC reactions: acetone, dimethyl ether, oxirane, methyl acetate, azomethane, dimethylamine, and trimethylamine, which we pre-optimized at B3LYP 35-37 /6-31+G** using Gaussian 16 (g16). 38 This selection was based on their similarity to organic substrates involved in "on water" reactions. 39 We used these organic molecules as probes to study the HB donating capabilities of −OH groups in small water clusters. For the sake of simplicity, we will often refer to these organic molecules as probes in this work. Table 1: Names of the water clusters considered in this work. The structures of all oligomers of water shown here can be seen in Figures 1 and 2  The complexes formed by the interactions between the probes and water clusters of Step II employs atom-centered potentials (ACPs) that are designed to mitigate the effects of incomplete correlation and incomplete basis sets associated with the HF-D3(BJ)/MINIs approach. 42 Step III uses a form of ACPs called basis set incompleteness potentials (BSIPs) designed to reduce the effects of incomplete basis sets used with conventional density-functional theory (DFT) based methods. 43 Steps I-III employed g16, and the process is summarized in Figure 2. Figure 2: Summary of the screening process followed in selecting the most stable complexes formed between water clusters and the probes. We applied this procedure to each of the possible 273 complexes (7 probes interacting with 39 water clusters).
Following the foregoing screening, a set of 250 configurations per water cluster-probe complex was obtained. These energy-ranked configurations were divided into ten bins containing 25 conformers each. We then selected the most stable complex from each bin to get 10 structures per complex and re-optimized the geometries of the probes within the complexes at B3LYP-D3(BJ) 35-37,44-47 /6-31+G**-BSIP 43 using g16; in all cases, the structures of the water clusters were kept fixed at the Temelso geometries. 31 This approach has the advantage of using high-level geometries; however, since they are fixed structures, the donating −OH group in the water clusters will not relax and, hence, the binding energies we calculate will all be upper bounds. The resulting structures constitute the set of complexes for our benchmark calculations. We obtained 2730 complex structures in total.
From our 2730 geometries, we removed the redundant geometries that arose from config-urations that converged to identical geometries when re-optimizing the probes in the complexes. For doing so, the COMPARE feature, as implemented in the program Critic2, 48,49 was employed. The final dataset contains 2376 non-redundant geometries. Our screening method ensured that strong interactions from the dangling −OH to the organic probes were present and other weaker interactions such as those involving C−H HB donors and dispersion interactions. The molecular systems were visualized using the software Chemcraft. 50 Reference binding energy calculations were calculated using DLPNO-CCSD(T) 51 with aug-cc-pVDZ and aug-cc-pVTZ basis sets, 34 (represented as DZ and TZ, respectively). A strategy similar to the one presented in reference 52 was applied for the extrapolation to the complete basis set (CBS) limit: we calculated the counterpoise (CP) corrected binding energies (BE CP ) 53,54 of the complexes as where the general nomenclature, E basis geometry (X), represents the electronic energy of species X (R for probe, (H 2 O) n for a water cluster of size n, and complex for the water cluster-probe aggregates), at the geometry specified by the subscript and using the basis set identified by the superscript. For example, the term E R Complex (R) represents the electronic energy of a probe at its geometry in a complex using the probe basis sets. Note that Equation 1 can be written alternatively as , and the deformation energy of the organic probe, The ORCA 55 package was used for the DLPNO-CCSD(T) calculations. Equations 2 and 3 were evaluated for both basis sets, DZ and TZ. A two-point extrapolation (DZ-TZ) to the CBS limit was applied, as described in reference 56. We also computed the average of the CP and non-CP approaches as BE i ave = BE i CP + BE i non−CP /2, with i = DZ, TZ, and DZ-TZ.
We defined our benchmark binding energies as the average of the CP and non-CP energies at the CBS extrapolation limit (i = DZ-TZ). Mackie and DiLabio showed that CP BEs tend to converge to the complete basis set limit from above, while the non-CP BEs converge from below. 52 Therefore, averaging the two quantities generally results in quicker convergence to the CBS limit. These average binding energies show a quick convergence to the CBS limit for most noncovalent interactions studied herein. A selected example of the convergence of BE i ave is depicted in Figure 3; since BE CP and BE non−CP converge to the same value at the CBS limit (see right-hand sides of equations 1 and 3), BE i ave will more rapidly converge to the complete basis set limit than the individual counterpoise and non-counterpoise binding energies. However, the convergence of the average of the CP and non-CP BEs is not always ideal, i.e., when the CP and non-CP BEs do not converge from above or below, respectively.

Results and Discussion
A diverse set of water cluster-probe configurations was obtained from the process described in the previous section. Figure 6 shows all benchmark binding energies (BE) obtained for the dimethyl ether molecule and all of the different water clusters. Due to its relatively simple structure, the discussion will be focused on this probe. The BE data cover an energy range from -0.19 kcal/mol to -8.42 kcal/mol. A careful inspection of the geometries of the complexes associated with Figure 6   to the organic probe is also donating an HB to a neighbouring molecule of water: accepting electronic density from both molecules and reducing the HB donating capabilities of the central water in DD. For SDSA and A+SDSA, the water molecule donating an HB to the probe is also acting as an HB acceptor of another water molecule. This is the well-known "cooperative" effect in hydrogen bonding: The water bound to the probe is electron deficient because it is donating electron density to a neighbouring water molecule, making the former water better able to accept electron density from the probe. This effect is also at play in the A+SDSA motif, which interactions are further enhanced by the secondary A interactions.
The monomer, which lacks neighbour water molecules acting as HB acceptors or donors, falls between the two dimer cases.
A+DDDA A+DDSA A+SDSA Figure 8: Some of the noncovalent bonds defined in Figure 7 can gain additional stability when combined with type A interactions. Colour code, red: oxygen, blue: hydrogen, and magenta: carbon. Dashed lines indicate either HB or C − H · · · (H 2 O) n interactions.
For the water clusters ranging in size from the trimer to the decamer, Figure 8 shows the following order of BE strength: DDSA < DDDA < A+DDSA < A+DDDA. In these cases, all interactions involve a water acting as a double hydrogen bond donor to a different water and the organic probe, and the water accepting a single HB from a neighbouring water (DDSA) or the water accepting two HBs from two neighbouring waters (DDDA).
For the same arguments given in the dimer example, a water molecule accepting more HBs from neighbouring waters is a better HB donor than one accepting fewer HBs, therefore DDDA motifs result in stronger BEs than DDSA. Likewise, A+DDDA represents a more Examining the most stable complex of each water cluster with each of the organic probes shows that the profile of benchmark binding energies ( Figure 10) has an irregular trend with increasing n in (H 2 O) n without appearing to reach an asymptote by n − − 10. To explain this trend in the curves, we noted that the water cluster displaying the strongest HB interactions belongs to the water heptamer -this it is the largest water cluster presenting book-like geometries, which allow for A+DDDA type interactions. Of the 11 heptamer structures we obtained from the literature, 31 6 are open structures (i.e. they do not form cages or prisms) which allow for A type interactions in addition to DDDA interactions. In contrast, the  Figure 11 illustrate this using complexes formed between (H 2 O) n and dimethylamine.
These findings suggest to us that the magnitude of the BEs obtained with the heptamer are likely to be quite close to the maximum HB strengths between the probes and water clusters.
The nature of the organic probe is also an important determinator of HB binding strengths.   Binding energy (kcal/mol) Figure 12: Benchmark binding energies for dimethyl ether and individual clusters as labelled in Table 1. show a similar distribution of binding energies (see also the additional data presented in the SI). We can expect that AOC reactions involving the moieties present in our probes would behave similarly when interacting with water molecules. We imagine real water-organic reactions involving significant dynamics in HB formation between the probes and the water systems with which they interact resulting in fluctuations in BE over time. We suspect that the maximum strength of the BEs achievable will likely be quite close those calculated for the probe-(H 2 O) 7 . Of course, entropy effects will impact the nature of the complexes formed in real systems.

Conclusions
Understanding noncovalent interactions between organic molecules and water clusters is required in order to understand the chemistry of aqueous organic phenomena. Secondary interactions play a critical role in forming organic molecule-water cluster complexes. While contact between a heteroatom in the organic probe and a −OH group from the cluster is required for a strong interaction, additional interactions from aliphatic hydrogens can provide extra stabilization to the complexes. We showed that the nature of the water molecule interacting with the organic species, whether it is acting as a single or double HB acceptor or donor, determines the strength of binding energies of the complexes. The nature of the organic probe is also an important determinant of HB strengths. Therefore, the computational study of AOC reactions requires consideration of the size of the water cluster and the nature of the secondary interactions, the HB donors and acceptors, and that of the organic molecules. The presence of different moieties in the probes can lead to a more complex distribution of the BEs; for instance, methyl acetate can engage in many of the classes of interactions identified in this work. Overall, our work illustrates the intricate nature of reactions in aqueous phase. This work also provides a novel benchmark data set for organic probe-water cluster binding energies that can be used, for example, in the assessment of computational methods that describe noncovalent interactions.

Acknowledgements
GAD acknowledges the Canada Foundation for Innovation, the National Science and En-