Tiling in the geometric model for water

Tiling in the geometric model for water

M. Girardi^I; W. Figueiredo^I; N. Guisoni^II; V. B. Henriques^II

^IDepartamento de Física, Universidade Federal de Santa Catarina, 88040-900, Florianópolis, Santa Catarina, Brasil

^IIInstituto de Física, Universidade de São Paulo, C.P. 66318, CEP 05315-970, São Paulo, SP, Brasil

ABSTRACT

Hydrogen bonded liquids like water present a rich thermodynamic behaviour due to the strength and directionality of the bonds. In a recent paper a geometric model based on Bernal's model for liquids was proposed to study the effects of the hydrogen bonds on the phase transitions of water, under pressure and temperature variations. Water molecules were assumed to stay at the vertices of coordination r (r = 4; 5; 6) of perfectly tiled polygons, and to have four links which allow up to four hydrogen bonds per particle. Mean field calculations yielded a phase diagram with three phases of different densities and a critical point at the end of the coexistence line between the high and low density phases. The three phases were considered to be liquids of different densities. In the present work we have shown that applying some geometric constraints to particle arrangements (thus correcting the system entropy, which was overestimated in the previous work), and allowing a variable number of links per molecule, leads to substantial alteration of the phase diagram. Three phases of different densities are still present, but no critical point appears. Two of the phases are solid, and one phase is amorphous.

1 Introduction

Water and other hydrogen bonded liquids are known to exhibit many anomalous properties [1-10], which seem to arise from the strength and directionality of hydrogen bonds (HB), and ensuing low density. The HB network, which percolates in liquid water, is thought to be responsible, among other things, for the increase of the isothermal compressibility and constant pressure specific heat upon cooling, and for the isothermal compressibility minimum at 46ºC.

In the context of molecular dynamics simulations, several microscopic models have been studied in an attempt to describe the peculiar properties of the associated liquid [1, 4, 5, 6]. These models focus attention mainly on the charge distribution on the rigid molecule and the parameters of the electrostatic and Lennard-Jones potential are adjusted to fit experimental data. Definition of hydrogen bonds depends on still controversial energetic and angular criteria, and are also parameter dependent.

A different approach is that of simplified models. In this case one hopes to include microscopic properties able to reproduce qualitatively the main features of liquid water behaviour, independently of adjustable parameters. One example is the square water model, which incorporates the directional character of the HB[10, 8, 11]. Square water is exactly solvable at T = 0 and gives an excellent estimation of the residual entropy of the ice [12]. The HB number is a decreasing function of temperature, but the model does not present an order-disorder transition driven by temperature[10, 8, 11]. Under an external electric field, the system displays a structural phase transition[10] at T=0. Used as a solvent, square water presents hydrophobic hydration [8], albeit insufficient to explain water solubility properties.

Density fluctuations seem to be an essential ingredient. The geometric model for water, proposed by two of us [9], allows for changes in the local environment of the molecules. It is based on the geometric description of liquids given by Bernal [13], and developed on the plane by Collins[17] a few decades ago. In two dimensions, the water molecules are disposed at the vertices of squares and triangles of equal sides, perfectly tiling the plane. Such tiling yields coordination numbers for sites (vertex) given by r = 4, 5 or 6. Four open links are allowed for each molecule, so that it can form up to four hydrogen bonds with its r neighbors. A mean-field treatment of the model, adapted from a previous study[17] for the inclusion of HBs, yielded coexistence of three phases of different densities, and a critical point at the end of a high-low densities coexistence line. The three phases present no regularity, and therefore were interpreted as liquid phases.

To obtain the above results, a tiling constraint on entropy was left aside. It was assumed that the entropic term in the free energy has two contributions: one coming from the distribution of four bonds among r neighbors, and another one, coming from the number of geometric arrangements of the sites. Particles were treated as independent, therefore the geometric constraint of perfect tiling was not taken into account in the calculation of the degeneracy of the spatial arrangements.

In the present work, we have modified the treatment of the geometric model for water imposing a geometric constraint to the neighborhood of each site, in order to have a perfect tiling of squares and triangles at least at the nearest-neighboring sites. We have also allowed for a variable number of links (n_c = 2, 3, 4), in order to represent possible misalignment under temperature variations. The geometric constraint on entropy was first calculated by Kawamura [14], for a purely hard core semi-infinite system and a solid-liquid transition was suggested. Similar corrections to entropy were found by Do et al. [15] and Yi et al. [16], in the treatment of the geometric model with a Lennard-Jones potential, for which they also proposed a melting line. We will show that the geometric correction for the water model produces a substantial change of the phase diagram: two of the phases of different densities are crystalline and no critical point is present.

This article is organized as follows: in Sec. 2 we review the geometric model for water and introduce a variable number of bonds per molecule. In Sec. 3 we perform a mean field calculation of the Gibbs free energy of the model, in which the constraint of tiling is introduced in the calculation of entropy. Finally, in Sec. 4, we present our discussions and conclusions.

2 The model

In 1960, Bernal proposed a geometric model for liquids, in which the fluid particles were placed at the vertices of regular or quasi-regular polyhedra. The choice of these randomly distributed three-dimensional objects was motivated by the presence of local ordering, together with the absence of long range order in the liquid phase. Collins [17] considered a two-dimensional version of Bernal's model that consists of a perfect tiling of triangles and squares of equal sides. He found a discontinuous phase change for a special energy condition, from a mean-field calculation, which he suggested to be analogous to a change of association number from one liquid phase to another. The noninteracting version of the model was employed later in the study of quasi-crystals [18, 19, 20], observed in alloy systems such as V-Ni and V-Ni-Si, in an attempt to explain their mechanism of formation.

In Fig. 1, we show a typical configuration of the system in an amorphous phase. The ordered (crystalline) phases are the square and triangular lattices, where the plane is completely filled by triangles or squares. The possible local structures around each site are shown in fig. 2. The sites can have three different coordination numbers¹ (r = 4, 5 and 6), whose specific volumes v_r are

where b is the intermolecular distance (polygon side), which is kept fixed (b = 1 in this work). Each lattice site is occupied by a molecule, and we define N_r as being the number of molecules with r neighbors. Volume and particle number conservation give

where V is the total volume and N is the total number of particles.

Hydrogen bonds will be defined as follows. Each molecule is connected to its r neighbors by r lines. It may have n_c = 2, 3 or 4 open links with its r neighbors. These links may or may not be available due to, for example, orientational restraints and, if available, are distributed randomly over the r lines. An HB will be present if two neighboring molecules point one of their links towards each other. In this way, the maximum number of hydrogen bonds for a given molecule is n_c. In Fig. 2, (S1) and (S2) represent pairs of molecules at sites of the type A and B, with coordination numbers r = 4 and s = 5, respectively. In the case of (S1), the A molecule has n_c = 4 links (represented by the full lines), and the B molecule, has n_c = 2. Only in (S2), where two links are in the line joining A and B, do we have a hydrogen bond between the two molecules.

Only hydrogen-bond energies shall be considered, since van der Waals interactions are one order of magnitude smaller than the hydrogen bond ones.

3 Mean field calculations

In this section, we consider a mean field approximation in order to find the Gibbs free energy of the present geometric model. The first step is to obtain the total energy of the system, which results solely from HBs. Following the definition of the previous section (see Fig. 2), we may write, for the probability of an HB between two molecules (assuming they are independent) whose sites have coordination numbers r and s

where n (n) is the number of open links of a molecule with coordination number r (s), and is the fraction of links per particle, which is independent of site-coordination. A natural choice for the potential energy between two sites is then f_rs = -eP_HB (here, e= 1).

The energy per particle is written as

where C_rs is the number of (rs) nearest pairs, given by C_rs = [rN_rp_s(r)+sN_sp_r(s)(1-d_r,s)]. The function p_s(r) is the probability of a pair of neighboring sites with coordination numbers r and s, which is given by

where

Here, N_5B and N_5C are the number of vertices of the type B and C, respectively (fig. 2). N_tot(r) is the number of neighbors to r particles, and the first neighbor geometric constraint, which precludes 4 and 6-particles from being neighbors is taken into account.

Now we must write the entropy of the system. Assuming N independent sites distributed among {N_r} sites, the number of spatial arrangements is simply

The last expression overestimates the number of possible states of the system, since tiling constraints were not taken into account. To correct the number of states given by eq. 12, we consider the factor[14, 16]

where n_r = N_r/N. The correcting factor F guarantees a perfect arrangement of polygons at least for the nearest-neighboring sites. For a further explanation concerning this factor, see references 14 and 16. In this way, the number of states due to this geometric constraint becomes

The newly calculated geometrically constrained number of spatial arrangements is significantly improved with respect to the non-constrained case. This can be seen from comparison with the known entropy of the athermal version of the square-triangle tiling model, given the area fraction occupied by each component. The exact partition function for the latter was obtained using the Bethe ansatz and, for the case of equal area fractions, the entropy per vertex [19, 21, 22] is

Maximization of W_o and W_g (eqs. 12 and 14) for this special case of area fractions yields, respectively, s_o ~ 1.38 and s_g ~ 0.55, with s = lnW, indicating a significant improvement of the estimated entropy by taking into account the constraint factor F.

We must yet compute the entropic contribution due to the distribution of the molecular links. Let us define as the number of sites with n_c links. Note that satisfies the relations M₂+M₃+M₄ = N and 2M₂+3M₃+4M₄ = N_c, where N_c = N . For a given set of values of , the number of arrangements of sites having n_c links is

The number of states for a given site with r neighbors and n_c links is the number of combinations of r over n_c and, therefore, the total number of states accessible to the system is

where is the number of sites having r neighbors and n_c links, which in this mean field approximation we take as = N_r/N. The entropy per particle is given by s = ln W and the partition function of the model is given by:

where P is the pressure, b = (k_BT)^-1, T is the temperature, and the prime in the summation indicates the constraints å_iN_i = N and å_iM_i = N. Writing z({ N_i,M_j} ) = e-Ts+Pv, we can identify the Gibbs free energy per particle as g(T,P) = z(T,P,{ N_i,M_j}_min), where {N_i,M_j}_min is the set of values of N_i and M_j which minimizes z. Solutions were sought for in the ranges 0 < N_i/N < 1, 0 < M_i/N < 1 and 2 < < 4 through the simulated annealing [23] algorithm.

4 Results and conclusions

In Fig. 3 we exhibit the phase diagram of the model in the plane P versus T for both fixed (four) and variable number of links. As can be seen in the figure, the diagram presents three phases of different densities, a triple point (tp), and no critical point. Note that allowing for a variable number of links (continuous lines in Fig. 3) stabilizes the intermediate density phase at higher temperatures, in relation to the case in that the number of links is fixed, which may be attributed to the additional source of disorder that comes from the variable number of links.

In Fig. 4 we have plotted the fractions n_r of r-coordinated sites, as a function of temperature at P = 0. It may be observed that the fractions n_r change abruptly at the transition temperatures. At low temperatures, the energetic term dominates, and we have a low density (LD) phase with n₄ = 1. As the temperature increases, the entropic term becomes important and we have two transitions to higher density phases (ID and HD phases). Note that the HD phase has n₆ = 1. In Fig. 5, we show the fractions n_r as a function of pressure at T = 0. At low pressures, the LD phase has n₄ = 1. Increasing the pressure, for P > 4.6 the LD phase also presents sites with coordination numbers r = 5 and r = 6, indicating the appearance of defects in this crystalline phase. Around P ~ 5 a transition to the HD phase is accomplished. We have a competition between the Pv term, which favours the HD phase, and the energy e contribution, which favours the LD phase. The behaviour of the n_r densities at temperatures and pressure away from the axes can be read from the table 1.

From the two figures and data in the table 1 it can be seen that the HD phase is the ordered triangular lattice, whereas the LD phase is the square lattice, except near the phase boundary, where defects appear. These two phases may be associated with crystalline solid phases. The ID phase has no long range order and can be described as amorphous.

Let us compare our results with those of the previous work[9], in which entropy of spatial arrangement of molecules was not geometrically constrained. At T = 0 the effects of entropy (of the spatial arrangements of molecules and of links) are absent and both phase diagrams coincide, presenting a LD-HD transition at P » 5 . However, away from this axis, the absence of a liquid-liquid phase transition line and a corresponding critical point, features present in our earlier model [9], are remarkable. The differences in the thermodynamic behaviour of the two systems arise mainly from the geometric constraint imposed on entropy, once variations on the number of links do not change the phase diagram qualitatively (see Fig. 3). Three points deserve our attention: i) The strongly reduced role of spatial entropy (see eqs. 12 and 14), associated with a restriction in the number of possible geometric arrangements, prevents the system from going continuously from the LD phase to HD phase; ii) The negative slope of the coexistence lines indicates, in accordance with Le Chatelier's principle, that entropy increases while volume decreases, on transitions to higher temperature. The energetic term favours the LD phase (f_rs is a minimum for r = s = 4), while the total entropy (geometric entropy plus entropy from the links) favours the phases of higher densities, leading to a competition between energy and entropy, at constant pressure. In our earlier work [9] we have shown that the unconstrained geometrical entropy favours the low density phase, while the link entropy favours the phases of higher density. These facts arise as interpretations of the phase diagram, which presents a coexistence line of positive inclination for the first case, of pure spatial disorder, and negative slope in case link disorder is included. In this work, spatial disorder is reduced, and link entropy is increased (through the introduction of variable link numbers). As a result, coexistence lines are all of negative slope. iii)The reduced role of spatial disorder is also apparent in the fact that the LD and HD phases are crystalline, for the variable link case. In these phases entropy increase with temperature is due solely to the link disorder, while the system remains quasi-crystalline. In the previous model the absence of the geometrical constraint allows smooth variation on each r-coordinated site fraction, implying in the presence of more heterogeneous phases, with no predominant n_r.

Figures 6 and 7 exhibit the fractions of hydrogen bonds and links per molecule as a function of temperature and pressure, respectively, at P = 0 and T = 0. As expected, in both figures, the fraction of hydrogen bonds, which is proportional to energy, is a decreasing function of temperature, and presents discontinuities at the transition temperatures. In the region of constant pressure (Fig. 6) the fraction of links also decreases with temperature, except at the ID-HD phase transition, at which it shows an evident increase. This behavior can be seen from the fact that total entropy must increase with temperature: while the HD phase has a low spatial entropy (n₆ = 1 as can be seen in Fig. 4) the number of links contributes to raise the entropy (eq. 16). For T = 0 (Fig. 7) the fraction of links remains constant, with = 4, and the number of hydrogen bonds is maximum (4 in the LD phase and 8/3 in the HD phase), corresponding to the minimum energy configurations.

In conclusion, we have improved a previous geometric model for water, including a variable number of links, and a first order correction to the geometric entropy, with the purpose of achieving a good tiling of the polygons. The phase diagram presents three phases of different densities, two crystalline (LD and HD) and one amorphous phase (ID). No critical point or liquid-liquid phase transition is present. However, the diagram has some resemblance to that near the ice-I_h, where, for some values of temperature, the system goes from the ice-I_h to the liquid water and from this to the ice-V as the pressure increases[24]. These qualitative features are seen in our Fig. 3, at intermediate temperatures. The model cannot account for the liquid-vapor transition since the distance between the molecules is kept fixed. However, the increase in the density as the temperature raises is akin to the behaviour of water at low temperatures.

Acknowledgements

MG and WF acknowledge financial support from the Brazilian agencies CNPq and FUNCITEC. NG and VBH acknowledge financial support from the Brazilian agency Fapesp (NG 03/00831-8 and VBH 01/11721-3).

References

[1] F. Franks, Water, A Comprehensive Treatise (Plenum, New York, 1972).

[2] P. H. Poole, F. Sciortino, T. Grande, H. E. Stanley, and C. A. Angell, Phys. Rev. Lett. 73, 1632 (1994).

[3] P. H. Poole, U. Essmann, F. Sciortino, and H. E. Stanley, Phys. Rev E 48, 4605 (1993).

[4] W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, J. Chem. Phys. 79, 926 (1983); M. W. Mahoney, J. Chem. Phys. 112, 8910 (2000).

[5] H. E. Stanley, M. C. Barbosa, S. Mossa, P. A. Netz, F. Sciortino, F. W. Starr, and M. Yamada, Physica A 315, 281 (2002).

[6] P. A. Netz, F. W. Starr, M. C. Barbosa, and H. E. Stanley, Physica A 314, 470 (2002).

[7] D. Chandler, Nature 417, 491 (2002).

[8] N. Guisoni and V. B. Henriques, Braz. J. Phys. 30, 736 (2000).

[9] N. Guisoni and V. B. Henriques, J. Chem. Phys. 115, 5238 (2001).

[10] M. Girardi and W. Figueiredo, J. Chem. Phys. 117, 8926 (2002).

[11] W. Nadler and T. Krausche, Phys. Rev. A 44, R7888 (1991); T. Krausche and W. Nadler, Z. Phys. B: Condens. Matter 86, 433 (1992).

[12] E. H. Lieb, Phys. Rev. 162, 162 (1967); E. H. Lieb, Phys. Rev. Lett. 18, 692 (1967).

[13] J. D. Bernal, Nature 185, 68 (1960); J. D. Bernal, Nature 183, 141 (1959).

[14] H. Kawamura, Prog. Theor. Phys. 70, 352 (1983).

[15] D. Yi-Jing, C. Li-Rong, and Y. Tzu-Tung, J. Phys. C: Solid State Phys. 15, 3059 (1982).

[16] Y. M. Yi and Z. C. Guo, J. Phys. Condens. Matter 1, 1731 (1989); Z. Guo and Y. Yi, Commun. in Theor. Phys. 8, 17 (1987).

[17] R. Collins, Proc. Phys. Soc. 83, 553 (1964).

[18] B. Rubinstein and S. I. Ben-Abraham, Mat. Sci. & Eng. A 294, 418 (2000).

[19] D. Joseph and V. Elser, Phys. Rev. Lett. 79, 1066 (1997).

[20] M. Oxborrow and C. L. Henley, Phys. Rev. B 48, 6966 (1993).

[21] M. Widom, Phys. Rev. Lett. 70, 2094 (1993).

[22] H. Kawamura, Physica A 177, 73 (1991).

[23] W. L. Goffe, G. D. Ferrier, and J. Rogers, J. Econom. 60, 65 (1994).

[24] A. M. Saitta and F. Datchi, Phys. Rev. E 67, 020201 (2003).

Received on 26 August, 2003

1

There are two possible geometries for the molecule with coordination number

r = 5 but in the previous work [9], it was considered just one.

Publication Dates

• Publication in this collection
  27 Apr 2004
• Date of issue
  Mar 2004

History

• Received
  26 Aug 2003