Diffusion in multicomponent aqueous alcoholic mixtures

The Fick diffusion coefficient matrix of the highly associating quaternary mixture water + methanol + ethanol + 2-propanol as well as its ternary and binary subsystems is analyzed with molecular dynamics simulation techniques. Three of the ternary subsystems are studied in this sense for the first time. The predictive capability of the employed force fields, which were sampled with the Green–Kubo formalism and Kirkwood–Buff integration, is confirmed by comparison with experimental literature data on vapor-liquid equilibrium, shear viscosity and Fick diffusion coefficient, wherever possible. A thorough analysis of the finite size effects on the simulative calculation of diffusion coefficients of multicomponent systems is carried out. Moreover, the dependence of the Fick diffusion coefficient matrix on the velocity reference frame and component order is analyzed. Their influence is found to be less significant for the main matrix elements, reaching a maximum variation of 19%. The large differences found for the cross elements upon variation of the reference frame hinder a straightforward interpretation of the Fick diffusion coefficient matrix with respect to the presence of diffusive coupling effects.

. Elements of the Fick diffusion coefficient matrix in the molar (black circles), volume (blue squares) and mass reference frame (red triangles) of the ternary subsystem water (1) + methanol (2) + ethanol (3) with x3 = 0.25 mol·mol −1 at 298.15 K and 0.1 MPa. The green crosses represent the expected value in the binary limit x1 → 0 for the molar reference frame.
Numerical results TABLE S1. Fick diffusion coefficient matrix in the molar, volume and mass reference frame of the ternary subsystem methanol (1) + ethanol (2) + 2-propanol (3) at 298.15 K and 0.1 MPa. The numbers in parentheses indicate the uncertainty in the last given digit.
No.  No.  No.   The asymptotic behavior towards the binary limits of the ternary Fick diffusion coefficient matrix can be obtained by comparing the diffusive flux equations of the ternary and binary mixtures. If one of the first two components is vanished, the following limits are given: Therein, D ij is the element in the ith row and the jth column of the ternary Fick diffusion coefficient matrix and D bin ij is the Fick diffusion coefficient of the binary subsystem i + j.

Frame of reference
The Fick diffusion coefficient matrix depends on the velocity reference frame. The molar reference frame is the natural choice for molecular simulations, but experimental data are typically evaluated on the volume reference frame. On the other hand, in non-equilibrium thermodynamics, usually the mass reference frame is assumed. To compare among different approaches, it is important to be aware how the according Fick diffusion coefficients relate to each other. The Fick diffusion coefficient matrix in the molar reference frame D M can be transformed into its form in the mass reference frame D w employing [2] where [x] and [w] are the diagonal matrices of the mole and mass fractions x i and w i , respectively. The elements of the matrices B wu and B uw are given by where δ ik is the Kronecker delta and the lower index n refers to the solvent. Similarly, D M can be transformed to the volume reference frame D V with and where the total molar volume v is the weighted sum of the partial molar volumes

Choice of the solvent
In a multicomponent mixture, there are several ways to order the components. Usually, the species with the highest concentration is chosen as the solvent (s) due to accuracy concerns. However, the choice of component numbering is arbitrary.
The Fick approach for quaternary mixtures involves three independent molar diffusion fluxes where the solvent component is s=4, c t is the mixture molar density and D ij are the elements of the Fick diffusion matrix in the molar reference frame. Varying component order changes the values of the Fick diffusion coefficient matrix, but not the actual molar fluxes. Since − (D 13 + D 23 + D 33 )∇x 4 .

Force fields
Molecular dynamics simulations rely on force fields that mimic the intermolecular interactions adequately. In this work, rigid and non-polarizable force fields of united-atom type were employed, which account for these interactions by a set of Lennard-Jones sites and point charges which may or may not coincide with respect to their site positions. The force fields for the three alcohols were developed by our group based on quantum chemical calculations, parameter optimization to experimental vapor-liquid equilibrium data and, in the case of 2-propanol, also to experimental selfdiffusion coefficient data [3][4][5][6]. For water, the TIP4P/2005 force field by Abascal and Vega [3] was employed. This force field was found to predict the transport properties of water and aqueous alcoholic mixtures with a better accuracy than other commonly used non-polarizable force fields [7]. The employed force fields have been tested with respect to their ability to predict structural, thermodynamic and transport properties of the corresponding pure substances [6][7][8], some binary [6,8,9], ternary [9,10] as well as their quaternary mixture [11]. Detailed information about these force fields can be found in Table S7 and the original publications [3][4][5][6].
To define a molecular model for a mixture on the basis of pairwise additive pure substance force fields, only the unlike interactions have to be specified. In case of the point charges, this can straightforwardly be done with Coulomb's law. However, for the unlike Lennard-Jones parameters, there is no physically sound approach so that combining rules have to be employed. The Lorentz-Berthelot combining rules were chosen here σ ab = (σ aa + σ bb )/2, and ε ab = √ aa ε bb , so that all present mixture data are strictly predictive. Transport data were sampled by equilibrium molecular dynamics simulation and the Green-Kubo formalism based on the net velocity auto-correlation function to obtain the phenomenological coefficient matrix [12] Here, N is the total number of molecules. N i and N j are the number of molecules of components i and j, respectively. v i,k (t) is the center of mass velocity vector of the k-th molecule of component i at time t. The brackets <...> denote the canonical (N V T ) ensemble average and equation (S18) corresponds to a reference frame in which the mass-averaged velocity of the mixture is zero [12].
With the phenomenological coefficients L ij , the elements of a matrix ∆ can be defined [12] which is related to the matrix B by its inverse, B = ∆ −1 .
In the case of a binary mixture, the Maxwell-Stefan diffusion coefficient Ð can be calculated with In the ternary case, the three Maxwell-Stefan diffusion coefficients are given by , . (S21)

Shear viscosity
The shear viscosity η was calculated concurrently with the phenomenological coefficients employing the Green-Kubo formalism. It is associated with the off-diagonal elements of the microscopic stress tensor J xy where V is the volume, k B is the Boltzmann constant and T die temperature.
Here, k and l denote different molecules of any species. The upper indices x and y stand for the spatial vector components, e.g. for velocity v x k or site-site distance r x kl . u(r kl ) is intramolecular the potential energy. Equations (S22) and (S23) may directly be applied to mixtures. To improve statistics, five independent terms of the stress tensor J xy p , J xz p , J yz p , (J xx p − J yy p )/2 and (J yy p − J zz p )/2 were considered [13].

Technical details
Molecular dynamics simulations were performed with the program ms2 [14,15] in two steps: First, a simulation in the isobaric-isothermal (N pT ) ensemble was carried out to calculate the density at the desired temperature, pressure and composition. In the second step, a canonic (N V T ) ensemble simulation was performed under the corresponding thermodynamic conditions to simultaneously determine the phenomenological coefficient and thermodynamic factor matrices. Newton's equations of motion were solved with a fifth-order Gear predictor-corrector numerical integrator and the temperature was controlled by velocity scaling. Throughout, the integration time step was 0.925 fs. The simulations contained 6000 molecules and were carried out in a cubic volume with periodic boundary conditions, where the cut-off radius was set to r c = 24.5 Å. Lennard-Jones long range interactions were considered using angle averaging [16]. Electrostatic long-range corrections were approximated by the reaction field technique with conducting boundary conditions ( RF = ∞).
The simulations in the N pT ensemble were equilibrated over 4 × 10 5 time steps, followed by a production run over 3 × 10 6 time steps. In the N V T ensemble, the simulations were equilibrated over 5 × 10 5 time steps, followed by production runs of 8 ×10 7 time steps. The phenomenological coefficients were calculated for up to 8 × 10 5 independent time origins of the correlation functions. The sampling length of the correlation functions was of 20 000 time steps or 18.5 ps throughout. The separation between the time origins was chosen such that all autocorrelation functions have decayed at least to 1/e of their normalized value to achieve their time independence [17]. The uncertainties of the predicted values were estimated with a block averaging method [18].

IV. FINITE SIZE EFFECTS
To asses finite size correction methods, simulations were performed for the three ternary subsystems methanol + ethanol + 2-propanol, water + methanol + 2-propanol and water + ethanol + 2-propanol with system sizes containing between 512 and 8000 molecules with composition x 1 = 0.125, x 2 = 0.625 and x 3 = 0.25 mol·mol −1 for all mixtures. Additionally, the ternary mixture water + ethanol + 2-propanol was regarded for the composition x 1 = 0.25, x 2 = 0.5 and x 3 = 0.25 mol·mol −1 . System size effects of the intra-, phenomenological, Maxwell-Stefan and Fick diffusion coefficients were evaluated by plotting the simulation results over the inverse of the edge length of the simulation volume L −1 and fitting a straight line to the data. The resulting intercept for L −1 → 0 corresponds to the infinite size value.  Coefficients subject to the fast correction procedure [11] (crosses) are compared with those corrected with the procedure by Jamali et al. [19] (cyan triangles). The green squares represent the Fick diffusion coefficients based on the individually extrapolated phenomenological coefficients.

Intra-diffusion coefficients
FIG. S16. Elements of the Fick diffusion coefficient matrix of the ternary subsystem water (1) + methanol (2) + 2-propanol (3) (xH2O = 0.25, xCH4O = 0.5 and xC3H8O = 0.25 mol·mol −1 ) as a function of the inverse edge length of the simulation volume L −1 at 298.15 K and 0.1 MPa. The blue dashed line is a linear fit to the uncorrected simulation results (blue bullets). Coefficients subject to the fast correction procedure [11] (crosses) are compared with those corrected with the procedure by Jamali et al. [19] (cyan triangles). The green squares represent the Fick diffusion coefficients based on the individually extrapolated phenomenological coefficients.