A new boundary driven NEMD scheme for heat and particle diffusion in binary mixtures

Abstract

A new method for Boundary Driven Non-Equilibrium Molecular Dynamics (BD-NEMD) simulation is presented. It allows the simultaneous imposition of both a constant temperature and concentration gradient. By varying the strength of the imposed gradients and obtaining stationary values of the mass and energy flux, this new technique can be used to extract the Onsager transport coefficients directly. We demonstrate the setup and measure explicitly the Soret and Dufour effects for a binary mixture of Argon and Krypton. We compare the Soret coefficient computed using this new scheme against the experimentally inspired protocol where the system has no mass current. However, this method is limited to only an estimate of the Soret coefficient rendered inefficient due to the poor signal-to-noise ratio and cannot give the transport coefficients explicitly. The new technique allows for accurate measurements of non-zero flux values and, if run in the linear response domain, the full set of transport coefficients can be extracted directly. We also speculate about how this scheme could be applied to more complex systems with applications for polymer separations and biological activities.

GRAPHICAL ABSTRACT

Keywords:

• Nonequilibrium molecular dynamics
• transport properties
• Onsager coefficients
• binary mixtures
• concentration gradient

1. Introduction

Thermophoresis or thermal diffusion is the process by which mixtures of multiple particle types separate due to an imposed temperature gradient, oftentimes referred to as the Soret effect. Thermophoresis is relevant in a variety of processes, ranging from large solar pond systems where a high salt gradient prevents natural convection amongst the layers helping to trap heat at the bottom of the pond; all the way to small biopolymers in which it is used to help characterise protein interactions in fluids [1]. Experimental estimation of this effect has proved to be challenging as flow induced by thermal gradients are typically of smaller order than those from concentration gradients. This is where molecular dynamics (MD) experiments have risen as a means to study this phenomenon, aiming initially at the calculation of transport coefficients by means of Green–Kubo formulas or by applying indirect (mechanical) perturbations [2, 3]. Prior work has examined the impacts of subjecting the system to a thermal gradient and waiting until mass equilibrium (${\mathbf{J}}_1=0$) is reached and the particles have naturally separated [4–7]. When this occurs the ratio of the resultant concentration profile over the thermal gradient gives the Soret coefficient ($S_T$), a measure of the ratio between the thermal diffusion coefficient ($D^{\mathrm{\prime }}$) to the binary diffusion coefficient (D) and allows insight into the key driving force behind particle separation. Given similar species, the lighter one will typically conglomerate in higher temperature regions and exhibit negative thermophoretic behaviour as it moves from the cold to hot region. One study has attempted to measure the rate at which this separation occurs by the sudden imposition of a temperature gradient and measurement of the transient response [8]. A larger value of $S_T$ indicates regimes in which thermal diffusion becomes more dominant, particularly relevant when liquids are near a glassy regime [9].

Conversely, the Dufour effect is the energy flux of the system created by the imposition of a concentration gradient. It has not been as well characterised in the literature, since it is typically negligible for liquid mixtures [10], but there is experimental evidence for its effect in carbon tetrachloride and cyclohexane mixtures [11]. In addition, the timescale for the establishment of a stationary thermal gradient is generally much quicker than that of a stationary concentration gradient, so studying this effect experimentally is constrained by the establishment of a constant concentration gradient during which time other effects such as buoyancy can corrupt precise measurements.

In this paper, we introduce a new schema to produce Boundary-Driven Nonequilibrium Molecular Dynamics (BD-NEMD). We attempt to study the impacts of imposing both thermal and concentration gradients, $\mathrm{\nabla }T$ and $\mathrm{\nabla }{x}_1,$ respectively, by means of fixing the kinetic energy and mole fraction at the edges and aim to answer two main questions.

1. Can a stationary state be established by imposing that only one gradient be different from zero, and from it, by monitoring the flux values, can an accurate and direct measurement of the Onsager transport coefficients (i.e. the thermal diffusion $D^{\mathrm{\prime }}$, the binary diffusion coefficient D and the thermal conductivity κ) be obtained?

2. Can we identify the regime under which linearity is observed and so provide a reliable estimate for the transport coefficients to be compared with prior simulation protocols?

In Section 2, we describe in detail the model and the simulation setup needed to drive the system out of equilibrium to enforce a stationary condition with the presence of linear temperature and concentration gradients. In Section 3, we recall from de Groot and Mazur [12] the constitutive relationships between the energy and mass fluxes and the corresponding temperature and chemical potential gradients in the linear regime for heat and mass transport in binary mixtures, which define the desired transport coefficients. In Section 4, we describe the technique to extract their precise values from the simulations and report our estimate of the of these values with comparison to results obtained using different protocols. Lastly, in Section 5, we draw some conclusions from the present study and discuss how this simulation setup can be applied to more complex molecular systems.

2. Model, method, and simulation setup

For a large class of irreversible thermal-diffusive phenomena produced by a range of experimental conditions the flows in a purely phenomenological theory are expressed as a linear function of the thermodynamic forces (1) ${\mathbf{J}}_i=\sum _j{L}_{ij}{\mathbf{X}}_{\mathbf{j}}$(1) where $\mathbf{J}$ and $\mathbf{X}$ are any of the Cartesian components of the independent fluxes and thermodynamic forces appearing in an expression for entropy production [12].

The matrix of coefficients, $\mathbf{L}$, is comprised of quantities called the Onsager transport coefficients and the relations written above are the corresponding phenomenological equations. In this work, restricted to heat and particle diffusion phenomena in binary, non-reactive, mixtures, we will model the external condition producing diffusion and transport of heat. In this case the relevant thermodynamic forces can be defined as the gradients of $\frac{1}{T}$ and $\frac{\mu }{T}$, where T is the temperature and μ an appropriate chemical potential, function of temperature, pressure and concentration, to be defined later on.

In this section, we introduce a procedure that allows for simultaneous imposition of a temperature and concentration gradient through BD-NEMD. The system we choose to study is a binary mixture of Argon (Ar) and Krypton (Kr) atoms which have been well characterised in the literature previously and allows for more accessible comparisons.

The system is modelled in reduced units using the Ar atom as reference, as characterised in Reith and Muller-Plathe [5]. The atomic interactions are described through 6-12 LJ pair potentials with $m_{\mathrm{Ar}}=1,{ϵ}_{\mathrm{Ar}}=1,{\sigma }_{\mathrm{Ar}}=1$ and $m_{\mathrm{Kr}}=2.1,{ϵ}_{\mathrm{Kr}}=1.39,{\sigma }_{\mathrm{Kr}}=1.07$. Mixing interactions between the two are constructed by means of the Lorentz–Berthelot convention (${\sigma }_{12}=\frac{{\sigma }_1+{\sigma }_2}{2}$ and ${ϵ}_{12}=\sqrt{{ϵ}_1{ϵ}_2}$). We use a timestep (ts) of 0.005 units (corresponding for real Ar to $1.09\times {10}^{-14}$  s) and a cut-off distance for energy calculation of $R_c=3{\sigma }_{\mathrm{Kr}}=3.21$.

The box size is $22.47\times 22.47\times 40.28$ with periodic boundary conditions are applied only in the x and y directions. In the z-dimension, we impose a confining external ‘wall’ potential that mimics the interaction of a molecule with the planar surface of a semi-infinite crystal made of LJ particles [13]. The z-dimension of the box is chosen so to achieve, on average, a number density n = 0.725 in the bulk. We choose to employ the same potential of the wall for both Ar and Kr particles, $U^{(w)}(d)=(\frac{C}{45d^9}-\frac{C}{6d^3})+U_{\mathrm{shift}}$ for $0<d\le d_c$ and $U^{(w)}(d)=0$ for $d>d_c$, where d is the particle distance from the wall, $C=\pi n_w$ and we set the wall number density $n_w$ equal to the average density of the fluid mixture. The values of the potential shift $U_{\mathrm{shift}}$ are chosen in such a way that $U^{(w)}(d_c)=0$. Setting the cut-off distance equal to the position of the minimum of the un-shifted potential $d_c=\sqrt[6]{2/5}$, the interaction of the particles with the wall becomes purely repulsive with both potential and force continuously going to zero with increasing distance $d\to d_c$. Shown in Figure 1 is a schematic diagram of the simulation setup together with a snapshot of the system with Kr atoms, in purple, larger in size than  Ar atoms, in green.

Figure 1. (a) The subdivision of the simulation box into the different regions, the inner zone, white background, is where the fields and currents are measured, while the regions where the thermostats act on particle velocities are highlighted by different shades of blue for the lower temperature $T_C$, on the left, and red for the higher temperature $T_H$, on the right, with particle exchanges taking place in the darker zones to fix the concentrations to the values $x_{1C}$ and $x_{1H}$, respectively, as described in the text. Periodic boundary conditions are only applied in the x, y-directions, while in the z-direction the system is confined by repulsive walls, dark grey area, placed at the edges of the box. (b) A snapshot is rendered with VMD [14] from the simulation data with temperature gradient $\mathrm{\nabla }T>0$ and molar fraction gradient $\mathrm{\nabla }{x}_1>0$ between the two reservoirs along the z-axis (blue arrow on-axis compass). Note a separation of Kr species (larger size, purple) particles gathering in the higher temperature region, $T_H$, as a result of the imposed concentration gradient. The empty spaces in the z-direction are the regions rarely explored by the particles where the wall potential is highly repulsive.

All simulations are performed with a fixed total number $N=\mathrm{14,112}$ of particles, on average with a 50/50 split of Ar to Kr atoms. The initial configuration is created by randomly assigning 7056 Ar and 7056 Kr particles with equal probability to the sites of an fcc lattice. After the particles have been placed and their identities assigned, the system is equilibrated by running a 200,000 ts constant temperature dynamics by means of the standard velocity Verlet algorithm coupled to an Andersen thermostat, where random velocity sampling takes place with a 10 ts frequency.

Thermal and concentration gradients are applied to the system in the BD-NEMD production runs. This is accomplished by first removing for convenience two small regions of width $0.88$ at the edges, most of the time empty of particles, and splitting the remaining volume, which is of dimension $21{\sigma }_{\mathrm{Kr}}\times 21{\sigma }_{\mathrm{Kr}}\times 36{\sigma }_{\mathrm{Kr}}$, along the z-axis into two reservoir regions at the edges and a central ‘bulk’ region where particles dynamics is not perturbed, but simply determined by the forces arising from intermolecular interactions and following Newton's equations of motion. Each of these regions is split into slices of width $1{\sigma }_{\mathrm{Kr}}$, with constant volume $V^{(s)}=540.244$ in reduced units. The two reservoir regions consist of 6 slices at both edges while the bulk region consists of 24 such slices in total, and in each slice the temperature, partial density and energy fields, together with their related currents are calculated at each time step.

In the reservoir regions, temperature and concentration can be imposed individually in each of the 6 slices on both sides of the simulation box. Temperature control is achieved by rescaling the kinetic energy separately for each atomic type.

Concentration control is achieved by a switch of particle identity at random in the reservoir, from Ar to Kr or vice versa, with mass changing while keeping fixed the momentum value, in such a way to attain within ±1 particle the desired concentration value. The origin of the latter limitation in the ability of fixing the concentration is simply due to the fact that the particle count in the reservoir slices cannot be a fractional number. Consider in particular the case of equal concentrations: it can be enforced only if the total number of particles in the slice is even. If this number happens to be odd choosing the type of the last unpaired particle will introduce an excess of one particle type. In the general case of arbitrary concentration, the chosen value to enforce will always fall between two fractional numbers differing at most by one unit in the integer in the numerator: the algorithm will then assign the unknown type for such last remaining particle at random with a probability equal to the desired concentration value. In this way, the concentration in each reservoir slice will asymptotically attain the desired imposed value on average.

In the production run, the 6 slices of each reservoir region are further split into two groups. While temperature control is applied in all the 6 slices, concentration control is only applied to the outermost, with respect to the bulk, 3 slices. This permits having a thermostatted buffer region of width equal to the potential cut-off distance $R_c$ between the fixed concentration reservoirs and the bulk region in which no particle swap takes place. The single-step switch of the particle identity does, in fact, introduce a numerical noise in the integration algorithm, particularly in the case of the smaller Ar to larger Kr transformation where a spike in the force value can result in poor integration of the dynamics at fixed timestep value. Fortunately for the Ar/Kr swap, the differences between the LJ potential parameters are small and with the imposed constant kinetic energy temperature control this problem is minimised. We found that avoiding direct contact between the concentration reservoirs and the bulk region is sufficient to obtain the desired accuracy in the dynamics of the total system. We did observe that without temperature control numerical errors introduced by the sudden particle identity switches, although very small for Ar/Kr on the local scale, results in an unwanted and sizeable overall temperature rise on the very long time scales required to get accurate averages, persisting even over several millions timesteps.

The thermostat regions are dubbed $T_H/T_C$ and are outlined in red(right)/blue(left) in Figure 1. To impose a temperature gradient we set $T_H=T+\mathrm{\Delta }T/2$ and $T_C=T-\mathrm{\Delta }T/2$. This establishes a constant thermal gradient with a value very close to $\mathrm{\Delta }T/\mathrm{\Delta }z$ around the mean system temperature $(T=1)$ across the inner 24 slices of width $\mathrm{\Delta }z=25.68$.

Velocity rescaling is performed following a similar procedure to the one outlined in a previous paper [15], where at each timestep, in every thermostatted slice, the total kinetic energy of the individual particle types are calculated after removing the group centre of mass motion from the individual particle velocity. Temperature control is performed in such a way that, individually for each particle type and as a result also globally, the total momentum in each slice is kept unchanged by the velocity rescaling procedure. The scaling factor is determined according to the instantaneous temperature $T_{\alpha }$ for species $\alpha$ computed from the kinetic energy ($E_{K\alpha }$) of the $N_{\alpha }$ particles of species $\alpha$, with $\alpha =Kr$ for species 1 or $\alpha =Ar$ for species 2. Locally for each slice: (2) $T_{\alpha }=\frac{2E_{K\alpha }}{3k_B(N_{\alpha }-1)}$(2) with (3) $E_{K\alpha }=\sum _{i=1}^{N_{\alpha }}\frac{{({\mathbf{p}}_{\alpha ,i}-m_{\alpha }{\mathbf{v}}_{\alpha })}^2}{2m_{\alpha }}$(3) where ${\mathbf{v}}_{\alpha }=\frac{1}{N_{\alpha }{m}_{\alpha }}\sum _{i=1}^{N_{\alpha }}{\mathbf{p}}_{\alpha ,i}$ is the instantaneous centre of mass velocity of particles of species $\alpha ,{\mathbf{p}}_{\alpha ,i}$ is the instantaneous individual particle momentum and the index i runs over the $N_{\alpha }$ particles within the inspected slice. We observed on the very long time scale of our simulations, particularly for the $\mathrm{\Delta }T=0$ cases, that thermostats using a velocity rescaling algorithm have a tendency to return a local temperature slightly above the desired value in the bulk region. The difference remains indeed small, below ${10}^{-4}$ parts, and it is hardly visible with respect to the statistical noise in the local averages even after averaging over several tens of million timesteps. This shift in temperature effect does not change the temperature gradient as it is constant across the bulk and we ignore its effect. All the values reported for temperature and concentration gradients, as well used in the calculation of transport properties, are those effectively measured by averaging over the particles in the bulk region.

All calculations were performed at fixed zero total momentum for the simulated system in the x- and y-directions, while the component in the z-direction fluctuates around zero only due to the interaction with the containing walls, since the two combined procedures leave the momentum unchanged at each time step. We choose to impose temperature gradients and concentration gradients in the ranges $\mathrm{\Delta }T=[0.,0.1]$ and $\mathrm{\Delta }{x}_1=[0.,0.1]$, where each simulation is run for several tens of million timesteps. Long equilibration times are needed to reach stationary flux values and are on the order of $t=30\times {10}^3$ LJ time units as shown in Figure 2 by plotting the running average of the mass flux with ($\mathrm{\Delta }T,\mathrm{\Delta }{x}_1)=(0.1,0.1$).

Figure 2. Running average of the mass flux for $(\mathrm{\Delta }T,\mathrm{\Delta }{x}_1)=(0.1,0.1)$ with error shown by the shaded region. Relaxation to the stationary value occurs after $t=30\times {10}^3$ units with minor long-living fluctuations around the mean value afterwards as shown by the horizontal green line. Shown in the inset is the signal-to-noise ($\frac{{\mu }^2}{{\sigma }^2}$) ratio which encouragingly continues to steadily climb with increasing statistics.

From the simulations, we extract the particle count, mass flux, heat flux, and temperature at each timestep and coarse grain them by averaging in each of the slices over 200 timesteps. The typical coarse-grained result along the z-dimension from a single simulation run are shown in Figure 3 for the case $(\mathrm{\Delta }T,\mathrm{\Delta }{x}_1)=(0.1,-0.045)$, corresponding (see Section 4) to ${\mathbf{J}}_1\to 0$ in the bulk region. The darker colours at the edges in the panel highlight the different reservoir regions with the adjoining slightly lighter colours indicating the additional thermostat slices.

Figure 3. A plot of the temperature (a), mol fraction (b), mass flux of Kr (c), and energy flux (d) as a function of position along the z-direction measured in units of slice width, i.e. ${\sigma }_{\mathrm{Kr}}$, with $(\mathrm{\Delta }T,\mathrm{\Delta }{x}_1)=(0.1,-0.045)$. Each point represents the average over a slice and $T_H/T_C$ are indicated by the red(right)/blue(left) shaded regions. The reservoir region in which the mol fraction is controlled consists of only the first and last 3 slices in $T_H/T_C$ as described in the text and are shaded darker. The values reported in (c) and (d) are for the flux averaged across the entire bulk region, neglecting the slices in direct contact with the thermostats and where the faint green line shows the zero level. The vertical green lines in panels (a) and (b) show the equidistant point in the bulk from the reservoirs and it may be observed that the concentration profile (b) is not symmetric in shape for the two species across the bulk and indicates a preference of the heavier Kr towards the colder regions.

3. Theory

The mass and energy flux drive the evolution of a system and the Onsager transport coefficients, notably $D^{\mathrm{\prime }}$, $D$, and $\kappa$, measure the strength of the response to imposed gradients. Note that, due to the Curie principle, viscosity belongs to independent transport equations and we will not be concerned with it in the present manuscript. It should be highlighted that these coefficients are only meaningful in the linear regime for which they are derived, where it has been assumed that the system is only at the first-order expansion away from equilibrium. Outside of this regime, the higher orders of the gradients become relevant and need to be considered explicitly.

For a binary system, one can write the constitutive relationships, Equation (1(1) ${\mathbf{J}}_i=\sum _j{L}_{ij}{\mathbf{X}}_{\mathbf{j}}$(1) ), relating these coefficients and the gradients to the energy and mass flux. Microscopically, spatial variations in the chemical potential and temperature are what drive this evolution and give rise to non-zero fluxes. The relevant chemical potential of the system ($\mu ={\mu }_1-{\mu }_2$ with ${\mu }_1,{\mu }_2$, respectively, the chemical potential of Kr and Ar, defined per unit mass) is taken to be a function of temperature, pressure, and the molar fraction of component $1x_1$, i.e. $\mu =\mu (T,p,x_1)$. We can exploit the fact that in the linear regime the average temperature and pressure remain constant and use Gibbs–Duhem and chain rule (see Appendix) to write $\mathrm{\nabla }\mu =\frac{\rho }{m_2{x}_{2}n}{\left(\frac{\mathrm{\partial }{\mu }_1}{\mathrm{\partial }{x}_1}\right)}_{T,p}\mathrm{\nabla }{x}_1$. We choose to take the specific constitutive relations from equation (229) and (230) on page 276 of chapter XI of de Groot–Mazur book [12] expressed in terms of the molar fraction $x_1$ in the following format: (4) $\begin{array}{rl}-{\mathbf{J}}_q& =[\kappa \mathrm{\nabla }T+x_{1}n{m}_1{\mu }_{11}T{D}^{\mathrm{\prime }\mathrm{\prime }}\mathrm{\nabla }{x}_1]\end{array}$(4) (5) $\begin{array}{rl}-{\mathbf{J}}_1& =\frac{m_1{m}_{2}n}{\rho }[x_1{x}_{2}n{D}^{\mathrm{\prime }}\mathrm{\nabla }T+nD\mathrm{\nabla }{x}_1]\end{array}$(5) where ${\mathbf{J}}_1$ is the mass flux of Kr, ${\mathbf{J}}_q$ the heat flux, ${\mu }_{11}=(\frac{\mathrm{\partial }{\mu }_1}{\mathrm{\partial }{x}_1}{)}_{T,p}$ the partial derivative of the Kr chemical potential at constant temperature and pressure, $\rho =m_1{x}_{1}n+m_2{x}_{2}n$ the mass density with $n$ the number density, $x_1=(\frac{n_{\mathrm{Kr}}}{n_{\mathrm{Ar}}+n_{\mathrm{Kr}}})$ the mole fraction of Kr, and $x_2=(\frac{n_{\mathrm{Ar}}}{n_{\mathrm{Ar}}+n_{\mathrm{Kr}}})$ the mole fraction of Ar. Onsager reciprocal relations express the equality of the cross coefficients and relate the mass flow induced by a temperature gradient to the energy flux coming from a concentration gradient and manifest itself in the above equations by setting $D^{\mathrm{\prime }}=D^{\mathrm{\prime }\mathrm{\prime }}$.

To estimate the heat flux in Equation (4(4) $\begin{array}{rl}-{\mathbf{J}}_q& =[\kappa \mathrm{\nabla }T+x_{1}n{m}_1{\mu }_{11}T{D}^{\mathrm{\prime }\mathrm{\prime }}\mathrm{\nabla }{x}_1]\end{array}$(4) ) for a given slice we can use equivalently, as shown by Ercole et al. [16], the energy flux defined by Irving and Kirkwood [17] as (6) $\begin{array}{rl}{\mathbf{J}}_e^{(s)}& =\frac{1}{V^{(s)}}[\sum _{i\in \mathrm{slice}}{e}_i({\mathbf{v}}_i-{\mathbf{v}}^{(s)})\\ & +\frac{1}{2}\sum _{i\in \mathrm{slice}}\sum _{j\ne i}({\mathbf{F}}_{ij}\cdot {\mathbf{r}}_{ij})({\mathbf{v}}_i-{\mathbf{v}}^{(s)})]\end{array}$(6) where $e_i$ is the sum of potential and kinetic energy, ${\mathbf{v}}_i$ the velocity for atom $i$, ${\mathbf{v}}^{(s)}$ the barycentric velocity of the tagged slice, ${\mathbf{F}}_{ij}$ the force on atom $i$ from atom $j$, and ${\mathbf{r}}_{ij}$ the vector joining the positions of atom $i$ and atom $j$. Note that index $j$ extends to particles outside the slice with the slice volume given by $V^{(s)}=21{\sigma }_{\mathrm{Kr}}\times 21{\sigma }_{\mathrm{Kr}}\times {\sigma }_{\mathrm{Kr}}$.

We estimate the mass flux in Equation (5(5) $\begin{array}{rl}-{\mathbf{J}}_1& =\frac{m_1{m}_{2}n}{\rho }[x_1{x}_{2}n{D}^{\mathrm{\prime }}\mathrm{\nabla }T+nD\mathrm{\nabla }{x}_1]\end{array}$(5) ) by summation of the species momentum in each slice (7) ${\mathbf{J}}_1^{(s)}=\frac{1}{V^{(s)}}\sum _{i\in \mathrm{slice}}({\mathbf{p}}_i-m_i{\mathbf{v}}^{(s)})$(7) where ${\mathbf{p}}_{\mathbf{i}}$ is the momentum of Kr atoms in slice $i$. By definition ${\mathbf{J}}_1^{(s)}+{\mathbf{J}}_2^{(s)}=0$. The choice of the fixed repulsive wall at the edges permits the system to remain immobile and allows the total momentum to fluctuate, throughout the simulation run, around the initial zero value. We specifically choose to impose along the z-direction confining walls rather than periodic boundary conditions as with the latter, we observed a permanent non-zero barycentric mass flow in the central bulk region compensated by an overall opposite flux created in the two reservoir regions since they are in direct contact through the periodic boundaries. This would contradict the hypothesis of zero bulk barycentric mass flow under which de Groot–Mazur's equations were derived. Finally, we choose to present the results in terms of the heavier particle 1 (Kr), note that only the sign of the Soret coefficient $S_T$ would change if referring to Ar.

We make in passing the remark that an additional way to write the constitutive relations has been formulated by Trimble and Deutch [18] and followed, not completely consistently, by other studies for transport coefficients [7, 15]. This results in defining the coefficients in a slightly different notation. However, conversion between the two can be found and is presented in the Appendix. Note that for the purpose of this manuscript, we present the results using Equations (4(4) $\begin{array}{rl}-{\mathbf{J}}_q& =[\kappa \mathrm{\nabla }T+x_{1}n{m}_1{\mu }_{11}T{D}^{\mathrm{\prime }\mathrm{\prime }}\mathrm{\nabla }{x}_1]\end{array}$(4) ) and (5(5) $\begin{array}{rl}-{\mathbf{J}}_1& =\frac{m_1{m}_{2}n}{\rho }[x_1{x}_{2}n{D}^{\mathrm{\prime }}\mathrm{\nabla }T+nD\mathrm{\nabla }{x}_1]\end{array}$(5) ), following de Groot–Mazur formalism.

In order for the theory to be correctly applied, there are two requirements that must be fulfilled, first that the system reaches an asymptotic state in which the fluxes reach a uniform stationary value across the central bulk region and second that the strength of the imposed gradients remains in the linear regime region. Once that is assessed, we also need to make a proper estimate of the errors associated with the statistical averages, which provide the values of the fluxes. We can postpone to the next section the checks for both requirements, while we conclude this section discussing the algorithm we will use to evaluate the errors.

We remind that values sampled along a trajectory generated by Molecular Dynamics are often correlated in time and for error analysis we need to estimate correctly the number of independent contributions to the averages. The error associated to the flux value is computed via a block averaging approach which allows for independent samples to be obtained by uniformly partitioning the trajectory into different blocks, giving $N_b$ blocks of equal length ${\tau }_b$, as described in Flyvbjerg and Petersen [19]. The $N_b$ necessary to determine when the independence of block averages occurs is taken through observation of when the estimated variance of the block averages of length ${\tau }_b$ divided by $N_b$ reaches a plateau. An example of this procedure is shown in Figure 4. The error of the overall average is then computed through the following set of equations: (8) $\begin{array}{rl}{J}_b& =\frac{1}{{\tau }_b}\sum _{i=1}^{{\tau }_b}{J}_{b,i},b=1,2,\dots ,N_b\end{array}$(8) (9) $\begin{array}{rl}\overline{J}& =\frac{1}{N_b{\tau }_b}\sum _{b=1}^{N_b}\sum _{i=1}^{{\tau }_b}{J}_{b,i}\end{array}$(9) (10) $\begin{array}{rl}{s}_{{\tau }_b}^2& =\frac{1}{N_b-1}\sum _{b=1}^{N_b}(J_b-\overline{J}{)}^2\end{array}$(10) (11) $\begin{array}{rl}{s}_J^2& =\frac{s_{{\tau }_b}^2}{N_b}\end{array}$(11) where $J_{b,i}$ is either the mass or energy flux value at an instant i in time within the particular time block b, $s_{{\tau }_b}^2$ is the, generally biased, estimate of the variance of $J_b$ and is a function of the time length ${\tau }_b$ of the blocks. For high frequency sampling rates, the estimate is biased due to the values being correlated in time. The unbiased estimated variance of $\overline{J}$ is taken as the value of $s_{{\tau }_b}^2/N_b$ as a function of ${\tau }_b$ in the range where it remains approximately constant.

Figure 4. The estimated standard deviation for $J_1$ from the block analysis of three independent trajectory segments each of ${10}^7$ timesteps for the ($\mathrm{\Delta }T=0.1,\mathrm{\Delta }{x}_1=0$) condition. Values for ($\sqrt{s_{{\tau }_b}^2/N_b}$) are plotted against the number of timesteps in each block ${\tau }_b$, where $N_b$ is the number of blocks and $N_b{\tau }_b={10}^7$. The region between 100 and 10,000 steps shows a relatively flat profile, which implies that there the block averages can be considered as independent samples. Going to shorter intervals of time, the decrease can be traced back to increasing correlations while at the other end the number of blocks $N_b$ becomes too small to provide a reliable estimate of the variance of the sample.

4. Results

We start by showing that both requirements mentioned in Section 3 are met, namely that the flux values are relaxing to stationary values and that the gradients applied are kept in the linear regime. The relaxation of the mass flux of particle 1 is shown in Figure 2 for the case of $(\mathrm{\Delta }T,\mathrm{\Delta }{x}_1)=(0.1,0.1)$. We see that the presence of long-lived fluctuations after an initial ‘fast’ relaxation requires averages over time intervals lasting longer than O(${10}^3$), but that the value does plateau and promisingly that the signal-to-noise ratio climbs steadily with time, testifying the achievement of asymptotic stationary conditions.

The first check for linearity can be quickly performed by comparing the values of the fluxes with the corresponding gradients when, while keeping $\mathrm{\Delta }T=0$, we double the imposed concentration boundary value from $\mathrm{\Delta }{x}_1=0.05\text{ to }\mathrm{\Delta }{x}_1=0.10$. Resulting in $\mathrm{\nabla }{x}_1^{(0.05)}=(1.540\pm 0.004)\times {10}^{-3}$ with ${\mathbf{J}}_1^{(0.05)}=(-7.73\pm 0.25)\times {10}^{-5}$ and $\mathrm{\nabla }{x}_1^{(0.1)}=(2.980\pm 0.004)\times {10}^{-3}$ with ${\mathbf{J}}_1^{(0.1)}=(-14.67\pm 0.29)\times {10}^{-5}$. Thus giving a ratio of $\mathrm{\nabla }{x}_1^{(0.1)}/\mathrm{\nabla }{x}_1^{(0.05)}={\mathbf{J}}_1^{(0.1)}/{\mathbf{J}}_1^{(0.05)}=1.9$. We recall that while for the temperature gradient $\mathrm{\nabla }T\approx \mathrm{\Delta }T/\mathrm{\Delta }z$ holds to a high degree of accuracy, this is not the case for the concentration gradient $\mathrm{\nabla }{x}_1$ in the presence of a non-zero temperature gradient in the bulk region. Therefore, the general test for linearity consists in showing that ${\mathbf{J}}_1(\mathrm{\nabla }{T}_a+\mathrm{\nabla }{T}_b,\mathrm{\nabla }{x}_{1,a}+\mathrm{\nabla }{x}_{1,b})={\mathbf{J}}_1(\mathrm{\nabla }{T}_a,\mathrm{\nabla }{x}_{1,a})+{\mathbf{J}}_1(\mathrm{\nabla }{T}_b,\mathrm{\nabla }{x}_{1,b})$ for any combination of the individual imposed gradients. Results from combining the imposed boundary conditions $\mathrm{\Delta }T=0.0,0.05,\mathrm{\Delta }{x}_1=0.0,0.05$, are analysed in Table 1 and show that, with small relative errors, the condition for linearity is well satisfied.

Table 1. Check of linearity using the measured values of the temperature gradient, the energy flux, the molar fraction gradient and the mass flux for component 1, reported in LJ units.

	${10}^3\cdot \mathrm{\nabla }T$	${10}^3\cdot {\mathbf{J}}_e$	${10}^3\cdot \mathrm{\nabla }{x}_1$	${10}^5\cdot {\mathbf{J}}_1$
(a)	1.861 (0.004)	−8.50 (0.03)	−0.121 (0.004)	−3.56 (0.25)
(b)	0.003 (0.004)	0.07 (0.03)	1.540 (0.004)	−7.73 (0.25)
(sum)	1.864 (0.006)	−8.32 (0.4)	1.419 (0.006)	−11.29 (0.35)
(a+b)	1.853 (0.004)	−8.39 (0.03)	1.395 (0.004)	−10.73 (0.25)
(rel.err)	0.6	0.8	1.7	5

The value obtained summing the fluxes for case (a) and case (b) are compared with the values measured with the case (a+b) with boundary values chosen to impose simultaneously both the temperature and the molar fraction gradients imposed in (a) and (b) separately. The fluxes satisfy the linear regime hypothesis within the numerical accuracy, with the resulting gradients adding up consistently. The relative error $\frac{(\mathrm{a}+\mathrm{b})-(\mathrm{sum})}{(\mathrm{sum})}$ is also included. Note that, since we can impose the boundary values but not the gradients, the check cannot be performed with exactly the same total gradients.

With the requirements for linearity and stationarity of the fluxes and gradients fulfilled, we next move to the calculation of the transport coefficients with the BD-NEMD method. We can see from the constitutive relations, Equations (1(1) ${\mathbf{J}}_i=\sum _j{L}_{ij}{\mathbf{X}}_{\mathbf{j}}$(1) ), that if we impose all but one gradient to be zero, we can measure one by one the elements of the matrix of coefficients while varying the strength of the non-zero gradient. So, at least in principle, in our procedure we have to set one of the two gradients equal to zero to extract directly the two transport relevant coefficients from the ratio of the resultant measured fluxes and the strength of the other gradient. As said, we only impose boundary values, so this is possible only after ensuring that the observed value of the first gradient is actually close to zero in the central bulk region ; otherwise, the residual value of the non-vanishing gradient needs to be taken into account explicitly. However, while we can impose the condition $\mathrm{\nabla }T=0$ with any desired degree of accuracy, in the presence of a non-zero temperature gradient a nonnegligible molar fraction gradient $\mathrm{\nabla }{x}_1\ne 0$ is unavoidable in the bulk region. Therefore we start by extracting the transport coefficients from Equations (4(4) $\begin{array}{rl}-{\mathbf{J}}_q& =[\kappa \mathrm{\nabla }T+x_{1}n{m}_1{\mu }_{11}T{D}^{\mathrm{\prime }\mathrm{\prime }}\mathrm{\nabla }{x}_1]\end{array}$(4) ) and (5(5) $\begin{array}{rl}-{\mathbf{J}}_1& =\frac{m_1{m}_{2}n}{\rho }[x_1{x}_{2}n{D}^{\mathrm{\prime }}\mathrm{\nabla }T+nD\mathrm{\nabla }{x}_1]\end{array}$(5) ) in the case when $\mathrm{\nabla }T=0$ and $\mathrm{\nabla }{x}_1\ne 0$. With both ${\mathbf{J}}_1$ and ${\mathbf{J}}_e$ known, the binary diffusion coefficient D and the product ${\mu }_{11}{D}^{\mathrm{\prime }}$ can be computed by rearranging the two equations to give (12) $\begin{array}{rl}{\mu }_{11}{D}^{\mathrm{\prime }}& =\frac{-{\mathbf{J}}_e}{m_1{x}_{1}nT\mathrm{\nabla }{x}_1}\end{array}$(12) (13) $\begin{array}{rl}D& =\frac{-{\mathbf{J}}_1\rho }{m_1{m}_2{n}^2\mathrm{\nabla }{x}_1}\end{array}$(13) The gradients $\mathrm{\nabla }T$ and $\mathrm{\nabla }{x}_1$ are the ones calculated from the temperature and molar fraction profiles in the bulk region between $(z=9{\sigma }_{\mathrm{Kr}}$ to $z=33{\sigma }_{\mathrm{Kr}})$, as shown in Figure 5(c) for the case $\mathrm{\Delta }T=0$ and $\mathrm{\Delta }{x}_1=0.1$. We can observe that the corresponding temperature profile is indeed flat, satisfying the condition $\mathrm{\nabla }T=0$. One can see that in the bulk region, there is a slight temperature increase with respect to the reservoir regions, but it remains uniform and therefore gives no contribution to the flux so that D and the product $D^{\mathrm{\prime }}{\mu }_{11}$ can be readily estimated from the above equations.

Figure 5. The temperature (top) and concentration (bottom) profiles where panels (a), (b) and (c) are for conditions $\mathrm{\Delta }T,\mathrm{\Delta }{x}_1=\{(0.1,0.1),(0.1,0),(0,0.1)\},$ respectively. Each point represents the average over a slice and regions $T_H/T_C$ are indicated by the red(right)/blue(left) shaded areas. The reservoir region in which the mol fraction is controlled consists of only the first and last 3 slices in $T_H/T_C$ as described in the text and are shaded darker. From (b), one can see that while the concentration is equal at both edges, it is not enough to force the profile in the bulk to be uniform and a small gradient emerges, with opposite sign to the temperature gradient.

We next impose the conditions $\mathrm{\Delta }{x}_1=0$ in order to obtain the heat conductivity κ and the thermal diffusion coefficient $D^{\mathrm{\prime }}$. However, in this case, as it can be observed from Figure 5(b) for $\mathrm{\Delta }T=0.1$, a residual, stationary, non-zero gradient in concentration is present in bulk region, despite the imposed equimolar ratio in the reservoirs, and its contribution to the mass flux of component 1 cannot be ignored. Thus in order to compute the latter two transport coefficients correctly we need to account for both the temperature and the concentration gradients from the observed profiles Figure 5(b) and plug their values in the following expressions for κ and $D^{\mathrm{\prime }}$, derived again from Equations (4(4) $\begin{array}{rl}-{\mathbf{J}}_q& =[\kappa \mathrm{\nabla }T+x_{1}n{m}_1{\mu }_{11}T{D}^{\mathrm{\prime }\mathrm{\prime }}\mathrm{\nabla }{x}_1]\end{array}$(4) ) and (5(5) $\begin{array}{rl}-{\mathbf{J}}_1& =\frac{m_1{m}_{2}n}{\rho }[x_1{x}_{2}n{D}^{\mathrm{\prime }}\mathrm{\nabla }T+nD\mathrm{\nabla }{x}_1]\end{array}$(5) ): (14) $\begin{array}{rl}\kappa & =\frac{-1}{\mathrm{\nabla }T}[{\mathbf{J}}_e+x_{1}n{m}_1{\mu }_{11}T{D}^{\mathrm{\prime }}\mathrm{\nabla }{x}_1]\end{array}$(14) (15) $\begin{array}{rl}{D}^{\mathrm{\prime }}& =\frac{-1}{x_1{x}_{2}n\mathrm{\nabla }T}[\frac{{\mathbf{J}}_1\rho }{m_1{m}_{2}n}+nD\mathrm{\nabla }{x}_1]\end{array}$(15) Summarising, from the simulation with $\mathrm{\nabla }T=0$ an accurate computation of $D$ and the product $D^{\mathrm{\prime }}{\mu }_{11}$ can be performed. Then these results can be used in a second simulation with $\mathrm{\nabla }T\ne 0$ to complete the calculation of all coefficients for the remaining values of $D^{\mathrm{\prime }}$ and κ. Note that also ${\mu }_{11}$ can now be computed.

With the flux values and necessary profiles in hand and following the procedure outlined above, we recollect the transport coefficients from the present study in Table 2.

Table 2. Heat and mass transport coefficients from the BD-NEMD protocol.

$D$	$D^{\mathrm{\prime }}$	$\kappa$	${\mu }_{11}$	$S_T=D^{\mathrm{\prime }}/D$
0.0501 (0.001)	0.087 (0.003)	4.603 (0.009)	−0.83 (0.34)	1.7 (0.1)

All values reported are in LJ units with error given in the parenthesis.

This numbers, reported to the Onsager coefficients, are compared in Table 3 with those calculated by Sarman and Evans [2], where the coefficients were extracted from Green-Kubo relations for a system of Ar–Kr using the Evans–Cummings scheme [20], an algorithm to drive NEMD simulations by equating the time derivative of the Hamiltonian to a product of the flux and a mechanical force which gives the exact Green-Kubo relations in the linear limit.

Table 3. Comparison of the Onsager coefficients between our results (BD-NEMD) and those from the mechanical perturbation scheme by Sarman and Evans [2] for a system of 1024 particles of Ar-Kr.

Onsager coefficient	BD-NEMD ($T=1.0,n=0.725$)	Ref. [2] ($T=0.965,n=0.7137$)
$L_{qq}$	4.60 (0.01)	4.24 (0.01)
$L_{1q}$	0.0214 (0.0007)	0.0182 (0.001)
$L_{11}$	0.0191 (0.0075)	0.0174 (0.0005)

Rather than the values of the transport coefficients $\kappa ,D^{\mathrm{\prime }},D$ and of ${\mu }_{11}$, which are not calculated directly in the mechanical perturbation scheme, we compare with the values of the Onsager coefficients $L_{qq},L_{1q}=L_{q1}$ and $L_{11}$ calculated from our data (see Appendix for the details of definitions).

We see relative agreement between the two methods despite a small difference in our thermodynamic points and it should also be considered that the Sarman and Evans results were obtained with less statistics and for smaller system size.

As a further check of the new method, we also compare the Soret coefficient defined as $S_T=D^{\mathrm{\prime }}/D$ [12] with the one computed using the protocol following studies used in [5, 15], which consists of imposing a temperature gradient at the edges of the simulation box and allowing the system to relax to mass equilibrium by the establishment of a natural concentration gradient. When ${\mathbf{J}}_1=0$, the Soret coefficient can be taken as the ratio of the observed concentration gradient over the observed thermal gradient [12], i.e. $S_T=-\frac{1}{x_1{x}_2}\frac{\mathrm{\nabla }{x}_1}{\mathrm{\nabla }T}$.

We explored the behaviour of our system resulting from the imposition of only the difference $\mathrm{\Delta }T=0.1$ in temperature between the two reservoir regions in order to estimate the naturally occurring concentration gradient. In practice, we found that the relaxation time needed to establish the asymptotic concentration and temperature profiles, i.e. to achieve mass equilibrium in the central bulk region, can be exceedingly long and that, correspondingly, the average mass flux remains small but non-zero even after very long times, reaching a final value that, divided by a number of the same order of magnitude (the product of the binary diffusion coefficient and thermal gradient) lead to an important correction for $S_T$, as we now show.

Making use of the binary diffusion coefficient D measured with the BD-NEMD scheme, we take the residual value of ${\mathbf{J}}_1$ into account to obtain the improved expression (16) $S_T=-\frac{1}{x_1{x}_2}\left(\frac{\mathrm{\nabla }{x}_1}{\mathrm{\nabla }T}\right)-\frac{\rho }{m_1{m}_2{x}_1{x}_2{n}^2}\left(\frac{{\mathbf{J}}_1}{D\mathrm{\nabla }T}\right)$(16) and use it to compute the Soret coefficient using the gradients in two equivalent ${\mathbf{J}}_1\approx 0$ conditions. In fact, thanks to the ability to fix both the concentration and the temperature in the reservoirs to our liking, we can pursue an alternative route to establish the same stationary flux ${\mathbf{J}}_1=0$ in a driven condition by imposing at the same time as $\mathrm{\Delta }T=0.1$ an appropriate boundary value for the concentration, $\mathrm{\Delta }{x}_1=-0.045$. Comparison of the gradients and fluxes obtained with these two alternative routes is summarised in Table 4. Both setups show a similar deviation from zero mass flux, ${\mathbf{J}}_1=-0.2$ in one case and ${\mathbf{J}}_1=+0.2$ in the other, so that, after treatment, both give $S_T\approx 1.7$.

Table 4. Results of the measured fluxes and gradients in the bulk for the experimentally inspired protocol where $J_1\approx 0$.

Condition	${10}^5\times J_1$	${10}^3\times J_e$	${10}^3\times \mathrm{\nabla }T$	${10}^3\times \mathrm{\nabla }{x}_1$	$S_T$
Natural	−0.2 (0.3)	−16.85 (0.03)	3.678 (0.004)	−1.512 (0.004)	1.69 (0.07)
Driven	+0.2 (0.1)	−16.89 (0.02)	3.691 (0.002)	−1.612 (0.002)	1.70 (0.04)

The results dubbed Driven were obtained imposing $\mathrm{\Delta }T=0.1$ and $\mathrm{\Delta }{x}_1=-0.045$, fixing both the thermal and molar fraction gradients, while Natural refers to values with only $\mathrm{\Delta }T=0.1$ imposed, fixing the temperature gradient and allowing the system to relax naturally until mass equilibrium is reached and a stationary concentration gradient is created.

Indeed, the values from Table 4 give for $\frac{-\mathrm{\nabla }{x}_1}{x_1{x}_2\mathrm{\nabla }T}$ two nonconsistent values, that is 1.747 (0.003) when we drive the system imposing both gradients and 1.644 (0.003) when we allow the system to reach the mass equilibrium state in the presence of a temperature gradient alone. Consistency is regained with Equation (16(16) $S_T=-\frac{1}{x_1{x}_2}\left(\frac{\mathrm{\nabla }{x}_1}{\mathrm{\nabla }T}\right)-\frac{\rho }{m_1{m}_2{x}_1{x}_2{n}^2}\left(\frac{{\mathbf{J}}_1}{D\mathrm{\nabla }T}\right)$(16) ) including the correction due to the non-zero mass flux, as shown in Table 4. In this way, we do also get a consistent agreement with the value $S_T=D^{\mathrm{\prime }}/D=1.7(0.1)$ previously calculated from the transport coefficients from BD-NEMD.

In Figure 6, we show the behaviour of the running average of ${\mathbf{J}}_1$ in the driven condition. We can see that even after ${10}^5$ LJ-units, the average mass flux is small but non-zero, reaching a final value ($J_1=0.2\times {10}^{-5}$) and that frustratingly the error associated with it is not insignificant in comparison ($s_J=0.1\times {10}^{-5}$ ). The increased uncertainty in the value of the Soret coefficient calculated from Equation (16(16) $S_T=-\frac{1}{x_1{x}_2}\left(\frac{\mathrm{\nabla }{x}_1}{\mathrm{\nabla }T}\right)-\frac{\rho }{m_1{m}_2{x}_1{x}_2{n}^2}\left(\frac{{\mathbf{J}}_1}{D\mathrm{\nabla }T}\right)$(16) ) stems from such non negligible error associated to the evaluation of ${\mathbf{J}}_1$ as highlighted in the inset of Figure 6.

Figure 6. Running average of the mass flux with the shaded region showcasing the error with $\mathrm{\Delta }T=0.1$ and $\mathrm{\Delta }{x}_1=-0.045$ to achieve mass equilibrium. Note that even after $t=40\times {10}^3$, there is still a measurable mass flux value and that at $t=100\times {10}^3$ the uncertainty is as large as the signal. This is shown more prominently in the inset where the signal-to-noise is given over time and contrasting to the inset of Figure 2 where the SNR rose monotonically with simulation time.

Note that we showed that not considering the second term with ${\mathbf{J}}_1\approx 0$ results in underestimating the error on the value of $S_T$. While it has been reported that the Soret coefficient has an extreme sensitivity to even a minor perturbation in the imposed gradient, our results point rather to the large uncertainties on $S_T$ coming from the difficulty to achieve the condition ${\mathbf{J}}_1=0$ with sufficient accuracy to justify neglecting its value. If this factor cannot be safely ignored, then knowledge of the binary diffusion coefficient must be obtained. In this respect, the benefit of the new BD-NEMD technique becomes apparent in that we can derive the value of the Soret coefficient directly from $D^{\mathrm{\prime }}/D$ computing the two transport coefficients in stationary non-equilibrium conditions, with an accuracy that can be improved at will by increasing the statistics or exploring in the linear regime a range of gradients conditions that can be plotted against the fluxes and the coefficients found through a simple linear fitting procedure.

5. Discussion

In the present manuscript, we have introduced a method that permits to set at will both a concentration and a temperature linear gradient by imposing simultaneously appropriate conditions at the edges of the simulation box and we used it to make a direct simulation of transport processes. The method has been applied here to binary mixtures but it can be straightforwardly generalised to multicomponent systems.

These are moreover real thermodynamics forces as opposed to the mechanical perturbations used by Evans et al. [2, 3]. The bulk region in between the edges is run with Newtonian dynamics and it is from this bulk that we can obtain stationary mass and energy flux values, varying the imposed temperature or concentration gradients while keeping the other gradient fixed to zero. This gives the desired transport coefficients using the constitutive relations presented by de Groot–Mazur and valid in the linear regime. We compare the results for $S_T$ to the one obtained by running with a protocol where the system establishes a mass equilibrium by creating a concentration gradient countering the temperature polarisation of the system so that $J_1\approx 0$. The two methods give consistent results, but to get complete alignment we need to consider the explicit residual value of $J_1$ from this setup even for apparently very small order of magnitudes in order to correctly determine the uncertainty. Because of this, two other benefits of running the new BD-NEMD technique over this protocol are that (i) because no flux signal is perfectly zero, the stationary values give a higher signal-to-noise ratio and (ii) we can explicitly compute all the various transport coefficients and not be constrained to only extract their relative ratios or have to rely on equilibrium sampling techniques. This opens the possibility for exploration into the driving forces for glassy systems or purification where the thermal and concentration gradients are of similar magnitude and it is desired to compute the specific contributions to the magnitude of the flux evolution. Moreover, by fixing one species type and varying the other, this technique allows for a range of transport properties to be determined and thus allows for creation of more suitable mixtures. This is particularly useful in areas involving nanofluids [21] in specific applications like refrigerant cooling where thermal diffusion is thought to play a more significant role and screening possible candidates based on observed properties can be difficult.

We highlight that while the technique is broad and applicable to any system, the results themselves can have a dependence on the thermodynamic point making comparison and estimation for different systems challenging. The results we present are only for a system of moderate packing density ($n\approx 0.725$) and only for the mild change of Ar to Kr and vice versa. In the future, we plan to explore more dramatic changes such as Na${}^{+}$ to K${}^{+}$, as the depolarisation of this specific type of concentration gradient is responsible in biological system signalling [22]. Other avenues include dramatically varying the particle size such as in the separation of polymers. This may require a modification to the current simulation setup as a sudden particle type change in the reservoir regions can lead to instabilities caused by large steric overlap. This could potentially be circumvented by introducing a different time integration scheme in the reservoir regions or by the gradual growing or shrinking as is commonly done in free energy perturbation methods [23–25]. We hope to apply this technique in a wide variety of applications where non-equilibrium states due to imposed gradients can be found and help elucidate the strength of the elusive cross-coefficient terms.

Acknowledgments

We are delighted to dedicate this small contribution in non-equilibrium molecular simulation to Mike Klein who has just achieved his 80 years. We wish he will continue to share with us his human and scientific wisdom for many years to come. This work is supported by grants from the Welch Foundation F-1896 and from the NIH GM059796.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work is supported by grants from the Welch Foundation F-1896 and from the Foundation for the National Institutes of Health (NIH) GM059796.

 

Appendix. Reconciling Trimble–Deutch to de Groot–Mazur

In the present manuscript, we attempted to make comparisons primarily to the results in  [3, 7, 15]. However, we found that the coefficients defined in the constitutive equations in  [18] relating the flux to the transport equations do not match properly the ones given by de Groot and Mazur [12]. We know that Onsager relations state that the mass flow per unit of imposed temperature difference should be equivalent to the energy flow per unit of imposed concentration gradient [26]. Under the assumption of microscopic reversibility this means that the cross coefficients $L_{q1},L_{1q}$, when defined with respect to the thermodynamic forces $\mathrm{\nabla }\frac{1}{T},\mathrm{\nabla }\frac{\mu }{T}$, should be identical. However, with the introduction of additional coefficients in the equations that is not immediately evident.

We first start by recalling the constitutive relations for a binary system in equations (212 and 213) from page 274 in Chapter XI of Non-Equilibrium Thermodynamics [12]. In order to simplify the comparison with the corresponding expressions of Trimble and Deutch [18], without having to include too many additional constant conversion factors that can ultimately hide the important part of the result, we choose here to follow de Groot–Mazur to work, rather than with the previously used molar fractions, in terms of the mass fractions $c_1=\frac{m_1{x}_1}{m_1{x}_1+m_2{x}_2}$ and $c_2=\frac{m_2{x}_2}{m_1{x}_1+m_2{x}_2}$. The relation within the two alternative descriptions if fully given in Chapter XI of  [12]. For a binary system the relevant chemical potential appearing in the expression of the thermodynamic force $\mathrm{\nabla }(\mu /T)$ driving mutual diffusion is given by $\mu ={\mu }_1-{\mu }_2$ where ${\mu }_1$ and ${\mu }_2$ are the specific chemical potential of components 1 and 2, respectively, and at constant $T$ and$p$ satisfy the Gibbs–Duhem relation $c_1(\mathrm{\nabla }{\mu }_1{)}_{T,p}+c_2(\mathrm{\nabla }{\mu }_2{)}_{T,p}=0$. The following expression for the gradient of μ at constant temperature in terms of the gradient of the mass fraction $c_1$, $\mathrm{\nabla }({\mu }_1-{\mu }_2)=\mathrm{\nabla }{\mu }_1/c_2={\mu }_{11}^c\mathrm{\nabla }{c}_1/c_2$ can then be derived, introducing here for brevity the quantity ${\mu }_{11}^c=(\frac{\mathrm{\partial }{\mu }_1}{\mathrm{\partial }{c}_1}{)}_{T,p}$. Note that the quantity ${\mu }_{11}=(\frac{\mathrm{\partial }{\mu }_1}{\mathrm{\partial }{x}_1}{)}_{T,p}$ introduced in Section 3 is the equivalent abbreviation when using ${\mu }_1$ as a function of the molar fraction $x_1$ rather than of the mass fraction $c_1$. We can transform from one to the other by means of the relation ${\mu }_{11}^c\mathrm{\nabla }{c}_1={\mu }_{11}\mathrm{\nabla }{x}_1$, that can be easily derived from the above definition of $c1$ and $c2$.

The transport coefficients $D^{\mathrm{\prime }}$ and $D^{\mathrm{\prime }\mathrm{\prime }}$ are then defined (A1) $\begin{array}{rl}{J}_e{}^{(\mathrm{dm})}& =-[\frac{L_{\mathrm{qq}}^{(\mathrm{dm})}}{T^2}\mathrm{\nabla }T+\frac{L_{\mathrm{q}1}^{(\mathrm{dm})}}{c_{2}T}{\mu }_{11}^c\mathrm{\nabla }{c}_1]\\ & =-\kappa {}^{(\mathrm{dm})}\mathrm{\nabla }T-\rho c_{1}T{D}^{\mathrm{\prime }\mathrm{\prime }}{\mu }_{11}^c\mathrm{\nabla }{c}_1\end{array}$(A1) (A2) $\begin{array}{rl}{J}_1{}^{(\mathrm{dm})}& =-[\frac{L_{1\mathrm{q}}^{(\mathrm{dm})}}{T^2}\mathrm{\nabla }T+\frac{L_{11}^{(\mathrm{dm})}}{c_{2}T}{\mu }_{11}^c\mathrm{\nabla }{c}_1]\\ & =-\rho c_1{c}_2{D}^{\mathrm{\prime }}\mathrm{\nabla }T-\rho D{}^{(\mathrm{dm})}\mathrm{\nabla }{c}_1\end{array}$(A2) and showed to satisfy the equality $D^{\mathrm{\prime }\mathrm{\prime }}=D^{\mathrm{\prime }}$ as a result of Onsager relation $L_{q1}^{(\mathrm{dm})}=L_{1q}^{(\mathrm{dm})}$. The superscripts $(\mathrm{dm})$ has been used here to identify the quantities defined using the formalism of de Groot-Mazur, while the superscript $(\mathrm{td})$ will be used in the following to distinguish similar, but somehow differing, quantities defined in Trimble and Deutch formalism. In equations (90) and (91) on page 167 of Trimble and Deutch [18] the Onsager relation is imposed in terms of somewhat different phenomenological coefficients $L_{\mathrm{q}1}^{(\mathrm{td})}=L_{1\mathrm{q}}^{(\mathrm{td})}$ and the equations are then recast defining different transport coefficients ${\kappa }^{(\mathrm{td})},D_T,D^{(\mathrm{td})}$: (A3) $\begin{array}{rl}{J}_e{}^{(\mathrm{td})}& =-[\frac{L_{\mathrm{qq}}^{(\mathrm{td})}}{T}\mathrm{\nabla }T+\frac{L_{\mathrm{q}1}^{\mathrm{td}}}{c_2}{\mu }_{11}^c\mathrm{\nabla }{c}_1]\\ & =-\kappa {}^{(\mathrm{td})}\mathrm{\nabla }T-\frac{D_T}{c_2}{\mu }_{11}^c\mathrm{\nabla }{c}_1\end{array}$(A3) (A4) $\begin{array}{rl}{J}_1{}^{(\mathrm{td})}& =-[\frac{L_{1\mathrm{q}}^{(\mathrm{td})}}{T}\mathrm{\nabla }T+\frac{L_{11}^{(\mathrm{td})}}{c_2}{\mu }_{11}^c\mathrm{\nabla }{c}_1]\\ & =-\frac{D_T}{T}\mathrm{\nabla }T-\frac{D^{(\mathrm{td})}}{c_2}{\mu }_{11}^c\mathrm{\nabla }{c}_1\end{array}$(A4) Finally, by comparing, term by term, the expressions in front of the imposed gradients in Equations (A1(A1) $\begin{array}{rl}{J}_e{}^{(\mathrm{dm})}& =-[\frac{L_{\mathrm{qq}}^{(\mathrm{dm})}}{T^2}\mathrm{\nabla }T+\frac{L_{\mathrm{q}1}^{(\mathrm{dm})}}{c_{2}T}{\mu }_{11}^c\mathrm{\nabla }{c}_1]\\ & =-\kappa {}^{(\mathrm{dm})}\mathrm{\nabla }T-\rho c_{1}T{D}^{\mathrm{\prime }\mathrm{\prime }}{\mu }_{11}^c\mathrm{\nabla }{c}_1\end{array}$(A1) ) and (A2(A2) $\begin{array}{rl}{J}_1{}^{(\mathrm{dm})}& =-[\frac{L_{1\mathrm{q}}^{(\mathrm{dm})}}{T^2}\mathrm{\nabla }T+\frac{L_{11}^{(\mathrm{dm})}}{c_{2}T}{\mu }_{11}^c\mathrm{\nabla }{c}_1]\\ & =-\rho c_1{c}_2{D}^{\mathrm{\prime }}\mathrm{\nabla }T-\rho D{}^{(\mathrm{dm})}\mathrm{\nabla }{c}_1\end{array}$(A2) ) with those in Equations (A3(A3) $\begin{array}{rl}{J}_e{}^{(\mathrm{td})}& =-[\frac{L_{\mathrm{qq}}^{(\mathrm{td})}}{T}\mathrm{\nabla }T+\frac{L_{\mathrm{q}1}^{\mathrm{td}}}{c_2}{\mu }_{11}^c\mathrm{\nabla }{c}_1]\\ & =-\kappa {}^{(\mathrm{td})}\mathrm{\nabla }T-\frac{D_T}{c_2}{\mu }_{11}^c\mathrm{\nabla }{c}_1\end{array}$(A3) ) and (A4(A4) $\begin{array}{rl}{J}_1{}^{(\mathrm{td})}& =-[\frac{L_{1\mathrm{q}}^{(\mathrm{td})}}{T}\mathrm{\nabla }T+\frac{L_{11}^{(\mathrm{td})}}{c_2}{\mu }_{11}^c\mathrm{\nabla }{c}_1]\\ & =-\frac{D_T}{T}\mathrm{\nabla }T-\frac{D^{(\mathrm{td})}}{c_2}{\mu }_{11}^c\mathrm{\nabla }{c}_1\end{array}$(A4) ) the conversion between the two formalisms can be fully expressed as (A5) $\begin{array}{rl}\kappa {}^{(\mathrm{dm})}& =\kappa {}^{(\mathrm{td})}(\mathrm{heat}\mathrm{conductivity})\end{array}$(A5) (A6) $\begin{array}{rl}D{}^{(\mathrm{dm})}& =\frac{{\mu }_{11}^c}{\rho c_2}D{}^{(\mathrm{td})}(\mathrm{binary}\mathrm{diffusion}\mathrm{coefficient})\end{array}$(A6) (A7) $\begin{array}{rl}{D}^{\mathrm{\prime }}& =\frac{D_T}{\rho c_1{c}_{2}T}(\mathrm{thermal}\mathrm{diffusion}/\mathrm{Dufour}\mathrm{coefficient})\end{array}$(A7) resulting in the following alternative expressions for the Soret Coefficient $S_T$ (A8) $S_T=\frac{D^{\mathrm{\prime }}}{D^{(\mathrm{dm})}}=\frac{D_T}{c_1{\mu }_{11}^{c}T{D}^{(\mathrm{td})}}$(A8) Note finally that, following de Groot and Mazur [12], for the case ${\mathbf{J}}_1=0$ the more commonly used expression for $S_T$ can be derived from Equation (A2(A2) $\begin{array}{rl}{J}_1{}^{(\mathrm{dm})}& =-[\frac{L_{1\mathrm{q}}^{(\mathrm{dm})}}{T^2}\mathrm{\nabla }T+\frac{L_{11}^{(\mathrm{dm})}}{c_{2}T}{\mu }_{11}^c\mathrm{\nabla }{c}_1]\\ & =-\rho c_1{c}_2{D}^{\mathrm{\prime }}\mathrm{\nabla }T-\rho D{}^{(\mathrm{dm})}\mathrm{\nabla }{c}_1\end{array}$(A2) ) (A9) $S_T=-\frac{1}{c_1{c}_2}\frac{\mathrm{\nabla }{c}_1}{\mathrm{\nabla }T}=-\frac{1}{x_1{x}_2}\frac{\mathrm{\nabla }{x}_1}{\mathrm{\nabla }T}.$(A9) where the right-hand side, in terms of molar fractions, is the actual expression we used when discussing our results in Section 4 and follows immediately by means of the equality $\mathrm{\nabla }{c}_1/c_1{c}_2=\mathrm{\nabla }{x}_1/x_1{x}_2$.