Wavevector- and frequency-dependent shear viscosity of water: the modified collective mode approach and molecular dynamics calculations

The transverse momentum time autocorrelation functions and wavevector-and frequency-dependent shear viscosity are calculated for an interaction site model of water using a modified collective mode approach and molecular dynamics simulations. The modified mode approach is based on a formulation which consistently takes into account non-Markovian effects into the kinetic memory kernels. As is demonstrated by comparing the theory results with the molecular dynamics data, the entire frequency dependence of the shear viscosity can be reproduced quantitatively over the whole wavelength range in terms of six generalized collective modes employing the kinetic memory kernel in the non-Markovian approximation of the third order. It is also shown that the results corresponding to the exact atomic and abbreviated molecular descriptions may differ considerably. In the infinite wavevector regime the dynamic correlations are completely determined by a single free motion of the molecules. It is a great pleasure and a great honour for us to dedicate this paper to our good colleague and friend Prof. Reinhard Folk who celebrates this year his 60th birthday.

Despite the progress in qualitative understanding the origin of collective excitations in water and other molecular liquids, very little is done in developing a consistent theory enabling us to present quantitatively the time correlation functions and generalized transport coefficients, such as the wavevector-and frequency-dependent shear and bulk viscosities as well as the thermal conductivity.Previous attempts in this direction have been made by Bertolini and Tani [10,11,14] who applied a generalized hydrodynamics description to the TIP4P and SPC/E models of water.In this description, however, the non-Markovian contributions to the kinetic memory functions are modelled in a semiphenomenological way by introducing adjustable parameters.These parameters are determined by fitting the theory results to MD data.Similar procedures have been used when investigating dielectric properties [4] and longitudinal momentum fluctuations [21] in the TIP4P and SPC/E water.
Recently [19], a modified formulation of the generalized collective mode (GCM) approach [23][24][25] has been proposed to describe the processes of dynamical polarization in interaction site models of molecular liquids.The modified formulation consistently takes into account non-Markovian terms in the kinetic memory kernels by involving additional dynamical quasivariables leading to the correlation times of higher orders.This has allowed us to accurately reproduce the entire frequency and wavevector dependence of dielectric functions for the TIP4P water without applying any fitting procedure.The spectrum of time correlation functions corresponding to the polarization fluctuations has been associated with diffusive Debye-like and propagating librational modes.
In [19] it has been concluded that for the case of molecular liquids, the standard formulation [23][24][25] of the GCM approach (which is based on the Markovian route for kinetic memory kernels) does not necessarily provide a rapidly converging approximation, even when large numbers of high-order dynamical variables are included into the kinetic set.This is contrary to the case of simple Lennard-Jones fluids, where the Markovian approximation has proved to be a reasonable route for a quantitative description of transport coefficients [26].This approximation has also been successfully used to study extended hydrodynamic, dielectric and magnetic properties in simple liquid mixtures [27,28], dipolar systems [29,30], and spin fluids [31][32][33].Due to a very complicated interplay between the hydrodynamic variables and those associated with internal degrees of freedom, the description of molecular liquids requires a more sophisticated approximation.
In the present study we develop the modified GCM approach to describe the processes concerning the propagation and dissipation of transverse momentum fluctuations in molecular liquids.The paper is organized as follows.An outline of the theory is given in section 2. The MD simulations are briefly described in section 3. The obtained theoretical results for the transverse momentum and its current time autocorrelation functions as well as the generalized shear viscosity of the TIP4P water are presented and compared with the MD data in section 4. The concluding remarks are summarized in section 5.

Theory
Let us consider a molecular liquid consisting of N rigid molecules each of which is composed of M interaction sites.Applying the spatial Fourier transform, one obtains the following expression for the microscopic operator of the total momentum of the system in the atomic description where r ia (t) and v ia (t) = (d/dt)r ia (t) are correspondingly the position and velocity of site a carrying mass m a within the ith molecule at time t.The extended hydrodynamic variable (1) can be split as J = J L + J T into the longitudinal J L = k J• k and transverse J T = k×[J× k] components with respect to the wavevector k, where k = k/k.The component J T is orthogonal to J L and other hydrodynamic variables such as the molecule number density n(k, t) = N,M i,a e −ik•r ia (t) and the total energy e(k, t) = N,M i,a e ia (t) e −ik•r ia (t) , where e ia (t with ϕ ab ij being the site-site intermolecular potentials.The transverse fluctuations can thus be studied separately from the density, longitudinal momentum, and energy correlations. The collective variable (1) satisfies the equation of motion where L is the Liouville operator.For rigid models it is convenient to express the atomic phase coordinates in terms of translational and rotational motions of the molecules as where δ ia is the position of site a belonging to molecule i relatively to its center of mass r i moving with the translational velocity v i , while w i denotes the molecule angular velocity.Then the Liouville operator takes the form where denotes the matrix of moments of inertia, m = a m a is the molecule mass, and 1 the unit matrix.The formal solution of equation ( 2) is where J T (k) ≡ J T (k, t = 0).The transverse momentum time autocorrelation function can thus be introduced in the form where denotes the equilibrium averaging.In view of (5), this function can be approximately presented at short times using the truncated Taylor series where f 2s (k) ≡ f 2s (k, t = 0) are nonzero static correlation functions (SCFs) with s = 0, 1, 2, . . ., S (f 0 (k) ≡ f (k)) and being the higher-order (HO) (s 1) time correlation functions (TCFs).The SCFs f 2s (k) are defined on the set {J T (k), LJ T (k), . . ., LS J T (k)} with Ls J T (k) being the HO time derivatives of J T (k, t) at t = 0.The main idea of the modified GCM approach [19] is based on the introduction of the extended set L−S J T (k), . . ., L−1 J T (k), J T (k), LJ T (k), . . ., LS J T (k) (8) which includes apart from the HO array Ls J T (k) the new lower-order (LO) components L−s J T (k) of the basic dynamical variable J T (k), where s = 1, 2, . . ., S .The LO components are evaluated in view of the rule L L−1 = 1 with 1 being the unit operator.Then by definition L−s where C s are some constants.Unlike the HO variables which can be written explicitly acting by the Liouville operator, the LO quantities require the direct integration of motion and cannot be presented in an explicit form.For this reason, these quantities should be treated as quasivariables.
Using set (8), one constructs the matrix of TCFs, where α, β = −S , −S +1, Since the usual (α, β 0) TCFs vanish at t → ∞, the constants C s can be determined by fulfilling the condition lim t→∞ f αβ (k, t) = 0 for all the rest elements of F(k, t).Then integrating by parts, the LO (α+β < 0) SCFs can be expressed uniquely in terms of the basic function It can be shown within the memory function formalism [19,24,34] that the TCFs satisfy the generalized Langevin equation where Ω(k) and Γ(k, ω) denote the frequency and memory matrices, respectively, is the Laplace transform.Assuming that at given values of S and S the matrix F(k, t) includes an almost complete set of slow TCFs, we can apply the Markovian approximation Γ(k, ω) ≈ Γ(k, ω = 0).Then equations (10) simplifies, [1iω + T(k)]F (S S) (k, ω) = F(k), and can be solved analytically.The result in time representation is f , where z γ is an eigenvalue corresponding to an eigenvector γl are the elements of matrices T and x −1 , respectively, and F (S S) denotes the TCFs in the Markovian approximation.Thus, each element of F (S S) (k, t) can be expressed as a weighted sum of S + S + 1 exponentials which are connected with the generalized collective modes z γ (k).More explicit expressions for F(k) and F(k, ω = 0) ≡ F(k) are (S = 3 and S = 3): ).It is necessary to underline that the Markovian approximation applied to the extended set (8) of dynamical variables does not concern the usual memory kernels Γ S (k, ω) corresponding to the HO array Ls J T (k).Indeed, solving the system of linear equations (10) with respect to f (k, ω), the result can be cast in the form of the S-th order continued fraction where ] and so on [35] are HO elements of the frequency matrix, and is the representation of the exact kinetic memory kernel Γ S (k, ω) in the S -th approximation [19].We see that at S = 0 the function Γ S S depends on frequency in a characteristic way and, therefore, may describe essential non-Markovian effects.For S = 0, the modified mode approach reduces to its standard formulation [23][24][25].
The coefficients a (s) s and b s appearing in the Pade-like approximation ( 12) are directly connected with the SCFs and correlation times τ ν (k).In particular, for S = 0 and S = 0, 1, . . .one finds that a For any S > 0 they are calculated using the recurrent relations a 0 ) 2 and so on for S = 2, 3, . .., where n = 0, 1, . . ., S − 1, and From the above it becomes clear that first 2S-fold derivatives at t = 0 and 2S time integral moments for the genuine f (k, t) and approximated f (S S) (k, t) functions coincide.In other words, the sum rules are fulfilled exactly up to order 2S in time space and to order 2S in the frequency domain.In such a way, the time (or frequency) dependence of f (k, t) (or f (k, ω)) is reproduced for arbitrary t (or ω).It is also worth emphasizing that the non-Markovian terms cannot be handled correctly by a simple expansion of the exact memory kernel Γ S (k, ω) at ω → 0 into the truncated Taylor series Γ S S (k, ω) = Γ S + ω∂Γ S /∂ω Despite the fact that both Γ S S and Γ S S lead to the exact sum rules up to order 2S in the frequency space, the function Γ S S distorts the true short-time behaviour of f (k, t), because the condition lim ω→∞ Γ S (k, ω)/ω = 0 = lim ω→∞ Γ S S /ω is broken for Γ S S .Only the Pade-like presentation Γ S S of Γ S is capable of satisfying the sum rules simultaneously in the time and frequency domains.
The Langevin equation (10) being rewritten directly for f (k, ω) yields the equali- where Γ(k, ω) (≡ Γ 1 (k, ω)) is the hydrodynamic memory function.The wavevectorand frequency-dependent shear viscosity η(k, ω) is simply connected with Γ(k, ω) by the relation [11,24,26,36] where ρ = N V m ≡ nm is the molecular mass density with V being the volume of the system.In view of equations ( 11), (13), and ( 14) one obtains the following analytical expression for the shear viscosity in the modified mode approximation Equivalently, expression (15) can also be presented in terms of generalized modes as

MD simulations
The shear viscosity of water was calculated by means of MD simulations in a lot of previous works [9,11,13,17,37,38].The most of them were devoted to the investigation of the usual viscosity η corresponding to η(k, ω) in the limits k → 0 and ω → 0. Only a few works dealt with the calculations of the generalized viscosity η(k, ω) with k = 0 or ω = 0 or both [9,11,17].However, the computations were restricted to rather narrow regions in the wavevector k and frequency ω domains.Note that there are several schemes enabling us to evaluate η at k = 0 and ω = 0.Among them it is worth pointing out the Green-Kubo method [39][40][41], the Einstein scheme [42], the approach by Hess and Evans [38,43], as well as the nonequilibrium [44] and reverse nonequilibrium [45] molecular dynamics.For finite k = 0 and ω = 0 (where the conventional methods are unapplicable), the generalized shear viscosity η(k, ω) is calculated in terms of the transverse momentum time autocorrelation function f (k, t) (see equation ( 6)) using relation (14).
In the case k → 0 and ω → 0, relation (14) reduces to the well-known Green-Kubo formula [39][40][41] where denotes the transverse component of the stress vector, T and k B are the temperature and Boltzmann's constant, respectively, and f ia (t) = m a ( Lv i + Lw i ×δ ia ) is the total force (intermolecular plus constraint) acting on atom a of molecule i.Indeed, explicitly taking the derivative of J T (k, t) with respect to time one finds so that the dynamical variable JT (k, t) at k → 0 tends to zero (because the total momentum is a conserved quantity) like lim k→0 JT (k, t) = −ik σ T (t), where the third Newton law i,a f ia = 0 has been applied.Therefore, the transverse current momentum autocorrelation function will behave at k → 0 as lim k→0 f 2 (k, t) = k 2 f σ (t).On the other hand, in view of the equality , so that taking into account equation ( 14) leads to the frequency-dependent shear viscosity where the relation lim k→0 f (k) = mk B T (see Appendix) has been used.For ω = 0, the right-hand side of equation ( 20) just reproduces the Green-Kubo formula (16).
For k = 0 and ω = 0 we obtain from equation ( 14) the wavevector-dependent shear viscosity Our MD simulations were carried out for the TIP4P (transferable intermolecular potential with four points, M = 4) model [46] of water in the microcanonical ensemble at a density of ρ = 1 g/cm 3 and at a temperature of T = 292 K.We considered N = 256 molecules in the cubic volume to which the toroidal boundary conditions have been applied.The interaction site reaction field geometry [47] has been utilized to reduce the finite-size effects.The equations of motion were integrated based on the matrix method [48] with a time step of 2 fs.The time of observation over the equilibrium state was 1 000 000 steps.The TCFs functions were evaluated for the wide wavenumber set k = [0, 1, . . ., 300, . ..] k min up to k max = 100 Å−1 , where k min = 2π/ 3 √ V = 0.319 Å−1 and 3 √ V is the length of the simulation box edge.The wavevector-and frequency-dependent shear viscosity η(k, ω) was calculated in view of equation ( 14) by applying the Laplace transform f (k, ω) = L iω f (k, t) to the transverse momentum time autocorrelation function (6) up to frequencies of ω max = 10 000 THz.The error of the MD calculations was estimated to be of order 0.1 ÷ 0.5%.
When generalized hydrodynamics is considered for molecular liquids, the atomic or a molecular description can be used.In the exact atomic description each atom yields its own contribution to the dynamical variables, modulated by its own phase factor.For instance, the above definitions for the dynamical variables J(k, t) (equation ( 1)) and J(k, t) (equation ( 18)) have been written in the atomic phase coordinates.In the approximate molecular description a single phase factor is employed for each molecule using the coordinate of its center of mass.In particular, the molecular counterparts of equations ( 1) and (18) are In the limit k → 0, the atomic and molecular descriptions lead to the same results for f (k, t) and η(k, w), because then lim k→0 J T (k, t) = lim k→0 J(k, t) (due to the property a m a δ ia (t) = 0 which follows from the definition of the center of mass).This is, however, not the case for finite wavenumbers k = 0 (and higher-order dynamical variables even at k → 0), where the details of the distribution of atoms within the molecule become important, and the k-dependence of the results can be significantly different.In our MD calculations we have used these both descriptions in order to explicitly demonstrate this difference.

Results and discussion
The reduced transverse momentum static autocorrelation function C(k) = f (k)/k B T and shear rigidity modulus G(k) = ρ/k 2 • f 2 (k)/f (k) of the TIP4P water are shown in subsets (a) and (b) of figure 1, respectively.Note that the function f (k) can be calculated exactly in both the atomic and molecular descriptions using analytical expressions (see Appendix).We also present the corresponding MD data in order to estimate the error of the MD calculations.As can be seen the deviations between the MD data and the exact results are very small (of order 0.5% or even less within the atomic description).On the other hand, the differences between the atomic and molecular presentations of C(k) and G(k) are significant and constitute about 10 ÷ 20%.These differences increase with the increasing k but they vanish for C(k) at k → 0. However, the rigidity modulus G(k) obtained in the exact atomic description considerably differs from that corresponding to the abbreviated molecular description even in the hydrodynamics limit k → 0. Thus the latter description cannot be used for interaction site models when a high accuracy of the calculations is required.
The wavevector-dependent shear viscosity η(k) is plotted in figure 2. According to definition (21), the function η(k) is inverse proportional to the basic correlation time τ 0 (k).In the hydrodynamic regime k → 0, the correlation time increases with the increasing k like τ 0 (k) ∼ 1/(ηk 2 ), where η = lim k→0 η(k).In our calculations we have obtained η = 6.02•10 −4 Pa s, while the experimental value for the shear viscosity of real water at normal conditions is about 8.9 •10 −4 Pa s [49].The TIP4P model of water thus underestimates η.In the infinite wavenumber regime k → ∞, the shear viscosity tends to zero.It is interesting to remark that the atomic and molecular descriptions lead almost to the same results for η(k) over the whole wavevector range with only slight deviations at intermediate and large values of k.Samples of the normalized, transverse momentum Φ(k, t) = f (k, t)/f (k) and current momentum Ψ(k, t) = f 2 (k, t)/f 2 (k) time autocorrelation functions are shown in figures 3 and 4, respectively, for various wavevector values.Within the modified GCM approach we have used the non-Markovian approximation at S = 2 with S = 3 that corresponds to six (S + S + 1) modes.Further increase of S and S practically does not affect the results.It is worth emphasizing that within the usual Markovian approximation when S = 0 we could not obtain a satisfactory description at any reasonable values of S. The convergence of the theoretical results with the increasing S was very slow at S = 0 even for extremely large values of S. The characteristic intervals of varying the functions Φ(k, t) and Ψ(k, t) at short and long times differs significantly (up to several orders) due to the presence of fast rotational and relatively slow translational motions of water molecules.This leads to an essential non-Markovian behaviour which cannot be handled by the standard approximation.On the other hand, choosing S = 3 allows us to speed up the convergence drastically, so that an excellent agreement between the theory and MD approach has been achieved already at S = 2.This can be seen in figure 3 and 4, where the theory and MD results for Φ(k, t) and Ψ(k, t) are almost undistinguishable.The functions Φ(k, t) and Ψ(k, t) obtained within the atomic and molecular descriptions differ slightly at small wavevectors.However, the difference increases considerably with rising k especially in the case of function Ψ(k, t) and can be observed at small and intermediate times t.With the increasing t all the functions tend to zero, where a long-time tail behaviour can be indicated at small and intermediate values of k.
In the large wavevector regime, the functions Φ(k, t) and Ψ(k, t) approach their free-motion counterparts.Such counterparts have been calculated by us in separate stochastic simulations.In these simulations the equations of translational and rotational motion were integrated in the absence of any intermolecular interactions, while the initial translational and angular velocities of the molecules were generated according to the Maxwell distribution corresponding to an equilibrium state with given temperature T .
We mention that the static wavevector-dependent correlation functions f 2s (k) (s = 1, . . ., S) as well as the wavevector-dependent correlation times τ ν (k) (ν = 0, 1, . . ., 2S ) should be considered as input quantities to the theory.In other words this theory predicts dynamical behaviour of the system provided the static properties are already known.Such quantities might be obtained within a microscopic theory.Note that the function f (k) is calculated analytically (see Appendix).The higher-order quantities f 2 (k), f 4 (k), and so on are expressed in terms of site-site distribution functions.The latter functions can be evaluated within the RISM integral equation theory [50,51]), where the hypernetted chain or other approximations is used for the closure relation.The high-order static correlation functions can also be expressed in terms of lower-order ones using various decoupling approximations.The basic correlation time τ 0 (k) is directly related to the static shear viscosity η(k) (see equation ( 21)).Higher-order (ν 1) correlation times τ ν (k) define the non-Markovian behaviour of kinetic memory kernel (12).With the increasing ν the chief contribution to τ ν (k) will be defined mainly by the shape of f (k, t) at long times, which in turn is formed by diffusion processes.Therefore, the quantities τ ν (k) can be calculated by solving the corresponding diffusion equations (like those presented in [52]).Alternatively, the correlation times τ ν (k) can be expressed in terms of f 2s (k) using the relations a 1 and so on for sufficiently large S, where the convergence in a s s and b s s is expected.These relations constitute a system of nonlinear equations for τ ν (k) which should be solved numerically.The above problems are beyond the scope of the present paper and will be considered elsewhere.Here, in order to avoid additional uncertainties, the required quantities f 2s (k) and τ ν (k) were evaluated exactly in the MD simulations.
To visualize the collective mode spectrum, the Laplace transforms Φ(k, ω) and Ψ(k, ω) of Φ(k, t) and Ψ(k, t) are plotted in figure 6.On the existing local maxima in functions Φ(k, ω) Ψ(k, ω) at nonzero frequencies ω = 0 and on the form of the shape of these functions at ω = 0 and near the maxima, we can clearly identify at least  The real part of the wavevector-and frequency-dependent shear viscosity of the TIP4P water obtained within the modified GCM approach (solid curves) versus the MD data (circles) related to the atomic description.The MD results corresponding to the molecular description (at k = k min , dashed curve) as well as to frequency-dependent (k = 0, atomic description) shear viscosity (dots) are included in subset (a).three primary modes.They should be related to the extended diffusive hydrodynamic mode z h (k) = σ h (k) and two propagating kinetic optical-like excitations z O1,2 (k) = σ O1,2 (k) + iω O1,2 (k).Other auxiliary modes are not so visible on the scale of figure 3 due to the smallness of their weights.In the hydrodynamic limit k → 0, the extended hydrodynamic mode can be reproduced even by the one-mode approximation (S = S = 0), where Γ(k, ω) ≈ Γ(k) = k 2 /ρ • η (see equation ( 14)).Then according to equation (13), one obtains Φ(k, ω) ≈ 1/(iω + k 2 η/ρ) that corresponds in time space to Φ(k, t) ≈ exp(−ηk 2 t/ρ) = exp(−z h t).So that the extended hydrodynamic mode behaves at small k like σ h = ηk 2 /ρ and appears to be purely dissipative.The above single-exponential decay of Φ(k, t) is well-known in the theory of simple liquids [39,40].In the limits k → 0 and t → ∞ the exponential presentation exp(−ηk 2 t/ρ) becomes exact.This can be used for an alternative calculation of the shear viscosity in MD simulations by fitting Φ(k, t) to an exponential decay at long times and calculating η from the decay constant [9].
With the rising k, the role of the extended hydrodynamic mode decreases and the contributions of kinetic optic-like excitations to Φ(k, t) and Ψ(k, t) cannot be neglected.The lower z o1 (k) = σ O1 (k) + iω O1 (k) and higher z o2 (k) = σ O2 (k) + iω O2 (k) frequency lying, ω O1 (k) < ω O2 (k), optical-like modes can clearly be observed in the momentum current correlation functions (see figure 5).The propagating frequencies ω O1 (k) and ω O2 (k) can be approximately related to the positions of two peaks in Ψ(k, ω), while the sharpness of these peaks can be related to the dissipation coefficients σ O1 (k) and σ O2 (k).The lower-and higher-frequency optical modes should be associated with roll and pitch (in the terminology of [21,53]) librational motions, respectively.The pitch librational oscillations are caused by the fast rotations of molecules (due to the large torques and low moment of inertia of the water molecule) in the averaged field created by their neighbours due to the interactions.The pitch librational frequency takes the values ω O2 (k) ∼ 50 ÷ 200 THz in depending on the wavevector.As can be seen, the higher-frequency optical mode vanishes at large wavenumbers k > 20k min , where the free-motion regime begins to dominate and the contribution from the molecular interactions can be neglected.On the other hand, the lower-frequency optical mode, with ω O1 (k) ∼ 10 ÷ 200 for k ∈ [k min , 100k min ], does not disappear even in the large wavenumber regime, reflecting a single-molecule character of roll librational excitations.
Contrary to the extended hydrodynamic mode, where σ h → 0 with k → 0, the propagating frequencies of kinetic optical-like excitations remain finite even in the hydrodynamic limit k → 0. The extended hydrodynamic mode z h = ηk 2 /ρ is also inherent in the transverse momentum spectrum of simple Lennard-Jones fluids [24,26].Thus, in the case of molecular liquids, such a mode originates from the center-of-mass motions of constituent molecules.The optical-like excitations arise as the result of the rotational motion of molecules with finite size mass distribution.The latter excitations are thus absent in simple fluids.
The real η (k, ω) and imaginary η (k, ω) parts of the wavevector-and frequencydependent shear viscosity η(k, ω) = η (k, ω) − iη (k, ω) of the TIP4P water obtained in the MD simulations within the atomic description are plotted in figures 6 and 7, respectively, in comparison with the results of the modified GCM approach.The MD data related to frequency-dependent (k = 0) shear viscosity as well as to the molecular description are included too (subset (a)).Both the descriptions lead nearly to the same values of η(k, ω) except some difference in the libration frequency region ω ≈ 200 THz.It can be also seen that despite a rather complicated behaviour of functions η (k, ω) and η (k, ω) with ω, especially at small and intermediate k, the theory presented is capable of describing the generalized shear viscosity not only qualitatively but also quantitatively in the whole wavevector and frequency range.

Conclusion
In the present study we have developed the collective mode approach in order to quantitatively describe the processes concerning the propagation and dissipation of transverse momentum and current momentum fluctuations in molecular liquids.As is now well established and has been confirmed in our MD simulations for the TIP4P water, such processes exhibit a highly non-Markovian behaviour and thus cannot be handled within the standard approaches where the Markovian convention is assumed.The theory developed has taken into account the non-Markovian effects by consistently including new terms into the kinetic memory kernels.This has allowed us to precisely reproduce the wavevector-and frequency-dependent shear viscosity of the TIP4P water using a relatively small number of modes.The main advantage of our approach with respect to other theories is that it contains no adjustable parameters and thus no semiphenomenological fitting procedures are required.It has been also pointed out that the molecular description cannot be recommended for accurate calculations of the shear viscosity in MD simulations, and the preference should be given to the exact atomic description.
The approach proposed can also be adapted to the investigations of fluctuations related to the number density, longitudinal momentum and total energy.These and related problems, as well as the calculation of the wavevector-and frequencydependent bulk viscosity and thermal conductivity will be considered in further studies.
Note that equilibrium distribution functions are factorized into the coordinate and velocity parts.In its turn, the translational and angular velocities are distributed independently of one another for each molecule.As a result, the nonzero contributions to equilibrium averaging in (A2) give only the terms with coincident molecular indexes (i = j) of summation.It is more convenient to consider each molecule in its own principal coordinate system in which δ ia ≡ ∆ a .Then taking into account equation (3), expression (A2) transforms into where j 0 (z) = sin(z)/z, j 1 (z), and j 2 (z) = 3j 1 (z)/z − j 0 (z) are the spherical Bessel functions of order zero, one, and two, respectively, ρ = ρ/ρ and averaging in (A4) is performed over orientations of k-vector with respect to the molecule.Owing to the identity of molecules and isotropy of the system we have that v 2

Figure 1 .
Figure 1.The reduced transverse momentum static autocorrelation function (subset (a)) and shear rigidity modulus (subset (b)) of the TIP4P water.The MD results in the atomic and molecular descriptions are shown as full and open circles, respectively.The analytical calculations are presented in subset (a) by the solid (atomic) and dashed (molecular) curves.

Figure 2 .
Figure 2. The wavevector-dependent shear viscosity obtained in the MD simulations for the TIP4P water using the atomic (full circles connected by the solid curve) and molecular (open circles connected by the dashed curve) descriptions.

Figure 3 .
Figure3.The normalized, transverse momentum-momentum time correlation functions of the TIP4P water obtained within the modified GCM approach (solid curves) versus the MD data related to the atomic (cirlces) and molecular (dashed curves) descriptions.The result corresponding to the free-motion regime (in the atomic description) is also included (crosses, subsets (g) and (h)).

Figure 4 .
Figure 4.The normalized, transverse current momentum time autocorrelation functions of the TIP4P water.The stress autocorrelation function (in the atomic description) corresponding to the low wavenumber limit is also plotted (squares, subset (a)).Other notations are the same as for figure 3.

Figure 5 .
Figure 5.The Laplace transforms of the normalized transverse momentum (solid curves) and current momentum (dashed curves) time autocorrelation functions of the TIP4P water obtained within the modified GCM approach.The corresponding MD data related to the atomic description are shown as full and open circles.

Figure 6 .
Figure6.The real part of the wavevector-and frequency-dependent shear viscosity of the TIP4P water obtained within the modified GCM approach (solid curves) versus the MD data (circles) related to the atomic description.The MD results corresponding to the molecular description (at k = k min , dashed curve) as well as to frequency-dependent (k = 0, atomic description) shear viscosity (dots) are included in subset (a).

Figure 7 .
Figure 7.The imaginary part of the wavevector-and frequency-dependent shear viscosity of the TIP4P water.Other notations are the same as for figure 6.