The Adam-Gibbs relation and the TIP4P/2005 model of water.

We report a numerical test of the Adam-Gibbs relation for the TIP4P/2005 model of water. The configurational entropy is here evaluated as the logarithm of the number of different basins in the potential energy landscape sampled in equilibrium conditions. Despite the non-monotonic behaviour which characterise the density dependence of the diffusion coefficient, the Adam-Gibbs relation is satisfied within the numerical precision in a wide range of densities and temperatures. We also show that expressions based on the excess entropy (the logarithm of the number of sampled microstates in phase space) fail in the region of densities where a tetrahedral hydrogen bond network develops.

Introduction

One of the most intriguing phenomena in condensed matter physics is the glass transition. The temperature (T) at which it occurs is conventionally defined at the T at which viscosity (η) reaches values of the order of  Pa s. Such high viscosities render flow impossible on reasonable observation timescales and hence the liquid turns into a disordered solid. The tendency of liquids to form disordered solids, i.e., glasses, on cooling varies. Protein solutions typically form disordered arrested states readily [1], making crystallisation the hard problem [2]. Others such as silica are both fairly easily obtained as a crystalline solid or as a glass [3]. Finally, there are substances that tend towards crystallisation, challenging the creativity of experimentalists if a glass is desired. An example of such a bad glass former is water. Several methods had to be devised to obtain it in its glassy form [4-7]. Besides the tendency of a material to form a glass, an interesting feature is the change of η or of the diffusion coefficient D upon cooling before reaching the glassy state. The way η or D change with T conveys information not only on the possibility to shape the supercooled liquid on its way to the glass, but also on the physical process controlling the slowing down of the dynamics. Typically two types of liquids are discerned. Strong liquids show an Arrhenius-type behaviour, i.e., vs. 1/T is linear, whereas fragile liquids display Super-Arrhenius behaviour, described by the Vogel–Fulcher–Tammann (VFT) relation [8]: Here B, and are (path dependent) constants.

An intriguing connection between thermodynamics and the VFT relation for dynamical quantities is offered by the Adam–Gibbs (AG) relation [9]: Here indicates the configurational entropy, proportional to the logarithm of the number of distinct relevant liquid configurations. The idea behind Equation (2) is based on the concept of ‘cooperatively rearranging regions’ (CRR), spatial regions in which relaxation processes take place cooperatively. The activation energy for such relaxations is assumed to be extensive in the number of atoms or molecules that make up the CRR. The increase in the size of the CRR becomes thus associated to a decrease of the configurational entropy and responsible for the slowing down of the dynamics on supercooling. Numerical studies have attempted to identify the CRR from atomistic configurations. For the case of water see for example Ref. [10]. If vanishes at finite T, the AG model is consistent with approaches predicting an underlying (thermodynamic) phase transition as origin of the divergence of relaxation time. If can be linearly expanded around the temperature at which it vanishes as , (where now is the coefficient of the linear expansion) Equation (2) coincides with Equation  (1). In this case is called the Kauzmann temperature, being the T where . Equations (1) and (2) can also be written in terms of the diffusion constant D or other structural relaxation times τ. Due to the possible decoupling between D and η the parameters in the equations are in some cases sensitive to the chosen observable. For a critical exam of the AG model see for example [11].

The potential energy landscape (PEL) approach identifies in the liquid state as the logarithm of the number of distinct basins sampled on the potential energy surface [12], each of them conventionally associated to the minimum potential energy of the basin, the so-called inherent structure (IS). The PEL approach offers also a consistent way to numerically evaluate [13-16], which avoids the approximation of as the difference between the liquid and the solid entropy commonly adopted in experiments [17-19]. PEL-based studies have offered the possibility to check the validity of the AG relation in several model potentials [20-27]. In most cases, it has been shown that a representation of or vs. is consistent with the numerical results.

Very recently we have reported a numerical study of the statistical properties of the PEL for the TIP4P/2005 model of water [28]. We have shown that a Gaussian PEL properly describes the simulation data, reproducing the thermodynamic anomalies characteristic of water and predicting the existence of a liquid-liquid critical point. In this contribution we expand the landscape analysis to dynamics, testing the validity of the Adam–Gibbs relation for TIP4P/2005, the most accurate classical water model to date [29,30]. For completeness, we also compare the T and density (ρ) dependence of the diffusion coefficient with other propositions relating the excess and the two-point entropy with dynamics [31].

Simulation details

NVT simulations

We perform NVT simulations of 1000 TIP4P/2005 molecules in a cubic box utilising GROMACS 5.1.2 [32] with a leap-frog integrator using a timestep of 1 fs. The temperature is controlled using a Nosé-Hoover thermostat [33,34] with a time constant of 0.2 ps. For the coulombic interactions we use a particle mesh Ewald treatment [35] with a Fourier spacing of 0.1 nm. For both the Lennard-Jones and the real space Coulomb interactions, a cut-off of 0.85 nm is used. Lennard-Jones interactions beyond 0.85 nm have been included assuming a uniform fluid density. Finally, we maintain the bond constraints using the LINCS (Linear Constraint Solver) algorithm [36] of 6th order with one iteration to correct for rotational lengthening. We investigate 14 different densities from 0.9 to 1.42 g/cm3 and seven different Ts between 200 and 270 K. Very long equilibration runs (up to 100 ns) followed by equally long production runs have been performed.

Diffusivity

To evaluate the diffusion coefficient D we use the standard approach employed in molecular simulations [37]. We calculated the mean-square displacement for all MD runs. Here the average is over all particles in the system and over different initial times t, while indicates the oxygen position of the generic molecule. Since the oxygen is significantly heavier than the hydrogen atoms, the oxygen position provides a good characterisation of the centre of mass. From the long-time limit of the MSD we evaluate D via the Einstein relation

Entropy

In this study we evaluate five different entropies. The entropy of the liquid , the vibrational entropy , the configurational entropy , the excess entropy of the liquid with respect to the ideal gas and a two body approximation of the latter .

at each state point is calculated from the total energy and the free energy of the liquid via In Ref. [28] we evaluated of TIP4P/2005 using a series of thermodynamic and Hamiltonian integration steps [37-39] as well as . The vibrational part of the entropy () is split into two components. A harmonic component and an anharmonic component . For the evaluation of at least 30 equally spaced configurations were extracted from the MD production runs. These configurations were minimised using a conjugate-gradient algorithm. The configurations obtained in this way are the inherent structures (IS). For the IS we calculated the density of states by diagonalising the Hessian matrix, thereby obtaining the eigenfrequencies . From these frequencies we calculate as Here ℏ denotes Planck's constant in its reduced form. To evaluate the entropic contribution arising form anharmonicities [28] we fit the difference between the potential energy and the IS energy along isochores as a polynomial in T (with coefficients ), corresponding to Then and are calculated via and The excess entropy with respect to the ideal gas is computed from is calculated from the canonical partition function of a system of non-interacting water-shaped molecules: This partition function can be split in a translational part and a rotational part where m is the mass of the water molecule, , and are the moments of inertia along the three principal axes, is Boltzmann's constant and h is Planck's constant. The factor in front of accounts for the water molecule's symmetry [40]. From this partition function we calculate the free energy of the ideal gas and the ideal gas entropy

Finally we compute a two-body approximation of the translational component of the excess entropy from the O-O-pair correlation function [31]: where ρ denotes the molecule number density.

Results

Figure 1 shows the self-diffusion constant as well as the calculated entropies for all studied T and ρ. All quantities show maxima in their respective ρ dependence, consistent with the unconventional dynamic and thermodynamic behaviour of water. As already pointed out [20] the maximum in can be explained as a balance of two mechanisms. At low ρ the structure of the liquid resembles the one of a tetrahedral network which is expected to be characterised by a small number of potential energy minima and therefore a small . On increasing ρ, the network is progressively distorted increasing the number of potential energy minima and thereby . At even higher ρ the Lennard-Jones repulsion becomes dominant, decreasing the number of minima and . This non monotonic behaviour of is already evident at high T. Indeed, the estimated total number of PEL basins does already show such non monotonic behaviour [28]. Also the maxima in and D are consistent with the previous results on the SPC/E model of water [20]. The observed maximum in in Figure 1 however contrasts the result for SPC/E, where no maximum was reported.

Figure 1.

Density dependence of the diffusion constant (a), the configurational entropy (b), the vibrational entropy (c) and the liquid entropy (d). The data are shown for all studied densities and temperatures.

From the T dependence of both and D it is possible to check the validity of the Adam–Gibbs relation (Equation  (2)). Figure 2(a) shows vs. . The same data are also reported in Figure 2(b) arbitrarily shifted in the y direction to highlight the linear dependence, consistent with the theoretical prediction of the Adam–Gibbs relation. The density dependence of the parameters A and are shown in Figure 2(c). Both A and appear to approach a density independent value at large ρ. Significant deviations are only visible at low ρ, in the region where a well connected network of hydrogen bonded molecules develops. We note on passing that if the thermal velocity is explicitly accounted for by normalising D, the quality of the AG fit remains identical.

Figure 2.

Semi-log plot of the diffusion constant D vs. 1/ for all studied state points. Part (a) shows the data as obtained and part (b) shows vertically shifted data. Here the solid lines represent best fits with the Adam–Gibbs relation (Equation (2)). Part (c) shows the parameters obtained by fitting the data shown in (a) with Equation (2). The main axis shows the parameter A and the alternative axis the pre-exponential constant .

For completeness, we compare the TIP4P/2005 data also with the excess entropy scaling [31,41-47] proposition which relates D with the excess entropy . The rationale behind this hypothesis is that dynamics is controlled by the total number of accessible microstates in phase space. The validity of such a scaling can be tested if an adimensional diffusion constant is plotted against . is commonly defined as [42] to scale out the the trivial thermal velocity contribution . Here ρ is the molecule number density and M the molecular mass.

We summarise the corresponding results in Figure 3(a). The data clearly show that there is no data collapse for densities lower than 1.1 g/cm3. A reasonable collapse is observed for only, i.e., in the region where the hydrogen bond network is significantly distorted. Finally we note that if one approximates the excess entropy with the translational contribution (see Figure3(b)) an even worse scaling is observed, in accordance with the results of Chopra et al. [31] for SPC/E.

Figure 3.

Plots of the reduced diffusion constant versus (a) and (b).

Conclusions

We have recently reported [28] an exhaustive investigation of the potential energy landscape of TIP4P/2005, one of the most accurate classic water models. This model is able to reproduce the complex pattern of thermodynamic anomalies which characterise liquid and supercooled water [48]. In that study an accurate evaluation of the number of PEL basins sampled at each T and ρ has been carried out, providing access to the configurational entropy. Here we combine the evaluated diffusion coefficient with the associated configurational and excess entropy to test the validity of the Adam–Gibbs relation and of the excess entropy scaling. Despite the complexity of the system (which is characterised by non-monotonic dependence of the diffusion coefficient on the density), is linear in for all densities. The present data, and the previously available data for the SPC/E model [20], provide strong support to the hypothesis that the number of PEL basins is a key quantity in controlling the dynamics in supercooled water.