{'fraction_non_alphanumeric': 0.03264497288277013, 'fraction_numerical': 0.017417605340008344, 'mean_word_length': 4.697563874034462, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Shear viscosity of n-pentane from diffusivity based molecular dynamics methods

The molecular dynamics methods for calculation of shear viscosity based on liquid diffusivity are tested against the classical Green–Kubo relation for n-pentane at 330 K and 0.601 g/cm3. The D-based method is shown to be as accurate as calculation of viscosity from the Green–Kubo formalism for pentane liquid. Stokes–Einstein relation is also in agreement with simulation results. However, it has much bigger uncertainty. The results of the simulation are in agreement with experimental data.


Introduction
The shear viscosity of liquid is a property that is important in many technical applications. There are conditions, when the experimental observation of this property becomes difficult [1][2][3][4]. In these cases, the molecular dynamics methods allow to predict this quantity quite accurately.
The equilibrium method of the shear viscosity calculation from molecular dynamics is usually based on the Green-Kubo (G-K) formalism. Viscosity is calculated as an integral of the stress auto-correlation function (SACF). For the atomic and simple liquids, the integral converges without problems [5][6][7][8]. This becomes a challenge for the molecular systems due to the high SACF decay time. In 2015, Maginn and coauthors proposed a technique to solve the convergence problems in the G-K method for highly viscous liquids [9]. This technique is called time decomposition method (TDM) and has been successfully applied in recent works for ionic liquids [10] and hydrocarbons [11]. The accurate error estimation of the G-K method for Lennard-Jones liquid has been done by Kim and coauthors in 2018 [12]. After all, the G-K approach is limited by the accuracy of SACF calculation and requires a lot of statistics.
On the other hand, there is an approach to the calculation of viscosity based on a change in the diffusion coefficient with the size of the modeling box. The analytical correlation for the correction of the diffusion coefficient by size was obtained theoretically and proved using molecular dynamics in [13]. This is also applied for the mutual diffusion coefficients [14] and rotational diffusion [15]. It was originally derived from the Navier-Stokes equations to provide accurate estimation of diffusivity that takes into account effect of boundary conditions in the simulation. Then it was used to calculate viscosity of some simple liquids [16]. Recently, D-  [17] and to use as a universal technique for viscosity calculation.
There is an approach of viscosity calculation based on the Stokes-Einstein relation, which is used in many theoretical models of diffusivity [18][19][20][21]. It connects liquid viscosity, diffusivity and size of liquid particles. It is used for rough estimation of the viscosity when accurate Green-Kubo simulations are not possible. The example of such situation is screening of ionic liquids [22].
Here, we test these equilibrium approaches against each other to get the estimation of their accuracy and obtain self-consistent value of n-pentane viscosity.
The article is structured as follows. In section 2, the equilibrium methods of shear viscosity calculation are described. The details of molecular dynamics simulation are presented in section 3. Section 4 contains the comparison of different viscosity computational methods.

Equilibrium methods for viscosity calculation 2.1. Green-Kubo integral
The shear viscosity η αβ can be calculated using the Green-Kubo formula [23]: where C σ (t) is an autocorrelation function; σ αβ are off-diagonal components of the stress tensor; V and T are system volume and temperature; k B is the Boltsmann constant; . . . in equation (2) is an average over the canonical ensemble; t ′ in (1) is a value of time when the autocorrelation function asymptotically decays to zero. The shear viscosity η is found as an average of the η xy , η xz and η yz . The stress tensor σ αβ is calculated from the next equation: where N is a number of atoms; r iα and v iα are α-components of coordinate and velocity of the i-th atom, and f iα is an α-component of the force that acts on the i-th atom.

D-based method
Liquid diffusivity depends on the simulation box size in molecular dynamics even when periodic boundary conditions are applied. This dependence is obtained analytically from hydrodynamics [13]: where D ∞ is the diffusivity at infinite system limit; D L is the diffusivity at simulation box size L; T is a temperature; η is a shear viscosity; k is the Boltzmann constant; ε ≈ 2.837 298 is a numerical constant. Increase in diffusivity with system size is also observed in nanopores [24] and hydrocarbons [25]. Deviations from this relation are only observed for rather specific systems [26][27][28]. Viscosity could be estimated from the slope of the dependence of D on system size L using relation (4).
The diffusivity D L at the system size L is computed from the Einstein relation where r 2 is a mean squared displacement of particles.

Stokes-Einstein relation
Viscosity is calculated directly from the diffusivity where R is the radius of a molecule; C is a constant that depends on the boundary conditions between the diffusing molecule and surrounding liquid; C = 6π for the stick boundary and C = 4π for the slip boundary conditions. Radius of gyration (R g ) [29,30] is used instead of R because pentane molecule is not spherical. Radius of gyration is defined by the relation where M is a mass of the molecule, r cm is a center of mass coordinate of the molecule. Size corrected (4) value of diffusivity is used.

Details of molecular dynamics simulation
The OPLS-AA (optimized potentials for liquid simulations-all atom) [31] all-atom force field is used to parameterize energy equations. Authors have compared another force fields by the predictive power of equation of state and diffusivity in previous work [20,32]. Since here relatively short n-alkane is studied, there is no need to use modified L-OPLS-AA [33] version developed for long hydrocarbons. All viscosity calculations are performed for balanced systems. Non-bonded 1-4 interactions are scaled with factor of 0.5. The cut off radius is 12Å. Long-range Coulomb interactions are calculated using the particle-particle particle-mash method [34]. The equations of motion are integrated using velocity-Verlet integrator, the timestep is 1 fs. The simulations are performed in periodic boundary conditions using the LAMMPS package [35][36][37].
The initial configuration is created with randomly distributed molecules. The density is 0.601 g/cm 3 . Then, the energy minimization is performed to avoid overlapping atom cores and N V T relaxation with Langevin thermostat is applied for 50 ps and then N V T relaxation with Nose-Hoover thermostat for 200 ps. The temperature is 330 K. After that, the calculations using G-K, D-based and S-E methods are performed.
For the G-K case, 1000 molecules are considered in the box. In the D-based method, we run simulations for 125, 512, 1300 and 35 000 molecules in the box.

Results
The results of the molecular dynamics simulation with different methods of viscosity estimation are presented in table 1. The results are in overall agreement with each other and experimental data.
The Green-Kubo method is the most well established approach for viscosity calculation. In this work, it is used as a reference value together with the experimental data [38] since the experimental value may differ from the simulation result due to the approximate nature of the interaction potential. The Green-Kubo value of viscosity is obtained using TDM method [9] via 100 statistically independent 1 ns trajectories.
The Stokes-Einstein relation has a rather small statistical uncertainty. On the other hand, it is important that this relation is approximate and there are a lot of different approximations for it. The simplest change frequently applied in the literature is a variation of boundary conditions between the diffusing molecule and surrounding liquid. This changes the value of C in (6). This leads to a wide range of possible values of the diffusivity with this method as shown in table 1. It makes the error of S-E value at least 50% due to the limitations of the method. Another problem is a determination of hydrodynamic radius of the molecule that might significantly differ from the radius of gyration.   190 - The D-based method is rarely used method for the viscosity calculation. However, it gives values of diffusivity that are in agreement with the Green-Kubo integral within the statistical uncertainty. The results of this method are taken from the slope of the dependence of diffusivity on the inverse simulation box length (figure 1).

Conclusions
The shear viscosity of liquid n-pentane at T = 330 K and ρ = 0.601 g/cm 3 is calculated from molecular dynamics via different approaches. The OPLS-AA force field is used. The results obtained from different methods are in agreement within statistical error.
The D-based method shows the same accuracy as well established Green-Kubo methods of viscosity calculation. Difference between these results is close to the statistical uncertainty of the simulation.
Stokes-Einstein relation gives accurate result when stick boundary condition is applied. However uncertainties that arise from the theoretical formulation of this relation are more significant than uncertainties of G-K and D-based methods. The results of the simulation are in agreement with the experimental data. This fact verifies the OPLS-AA force field and the used approaches.