E water model

Accepted Manuscript Rotational and translational dynamics of the SPC/E water model

Nadège Meyer, Vincent Piquet, Jean-François Wax, Hong Xu, Claude Millot PII: DOI: Reference:

S0167-7322(18)33112-X doi:10.1016/j.molliq.2018.08.024 MOLLIQ 9469

To appear in:

Journal of Molecular Liquids

Received date: Revised date: Accepted date:

15 June 2018 30 July 2018 5 August 2018

Please cite this article as: Nadège Meyer, Vincent Piquet, Jean-François Wax, Hong Xu, Claude Millot , Rotational and translational dynamics of the SPC/E water model. Molliq (2018), doi:10.1016/j.molliq.2018.08.024

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Rotational and Translational Dynamics of the SPC/E Water Model Nad`ege Meyera , Vincent Piquetb , Jean-Fran¸cois Waxa , Hong Xua , Claude Millotb,∗ a Universit´ e

de Lorraine, LCP-A2MC, F-57000 Metz, France de Lorraine, CNRS, LPCT, F-54000 Nancy, France

RI PT

b Universit´ e

SC

Abstract

Translational and rotational diffusion coefficients and the translation-rotation couplings of the SPC/E water model are investigated in liquid phase by molec-

NU

ular dynamics simulations in the micro-canonical ensemble using a sample of 512 molecules in a cubic box with periodic boundary conditions. The diffusion

MA

coefficients and cross-translation-rotation coupling coefficients are computed in the laboratory and/or molecular systems of axes using a Green-Kubo approach. The studied thermodynamic states explore pressures and temperatures in the following ranges: 210-373 K (1 bar); 300-700 K (400, 1000, and 3000 bar), 277 K

D

(1-6000 bar). Translation diffusion is also analyzed using the Einstein approach.

PT E

The rotational diffusion is also characterized by computing the first and second order correlation times of the OH bond vectors and of the axes of inertia of the water molecule.

Keywords: Rotational diffusion. Translational diffusion. Translation-rotation

CE

coupling. Water. Molecular dynamics simulation. SPC/E model.

AC

1. Introduction

The recently reported question of the existence of plastic phase in liquid wa-

ter under very high pressure has rekindled the interest for diffusion in molecular ∗ Corresponding

author Email address: [email protected] (Claude Millot)

Preprint submitted to Journal of LATEX Templates

August 6, 2018

ACCEPTED MANUSCRIPT

liquids[1, 2]. Indeed, in such phases, the translational dynamics stops while the rotational dynamics still survives. On the experimental point of view, the translational diffusion coefficient can be measured quite easily and accurately in nuclear magnetic resonance (NMR) experiments for example. On the other hand, the situation is not so clear for

RI PT

the rotational diffusion as the corresponding diffusion coefficients can only be obtained indirectly, for example, from nuclear magnetic relaxation, but relying on the assumption that a given diffusion model is valid.

On the simulation point of view, both diffusion properties can be computed

SC

very easily and accurately from molecular dynamics (MD) simulations. The translational diffusion coefficient can either be obtained from the velocity auto-

NU

correlation function or from the mean squared displacement while the rotational diffusion coefficient can be computed from the angular velocity autocorrelation function or from mean square angular displacement without resorting to any

MA

approximate model. Moreover, as the motion of the molecule can be described either in the laboratory frame or in the molecule frame, possible anisotropy of diffusion properties can also be studied.

D

Though the translational diffusion is commonly analyzed from MD simulation, relatively few studies are devoted to the explicit calculation of the rota-

PT E

tional diffusion coefficients and particularly to the anisotropy of the rotational diffusion tensor. Most of the time, the anisotropy of the rotational diffusion is interpreted through the analysis of correlation times of vectors linked to the

CE

molecule, such as the electric dipole vector, in-plane and out-of-plane vectors or the OH-bond vector of a water molecule, for example. Beyond the standard determination of translational diffusion coefficients, MD simulations, generally

AC

in connection with NMR relaxation measurements, have been used to analyze directly the rotational diffusion coefficient, and particularly its anisotropy in the molecular frame, for several molecular liquids and solutions. For example, simulation studies have been devoted to water[3, 4, 5], benzene[6, 7], anthracene[8], perylene[9], ammonium ion[10], quinoline and isoquinoline[11, 12, 13] in neat liquid phases and/or solutions. The coupling between translation and rotation 2

ACCEPTED MANUSCRIPT

has been the subject of numerous works and can be investigated by computing correlation functions involving translation and rotational dynamical properties [14, 15, 16, 17, 18]. The analysis of translational diffusion in water by MD simulation has been reported many times in the literature with numerous intermolecular potentials.

RI PT

Surprisingly, it is no so common to find a simultaneous extensive study of the translational and rotational diffusion. Simple intermolecular potentials such as TIP3P, TIP4P, SPC, SPC/E, TIP4P-2005 are popular candidates which have been the subject of many studies [19, 20, 21, 22, 23, 24, 25], but at selected

SC

thermodynamic conditions. The TIP4P-2005 model [26], designed in order to improve the modeling of low temperature water with respect to the predictions

NU

of the popular SPC/E model [27], has been the subject of an extensive study of the translation and rotational diffusion properties [4], including the analysis of the coupling between translation and rotational degrees of freedom. However,

MA

this study was limited to 1 bar over a temperature range between 210 and 310 K. Moreover, this potential is not so efficient at higher temperatures and pressures. The purpose of this study is to provide a consistent systematic investigation

D

with the SPC/E model of the rotational diffusion properties together with the translational properties for completeness and emphasizing also the translation-

PT E

rotation couplings. Different values of the components of the rotational diffusion tensor at ambient conditions have been published with this model by different groups [3, 5, 28]. Early calculations of the coupling between translational veloc-

CE

ity and angular velocity of the SPC/E water model [29] and more recent ones [5] have been obtained at a few state points. In this study, these data will be completed by systematic calculations as a function of temperature and pressure.

AC

When available, we also will compare our results with the TIP4P-2005 model. We shall also discuss several points as the anisotropy of these properties, the

accuracy of models used to estimate the rotational diffusion coefficient and the pressure effect. In Section 2, the details of the MD simulations are presented as well as the calculation of the quantities of interest (diffusion coefficients, translation3

ACCEPTED MANUSCRIPT

rotation couplings, reorientational correlation times). Then the results are presented and discussed in the light of other simulation or experimental results, emphasizing the temperature effect at 1 bar and the pressure effect. The translation properties are described in Section 3, the rotational diffusion coefficients and the reorientational correlation times are presented in Sections 4 and 5, re-

RI PT

spectively, and the translation-rotation diffusion constants is the topic of Section 6.

SC

2. Technical details 2.1. Molecular dynamics simulations

Molecular dynamics simulations are performed using periodic boundary con-

NU

ditions (PBC) with N =512 rigid molecules in the central cubic box interacting through the non polarizable SPC/E intermolecular potential [27]. The time step

MA

is 1 fs.

Equations of motion are solved using center of mass degrees of freedom and a leapfrog quaternion algorithm due to Svanberg [30]. It is derived from an algorithm proposed by Fincham [31] and improves the conservation of energy

D

significantly. In this version, the quaternion vector Q(t+Δt/2) for each molecule

PT E

is obtained iteratively. In our implementation of the algorithm, iterations are stopped when the norm of its variation between iterations i and i + 1 is lower than 10−9 , leading to convergence after about 4 iterations on average.

CE

Electrostatic interactions are computed using lattice summations based on the Ladd approach, initially proposed for assemblies of point dipoles [32, 33, 34]. The generalization to systems composed of point charges [35, 36] leads to

AC

expressions equivalent to the ones proposed by different groups [37, 38]. Both Ladd and Ewald approaches are physically equivalent when the same boundary conditions are applied. In our simulations, a reaction field contribution is added with the dielectric constant of the surrounding continuum equal to infinity. Thus, the Ladd approach becomes equivalent to the standard Ewald approach [34]. The expansion of the electrostatic interaction energy in cubic harmonics

4

ACCEPTED MANUSCRIPT

has been truncated at the rank 10, which has been proven to be sufficient to get good accuracy [34, 39]. The van der Waals interactions are truncated beyond a molecular-based cutoff of 11.3 ˚ A. Usual corrections on potential energy and pressure due to this cutoff have been applied [40]. The system has been studied at ambient pressure and 10 temperatures in the

RI PT

range 210-373 K. Pressure effects have been addressed by modeling the system at 400, 1000, and 3000 bar and nine temperatures in the range 300-700 K and in the pressure range 1-6000 bar at 277 K. The ambient pressure isobar and the isotherm at 277 K have already been considered in the literature [41] for the

SC

SPC/E water model. These new simulations at these thermodynamical states are simply intended to focus on the rotational diffusion coefficients and the

NU

translation-rotation coupling. Nevertheless, translational diffusion will also be discussed in that cases for completeness and comparison with previous studies. Each simulation, labeled by pressure P and temperature T , is obtained in

MA

three stages. During the first one, the system is simulated during 10 5 time steps at 10 densities with scaling of velocities and angular momenta every 100 steps in order to impose the chosen temperature. The pressure is computed and used

D

to guide the changes of density. Then, the density is fitted as a function of pressure by a polynomial of rank 3 to get the value of the density at the desired

PT E

pressure. The second stage consists of 5 ∙ 105 thermalization time steps under the same conditions, using the optimized value of the density. Figure S1 of the Supplementary Material presents our results for the density of liquid water at

CE

1 bar in the temperature range 210-373 K compared to results of the literature. For each simulation, the last stage is a production run of 2 ∙ 106 time steps at constant energy. The imposed energy is the average energy obtained in the last

AC

quarter of the previous thermalization stage. In order to correct the small drift in energy due to finite time step and cutoff of Lennard-Jones interactions, a rescaling of velocities and angular momenta is performed every 10 4 time steps

in order to impose energy conservation. Thus, along the 2 ∙ 106 time steps, the remaining energy drift with respect to the prescribed energy does not exceed 0.35% and 0.65% for simulations below and beyond 400 K, respectively. Tables 5

ACCEPTED MANUSCRIPT

S1-S2 of the Supplementary Material contain the values of density and imposed energy E of each simulation run.

2.2. Diffusion tensors

RI PT

In our study, the molecular frame xyz is defined with the z-axis in the opposite direction of the molecular electric dipole moment and the x-axis parallel to the H-H direction. The laboratory frame XY Z is bound to the simulation box. When necessary, we will distinguish between properties expressed in the

SC

molecular and laboratory frames by using ’mol’ and ’lab’ labels, respectively. The correlation function approach is used to calculate the translational [42,

NU

43, 44] and rotational [45, 46, 47] self-diffusion coefficients. Components of the corresponding tensors are defined by: Z ∞ Dtrans,ij = hvi (0).vj (t)i dt

(1)

MA

0

and

Drot,ij =

Z

∞

0

hωi (0).ωj (t)i dt

(2)

D

where vi and ωj are i and j components of the translational and rotational ve-

PT E

locity of a given molecule, respectively. In the laboratory frame, for an isotropic liquid, both tensors are diagonal with three equal components. In the molecular frame, anisotropy of the diffusion with respect to molecular axes can be revealed.

AC

CE

In the laboratory frame, the diffusion coefficients Z 1 ∞ Dtrans = hv(0).v(t)i dt 3 0

and

Drot

1 = 3

Z

∞ 0

hω(0).ω(t)i dt

are the third of the trace of these diffusion tensors. Brackets denote a statistical average over the molecules and the time origins. In the case of water, in the

6

ACCEPTED MANUSCRIPT

molecular system of axes, one obtains three different translation diffusion coefmol mol ficients Dtrans,αα and three different rotational diffusion coefficients Drot,αα (α

being one of the x, y and z axes of the molecular frame). Coupling between translation and rotation can be defined by computing the

D trans,i−rot,j =

Z

∞ 0

RI PT

cross-correlation self-diffusion tensor: hvi (0).ωj (t)i dt.

(3)

In the laboratory frame, all the elements of this tensor are zero. However, in the molecular frame, according to the C2v symmetry of the water molecule and

SC

to our choice of axis, both components D trans,y−rot,x and D trans,x−rot,y are non zero [48].

In our simulations, the correlation functions are integrated over 6.1 ps. Cal-

NU

culation of diffusion coefficients in laboratory and molecular systems of axes is discussed in the Supplementary Material (Section 2 including Fugure S2).

MA

In addition to this approach relying on velocity and cross-velocity autocorrelation functions (VAC), Dtrans has also been determined using the mean square

(4)

D

displacement (MSD) and the Einstein relationship:

 |ri (0) − ri (t)|2 Dtrans = lim t→∞ 6t

PT E

where ri (t) is the position of molecule i at time t. Practically the slope of this curve has been obtained using the time interval 100-200 ps of the MSD averaged over the molecules and 1.8 105 time origins.

CE

2.3. Calculation of reorientational correlation times In simulation, the molecular reorientational correlation times are obtained

AC

from the single molecule reorientational correlation functions Cl (t) = hPl (u(0).u(t))i

(5)

where Pl (x) is the Legendre polynomial of order l and u(t) is a unit vector

bound to the molecule at time t. The corresponding correlation times are given by τl =

Z

∞ 0

7

Cl (t)dt.

(6)

ACCEPTED MANUSCRIPT

The reorientational correlation times τ1 and τ2 have been computed for four different vectors: the OH-bond vector, the molecular dipole vector (along the z axis), the HH-vector (x axis) and the y axis vector (perpendicular to the molecular plane). To analyze rotational diffusion, these quantities are generally computed by simulation rather than the rotational diffusion coefficients

RI PT

themselves. We will discuss the validity of this substitution in the light of our results.

SC

3. Translational diffusion coefficients 3.1. Temperature effect at 1 bar

Figure 1a presents the evolution of the translational diffusion coefficients

NU

with temperature and a comparison with experimental results [49, 50, 51, 52, 53, 54, 55] and with the TIP4P-2005 model from simulations with N = 1000

MA

molecules in PBC[4]. Note that this simulation study of the TIP4P-2005 model has been performed between 210 and 310 K and the authors have published analytical formula fitting the various diffusion coefficients within this temperature range. For comparison with our simulations in the range 210-373 K, we have

D

extrapolated the diffusion coefficients of the TIP4P-2005 model in the range

PT E

310-373 K using the published fitted functions. Thus, our comparison between SPC/E and TIP4P-2005 between 310 and 373 K is only indicative because it is related to the robustness of the fit and does not necessarily reflect the exact be-

CE

haviour of the TIP4P-2005 model. This remark applies also to the comparison of rotational diffusion coefficients and of translation-rotation couplings discussed in the next sections.

AC

The translational diffusion coefficient for a system of infinite size has been

estimated from the equation[56, 57]: lab lab D∞ = DN +

2.837297kB T 6πηL

(7)

lab where DN and η are the translational diffusion coefficient (in the laboratory

frame) and the shear viscosity for the system of N molecules in PBC, kB the 8

ACCEPTED MANUSCRIPT

Boltzmann’s constant, T the temperature, and L the edge of the cubic box for lab the system of N molecules. This correction has been applied to DM SD . The

shear viscosity has been obtained from the Green-Kubo approach, by integrating the six off-diagonal elements of the stress-autocorrelation function[44] over 6.1 lab ps. In the investigated temperature range, the corrective term added to DN to

RI PT

lab get D∞ amounts to 10-11%, except at 210 K (40%) and 230 K (14%). Note that

our diffusion coefficients obtained for N = 512 from VAC and MSD methods coincide nearly exactly as can be seen from Figure 1a.

˚2 .ps−1 for N = 512 and 0.276 A ˚2 .ps−1 for At 298 K, our result is 0.250 A

SC

an infinite size system. Comparison with a compilation of literature results [3, 5, 19, 20, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70], for the SPC/E

NU

model at 298-300 K as a function of the system size (N −1/3 ) is presented in Figure S3 of the Supplementary Material. One can point out that the published values are somewhat scattered and dependent on simulation features.

MA

Globally, compared to the experiment, the SPC/E model slightly overestimates the diffusion, especially at the lowest temperatures and is rather good beyond ≈325 K. On the contrary, the TIP4P-2005 model is very good at low

D

temperature, but seems to depart significantly from the experimental values beyond ≈325 K. This is not surprising remembering that this model has specially

PT E

been designed to improve the description of the low temperature region of the phase diagram of water.

CE

Fitting the diffusion coefficients by a Speedy-Angell law[71]  γ T −1 D = D0 Ts

(8)

gives (D0 /˚ A2 .ps−1 , Ts /K, γ) parameters equal to (1.0694, 201.401, 1.9744),

AC

(1.0833, 201.877, 1.97093) and (1.16191, 199.896, 2.01737) for Dtrans,V AC ,

Dtrans,M SD , and Dtrans,∞ in the laboratory frame, respectively. By comparison, the corresponding parameters are (1.7176, 218.341, 1.99572) using sets of experimental values of Figure 1a and (1.0074, 220.76, 1.50306) for the TIP4P2005 model using data from Ref. [4]. These fits consistently show that the arrest temperature Ts for TIP4P-2005 (220.8 K) is in better agreement with the ex9

ACCEPTED MANUSCRIPT

perimental value (218.3 K) than for SPC/E (199.9 K). However, the evolution with temperature governed by the power γ (1.996 experimentally) is better with SPC/E (2.017) than with TIP4P-2005 (1.503). In order to allow a direct comparison with the TIP4P-2005 model, the temperature dependence of the diffusion coefficients has also been fitted to the

RI PT

following functional form used by Rozmanov and Kusalik[4]:   B D(T ) = A exp + C. T − T0

(9)

We used the same fitting strategy as these authors. In order to improve the

SC

quality of the fits by preserving balance between low-temperature and hightemperature points, the fitted function has been chosen to be 20D(T )+ln(D(T ))

NU

rather than D(T ). The results are gathered in Table 1. This functional form describes the diffusion of a glass forming liquid, interpreting T0 as the temperature of diffusion quenching. Again, the comparison with TIP4P-2005 shows that the

MA

dynamical arrest temperature T0 , which is around 175 K for TIP4P-2005 (in

AC

CE

PT E

D

laboratory frame) [4], takes lower values (around 100 K) for the SPC/E model.

10

ACCEPTED MANUSCRIPT

B

C

T0

(˚ A2 .ps−1 )

(K)

(˚ A2 .ps−1 )

(K)

lab Dtrans,V AC

14.8379

−0.003676

102.30

0.003911(1.2)

lab Dtrans,M SD

−797.29

18.9024

−894.18

−0.006659

90.52

0.004790(4.7)

22.358

mol Dtrans,V AC

−920.36

6.027

mol Dtrans,xx,V AC

−524.75

6.010

mol Dtrans,yy,V AC

−548.94

7.420

mol Dtrans,zz,V AC

−504.88

4.694

−530.17

A

B

0.005182(4.5)

0.004856

130.94

0.001898(1.1)

0.003418

128.82

0.001741(1.4)

0.006259

132.32

0.003147(1.0)

0.004867

131.29

0.001122(1.0)

C

T0

RMSD(RMSD%)

(K)

(rad2 .ps−1 )

(K)

−1274.74

0.013922

38.92

0.002554(1.4)

−1182.37

0.014464

47.43

0.003262(1.5)

17.299

−1009.44

0.014759

64.69

0.002724(1.4)

12.564

−1245.59

0.019037

33.99

0.001294(1.1)

47.768

−1313.12

0.009150

37.39

0.006031(1.9)

A

B

C

T0

RMSD(RMSD%)

30.331

mol Drot,V AC

24.272

MA

lab Drot,V AC

mol Drot,yy,V AC

PT E

D

mol Drot,zz,V AC

(˚ A.rad.ps

mol Dtrans,x−rot,y,V AC

CE

mol Dtrans,y−rot,x,V AC

SC

88.12

(rad2 .ps−1 )

mol Drot,xx,V AC

RMSD(RMSD%)

−0.006067

NU

lab Dtrans,∞

RI PT

A

−1

)

(K)

−0.1846 −1.6654

(˚ A.rad.ps

−1

)

(K)

−269.98

−0.0000487

195.31

0.002612(10.2)

−505.18

−0.0004578

137.94

0.001499(2.1)

Table 1: Fits of diffusion coefficients of the SPC/E water model at 1 bar using Eq. (9)

AC

in the temperature range 210-373 K. Simulation results with N = 512 molecules. RMSD and RMSD% are the root mean square deviation and root mean square percentage deviation between the simulation results and the fitted values, averaged over the 9 studied temperatures (RMSD has the same unit as A).

mol Figure 1b presents an Arrhenius plot of the three components of Dtrans , and lab mol the average in both frames, Dtrans and Dtrans . As it has already been observed,

11

ACCEPTED MANUSCRIPT

there is a marked non-Arrhenius behavior at low temperature (below ≈298 K)

mol mol and a marked anisotropy of translational diffusion with Dtrans,yy > Dtrans,xx >

mol Dtrans,zz at all temperatures. The behavior is similar to what has been ob-

served for the TIP4P-2005 model (see Figure 1c, where results from Ref. [4] are reported for comparison). Above ≈298 K, however, the temperature evolution

RI PT

of Dtrans follows relatively well the Arrhenius behavior D = D0 e−ΔE/kB T and results of such fits are given in Table 2. We find an activation energy ΔE of 3.3 kcal/mol for the SPC/E model at 1 bar. The increasing divergence at the lab mol lowest temperatures between Dtrans and Dtrans in Figure 1b should be consid-

SC

ered carefully. The absolute values are really small and should not be given any physical meaning. Indeed, logarithmic scale overemphasizes the differences be-

NU

tween both curves. The same behavior has also been reported for TIP4P-2005

AC

CE

PT E

D

MA

(see Figure 1c).

12

ACCEPTED MANUSCRIPT

N→∞

0.8

0.6

2

Dtrans / Å .ps

-1

N=512

0.4

0.2

250

300

350

T/K

1

(b)

Dtrans, mol yy xx

1

aver. yy Dtrans, mol xx

2

Dtrans / Å .ps

-1

aver.

-1

0.1

0.1

SC

2

Dtrans / Å .ps

400

(c)

zz

zz

0.01

0.01

Dtrans, lab

2.5

3

4

3.5

1000/T / K

4.5

-1

5

2

NU

2

RI PT

(a) 200

2.5

3

Dtrans, lab

4

3.5

1000/T / K

4.5

5

-1

MA

Figure 1: Temperature dependence of the translational diffusion coefficient at 1 bar. Panel (a): experimental values from Refs.[49, 50, 51, 52, 53, 54, 55] (green dotted line). Simulation results with the SPC/E model from MSD with N =512 (magenta dotted line) and from VAC (black lines) with N =512 molecules and infinite size system corrected values N → ∞. Simulation results with the TIP4P-2005 model from VAC with N =1000 molecules[4] (blue lines). Contin-

D

uous lines and dashed lines are obtained in the laboratory and molecular frames, respectively.

PT E

Panel (b): Arrhenius plots of the non-zero components of the translational diffusion tensor in the molecular frame obtained with the SPC/E model and N =512: xx component (black continuous line), yy component (blue dot-dashed line), and zz component (dark-green dashed line). The violet dotted line and the red continuous line represent the average in molecular frame and in laboratory frame, respectively. Panel (c): same as panel (b) for the TIP4P-2005

CE

model with N =1000 (from Ref. [4]). Portions of curves for 1000/T < 3.22 (T >310 K) are extrapolated using analytical functions fitted from simulations results in the range 210-310 K

AC

(Ref. [4] and see Section 3.1).

3.2. Pressure effect Figure 2 presents the evolution with temperature at different pressure of

Dtrans obtained for VAC with N = 512. At moderate temperatures, Dtrans follows reasonably an Arrhenius evolution with temperature. As mentioned before, Table 2 reports such Arrhenius fits of the different determinations of 13

ACCEPTED MANUSCRIPT

translational diffusion coefficients. Up to 3000 bar, one observes a decrease of the activation energy as the pressure increases. This quantity is the smallest for the diffusional translation along the y molecular axis and the largest along the x molecular axis. lab Table 3 gathers the numerical values for DVlabAC and DM SD at 1, 400, 1000

RI PT

and 3000 bar with the error bars. At 300 K, our simulation results and especially DVlabAC present a maximum somewhere between 400 and 3000 bar with sufficiently small error bars to clearly show the presence of a maximum. The small number of pressures investigated does not allow to precisely locate the

SC

maximum. At 277 K, the translation diffusion coefficient seems to present also a maximum somewhere between 1000 and 3000 bar. The effect is subtle but

NU

noticeable. However, the conclusion must be tempered by the importance of the error bars. Figure 3 presents Dtrans at 277 K for pressure up to 6000 bar compared to the experiment; a maximum is located around 1750 bar. Although

MA

the error bars are too large to be absolutely conclusive, the trend corresponds to what has been observed by Bagchi et al.[41]. Though in that previous work, the error bars are also too large to say that the existence of a maximum is

D

reproduced, the trend is the same and their simulation data seem to present a maximum between 1000 and 2000 bar. Noticeably, this maximum is also recov-

PT E

ered in the molecular frame, whatever the axis considered. At 350 K and higher temperatures, this effect has disappeared. To compare these findings with experimental results, one can define at a

CE

given temperature the amplitude of the diffusion increase by the ratio between the maximum diffusion coefficient and the value at 1 bar. From experimental results, this amplitude is in the range 1.03-1.10 at 277 (or 273) K (with a

AC

pressure of the maximum ranging from '500 to '2000 bar) and in the range

1.02-1.08 at 298 K (with a pressure of the maximum in the range '900-'1300

bar) [51, 53, 72, 73, 74, 75, 76]. Our results are around '1.03 at 277 K (with a

maximum at pressure around 2000 bar) and '1.04 at 300 K (with a maximum at pressure around 1000 bar). Qualitatively, the SPC/E model reproduces the increase of translational diffusion induced by pressure. 14

ACCEPTED MANUSCRIPT

DVlabAC

lab DM SD

lab D∞

DVmol AC

DVmol AC,xx

DVmol AC,yy

DVmol AC,zz

65.596

67.334

77.411

29.175

29.683

35.141

23.078

3.284

3.300

3.326

2.765

2.840

2.697

2.794

53.408

56.355

65.385

20.768

20.703

24.978

16.870

3.179

3.221

3.242

2.554

2.613

2.489

2.597

39.461

48.561

55.405

17.775

19.050

20.603

14.298

2.979

3.168

3.175

2.465

2.573

2.381

2.493

24.125

23.722

2.692

2.667

D0 /˚ A2 .ps−1 ΔE/kcal.mol−1

RI PT

1 bar

400 bar D0 /˚ A2 .ps−1 ΔE/kcal.mol−1 D0 /˚ A2 .ps−1 ΔE/kcal.mol−1

SC

1000 bar

D0 /˚ A2 .ps−1

27.626

14.439

15.147

16.828

11.548

2.694

2.379

2.461

2.316

2.388

MA

ΔE/kcal.mol

−1

NU

3000 bar

Table 2: Arrhenius fits of the translational diffusion coefficients of the simulation results with N = 512 molecules with the SPC/E model: at 1 bar in the temperature range 298-373 K; at 400, 1000, and 3000 bar in the range 300-550 K for D lab and in the range 300-400 K for

AC

CE

PT E

D

D mol .

15

ACCEPTED MANUSCRIPT

1 bar

400 bar

1000 bar

3000 bar

277 K

0.156(3)

0.157(1)

0.158(2)

0.156(2)

300 K or 298 K

0.250(3)

0.263(5)

0.259(2)

0.243(2)

350 K or 353 K

0.603(3)

0.574(4)

0.555(4)

0.500(4)

0.998(8)

0.939(6)

0.827(3)

400 K

RI PT

A2 .ps−1 DVlabAC /˚

lab ˚2 −1 DM SD /A .ps

0.157(3)

0.157(3)

0.156(3)

0.152(4)

300 K or 298 K

0.250(5)

0.248(5)

0.263(2)

0.243(3)

350 K or 353 K

0.612(6)

0.558(14)

0.557(4)

0.482(5)

0.993(6)

0.963(23)

0.821(8)

NU

400 K

SC

277 K

Table 3: Pressure effect on translational diffusion coefficient. Simulation results with the

AC

CE

PT E

D

MA

SPC/E model and N = 512 molecules. Absolute error on last digits are given in parentheses.

16

ACCEPTED MANUSCRIPT

-1

8

2

Dtrans / Å .ps

6

2

RI PT

4

(a)

300

400

600

500

T/K

SC

10

700

(b)

1

NU

1

yy xx zz Dtrans, mol

MA

2

Dtrans / Å .ps

-1

Dtrans, lab

1.5

2

3

2.5

1000/T / K

3.5

4

-1

D

Figure 2: Pressure and temperature dependence of the translational diffusion coefficient of the SPC/E water model from VAC and N =512 molecules. Panel (a): diffusion coefficient in

PT E

laboratory frame vs. temperature : 400 bar (black continuous line), 1000 bar (blue dot-dashed line), and 3000 bar (dark-green dashed line). Panel (b): Arrhenius plots of components in molecular frame. Thin lines : 400 bar; medium lines : 1000 bar; thick lines : 3000 bar. Black continuous line : xx; blue dot-dashed line : yy; dark-green dashed line : zz. The continuous

AC

CE

red line represents the average in laboratory frame.

17

ACCEPTED MANUSCRIPT

0.24

0.2

2

Dtrans / Å .ps

-1

0.22

0.18 0.16

0

(a)

1000

2000

3000

4000

P / bar

-1

0.2

2

Dtrans / Å .ps

6000

(b)

0.22

NU

0.18 0.16 0.14

MA

0.12

0

5000

SC

0.24

0.1

RI PT

0.14 0.12

1000

2000

3000

4000

5000

6000

P / bar

Figure 3: Pressure dependence of the translational diffusion coefficient of the SPC/E model

D

at 277 K. Panel (a) : black continuous line : from VAC, N =512, in the laboratory frame.

PT E

Black dashed line : from VAC, N =512, in the molecular frame. Blue continuous line : from MSD, N =512. Blue dashed line: from MSD, N→ ∞. Green dotted line : experiment at 278 K from Ref. [76]. Panel (b) : components in the molecular frame. Black continuous line : xx;

CE

blue dot-dashed line : yy; dark-green dashed line : zz; violet dotted line : average.

4. Rotational diffusion coefficients

AC

4.1. Temperature effect at 1 bar Figure 4 presents the rotational diffusion coefficients obtained from angu-

lar velocity autocorrelation functions. Figure 4a shows a comparison with the TIP4P-2005 model which gives rotational diffusion coefficients smaller than the SPC/E ones. Figure 4b shows the evolution with temperature of the ratio Drot /Dtrans , evaluated in the laboratory system of axes. This ratio remains

18

ACCEPTED MANUSCRIPT

roughly constant at high temperature, but increases significantly below 250 K, sign of a decoupling between translation and rotation. This curve is qualitatively similar to the result obtained by Mazza et al. [28] with the same potential, but the numerical values are quite different. However, this curve is in contradiction with the curve obtained from experimental data using an empirical approximate

RI PT

expression of Drot [77].

The anisotropy of the rotational diffusion coefficient is substantial (Figmol ure 4c), whatever the temperature. At 298 K, one finds components (Drot,xx , mol mol Drot,yy , Drot,zz ) equal to (0.244, 0.132, 0.322) in rad 2 .ps−1 . In the literature,

SC

different sets of rotational diffusion coefficient components have been proposed for the SPC/E model near ambient conditions. For example, the values ob-

NU

tained at 298 K by Svishchev and Kusalik[3] are equal to (0.46, 0.22, 0.45) and (0.211, 0.114, 0.272) as obtained by Chevrot et al.[5]. At 300 K and a density of '0.9756 g/cm3 (2.3% lower than our density), Zasetsky et al.[64] obtained

MA

three rotational diffusion coefficients equal to 0.30, 0.45 and 0.82 rad 2 .ps−1 .

Our values are in good agreement with the values obtained by Chevrot et al.[5] The difference of '15%, similar to that observed for the translational

D

diffusion coefficient in Figure S3 of the Supplementary Material, remains unexplained and can perhaps be attributed to the differences in simulation details

PT E

(integration of atomic equations of motion with holonomic constraints to rigidify the molecules vs. rigid body quaternion algorithm here, Particle Mesh Ewald vs. Ladd summation for the treatment of electrostatic interactions, notably).

CE

Figures 4a and 4c show that the rotational coefficients computed in the laboratory frame and the average in the molecular frame remain very similar in the temperature range 210-373 K. Moreover, in that range, the average values

AC

and the components in the molecular frame obey an Arrhenius law from the highest temperature down to 250 K. Table 4 contains results of such fits between 273-373 K with an activation energy ΔE equal to 3.1 kcal/mol for the average

Drot and 3.0, 2.8, and 3.2 kcal/mol for the diffusion coefficients around the x, y and z molecular axes, respectively.

19

ACCEPTED MANUSCRIPT

DVlabAC

DVmol AC

DVmol AC,xx

DVmol AC,yy

DVmol AC,zz

45.995

41.007

36.973

13.988

76.891

3.121

3.064

2.971

2.762

3.245

41.807

36.454

33.476

3.058

2.978

2.898

37.656

32.795

30.756

2.994

2.917

2.847

29.674

28.819

2.829

2.826

D0 /rad2 .ps−1 ΔE/kcal.mol−1 400 bar D0 /rad2 .ps−1 ΔE/kcal.mol−1

15.864

67.171

2.843

3.157

12.832

56.625

2.698

3.054

28.781

12.266

46.420

2.773

2.653

2.936

ΔE/kcal.mol−1

SC

1000 bar D0 /rad2 .ps−1

MA

ΔE/kcal.mol

−1

NU

3000 bar D0 /rad2 .ps−1

RI PT

1 bar

Table 4: Arrhenius fits of the rotational diffusion coefficients of the SPC/E water model in the temperature ranges 273-373 K at 1 bar and 300-400 K at 400, 1000 and 3000 bar. Simulation

AC

CE

PT E

D

results with N = 512 molecules.

20

ACCEPTED MANUSCRIPT

8 7

0.6 -2 2

0.4 0.3 0.2

(a) 300

250

3 2 1

400

350

5 4

RI PT

200

(b)

6

Drot/Dtrans / rad .Å

2

Drot / rad .ps

-1

0.5

200

T/K

(c) zz

1

2

0.1

2.5

3

4

3.5

1000/T / K

4.5

-1

5

yy

0.1

NU

2

Drot, lab

(d)

SC

Drot / rad .ps

yy

400

350

zz Drot, mol xx aver.

-1

Drot, mol xx aver.

2

Drot / rad .ps

-1

1

300 T/K

250

2

2.5

3

Drot, lab

4

3.5

1000/T / K

4.5

5

-1

MA

Figure 4: Temperature dependence of the rotational diffusion coefficient and of its ratio with the translational diffusion coefficient at 1 bar. Panel (a): SPC/E model (N =512) in black, TIP4P-2005 model (N =1000) in blue (from Ref. [4]). Continuous lines with squares and dashed lines with crosses correspond to results from VAC in the laboratory and molecular

D

frames, respectively. Panel (b): rotational to translational diffusion coefficient ratio, SPC/E model. Panel (c): Arrhenius plots of the non-zero components of the rotational diffusion tensor

PT E

of the SPC/E model in the molecular frame with N =512: xx component (black continuous line), yy component (blue dot-dashed line), and zz component (dark-green dashed line). The violet dotted line and the red continuous line represent the average in molecular and laboratory frames, respectively. Panel (d): same as panel (c) for the TIP4P-2005 model (from Ref. [4]).

CE

For TIP4P-2005 model in panels (a) and (d), the portions of curves for 1000/T < 3.22 (T >310 K) are extrapolated using analytical functions fitted from simulations results in the range 210-

AC

310 K (Ref. [4] and see Section 3.1).

Comparison between Figures 4c (SPC/E) and 6d (TIP4P-2005[4]) shows

mol that in the temperature range 210-373 K, the SPC/E model leads to Drot,zz >

mol mol Drot,xx > Drot,yy contrary to the TIP4P-2005 model for which an inversion

mol mol mol occurs below ' 240 K where the order becomes Drot,yy > Drot,zz > Drot,xx .

As in the case of translational diffusion, the rotational diffusion coefficients

21

ACCEPTED MANUSCRIPT

have been fitted using the functional form of Eq. (9). The results are given in Table 1. For rotational diffusion coefficients too, the SPC/E model presents an arrest temperature T0 lower than the TIP4P-2005 model (' 40 K vs. ' 110 K lab for Drot , for example).

The temperature dependence of the rotational diffusion with the SPC/E

RI PT

model has been studied previously by MD simulations with PBC at a density of 1 g/cm3 between 200 and 350 K by Mazza et al.[28, 78] for a system of N =1728 molecules and at densities corresponding to ambient pressure between 222 and 356 K for a system of N =256 molecules by Galamba[58]. This last

SC

author characterizes the rotational diffusion by correlation times τl associated with the Legendre polynomial of rank l = 1, 2, 3 of the reorientational time-

NU

correlation function of the OH bond normalized vector and makes correlation with the diffusion coefficient using the Debye model (Eq. (11), Section 5.3), emphasizing, however, the limitations of this relation. We will come back later

Einstein-type approach as

MA

to this statement. Mazza et al. define rotational diffusion coefficients by an

 1

Δϕ(t)2 t→∞ 4t

D

Drot = lim

(10)

PT E

with Δϕ(t) being an unbounded angular displacement of one axis bound to the molecule during a time interval t. With this definition, the rotational diffusion coefficient obtained from the reorientation of a given molecular axis is related to the components of the rotational diffusion tensor corresponding to rotations

CE

around both other perpendicular axes. At 300 K, these authors obtain three rotational diffusion coefficients with numerical values '0.090, '0.115, and '0.080

AC

rad2 /ps from the reorientation of molecular x, y and z axes, respectively (using

our definition of inertia axes). Our values for Drot,αα (α = x, y, z) are significantly larger. A possible explanation of the discrepancies can be related to the fact that the use of Eq. (10) to get rotation is particularly delicate, as noted by Laage and Hynes[79]. By contrast, the use of Eq. (2) is straightforward.

22

ACCEPTED MANUSCRIPT

4.2. Pressure effect Figure 5 shows the pressure effect on Drot at 277 K up to 6000 bar and completes the study of Bagchi et al.[41] who analyzed the rotational diffusion

RI PT

0.2

2

Drot / rad .ps

-1

0.25

0.15

0

1000

2000

3000

5000

6000

MA

-1 2

Drot / rad .ps

4000

NU

P / bar

SC

(a)

0.2

0.1

D

(b) 1000

PT E

0

2000

3000

4000

5000

6000

P / bar

Figure 5: Pressure dependence of the rotational diffusion coefficient of the SPC/E model at 277 K from VAC with N =512. Panel (a): in the laboratory frame (black continuous line with

CE

squares) and in the molecular frame (black dashed line with crosses). Panel (b): components in the molecular frame, xx (black continuous line), yy (blue dot-dashed line), zz (dark-green

AC

dashed line), and average (violet dotted line).

in terms of reorientational correlation times only. It is found that Drot

slightly increases from 0.159 to 0.182 rad 2 .ps−1 between 1 and 6000 bar (in

the laboratory frame). Bagchi et al. interpreted this behaviour by an analysis of the structure showing that the pressure increases the presence of interstitial water molecule in the first solvation shell of the molecules. Interstitial molecules

23

ACCEPTED MANUSCRIPT

have been previously observed in water with the SPC/E model by Svishchev and Kusalik [80]. This pressure induced structural change leads to a weakening of hydrogen bonding and an increase of the diffusion. We also observed it considering Dtrans . At 277 K, when the pressure increases up to '4000 bar, mol mol the components of Drot increase. Between 4000 and 6000 bar, Drot,zz remains

RI PT

nearly constant and Drot,yy increase very little. This is probably the first sign that the broad maximum is nearly reached and Drot will soon start to decrease at additional pressure increase.

lab Figure 6a presents the evolution of the rotational diffusion coefficients Drot

SC

with temperature at 400, 1000 and 3000 bar and, at the scale of the plot, illustrates clearly that, at sufficiently high temperature (for example, 400 K), the

NU

rotational diffusion decreases with an increase of pressure. This is also apparent from Figure 6b presenting Arrhenius plots of the three rotational diffusion comol efficients Drot,αα (α = x, y and z). Figure 6c compares the rotational diffusion

MA

coefficients obtained in the laboratory and molecular frames. Up to ' 450 K, both determinations are close to each other. Figure 6d shows an Arrhenius plot lab of Drot below 400 K at the four considered pressures. This plot shows that,

D

at the lowest temperatures, the rotational diffusion coefficient increases as the pressure increases and that the inversion temperature is around 338 K. Looking

PT E

mol at Drot , a very similar plot is obtained (not shown here) and the inversion tem-

perature is around 336 K. One can conclude that, starting from atmospheric pressure, the rotational diffusion coefficient of the SPC/E model increases if

CE

pressure increases below 340±5 K, and the opposite behavior occurs beyond that temperature.

Concerning the anisotropy of the rotational diffusion, one can observe in Fig-

AC

ure 6b that, at all investigated pressures and temperatures, this model predicts mol mol mol that Drot,zz > Drot,xx > Drot,yy . mol mol mol mol Figure 6e shows the ratio Drot,xx /Drot,zz and Drot,yy /Drot,zz at different

pressures as a function of temperature. At 300 K and above, the anisotropy is significant. Above 300 K, the anisotropy of rotational diffusion is only slightly mol mol temperature dependent, especially Drot,yy /Drot,zz that remains nearly constant.

24

ACCEPTED MANUSCRIPT

mol mol The ratio Drot,xx /Drot,zz slightly decreases by about 10-15 % between 300 and

700 K. In this temperature range, a pressure increase leads to a small increase of these ratios, meaning a lowering of the anisotropy. Below 300 K, the Hbonds become more and more important. The involvement of a molecule in a strong rigid three-dimensional H-bond network at very low temperature removes

RI PT

the existence of facilitated rotation around a given axis, the rotation becomes

AC

CE

PT E

D

MA

NU

SC

equally difficult around any molecular axis.

25

ACCEPTED MANUSCRIPT

100

20

(b)

(a) zz

10

xx Drot, mol yy

2

Drot / rad .ps

2

Drot / rad .ps

-1

-1

15

10

1

300

400

1

700

600

500

1

(c) Drot, lab -1

Drot, mol

3

2.5

1000/T / K

4

3.5

-1

0.7 0.6 0.5

3

2.5

MA

0.8

mol mol Drot,xx /Drot,zz

(d)

1000/T / K

3.5 -1

(e)

mol mol Drot,yy /Drot,zz

0.4

D

anisotropy of rotational diffusion

1 0.9

NU

1

4

3.5

-1

SC

2

Drot / rad .ps

2

Drot / rad .ps

-1

Drot, lab

2

3

2.5

1000/T / K

10

1.5

2

1.5

T/K

1

RI PT

5

0.3

PT E

200

300

400

500

600

700

T/K

Figure 6: Pressure and temperature dependence of the rotational diffusion coefficient of SPC/E water model from VAC and N =512 molecules. Panel (a): diffusion coefficient in laboratory

CE

frame vs. temperature at 400 bar (black continuous line), 1000 bar (blue dot-dashed line), and 3000 bar (dark-green dashed line). Panel (b): Arrhenius plots of rotational diffusion coefficients in molecular frame: xx component (black continuous lines), yy component (blue

AC

dot-dashed lines), zz component (dark-green dashed lines), and average in the laboratory

frame (red continuous lines). Thin lines : 400 bar; medium lines : 1000 bar; thick lines : 3000 bar. Panel (c): Arrhenius plots of average in the laboratory frame (red continuous lines) and in the molecular frame (violet dotted lines). Thin line : 400 bar; medium : 1000 bar; thick : 3000 bar. Panel (d): Arrhenius plots of average in the laboratory frame for T ≤ 400 K at

1 bar (green thin continuous line), 400 bar (magenta dotted line), 1000 bar (blue dot-dashed line), and 3000 bar (red thick continuous line). Panel (e): anisotropy of rotational diffusion mol /D mol mol mol represented by the ratios Drot,xx rot,zz and Drot,yy /Drot,zz . Temperature dependence at

1 bar (green thin continuous line), 400 bar (magenta dotted line), 1000 bar (blue dot-dashed

26 line), and 3000 bar (red thick continuous line).

ACCEPTED MANUSCRIPT

5. Reorientational correlation times 5.1. Temperature effect at 1 bar Figure 7 presents our results for τ1 and τ2 for the inertia axes and the OH bond vector at 1 bar. Our values at 298 K are 4.31, 2.77, 4.42, 4.34 ps (τ1 ) and

RI PT

1.90, 1.10, 1.47, 1.71 ps (τ2 ) for the x, y, z axes and OH bond vector, respectively. These values are consistent with some determinations from literature for the SPC/E model[3, 64, 79, 81, 82] and differences larger than ± 10 % can be found with others[19, 83]. Compared to experimental results (τ2,x = 2.0 ps[84],

SC

τ1,z = 4.76[85], τ2,z = 1.92 ps[85, 86], τ2,OH = 1.95 ps[87, 88, 89]), our results are in reasonable agreement, with rotational motions slightly too fast with this intermolecular potential.

NU

The temperature dependence of τ2,OH and τ2,y between 275 K and 350 K obtained from NMR relaxation data and from MD simulations with the SPC/E

MA

model for D2 −16 O liquid water have been reported before [90]. Our values for H2 −16 O are consistently smaller (about 10-20 %) than the MD values of that study due to the difference of inertia momenta. Compared with the simulation of Galamba[58] at 1 bar in the temperature range 222-356 K with N = 256

D

molecules, our values are in good agreement except at the lowest temperatures.

PT E

At 230 K, our values for τ1,OH and τ2,OH are smaller by about 20 %. Though our densities are slightly different, the discrepancy reflects mainly the statistical

AC

CE

error.

27

ACCEPTED MANUSCRIPT

10

τ1

RI PT

τ / ps

100

z x OH y x OH z y

1

τ2 200

300 T/K

350

400

SC

250

Figure 7: Reorientational correlation times τ1 (thin lines, upper group) and τ2 (thick lines, lower group) for the inertia axes and the OH bond vector at 1 bar: x-axis (black continuous

NU

line), y-axis (blue dot-dashed line), z-axis (green dashed line), and OH-bond vector (red dotted line) for the SPC/E model.

MA

5.2. Pressure effect

The pressure dependence of correlation times τ1 and τ2 of the inertia axes and −1 OH-bond vector is presented in Figure 8. Paralleling the behavior of Drot , the

D

correlation times decrease when pressure increases at low-enough temperature

PT E

and has the opposite behavior at higher temperature. The turning temperatures are 405±10, 410±15, 465±10, 420±10 K for τ1 and τ2 of the x, y, z and OH−1 bond axis, respectively. As the relationship τl ∝ Drot is not strictly valid, the

turning temperature obtained from the correlations times does not necessarily

AC

CE

have to be identical to the one obtained from Drot .

28

ACCEPTED MANUSCRIPT

100

τ1

1

τ2

0.1

4

100

10

τ1

1

τ2

4

2

3 -1 1000/T / K

4

5

5

τ1

τ2

1

0.1

(c) 3 -1 1000/T / K

(b)

OH bond

10

0.1

2

τ2

1

z axis

1

1

5

τ / ps

τ / ps

100

3 -1 1000/T / K

τ1

0.1

(a) 2

10

SC

1

y axis

RI PT

10

τ / ps

x axis

NU

τ / ps

100

1

2

(d) 3 -1 1000/T / K

4

5

MA

Figure 8: Pressure and temperature dependence of the reorientational correlation times τ1 and τ2 for the inertia axes and the OH-bond vector at 1 bar (green thin continuous line), 400 bar (magenta dotted line with circles), 1000 bar (blue dashed line), and 3000 bar (red thick continuous line) for the SPC/E model. Panel (a): x axis; panel (b): y axis; panel (c): z axis;

D

panel (d): OH-bond vector.

PT E

Though we did not explore sufficiently high pressures to locate the minimum in the curves of correlation times vs. pressure (or the maximum for rotational diffusion coefficient vs. pressure), at 277 K and 6000 bar the system is close

CE

to the broad minimum that Bagchi et al. have found. Defining the amplitude of the pressure effect by the ratio between a correlation time at 1 bar and the correlation time at 6000 bar, one obtains an amplitude of 1.33 (using τ1 ) or

AC

1.34-1.37 (using τ2 ) at 277 K. At 300 K, as the pressure has been limited to

3000 bar, one can just say that the amplitude is larger than 1.16-1.18 from τ2 , for example. Examining the amplitude of the pressure effect by the ratio

between the maximum rotational diffusion coefficient and the value at 1 bar, one gets quantitatively different results: ≈1.15 at 277 K and >1.09 at 300 K. From experiments, defining the amplitude by ratio of τ2 correlation times 29

ACCEPTED MANUSCRIPT

or inverse ratio of longitudinal relaxation times, values of 1.29 [75], 1.26 [91], 1.16 [92] at 272-273 K, and 1.12 [93], 1.10 [94], 1.26 [91] at 283 K have been published. The pressure corresponding to the extremum for rotational diffusion is generally found experimentally around 2000 bar at these temperatures. The SPC/E model reproduces the amplitude of this pressure effect for rotational

RI PT

diffusion reasonably well.

Wakai and Nakahara[95] have investigated the NMR quadrupolar relaxation of deuterium in liquid D 2 O from 283 K to 323 K up to a pressure of 3000 bar. At each pressure, they observe an obvious decrease of the correlation time τ2

SC

of the OD-bond vector when the temperature increases. At each temperature, they obtained a monotonous decrease of the correlation time when the pres-

NU

sure increases. The SPC/E model agrees with the experimental observation. Correcting our values for isotopic effects, one can conclude that in these temperature and pressure ranges the SPC/E model predicts correlation times about

MA

10% smaller than experiment. This too fast rotational diffusion can be correlated with the observation that this model also overestimates the translational diffusion, compared to the experiment. The too fast dynamics of this model

D

had already been pointed out[19, 20].

PT E

5.3. Test of the Debye model

Within the small-step rotational diffusion model (Debye model), the correlation function for the reorientation of a vector bound to a molecule Cl (t) =

CE

hPl (u(0).u(t))i decreases exponentially and τl =

1 . l(l + 1)Drot

(11)

AC

Though this model is known not to be quantitative, because large-angle reorientation can occur and inertial effects have to be taken into account, this simple model is often used to analyze experimental or simulation results. In order to test this model, Figure 9 presents plots of the ratio of the rotational coefficient Drot (Debye) = 1/(6τ2,axis ) and the ’exact’ rotational diffusion

30

ACCEPTED MANUSCRIPT

coefficient obtained from our simulations. As it is well-known, this clearly illustrates that the model is not satisfactory, even under ambient conditions. As expected, at low temperature, the importance of large-jump reorientation increases and the applicability of the model deteriorates. At high temperature, inertial effects are important and not taken into account by the model. In addi-

RI PT

tion, an increase of pressure plays a role opposite to an increase of temperature. This Figure shows that the error on Drot resulting from the application of the Debye model would be at least 25-30 % (using the y axis reorientation) and

SC

50-55 % (using the OH-bond vector) in the range 300-400 K.

NU

1.4

Debye model valid

1

y axis

0.8

MA

1/(6 τ2,axis Drot)

1.2

0.6

OH axis

0.4

300

400

500

600

700

T/K

PT E

200

D

0.2

Figure 9: Test of the Debye model for the reorientation of OH-bond vector and the y axis perpendicular to the molecular plane. If the Debye model is valid, Drot (Debye) = 1/(6τ2,axis ).

CE

lab at different pressures: 1 bar (green thin conTemperature dependence of Drot (Debye)/Drot

tinuous line), 400 bar (magenta dotted line), 1000 bar (blue dashed line), and 3000 bar (red thick continuous line).

AC

Another possible test of the model is the ratio τ1 /τ2 (not shown here), which

should be equal to 3 if Debye model is valid. From our simulation data, for any axis, this ratio takes a maximum value of about 2.55 around 350 K and its evolution with temperature and pressure follows a trend similar to the plots of Figure 10. At the lowest temperatures, the ratio reaches values around 2.252.30. At 700 K, τ1 /τ2 for y axis decreases to '2.1 at 400 bar and '2.25 at 3000 31

ACCEPTED MANUSCRIPT

bar. For the OH-bond vector, this ratio is '2.3 at 400 bar and '2.35 at 3000 bar. These results confirm if necessary that the Debye model is not accurate to determine rotational diffusion coefficients from correlation times. This is important to remind when comparing experimental results obtained with this

RI PT

model and simulation results, for instance. Beyond the value of the rotational diffusion coefficient itself, it could lead to wrong tendencies vs. temperature or pressure.

SC

5.4. Anisotropy of reorientational motions

This property is illustrated in Figure 10 presenting the ratios τl,x /τl,z and

NU

τl,y /τl,z for l = 1 and 2. A feature worth noticing is that τl,x and τl,z remain markedly larger than τl,y at any temperature and pressure. In addition to Figure mol mol mol mol 6e (showing the ratio Drot,xx /Drot,zz and Drot,yy /Drot,zz ), these plots illustrate

MA

the importance of anisotropy in rotational diffusion of water, well known from experiment [90]. The anisotropic diffusion model[96, 97] taking into account the

AC

CE

PT E

D

correlated motions of the three inertia axes would tackle these features.

32

1

(a)

τ1,x / τ1,z

0.8

τ1,y / τ1,z 0.6 200

300

400

500

600

T/K

700

SC

1.4

(b)

τ2,x / τ2,z

NU

1.2 1

τ2,y / τ2,z

0.8 0.6 200

MA

anisotropy of reorientational motion

1.6

RI PT

anisotropy of reorientational motion

ACCEPTED MANUSCRIPT

300

400

500

600

700

T/K

Figure 10: Anisotropy of reorientational motion of the SPC/E model represented by the

D

ratios τl,x /τl,z and τl,y /τl,z . Temperature dependence at different pressures: 1 bar (green

PT E

thin continuous line), 400 bar (magenta dotted line), 1000 bar (blue dashed line), and 3000 bar (red thick continuous line). Panel (a): l = 1; panel (b): l = 2.

CE

6. Translation-rotation coupling As already mentioned, for water with our choice of molecular axes, the only non-zero components of the translation-rotation diffusion tensor (expressed

AC

in the molecular system of axes) are the translation-x / rotation-y and the

translation-y / rotation-x ones, if the C2 symmetry axis is the z-axis. We will label trans, x − rot, y and trans, y − rot, x these translation-rotation cou-

plings. Such components are sensitive to the choice of molecular system of axes. Chevrot et al.[5] put the molecule in the (y, z)-plane and the correspondence with our choice requires to make the substitution (x → y, y → −x, z → z). 33

ACCEPTED MANUSCRIPT

Kusalik and coworkers[4, 29] put the molecule in the (x, z)-plane with the z-axis having a direction opposite to our choice. The correspondence with our choice requires the substitution (x → x, y → −y, z → −z). Consequently, our values are all negative when Chevrot et al. and Kusalik et al. obtain positive values and our trans, x − rot, y and trans, y − rot, x values are exchanged with respect

RI PT

to Chevrot et al.’s results. In Figure S4 of the Supplementary Material, the translation-rotation cross-correlation functions at ambient temperature and 1 bar and 3000 bar are presented.

SC

6.1. Results at 1 bar

mol mol At 298 K, our values for Dtrans,x−rot,y and Dtrans,y−rot,x are −0.0126 and

−0.0719 ˚ A.rad.ps−1 . Chevrot et al. obtain (with our axis definition) −0.0104

NU

and −0.0630 ˚ A.rad.ps−1 , respectively. Figure 11a shows the evolution of these

mol coefficients with the temperature. Dtrans,y−rot,x has the largest value and a

MA

stronger temperature dependence. The coupling diminishes when the temperature decreases and tends to vanish around 200 K. This is consistent with the results of the ratio Drot /Dtrans . The behavior at low temperature is consistent

D

with the decoupling of translation and rotation. mol The cross-correlation translation-rotation diffusion coefficients Dtrans,α−rot,β

PT E

have been fitted as a function of the temperature using Eq. (9) and the results are given in Table 1. Note that for the these transport coefficients D(T ) which are negative, the strategy described in subsection 3.1 had to be modified and

CE

the fitted function was −20D(T ) + ln (−D(T )). The comparison with the TIP4P-2005 model shows very similar numerical mol mol values for Dtrans,y−rot,x (T ), but values of Dtrans,x−rot,y (T ) about 40-60 % larger

AC

for the SPC/E model. 6.2. Pressure effect The effect of pressure on the translation-rotation coupling is illustrated in Figures 11b and 11c. This effect is moderate. At 300 K or at lower temperatures, there is nearly no pressure effects up to a few thousands of bar. The isotherm

34

ACCEPTED MANUSCRIPT

at 277 K, where pressures up to 6000 bar have been investigated, confirms this trend. Up to 550 K, a pressure increase of 3000 bar diminishes (in absolute mol value) the effect of the coupling by 5-10 % for Dtrans,y−rot,x (T ) and by 20-30 % mol for Dtrans,x−rot,y (T ). This observation, linked to the alteration of the hygrogen

[93, 94].

0 0 -0.1

SC

-1

Dtran,α-rot,β / Å.rad.ps

-0.1

trans,y-rot,x

-0.2

(a) 300 T/K

250

0

(b)

300

trans,x-rot,y

D

-0.02

350

400

450

500

550

T/K

(c)

-0.03

-0.04

PT E

Dtrans,α-rot,β / Å.rad.ps

-1

-0.01

trans,y-rot,x

-0.5

-0.6

400

350

-0.3

-0.4

MA

200

trans,x-rot,y

-0.2

NU

Dtrans,α-rot, β / Å.rad.ps

-1

trans,x-rot,y

RI PT

bond network due to pressure, is consistent with previous experimental works

trans,y-rot,x

-0.05

0

1000

2000

3000 P /bar

4000

5000

6000

CE

Figure 11: Translation-rotation diffusion coefficient components of the SPC/E model in the molecular system of axes: translation-x / rotation-y (green lines) and translation-y / rotationx (blue lines) components. Panel (a): at 1 bar and 210 ≤ T ≤ 373 K. Panel (b): at 400 bar

AC

(thin red continuous lines), 1000 bar (medium blue dot-dashed lines) and 3000 bar (thick orange continuous lines) and 300 ≤ T ≤ 550 K) Panel (c): at T = 277 K and 1 ≤ P ≤ 6000

bar.

35

ACCEPTED MANUSCRIPT

7. Conclusion This work has presented an extensive study of the diffusional properties of the popular SPC/E model of water in the liquid phase at 1 bar (210-373 K), 400, 1000, 3000 bar (300-700 K), and from 1 to 6000 bar at 277 K. In spite

RI PT

of numerous studies already devoted to this water model, this work appears to give for the first time an extensive accurate estimation of the rotational diffusion coefficients and of the translation-rotation coupling over a relatively wide range of temperatures and pressures.

SC

When available, we compared our results to experimental data and TIP4P2005 model predictions. We observed that this last potential is more accurate at low temperature, but SPC/E remains competitive at higher ones.

NU

Our results analyzed in detail the anisotropy of translational and rotational dynamics in the molecular frame, what is especially interesting for the rotational

MA

diffusion because it can be obtained experimentally whereas no experiment can analyze the anisotropy of translational diffusion. We have analyzed the anomalous behaviour of increased diffusion induced by pressure, due to the weakening of hydrogen bonds, and found that the SPC/E

D

model can describe this effect. The amplitude of this effect is larger for the

PT E

rotational diffusion than for translational diffusion and thus is more accurately observable from simulation of the rotational diffusion coefficient than from the translational diffusion coefficient. With the SPC/E model, the anomalous behaviour disappears for the translational diffusion coefficient at a temperature

CE

between 300 and 350 K and around 340 K for the rotational diffusion coefficient. Our results also illustrate in a systematic manner that using the Debye model

AC

to determine the rotational diffusion coefficient from correlation times leads to wrong results and wrong tendencies as well. It has also been observed that earlier determination of the rotational diffusion coefficients of the literature are not absolutely reliable and the present values improve these former values. The translation-rotation coupling coefficients have also been computed for the first time with this model over a wide temperature and pressure range. Our

36

ACCEPTED MANUSCRIPT

results show a decoupling between rotational and translational diffusion as the temperature decreases. Such a behavior is also expected at very high pressures in plastic phases.

RI PT

Supplementary Material This work is accompained by a Supplementary Material file containing five Sections including four Figures (Figure S1 to S4) and two Tables (Table S1 and

SC

S2).

Acknowledgments

NU

We are grateful to the Universit´e de Lorraine and the CNRS for their support. We particularly thank the Institut Jean Barriol de Chimie et Physique Mol´eculaire et Biomol´eculaire (FR 2843) for a specific support through a gen-

MA

erous special funding. We also thank the Lorraine M´esocentre EXPLOR for providing us with computer time.

D

References

PT E

[1] Y. Takii, K. Koga, H. Tanaka, A plastic phase of water from computer simulation, J. Chem. Phys. 128 (2008) 204501. [2] J. L. Aragones, C. Vega, Plastic crystal phases of simple water models, J.

CE

Chem. Phys. 130 (2009) 244504. [3] I. M. Svishchev, P. G. Kusalik, Dynamics in liquid h 2 o, d2 o, and t2 o: A

AC

comparative simulation study, J. Phys. Chem. 98 (1994) 728–733. [4] D. Rozmanov, P. G. Kusalik, Transport coefficients of the tip4p-2005 water model, J. Chem. Phys. 136 (2012) 044507.

[5] G. Chevrot, K. Hinsen, G. R. Kneller, Model-free simulation approach to molecular diffusion tensors, J. Chem. Phys. 139 (2013) 154110.

37

ACCEPTED MANUSCRIPT

[6] A. Laaksonen, P. Stilbs, R. E. Wasylishen, Molecular motion and solvation of benzene in water, carbon tetrachloride, carbon disulfide and benzene: A combined molecular dynamics simulation and nuclear magnetic resonance study, J. Chem. Phys. 108 (1998) 455–468.

RI PT

[7] R. Witt, L. Sturz, A. D¨ olle, F. M¨ uller-Plathe, Molecular dynamics of benzene in neat liquid and a solution containing polystyrene.

13

c nuclear mag-

netic relaxation and molecular dynamics simulation results, J. Phys. Chem. A 104 (2000) 5716–5725.

SC

[8] G. S. Jas, Y. Wang, S. W. Pauls, C. K. Johnson, K. Kuczera, Influence of temperature and viscosity on anthracene rotational diffusion in organic sol-

NU

vents: Molecular dynamics simulations and fluorescence anisotropy studies, J. Chem. Phys. 107 (1997) 8800–8812.

[9] G. S. Jas, E. J. Larson, C. K. Johnson, K. Kuczera, Microscopic details of

MA

rotational diffusion of perylene in organic solvents : Molecular dynamics simulation and experiment vs debye-stokes-einstein theory, J. Phys. Chem. A 104 (2000) 9841–9852.

D

[10] T.-M. Chang, L. X. Dang, On rotational dynamics of an nh + 4 ion in water,

PT E

J. Chem. Phys. 118 (2003) 8813–8820. [11] C. Millot, J.-C. Soetens, N. Ahmad, R. Adnan, Molecular simulation of unusual dynamical properties of quinoline in liquid phase, Euro. Phys. Lett.

CE

96 (2011) 43002.

[12] J.-C. Soetens, N. Ahmad, R. Adnan, C. Millot, Molecular dynamics sim-

AC

ulations of quinoline in the liquid phase, J. Phys. Chem. B 116 (2012) 5719–5728.

[13] N. Ahmad, R. Adnan, J.-C. Soetens, C. Millot, Molecular dynamics simulations of liquid isoquinoline as a function of temperature, Chem. Phys. 407 (2012) 29–38.

38

ACCEPTED MANUSCRIPT

[14] B. J. Berne, J. A. Montgomery, Jr, The coupling between translational and rotational motions, Mol. Phys. 32 (1976) 363–378. [15] A. R. Davies, G. J. Evans, M. W. Evans, Coupling of rotation and translation in molecular fluids, Faraday Disc. Chem. Soc. 66 (1978) 231–238.

RI PT

[16] A. Idrissi, P. Damay, S. Krishtal, M. Kiselev, Analysis of the effect of translation-rotation coupling on diffusion along the molecular axes, J. Phys. Chem. B 110 (2006) 18560–18565.

SC

[17] A. Idrissi, S. Longelin, P. Damay, S. Krishtal, M. Kiselev, Analysis of the transverse and the longitudinal pseudodiffusion of co 2 in sub- and supercritical states: A molecular-dynamics analysis, J. Chem. Phys. 125 (2006)

NU

224501.

[18] A. Faraone, L. Liu, S.-H. Chen, Model for the translation-rotation coupling

MA

of molecular motion in water, J. Chem. Phys. 119 (2003) 6302–6313. [19] D. van der Spoel, P. J. van Maaren, H. J. C. Berendsen, A systematic study of water models for molecular simulation: Derivation of water models

D

optimized for use with a reaction field, J. Chem. Phys. 108 (1998) 10220–

PT E

10230.

[20] P. Mark, L. Nilsson, Structure and dynamics of the tip3p, spc, and spc/e water models at 298 k, J. Phys. Chem. A 105 (2001) 9954–9960.

CE

[21] B. Guillot, A reappraisal of what we have learnt during three decades of computer simulations on water, J. Mol. Liq. 101 (2002) 219–260.

AC

[22] D. Paschek, Temperature dependence of the hydrophobic hydration and interaction of simple solutes: An examination of five popular water models, J. Chem. Phys. 120 (2004) 6674–6690.

[23] C. Vega, C. McBride, E. Sanz, J. L. F. Abascal, Radial distribution functions and densities for the spc/e, tip4p and tip5p models for liquid water

39

ACCEPTED MANUSCRIPT

and ices ih , ic , ii, iii, iv, v, vi, vii, viii, ix, xi and xii, Phys. Chem. Chem. Phys. 7 (2005) 1450–1456. [24] C. Vega, J. L. F. Abascal, E. Sanz, L. G. MacDowell, C. McBride, Can simple models describe the phase diagram of water ?, J. Phys.: Condens.

RI PT

Matter 17 (2005) S3283–S3288. [25] C. Vega, J. L. F. Abascal, Relation between the melting temperature and the temperature of maximum density for the most common models of water,

SC

J. Chem. Phys. 123 (2005) 144504.

[26] J. L. F. Abascal, C. Vega, A general purpose model for the condensed phases of water: Tip4p/2005, J. Chem. Phys. 123 (2005) 234505.

NU

[27] H. J. C. Berendsen, J. R. Grigera, T. P. Straatsma, The missing term in effective pair potentials, J. Phys. Chem. 91 (1987) 6269–6271.

MA

[28] M. G. Mazza, N. Giovambattista, H. E. Stanley, F. W. Starr, Connection of translational and rotational dynamical heterogeneities with the breakdown of the stokes-einstein and stokes-einstein-debye relations in water, Phys.

D

Rev. E 76 (2007) 031203.

PT E

[29] I. M. Svishchev, P. G. Kusalik, Roto-translational motion in liquid water and its structural implication, Chem. Phys. Lett. 215 (1993) 596–600. [30] M. Svanberg, An improved leap-frog rotational algorithm, Mol. Phys. 92

CE

(1997) 1085–1088.

[31] D. Fincham, Leapfrog rotational algorithms, Mol. Simul. 8 (1992) 165–178.

AC

[32] A. J. C. Ladd, Monte carlo simulation of water, Mol. Phys. 33 (1977) 1039– 1050.

[33] A. J. C. Ladd, Long-range dipolar interactions in computer simulations of polar liquids, Mol. Phys. 36 (1978) 463–474.

40

ACCEPTED MANUSCRIPT

[34] M. Neumann, Dielectric properties and the convergence of multipolar lattice sums, Mol. Phys. 60 (1987) 225–235. [35] C. Millot, J.-L. Rivail, R. Diguet, Static dielectric constant density and temperature dependence for the tips model of liquid methyl chloride, Chem.

RI PT

Phys. Lett. 160 (1989) 228–232. [36] F. Dehez, M. T. C. Martins Costa, D. Rinaldi, C. Millot, Long-range electrostatic interactions in hybrid quantum and molecular mechanical dynam-

SC

ics using a lattice summation approach, J. Chem. Phys. 122 (2005) 234503. [37] D. J. Adams, G. S. Dubey, Taming the ewald sum in the computer simulation of charged systems, J. Comput. Phys. 72 (1987) 156–176.

NU

[38] B. Cichocki, B. U. Felderhof, K. Hinsen, Electrostatic interactions in periodic coulomb and dipolar systems, Phys. Rev. A 39 (1989) 5350–5358.

MA

[39] C. Millot, J.-C. Soetens, M. T. C. Martins Costa, Static dielectric constant of the polarizable stockmayer fluid. comparison of the lattice summation and reaction field methods, Mol. Simul. 18 (1997) 367–383.

D

[40] M. P. Allen, D. J. Tildesley, Computer simulation of liquids, Oxford Uni-

PT E

versity Press, Oxford, 1987. [41] K. Bagchi, S. Balasubramanian, M. L. Klein, The effects of pressure on structural and dynamical properties of associated liquids: Molecular dy-

CE

namics calculations for the extended simple point charge model of water, J. Chem. Phys. 107 (1997) 8561–8567.

AC

[42] R. Kubo, in: U. of Colorado (Ed.), Lectures in Theoretical Physics, Vol. 1, Wiley (Interscience), New York, 1959.

[43] H. Mori, I. Oppenheim, J. Ross, in: J. de Boer, G. E. Uhlenbeck (Eds.), Studies in Statistical Mechanics, North-Holland Publishing Company, Amsterdam, 1962.

41

ACCEPTED MANUSCRIPT

[44] R. Zwanzig, Time-correlation functions and transport coefficients in statistical mechanics, Annu. Rev. Phys. Chem. 16 (1965) 67–102. [45] W. A. Steele, Molecular reorientation in liquids. i. distribution functions and friction constants, J. Chem. Phys. 38 (1963) 2404–2410.

RI PT

[46] W. A. Steele, Molecular reorientation in liquids. ii. angular autocorrelation functions, J. Chem. Phys. 38 (1963) 2411–2418.

[47] W. T. Huntress Jr, The study of anisotropic rotation of molecules in liquids

SC

by nmr quadrupolar relaxation, Adv. Magn. Reson. 4 (1970) 1–37. [48] M. W. Evans, Computer simulation of liquid anisotropy. v. nonlinear molec-

NU

ular dynamics at high field strengths, J. Chem. Phys. 78 (1983) 925–930. [49] K. T. Gillen, D. C. Douglass, M. J. R. Hoch, Self-diffusion in liquid water

MA

to −31o c, J. Chem. Phys. 57 (1972) 5117–5119. [50] R. Mills, Self-diffusion in normal and heavy water in the range 1−45o , J. Phys. Chem. 77 (1973) 685–688.

D

[51] K. Krynicki, C. D. Green, D. W. Sawyer, Pressure and temperature de-

PT E

pendence of self-diffusion in water, Faraday Discuss. Chem. Soc. 66 (1978) 199–208.

[52] K. R. Harris, L. A. Woolf, Pressure and temperature dependence of the self

CE

diffusion coefficient of water and oxygen-18 water, J. Chem. Soc. Faraday Trans. I 76 (1980) 377–385. [53] F. X. Prielmeier, E. W. Lang, R. J. Speedy, H.-D. L¨ udemann, The pressure

AC

dependence of self diffusion in supercooled light and heavy water, Ber. Bunsenges. Phys. Chem. 92 (1988) 1111–1117.

[54] A. J. Easteal, W. E. Price, L. A. Woolf, Diaphragm cell for hightemperature diffusion measurements. tracer diffusion coefficients for water to 363 k, J. Chem. Soc. Faraday Trans. 1 85 (1989) 1091–1097.

42

ACCEPTED MANUSCRIPT

[55] M. Holz, S. R. Heil, A. Sacco, Temperature-dependent self-diffusion coefficients of water and six selected molecular liquids for calibration in accurate 1

h nmr pfg measurements, Phys. Chem. Chem. Phys. 2 (2000) 4740–4742.

[56] B. D¨ unweg, K. Kremer, Molecular dynamics simulation of a polymer chain

RI PT

in solution, J. Chem. Phys. 99 (1993) 6983–6997. [57] I.-C. Yeh, G. Hummer, System-size dependence of diffusion coefficients and viscosities from molecular dynamics simulations with periodic boundary

SC

conditions, J. Phys. Chem. B 108 (2004) 15873–15879.

[58] N. Galamba, On the hydrogen-bond network and the non-arrhenius transport properties of water, J. Phys.: Condens. Matter 29 (2017) 015101.

NU

[59] K. Watanabe, M. L. Klein, Effective pair potentials and the properties of water, Chem. Phys. 131 (1989) 157–167.

MA

[60] I. I. Vaisman, L. Perera, M. L. Berkowitz, Mobility of stretched water, J. Chem. Phys. 98 (1993) 9859–9862.

[61] A. Chandra, T. Ichiye, Dynamical properties of the soft sticky dipole model

PT E

2709.

D

of water: Molecular dynamics simulations, J. Chem. Phys. 111 (1999) 2701–

[62] M. W. Mahoney, W. L. Jorgensen, Diffusion constant of the tip5p model

CE

of liquid water, J. Chem. Phys. 114 (2001) 363–366. [63] Y. Wu, H. L. Tepper, G. A. Voth, Flexible simple point-charge water model with improved liquid-state properties, J. Chem. Phys. 124 (2006) 024503.

AC

[64] A. Y. Zasetsky, S. V. Petelina, A. K. Lyashchenko, A. S. Lileev, Computer simulation study of rotational diffusion in polar liquids of different types, J. Chem. Phys. 133 (2010) 134502.

[65] O. Gereben, L. Pusztai, System size and trajectory length dependence of the static structure factor and the diffusion coefficient as calculated from

43

ACCEPTED MANUSCRIPT

molecular dynamics simulations: The case of spc/e water, J. Mol. Liq. 161 (2011) 36–40. [66] G. Guevara-Carrion, J. Vrabec, H. Hasse, Prediction of self-diffusion coefficient and shear viscosity of water and its binary mixtures with methanol

RI PT

and ethanol by molecular simulation, J. Chem. Phys. 134 (2011) 074508. [67] U. Dahal, N. P. Adhikari, Molecular dynamics study of diffusion of heavy water in normal water at different temperatures, J. Mol. Liq. 167 (2012)

SC

34–39.

[68] S. Tazi, A. Botan, M. Salanne, V. Marry, P. Turq, B. Rotenberg, Diffusion

Matter 24 (2012) 284117.

NU

coefficient and shear viscosity of rigid water models, J. Phys.: Condens.

[69] S. H. Lee, Temperature dependence on structure and self-diffusion of water:

MA

a molecular dynamics simulation study using spc/e model, Bull. Korean Chem. Soc. 34 (2013) 3800–3804.

[70] E. Galicia-Andr´es, H. Dominguez, L. Pusztai, O. Pizio, On the composition

D

dependence of thermodynamic, dynamic and dielectric properties of water-

PT E

methanol model mixtures. molecular dynamics simulation results, Condens. Matter Phys. 18 (2015) 43602. [71] R. J. Speedy, C. A. Angell, Isothermal compressibility of supercooled water

CE

and evidence for a thermodynamic singularity at −45o c, J. Chem. Phys. 65 (1976) 851–858.

AC

[72] V. V. Kisel’nik, N. G. Malyuk, A. I. Toryanik, V. M. Toryanik, Effect of pressure and temperature on the self-diffusion of water, Zh. Strukt. Khim. 14 (1973) 963–967.

[73] L. A. Woolf, Tracer diffusion of tritiated water (tho) in ordinary water (h2 o) under pressure, J. Chem. Soc. Faraday Trans. I 71 (1975) 784–796.

44

ACCEPTED MANUSCRIPT

[74] C. A. Angell, E. D. Finch, L. A. Woolf, P. Bach, Spin-echo diffusion coefficients of water to 2380 bar and −20o c, J. Chem. Phys. 65 (1976) 3063–3066. [75] F. X. Prielmeier, E. W. Lang, R. J. Speedy, H.-D. L¨ udemann, Diffusion in supercooled water to 300 mpa, Phys. Rev. Lett. 59 (1987) 1128–1131.

RI PT

[76] K. R. Harris, P. J. Newitt, Self-diffusion of water at low temperatures and high pressure, J. Chem. Eng. Data 42 (1997) 346–348.

[77] J. Qvist, C. Mattea, E. P. Sunde, B. Halle, Rotational dynamics in supercooled water from nuclear spin relaxation and molecular simulations, J.

SC

Chem. Phys. 136 (2012) 204505.

[78] M. G. Mazza, N. Giovambattista, F. W. Starr, H. E. Stanley, Relation be-

Rev. Lett. 96 (2006) 057803.

NU

tween rotational and translational dynamic heterogeneities in water, Phys.

MA

[79] D. Laage, J. T. Hynes, On the molecular mechanism of water reorientation, J. Phys. Chem. B 112 (2008) 14230–14242. [80] I. M. Svishchev, P. G. Kusalik, Structure in liquid water: A study of spatial

D

distribution functions, J. Chem. Phys. 99 (1993) 3049–3058.

PT E

[81] C. Vega, J. F. L. Abascal, Simulating water with rigid non-polarizable models: A general perspective, Phys. Chem. Chem. Phys. 13 (2011) 19663– 19688.

CE

[82] D. Braun, S. Boresch, O. Steinhauser, Transport and dielectric properties of water and the influence of coarse-graining: Comparing bmw, spc/e and

AC

tip3p models, J. Chem. Phys. 140 (2014) 064107. [83] S. Chowdhuri, M.-L. Tan, T. Ichiye, Dynamical properties of the soft sticky dipole-quadrupole-octupole water model: A molecular dynamics study, J. Chem. Phys. 125 (2006) 144513.

[84] B. Halle, H. Wennerstr¨ om, Interpretation of magnetic resonance data from water nuclei in heterogeneous systems, J. Chem. Phys. 75 (1981) 1928–1943. 45

ACCEPTED MANUSCRIPT

[85] M. S. P. Sansom, I. D. Kerr, J. Breed, R. Sankararamakrishnan, Water in channel-like cavities: Structure and dynamics, Biophys. J. 70 (1996) 693–702. [86] K. Krynicki, Proton spin-lattice relaxation in pure water between 0 o c and

RI PT

100 o c, Physica 32 (1966) 167–178. [87] J. R. C. van der Maarel, D. Lankhorst, J. de Bleijser, J. C. Leyte, On the single-molecule dynamics of water from proton, deuterium and oxygen-17

SC

nuclear magnetic relaxation, Chem. Phys. Lett. 122 (1985) 541–544. [88] R. P. W. J. Struis, J. de Bleijser, J. C. Leyte, Dynamic behaviour and some of the molecular properties of water molecules in pure water and in mgcl 2

NU

solutions, J. Phys. Chem. 91 (1987) 1639–1645.

[89] R. Ludwig, Nmr relaxation studies in water-alcohol mixtures: The water-

MA

rich region, Chem. Phys. 195 (1995) 329–337. [90] J. Ropp, C. Lawrence, L. T. C. Farrar, J. L. Skinner, Rotational motion in liquid water is anisotropic: A nuclear magnetic resonance and molecular

D

dynamics simulation study, J. Am. Chem. Soc. 123 (2001) 8047–8052.

PT E

[91] E. W. Lang, H.-D. L¨ udemann, High pressure o-17 longitudinal relaxation time studies in supercooled h 2 o and d2 o, Ber. Bunsenges. Phys. Chem. 85 (1981) 603–611.

CE

[92] E. Lang, H.-D. L¨ udemann, Pressure and temperature dependence of the longitudinal proton relaxation times in supercooled water to −87o c and

AC

2500 bar, J. Chem. Phys. 67 (1977) 718–723. [93] J. Jonas, T. DeFries, D. J. Wilbur, Molecular motions in compressed liquid water, J. Chem. Phys. 65 (1976) 582–588.

[94] T. DeFries, J. Jonas, Pressure dependence of nmr proton spin-lattice relaxation times and shear viscosity in liquid water in the temperature range −15-10o c, J. Chem. Phys. 66 (1977) 896–901. 46

ACCEPTED MANUSCRIPT

[95] C. Wakai, M. Nakahara, Pressure dependencies of rotational, translational, and viscous friction coefficients in water-d 2 , acetonitrile-d3 , acetonitrile, chloroform, and benzene, J. Chem. Phys. 100 (1994) 8347–8458. [96] L. D. Favro, Theory of the rotational brownian motion of a free rigid body,

RI PT

Phys. Rev. 119 (1960) 53–62. [97] W. T. Huntress Jr, Effects of anisotropic molecular rotational diffusion on

AC

CE

PT E

D

MA

NU

SC

nuclear magnetic relaxation in liquids, J. Chem. Phys. 48 (1968) 3524–3533.

47

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN

US

CR

IP

T

Graphical abstract

ACCEPTED MANUSCRIPT Highlights

CE

PT

ED

M

AN

US

CR

IP

T

SPC/E water model predicts increase of diffusion with pressure at low temperature Translation-rotation coupling in liquid water decreases as temperature decreases Translation-rotation coupling diminishes only slightly with pressure in kbar range Quantitative limitation of rotational diffusion Debye model is assessed

AC

   

Figure 1

Figure 2

Figure 3

Figure 4

Figure 5

Figure 6

Figure 7

Figure 8

Figure 9

Figure 10

Figure 11