Hydrogen-bond pattern to characterize water network

Graph and digraph

A graph is a set of vertices and edges. Vertices are often called nodes or points, and edges are also called lines. The edge donates the connectivity between two vertices. Vertices and edges may be distinguishable by labels. Graphs are widely used for schematic representations (Fig. 1). In chemistry, a chemical formula is described as a graph whose vertices represent atoms and edges represent chemical bonds: labels of vertices donate elements symbols and labels of edges mean kinds of bonds such as single bonds, double bonds, aromatic bonds, and so on.

Fig. 1:

A graph representation of cyclopropane and the corresponding adjacency matrix.

Graphs can be represented by various ways. One convenient way is to use an adjacency matrix. A graph in which edges are not directed is called an undirected graph. The adjacency matrix of an undirected graph with n vertices is the n^th order square matrix, whose element a_ij is equal to 1, if there is an edge between vertex i and vertex jd, and 0 otherwise. Hence the matrix is symmetric.

A graph is a mathematical structure which represents the topology of a molecule. Molecular structures can be represented with graphs, and each of different graphs corresponds to a different adjacency matrix. By imposing the conditions which need to be satisfied, all possible adjacency matrices can be created, and then all possible graphs can be obtained, thus all possible topologically distinct structures can be obtained. In the case of hydrocarbon, the graph must satisfy the following two conditions: (i) the maximum number of the edges connected to a vertex is 4, (ii) the graph is connected. All possible isomers with three carbon atoms can be obtained by counting up all possible graphs of n=3 with the imposed conditions satisfied.

A graph invariant is a property of graph. The number of vertices and the number of edges are two simple graph invariants. The Hosoya index [4] is also a kind of graph invariant. If two graphs are equivalent, the graphs are called isomorphic. Isomorphic graphs have the same graph invariants.

The most important feature of hydrogen bond is that it possesses the direction. An HB network can be represented by a digraph. A digraph (or directed graph) is a graph whose edges are directed (Fig. 2). In this case, vertices represent components forming hydrogen bonds, and arcs stand for hydrogen bonds from hydrogen donors to accepters.

Fig. 2:

A graph representation of water trimer and the corresponding hydrogen bond matrix.

A digraph contains more information than a graph, concerning the relation between two vertices, i.e. the direction from a vertex toward another vertex. A directed edge is called an arc. An HB network indicates how water molecules are connected. Since a water molecule possesses two hydrogen donors and two hydrogen acceptors, it can form up to four hydrogen bonds.

Radhakrishnan and Herndon [12] in 1991 used digraphs to represent hydrogen-bonded water clusters, and the interaction energy of water clusters were fitted as a function of the graph invariants. The correlation was very good if proper graph invariants were selected for the analyses. The selected graph invariants include not only the number of arcs, which represents the number of the hydrogen bonds formed in the cluster, but also information about the network, e.g. the numbers of paths of different lengths, the number of three- and four-membered rings, and the total number of double donation or acceptance of arcs at each vertex. This result indicates that the cooperativity of hydrogen bonds is inherently associated with the HB network.

Graph theoretical techniques were used to generate stable arrangements of water clusters by McDonald et al. [13] in 1998. Graph invariants were used to present the features of water clusters by Kuo et al. [14] in 2001, based on the bond variables b_ij, in which b_ij is equal to 1, if water i donates to water j, −1 if water j donates to water i, and 0 if there is no hydrogen bond between i and j. Physical properties of water clusters are fit to some graph invariants [14].

Digraph and HB matrix, and enumeration of HB patterns

An important aspect of the matrix representation of digraph is that it allows efficient enumeration of water clusters. In order to enumerate all possible digraphs corresponding to water clusters and also to distinguish topologically distinct isomers, a special adjacency matrix was defined by Miyake and Aida [15], where the element h_ij of the adjacency matrix is 1 if, and only if, the arc from vertex i toward vertex j exists. This matrix for digraph was named “hydrogen bond (HB) matrix” (Fig. 2).

A digraph satisfying the following three conditions is considered as representing a water cluster: (i) maximum number of arcs directed toward a vertex is 2 and maximum number of arcs directed from a vertex is 2, (ii) up to one arc is allowed between any two vertices, (iii) the digraph must be connected. A digraph which satisfies these conditions corresponds to HB network. Miyake and Aida [15], [16] enumerated the number of the digraph satisfying the above conditions by means of adjacency and HB matrices up to eight vertices. Vukičević et al. [17] extended up to 12 vertices using a more sophisticated algorithm. Then, a different algorithm was created to generate up to 16 vertices by Brinkmann [18]. The numbers guarantee the upper limits of the HB patterns of water clusters; in other words, there are no other topologically distinguishable isomers.

All possible graphs and digraphs which correspond to water clusters, for n=3–5, are given in [16]. Here, all possible graphs and all possible digraphs with three vertices are shown in Fig. 3. The number of the graphs is two (one linear and one cyclic). The graphs with n=3 correspond to all possible isomers of saturated hydrocarbon with three carbon atoms. Actually, the linear isomer is propane and the cyclic isomer is cyclopropane. The digraphs with n=3 correspond to water trimer, and the number of those digraphs is 5 (3 linear and 2 cyclic). The linear digraphs are created from the linear graph. The cyclic digraphs are created from the cyclic graph. The important point here is that we know all the possible topologically distinct HB patterns for water trimer: there are no HB patterns for water trimer other than the digraphs shown in Fig. 3.

Fig. 3:

All possible graphs and digraphs with three vertices.

Rooted digraph and HB matrix

A rooted graph is a graph in which one vertex is different from the other vertices. Molecular structures with a heteroatom can be represented with rooted graphs, where vertices correspond to a heteroatom or carbon atoms, and edges correspond to bonds. Rooted digraph is used to represent the feature of protonated water (PW) cluster, since PW is distinguished from the other water molecules in the PW cluster [19], [20]. A rooted digraph is a set of vertices and arrows (Fig. 4). The HB matrix can be applied to a rooted digraph: for a rooted digraph with n vertices, the HB matrix H is the n^th order square matrix, whose element h_ij is equal to 1 for an arrow directed from vertex i to vertex j, and 0 otherwise.

Fig. 4:

A graph representation of protonated water trimer and the corresponding hydrogen bond matrix.

The conditions for a rooted digraph to be equivalent to a PW cluster are the following. In a PW cluster, one PW can accept up to one proton from other water molecules, and donate up to three protons to other water molecules. One water molecule can accept up to two protons from other water molecules or PW, and donate up to two protons to other water molecules or PW. In addition, in a PW cluster, all of the neutral water molecules and PW should be connected by hydrogen bonds.

Generation of water cluster structures

TTM2-R potential energy function [23] was used to calculate interaction energies of water clusters. The permanent dipole moment of this model originating from partial charges is 1.853 D, which reproduces the experimental value of 1.8546 D [24] in the gas phase. In addition, atomic polarizabilities placed on O- and H-sites produce induced dipole moments influenced by surrounding molecules. The repulsive exponential term acting between the O-sites to the two-body pairwise terms was added [22] in order to correct the unphysical stabilization at a short O···O distance of the original TTM2-R function. The modified two-body pairwise terms are

Vpair=∑i=1N−1∑j>iN[Ar12+Br10+Cr6+De−Er].

The values of parameters A, B, and C are taken from the original paper [23]. The parameters for the added terms are D=1.0×10¹³ kcal/mol and E=13.2 Å⁻¹ [22].

The MC simulations of water clusters, (H₂O)_n, n=2–10, were performed with the standard Metropolis algorithm [25]. The details concerning the simulation procedure are given in [22]. Configurations were sampled at every MC step, and 10⁹ configurations were classified into HB patterns by solving the graph isomorphism problem for each system. We defined the formation of hydrogen bonding, if the hydrogen-bond distance r(OH···O) is less than 2.45 Å and at the same time the hydrogen-bond angle θ(O–H···O) is wider than 120°.

The relative Helmholtz energy of each HB pattern is evaluated from the population of HB patterns, because the MC simulations were performed in the canonical ensemble.

ΔAIJ=−RTlnNINJ.

Here, ΔA_IJ is the Helmholtz energy of HB pattern I relative to HB pattern J. R is the gas constant and T is temperature. N_I denotes the number of the configurations classified into the HB pattern I.

For each HB pattern, the energy average 〈E〉_I is given by

〈E〉I=∑i∈IViNI,

where the sum is over the configurations belonging to the HB pattern I, and V_i is the potential energy associated with a configuration i. The internal energy of HB pattern I relative to HB pattern J is obtained by

ΔUIJ=〈E〉I−〈E〉J.

The relative Helmholtz energy between HB patterns I and J is related to the relative internal energy and entropy,

ΔAIJ=ΔUIJ−TΔSIJ,

where ΔS_IJ is the relative entropy between the HB patterns I and J. The intramolecular vibrational entropic contributions of water molecules are not included, because the rigid water model is used. Note that the internal energy and the entropy of an HB pattern are implicitly temperature-dependent.

Distribution of HB patterns of water trimer

From the MC simulations of water clusters (H₂O)n, n=3–8, the relative values of Helmholtz energy (ΔA), internal energy (ΔU), and entropic term (−TΔS) for several HB patterns of each n-water system were evaluated, and the results of some of low Helmholtz energy HB patterns for each system are given in Ref. [22].

Here, the results for water trimer are summarized. As shown in Fig. 3, the number of topologically distinct patterns of the system with three vertices is 5; i.e. it is certain that there are only five digraphs (HB patterns) for water trimer. The relative and cumulative populations of HB patterns of water trimer are shown in Ref. [22]. At 200 K, a cyclic HB pattern 3-3A is the most populated, namely, the lowest in Helmholtz energy. The Helmholtz energy (blue bar), the internal energy (red bar) and the entropic term (green bar) of other HB patterns relative to the corresponding values of 3-3A at 200 K are shown in Fig. 6. The Helmholtz energy (in blue) of 3-2A is 2.9 kJ/mol higher than that of 3-3A, followed by 3-2B, 3-2C and 3-3B.

Fig. 6:

Helmholtz energies, internal energies, entropic terms relative to the most populated HB pattern for water trimer at 200 and 300 K.

At 300 K, on the other hand, the linear pattern 3-2A is the most populated, lying 1.0 kJ/mol lower than the cyclic pattern 3-3A in Helmholtz energy. The Helmholtz energy, the internal energy and the entropic term of each of other HB patterns relative to the corresponding values of 3-2A at 300 K are shown also in Fig. 6. The internal energy of 3-2A is 12.0 kJ/mol higher than that of 3-3A which remains the lowest among the five patterns even at 300 K. The entropic term demonstrates that the cyclic patterns (3-3A and 3-3B) are entropically less favorable than the linear patterns (3-2A, 3-2B, and 3-2C).

All of the global minimum structure and several stationary points of water trimer, which were obtained by means of high-level ab initio molecular orbital (MO) calculations [26], [27], are identified in the configuration space of 3-3A. We have calculated the local minimum geometries of water trimer at the theoretical level of MP2/aug-cc-pVDZ using Gaussian 09 program package [28], and have found four local minima. They are shown in Fig. 7. The lowest two geometries correspond to 3-3A. In the HB pattern, the orientation of dangling OH groups is not distinguished. The other two geometries correspond to 3-2C and 3-2B, respectively. Local minima do not exist corresponding to either 3-2A or 3-3B. It is worthy to note that there is no local minimum on the potential energy surface of water trimer which corresponds to the linear 3-2A, which is the most popular HB pattern at 300 K. This is a very important finding; local minimum geometries cannot be always representative.

Fig. 7:

Local minimum structures of water trimer with potential energy (ΔE) and Gibbs energy (ΔG) at 298.15 K relative to those of the global minimum structure, and the HB patterns of water trimer with internal energy (ΔU) and Helmholtz energy (ΔA) relative to those of 3-3A at 300 K.

The internal energies of the HB patterns 3-2C and 3-2B relative to that of 3-3A at 300 K are 23.5 and 23.4 kJ/mol, respectively, and they are similar to but lower than the relative energies of the corresponding local minima (28.8 and 29.1 kJ/mol, respectively, at MP2/aug-cc-pVDZ). This is fairly reasonable and it may indicate the reliability of the HB pattern analysis based on the MC simulation using TTM2-R potential energy function [23].

Relative Gibbs energy ΔG of each local minimum structure in Fig. 7 is a value relative to Gibbs energy of the global minimum geometry. Gibbs energy of each local minimum structure has been calculated at 298.15 K, at the theoretical level of MP2/aug-cc-pVDZ using Gaussian 09 program package [28]. On the other hand, relative Helmholtz energy ΔA between each HB pattern and the most populated HB pattern in Fig. 7 is a value evaluated from MC simulation at 300 K [22]. Since no local minima are found corresponding to either the linear HB pattern 3-2A or the cyclic HB pattern 3-3B, we cannot estimate ΔG of a structure corresponding to either HB pattern.

Quantum chemical evaluation of ΔG at 298.15 K of a molecular system is often performed based on the left method in Fig. 5. It must be noticed that other isomers which do not correspond to local minima on the potential energy surface of a molecular system can be observed at a finite temperature. On the other hand, the results of molecular simulation depend severely on the quality of the potential energy function used in the simulation. In the current case [22], the quality of the TTM2-R potential energy function [23] is reasonably high as described above and also later.

Distribution of HB patterns of water octamer

Here we show another example that an isomer which corresponds to the global minimum on the potential energy surface may not be observed at a finite temperature.

The number of topologically distinct patterns of the system with eight vertices is 164 461; i.e. it is certain that there are 164 461 digraphs (HB patterns) for water octamer at the most [15]. From the MC simulations of water clusters, the results for water octamer are shown here. The relative and cumulative populations of HB patterns of water octamer are given in ref [22].

At 200 K, 18 589 HB patterns were sampled in the MC simulation. At 200 K, one of the cubic HB pattern (S₄ symmetry) is the most populated, namely, the lowest in Helmholtz energy. The Helmholtz energy (blue bar), the internal energy (red bar) and the entropic term (green bar) of other HB patterns relative to the corresponding values of the S₄ cubic HB pattern at 200 K are shown in Fig. 8, where two cubic HB patterns and three pseudo-cubic HB patterns in which one hydrogen bond is broken are included. In each of these pseudo-cubic HB patterns, the internal energy becomes higher, while the entropic term becomes lower than the corresponding cubic pattern. In Fig. 8, various HB patterns other than cubic pattern are included. In those patterns, the internal energy becomes higher, while the entropic term becomes lower than the cubic pattern. The spirocyclic HB pattern shows the lowest entropic term in the HB patterns in Fig. 8, while it is the seventh most stable in Helmholtz energy at 200 K.

Fig. 8:

Helmholtz energies, internal energies, entropic terms relative to the most populated HB pattern for the 10 most populated HB patterns of water octarmer at 200 K.

At 300 K, 37 978 HB patterns were obtained in the MC simulation. At 300 K, the spirocyclic HB pattern, which is the lowest in entropy term at 200 K in Fig. 8, is the most populated, namely, the lowest in Helmholtz energy. The Helmholtz energy, the internal energy and the entropic term of each of other HB patterns relative to the corresponding values of the spirocyclic isomer at 300 K are shown in Fig. 9, where many other patterns are included. In those patterns, the entropic terms are further lower. It is worthy to note that the populations of cubic isomers are very small at 300 K.

Fig. 9:

Helmholtz energies, internal energies, entropic terms relative to the most populated HB pattern for the 10 most populated HB patterns of water octarmer at 300 K.

For the cubic water octamers, there are 14 HB patterns, and 14 local minimum structures are found at the theoretical level of MP2/aug-cc-pVDZ. One local minimum is found for each HB pattern, and they are shown in Fig. 10. The symmetry of each structure (1–14) in Fig. 10 is as follows: 1 S₄, 2 D_2d, 3 C_i, 4 C₂, 8 C_s, 13 C₄ and 14 C₄ symmetry. Other structures in Fig. 10 are in C₁ symmetry. Relative Gibbs energy ΔG of each local minimum structure is plotted in Fig. 11, which is a value relative to Gibbs energy of the structure 1. Gibbs energy of each local minimum structure has been calculated at 298.15 K, at the theoretical level of MP2/aug-cc-pVDZ using Gaussian 09 program package [28].

Fig. 10:

Cubic water octamers.

Fig. 11:

Relative free energies of HB patterns and relative free energies of local minimum structures of cubic water octamers.

Relative Helmholtz energy ΔA between each HB pattern and the most populated HB pattern 1 is plotted in Fig. 11. We have evaluated the relative Helmholtz energies of the cubic HB patterns from the MC simulation at 300 K [22] using TTM2-R potential energy function [23]. Note that the structures corresponding to HB patterns 13 and 14 were not sampled in the MC simulation, which indicates that the energies of those HB patterns are too high to be observed at 300 K.

As shown in Fig. 11, the order of ΔG of the cubic structures and the order of ΔA of the cubic HB patterns are similar: it may indicate the reliability of the HB pattern analysis based on the MC simulation using TTM2-R potential energy function [23].

The cubic structure 2 (D_2d symmetry) is the global minimum on the potential energy surface of water octamer. The cubic structure 1 (S₄ symmetry) is the most populated in a low temperature. At higher temperature, entropically unfavorable cubic structures are less populated. The populations of cubic HB patterns are very small at 300 K: any cubic structures may not be observed at 300 K. Even though ΔG of the cubic structure 1 is the lowest among the 14 cubic structures at 298.15 K, it does not mean the cubic structure 1 is the most populated among water octamer at 298.15 K.

The relative stability in the 14 cubic water octamers is related to the HB patterns, and is attributed to the strong dependence of the HB strength on the local HB types (see Section “HB type of water molecule” and Fig. 12) of the contributing water molecules [9]. The hydrogen bond between the water molecules of the type dda and the type daa (i.e. dda water molecule is proton acceptor and daa water molecule is proton donor, dda←daa) is the strongest; the origin of the strong HB energy is the large charge transfer energy as well as the large dispersion energy [9]. The number of the dda←daa pair in each HB pattern (1–14) in Fig. 10 is as follows: 4 for HB patterns 1 and 2, 2 for 3–8, 1 for 9–12, 0 for 13 and 14.

Fig. 12:

Molecular local HB type and dipole moment of water molecule estimated from MC simulation of water decamer.

Local HB type of water molecule

The constituent water molecules in water cluster can be classified into the local HB types (free, d, a, dd, aa, da, dda, daa, and ddaa) (Fig. 12) according to the way how the water molecule participates in the HB network [22].

Although the TTM2-R model [23] has the fixed geometry of a water molecule, a dipole moment of a constituent water molecule of a water cluster changes because of the induced dipole moments which are generated by the interaction with surrounding molecules.

Molecular dipole moments of constituent water molecules

The ensemble averages of molecular dipole moments of water molecules in water clusters (in the MC simulation [22] at 300 K) for each system size (the number of water molecules in the system) are given in [22]. The dipole moments of hydrogen-bonded water molecules increase with the size of the water cluster. The molecular dipole moment depends on the local HB network around the water molecule [22]. The magnitude of the molecular dipole moment increases in the following order: free<d<a<dd<aa<da<dda<daa<ddaa. The average dipole moment for each local HB type of water decamer is shown in Fig. 12. The dipole moments of not-hydrogen-bonded (free) molecules are almost constant independent of the number of water molecules in the system, whereas the dipole moments of hydrogen-bonded molecules increase with the number of water molecules in the system.

The dipole moments of free molecules are close to that of an isolated water molecule (1.853 D). The dipole moment of a water molecule is enhanced by forming hydrogen bonds, and the enhancement depends on the type of the local HB environment. The dipole moments of da type molecules are larger than those of aa and dd types. It may indicate that a water molecule which accepts a hydrogen bond and at the same time donates a hydrogen bond is polarized by the continuous hydrogen bonds. The dipole moments of water molecules of dda and daa types are larger than those of da type molecules. The dipole moments of ddaa type water molecules are larger than those of dda or daa type molecules. The ddaa type water molecule which accepts two hydrogen bonds and at the same time donates two hydrogen bonds is highly polarized.

As was indicated in [22], it is noteworthy that the dipole moments of ddaa type water molecules in water decamer is close to that in bulk water calculated in the MD simulation with the TTM2-R water model (2.65 D) [23]. Thus, the large dipole moment of a water molecule in bulk water is attributed to the large population of the ddaa type water molecules, which is highly polarized inherently. A dipole moment of a water molecule depends on the local type of hydrogen bonding. The similar results have been obtained from DFT calculations of water clusters, in which the importance is stressed of the polarization and delocalization between water molecules [39]. In bulk water, ddaa type water molecules are dominant in average. The average number of hydrogen bonds per molecule in liquid water at 300 K was reported [40] to be 3.6 from the MD simulation with the SPC/E model.

Classification of O–H bond in PW cluster

All the O–H bonds in PW cluster are classified into 10 topological types (Fig. 13), utilizing the combination of the local HB types of the contributing water molecules in the O–H bond. Six of them are hydrogen-bonded O–H bonds. Four of them are free O–H bonds [43]. The free O–H bonds are explicitly represented by small vertices, each of which links to a larger vertex in the digraph representation in Fig. 13.

Fig. 13:

O–H topological types in protonated water clusters. Local minimum structures of representative small protonated water clusters containing two Zundel and three Eigen isomers are shown. In their digraph representations, the O–H topological types are illustrated. Each of small vertex with undirected edge in digraphs represents free O–H. The O–H topological types contain six H-boned O–H and four free O–H. Free O–H of acceptor water A has two stretching modes: symmetric stretching A(s) and anti-symmetric stretching A(a).

Distribution of PW octamer

Rooted digraph is used to represent topologically distinct isomers of PW cluster [19], [20]. The numbers of enumerated rooted digraphs of PW clusters H₂O⁺(H₂O)_n−1 (n=2–8) are given in [19]. The number of rooted digraphs for PW octamer (n=8) is 896 063: this is the maximum number of its topologically distinct structures. It is not trivial to search full potential energy surface of PW octamer using a reliable level of calculation. Making use of different computational methods [41], [42], 134 topologically distinct local minima of PW octamers were obtained [43] at the theoretical level of MP2/aug-cc-pVDZ. The lowest energy structure is TP1, which is Eigen type. TP26 is the 26^th stable isomer, and the most stable isomer of the Zundel type. To gain insights into the thermodynamic stability of topologically distinct isomers, temperature-dependent populations of the optimized local minimum structures were estimated on the basis of the molecular partition functions [43] using the left method of Fig. 5. At a temperature lower than ~140 K, TP1 is dominant. In the range between ~150 and ~250 K, TP26 (Zundel isomer) is dominant being followed by TP33. At higher temperature (>270 K), a chain structure TP94 is dominant. The optimized structures and the corresponding digraphs of these four PW octamers are shown in Fig. 14. The O–H bonds are drawn with coloring same as in Fig. 13.

Fig. 14:

Optimized structures and digraph representations of representative protonated water octamers. The O–H topological types in each structure are indicated in the corresponding rooted digraph.

Vibrational frequency of O–H topological type

All O–H vibrational frequencies of 25 selected isomers of PW octamer are plotted in Fig. 15. The O–H vibrational frequencies spread from ~1200 cm⁻¹ to ~3800 cm⁻¹: they can be classified according to the topological type of the O–H bond [43]. The vibrational frequency for free O–H of acceptor water molecule A is further classified into symmetric stretching A(s) and anti-symmetric stretching A(a) according to the vibration modes (Fig. 13).

Fig. 15:

Harmonic frequencies of O–H stretching modes of protonated water octamers. The wavenumber was calculated at the MP2/aug-cc-pVDZ level of theory and scaled by 0.955.

If the O–H topological type (Fig. 13) is the same, the corresponding O–H vibrational frequency is similar even in different isomers. The O–H topological types can be used to distinguish spectral signatures of O–H stretching modes. Assign the O–H bond the topological type, and you can make a rough estimate about the vibrational frequency.

Calculated IR spectra and assignment of coexisting isomers

The calculated harmonic IR spectra of the representative structures, TP1, TP26, TP33 and TP94 in the range of 2000–3800 cm⁻¹ are given in [43]. The range of 3600–3800 cm⁻¹ corresponds to the free O–H stretching modes of the topological types including A(s) (~3630 cm⁻¹), D1A2 (~3670 cm⁻¹), D1A1 (~3700 cm⁻¹) and A(a) (~3750 cm⁻¹). The number of the bands of TP26 and TP33 is four in this range and it corresponds with the number of the bands of the observed spectra of PW octamer [44], [45]. Since there is a spectral signature of Zundel core (the shared proton oscillation) in the observed spectra [44], it is certain that TP26 exists. The coexistence of either TP33 or TP1, or both of them, with TP26 does not conflict with the observed spectra [44].