Shear Viscosity Computed from the Finite-Size Effects of Self-Diffusivity in Equilibrium Molecular Dynamics.

A method is proposed for calculating the shear viscosity of a liquid from finite-size effects of self-diffusion coefficients in Molecular Dynamics simulations. This method uses the difference in the self-diffusivities, computed from at least two system sizes, and an analytic equation to calculate the shear viscosity. To enable the efficient use of this method, a set of guidelines is developed. The most efficient number of system sizes is two and the large system is at least four times the small system. The number of independent simulations for each system size should be assigned in such a way that 50%-70% of the total available computational resources are allocated to the large system. We verified the method for 250 binary and 26 ternary Lennard-Jones systems, pure water, and an ionic liquid ([Bmim][Tf2N]). The computed shear viscosities are in good agreement with viscosities obtained from equilibrium Molecular Dynamics simulations for all liquid systems far from the critical point. Our results indicate that the proposed method is suitable for multicomponent mixtures and highly viscous liquids. This may enable the systematic screening of the viscosities of ionic liquids and deep eutectic solvents.

Introduction

The shear viscosity plays an important role in quantifying the required energy for mechanical and chemical processes, and it is essential for solving the Navier–Stokes equations. Various empirical/semiempirical models for predicting the shear viscosity have been developed, including equations of state combined with scaling relations,[1−11] free-volume theory,[12] friction theory,[13] as well as fitting to functional forms.[14,15] Alternatively, the shear viscosity can be computed via equilibrium or nonequilibrium molecular dynamics (MD) simulations.[7,8,16−24] In nonequilibrium MD (NEMD), the viscosity is calculated from the response of the system to an external shear field.[18,25−29] Although NEMD is computationally efficient at large external fields,[17,30] the computed viscosity can depend on the applied shear rate.[20,21] In equilibrium MD (EMD), the shear viscosity can be computed from the Einstein relation:[31−34]where Pαβ represents the off-diagonal components of the stress tensor (i.e., P, P, and P), kB is the Boltzmann constant, t is the correlation time, and V and T are the volume and temperature of the system, respectively. The angle brackets denote an ensemble average. Since the stress tensor is defined for the entire simulation box, the shear viscosity is a property of the system as a whole. This means that an increase in the system size does not improve the statistical uncertainty of the computed shear viscosity.[35] Due to the slow relaxation of highly viscous fluids such as ionic liquids[14,36] and deep eutectic solvents,[37] and due to the large fluctuations in the components of the stress tensor[17] (see Figure S1 in the Supporting Information), very long MD simulations are required to sufficiently sample the stress tensor components.

Unlike the shear viscosity, the self-diffusivity is a single-molecule property and can be calculated from the mean-squared displacement (MSD) of all individual molecules of the same species:[34,38−40]where r(t) is the position of the jth molecule of species i at time t, and N is the number of molecules of species i in the system. The statistical uncertainties of self-diffusivities decrease as the number of molecules in the system increases.[35,41,42] The simulation length needed to obtain a linear relation between the MSD and time is much smaller than the length of an MD simulation required for computing the shear viscosity. Therefore, accurate self-diffusivities can be computed from short MD simulations with a large number of molecules. In the same short simulation, such a high accuracy cannot be achieved for the shear viscosity due to the smaller number of samples for the stress tensor components.[35]

Typical MD simulations for computing transport properties use hundreds to thousands of molecules, which is orders of magnitude smaller than the thermodynamic limit. To overcome this issue, simulation boxes with periodic boundary conditions are used.[38] While this approach mimics the presence of an infinite bulk system surrounding the simulated system, computed properties may depend on the size of the simulation box. Large finite-size effects for thermodynamic properties such as activity coefficients,[43] chemical potentials,[44] and Kirkwood–Buff coefficients[45,46] have been reported in literature. In the case of transport properties, self-diffusion coefficients computed from MD simulations depend strongly on the system size while the shear viscosity does not show any system-size dependency.[47−49] Dünweg and Kremer[50] showed a linear increase of the self-diffusivity of polymers in a good solvent with the inverse of the simulation box length. Yeh and Hummer[47] investigated the system-size dependency of computed self-diffusivities of pure liquids, consisting of spherical molecules, and derived a relation between the infinite (Dself∞) and finite-size (DselfMD) self-diffusivity:[47]where DYH is the finite-size correction, here referred to as the Yeh–Hummer (YH) correction, L is the length of the simulation box, which is proportional to N1/3, where N is the number of molecules, η is the shear viscosity of the liquid, and ξ is a dimensionless constant equal to 2.837297 for a cubic simulation box with periodic boundary conditions.[47,51−53] For the derivation of this analytic correction, Yeh and Hummer[47] used a hydrodynamic model of a spherical particle diffusing in a medium of viscosity η in a box with periodic boundary conditions imposed in every direction. In this model, the hydrodynamic interactions of the particle with the surrounding molecules and the periodic images affect the self-diffusion coefficient. The finite-size effect of the self-diffusivity is caused from the latter, i.e., the self-interaction of the particle through the periodic boundaries. Based on this hydrodynamic model, Yeh and Hummer[47] derived eq for spherical particles and verified this correction for a single-component Lennard-Jones (LJ) fluid and pure water. Moultos et al.[48] showed that the YH correction is also applicable to molecules of varying size and shape, such as n-alkanes and glymes, provided that the system size exceeds 250 molecules.

In the studies by Spångberg et al.,[54] Kühne et al.,[55] and others,[56−58] the shear viscosities and self-diffusivities of water in the thermodynamic limit were calculated from eq . In these studies, several systems sizes were used to compute finite-size self-diffusivities, which were then fitted with a linear regression. To the best of our knowledge, the computation of shear viscosities from finite-size self-diffusivities has not been considered previously for multicomponent mixtures or highly viscous liquids. In this study, we use weighted least-squares linear regression analysis to develop a well-structured methodology for computing shear viscosities from finite-size effects of self-diffusivities. To allocate the available computational resources efficiently, a set of guidelines for choosing simulation parameters, such as the optimum number of system sizes and their size differences is provided. The application of the proposed method is verified for pure water, a large number of binary and ternary LJ systems, and the ionic liquid [Bmim][Tf2N].

This paper is organized in five sections. The proposed method is described in Section . In Section , details of MD simulations are briefly explained. The results of the MD simulations for pure water, binary and ternary LJ systems, and [Bmim][Tf2N] are discussed in Section , along with a set of guidelines for the efficient use of the proposed method. The conclusions are provided in Section .

Method

To develop a systematic method to compute shear viscosities from the finite-size effects of self-diffusivities, we write eq in a linear form, y = ax + b:where −ξkBT/6πL and DselfMD are the independent and dependent variables, respectively. The intercept of this line with the vertical axis (L → ∞) is the self-diffusivity in the thermodynamic limit, Dself∞. The inverse of the slope is the shear viscosity of the fluid, η.

To compute the shear viscosity, this method uses self-diffusivities of at least two system sizes. For each system size, the average self-diffusivity and its variance can be estimated from the mean (D̅) and sample variance (S2) of the self-diffusivities computed from several independent simulations:where Nsim, is the number of independent simulations for the jth system size and D indicates the self-diffusivity computed from the kth independent simulation for system size j. The parameters of interest (1/η and DselfMD) are then fitted to eq with weighted least-squares linear regression.[59] The weighted least-squares linear regression is briefly explained in the Supporting Information. The linear regression analysis requires the standard errors (SEs) of the average self-diffusivities for all system sizes. The inverse of squared SEs are used as the weighting factors for each data point. Since no prior knowledge of these SEs is available, the SE of the self-diffusivity of each system size can be estimated from the sample variance:According to the work of Pranami and Lamm,[41] multiple independent simulations are required to correctly compute the sample variance of the mean self-diffusivity. For a single MD simulation, the computed MSD depends on the initial configuration regardless of the simulation length. Weighted least-squares linear regression analysis yields both the averages and the variances of the parameters in eq . The self-diffusivity in the thermodynamic limit is a direct outcome of this analysis. The average shear viscosity is equal to the inverse of the average slope of the fitted line, a̅. If this slope has a statistical uncertainty of δa, the statistical uncertainty of the shear viscosity, δη, can be calculated from error propagation:

In a recent study, we showed that finite-size corrections to self-diffusivities of different species in a mixture are identical and equal to the YH correction (eq ).[49] This conclusion is based on a detailed study of 250 binary LJ systems with a wide range of LJ energy (ε), size (σ), and mass (m) ratios. The results show that the viscosity of a mixture can be predicted from the finite-size effects of the self-diffusivity of each species regardless of the mass or size ratios. To maximize the statistical information on a single simulation, we introduce a new quantity, the average self-diffusivity (Davg), which is the arithmetic mean of the self-diffusion coefficients of all species, weighted by their corresponding mole fractions. By using the definition of the MSD in multicomponent mixtures (eq ),[60−62] it can be shown that Davg is constructed from the self-diffusion of all molecules in the mixture:where n and N are the total number of species and molecules in the mixture, respectively, and x is the mole fraction of species i. Since the self-diffusivities of all species experience an identical finite-size effect,[49] the same YH correction (eq ) can be applied to Davg:By combining eqs and 10, the shear viscosity of a mixture can be obtained from Davg, similar to the approach used for pure liquids. Hereafter, the proposed method will be called the D-based method, and η denotes the corresponding computed shear viscosity.

The D-based method can be used to compute shear viscosities of highly viscous systems such as ionic liquids and deep eutectic solvents. The shear viscosity of these systems can be as large as several hundred cP at room temperature.[37,63−65] The length of MD simulations for computing self-diffusivities depends directly on how fast the constituent ions diffuse in the bulk liquid. In MD simulations, Fickian diffusion is observed at time scales in which a linear relation between the MSD and time (i.e., a slope of 1 on a log–log plot) is established.[66] However, this criterion does not ensure Gaussian diffusion, which corresponds to a Gaussian distribution for the displacement probability of ions.[67] Ionic liquids and deep eutectic solvents consist of highly associated pairs of ions and temporarily form cages.[68] While cage effects can also be present in simple fluids,[69] it plays an important role in determining the minimum length of an MD simulation for ionic systems.[68,70] As discussed in detail in the work of Casalegno et al.,[68] for short time scales, each ion fluctuates around a certain position in the cage. Due to dynamical heterogeneity,[67,68,71] a non-Gaussian distribution is observed for the diffusion probability of ions trapped in the cage. For longer time scales, ions jump from a cage to another. The displacement probability forms a Gaussian distribution corresponding to Gaussian diffusion.[68] According to Casalegno et al.,[68] for room-temperature ionic liquids, a rough estimate of the time scale corresponding to Gaussian diffusion can be made based on a minimum average displacement of 1.5 nm for all constituent ions. This criterion ensures that the simulation time is sufficiently long for all ions to break the local ion cages and diffuse in the bulk liquid.[70,72] This criterion can be used for room-temperature ionic liquids for performing long enough MD simulations to compute finite-size self-diffusivities with the D-based method.

Simulation Details

To validate the D-based method, three different systems are considered: binary and ternary LJ systems, pure water, and the ionic liquid [Bmim][Tf2N]. All MD simulations are carried out with LAMMPS (version 16, Feb 2016).[73] The order-n algorithm is used for an efficient sampling of time correlations for the calculation of self-diffusivities and shear viscosities.[66]

MD simulations of 26 ternary LJ systems were carried out. The results for 250 binary LJ systems are obtained from the Supporting Information of our previous work, ref (49). For the interaction between dissimilar LJ particles (i and j), we use the Lorentz–Berthelot mixing rules with a modification factor (k) to include nonideality in the systems:[34]where ε and σ are the parameters of the LJ potential. The simulation parameters of all binary and ternary systems are provided in Tables and 2, respectively. All parameters are reported in reduced units, where σ1 = σ = 1, ε1 = ε = 1, and m1 = m = 1 (mass) are the basis units.[34] All LJ interactions are truncated and shifted at a cutoff radius of 4σ. Simulations for two system sizes (500 and 4000 LJ particles) were performed at a reduced temperature of 0.65 and a reduced pressure of 0.05. The number densities of the binary[49] and ternary systems are in the range of 0.14–0.89 and 0.78–0.88, respectively. An integral time step of 0.001 is used, and the simulation lengths of the binary and ternary systems are 200 million and 100 million time steps, respectively.

Table 1

Specifications of 250 Binary LJ Systems at a Temperature of 0.65 and a Pressure of 0.05[49]a

specification	values
x1	0.1, 0.3, 0.5, 0.7, 0.9
ε2/ε1	1.0, 0.8, 0.6, 0.5
σ2/σ1	1.0, 1.2, 1.4, 1.6
m2/m1	(σ2/σ1)3
kij	0.05, 0.0, −0.3, −0.6

LJ species 1 has σ1 = σ = 1.0, ε1 = ε = 1.0, and mass = m1 = 1.0 in reduced units.[34]k is an adjustable parameter in the Lorentz–Berthelot mixing rules (eq ).

Table 2

Specifications of 26 Ternary LJ Systems at a Temperature of 0.65 and a Pressure of 0.05a

specification	values
k12	0.05, −0.3, −0.6
k13	0.05, −0.3, −0.6
k23	0.05, 0.0, −0.3, −0.6

The LJ interaction parameters of species 1, 2, and 3 are ε1 = ε = 1.0, ε2 = 0.8, and ε3 = 0.6, respectively. All species have the same size (σ1 = σ2 = σ3 = σ = 1.0) and mass (m1 = m2 = m3 = 1.0). All values are reported in reduced units.[34]x indicates the mole fraction of species i, where x1 = 0.4, x2 = 0.3, and x3 = 0.3.

For the molecular systems, the initial configurations were made in Packmol,[74] and the LAMMPS input files were created with VMD.[75] The three-site SPC/E model is used for water.[76] The force field parameters of [Bmim][Tf2N] are obtained from the work of Zhang et al.[77] The LJ interactions for water and the ionic liquid are truncated at 9 and 12 Å, respectively, with analytic tail corrections considered for energy and pressure. The particle–particle particle–mesh (PPPM) method with a relative precision of 10–6 is used for the long-range electrostatic interactions. The Verlet algorithm is used to integrate Newton’s equations of motion with a time step of 1 fs. MD simulations for water performed at 298 K and 1 atm. MD simulations of [Bmim][Tf2N] were performed at three temperatures of 300, 400, and 500 K and a pressure of 1 atm. The simulation times used for computing shear viscosities of water and [Bmim][Tf2N] are 50 and 200 ns, respectively. Self-diffusivities of water and [Bmim][Tf2N] are computed from independent simulations of 0.5 and 20 ns, respectively. Due to the slow relaxation of [Bmim][Tf2N] at 300 K, MD simulations of 50 and 450 ns are needed to obtain the self-diffusivities and shear viscosity at this low temperature, respectively. The choice of these simulation lengths is made so that a minimum average displacement of 1.5 nm is obtained for all systems (discussed thoroughly in Section ).

Results and Discussions

Pure Water

A set of simulations consisting of seven system sizes—250, 500, 1000, 2000, 4000, 8000, and 16 000 water molecules—was carried out. The average self-diffusivities were obtained from 100 independent simulations of 0.5 ns for each system size. In Figure , the average self-diffusivities of water are shown as a function of the system size. These finite-size self-diffusivities lie on the line fitted to eq :where N equals the number of water molecules in the simulation box. From the intercept with the vertical axis, the self-diffusivity of the SPC/E water model in the thermodynamic limit at 298 K and 1 atm is estimated to be (2.870 ± 0.004) × 10–9 m2 s–1 and the shear viscosity 0.708 ± 0.014 cP. The shear viscosity of water computed from the Einstein relation (eq ) is 0.694 ± 0.010 cP. These two values are in excellent agreement, considering the wide range of shear viscosities reported in the literature for the SPC/E water model at the same conditions: 0.68,[56] 0.71,[54] 0.729,[78,79] and 0.82 cP.[80] This agreement confirms the applicability of the D-based method for pure water.

Figure 1

Computed self-diffusion coefficients of the SPC/E water model at 298 K and 1 atm for seven system sizes: N = 250, 500, 1000, 2000, 4000, 8000, and 16 000 water molecules. A total of 100 independent simulations of 0.5 ns were performed for each system size. The dashed line is fitted with the weighted least-squares linear regression (eq ) to the YH equation (eq ). Error bars are smaller than the symbol sizes. The self-diffusivities are tabulated in Table S1 of the Supporting Information.

The extensive data set of finite-size water self-diffusivities provides a suitable estimation of the variances (S2) and the standard deviations (S) of self-diffusivities as a function of the system size. These standard deviations will be used in the next section for optimizing the D-based method. Figure shows the estimated standard deviations as a function of the number of molecules (N–1/3). Since no prior knowledge of the functional form is available, an initial guess for the functional form would be a power-law function (aN) with two fitting parameters:The exponent −0.31 and the linear arrangement of the data points in Figure suggest that the standard deviation can also be a linear function in N–1/3, which decreases the number of fitting parameters to one (aN–1/3):

Figure 2

Estimated standard deviations of the self-diffusion coefficients of the SPC/E water model at 298 K and 1 atm, computed for seven system sizes: N = 250, 500, 1000, 2000, 4000, 8000, and 16 000 water molecules. A total of 100 independent simulations of 0.5 ns were performed for each system size. The red and blue dashed lines are fits to a power-law (eq ) and a linear function (eq ), respectively.

In the next section, eq combined with eq will be used as a model for the normal distribution of finite-size self-diffusivities of water as a function of system size.

Optimization

In this section, a set of optimum simulation parameters for which the computed shear viscosities have a minimum statistical uncertainty is proposed. The following parameters are considered: the number of system sizes, the size difference between systems, and the allocation of computational resources to each system size. The size difference between system i and system j is normalized by the size of the smallest system (system 1): (N – N)/N1. According to the work of Moultos et al.,[48] the smallest system should contain at least 250 molecules. This criterion ensures that the YH correction (eq ) provides an accurate prediction for the finite-size effects of self-diffusivities. As a constraint on the optimization problem, the total computational resources are fixed. The computational resources scale linearly with the number of independent simulations (Nsim,) and polynomially with the number of molecules (N) in the simulation box, depending on computer hardware, the employed computational methods, and the scalability of MD simulations.[34] Thus, the ratio, α, between the computational resources allocated to system i and system j iswhere γ indicates the scalability of MD simulations. Two values of γ are considered here for the scalability of MD simulations: γ = 1 and 2. For γ = 1, computational requirements grow linearly with the system size. For γ = 2, the growth in computational requirements is quadratic and thus faster than the linear growth. Current state-of-the-art MD packages[73,81,82] have a good computational scalability, γ, close to 1.

The simulation results presented for the SPC/E water model can be used as a basis for finding the optimum combination of the simulation parameters for the D-based method. eqs and 14 are used to model the average and standard deviation of finite-size self-diffusivity of water as a function of system size. This model enables us to predict the self-diffusivity of water for a hypothetical MD simulation with a specified number of molecules. This can be achieved by generating a random number from a normal distribution with a mean self-diffusivity and a standard deviation determined by eqs and 14. For a set of simulations with a number of system sizes and independent simulations per system size, a set of finite-size self-diffusivities is constructed and the corresponding shear viscosity is calculated from the D-based method. Since this shear viscosity depends on the set of finite-size self-diffusivities, the procedure should be repeated for many times and all shear viscosities are recorded in a histogram. For a specific set of simulation parameters, the data stored in the histogram yield an estimate for the variance (S2) and standard deviation (S) of the shear viscosity. The aim of this optimization procedure is to find the simulation parameters that minimize the standard deviation (S) of the shear viscosity.

The first scenario considered here is the optimization for two system sizes. The simulation parameters studied are (i) the normalized size difference ((N2 – N1)/N1), and (ii) the ratio between the computational resources allocated to system 1 and 2 (α). The objective function is the estimated standard deviation of shear viscosities normalized by Smin, the global minimum estimated standard deviation of the shear viscosity for a specified value of γ and all values of α and (N2 – N1)/N1. In Figure , the normalized estimated standard deviations (S/Smin) are shown for several values of α, ranging from 0.2 to 5.0, and the two values of γ. As can be seen in Figure , a range of α between 1 and 2 yields the smallest value of S/Smin. This means that the number of independent simulations for each system size should be distributed in such a way that 50%–70% of the computational resources is allocated to the large system. Furthermore, it can be observed that the large system should be at least 4 times the size of the small system. Depending on the scalability of MD simulations, the optimum normalized size difference ((N2 – N1)/N1) ranges from 3 (γ = 2) to 40 (γ = 1).

Figure 3

Normalized estimated standard deviation (S/Smin) of the shear viscosity as a function of the normalized size difference between two systems ((N2 – N1)/N1). The total amount of computational resources is fixed. Different colors indicate various ratios α of the computational resources allocated to large and small system sizes (eq ): 0.2 (black), 0.5 (magenta), 1.0 (green), 2.0 (red), and 5.0 (blue). Two types of scalability for MD simulations are considered: (a) high scalability (γ = 1) and (b) low scalability (γ = 2).

A similar investigation can be carried out for three system sizes. Here, α is set to 1, indicating that the computational resources are equally distributed between the three system sizes. In Figure , normalized estimated standard deviations (S/Smin) are shown as a function of the normalized size difference between the small and medium systems ((N2 – N1)/N1), and the medium and large systems ((N3 – N2)/N1). Figure a (γ = 1) shows that a minimum cost function (S/Smin) is achieved for (N2 – N1)/N1 > 10. This suggests that, while the choice of the size difference between the small and medium systems is important, the size difference between the medium and large systems does not play a significant role. For very small values of the normalized size difference between the medium and large systems ((N3 – N2)/N1 → 0), it can be deduced that the optimization problem reduces to a problem of two system sizes, which has already been discussed. In Figure b, this is also observed for γ = 2. The optimum condition is achieved at (N2 – N1)/N1 = 3, which is in agreement with what is shown in Figure b for two system sizes. This suggests that the optimum condition observed for three system sizes can also be achieved with two systems.

Figure 4

Normalized estimated standard deviations (S/Smin) of the shear viscosities as a function of the normalized size difference between small and medium systems ((N2 – N1)/N1), and between medium and large system sizes ((N3 – N2)/N1). A fixed amount of computational resources is equally distributed between the three system sizes. Two types of scalability for the MD simulations are considered: (a) high scalability (γ = 1) and (b) low scalability (γ = 2).

Based on the optimization results, the use of either two or three system sizes yields optimum shear viscosities. In the case of a limited amount of computational resources, adding more system sizes leads to a smaller number of independent simulations per system size. As a consequence, the limited number of independent simulations leads to poor sampling that may not yield an accurate average and standard deviation for self-diffusivities. This adversely affects the accuracy of average finite-size self-diffusivities and consequently the computed shear viscosity. Hence, the use of more than two system sizes is not justified. The choice of the optimum conditions may vary depending on the MD software[73,81,82] as well as the algorithms used in the software (e.g., handling of long-range electrostatic interactions, details of the neighbor lists[34]). These specifications determine the optimum size of the system and number of independent simulations according to the results shown in Figure , and consequently affect the computational requirements of the D-based method.

Lennard-Jones Systems

To examine the accuracy of the D-based method for multicomponent fluid mixtures, the shear viscosities of 250 binary[49] and 26 ternary LJ systems were computed. The comparison between the D-based shear viscosities and those computed from the Einstein relation (ηEMD, eq ) are shown in Figure . By varying the characteristics of the studied LJ systems as mentioned in Tables and 2, a wide range of shear viscosities is covered. While all data points show a good agreement between the two methods, the smallest deviation from the diagonal dashed lines is observed for the quantity Davg, for both binary and ternary systems. This is expected since Davg is constructed from the self-diffusivities of all species present in the mixture. Therefore, for multicomponent mixtures, Davg should be used for calculating the shear viscosity rather than the self-diffusivity of only a single species.

Figure 5

Comparison between the shear viscosities of (a) 250 binary and (b) 26 ternary LJ systems computed from the Einstein relation (ηEMD, eq ) and the D-based method (η, eq ) at a reduced temperature of 0.65 and a reduced pressure of 0.05. The D-based shear viscosities are computed from the self-diffusion coefficients of species 1 (blue circles), species 2 (green squares), and species 3 (cyan diamonds; only for ternary systems) as well as the average self-diffusivity (eq , red crosses). Error bars are omitted for clarity. All computed shear viscosities and their statistical uncertainties are listed in the Supporting Information, Tables S4 and S7.

In Figure , the maximum deviation from the diagonal line, which represent perfect agreement between the D-based method and the Einstein relation, are observed for binary LJ systems with very low densities (inset of Figure a). As discussed in our previous study,[49] these outliers have very dissimilar size (σ2/σ1 = 1.4 or 1.6) and interaction (ε2/ε1 = 0.5 or 0.6) parameters. According to the work of Heyes et al.[83] for a hard-sphere fluid, the exponent −1/3 in N–1/3 (i.e., L–1) in eq is valid only for a range of packing fractions. Since the validity of the YH correction have been shown for many real molecular liquids,[47−50] it is expected that the range of packing fractions for which the YH correction holds correspond to a liquid phase. From Figure , a similar observation to the work of Heyes et al.[83] can be made for LJ systems. In this figure, the normalized difference between the shear viscosities computed from the D-based method (using Davg) and the Einstein relation is shown as a function of (a) the shear viscosity and (b) the density for all binary and ternary LJ systems. For LJ systems at high densities, a good agreement between the two methods is observed. For the systems with low densities/viscosities, shown also in the inset of Figure , substantial deviations (up to 30%) between the D-based method and the Einstein relation are observed. Figure shows that as the density of a system decreases and thus hydrodynamic interactions become weaker, the scaling proposed by Yeh and Hummer does not hold anymore. At low densities, the finite-size effects of self-diffusivities do not scale as N–1/3, and the exponent varies according to the packing fraction of the fluid (e.g., see the work of Heyes et al.[83]). As discussed briefly in the Supporting Information, the outliers shown in Figures and 6 correspond to conditions close to the critical point. This suggests that neither the YH correction nor the D-based method is applicable at conditions close to the critical point. Nevertheless, our findings show that the D-based method is able to predict accurate shear viscosities for all real molecular systems in a liquid phase far from the critical point.

Figure 6

Normalized absolute difference between shear viscosities computed from the Einstein relation and the D-based method as a function of (a) the shear viscosity and (b) the density. Data are shown for binary (blue diamonds) and ternary (green squares) LJ systems at a temperature of 0.65 and a pressure of 0.05. The D-based method shear viscosities are computed from average self-diffusivities (Davg).

[Bmim][Tf2N]

As a representative test case, shear viscosities of the ionic liquid [Bmim][Tf2N] are calculated from the D-based method and the Einstein relation (eq ) at a pressure of 1 atm and three temperatures: 300, 400, and 500 K. MD simulations of 20 ns were performed for computing self-diffusivities of [Bmim][Tf2N] at 400 and 500 K. To fulfill the minimum average displacement criterion of 1.5 nm (explained in Section ), the MD simulations at 300 K were carried out for 50 ns. According to the set of guidelines proposed for the D-based method, 40 and 8 independent MD simulations were performed for two systems of 150 and 1200 pairs of ions, respectively.

In Figure , the computed self-diffusivities for both [Bmim] and [Tf2N] ions as well as Davg are shown. It can be seen that the finite-size self-diffusivities of [Bmim] and [Tf2N] greatly differ from each other, but the finite-size effects remain equal for both ions and Davg. Therefore, shear viscosities computed from the slopes of the lines connecting finite-size self-diffusivities are equal within the error bars. The computed shear viscosities from the D-based method and the Einstein relation, and their 95% confidence intervals are reported in Table S11 in the Supporting Information. Figure shows that the use of Davg for the D-based method is not limited to molecular mixtures and that this quantity can be used for ionic systems as well.

Figure 7

Computed self-diffusion coefficients of [Bmim][Tf2N] at 300, 400, and 500 K and 1 atm. In all, 40 and 8 independent simulations were performed for two system sizes of 150 and 1200 ion pairs, respectively. Self-diffusivities are shown for [Bmim] (blue circles), [Tf2N] (green diamonds), and Davg (red squares). The slope of the line connecting each two points yields the shear viscosity. The simulation length at 400 and 500 K is 20 ns. Simulations at 300 K were performed for 50 ns. All self-diffusivities are listed in Tables S8–S10 of the Supporting Information. Error bars are smaller than the symbol sizes.

Figure shows the computed shear viscosities of the ionic liquid as a function of temperature. The lines are the fits to the Vogel equation:[84]where A, B, and C are the coefficients of this equation. For all temperatures, a good agreement is observed between the D-based method, the Einstein relation (eq ), and the estimates of shear viscosities from the work of Zhang et al.[77] This agreement confirms the applicability of the D-based method to complex and highly viscous mixtures of nonspherical molecules/ions. To compute the shear viscosity of any ionic liquid or deep eutectic solvent, the set of guidelines for the D-based method along with the criterion on the minimum average displacement can be used to specify an optimum set of MD simulations.

Figure 8

Shear viscosity of [Bmim][Tf2N] as a function of temperature at 1 atm, computed from the D-based method (blue circles, eq ) and the Einstein relation (red squares, eq ). The lines are fits to the Vogel equation[84] (eq ) for the D-based method (blue dashed, eq ), Einstein relation (red dashed), and Green–Kubo relation (green solid; data extracted from the work of Zhang et al.[77]). The shear viscosities and coefficients of the Vogel equation are provided in Tables S11 and S12 of the Supporting Information.

Conclusion

A systematic methodology, called the D-based method, is proposed for accurately computing the shear viscosity of a liquid from the finite-size effects of self-diffusivities. The computational requirements of this method and the statistical uncertainty of the computed viscosity are comparable to the conventional methods, e.g., the Einstein relation. By performing weighted least-squares linear regression analysis, the shear viscosity can be computed from the slope of a line fitted to computed finite-size self-diffusivities. To obtain accurate shear viscosities at a minimum computational requirement, a set of guidelines for this method was proposed. The optimum number of system sizes is two, and depending on the available computational resources and the scalability of the MD simulations, the large system size should be 4–40 times the small system size. The number of independent simulations per system size should be assigned in such a way that 50%–70% of the computational resources is allocated to the MD simulations of the large system. For multicomponent mixtures, the D-based method performs best when the average self-diffusivity of all species (Davg) is used instead of using the self-diffusivity of a single species. The D-based method was verified for pure water, a large number of binary and ternary Lennard-Jones systems, and an ionic liquid ([Bmim][Tf2N]). The results of the D-based method were in good agreement with those obtained from the Green–Kubo and Einstein relations for all molecular systems. These results suggest that the D-based method can be a potential method for computing shear viscosities of highly viscous liquids and multicomponent mixtures.