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

ABSTRACT 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 10 12 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][5][6][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 CONTACT  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., ln η vs. 1/T is linear, whereas fragile liquids display Super-Arrhenius behaviour, described by the Vogel-Fulcher-Tammann (VFT) relation [8]: Here B, T 0 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 S conf 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 S conf 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 TS conf can be linearly expanded around the temperature T 0 at which it vanishes as TS conf = (A/B) · (T − T 0 ) /T 0 , (where now A/BT 0 is the coefficient of the linear expansion) Equation (2) coincides with Equation (1). In this case T 0 is called the Kauzmann temperature, being the T where S conf = 0.
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 S conf 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 S conf [13][14][15][16], which avoids the approximation of S conf as the difference between the liquid and the solid entropy commonly adopted in experiments [17][18][19]. PEL-based studies have offered the possibility to check the validity of the AG relation in several model potentials [20][21][22][23][24][25][26][27]. In most cases, it has been shown that a representation of ln D or ln η vs. (TS conf ) −1 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 S exc and the two-point entropy with dynamics [31].

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/cm 3 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 r 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 lim t →∞

Entropy
In this study we evaluate five different entropies. The entropy of the liquid S liq , the vibrational entropy S vib , the configurational entropy S conf , the excess entropy of the liquid with respect to the ideal gas S exc and a two body approximation of the latter S (2) . S liq at each state point is calculated from the total energy E liq and the free energy of the liquid F liq via In Ref. [28] we evaluated F liq of TIP4P/2005 using a series of thermodynamic and Hamiltonian integration steps [37][38][39] as well as E liq . The vibrational part of the entropy (S vib ) is split into two components. A harmonic component S harm and an anharmonic component S anh .
For the evaluation of S harm at least 30 equally spaced configurations were extracted from the MD production runs. These configurations were minimised using a conjugategradient algorithm. The configurations obtained in this way are the inherent structures (IS). For the IS we calculated the density of states by diagonalising the 6N × 6N Hessian matrix, thereby obtaining the eigenfrequencies ω i . From these frequencies we calculate S harm as Here denotes Planck's constant in its reduced form. To evaluate the entropic contribution arising form anharmonicities [28]S anh we fit the difference between the potential energy and the IS energy along isochores as a polynomial in T (with coefficients c i ), corresponding to Then S vib and S conf are calculated via and The excess entropy with respect to the ideal gas is computed from S id is calculated from the canonical partition function Z id 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, I x , I y and I z are the moments of inertia along the three principal axes, k B is Boltzmann's constant and h is Planck's constant. The factor 1/2 in front of Z R accounts for the water molecule's C 2v symmetry [40]. From this partition function we calculate the free energy of the ideal gas (14) and the ideal gas entropy Finally we compute a two-body approximation of the translational component of the excess entropy S (2) from the O-O-pair correlation function g(r) [31]: (16) where ρ denotes the molecule number density. 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 S conf 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 S conf . On increasing ρ, the network is progressively distorted increasing the number of potential energy minima and thereby S conf . At even higher ρ the Lennard-Jones repulsion becomes dominant, decreasing the number of minima and S conf . This non monotonic behaviour of S conf 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 S liq and D are consistent with the previous results on the SPC/E model of water [20]. The observed maximum in S vib in Figure 1 however contrasts the result for SPC/E, where no maximum was reported. From the T dependence of both S conf and D it is possible to check the validity of the Adam-Gibbs relation (Equation (2)). Figure 2(a) shows ln(D) vs. 1/TS conf . 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 D ∞ are shown in Figure 2(c). Both A and D ∞ 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 k B T/M is explicitly accounted for by normalising D, the quality of the AG fit remains identical.

Results
For completeness, we compare the TIP4P/2005 data also with the excess entropy scaling [31,[41][42][43][44][45][46][47] proposition which relates D with the excess entropy S exc . The  (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 D ∞ . 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 D * is plotted against S exc . D * is commonly defined as [42] to scale out the the trivial thermal velocity contribution k B T/M. 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/cm 3 . A reasonable collapse is observed for 1.1 < ρ < 1.3 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 S (2) (see Figure3(b)) an even worse scaling is observed, in accordance with the results of Chopra et al. [31] for SPC/E.

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), ln D is linear in 1/TS conf 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.