Fourier Correlation Method for Simulating Mutual Diffusion

Article pubs.acs.org/IECR

Fourier Correlation Method for Simulating Mutual Diffusion Coefficients in Condensed Systems at Equilibrium Joseph W. Nichols and Dean R. Wheeler* Department of Chemical Engineering, Brigham Young University, Provo, Utah 84602, United States ABSTRACT: We introduce here the Fourier correlation (FC) method for obtaining bulk diffusion coefficients in simulations of condensed, nonideal mixtures. The method in particular allows a Fickian mutual diffusion coefficient to be calculated directly in a molecular dynamics simulation and provides for a correction of errors due to finite system sizes. It is an equilibrium method that may be performed at the same time as standard thermodynamic computations. The method is tested by simulating a binary system of LJ particles meant to approximate a nonideal liquid mixture of methane and carbon tetrafluoride. Results are validated by comparison to Green−Kubo and nonequilibrium methods and simulation at multiple system sizes.

■

INTRODUCTION The diffusion coefficient is necessary to describe mass transport in nonequilibrium systems. Diffusion coefficients for many systems of interest have been measured experimentally.1−5 Because diffusivity is an expression of an empirical transport law, it potentially depends on all local conditions including composition and thermodynamic state point. Therefore, the applicability of diffusivity values is generally confined to a relatively narrow range of conditions close to the experimental conditions. Moreover, for some systems, measurement of diffusion coefficients is difficult or impractical. For these reasons it is desirable to find a reliable theoretical method for predicting diffusion coefficients and associated trends. There are systems, such as gas mixtures and point defects in crystal lattices, for which statistical mechanical theories can be used to determine diffusivities with reasonable accuracy.6−8 Generally these approaches require a system for which entropic effects (i.e., free energy, free volume, and many-body correlations) are tractable due to the relatively low concentration of the mobile or diffusing species. However, there are many condensed-matter systems of interest for which direct theoretical derivation of mass-transport properties is not likely to be accurate. It is natural in these cases to turn to molecular dynamics (MD) and similar brute-force approaches to simulate diffusion. In this work we introduce the Fourier correlation (FC) method for obtaining bulk diffusion coefficients in molecular dynamics simulations of mixtures. The method in particular allows a Fickian mutual diffusion coefficient to be calculated directly in an MD simulation and provides for a correction of errors due to finite system sizes. It is an equilibrium method that may be performed at the same time as standard thermodynamic computations. The method is tested by simulating a binary system of LJ particles meant to approximate a nonideal liquid-phase mixture of methane and carbon tetrafluoride. Results are validated by comparison to concurrent Green−Kubo (GK) calculations as well as separate nonequilibrium diffusion calculations. Indeed, the FC method can be considered complementary with the GK method and both can be computed simultaneously. The results show that the FC © XXXX American Chemical Society

method is an accurate and efficient route to obtaining Fickian diffusivities.

■

DIFFUSION IN A BINARY MIXTURE The term diffusion in this paper refers specifically to bulk molecular diffusion, or the diffusion of one species through another as a result of species chemical potential gradients, as opposed to self-diffusion or tracer diffusion, which refers to the diffusion of a single particle due purely to Brownian motion. Mutual diffusion is the special case of bulk diffusion in which only two species are present in the system. To describe mass transport in a mixture, one must choose a reference velocity. The most common choices are mass-, molar-, and volume-average reference velocities. As described by Cussler, the volume-average one is particularly useful as it can easily be reduced to simpler cases of constant-density or constant-total-concentration mixtures.9 The volume-average case will serve as our starting point for a macroscopic description of mass transport. Ni, species i total molar flux, can be expressed as a combination of diffusive and advective fluxes: Ni = J□ + ci v □ i

(1)

where ci is species molar concentration. The diffusive flux relative to v□ for either species 1 or 2 in a binary mixture is9

J□ = −D12∇ci i

(2)

where D12 is a Fickian mutual diffusivity. Finally, the volumeaverage velocity is v □ = V1̅ N1 + V2̅ N2

(3)

where V̅ i is species partial volume. Note that the thermodynamics of mixtures stipulates that c1V1̅ + c 2V2̅ = 1

(4)

Received: August 3, 2015 Revised: November 12, 2015 Accepted: November 17, 2015

A

DOI: 10.1021/acs.iecr.5b02849 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

use MD simulations to predict bulk diffusion coefficients have been ongoing since the 1970s.15 Early efforts were generally limited to binary mixtures of Lennard-Jones fluids.16−22 Simulations of penetrant diffusion in membranes have attracted some attention.23−27 While there has been some work simulating systems of three or more species, relatively few have attempted to find a full set of diffusion coefficients for these systems.28−34 Figure 1 gives an overview of conceivable methods for determining microscopic transport coefficients using molecular

For the alternative case of molar-average reference velocity v*, the corresponding flux expressions are10 Ni = J*i + ci v*

(5)

J*i = −D12c∇yi

(6)

v* =

1 (N1 + N2) c

(7)

where c is the total concentration and yi = ci/c is the number or mole fraction. Both eqs 2 and 6 have been called Fick’s first law of diffusion. For the binary system the two types of diffusive fluxes are related by J1□ = cV2̅ J1* J□ = cV1̅ J*2 2

(8)

which can be proved with a combination of the above equations.9 The two types of fluxes are exactly equivalent (J□ i = Ji*) when c is constant, which is true for mixtures of ideal gases or in the limit of an infinitely dilute solute (e.g., c1 ≪ c2). Likewise under this condition, v□ = v*. Mass transport is alternatively described by the Maxwell− Stefan transport equation, generalized to use a thermodynamic driving force.5,9−12 In the case of a binary mixture, the Maxwell−Stefan equation reduces to eq 5, provided we define diffusive flux as c J*i = −+12 i ∇μi kBT (9)

Figure 1. Matrix of microscopic transport property simulation methods. Asterisks indicate methods compared in this work.

dynamics. Grouped to the left are nonequilibrium molecular dynamics (NEMD) simulation methods, which depend on the application of an external driving force or boundary condition to create transport events that then can be observed.27,35 The rightmost column gives the equilibrium methods that require no external fields or boundary conditions but rather track molecular transport under spontaneous thermal fluctuations. Equilibrium methods have the advantage of allowing the calculation of multiple system properties in a single simulation as well as eliminating the need for the exceptionally large driving forces which are generally required for NEMD systems. On the other hand, NEMD methods frequently have greater efficiency in calculating a single transport property of interest.31 The rows in Figure 1 further categorize methods: methods on the first row are derived only from a linear constitutive equation. Methods on the second row are additionally and explicitly built upon an equation of change,10 that is, a transient-type macroscopic conservation law. First-row methods are much more often realized.27,30,35−39 Indeed, it is the primary purpose of this work to introduce the Fourier correlation equilibrium method in the second row. While Figure 1 is meant to illustrate general principles, the focus of this work is mass transfer. Accordingly, when applied to mass transport, the methods where concentration variations are explicitly generated or tracked generally yield Fickian diffusivity (D12), while the remaining methods generally yield activity diffusivity (+12 ). Experimental work to determine diffusivity generally uses a concentration gradient as a driving force and therefore reports Fickian diffusivity rather than activity diffusivity.10 Therefore, a simulation method such as the Green−Kubo (GK) method that gives activity diffusivity can only be compared to experiment after applying the thermodynamic factor.14,20 In contrast, the proposed FC method calculates Fickian diffusivity directly. We note that a matrix of mass-transport simulation methods similar to Figure 1 was produced by Liu et al., though they did not anticipate the FC equilibrium method.34

and recognize that Ni ≡ civi. +12 is a binary interaction parameter or activity-based diffusivity, kB is Boltzmann’s constant, T is the absolute temperature, and μi is the chemical potential. Because eq 9 has a more thermodynamically rigorous driving force than Fick’s law, it is conjectured that +12 will be less composition-dependent than D12, though in practice this is not necessarily true.9,13 Note that, because we are dealing with microscopic systems, Ni, Ji, ci, and μi are given here in terms of molecules, rather than moles, with the understanding that conversion between the two unit systems involves the trivial application of Avogadro’s number. By enforcing the equality between eqs 6 and 9, one can show that

D12 = Q 11+12

(10)

where Q11 is a thermodynamic factor given by Q 11 =

⎛ ∂ln γ ⎞ y1 ⎛ ∂μ1 ⎞ 1⎟ ⎜⎜ ⎟⎟ = 1 + ⎜⎜ ⎟ kBT ⎝ ∂y1 ⎠ ln y ∂ ⎝ 1 ⎠T , P T ,P

(11)

where γ1 is an activity coefficient and the derivatives are taken at constant temperature and pressure.14 From this it is obvious that, for an ideal (Lewis) binary mixture in which Q11 = 1 at all compositions, the Fickian and activity diffusivities are the same. For a nonideal mixture Q11 = 1 at the composition end points (y1 = 0,1), so that +12 → D12 in the limit of infinite dilution of one species.

■

OVERVIEW OF SIMULATION METHODS Molecular dynamics simulations are an important tool for calculating diffusivity for model systems and relating macroscopic transport parameters to microscopic behavior. Efforts to B


Article


■

THEORY In this section we develop the mathematical basis of the Fourier correlation method for a binary mixture. In addition we summarize some previously developed simulation methods, including the Green−Kubo formula for mutual diffusivity, because they serve as a means of comparison and validation. Both FC and GK methods depend on the fact that, at equilibrium, random thermal fluctuations create local concentration or chemical potential gradients. These local fluctuations decay toward an equilibrium average by the same processes as do macroscopic nonequilibrium states, in the limit of small perturbations (i.e., linear processes). This important idea is the essence of Onsager’s regression hypothesis and is the basis for the fluctuation−dissipation theorem.40−42 Despite the long history and now-routine use of this idea, it is rather remarkable to consider that one can extract a multitude of macroscopic nonequilibrium properties, including diffusion coefficients, from an equilibrium simulation of around a thousand molecules. Development of FC Method. The origin of the Fourier correlation method can be traced to a series of two papers intended to prove Onsager’s reciprocal relation for multicomponent mass transfer.43,44 In developing the proof, it became apparent that similar analysis could be used to compute a Fick’s law diffusivity in an equilibrium simulation. As indicated above, the FC method requires the use of a macroscopic conservation law, namely species continuity:

composition. Next we substitute these expansions into eq 12 and keep only first-order terms as part of the linearization process. For instance, the terms like ∇ci·v□ and ∇ci·∇V̅ i are eliminated because they are products of first-order terms and are therefore second-order. This gives the following governing equation: ∂ci(1) = D12∇2 ci(1) ∂t

also known as Fick’s second law. To simplify the notation here and below, we take D12 to be the value at the equilibrium mixture composition. Nevertheless, for the microscopic systems simulated here, D12 remains a function of wavenumber q as noted below. (The same is true for mixture quantities Q11 and +12 .) A finite spatial Fourier transform of the governing equation turns it into a simpler ordinary differential equation (ODE). Anticipating this course, we generate the structure factor for species i, which is the Fourier transform of the species concentration:

(12)

ci(1)(r, t ) =

(13)

(14)

A combination of the above equations results in ∇·v □ = −D12(∇c1·∇V1̅ + ∇c 2·∇V2̅ )

(15)

We now consider a system at equilibrium, meaning on average that species fluxes are zero. In the thermodynamic (large-system) limit, thermal fluctuations from the average state are vanishingly small. Linearization of the mixture properties for a finite system about this equilibrium state yields

ci =

Ni + ci(1) + O(2) V

v □ = 0 + v □(1) + O(2)

a

a∈i

(19)

N V

∑ ψi(q, t )eiq·r q≠0

(20)

Vector q is known as the reciprocal lattice vector or wave vector and functions as an eigenvalue for the governing equation. In order to satisfy periodic boundary conditions in a cubic system, q = 2πn/L, where L is the box length and each of the three elements of n take integer values. Note that n2 = n·n, also an integer, gets used below as a shorthand for indicating the discrete values of q = |q| used in the FC algorithm. In order for a macroscopic quantity to be obtained in a molecular simulation, it must be averaged over a long trajectory, here indicated by angle brackets. For the assumed case of a single-phase homogeneous equilibrium simulation, ⟨ψi(q, t)⟩ = 0 for q ≠ 0, and ⟨c(1) i (r, t)⟩ = 0. These expressions follow from the fact that thermal disorder causes each location, over sufficient time, to revert to the macroscopic equilibrium or average state. This is a manifestation of the ergodic hypothesis of statistical mechanics. We are now ready to solve the governing diffusion equation (eq 18) by the substitution of a Fourier series for the dependent variable (eq 20), a traditional method for solving partial differential equations. This yields an equation in which each of the terms is multiplied by an eigenfunction, namely eiq·r. By multiplying by an orthogonal eigenfunction and integrating over the volume, the sums over q are eliminated and one is left with an ODE for the structure factors at any particular nonzero q value:

Combining this with the total differential of eq 4, one can show that V1̅ dc1 + V2̅ dc 2 = 0

∑ e−iq·r (t)

where i is the imaginary number and N is the total number of molecules in the system. The sum is over all molecules a belonging to species i and requires molecule center-of-mass positions ra(t). The linear concentration field is then constructed as a Fourier series:

where t is time. The second line comes from substituting eqs 1 and 2. (Because eqs 1 and 5 are mathematically equivalent expressions for Ni, it is ultimately immaterial which expression one chooses to substitute.) In order to get a usable expression for ∇·v□ we first introduce the Gibbs−Duhem equation for volume (at constant T and P) c1dV1̅ + c 2dV2̅ = 0

1 N

ψi(q, t ) =

∂ci = −∇·Ni ∂t = ∇·(D12∇ci) − ∇ci·v □ − ci∇·v □

(18)

(16) (17)

where Ni and V are species number and system volume, respectively. The first term in each series is the equilibrium average (a constant), and the second term is a first-order fluctuation that depends on time and position. We can likewise expand D12 and V̅ i into series with constant leading terms, followed by fluctuation terms that depend linearly on the local

dψi(q, t ) dt

= −D12(q)q2ψi(q, t )

(21)

Integrating this ODE gives C


Article

Industrial & Engineering Chemistry Research ψi(q, t ) = ψi(q, 0) exp[−D12(q)q2t ]

One can also think of eq 24 as describing a virtual fluctuation experiment in which the simulation box is divided into two regions, each of which is not necessarily contiguous. Diffusivity is related to the fluctuating rate of spontaneous mass transfer between the two regions. As shown in Figure 2, vector q

(22)

This solution would be appropriate for a transient molecular dynamics simulation in which the system was prepared with a concentration inhomogeneity and allowed to relax toward equilibrium. However, an ensemble-average of the two sides of eq 22 in a purely equilibrium simulation leads to a trivial result, 0 = 0. Instead, we transform each side into an autocorrelation quantity through multiplication by the complex conjugate of ψi(q, 0) before performing the average:

Figure 2. Representations of two reciprocal lattice vectors, showing division of the simulation cell into alternating regions that exchange mass.

⟨ψi(q, t )ψi( −q, 0)⟩ = ⟨ψi(q, 0)ψi( −q, 0)⟩ exp[−D12(q)q2t ]

(23)

This can be rewritten as 2

Sii(q, t ) = Sii(q, 0) exp[−D12(q)q t ]

represents how the box is divided in each coordinate direction into rectangular slabs that alternate between the two regions. At larger q values the alternating regions are more closely spaced and diffusive relaxation takes place on shorter, more accessible time scales. On the other hand, smaller q values better represent desired macroscopic diffusion. A nonlinear least-squares fit was used to regress the best value of D12 at each discrete q value by minimizing the sum of squared residuals between a set of simulated S11 relaxation points and eq 24. Figure 3 illustrates the results for one such fit,

(24)

where, to generalize for arbitrary species i and j, Sij is defined as Sij(q, t ) = N ⟨ψi(q, t )ψj( −q, 0)⟩

(25)

Sij(q, t), which we term a Fourier pairwise time-correlation function, plays a central role in the FC method. Sij as defined is real-valued, intensive (not explicitly dependent on system size), and readily calculable in a molecular simulation from knowledge of molecular positions in time. Equation 24 describes how microscopic thermal fluctuations decay, on average, toward zero just as would macroscopic nonequilibrium disturbances, in accord with Onsager’s regression hypothesis. It should be obvious that one can use eq 24 as the basis for regressing D12(q) for any nonzero q value. In practice one seeks a macroscopic value of D12, which corresponds to the thermodynamic limit, namely q → 0 or L → ∞. This means that a range of finite q values are used and one must extrapolate to the q = 0 result, as discussed more fully below. All other things being equal, the larger the box size L, the smaller the nonzero values of q that can be simulated, leading to a more reliable extrapolation. The algorithmic details for collecting Sij(q, t) in an MD computer code are virtually identical to how other timecorrelation functions (see eqs 26 and 27) are collected.32 For instance, a routine technique to maximize data collection is the use of multiple time origins.45 In this case, multiple instances of eq 25 are collected concurrently at discrete times t and with staggered time origins. These multiple instances are then averaged. All times are necessarily integer multiples of the simulation time step Δt. We used two other ideas in this work to maximize data collection. The first relies on the fact that the liquid is isotropic, and so Sij(q, t) should be the same for multiple q values that have the same scalar magnitude q. Furthermore, Sij(q, t) and Sij(−q, t) are trivially the same and so there is no need to collect both. Thus, we collected Sij values in a hemispherical fashion in q space and reduced them to a radially symmetric function Sij(q, t) evaluated at discrete q and t values. The discrete eigenvalues q chosen depend on box size L and a set of integer n2 values (see discussion following eq 20). Here we collected Sij for q values in the interval corresponding to n2 = 1−17. The second idea relies on the fact that the S22 signal is equivalent to S11 in terms of time decay (see eq 24). The final diffusivities for the FC method were therefore determined from an average of results for both the S11 and S22 signals, though for convenience we focus on S11 in the discussion below.

Figure 3. Representative Fourier correlation function versus time. The semilog plot linearizes the simulation and fitting curves (eq 24) to illustrate the adequacy of the method, though the fitting is not done in this semilog space. The fitting process produces a Fickian diffusion coefficient for one simulation run at one q value.

showing good agreement. The fitting process is repeated for other simulation runs and q values in order to produce D12 averages and error bars shown in Figures 8 and 9. The obvious departure of the simulation from the fitting curve at large times is an indication of the relative uncertainty or noise of the S11 signal, which in this case is about 1% of the zero-time S11 value. In any case, the relative noisiness of the signal at large times does not substantially affect the nonlinear fit because the exponentially decaying S11 signal is already close to zero at that point. We note that the fitting of D12 can in theory be complicated by a Gaussian-shaped “inertial” or “ballistic” or “Knudsen” region in S11 near the time origin. This region represents the brief period before a sufficient number of collisions between molecules leads to proper diffusive dynamics. For condensedphase systems at reasonably low q values, we observed that this effect concludes at such short times that it has essentially no effect on fitting results. Indeed, the ballistic region, which in Figure 3 would appear as an initial portion of the curve shaped like a parabola and having zero slope at the time origin, is not visible on this time scale. On the other hand, for low-density gas systems (not shown) the ballistic region extends to relatively long times and can confound fitting of results, as eq D


Article


extrapolation to infinite time. In this work we used the GK form and did not consider the Einstein formula separately. As a practical matter, one can carry out the integration in eq 26 only to a finite time. Therefore, an inherent difficulty with any Green−Kubo relation is how to deal with the neglected portion of the integral, known as the long-time tail. Even if longer times are accessible in a simulation, it is not desirable to carry the integral past the point where random noise in the integrand (due to finite-size simulation sampling) is larger in magnitude than the underlying signal. Integrating too long over a noisy correlation function will simply produce a “random walk” that increases uncertainty in the resulting transport property.32 On the other hand, simple truncation of the longtime tail results in a systematic bias to the transport property.46 Alternatively, if the long-time tail can be represented by a simple semiempirical function, the resulting expression can be used to extrapolate the GK integral out to infinite time. In the case of viscosity, mode coupling theory predicts a positive longtime tail with the analytic form of At−3/2 .47 By analogy, we apply the same form here, where A is an adjustable fitting parameter. An example of the application of the long-time tail correction is shown in Figure 4. This is the result for one representative simulation run; to get the averages and error bars on +12 in Figures 8 and 9, the process is repeated. For our system the extrapolation augments the truncated GK integral by around 5−10% and brings the final average value of the activity diffusivity much closer to the value given by a nonequilibrium method. Finally, we note that Maginn and co-workers have recently developed an improved GK fitting and time-extrapolation procedure in the context of viscosity calculations.48 Thermodynamic Factor. As eq 10 shows, a thermodynamic factor is required to convert between Fickian and activity mutual diffusivities. We previously showed how to obtain efficiently the thermodynamic factor from a periodic molecular simulation using an improved implementation of Kirkwood− Buff solution theory.14 For a binary mixture the relevant formula is

24 accounts only for long-time diffusive behavior. Therefore, the FC method may be one of the few methods where it is more convenient to predict a transport property for a liquid rather than for a gas. Green−Kubo Method. It is instructive to contrast the FC method with the traditional means of calculating transport properties in an equilibrium simulation, namely the Green− Kubo method. The GK formula for mutual diffusion is32 +12 =

N1N2 3N

∫0

∞

dt ⟨[v2(t ) − v1(t )]·[v2(0) − v1(0)]⟩ (26)

where the integrand is a species velocity correlation function (VCF) that depends on species-average velocities vi(t). This formula can be shown to be identical to Jaccuci and McDonald’s original derivation when combined with a statement of system momentum conservation, m1N1v1 + m2N2v2 = 0, where mi is molecular mass.18 Figure 4 shows an

Figure 4. Example of a running GK integral for activity diffusivity, including the application of a long-time-tail correction. The equation At−3/2 was fit to the integrand between 2 and 8 ps (shown by the overlap between the two curves). After 8 ps variation in the running integral was determined solely by the fitted equation.

Q 11(q) =

example of a running GK integral for activity diffusivity obtained from a single simulation. As with the FC method, multiple time origins are used here to increase sampling efficiency. The so-called Einstein or mean-square-displacement formula for mutual diffusion may be readily obtained from eq 26: NN d +12 = 1 2 lim ⟨|[R 2(t ) − R1(t )] − [R 2(0) − R1(0)]|2 ⟩ 6N t →∞ dt

N1N2 N22S11(q ,

0) − 2N1N2S12(q , 0) + N12S22(q , 0) (28)

This means that quantities already required by the FC method, namely radially averaged Sij(q, t) functions, contain all the information required to make the diffusivity conversion. For a nonideal mixture, experience suggests Q11 is strongly dependent on q for the typically accessible range of q values, and so some care must be used in performing the extrapolation to the thermodynamic limit. This matter receives considerable attention below. Nonequilibrium External-Field Method. As a check on the FC and GK equilibrium methods employed here, we also performed NEMD simulations of activity diffusivity using the method of a homogeneous external field. The general approach and algorithmic details are given in prior work.31 As NEMD methods generally require extremely large perturbations in order to measure a transport property, nonlinear effects become important. Therefore, one must extrapolate results back to the linear regime, effectively the zero-perturbation limit. In the present work, each NEMD diffusivity point in Figures 8 and 9 was based on an extrapolation from multiple simulations at different field strengths as illustrated in Figure 5. As shown,

(27)

where Ri(t) is the time-varying aggregate center-of-mass position for species i, and must be based on a separate set of “unfolded” molecule positions where periodic boundaries have not been applied.32 An Einstein formula merely uses particle equations of motion (built into the MD algorithm) to twiceintegrate the GK integrand, thus requiring a differentiation in the end to recover the transport property. As long as such integrations and differentiations are performed with suitable accuracy, the Einstein and GK representations contain the same information and should give identical transport-property results. Nevertheless some researchers prefer one form over the other for descriptive or analysis purposes, or ease of E


Article


The code automatically adjusts the size of the time step with changing composition to maintain numerical stability while maximizing the sampling of phase space. To this end, Δt is chosen to maintain a particular root-mean-square displacement per time step for the aggregate population of molecules in the simulation, with this average displacement empirically selected to be 1.5 fm. For y2 = 0.6 this corresponds to Δt = 3.8 fs. All simulations were equilibrated for 2 × 106 time steps, followed by sequential production runs of 5 × 106 time steps. Error bars in the figures below represent 95% confidence intervals based on 3−6 sequential NVT runs at the same point. As shown below, Green−Kubo results were confirmed using steady-state field-driven NEMD at the composition point y2 = 0.6 and multiple system sizes. As with the equilibrium simulations, NEMD simulations were run for 5 × 106 time steps following equilibration. Body forces of 1, 2, 3, 5, and 7 MeV/m were applied to the methane particles, and an opposing force was applied to the CF4 particles to produce net zero system acceleration. The activity diffusivity +12 was calculated for each field value, and extrapolation was performed as described for Figure 5.

Figure 5. Extrapolation on a set of NEMD simulation results to find the activity diffusivity at zero field strength. Error bars are 95% confidence intervals for a set of three simulations of N = 500 particles at each field strength.

uncertainties increase as field strengths are decreased due to a smaller signal-to-noise ratio. The diffusivity extrapolation was based on a weighted linear least-squares fit that takes into account the uncertainty at each field strength.

■

SIMULATION RESULTS AND DISCUSSION This section compares the results for the respective simulation methods. Figures 6 through 8 are detailed looks at the leastideal system composition, y2 = 0.6, which provides a robust measure of the suitability of the proposed FC method and extrapolation techniques. In fact, this state point is at a temperature just above a liquid−liquid miscibility gap, which

■

SIMULATION DETAILS The test platform was a binary liquid mixture of Lennard-Jones particles. The LJ parameters, given in Table 1, were chosen to Table 1. Lennard-Jones Parameters for a Binary Mixture of CF4 (Species 1) and CH4 (Species 2) Taken from Nichols et al.14 species pair σ (nm) ϵ/kB (K)

11

22

12

0.4150 175.0

0.3728 149.0

0.3939 142.1

approximate carbon tetrafluoride (CF4) and methane (CH4). These parameters, including the cross interactions, have been shown to approximate reasonably well the nonideal liquid-state behavior of this mixture, including the presence of an upper solution critical point around Tcrit = 94.5 K and y2,crit = 0.57.14,49 The molecular dynamics computer code used to carry out the simulations has been described in prior work.14,31,32 The GK activity diffusivity, FC Fickian diffusivity, and thermodynamic factor were each calculated from the same set of equilibrium simulations. Simulation conditions were chosen to obtain a homogeneous but nonideal liquid mixture over the entire composition range. All simulations were run at a temperature of 150 K. Preliminary constant-NPT runs at each composition point were performed at pressure 5 MPa in order to obtain appropriate density values. Subsequent constant-NVT production runs were then based on these densities. Compositions were tested at the following mole fractions of methane: y2 = 0.1, 0.2, 0.4, 0.6, 0.8, and 0.9. System size N = 1200 was generally used for the six composition points. At the composition point y2 = 0.6, the least ideal of those tested, additional simulations at system sizes N = 500, 700, 1000, 1400, 2000, 4000 particles were used to explore system-size effects.

Figure 6. Representative subset of FC results as a function of wavenumber for four system sizes: N = 500, 1200, 2000, and 4000. Additional runs at these and other system sizes were performed (see Figures 7 and 8) but are not shown here. The activity diffusivity points (+12 ) were generated by dividing the Fickian diffusivity points (D12) by the corresponding thermodynamic factors (Q11). The two types of curves correspond to linear (Lin) and quasilinear (QLin) fitting methods, applied to N = 2000 data only. F


Article

Industrial & Engineering Chemistry Research occurs when Q11 ≤ 0.14 For the more-ideal simulation points in this work, and for more-ideal systems generally, the differences between the extrapolation techniques, and between D12 and +12 , are accordingly diminished. Figure 6 shows results for the FC method for four different system sizes at y2 = 0.6 in order to illustrate the process of extrapolating to zero wavenumber. For each simulation we calculated average Fickian diffusivity D12 and thermodynamic factor Q11 values for the 15 discrete q values in the range 1 ≤ n2 ≤ 17. The corresponding activity diffusivity, +12 = D12 /Q 11, is also shown for each point. As previously indicated, the calculation of the thermodynamic factor is embedded in the Fourier correlation method and can be obtained at essentially no additional cost. The data for the three indicated properties show two notable features: nearly linear behavior for q < 5 nm−1 and significant overlap between data sets from differentsize simulations. These features give confidence that it is possible to determine accurate macroscopic diffusivities from relatively small simulations, by extrapolating to the thermodynamic (q → 0 or L → ∞) limit. The first extrapolation procedure was a simple linear leastsquares fit of points for q < 5 nm−1. The extrapolations to q = 0 are shown by the solid lines in Figure 6, based only on data from the N = 2000 simulation. Theoretical considerations show that structure-factor-based functions of q must be even, meaning that these functions have zero slope at q = 0.14 This suggests that a simple linear extrapolation of Q11 or D12 could underpredict the true intercepts (macroscopic values). We therefore developed a second extrapolation technique termed quasilinear or QLin. The Q11 values for wave numbers below 5 nm−1 were leastsquares-fit to the following empirical model that both follows the observed linear trend and satisfies the zero-slope condition as q → 0: Q 11 =

Aq2 q + 1nm−1

Figure 7. Extrapolated values of the thermodynamic factor at multiple simulation system sizes for y2 = 0.6, the least ideal composition point. Error bars represent 95% confidence intervals.

other. Rather, the QLin method in this study generated Q11 values that made for modestly better agreement between the GK and FC D12 values shown in Figures 8 and 9. Figure 7 shows that the preferred QLin method is relatively sizeindependent for N ≥ 1000, again suggesting that with an appropriate extrapolation to q = 0 one can obtain reliable results with a relatively small simulation size. Figure 8 compares the NEMD, Green−Kubo, and Fourier correlation methods in calculating diffusivity at different system

+B (29)

where A and B are adjustable fitting parameters as one would encounter in a linear fit. Corresponding QLin fits were used for D12 and +12 as shown in Figure 6. The implied constant of 1 nm−1 in eq 29 represents the wavenumber below which the zero-slope condition begins to have a noticeable effect. It can be adjusted empirically, but the present value worked well over the full range of system points examined here. As implied in Figure 6, the Fourier correlation method can be used to calculate the activity diffusivity by dividing the directly calculated Fickian diffusivity by the thermodynamic factor. While this method gives reasonable values for the activity diffusivity, it appears to have a greater system-size bias, which could be an indicator of the inadequacy of a linear extrapolation on the D12/Q11 curve. Figure 7 shows the results for extrapolated Q11 values at y2 = 0.6 for the two extrapolation techniques with changing system size. The notable offset between the two methods is caused by the steep slope of Q11 vs q for this nonideal state point. The generally increased uncertainty in extrapolated Q11 values with system size is counterintuitive as greater computational cost has been expended at larger system sizes. This may be due to the fact that smaller periodic systems tend to suppress persistent spatial concentration fluctuations, due to the free energy penalty of concentration gradients.50,51 There is little in Figure 7 itself to suggest one extrapolation method is superior to the

Figure 8. Comparison of diffusivities obtained by different methods and at different system sizes for y2 = 0.6. The upper data set represents activity diffusivities. The lower data set represents Fickian diffusivities. Points are offset horizontally for clarity. Error bars represent 95% confidence intervals.

sizes. Both Fickian and activity diffusivity data are shown for comparison. For the upper data set (open symbols), FC Fickian diffusivities are converted to activity diffusivities with Q11 values using the two extrapolation techniques. NEMD results are included as well. The excellent agreement between NEMD and corresponding GK results gives us confidence that the GK results, including associated long-time-tail corrections, are reliable. For the lower data set of Figure 8 (filled symbols), the Green−Kubo activity diffusivity is modified by an extrapolated Q11 value in order to compare to the FC result. Even though the GK value itself is not explicitly corrected to the thermodynamic limit, the Q11 value is. The resulting product for the QLin extrapolation is in close agreement with the Fickian diffusivity directly calculated by the Fourier correlation method. Which method is preferred to obtain Fickian diffusivities? The observed statistical uncertainties for each method are comparable. Both methods require an extrapolation to q = 0 (in the case of the GK result this is because it must be multiplied G


Article


thus the extrapolation naturally corrects for error due to finite simulation size. In this work preliminary NPT simulations were used to calculate appropriate system densities. This is a relatively easy calculation as densities are one of the more convergent and stable properties to obtain in molecular simulations. Subsequent NVT simulations at these densities were used to perform FC, GK, and NEMD analysis. In the case of the FC method, constant-volume simulations have the advantage that timedependent structure factors need only be collected for a small number of discrete wave numbers. A constant-pressure simulation, in which volume fluctuates, would produce properties for a quasidiscrete or smeared distribution of wave numbers. The wavenumber could still be discretized into bins (i.e., as for a histogram) to determine Sij(q, t) from an NPT simulation. In this manner we believe the FC method is easily adaptable to calculating diffusivity directly in NPT simulations. In subsequent work we intend to show the extension of the FC method to systems of greater than two components in which a Maxwell−Stefan framework is used, 10 as was anticipated by prior theoretical development.44

by Q11). The GK method has the additional uncertainty of a necessary long-time-tail correction. For this reason we believe the FC method is likely to be more reliable in general. Nevertheless, the FC and GK results in Figure 8 were obtained simultaneously from the same set of simulations, and ultimately there is no need to choose one method at the expense of the other. Figure 9 provides a summary of Fickian diffusivities for the simulated CF4−CH4 system at six different compositions. Six

■

AUTHOR INFORMATION

Corresponding Author

Figure 9. Dependence of the Fickian mutual diffusion coefficient on methane mole fraction for three different methods. Error bars represent 95% confidence intervals. Points are offset horizontally for clarity. Lines are guides for the eye.

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■ ■

independent equilibrium simulations (N = 1200, 5 × 10 time steps each) were run for each composition to get the indicated averages and uncertainties. Additional NEMD simulations were run at three compositions. The diffusion coefficients calculated by the Fourier correlation method are in excellent agreement with Green−Kubo and NEMD values that were modified by extrapolated Q11 values. The maximum difference between the FC diffusivity and the modified GK diffusivity is less than 3%. We are not aware of experimental values to which these diffusivities may be compared. 6

ACKNOWLEDGMENTS This work was supported by the U.S. National Science Foundation through Grant CTS-0547610. REFERENCES

(1) McCall, D. W.; Douglass, D. C.; Anderson, E. W. J. Chem. Phys. 1959, 31, 1555. (2) Reis, R. A.; Nobrega, R.; Oliveira, J. V.; Tavares, F. W. Chem. Eng. Sci. 2005, 60, 4581−4592. (3) Verbrugge, M. W.; Koch, B. J.; Schneider, E. W. J. Appl. Electrochem. 2000, 30, 269−275. (4) Tyrrell, H. J. V.; Harris, K. R. Diffusion in Liquids; Butterworths: London, 1984. (5) Rutten, P. W. M. Diffusion in Liquids; Delft University Press: Delft, The Netherlands, 1992. (6) Beijeren, H. V.; Ernst, M. H. Physica 1973, 70, 225−242. (7) Kurochkin, V. I.; Makarenko, S. F.; Tirskii, G. A. J. Appl. Mech. Tech. Phys. 1984, 25, 218−225. (8) Heitjans, P., Kärger, J., Eds.; Diffusion in Condensed Matter: Methods, Materials, Models, 2nd ed.; Springer: Berlin, 2005. (9) Cussler, E. L. Diffusion: Mass Transfer in Fluid Systems, 3rd ed.; Cambridge University Press: New York, 2009. (10) Bird, R. B.; Stewart, W. E.; Lightfoot, E. N. Transport Phenomena, 2nd ed.; Wiley: New York, 2002. (11) Newman, J. Chem. Eng. Sci. 2009, 64, 4796−4803. (12) Hirschfelder, J. O.; Curtiss, C. F.; Bird, R. B. Molecular Theory of Gases and Liquids; Wiley: New York, 1954. (13) Newman, J.; Thomas-Alyea, K. E. Electrochemical Systems, 3rd ed.; Wiley: New York, 2004. (14) Nichols, J. W.; Moore, S. G.; Wheeler, D. R. Phys. Rev. E 2009, 80, 51203. (15) Erdey-Gruz, T. Transport Phenomena in Aqueous Solutions; Akademiai Kiado: Budapest, 1974. (16) Hansen, J.-P.; McDonald, I. R. Theory of Simple Liquids, 2nd ed.; Academic Press: New York, 1986. (17) Hansen, J.-P.; McDonald, I. R. J. Phys. C: Solid State Phys. 1974, 7, L384−L386.

■

CONCLUSION The Fourier correlation method provides an efficient alternative to the Green−Kubo or Einstein methods for calculating diffusion coefficients from equilibrium simulations of condensed-phase systems. The method is based on the use of a finite Fourier transform that is fully consistent with the periodic boundaries typically used in molecular simulations of bulk materials. Agreement between FC, GK, and NEMD methods was satisfactory for a nonideal LJ liquid mixture, once appropriate corrections and extrapolations were used. The FC method generates Fickian diffusivity, while the GK and constant-field NEMD methods generate activity diffusivity. The use of a thermodynamic factor Q11, naturally included in the FC method, allows conversion and comparison between the two diffusivity representations. Some might criticize the FC method because it requires an extrapolation to zero wavenumber. However, it should be noted that other methods likewise require extrapolations. In the case of constant-field NEMD, the extrapolation is to zero field. In the case of the Green−Kubo and Einstein methods, the extrapolation is to infinite time. One practical benefit of the extrapolation in the FC method is that the end point is equivalent to the thermodynamic or macroscopic limit, and H


Article

Industrial & Engineering Chemistry Research (18) Jacucci, G.; McDonald, I. R. Phys. A 1975, 80, 607−625. (19) Jolly, D. L.; Bearman, R. J. Mol. Phys. 1980, 41, 137−147. (20) Hoheisel, C.; Vogelsang, R. Comput. Phys. Rep. 1988, 8, 1−70. (21) Morriss, G. P.; Evans, D. J. Mol. Phys. 1985, 54, 629. (22) Rowley, R. L.; Stoker, J. M.; Giles, N. F. Int. J. Thermophys. 1991, 12, 501−513. (23) Müller-Plathe, F.; Rogers, S. C.; van Gunsteren, W. F. J. Chem. Phys. 1993, 98, 9895−9904. (24) MacElroy, J. M. D. J. Chem. Phys. 1994, 101, 5274. (25) Sunderrajan, S.; Hall, C. K.; Freeman, B. D. J. Chem. Phys. 1997, 107, 10714−10722. (26) Thompson, A. P.; Heffelfinger, G. S. J. Chem. Phys. 1999, 110, 10693−10706. (27) Arya, G.; Chang, H.-C.; Maginn, E. J. J. Chem. Phys. 2001, 115, 8112. (28) Payne, R.; Theodorou, I. E. J. Phys. Chem. 1972, 76, 2892−2900. (29) Soetens, J.-C.; Millot, C.; Maigret, B. J. Phys. Chem. A 1998, 102, 1055−1061. (30) Newman, J.; Thomas, K. E.; Hafezi, H.; Wheeler, D. R. J. Power Sources 2003, 119, 838−843. (31) Wheeler, D. R.; Newman, J. J. Phys. Chem. B 2004, 108, 18362− 18367. (32) Wheeler, D. R.; Newman, J. J. Phys. Chem. B 2004, 108, 18353− 18361. (33) Liu, X.; Martin-Calvo, A.; McGarrity, E.; Schnell, S. K.; Calero, S.; Simon, J.-M.; Bedeaux, D.; Kjelstrup, S.; Bardow, A.; Vlugt, T. J. H. Ind. Eng. Chem. Res. 2012, 51, 10247−10258. (34) Liu, X.; Schnell, S. K.; Simon, J.-M.; Krüger, P.; Bedeaux, D.; Kjelstrup, S.; Bardow, A.; Vlugt, T. J. H. Int. J. Thermophys. 2013, 34, 1169−1196. (35) Morriss, G. P.; Evans, D. J. Phys. Rev. A: At., Mol., Opt. Phys. 1987, 35, 792−797. (36) Thomas, J. C.; Rowley, R. L. J. Chem. Phys. 2007, 127, 174510. (37) Wheeler, D. R.; Fuller, N. G.; Rowley, R. L. Mol. Phys. 1997, 92, 55−62. (38) Stoker, J. M.; Rowley, R. L. J. Chem. Phys. 1989, 91, 3670−3676. (39) Sarman, S.; Evans, D. J. Phys. Rev. A: At., Mol., Opt. Phys. 1992, 45, 2370−2379. (40) Onsager, L. Phys. Rev. 1931, 37, 405−426. (41) Onsager, L. Phys. Rev. 1931, 38, 2265−2279. (42) Callen, H. B.; Welton, T. A. Phys. Rev. 1951, 83, 34−40. (43) Monroe, C. W.; Newman, J. Ind. Eng. Chem. Res. 2006, 45, 5361−5367. (44) Monroe, C. W.; Wheeler, D. R.; Newman, J. Ind. Eng. Chem. Res. 2015, 54, 4460−4467. (45) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, U.K., 1990. (46) Fernandez, G. A.; Vrabec, J.; Hasse, H. Int. J. Thermophys. 2004, 25, 175−186. (47) Erpenbeck, J. J. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38, 6255−6266. (48) Zhang, Y.; Otani, A.; Maginn, E. J. J. Chem. Theory Comput. 2015, 11, 3537−3546. (49) Schoen, M.; Hoheisel, C.; Beyer, O. Mol. Phys. 1986, 58, 699− 709. (50) Moore, S. G.; Wheeler, D. R. J. Chem. Phys. 2011, 134, 114514. (51) Moore, S. G.; Wheeler, D. R. J. Chem. Phys. 2012, 136, 164503.

I


Fourier Correlation Method for Simulating Mutual Diffusion

Recommend Documents