Classification of hydrogen bond flips in small water polyhedra applied to concerted proton tunneling

Abstract

Recently a new mechanism of proton tunneling in a prism-like water hexamer was revealed [Richardson et al., Science, 2016, 351, 1310]. The tunneling motion involves the concerted breaking of two hydrogen bonds and rotations of two nearest water molecules. Eventually, this structural transformation means flipping one of the hydrogen bonds without the creation of defects in the hydrogen bond network. On the surface of polyhedral water clusters, there are five essentially different types of hydrogen bonds, and only two of them can be changed in this manner. In this article, the topological classification of such transformations for five small water polyhedra: triangular, pentagonal, and hexagonal prisms as well as cube and polyhedron 4⁴5⁴, consisting of four square and four pentagonal faces, is presented. Our classification includes the enumeration of all possible one-bond-flips with consideration of the types of hydrogen bonds on the polyhedral surface. Attention is paid to the most stable proton configurations which can be studied in experiments. It was established that a number of one-bond-flip transitions between the low energy configurations are possible in clusters in the shape of triangular and pentagonal prisms.

---

1. Introduction

The properties of water and ice are still far away from a complete understanding. In particular, it is due to the fact that the ice is not a crystal in the usual sense of the word. The crystal lattice of ice determines the position of the oxygen atoms only. The positions of the hydrogen atoms are disordered. And, the position of hydrogen atoms (protons) changes constantly.^1,2 Proton transfer along hydrogen (H–) bonds is a quantum process that has been closely studied for several decades. A considerable attention is paid to quantum tunneling in small water clusters,^3–6 taking into account the presence of powerful experimental methods of terahertz spectroscopy and high accuracy of quantum chemical calculations. Recently, the ring-polymer instanton method has been applied with much success for calculating tunneling splitting in water clusters.^7,8

A large number of studies have been devoted to the investigation of rotational transitions of individual water molecules in clusters.^4–6 But recently it has been reported that the tunneling motion in the water hexamer prism involves the concerted breaking of two H-bonds with the rotation of two hydrogen-bonded molecules.⁹ The cluster in the form of a triangular prism is the smallest polyhedral water cluster (PWC). Each of the molecules takes part in three hydrogen bonds. Eventually, the concerted H-bond breaking means changing the direction of one H-bond (one bond flipping, see Fig. 1) without the formation of defects, that is, according to the Bernal–Fowler ice rules.¹⁰ In graph theoretical representation, according to these rules for PWCs one or two arrows (not 0 and not 3) point in and out in each node.

Fig. 1 Transition between the tunneling states of the lowest-energy structure of the water hexamer (a) and corresponding one bond flipping in the graph-theoretical representation (b).

In spite of the complex nature of the intermolecular interactions, there are effective approaches to determine the correlation between the energy of PWCs and topological features of the H-bond network.^11–21 One of these approaches is the Strong–Weak-Effective-Bond (SWEB) discrete model of the molecular interaction.¹⁷ This physically based model, that takes into account the peculiarities of the interaction between the first, second and third neighbors in PWCs, allows the understanding of the statistical regularities. All proton configurations of PWCs can be separated into distinct energy levels using only local topological features of the H-bond network. At first we consider the transitions between the lowest energy configurations which can be observed in experiments. Note that comprehensive statistics of all proton configurations for PWCs were obtained in ref. 22–24.

The present article is devoted to the topological classification of all one-bond-flip transitions in small water polyhedra including the triangular, cubic, pentagonal and hexagonal prisms as well as octahedron 4⁴5⁴ that consists of four square and four pentagonal faces.

2. Methodology

A. Statistics of proton disorder

PWCs can be assigned to two-dimensional systems by the degree of correlation. All molecules are located on one surface. Each of the molecules forms three H-bonds with the neighbor molecules as well as in the hexagonal ice monolayer. The specific residual entropy of these systems (in the dimensionless form) S₀ approximately equals ,²² in contrast to well-known Pauling's formula for three-dimensional systems: S₀ = ln(3/2).²⁵ An exact statistic of proton disorder in PWCs is a starting point for our analysis. It is important to emphasize that for the study of proton tunneling not only the set of symmetry distinct configurations is of interest but also a complete set of all proton configurations, corresponding to the notion of the residual entropy. Transitions also occur between different symmetry equivalent configurations. For many of the PWCs, the numbers of all configurations were obtained earlier.²³ An exact formula for prism-shaped water clusters is given in ref. 22:where X_n is a total number (integer) of all proton configurations for the clusters in the form of an n-gonal prism. For triangular, rectangular (cube), pentagonal, and hexagonal prisms, X_n equals 104, 450, 2008, and 9074. The numbers of symmetry distinct configurations for these clusters are 10, 14, 102, and 408, respectively.^11,16 For water octahedron 4⁴5⁴ these numbers are equal to 8976²³ and 1122.¹⁶ The listed clusters possess the following symmetry: D_3h (12), O_h (48), D_5h (20), D_6h (24) and D_2d (8). Here, the order of symmetry groups is given in brackets.

B. Energetics of proton configurations

PWCs and gas hydrate frameworks are formed only by eclipsed (mirror-symmetrical) water–water bonds, whereas ordinary hexagonal ice is formed by both eclipsed and staggered (center-symmetrical) bonds.^1,2 These bond types are shown in Fig. 2. As it was noted in the early works of Bjerrum, each of these H-bonded pairs can be divided into two subtypes depending on the mutual orientations: more and less energetically favorable (strong/weak).^26,27 For the eclipsed bonds, the strong H-bonds correspond to trans-configurations of water dimers (Fig. 2a) and the weak H-bonds correspond to cis-configurations (Fig. 2b). This was first noted by Radhakrishnan and Herndon and applied to ring-like water clusters.²⁸

Fig. 2 Eclipsed (a and b) and staggered (c) H-bonds. Classification of the hydrogen bonds on the surface of polyhedral water clusters (d). Arrows indicate the direction of the H-bonds: from the proton donor to the acceptor.

In PWCs, the same types of H-bonds can be differently oriented to the polyhedral surface. Therefore, there are five distinct types of H-bonds in this case.^12,15,17,23 Two of them (t0, t1) correspond to the trans-configuration of H-bonded pairs (Fig. 2d), and three types correspond to the cis-configuration (c0, c1, and c2). This graph theoretical classification of H-bonds by the edge node forms the basis of the Strong–Weak-Bond (SWB) discrete model.²³ Subsequently, it was established that only one of these bond types (t1) is preferable taking into account the interaction between the second and third neighbors. For large PWCs, the high predictive ability of a more accurate SWEB¹⁷ model was examined by different methods. Correlation between the energy and the number of t1-bonds is also high even at relatively considerable bending of H-bonds in small PWCs (Fig. 3). Geometrically optimized configurations in the TIP4P water model²⁹ were obtained using TINKER Molecular Modeling package.³⁰ The numbers of symmetry distinct configurations are shown in Fig. 3 for each value of t1. One configuration of the triangular prism-like cluster has changed its shape during optimization (configuration 1, see below). The lowest energy and near-lowest energy configurations are of the most interest for our analysis. In the SWEB model this means the maximum and the next-maximum number of t1-bonds.

Fig. 3 The binding energy per molecule (kJ mol⁻¹) for proton configurations of PWCs in the TIP4P water model as a function of the number of t1-bonds which determines the energy levels. The numbers at points correspond to the level populations. The configuration denoted by square (Prism_3) has changed its shape during optimization (unstable).

Note that all configurations of the lowest energy levels have a minimal total dipole moment (zero for cubic and hexagonal prisms).

C. Transfer matrices

First of all, let us draw attention to the fact that the one-bond-flip process is possible only for t0- and c0-bonds. Indeed, for the other three types of H-bonds, this leads to the violation of the Bernal–Fowler ice rules. This allows us to suggest the following algorithm for the calculation of the one-bond-flips:

1. Find all proton configurations of the PWC that satisfy the Bernal–Fowler ice rules (N). Note that in a topological representation, each proton configuration is defined by a row of the H-bond directions with respect to somewhat canonical orientation (numbers −1 and +1). For the water hexamer prism, it means 104 rows by 9 numbers. This is the complete list of configurations.

2. Separate all symmetry distinct configurations (n). For the triangular prism, there are ten rows. This is the short list of proton configurations. Each of these configurations has a number of physically equivalent configurations that differ by a rotation or a reflection only. The complete list is formed by combining the classes of equivalent configurations. Each of the symmetry distinct configurations is a representative of a whole class of symmetry equivalent configurations. The maximum number of equivalent configurations does not exceed the order of the symmetry group ν for the given polyhedron. Therefore, n ≥ N/ν. For large PWCs, n ≈ N/ν, since the majority of proton configurations are unsymmetrical (symmetry C1). This relation becomes exact if there are no symmetrical configurations (octahedron 4⁴5⁴, see above).

3. Turn sequentially all t0- and c0-bonds in configurations of the short list and identify the derived configurations in the complete list. Then find the equivalent configuration in the short list. The result can be represented in two square matrices by the size n × n. The elements of the first matrix are the numbers of one-bond-flips from i to j configurations. The elements of the second matrix are the series of numbers corresponding to the serial numbers of changeable H-bonds.

3. Results and discussion

A. Triangular prism

All symmetry distinct configurations of the water cluster in the form of a triangular prism are shown in Fig. 4. The arrows indicate the direction of the H-bonds: from the proton donor to the acceptor. The bold grey arrows correspond to t1-bonds, the number of which determines the energy of the cluster in the SWEB model. The double and dashed arrows correspond to t0- and c0-bonds, which allow the change in the H-bond direction without violating the Bernal–Fowler rules. The matrix of possible one-bond flips has the following form:Each non-zero element a_ij of this matrix corresponds to the number of one-bond flips that transform configuration i into configuration j. The numbers of t1-bonds of each configuration are shown at the left and above the matrix. The serial numbers of the flipping H-bonds are given in the following matrix:

Fig. 4 Symmetry distinct configurations for the triangular prism-like cluster. All possible one-bond-flips in the lowest energy configurations (t1 = 2) are indicated by ovals. The grey, double, and dashed arrows correspond to t1-, t0- and c0-bonds.

The H-bond numbering is shown in configuration 1 (Fig. 4).

In the TIP4P water model, configuration 3 is the most stable (Fig. 3). Note that this prism-like configuration is the most energetically favorable structure according to ab initio calculations.^31,32 The other two configurations (7 and 8) of the same level (t1 = 2) are less favourable (see Fig. 3). According to matrix (2), in each of these configurations there are H-bonds that allow the one-bond-flip transformation to itself (automorphic transformations). Each of the configurations (7 and 8) has two changeable H-bonds. In Fig. 4 these H-bonds are indicated by ovals.

The proton configurations before and after these transformations are physically equivalent. One can be obtained from another by a rotation or reflection. In the short list of configurations, they are represented by the same structure. At the same time, they are different. Transition between them is observed in the experiment. The automorphic tunneling transition in the lowest energy configuration 3 was considered by Richardson et al.⁹ Note that the energy barriers can be symmetrical for the automorphic transitions.

One can see in Fig. 4 that automorphic transitions in configurations 3 and 8 go through t0-bonds (trans-configurations), whereas in configuration 7 this transition goes through the c0-bond (cis-configuration). Although, the notions of trans- and cis-configurations are considered here from the topological point of view, all H-bonds are bent, especially in triangles. It is not difficult to verify that the t0-bond transitions connect the proton configurations that differ by the rotation around the 2-fold axis passing through the middle of the changing H-bond, whereas the c0-bond transitions connect configurations that differ by 3-fold rotatory reflection, that is, which are chiral isomers. Also note that the transfer matrix (2) is not symmetrical. This is related to the fact that the configurations themselves have different symmetries. For example (see matrix (2)), in high-symmetrical configuration 10 each of the transitions through the 3th, 5th, and 6th H-bonds leads to the same configuration 8. But the reverse transition from configuration 8 to configuration 10 is possible through the 3th H-bond only (Fig. 4).

In Fig. 5, all transitions between the cluster configurations are presented in the form of connectivity graphs. The vertical positions of the graph's nodes correspond to the energy levels within the framework of the SWEB discrete model (t1 = 0, 1, 2). Here, along with the automorphic transition in configuration 3 (see above), there are two pairs of analogous transitions in configurations 7 and 8 (t1 = 2). But there is no transition between these lowest energy configurations. Note one more peculiarity of this graph. The pair of configurations 2 and 6 is a duplicate of pairs 1 and 5 and vice versa. The topological connectivity of each configuration 1 and 2 with the others is identical, the same situation with configurations 5 and 6. This is discussed in Section 3E.

Fig. 5 Connectivity graph of the water cluster in the form of the triangular prism.

B. Cubic cluster

It is well known that the cubic water cluster has two lowest energy isomers of D_2d and S₄ symmetries (t1 = 4).³³ In the SWEB model, the next energy level (t1 = 3) is not occupied. The connectivity graph is shown in Fig. 6. Here, configuration 5 is the duplicate of configuration 3 (see Section 3E). Unfortunately, for the lowest energy configurations 12 and 13 there is neither automorphic transition nor transition between them. Only the transitions that significantly increase the energy are possible. Symmetry distinct configurations of the cubic cluster are shown in Fig. 7.¹¹ It was established that for the lowest energy configurations, only the transitions 12 ⇄ 2, 12 ⇄ 9, and 13 ⇄ 7 are possible. The changeable H-bonds are indicated by ovals (Fig. 7). Configurations in these pairs differ from each other by marked bonds only. Transition 12 ⇄ 9 goes through the c0-bond (cis-configuration), whereas the other two transitions go through the t0-bond (trans-configuration). Note that the lack of the lowest energy transitions makes the cubic cluster less interesting for our analysis.

Fig. 6 Connectivity graph for the cubic cluster.

Fig. 7 Symmetry distinct configurations of the cubic cluster. One-bond-flip transitions 12 ⇄ 2, 12 ⇄ 9, and 13 ⇄ 7 with the participation of the lowest energy configurations 12 and 13 are indicated by ovals.

C. Pentagonal prism

For the pentagonal prism-like cluster, there are six configurations with the maximum number of t1-bonds (t1 = 4, see Fig. 3). Note that these six configurations were also constructed by the insertion method, and four of them were found experimentally.³⁴ These configurations and also their connectivity with the second level configurations (t1 = 3) are shown in Fig. 8. The numbers with the asterisk correspond to pairs of configurations with the opposite direction of all H-bonds. Therefore, the total number of second level configurations equals 10 according to Fig. 3.

Fig. 8 Symmetry distinct configurations of the lowest energy level for the pentagonal prism-like cluster (a). The automorphic one-bond-flip transitions are indicated by ovals. Dashed ovals correspond to the transitions between different configurations. Connectivity graph (b). The numbers with asterisk correspond to pairs of configurations with the opposite direction of all H-bonds.

As it was noted by Richardson et al.⁹ this cluster can be of interest for the study of the concerted rotational transitions in water clusters. For the lowest energy configurations, there are both automorphic transitions and transitions between them. The transitions in configuration 1 go through the c0-bond (cis-configuration) and the residual transitions go through the t0-bond (trans-configuration) (Fig. 8). The most stable configuration in the model TIP4P is configuration 1,³¹ although configuration 2 is only slightly less stable. In Fig. 3 these configurations are indistinguishable. Five of these six configurations (Fig. 8a), with the exception of configuration 4, were found using the high-level method MP2/6-31G*.³² In this case, configuration 2 is the most stable. The same result was obtained by the method RHF/6-31G(d,p).³¹ It should be emphasized that the effects of tunneling were detected in three isomers of this cluster.³⁴ Also note that for the water decamer, the prismatic form is more energetically favorable.

D. Hexagonal prism and octahedron 4⁴5⁴

The lowest energy configurations of clusters in the form of the hexagonal prism and octahedron 4⁴5⁴ with the maximum number of t1-bonds are shown in Fig. 9. Additionally, for the last cluster, the lowest energy configuration of the second energy level is shown. Note that all such figures were drawn with the help of a special program. Strictly speaking, the number of optimal configurations for octahedron 4⁴5⁴ equals 8 (not 6), see Fig. 3. Two additional configurations correspond to changing the directions of all H-bonds in configurations 1 and 3. It was established that there is no on-bond-flip transition between configurations of the lowest energy levels for each of these clusters. In this connection, the lowest energy configuration of the second level for the octahedron-like cluster 4⁴5⁴ is of great interest. Here, there are two automorphic transitions through t0-bonds (Fig. 9). It should be noted that this configuration essentially exceeds all other configurations of the second energy levels for both cluster forms in the binding energy for the TIP4P model (lowest point in Fig. 3). According to the Cambridge Cluster Database,³¹ this configuration is the lowest in energy in the model TIP5P.³⁵

Fig. 9 Symmetry distinct configurations of the lowest energy levels for the clusters in the form of hexagonal prism (a) and octahedron 4⁴5⁴ (b). In the low row, an additional configuration of the second level is shown. The asterisk corresponds to the pair of configurations with the opposite direction of all H-bonds.

E. Antisymmetry

We consider the reversal of all H-bonds as an additional nongeometric operation of symmetry, more precisely antisymmetry.^36–38 This notion makes the enumeration and classification of the proton configurations easier, especially for the large systems where the majority of proton configurations are nonantisymmetrical (ordinary). The ordinary configurations form pairs which are distinguished by flips of all H-bonds and possibly by a rotation or reflection. Therefore, the list of essentially different proton configurations can be almost twice shortened. The results obtained in this study corroborate the utility of the new concept of the hydrogen-bonding antisymmetry. As it is easy to see (Fig. 5), the antipodal configurations with the opposite direction of all H-bonds can be symmetrically located in the connectivity graph. Configurations 2 and 6 (Fig. 5) are the antipodes of 1 and 5 (respectively), though the antipodal configurations are distinguished by the front H-bond only (bond #6, Fig. 4). For the cubic cluster (Fig. 6), the antipodes 3 and 5 are also symmetrically located in the connectivity graph. Under the extended symmetry (antisymmetry) the splitting disappears, since the antipodal configurations have become equivalent. In the small clusters, the antipodal configurations can be essentially different. So, one of the antipodes (configuration 1) has changed its shape during the geometrical optimization using the TIP4P water model (unstable, see Fig. 3). The majority of the lowest energy configurations of the considered clusters are antisymmetrical except configurations 1 and 3 for octahedron 4⁴5⁴. The antisymmetrical configurations are transformed into symmetrically equivalent ones under flipping of all H-bonds. We can say that the automorphic configurations considered are locally antisymmetrical.

4. Concluding remarks

In this article we present a combinatorial-topological classification of the concerted H-bond breaking by quantum tunneling for five smallest water polyhedra. We show that such structural transformations mean the change of the direction of one H-bond (one bond flipping) without the creation of defects. Our classification includes enumeration of all possible one-bond-flips in the lowest energy proton configurations considering the types of H-bonds on the polyhedral surface. It would seem natural that the proton transitions through t0- and c0-bonds are distinguished.

It was established that in the energy lowest configurations, the one-bond-flip transitions are possible only in the clusters in the shape of triangular and pentagonal prisms and also of octahedron 4⁴5⁴. For the pentagonal prism-like cluster, not only the automorphic transitions but also the transitions between near-in-energy configurations are possible. We mention in passing that the large PWC in the form of a pentagonal dodecahedron that consists of twenty molecules can also be of interest as a model system. The lowest energy level (t1 = 7) involves 106 symmetry distinct configurations, which are divided into three groups according to the total dipole moment.¹⁷ But within the lowest energy group (28 configurations with the smallest dipole moment) there are neither automorphic transitions nor transitions between them.

Finally note that the quantum description of the proton tunneling can help to parameterize the kinetics of the proton transfer. This gives the opportunity of more realistic modeling of ice-water systems. In such an event, the transfer matrices like (2) can be used to describe the proton dynamics in PWCs, taking into account the H-bond type (t0, c0), the energy difference as well as the numbers of symmetrically equivalent structures before and after the transitions.

Acknowledgements

The reported study was partially supported by RFBR, research project No. 15-03-04274.

References

1. V. F. Petrenko and R. W. Whitworth, Physics of Ice, Oxford University Press, 1999.
2. E. B. Moore and V. Molinero, Phys. Chem. Chem. Phys., 2011, 13, 20008–20016.
3. A. C. Legon and D. J. Millen, Chem. Soc. Rev., 1992, 21, 71–78.
4. N. Pugliano and R. J. Saykally, Science, 1992, 257, 1937–1940.
5. R. J. Saykally and G. A. Blake, Science, 1993, 259, 1570–1575.
6. C. Pérez, S. Lobsiger, N. A. Seifert, D. P. Zaleski, B. Temelso, G. C. Shields, Z. Kisiel and B. H. Pate, Chem. Phys. Lett., 2013, 571, 1–15.
7. J. O. Richardson, S. C. Althorpe and D. J. Wales, J. Chem. Phys., 2011, 135, 124109.
8. J. O. Richardson, D. J. Wales, S. C. Althorpe, R. P. McLaughlin, M. R. Viant, O. Shih and R. J. Saykally, J. Phys. Chem. A, 2013, 117, 6960–6966.
9. J. O. Richardson, C. Pérez, S. Lobsiger, A. A. Reid, B. Temelso, G. C. Shields, Z. Kisiel, D. J. Wales, B. H. Pate and S. C. Althorpe, Science, 2016, 351, 1310–1313.
10. J. D. Bernal and R. H. Fowler, J. Chem. Phys., 1933, 1, 515–548.
11. S. McDonald, L. Ojamae and S. J. Singer, J. Phys. Chem., 1998, 102, 2824–2832.
12. D. J. Anick, THEOCHEM, 2002, 587, 97–110.
13. D. J. Anick, J. Chem. Phys., 2003, 119, 12442–12456.
14. S. D. Belair and J. S. Francisco, Phys. Rev. A: At., Mol., Opt. Phys., 2003, 67, 063206.
15. V. Chihaia, S. Adams and W. F. Kuhs, Chem. Phys., 2004, 297, 271–287.
16. A. Lenz and L. Ojamäe, Phys. Chem. Chem. Phys., 2005, 7, 1905–1911.
17. M. V. Kirov, G. S. Fanourgakis and S. S. Xantheas, Chem. Phys. Lett., 2008, 461, 180–188.
18. P. Parkkinen, S. Riikonen and L. Halonen, J. Phys. Chem. A, 2013, 117, 9985–9998.
19. A. M. Tokmachev, A. L. Tchougréeff and R. Dronskowski, Theor. Chem. Acc., 2015, 134, 115.
20. V. Shilpi, S. P. Kaur and C. N. Ramachandran, RSC Adv., 2015, 5, 74270–74273.
21. S. Iwata, D. Akase, M. Aida and S. S. Xantheas, Phys. Chem. Chem. Phys., 2016, 18, 19746–19756.
22. M. V. Kirov, J. Struct. Chem., 1994, 35, 126–128.
23. M. V. Kirov, J. Struct. Chem., 1996, 37, 84–91.
24. M. V. Kirov, Physica A, 2009, 388, 1431–1445.
25. L. Pauling, J. Am. Chem. Soc., 1935, 57, 2680–2684.
26. N. Bjerrum, Kgl. Danske Videnskab. Selskab, Math-fys. Medd., 1951, 27, 1–56.
27. N. Bjerrum, Science, 1952, 115, 385–390.
28. T. P. Radhakrishnan and W. C. Herndon, J. Phys. Chem., 1991, 95, 10609–10617.
29. W. L. Jorgensen, J. Chandresekhar, J. D. Madura, R. W. Impey and M. L. Klein, J. Chem. Phys., 1983, 79, 926–935.
30. J. W. Ponder, TINKER, Software Tools for Molecular Design, Version 6.2, Washington University School of Medicine, Saint Louis, MO, 2012http://dasher.wustl.edu/tinker.
31. D. J. Wales, J. P. K. Doye, A. Dullweber, M. P. Hodges, F. Y. Naumkin, F. Calvo, J. Hernaández-Rojas and T. F. Middleton, The Cambridge Cluster Database, 2001, http://www-wales.ch.cam.ac.uk/CCD.html.
32. R. M. Shields, B. Temelso, K. A. Archer and G. C. Shields, J. Phys. Chem. A, 2010, 114, 11725–11737.
33. C. J. Tsai and K. D. Jordan, J. Chem. Phys., 1991, 95, 3850–3853.
34. C. Pérez, D. P. Zaleski, N. A. Seifert, B. Temelso, G. C. Shields, Z. Kisiel and B. H. Pate, Angew. Chem., Int. Ed., 2014, 53, 14368–14372.
35. M. W. Mahoney and W. L. Jorgensen, J. Chem. Phys., 2000, 112, 8910–8922.
36. M. V. Kirov, J. Struct. Chem., 2002, 43, 790–797.
37. M. V. Kirov, J. Phys. Chem. B, 2014, 118, 13341–13348.
38. M. V. Kirov, Acta Crystallogr., Sect. A: Found. Adv., 2016, 72, 395–405.

---