System-size corrections for self-diffusion coefficients calculated from molecular dynamics simulations: The case of CO₂, n-alkanes, and poly(ethylene glycol) dimethyl ethers

Self-diffusion coefficient is a fundamental transport property, essential for understanding the structure and thermophysical behavior of a fluid, as well as for the design of industrial processes.^1–3 A comprehensive review of the available experimental measurements of self-diffusion coefficients for a wide range of substances (from noble gases to complex organic molecules), along with a detailed discussion of the different experimental measurement techniques, is provided by Suárez-Iglesias et al.⁴ 

Employing molecular dynamics (MD) simulations for the calculation of self-diffusion coefficients is an appealing alternative to experiments,⁵ which very often include high cost or technical difficulties when high pressures and/or temperatures are involved. The first such simulation study calculating self-diffusion coefficients was performed approximately fifty years ago for the hard sphere fluid by Alder and Wainwright.⁶ Liu et al.,⁷ Vaz et al.,⁸ Ohtori and Ishii,^9,10 and Pranami and Lamm¹¹ reported detailed MD studies for the self-diffusion coefficients of model fluids (e.g., hard sphere, Lennard-Jones, and square-well). Additionally, numerous MD studies examining real fluids, such as CH₄, CO₂, and H₂O, can be found in the recent literature.^12–18 

In typical MD simulations (using up to a few thousand molecules under periodic boundary conditions), long-range interactions can lead to significant system-size effects and thus to systematic errors in the calculation of self-diffusion coefficients, if no corrections are applied.^19–22 The approach of Yeh and Hummer (YH)²¹ has been used in the literature^16,23–25 for correcting results obtained with finite systems to estimate infinite-system (thermodynamic limit) extrapolated values. According to YH, valid for cubic boxes, the following correction, originating from hydrodynamic interactions, needs to be applied in order to account for system-size effects:

where D_MD is the self-diffusion coefficient obtained from the MD simulations, D_∞ is the self-diffusion coefficient corrected for the system-size effects, L is the length of the cubic simulation box, k_B is the Boltzmann constant, T is the absolute temperature, η is the shear viscosity, and ξ ≈ 2.837 298 is a dimensionless constant determined by an Ewald-like summation of a periodic lattice.^13,21,26 Eq. (1) suggests that D_MD ∝ N^−1/3, since L ∝ N^1/3, where N is the number of molecules.

Recently, however, Heyes et al.¹² performed a series of simulations for the hard sphere fluid over a broad range of conditions and found that Eq. (1) is, in many cases, not adequate to represent the system size dependence of D_MD. Heyes et al. showed that a more complex density dependency is needed to describe data over a wide range of conditions. Consequently, they proposed the following empirical expression for the self-diffusion coefficient:

where A and α are calculated through polynomial expressions which are functions of the density. When α = 1/3, Eq. (2) reduces to a dependence on N similar to the YH expression. Heyes et al. observed that for the hard sphere fluid, α approaches 1 at low densities, goes through a minimum (of α ≈ 1/3) at intermediate densities, and increases to 1 (or higher) for high densities.

Despite the fact that system size corrections must be applied to the calculation of self-diffusion coefficients for all fluids, relatively few MD studies take them into account. Popular alternatives are the use of relative large system sizes (1000 molecules or more) or bigger cutoff radii of the potential.^17,22,27 Such alternatives are shown to mitigate the system size dependency at ambient conditions in some cases (i.e., H₂O).^16,22 In the current study, we examine the validity of the YH approach for a number of different fluids. Of particular interest are real fluids belonging to three characteristic groups of fluids that are all commonly used as solvents. The fluids are carbon dioxide, n-alkanes (methane, propane, n-hexane, and n-hexadecane), and poly(ethylene glycol) dimethyl ethers (glymes in short) CH₃O–(CH₂CH₂O)_n–CH₃ with n = 1, 2, 3, and 4 (G1, G2, G3, and G4, respectively).

Equilibrium MD simulations for carbon dioxide and n-alkanes (methane, propane, n-hexane, and n-hexadecane) were conducted with the open source package GROMACS (version 4.6.5),²⁸ whereas the four different glymes were simulated using LAMMPS (version 18 February 11).²⁹ Intermolecular and intramolecular interactions for all systems studied were described using parameters taken from the TraPPE force field.^30–32 Bond stretching parameters for glymes were adopted from the general Amber force field (GAFF)³³ so that the molecules were fully flexible. The 1–4 interactions for n-alkanes were ignored, while for the case of glymes, they were scaled by a factor of 0.5, which was found to generate good density and dynamic properties.³⁴ A cutoff distance of 14 Å was applied for both Coulombic and van der Waals interactions. The long-range electrostatic interactions were calculated using the particle mesh Ewald (PME)^35,36 in GROMACS and the particle-particle particle-mesh (PPPM) method³⁷ in LAMMPS.

All simulations of carbon dioxide and n-alkanes were performed in the isothermal-isobaric (NPT) ensemble with each state point being calculated by averaging the results of 10 independent production runs of 2–10 ns long, each one starting from a different initial configuration. The simulation temperature and pressure were maintained constant using the Nosé-Hoover^38,39 thermostat and Parrinello-Rahman⁴⁰ barostat, respectively, with coupling constants equal to 1 ps.

Simulation boxes containing 125, 250, 500, 1000, 2000, and 4000 glyme molecules were constructed for each glyme by putting the molecules randomly in a cubic box using the package Packmol.^41,42 Energy minimization was then carried out to remove overlaps. The minimization was followed by a 2 ns equilibration run in the isothermal-isobaric (NPT) ensemble at a pressure of 1 atm (except for G1, where the pressure was set at 30 atm to maintain a liquid phase) and a temperature of 400 K. The densities obtained from these simulations were then used to carry out three independent production runs at the same temperature in the canonical (NVT) ensemble. Each NVT production run was 10 ns long; the first 1 ns of which was considered equilibration and the last 9 ns was used to calculate the self-diffusion coefficient. The Nosé-Hoover thermostat^38,39 and the extended Lagrangian approach⁴³ were used to control the temperature and pressure, respectively.

Periodic boundary conditions were imposed in all simulations, and a time step of 1 fs was used. Tail corrections were applied for energy and pressure. For each system, the self-diffusion coefficients were calculated using the Einstein relation,

where r_i(t) is the position of the center of mass of species i at time t, and the angle brackets indicate an ensemble average. To improve the statistics, different time origins were taken along the simulated trajectory in the calculation of the mean square displacements $rit−ri02$⁠. The calculated mean square displacements as a function of time for the glymes, as the most complex fluid studied, together with the fit to the Einstein relation are shown in Fig. S1 of the supplementary material to demonstrate the linearity. A representative snapshot of the simulation box containing 125 and 4000 G4 molecules, respectively, is provided in Fig. S2 of the supplementary material.

The shear viscosity of glymes was calculated using the time-decomposition method.⁴⁴ For each system, 30 independent 5 ns NVT trajectories were generated. For each trajectory, the stress tensor was saved every 5 fs. The first 1 ns of each trajectory was considered equilibration and the remaining 4 ns was used to calculate the shear viscosity following the Green-Kubo relation,

where V is the system volume, and τ_αβ denotes the αβ-component of the stress tensor. Similar to the case of the self-diffusion coefficient, different time origins were used in the calculation of the correlation function to decrease the uncertainty. The average running integral over 30 trajectories was fit to a double exponential function of the form

with a weighting function derived from the corresponding time-dependent standard deviation of the running integrals

In the above equations, N is the number of independent trajectories, and A, α, τ₁, and τ₂ are fitting parameters. α in Eq. (5) should not be confused with α in Eq. (2). The value of the fitting results at the long time limit was taken as the calculated viscosity of each system. Viscosity at each state point for the CO₂ and n-alkanes was calculated by averaging 5 independent 10 ns simulation results using the Green-Kubo method (Eq. (4)). These computed viscosities were then used in the YH correction.

The self-diffusion coefficients of carbon dioxide and methane were calculated for systems ranging from 250 to 2000 molecules at 323.15 K and 200 bars and 150.15 K and 500 bars, respectively. The results are shown in Fig. 1 and in tabular form in Table S1 of the supplementary material. Available experimental data are also provided in Table S3 of the supplementary material for comparison. The calculated viscosities agree with the experimental results very good. The self-diffusion coefficients for both fluids scale linearly with 1/N^1/3. In Fig. 1, the extrapolated value of the self-diffusion coefficient for infinite system size, D_∞,MD, is shown as a blue dashed line. The values are also reported in Table I. Using Eq. (1), the YH correction was applied, and the corrected results are shown as red diamonds, with a mean value represented as a red dashed line. The viscosity used in Eq. (1) was calculated from MD simulations for both fluids. The calculated viscosity values and available experimental data are given in the supplementary material. As indicated clearly in Fig. 1, the corrected values, D_{∞,Y H}, are in excellent agreement with the extrapolated diffusivity D_∞,MD, except for the smallest system size (250 molecules) for the case of carbon dioxide. This difference can be attributed to the poor statistics in the calculation of self-diffusion coefficient due to the small number of molecules. We have to notice here that the actual magnitude of the correction is significant. For the case of carbon dioxide, it ranges from approximately 15% (for the smallest size of 250 molecules) to approximately 8% (for the largest size considered, i.e., 2000 molecules), while for methane the range is 6%–11%. This finding suggests that even the use of relatively large system sizes in the calculation of self-diffusion coefficient of fluids cannot guarantee an accurate prediction and in most cases a correction, like the form suggested by YH, must be applied. It should be noted here that the densities obtained from the different system sizes were practically identical. This finding is in line with the work by Yeh and Hummer²¹ and other studies in the literature.^16,23–25

As mentioned in the Introduction, for the hard sphere fluid, Heyes et al. showed that the value of the exponent α in Eq. (2) can have different values depending on the density.¹² In the current work, in addition to α = 1/3 as used in the YH method,²¹ the calculated self-diffusion coefficients for CO₂ and methane were also fitted to a linear function of the form of Eq. (2) with α = 1/2 and α = 2/3, respectively. Table I shows the calculated values of D_∞,MD, for the three values of α. The fitting results are provided in Figs. S3 and S4 of the supplementary material. The D_∞,MD values of CO₂ and methane using the three different exponents are very close to each other (deviation approximately 4%), suggesting that, for the densities studied in the current work for CO₂ and methane, the YH method corrects the system size effect reasonably well.

In Fig. 2, the self-diffusion coefficient of propane (Fig. 2 top), n-hexane (Fig. 2 middle), and n-hexadecane (Fig. 2 bottom) is presented as a function of 1/N^1/3 for system sizes ranging from 500 to 4000 molecules for propane and n-hexane and from 250 to 2000 for n-hexadecane. Numerical values and available experimental data are listed in Tables S1 and S3 of the supplementary material. Similar to the cases of carbon dioxide and methane, the diffusivities obtained by the different system sizes scale linearly with 1/N^1/3. For the case of propane, the MD results for the different system sizes deviate from D_∞,MD by approximately 2.4%–4.9% (4000 to 500 molecules). The corrected values D_{∞,Y H}, using the YH formula (Eq. (1)), were found to deviate from D_∞,MD by approximately 0.5%. This value is, in most cases, within the statistical uncertainty of the MD simulations. For the case of n-hexane (Fig. 2 middle panel), the magnitude of the YH correction ranges from 5.3% to 9.9%. YH corrected values D_∞,YH overshoot D_∞,MD by 1.3%. The correction in n-hexadecane diffusion coefficients for the different system sizes is in the range of 9%–18% (2000 to 250 molecules), while the value D_{∞,Y H} is lower compared to D_∞,MD by 1.7%. From these findings, it becomes obvious that even the use of 2000 molecules, which is a relatively large system size, leads in self-diffusion coefficients that significantly deviate (3%–9%) from the ones extrapolated to infinite system size. For some cases, like n-hexane and n-hexadecane shown in Fig. 2, even the use of a very large number of molecules (i.e., 10 000 or 1/N^1/3 = 0.046) results in errors in the range of 2%–5%. These error values, most probably, will be higher than the uncertainties of the calculations, given that as the number of molecules increases, the statistical uncertainties in the MD simulation decrease substantially. At the same time, the computational cost becomes very high.

As for the case of CO₂ and methane, Table I includes the D_∞,MD values of propane, n-hexane, and n-hexadecane using the three different values of α and the fitting results are provided in Figs. S5–S7 of the supplementary material. The largest deviation between the results for different α values was found for n-hexadecane (6.4%) and the smallest for propane (1.8%), while for n-hexane, the deviation is approximately 4%.

As described in Sec. II, the equilibrium box size was determined for each glyme through separate NPT ensemble simulations. For each glyme, the deviation in the calculated density for systems with different numbers of molecules varied by no more than 0.2%, indicating that density is independent of system size. Therefore any differences in the self-diffusivities or viscosities as a function of system size are not attributable to density differences. In addition, the calculated density of each glyme agrees with available experimental results^45,46 very well, with a maximum deviation of 1.6%. All density results and experimental values are provided in Tables S2 and S3, respectively, of the supplementary material.

Using the time-decomposition method, the viscosities of the glymes were calculated for each system size. The results are shown in Fig. 3, and values derived from available experimental data are provided in the supplementary material. The calculated viscosities show some differences at different system sizes, but there is no clear system size dependence, consistent with previous studies.^21,47,48 Additionally, a systematic discussion on the effect of system size on MD calculated viscosities is presented in the recent study by Bernardi et al.⁴⁹ 

Based on NVT ensemble simulations, the center of mass self-diffusion coefficients was calculated using the Einstein relation (Eq. (3)). The results are shown in Fig. 4, with the numerical values listed in Table S2 of the supplementary material. The calculated self-diffusion coefficients were found to be larger in the boxes having more glyme molecules. Similar to the case of carbon dioxide, methane, and the n-alkanes, a linear dependence on 1/N^1/3 was observed for each glyme except for the smallest boxes (125 molecules) in the case of G2 and G4. Similar behavior was previously reported by Gabl and co-workers.⁴⁸ Relatively large uncertainty in the calculated self-diffusion coefficients in the small boxes was also observed for G1 and G3. Considering the limited number of molecules in these systems and the finite time scale of the simulations, such large uncertainty and deviation from the linear behavior can be the result of insufficient statistics. Therefore we did not include the results with 125 molecules in the fitting. The self-diffusion coefficients calculated for the five largest boxes were fit to a linear function and extrapolated to infinite box size to obtain an estimate of D_∞,MD. The D_∞,MD values obtained this way are shown in Fig. 4 as blue dashed lines. For each glyme, the calculated self-diffusion coefficient for the largest box (4000 molecules) was found to be 4%–7% lower than the corresponding D_∞,MD, whereas that calculated for the smallest box (250 molecules) was found to be up to 17% lower. The finite size effect of the calculated self-diffusion coefficient was also corrected using Eq. (1) for each system, and the resultant D_∞,YH values are shown in Fig. 4. The average viscosity calculated over the five largest box sizes was used in the calculation for each glyme. Except for the smallest box of G2 and G4, the calculated D_{∞,Y H} values all agree with one another within the statistical uncertainties and are independent of the system size. This observation is consistent with that reported previously for water, Lennard-Jones, and hard-sphere fluid²¹ as well as the results presented here for carbon dioxide, methane, and the n-alkanes.

The average of D_{∞,Y H} over the five largest boxes was calculated for each glyme, and the results are shown as red dashed lines in Fig. 4. In each case, the calculated self-diffusion coefficient was found to be 3%–6% lower than the average D_{∞,Y H} for the largest box (4000 molecules), whereas that calculated from the smallest box (250 molecules) was found to be around 12%–14% lower. Using polarizable force fields, Borodin studied the finite size effect on self-diffusion coefficients of several ionic liquid systems and reported that the corrections are typically 10%–20%.²¹ A correction of 15% was also reported for an ionic liquid 1-ethyl-3-methylimidazolium triflate by Gabl and co-workers.⁴⁸ These results fall in the same range found for the small boxes studied in the current work. In Table I, the self-diffusion coefficient values D_∞,MD using the three different α values are reported. As discussed previously for CO₂ and alkanes, the values obtained by different α values are almost identical (deviation of less than 4%). Figures of diffusion coefficient as a function of 1/N^1/3, 1/N^1/2, and 1/N^2/3 are given in Figs. S8–S11 of the supplementary material.

In a previous study, the corrected self-diffusion coefficients using Eq. (1) were found to agree with extrapolated MD results almost perfectly for the Lennard-Jones and hard-sphere fluids.²¹ A slightly over-correction was observed for water, which, according to the authors, was caused by the long-range electrostatic interaction.²¹ In the current work, relative to D_∞,MD, the finite size effect was found to be slightly over-corrected by Eq. (1) for G1 and under-corrected for the other three glymes. However, for all the glymes studied here, the deviations between D_{∞,Y H} and D_∞,MD are smaller than 2% and are within the statistical errors of the simulations. These results suggest that the correction (Eq. (1)) proposed by Yeh and Hummer²¹ works well even for molecular fluids such as glymes. A more thorough study is needed to describe the precise mechanisms behind the system size dependencies for the self-diffusion coefficient of real fluids and mixtures.

This study reveals that for MD simulation of fluids, the use of “large” systems of several thousand molecules is not always adequate for eliminating the dependence of self-diffusion coefficients on system size, and that explicit corrections must be applied to account for systematic errors introduced by finite size effects. To this extent, YH corrections are shown to be adequate for all systems examined, namely, carbon dioxide, methane, propane, n-hexane, n-hexadecane, and various poly(ethylene glycol) dimethyl ethers (G1–G4). The magnitude of the correction is significant, with the smallest systems examined (N = 250) giving for some cases a self-diffusion coefficient approximately 15% lower than the infinite system-size extrapolated value. Even a fairly large system of N = 2000 shows deviations of 3%–8% from the infinite size extrapolated values, indicating that proper corrections must be applied. These results suggest that finite size corrections to computed self-diffusivities must be used in order to obtain accurate results.

See supplementary material for tables and figures with additional results from MD simulations and experiments.

This publication was made possible by NPRP Award No. NPRP 6-1157-2-471 from the Qatar National Research Fund (a member of The Qatar Foundation) to the Texas A&M University at Qatar team. The statements made herein are solely the responsibility of the authors. Y.Z. and E.M. are supported by the U.S. Department of Energy, Basic Energy Science, Joint Center for Energy Storage Research under Contract No. DE-AC02-06CH11357. We thank Professor A. Z. Panagiotopoulos for very useful discussions. Computational resources were provided by the High Performance Computing Center of Texas A&M University at Qatar and the Center for Research Computing (CRC) at the University of Notre Dame.