Electronic origin of the dependence of hydrogen bond strengths on nearest-neighbor and next-nearest-neighbor hydrogen bonds in polyhedral water clusters (H₂O)_n, n = 8, 20 and 24

Abstract

The influence of the nearest neighbor and next-nearest neighbor water molecules on the strength of the hydrogen (H) bonds was examined for the polyhedral clusters of cubic (H₂O)₈, dodecahedral (H₂O)₂₀ and tetrakaidecahedral (H₂O)₂₄ cages. The relative stability and the characteristics of the H bond networks are also studied. The charge-transfer (CT) and dispersion interaction terms of every pair of H bonds are evaluated using perturbation theory based on the locally-projected molecular orbitals (LPMO PT). Every water molecule and every H-bonded pair in these polyhedral clusters are classified by the types of the neighbor molecules and H bonds. The relative binding energies among the polyhedral clusters are grouped by these classifications. The optimized O⋯O distances, which are strongly correlated with the calculated pairwise CT terms, are dependent on the 49 sub-groups of the H bonds determined by the type of the neighbor molecules. The electronic origin of this dependence is analyzed using Mulliken's charge-transfer theory, and employing a few assumptions, the analytical formulas for the contribution of the CT terms to the H bond energy are derived.

---

Introduction

Among the finite size of pure water clusters, the polyhedral forms exhibit some unique features. The cubic cage of (H₂O)₈ consists of 6 four-membered rings (4⁶),¹ the dodecahedral cage of (H₂O)₂₀ consists of 12 five-membered rings (5¹²),^1,2 and the tetrakaidecahedral cage of (H₂O)₂₄ consists of 12 five-membered rings and 2 six-membered rings (5¹²6²).^1,3 There are several graph theoretical studies for water clusters,^1,4–9 and some are specific for polyhedral forms, and the others are more general. These studies showed that there are 14 distinct isomers of the cubic cage of (H₂O)₈, while there are 30 026 symmetry-distinct configurations for the dodecahedral cage of (H₂O)₂₀,^1,5 and 3 043 836 for the tetrakaidecahedral cage of (H₂O)₂₄.^3,6 Every water molecule in these polyhedral clusters participates in three hydrogen bonds (H bonds) to neighbors; one hydrogen-donating OH and two hydrogen-accepting O⋯H (hereafter, this type of water is called daa),^2,10–15 or two hydrogen-donating OH and one hydrogen-accepting O⋯H (dda). Anick called them F and L, respectively.^7,8 Water dda (F) has a “free” OH, often called a dangling OH bond. The polyhedral cages consist of n/2 dda and n/2 daa water molecules.⁷ Therefore, there are only four types of hydrogen bonds in the polyhedral clusters, viz. daa⇐daa, daa⇐dda, dda⇐daa, and dda⇐dda. Among them, the bond dda⇐daa is the strongest.^16–18 Smith and Dang called the bonds “defect” if both of the neighbor pair are either dda or daa, suggesting that dda⇐dda and daa⇐daa are weaker than the other types of pairs.^1,2 Because of the fixed arrangement of the oxygen atoms among numerous distinct isomers and the restricted types of the hydrogen bonds in the clusters, the polyhedral water clusters are ideal to probe the characteristics of the hydrogen bonded networks and the factors determining the strength of the H bonds. The above four types of H bonds are further sub-grouped by the types of the neighboring water molecules. One of the ways is to use the “topological index” of an H bond ξ,¹⁹ which is defined as a sum of the number of dd neighbors at the H donor and the number of aa neighbors at the H acceptor. The index ranges from 0 to 4. Thus, there are 20 subsets of the types of the H bonds. Anick used the index in his analysis of the calculated OH harmonic frequencies, and further classified the H bonds by introducing the descriptors (t_1P,t_1Q) for the H acceptor and (t_2P,t_2Q) for the H donor; their values are 0 if the neighbor (secondary) water is daa, and 1 if the neighbor water is dda.^8,20,21 With this classification, there are 49 sub-groups of the H bonds (see below). Using these descriptors, he demonstrated that the H bond distances can be fitted linearly. He further extended the models to include the influence from the ternary water molecules.

Recently, one of the present authors (S.I.) has developed the perturbation theory based on the locally projected molecular orbital (LPMO PT) and applied it to several molecular clusters.^22,23 The zero-order energy is exactly free of Basis Set Superposition Error (BSSE) with any basis set, and because of the proper definition of the excited orbitals, the perturbation correction is approximately BSSE-free.²⁴ By using the augmented basis functions, and by adding the dispersion (Disp) terms properly, the method provides reasonably accurate binding energies for H bonded clusters, halogen bonded dimers and weak charge-transfer (CT) complexes.^23,25 One of the advantages of this method is that the contribution to the interaction energy from the CT terms can be expressed by a sum of the terms of every pair of H bonded molecules as the contributions from the Disp terms are. The pairwise terms of the CT and Disp terms in the water clusters have been shown to be strongly correlated with the H bond distance in water clusters,^26,27 and they have been used as an alternative measure of the strength of the H bonds in analyzing the H bonded networks in (H₂O)₂₀ and (H₂O)₂₅.²⁷ The present paper is an extension of these previous studies to representative polyhedral water clusters.

With the aid of the CT term of each H bond, the electronic origins of the dependence of the H bond strength on the types of the neighboring waters and on the H bond networks are elucidated. To analyze them, we adopt the classical model theory of the charge-transfer (electron-donor–aceptor) interaction proposed by Mulliken.²⁸ The wave function for an H bond is given by

where and stand for the hydrogen acceptor and donor.^29–31 The electron donor orbital is the non-bonding orbital of the hydrogen acceptor water molecule , and the electron acceptor orbital is the anti-bonding σ_OH* of the hydrogen donor water molecule . These frontier molecular orbitals, which participate the CT electronic configurations, undergo changes by the neighboring H bonds. Based on a few simple assumptions, the analytical forms of the CT contributions to the H bond energy of 49 sub-groups defined by Anick^8,21 are derived, and they are well correlated with the H bond distances and with the calculated CT terms, averaged over all of cube isomers of (H₂O)₈ and over 13 most stable dodecahedral (H₂O)₂₀. The model theory can be extended to any H bonded networks, which are found in more complicated caged clusters than in the polyhedral clusters.

The goals of the present study are (1) to demonstrate that the strength of H bonds is influenced by the types of the other two neighbor H bonds in the polyhedral water clusters and that accordingly, there are 49 sub-grouped H bonds, whose strengths are distinct from those of the others, and (2) to elucidate the electronic origin of this dependence in terms of the Mulliken's CT theory. (3) Besides, the relative stability among the isomers is related to the different numbers of the sub-grouped H bonds in the isomers.

Theoretical and computational procedure

The binding energy and pairwise energy

The perturbation expansion theory based on the locally projected molecular orbital (LPMO PT) has been previously developed by one of the authors.^{22,23,25,32–34} Here we describe only a few salient features that pertain to the current discussion. The zero order wave function Ψ_LPMO for a molecular cluster is a single Slater determinant constructed from sets of strictly local MOs variationally determined by a coupled set of Hartree–Fock equations. The equations are derived under the strong constraint that the occupied MOs for each molecule are expanded only with the basis sets placed on that molecule. The zero and first order wave function iswhere |LE_X〉 stands for the single excitations within molecule X and |CT_X→Y〉 stands for the single excitations from molecule X to Y, whereas |Disp_X−Y〉 denotes for the double excitations of the dispersion type. When the number of molecules in the cluster is larger than 3, the required computation time up to the single excitation PT calculation is much shorter than for the supermolecule HF calculation. This expansion is made possible by defining the locally projected excited MOs, most of which are expanded in terms of the basis sets on each molecule X.²⁴ The binding energy of the clusters is defined as The first parenthesis in eqn (3) represents the binding energy evaluated by LPMO, and defines the energy E^LPMO_BindE, which contains the electrostatic and exchange-repulsion terms and a part of induction terms. The second parenthesis in eqn (3) is the second and third order correction PT contribution from the single excitations, which is expressed by the sums The contribution from the local excitations is non-zero only at the third order and always much smaller than the CT terms, because the occupied LPMOs are canonical within each molecule. In other words, E^LPMO_BindE does contain most of the induction terms, which are defined by the second order perturbation theory expanded in terms of the products of the wave functions of the isolated monomers. According to the previous literature,^35–37 the induction terms are the contribution from the products in which one of the wave functions of the monomers is for the excited state, and the dispersion terms are that from the products in which two of the wave functions are for the excited states. The polarization terms are the sum of both. Under this formalism, the CT terms are a part of the induction terms, but note that E^LPMO_BindE in our theory does not contain the CT terms, which are separately evaluated by the second and third order PT as E_CT+LE. By using the MC⁺BS (monomer-centered plus basis sets), which allows that the excited orbitals are delocalized over molecules, the SAPT could evaluate the CT contribution.³⁷ In the polarizable model functions such as TTM2-R,³⁸ the atomic (not molecular) centered polarizability might be able to effectively include the CT contributions. The term E^LPMO_BindE in eqn (3) is the dominant term in the binding energy. The LPMO wave function Ψ_LPMO is a single Slater determinant, and in that sense, it is a kind of the supermolecule wave function. Therefore, the decomposition of E^LPMO_BindE to the sum of electrostatic, exchange-repulsion, and induction terms requires further model calculations, and in particular, for the clusters of many molecules, the analysis is not straightforward.

The CT term for a pair of water molecules is used to identify the H bond. Because the clusters studied in the present study are of polyhedral form, the O–O distance or O–HO distance can be used as the criteria of the H bond, as often used. In the previous studies,²⁶ the CT criteria helped to find the weak H bonds in some complex cage structures of water clusters and to locate the cyclic networks. Subsequently, the hydrogen bonding matrix H_B, defined by Miyake and Aida,⁹ is constructed from the digraph for the cluster.⁴ By manipulating the matrix, the types of the water molecules and of the hydrogen bonds are identified. To analyze the clusters, a code was developed with python3/numpy.

Computational details

The geometries of cubic water clusters (H₂O)₈ were optimized with the MP2/aug-cc-pVDZ level in the present study. For (H₂O)₂₀ we used the thirteen isomers optimized by Xantheas at the MP2/aug-cc-pVDZ level,³⁹ which are among a larger set of isomers determined by empirical potential energy functions. The isomers of tetrakaidecahedral (H₂O)₂₄ are the 20 most stable ones, optimized at the DFT/B3LYP/aug-cc-pVDZ level.⁴⁰ A few additional higher lying stable isomers are selected for comparison.

The calculations were carried out using the in-house ab initio MO package, MOLYX, on a local Mac Pro and on the SGI UV2000 at the RCCS, Okazaki Research facilities of National Institutes of Natural Science (NINS). The required computer resources for (H₂O)₂₄/aug-cc-pVDZ on the SGI UV₂₀₀₀ were 32 cores (16G memory each) and 236 h 42 m cpu (11 h 05 m wall) time for LPMO 3SPT + Disp. The present code is partially parallelized using OpenMP.

Results and discussion

Relative energies

The relative energies of the polyhedral cages of (H₂O)_n isomers for n = 8, 20 and 24 studied in the present work are shown in Fig. 1. Table 1 compares the relative (with respect to the most stable isomer) binding energies of all of the cube (H₂O)₈ isomers evaluated by eqn (4) of LPMO PT, and by the MP2 method. The energy in the third column is evaluated by replacing the HF energy with E_HF(Ψ_LPMO) + E^2SPT + E^3SPT in eqn (3). In previous studies,^23,34,41 the evaluated energy is shown to be close to the counterpoise (CP) corrected HF binding energy if the augmented basis functions are used. By adding the dispersion term E_Disp as in eqn (4), E^3SPT+Disp_BindE comes close to the approximately CP corrected HF plus MP2 binding energy E^{apprCP_HF+MP2}. For instance, E^3SPT+Disp_BindE/aug-cc-pVDZ for Cl⁻(H₂O) and Cl⁻(H₂O)₂ is −62.0 and −125.5 kJ mol⁻¹, respectively.⁴¹ The corresponding E^{apprCP_HF+MP2}/aug-cc-pVxZ binding energies are −57.6 (x = D), −61.2 (x = T), −62.1 (x = Q) and −62.5 (x = 5) kJ mol⁻¹ for Cl⁻(H₂O),⁴² and −112.7 (x = D), −125.3 (x = T), −114.3 (x = Q) kJ mol⁻¹ for Cl⁻(H₂O)₂.⁴¹

Fig. 1 Relative binding energy of polyhedral water clusters of (H₂O)_n. (a) n = 8, (b) n = 20, and (c) n = 24. See text and Table 4 for A, B, C and D in (c).

Table 1 Relative energies (kJ mol⁻¹) of the cubic isomers of (H₂O)₈

Isomer^a	E ^3SPT+DispBindE (eqn (4))	E ^apprCP_HF+MP2BindE	E ^HF+MP2BindE	dda(s,t)^c
a See Fig. 2 for the labelling of isomers. The first number is the number of the dda⇐daa pairs, Mdda⇐daa (≡bFL), and the second is the ordering of the binding energy. b The HF binding energy is replaced with E^LPMOBindE + ECT+LE in eqn (4). c The number of dda(s,t) type waters; (s,t) = (0,3), (1,2), (2,1), (3,0). daa(s,t) = dda(t,s). See text.
4-1	0.0	0.0	0.0	4,0,0,0
4-2	0.32	0.07	0.16	4,0,0,0
2-1	12.51	12.09	12.42	0,4,0,0
2-2	12.76	12.27	12.71	0,4,0,0
2-3	16.65	16.84	17.43	1,2,1,0
2-4	16.66	16.50	17.04	1,2,1,0
2-5	16.89	17.07	17.71	1,2,1,0
2-6	16.99	16.12	16.44	1,2,1,0
1-1	26.37	26.49	27.44	0,2,2,0
1-2	27.01	27.03	28.20	0,2,2,0
1-3	27.10	26.68	27.48	0,2,2,0
1-4	30.93	30.83	31.73	0,3,0,1
0-1	42.07	44.19	46.02	0,0,4,0
0-2	42.10	44.40	46.26	0,0,4,0

---

Table 1 and Fig. 1a clearly show that in terms of the stability there are four groups of the cubic (H₂O)₈ isomers, whose Schlegel diagrams are given in Fig. 2.⁴³ The corresponding ball-and-stick figures are given in ESI.†Fig. 3 shows the correlation between the CT and dispersion (Disp) terms and the O⋯O distance for every pair of the H bonds in all cubic (H₂O)₈ isomers and in the thirteen dodecahedral (H₂O)₂₀ isomers studied. The plots distinguish four types of the H bonds, viz. daa⇐daa, daa⇐dda, dda⇐daa, and dda⇐dda; for each type. The strongest H bonds are of the dda⇐daa type, as previously reported.^{8,10,12,17,44} By proposing the ‘Strong–Weak-Effective-Bond’ (SWEB), Kirov, Fanourgakis and Xantheas classified the H bonds on the surface of polyhedral water clusters to five groups, namely t1d, t1a, c2, c0 and c1a.⁴⁴ Their classification was based on the relative orientations as well as the connectivity to neighbors. As shown in Fig. 3, the present four types correspond to theirs except that our daa⇐dda includes both t1a and c1a. Similarly, Chihaia et al. classified the H bonds in the buckyball water clusters to five types.³ By explicitly counting the types of the neighboring water molecules, Kuo et al. defined the “topological index” ξ for H bonds as a sum of the number of dd neighbors at the H donor and the number of aa neighbors at the H acceptor.^19,21 The index was shown to be related to the H bond distance. Anick further classified the H bonds by introducing the descriptors (t_1P,t_1Q) for the H acceptor and (t_2P,t_2Q) for the H donor; there are 49 sub-groups of H bonds,^8,20 which are more explicitly analyzed below in detail. The binding energy of the clusters is correlated with the types of the H bonds in the cluster.

Fig. 2 Schlegel diagrams of all of cubic (H₂O)₈ isomers. The triangles are dda waters, and the squares are daa waters. See the corresponding ball-and-stick figures in ESI.†

Fig. 3 Correlation of the charge-transfer energy E^X,Y_CT and dispersion E^X,Y_Disp energies per H bond for (a) cubic (H₂O)₈ and (b) dodecahedral (H₂O)₂₀ with the O⋯O bond distance/Å.

The first column of Table 1 is the label of the cubic (H₂O)₈ isomers; the first number of the label is the number of the pairs, M_dda⇐aad, of dda⇐daa, which is the strongest among the four types of the H bonds. The total binding energies are grouped by M_dda⇐daa, which is equivalent to parameter b_FL defined by Anick.⁸ He demonstrated that the total energy of the cube isomers is linear to b_FL.

Among the isomers of M_dda⇐daa = 2, the first two isomers, which have the mirror-imaged isomers, are more stable than the next four isomers. To examine the cause of the difference among the M = 2 isomers, the dda water molecules are further sub-grouped by the types of three nearest neighbor water molecules as dda(s,t) with (s,t) = (0,3), (1,2), (2,1), or (3,0), where (s,t) implies that the water molecule bonds to s water molecules of the dda type and t water molecules of the daa type. The list (0,4,0,0) implies that isomers 2-1 and 2-2 have four dda(1,2). Their Schlegel diagrams show that every (triangle) dda molecule has one triangle dda and two squares daa as the neighbor molecules. The number of dda(2,1) waters is equivalent to a_Th of Anick,⁸ who showed that larger a_Th clusters are less stable. The next four isomers share the common list (1,2,1,0). Among the isomers of M_dda⇐daa = 1, the first three of them have the common list (0,2,2,0), and their binding energies are quite close to each other.

This dependence of the relative stability on the list of the sub-groups of water molecules in the cluster is the consequence of the strong dependence of the H bond strength on the pairs of neighboring molecules. To demonstrate it clearly, some of the pairs of the H bonds are explicitly indicated in Fig. 3.

Interestingly, the fitted function in the early work by McDonald et al.¹ successfully distinguished M_dda⇐daa from that of the other types of pairs, but failed to differentiate among the isomers of the same M_dda⇐daa.

All dodecahedral (H₂O)₂₀ and tetrakaidecahedral (H₂O)₂₄ cage isomers studied are a family of M_dda⇐daa = 7 and 9, respectively, as reported earlier. They are only those of the most stable isomers,^2,20,45,46 therefore, the lists of the sub-groups of water molecules are (4,6,0,0) for (H₂O)₂₀ and (6,6,0,0) for (H₂O)₂₄ in the present study. As mentioned in the Introduction, there is a large number of isomers, which belong to the other lists.

Effect of the next-nearest-neighbor hydrogen bonds

In polyhedral water clusters, there are seven types of H donors and of H acceptors, depending on the types of two neighboring waters.⁸Fig. 4 depicts these seven types of H donors and H acceptors. When a dda water (triangle in Fig. 2 and 4) is a hydrogen donor for a hydrogen bond, this water has an H donating to a neighbor water and an HO accepting from a neighbor water. Because the neighbor waters are either dda or daa, as shown in the first row of Fig. 4, there are four types of the hydrogen donors of dda depending on the types of waters at the other end of the neighbor hydrogen bonds. Symbol d̲d_ddaa_daa implies that the hydrogen donor, d̲da, donates a H atom to a neighbor dda water and accepts a H atom from a neighbor daa water. As shown in the second row of Fig. 4, when daa (square) is the hydrogen donor for a hydrogen bond, this water should have two H atoms accepting from the neighboring waters, and there are three types of the hydrogen donor daa in polyhedral water clusters. Similarly, when a dda water is the hydrogen acceptor in a hydrogen bond, as shown in the third row of Fig. 4, there are three types of hydrogen acceptors, depending on the neighbor waters. On the other hand, when a daa water is a hydrogen acceptor, the fourth row of Fig. 4 shows that there are four types of the hydrogen acceptors. Symbol d_ddad_ddaa̲ implies that the hydrogen acceptor is dda, donating a H atom to two dda waters. Thus, if two neighbor waters of both of the hydrogen donor and acceptor waters are taken into account, there are 49 types of H bonds as shown in Fig. 4.⁸

Fig. 4 Seven types of the hydrogen donors and acceptors. The upper seven hydrogen donors correspond to the column in Tables 2 and 3, and the lower seven hydrogen acceptors correspond to the row.

These 49 sub-groups of the hydrogen bonds are same with those characterized by the descriptors (t_1P,t_1Q) for the H acceptor and (t_2P,t_2Q) for the H donor by Anick.²⁰ The average hydrogen bond length R(O⋯O) and the CT term for each type of the hydrogen bonds are given in Table 2 for (H₂O)₈ and in Table 3 for (H₂O)₂₀. In Table 2, the “topological index”, ξ, of Kuo et al. is given in square brackets.¹⁹ The values in the tables are the averaged ones over all isomers studied; note that for (H₂O)₈ all of the possible polyhedral isomers are studied, but for (H₂O)₂₀ only the most stable thirteen isomers are examined. Therefore, the averaged value for the H bond types for (H₂O)₂₀ is favored for the strong H bonds, and therefore, Table 3 lacks some of the H bond sub-groups found in Table 2. The corresponding table for (H₂O)₂₄ is very similar to Table 3 for (H₂O)₂₀. Because the clusters of (H₂O)₂₄ were optimized at the DFT/B3LYP/aug-cc-pVDZ level, while the clusters of (H₂O)₈ and (H₂O)₂₀ were optimized at the MP2/aug-cc-pVDZ level, the direct numerical comparison is not meaningful.

Table 2 The CT energy (kJ mol⁻¹) and O⋯O distance (Å) for the hydrogen bond types in cube (H₂O)₈ isomers.^a In each pair of H-bond types, the first line is the number of cases of the pairs and the topological index ξ in square brackets. The second and third lines are the averaged CT energy and O⋯O distance with the standard deviation in parentheses

	d̲d daa a dda	d̲d dda daa	d̲d daa a dda	d̲d daa a daa	d̲a dda a dda	d̲a dda a daa	d̲a daa a daa
a See Fig. 4 for the types of hydrogen donors and of acceptors.
d dda d dda a̲	0	0	0	1, [0]	0	1, [1]	1, [0]
				−11.45(0.00)		−19.22(0.00)	−21.18(0.00)
				2.710(0.00)		2.605(0.00)	2.587(0.00)
d daa d dda a̲	0	1, [2]	10, [2]	5, [1]	2, [3]	7, [2]	1, [1]
		−8.21(0.00)	−7.93(0.05)	−9.77(0.86)	−13.93(0.01)	−16.20(1.20)	−19.84(0.00)
		2.778(0.00)	2.799(0.003)	2.744(0.020)	2.665(0.001)	2.635(0.011)	2.597(0.00)
d daa d daa a̲	2, [4]	3, [3]	6, [3]	4, [2]	10, [4]	2, [3]	0
	−5.65(0.00)	−7.09(0.61)	−6.91(0.51)	−8.31(0.02)	−12.54(0.62)	−14.29(0.05)
	2.868(0.000)	2.803(0.020)	2.826(0.017)	2.776(0.001)	2.685(0.010)	2.658(0.002)
d dda a dda a̲	0	1, [1]	5, [1]	18, [0]	4, [2]	5, [1]	1, [0]
		−4.61(0.00)	−4.76(0.05)	−5.46(0.40)	−8.31(0.03)	−9.93(0.88)	−11.99(0.00)
		2.869(0.000)	2.890(0.003)	2.844(0.016)	2.758(0.005)	2.718(0.020)	2.678(0.000)
d dda a daa a̲	1, [3]	3, [2]	5, [2]	1, [1]	3, [3]	1, [2]	0
	−3.14(0.00)	−4.02(0.41)	−4.16(0.35)	−4.99(0.00)	−7.59(0.54)	−8.77(0.00)
	2.985(0.000)	2.910(0.024)	2.921(0.021)	2.868(0.000)	2.778(0.016)	2.745(0.000)
d daa a dda a̲	4, [3]	5, [2]	16, [2]	5, [1]	6, [3]	10, [2]	0
	−3.16(0.02)	−3.93(0.47)	−4.10(0.34)	−4.82(0.15)	−6.87(0.50)	−7.92(0.23)
	2.976(0.007)	2.908(0.025)	2.923(0.021)	2.871(0.005)	2.796(0.016)	2.762(0.006)
d daa a daa a̲	11, [4]	1, [3]	4, [3]	0	2, [4]	0	0
	−2.63(0.22)	−3.09(0.00)	−3.55(0.07)		−6.12(0.00)
	3.028(0.019)	2.965(0.000)	2.965(0.005)		2.826(0.000)

---

Table 3 The CT energy (kJ mol⁻¹) and O⋯O distance (Å) for the hydrogen bond types in dodecahedral (H₂O)₂₀ isomers

H-donor^a	d̲d dda a daa	d̲d daa a dda	d̲d daa adaa	d̲a dda a dda	d̲a dda a daa
H-acceptor^a
a See Fig. 4 for the types of hydrogen donors and of acceptors.
d daa d dda a̲				−16.19	−19.13
				2.627	2.594
d daa d daa a̲			−9.59	−14.14	−16.01
			2.735	2.654	2.627
d dda a dda a̲	−6.01	−6.24	−7.07	−9.98
	2.833	2.838	2.800	2.725
d dda a daa a̲	−5.16	−5.15	−6.16
	2.877	2.885	2.838
d daa a dda a̲	−5.01	−5.05	−5.92
	2.877	2.888	2.841

---

As shown in Table 2, the standard deviations both for the CT terms and O⋯O distances are small. The largest ones in cubic isomers (H₂O)₈ are 0.47 kJ mol⁻¹ and 0.025 Å for d_daaa_ddaa̲⇐d̲d_ddaa_daa (the number of cases is 5). For (H₂O)₂₀, the largest standard deviations are 0.46 kJ mol⁻¹ and 0.008 Å, respectively, for d_daad_ddaa̲⇐d̲d_ddaa_daa (the total number of cases is 26).

Tables 2 and 3 are split to four blocks of dda⇐daa, dda⇐dda, daa⇐daa, and daa⇐dda; which correspond to four symbols in Fig. 3. In each block, smaller ξs (the values in square brackets) are for larger CT terms and shorter H bond distances. The “defect” H bonds, as defined by Smith and Dang,² are those in the upper left and lower right blocks in tables.

The large distribution of the plots for the same types of the H bonds found in Fig. 3 results from the dependence of the H bond strength on the types of two neighbor water molecules.

Table 4 shows the average CT energy and O⋯O distance for four of the strongest types of the H bonded pairs in the (H₂O)₂₄ isomers in Fig. 1(c), where the first 20 isomers are labeled A, B, C, C′ and D. Their stability is grouped by the number of the H bonded pair sub-groups. The most stable isomers (A) have three d_daad_ddaa̲⇐d̲a_ddaa_dda, and the next group (B) has two.

Table 4 Characteristics of tetrakaidecahedral (H₂O)₂₄ isomers

H acceptor H donor			d daa d dda a̲ d̲a dda a daa	d daa d dda a̲ d̲a dda a dda	d daa d daa a̲ d̲a dda a daa	d daa d daa a̲ d̲a dda a dda
CT/kJ mol^−1			−18.21	−21.46	−15.28	−18.23
R(O⋯O)/Å			2.602	2.571	2.638	2.597
(H2O)24	Isomers	Relative E/kJ mol^−1
A	w01–w04	−0.50–0.40	0	3	6	0
B	w05–w10	0.85–1.23	1	2	5	1
C	w11–w13	1.37–1.84	2	1	4	2
C′	w14–w18	1.69–2.14	2	1	4	2
D	w19–w20	2.55–2.64	3	0	3	3

---

Charge-transfer theory of Mulliken

Tables 2–4 clearly show that the strength of an H bond is influenced from the types of the two other neighbor waters of the H donor and H acceptor molecules. To further look for the electronic structure origin of this finding, we adopt the classical model theory of the charge-transfer (electron-donor–acceptor) interaction by Mulliken.^29–31

The wave function for an H bond is given by eqn (1), where, for instance, and . By assuming a weak interaction, the mixing coefficient is approximately given [Assumption A] as

where is the resonance integral and the Coulomb attraction term. The effective ionization energy of the donor orbital, , and the electron affinity of the acceptor orbital, , are assumed to be functions of the net charges, and , of the molecules [Assumption B]. Recently Xu and coworkers computationally demonstrated the linear dependence of the ionization energy and electron affinity on a small fractional charge,⁴⁷ using their double hybrid DFT.⁴⁸ The following linear forms

I^A_n(δ) = I^W₀ + p̃^Wδ (7)

A^D_σ(δ) = A^W₀ + ã^Wδ (8)

are used in Assumption B. The linear coefficients p̃^W and ã^W are positive, because the electron is more tightly bound as the system is more positive. By inserting eqn (7) and (8) into eqn (6) and by assuming the change is small, the coefficient can be expanded as The second line defines the positive parameters, and , whereas the third line defines the characteristic factors of the H acceptor and of the H donor.

Effective net charges for the Frontier orbitals

The simplest hydrogen bond is a water dimer, W⇐W, which can be thought of as a reference system for the chains of H bonds. The net charge δ of both molecules is zero, and thus Ω = 0. The coefficient is This simplest H bond induces the net charge −ζ_W2 on the H donor water molecule and +ζ_W2 on the H acceptor water molecule,

ζ_W2 ≡ (−λ_W–W)² (13)

This change of the net charge influences a new H bond. For instance, suppose that a water molecule is H-bonded to the H donor side of the water dimer as (W_a⇐W_da)⇐W_d. To estimate the strength of the new H bond W_da⇐W_d, the change of the ionization energy of the non-bonding orbital in W_da should be taken into account. Here, the Effective Net Charge for the Frontier Orbital (ENCFO), δ(W_da), on water W_da is assumed to be equal to −ζ_W2, because molecule W_da is the hydrogen donor and the electron acceptor in the water dimer. Accordingly, the ionization energy I^A_n(δ) in eqn (6) for the new hydrogen bond, W_da⇐W_d, is changed as in eqn (7). Thus, the CT coefficient for the new H bond, eqn (6), iswhere δ(W_da) = −ζ_W2. By using eqn (10) and (13) Thus, the ENCFO of molecule W̲_d in the chain is

δ((W_a⇐W_da)⇐W̲_d) = −(−λ_{(Wa⇐Wda)⇐Wd})² (16)

For further analysis, the formulas are simplified by reducing the subscripts. The resonance integral is an exponentially decaying function of R(O⋯O), and the Coulomb term is proportional to the inverse of R(O–O). As Tables 2 and 3 show, R(O–O) varies depending on the pairs of the H donor and acceptor. But for the time being, we adopt an ambitious Assumption C that and By neglecting the higher order of ζ_W2, we have

δ((W_a⇐W_da)⇐W̲_d) = −(ζ_W2 + 2p^Wζ_W2^3/2) (17)

So when molecule W̲_d becomes an H acceptor to a new H bond, its effective net charge for the frontier orbital (ENCFO) is negative, and its ionization energy decreases, and thus it becomes a better electron donor (H acceptor).

The next H bond arrangement examined is ((W_a⇐)²W_dda)⇐W̲_d, which is a part of d_ddad_ddaa̲ as shown in Fig. 4. Because molecule W_dda has two donating OH, two of the σ_OH* orbitals accept the electron from the H acceptor molecules, and δ(W_dda) = −2ζ_W2, if the two σ_OH* orbitals do not interfere with each other (Assumption D). Thus we have

δ(((W_a⇐)²W_dda)⇐W̲_d) ≃ −(ζ_W2 + 4p^Wζ_W2^3/2) (18)

The H acceptor d_ddad_ddaa̲ has two of the above chain as ((W_a⇐)²W_dda)²⇐W̲_d. Again by assuming no interference of the H bonds, we obtain the ENCFO and the characteristic factor for the hydrogen acceptor water of row-1 in Table 2 as

δ(d_ddad_ddaa̲) = δ(((W_a⇐)²W_dda)²⇐W̲_d) ≃ −2(ζ_W2 + 4p^Wζ_W2^3/2)

Ω^A(d_ddad_ddaa̲) = −p^Wδ(d_ddad_ddaa̲) ≃ +2p^W(ζ_W2 + 4p^Wζ_W2^3/2) (row-1)

Another chain of the type W̲_a⇐(W_da⇐W_d) is examined, where the H acceptor of a water dimer, W_da⇐W_d, is an H donor in the new H bond. In eqn (6), and . Therefore, the CT coefficient is

where δ(W_da) = ζ_W2 is used. The ENCFO of the end water of the chain under the above assumptions is

δ(W̲_a⇐(W_da⇐W_d)) ≃ +(ζ_W2 + 2a^Wζ_W2^3/2) (21)

and this chain is found in the rightest chain of four types of the hydrogen acceptors at the bottom of Fig. 4. From eqn (17) and (21), the ENCFO cannot be the net charges of molecules in the chain water trimer, W⇐W⇐W, because the sum of the two terms are not zero, unless p^W = a^W. The ENCFO, therefore, is the hypothetical quantity to estimate the changes of the ionization energy and the electron affinity in Mulliken's CT theory.

For chain W̲_a⇐(W_dd⇒W_a), and . Therefore, δ(W_dd) = −ζ_W2. Thus, for a branched chain W̲_a⇐(W_dda(⇐W_d)(⇒W_a)), which is found at the right half of d_ddaa_ddaa̲ in Fig. 4,

δ(W_dda) ≃ +ζ_W2 − ζ_W2 = 0 (22)

Thus,and

δ(W̲_a⇐(W_dda(⇐W_d)(⇒W_a))) (24)

= (−λ_{W̲a−(Wdda(⇐Wd)(⇒Wd))}) = ζ_W2 (25)

By combining eqn (25) with eqn (18), we have for d_ddaa_ddaa̲,

δ(d_ddaa_ddaa̲) = −(ζ_W2 + 4p^Wζ_W2^3/2) + ζ_W2 (26)

  ≃−4p^Wζ_W2^3/2 (27)

Ω^A(d_ddaa_ddaa̲) ≃ +4(p^W)²ζ_W2^3/2 (row-4) (28)

By repeating the above procedure, we have

Ω^A(d_daad_ddaa̲) ≃ +2p^W(ζ_W2 + 2p^Wζ_W2^3/2) (row-2)

Ω^A(d_daad_daaa̲) ≃ 2p^Wζ_W2 (row-3)

Ω^A(d_ddaa_daaa̲) ≃ 4p^W(p^W − a^W)ζ_W2^3/2 (row-5)

Ω^A(d_daaa_ddaa̲) ≃ 0 (row-6)

Ω^A(d_daaa_daaa̲) ≃ −4p^Wa^Wζ_W2^3/2 (row-7)

In the study on the cubic isomers of X⁻(H₂O)₇ (X = F and Cl),⁴¹ the ENCFO of the ligand water molecule (L), which is also induced by the chain of the hydrogen bonds, was derived. The procedure to obtain the formulas for was same with the above one, by replacing “L” in X⁻(H₂O)₇⁴¹ with “D”. The characteristic factors of the H donor in Tables 2 and 3 are

Ω^D(d̲a_daaa_daa) ≃ 2a^W(ζ_W2 + 4a^Wζ_W2^3/2) (column-7)

Ω^D(d̲a_ddaa_daa) ≃ 2a^W(ζ_W2 + 2a^Wζ_W2^3/2) (column-6)

Ω^D(d̲a_ddaa_dda) ≃ 2a^Wζ_W2 (column-5)

Ω^D(d̲d_daaa_daa) ≃ 4(a^W)²ζ_W2^3/2 (column-4)

Ω^D(d̲d_daaa_dda) ≃ 0 (column-3)

Ω^D(d̲d_ddaa_daa) ≃ −4a^W(p^W − a^W)ζ_W2^3/2 (column-2)

Ω^D(d̲d_ddaa_dda) ≃ −4a^Wp^Wζ_W2^3/2 (column-1)

The charge-transfer term in the Mulliken theory is

The relations

Ω^A for upper row > Ω^A for lower row (31)

Ω^D for right column > Ω^D for left column (32)

hold if p^W ≳ a^W > 0. The CT terms shown in Tables 2 and 3 are consistent with these relations for Ω^A and Ω^D with a very few exceptions between columns 3 and 2, and between rows 5 and 6 in Table 2.

The large step between the fifth column (d̲a_ddaa_dda) and the fourth column (d̲d_daaa_daa) corresponds to the well-known empirical rule that the daa water is a better H donor than the dda molecule. The step between the third row (d_daad_ddaa̲) and the fourth row (d_ddaa_ddaa̲) also corresponds to the rule that dda is a better H acceptor. The four blocks in Tables 2 and 3 result from the linear term ζ_W2 in Ω^A and Ω^D.

The sum is a measure of the strength of the hydrogen bonds. For instance, among the dda⇐daa bonds, the sum for the strongest bond is

Ω^A(d_ddad_ddaa̲) + Ω^D(d̲a_daaa_daa) = 2(p^W + a^W)ζ_W2 + 8((p^W)² + (a^W)²)ζ_W2^3/2

while the sum for the weakest one is

Ω^A(d_daad_daaa̲) + Ω^D(d̲a_ddaa_dda) = 2(p^W + a^W)ζ_W2

As is seen in Table 2, the hydrogen bond length of d_ddad_ddaa̲⇐d̲a_daaa_daa is shorter by 0.1 Å than that of d_daad_daaa̲⇐d̲a_ddaa_dda, and the difference in the CT term is as large as 8.64 kJ mol⁻¹. Similarly, the dependence of the H bond strength on Kuo's topological index ξ within each of four blocks is elucidated by the sum ; the coefficient of ζ_W2^3/2 is related with the topological index ξ.

To further demonstrate the relation between the sum of the characteristic factors, , and the CT terms and O⋯O distances, a matrix is constructed by extracting the columns and rows of Tables 2 and 3, and is shown in Fig. 5. Because all the parameters are safely assumed to be positive and ζ < 1, the large CT terms and short O⋯O distances strongly are correlated with the sum. If p^W ≈ a^W, the matrix of becomes approximately symmetric, and Fig. 5 numerically shows that it is nearly symmetric.

Fig. 5 Comparison of the sum of the characteristic factors, , for some of the typical hydrogen donor and acceptor pairs with the CT term/kJ mol⁻¹ (the upper value) and the O⋯O distance per Å. The left and right values are for (H₂O)₈ and for (H₂O)₂₀, respectively. The subscripts and superscripts in the equations in text are deleted.

The differences of the characteristic factors for columns 7 and 6, and for columns 6 and 5 in Table 2 are . The corresponding numerical differences are [1.95, 3.64, *, 2.06, *, *] and [*, 2.27, 1.75, 1.62, 1.18, 1.05] kJ mol⁻¹, respectively, from the top row, where “*” implies that no data are available in Table 2. The averaged value (the standard deviation) of these two differences is 1.94 (0.80) kJ mol⁻¹. Similarly, the difference of the characteristic factors for rows 2 and 1 and for rows 3 and 2 is . The corresponding numerical differences in Table 2 are [1.34, 3.02, *, 1.68, *, *] and [*, 1.91, 1.35, 1.46, 1.03, 1.12] kJ mol⁻¹, respectively, from the rightest of rows. The average value (standard deviation) of two differences is 1.61 (0.63) kJ mol⁻¹. The other examples are seen in the differences of the characteristic factors for columns 3 and 1 and for rows 7 and 6; both of them are . The corresponding numerical difference values in Table 2 are [*, *, 1.25, *, 0.88, 0.94, 0.92] kJ mol⁻¹ between columns 3 and 1 from the top, and [*, *, 0.75, *, 0.55, 0.84, 0.53] kJ mol⁻¹ between rows 7 and 6 from the right. The average value (standard deviation) of two differences is 0.83 (0.23) kJ mol⁻¹. Although the standard deviations are not small, we may safely claim that these numerical tests qualitatively support the model derived from Mulliken's CT theory. But, if the product are assumed to be common for all of the above cases, the ratio a^W/p^W could be estimated as

which shows that the simple assumption is broken down. The refinement of the model is required to include the influence from the change of the H bond distance.

In Fig. 5, the values for (H₂O)₈ (left numbers in each box) and (H₂O)₂₀ (right) are compared with each other. The strengthening of the hydrogen bonds in (H₂O)₂₀ with respect to (H₂O)₈ is evident. One possible cause is the size of the ring structure; five membered rings in (H₂O)₂₀ and four membered rings in (H₂O)₈, resulting in less strain in the ring structure for (H₂O)₂₀ than for (H₂O)₈. However, the above analysis based on Mulliken's CT theory strongly supports an alternative cause, namely that the longer lengths of the hydrogen bond chains in (H₂O)₂₀ influence the hydrogen bond strength. As shown above, the change in ENCFO propagates positively or negatively through the long hydrogen bonded chain. In his regression analysis,²¹ Anick added the parameters of the contribution from the ternary water molecules (next-nearest neighbors), and Recently Akase and Aida examined the type and size dependence of the dipole moment of a water molecule in various geometric configurations of water clusters (H₂O)_n, n = 2–10,¹⁵ using the all-atom polarizable model potential, TTM2-R.⁴⁹ The dipole moments of dda and daa waters do slowly converge up to n = 8. This finding supports that the difference between (H₂O)₈ and (H₂O)₂₀ in Tables 2 and 3 and in Fig. 5 results from the size dependence of the CT interaction and of the hydrogen bond strength. Similarly Bakó and Mayer reported the collective effects on the dipole moments caused by “non-immediate environment” in water clusters.⁵⁰ It suggests that the regression analysis for a set of the water clusters should have better carried out with a constant size, to avoid the unintentional statistical errors,^1,8,12,20,51 unless the lengths of the H bonded networks are taken into account.

Concluding remarks

The influence of the neighbor and next-neighbor water molecules on the strength of the H bonds was examined for the polyhedral clusters of (H₂O)₈, (H₂O)₂₀ and (H₂O)₂₄. The polyhedral water clusters are suitable for the study, because they have only two types of molecules, either dda type or daa type, and they are further sub-grouped by the types of three neighbor (H-bonded) molecules as dda(s,t) and daa(s,t), s + t = 3. The stability of the isomers is related with the numbers of the eight sub-grouped molecules. In these clusters there are four types of the H bonds, and they are further divided to 49 sub-types of the H bonds when the neighbors of each of the H donor and H acceptor water molecules are taken into account. The H bond strengths are classified by these 49 sub-types. Based on the Mulliken's CT theory, the formulas for the CT contribution to the H bond energy are derived, and they explain the sub-type dependence of the H bond strength. This model theory can be extended to any H bonded networks of water molecules, and the works are in progress for hundreds of isomers of (H₂O)₈ and (H₂O)₁₂, which include various H bonded networks.⁵² The model theory with properly determined parameters could help to improve the empirical potential functions containing the influence of the H bonded chains and networks.

Acknowledgements

The work is partially supported by the Grant-in-Aid for Science Research of JSPS (S.I.). Part of this work was supported by the US Department of Energy, Office of Science, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences and Biosciences at Pacific Northwest National Laboratory. Battelle operates the Pacific Northwest National Laboratory for the US Department of Energy (S. S. X.). S. I. acknowledges Prof. H. Matsuzawa and Dr C. Ishibashi for their helpful discussion on the H bond analysis.

References

1. S. McDonald, L. Ojamae and S. J. Singer, J. Phys. Chem. A, 1998, 102, 2824–2832.
2. D. E. Smith and L. X. Dang, J. Chem. Phys., 1994, 101, 7873–7881.
3. V. Chihaia, S. Adams and W. F. Kuhs, Chem. Phys., 2004, 297, 271–287.
4. W. C. Radhakrishnan and T. P. Herndon, J. Phys. Chem, 1991, 95, 10609–10617.
5. J.-L. Kuo, J. V. Coe, S. J. Singer, Y. B. Band and O. Lars, J. Chem. Phys., 2001, 114, 2527–2540.
6. M. V. Kirov, J. Struct. Chem., 2002, 43, 790–797.
7. D. J. Anick, THEOCHEM, 2002, 587, 87–96.
8. D. J. Anick, THEOCHEM, 2002, 587, 97–110.
9. T. Miyake and M. Aida, Chem. Phys. Lett., 2006, 427, 215–220.
10. U. Buck, I. Ettischer, M. Melzer, V. Buck and J. Sadlej, Phys. Rev. Lett., 1998, 80, 2578–2581.
11. J. Brudermann, M. Melzer, U. Buch, K. Kazimirski, J. Sadlej and V. Bush, J. Chem. Phys., 1999, 110, 10649–10652.
12. J. Sadlej, Chem. Phys. Lett., 2001, 333, 485–492.
13. M. Miyazaki, A. Fujii, T. Ebata and N. Mikami, Science, 2004, 304, 1134–1137.
14. J.-W. Shin, N. I. Hammer, E. Diken, M. A. Johnson, R. S. Walters, T. D. Jaeger, M. Duncan, R. A. Christle and K. D. Jordan, Science, 2004, 304, 1137–1140.
15. D. Akase and M. Aida, J. Phys. Chem. A, 2014, 118, 7911–7924.
16. S. S. Xantheas, Chem. Phys., 2000, 258, 225–231.
17. K. Ohno, M. Okimura, N. Akai and Y. Katsumoto, Phys. Chem. Chem. Phys., 2005, 7, 3005–3014.
18. S. Yoo and S. S. Xantheas, in Handbook of Computational Chemistry, ed. J. Leszczynski, Springer, 2012, pp. 761–792.
19. J.-L. Kuo, C. V. Ciobanau, L. Ojamae, I. Shavitt and S. J. Singer, J. Chem. Phys., 2003, 118, 3583–3588.
20. D. J. Anick, J. Chem. Phys., 2003, 119, 12442–12455.
21. D. J. Anick, J. Phys. Chem. A, 2006, 110, 5135–5143.
22. S. Iwata and T. Nagata, Theor. Chem. Acc., 2006, 116, 137–144.
23. S. Iwata, J. Chem. Phys., 2011, 135, 0941011.
24. T. Nagata and S. Iwata, J. Chem. Phys., 2004, 120, 3555–3563.
25. S. Iwata, Phys. Chem. Chem. Phys., 2012, 14, 7787–7794.
26. S. Iwata, P. Bandyopadhyay and S. S. Xantheas, J. Phys. Chem. A, 2013, 117, 6641–6651.
27. S. Iwata, Phys. Chem. Chem. Phys., 2014, 16, 11310–11317.
28. R. S. Mulliken and W. B. Person, Molecular Complexes, Wiley Interscience, NY, 1969.
29. R. S. Mulliken, J. Am. Chem. Soc., 1952, 74, 811–824.
30. K. Nukazawa, J. Tanaka and S. Nagakura, J. Phys. Soc. Jpn., 1953, 8, 792–793.
31. H. Tsubomura, Bull. Chem. Soc. Jpn., 1954, 27, 445–450.
32. S. Iwata, Chem. Phys. Lett., 2006, 431, 204–209.
33. S. Iwata, J. Phys. Chem. B, 2008, 112, 16104–16109.
34. S. Iwata, J. Phys. Chem. A, 2010, 114, 8697–8704.
35. I. G. Kaplan, in Types of Intermolecular Interactions: Qualitative Picture, ed. S. Wilson, Wiley Interscience, NY, 2003, ch. 6, vol. 3, pp. 161–181.
36. I. G. Kaplan, Intermolecular Interactions, John Wiley & Sons, Ltd., 2006.
37. K. Jeziorski and B. Szalewicz, in Symmetry-adaptive Perturbation Theory, ed. S. Wilson, Wiley Interscience, NY, 2003, ch. 9, vol. 3, pp. 232–277.
38. C. J. Burnham and S. S. Xantheas, J. Chem. Phys., 2002, 116, 1500–1510.
39. S. S. Xantheas, Can. J. Chem. Eng., 2012, 90, 843–853.
40. S. Yoo, M. V. Kilov and S. S. Xantheas, J. Am. Chem. Soc., 2009, 131, 7564–7566.
41. C. Ishibashi, S. Iwata, K. Onoe and H. Matsuzawa, J. Phys. Chem. A, 2015, 119, 10241–10253,  DOI:10.1021/acs.jpca.5b07244.
42. S. S. Xantheas, J. Chem. Phys., 1996, 104, 8821–8824.
43. V. Schlegel, Verhandlungen der Kaiserlichen Leopoldinisch-Carolinischen Deutschen Akademie der Naturforscher, 1883, 44, 343.
44. M. V. Kirov, G. S. Fanourgakis and S. S. Xantheas, Chem. Phys. Lett., 2008, 461, 180–188.
45. L. Lenz and A. Ojamae, Chem. Phys. Lett., 2006, 418, 361–367.
46. D. J. Anick, J. Chem. Phys., 2010, 132, 164311.
47. I. Y. Zhang and X. Xu, Phys. Chem. Chem. Phys., 2012, 14, 12554–12570.
48. I. Y. Zhang and X. Xu, Int. Rev. Phys. Chem., 2011, 30, 115–160.
49. S. S. Xantheas, C. J. Burnham and R. J. Harrison, J. Chem. Phys., 2002, 116, 1493–1499.
50. I. Bakó and I. Mayer, J. Phys. Chem. A, 2016, 120, 631–638.
51. L. Lenz and A. Ojamae, Phys. Chem. Chem. Phys., 2005, 07, 1905–1911.
52. D. Akase, M. Aida, K. Ohno and S. Iwata, PacifiChem, 2015.

---

Footnotes

† Electronic supplementary information (ESI) available: The xyz files (cubic_W8.xyz.zip) the ball-and-stick figures for the cubic (H₂O)₈ are available. See DOI: 10.1039/c6cp02487d

‡ Visiting Professor, 2012–2016.

---