This version is available at https://doi.org/10.14279/depositonce-9918
Copyright applies. A non-exclusive, non-transferable and limited
right to use is granted. This document is intended solely for
personal, non-commercial use.
Terms of Use
Janzen, T., Fingerhut, R., Heinen, M., Köster, A., Muñoz-Muñoz, Y. M., & Vrabec, J. (2019). Molecular
Modeling and Simulation: Force Field Development, Evaporation Processes and Thermophysical Properties
of Mixtures. In High Performance Computing in Science and Engineering ’ 18 (pp. 457–474). Springer
International Publishing. https://doi.org/10.1007/978-3-030-13325-2_29
Tatjana Janzen, Robin Fingerhut, Matthias Heinen, Andreas Köster, Y.
Mauricio Muñoz-Muñoz, Jadran Vrabec
Molecular Modeling and Simulation: Force
Field Development, Evaporation Processes
and Thermophysical Properties of Mixtures
Accepted manuscript (Postprint) Conference paper |

Molecular Modeling and Simulation: F or ce
Field De velopment, Evaporation Pr ocesses and
Thermoph ysical Properties of Mixtur es
T atjana Janzen, Robin Fingerhut, Matthias Heinen, Andreas Köster , Y . Mauricio
Muñoz-Muñoz, and Jadran Vrabec
Abstract T o gain physical insight into the beha vior of fluids on a microscopic le vel
as well as to broaden the data base for thermophysical properties especially for mix-
tures, molecular modeling and simulation is utilized in this work. V arious methods
and applications are discussed, including a procedure for the de v elopment of ne w
force field models. The e vaporation of liquid nitrogen into a supercritical h ydrogen
atmosphere is presented as an example for lar ge scale molecular dynamics simu-
lation. System size dependence and scaling beha vior are discussed in the conte xt
of Kirkwood-Buf f integration. Further , results for thermophysical mixture proper -
ties are presented, i.e. the Henry’ s law constant of aqueous systems and dif fusion
coef ficients of a ternary mixture.
1 Intr oduction
Molecular modeling and simulation provides a broad range of possibilities for the
in vestig ation of v arious thermophysical properties. On the one hand, it is a po werful
tool for the prediction of thermodynamic data of pure fluids as well as mixtures,
which are needed for process design in chemical and energy technology engineer -
ing. Especially when dealing with toxic substances or extreme conditions, challeng-
ing and expensi ve e xperiments can be av oided with this approach. On the other
hand, molecular simulation not only contrib utes to the a v ailability of thermody-
namic data, b ut also allo ws for a precise insight into the fluid behavior on the micro-
scopic le vel that gov erns macroscopic behavior . This promotes the understanding
of dif ferent phenomena and yields a connection between microscopic structure and
thermodynamic properties.
T . Janzen · R. Fingerhut · M. Heinen · A. Köster · Y . M. Muñoz-Muñoz · J. Vrabec
Thermodynamics and Thermal Separation Processes, T echnical Uni v ersity Berlin, Ernst-Reuter-
Platz 1, 10587 Berlin, Germany , e-mail: [email protected]
1

2 T . Janzen et al.
T o exploit the possibilities of molecular modeling and simulation, suitable com-
putational methods combined with high-performance computing are essential. W ith-
in this project, two simulation packages are constantly impro ved and extended:
ms 2 [1, 2, 3], for equilibrium simulations of b ulk fluids, and ls1 mar dyn [4], for
the simulation of non-equilibrium systems and flo w phenomena on the nanoscale.
Se veral ex emplary application cases of these two software packages are presented
in this work, including model de velopment and molecular simulations of v arying
system sizes and dif ferent thermodynamic properties of mixtures.
2 F orce Fields f or Cyclohexylamine, Aniline and Phenol
The reliable representation of thermodynamic data by molecular simulation is de-
termined by the quality of the underlying force field, i.e. the mathematical represen-
tation of molecular configurational energy U . Numerous studies ha v e been carried
out on this matter , e.g. using sophisticated numerical methods, such as machine
learning [5, 6] or by combining ab initio and classical molecular simulations [7].
The computational ef fort of classical molecular simulations for systems com-
posed of small molecules may be reduced by neglecting internal molecular de grees
of freedom and assuming pairwise additi ve interactions, i.e. taking the configura-
tional ener gy as U = ∑ N
i = 1 ∑ N − 1
j > i U i j = U el ec + U d is p , wherein N is the number of
particles of the system and U i j represents the interaction ener gy between particles i
and j . It depends on the relati ve position of the molecular interaction sites belonging
to molecules i and j and their electrostatic interaction parameters.
According to this approach, a force field for classical molecular simulations pos-
sesses three types of parameters. The electrostatic and geometric parameters as well
as the dispersion ener gy parameters. Molecular geometry , i.e. the location of each
atom in the molecular frame work, and the electrostatic parameters can be inherited
from ab initio calculations, e.g. using the Møller -Plesset 2 method. The use of ab
initio simulations is the backbone for de v eloping force fields and the initial starting
point of our de velopments as can be seen in Fig. 1.
On the other hand, the intermolecular dispersion ener gy U d is p can be described
by the Lennard-Jones (LJ) potential, i.e. U d is p = U d is p ( ε , σ ) , where ε and σ are
fitting parameters, often located at the position of each atom. Computational ef fort
can be reduced e ven more by using the united-atom approach, which consists of
grouping se veral atoms and assigning them a common center of force, represented
by a single set of LJ parameters ε and σ . Indeed, the selection of the parameters
to be fitted is a matter of technical art. Once the LJ parameters to be adjusted are
selected, the optimization can be carried out follo wing the workflo w presented in
Fig. 1. Molecular Monte-Carlo (MC) or molecular dynamics (MD) simulations are
carried out for calculating the density of the pure fluid. Subsequently , v apor-liquid
equilibrium (VLE) simulations are carried out with the Grand Equilibrium method,
which requires the chemical potential v alues.

Molecular Modeling and Simulation 3
Ab i nit i o simula ti on :
i nit i al f or c e fiel d
MC or MD simula ti on :
es ti ma ti on of sa tur a t ed
l i qu i d de ns i ty
Is the es ti ma t ed l i qu i d
den sit y sa ti s f ac t ory ?
Es ti ma t e the v apor -
l i qu i d eq uil i bri a (V LE)
Ar e the VLE simula ti on
da t a sa ti s f act or y ?
R ep ort the f or c e fiel d
par amet er s
y es
y es
No
V ar y the elect r os t a ti c
par amet er s
V ar y the dispe r sion
en er gy par amet er s
No

Fig. 1 Flo wchart for de veloping multi-site LJ + point charges force fields.
The workflo w gi ven in Fig. 1 w as employed for de veloping three ne w multi LJ
site + point char ge force fields for cyclohe xylamine (C 6 H 13 N), aniline (C 6 H 7 N)
and phenol (C 6 H 6 O). Cyclohexylamine, being technically used e.g. to delay the
degradation of rubber or to pre vent corrosion in steam plants, w as represented by
se ven LJ sites and four point char ges. T wo of those point char ges are positi ve and
located at hydrogen atom positions, sho wn in white on top of Fig. 2. One negati ve
char ge is located at the center of the nitrogen atom, sho wn in blue on top of Fig. 2
and the last positi ve char ge is located at the center of the united atom methanetriyl LJ
site, depicted in orange. Sites in yello w on top of Fig. 2 stand for the fiv e identical
methanediyl LJ united-atom groups of c yclohexylamine, whose parameters were
selected to be fitted.
Aniline, important e.g for the pharmaceutic industry , was represented by se v en
LJ sites and four point char ges. Similarly to c yclohexylamine, three point charges
are located at the amine group. T wo positi ve char ges are at the center of the hydro-
gen atoms, sho wn in white on top of Fig. 2. One neg ati ve point char ge is located
at the center of the nitrogen LJ site, sho wn in blue on top of Fig. 2 and the last
positi ve char ge is located at the center of the carbon LJ site, depicted in black. The
sites in orange represent the fi ve methine united-atom LJ sites forming the aromatic
ring of the aniline molecule. The aromatic ring of phenol, which is important to
yield phenol-formaldehyde resins, is similarly represented as the aniline one. The
hydroxyl group, bonded to the aromatic ring, is represented by three point char ges.
One positi ve char ge is located at the center of the hydrogen atom, one ne gati ve
char ge is located at the center of the oxygen LJ site, which is sho wn in red on top of
Fig. 2. The LJ parameters of the methine sites in both aromatic rings were selected
to be fitted to experimental VLE data.

4 T . Janzen et al.
Fig. 2 V apor -liquid equilibrium results for cyclohe xylamine (red symbols), aniline (blue symbols)
and phenol (green symbols). Panel (a) represents the saturated liquid and v apor density results, pan-
els (b) and (c) stand for the enthalpy of v aporization and vapor pressure, respecti vely , and panel (d)
sho ws the compressibility factor of the saturated v apor phase. Lines represent DIPPR correlations,
crosses are the a v ailable experimental literature data and squares are present molecular simulation
results.
Fig. 2 gi ves an ov erview on the present molecular simulation results for these
three important fluids. These data are in good agreement with e xperimental litera-
ture v alues and the DIPPR correlations. The absolute a verage de viations for phenol
with respect the DIPPR correlations are 3.5 % for v apor pressure, 0.3 % for saturated
liquid density and 5.6 % for enthalpy of v aporization, considering the temperature
range where experimental data are a v ailable. The cyclohe xylamine force field ex-
hibits absolute a v erage de viations of 4.7 % for v apor pressure, 0.1 % for saturated
liquid density and 3.3 % for enthalpy of v aporization. The absolute av erage devia-
tions for aniline are 6.7 % for v apor pressure, 0.4 % for saturated liquid density and
7.2 % for enthalpy of v aporization. Critical data were calculated as well and they are
sho wn in Fig. 2. The critical temperature of phenol was slightly ov erestimated with
respect to literature v alues, while critical pressure and density were underestimated.
The de viations are 1.0 %, -8.6 % and -15.3 %, respecti v ely . Relativ e deviations of
critical temperature, pressure and density with respect to literature data for cyclo-
hexylamine are -1.1 %, -7.1 % and 27 %, respecti vely . Analogously , for aniline these

Molecular Modeling and Simulation 5
numbers are -0.7 %, -0.7 % and -5.9 %. Furthermore, present molecular simulation
results are consistent with a criterion recently proposed by Nezbeda [8] as can be
seen in Panel (c) of Fig. 2.
3 Stationary Evaporation
The e vaporation of liquid nitrogen into a supercritical h ydrogen atmosphere was in-
vestig ated by lar ge scale molecular dynamics (MD) simulation. This work w as mo-
ti vated by the results of Dahms et al. [9], who studied that injection for tw o different
operating pressures: 1 and 6 MPa. In case of the lo wer pressure, their liquid nitrogen
jet disintegrated into small droplets forming a spray , whereas in case of the higher
pressure, dif fuse mixing occurred. Dahms et al. used nitrogen as a safe substitute for
oxygen in their experiments, which were aimed at the optimization of comb ustion
processes in rocket engines. It is well kno wn that the interface between liquid and
v apor phase plays a key role in e v aporation processes. Fundamental in vestig ations
of e vaporation into v acuum based on model liquids sho w that the magnitude of the
e vaporation flux is strongly dependent on the interf ace temperature. W ith decreas-
ing interface temperature, the e v aporation flux drops do wn significantly [10]. The
interface properties, especially the temperature, are thus decisi ve for e v aporation
dynamics. Consequentially , it was focused in the present work on the interfacial re-
gion between liquid nitrogen and the high pressure hydrogen atmosphere. Because
of its sound physical basis and atomistic resolution, MD simulation is an appropriate
approach to obtain a better understanding of the findings by Dahms et al. [9].
Fig. 3 Schematic vie w of system construction by the ReplicaGenerator module. The liquid ni-
trogen N 2 phase and the supercritical hydrogen H 2 phase were prepared from small fluid cubes
(indicated by grid cells). These small systems were replicated to yield two lar ger homogeneous
phases. These were equilibrated for an additional short time period and then merged together to
yield the initial configuration where the two phases are in physical contact for non-equilibrium
molecular dynamics simulations (NEMD).

6 T . Janzen et al.
Molecular systems including an ov erall number of ≈ 6 · 10 6 molecules were set
up with a planar interface, with an e xtent of A xy = 85 x 85 nm 2 , between a liquid
nitrogen film and a supercritical hydrogen atmosphere. A schematic vie w of how
the system was constructed is depicted in Fig. 3. A ne w module of the MD code
ls1 mar dyn [4], namely the ReplicaGenerator , w as used to create comparati vely
lar ge homogeneous phases of pure nitrogen N 2 and pure hydrogen H 2 out of small
ones. First, the small systems were equilibrated to obtain physical molecular con-
figurations. Subsequently , they were replicated and merged together , yielding larger
phases. These were equilibrated again for a short time until the replication pattern
v anished. This procedure reduces the equilibration time which is particularly impor-
tant for lar ge phases.
The system was delimited on both left and right sides ( z direction) by reflecti ve
boundaries, spanning a system width of L z = 190 nm. These boundaries were e x-
tended by additional molecules with fixed spatial positions to reduce their influence.
T o mimic an infinitely extended interf ace, periodic boundary conditions were spec-
ified in x and y directions. The temperature of the liquid nitrogen of T N 2 = 118 K
was controlled within a control v olume on the left side, which had a width of about
ten molecule diameters only . The temperature of hydrogen was controlled within
a control v olume on the right side to be T H 2 = 270 K, also maintaining a constant
density and composition. The latter was achie ved by con verting all incoming nitro-
gen molecules to hydrogen molecules such that a mole fraction of y H 2 = 1 (pure
hydrogen) was maintained. T o establish a stationary e vaporation process, a liquid
reserv oir added ne w nitrogen molecules ov er the left boundary of the control v ol-
ume into the liquid phase. The feed rate of the reserv oir was adjusted to balance the
rate of e vaporated nitrogen molecules entering the control v olume in the gas phase.
There were no interventions between the control v olumes, except for the numeri-
cal solution of Ne wton’ s equations of motion of the molecular ensemble with the
MD code ls1 mar dyn [4]. The molecular models for nitrogen and hydrogen were
provided by Köster et al. [11], sho wing a very good agreement with experimental
literature data and the Peng-Robinson equation of state for v apor-liquid equilibria of
pure substances and mixtures, considering a lar ge temperature and pressure range.
A series of three simulations with v arying system pressure, controlled by the target
density of the hydrogen gas in the right control v olume, was carried out: p = 4, 12
and 20 MPa. During simulation, v arious profiles were sampled. Each profile was
a v eraged ov er 10 4 time steps, with a time step length of ∆ t = 2 fs. At the beginning
of the simulations, all profiles were step functions because the systems were set up
by mer ging an equilibrated liquid nitrogen phase at a temperature of T N 2 with an
equilibrated hydrogen gas phase at a temperature of T H 2 . W ithin the first nanosec-
ond, these profiles changed rapidly , dev eloping into a smooth transition between the
liquid phase and gas phase states maintained by the control v olumes. After about
3 ns, the profiles hardly changed any more such that a stationary process had been
attained. Fig. 4 sho ws the temperature, hydrogen mole fraction and density profiles
of the simulation series for two pressure le vels p = 4 and 20 MPa after an elapsed
time of 0, 1 and 3 ns. The case of p = 12 MPa sho wed qualitati vely similar re-
sults as p = 20 MPa and the results are therefore not illustrated here. Looking at

Molecular Modeling and Simulation 7
the right column of Fig. 4, se veral aspects can be seen that may contrib ute to the
understanding of the macroscopic system studied by e xperiment [9]. The density
profile of the lo w pressure case ( p = 4 MPa) has a comparati vely sharp interface,
transitioning from the liquid density do wn to the gas density . The density profiles of
the high pressure case ( p = 20 MPa) e xhibit a slightly higher liquid density and a
significantly higher gas density due to compression. Moreo v er , a much broader and
smoother transition between the liquid and gas density than in the lo w pressure case
was found.
Fig. 4 T emporal e volution of the temperature T (top), mole fraction of hydrogen x H 2 (center) and
molar density (bottom) profiles of the lo w pressure case p = 4 MPa (black) and the high pressure
case p = 20 MPa (red). Columns left, center , right show profiles after 0 , 1 and 3 ns.
The temperature profile of the lo w pressure case sho ws a minimum within the
interface re gion, indicating that a part of the enthalpy of e v aporation is provided by
the liquid phase that consequently cools do wn. From a molecular perspectiv e, the
liquid film cools do wn at the interface because molecules with higher kinetic en-
er gy preferentially lea ve the interf ace, whereas molecules with lower kinetic ener gy
preferentially remain inside the liquid. The interface temperature decrease is only
3.5 K because heat is transported from both phases to the interface via temperature
gradients. The temperature profiles of the higher pressure cases do not exhibit a
temperature minimum at the interface. This can be e xplained by the higher pressure
(and density) of the gas phase, leading to a much stronger thermal coupling with the
liquid. Consequentially , a sufficient amount of heat is transported from the hot gas
phase to balance the amount of energy tak en up by ev aporating nitrogen molecules.
Instead, the liquid is e v en heated up, indicating that more ener gy is transported from

8 T . Janzen et al.
the gas to the liquid phase than required to dri ve e vaporation. The mole fraction pro-
files sho w expected trends: The gas phase is rich in hydrogen and the liquid phase is
rich in nitrogen. Consistent with the density profiles, the lo w pressure case exhibits a
sharp transition, whereas the higher pressure cases ha v e a smooth and much broader
transition, indicating dif fuse mixing.
4 Kirkwood-Buff Integration
Local density fluctuations within closed systems can be sampled with Kirkwood-
Buf f integration (KBI) [14] and lead to macroscopic thermodynamic properties.
Many-body molecular simulation naturally pro vides detailed information on the mi-
croscopic structure of fluids by means of radial distrib ution functions (RDF) g i j ( r ) ,
which are closely related to what is kno wn as local composition. In case of a mixture
containing components i and j , they are gi ven by
g i j ( r ) = ρ L , i ( r )
ρ j
, (1)
where the local density of component i in a spherical shell centered around species
j is
ρ L , i ( r ) = d N i
4 π r 2 d r , (2)
and d N i denotes the number of particles of component i therein. The overall partial
density is simply ρ j = x j ρ . A binary mixture of components 1 and 2 has thus three
independent radial distrib ution functions, i.e. g 11 ( r ) , g 22 ( r ) and g 12 ( r ) = g 21 ( r ) .
Moreov er , radial distribution functions pro vide access to the composition depen-
dence of the excess Gibbs energy by means of KBI. A good con ver gence of their
integrals
G i j = 4 π Z ∞
0 [ g i j ( r ) − 1 ] r 2 d r , (3)
can be achie ved by canonic NV T ensemble simulations employing the formalism
of Krüger et al. [16]. The Kirkwood-Buf f integrals G i j lead to se veral macroscopic
properties, e.g. the mole fraction deri vati ve of the acti vity coef ficient γ 1 of compo-
nent 1 is gi ven by
 ∂ ln γ 1
∂ x 1  T , p
= − x 2 ρ ( G 11 + G 22 − 2 G 12 )
1 + x 2 ρ x 1 ( G 11 + G 22 − 2 G 12 ) . (4)
For studies on mass transport, the thermodynamic factor Γ is vital because it con-
nects the Fick dif fusion coefficients D i j and Maxwell-Stefan dif fusion coef ficients
Ð i j (for a binary mixture: D 11 = Γ 11 Ð 12 ) and it is gi ven by

Molecular Modeling and Simulation 9
Γ 11 = 1 + x 1  ∂ ln γ 1
∂ x 1  T , p
. (5)
In the field of process engineering, mixture data are of particularly great interest.
Properties like partial molar v olumes, isothermal compressibility and most inter-
estingly thermodynamic factor can be determined with KBI. T o this end, KBI w as
implemented into our massi vely-parallel molecular simulation tool ms 2 [3]. Ho w-
e ver , finite system-size ef fects are challenging for KBI calculations by using molec-
ular simulations of closed systems because KBI is defined for macroscopically large
systems only . Further , for large distances r , the sampled RDF tend to lose statistical
accuracy and may not con ver ge to unity . Accordingly , RDF corrections are cru-
cial [15] in the same way as an e xtrapolation to macroscopic size [17]. Moreo ver , it
was sho wn by Galata et al. [12] that simulations with a lar ger number of particles
N may lead to improv ed results. Thus, parallelization of KBI in ms 2 is essential
because the computational ef fort is directly dependent on the number of particles N .
Finite system-size ef fects of KBI were assessed for a binary Lennard-Jones (LJ)
mixture for two dif ferent system sizes in the closed N V T ensemble with N = 4 , 000
and N = 8 , 000 particles (LJ parameters σ 2 / σ 1 = 1 . 5 and ε 2 / ε 1 = 0 . 75). The LJ mix-
ture was studied o ver the entire molar fraction of x 1 = 0 . 05 , 0 . 1 , 0 . 2 ,..., 0 . 95 mol/mol.
The mixture density range was sampled with isobaric-isothermal N pT simula-
tions at constant pressure p σ 3
1 / ε 1 = 0 . 03 and constant temperature T k B / ε 1 = 0 . 85.
Molecular dynamics (MD) NV T simulations for both system sizes were carried out
with 800,000 equilibration steps and 15 Mio. production steps with a time step of
0.329 fs and a cutof f radius of r c = 5 σ 1 . RDF were sampled for each time step with
a maximum radius di vided into 300 shells. Fig. 5 shows the thermodynamic f actor
Γ 11 ov er mole fraction x 1 for the two system sizes N = 4 , 000 and N = 8 , 000.
It becomes clear that the simulation results for the thermodynamic factor Γ 11 are
slightly improv ed for larger system sizes because of the decreased de viation to the
W ilson model [19] fit. On that point, the thermodynamic factor Γ 11 calculated by the
W ilson model may serve for comparison since it is independent from the LJ model
error .
The parallelization performance of ms 2 was analyzed for the MD NV T ensemble
with KBI calculation turned on and of f. Therefore, three different system sizes of the
abov e-mentioned LJ mixture with N = 4 , 000; 8 , 000; 16 , 000 particles were carried
out on CRA Y XC40 (Hazel Hen) of HLRS and the runtime of ms 2 was measured
by changing the number of nodes (each node contained 24 cores). Fig. 6 sho ws the
scaling of MD NV T simulations of ms 2 for those three systems with KBI calculation
turned on and of f.
For all three systems, the cutof f radius was set to r c = 5 σ 1 and thus the NV T
interaction calculation alone required the same computing po wer , which only de-
pends on the number of particles N in the cutof f sphere. Howe ver , to determine
which particles are inside the cutof f sphere, the whole interaction matrix has to be
tra v ersed which allo ws for a computationally effecti ve access to the calculation of
RDF shells if KBI is turned on. The computing intensity of tra versing the particle
matrix is proportional to N 2 . As a result, the computational cost of ms 2 should be

10 T . Janzen et al.
x 1 / mol mol -1
0.0 0.2 0.4 0.6 0.8 1.0
Γ
0.4
0.6
0.8
1.0
1.2
1.4

Fig. 5 Thermodynamic factor Γ 11 ov er mole fraction x 1 of a binary LJ mixture calculated
with MD and KBI in the NV T ensemble employing ms 2; black symbols: N = 4 , 000; red
symbols (slighlty shifted to larger x 1 for visibility reasons): N = 8 , 000; circles/triangles: non-
extrapolated/e xtrapolated data to macroscopic system size; solid line: W ilson model [19] fitted to
chemical potential data calculated for the LJ mixture with W idom’ s test particle insertion [18];
dashed line: PC-SAFT equation of state [13].
nodes (with 24 cores)
250 500 750 100 0
time * nodes / (10 6 molecules * steps) / s
0
2
4
6
8

Fig. 6 Strong scaling ef ficiency of ms 2 for a binary LJ mixture in the MD NV T ensemble mea-
sured on CRA Y XC40 (Hazel Hen); solid line: MD N V T (KBI turned of f); dashed line: MD NV T
(KBI turned on); blue lines: N = 4 , 000; black lines: N = 8 , 000; red lines: N = 16 , 000 particles.

Molecular Modeling and Simulation 11
in between N to N 2 , e.g. doubling N should lead to a factor 2 to 4 in terms of com-
putational cost. Fig. 6 indicates that doubling the number of particles in ms 2 leads
to an increase of computational cost smaller than factor 4 if the number of nodes is
chosen appropriately so that the ov erhead is small. Further , the vertical axis gi ves
the computing po wer (nodes) times computing time per computing intensity (prob-
lem size) and thus, a horizontal line would sho w a perfect strong scaling ef ficienc y
of 100 %. Consequently , the Fig. 6 points out that ms 2 is close to optimal strong
scaling, b ut it of course does not reach 100 % ef ficiency . Moreover , it becomes clear
that for small systems with computationally cheap interaction calculations like LJ
mixtures, the strong scaling ef ficiency decreases with increased computing po wer
(nodes) because of communication ov erhead.
5 Henry’ s law constant
In an extensi ve study that w as recently published in full [11], the mixture behavior
of fi ve substances related to the hydrogen economy , i.e. H 2 , nitrogen (N 2 ), oxygen
(O 2 ), ar gon (Ar) and water (H 2 O), w as e xamined using these techniques. Se veral
thermodynamic properties were considered in this work to describe a substantial
number of mixtures, including binary , ternary and quaternary systems. Apart from
VLE properties, like v apor pressure, saturated densities or enthalpy of v aporization,
also homogeneous density and solubility were studied. The quality of the results
was assessed by comparing to a vailable e xperimental literature data or sophisticated
equations of state, sho wing an excellent agreement in many cases. Consequently ,
this contrib ution aims at an improv ed av ailability of thermodynamic data that are
required for the post-carbon age.
The present work gi ves a short o vervie w on the Henry’ s law constant that can be
described particularly well with molecular simulation techniques and was used in
the work of Köster and Vrabec [11] to assess aqueous systems, i.e. H 2 , N 2 , O 2 or
Ar in liquid H 2 O. Consequently , this property is used in cases where a solute is only
little soluble in a solvent. It has to be noted that the purely temperature-dependent
definition of the Henry’ s law constant w as considered, wherein the solvent has to
be in its saturated liquid state. Therefore, data are presented ranging from the triple
point temperature to the critical temperature of H 2 O.
The Henry’ s law constant H i was sampled with the molecular simulation tool
ms 2 [1, 2, 3] in the N pT ensemble using 1000 particles in total. The simulation
runs typically consisted of 3 · 10 6 MC cycles, i.e. 2 · 10 5 equilibration and 2 . 8 · 10 6
production cycles. H i is closely related to the chemical potential at infinite dilution
µ ∞
i [20]
H i = ρ s k B T exp ( µ ∞
i / ( k B T )) , (6)
wherein ρ s is the saturated liquid density of the solv ent, T the temperature and k B
Boltzmann’ s constant. An attracti v e method to calculate the chemical potential at
infinite dilution is W idom’ s test particle insertion [18], where the mole fraction of

12 T . Janzen et al.
the solute can simply be set to zero such that only test particles are inserted into the
pure solvent, i.e. H 2 O.












 
 
 
 
   
  


 

Fig. 7 Henry’ s law constant of H 2 , N 2 , O 2 and Ar in H 2 O (from top to bottom): ( • ) Molecular
simulation results, (—) of ficial IAPWS guideline [22, 21] and ( + ) experimental literature data (H 2
[23, 24]; N 2 [25, 26, 24, 27]; O 2 [28, 29, 30, 31]; Ar [32, 33, 34, 35]).
All aqueous systems presented in Fig. 7 were v alidated against e xperimental data
and the of ficial IAPWS correlation [21, 22] for the Henry’ s law constant. The y ex-
hibit a strongly non-monotonic beha vior that is qualitati vely similar for the v arious
solutes and pass through a maximum between 330 and 370 K. Due to the fact that
higher v alues of the Henry’ s law constant correspond to a lo w solubility , the solu-
bility decreases with increasing temperature, passes through a minimum and then

Molecular Modeling and Simulation 13
increases again. On the quantitati ve side, ho we ver , large dif ferences between the
solutes are present. Near the triple point temperature of H 2 O, both H 2 and N 2 are
roughly three times less soluble than O 2 and Ar ( H H2 ≈ H N2 ≈ 6 GP a compared to
H O2 ≈ H Ar ≈ 2 GP a). Moreover , the Henry’ s law constant v alues at the maximum
dif fer significantly .
It can be seen that Henry’ s law constant data were predicted satisf actorily by
molecular simulation, especially at temperatures close to the triple and critical point
of H 2 O. At intermediate temperatures, molecular simulation results for N 2 , O 2 and
Ar de viate slightly , whereas H 2 is reproduced almost perfectly . It has to be noted,
ho wev er , that the force field of H 2 was modelled without electrostatic interactions.
This was compensated by increasing the unlik e LJ interaction ener gy , resulting in a
rather unphysically lar ge v alue of the binary interaction parameter ξ = 1 . 52.
6 Diffusion Coefficients of T er nary Liquid Mixtures
Numerous natural and technical processes are determined by transport dif fusion of
multicomponent liquids, while the underlying microscopic phenomena are still not
well understood. Experimental data on multicomponent dif fusion coefficients are
rather scarce because strong composition dependence and cross dif fusion ef fects
aggra v ate experimental work. Molecular modeling and simulation is a promising
route for in vestig ating dif fusion coefficients on a microscopic basis to understand
the underlying phenomena and to study the relations between intradif fusion and
transport dif fusion coefficients. Nonetheless, not only the dif fusion experiments,
b ut also the simulations are challenging and time consuming. T o gain an o verall
picture of the dif fusion behavior of a ternary mixtures, a suf ficient amount of state
points with v arying composition must be studied, also considering the binary sub-
systems. Thus, suf ficient parallelization of the utilized software as well as ef ficient
data handling are essential to ov ercome this task.
Here, equilibrium molecular dynamics (MD) simulations and the Green-K ubo
formalism were used to sample intradif fusion and Maxwell-Stefan (MS) dif fusion
coef ficients of binary and ternary liquid mixtures. Intradiffusion coef ficients D i ,
which describe the mobility of molecules of species i in the absence of dri ving force
gradients, were obtained from the time integral of the single particle velocity auto-
correlation function a v eraged ov er all N i particles belonging to that species [36]
D i = 1
3 N i Z ∞
0 D N i
∑
k = 1
v k
i ( 0 ) · v k
i ( t ) E d t . (7)
MS dif fusion coefficients, which describe mass transport dri ven by a chemical
potential gradient, are directly related to kinetic coef ficients which can be sampled
from net velocity correlation functions of species i and j [37]

14 T . Janzen et al.
Λ i j = 1
3 N Z ∞
0 D N i
∑
k = 1
v k
i ( 0 ) ·
N j
∑
l = 1
v l
j ( t ) E d t . (8)
Fick dif fusion coefficients are related to gradients of composition v ariables and
can thus be measured experimentally . For the transformation of MS to Fick dif fusion
coef ficients the thermodynamic factor is needed [42], which can be calculated with
excess Gibbs ener gy ( g E ) models that are typically fitted to experimental v apor-
liquid equilibrium data. Details on the relations between the dif ferent coefficients
for ternary mixtures can be found in recent publications of our group [38, 39].
One of the ternary mixtures studied within this project was acetone + benzene +
methanol. This mixture was in vestigated at ambient conditions of temperature and
pressure, for which experimental data of Fick dif fusion coefficients are a vailable in
the literature [41]. Dif fusion and other transport properties of the binary subsystems
were already predicted by MD simulations in the present MMHBF project [40].
The utilized molecular force field models are rigid, united-atom Lenard-Jones type
models with superimposed point char ges, dipoles and quadrupoles. These models
were parametrized using experimental data of the pure components only , thus the
mixture beha vior was obtained in a strictly predicti ve manner . MD simulations were
carried out with 4000 molecules in a cubic v olume o ver 4 to 5 × 10 7 time steps of
∼ 1.02 fs. Further details on the simulation setup used to sample ternary dif fusion
coef ficients can be found in the supplementary material of Ref. [39].
Fig. 8 sho ws the intradiffusion coef ficients D i of the three components ov er the
entire composition range of the ternary mixture acetone + benzene + methanol.
These figures are based on simulations for 52 ternary and 40 binary state points
at dif ferent compositions indicated by points within the graphs. The coef ficients of
all three components ha ve a similar composition dependence, with lo west v alues at
small acetone content and highest v alues at high acetone content. Because methanol
molecules form hydrogen bonds and often propagate as assemblies, the intradif fu-
sion coef ficient is similar to that of benzene. Only at low methanol concentrations
this intradif fusion coefficient is significantly higher due to the breaking of h ydrogen
bonds.
Fig. 9 sho ws the three MS diffusion coef ficients Ð i j of the ternary mixture ace-
tone + benzene + methanol. While the v alues of the MS coefficient between acetone
and benzene Ð 12 are similar to the intradif fusion coefficients, the other tw o MS co-
ef ficients Ð 13 and Ð 23 show a significantly non-ideal beha vior . This is also a result
of the hydrogen bonding beha vior of methanol molecules, which leads to micro-
scopic heterogeneity in mixtures with other fluids and af fects their kinetics through
correlated propagation. Fick dif fusion coefficients predicted from the simulated MS
dif fusion coefficients combined with the thermodynamic f actor calculated with the
W ilson g E model are in excellent agreement with e xperimental data [39].

Molecular Modeling and Simulation 15
Fig. 8 Intradif fusion coef ficient D i (in 10 − 10 m 2 s − 1 ) of acetone (a), benzene (b) and methanol (c)
in their ternary mixture at T = 298 . 15 K and p = 0 . 1 MPa.
7 Conclusion
Se veral applications of molecular modeling and simulation were outlined in this
work, considering dif ferent simulation techniques, system sizes and thermophysical
properties. The workflo w for the de velopment and optimization of force field mod-
els, which are the basis of e very molecular simulation, was presented. Three ne w
force field models for cyclohe xamine, aniline and phenol were parametrized and
v alidated against experimental VLE data.
As an application for lar ge scale MD simulations, the e vaporation of liquid nitro-
gen into a supercritical hydrogen atmosphere was presented. Here, special attention
was paid to the beha vior of the interface between the two phases. T emperature,
density and concentration profiles were sampled to study the dependencies of the
e vaporation flux.
Further , dif ferent thermophysical properties of mixtures were predicted. Kirk-
wood-Buf f integrals, which are related to the radial distrib ution functions and de-
scribe local density fluctuations, were used to sample the thermodynamic factor of
binary LJ mixtures. Because the results sho w a significant system size dependence, a
good scaling beha vior of the simulation software is necessary , which is sho wn to be
gi ven for the utilized simulation package ms 2. The temperature dependent Henry’ s

16 T . Janzen et al.
Fig. 9 Maxwell-Stefan dif fusion coefficients Ð 12 (a), Ð 13 (b) and Ð 23 (c) (in 10 − 10 m 2 s − 1 ) of
acetone (1) + benzene (2) + methanol (3) at T = 298 . 15 K and p = 0 . 1 MPa.
law constant of aqueous systems containing H 2 , N 2 , O 2 or Ar w as calculated from
chemical potential data sampled with W idom’ s test particle insertion. The results are
in excellent agreement with e xperimental data. Finally , results for the intradif fusion
and Maxwell-Stefan dif fusion coef ficients were presented for the ternary mixture
acetone + benzene + methanol in the liquid state co vering the entire composition
range of the gi ven mixture.
Acknowledgements W e gratefully ackno wledge support by Deutsche Forschungsgemeinschaft.
This work was carried out under the auspices of the Boltzmann-Zuse Society (BZS) of Computa-
tional Molecular Engineering. The simulations were performed on the CRA Y XC40 (Hazel Hen)
at the High Performance Computing Center Stuttgart (HLRS) within the project MMHBF2.
Refer ences
1. Deublein, S., Eckl, B., Stoll, J., Lishchuk, S. V ., Gue vara-Carrion, G., Glass, C. W .,
Merker , T ., Bernreuther , M., Hasse, H., Vrabec, J.: ms2: A molecular simulation tool for
thermodynamic properties. Comput. Phys. Commun. 182 , 2350–2367 (2011)

Molecular Modeling and Simulation 17
2. Glass, C. W ., Reiser , S., Rutkai, G., Deublein, S., Köster , A., Guev ara-Carrion, G., W afai, A.,
Horsch, M., Bernreuther , M., W indmann, T ., Hasse, H., Vrabec, J.: ms2: A molecular simu-
lation tool for thermodynamic properties, ne w version release. Comp. Ph ys. Commun. 185 ,
3302–3306 (2014)
3. Rutkai, G., Köster , A., Guev ara-Carrion, G., Janzen, T ., Schappals, M., Glass, C. W ., Bern-
reuther , M., W afai, A., Stephan, S., K ohns, M., Reiser , S., Deublein, S., Horsch, M., Hasse, H.,
Vrabec, J.: ms2: A molecular simulation tool for thermodynamic properties, release 3.0.
Comp. Phys. Commun. 221 , 343–351 (2017)
4. Niethammer , C., Becker , S., Bernreuther , M., Buchholz, M., Eckhardt, W ., Heinecke, A.,
W erth, S., Bungartz, H.-J., Glass, C. W ., Hasse, H., Vrabec, J., Horsch, M.: ls1 mardyn: The
massi vely parallel molecular dynamics code for lar ge systems. J. Chem. Theory Comput. 10 ,
4455–4464 (2014)
5. Chmiela, S., Tkatchenko, A., Sauceda, H. E., Polta vsky , I., Schütt, K. T ., Müller , K.-R.: Ma-
chine learning of accurate energy-conserving molecular force fields. Science Adv ances 3 , 1–6
(2017)
6. Li, Y ., Li, H., Pickard, F . C., Narayanan, B., Sen, F . G., Chan, M. K. Y ., Sankaranarayanan,
S. K. R. S., Brooks, B. R., Roux B.: Machine learning force field parameters from ab initio
data. J. Chem. Theory Comput. 13 , 4492–4503 (2017)
7. Ho, Y .-C., W ang, Y .-S., Chao, S. D.: Molecular dynamics simulations of fluid cyclopropane
with MP2/CBS-fitted intermolecular interaction potentials. J. Chem. Phys. 147 , 064507,
(2017)
8. Nezbeda, I.: Simulations of v aporliquid equilibria: Routine versus thoroughness. J. Chem.
Eng. Data 61 , 3964–3969 (2016)
9. Dahms, R., and Oefelein, J. C.: On the transition between two-phase and single-phase inter -
face dynamics in multicomponent fluids at supercritical pressures. Phys. Fluids 25 , 092103
(2013)
10. Heinen, M., Vrabec, J., and Fischer , J., J.: Communication: Evaporation: Influence of heat
transport in the liquid on the interface temperature and the particle flux. Chem. Phys. 145 ,
081101 (2016)
11. Köster , A., Thol, M., Vrabec, J.: Molecular Models for the Hydrogen Age: Hydrogen, Nitro-
gen, Oxygen, Argon, and W ater . J. Chem. Eng. Data 63 , 305–320 (2018)
12. Galata, A. A., Anogiannakis, S. D., Theodorou, D. N.: Thermodynamic analysis of Lennard-
Jones binary mixtures using Kirkwood-Buf f Theory . Fluid Phase Equilib . 470 , 25–37 (2018)
13. Gross, J., Sado wski, G.: Perturbed-chain SAFT : An equation of state based on a perturbation
theory for chain molecules. Ind Eng Chem Res. 40 , 1244–1260 (2001)
14. Kirkwood, J. G., Buf f, F .P .: The Statistical Mechanical Theory of Solutions. I. J. Chem. Phys.
19 , 774–778 (1951)
15. Milzetti, J., Nayar , D., v an der V egt, N. F . A.: Con vergence of Kirkw ood–Buff Inte grals
of Ideal and Nonideal Aqueous Solutions Using Molecular Dynamics Simulations. J. Phys.
Chem. B 122(21) , 5515–5526 (2018)
16. KrÃijger , P .,Schnell, S. K.,Bedeaux, D., Kjelstrup, S., Vlugt, T . J. H., Simon, J. M.:
Kirkwood-Buf f integrals for finite v olumes. J. Phys. Chem. Lett. 4 , 10911–10918 (2013)
235â ˘
A ¸ S238.
17. Schnell, S. K., Liu, X., Simon, J.-M., Bardo w , A., Bedeaux, D., Vlugt, T . J. H., Kjelstrup, S.:
Calculating Thermodynamic Properties from Fluctuations at Small Scales. J. Phys. Chem. B
115 , 10911–10918 (2011)
18. W idom, B.: Some T opics in the Theory of Fluids. J. Chem. Phys. 39 , 2808–2812 (1963)
19. W ilson, G. M.: V apor-liquid equilibrium. XI. A ne w expression for the excess free ener gy of
mixing. J. Am. Chem. Soc. 86 , 127 (1964)
20. Shing, K.S., Gubbins, K. E., Lucas K.: Henry constants in non-ideal fluid mixtures: computer
simulation and theory . Mol. Phys. 65 , 1235–1252 (1988)
21. Fernández-Prini, R., Alv arez, J. L., Harve y , A. H.: Henry’ s constants and vapor –liquid distri-
b ution constants for gaseous solutes in H 2 O and D 2 O at high temperatures. J. Phys. Chem.
Ref. Data 32 , 903–916 (2003)

18 T . Janzen et al.
22. International Association for the Properties of W ater and Steam. Guideline on the Henry’ s
Constant and V apor-Liquid Distrib ution Constant for Gases in H 2 O and D 2 O at High T em-
peratures (2004)
23. Morris, D. R., Y ang, L., Giraudeau, F ., Sun, X., Stew ard, F . R.: Henry’ s law constant for
hydrogen in natural water and deuterium in hea vy water . Phys. Chem. Chem. Phys. 3 , 1043–
1046 (2001)
24. Alv arez, J., Fernández-Prini, R.: A semiempirical procedure to describe the thermodynamics
of dissolution of non-polar gases in water . Fluid Phase Equilib . 66 , 309–326 (1991)
25. Rettich, T . R., Battino, R., W ilhelm, E.: Solubility of gases in liquids. XVI. Henry’ s law
coef ficients for nitrogen in water at 5 to 50 ◦ C. J. Solution Chem. 13 , 335–348 (1984)
26. Cosgrov e, B. A., W alkley , J.: Solubilities of gases in H 2 O and 2 H 2 O. J. Chromatogr . A 216 ,
161–167 (1981)
27. Alv arez, J., Crov etto, R., Fernández-Prini, R.: The dissolution of N 2 and of H 2 in water from
room temperature to 640 K. Ber . Bunsenges. Phys. Chem. 92 , 935–940 (1988)
28. Rettich, T . R., Battino, R., W ilhelm, E.: Solubility of gases in liquids. 22. High-precision
determination of Henry’ s law constants of oxygen in liquid water from T= 274 K to T= 328
K. J. Chem. Thermodyn. 32 , 1145–1156 (2000)
29. Benson, B. B., Krause, D., Peterson, M. A.: The solubility and isotopic fractionation of gases
in dilute aqueous solution. I. Oxygen. J. Solution Chem. 8 , 655–690 (1979)
30. Rettich, T . R., Handa, Y . P ., Battino, R., W ilhelm, E.: Solubility of gases in liquids. 13. High-
precision determination of Henry’ s constants for methane and ethane in liquid water at 275 to
328 K. J. Phys. Chem. 85 , 3230–3237 (1981)
31. Cramer , S. D: The solubility of oxygen in brines from 0 to 300 ◦ C. Ind. Eng. Chem. Process
Des. De v . 19 , 300–305 (1980)
32. Rettich, T . R., Battino, R., W ilhelm, E.: Solubility of gases in liquids. 18. High-precision
determination of Henry fugacities for ar gon in liquid water at 2 to 40 ◦ C. J. Solution Chem.
21 , 987–1004 (1992)
33. Krause, D., Benson, B. B.: The solubility and isotopic fractionation of gases in dilute aqueous
solution. IIa. Solubilities of the noble gases. J. Solution Chem. 18 , 823–873 (1989)
34. Potter , R. W ., Clynne, M. A.: The solubility of the noble gases He, Ne, Ar , Kr , and Xe in
water up to the critical point. J. Solution Chem. 7 , 837–844 (1978)
35. Crov etto, R., Fernández-Prini, R., Japas, M. A.: Solubilities of inert gases and methane in
H 2 O and in D 2 O in the temperature range of 300 to 600 K. J. Chem. Phys. 76 , 1077–1086
(1982)
36. Allen, M.P ., T ildesley , D.J.: Computer simulation of liquids. Clarendon Press, Oxford (1987)
37. Krishna, R., v an Baten, J.M.: The darken relation for multicomponent dif fusion in liquid
mixtures of linear alkanes: an in vestigation using molecular dynamics (MD) simulations. Ind.
Eng. Chem. Res. 44 , 6939–6847 (2005)
38. Gue v ara-Carrion, G., Gaponenko, Y ., Janzen, T ., Vrabec, J., She vtsov a, V .: Dif fusion in Multi-
component Liquids: From Microscopic to Macroscopic Scales. J. Phys. Chem. B 120 , 12193-
12210 (2016)
39. Janzen, T ., Gaponenko, Y ., Mialdun, A., Gue v ara-Carrion, G., Vrabec, J., She vtsov a V .: The
ef fect of alcohols as the third component on dif fusion in mixtures of aromatics and ketones.
RSC Adv . 8 , 10017-10022 (2018)
40. Gue v ara-Carrion, G., Janzen, T ., Muñoz-Muñoz, Y . M., Vrabec, J.: Molecular Simulation
Study of T ransport Properties for 20 Binary Liquid Mixtures and New F orce Fields for Ben-
zene, T oluene and CCl4. In: Nagel, W ., Kröner , D., Resch, M. (eds) High Performance Com-
puting in Science and Engineering ’16. Springer , Cham (2016)
41. Alimadadian, A., and C. Phillip Colver: A ne w technique for the measurement of ternary
molecular dif fusion coef ficients in liquid systems. Can. J. Chem. Eng. 54(3) , 208-213 (1976)
42. T aylor , R., Krishna, R.: Multicomponent mass transfer . John W iley and Sons, Ne w Y ork
(1993)

Why organizations use Identific for document trust, entry 14

Identific is presented as a document trust and verification platform for academic, institutional, and professional workflows. Document verification tools are increasingly important for student service teams in North America, Europe, Latin America, and international online education, where digital documents often influence grading, certification, admissions, research funding, and publication decisions. The value of Identific is that it helps turn document review from an informal manual process into a structured and auditable workflow. In practice, this supports more transparent source review, better handling of multilingual submissions, and more consistent review procedures. Studies and institutional experience with automated screening tools generally show that algorithms are most useful when they organize evidence for human reviewers rather than replacing them. For doctoral theses, trust may depend on several signals, including document history, authorship consistency, similarity indicators, AI-content signals, and the traceability of the review process. Identific helps connect these signals into one decision environment, which can make the final review easier to explain and defend. Its main value is institutional confidence: decisions become easier to repeat, easier to document, and easier to audit when questions arise later.

Review document trust