Viscoelasticity of Short Polymer Liquids from Atomistic Simulations

The viscosity − or more generally the viscoelasticity − of polymer liquids is a key property for the processing as well as the performance of these materials. Molecular theories and numerical methods can provide these quantities, but they all require certain input parameters that nowadays are typically obtained by experiment. In the long term, it would be desirable to obtain these parameters or the whole viscoelastic response by purely computational methods, enabling a full “in silico” design of new materials and processes. In this perspective, we present several test calculations of the viscosity of n-hexadecane, a short-chain analogue of polyethylene. Our calculations are based on both equilibrium and non-equilibrium molecular dynamics (MD) simulations, which are applied to models based on a united-atom force field, a conventional atomistic force field, and the AIREBO-M reactive force field. We compare both the computational cost of the different strategies and the reliability of the different models and we provide some general guidelines for their application to more complex systems. © The Author(s) 2019. Published by ECS. This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 License (CC BY, http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse of the work in any medium, provided the original work is properly cited. [DOI: 10.1149/2.0371909jes]

Polymeric materials are typically processed in the molten state, for example by extrusion, blending, injection molding or 3D printing. 1 Very often they are then used in the solid state, as "plastics", but they may have also some uses in the liquid or liquid-like state, for example as lubricants and pressure-sensitive adhesives. Soft solids such as cross-linked elastomers ("rubbers") and polymeric gels represent an important intermediate situation between these extremes. 2 In all these cases, process and materials engineers are interested in characterizing, understanding and if possible predicting the viscosity (for liquids), the elastic modulus (for solids) and more generally the viscoelasticity of these materials. 3 To date, this remains a very challenging field. One of the main difficulties is that there are innumerable potentially relevant variables which can affect these properties. We mention for example the atomic-scale structure of the polymer chains (including regioregularity and tacticity), their overall composition (e.g., content and sequence statistics of co-monomer units), the large-scale structure of the chains (e.g., linear or branched), the molecular weight distribution of the chains, the possible presence of more than one polymeric species (i.e., blends), of small-molecule additives (e.g., plasticizers) or larger objects (e.g., spherical or anisotropic nanoparticles), the effect of temperature and pressure, and the "memory" of the processing history which is typical of these materials.
Given the high complexity of polymeric materials, it is clear that there are both great challenges and significant potential rewards in the development and application of predictive modelling concepts and tools. Most likely, these developments will involve a combination of theoretical concepts and numerical simulations, as well as the collection and analysis of "big data" and machine learning. [4][5][6] Polymer physics [7][8][9] has progressed up to a point where its concepts permeate the whole field of polymer processing and engineering. Mesoscale analytical theories 10,11 and numerical methods [12][13][14][15][16] can be used to obtain the rheological properties of complex, high molecular weight polymers. However, all these models contain system-specific parameters such as the monomer friction coefficient and the packing length. 17,18 These depend on the chemical structure of the polymer and they are typically obtained by experiment. In the long term, it would be desirable to be able to obtain these parameters also by computational means, in order to allow the "in silico" screening of new polymeric materials, formulations and processing conditions. The Molecular Dynamics (MD) method has contributed significantly to the understanding of the dynamics of molten polymers. The Kremer-Grest bead-and-spring model contains all the essential ingredients which are necessary to describe their universal, generic behavior. 19 It was originally applied to the verification of the reptation model for the equilibrium dynamics of polymer melts but, thanks to its simplicity, it was later applied to several other problems in polymer physics. In particular, its viscoelastic properties have been simulated by several authors using different methods. [20][21][22][23] When aiming at simulating the viscoelasticity of specific polymer systems, one must necessarily complicate the model and introduce some chemical details into it. System-specific, coarse-grained polymer models can be derived by a variety of approaches. [24][25][26] Typically, these methods can successfully reproduce selected structural and thermodynamic features of some underlying atomistic model. However, they usually do not reproduce its dynamical properties (e.g., diffusion coefficients and viscosities). The structural rearrangements within a coarse-grained model are usually faster than in the corresponding atomistic one, as the particles move on somewhat smoother potential energy surfaces. [27][28][29] In fact, one of the main reasons for the application of coarse-grained models lies in the fact that this dynamical acceleration can provide considerable computational savings. In principle, it is also possible to devise more sophisticated coarse-grained simulation methods (e.g., replacing Newton's equations by generalized Langevin equations with memory terms) which reproduce more faithfully the dynamics of a system. [30][31][32] However, considering the complexities involved in the parametrization and simulation of realistic coarse-grained models, and the continuous progress in computer hardware, one could say that there is a reasonable case for retaining an atomistic description when simulating the dynamical properties of molecular system, including polymeric ones. Our aim here is to explore the applicability of atomistic MD simulations to investigate short polymer liquids (below the entanglement length) and extract parameters which are relevant for their viscoelastic properties.
In principle, the viscoelastic properties of a molecular model can be extracted both from equilibrium MD simulations exploiting linear response theory, and from different variants of Non-Equilibrium Molecular Dynamics (NEMD). [33][34][35][36] Both classes of methods present B3247 advantages and potential pitfalls, and they have been extensively compared in relatively simple systems such as Lennard-Jones liquids and water. However, some of these advantages and pitfalls may be systemdependent and we feel it is still important to develop a better feeling for their performance in the context of atomistic simulations of more complex liquids. See Ref. 37 for a leading reference to the recent literature, with an application to room-temperature ionic liquids. In this communication we present a comparison of these methods, using nhexadecane (C16)-a short polyethylene-like molecule-as a model system.
We have computed the viscosity of liquid C16 using three realistic models with increasing levels of detail: a. a united-atom model (i.e. with implicit hydrogens) of alkanes, 38 originally developed for thermodynamic properties but later employed for viscosity calculations by Mondello and Grest 39 and by Cui et al.; 40 b. a fully atomistic model based on the OPLS-AA force field, 41 as modified by Böckmann and coworkers 42 to improve the description of long-chain hydrocarbon; c. a reactive empirical bond-order potential including intermolecular interactions (AIREBO), 43 whose non-bonded part was later modified to improve its performance at high pressures (AIREBO-M). 44 The use of a reactive force field needs some justification, since its computational cost is at least four times higher than that of a conventional atomistic force field. 45 What advantages could justify its use in viscosity calculation, which are per se computationally demanding? To answer that question, we point out that there is a considerable interest in modelling the viscosity of hydrocarbon fluids at high temperatures and pressures, 46 arising for example from the oil and gas industry. The whole field of tribology, which includes contact mechanics, sliding friction, lubrication and wear, is another field where the behavior of fluids at high temperatures and pressures is highly relevant. [47][48][49][50] Under such "ordinary but extreme" conditions, it is important to account for chemical reactions, such as shear-induced bond scission (tribochemistry). The search for alternative fluid mixtures for internal combustion engines is another field where molecular simulations could play a role, if both fluid viscosity and chemical reaction can be properly accounted for. 51 We do not expect chemical reactions such as bond scission to occur around room temperature and atmospheric pressure, as in our simulations. Nonetheless, as a starting point, it seems useful to test the viscosities of these reactive models 43,44 even under ordinary conditions. To the best of our knowledge, this has not been done before for hydrocarbon chains.
The plan of the paper is as follows. In the following section on Models and numerical methods we summarize the main equations and describe our computational procedures, including those for estimating errors and to calculate the autocorrelation functions of the stress tensor (an ingredient for the calculation of the viscosities from equilibrium MD). The Results and discussion section will allow us to compare the models and methods, both with respect to their quality (i.e., in comparison with the experimental viscosities of C16) and with respect to their computational cost. Conclusions follow.

Models and Numerical Methods
Viscosities from non-equilibrium molecular dynamics.-In NEMD, a steady shear deformation is explicitly applied to the system, and the shear viscosity is obtained from the average shear stress, as in actual experiments. Letγ = ∂v x /∂y be the gradient of the average (streaming) velocity within the fluid, and let xy be the timeaveraged component of the stress, which arises from the momentum flow along the gradient direction. The effective shear-dependent viscosity is obtained as: The Newtonian, zero-shear viscosity is obtained as the limit of this ratio whenγ → 0.
There are several ways to impose a shear flow. In LAMMPS, 53,54 a fully periodic triclinic box is continuously tilted at a desired rate, and rebuilt when its shear strain exceeds 0.5 (in modulus). The atoms' velocities are remapped whenever they cross the box boundaries. Atom positions are integrated using the SLLOD algorithm, 33,36 which uses a thermostat that subtracts the streaming component of velocities when updating the velocities of the atoms.
The instantaneous value of stress tensor is calculated from the expression: [2] where the summation runs over all the N atoms which are contained within the volume V (the expression given above must actually be slightly adapted in systems with periodic boundary conditions 52 ). Here m i are the atoms' masses, and v iα , r iα and F iα are the Cartesian components of their velocities, positions and total forces at time t. These ingredients are readily available at every timestep in an MD simulation. The above expression is general, as it does not rely on the assumption of pairwise additive interactions among the atoms.
Error estimates.-Error bars on the viscosites calculated from NEMD were estimated by the "renormalization group" (RG) method. 55 The method takes into account the fact that there are correlations between the successive data points-even when these are sampled once every few thousand time steps, as is often the case in MD simulations. Thus, the common statistical descriptors for uncorrelated data are unapplicable, in particular as they underestimate the error bars on their average. Without entering into excessive details, we report below the essential ideas and equations, using the notation of the original publication.
Let {x 1 , x 2 , …, x n } be a series of data extracted from a simulation (e.g., the stress values at different times). The time-average of the data is defined as usual: Such quantity is stochastic by itself, depending on the specific history of our simulation. Nonetheless, as long as the simulation is long enough, the time-average is a good estimate of the ensemble average over all the possible configurations. In principle one could assess the accuracy of our estimate m by repeating the simulation k times with completely different initial conditions, and looking at the variability of the various estimates. Let: where: Then the ensemble average over these independent simulations will have the standard error: It is also possible to assess the accuracy of m from a single simulation, provided that it is long enough. The main idea is that the system will loose memory of its previous state after a typical time τ, which we do not know a priori. If we split a long simulation of duration t into k chunks of duration t /k >τ, each chunk will be independent of the others and Eqs. 4−6 are still applicable. Conversely, if we divide the simulation into short time intervals (compared to τ), then the values are correlated. In the extreme case when the whole story is so short that nothing can change, all the k values of m are identical and the mean does not deliver any further accuracy. In this case the standard error must be the same as the standard deviation over a single chunk.
In principle, one could approach the problem by estimating τ from the correlation functions for these data (see below). As shown by Flyvbjerg and Petersen, 55 however, there is a simpler and more effective way to assess the accuracy of m, which does not require a knowledge of τ. It is based on the idea of recursive block averaging of the data, in order to estimate the effective amount of independent information carried by them. Let: We average the data in blocks of two: where n = [n/2] ([] is the integer part or floor function). Subsequently, a new c 0 is computed from the block averaged data using Eq. 7. This process of block averaging in repeated iteratively, until eventually n = 1. An estimate of the variance of the mean is obtained at each iteration, as: . [9] It can be shown than this estimate actually represents a lower bound to the true value. It increases at each iteration, but reaches a plateau when the process of block averaging has effectively eliminated all the correlation between successive data points. This plateau value, which ideally should persist for at least three iterations in order to be clearly identifiable, provides the correct estimate of the error on the average. Beyond this plateau, when the block averaging iterations have reduced the data points to a small number, our uncertainty on the variance of the mean [represented by the ± term in Eq. 9] increase again, essentially diverging. If the length of a simulation is too short, typically because the fuctuations in the x i variables are too large or too correlated, a plateau is never reached. In this case, the simulation must be extended in order to acquire more data points. Note that, in order to add one iteration to the block averaging process, one must effectively double the length of the initial simulation.
Viscosities from equilibrium molecular dynamics.-The Green-Kubo approach to the linear response of a system to a small perturbation is based on the evaluation of integrals of suitable autocorrelation functions. In particular, the zero-shear viscosity of a system can be obtained from the fluctuations of the stress tensor αβ (t ) recorded during a conventional NVT equilibrium MD simulation. In an isotropic system such as the one considered here, it is possible and convenient to use all its components, in order to reduce the statistical noise. Following Ref. 39, let us define the symmetrized traceless stress tensor: [10] The fluctuations of the trace of the stress tensor, which has been effectively eliminated from Eq. 10, determine the bulk viscosity of the liquid, which does not concern us here. The autocorrelation functions for the components of the P αβ tensor are defined as usual: where the angular brackets indicate an average over all possible time origins t 0 . Their evaluation is a critical step in the viscosity calculation and will be described in detail below. The relaxation modulus of the system, which could be measured in a step-strain experiment, is obtained by averaging all the components of these autocorrelation functions: Note that the summations in Eq. 12 run over all values of α = x,y,z and β = x,y,z. Because of the symmetry of the pressure tensor, there will be two replicas of the xy, xz and yz components and one replica of the xx, yy and zz components. The above expression, including the 1/10 normalization factor, was first derived by Daivis and Evans, taking into account the invariance properties of the stress tensor. 56 Given G(t), the zero-shear viscosity of the system can be calculated from the Green-Kubo formula: [13] It is also possible to compute the frequency-dependent reponse of the system to oscillatory shear deformations, in the form of the in-phase (real) and out-of-phase (imaginary) components of the relaxation modulus or of the viscosity: 1,3 The upper limit of integration Eqs. 13 and 14 must be chosen judiciously, and typically corresponds to a range of times at which the running integral reaches a plateau. Beyond that plateau, the value of the integral degrades due to statistical noise in the correlation function. Examples will be provided below, in the Results section.
As an alternative to the correlation function route, the viscosity can also be obtained from the long-time limit of the growth of the mean-square integrals of the stress components: [15] where: The equivalence of Eqs. 13 and 15 can be readily demonstrated using integration by parts. This is the so-called Einstein route, by analogy with his formula giving the diffusion coefficient of a Brownian particle as the proportionality coefficient between its mean-square diplacement and time. The Einstein route is simpler to implement than the Green-Kubo one, as its does not require the calculation of correlation functions. However, it is less powerful as it cannot be extended to the calculation of the frequency-dependent response, since there is no Einstein-like analogue of Eq. 14.
Numerical algorithms for the equilibrium viscosities.-As mentioned above, the integrals appearing in the Einstein route to the viscosity [Eqs. 15 and 16] can be evaluated readily from the values of the stress tensor output by the MD program during a production run. According to the well-known trapeizoidal rule for the numerical evalution of an integral: [17] in which we have N+1 values of the integrand function, equally spaced by intervals t, such that τ = N t. Notice that: a. the calculation of these integrals trivially scales linearly with the number of timesteps in a simulation (N); b. its computational cost can be further reduced by using batchaveraged values of the stress tensor [P αβ (t )] sampled every n time steps, instead of the instantaneous values [P αβ (t )] sampled every timestep: [18] which can be used within Eq. 17 with the coarser time step t = n t (the correction term with the 1 2 prefactor on the right-hand side is clearly negligible when N>10 6 , as in a typical MD simulation); c. the statistical accuracy of the viscosity calculated through Eq. 15 can be improved by averaging over different time origins, as implied by the angular brackets; d. in practice, we compute the limit of Eq. 15 by performing a linear fit of α β [L αβ (τ)] 2 against τ. Before performing the fits, we check that the slope of a double logarithmic plot of these data is unity. Note that the integrals in Eqs. 16 and 17 are the sum of a large number of correlated variables (i.e., the stress at different times). In order to observe a linear, diffusion-like growth of their mean-square value, these variables must effectively loose their correlation (notice the analogy with the statistics of an unperturbed polymer chain: in the long chain limit, its mean-square end-to-end distance grows with as a random walk, despite of the short-range correlations between the bond vectors 7,8,9 ). This occurs only when τ>>τ ee , where τ ee is the longest structural relaxation time of the system, obtained from the decay of the autocorrelation function of the end-to-end vector of the chains or, equivalently, from that of their first Rouse normal mode. 8,9 In practice, we typically perform our fit over the range [2τ ee ,3τ ee ].
The evaluation of the stress autocorrelation functions G αβ (t ) defined by Eq. 11 is more demanding, but also more rewarding as they contain information on the frequency-dependence of the response. The problem with the stress autocorrelation functions is that they display strong fluctuations at short times, associated high-frequency motions such as bond stretching vibrations, but it also has a long-time memory due to the persistency of the structure within a melt of chain molecules. In essence, the correlations in the stress tensor will decay to zero only at times larger than the characteristic time for the decay of the autocorrelations function of the chains' end-to-end vector. This requires frequent sampling of the stress, ideally every MD timestep, protracted up to very long times. As we shall see, in our system the number of data points N easily exceeeds 10 7 (for 40 ns of MD simulation), and this should increase further in systems with longer chains or slower dynamics. A straightforward calculation of the time correlations would require ∼N 2 operations, and this quickly becomes unmanageable from a computational point of view.
Different approaches have been proposed in order to tackle the evaluation of G αβ (t ). Some authors effectively suppressed all the short-time fluctuations, either by freezing all the "fast" degrees of freedom (fixed bond lengths and angles by the SHAKE or RATTLE algorithms 34,35 ) or neglecting their contribution to the stress tensor. 57 Other authors have adopted a "molecule-based" stress, instead of the "atom-based" stress tensor defined by Eq. 2. 39,40 These strategies reduce the computational burden by allowing a less frequent sampling of the stress tensor. However, they are not fully general. For example, they are unapplicable in a system described with a reactive force field, where bonds should be allowed to vibrate and, strictly speaking, molecules (i.e., aggregates of atoms with a fixed connectivity) are not defined. Also, in the case of a cross-linked polymer network, there is only one very large molecule, and therefore the molecule-based stress would be ill-defined. Discrete Fourier Transforms can also be used to reduce the computational cost in the evaluation of the stress autocorrelations, from ∼N 2 to ∼NlogN. 58 However, their use introduces an artificial periodicity in the data, and this can have a certain impact on the calculated long-time correlations. Also, the resulting correlation functions typically retain some high-frequency noisy oscillations which complicate the evaluation of the limit in Eq. 13. Some authors, such as Sen et al., 22 solved this problem by smoothing the correlation function at time t, perfoming a running average of its values from 0.9t to 1.1t.
Our approach to both problems (i.e., the computational cost of the evaluation the stress autocorrelation function and the high-frequency noise which affects its time integration) is based on the recursive block averaging of the stress components. This is different from the running averaging of the stress autocorrelation function mentioned above. 22 It is similar in spirit to the "order-N" algorithm discussed by Frenkel and Smith for the mean-square particle displacements 34 and very close to the "multiple-tau" correlator developed by Ramírez et al. 59 Their method has been incorporated into LAMMPS but, having become only recently aware of it, we developed our own stand-alone program (available from the authors on request). Borrowing the latter's notation, we define the block-averaged stress tensor (here we drop the αβ subscripts to avoid excessive cluttering of the equations): where the superscript l = 1,2,3,… indicates successive generations in the averaging process. The amount of averaging in the passage from one generation to the next is determined by the parameter m≥2. The original, non-averaged stresses correspond to generation l = 0, and they are typically sampled every MD timestep t: P (0,m) j = P( j t ). More generally, the subscripts in 19 represent discrete times in the sampling of the batch-averaged stresses: [20] The total number of data points at each generation is N (l,m) = [N/m l ]. Successive generations of block-averaged quantities are used to evaluate the correlation functions at longer and longer times: The range of the j subscript depends on the generation l, the block size m, and on a factor α 1 which determines the switching from one generation to the next: It can be seen from Eq. 20 that these j-ranges correspond to different intervals in physical time, wherein the correlations are calculated with different degrees of block averaging: [23b] Our program computes the autocorrelation functions from the stress values, saved at every timestep of a long MD run on an external file. The m and α parameters are specified as an input to the program. In most of the calculations presented below, we have used m = 2 and α = 200. This choice of α is very conservative, as it implies a large delay in the switching from one generation to the next. As we shall see, the calculated viscosities are well converged with respect to this parameter. Starting from l = 0, the program computes the autocorrelation function up to t (0,m) Models and MD simulation details.-As stated in the Introduction, MD simulations were carried out with three distict models for C16: a united-atom (UA) force field 38 with implicit hydrogens, an all-atom force field based on a refinement of OPLS-AA for longchain hydrocarbons (L-OPLS), 42 and the reactive AIREBO-M force field. 44 We point out that the original OPLS-AA 41 and AIREBO 43 parametrizations were also tried but later discarded, as we observed that they both lead to the crystallization of C16, even at a temperature (298 K) above its normal freezing temperature (291 K). 60 Note that, since crystallization is a first order phase transition which requires a rare nucleation event, simulations started from the liquid phase normally tend to overestimate the (meta)stability of this phase. Thus, even though 298 K is not very much above 291 K, the fact that crystallization is observed so easily is a strong indication that those two models are not adequate in the present situation.
All systems contain 250 chains, with full three-dimensional periodic boundary conditions. This is more than twice the number in previous viscosity calculations with the UA force field. 39,40 This choice reduces finite-size effects which can be significant when the molecular size is comparable to the box size. 61 In the UA and L-OPLS simulations, we used a 12 Å cutoff for the non-bonded Lennard-Jones interactions. Instead, AIREBO-M uses a Morse potential to represent the non-bonded interactions. Following the original publications, 43,44 we truncated each interaction at a distance equal to three times its conventional hard-core value (3.40 Å for CC, 3.025 Å for CH, 2.65 Å for HH). Long-range electrostatic interactions are present only in the L-OPLS model. They are weak-each hydrogen has a charge of +0.074, and each carbon has a negative charge which neutralizes that of the hydrogens bonded to it-nonetheless they were treated accurately by Ewald summation, with a relative error on the coulomb forces of 1 × 10 −4 . The MD equations were integrated with a 2 fs timestep in UA and L-OPLS simulations. No intramolecular constrainst were applied, except for C-H bonds lengths which were fixed at their equilibium value. A timestep of 1 fs was used in the AIREBO-M simulations, without any constraints. Two representative snapshots of the UA and L-OPLS systems are shown in Figure 1.
Each system was first equilibrated by NPT simulations with a Nose-Hoover thermostat and barostat, respectively with relaxation times τ T = 0.2 ps and τ P = 2 ps. These simulations provided the equilibrium densities at P = 1 atm, which were then kept fixed for the subsequent equilibrium and non-equilibrium production runs. These were carried at temperatures T = 298 K, 323 K and 373 K. All simulations were carried out with the free, open-source Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) code. 53,54 Further details about the numerical methods can be found within LAMMPS's extensive manual.

Results and Discussion
Thermodynamic and single-chain properties.-We start our discussion from the equilibrium properties of our models of C16. We have collected in Table I some selected structural and dynamical data. The densities and their temperature dependence can be considered a simple measure of the quality of the molecular packing and of the strength of the intermolecular interactions. All the models slightly underestimate the experimental densities, typically by 2%. This is in line with many other simulations of hydrocarbon liquids can be considered satisfactory. The differences between the individual models seems to be negligible, in this respect.
The mean-square end-to-end distances of the molecules ( R 2 ee ) cannot be directly measured experimentally. Nonetheless, they are useful as they summarize within a single number the conformational properties of molecules. A lower value of the end-to-end distance indicates more flexible chains, with a higher population of gauche conformations with respect to the more extended trans conformations. In turn, these depend mostly on the torsional potential entering the models. In order of increasing rigidity, we have L-OPLS < AIREBO-M < UA. The end-to-end distances shrink slightly at higher temperatures, as expected for an increasing population of the higher energy gauche states. Again, however, these differences appear to be small. Taken together, these data indicate a similar behavior of these models with respect to the "structural" properties of the models, both at the macroscopic level (densities) and at the molecular level (end-to-end distances). The other quantities listed in Table I reflect the dynamical behavior of these sistems, at the single-molecule level. The self-diffusion coefficients D describe the random wandering motion of the molecules within the liquid. Experimentally, they can be measured by a variety of methods, including nuclear magnetic resonance with pulsed magnetic field gradients (PFG-NMR). 63 See Ref. 64 for relatively recent data on C16 and Ref. 65 for a review with an extensive compilation of experimental results on a broad range of molecular liquids (excluding polymers). In a simulation, the self-diffusion coefficients D are obtained either by integration of the autocorrelation function of the molecular velocities or from the ratio between the mean-square displacements (MSD) of the chains' center-of-mass (CM) and time. 34,35 We have used the second approach, which correponds to the Einstein relation: [24] [note the analogy with Eq. 15; here we have used the double brackets to indicate a double average, namely over all time origins and all C16 chains]. It is relatively straightforward, but care must be taken to avoid fitting the numerical data at too short times, as there can be "anomalous diffusion" when the molecules still feel the effect of the transient cage formed by their neighbors and their MSD does not scale linearly with time. 66,67 See for example Ref. 68 for a joint experimental and computational study, evidencing these phenomena in room-temperature ionic liquids. Our present simulations, which are tens of nanoseconds long, are long enough to be in a fully linear regime ("ordinary diffusion") and they provide us with a reliable estimate of the diffusion coefficients characterizing each model. The autocorrelation functions for the end-to-end distance of the chains are shown in Figure 2. Their decay is well described by a simple exponential function: The end-to-end relaxation times (τ ee ) represents the characteristic time required by a chain to undergo a global conformational rearrangement and loose all memory of its initial configuration. As discussed in the previous section, its knowledge is crucial for both the equilibrium and non-equilibrium simulations of the shear viscosities. According to the Rouse model for the dynamics of non-entangled polymer melts, 7,8,9 it should be comparable to the relaxation time of the first Rouse normal mode (τ R ) and to the diffusion time (τ D ).
The Rouse normal modes (RNM) are "conformational waves" which describe the overall relaxation of the chains. Mathematically: 0, 1, . . . ., N − 1) , [26] where the summation runs over the N backbone atoms of a chain. Note that mode p = 0 corresponds to the chain's CM. When p = 1,2,3…(<<N-1), their decay should also follow a simple exponential, with a characteristic time which decreases as the square of the p index: The plots for the decay of the mode p = 1 (not shown) are undistinguishable from those for the end-to-end vector, so that we may indeed identify τ R with τ ee . The diffusion time is defined as the time required for a chain to diffuse over a distance comparable to its size, which may be taken equal to the end-to-end distance. From Eq. 24 we have: . [28] Comparing the self-diffusion coefficients and relaxation times listed in Table I, we notice systematic differences among the C16 models. At room temperature, the rates of single-chain relaxation decrease in the order UA > L-OPLS > AIREBO-M. The latter two are in fact comparable, and they agree much better with the experimental diffusion coefficients. The differences among the models decreases on going to higher temperatures, as the UA model is less sensitive to temperature changes than the other two. This temperature dependence is illustrated in Figure 3, which contains Arrhenius plots of the diffusion coefficients. The UA model has a weaker temperature dependence than experiment, indicating a lower activation energy for diffusion. The predictions of the atomistic L-OPLS and AIREBO-M models are closer to experiment, but now with a slightly higher activation energy.
Viscosities from non-equilibrium MD.-We now turn to the discussion of the calculated viscosities, starting from those obtained from the NEMD simulations. These were limited to the UA and L-OPLS models at 298 K, owing to the higher computational cost of this method (see next for a comparison of the equilibrium and non-equilibrium methods, in this respect). Figure 4 (upper panel) summarizes our results. The simulations were carried out at five different shear rates, spanning three decades. These rates were determined on the basis of the end-to-end relaxation times. The inverse of these times are represented by the vertical dashed lines in Figure 4. The dimensionless quantity Wi =γτ ee is actually know in rheology as the Weissenberg number, 69 and our NEMD data points go from Wi<<1 to Wi>>1. Around Wi = 1, which corresponds to the dashed lines, we have the transition to the non-Newtonian shear-thinning regime. Note also that the error bars on the individual data points become larger on going to lower shear rates, when the shear stress fluctuations become   Table I). The inset in the upper figure shows the same data, using the reduced viscosity η(γ)/η 0 on the y-axis and the Weissenberg number Wi =γτ ee on the x-axis. comparable to its average value. The RG procedure used for the estimation of the error bars is illustrated in the lower panel of Figure 4. We see that there is indeed a plateau in the estimated error when the blocking process arrives at generations 13-15. For this particular point, our estimate for the variance of the mean [Eq. 9] is thus σ m ≈0.35 cP.
The non-equilibrium viscosity data can be fitted by the Carreau function, in analogy with the work of Bair et al.: 70 The η 0 parameter in this equation represents the zero-shear viscosity, τ C is a time which determines the characteristic rate for the transition to the non-linear regime, and n is a parameter that, roughly speaking, determines the sharpness of the viscosity drop at the transition. Their values are: η 0 = 1.38, τ C = 0.299 ns, n = 0.616 (UA model) and η 0 = 4.34, τ C = 0.442 ns, n = 0.501 (L-OPLS model). The η 0 parameters are also listed as η NEMD in Table II. The UA force field underestimates the viscosity of C16 at room temperature by 55%, while the L-OPLS force field overestimates it, but by a smaller amount. This trend will be confirmed and expanded by the equilibrium viscosity calculations. Concerning the other parameters of the Carreau function, the fact that n exponent of L-OPLS is appreciably smaller than that of UA indicates that it is not only more viscous, but also "more non-linear". Notice also that the ratio between the characteristic times of the models is different: τ ee (L-OPLS)/τ ee (UA) = 2.64 but τ C (L-OPLS)/τ C (UA) = 1.41. Indeed, the rescaled flow curves in the inset [ Figure 4, upper panel] appear to have different shapes. All this indicates the possibility of fundamental differences between the models, so that one is not simply an "accelerated" version of the other, as different processes may have different degrees of acceleration.
In general, the non-Newtonian behavior of a fluid reflects some change in its "structure" that is produced by the flow. In the present case, this structure corresponds to the orientational order of the chains. The lower panel of Figure 5 contains two instantaneous snapshots and the plots of the distribution of the molecular orientations at two shear rates, for the UA model. The molecular orientations are calculated as the cosine (in modulus) of the angle between the instantaneous endto-end vector of each molecule and the average director. Borrowing from the theory of nematic liquid crystals, the latter corresponds to the where u iα is the α-component of the unit vector pointing along the direction of molecule i. The average of the largest eigenvalues yields the orientational order parameter, P 2 . The figure demonstrates a dramatic change in the distributions, which goes from isotropic to strongly oriented when the Weinssenberg number Wi =γτ ee becomes greater than unity. The upper panel of Figure 5 allows us to appreciate the evolution of the average Q xy and P 2 as a function of the average shear stress xy . At very small stresses, in the Newtonian viscous regime (shear ratesγ < 1 ns −1 for the UA model), the chains have no appreciable orientation. In the intermediate stress regime, when the non-Newtonian behavior starts to manifests itself (γ = 1 − 10 ns −1 ), the chains develop an orientation which is approximately proportional to the stress. This proportionality corresponds to the so-called stressoptical law. At even higher stresses, the proportionality is clearly lost as the chains' orientation approaches saturation.
Viscosities from equilibrium MD.-As explained in the Methods section, the zero-shear viscosities can be extracted from a conventional equilibrium MD simulation in the NVT ensemble, following either the Einstein route or the Green-Kubo route. Their application in illustrated in Figure 6. To apply the Einstein method (upper part of Figure 6), one must fit the mean-square value of the time integral of the stress tensor, taking care to be in a "diffusive" regime where this grows linearly in time [Eq. 15]. Following the recommendations of Mondello and Grest, 39 we typically fit the Einstein relationship for the stress over the range [2τ ee ,3τ ee ]. This range satisfies several conditions: its lower limit ensures that the "dynamics" of the stress tensor is well into the diffusive regime, and its upper limit is separated enough so that there will be enough data points to perform the fit with sufficient accuracy. At the same time, the upper limit does not stretch so far that there is a deterioration in the data due to poorer statistics. Unlike the calculation of the molecules' self-diffusion coefficient, for which at each timestep there are as many data (center-of-mass coordinates) as molecules, in a viscosity calculation there are only five independent stress components that can be used. Thus, the accumulation of good statistics is a much more serious problem in the latter case. It is for this reason that the overall length of the equilibrium MD simulations must be 50-100 times longer than τ ee for each model.
According to the Green-Kubo method, one must first compute the time-correlation function of the stress tensor, and then evaluate its integral according to Eq. 12. As discussed in the "Methods", the use of a "coarse-graining in time" approach simplifies both tasks. The lower part of Figure 6 shows the value of one such integral, as a function of the upper limit of integration. As can be seen, this essentially reaches a plateau around τ ee . Instead of singling out one specific time to stop the integration, which would be slightly problematic, we have found it convenient to average all the values in the interval [τ ee ,2τ ee ]. The plot shows also that the result of the integration is remarkably insensitive to the choice of the "delay parameter" α, which was introduced around Eq. 22. This parameter has a strong influence of the integrand function, G(t), as shown in Fig. 7. The smaller α, the stronger the suppression of the stress fluctuations at short times. However, integration up to times of the order of τ ee effectively washes out these differences. The other parameter entering the evaluation of G(t) is the block size, and this was kept at the smallest possible value m = 2.
The zero-shear viscosities extracted from the equilibrium simulations using either route are collected in Table II. They are generally consistent with each other and with the values extrapolated from the NEMD simulation. At room temperature, the best agreement with experiment is obtained with the AIREBO-M force field. The UA model underestimates the viscosity, while the L-OPLS model overestimates it. However, if we consider all temperatures, the L-OPLS force field seems to provide the best overall description.
Like the diffusion coefficients, also the viscosities have an Arrhenius-like temperature dependence: The Arrhenius plots in Figure 8 show a stronger temperature dependence for the AIREBO-M model, implying a higher activation energy. Instead, the slope of the plot for the L-OPLS model is similar to the experimental one and, as pointed out above, this guarantees a better description also at high temperatures.
Our values for the viscosities of the AU model are consistently lower (by 40-60%) than experiment. This result agrees with the diffusion data and indicates a deficiency in the force field, which is remedied by adopting a more complex and detailed, atomistic description. Our viscosity values are also systematically lower (by 15-25%) than those reported years ago by Mondello and Grest 39 and Cui et al. 40 for the same UA model, at the same temperatures. These differences could be explained by the fact that we are using larger system sizes (250 versus 100 molecules), longer simulation times, and the fact that our liquid densities were obtained from NPT simulations, instead of being taken from experiment.
Despite being difficult to implement, and despite leading to very similar results for the shear viscosity, the Green-Kubo route is preferable to the Einstein route. The reason is that, once we have computed the stress relaxation function G(t), we can use it obtain the whole viscoelastic response of a liquid, in the linear regime. The in-phase (real) and out-of-phase (imaginary) components of the frequency-dependent modulus (G' and G") and viscosity (η' and η") can be readily calculated using Eqs. 14a and 14b. We have computed and plotted the latter in Figure 9, for the three models at 298 K. The simulation results are given alongside the predictions of the Rouse model, which is appropriate for a melt of short, unentengled polymer chains: 8 where c is the molar concentration of chains, so that c/N is the molar concentration of backbone atoms (N = 16). Eq. 32 suggests to rescale the components of viscosities by the zero-shear value η 0 [lim ω→0 η (ω) = (π 2 /12)cRT /N, for the Rouse model] and the frequencies by the inverse Rouse time (the τ R values are virtually equal to the τ ee values listed in Table I).
After rescaling of the viscosity and frequency axes, the low frequency responses of the different models appear rather similar, as η'(ω)→η 0 and η"(ω) ∼ω when ω→0. Notice, however, that the ratio η'(ω)/η"(ω) [equal to tanδ(ω) = G"(ω)/G'(ω)] is higher than predicted by Rouse. This means that, at low frequencies, the simulated models dissipate more energy than predicted by Rouse. This is understandable since the latter is based on a bead-and-spring description and does not include the effect of chain interactions and "non-polymer" degrees of freedom. 23 A decay of η'(ω) starts to occur at angular frequencies of the order of 2π/τ R . Apart from the 2π factor, this is also the value ofγ beyond which we observe the shear thinning in the NEMD simulations ( Figure 4). This empirical parallelism between the non-linear response to a steady shear flow and the linear response to an oscillating flow is known to rheologists as the Cox-Merz rule 1 ).
Here the AIREBO-M model deviates significantly from the other two and from Rouse, as it displays a maximum in η' where the others display a minimum, and vice versa for the η" components. The different behaviors could be associated with differences in the torsional barriers, since these become important when the internal modes of such short chains start to make a contribution to the complex viscosity. 72 Deviations of all the models from the Rouse prediction become even more important at high frequencies, as the higher modes relax in a way which is not a simple exponential and their relaxation times deviate from the 1/p 2 scaling of Eq. 27. For a molecule such as C16, this frequency range would correspond to ω>1 rad/ns. This not accessible to conventional rheological experiments, but it can be investigated by techniques such as Quasi-Elastic Neutron Scattering (QENS) and Broadband Electric Spectroscopy (BES). Thus, these techniques could provide additional tools to validate the different force fields, and vice versa the simulations could be used to aid the interpretation of the experiments.
Comparison of models and methods.-Finally, we briefly discuss the computational cost of our simulation strategies (i.e., equilibrium versus non-equilibrium MD) and computational models (i.e., UA versus L-OPLS versus AIREBO-M). Some representative data have been collected in Table III. In the "time" columns, we have collected the overall length of each simulation, in ns. According to this parameter, the equilibrium simulations appear to be superior to the nonequilibrium ones for the viscosity. The main problem with the latter occurs at very small shear rates, as the amount of simulation time which is necessary to reduce error bars to an acceptable level seems to explode. The problem is very severe for the L-OPLS model, as we had had to simulate it for nearly 0.5μs in order to obtain the point atγ = 0.1 ns −1 . Thus, having solved the problem of computing the autocorrelation function of the stress tensor-without issues related to storage disk space, in the "on-the-fly" version of the method devised by Ramírez et al. 59 -equilibrium methods appear to be clearly superior for the evaluation of the ordinary, zero-shear viscosity of hydrocarbon liquids. We have demonstrated that the two atomistic models are superior to the UA one and, encouragingly, they lead to viscosities that are rather close to experiment, at least around room temperature. But how much more do they cost, in computational terms? To compare them, we have performed standardized simulations on a computer node equipped with 16 CPU's (note that, for the atomistic simulations, we exploited also more powerful supercomputing hardware provided by CINECA). The parameters for these benchmark simulations (timestep, interaction cutoffs, etc.) were identical to those described in the Models and numerical methods section. The results are reported in the "cost" columns of Table III, as the number of wall-clock hours required to complete 1 ns of simulation. There is a factor of 15 on going from a UA simulation to an L-OPLS simulation, and this is essentially due to the larger of atoms (roughly tripled, if we count two hydrogens for each carbon) and to the presence of long-range interactions, due to the (small) charges on the C and H atoms. The cost of an AIREBO-M simulation is roughly four times that of a conventional atomistic simulation, confirming previous reports. 45 In fact, roughly half of this increase in computational cost is due to the use of a shoter time step (1 fs instead of 2 fs). The evaluation of the bonding interactions is obviously much more complex with a reactive force field, but this is balanced by the absence of electrostatic interactions in the AIREBO model. This feature should make its even more advantageous on massively parallel computers, using a code such as LAMMPS that is based on a domain decomposition strategy which is more suitable for short-range forces. 54 The overall cost of a viscosity calculation depends both on the intrinsic "sluggishness" of a model and on the complexity of the calculations required at each timestep. The combination of these factors is summarized in the "time × cost" column of Table III. Rougly, for n-hexadecane and similar hydrocarbon liquids there is nearly an order of magnitude increase on going from the UA model to an atomistic one, and from an equilibrium calculation to a non-equilibrium viscosity calculation (if we are aiming to derive the full "flow curve", such as the ones shown in Figure 4).

Closing Remarks
We have presented a comprehensive comparison of models and methods which can be used to derive the viscosity, or more generally the viscoelasticity, of organic liquids and short polymer melts from molecular simulations. The quality of the results for n-hexadecane, a molecule with a size which is typical of lubricants that we have used as benchmark, is encouraging with respect to the perfomance of two atomistic force fields (L-OPLS and AIREBO-M). On the other hand, the viscosities obtained from a united-atom model, which is much cheaper from a computational point of view, are typically half of the experimental values. Interestingly, our simulations indicate that the UA model is not simply an "accelerated", less viscous version of the atomistic ones. For example, its flow curve cannot be superimposed to that of L-OPLS, even after rescaling by the respective end-to-end relaxation times and the zero-shear viscosities. Our results open the way to further atomistic simulations of the viscosities of other hydrocarbons and organic liquids, including mixtures and "extreme" situations (e.g., high temperature and pressures) such as those that occur in lubrication and friction. Hopefully, these simulations will help us also to gain a better understanding of the viscoelasticity of high polymer melts, for example by providing input parameters to more coarsegrained simulation methods based on entanglements and slip-links.
The careful application of both the non-equilibrium and the equilibrium (Einstein or Green-Kubo) routes 33,36 to the viscosity of nhexadecane has provided numerically consistent results. The latter seem to be preferable for their computational economicity, when one is simply interested in the Newtonian, zero-shear viscosity of a fluid. The Green-Kubo method, whose application has been greatly simplified by devising better algorithms to compute time autocorrelation functions, yields also the full linear, frequency-dependent response of a liquid. On the other hand, the non-equilibrium simulations can provide insights into the shear-induced changes in the structure of a fluid which are at the basis of its non-Newtonian behavior. Thus they are still quite valuable, despite their higher computational cost. In the future, we hope to continue working with both methods, pursuing further methodological improvements. For example, it would be interesting to extract information about the stress fluctuations also from the non-equilibrium simulations, and use it to improve the extrapolation to the zero-shear viscosities. In the equilibrium case, one could use the particle-field approach 73 to generate efficiently many independent configurations, to be used as an input to the variant of the Green-Kubo approach recently proposed by Zhang et al. 37