Predicting and Rationalizing the Soret Coefficient of Binary Lennard-Jones Mixtures in the Liquid State [original]

RESEARCH ARTICLE

www.advtheorysimul.com

Predicting and Rationalizing the Soret Coeﬃcient of Binary

Lennard-Jones Mixtures in the Liquid State

Nils E. R. Zimmermann,* Gabriela Guevara-Carrion, Jadran Vrabec,* and Niels Hansen*

The thermodiﬀusion behavior of binary Lennard-Jones mixtures in the liquid

state is investigated by combining the individual strengths of non-equilibrium

molecular dynamics (NEMD) and equilibrium molecular dynamics (EMD)

simulations. On the one hand, boundary-driven NEMD simulations are useful

to quickly predict Soret coeﬃcients because they are easy to set up and

straightforward to analyze. However, careful interpolation is required because

the mean temperature in the measurement region does not exactly reach the

target temperature. On the other hand, EMD simulations attain the target

temperature precisely and yield a multitude of properties that clarify the

microscopic origins of Soret coeﬃcient trends. An analysis of the Soret

coeﬃcient suggests a straightforward dependence on the thermodynamic

properties, whereas its dependence on dynamic properties is far more

complex. Furthermore, a limit of applicability of a popular theoretical model,

which mainly relies on thermodynamic data, was identiﬁed by virtue of an

uncertainty analysis in conjunction with eﬃcient empirical Soret coeﬃcient

predictions, which rely on model parameters instead of simulation output.

Finally, the present study underscores that a combination of predictive

models and EMD and NEMD simulations is a powerful approach to shed light

onto the thermodiﬀusion behavior of liquid mixtures.

1. Introduction

In the nineteenth century, Ludwig[1] and Soret[2] observed in-

dependently from each other that a temperature gradient in

an isotropic mixture can induce a concentration gradient (Fig-

ure 1). This coupling phenomenon between heat and mass ﬂow

N. E. R. Zimmermann, N. Hansen

Institute of Thermodynamics and Thermal Process Engineering

University of Stuttgart

Pfaﬀenwaldring 9, 70569 Stuttgart, Germany

E-mail: [email protected]; [email protected]

G. Guevara-Carrion, J. Vrabec

Thermodynamics and Process Engineering

Technische Universität Berlin

Ernst-Reuter-Platz 1, 10587 Berlin, Germany

E-mail: [email protected]

The ORCID identiﬁcation number(s) for the author(s) of this article

can be found under https://doi.org/10.1002/adts.202200311

Wiley-VCH GmbH. This is an open access article under the terms of the

Creative Commons Attribution License, which permits use, distribution

and reproduction in any medium, provided the original work is properly

cited.

DOI: 10.1002/adts.202200311

is known as thermal diﬀusion, thermodiﬀu-

sion, Ludwig–Soret eﬀect, or simply Soret

eﬀect. Coupled heat and mass transfer oc-

curs in many industrially relevant processes

such as absorption, distillation, extraction,

drying, crystallization, or condensation.[3]

As a second-order eﬀect, the mass ﬂux in-

duced by the Soret eﬀect is usually rel-

atively small, but it can play an impor-

tant role in many natural[4,5] and technical

processes.[6] Thermal diﬀusion can be em-

ployed to separate mixtures that are hard

to handle by conventional methods,[3] for

example, the fractionation of isomers,[7]

isotopes,[8] polymers,[9] and colloids.[10] Fur-

ther applications include surface coating,[11]

geological carbon dioxide storage,[12] mod-

eling of hydrocarbon reservoirs,[13] and hy-

drogen enriched ﬂames.[6]

Chemical systems that undergo ther-

modiﬀusion in practical settings are fre-

quently binary systems in the liquid state.

In such cases, the separation of the two

components can be quantiﬁed by the Soret

coeﬃcient ST, which is a preferred quantity

in the study of liquids.[14] In case of multi-

component mixtures, the Soret coeﬃcient becomes a vectorial

quantity.[15] The Soret coeﬃcient relates the thermal diﬀusion co-

eﬃcient DTto the Fick diﬀusion coeﬃcient D

ST=DT

D(1)

From a technological perspective, the Soret coeﬃcient provides a

compact metric to rate the separation eﬃciency of a process that

is based on thermodiﬀusion (Figure 1).

Two factors render the measurement of the Soret coeﬃ-

cient usually challenging: ﬁrst, the small magnitude of ST

and, second, the presence of undesirable convective currents

due to gravitation.[16,17] Despite eﬀorts to produce benchmark

data,[12,18–21] some even in microgravity environments,[22,23] the

measurement of the Soret coeﬃcient is still challenging, in par-

ticular, for liquid mixtures with more than two components.[16,24]

Consequently, several theoretical models were proposed to pre-

dict the Soret coeﬃcient,[25–27] which were thoroughly reviewed,

for example, in refs. [16, 28]. However, even the most established

models based on non-equilibrium thermodynamics[16] often fail

when dealing with strongly interacting liquid mixtures.

Molecular modeling and simulation oﬀer a complimentary ap-

proach to gain a deeper understanding of the microscopic mech-

25130390, 2022, 11, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/adts.202200311 by Cochrane Germany, Wiley Online Library on [30/11/2022]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

www.advancedsciencenews.com www.advtheorysimul.com

Figure 1. An increasing Soret coeﬃcient STyields stronger separation of

components if the temperature gradient, the Fick diﬀusion coeﬃcient D,

and the mean composition remain constant. The displayed behavior is

referred to as positive thermodiﬀusion because the present component

(typically the heavier one) accumulates in the low temperature region. By

contrast, the component that accumulates in the high temperature region

(typically the lighter one) undergoes negative thermodiﬀusion.

anisms of thermal diﬀusion and also allow for the prediction of

Soret and other transport coeﬃcients. The Lennard-Jones inter-

action potential is particularly attractive because it is a computa-

tionally eﬃcient and yet realistic molecular model, thus helping

to rapidly generate a fundamental understanding of molecular

processes. Furthermore, Lennard-Jones mixtures can exhibit a

rich vapor–liquid equilibrium behavior,[29–31] which can facilitate

constructing new theoretical models and testing application lim-

its of existing ones. Knowing the phase diagram is in itself central

to studying thermodiﬀusion because phase changes are typically

to be avoided, that is, gas bubbles and solid precipitates are un-

wanted when investigating thermal diﬀusion in the liquid state.

Thermodiﬀusion can be simulated in various ways by molec-

ular simulations, where non-equilibrium molecular dynamics

(NEMD) simulations are a well-established standard approach.

These simulations closely mimic experimental scenarios for ther-

modiﬀusion measurements, making them the ideal candidates

to construct benchmark data. NEMD simulations monitor the

response to system perturbations out of equilibrium. Synthetic

NEMD simulations[32–34] add external forces on molecules, thus

altering their dynamics in an artiﬁcial way. On the contrary,

boundary-driven NEMD methods[35–42] aim at paralleling experi-

mental situations much more directly by their conceptual set-up.

These methods impose diﬀerent temperatures within the simu-

lation volume that give rise to temperature gradients, which, in

turn, exert thermodynamic forces on the molecules in an indirect

and natural way. Typically, a narrow region along one spatial di-

rection is assigned as a heat source volume, and a corresponding

heat sink region is placed half the simulation volume apart. The

diﬀerence between the various boundary-driven methods is the

way in which the temperature is controlled in the two regions.

Another approach to simulate thermal diﬀusion and to calcu-

late the Soret coeﬃcient is based on the spontaneous ﬂuctuations

of the microscopic ﬂuxes that are present in equilibrium molec-

ular dynamics (EMD) simulations through the Green–Kubo for-

malism. Although the corresponding ﬂux expressions for mix-

tures were thoroughly derived in the 1980s,[43,44] only a few

molecular simulation studies on thermal diﬀusion employing

EMD techniques can be found in the literature.[32,33,45–50] Further-

more, EMD methods have been mostly used for validation pur-

poses, that is, together with NEMD simulations.[33,45–48,50] Pre-

sumably, this originates from the diﬃculties inherent to EMD

simulations. First, cross-correlation functions are weaker and

thus more challenging to sample than auto-correlation functions.

Second, reliable information on the thermodynamic factor and

the partial molar enthalpy is required to calculate the Soret coef-

ﬁcient. However, EMD simulations are valuable to gain insight

into the sign and magnitude of the individual contributions to

thermal diﬀusion processes, as recognized by Miller et al.[49]

In this work, the strengths of NEMD and EMD simulations

were combined to reliably predict and understand the Soret co-

eﬃcient of three diﬀerent Lennard-Jones mixtures in the liquid

state.[29,30] Furthermore, a popular theoretical model,[27] which

mainly relies on thermodynamic input data, and a regression

model to predict the Soret coeﬃcient were tested, investigating

their applicability limits. Our results and insights underscore

that NEMD and EMD simulations in conjunction with theoret-

ical models build a strong framework to foster understanding of

thermodiﬀusion in liquid mixtures.

2. Methodology

When two or more processes occur simultaneously in a sys-

tem, they may couple, causing new eﬀects.[51] Examples of such

cross-phenomena are the Soret and Dufour eﬀects. The ther-

modiﬀusion or Soret eﬀect addresses the macroscopic motion

of molecules that results in a concentration gradient caused by

a temperature gradient. The reciprocal phenomenon, that is,

a temperature gradient caused by a concentration gradient, is

known as Dufour eﬀect.

2.1. Non-Equilibrium Thermodynamics

Linear non-equilibrium thermodynamics (NET) provides a com-

pact framework to describe the behavior of systems subject

to various driving forces and is fundamentally based on four

postulates:[14]

1) the quasi-equilibrium postulate, requiring that gradients are

not too large;

2) the linearity postulate, stating that all ﬂuxes are linear func-

tions of all relevant driving forces;

3) Curie’s postulate, constraining which ﬂuxes and forces are

coupled; and

4) Onsager’s reciprocal relations, requiring symmetric coeﬃ-

cient matrices in the ﬂux-force relations.

Considering the Soret and Dufour eﬀects only, the equa-

tions for heat ﬂux JQand mass ﬂux of component 1 Jm

1in a binary

mixture can be written in the NET framework as[52]

JQ=−LQQ

∇T

T2−LQ1(𝜕𝜇1

𝜕w1)T,p

∇w1

(1 −w1)T(2)

1=−L1Q

∇T

T2−L11(𝜕𝜇1

𝜕w1)T,p

∇w1

(1 −w1)T(3)

www.advancedsciencenews.com www.advtheorysimul.com

where Tis the temperature, pis the pressure, and 𝜇1and w1

are the chemical potential and mass fraction of component 1,

respectively. Lij are phenomenological Onsager coeﬃcients that

describe the proportionality between the thermodynamic forces

and the ﬂuxes that they induce.[53] The auto-correlated properties

L11 and LQQ are the mass–mass and heat–heat phenomenological

coeﬃcients. The cross-correlated phenomenological coeﬃcients

L1Q=LQ1follow Onsager’s reciprocity relations.[54]

On the other hand, experimentalists usually employ the follow-

ing constitutive equations[46,52] to describe heat and mass ﬂuxes

JQ=−𝜆∇T−𝜌mw1(𝜕𝜇1

𝜕w1)T,p

TDD

T∇w1(4)

1=−𝜌mw1w2DS

T∇T−𝜌D∇w1(5)

where 𝜆is the thermal conductivity, 𝜌mis the mass density of the

mixture, and DS

Tand DD

Trepresent the thermal diﬀusion coeﬃ-

cients of Soret- and Dufour-type, respectively.

When the steady state is reached (i.e., Jm

1=0) and a constant

mass density can be assumed, Equation (5) can be rearranged to

ST=DS

D=− 1

w1w2(∇w1

∇T)Jm

1=0

=− 1

x1x2(∇x1

∇T)J1=0

(6)

where xiis the mole fraction of component iand J1the molar

ﬂux of component 1. Another frequently determined and closely

related quantity is the thermal diﬀusion factor[52] 𝛼T

𝛼T=TST=TDT

D(7)

The relations between phenomenological coeﬃcients and

transport coeﬃcients can be determined by comparing the NET

phenomenological Equations (2) and (3) with the constitutive

Equations (4) and (5). Consequently, the thermal diﬀusion coef-

ﬁcients of Soret- and Dufour-type are thus given by[46]

T=L1Q

𝜌mw1w2T2(8)

T=LQ1

𝜌mw1w2T2(9)

From Onsager’s reciprocity relations, it follows that DD

T=DS

DT.

The Fick diﬀusion coeﬃcient is related to the mass–mass phe-

nomenological coeﬃcient L11 by

D=L11

𝜌mw2T(𝜕𝜇1

𝜕w1)T,p

=L11

𝜌w1w2M1M2

Γ(10)

where 𝜌is the molar density of the mixture. Miis the molar mass

of component i,kBis the Boltzmann constant, and Γthe thermo-

dynamic factor. The Soret coeﬃcient is thus[46]

ST=1

w1T

L1Q

L11 (𝜕w1

𝜕𝜇1)T,p

=w1M2

x1kBT2Γ

L1Q

L11

(11)

Figure 2. a) The key elements of the Langevin NEMD simulations are two

temperature control regions. One of those regions (center) is maintained

at a higher temperature Thigh, whereas the outermost regions that are con-

nected via periodic boundary conditions are kept at a lower temperature

Tlow. The temperature Tand the mole fractions xias functions of the Carte-

sian z-direction are obtained from histograms using 50 bins. b) The pro-

ﬁles of Tand x1obtained from Langevin NEMD simulations are linear.

Two regions (green) can be used to determine slopes for further analysis

to calculate the Soret coeﬃcient.

2.2. Non-Equilibrium Molecular Dynamics

NEMD simulations closely mimic experimental scenarios to

study thermodiﬀusion because these simulations exhibit regions

of high and low temperature within the simulation volume. Fur-

thermore, linear concentration proﬁles and linear temperature

proﬁles in between these regions are obtained at steady state, the

latter implying a constant thermal conductivity.[35] The thermal

conductivity of pure Lennard-Jones ﬂuids and binary mixtures

increases with temperature and density,[35,55–57] the two eﬀects

possibly compensating each other. Thus, a narrow region in the

center of the simulation volume is deﬁned as a heat source and a

corresponding region as a heat sink at the edges of the simulation

volume in z-direction (Figure 2a).

ThetemperatureproﬁleT(z) and the mole fraction proﬁles

xi(z) were determined on the basis of particle trajectories ob-

tained with LAMMPS[58] (stable version: October 29, 2020). Both

temperature-control regions were of equal size, having a width

of 0.3𝜎1, unless stated otherwise, where 𝜎1is the Lennard-Jones

size parameter of component 1. The resulting width was typically

≈1% of the simulation volume. Two Langevin thermostats[59]

www.advancedsciencenews.com www.advtheorysimul.com

were used, where one thermostat heated one temperature-control

region and the other thermostat cooled the other T-control re-

gion, so that we refer to this method as Langevin NEMD.

The Langevin NEMD simulations were typically carried out

with 1200 molecules. The z-direction of the simulation volume

was divided into 50 bins of constant width Δz=lz∕50 (Figure 2a).

Each molecule iwas assigned to one of the bins j,accordingtoits

position zi:j=int(zi∕Δz). For each bin j, the kinetic energy of all

molecules located in the bin Ekin(zj) was determined

Ekin(zj)=1

N(zj)

∑

i=1

mi|vi−v|2(12)

where N(zj) is the number of molecules in bin j,mi,andvithe

mass and the velocity of molecule i, respectively, and vis the

center-of-mass velocity of the system. Since the total momentum

of the system was set to zero, vvanishes. The local temperature

in bin jT(zj) was calculated by[35]

T(zj)=2Ekin(zj)

3N(zj)kB

(13)

The number of degrees of freedom 3 N(zj) should be adjusted

to account for ﬁxing the total momentum of the entire system,

by subtracting three in total, but distributed over the entire tem-

perature proﬁle. Since 50 bins were typically used, we would

need to subtract approximately 3∕50 =0.06 from the number

of degrees of freedom in each bin. As the number of degrees

of freedom was typically (3 ×1200)∕50 =72 in each bin with

N(zj)≈N∕Nbins =1200∕50, the adjustment would have been

on the order of 0.1% and thus negligible. In this context, we

investigated the inﬂuence of system size (i.e., the number of

molecules) on the Soret coeﬃcient—and thus on the tempera-

ture proﬁle—while keeping the number of bins constant, and

found no measurable inﬂuence when using larger systems, in

agreement with previous work.[60] Usually, four diﬀerent temper-

ature diﬀerences ΔTbetween the region of high Thigh and low

temperature Tlow were considered. The relative imposed temper-

ature diﬀerences (Thigh −Tlow)∕(0.5⋅(Thigh +Tlow)) were set to

13.3%, 26.6%, 40%, and 53.3%, respectively. The largest temper-

ature gradient observed in the simulations was 0.026 such that

violation of local equilibrium is not expected.[61,62]

After a certain response time, a steady state was observed,

where the temperature gradient induced a mole fraction gradi-

ent. The simulation time of 3500 to 5000 Lennard-Jones units

spent to establish the steady state was selected based on transient

density proﬁles reported by Bonella et al.[63] The mole fraction

of component kin bin jx

k(zj) was calculated as the ratio of the

number of molecules of component kin bin jN

k(zj) over the total

number of molecules in bin jN(zj)

xk(zj)=Nk(zj)

N(zj)(14)

The slopes of the mole fraction dxk∕dzand temperature proﬁles

dT∕dzas well as the average mole fraction of component k

provide all required input data to calculate the Soret coeﬃcient of

component kusing Equation (6). For each ΔT, ten simulations

Figure 3. a) Parity plot of mean temperature in the proﬁle ﬁt region of the

Langevin NEMD simulations T(z)|Vﬁt and the global mean temperature

⟨T⟩for all systems and state points studied in this work (cf. Table 1). b)

Mean mole fraction of component 1 in the proﬁle ﬁt region of the Langevin

NEMD simulations x1(z)|Vﬁt for all systems and state points.

with independent initial conﬁgurations were performed in order

to average the individual Soret coeﬃcient data and obtain reliable

uncertainty estimates via the sample standard deviation.[64,65]

A transition regime was observed, in which the tempera-

ture and mole fraction proﬁles were still inﬂuenced by the

temperature-control regions (Figure 2b). Therefore, all data

points that were closer than ≈0.05 ×lzto the center of each ther-

mostat layer were discarded. Thus, roughly 20% of the proﬁle

data were discarded in order to maximize the likelihood to ana-

lyze linear proﬁles only.

Unless it is speciﬁed otherwise, the average temperature in

the proﬁle ﬁt region Vﬁt was identiﬁed as the temperature of the

NEMD simulations, that is, T=T(z)|Vﬁt . The average values from

proﬁles agree very well with the corresponding global mean tem-

perature obtained from all molecules in the entire NEMD sim-

ulation volume (Figure 3a). Similarly, the mean mole fraction

of component 1 in the proﬁle ﬁt region of the NEMD simula-

tions x1(z)|Vﬁt was close to the globally imposed value of x1=

0.5molmol

−1(Figure 3b).

www.advancedsciencenews.com www.advtheorysimul.com

Newton’s equations of motion were integrated with a velocity-

Verlet algorithm[66,67] during the production phase. More details

about the NEMD simulations and the corresponding computa-

tional workﬂow are provided in Section S1.1, Supporting Infor-

mation.

2.3. Equilibrium Molecular Dynamics

All EMD simulations were carried out with the ms2 package.[68]

The mass–mass, heat–heat, and heat–mass phenomenolog-

ical coeﬃcients deﬁned by the NET framework[52] can be

sampled with EMD simulations, employing the Green–Kubo

formalism.[69,70] In this framework, the phenomenological mass–

mass phenomenological coeﬃcient L11 is given by[44]

L11 =V

3kB

∞

∫

0⟨Jm

1(0)Jm

1(t)⟩dt(15)

where tis the time. The microscopic mass ﬂux Jm

iis deﬁned as

i(t)=1

∑

mi(vk

i(t)−⟨v⟩) (16)

Here, miis the mass of a molecule of component i,vk

i(t) the cen-

ter of mass velocity vector of molecule kat some time t,Niis the

number of molecules of component i, and the angular brackets

indicate the canonical (NVT) ensemble average. After setting the

total momentum to zero, ⟨v⟩deviated from zero during simula-

tion only within machine error such that Jm

1+Jm

2=0.

The Green–Kubo expression for the heat–mass phenomeno-

logical coeﬃcient LQ1is

LQ1=V

3kB

∞

∫

0⟨JQ(t)Jm

1(0)⟩dt(17)

In EMD simulations, the heat ﬂux JQcannot be directly

determined—in contrast to the internal energy ﬂux JEwhich can

be directly computed. Thus, the associated internal energy–mass

phenomenological coeﬃcient LE1was calculated by

LE1=V

3kB

∞

∫

0⟨JE(t)Jm

1(0)⟩dt(18)

where the internal energy ﬂux JEis given by

JE=1

∑

i=1

∑

k=1

i(vk

i−⟨v⟩)2

⋅(vk

i−⟨v⟩)

−1

∑

i=1

∑

j=1

∑

k=1

∑

l≠k[rkl

ij :

𝜕ukl

𝜕rkl

−I⋅ukl

ij ]⋅(vk

i−⟨v⟩) (19)

Here, ukl

ij is the inter-molecular potential energy and rkl

ij is the dis-

tance vector between molecules kand l. The indices iand jde-

note the components of the mixture. The dyadic product is rep-

resented by : and Iis the unitary tensor. In a binary mixture, LQ1

can be determined from LE1if the partial molar enthalpies hiof

both components are known[71]

L1Q=L1E−L11(h1

−h2

M2)(20)

Because of Onsager’s reciprocity LE1=L1E. Thus, both cross-

correlated coeﬃcients were independently sampled during sim-

ulation and averaged to improve statistics.

2.3.1. Partial Molar Enthalpy

In this work, the partial molar enthalpy was calculated from EMD

simulations in two steps. First, the residual molar enthalpy of the

binary mixture hres was calculated for ten diﬀerent mole fractions

in the isobaric–isothermal (NpT) ensemble over the composition

range x1=0.5±0.05 mol mol−1. Subsequently, an appropriate

function h=f(x1) was ﬁtted by a least squares optimization to the

simulation data to generate an expression for the molar enthalpy

h=hres +hid.

The partial molar enthalpy hiwas then obtained from

hi=h+xj(𝜕h

𝜕xi)T,p

(21)

2.3.2. Thermodynamic Factor

To determine the Fick diﬀusion coeﬃcient from the mass–

mass phenomenological Onsager coeﬃcients according to Equa-

tion (10), the thermodynamic factor Γis required

Γ= x1

kBT(𝜕𝜇1

𝜕x1)T,p

(22)

It describes the thermodynamic non-ideality of a mixture and can

be obtained from Kirkwood–Buﬀ integrals[72] Gij

Gij =4𝜋

∞

∫

0(gij(r)−1)r2dr(23)

employing

Γ=1−𝜌1𝜌2(G11 +G22 −2G12)

𝜌1+𝜌2+𝜌1𝜌2(G11 +G22 −2G12)(24)

where gij (r) is the radial distribution function, while 𝜌1and

𝜌2are the partial densities 𝜌i=xi𝜌. The truncation method by

Krüger et al.[73] and the corrections of the radial distribution func-

tion by Ganguly and van der Vegt[74] were employed because of

convergence issues.[75]

www.advancedsciencenews.com www.advtheorysimul.com

Table 1. Lennard-Jones parameters and state points.

System, Mass ratio Size ratio Energy ratio Temperature Pressure Cutoﬀ radius Time step size

state point m2∕m1𝜎2∕𝜎1𝜀2∕𝜀1kBT∕𝜀1p𝜎3

1∕𝜀1rcut∕𝜎1Δt∕103

A1 0.5 1.5 0.75 0.75 0.09 3.75 0.7071

A2 0.5 1.5 0.75 1 0.09 3.75 0.7071

A3 0.5 1.5 0.75 1 1.8 3.75 0.7071

B1 3.375 1.5 1 0.75 0.09 3.75 1

B2 3.375 1.5 1 1 0.09 3.75 1

B3 3.375 1.5 1 1 1.8 3.75 1

C1 1 1 0.75 0.75 0.09 2.5 1

C2 1 1 0.75 1 0.09 2.5 1

C3 1 1 0.75 1 1.8 2.5 1

2.4. Models

Equimolar (i.e., x1=x2=0.5molmol

−1) binary liquid mixtures

of molecules that interact via the Lennard-Jones potential were

considered

ukl

ij =4⋅𝜀ij[(𝜎ij

rkl

ij )12

−(𝜎ij

rkl

ij )6](25)

where 𝜀ij and 𝜎ij are the Lennard-Jones energy and size parameter

for the interaction between molecule pairs belonging to compo-

nents iand j,andrkl

ij denotes the distance between molecules k

and l. The interaction parameters if i=j𝜎iand 𝜀iare provided

in Table 1. The parameters for the interaction between diﬀer-

ent components were determined with the Lorentz–Berthelot

combining rules[76] unless stated otherwise: 𝜎ij =0.5(𝜎i+𝜎j)and

𝜀ij =√𝜀i𝜀j. Furthermore, the potential was explicitly evaluated

up to a cutoﬀ radius rcut =2.5×max(𝜎1,𝜎2),[39] and analytic

tail corrections were applied to pressure and energy terms.[77,78]

Small time step sizes Δtwere chosen to integrate Newton’s equa-

tions of motion[79] (Table 1).

All results are presented in Lennard-Jones units, where the

mass m1, the Lennard-Jones energy 𝜀1, and Lennard-Jones size

parameters of component 1 𝜎1were used for reducing all quanti-

ties: m∗=m∕m1,E∗=E∕𝜀1,andl∗=l∕𝜎1,whereEdenotes en-

ergy and llength or distance. The reduced time and temperature

are t∗=t√𝜀1∕(m1𝜎2

1)andT∗=TkB∕𝜀1, respectively. For the sake

of brevity, the asterisks on the reduced quantities are omitted in

the remainder of this work.

3. Results

First, we present NEMD results because the corresponding Soret

coeﬃcients serve as a benchmark for the subsequent EMD sim-

ulations. For both approaches, speciﬁc issues, pitfalls, and nec-

essary procedures are highlighted to obtain reliable Soret coeﬃ-

cient data. This is followed by a comparison and discussion of the

data for all Lennard-Jones models considered in this work, where

two lean prediction models are also included. The ﬁrst model

is based on thermodynamic and approximate kinetic data,[27]

Figure 4. Soret coeﬃcient of krypton STin its mixture with argon as a func-

tion of temperature Tusing the present Langevin NEMD method (circles).

The color from blue to orange indicates an increasing relative temperature

diﬀerence and thus a rising Tgradient in the simulation volume. Litera-

ture values[32,33,35,41,49] (gray diamonds) and their average value[41] (black

square) are provided for comparison.

whereas the second one is a regression model that relies on

Lennard-Jones and mass parameters of the force ﬁeld.[39]

3.1. NEMD: Creating a Benchmark

3.1.1. Validation

The Langevin NEMD method and the corresponding analysis

workﬂow was validated on the basis of literature data for the

mixture argon/krypton.[32,33,35,41,49] The present Soret coeﬃcient

data (Figure 4) agree well with data from literature, which were

obtained using a broad variety of approaches: NEMD with syn-

thetic NEMD simulations,[32,33] the heat exchange method,[35] as

well as the modiﬁed heat exchange method and reverse NEMD

simulations[41] and Green–Kubo-based EMD simulations.[49] The

average temperature that was observed in the central proﬁle-ﬁt

region of our simulation volume varied with the imposed tem-

perature diﬀerence. The colors on the depicted data points cor-

respond to the observed relative temperature diﬀerence ΔT∕T.It

can be concluded that the present workﬂow is reliable but that it

www.advancedsciencenews.com www.advtheorysimul.com

Figure 5. Soret coeﬃcient of component 1 STfor diﬀerent temperatures T

and target pressures pas obtained from Langevin NEMD simulations for

the Lennard-Jones mixture A studied in this work (cf. Table 1). The color

from blue to orange indicates an increasing relative temperature diﬀerence

ΔT∕Tin the Langevin NEMD simulations. Open symbols represent linear

interpolations of STand its standard deviation around the desired target

temperature Ttarget. Corresponding plots for the other two Lennard-Jones

mixtures are provided in the Supporting Information (Figure S5).

leaves some imprecision regarding the target temperature, which

is on the order of 1%.

The Langevin NEMD procedure was also validated based on se-

lected Lennard-Jones systems,[39,71] and the agreement was usu-

ally excellent (Figure S2, Supporting Information). Furthermore,

the Langevin NEMD method and the enhanced heat exchange

(eHEX) algorithm[40] as implemented in LAMMPS[58] yield con-

sistent results (Figure S4, Supporting Information). Therefore,

it can be concluded that the Langevin NEMD procedure reliably

yields Soret coeﬃcient data for Lennard-Jones systems.

3.1.2. Interpolation

The Soret coeﬃcient from a given Langevin NEMD simulation

is associated with the mean temperature in the ﬁt region, which

is the common practice when handling the rather complicated

temperature-dependency of thermodiﬀusion quantities.[14] Typ-

ically, Soret coeﬃcient data from NEMD simulations should

be extrapolated to the zero-force regime,[32,39,80] where a linear

function can be used.[39] However, MacGowan and Evans esti-

mated zero-force transport coeﬃcients by weighted averaging if

a distinct force dependence of a given coeﬃcient could not be

detected.[32] On the one hand, the variation of the Soret coeﬃ-

cient is typically small for diﬀerent thermodynamic driving forces

around a given state point in this work (Figure 5). On the other

hand, the average value of the temperature in the gradient ﬁt re-

gion does not exactly match the target temperature (i.e., mean

value of imposed high and imposed low temperature), but the

deviations are small: a few percent. Furthermore, the Soret coef-

ﬁcient data determined with diﬀerent gradients and at diﬀerent

temperatures usually exhibit consistent temperature trends (Fig-

ure 5). Therefore, comparison to EMD predictions, in which the

target temperature is reached exactly, was facilitated by interpo-

lation. Speciﬁcally, STand its standard deviation were linearly in-

terpolated in the small window around the target temperature to

provide the most likely values at the target temperature.

3.1.3. Screening Other Inﬂuences

Validation on the basis of literature data and solving the interpola-

tion issue were central to verifying the reliability of our Langevin

NEMD method. However, other technical and model-related fac-

tors that can possibly inﬂuence the Soret coeﬃcient were also

sounded out. Speciﬁcally, it was reassured that the Langevin

damping constant, the width of the temperature control regions,

the system size, and the cutoﬀ radius have no systematic impact

on ST(Section S2.4, Supporting Information), which leads us to

the conclusion that the Langevin NEMD method proposed here

yields reliable and robust results.

3.2. EMD: Obtaining Causal Data

The EMD approach based on correlation functions and the

Green–Kubo formalism is challenging for dense liquids because

the cross-correlations required to calculate thermal diﬀusion are

at least one order of magnitude smaller than the auto-correlation

functions[44] required to determine the self-diﬀusion coeﬃcient.

In a ﬁrst step, EMD simulations were performed using the same

simulation parameters as for the NEMD simulations (Table 1).

The resulting Soret coeﬃcient data are provided in Figure 8 and

Table S1, Supporting Information. In a second step, the simula-

tion parameters were optimized to yield consistent data with the

NEMD simulations.

3.2.1. Finite Size Eﬀects and Cutoﬀ Radius

It is well known that small molecular systems under periodic

boundary conditions are associated with systematic errors when

diﬀusion coeﬃcients are calculated. In this context, the size de-

pendence of the mass–mass L11 and internal energy–mass LE1

phenomenological coeﬃcients was studied for one of the con-

sidered systems (i.e., mixture B2). A series of simulations with

varying system size containing 600, 1200, 2400, 4800, and 8000

molecules was performed. As expected, a clear system size de-

pendence was observed for the mass–mass phenomenological

coeﬃcient L11 (Figure 6a). Its value in the thermodynamic limit,

obtained by extrapolation to inﬁnite size, was found to be 5.5%

higher than its sampled value for 1200 molecules. Therefore, L11

had to be corrected to account for ﬁnite size eﬀects when working

with small systems, for example, as proposed in ref. [81]. A ﬁnite

size correction of L11 based on a modiﬁed version of the correc-

tion factor by Yeh and Hummer [82,83] for binary mixtures yields

an improvement, but L11 remains underestimated by 4%. On the

other hand, no clear system size dependence could be inferred

for the internal energy–mass coeﬃcient LE1despite its diﬀusive

nature (Figure 6b).

The prediction of the thermodynamic factor by Kirkwood–Buﬀ

integration is also subject to system size eﬀects, as can be seen in

Figure 7. Because this formalism was originally derived for the

grand canonical (𝜇VT) ensemble, it has to be extrapolated to the

www.advancedsciencenews.com www.advtheorysimul.com

Figure 6. a) Mass–mass and b) internal energy–mass phenomenological

coeﬃcients of mixture B2 as functions of the inverse edge length of the

simulation volume l−1. The dashed lines are a linear ﬁt and a ﬁt of a con-

stant (average) to the simulation results, respectively.

thermodynamic limit when the NVT ensemble is used. The cor-

rection method by Krüger et al.[73] was successfully employed in a

previous work[75] to determine extrapolated values when dealing

with small systems. However, when suﬃciently large systems are

employed, this methodology can lead to an over-correction of the

thermodynamic factor (Figure 7).

Because of ﬁnite size-related issues, additional EMD sim-

ulations were performed with larger volumes for all studied

mixtures. For this purpose, relatively large systems with 8000

molecules were employed and the cutoﬀ radius was extended to

half of the edge length of the simulation volume l∕2. For mixture

A2, the diﬀerence between the extrapolated and sampled mass

transfer coeﬃcient was within its statistical uncertainty (i.e.,

≈2%). Further, the sampled thermodynamic factor diﬀered

by +2.5%from its extrapolated value. Therefore, the sampled

values of the mass–mass phenomenological coeﬃcient and the

thermodynamic factor were employed directly without further

ﬁnite size corrections for the calculation of the Soret coeﬃcient.

On the other hand, the coeﬃcients sampled in a smaller sim-

ulation volume, containing 1200 molecules, do require ﬁnite

size corrections.

Figure 8 shows the Soret coeﬃcient of the studied Lennard-

Jones mixtures for diﬀerent system sizes (i.e., 1200 and 8000

molecules). As can be seen, there is a very good agreement be-

Figure 7. a) Thermodynamic factor of the mixture B2 as a function of

the inverse edge length of the simulation volume l−1, including a linear

ﬁt (dashed line) and the resulting extrapolated value in the thermody-

namic limit (ochre cross). b) Thermodynamic factor sampled for systems

with 1200 (solid circles) and 8000 (solid squares) molecules, respectively,

and their corrected values (empty symbols) following Krüger et al.[73] as

functions of the calculated Soret coeﬃcient. The ochre cross represents

the linearly extrapolated thermodynamic factor of mixture B2, as obtained

from (a).

Figure 8. Soret coeﬃcient calculated with 1200 (solid circles) and 8000

(solid squares) molecules. The open symbols represent coeﬃcients cal-

culated with the corrected mass–mass phenomenological coeﬃcients[81]

L11 and the corrected thermodynamic factor[73] Γ. The ochre symbols rep-

resent data calculated with 1200 molecules and a cutoﬀ radius of l∕2.

www.advancedsciencenews.com www.advtheorysimul.com

Figure 9. Range of the calculated Soret coeﬃcient when the magnitude of

the mass–mass phenomenological coeﬃcient (ochre) changes by ±10%

and when the partial molar enthalpy diﬀerence (h1∕M1−h2∕M2)(blue)

changes by ±10%. The original Soret coeﬃcient data with uncertainties

(solid squares) were added for comparison.

tween the Soret coeﬃcient based on the simulations with 1200

molecules and those obtained from the simulations with 8000

molecules for the mixtures B and C. The peculiarly strong diﬀer-

ences found for mixture A could be traced back to the use of a rel-

atively short cutoﬀ in the simulations with 1200 molecules. When

the cutoﬀ radius was changed to l∕2, the agreement with the 8000

molecules simulation was signiﬁcantly improved. A strong de-

pendence of the cross-correlated phenomenological coeﬃcient

LE1on the cutoﬀ radius was found especially for mixture A1.

Changes of other transport properties due to the cutoﬀ radius

(e.g., diﬀusion coeﬃcients, shear viscosity, and thermal conduc-

tivity) were found to be within their statistical uncertainties. Neg-

ligible changes were also found for the thermodynamic equilib-

rium properties. Solely, the average pressure was 15% lower for

the simulations with a smaller cutoﬀ radius.

3.2.2. Sensitivity

Two factors were identiﬁed to highly aﬀect the Soret coeﬃcient

value when predicted from EMD simulations: the mass–mass

phenomenological coeﬃcient and the diﬀerence between the par-

tial molar enthalpies of both components (h1∕M1−h2∕M2). For

the majority of considered mixtures, a change of a few percent of

either of these quantities leads to a large change of the Soret coef-

ﬁcient. For example, Figure 9 shows the resulting changes of the

Soret coeﬃcient when the mass–mass phenomenological coeﬃ-

cient L11 and the enthalpy diﬀerence (h1∕M1−h2∕M2) are varied

by ±10%. The strong variation of the Soret coeﬃcient is related to

LQ1. As can be seen in Figure 10, the contribution of the enthalpy

on the heat–mass phenomenological coeﬃcient is large for most

of the studied mixtures. Furthermore, in cases where terms of

Equation (20) have a very similar magnitude a small change of

the second term may even lead to a change of sign of the Soret co-

eﬃcient.

Because of the strong sensitivity of the EMD simulations for

the studied mixtures, the results obtained from the larger sys-

tems featuring 8000 molecules are regarded as the more ade-

Figure 10. Soret coeﬃcient of component 1 as a function of the heat–

mass phenomenological coeﬃcient LQ1and the internal energy–mass

phenomenological coeﬃcient LE1, respectively, obtained from EMD simu-

lations with 8000 molecules.

quate values for comparison with NEMD results. The numerical

results are listed in Table S2, Supporting Information.

3.2.3. Other Properties

EMD simulations for the calculation of the Soret coeﬃcient are

challenging and have to be performed with carefully chosen pa-

rameters. However, all relevant quantities for gaining a micro-

scopic understanding of the heat and mass transfer processes can

be calculated directly with EMD: the Fick and the thermal diﬀu-

sion coeﬃcients as well as the thermal conductivity or the degree

of coupling q=LQ1∕(L11LQQ )0.5. These quantities are listed in Ta-

ble S3, Supporting Information.

3.3. NEMD versus EMD: Rationalizing ST

Figure 11 shows the ﬁnal data set with Soret coeﬃcients of

Lennard-Jones component 1 STfrom NEMD and EMD simula-

tions as functions of temperature for two diﬀerent average pres-

sures and three Lennard-Jones mixtures (A, B and C; cf. Table 1).

None of the systems studied gives rise to a sign change in the

Soret coeﬃcient and thus to a switch in thermophoretic behav-

ior when the temperature or pressure is altered. The temperature

and pressure ranges investigated are small and the data sparse.

The observation is yet important because the change from posi-

tive to negative thermodiﬀusion behavior is subject of many stud-

ies due to the technological signiﬁcance of this phenomenon.

The pressure impact on STis consistent across the diﬀerent

systems; however, the temperature inﬂuence can be diﬀerent. In

all cases, an increase in pressure from 0.09 to 1.8 yields a decrease

of the Soret coeﬃcient of component 1: by 30%, 417%, and 53%

for mixtures A, B, and C, respectively. By contrast, STincreases

with rising temperature from 0.75 to 1 at p=0.09 for mixture B

by 50%, whereas it decreases by 44% and 6% for mixtures A and

C, respectively. However, half of these relative change ﬁgures

dST=100 ×ST(T2or p2)−ST(T1or p1)

|ST(T1or p1)|%(26)

www.advancedsciencenews.com www.advtheorysimul.com

Figure 11. Soret coeﬃcient of component 1 STfor diﬀerent temperatures

Tand pressures p(squares: p=0.09; circles: p=1.8) of three diﬀerent

Lennard-Jones mixtures A (top), B (center), and C (bottom) as obtained

from NEMD (black symbols) and EMD simulations (blue symbols). The-

oretical predictions according to the model by Shukla and Firoozabadi[27]

(small purple symbols) are also included as well as predictions on the ba-

sis of Reith’s and Müller-Plathe’s regression model (ochre solid line).[39]

exhibit large uncertainties when relying on a 2𝜎dSTthreshold

to achieve 95% conﬁdence. For the temperature-induced ST

changes, the uncertainties are 11%, 67%, and 37% and, for the

pressure-induced STchanges, they are 15%, 622%, and 23% for

mixtures A, B, and C, respectively. We used the STuncertain-

ties 𝜎STfrom the NEMD simulations to obtain the uncertain-

ties 𝜎dSTand applied conventional uncertainty propagation rules,

which yields, for example, for the temperature-induced relative

STchange

𝜎dST=100 ×[(𝜎ST(T2)

|ST(T1)|)2

+(𝜎ST(T1)×[−|ST(T1)|−(±1)(ST(T2)−ST(T1))]

[ST(T1)]2)2]0.5

%(27)

where the factor ±1 depends on whether ST(T1)≷0.

NEMD and EMD simulations usually yield consistent results

because the 1𝜎STuncertainties always overlap—with exception of

mixture B, state point 2. The statistical uncertainties of STfrom

EMD simulations are larger than those from NEMD simulations,

which can be explained by the relatively small signal-to-noise ra-

tio inherent to the cross-correlation functions. For mixtures A, all

1𝜎STuncertainties are small, whereas those of mixtures B and C

might seem moderate or large. This, however, is only a plotting

eﬀect because the vertical axis range in Figure 11 is more narrow

for mixture B and much more narrow for mixture C.

3.4. Simulations versus Theoretical Models: Eﬃciently Predicting

Several theoretical models[25,27,84,85] for predicting STleverage

readily accessible thermodynamic information. In particular, the

model by Shukla and Firoozabadi[27] is successful for gaseous and

liquid mixtures.[86] The Soret coeﬃcient is given by[27,86]

SSF

T=− u1∕𝜏1−u2∕𝜏2

x1(𝜕𝜇1∕𝜕x1)T,p T−(v2−v1)(x1u1∕𝜏1+x2u2∕𝜏2)

(x1v1+x2v2)x1(𝜕𝜇1∕𝜕x1)T,p T(28)

where uiand viare the partial molar internal energy and the par-

tial molar volume of component i, respectively. The component-

speciﬁc factor 𝜏i, which is the cohesive energy of the liquid di-

vided by the activation energy supplied for molecular motion,[27]

ranges typically between three and ﬁve for liquids.[27,86] The pre-

dictions are included in Figure 11 using 𝜏i=3.5 (cf. ref. [86] and

Section S2.6, Supporting Information). Since none of the state

points chosen for the three considered Lennard-Jones mixtures

is in the proximity of the critical point, the model should work

well.[27]

We performed a detailed uncertainty propagation analysis[87]

for the Soret coeﬃcient given by the Shukla–Firoozabadi model

(Section S2.6, Supporting Information). This could explain some

of the discrepancies observed between NEMD, EMD, and pre-

dictive model results. For mixtures A and C, the model predicts

the correct sign and most of the temperature and pressure de-

pendences of ST. The diﬀerence between NEMD simulation and

model prediction relative to the model uncertainty is, on an av-

erage, however as large as 285% for these two mixtures. While

uncertainties cannot be made fully responsible, they could con-

tribute signiﬁcantly to the deviations. In order to better under-

stand the precise source of the model uncertainty, its relative con-

tributions were calculated (Section S2.6 and Table S4, Supporting

Information). The by-far largest contribution to the uncertainty

in SSF

Tof component 1 stems from the 𝜏1parameter 𝜎𝜏1.Thenext

largest contribution originates usually from 𝜎𝜏2. The critical role

of the 𝜏parameters is also reﬂected in the inability of the model

by Shukla and Firoozabadi to predict the correct sign of STfor

highly asymmetric mixture B. Here, the regression model by Re-

www.advancedsciencenews.com www.advtheorysimul.com

ith and Müller-Plathe[39] might help shed light onto the reason

for the prediction failure.

The Reith–Müller-Plathe model consists of three additive

quadratic functions, each of which was separately ﬁtted to a se-

ries of NEMD-derived Soret coeﬃcients. In each series, a single

model parameter of component 1, m1,𝜎1,or𝜀1, was systemati-

cally varied so that the impact of one model parameter was en-

coded in one quadratic function. The resulting overall Soret coef-

ﬁcient, which was successfully validated on realistic systems such

as isotopic argon–argon, argon–krypton, and xenon–methane, is

given by[39]

SRMP

T=−0.7(m1

m2)2

+9.5(m1

m2)−8.8

+67.4(𝜎1

𝜎2)2

−179.3(𝜎1

𝜎2)+111.9

+4.4(𝜀1

𝜀2)2

+3.5(𝜀1

𝜀2)−7.9 (29)

The applicability ranges for the diﬀerent contributions are lim-

ited by originally ﬁtting within m1∕m2≤8, 𝜎1∕𝜎2≤1.25, and

𝜀1∕𝜀2≤1.75. Moreover, the Soret coeﬃcient is not given in

Lennard-Jones units but in (10−3K−1). Hence, SRMP

Thas to be ad-

justed to the present unit system.

The regression model by Reith and Müller-Plathe can reliably

predict the sign of the Soret coeﬃcient. Furthermore, the magni-

tude of SRMP

Tis reasonable, considering the fact that this predic-

tion does not incorporate any temperature or pressure/density

dependence. Another drawback of this model is the assumption

that there are no coupling eﬀects among the Lennard-Jones mass,

size, and energy parameters. However, it has been found that

such eﬀects play a non-negligible role for the prediction of ther-

modiﬀusion of equimolar Lennard-Jones mixtures.[88]

The insight that was gained from these predictions is that

mass and Lennard-Jones size parameter have opposing eﬀects

on the Soret coeﬃcient. The lower mass of component 1 alone

would result in negative thermodiﬀusion, whereas the smaller

size parameter would trigger positive thermodiﬀusion. It thus

depends on which eﬀect is stronger and hence which type of

model parameter ratio is larger. In the case of mixture B, the

mass ratio ultimately dictates the negative sign. The 𝜏parameters

used in the model by Shukla and Firoozabadi[27] were originally

determined for real systems, for which size and mass typically

correlate. It is therefore possible that these parameters require

a re-optimization if size and mass strongly diverge, as seems

to hold for mixture B. Very recently, Hoang and Galliero [89]

came to a similar conclusion on the lack of reliability of the

Shukla–Firoozabadi model for the prediction of thermal diﬀu-

sion factors for Lennard-Jones mixtures that are asymmetric in

size and mass.

Finally, note that another possibility to describe the eﬀects

of mass, size, and energy parameter ratios on transport prop-

erties of Lennard-Jones mixtures is the van der Waals one-ﬂuid

model.[76,90] This model considers the mixture to be a hypothet-

ical pseudo-compound that approximates its equilibrium and

transport properties. [91] However, this approach has signiﬁ-

cant limitations when dealing with strongly asymmetric mix-

tures [56,91] so that is not applied to the present systems.

4. Conclusion

We used Langevin NEMD alongside EMD simulations to in-

vestigate the Soret coeﬃcient of binary Lennard-Jones mixtures.

The Langevin NEMD simulations served as benchmarks because

EMD predictions based on the Green–Kubo formalism are chal-

lenging due to a weak signal-to-noise ratio of the cross-correlation

functions, the strong system-size dependence of the mass–mass

phenomenological coeﬃcient, and the uncertainties related to

the calculation of the thermodynamic factor. Once optimal set-

tings for the EMD simulations were determined, the data-set

from the EMD simulations was successfully used to rational-

ize the behavior of the Soret coeﬃcient. Furthermore, the EMD

data were leveraged to test the applicability of a lean theoretical

model,[27] which allows for predictions of the Soret coeﬃcient

mainly based on thermodynamic data.

We could identify a limitation of the theoretical model by

Shukla and Firoozabadi.[27] Importantly, a detailed uncertainty

analysis, which has, to the best of our knowledge, not been per-

formed for the model yet, and complementing insights from a re-

gression model,[39] which leverages model parameters of the mix-

ture chosen, helped us pinpointing the probable reason. The the-

oretical model is unable to account for opposing eﬀects of mass

and size on the Soret coeﬃcient. A parameter that is mainly ki-

netic in nature plays here a major role in dictating the sign of

the Soret coeﬃcient and thus in the applicability of the analyti-

cal model. More elaborate models such as those by Artola et al.[92]

could be used instead, but this model requires, for example, more

transport data.

Extending the combination of NEMD and EMD simulations to

more complex systems is expected to provide insight into prac-

tical applications of thermodiﬀusion. The NEMD method was

successfully applied to study aqueous alkali halide solutions[93]

as well as to ternary[94,95] and quaternary[94] liquid mixtures of or-

ganic molecules. For both methods, the sampling of species with

low concentration poses a challenge in these applications.

Supporting Information

Supporting Information is available from the Wiley Online Library or from

the author.

Acknowledgements

The authors N.E.R.Z. and N.H. acknowledge funding by the Deutsche

Forschungsgemeinschaft (DFG, German Research Foundation) under

Germany’s Excellence Strategy - EXC 2075 – 390740016. The authors

appreciate the support by the Stuttgart Center for Simulation Science

(SimTech) and by the High Performance and Cloud Computing Group

at the Zentrum für Datenverarbeitung of the University of Tübingen, the

state of Baden-Württemberg through bwHPC, and the DFG through grant

no INST 37/935-1 FUGG. The authors G.G.-C. and J.V. acknowledge the

ﬁnancial support of the DFG under the grant VR 6/11 and want to give

special thanks to Jean-Marc Simon for his valuable help on the implemen-

tation and validation EMD Soret calculations. EMD simulations were per-

formed on the HPE Apollo system Hawk at the High Performance Com-

puting Centre Stuttgart (HLRS) contributing to the project MMHBF2. The