Application of Reverse Nonequilibrium Molecular Dynamics to the

Article Cite This: J. Phys. Chem. B XXXX, XXX, XXX−XXX

pubs.acs.org/JPCB

Application of Reverse Nonequilibrium Molecular Dynamics to the Calculation of the Mutual Diffusion Coefficient of Alkane Mixtures Hari Krishna Chilukoti,*,† Florian Müller-Plathe,† and Hua Yang‡ †

Downloaded via UNIV OF LOUISIANA AT LAFAYETTE on September 25, 2018 at 01:42:57 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

Technische Universität Darmstadt, Eduard-Zintl-Institut für Anorganische und Physikalische Chemie, Alarich-Weiss-Strasse 8, 64287 Darmsadt, Germany ‡ Key Laboratory of Inorganic−Organic Hybrid Functional Material Chemistry, Ministry of Education, College of Chemistry, Tianjin Normal University, Tianjin, 300387, China S Supporting Information *

ABSTRACT: In a recent publication, a reverse nonequilibrium molecular dynamics (RNEMD) method was presented for computing the mutual diffusion coefficient of liquid mixtures. A concentration gradient and a subsequent mass flux are induced in the system by suitably exchanging molecules in different regions. The algorithm has been successfully tested on Lennard-Jones mixtures and molecular fluid mixtures with molecules having the same number of particles. In this work, a modification is made to the RNEMD method to determine the mutual diffusion coefficient of binary liquid mixtures with molecules having different sizes and masses. To migrate molecules of a different type, the splitting method has been used in this work. Investigation of the resulting steady-state mass fraction profile allows the evaluation of the mutual diffusion coefficient. For validation, the mutual diffusion coefficients of ethane−propane and ethane−pentane liquid mixtures at different compositions and temperatures have been obtained using this method. The mutual diffusion coefficients obtained from the RNEMD simulations are within the error bars of values obtained by equilibrium molecular dynamics for the identical model and conditions. The excess energy released due to the exchange of molecules is efficiently removed by strongly coupling a local thermostat in the region around the insertion point. There is no heating of the analysis region.

1. INTRODUCTION Studies on mass transport in liquids by mutual diffusion are of fundamental importance for many applications in biology, chemistry, and chemical engineering.1−3 Accurate measurement of mutual diffusion coefficients over wide ranges of pressure and temperature is necessary because of its importance in the design of extractors, reactors, separators, and so on.4 Fick’s law and Maxwell−Stefan (MS) theory are commonly used to describe mutual diffusion in mixtures.5 Fick diffusivities (D12) are readily accessed from experiments since they are defined in terms of measurable concentration gradients. In contrast, MS diffusivities (Đ12) are defined by the chemical potential gradient and more cumbersome to obtain from experiments, because chemical potentials cannot be measured directly.6 Mass diffusion coefficients can be obtained from experiments, by analytical models, or by performing molecular dynamics (MD) simulations. Simulation has emerged as the preferred choice to determine mass transport properties when experiments are difficult or expensive to perform. Maxwell− Stefan diffusivities can be directly obtained from MD simulations. However, if we want to determine the Fick diffusion coefficient from MD simulations, the thermodynamic © XXXX American Chemical Society

factor, which links concentration and chemical potential, needs to be evaluated. Then, the thermodynamic factor is multiplied by the MS diffusion coefficient to evaluate Fick diffusion coefficient.5 The methods for the calculation of mutual diffusion coefficient from MD simulations have been reviewed by Liu et al.7 Numerous equilibrium MD (EMD) simulations have been performed to evaluate mutual diffusion coefficients of binary and ternary systems; for examples, see references.8−16 A few researchers have also performed nonequilibrium molecular dynamics (NEMD) simulations to determine the mutual diffusion coefficient of binary systems.4,17−23 In a recent work, we proposed a reverse nonequilibrium molecular dynamics (RNEMD) method to determine the mutual diffusion coefficient of binary liquid mixtures. This method has been tested on Lennard-Jones mixtures (Ar/Kr) and molecular liquid mixtures with molecules having the equal number of (united) atoms (water/deuterium oxide (H2O/ D 2 O) and united-atom benzene−cyclohexane (C 6 H 6 / C6H12)).22,23 This method follows the algorithm used by the Received: July 18, 2018 Revised: September 13, 2018 Published: September 13, 2018 A

DOI: 10.1021/acs.jpcb.8b06886 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B Table 1. Details of the Systemsa liquid mixture

molar fraction (xC2)

tempe (K)

no. of molecules (NC2)

no. of molecules (NC3/NC5)

C2H6−C3H8

0.0 0.32 0.50

0 450 700

1400 950 700

0.68 1.0

250 250 200 225 250 275 250 250

950 1400

450 0

0.0

250

0

900

0.32 0.50

250 225 250 275 300 250

450 700

950 700

950

450

C2H6−C5H12

0.68

density (kg/m3) at pressure of 30 MPa 563.2 ± 534.6 ± 579.8 ± 549.5 ± 516.2 ± 477.6 ± 496.0 ± 450.2 ± 500.3# 672.3 ± 688.5# 631.6 ± 627.6 ± 600.9 ± 572.4 ± 541.2 ± 561.4 ±

0.8 591.3# 0.8 0.7 0.9 1.1 1.0 1.0 1.1 1.2 0.9 0.9 1.0 0.9 1.0 1.0

The values with # are from experiments.31

a

original RNEMD for thermal conductivity24 and shear viscosity.25 A constant unphysical mass transport is imposed on the system by artificially exchanging positions and velocities of molecules periodically. Molecules are thus moved against their concentration gradient. As no molecules are lost, they will diffuse back by physical diffusion, trying to equalize their concentrations. The imposed unphysical transport and the diffusive flux are, thus, equal, once the stationary state is reached. This results in a stationary concentration profile between the exchange regions, and the concentration gradient is measured to obtain the mutual diffusion coefficient. The main challenge in this method is the exchange of molecules between different regions without generating excessive energy, which would increase the system temperature. For liquid mixtures of similar molecules (Ar/Kr, H2O/D2O, and C6H6/ C6H12) a thermostat with typical settings for liquids is sufficient to remove the excess energy created in the system due to the exchange of molecules. However, moving molecules of dissimilar sizes or masses or interaction strengths between different regions is challenging. This is the main obstacle for this method of calculating mutual diffusion of a liquid mixture with molecules of different sizes or masses. In the present paper, the RNEMD method is modified to extend its application range to the mutual diffusion coefficients of binary liquid mixtures with molecules of different sizes and masses. In order to exchange molecules of unlike size and mass between different regions, the splitting method is used, by which an existing molecule is split in a number of time steps into two daughter molecules. The excess energy generated due to the insertion of molecules in the splitting method is removed by strongly coupling a local thermostat to the insertion and neighboring regions. Thus, the temperature in the system is maintained close to the target value. The new version of the RNEMD method is validated by evaluating mutual diffusion coefficients of ethane−propane and ethane− pentane liquid mixtures.

test the applicability of the modified RNEMD method, two types of liquid alkane mixtures, ethane (C2H6)−propane (C3H8) and ethane−pentane (C5H12), are simulated in this work. The propane molecule is slightly dissimilar in mass (1.5 times heavier) and length (one C−C bond length longer) as compared to the ethane molecule. However, the pentane molecule is about 2.4 times heavier and two C−C bond lengths longer than ethane molecule. The ethane−propane and ethane−pentane liquid mixtures are simulated at three molar fractions of ethane, xC2 = 0.32, 0.5, and 0.68. Additionally, both liquids are also simulated at four different temperatures. The computational details of the examined systems are given in Table 1. This study is a test of the RNEMD method. Thus, comparison to the other method using the same model is done. The match to experiment is less important. Therefore, a united atom model is used to describe alkane molecules because of its simplicity and computational advantage. In the united atom representation of an alkane molecule, CH3 and CH2 groups are treated as single interaction sites. The molecular interactions between united atoms are treated using the TraPPE force field.26 In this force field, intramolecular interactions consist of bond bending and torsional motion. The harmonic potential is used to model the bending among any three united atoms in a molecule. V (θ ) =

1 kθ(θ − θ0)2 2

(1) −1

−2

where kθ is 519.65 kJ mol rad and θ0 is 114°. The bond angle is θ. The torsion between two united atoms in a molecule is considered as follows V (ϕ) = c0 + c1(1 + cos ϕ) + c 2(1 − cos 2ϕ) + c3(1 + cos 3ϕ)

(2)

where ϕ is dihedral angle, and c0 = 0 kJ/mol, c1 = 5.90 kJ/mol, c2 = −1.13 kJ/mol, c3 = 13.16 kJ/mol. The bond distances (CH3−CH2 and CH2−CH2) are maintained at a constant value of 1.54 Å by the SHAKE method.27 Intermolecular interaction between united atoms that are located on different molecules and united atoms situated on a molecule that are

2. COMPUTATIONAL DETAILS The objective of this paper is to extend the RNEMD method to evaluate the mutual diffusion coefficient of liquid mixtures with molecules having dissimilar masses and sizes. In order to B


Article

The Journal of Physical Chemistry B

temperature during the RNEMD simulations, respectively. The system reached a steady state after ≥24 ns in RNEMD simulations. For measuring the mutual diffusion coefficient, the data were collected for 60 ns after the steady state had been reached. The RNEMD procedure for the calculation of D12 of an ethane−pentane mixture is shown in Figure 1. The simulation system is divided into 20 slabs of equal thickness in the zdirection. The slab at the center of the box is defined as slab B, where an excess of ethane is desired. The slabs A1 and A2, appear at the end of the system, are the two symmetric parts of slab A where an excess of pentane or propane is desired. The mass flux is introduced into the system along the z-direction by moving ethane molecules from slab A to slab B and pentane molecules from slab B to slab A at periods of several thousand time steps (every 30−100 ps). We have demonstrated that a molecule can be directly exchanged without creating much excess energy when the liquid mixtures have molecules with an equal number of atoms and are similar in volume.22,23 In the ethane−pentane system, pentane molecule is considerably bigger in length and mass than ethane molecule. It is not possible to exchange ethane and pentane molecules directly between slabs A and B without creating overlaps and, hence, colossal forces. Therefore, we have used the splitting method,33 which was originally developed to examine steady-state evaporation of nanodroplets, to move ethane and pentane molecules between the exchange regions. The splitting method is briefly discussed below. For details, readers are referred to ref 33. Artificially moving a molecule means it must be destroyed in the from-region and recreated in the to-region. Deletion in the from-region is done without further ado and it is left to the surrounding molecules to fill the void that has been created. The recreation of the molecule in the to-region is performed by the splitting method, which is much more benign than a straight random insertion. First, a molecule of the same species, which resides in the to-region, is selected. Then an identical copy of this molecule with identical atomic coordinates is created. Thus, the parent molecule has been split into two molecules. In the implementation, coordinates of the deleted molecule are overwritten with the coordinates of the parent molecule in the insertion region (two molecules overlap completely). To avert the infinite large force between the overlapping molecules, the LJ interactions between them are switched off. However, the daughter molecules interact with other molecules in a normal manner. The atoms in the overlapping molecules interact with each other via Gaussian potential to push them apart gradually, which is softcore at a short distance

three united atoms apart is considered using the Lennard-Jones (LJ) potential. The interaction parameters between CH3 and CH2 groups are obtained using the Lorentz−Berthelot combining rules. σii + σjj εij = εiiεjj and σij = (3) 2 The LJ interactions beyond a cutoff radius (rc) of 14 Å are neglected and long-range corrections for LJ interactions are not applied in this work. All simulations in this work were performed using a modified version of the YASP package.28 The leapfrog algorithm29 was used to integrate equations of motion for united atoms and an integration time step 2 fs is used for all simulations. Rectangular parallelepiped simulation systems have been constructed as shown in Figure 1. The lengths of the systems

Figure 1. Schematic of the RNEMD method for calculating the mutual diffusion of an ethane−pentane liquid mixture. The green dashed line boxes indicate the observation volumes in which the mass flux (J1z) is measured relative to the mixture velocity (vz). Division of slabs in the computational system is also illustrated in the figure.

studied in this work are Lx = Ly ≈ 30 Å and Lz is in the range 175−225 Å depending on the type of liquid mixture, composition, and temperature. Periodic boundary conditions are applied in all three directions. The pressure and temperature of the system were maintained at a target value by connecting it to the Berendsen thermostat and barostat.30 Temperature coupling times of 0.2 and 0.002 ps are used in the analysis regions and in the insertion regions, respectively. A pressure coupling time of 5 ps is used during equilibration at the pressure of 30 MPa. Initial configurations were constructed by randomly positioning two types of alkane molecules in the simulation box at a desired composition using the Packmol package.32 After arranging the molecules, the system was equilibrated in NPT and NVT ensembles. After the system has reached the equilibrium state, the liquid density was calculated and the diffusion coefficients were calculated using Einstein relations for reference. Data were collected over a time of 50 ns for ethane−pentane system and 30 ns for ethane−propane system. The equilibrated states were used as starting configurations for subsequent RNEMD simulations. The RNEMD simulations were performed in NVT ensemble. It should be noted that the atoms in the insertion regions are only coupled (coupling time 0.002 ps) to a heat bath that is at the target temperature and atoms in the analysis regions are only coupled (coupling time 0.2 ps) to a heat bath at the target

2

Vg = −fk e−(x − xc)

/2w 2

(4)

where f k, xc, and w are amplitude, the center, and the width of the trough. The values of f k, xc, and w used in the present simulations are 75 kJ/mol, 0.5 nm, and 0.3 nm, respectively. The magnitude of f k, influences the speed at which the overlapped molecules are pushed away from each other. The interaction between the daughter molecules is switched back to the full LJ potential when the longest distance between any two atoms is over 0.7, 1.2, and 1.2 nm for ethane, propane, and pentane molecules, respectively (typically, it takes about 10 ps to become normal molecules). At this point, two daughter molecules have become normal molecules. The main technical aspect of this method is to reduce and to limit the insertion C


Article


similar at different exchange periods (D12 = 6.09 × 10−9, 5.94 × 10−9, 6.12 × 10−9 m2/s for texch = 90, 100, 120 ps, respectively). This also suggests that the calculated D12 is independent of exchange period in the present system. However, it should be noted that extremely low exchange periods might result in some error in the obtained mutual diffusion coefficients due to the influence of a strong concentration gradient. An exchange period of 100 ps is used for ethane−pentane mixtures, and 30 ps is used for ethane− propane mixtures in this study.

effects on local properties. The local excess energy generated due to the creation of molecules is removed by a strongly coupled thermostat which acts in the insertion slab and the two neighboring slabs (coupling time 0.002 ps). The other advantage of this method is that it can handle multiple insertions. Note that a different set of parameters might be needed in the Gaussian potential to push the daughter molecules away from each other when the lengths and masses of the molecules are different. A good criterion to select the parameters in the Gaussian potential is that the energy released during the insertion of molecules is small enough to be removed by the local thermostat and keep the temperature of the system close to the target temperature. This method is used to move both ethane and pentane molecules between slabs A and B, or ethane and propane molecules. It is to note if we exchange an equal number of molecules that are different in size and mass between slabs A and B, in addition to the diffusive mass flux due to the concentration gradient, there will be bulk fluid motion due to the pressure difference. Therefore, we exchanged on average 2.4 ethane molecules for one pentane molecule in the ethane−pentane system and on an average 1.5 ethane molecules for one propane molecule in the ethane− propane system to minimize the bulk fluid motion due to the pressure difference. In a binary mixture, the one-dimensional form of Fick’s first law of diffusion is J1z = −ρD12

dω1 dz

3. RESULTS AND DISCUSSION As discussed previously, maintaining the total energy in the presence of molecular exchange is a key technical necessity of the RNEMD method. If the total energy of the system increases, the system temperature will rise. The surge in energy mainly originates from inserting the larger of the two types of exchanged molecules. This could be particularly severe when positions of molecules with dissimilar masses and sizes are exchanged directly to move them to different regions. For example, in an ethane−pentane mixture, it is relatively easy to position an ethane molecule in place of a pentane molecule without creating a lot of excess energy. Conversely, it is extremely difficult to put a pentane molecule in place of an ethane molecule without causing a major overlap and a surge of energy. One part of the solution, to insert both molecules in already a dense liquid region, is the splitting method. The perturbations are local and short-lived.33 The other ingredient is a localized removal of the excess energy. The atoms in the splitting (inserting) slabs and the two slabs next to it are strongly coupled to a temperature bath (coupling time 0.002 ps) by a Berendsen thermostat30 to remove the excess energy and to limit the temperature rise in the system. The thermostat removes the excess energy quickly and prevents it from propagating into the other parts of the simulation box where the analysis is performed. Any artifacts would be limited to the insertion regions. It should be noted that the usage of the splitting method in the present calculations also confirms that it can be used to insert complicated molecules in a dense liquid region without altering properties far away from the insertion regions. To know how efficiently the excess energy is removed by the local thermostat when the molecules are exchanged with the splitting method, the average temperature in each slab for an equimolar ethane−pentane system at 250 K (global temperature with coupling time of 0.2 ps) and 30 MPa (global pressure with coupling time of 5 ps) is obtained over the production run (60 ns) when the molecules are exchanged for a period of 100 ps and shown in Figure 2. The average temperature of the system is 250.29 K, which is very close to the target temperature. The maximum variation in temperature in the system is about 2 K, which is considered acceptable in the calculation of D12. This indicates that the strongly coupled thermostat is able to efficiently remove the excess energy released during the insertion of molecules. The cumulative average of temperature against time in slabs A, 3, 8, and B are shown in Figure 3 for equimolar ethane−pentane mixture at 300 K and 30 MPa with particle exchange period 100 ps. The figure suggests that the local temperature in every slab fluctuates around the mean, which is close to the target temperature and there is no rise in the temperature. This indicates that the split molecules are gradually moved away from each without generating colossal forces and the excess energy generated during the splitting of molecules is efficiently

(5)

where J1z is the molecular mass flux of species “1” along zdirection, which should be defined relative to the mass average velocity of a mixture.34 ρ is the density and ω1 is the mass fraction of species 1. The mass flux J1z is defined as J1z = ρω1(v1z − vz)

(6)

where v1z is the instantaneous arithmetic average of the molecular velocity of species 1 in the observation volume. The mass average velocity vz in the z-direction is given as vz = ω1v1z + ω2v2z (7) By substituting eq 6 into eq 5, D12 can be evaluated as D12 = −

ω1(v1z − vz) dω1 dz

(8)

The mass flux introduced in the system relative to the mass average velocity was calculated using eq 7 in right and left observation volumes which are shown in Figure 1. The average value is used in the calculation of mutual diffusion coefficient. Due to the imposed molecule transfer, a concentration profile is established in the system. After a long enough time, the system reaches a steady state. The average mass fraction profile between the exchange regions is fitted to a straight line to evaluate the mass fraction gradient

dω1 dz

( ). Then the mutual

diffusion coefficient is calculated using eq 8. The magnitude of mass flux and the resulting mass fraction gradient are dependent on the time between molecule exchanges which is referred to as exchange period texch. It has been demonstrated that the mutual diffusion coefficient calculated using the RNEMD does not strongly depend on the exchange period.22 We have also checked that the D12 values calculated for ethane−pentane mixture at xC2 = 0.32, 250 K, and 30 MPa are D


Article


concentration of ethane is evaluated at different simulation times. The cumulative average of the concentration of ethane against time in slabs A, 3, 8, and B over the production run are shown in Figure 4 for equimolar ethane−pentane mixture at

Figure 2. Temperature profile of slabs averaged over the production run of 60 ns for an ethane−pentane system at xC2 = 0.5, 250 K, and 30 MPa with molecule exchange period texch = 100 ps. The production run of 60 ns is divided into five subsets and from each subset the temperature of the slabs is evaluated, and these values are used to obtain the standard deviation, which is shown as error bars.

Figure 4. Cumulative average of concentration of ethane molecules versus time in slabs A, 3, 8, B for an ethane−pentane system at xC2 = 0.5, 250 K, and 30 MPa with molecule exchange period texch = 100 ps.

300 K and 30 MPa with particle exchange period 100 ps. It is noticed that concentration exhibit no systematic trend over the production time 60 ns. The figure also suggests that the concentration of ethane in slabs only fluctuate around a mean. These tendencies indicate that our systems have reached a steady state during last 60 ns of the simulation time. The gradient of the concentration profile is related to the mutual diffusion. The mass fraction of ethane along z-direction is illustrated in Figure 5 for the ethane−pentane system at

Figure 3. Cumulative average of temperature versus time in slabs A, 3, 8, B for an ethane−pentane system at xC2 = 0.5, 250 K, and 30 MPa with molecule exchange period texch = 100 ps.

Figure 5. Mass fraction of ethane as a function of distance in an ethane−pentane mixture at 250 K and 30 MPa with molecule exchange period texch = 100 ps. The production run of 60 ns is divided into five subsets and from each subset the mass fraction of ethane is evaluated, and these values are used to obtain the standard deviation, which is shown as error bars.

removed by the local thermostat. It is noticed that temperature exhibits no systematic trend over the production time 60 ns. The temperature fluctuations become smaller when the dissimilarity between the exchanged molecules is relatively small. The systems attain a steady state after molecules are exchanged for a time of over 24 ns in RNEMD simulations. The data collected after this time is used for the calculation of the mutual diffusion coefficient. In order to demonstrate that our systems have reached a steady state, the mass

different compositions, when the molecules are exchanged at a period of 100 ps. Different mass fraction gradients of ethane are observed for different compositions in the ethane−pentane mixture at the same exchange period. The average mass fraction of ethane in slabs related by symmetry (slabs 3 and 18, 4 and 17, etc.) is within the error bars. This is also an E


Article

The Journal of Physical Chemistry B indication that the systems have reached a steady state in spite of their natural fluctuations. The mass fraction of pentane changes in opposite direction to that of mass fraction of ethane in all systems (not shown). The mass fraction profiles are linear in between exchange regions. The mass fractions from slab 3 to slab 9 and slab 13 to slab 19 are fitted to a line to get the slope values ( dω1 ). The average of two slope values is used dz in the calculation of the mutual diffusion coefficient by eq 8. The mass flux distribution of ethane in an ethane−propane system also exhibits tendencies similar to those of the ethane− pentane system (not shown) at different compositions. As described previously, the mass flux should be measured with respect to the mass average velocity. Figure 6 shows the

Table 2. Mutual Diffusion Coefficients for Ethane−Propane and Ethane−Pentane Liquid Mixtures at Different Temperatures and Concentrations at 30 MPaa liquid mixture C2H6−C3H8

C2H6−C5H12

molar fraction (xC2) 0.0 0.32 0.50

0.68 1.0 0.0 0.32 0.5

0.68 a

Figure 6. Cumulative average of

J1z ρ

250 250 200 225 250 275 250 250 250 250 225 250 275 300 250

mutual diffusion coeff (m2/s × 10−9)

self-diffusion coeff (m2/s × 10−9) 10.22 ± 0.16

9.40 5.45 7.61 10.31 11.64 15.09 11.07

± ± ± ± ± ± ±

0.45 0.33 0.69 0.59 0.26* 1.04 0.45 18.60 ± 0.27 5.17 ± 0.14

5.91 6.05 7.59 7.82 10.48 12.48 6.38

± ± ± ± ± ± ±

0.46 0.48 1.26 0.13* 1.66 2.76 1.41

The values with * are obtained by EMD simulations.

Relatively large error bars are observed for ethane−pentane mixtures, which is due to a lower signal-to-noise ratio, a known phenomenon in RNEMD calculations. It is to note that the uncertainty of the Fick diffusion coefficient calculated for each subset is primarily due to the uncertainty of the concentration gradient. The relation between diffusion coefficient and temperature can be described by Arrhenius equation

versus time in an ethane−pentane

mixture at xC2 = 0.5, 250 K, and 30 MPa with molecule exchange period texch = 100 ps.

D = D0 exp( −E /RT )

J1z

cumulative average of the quantity mass flux by density ( )

(9)

where D0 is a constant, E is the activation energy, R is the universal gas constant, and T is the temperature of the examined liquid. By fitting a linear line to ln(D) versus inverse of temperature (Figure S1 in the Supporting Information), we can get the activation energy. The activation energies obtained from our simulations for ethane−propane and ethane−pentane mixtures are 6.15 ± 0.16 and 5.59 ± 0.56 kJ/mol, respectively. All the obtained D12 values are in the typical range of liquids’ diffusion coefficients 10−9.2 With an increase in the molar fraction of ethane (xC2), the D12 of ethane−propane liquids increases in the examined range of xC2. A similar trend is observed by Krishna and van Baten35 using MD simulations for methane−propane at 30 MPa and 333 K. However, the mutual diffusion coefficient of the ethane−pentane mixture for xC2 = 0.68 is lower than for xC2 = 0.33 and 0.5. This kind of tendency was observed for a methane−hexane mixture at 333 K and 30 MPa.35 Krishna and van Baten35 also reported that the Fick diffusion coefficient reaches a minimum at around xC1 = 0.9, which is due to the proximity to the spinodal composition. To benchmark our mutual diffusion coefficient values obtained from the modified RNEMD, diffusion coefficients of equimolar ethane−pentane liquid mixture at 250 K and 30 MPa have also been determined by equilibrium molecular dynamics simulations. The self-diffusion coefficients of ethane and pentane were calculated by evaluating the slope of the mean square displacement (MSD) of ethane or pentane molecules against time line in its linear region by the Einstein

ρ

against time over the last 60 ns in the left and right observation J volumes, which are depicted in Figure 1. The quantity 1z in an ρ

observation volume is evaluated using eq 6. The figure suggests that the mass flux is constant during the production time in both observation volumes as it should be in a steady state. The mass flux in the left observation volume shows a negative value J due to directionality and the absolute magnitude average of 1z ρ

is used in the calculation of the mutual diffusion coefficient. Using eq 8, the mutual diffusion coefficient D12 is determined J with the averaged values of 1z and dω1 . ρ

temp (K)

dz

To examine the influence of temperature on mutual diffusion, ethane−propane and ethane−pentane mixtures are simulated at four different temperatures for xC2 = 0.5 and p = 30 MPa. Additionally, the effect of concentration on mutual diffusion is examined by simulating both mixtures at three compositions, xC2 = 0.32, 0.5, 0.68. The obtained mutual diffusion coefficient values of ethane−propane and ethane− pentane mixtures for different molar fractions and temperatures are given in Table 2, which were calculated with exchange periods of 30 and 100 ps, respectively. To obtain the error bars, the production run of 60 ns was divided into five equal subsets. From each subset, a mutual diffusion coefficient is obtained and these values are used to obtain the standard deviation, which is shown as the error value in the table. F


Article


conditions. This indicates the new RNEMD method is a feasible alternative to calculate mutual diffusion coefficient of liquid mixtures. Moreover, the present results also demonstrate that the splitting method can be used to insert molecules at high frequency into dense liquid phase without creating artifacts in computational systems.

relation (eq 10). The mutual diffusion coefficient was obtained from the slope of the MSD of center of mass of all ethane or pentane molecules versus time in its linear region (Figure S2 in the Supporting Information) with an additional prefactor.16 Di =

1 d lim ⟨|ri(t ) − ri(0)i|2 ⟩ 6 t →∞ dt

■

(i: molecule index) (10)

ij 1 1 yzz N d Đ12 = x1x 2jjj + zzz xiM i)2 jx M x M 6 dt 2 2{ k 1 1

S Supporting Information *

2

× ⟨|rAcm(t ) − rAcm(0)|2 ⟩

(A = C2H6orC5H10)

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.8b06886. Fick diffusion coefficient versus inverse of temperature (Figure S1), MSD versus time plots for MS diffusivity calculation (Figure S2), running Kirkwood−Buff integrals (Figure S3), and thermodynamic factor calculation procedure (PDF)

(11)

where ri(t) is the position vector of center of mass of molecule i (C2H6 or C5H12) at time t. rcm A (t) =

N

∑i ∈AA ri(t )·M i ∑i ∈ A M i

ASSOCIATED CONTENT

is the position

■

vector of center of mass of all particles of species A (C2H6 or C5H12) at time t. N is the total number of molecules in the system. xA and MA are molar fraction and molecular weight of species A (C2H6 or C5H12). It should be noted that the value obtained from eq 11 needs to be multiplied by a thermodynamic factor Γ in order to consider the nonideality of the liquid mixture.36,37 We have adopted the Kirkwood− Buff coefficient method6,7,38 (details are given in the Supporting Information) to calculate the thermodynamic factor for ethane−pentane mixture, which is equal to 1.02 ± 0.01. The diffusion coefficients DC2H6, DC5H12, and D12 obtained from EMD using eqs 10 and 11 for ethane−pentane mixture are (10.00 ± 0.04) × 10−9, (7.47 ± 0.35) × 10−9, and (7.82 ± 0.12) × 10−9 m2/s, respectively. The mutual diffusion coefficient D12 obtained from the RNEMD is (7.59 ± 1.26) × 10−9 m2/s, which agrees well with the value determined from EMD. The present results suggest that the RNEMD method is a viable alternative method for calculating reliable mutual diffusion coefficient for liquid mixtures with molecules of dissimilar masses or sizes.

AUTHOR INFORMATION

Corresponding Author

*H. K. Chilukoti. E-mail: [email protected]. ORCID

Hari Krishna Chilukoti: 0000-0002-7882-0805 Florian Müller-Plathe: 0000-0002-9111-7786 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by the Collaborative Research Centre Transregio 75 Droplet Dynamics under Extreme Ambient Conditions of the Deutsche Forschungsgemeinschaft (DFG). H.K. acknowledges correspondence with Tatjana Janzen, Jianguo Zhang, Prof. Higashi Hidenori, and Prof. Rajamani Krishna.

■

REFERENCES

(1) Krishna, R.; Wesselingh, J. A. The Maxwell-Stefan approach to mass transfer. Chem. Eng. Sci. 1997, 52, 861−911. (2) Cussler, E. L. Diffusion mass transfer in fluid systems; Cambridge University Press: Cambridge, U.K., 2009. (3) Hendriks, E.; Kontogeorgis, G. M.; Dohrn, R.; de Hemptinne, J.C.; Economou, I. G.; Ž ilnik, L. F.; Vesovic, V. Industrial requirements for thermodynamics and transport properties. Ind. Eng. Chem. Res. 2010, 49, 11131−11141. (4) Higashi, H.; Tamura, K.; Seto, T.; Otani, Y. Direct calculation of mutual diffusion coefficients of binary system using non-equilibrium molecular dynamics simulation. Fluid Phase Equilib. 2015, 402, 83− 88. (5) Taylor, R.; Krishna, R. Multicomponent mass transfer; Wiley, 1993. (6) Liu, X.; Schnell, S. K.; Simon, J.-M.; Bedeaux, D.; Kjelstrup, S.; Bardow, A.; Vlugt, T. J. H. Fick diffusion coefficients of liquid mixtures directly obtained from equilibrium molecular dynamics. J. Phys. Chem. B 2011, 115, 12921−12929. (7) Liu, X.; Schnell, S. K.; Simon, J.-M.; Krüger, P.; Bedeaux, D.; Kjelstrup, S.; Bardow, A.; Vlugt, T. J. H. Diffusion coefficients from molecular dynamics simulations in binary and ternary mixtures. Int. J. Thermophys. 2013, 34, 1169−1196. (8) 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−6947. (9) Guevara-Carrion, G.; Janzen, T.; Muñoz-Muñoz, Y. M.; Vrabec, J. Mutual diffusion of binary liquid mixtures containing methanol,

4. CONCLUSIONS The RNEMD method has been extended to the calculation of the mutual diffusion coefficient of liquid mixtures from molecules with the same number of particles to molecules having dissimilar masses and sizes. It has been tested on ethane−propane and ethane−pentane liquid mixtures. Pentane and propane molecules are about 2.4 and 1.5 times more voluminous than ethane molecule, respectively. A mass flux is introduced in the computational systems by periodically exchanging molecules weighted with their molecular weight between different regions. Since the molecules in the tested liquids are considerably different in size and mass, a direct exchange of molecules is not possible without generating colossal excess energy, which in turn leads to an increase in temperature of the system. Therefore, the splitting method has been used to transfer molecules between different regions. The excess energy generated due to the insertion of molecules in the splitting method is efficiently removed by a strongly coupled local thermostat acting in the exchange and neighboring regions. It prevents heat from spreading into other parts of the system. Thus, the temperature in the system is maintained at the target value. The obtained mutual diffusion coefficients of alkane mixtures are very close to those determined using Einstein relation in equilibrium molecular dynamics simulations for the same model and G


Article

The Journal of Physical Chemistry B ethanol, acetone, benzene, cyclohexane, toluene, and carbon tetrachloride. J. Chem. Phys. 2016, 144, 124501. (10) Harmandaris, V. A.; Angelopoulou, D.; Mavrantzas, V. G.; Theodorou, D. N. Atomistic molecular dynamics simulation of diffusion in binary liquid n -alkane mixtures. J. Chem. Phys. 2002, 116, 7656−7665. (11) Keffer, D. J.; Adhangale, P. The composition dependence of self and transport diffusivities from molecular dynamics simulations. Chem. Eng. J. 2004, 100, 51−69. (12) Rowley, R. L.; Stoker, J. M.; Giles, N. F. Molecular-dynamics simulation of mutual diffusion in nonideal liquid mixtures. Int. J. Thermophys. 1991, 12, 501−513. (13) Schoen, M.; Hoheisel, C. The mutual diffusion coefficient D 12 in liquid model mixtures A molecular dynamics study based on Lennard-Jones (12−6) potentials. Mol. Phys. 1984, 52, 1029−1042. (14) Par̆ez, 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. (15) Wheeler, D. R.; Newman, J. Molecular dynamics simulations of multicomponent diffusion. 2. Nonequilibrium method. J. Phys. Chem. B 2004, 108, 18362−18367. (16) Schmitz, H.; Faller, R.; Müller-Plathe, F. Molecular mobility in cyclic hydrocarbons: A simulation study. J. Phys. Chem. B 1999, 103, 9731−9737. (17) Thompson, A. P.; Ford, D. M.; Heffelfinger, G. S. Direct molecular simulation of gradient-driven diffusion. J. Chem. Phys. 1998, 109, 6406. (18) Tsige, M.; Grest, G. S. Molecular dynamics simulation of solvent−polymer interdiffusion: Fickian diffusion. J. Chem. Phys. 2004, 120, 2989−2995. (19) Hafskjold, B.; Ratkje, S. K. Criteria for local equilibrium in a system with transport of heat and mass. J. Stat. Phys. 1995, 78, 463− 494. (20) Frentrup, H.; Avendaño, C.; Horsch, M.; Salih, A.; Müller, E. A. Transport diffusivities of fluids in nanopores by non-equilibrium molecular dynamics simulation. Mol. Simul. 2012, 38, 540−553. (21) Erpenbeck, J. J.; Kincaid, J. M. Calculation of the mutual diffusion coefficient by equilibrium and nonequilibrium molecular dynamics. Int. J. Thermophys. 1986, 7, 305−317. (22) Yang, H.; Zhang, J.; Müller-Plathe, F.; Yang, Y. A reverse nonequilibrium molecular dynamics method for calculating the mutual diffusion coefficient for binary fluids. Chem. Eng. Sci. 2015, 130, 1−7. (23) Yang, H.; Zhang, J.; Müller-Plathe, F. Extending reverse nonequilibrium molecular dynamics to the calculation of mutual diffusion coefficients in molecular fluid mixtures. Mol. Simul. 2016, 42, 1379−1384. (24) Müller-Plathe, F. A simple nonequilibrium molecular dynamics method for calculating the thermal conductivity. J. Chem. Phys. 1997, 106, 6082. (25) Müller-Plathe, F. Reversing the perturbation in nonequilibrium molecular dynamics: An easy way to calculate the shear viscosity of fluids. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1999, 59, 4894−4898. (26) Martin, M. G.; Siepmann, J. I. Transferable potentials for phase equilibria. 1. united-atom description of n-alkanes. J. Phys. Chem. B 1998, 102, 2569−2577. (27) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 1977, 23, 327−341. (28) Müller-Plathe, F. YASP: A molecular simulation package. Comput. Phys. Commun. 1993, 78, 77−94. (29) Allen, M. P.; Tildesley, D. J. Computer simulation of liquids; Clarendon Press, 1989. (30) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 1984, 81, 3684−3690.

(31) Afeefy, H. Y.; Liebman, J. F.; Stein, S. E. Neutral Thermochemical Data. In NIST Chemistry WebBook; Linstrom, P. J., Mallard, W. G., Eds; NIST Standard Reference Database Number 69; National Institute of Standards and Technology: Gaithersburg, MD, https://www.nist.gov/ (retrieved August 27, 2018). (32) Martínez, L.; Andrade, R.; Birgin, E. G.; Martínez, J. M. PACKMOL: A package for building initial configurations for molecular dynamics simulations. J. Comput. Chem. 2009, 30, 2157− 2164. (33) Zhang, J.; Müller-Plathe, F.; Yahia-Ouahmed, M.; Leroy, F. A steady-state non-equilibrium molecular dynamics approach for the study of evaporation processes. J. Chem. Phys. 2013, 139, 134701. (34) Bird, R. B.; Stewart, W. E.; Lightfoot, E. N. Transport phenomena; Wiley, 2007. (35) Krishna, R.; van Baten, J. M. Describing diffusion in fluid mixtures at elevated pressures by combining the Maxwell−Stefan formulation with an equation of state. Chem. Eng. Sci. 2016, 153, 174−187. (36) Jolly, D. L.; Bearman, R. J. Molecular dynamics simulation of the mutual and self diffusion coefficients in Lennard-Jones liquid mixtures. Mol. Phys. 1980, 41, 137−147. (37) Galliero, G.; Volz, S. Thermodiffusion in model nanofluids by molecular dynamics simulations. J. Chem. Phys. 2008, 128, 064505. (38) Schnell, S. K.; Liu, X.; Simon, J.-M.; Bardow, A.; Bedeaux, D.; Vlugt, T. J. H.; Kjelstrup, S. Calculating thermodynamic properties from fluctuations at small scales. J. Phys. Chem. B 2011, 115, 10911− 10918.

H


Application of Reverse Nonequilibrium Molecular Dynamics to the

Recommend Documents