INTRODUCTION
Transport properties are of great importance for theoretical and industrial purposes. In industry, accurate values of transport properties are essential for the correct design of equipment and processes to perform tasks such as the transport of complex fluids (e.g., polymer melts or solutions) and efficient heat exchange (^{Dysthe et al., 1998}). For theoretical and academic purposes, the development of new models and methodologies capable of predicting these properties can decrease the need of experiments, mainly in extreme (and thus expensive) experimental conditions. Along with the well-known energy, mass, and momentum transport, coupled transport phenomena such as Soret and Dufour effects may take place in a relevant way. The Soret effect, also known as thermal diffusion, is one of the most studied coupled transport phenomena (^{Demirel, 2007}) and is the one addressed here.
Thermal diffusion is observed when a mixture is subjected to a temperature gradient which induces not only heat flux, but also a material flux that causes segregation of chemical species in the mixture. This distribution is consistent with the chemical potential gradient for each species, which is the driving force for diffusive flux. In a closed system, a steady state is reached, resulting in a well-developed stationary concentration gradient (^{Leahy-Dios, 2008}; ^{De Groot and Mazur, 1984}).
Even around 160 years after the first observations of thermal diffusion made by Ludwig in 1856, a complete knowledge of the behavior of mixtures under temperature gradients is still lacking (^{Kincaid and Hafskjold, 1994}; ^{Leahy-Dios, 2008}; ^{Würger, 2013}). Nevertheless, thermal diffusion plays a key role in several natural and industrial processes. Currently, the most discussed application is in the petroleum industry. Because of geothermal gradients, diffusive fluxes can occur in the rock matrix of an oil reservoir, interfering in the distribution of components along its height. Therefore, thermal diffusion might be the reason for the unexpected compositional gradient observed in some oil reservoirs. A typical and famous example is the Yufutsu field in Japan, where heavy compounds are concentrated in the top of the reservoir, and vice-versa (^{Gorayeb et al., 2000}).
Over the years, several researchers have focused on measurement and prediction of the behavior of mixtures under a thermal gradient. ^{Chapman (1916)} and ^{Enskog (1911)} were among the first authors to theoretically describe segregation due to thermal diffusion in dilute gas mixtures. More recently, Molecular Dynamics (MD) has emerged as a useful tool to better understand the phenomenon on the microscopic scale. Several methodologies based on equilibrium and non-equilibrium MD were developed to predict thermodynamic and transport properties, including thermal diffusion coefficients (^{Artola and Rousseau, 2013}). Transport properties such as thermal conductivity, viscosity and diffusion coefficient can be computed by equilibrium MD using the well-known Green-Kubo relations (^{Dysthe et al., 1998}; ^{Fernández et al., 2004}; ^{Liang and Tsai, 2010}). For coupled transport properties, the most successful methods are the Synthetic Non-Equilibrium Molecular Dynamics (S-NEMD) and the Boundary-driven Non-Equilibrium Molecular Dynamics (BD-NEMD) methods. In the S-NEMD method, an external field, which must satisfy the adiabatic incompressibility of phase space, is applied to drive the system out of equilibrium (^{Morris and Evans, 1990}). The S-NEMD method, together with Green-Kubo relations, can be used to determine the phenomenological coefficients. The main drawback of this method consists of an indirect measurement of the heat flux by molecular simulation for the case of non-ideal mixtures, since only the internal energy flux can be expressed in terms of microscopic quantities (^{Perronace et al., 2002}; ^{Artola and Rousseau, 2013}). In the BD-NEMD method, a simulation box with periodic boundary conditions is divided along one direction (e.g., the z-axis) into slabs of equal thickness. Then, two slabs far apart from each other are acted upon by adding energy to one slab (tagged as "hot") while removing energy from the other slab (tagged as "cold"). In this way, energy will naturally flow across the system in the opposite direction, that is, from the hot slab to the cold slab, leading to a steady state condition with a well-developed temperature gradient. There are different methods described in the literature to carry out the described energy transfer. The most cited algorithms are the one developed by ^{Hafskjold et al. (1993)}, referred to as the heat exchange or HeX algorithm, and the one developed by ^{Müller-Plathe (1997)}, referred to as the momentum exchange or PeX algorithm. In the HeX method, internal energy is added to or removed from a slab by rescaling local velocities according to the amount of energy to be transferred (^{Hafksjold et al., 1993}). In the PeX method, energy transfer is achieved by a velocity swapping between the coldest (i.e., slowest) particle of the hot slab and the hottest (i.e., fastest) particle of the cold slab (^{Müller-Plathe, 1997}). Advantages of BD-NEMD over S-NEMD are the use of ordinary Newtonian dynamics in intermediate regions and the exact knowledge of the amount of energy exchanged between the tagged slabs, which avoids the definition of a microscopic heat flux (^{Müller-Plathe and Reith, 1999}; ^{Müller-Plathe, 1997}). Despite these advantages, several characteristics are required for the method to be used properly. Because the simulation box is divided along one direction into several (commonly 20) slabs, the number of atoms inside each slab must be large enough for it to behave as a thermodynamic system (^{Hafskjold, et al., 1993}; ^{Hafskjold and Ratkje, 1995}; ^{Zhang and Müller-Plathe, 2005}; ^{Polyakov et al., 2008}). Another problem is related to the relatively large temperature gradients needed to obtain a suitable estimate of the Soret coefficient. Because cold slabs can reach very low temperatures, the dynamics of the system may be affected and the low energy regions may compromise the movement of particles in one direction due to formation of regions behaving like a glassy system. Another problem frequently observed is an artificial drift in the total energy of the system during the simulation run, which can interfere in the resulting transport property value. According to ^{Zhang et al. (2005)}, this problem can be solved by weakly coupling a thermostat to the system. However, the addition of a thermostat to the system causes the intermediate slabs not to follow Newtonian mechanics, which can artificially change the transport properties. ^{Zhang et al. (2005)} compared the influence of a Berendsen thermostat (^{Berendsen et al., 1984}) on the computed thermal conductivity of pure benzene, cyclohexane, n-hexane, and of mixtures of benzene and cyclohexane, using the PeX algorithm of ^{Müller-Plathe (1997)}. ^{Zhang et al. (2005)} observed a difference of around 10% for thermal conductivities obtained in thermostated simulations compared to thermal conductivities obtained in the absence of a thermostat. These deviations were considered to be small because common coupling times used in BD-NEMD simulations are one order of magnitude larger than the one used in the work of ^{Zhang et al. (2005)}. However, the authors carried out simulations with total times of 8000 ps. As we are going to show here, this time scale could be too small to compute thermal diffusion coefficients properly. Therefore, the conclusions drawn for thermal conductivity calculations may be not directly extendable to thermal diffusion coefficients. According to ^{Zhang et al. (2005)}, the energy drift can be explained by a round-off or cutoff noise, which is more prominent for charged molecules. According to ^{Zhang et al. (2005)}, the use of a Berendsen thermostat with a small coupling time showed no significant effect on the thermal conductivity value. According to these authors, it is because in the Berendsen thermostat the atom velocities are uniformly rescaled and not changed with respect to local friction, maintaining the velocity directions, but changing only their magnitude, which causes other thermostats to be inappropriate for the task.
Other authors employed the BD-NEMD method to determine thermal diffusion factors and/or Soret coefficients by using a weakly coupled thermostat to avoid the energy drift (^{Reith and Müller-Plathe, 2000}; ^{Zhang and Müller-Plathe, 2005}; ^{Polyakov et al., 2008}"). In addition, other methods that use two independent thermostats coupled to different regions of the simulation box may be used to generate the temperature gradient, also avoiding the energy drift problem (^{Maier et al., 2011}; ^{Maier et al., 2012}). ^{Simon et al. (1998}; ^{1999)} obtained thermal diffusion factors for mixtures of methane and n-decane; ^{Perronace et al. (2002)} computed Soret coefficients for mixtures of n-pentane and n-decane at three concentrations and obtained good results; ^{Nieto-Draghi and Ávalos (2005)} determined Soret coefficients for various aqueous systems; ^{Zhang and Müller-Plathe (2005)} determined Soret coefficients of benzenecyclohexane mixtures. Recently, the method was used to obtain Soret coefficients for fluids confined in slit pores (^{Galliéro et al., 2006}, ^{Hanaoui et al., 2011}, and ^{Hannaoui et al., 2013}). The BD-NEMD method was also applied for mixtures containing more than two components (^{Galliéro et al., 2010}; Galliéro et al., 2013).
Here, we used the BD-NEMD with the PeX algorithm proposed by ^{Nieto-Draghi and Ávalos (2003)} to obtain the Soret coefficient of an n-pentane/ndecane mixture at 298 K and 1 atm. The main goal of this work was the evaluation of how the parameters and conditions used in the simulations affect the observed energy drift and, consequently, the computed Soret coefficients and their fluctuations. Parameters and conditions analyzed are the size (number of molecules), time step size, frequency of momentum exchange, and the use of a Berendsen thermostat. The n-pentane/n-decane mixture was chosen so that the results obtained here could be compared to experimental data and to other molecular simulation results (^{Perronace et al., 2002}). The TraPPE-UA force field (^{Martin and Siepmann, 1998}) was used to model n-alkanes because it provides reliable results for density and diffusion coefficients, as pointed out by ^{Makrodimitri et al. (2011)}.
THEORETICAL BASIS
In a binary mixture under a thermal gradient, the mass flux of component 1 may be described by Equation (1).
where J_{1} is the mass flux of component 1, T is the temperature, λ is the heat conduction coefficient, D is the diffusion coefficient, DT is the thermal diffusion coefficient and w is the mass fraction of a given component. Once a steady state is reached, the mass flux due to thermal diffusion cancels out with the diffusive flux generated by the concentration gradient, leading to a zero net mass flux (J_{1} = 0), which results in
where x_{1} and x_{2} are mole fractions and S_{T} is the so-called Soret coefficient. This equation can be directly applied to determine the Soret coefficient from a molecular dynamics simulation, provided that linear temperature and composition gradients are obtained at the steady state condition. The theory of non-equilibrium thermodynamics is fully described in the literature (^{De Groot and Mazur, 1984}; ^{Demirel, 2007}).
METHODOLOGY
Simulation Details
In this work, a united atom force field was employed to model n-alkanes by considering CH_{2} or CH_{3} groups as single pseudo atoms. The potential energy is taken as a sum of dihedral torsion angle, bond angle bending, bond stretching, and inter- and intramolecular non-bonded contributions, as follows:
^{Makrodimnitri et al. (2011)} showed that the TraPPE-UA force field (^{Martin and Siepmann, 1998}), augmented with a harmonic bond stretching potential, can be used to compute densities and diffusion coefficients that agree well with experimental data for small molecules in heavy hydrocarbon. The individual contributions to the augmented TraPPE-UA force field used here are expressed by:
where k_{θ} and k_{l} are the force constants; c_{1}, c_{2} and c_{3} are constants of the dihedral model, φ is the angle formed by dihedral projection; ε_{ij} and σ_{ij} are the energy and size parameter of the Lennard-Jones model, respectively, rc is the cutoff distance which was used according to the TraPPE-UA force field. An united atoms force field was chosen to simulate systems containing a larger number of molecules with a smaller computational time. In this way, fluctuations in composition gradients are reduced. The force field parameters used are given in Table 1. Values of ε_{ij} and σ_{ij}, for unlike interactions were determined by the Lorentz-Berthelot combining rule, defined as follows:
All simulations were carried out on a GPU-based cluster using the software LAMMPS (^{Plimpton, 1995}; ^{Brown et al., 2011}; ^{Brown et al., 2012}). The velocity-Verlet algorithm (^{Verlet, 1967}) was used to integrate the equations of motion. The simulation box was considered with periodic boundary conditions in all directions. The time step size, number of molecules, and the statistical ensemble (NVE, NVT, or NPT) used depended on the simulation run (each case is discussed in the next sections). Standard long-range corrections were taken into account in all simulations (^{Allen and Tildesley, 1989}).
Interaction Type | Group | Parameters | ||||
---|---|---|---|---|---|---|
Bond stretching | k_{r} /k_{b} (K/Å^{2}) | r_{eq} (Å) | ||||
C - C | 96500 | 1.54 | ||||
Bond angle bending | C - C - C |
k_{θ} / k_{b} (K/ rad^{2}) 62500 |
θ_{eq}(º) 114.0 |
|||
Dihedral torsional angle | c_{0}/ k_{b} (K) | c_{1} / k_{b} (K) | c_{2} / k_{b} (K) | c_{3}/ k_{b} (K) | ||
C - C - C - C | 0 | 355.03 | -68.19 | 791.32 | ||
Lennard-Jones | ε/ k_{b} (K) | σ(Å) | r_{c}(Å) | |||
CH_{2} | 46.0 | 3.950 | 14.0 | |||
CH_{3} | 98.0 | 3.750 | 14.0 |
Determination of Density
Equilibrium MD simulations using NPT and NVT ensembles with a time step of 1fs were carried out in order to determine the density of pure components and mixtures. For each pure component (n-pentane and n-decane), a system consisting of 800 molecules was used. For binary mixtures, equimolar mixtures of 1200 molecules were used. To pack the molecules into a cubic simulation box, the Packmol (^{Martínez, et al., 2009}) software was used to reach a density between 0.500 and 0.600 g/cm^{3}. To avoid trapping in a local minimum of energy, the systems were first warmed up from 298 to 700 K with a 1 fs time step during 25 ps employing a ramped Nosé-Hoover chain thermostat using the formulation of ^{Shinoda et al. (2004)} with a coupling time of 100 fs. After reaching 700 K, the system was equilibrated in the NVT ensemble for 2 ns. Next, it was gradually cooled down to a temperature of 298 K during 300 ps, being thermally equilibrated at 298 K for 2 ns in the NVT ensemble. Once equilibrated, the NPT ensemble with the Nosé-Hoover barostat (following the formulation of ^{Shinoda et al., 2004}) with a coupling time of 1000 fs was employed to equilibrate the volume using a target pressure of 1 atm. The system volumes were allowed to equilibrate for 2 ns and a sampling time of at least 10 ns was used to determine the average volume.
Determination of Self-Diffusion Coefficients
Equilibrium MD simulations were carried out with pure component systems to determine their self-diffusion coefficients (D). The results were used to test the accuracy of the force field in the determination of this transport property by comparison with experimental data (^{Douglas and McCall, 1958}). A time step of 1 fs was used in all simulation runs. After the determination of density, a new system containing 800 molecules of the desired component was packed into a cubic simulation box with periodic boundary conditions and a volume equal to the previously determined equilibrium volume. The procedure described in the preceding section, excluding the volume equilibration step, was used again to avoid local minima. Next, an equilibration run in the NVE ensemble was carried out for 5 ns. The final configuration obtained was used as initial configuration in the simulations for determining D, which is also performed in the NVE ensemble. Einstein's equation was used directly, that is,
where Δr^{2} is the mean square displacement of the center of mass of molecules averaged over all centers of mass of each kind of molecule and multiple time origins and d is the dimensionality of the system. The simulation time (t) was at least 10 ns for all self-diffusion coefficient determinations, which was large enough to achieve a constant slope behavior, corresponding to D in Equation (10).
The BD-NEMD Algorithm
Boundary-driven non-equilibrium molecular dynamics was employed to determine the Soret Coefficient of an equimolar mixture of n-decane and npentane at 298 K. The PeX algorithm of ^{Nieto-Draghi and Avalos (2003)} was used in all simulations. This algorithm is similar to the original PeX algorithm of ^{Müller-Plathe (1997)}, but including the possibility of energy exchanges between atoms with different masses, changing their momenta through a virtual collision.
In this work, we used a tetragonal simulation box with periodic boundary conditions and dimensions L × L × 3L. The box was divided in the z-axis direction (length equal to 3L) into 20 slabs of equal thickness. With a predefined time interval, the atom with the smallest velocity in the hot slab (slab 10) and the atom with the highest velocity in the cold slab (slab 0) undergo a virtual elastic collision. In this way, kinetic energy is transferred from the fastest atom to the slowest atom, which most frequently yields a transfer from the cold slab to the hot slab as the system evolves towards a steady state. Depending on the frequency of virtual collisions, the system reaches steady state with a linear temperature profiles in both intermediate regions between the hot and cold slabs. The temperature gradient induces a thermo-diffusive mass flux throughout the system, thus inducing a compositional gradient. Once both temperature and composition gradients are established and show a linear behavior, Equation (2) can be applied to compute the Soret coefficient.
For a system containing 600 n-pentane and 600 ndecane molecules, exchange periods (N_{exch}) of 300 and 400 fs were tested with time steps of 1 and 2 fs. For the equimolar system with 1800 and 2400 molecules, exchange periods of 300 and 400 fs were tested. To take into account the influence of a thermostat on the Soret coefficient determination, the method was applied with a Berendsen thermostat with a weak-coupling time of 50 ps for equimolar systems with 1200 molecules and time step equal to 2 fs. This weak coupling factor was chosen following other studies in the literature (Zhang and Müller-Plathe, 2008; ^{Polyakov et al. 2008}) in order to minimize the influence of the thermostat on the system. We evaluated whether the addition of a thermostat, even with a weak coupling factor, can change the dynamics of the system in order to affect the determined Soret coefficient. The exchange periods of thermostated systems were 300 and 400 fs. The time step of 2 fs was chosen following studies published in the literature that used 2 fs or higher time steps (^{Simon et al., 1998}; ^{Perronace et al., 2002}; ^{Zhang and Müller-Plathe, 2005}; ^{Polyakov et al., 2008}). To avoid energy local minima, the systems were heated and cooled in the same way described for self-diffusion coefficient determinations. The final configurations obtained after the heating and cooling cycles were used as initial configurations for all systems. To evaluate the drift in total energy, the same initial configurations were used without applying the PeX algorithm. To obtain Soret coefficients, total simulation times of 48 ns were used in all cases. The local concentration in each slab was computed by counting the centers of mass of each species. Both concentration and temperature in each slab were sampled using intervals of 20 fs. It is important to point out that in BD-NEMD methods the assumption of linearity between fluxes and forces must hold, which also means that linear composition and temperature gradients must be obtained at the steady state condition.
RESULTS AND DISCUSSION
Densities of pure components and mixtures were determined via equilibrium NPT simulations as described previously and are compared to experimental results in Table 2, where the self-diffusion coefficients are also reported and compared with experimental and simulation data available in the literature. The standard deviations for experimental densities and self-diffusion coefficients were estimated by Equation (11) using uncorrelated data obtained by standard autocorrelation function analysis (Allen and Tildesley, 1987).
where is the standard deviation of the sample average with respect to the trajectory average, υ_{Y} is the standard deviation and M is the number of uncorrelated measurements of Y.
The results for density shown in Table 2 are in good agreement with data from the literature with deviations smaller than 1%. The self-diffusion coefficients obtained using the augmented TraPPE-UA force field were about 35% higher than the experimental values for n-pentane and 65% higher for n-decane. The discrepancies observed could be related to the fact that the TraPPE-UA force field was not parameterized taking into account transport properties. ^{Makrodimitri et al. (2011)} observed differences between simulated and experimental data ranging from 0.50 to 32%, depending on the condition, and considered this to be a good result. Comparing with results obtained with other force fields (anisotropic united atoms - AUA - ^{Toxvaerd, 1990}), the self-diffusion coefficients obtained here are 25% higher for n-pentane and 42% higher for n-decane. ^{Polyakov et al. (2008)} also observed values 20% higher or more than experimental data for self-diffusion coefficients using the TraPPE-UA force field. The same authors found a systematic deviation of 40% for the mutual diffusion coefficients for heptane-benzene mixtures. These results, when compared to simulations using different force fields, suggest that the TraPPE force field overestimates the mobility of the substances addressed above. Even though the results obtained for self-diffusion coefficients present considerably high deviations when compared to experimental data, it is important to evaluate the performance of the TraPPE-UA force field in the case of long alkanes.
Components | ρ (kg/m^{3}) | D (10^{-9} m^{2}/s) | ||||
---|---|---|---|---|---|---|
this work | exp. | this work | sim. (AUA1) |
sim. (AUA2) |
exp. | |
n-C5 | 621.5±0.1 | 621a | 7.59±0.17 | 612b | 618b | 5.62c |
n-C10 | 734.1±0.1 | 730a | 2.21±0.03 | 1.49b | 1.66b | 1.31c,1.38d |
n-C5 + n-C10 | 694.2±0.1 |
^{a}- ^{Dreisbach (1959)};
^{b}- ^{Mutoru et al. (2013)};
^{c}- ^{Douglas and McCall (1958)};
^{d}- ^{Ertl and Dullien (1973)}.
The Soret coefficients were determined for each condition using Equation (2). The first 3 ns were excluded from the calculations to avoid non-steady state perturbations. For each simulation, two values of composition and temperature gradients were obtained from the nine slabs above the hot slab and from the nine slabs below the hot slab. These values were used to evaluate and compare the methodology through the Soret coefficient determined for each region of the simulation box. Conditions used in each simulation and the Soret coefficients obtained are summarized in Table 3. The standard deviations were calculated using Equation (11).
Thermostat | Number of | Time step | N_{exch} | S_{T,upper} | S_{T,lower} | |ΔS^{T}| |
---|---|---|---|---|---|---|
molecules | (fs) | (fs) | (10. K^{-1}) | (10^{-3}. K^{-1}) | (10^{-3}. K^{-1}) | |
None | 1200 | 1 | 300 | 2.655±0.364 | 2.775±0.379 | 0.120 |
400 | 2.317±0.561 | 2.576±0.508 | 0.259 | |||
2 | 300 | 1.921±0.360 | 2.336±0.364 | 0.415 | ||
400 | 2.291±0.406 | 1.801±0.413 | 0.490 | |||
1800 | 1 | 300 | 2.631±0.324 | 2.820±0.362 | 0.189 | |
400 | 2.662±0.363 | 2.077±0.454 | 0.585 | |||
2400 | 1 | 300 | 2.677±0.383 | 2.335±0.309 | 0.342 | |
400 | 3.365±0.456 | 3.034±0.389 | 0.331 | |||
Berendsen | 1200 | 2 | 300 400 | 2.454±0.341 2.100±0.417 | 2.620±0.373 2.283±0.391 | 0.166 0.183 |
The differences between data using the upper and lower regions of the simulation box were higher for systems with 1200 molecules using a time step of 2 fs and for the system with 1800 molecules using an exchange period of 400 fs. Because the energy exchanged was nearly the same, simulations containing 1800 and 2400 molecules presented smaller gradients, which may explain why results with the larger number of molecules maintained similar fluctuations. However, the Soret coefficients determined with 2400 molecules are statistically similar (i.e., within the uncertainty) to those obtained with 1200 and 1800 molecules. This result confirmed that 1200 molecules are sufficient to describe the system. The system with 1800 molecules and N_{exch} equal to 300 presented almost the same gradient as the system containing 1200 molecules with N_{exch} equal to 400, showing Soret coefficient values next to each other, which reinforces the conclusion that a system with 1200 molecules is enough under the conditions studied. Furthermore, the maximum observed difference between all data is 1.564·10^{-3} K^{-1}.
According to ^{Zhang and Müller-Plathe (2005)}, a high imposed perturbation (i.e. frequent momentum exchanges) reduces the effect of noise in the calculations. However, exceedingly high perturbations must be avoided because they lead to very high temperature gradients in the simulation box. This can cause the fluid to behave like a supercooled liquid (i.e., be trapped in a free energy barrier) in the coldest slabs and, therefore, slow down the molecular motion in such slabs. Another implication of high perturbations is the risk of observing non-linear gradients caused by a too large departure from equilibrium. Another point that should be carefully analyzed is the energy drift observed during the simulations. A continuous increase in total energy causes an increase in the average temperature along the simulation run, which thus prevents the system from reaching a steady state. Finally, we note that the time needed for both temperature and composition profiles to reach a plateau will differ depending on the simulation conditions.
The effects of two different perturbation levels were studied, with momentum being exchanged after every 300 or 400 fs. According to the data reported in Table 3, higher deviations between statistically equivalent regions generally occur with N_{exch} equal to 400 fs for all studied systems. It suggests, in accordance with the literature (^{Zhang and Müller-Plathe, 2005}), that weaker perturbations imply a lower signal to noise ratio and thus differences between the Soret coefficients determined from both regions. This behavior can be partially explained by the small composition difference between the adjacent intermediate slabs, which can make the average composition in each slab to be of the same order of magnitude as the fluctuations. ^{Zhang and Müller-Plathe (2005)} observed similar trends. Nevertheless, the temperature and composition profiles obtained here showed linear behaviors for all systems, with determination coefficients (R^{2}) around 0.98 for composition gradients of the systems with time step equal to 1 fs and 0.93 to 0.97 for the systems with time step equal to 2 fs. The R^{2} for temperature gradients were higher than 0.99 for all studied systems.
Figure 1 shows the composition and temperature profiles for the system containing 1200 molecules, simulated with a time step of 1 fs and with an exchange period of 300 fs.
In Figure 1, composition and temperature of each slab were taken as an average between equivalent slabs from the upper and bottom regions. The error bars were computed through the evaluation of auto correlation functions of composition and temperature for each slab. The correlation time was around 2 ns for composition and 0.2 ns for temperature. The data used to calculate errors were averages of subsamples containing 2 ns of data for both composition and temperature, which are uncorrelated. The standard errors were estimated according to Equation (11).
As can be seen in Figure 1, there is a high temperature difference between hot and cold slabs, caused by the high driving forces imposed through unphysical momentum exchanges. The temperature differences were around 110 and 140 K for perturbations at every 400 and 300 fs, respectively. According to experimental data, the boiling points of n-pentane and n-decane are 309.2 and 447.2 K, respectively. The steady state temperatures reached values around 235 K in the cold slab and 370 K in the hot slab for the systems with N_{exch} equal to 300 fs. In the simulation with N_{exch} equal to 400 fs, the cold slab reached a temperature of 243 K and the hot slab of 360 K. Therefore, the hot region reached higher temperatures that in the vapor-liquid equilibrium region for an equimolar n-pentane/n-decane mixture. ^{Zhang and Müller-Plathe (2005)} observed similar behavior for a benzene/cyclohexane mixture at 324 K and mole fraction of 0.25. These authors observed a temperature of the hot slab around 360 K, which is higher than the boiling point of both substances (353.1 K for benzene and 353.7 K for cyclohexane), indicating a possible vapor-phase transition. To confirm the liquid behavior in the whole system ^{Zhang and Müller-Plathe (2005)} plotted the density profile in the simulation box and computed the mutual diffusion coefficient. For the cases studied here, the exchange period of 300 fs for 1200 molecules corresponds to the worse case, generating the highest temperature gradient. To be sure of the liquid-like behavior in the entire system, we plotted the steady state density profile in Figure 2.
Following ^{Zhang and Müller-Plathe (2005)}, we can affirm that the system behaved like a liquid over the entire domain of the system. To test the mobility of molecules, the intra-diffusion coefficients (also known as self-diffusion or tracer diffusion coefficients) for n-pentane and n-decane in the mixture were computed. To take this property into account, Equation (10) was used, but the average was taken only over the squared displacement of the center of mass of a given kind of molecule. The intra-diffusion coefficients for n-decane and n-pentane in the mixture were around respectively 3.2·10^{-9} and 5.2·10^{-9} m^{2}·s^{-1}, which agrees with the self-diffusion coefficients in the liquid state. Furthermore, the self diffusion coefficients of the components in the mixture are closer to each other than those of pure components, as pointed out by ^{Buhn et al. (2004)}, which reinforces that the TraPPE-UA force field correctly describes this trend of transport properties in mixtures. One possible explanation for why the system does not present vapor phase behavior in the hot slab is because we started the simulation with liquid using NVT or NVE ensemble. Since the density was fixed at a constant value and small differences in the box density (like a vapor phase) can cause an increase in pressure, this may be sufficient to prevent a phase transition.
In our simulations, very low temperatures were reached, i.e., temperatures next to the melting point of n-decane were observed. The experimental melting points of n-pentane and n-decane are respectively 143.4 and 243.3 K (^{Dreisbach, 1959}). In order to test whether the coldest slabs of the system were acting like a supercooled liquid and the corresponding energy barrier impeded molecular movements, preventing the motion of molecules in the z-axis, the self-diffusion coefficients on the x and y axes were compared to that in the z direction. Because the non-equilibrium system is subjected to a temperature gradient, some level of organization or preferential motion due to thermal diffusion may be observed. Therefore, because of the preferential motion addressed here, slight differences in self-diffusion coefficients of a given mixture component can occur in the z direction with respect to the x and y directions. For the systems of 1200 molecules and time step equal to 1 fs the differences varied from 15 to 25% for N_{exch} equal to 300 fs and 15% for N_{exch} equal to 400 fs. For systems with 1200 molecules and time step of 2 fs, differences were about 15 to 20% for N_{exch} equal to 300 and 400. To evaluate whether these differences are due to some level of organization or natural differences between the self-diffusion measured in the three directions, an equilibrium molecular dynamics simulation of a system with the same initial configuration of the system containing 1200 molecules and using a 1 fs time step was carried out during the same 48 ns in the NVE ensemble. The deviations between the self-diffusion coefficients calculated along the z-axis and the average of the x and y axes for the system not subjected to the driving force was between 10 and 15%. These differences are almost the same as that observed for the non-equilibrium system. In addition, self-diffusion coefficients determined for both kinds of systems (equilibrium and non-equilibrium) were of the same order of magnitude, which suggests that the molecules are not jamming together due to supercooled liquid-like behavior. Even with large perturbations and high temperature gradients, all systems simulated here behaved like a liquid in the sense of self-diffusion coefficients and densities. For the systems with a coupled thermostat, similar behaviors were observed.
The increase in total energy of the system during the simulation run was evaluated for all conditions studied. The relative energy drift is shown in Figure 3 for all simulations with 1200 molecules.
To compute this effect, we determined the time-dependent ratio of total energy, β_{E}(t), the overall ratio of total energy, β_{E}', and the temperature variation during the simulation (ΔT), defined as follows,
where ET and T are, respectively, the instantaneous total energy and temperature of the system, and the symbols <>_{0} and <>_{f} denote, respectively, averages taken at the first 2 ns and at the last 2 ns of the simulation run. The values of β_{E}' and ΔT are shown in Table 4, while β_{E}(t) is shown in Figure 3.
Results shown in Figure 3 suggest that the stronger the perturbations the higher the energy drift during the simulation. In addition, higher energy drifts were observed for time steps of 2 fs. When the total energy of the system rises, the overall temperature increases.
Thermostat | time-step | N_{exch} | β'_{E} | ΔT |
---|---|---|---|---|
(fs) | (fs) | (%) | (K) | |
None | 1 | 300 | 1.79 | 1.6 |
400 | 1.06 | 0.6 | ||
2 | 300 | 7.47 | 7.3 | |
400 | 4.68 | 4.6 | ||
Berendsen | 2 | 300 400 |
-0.44 -0.06 |
-0.1 0.0 |
According to results shown in Table 4, higher increases in overall temperature were found for systems with higher time steps and larger perturbations, except, as expected, for the system in which Berendsen's thermostat was used. According to ^{Zhang (2006)}, the increase in energy which leads to a warm up of the system is caused by round-off or cut-off noise. However, when a momentum exchange occurs, the next step will deviate from Newtonian mechanics, causing a perturbation in the system and thus not conserving the total energy as supposed and recommended. From Equations (3) to (7), it is possible to perceive that all terms will depend in the last instance on the position vectors of each atom. Thus, when a perturbation takes place, the directions of velocity vectors of the participating atoms will change and the Hamiltonian will not be conserved in the next integration step. When this occurs, stretching out of angles, bond lengths, and dihedrals and overlapping may take place, which can cause an unphysical addition of energy to the overall system. This information gains a little more strength if we compare the data obtained in this work for the systems described in Table 4, except for simulations using the thermostat. When the time step is higher, the next step of integration will displace more the atoms, favoring stretching and overlapping. Also, doing perturbations more frequently, we may favor overlaps and stretching, which can be the reason for higher energy drifts and higher overall temperatures. As already mentioned, this effect can input more fluctuations and may complicate the convergence of thermal and composition profiles. A reasonable solution to this problem could be the addition of a Berendsen thermostat with a gentle coupling so as not to be a source of instability. The Berendsen thermostat rescales the velocity of the atoms of the system, which maintain the directions of velocity vectors and, as a first approximation, do not cause large disturbances in the dynamic of the system. To evaluate the use of the thermostat, the convergence of composition and temperature profiles in each region of the simulation box and the standard errors in Soret coefficient must be made in conjunction. Figures 4 and 5 show the time evolution of thermal and composition profiles for systems with 1200 molecules. The lines are averages accumulated since the outset of the production run, while each circle is an average taken from the preceding 2 ns.
According to the results shown in Figure 4 and 5, it is possible to conclude that the temperature profiles converge on a very short time scale, before the time shown on the graph (the graphs start at 3 ns), while composition profiles require a long time scale, at least 20 ns. In addition, these calculations show that the composition gradient fluctuates much more than the temperature gradient. The high fluctuation in the composition gradient could be related to molecular size. Because the composition profiles are calculated by counting the number of centers of mass in a slab, a large portion of a molecule could be in the adjacent slab while the center of mass is in another. Therefore, due to the molecular size and its movement, the centers of mass could jump between two slabs, contributing to an increase in fluctuations. These results are consistent with results reported by ^{Perronace et al. (2002)}. The linearity of composition and temperature profiles were evaluated through the R^{2}, obtaining respectively values higher than 0.96 and 0.99 for composition and temperature. Another interesting observation is that the R^{2} of the composition profiles reached acceptable values (higher than 0.96; data not shown) before the profile reached the convergence value, underestimating the composition gradient and thus the Soret coefficient value. In this way, it is possible to conclude that one of the major source of errors in the Soret coefficient may be due to fluctuations of the composition gradient or to an inadequate simulation time to achieve the plateau of the composition gradient. In this way, it is necessary to perform large simulations and continuously check the convergence. A large simulation is going to decrease the uncertainties of the Soret coefficient, which is high due to the fluctuations of the composition profile. The addition of the Berendsen thermostat showed no visible change in the behavior of the curves of Figures 4 and 5 and the same observations stated above are applied to the thermostated simulations.
The standard errors in the Soret coefficient in each simulation were determined using the statistically equivalent slabs of both regions as a replicate to perform the linear regression of temperature and composition profiles. The errors in composition and temperature in each slab were taken into account to compute profiles and their uncertainties. The standard error of the Soret coefficient was computed with the composition and temperature profile uncertainties using standard error propagation methods (^{Bevington and Robinson, 2003}).
Soret coefficients and uncertainties are presented in Table 5 and compared with experimental data and other simulation results from the literature.
Thermostat | t-step/N_{exch} (fs) | S_{T}(10^{-3}. K^{-1}) | |||
---|---|---|---|---|---|
This worka | BD-NEMD (HeX) | S-NEMD | exp. | ||
None | 1/300 | 2.70±0.24 | |||
1/400 | 2.42±0.38 | 2.91±0.45b | 3.59±0.83b | 3.27±0.23c | |
2/300 | 2.18±0.26 | ||||
2/400 | 2.07±0.28 | ||||
Berendsen | 2/300 | 2.58±0.26 | 3.07d | ||
2/400 | 2.16±0.28 |
^{a}- R-NEMD, S_{T} ± 2·υ_{<Y>};
^{b}- ^{Perronace et al. (2002)}, 300 K;
^{c}- ^{Perronace et al. (2002)}, 300 K;
^{d}- ^{Mezquia et al. (2014)}, 298 K.
The Soret coefficient for n-decane in n-pentane determined here and compared with experimental data showed deviations around 17.5 to 37%. The standard errors determined in this work were in general smaller than those of other studies from the literature, which may be due to the higher simulation time and smaller time steps employed compared to the work of ^{Perronace et al. (2002)}, where the simulation time was 32 ns with a time step of 4 fs. It is possible to notice a difference between the absolute values of the Soret coefficient obtained by ^{Perronace et al. (2002)} and those obtained here. Those differences and deviations compared to the experimental data are possibly due to the different NEMD methodology employed and to the force field used to model the n-alkanes. While ^{Perronace et al. (2002)} used the SKS force field (^{Smit et al., 1995}) with the HeX (BD-NEMD) and the S-NEMD algorithms, we tested the R-NEMD algorithm (also a BD-NEMD algorithm) of ^{Müller-Plathe (1997)} with the TraPPE-UA force field (^{Martin and Siepmann, 1998}). Better agreement with experimental data was found for the system with N_{exch} equal to 300 fs and a time step of 1 fs, showing deviations of 17.5% from experimental data, which can be related to a smaller energy drift. In addition, a 1 fs time step will lead to smaller bond, angle, and dihedral stretching caused by the energy exchange steps. Results using a time step of 2 fs without the coupling of the Berendsen thermostat presented the worst agreement compared to literature data, showing deviations ranging from 25 to 42% compared with other simulations and 33 to 37% compared with experimental data. These deviations may be related to the observed energy drift, causing a continuous increase of perturbations and changing the dynamics of the system at intermediate slabs. Although a deviation of 17.5% could be considered high, the difficulty in obtaining cross coefficients, either by experiments or by molecular dynamics, means that the result can be considered to be in good agreement with the literature, consistent with what ^{Perronace et al. (2002)} have addressed. The use of the thermostat showed higher discrepancies than the system using a time step of 1 fs for different values of N_{exch}. In addition, the deviations between experimental data and S_{T} obtained from simulations using a 1 fs time step and N_{exch} equal to 300 were smaller than others, which suggest that the use of a small time step can provide better S_{T} values. In principal, the use of a thermostat could control the energy drift and improves results. However, this effect did not clearly occur. It suggests that the use of a thermostat may cause no expected deviations, which may be due to the change of dynamics of the system on intermediate slabs, since the velocity rescaling is applied on the entire domain.
CONCLUSION
The Soret coefficients of an n-decane/n-pentane equimolar mixture were determined using the Boundary Driven Non-Equilibrium Molecular Dynamics method. The main contribution of this work is to show a systematic analysis of the method, testing the effect of different time steps, energy exchange, and thermostat control over the steady state conditions, the Soret coefficient and its fluctuation.
A high amount of energy exchange produced a high driving force in the simulation box, where we observed temperature differences (i.e., the temperature difference between the hot and cold slabs) ranging from 110 to 140 K. Although this driving force generated temperatures in the vapor-liquid equilibrium region for the n-pentane/n-decane mixture, the system behaved as a liquid. The liquid-like behavior of the system was confirmed by evaluating the self-diffusion coefficients of each component and the density profile of the mixture over the simulation box. According to the results, the self-diffusion of both substances remained on the order of magnitude expected for the liquid phase in all directions of the simulation box. A small deviation was observed for self-diffusion along the z-axis, which may be due to the energy flux in this direction. However, when compared to a system with the same initial configuration and simulated using equilibrium molecular dynamics, the deviations were of the order of 10%. The density profile also confirms the liquid-like behavior in all slabs of the box.
We studied the use of a thermostat to calculate the Soret coefficients and the influence of the time step and the amount of energy exchange. Higher overall energy drifts were observed for simulations running without the thermostat with time steps equal to 2 fs. This energy drift led to an overall temperature variation of around 7 K. Although the Soret coefficient deviations obtained from these simulations were not smaller, the energy drift may not be the only cause of deviations. Compared with experimental data, Soret coefficients obtained from simulations with a time step of 1 fs presented better results. The energy drift for this system was not significant since variations around 1 K were observed. For simulations using the thermostat, the Soret coefficient for energy-exchange periods of 300 fs are in agreement with literature data. However, the worst results were observed for simulations running without the thermostat with exchange periods of 400 fs and a time step equal to 2 fs. Additionally, at least 20 ns are needed to acquire an ideal convergence of the composition profile. Althought smaller time steps generate reliable results without the use of a thermostat, a higher computational effort must be employed, which can be a problem.
We have shown that temperature profiles converge on a short time scale (around 2 ns for our system), while the composition profiles require a very long time scale (around 20 ns for our system). Therefore, to obtain reliable Soret coefficients, we need not only a good force field, but also large simulations in terms of the number of particles and simulation time. Large simulations smooth out composition-profile fluctuations and decrease uncertainties.