This version is available at https://doi.org/10.14279/depositonce-10426
Copyright applies. A non-exclusive, non-transferable and limited
right to use is granted. This document is intended solely for
personal, non-commercial use.
Terms of Use
Guevara-Carrion, G., Fingerhut, R., & Vrabec, J. (2020). Fick Diffusion Coefficient Matrix of a Quaternary
Liquid Mixture by Molecular Dynamics. The Journal of Physical Chemistry B, 124(22), 4527–4535.
https://doi.org/10.1021/acs.jpcb.0c01625
Gabriela Guevara-Carrion, Robin Fingerhut, Jadran Vrabec
Fick Diffusion Coefficient Matrix of a
Quaternary Liquid Mixture by Molecular
Dynamics
Accepted manuscript (Postprint)Journal article |
Fick Diffusion Coefficient Matrix of a
Quaternary Liquid Mixture by Molecular
Dynamics
Gabriela Guevara-Carrion, Robin Fingerhut, and Jadran Vrabec∗
Thermodynamics and Process Engineering, Technical University Berlin, 10587 Berlin, Germany
E-mail: vrabec@tu-berlin.de
Phone: +49 30 314 22755
1
Abstract
For the first time, the Fick diffusion coefficient matrix of a quaternary liquid mixture is
sampled consistently by means of molecular dynamics simulation. The required phenomeno-
logical diffusion coefficient and thermodynamic factor matrices of the mixture water + methanol
+ ethanol + 2-propanol are determined following the Green-Kubo formalism and Kirkwood-
Buff theory. Further, a system size correction methodology for the Fick diffusion coefficient of
multicomponent mixtures is proposed. Ten compositions are studied under ambient conditions
and validated by analyzing the ternary limits of the quaternary Fick diffusion matrix. Because
of complex intermolecular interactions due to the presence of hydrogen bonding, the elements
of the Fick diffusion coefficient matrix exhibit a significant composition dependence. The
magnitude of several cross coefficients indicate important coupling effects mainly affecting
the diffusive flux of water. These effects are explained in the light of the structural infor-
mation given by the radial distribution functions of the mixture. This work that solely rests
on molecular dynamics simulation techniques to predict the Fick diffusion coefficient matrix
of quaternary mixtures is expected to be a significant step forward for the understanding of
multicomponent diffusion.
Graphical TOC Entry
2
Introduction
Multicomponent solutions are involved in the majority of mass transfer processes occurring in
nature and technical applications1. Almost all separation processes in chemical engineering, such
as distillation, absorption or extraction, are affected by multicomponent diffusion in liquids. Rate-
based methods employed for modeling, design and control of these unit operations involve mass
and energy transfer models, which require diffusion coefficient data for mixtures2.
Mass transport in multicomponent systems is complex. For instance, describing isothermal-
isobaric diffusion of a quaternary mixture by Fick’s law requires a matrix with nine different diffu-
sion coefficient elements that are composition dependent. Related experimental work is challeng-
ing not only because of the elaborate equipment that is required to distinguish between different
components, but also because of the presence of coupling effects that hinder data processing and
interpretation3. Consequently, the availability of experimental data on diffusion coefficients for
mixtures containing three or more components is very limited4. In fact, they have experimentally
been measured only for 17 quaternary mixtures5–17, but in some of these cases, only the main ele-
ments of the diffusion matrix are reported. With diffusion having entered the scientific arena in the
1850s through the contributions of Graham18 and Fick19, experimental measurements alone are
obviously not able to satisfy the growing need for accurate mass transport properties20, particularly
for liquids that are constituted of many components.
Most predictive equations for diffusion coefficients of multicomponent liquids are extensions
of the Darken relation21–23, which is not valid for mixtures with strong intermolecular interactions.
The underlying physical phenomena are not well understood and the lack of experimental data im-
pedes the development and verification of new predictive equations. On the other hand, molecular
dynamics simulation is a compelling alternative for such predictions. In fact, such simulations
have the capability to accurately predict Fick diffusion coefficients of binary24,25 and ternary mix-
tures26,27. Krishna and van Baten21 were the first to deduce the required mathematical framework
and sampled the Maxwell-Stefan diffusion coefficient matrix of quaternary n-alkane mixtures with
molecular dynamics. Later, Liu et al.22 determined the Maxwell-Stefan diffusion coefficients of
3
quaternary model mixtures interacting with the Weeks-Chandler-Andersen potential. However,
none of these works predicted the Fick diffusion coefficient, perhaps because of the lack of a
mathematical framework to sample the required thermodynamic factor matrix. In a recent work
of our group28, expressions for the thermodynamic factor matrix of quaternary mixtures based on
Kirkwood-Buff theory were derived to close this gap. To the best of our knowledge, this is the
first work on the prediction of quaternary Fick diffusion coefficients solely based on molecular
dynamics simulation techniques.
The quaternary mixture water + methanol + ethanol + 2-propanol was chosen in this context
since the Fick diffusion coefficient of most of the involved subsystems has been predicted in pre-
vious work. In fact, the binary subsystems water + methanol29, water + ethanol29, methanol +
ethanol30, water + 2-propanol25, and the ternary subsystem water + methanol + ethanol26,27 have
successfully been studied in this sense.
Theory
In the framework of the generalized form of Fick’s law, the molar flux of component iin a mixture
of four components is written as a linear combination of concentration gradients ∇cj31
Ji=−
3
∑
j=1
Di j∇cj,(i=1,2,3),(1)
where Dii are the main diffusion coefficients that relate the molar flux of component ito its own
concentration gradient and Di j are the cross diffusion coefficients that relate the molar flux of
component ito the concentration gradient of component j32. The Fick approach involves three
independent diffusion fluxes and a 3×3 diffusion coefficient matrix, which is generally not sym-
metric, i.e. Di j 6=Dji. Further, the numerical values of Di j depend both on the reference frame for
velocity (molar-, mass- or volume-averaged) and on the order of the components. In this work, the
molar-averaged reference frame is employed throughout.
The main shortcoming of Fick’s law is the fact that concentration gradients are not the true
4
thermodynamic driving forces for diffusion, which are rather given by chemical potential gradi-
ents. Maxwell-Stefan theory follows this path, assuming that chemical potential gradients ∇µiare
balanced by friction forces between the components that are proportional to their mutual velocity
(u
u
ui−u
u
uj)33
4
∑
j6=i=1
xj(u
u
ui−u
u
uj)
Ði j
=−1
kBT∇µi(i,j=1,...,4),(2)
where xjis the mole fraction of component j,kBBoltzmann’s constant and Tthe temperature. The
Maxwell-Stefan diffusion coefficient Ði j plays the role of an inverse friction coefficient between
components iand j34. Its matrix does not depend on the component order and is symmetric so that
it has only six independent elements.
Maxwell-Stefan diffusion coefficients are associated with chemical potential gradients and thus
cannot directly be measured in the laboratory. However, they are accessible with equilibrium
molecular simulation techniques, i.e. the Green-Kubo formalism35,36 or the Einstein approach37.
For multicomponent mixtures, Fick and Maxwell-Stefan diffusion coefficients are related by33
D=B−1·Γ
Γ
Γ,(3)
in which all three symbols represent 3×3 matrices and the elements of Bare given by21
Bii =xi
Ði4+
4
∑
j6=i=1
xj
Ði j
,Bi j =−xi1
Ði j −1
Ði4.(4)
Fick diffusion coefficients can be calculated from the Maxwell-Stefan diffusion coefficients if the
thermodynamic factor matrix Γ
Γ
Γ
Γi j =δi j +xi
∂lnγi
∂xjT,p,xk,k6=j=1...3
,(5)
is known. Therein, δi j is the Kronecker delta function and γithe activity coefficient of component
i. The partial derivative has to be evaluated at constant temperature, pressure and mole fraction of
5
all other components.
Ternary limits of quaternary Fick diffusion coefficients
The consistency of the present simulation results can be assessed by analyzing the asymptotic
behavior of quaternary diffusion coefficients when approaching the ternary limits. This type of
analysis has been developed for ternary mixtures27,38 and is expanded here to quaternary mixtures.
Following the previously outlined procedure27, the asymptotic behavior of some elements of the
quaternary Fick diffusion coefficient matrix in the molar-averaged reference frame can be obtained.
When the water content of the studied quaternary mixture water (1) + methanol (2) + ethanol (3)
+ 2-propanol (4) vanishes, ∇x1→0, the mixture approaches its ternary subsystem methanol +
ethanol + 2-propanol. By comparing the diffusive flux equations of the quaternary and ternary
mixtures, it follows that
Dquat
22 →Dtern
11 ,Dquat
23 →Dtern
12 ,Dquat
32 →Dtern
21 ,Dquat
33 →Dtern
22 ,(6)
where Dquat
i j and Dtern
i j denote the respective elements of the quaternary and ternary Fick diffusion
coefficient matrices.
From a similar analysis for vanishing methanol content, x2→0, it follows for the ternary limit
water + ethanol + 2-propanol that
Dquat
11 →Dtern
11 ,Dquat
13 →Dtern
12 ,Dquat
31 →Dtern
21 ,Dquat
33 →Dtern
22 .(7)
Analogously, when ethanol disappears from the mixture, x3→0, the asymptotic behavior to-
wards the ternary limit water + methanol + 2-propanol is given by
6
Dquat
21 →Dtern
21 ,Dquat
22 →Dtern
22 ,Dquat
11 →Dtern
11 ,Dquat
12 →Dtern
12 .(8)
The behavior of quaternary diffusion when approaching the ternary subsystem water + methanol
+ ethanol, x4→0, requires some transformation of the expressions for the diffusive fluxes. Em-
ploying x1+x2+x3+x4=1 and ∇x3=−∇x1−∇x2−∇x4, it follows that
Dquat
11 −Dquat
13 →Dtern
11 ,Dquat
12 −Dquat
13 →Dtern
12 ,
Dquat
21 −Dquat
23 →Dtern
21 ,Dquat
22 −Dquat
23 →Dtern
22 .(9)
A detailed derivation of Eqs. (6) to (9) is given in the supplementary material.
Methods
The Fick diffusion coefficient matrix of quaternary and ternary mixtures was calculated from Eqs.
(3) to (5) and matrices Band Γ
Γ
Γthat were sampled exclusively by means of equilibrium molecular
dynamics simulation techniques.
The primary requirement for this task is the availability of molecular models that mimic the in-
termolecular 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 molecular
models for the three alcohols were developed by our group based on quantum chemical calcu-
lations and parameter optimization to experimental vapor-liquid equilibrium and, in the case of
2-propanol, also to self-diffusion data25,39–41. For water, the TIP4P/2005 model by Abascal and
Vega39 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
7
fields29. The interested reader is referred to the original publications25,39–41 for detailed informa-
tion about the four molecular pure substance models and their parameters. It has been shown that
all molecular models are suitable for the prediction of structural, thermodynamic and transport
properties of the corresponding pure substances25,29,30 as well as four of the binary25,26,30 and one
of the ternary26,27 subsystems of the regarded quaternary mixture.
To define a molecular model for a mixture on the basis of pairwise additive pure substance
models, 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 simple
Lorentz-Berthelot combining rules were chosen here, i.e., σab = (σaa +σbb)/2 and εab =√εaaεbb,
so that the present mixture data are strictly predictive.
Phenomenological coefficients
Transport data were sampled by equilibrium molecular dynamics simulation and the Green-Kubo
formalism35,36 based on the net velocity auto-correlation function to obtain the 4×4 phenomeno-
logical coefficient matrix21
Li j =1
3NZ∞
0dt DNi
∑
k=1
vi,k(0)·
Nj
∑
l=1
vj,l(t)E.(10)
Here, Nis the total number of molecules, Nithe number of molecules of component iand vi,k(t)
the center of mass velocity vector of the k-th molecule of component iat time t. The brackets <...>
denote the canonical (NV T) ensemble average and Eq. (10) corresponds to a reference frame in
which the mass-averaged velocity of the mixture is zero21.
With the phenomenological coefficients Li j, the elements of a 3×3 matrix ∆
∆
∆can be defined21
∆i j = (1−xi)Li j
xj−Li4
x4−xi
4
∑
k=16=iLk j
xj−Lk4
x4,(11)
This matrix can be directly employed in Eq. (3), since it is related to the matrix Bby its inverse,
8
B=∆
∆
∆−1.
System size corrections
It has been shown that under periodic boundary conditions, long-range interactions can lead to
an important dependence on system size inducing systematic errors on the calculation of the self-
diffusion or intra-diffusion coefficients42–45. Yeh and Hummer43 found that hydrodynamic self-
interactions in finite periodic systems are mainly responsible for system size effects and derived
a correction term based on the shear viscosity ηand the edge length of the simulation volume L,
i.e. 2.837297·kBT/(6πηL). Later on, Heyes et al.44 performed an investigation of hard sphere
fluids over a wide range of thermodynamic conditions and concluded that the Yeh and Hummer
correction term is not always adequate and that a more complex density dependence is needed.
Recently, Jamali et al.46 proposed a correction for the Maxwell-Stefan diffusion coefficient based
on the Yeh and Hummer43 term and the thermodynamic factor, but this correction is only valid for
binary mixtures. To the best of our knowledge, there has been no attempt to perform any system
size corrections for mutual diffusion coefficients of multicomponent mixtures.
Instead of applying the system size corrections directly to the Maxwell-Stefan or Fick diffusion
coefficients, they were rather applied here to the phenomenological coefficients Li j, cf. Eq. (11).
Once their system size effects were assessed, corrected values of the phenomenological coefficients
were employed in Eq. (3) to obtain the corresponding Fick diffusion coefficients. For this purpose,
simulations were performed for nine system sizes containing between 512 and 6000 molecules for
compositions of the present mixture with the lowest and highest molar density. The normalized
simulation results were then plotted over 1/N1/3to asses the system size dependence. For all
phenomenological cross coefficients Li j, no clear system size dependence could be inferred, which
might be due to their large statistical uncertainty, so that no corrections were made. In the case
of the four main coefficients Lii, a significant system size dependence was observed. A straight
line fitted to these data was extrapolated to infinite size, 1/N1/3→0. Very similar size corrections
for all Lii coefficients were found. Consequently, all main coefficients were normalized with the
9
results for 6000 molecules and plotted together in Figure 1. The resulting intercept of a single
linear fit was then employed to obtain the relative size correction for all sampled main coefficients
Lii. Note that no density effects were observed. The size correction calculated in this way led to
an increase of the phenomenological coefficients Lii by approximately 6%. A correction following
Yeh and Hummer would have led to around 15% larger corrections.
Figure 1: System size dependence of the main phenomenological coefficients Lii of the quaternary
mixture x1= 0.5 mol mol−1,x2= 0.125 mol mol−1and x3= 0.125 mol mol−1at 298.15 K and 0.1
MPa.
Thermodynamic factor
The thermodynamic factor of the quaternary mixture was estimated from information on the mi-
croscopic structure given by radial distribution functions gi j(r)based on Kirkwood-Buff theory.
Kirkwood-Buff integrals Gi j are defined in the grand canonical (µVT) ensemble47 by
Gi j =4πZ∞
0gi j(r)−1r2dr.(12)
10
When the canonical (NV T) ensemble is employed to sample this type of property, convergence
issues are possible48 so that corrections are required. In this context, the truncation and correc-
tion method developed by Krüger et al.49, that allows to obtain Kirkwood-Buff integrals from
simulations in the NV T ensemble, was applied here. Moreover, corrections of the radial distri-
bution functions are required. Therefore, Kirkwood-Buff integrals were calculated based on the
methodology proposed by Ganguly and van der Vegt50, which was found to be the most adequate
in previous work51. Extrapolation to the thermodynamic limit was not necessary because of the
rather large ensemble size N= 6000.
Expressions for the thermodynamic factor matrix of quaternary mixtures were derived in recent
work28 based on Ben-Naim’s formalism to determine particle number derivatives of the chemical
potential from Kirkwood-Buff integrals52. These lengthly equations have recently been published
elsewhere28 are not repeated here.
Results and Discussion
Predictive equilibrium molecular dynamics simulations of diffusion coefficients and the thermo-
dynamic factor of the quaternary mixture water (1) + methanol (2) + ethanol (3) + 2-propanol (4)
and its pure, binary and ternary and subsystems were carried out at 298.15 K and 0.1 MPa for the
compositions depicted in Figure 2. A total of ten quaternary compositions of this mixture, that lay
on the mole fraction plane x4= 0.25 mol mol−1, are discussed in this work.
To validate the simulation results from Kirkwood-Buff theory for the quaternary mixture, the
thermodynamic factor of all six binary subsystems was analyzed. For this purpose, additional
simulations were performed to sample the thermodynamic factor for all involved binary subsys-
tems and compared with the values obtained from a classical approach, i.e. a fit of the Wilson
excess Gibbs energy model53 to experimental vapor-liquid equilibrium data. A good agreement
was found between simulation and the classical approach, suggesting the trustworthiness of the
employed simulation methodology to access the chemical potential derivatives of the quaternary
11
Figure 2: Sampled compositions of the quaternary mixture water (1) + methanol (2) + ethanol (3)
+ 2-propanol (4) at 298.15 K and 0.1 MPa. The surface delimited by red lines indicates the plane
with a constant 2-propanol mole fraction x4= 0.25 mol mol−1.
12
mixture. Moreover, the Wilson model53 was fitted to data of the binary subsystems and the ob-
tained binary parameters were employed to predict the quaternary thermodynamic factor matrix for
the studied compositions. The resulting Wilson-based quaternary thermodynamic factor matrix is
shown in comparison with the one sampled directly with molecular dynamics for selected compo-
sitions in Figure 3. More comparisons of this type are provided in the supplementary material. A
good agreement was found especially for the main elements of this matrix with a relative averages
deviations between 2.3 and 3.6%.
Figure 3: Elements of the thermodynamic factor matrix Γi j for three selected quaternary compo-
sitions; left: x1= 0.125 mol mol−1,x2= 0.125 mol mol−1and x3= 0.5 mol mol−1; center: x1=
0.25 mol mol−1,x2= 0.125 mol mol−1and x3= 0.375 mol mol−1; right: x1= 0.5 mol mol−1,x2=
0.125 mol mol−1and x3= 0.125 mol mol−1. Present quaternary results based on Kirkwood-Buff
integration (red symbols) are compared with Γi j calculated with the Wilson excess Gibbs energy
model fitted to binary simulation results (blue symbols). The statistical uncertainties are within
symbol size.
The main element Γ11 of the thermodynamic factor matrix related to water evidences the high-
est non-ideality, yielding values that are considerably lower than unity, while the other two main
elements attain values close to unity. The cross elements of the thermodynamic factor matrix can
be of either sign and have values ranging from -0.2 to 0.2 for the studied compositions, cf. Figure 3.
The Fick diffusion matrix calculated with Eq. (3) for ten compositions of the quaternary mix-
13
ture is given in Table 1. Its main elements are positive and have values between 3.7 and 12.7×10−10
m2s−1. Further, the main element related to methanol is the highest and that related to water the
lowest, i.e. D22 >D33 >D11. Lower values for the main diffusion coefficients of water can be
associated with the presence of clustering, exerted by strong hydrogen bonding networks, which
hinders diffusion.
In order to visualize and test the consistency of the calculated Fick quaternary diffusion coeffi-
cient matrix, selected elements were plotted when the mole fraction of two species is kept constant
in the mixture. Exemplarly, Figure 4 shows the composition dependence of main and cross ele-
ments together with their expected asymptotic values, cf. Eqs. (6) to (8). Additional figures are
presented in the supplementary material. In most cases, the elements of the Fick diffusion coeffi-
cient matrix behave as expected, i.e. the sampled quaternary data converge to the projected value
when one species vanishes. This trend is clear for the main coefficients, but the cross coefficients
scatter, which may be a consequence to their inherently larger statistical uncertainties.
The main diffusion coefficient related to water D11 increases with increasing methanol con-
tent. This can be explained by the disturbance of the microscopic water-water structure through
the presence of water-methanol hydrogen bonding and is evidenced by changes in the magnitude
and location of the first and second peaks of the center of mass water-water radial distribution
function g11(r), cf. Figure 5. A similar observation was made for the somewhat weaker increase
of D11 with rising ethanol content; the corresponding radial distribution functions are shown in the
supplementary material.
The decrease of D22 and D33 with rising water content is consistent with the presence of water-
methanol and water-ethanol networks, which can be also deduced from the corresponding radial
distribution functions g12(r)and g13(r). When the water content is kept constant, the main coeffi-
cient related to ethanol D33 decreases with rising ethanol mole fraction due to a strengthening of
the ethanol-ethanol hydrogen bonding network, cf. Figure 4. A small decrease of D22 with rising
ethanol mole fraction was also observed for a constant water content, which suggests an increase
of the presence of methanol-ethanol networks when the number methanol molecules available for
14
Table 1: Fick diffusion coefficient matrix and its three eigenvalues of the quaternary mix-
ture water (1) + methanol (2) + ethanol (3) + 2-propanol (4) at 298.15 K and 0.1 MPa. The
numbers in parentheses indicate the uncertainty in the last given digit.
x1x2x3Dˆ
D
mol mol−110−10m2s−110−10m2s−1
0.125 0.125 0.5
7.1(2)−1.0(3)−0.7(1)
0.0(2)10.9(3)0.3(1)
0.2(2)−0.7(4)8.9(2)
7.1
10.8
8.9
0.125 0.25 0.375
7.9(2)−1.0(2)−0.4(2)
0.4(2)10.8(2)−0.2(2)
−0.4(2)−0.3(3)9.2(2)
7.8
10.7
9.3
0.125 0.375 0.25
8.6(2)−1.2(2)−0.7(2)
0.1(3)12.7(2)−0.1(2)
0.0(2)−0.7(2)11.0(2)
8.6
12.7
11.0
0.125 0.5 0.125
9.8(2)−0.9(2)−0.7(2)
0.0(2)12.6(2)−0.6(3)
−0.3(2)−0.6(1)11.3(2)
9.6
12.8
11.3
0.25 0.125 0.375
5.8(1)−0.9(3)−0.5(2)
1.3(1)10.3(2)−0.1(1)
−0.5(1)−1.5(3)7.9(1)
6.1
10.2
7.7
0.25 0.25 0.25
6.3(1)−1.5(2)−0.6(2)
0.7(1)10.4(2)−0.1(2)
0.0(1)−0.6(2)8.5(2)
6.6
10.1
8.4
0.25 0.375 0.125
7.2(1)−1.3(2)−0.6(3)
−0.4(1)11.2(2)−0.5(3)
0.9(1)−0.8(1)9.1(2)
7.4
11.5
8.6
0.375 0.125 0.25
4.3(1)−1.6(3)−1.1(2)
1.3(1)9.0(2)−0.1(1)
−0.28(8)−0.7(2)7.4(2)
4.8
8.7
7.2
0.375 0.25 0.125
5.1(1)−2.0(3)−1.1(3)
0.1(1)10.2(2)0.1(2)
0.97(7)−0.4(1)8.0(2)
5.5
10.1
7.6
0.5 0.125 0.125
3.7(1)−2.3(3)−0.8(3)
0.55(6)8.6(2)−0.3(2)
0.41(5)−0.3(1)6.6(1)
4.2
8.5
6.2
15
Figure 4: Main and cross elements of the Fick diffusion coefficient matrix as a function of methanol
mole fraction for selected compositions with x1= 0.125 mol mol−1and x4= 0.25 mol mol−1.
Simulation results for the quaternary mixture (blue symbols) are shown together with the results
of the diffusion matrix of the ternary subsystem where x2→0 (green symbol) and the predictive
equations by Allie-Ebrahim et al.23 (red symbols).
16
0
2
4
6
8
10
2 4 6 8
g11 (r)
r/ Å
x1
Figure 5: Water-water radial distribution function g11(r)for varying water mole fraction x1= 0.125
mol mol−1(black line), 0.25 mol mol−1(blue line), 0.375 mol mol−1(red line) and 0.5 mol mol−1
(green line) with constant x3= 0.125 mol mol−1and x4= 0.25 mol mol−1.
17
self-association is reduced. This explanation is in line with the changes observed for the radial
distribution function g23(r), i.e. an increase of the radii of the first and second coordination shells,
cf. supplementary material.
For most of the studied compositions, the cross elements Di j are not zero, i.e. there is an
influence of the concentration gradient of another species on the diffusion of component i. In fact,
both cross coefficients related to water have the same order of magnitude as the main coefficient
D11, which might lead to important coupling effects when a significant concentration gradient
of the corresponding alcohol is present. The ratio between cross and main elements of the Fick
diffusion coefficient matrix |Di j/Dii|provides a characterization of coupling effects. Figure 6
shows the composition dependence of this ratio for the water-related cross coefficients when the
mole fraction of two alcohols remains constant. In general, |D12/D11|increases with water content,
which suggests important water clustering around methanol by hydrogen bonding. In the case of
the cross coefficient D13, the tendency is not that clear, although a moderate increase with water
mole fraction can be inferred. Further, the influence of the concentration gradient of methanol on
water diffusion is greater than that of ethanol. This observation could be explained by the larger
hydrophobic tail of ethanol molecules, which reduce the occurrence and size of mixed water-
ethanol clusters.
In order to demonstrate the relationship between the microscopic structure and the coupling
effects, the average number of water(1) molecules in the first solvation shell of both alcohols ki1=
4πxiρRrc
0gi j(r)r2dr was calculated. Therein, istands for the methanol or ethanol site surrounded
by water, ρis the mixture density and rcis the radius of the first coordination shell, i.e. the location
of the first minimum of gi1at rc∼3.6 Å. These results are shown together with those from the
ratio |Di j/Dii|in Figure 6. It is clear that the coupling effects are related to the average number
of water molecules located in the first solvation shell of methanol or ethanol, i.e. water molecules
are transported together with methanol or ethanol molecules when a mass flux of either alcohol is
given.
On the other hand, most cross coefficients related to D22 and D33 are relatively small and the
18
Figure 6: Ratio between cross and main elements of the Fick diffusion coefficient matrix |Di j/Dii|
(top) and average number of water molecules in the first solvation shell of methanol or ethanol
(bottom) as a function of water mole fraction for selected quaternary compositions with x2= 0.125
mol mol−1and x4= 0.25 mol mol−1(left) as well as x3= 0.125 mol mol−1and x4= 0.25 mol
mol−1(right).
19
Fick matrix eigenvalues approach the values of the main coefficients, nearly implying an indepen-
dence of the diffusion fluxes of methanol and ethanol on the concentration gradients of the other
species, cf. Table 1.
The intra-diffusion coefficients of the four species in the mixture were sampled simultaneously
with the phenomenological diffusion coefficients. These data, together with the thermodynamic
factor, were employed to test the Darken-based predictive equations by Liu et al.22 and Allie-
Ebrahim et al.23 on their ability to reproduce the Fick diffusion coefficients calculated in this work.
While the predictive equation by Liu et al.22 was not able to return acceptable results, the equations
of Allie-Ebrahim et al.23 are in very good agreement with the present results for the main elements
of the Fick diffusion matrix with relative average deviations below 10%, cf. Figure 4. Further, in
most cases there is also an agreement for the cross elements within the statistical uncertainties of
the simulation data.
Conclusions
Diffusion processes in multicomponent liquid mixtures are of great importance in science and
engineering research. However, because of the serious difficulties associated with the measurement
of transport diffusion coefficients, the data availability for mixtures with three or more components
is very poor. Molecular dynamics simulation has been identified as an alternative to mitigate
data shortage for ternary mixtures. For the first time, as a powerful approach to Fick diffusion
coefficients of real quaternary liquid mixtures is presented here.
The Fick diffusion coefficient matrix of the strongly non-ideal liquid mixture water + methanol
+ ethanol + 2-propanol at 298.15 K and 0.1 MPa was sampled for ten compositions along the
plane with constant 2-propanol mole fraction, x4= 0.25 mol mol−1. The required phenomenologi-
cal coefficient matrix was calculated with the Green-Kubo formalism35,36 and the thermodynamic
factor matrix employing Kirkwood-Buff theory. Because of the lack of experimental data, only
consistency tests were made to verify the sampled Fick diffusion coefficients. However, present
20
predictions are expected to be accurate because convincing results were obtained in previous work
for four binary and one of the ternary subsystems on the basis of the same molecular models.
The asymptotic behavior of the elements of the quaternary Fick diffusion coefficient matrix was
analyzed and found to agree within the statistical uncertainties with the values of the ternary sub-
systems.
From the analysis of the Fick diffusion matrix, significant coupling effects were found mainly
for the diffusive flux for water, which were related to water-alcohol hydrogen bonding networks.
The lowest values for the main coefficient of water and its decrease with rising water content were
explained with the presence of water clusters.
Two Darken-based predictive models were tested and the predictive equations of Allie-Ebrahim
et al.23 were found to be in very good agreement with the Fick diffusion matrix predicted in this
work.
Acknowledgement
This work contributes to the Collaborative Research Center SFB-TRR 75 of Deutsche Forschungs-
gemeinschaft (DFG) and was funded under Grant No. VR 6/11. We gratefully acknowledge the
Paderborn Center for Parallel Computing (PC2) for the generous allocation of computer time on
the OCuLUS and Noctua clusters as well as the High Performance Computing Center Stuttgart
(HLRS) under the grant MMHBF2.
Supporting Information Available
Contains the detailed calculation of the ternary limits of quaternary Fick diffusion coefficients as
well the simulation procedure. The complete set of Figures for the main and cross elements of the
Fick diffusion coefficient matrix for the quaternary mixture as shown in Figure 4. The graphical
representation of selected radial distribution functions.
21
References
(1) Miller, D. G. Some Comments on Multicomponent Diffusion: Negative Main Term Diffusion
Coefficients, Second Law Constraints, Solvent Choices, and Reference Frame Transforma-
tions. J. Phys. Chem. 1986,90, 1509–1519.
(2) Kenig, E. Y.; Blagov, S. In Distillation: Fundamentals and Principles; Gòrak, A.,
Sorensen, E., Eds.; Elsevier: London, 2014.
(3) Guevara-Carrion, G.; Ancherbak, S.; Mialdun, A.; Vrabec, J.; Shetsova, V. Diffusion of
Methane in Supercritical Carbon Dioxide Across the Widom Line. Sci. Rep. 2019,9, 8466.
(4) Wambui Mutoru, J.; Firoozabadi, A. Form of Multicomponent Fickian Diffusion Coefficients
Matrix. J. Chem. Thermodyn. 2011,43, 1192–1203.
(5) Rai, G. P.; Cullinan, H. T. Diffusion Coefficients of Quaternary Liquid System Acetone-
Benzene-Carbon Tetrachloride-n-Hexane at 25 0C. J. Chem. Eng. Data 1973,18, 213–214.
(6) Noulty, R. A.; Leaist, D. G. Quaternary Diffusion in Aqueous Potassium Chloride-Potassium
Dihydrogen Phosphate-Phosphoric Acid Mixtures. J. Phys. Chem. 1987,91, 1655–1658.
(7) Leaist, D. G. Diaphragm-Cell Studies of Diffusion in the Four-Component System Hy-
drochloric Acid–Sodium Chloride–Sodium Iodide–Water. J. Chem. Soc. Faraday Trans. 1
1987,83, 829–839.
(8) Leaist, D. G. Diffusion Coefficients of Four-Component Systems from Taylor Dispersion
Profiles. Sucrose + LiCl + KCl + Water and NBu4Br + LiBr + KBr + Water. Ber. Bunsenges.
Phys. Chem. 1991,95, 117–122.
(9) Paduano, L.; Sartorio, R.; Vitagliano, V.; Albright, J. G.; Miller, D. G. Measurement of
the Mutual Diffusion Coefficients at One Composition of the Four-Component System
α-Cyclodextrln-L-Rhenylalanlne-Monobutylurea-H2O at 25 0C. J. Phys. Chem. 1992,96,
7478–7483.
22
(10) Leaist, D. G.; Hao, L. Quaternary Diffusion Coefficients of NaCl-MgCl2-Na2SO4-H2O Syn-
thetic Seawaters by Least-Squares Analysis of Taylor Dispersion Profiles. J. Solution Chem.
1993,22, 263–277.
(11) McEwan, K.; Leaist, D. G. Quaternary Mutual Diffusion Coefficients for Aqueous Solutions
of a Cationic–Anionic Mixed Surfactant from Moments Analysis of Taylor Dispersion Pro-
files. Phys. Chem. Chem. Phys. 2003,18, 3951–3958.
(12) Annunziata, O.; Vergara, A.; Paduano, L.; Sartorio, R.; Miller, D. G.; Albright, J. G. Precision
of Interferometric Diffusion Coefficients in a Four Component System Relevant to Protein
Crystal Growth: Lysozyme–Tetra(ethylene glycol)–NaCl–H2O. J. Phys. Chem. B 2003,107,
6590–6597.
(13) Vanag, V. K.; Rossi, F.; Cherkashin, A.; Epstein, I. R. Cross–Diffusion in a Water–in–Oil
Microemulsion Loaded with Malonic Acid or Ferroin. Taylor Dispersion Method for Four–
Component Systems. J. Phys. Chem. B 2008,112, 9058–9070.
(14) Annunziata, O.; Vergara, A.; Paduano, L.; Sartorio, R.; Miller, D. G.; Albright, J. G.
Quaternary Diffusion Coefficients in a Protein–Polymer–Salt–WaterSystem Determined by
Rayleigh Interferometry. J. Phys. Chem. B 2009,113, 13446–13456.
(15) Ribeiro, A. C. F.; Santos, C. I. A. V.; Lobo, V. M. M.; Esteso, M. A. Quaternary Diffusion Co-
efficients of β–Cyclodextrin + KCl + Caffeine + Water at 298.15 K Using a Taylor Dispersion
Method. J. Chem. Eng. Data 2010,55, 2610–2612.
(16) Rossi, F.; Vanag, V. K.; Tiezzi, E.; Epstein, I. R. Quaternary Cross-Diffusion in Water–in–Oil
Microemulsions Loaded with a Component of the Belousov–Zhabotinsky Reaction. J. Phys.
Chem. B 2010,114, 8140–8146.
(17) Peters, C.; Thien, J.; Wolff, L.; Koß, H.-J.; Bardow, A. Quaternary Diffusion Coefficients
in Liquids from Microfluidics and Raman Microspectroscopy: Cyclohexane + Toluene +
Acetone + Methanol. J. Chem. Eng. Data 2020,65, 1273–1288.
23
(18) Graham, T. The Bakerian Lecture.– On the Diffusion of Liquids. Philos. Trans. R. Soc. 1850,
140, 1–46.
(19) Fick, A. E. Über Diffusion. Pogg. Ann. Phys. Chem. 1855,94, 59–86.
(20) Wakeham, W. A. In Chemical Thermodynamics for Industry; Letcher, T. M., Ed.; The Royal
Society of Chemistry: London, 2004.
(21) Krishna, R.; van Baten, J. M. The Darken Relation for Multicomponent Diffusion in Liquid
Mixtures of Linear Alkanes: An Investigation Using Molecular Dynamics (MD) Simulations.
Ind. Eng. Chem. Res. 2005,44, 6939–6847.
(22) Liu, X.; Vlugt, T. J. H.; Bardow, A. Predictive Darken Equation for Maxwell-Stefan Diffu-
sivities in Multicomponent Mixtures. Ind. Eng. Chem. Res. 2011,50, 10350–10358.
(23) Allie-Ebrahim, T.; Russo, V.; Ortona, O.; Paduano, L.; Tesser, R.; Di Serio, M.; Singh, P.;
Zhu, Q.; Moggridge, G. D.; D’Agostino, C. A Predictive Model for the Diffusion of a Highly
Non-Ideal Ternary System. Phys. Chem. Chem. Phys. 2018,20, 18436–18446.
(24) Guevara-Carrion, G.; Janzen, T.; Munoz-Munoz, Y. M.; Vrabec, J. Mutual diffusion of Binary
Liquid Mixtures Containing Methanol, Ethanol, Acetone, Benzene, Cyclohexane, Toluene,
and Carbon Tetrachloride. J. Chem. Phys. 2016,144, 124501.
(25) Muñoz-Muñoz, Y. M.; Guevara-Carrion, G.; Vrabec, J. Molecular Insight into the Liquid
Propan–2–ol +Water Mixture. J. Phys. Chem. B 2018,122, 8718–8729.
(26) Parez, S.; Guevara-Carrion, G.; Hasse, H.; Vrabec, J. Mutual Diffusion in the Ternary Mixture
of Water + Methanol + Ethanol and its Binary Subsystems. Phys. Chem. Chem. Phys. 2013,
15, 3985–4001.
(27) Guevara-Carrion, G.; Gaponenko, Y. A.; Janzen, T.; Vrabec, J.; Shetsova, V. Diffusion in
Multicomponent Liquids: From Microscopic to Macroscopic Scales. J. Phys. Chem. B 2016,
120, 12193–12210.
24
(28) Fingerhut, R.; Herres, G.; Vrabec, J. Thermodynamic Factor of Quaternary Mixtures from
Kirkwood–Buff Integration. Mol. Phys. 2020,doi: 10.1080/00268976.2019.1643046.
(29) Guevara-Carrion, G.; Vrabec, J.; Hasse, H. Prediction of Self-diffusion Coefficient and Shear
Viscosity of Water and its Binary Mixtures with Methanol and Ethanol by Molecular Simu-
lation. J. Chem. Phys. 2011,134, 074508.
(30) Guevara-Carrion, G.; Nieto-Draghi, C.; Vrabec, J.; Hasse, H. Prediction of Transport Prop-
erties by Molecular Simulation: Methanol and Ethanol and their Mixture. J. Phys. Chem. B
2008,112, 16664–16674.
(31) Cussler, E. L. Mass Transfer in Fluid Systems, 2nd ed.; Cambridge University Press: Cam-
bridge, 1997.
(32) Buzatu, D.; Buzatu, F. D.; Paduano, L.; Sartorio, R. Diffusion Coefficients for the Ternary
System Water + Chloroform + Acetic Acid at 25 C. J. Solution Chem. 2007,36, 1373–1384.
(33) Taylor, R.; Krishna, R. Multicomponent Mass Transfer; John Wiley & Sons: New York, 1993.
(34) Krishna, R. Uphill diffusion in multicomponent mixtures. Chem. Soc. Rev. 2015,44, 2812–
2836.
(35) Green, M. S. Markoff Random Processes and the Statistical Mechanics of Time-Dependent
Phenomena. II. Irreversible Processes in Fluids. J. Chem. Phys. 1954,22, 398–414.
(36) Kubo, R. Statistical-Mechanical Theory of Irreversible Processes I. General Theory and Sim-
ple Applications to Magnetic and Conduction Problems. J. Phys. Soc. Jpn. 1957,12, 570–
586.
(37) Meier, K.; Laesecke, A.; Kabelac, S. Transport Coefficients of the Lennard-Jones Model
Fluid. II Self-Diffusion. J. Chem. Phys. 2004,121, 9526–9535.
25
(38) Sechenyh, V.; Legros, J. C.; Mialdun, A.; Ortiz de Zárate, J. M.; Shevtsova, V. Fickian Diffu-
sion in Ternary Mixtures Composed by 1,2,3,4-Tetrahydronaphthalene, Isobutylbenzene, and
n-Dodecane. J. Phys. Chem. B 2016,120, 535–548.
(39) Abascal, J. L. F.; Vega, C. A General Purpose Model for the Condensed Phases of Water:
TIP4P/2005. J. Chem. Phys. 2005,123, 234505.
(40) Schnabel, T.; Vrabec, J.; Hasse, H. Henry’s Law Constants of Methane, Nitrogen, Oxigen
and Carbon dioxide in Ethanol from 273 to 498 K: Prediction from Molecular Simulation.
Fluid Phase Equilib. 2005,233, 134–143.
(41) Schnabel, T.; Srivastava, A.; Vrabec, J.; Hasse, H. Hydrogen Bonding of Methanol in Super-
critical CO2: Comparison between 1H-NMR Spestroscopic Data and Molecular Simulation
Results. J. Phys. Chem. B 2007,111, 9871–9878.
(42) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford,
1987.
(43) Yeh, I. C.; Hummer, G. System-Size Dependence of Diffusion Coefficients and Viscosities
from Molecular Dynamics Simulations with Periodic Boundary Conditions. J. Phys. Chem.
B2004,108, 15873–15879.
(44) Heyes, D. M.; Cass, M. J.; Powles, J.; Evans, W. A. B. Self-Diffusion Coefficient of the Hard-
Sphere Fluid: System Size Dependence and Empirical Correlations. J. Phys. Chem. B 2007,
111, 1455–1464.
(45) Moultos, O. A.; Zhang, Y.; Tsimpanogiannis, I. N.; Economou, I. G.; Maginn, E. J. System–
Size Corrections for Self–Diffusion Coefficients Calculated from Molecular Dynamics Sim-
ulations: The Case of CO2, n–Alkanes, and Poly(ethylene Glycol) Dimethyl Ethers. J. Chem.
Phys. 2016,145, 074109.
26
(46) Jamali, S. H.; Wolff, L.; Becker, T. M.; Bardow, A.; Vlugt, T. J. H.; Moultos, O. A. Finite-
Size Effects of Binary Mutual Diffusion Coefficients from Molecular Dynamics. J. Chem.
Theory Comput. 2018,14, 2667–2677.
(47) Kirkwood, J. G.; Buff, F. P. The Statistical Mechanical Theory of Solutions. J. Chem. Phys.
1951, 774–777.
(48) Milzetti, J.; Nayar, D.; van der Vegt, N. F. A. Convergence of Kirkwood–Buff Integrals
of Ideal and Nonideal Aqueous Solutions Using Molecular Dynamics Simulations. J. Phys.
Chem. B 2018,122, 5515–5526.
(49) Krüger, P.; Schnell, S. K.; Bedeaux, D.; Kjelstrup, S.; Vlugt, T. J. H.; Simon, J.-M. Kirkwood-
Buff Integrals for Finite Volumes. J. Phys. Chem. Lett. 2013,4, 235–238.
(50) Ganguly, P.; van der Vegt, N. F. A. Convergence of Sampling Kirkwood–Buff Integrals of
Aqueous Solutions with Molecular Dynamics Simulations. J. Chem. Theory Comput. 2013,
9, 1347–1355.
(51) Fingerhut, R.; Vrabec, J. Kirkwood-Buff integration: A Promising Route to Entropic Proper-
ties? Fluid Phase Equilib. 2019,485, 270–281.
(52) Ben-Naim, A. Molecular Theory of Solutions; University Press: Oxford, 2006.
(53) Wilson, G. M. Vapor-Liquid Equilibrium. A New Expression for the Excess Free Energy of
Mixing. J. Am. Chem. Soc. 1964,86, 127–130.
27
Supplementary Material to:
Fick Diffusion Coefficient Matrix of a Quaternary
Liquid Mixture by Molecular Dynamics
Gabriela Guevara-Carrion, Robin Fingerhut, and Jadran Vrabec∗
Thermodynamics and Process Engineering, Technical University Berlin, 10587 Berlin,
Germany
E-mail: vrabec@tu-berlin.de
Phone: +49 30 314 22755
S1
Ternary limits of quaternary Fick diffusion coefficients
Here we present a complete analysis of the asymptotic behavior of diffusion coefficients for
quaternary mixtures. The three independent diffusive fluxes in the molar-averaged frame of
reference can be written as
−J1/ρ =Dquat
11 ∇x1+Dquat
12 ∇x2+Dquat
13 ∇x3,(1)
−J2/ρ =Dquat
21 ∇x1+Dquat
22 ∇x2+Dquat
23 ∇x3,(2)
−J3/ρ =Dquat
31 ∇x1+Dquat
32 ∇x2+Dquat
33 ∇x3.(3)
When the water content of the studied quaternary mixture water (1) + methanol (2) +
ethanol (3) + 2-propanol (4) vanishes, so does its as mass flux i.e.,∇x1→0and J1→0.
Then, from Eq. (1) it follows that Dquat
12 ∇x2+Dquat
13 ∇x3→0. Since ∇x2and ∇x3are
independent variables Dquat
12 →0and Dquat
12 →0. Further, Eqs. (2) and (3) can be written
as
−J2/ρ →Dquat
22 ∇x2+Dquat
23 ∇x3,(4)
−J3/ρ →Dquat
32 ∇x2+Dquat
33 ∇x3.(5)
The two independent diffusive fluxes for the ternary mixture methanol (2) + ethanol (3) +
2-propanol (4) are
−J2/ρ =Dtern
11 ∇x2+Dtern
12 ∇x3,(6)
−J3/ρ =Dtern
21 ∇x2+Dtern
22 ∇x3.(7)
S2
From the comparison of Eq. (4) with Eq. (6) as well as Eq. (5) with Eq. (7), it follows that
Dquat
22 →Dtern
11 , Dquat
23 →Dtern
12 , Dquat
32 →Dtern
21 , Dquat
33 →Dtern
22 ,(8)
A similar analysis can be done for x2→0and x3→0. However, the asymptotic behavior
for x4→0requires a slightly different approach. If x4→0, it follows from the condition
x1+x2+x3+x4= 1 that ∇x3→ −∇x1− ∇x2and the expressions for the diffusive fluxes
JV
1and JV
2take the form
−J1/ρ →Dquat
11 −Dquat
13 ∇x1−Dquat
11 −Dquat
13 ∇x2,(9)
−J2/ρ →Dquat
21 −Dquat
23 ∇x1−Dquat
22 −Dquat
23 ∇x2.(10)
Taking into account the two independent diffusive fluxes for the ternary mixture water (1)
+ methanol (2) + ethanol (3) leads to
−J1/ρ =Dtern
11 ∇x1+Dtern
12 ∇x2,(11)
−J1/ρ =Dtern
21 ∇x1+Dtern
22 ∇x2.. (12)
Applying the same logic as above, it follows that
Dquat
11 −Dquat
13 →Dtern
11 , Dquat
12 −Dquat
13 →Dtern
12 ,(13)
Dquat
21 −Dquat
23 →Dtern
21 , Dquat
22 −Dquat
23 →Dtern
22 .(14)
S3
Simulation details
Molecular dynamics simulations were performed with the program ms21in two steps: First,
a simulation in the isobaric-isothermal (NpT) ensemble was carried out to calculate the
density at the desired temperature, pressure and composition. In the second step, a canonic
(NV T) ensemble simulation was performed at the corresponding thermodynamic conditions
to simultaneously determine the phenomenological coefficient and thermodynamic factor ma-
trices. 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.93 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 rc= 24.5Å. Lennard-Jones long range interactions were considered using an-
gle averaging.2Electrostatic long-range corrections were approximated by the reaction field
technique with conducting boundary conditions (RF =∞).
The simulations in the NpT ensemble were equilibrated over 4×105time steps, followed
by a production run over 3×106time steps. In the NV T ensemble, the simulations were
equilibrated over 5×105time steps, followed by production runs of 1.4×108time steps. The
phenomenological coefficients were calculated for up to 1.1×106independent time origins
of the autocorrelation functions. The sampling length of the autocorrelation functions was
17.5ps throughout. The separation between the time origins was chosen such that all
autocorrelation functions have decayed at least to 1/eof their normalized value to achieve
their time independence.3The uncertainties of the predicted values were estimated with a
block averaging method.4
S4
Figure S1: Elements of the thermodynamic factor matrix Γij for six quaternary compositions;
a) x1= 0.125 mol mol−1,x2= 0.25 mol mol−1and x3= 0.375 mol mol−1; b) x1= 0.125
mol mol−1,x2= 0.375 mol mol−1and x3= 0.25 mol mol−1; c) x1= 0.25 mol mol−1,x2=
0.125 mol mol−1and x3= 0.375 mol mol−1; d) x1= 0.25 mol mol−1,x2= 0.25 mol mol−1
and x3= 0.25 mol mol−1; e) x1= 0.25 mol mol−1,x2= 0.375 mol mol−1and x3= 0.125 mol
mol−1; f) x1= 0.375 mol mol−1,x2= 0.125 mol mol−1and x3= 0.25 mol mol−1. Present
quaternary results based on Kirkwood-Buff integration (red symbols) are compared with Γij
calculated with the Wilson excess Gibbs energy model fitted to binary simulation results
(blue symbols). The statistical uncertainties are within symbol size.
S5
Figure S2: Main and cross elements of the Fick diffusion coefficient matrix as a function of
water mole fraction for selected compositions with x2= 0.125 mol mol−1and x4= 0.25 mol
mol−1. Simulation results for the quaternary mixture (blue symbols) are shown together with
the results of the diffusion matrix of the ternary subsystem where x1→0(green symbol)
and the predictive equations by Allie-Ebrahim et al.5(red symbols).
S6
Figure S3: Main and cross elements of the Fick diffusion coefficient matrix as a function of
water mole fraction for selected compositions with x3= 0.125 mol mol−1and x4= 0.25 mol
mol−1. Simulation results for the quaternary mixture (blue symbols) are shown together with
the results of the diffusion matrix of the ternary subsystem where x1→0(green symbol)
and the predictive equations by Allie-Ebrahim et al.5(red symbols).
S7
Figure S4: Main and cross elements of the Fick diffusion coefficient matrix as a function
of methanol mole fraction for selected compositions with x3= 0.125 mol mol−1and x4=
0.25 mol mol−1. Simulation results for the quaternary mixture (blue symbols) are shown
together with the results of the diffusion matrix of the ternary subsystem where x2→0
(green symbol) and the predictive equations by Allie-Ebrahim et al.5(red symbols).
S8
Figure S5: Main and cross elements of the Fick diffusion coefficient matrix as a function
of ethanol mole fraction for selected compositions with x1= 0.125 mol mol−1and x4=
0.25 mol mol−1. Simulation results for the quaternary mixture (blue symbols) are shown
together with the results of the diffusion matrix of the ternary subsystem where x3→0
(green symbol) and the predictive equations by Allie-Ebrahim et al.5(red symbols).
S9
Figure S6: Main and cross elements of the Fick diffusion coefficient matrix as a function
of ethanol mole fraction for selected compositions with x2= 0.125 mol mol−1and x4=
0.25 mol mol−1. Simulation results for the quaternary mixture (blue symbols) are shown
together with the results of the diffusion matrix of the ternary subsystem where x3→0
(green symbol) and the predictive equations by Allie-Ebrahim et al.5(red symbols).
S10
0
2
4
6
8
10
12
2 4 6 8
g11 (r)
r/ Å
x1
Figure S7: Water-water radial distribution function g11(r)for varying water mole fraction x1
= 0.125 mol mol−1(black line), 0.25 mol mol−1(blue line), 0.375 mol mol−1(red line) and
0.5 mol mol−1(green line) with constant x2= 0.125 mol mol−1and x4= 0.25 mol mol−1.
S11
0
1
2
3
4
5
2 4 6 8
g12 (r)
r/ Å
x1
Figure S8: Water-methanol radial distribution function g12(r)for varying water mole fraction
x1= 0.125 mol mol−1(black line), 0.25 mol mol−1(blue line), 0.375 mol mol−1(red line) and
0.5 mol mol−1(green line) with constant x2= 0.125 mol mol−1and x4= 0.25 mol mol−1.
S12
0
1
2
3
4
5
2 4 6 8
g12 (r)
r/ Å
x1
Figure S9: Water-methanol radial distribution function g12(r)for varying water mole fraction
x1= 0.125 mol mol−1(black line), 0.25 mol mol−1(blue line), 0.375 mol mol−1(red line) and
0.5 mol mol−1(green line) with constant x3= 0.125 mol mol−1and x4= 0.25 mol mol−1.
S13
0
0.5
1
1.5
2
2 4 6 8 10 12 14
g23 (r)
r/ Å
x2
Figure S10: Methanol-ethanol radial distribution function g23(r)for varying methanol mole
fraction x2= 0.125 mol mol−1(black line), 0.25 mol mol−1(blue line), 0.375 mol mol−1(red
line) and 0.5 mol mol−1(green line) with constant x1= 0.125 mol mol−1and x4= 0.25 mol
mol−1.
S14
References
(1) Glass, C. W.; Reiser, S.; Rutkai, G.; Deublein, S.; Köster, A.; Guevara-Carrion, G.;
Wafai, A.; Horsch, M.; Bernreuther, M.; Windmann, T. et al. ms2: A Molecular Simula-
tion Tool for Thermodynamic Properties, New Version Release. Comp. Phys. Commun.
2014,185, 3302–3306.
(2) Lustig, R. Angle-average for the powers of the distance between two separated vectors.
Mol. Phys. 1988,65, 175–179.
(3) Schoen, M.; Hoheisel, C. The mutual diffusion coefficient D_12 in binary liquid model
mixtures. Molecular dynamics calculations based on Lennard-Jones (12-6) potentials.
Mol. Phys. 1984,52, 33–56.
(4) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford,
1987.
(5) Allie-Ebrahim, T.; Russo, V.; Ortona, O.; Paduano, L.; Tesser, R.; Di Serio, M.;
Singh, P.; Zhu, Q.; Moggridge, G. D.; D’Agostino, C. A Predictive Model for the Dif-
fusion of a Highly Non-Ideal Ternary System. Phys. Chem. Chem. Phys. 2018,20,
18436–18446.
S15
