Generalized Form for Finite-Size Corrections in Mutual Diffusion Coefficients of Multicomponent Mixtures Obtained from Equilibrium Molecular Dynamics Simulation.

The system-size dependence of computed mutual diffusion coefficients of multicomponent mixtures is investigated, and a generalized correction term is derived. The generalized finite-size correction term was validated for the ternary molecular mixture chloroform/acetone/methanol as well as 28 ternary LJ systems. It is shown that only the diagonal elements of the Fick matrix show system-size dependency. The finite-size effects of these elements can be corrected by adding the term derived by Yeh and Hummer (J. Phys. Chem. B 2004, 108, 15873-15879). By performing an eigenvalue analysis of the finite-size effects of the matrix of Fick diffusivities we show that the eigenvector matrix of Fick diffusivities does not depend on the size of the simulation box. Only eigenvalues, which describe the speed of diffusion, depend on the size of the system. An analytic relation for finite-size effects of the matrix of Maxwell-Stefan diffusivities was developed. All Maxwell-Stefan diffusivities depend on the system size, and the required correction depends on the matrix of thermodynamic factors.

Introduction

The past few decades, Equilibrium Molecular Dynamics (EMD) simulation has emerged as a powerful tool for computing diffusion coefficients of pure components and multicomponent mixtures.[1−18] Typically, system sizes in the order of hundreds to a few thousands molecules are used, combined with periodic boundary conditions.[19,20] As shown in the pioneering work by Dünweg and Kremer almost 30 years ago,[21] the choice of the system size strongly affects computed self-diffusivities, which scale linearly with the inverse of the simulation box size. In early 2000s, Yeh and Hummer[22] (YH) derived an analytic hydrodynamic correction which should be added to self-diffusivities computed from EMD, to obtain diffusivities at the thermodynamic limit, i.e., the quantity measured in experiments. The practically more relevant mass transport due to concentration gradients depends on mutual diffusion coefficients (i.e., Fick and Maxwell–Stefan).[23−25] Very recently, Jamali et al.[26] proposed an empirical correction that should be applied to computed Maxwell–Stefan (MS) diffusivities in binary mixtures. To the best of our knowledge, it is still unknown if and how finite-size corrections should be applied to mutual diffusivities in multicomponent systems. In this study, we show that such corrections should indeed be applied, and we derive a general correction term for mutual diffusion coefficients of multicomponent mixtures computed from EMD. In this way, multicomponent diffusivities can be computed in a reliable way, allowing for quantitative comparison with experiments. The derivation of the new generalized correction presented in the present study proves the validity of the empirical correction term proposed by Jamali et al.[26] for binary mixtures and provides fundamental understanding of the underlying mechanisms of system size effects in multicomponent mutual diffusion.

The manuscript is organized as follows: Some essential theoretical definitions are provided in the following paragraphs. The generalized correction term for the mutual diffusivities in multicomponent mixtures is derived in the Results and Discussion section, where also the verification of the new correction is shown for 28 ternary Lennard-Jones (LJ) mixtures and for a ternary molecular mixture. In the Conclusions section, the main findings of this study are summarized.

Background: Mutual Diffusivity and Finite-Size Effects

Fick (mutual) diffusion coefficients are commonly used to describe transport diffusion in a mixture.[23,24,27] According to Fick’s law, the molar flux (J) of species i in an n-component mixture is proportional to the concentration gradients of (n – 1) constituent species (c):[23,28]in which D denotes Fick diffusion coefficients. For an n-component mixture, the matrix of Fick diffusivities consists of (n – 1)2 elements. The empirical nature of Fick’s law has been overcome by the MS formulation of diffusion, according to which the driving force for diffusion is the gradient in chemical potential (μ), which is in balance with frictional forces:[23,24,28,29]where x is the mole fraction of species i, c is the total concentration, R is the universal gas constant, and T is the temperature of the system. The friction force is due to the difference in the fluxes (velocities) of species i and j, (i.e, xJ – xJ). The inverse of the MS diffusion coefficient (i.e., ) can be considered as a friction coefficient. For an n-component mixture, n·(n – 1)/2 diffusion coefficients are defined. In a homogeneous mixture, D and Đ describe the same physical phenomenon. These diffusivities are related via the so-called matrix of thermodynamic factors, [Γ]:[23,24,28]where [DFick] is the matrix of Fick diffusivities. [Δ] is the symmetric phenomenological diffusion coefficient matrix which is related to MS diffusivities according to[11,28−30]where i and j take values from 1 to n – 1 and i ≠ j. [Γ] is an asymmetric matrix, whose elements can be computed from[28,31−33]where δ is the Kronecker delta and the subscript “T,p,Σ” denotes constant temperature, pressure, and mole fractions of all species, except for the nth species.[28,29,33] More information about the computation of [Δ] and [Γ] from EMD simulations can be found in the literature.[28−30,34−38]

The Brownian motion of molecules in a pure or a multicomponent fluid mixture can be described by the self- (or tracer) diffusivity.[39] In EMD, self-diffusivities are computed from ensemble averages of mean-squared displacements of individual molecules.[28,40] According to the studies of Dünweg and Kremer,[21] and Yeh and Hummer,[22] computed self-diffusivities from EMD (DMD) scale linearly with the inverse of the simulation box length (L). Based on hydrodynamic arguments, these authors derived an analytic correction (hereafter referred to as the YH correction, DYH) for the self-diffusivity, D, of species i:in which kB is the Boltzmann constant, η is the shear viscosity of the system computed from EMD, and D∞ is the self-diffusivity in the thermodynamic limit. η computed in EMD does not show finite-size effects.[22,26,41] ξ is a constant which depends on the shape of the simulation box[42−45] (for a cubic simulation box,[22] ξ = 2.837297). The validity of eq has been extensively verified for various conditions and types of molecules.[22,26,41,46] Moultos et al.[41] showed that the YH correction holds for nonspherical molecules when a minimum number of 250 molecules is used in the simulation. In a recent study,[46] we showed that the YH finite-size correction also holds for self-diffusivities of mixtures. This finding allowed us to develop a method to compute the shear viscosity of a mixture from the computed self-diffusivities of its constituent components.[46]

In our previous work,[26] finite-size effects of binary mutual diffusivities were studied. For a binary mixture, the matrix form of eq reduces to an algebraic equation as there is a single DFick, ĐMS, and Γ. By performing simulations for more than 200 distinct Lennard-Jones and 9 molecular systems, we derived the following phenomenological correction for the computed (i.e., finite-size) MS diffusivity, ĐMSMD:[26]By combining eq and eq for a binary mixture, it can be shown that finite-size effects of the binary Fick diffusion coefficient require the same correction as self-diffusivities:

Simulation Details

The open-source software package LAMMPS[47] (version 16, Feb. 2016) was used to perform MD simulations for two sets of simulations: ternary Lennard-Jones (LJ) systems and a ternary molecular mixture. To compute transport properties and thermodynamic factors of these systems, the OCTP plugin was used.[48] The scheme used for computing Maxwell–Stefan (MS) diffusivities from Onsager coefficients and thermodynamic factors from Kirkwood–Buff coefficients is described in detail in the main text and Supporting Information of our published papers.[26,28,48,49]

Molecular Mixture

MD simulations were performed for a ternary mixture of (1) chloroform, (2) acetone, and (3) methanol at 298 K and 1 atm, corresponding to a density of 1025 kg/m3. The mole fractions of these components in the mixtures are 0.3, 0.3, and 0.4, respectively. Similar to the work of Liu et al.,[38] force field parameters for methanol,[50] acetone,[51] and chloroform[52] are obtained from literature. Schematic representations of these molecules are shown in Figure S1 in the Supporting Information, and force field parameters are listed in Table S2 in the Supporting Information. All molecules are considered rigid. The bond lengths and angles are listed in Table S3 in the Supporting Information. LJ interactions are truncated at a cutoff radius of 12 Å, and analytic tail corrections are included for the energy and pressure of the system.[20] Long-range electrostatic interactions are considered by using the particle–particle particle–mesh (PPPM) method with a relative precision of 10–6.[20] Diffusion coefficients were computed for four system sizes of 250, 500, 1000, and 2000 molecules. Initial molecular configurations were constructed in PACKMOL,[53] and LAMMPS input files were generated in VMD.[54] A time step of 1 fs was used, and all simulations have a total length of 100 ns. A total of 100 independent simulations were performed to obtain low statistical uncertainties. Thermodynamic factors (Γ) for this system are obtained from the work of Liu et al.[38] Γ11, Γ12, Γ21, and Γ22 are 0.61, −0.40, −0.31, and 0.79, respectively.

Lennard-Jones Systems

In this study, 28 distinct LJ ternary systems are considered. Properties are reported in reduced units, where the LJ energy (ϵ), size (σ), and mass (m) are equal to 1. To create asymmetric systems, the LJ energy parameter of each species is different (ϵ1 = 1.0, ϵ2 = 0.8, and ϵ3 = 0.6), while the sizes and masses for all particles are equal (σ1 = σ2 = σ3 = σ and m1 = m2 = m3 = m). All 28 ternary LJ systems correspond to homogeneous mixtures away from demixing. Interactions between dissimilar LJ particles are calculated according to the Lorentz–Berthelot mixing rules with an adjustable parameter (k) for the energy parameters:[20]The k parameters for the 28 ternary systems are listed in Table S1 in the Supporting Information. LJ interactions are truncated at a cutoff radius of 4σ. Analytic tail corrections are considered for energy and pressure calculations.[20] The mole fraction of species 1 is 0.4, and species 2 and 3 have mole fractions equal to 0.3. Simulations were performed at a reduced temperature of 0.65 and a reduced pressure of 0.05. The corresponding densities of these systems range from 0.78 to 0.89. Transport properties were computed in the microcanonical (NVE) ensemble for four system sizes consisting of 500, 1000, 2000, and 4000 LJ particles. To achieve sufficient statistics, 100 independent simulations with different initial configurations were performed for each system. To calculate radial distribution functions (RDFs) and consequently thermodynamic factors, MD simulations were performed in the canonical (NVT) ensemble. The statistical uncertainties of thermodynamic factors were obtained by performing 5 independent simulations, each one consisting of 25000 LJ particles. A time step of 0.001 was used to integrate equations of motion by using the velocity Verlet algorithm.[20] The total simulation lengths for computing transport properties and RDFs are 100 and 10 million timesteps, respectively.

Results and Discussion

Generalized Correction for the Finite-Size Effects in Fick Diffusivities

Equation is used as the initial step for investigating finite-size effects in multicomponent mixtures. Fick diffusion coefficients were computed for 4 system sizes for the ternary mixture (1) chloroform/(2) acetone/(3) methanol. Simulation details can be found in the Supporting Information. The finite-size effects are shown in Figure for all elements of the Fick diffusivity matrix. The extrapolation to the thermodynamic limit is performed by fitting a line to the 4 finite-size diffusivities. Figure strongly suggests that the off-diagonal elements of the Fick diffusivity matrix do not show any system-size dependence. In sharp contrast, the diagonal elements depend on the size of the simulation box and experience a finite-size effect of the same magnitude as the YH correction. Hence, we will hypothesize that the YH correction should be applied to the diagonal elements of the Fick diffusivity matrix, while no finite-size correction is needed for off-diagonal elements:where [DFickMD] and [DFick∞] are the Fick diffusivity matrices for a finite-size simulation box and in the thermodynamic limit, respectively. [I] is the identity matrix.

Figure 1

Computed elements of the Fick diffusivity matrix for the ternary mixture of (1) chloroform, (2) acetone, and (3) methanol (xchloroform = xacetone = 0.3) as a function of the inverse simulation box length L: (a) Diagonal element DFick1-1, (b) off-diagonal element DFick1-2, (c) off-diagonal element DFick2-1, and (d) diagonal element DFick2-2. Blue circles are the computed diffusion coefficients in MD simulations. Red squares are corrected diffusivities using the YH correction (eq ). Dashed lines show extrapolation to the thermodynamic limit, and solid lines are the extrapolated values. All data related to Fick diffusion coefficient computations are provided in the Simulation Details section and in the Supporting Information.

To examine the validity of this hypothesis, EMD simulations were performed to obtain the mutual diffusion coefficients of 28 different ternary LJ mixtures. For each mixture, 4 system sizes were simulated. All simulation details are provided in the Supporting Information. The differences between Fick diffusivities for the finite size and for the extrapolation to the thermodynamic limit are shown in Figure . Clearly, the results for the LJ and the molecular systems are in-line with our hypothesis (eq ).

Figure 2

Finite-size corrections required for elements of the Fick diffusivity matrix (a) DFick,1-1, (b) DFick,1-2, (c) DFick,2-1, and DFick,2-2 as a function of the YH correction (DYH, eq ) for 28 LJ systems containing 500 particles (blue circles), 1000 particles (red squares), 2000 particles (green diamonds), and 4000 particles (magenta pentagons). Reduced temperature is 0.65 and reduced pressure 0.05. The diagonal dashed lines indicate perfect agreement between the YH correction and the required finite size corrections for Fick diffusivities. The vertical dashed lines indicate no finite-size effects of Fick diffusivities. Raw data are provided in the Supporting Information.

The system-size dependency of multicomponent Fick diffusivities can also be investigated using the eigenvalue [D̂Fick] and modal [P] matrices of the Fick diffusivity matrix [DFick]. [D̂Fick] is a diagonal matrix containing all eigenvalues of the Fick diffusivity matrix, while [P] consists of all eigenvectors.[23] By definition, the modal matrix [P] has the property that[23]Eigenvalue and modal matrices are important for the linearized theory of multicomponent diffusion.[23,55,56] Importantly, the linearized theory of multicomponent diffusion allows transforming the (n – 1) equations of motion to (n – 1) pseudobinary eigen-species representing (n – 1) hydrodynamic diffusion modes. The speed of these hydrodynamic diffusion modes is expressed by the eigenvalues of the Fick matrix. Thus, we expect that the hydrodynamic YH correction should affect the eigenvalues of the Fick matrix. Techniques like light dynamic light scattering (DLS), that can access hydrodynamic modes, determine the eigenvalues of the Fick matrix in multicomponent mixtures.[57,58]

By definition, the ith eigenvalue of the matrix of finite-size Fick diffusivities (D̂MD) is calculated from the following equation:A similar equation can be written for the ith eigenvalue of Fick diffusivities in the thermodynamic limit (D̂∞):By combining eq with eqs and 14 we obtainThe last equation holds only when the term in the parentheses equals zero. Hence, all eigenvalues of the Fick diffusivity matrix need a finite-size correction equal to DYH. This means that the eigenvalue matrix of finite-size Fick diffusivities, [D̂FickMD], whose diagonal elements are the eigenvalues of [DFickMD], should be corrected according toEquation is a general expression according to which one should correct the finite size effects of Fick diffusivity computed in EMD. Although the validity of this expression is shown in Figures and 2 for ternary mixtures, the theoretical derivation presented in this section makes it valid for any n-component mixture.

A similar finite-size effect investigation can be performed for the modal matrix, [P]. Based on eq , for the Fick diffusivity in the thermodynamic limit, we can writeBy using eq and simplifying the identity matrix (i.e., [I] = [P∞]−1 [P∞]), one can rewrite eq asAccording to eq , the left-hand side of this equation can be substituted by [D̂FickMD] + DYH[I]. Therefore, the terms DYH[I] on both sides of the equation cancel out, and the following relation can be obtained:Similarly to eq , for finite-size diffusivities computed by EMD, one can writeBy comparing eqs and 25, one can conclude that P∞ = PMD. This means that the modal (eigenvector) matrix of Fick diffusivities does not have any system-size dependency, while according to eq the eigenvalue matrix of Fick diffusivities should be corrected. This finding leads to the conclusion that while the size of the simulation box affects the speed of diffusion (i.e., the eigenvalues), the direction of the diffusion process (i.e., the eigenvectors) is unaffected.

Generalized Correction for the Finite-Size Effects in Maxwell–Stefan Diffusivities

As Fick and MS diffusivity are connected via the thermodynamic factors (see eq ), a correction term for the finite-size MS diffusivities can be derived as follows:and by using [DFick∞] = [Δ∞] [Γ] in the equation above leads tofrom which the MS diffusivities can be computed. Equation indicates that both the diagonal and off-diagonal elements of the matrix of MS diffusivities show system-size dependency, in sharp contrast to the Fick diffusivity matrix.

In Figure , the results of eq are shown for the finite-size effects of the computed MS diffusivities for the ternary molecular mixture (1) chloroform/(2) acetone/(3) methanol. All three MS diffusivities can be corrected using eq , while the magnitudes of the required corrections are different. The validity of this correction was further examined for 28 ternary LJ systems (see Figure ). For all systems, very good agreement between the computed correction and eq is observed. Figures and 4 clearly show that the correction (i.e., eq ) can accurately predict the finite-size effects of the matrix of MS diffusivities. It is important to note here that, similarly to eq , eq is a generalized expression that can be used for any multicomponent mixture, independently of the number of components.

Figure 3

MS diffusion coefficients for a mixture of (1) chloroform, (2) acetone, and (3) methanol (xchloroform = xacetone = 0.3) as a function of the simulation box length (L). Blue circles are the computed diffusion coefficients in MD simulations. Red squares are the corrected diffusivities using the proposed correction in eq . Dashed lines show extrapolation to the thermodynamic limit, and solid lines are the extrapolated values.

Figure 4

Comparison between finite-size corrections required for MS diffusivities (ĐMS∞ – ĐMSMD) and the proposed correction for MS diffusivities (ĐMScorrection, eq , for 28 LJ systems containing 500 particles (blue circles), 1000 particles (red squares), 2000 particles (green diamonds), and 4000 particles (magenta pentagons). Parts 1-2, 1-3, and 2-3 indicate Đ1-2, Đ1-3, and Đ2-3, respectively. Reduced temperature is 0.65 and reduced pressure 0.05. The dashed lines indicate perfect agreement between the proposed correction and the required finite size corrections for MS diffusivities. Raw data are provided in the Supporting Information.

Conclusions

In this study, we investigated the finite-size dependency of mutual diffusion coefficients of multicomponent mixtures, and a generalized correction is derived, i.e., eq . While the off-diagonal elements of the Fick diffusivities do not show any system-size dependency, the diagonal elements should be corrected with the term proposed by Yeh and Hummer for the finite-size self-diffusivities.[22] An eigenvalue analysis of the finite-size effects of the matrix of Fick diffusivities revealed that the eigenvector matrix of Fick diffusivities is unaffected by the size of the simulation box. Only eigenvalues, which describe the speed of diffusion, depend on the size of the system. This is in-line with the hydrodynamic nature of the finite-size effects. An analytic relation for finite-size effects of the matrix of Maxwell–Stefan diffusivities was also developed. The finite-size correction term for the mutual diffusivities was examined for the ternary molecular mixture chloroform/acetone/methanol as well as 28 ternary LJ systems. All simulation results are in good agreement with the proposed corrections for Fick and Maxwell–Stefan mutual diffusion coefficients.