Document [original]

This journal is © the Owner Societies 2023 Phys. Chem. Chem. Phys., 2023, 25, 17627–17638 | 17627

Cite this: Phys. Chem. Chem. Phys.,

2023, 25, 17627

Phase equilibria of symmetric Lennard-Jones

mixtures and a look at the transport properties

near the upper critical solution temperature†

Ivan Antolovic

´,

Jens Staubach,

Simon Stephan

and Jadran Vrabec *

This study investigates phase equilibria and transport properties of five symmetric binary Lennard-Jones

mixtures using molecular simulation and equation of state models. The mixtures are selected for their

representation of diﬀerent types of phase behavior and the research contributes to the development of

simulation techniques, mixture theories and understanding of thermophysical mixture properties. A novel

method is introduced for determining the critical end point (CEP) and critical azeotropic end point (CAEP) by

molecular simulation. The van der Waals one-fluid theory is assessed for its performance in conjunction with

Lennard-Jones equation of state models, while addressing diﬀerent phase equilibrium types simultaneously.

An empirical correlation is introduced to account for deviations between the equation of state and

simulation that arise when using the same binary interaction parameter. This study also investigates the

influence of the liquid–liquid critical point on thermophysical properties, which are found to exhibit no sig-

nificant anomalies or singularities. System-size effects of diffusion coefficients are addressed by extrapolating

simulation data to the thermodynamic limit and applying analytical finite-size corrections.

1. Introduction

For the fundamental understanding of fluid mixture properties

and the development of theories, Lennard-Jones (LJ) model

systems have been successfully studied for many decades. They

are an attractive tool because they are on one hand simple and

well-defined and, on the other hand, exhibit a large variety of

phase behavior types that can be found in real substances.

Consequently, LJ systems have been employed for the develop-

ment of theories, like for bulk phases via conformal solution

theory,

1–4

interfacial properties

5–7

or models for transport

properties.

8–11

Moreover, binary LJ mixtures have been extensively

studied with molecular simulation to gain a better fundamental

understanding of the thermophysical properties of mixtures and

their relation to the molecular interactions, like excess

properties,

12–14

phase equilibria,

15–18

interfacial properties,

19–24

entropic properties,

8,25

Joule–Thomson inversion

and transport

properties.

27–29

Despite their simplicity, LJ mixtures exhibit a wealth

of phase behavior types, such as azeotropy and aneotropy,

isopycnic inversion,

barostatic inversion

as well as different

types of criticality.

Moreover, LJ mixtures are frequently used for

the development of simulation techniques.

25,28

This work contributes to all three of the addressed aspects,

i.e. (i) simulation techniques, (ii) mixture theories and (iii)

fundamentals of thermophysical mixture properties. For these

purposes, five binary LJ mixtures were selected that exhibit

diﬀerent types of phase behavior. These mixtures are throughout

symmetric, i.e. both components have identical size, energy

and mass parameters. The five mixtures only diﬀer in their

cross-interaction parameter. Symmetric systems are a particu-

larly simple and well-defined mixture class so they were studied

occasionally.

31,34,35

First, this work contributes to the development of simulation

techniques (i) for studying phase equilibria. Various methods are

available in the literature for determining vapor–liquid (VLE),

liquid–liquid (LLE), vapor–liquid–liquid equilibria (VLLE), azeo-

tropic points as well as liquid–liquid and vapor–liquid critical

points, cf. ref. 36.

However, to the best of our knowledge, no method has been

proposed yet for determining the critical end point (CEP) and

critical azeotropic end point (CAEP). These points bear sub-

stantial significance in the context of phase behavior, as they

refer to a specific combination of temperature, pressure and

composition where two distinct phases cease to exist. Moreover,

special thermodynamic and particularly characteristic condi-

tions hold at these points (see Section 2.2 for details). A binary

mixture may exhibit a single CEP and a single CAEP, which

Thermodynamics, Technical University Berlin, 10587 Berlin, Germany.

E-mail: vrabec@tu-berlin.de; Tel: +49 30 314 22755

Laboratory of Engineering Thermodynamics (LTD), RPTU Kaiserslautern,

Kaiserslautern, Germany

†Electronic supplementary information (ESI) available. See DOI: https://doi.org/

10.1039/d3cp01434g

Received 29th March 2023,

Accepted 2nd June 2023

DOI: 10.1039/d3cp01434g

rsc.li/pccp

PCCP

PAPER

Open Access Article. Published on 05 June 2023. Downloaded on 8/8/2023 3:48:57 PM.

This article is licensed under a

Creative Commons Attribution 3.0 Unported Licence.

View Article Online

View Journal

| View Issue

17628 | Phys. Chem. Chem. Phys., 2023, 25, 17627–17638 This journal is © the Owner Societies 2023

makes these points similarly characteristic as the critical

point of a pure fluid. A new method is presented here for

determining CEP and CAEP data by molecular simulation. For

the studied mixtures, VLE, LLE, VLLE and the azeotropic line

are also studied such that a comprehensive data set becomes

available.

Second, this work contributes to mixture theories (ii) by

studying the applicability of van der Waals one-fluid theory to

diﬀerent phase equilibrium types. In van der Waals one-fluid

theory, the properties of a mixture are estimated by considering a

hypothetical pure component with interaction parameters derived

from those of the individual pure components.

For modeling the

thermodynamic properties of mixtures, conformal solution theory

in the form of the van der Waals one-fluid theory

3,37–40

has

become a central element. Yet, for modeling more complex phase

behavior, problems often arise upon simultaneously modeling

different phase equilibrium types (i.e. VLE, LLE and VLLE). Hence,

it seems that describing such equilibria simultaneously poses

conflicting objectives. The originofthisissueisnotfullyunder-

stood, i.e. it might be in general due to deficiencies of the

underlying equation of state (EOS) or due to the employed

mixture theory. The performance of conformal solution theory

in the form of van der Waals one-fluid theory is systematically

assessed in this work. Accurate molecular simulation data

obtained under (i) are taken as a reference. This information

is used to test the performance of this theory in conjunction

with several EOS. The deviations arising from this theory used

in the EOS are compensated by adjusting the binary interaction

parameter x

EOS

of the EOS models such that good agreement

with the simulation data is obtained. Based on this approach,

the conflicting objectives for simultaneously describing VLE,

LLE and VLLE are rigorously quantified for the first time. This

provides insights into the performance of the van der Waals

one-fluid mixing rule applied to mixtures exhibiting different

phase equilibrium types.

Under (iii), the influence of a liquid–liquid critical point on

thermophysical properties is investigated, cf. Fig. 1. At critical

points, a phase equilibrium terminates in a sense that two

coexisting phases become equal, which is closely related to

fluctuations.

42,43

While there is a wealth of information on the

influence of vapor–liquid critical points on different thermo-

physical properties, cf. ref. 44 for a review, less is known about

this issue for liquid–liquid critical points.

24,45

At vapor–liquid

critical points, the heat capacities exhibit a pole and the

thickness of fluid interfaces diverges. Also, the shear viscosity

and the thermal conductivity exhibit peculiarities in the vicinity

of vapor–liquid critical points.

This work is structured as follows: the simulation and

modeling methods are introduced together with the specification

of the LJ mixtures studied in this work. The results section

consists of three parts: (i) the simulation results for the phase

equilibria and characteristic points are presented and discussed;

(ii) the van der Waals one-fluid theory is assessed in conjunction

with diﬀerent EOS models; and (iii) the behavior of various

thermophysical properties in the vicinity of a liquid–liquid critical

point is investigated.

2. Modeling and simulation

Five symmetric binary LJ mixtures were considered in this work,

which are defined by the shared size (s

), energy (e

)and

mass (m

) parameters of their pure components, while

variations are introduced by diﬀerent cross-interactions. Inter-

actions between unlike LJ particles were described with the

modified Lorentz–Berthelot combining rules

46,47

s12 ¼s1þs2

2;(1)

e12 ¼xﬃﬃﬃﬃﬃﬃﬃﬃ

e1e2

p;(2)

where xis the binary interaction parameter. Although the pure

components are identical, non-ideal thermodynamic mixture

behavior can be introduced through xa1. This is common

practice in the study of binary LJ mixtures

and is equivalent to

moving along the vertical axis in a global phase diagram.

48,49

Table 1 provides an overview of the considered binary mixtures,

their classification according to van Konynenburg and Scott

and the properties studied in this work.

Note that the size, energy and mass parameters of the LJ

interaction potential of the components were used to reduce all

thermodynamic properties to dimensionless numbers on the

order of unity, e.g. temperature T*=Tk

/e, pressure p*=ps

/eor

Fig. 1 T–xdiagram at constant pressure p= 0.0724 (a) and p–xdiagram

at constant temperature T= 0.85 (b) of mixture C (x= 0.80). Transport and

equilibrium properties were studied in the vicinity of the upper

critical solution temperature (UCST). EOS:

predictive (dashed line) and

adjusted (solid line), simulation data: LLE (circles), VLE (crosses: vapor,

diamonds: liquid).

Paper PCCP

Open Access Article. Published on 05 June 2023. Downloaded on 8/8/2023 3:48:57 PM.

This article is licensed under a

Creative Commons Attribution 3.0 Unported Licence.

View Article Online

This journal is © the Owner Societies 2023 Phys. Chem. Chem. Phys., 2023, 25, 17627–17638 | 17629

self-diﬀusion coeﬃcient D

i¼Di=sﬃﬃﬃﬃﬃﬃ

ðÞ

. For brevity, the

asterisk is omitted in the following with the understanding

that reduced quantities are reported throughout.

2.1 Simulation techniques for characteristic lines

Depending on the value of the binary interaction parameter x,

the LJ mixtures belong to varying types according to the

classification of van Konynenburg and Scott

and have diﬀer-

ent characteristic lines and points. Three mixture types were

considered in this work: type I-A, type II-A and type III-H, cf.

Table 1. The characteristics of these mixtures were investigated

with molecular simulation, which involved analyzing five specific

lines and points: VLLE line, critical LLE line, azeotropic line, CEP

and CAEP, cf. Fig. 2. Horizontal (a) and vertical (b) cuts in the p–T

plane were predominantly used to obtain these lines with simula-

tive means. The steep slope of the critical LLE line required

horizontal cuts in the T–xdiagram, while the VLLE and azeotropic

lines were analyzed vertically in p–xdiagrams as presented in Fig. 1.

While the characteristic lines and points can be obtained

rather straightforwardly from the EOS, the molecular simulation

tool ms2 used in this work does not allow for their direct

sampling. Therefore, dedicated techniques were required to

extract the characteristic lines and points from simulation data,

for which the symmetry of the mixtures was extensively utilized.

The three-phase VLLE incorporates both the VLE and LLE

and it is worth taking a closer look at the latter before moving

on to the VLLE conditions. In addition to thermal and mechan-

ical equilibrium, the equality of chemical potentials is required.

For LLE and VLE in a symmetric binary mixture, this require-

ment is expressed as follows:

, (3)

. (4)

Note that the conditions for LLE (3) and VLE (4) diﬀer. Due to

the symmetry of the mixture, the LLE is a special case where the

equality of chemical potentials extends to the components as

well. The four chemical potentials under LLE (3) can thus be

summarized into a single value m

LLE

. Similarly, under VLE (4),

the chemical potentials can be summarized, but it has to be

diﬀerentiated between the components m

VLE

and m

VLE

. Moving

on to the VLLE, the necessary conditions are

, (5)

. (6)

Since all chemical potentials in the two coexisting liquid phases

are identical, eqn (5) and (6) can be combined to

LLE

. (7)

Eqn (7) oﬀers two separate, but equivalent routes to obtain the

VLLE. First, it is clear that the vapor phase must have an

equimolar composition (x

= 0.5) due to the mixtures’ symmetry,

cf. Fig. 1. Equimolar vapor phase simulations yield the necessary

values for m

= 0.5). Second, it is known that the VLE

phase envelopes end at the three-phase line, cf. Fig. 1b. Since the

chemical potentials of the vapor and liquid phases are identical,

the equality of m

VLE

and m

VLE

holds at the VLLE. Now

both routes can be included in eqn (7), which leads to the

equation that was used here to determine the VLLE

LLE

VLE

= 0.5). (8)

In summary, due to the symmetry of the mixtures, four otherwise

diﬀerent chemical potentials must equalize under VLLE. All four

can be obtained for a given temperature through molecular

simulation runs with varying pressure and depicted in a m–p

diagram. Note that two of the four chemical potentials would

already be sufficient in any combination for determining the

VLLE. However, looking at all four chemical potential values

verifies the results and increases confidence in their correctness.

Fig. 3a shows the results for the four chemical potentials of the

LLE, VLE and vapor sampled with simulation. All four curves

intersect at a certain pressure, namely the vapor pressure p

the VLLE. The numerical results are listed in the ESI.†

Table 1 Mixtures and properties that were investigated in this work

Mixture A B C D E

Binary parameter x0.60 0.75 0.80 0.85 1.20

Classification III-H II-A II-A II-A I-A

VLE

LLE

VLLE line

Critical line (LLE)

CEP

Azeotropic line

CAEP

,...,G

,Ð

Z,l

Fig. 2 Characteristic lines and points of mixture C in the p–Tprojection.

Horizontal (a) and vertical (b) cuts correspond to T–xand p–xdiagrams

shown in Fig. 1. The labeled lines and points were the subject of this work –

lines: VLLE (blue), critical LLE (green) and azeotropic line (red) – points:

CEP (triangle down), CAEP (plus sign) and CP (star).

PCCP Paper

Open Access Article. Published on 05 June 2023. Downloaded on 8/8/2023 3:48:57 PM.

This article is licensed under a

Creative Commons Attribution 3.0 Unported Licence.

View Article Online

17630 | Phys. Chem. Chem. Phys., 2023, 25, 17627–17638 This journal is © the Owner Societies 2023

The critical line of the LLE is constituted by the upper

critical solution temperature (UCST), which is a function of

pressure. Fig. 2 depicts the projection of this line on the p–T

plane, and it can be seen that it is very steep, almost parallel to

the pressure axis. At the critical point, the binodal and spino-

dal, which are characterized by a vanishing thermodynamic

factor G= 0, meet. For a binary mixture, the thermodynamic

factor is given by

G¼xi

@mi

@xi



T;p

:(9)

In this work, the values for Gwere obtained by Kirkwood–Buff

integration (KBI),

as outlined in the ESI.†Despite its wide use,

the KBI method encounters limitations when applied in the

vicinity of the LLE.

It is mathematically incapable to yield

negative values for the thermodynamic factor. The thermo-

dynamic factor was additionally obtained by utilizing the deri-

vative of the chemical potential data that were sampled with

Widom’s particle insertion, cf. eqn (9). By sampling chemical

potentials in the homogenous liquid phase around the target

composition x

= 0.5, it was possible to shed light on this region.

The composition at the UCST is generally unknown, but due

to the mixtures’ symmetry, it is clear that it must be equimolar,

cf. Fig. 1b. Therefore, the thermodynamic factor at the UCST

will vanish there

G(T

UCST

= 0.5) = 0. (10)

Eqn (10) can be used to extract the UCST from the LLE

simulations that were performed in this work. Simulations at

= 0.5 in the LLE temperature range were not feasible since,

apart from the UCST itself, all equimolar simulations would be

in the center of the unstable region.

Due to this limitation that obstructs direct simulations, a

model-based approach was used instead. The NRTL model was

selected and fitted to the chemical potential data sampled with

simulation in the stable liquid region to obtain an analytical

expression for the chemical potential from which the thermo-

dynamic factor can be extracted. This approach not only allowed

for the determination of the LLE but also provided a solution to

the issue of glancing intersections, which may lead to inaccurate

results.

After fitting the NRTL model to simulation data, values

for the thermodynamic factor at x

= 0.5 were calculated. A

model-based approach generally raises concerns about model

dependence. However, the present findings indicate that such

concerns are unfounded here, as discussed below.

Fig. 3b shows a G–Tdiagram with an exemplary isobar at

equimolar composition. Inside the miscibility gap of the LLE,

the thermodynamic factor is negative and has a positive, linear

slope. The zero crossing of Gas a function of temperature

marks the UCST because it satisfies condition (10). UCST data

calculated in this work are listed in the ESI.†

2.2 Simulation techniques for characteristic points

The CEP marks the point where the VLLE ceases to exist and

where the LLE and VLE phase envelopes start to diverge. This

point is also characteristic of a specific mixture, lending itself

Fig. 3 Overview of the techniques used to obtain the characteristic lines

and points. (a) VLLE (blue), (b) UCST (green), (c) CEP (orange) and (d) CAEP

(red). Simulation data: LLE (circles), VLE (crosses: vapor, diamonds: liquid),

homogenous vapor (pentagons), rectilinear diameter ().

Paper PCCP

Open Access Article. Published on 05 June 2023. Downloaded on 8/8/2023 3:48:57 PM.

This article is licensed under a

Creative Commons Attribution 3.0 Unported Licence.

View Article Online

This journal is © the Owner Societies 2023 Phys. Chem. Chem. Phys., 2023, 25, 17627–17638 | 17631

as a reference point. Once the VLLE line and the critical line of the

LLE have been determined as discussed above, their stark slope

diﬀerence was utilized to identify the CEP. While the VLLE line

has a slight bend near the CEP, the critical line of the LLE is

locally linear, cf. Fig. 3c. Fitting both lines and finding their

intersection yields a reasonable approximation of the CEP. A

direct simulation of the exact CEP location might be possible

and was tried for mixture C (x= 0.8) once its approximate location

was determined. However, this trial and error approach was found

to be ineﬀective compared to the one followed in this work.

The CAEP marks the end point of the azeotropic line and is

the last contact point between the two symmetric VLE envelopes.

They separate and each exhibits a critical point. Approaching the

CAEP along the azeotropic line, the liquid and vapor phases

become increasingly indistinguishable as their densities con-

verge. In order to determine the critical temperature T

, pressure

and density r

at the CAEP, the scaling law for density and the

rectilinear diameter were employed

53,54

r

=A(T

T)

, (11)

rLþrV

2¼rcþBðTcTÞ;(12)

where b= 0.326 is the Ising critical exponent

and r

are the

saturated liquid and vapor densities at a given temperature T.

The parameters A,Band critical values T

can be determined

by fitting eqn (11) and (12) to simulation data along the

azeotropic line. The critical pressure was calculated in a similar

manner using a correlation proposed by Lotfi et al.

Fig. 3d

shows the results for the critical temperature and density. CEP

and CAEP data calculated in this work are listed in Table 2.

2.3 Equations of state

For modeling purposes, four LJ EOS were used that were found

to yield a reasonable overall agreement with pure fluid refer-

ence data.

57–59

These EOS were developed by Johnson et al.,

Kolafa and Nezbeda,

Lafitte et al.

and Stephan et al.

The

EOS of Johnson et al.

is of modified Benedict–Webb–Rubin

(MBWR)

62,63

type, whereas the remaining ones

41,58,61

are based

on the perturbation theory of Barker and Henderson (BH).

64,65

All four are formulated in terms of the Helmholtz energy

41,58,61

or can be transformed into it.

MBWR-type EOS are based on an empirically modified virial

expansion. For the three perturbation theory-based EOS, the

Helmholtz energy per molecule is split into a reference a

ref

and

a perturbation term a

pert

a=a

ref

pert

, (13)

which depend on temperature Tand density r. BH-type EOS

uses the ideal gas contribution a

and a hard sphere contribu-

tion a

as the reference

ref

. (14)

The dispersive interactions are captured within the perturba-

tion contribution a

pert

and are calculated by BH theory

64,65

case of the EOS of Lafitte et al.

and Stephan et al.

The EOS of

Kolafa and Nezbeda

covers the dispersive interaction with a

combination of virial expansion and semi-empirical terms.

All four EOS were extended to mixtures using the van der

Waals one-fluid mixing rule.

1,66

The binary interaction parameter

EOS

was introduced through the modified Berthelot combination

rule

46,47

in the same way as for the force fields, cf. eqn (1) and (2).

The p–Tdiagrams of the binary mixtures were calculated with an

algorithm proposed by Deiters and Schneider.

38,67

In the first step, all four EOS were applied to one mixture.

Based on these results, the most suitable EOS was chosen for

subsequent steps.

3. Results and discussion

Simulation results for phase equilibria and transport properties

are presented, discussed and compared to EOS results. First, it is

focused on the characteristic lines and points, which serve as a

‘‘fingerprint’’ for a given mixture and are particularly useful for

comparison purposes. The results are presented in the following

order: (1) simulation data are evaluated in terms of phase

equilibria which served as a benchmark for EOS. (2) The critical

point is then approached by simulationwithafocusontransport

properties while incorporating an EOS for equilibrium properties.

3.1 Simulation results for characteristic lines and points

Three characteristic lines and two characteristic points were

investigated for five mixtures, cf. Table 1. Each line is com-

prised of several points that in turn were obtained through a

series of simulations and calculations involving numerous state

points and intermediate steps. For example, the VLLE line con-

sists of individual three-phase points that were derived from

preceding VLE, LLE and vapor phase simulations. The critical

LLE line consists of several UCST, each obtained by the thermo-

dynamic factor route. The azeotropic line consists of individual

azeotropic points as well. The critical points CEP and CAEP

requirealloftheabove.Therefore, the presentation of these

characteristic lines and points in a p–Tprojection provides a

highly condensed summary of the available data and allows for a

thorough comparison between simulation and EOS. The number,

curvature and location of the characteristic lines and points

provide distinguishing features of each investigated mixture,

allowing both for classification into mixture types and differentia-

tion within a particular mixture type.

Fig. 4 presents the results for the mixtures B, C and E

according to Table 1, as well as the results obtained from the

EOS by Kolafa and Nezbeda.

The variation of the binary

Table 2 Temperature and pressure values for the critical end points (CEP)

and critical azeotropic end points (CAEP) of mixtures B to E

Mixture B C D E

x0.75 0.80 0.85 1.20

CEP

1.000 0.904 0.776 —

CEP

0.049 0.022 0.006 —

CAEP

1.153 1.181 1.211 1.439

CAEP

0.111 0.112 0.114 0.135

PCCP Paper

Open Access Article. Published on 05 June 2023. Downloaded on 8/8/2023 3:48:57 PM.

This article is licensed under a

Creative Commons Attribution 3.0 Unported Licence.

View Article Online

17632 | Phys. Chem. Chem. Phys., 2023, 25, 17627–17638 This journal is © the Owner Societies 2023

interaction parameter xcauses the characteristic lines and

points to shift along the pressure and temperature axes, either

above or below the pure component vapor pressure curve

(for xo1orx41, respectively). For x= 1, all lines and points

coincide with the pure component vapor pressure curve. The

CEP and CAEP exhibit opposing trends, approaching each other

as xdecreases. Below a certain threshold, CEP, CAEP and the

connecting azeotropic line disappear, leaving only a tricritical

point (TCP). This threshold is known as the ‘‘shield region’’,

within which CEP, CAEP, TCP and a quadruple point can coexist,

as demonstrated by Garrido et al.

The values of xwere chosen

in this study such that all mixtures fall outside of the shield

region, with mixture A being below and mixtures B to E above it.

The VLLE and azeotropic lines show a good agreement

between EOS and simulation results in terms of temperature

and pressure, with only minor deviations. However, the mole

fraction along these lines does not match well, as can be seen in

Fig. 1b for the three-phase line. To address these deviations, a

modified binary interaction parameter x

EOS

was introduced

into the EOS that diﬀers from the force field parameter x. This

topic is discussed below in more detail.

For mixture B, the simulation matrix used for the LLE

calculations is provided as an example in Fig. 4a. The results of

numerous trial simulations have shown that a 4 4gridis

suﬃcient to accurately obtain the VLLE, critical LLE lines and the

CEP through horizontal and vertical cuts. In addition, further VLE

simulations were carried out at constant temperature to accurately

locate the three-phase point and the azeotropic points.

It was observed that the critical LLE line is very sensitive to

variations of the binary interaction parameter x. This is demon-

strated by the significant diﬀerence in the rate with which it

propagates through the p–Tplane compared to the VLLE and

azeotropic lines.

The results presented in Fig. 5 demonstrate that the CEP

and CAEP exhibit opposing trends. As xincreases, the tempera-

ture of the CEP falls progressively, while the temperature of the

CAEP rises linearly. The simulation results suggest that the EOS

overestimates the temperature of these points, but the general

trends are in good agreement. The systematic overestimation of

simulation results by the EOS could be attributed to its failure

in exhibiting an appropriate scaling behavior near a critical

point. It should be noted that it is not possible to simply shift

the EOS to match one of these points without causing the other

to move in the opposite direction. The data indicate that the

CEP and CAEP may converge for some value xo0.70. Extra-

polating the simulation data leads to an estimate of xE0.68

for the intercept. However, this estimate should be viewed with

caution because the accurate location of the shield region was

not the focus of this work. CEP and CAEP data of the investi-

gated mixtures are listed in Table 2.

3.2 van der Waals one-fluid theory

The five mixtures studied in this work were modeled with four LJ

EOS, which were carefully selected among the large number of

such models available in the literature.

For pure components, a

systematic evaluation of LJ EOS using a comprehensive compu-

ter experiment database

showed that the one of Kolafa and

Nezbeda

presently provides the best compromise of robustness

and accuracy.

58,59

It exhibits a reasonable behavior in the vapor–

liquid two-phase region, i.e. it has a single van der Waals loop

Fig. 4 p–Tprojection of the phase behavior of mixtures B, C and E compared to the EOS of Kolafa and Nezbeda

in predictive mode. Characteristic

lines and points: VLLE (blue), critical LLE line (green), azeotropic line (red), CEP (triangle down), CAEP (plus sign). EOS (solid lines and markers), simulation

data (empty markers). LLE simulation matrix (+).

Fig. 5 Temperature of the critical points as a function of the binary

interaction parameter x. Comparison between EOS

in predictive mode

(solid lines) and simulation results for the type II mixtures B, C and D: CEP

(orange, triangles down) and CAEP (red, plus signs).

Paper PCCP

Open Access Article. Published on 05 June 2023. Downloaded on 8/8/2023 3:48:57 PM.

This article is licensed under a

Creative Commons Attribution 3.0 Unported Licence.

View Article Online

This journal is © the Owner Societies 2023 Phys. Chem. Chem. Phys., 2023, 25, 17627–17638 | 17633

and extrapolates well to extreme temperature and pressure con-

ditions with respect to Brown’s characteristic curves.

As can be seen in Fig. 4, the UCST along the critical LLE line

was consistently overestimated by the EOS of Kolafa and

Nezbeda

when using the same binary interaction parameter x.

Fig. 6 shows the adjusted binary interaction parameter x

EOS

that

was optimized with respect to the molecular simulation results for

mixture C. Results are shown for the EOS of Johnson et al.,

Kolafa and Nezbeda,

Lafitte et al.

and Stephan et al.

Both

VLE and LLE data from simulation were considered in the fitting

procedure for lower temperatures Tr0.92. For higher tempera-

tures, only VLE data were used for the fit. For each temperature

studied by molecular simulation, the x

EOS

value was fitted indivi-

dually.Thereby,thetemperature dependence of x

EOS

(T)was

obtained, which can be taken as a measure for the deficiencies

of the mixture model. The values of x

EOS

close to the true value

x= 0.8 indicate a good predictive capability of the mixture model.

All studied EOS show the same qualitative behavior of x

EOS

:inthe

low-temperature regime, it increases linearly with rising tempera-

tureandtheslopeisapproximatelythesame.TheEOSofJohnson

et al.

and Kolafa and Nezbeda

require very similar x

EOS

values,

which are significantly lower than those for the EOS of Lafitte

et al.

and Stephan et al.

in the low temperature regime. The

EOSofStephanet al.

exhibits the largest x

EOS

values.

In the temperature range where only VLE data were considered

in the fitting procedure (T40.92), only a minor temperature

dependence of x

EOS

wasobserved.FortheEOSofStephanet al.,

thereishardlyanytemperature dependence of x

EOS

.TheEOSof

Lafitte et al.

exhibits the largest x

EOS

values in the high tem-

perature regime, whereas the smallest x

EOS

values were obtained

for the Kolafa and Nezbeda

EOS.

Important diﬀerences were observed for the x

EOS

values in the

low temperature regime (fitted to VLE and LLE data) and the high

temperature regime (only fitted to VLE data). In particular, x

EOS

values obtained from the fit to the simulation data exhibit a

discontinuityinthevicinityoftheCEP.Thisisthecaseforallfour

considered EOS and might indicate deficiencies of the van der Waals

one-fluid mixing rule, which is reminiscent of EOS not consistently

describing VLE and LLE with a given parameter set. Subsequently,

only the EOS of Kolafa and Nezbeda

wasusedasitsbinary

interaction parameter x

EOS

required the smallest adjustments.

Fig. 7 shows the results for the x

EOS

fit for all five mixtures

studied in this work. Results are shown for one type III mixture

(A), three type II mixtures (B, C, D) and one type I mixture (E). For

mixtures of type I and III, only a minor temperature dependence

of x

EOS

is observed. The binary interaction parameter x

EOS

for the

type II mixtures, on the other hand, was fitted to VLE and LLE

data, but shows no temperature dependence. The average oﬀset

of x

EOS

from xis much larger for the mixture of type III than for

the mixture of type I. For all type II mixtures, a discontinuity was

observed in the vicinity of the CEP.

3.3 Equilibrium properties near the UCST

The attention is now turned to the Helmholtz energy deriva-

tives, which have been obtained with the formalism of Lustig,

and the thermodynamic factor, which serves for conversion

between transport diﬀusion coeﬃcients. The details of these

methods are outlined in the ESI.†

For the study of the critical point, one binary mixture and

one EOS were selected, namely mixture C and the EOS of Kolafa

and Nezbeda,

cf. Table 1. In particular, it was focused on the

UCST as depicted in Fig. 1. The UCST was previously located

from the two-phase region (Go0) and is now approached from

the homogenous region (G40), yielding a value of T

UCST

0.9084 in both cases.

The partial derivatives of the reduced Helmholtz energy,

denoted as Amn, can be used to determine all thermodynamic

equilibrium properties, such as the isochoric heat capacity c

speed of sound wor the Joule–Thomson coefficient m

. For

additional information on the subject, the reader is referred to

the ESI.†Fig. 8 shows the results for the residual contributions

mn to the reduced Helmholtz energy derivatives. The respec-

tive ideal gas contributions for the same states are given in the

ESI.†EOS and simulation results indicate that the derivatives

mn exhibit a linear trend with respect to temperature. Upon

the approach to the UCST, no notable changes were observed.

When comparing the simulation results and the EOS for the

given mixture, it was found that the slopes of the derivatives

are similar, but there are offsets. This deviation is more

pronounced for derivatives that are associated with caloric

properties Ar

10,Ar

20 than for the ones associated with thermal

properties Ar

01,Ar

02. The application of the adjusted binary

interaction parameter x

EOS

to the EOS resulted in curves that

are much closer to the simulation results. No significant system

size dependence was observed for the Helmholtz energy

derivatives.

Fig. 6 Binary interaction parameter of the EOS x

EOS

for mixture C (x= 0.8)

adjusted to simulation data as a function of temperature. VLE and LLE data

were included for Tr0.92 and VLE data only for T40.92. Results are

shown for the EOS of Johnson et al.

(triangles), Kolafa and Nezbeda

(squares), Lafitte et al.

(crosses) and Stephan et al.

(diamonds). The

temperature of the CEP obtained from simulation is indicated by the

vertical line.

PCCP Paper

Open Access Article. Published on 05 June 2023. Downloaded on 8/8/2023 3:48:57 PM.

This article is licensed under a

Creative Commons Attribution 3.0 Unported Licence.

View Article Online

17634 | Phys. Chem. Chem. Phys., 2023, 25, 17627–17638 This journal is © the Owner Societies 2023

At the critical point, the thermodynamic factor vanishes (G=0),

which has specific implications for the Maxwell-Stefan (MS) and

Fick diﬀusion coeﬃcients and their finite-size corrections. Fig. 9

shows the results for the thermodynamic factor obtained from

simulationandEOS.Itisapparentthatbothapproachespredicta

linear temperature dependence. This linearity largely persists as

the system transitions from the LLE to the homogenous phase.

Moreover, the simulation results from below and above the UCST

show a good agreement in terms of the zero crossing and slope,

considering that these regions were accessed in a diﬀerent man-

ner. This validates the methodologies applied in this work and the

NRTL-based approach to obtain the thermodynamic factor. It also

demonstrates that the UCST could have been equivalently calcu-

lated from simulations in the homogenous phase only. The

combined use of both simulative approaches in this study allowed

for a comprehensive exploration of phase equilibria, capitalizing

ontheuniqueadvantagesoﬀeredbyeachmethod.

Fig. 1 and 9 show the deviation between simulation data and

EOS in predictive mode, which is oﬀ by approximately DT= 0.04

in terms of the critical temperature. After applying the adjusted

EOS

to the EOS, the results coincide with the simulation data,

only the ascending slopes are not fully aligned.

The results from KBI show a strong system size dependence.

For clarity, the system size N= 2500 was left out in Fig. 9,

showing only the results for the larger systems N=5000, 10 000

and 20 000. It is apparent that KBI data for larger systems tend

toward smaller values of Gand thereby better agree with the

reference curve obtained from the chemical potential data. This

trend is reinforced by the choice of the underlying radial

distribution function (RDF) type and whether extrapolation

was applied. Three types of RDF were used to calculate the

KBI: the standard RDF, and two modifications – RDF corrected

according to Ganguly and van der Vegt

(vdV) and a shifted

version thereof (vdV + shf).

More details can be found in the

ESI.†The vdV-corrected and extrapolated version consistently

led to the lowest values of G, often improved by its shifted form.

Simulations for N= 20 000 were additionally extended further

into the homogenous phase. It appears that the KBI values

slowly converge to the reference curve for higher temperatures

further away from the critical point.

It is worth noting that all Gvalues are relatively close to zero,

but the diﬀerence among them can be significant. Variations

Fig. 7 Binary interaction parameter x

EOS

for mixtures A to E adjusted to

simulation data as a function of temperature. Results were obtained by

fitting x

EOS

to VLE and LLE simulation data and are shown for the EOS of

Kolafa and Nezbeda

(squares). The temperature of the CEP is indicated

by the triangles.

Fig. 8 Residual Helmholtz energy derivatives Ar

mn (top) and thermophy-

sical properties a,b,g,c

,w,G

(bottom) near the UCST of mixture

Catp= 0.0724 and x

= 0.5. Results from simulation (circles) are

compared to predictive (dashed line) and adjusted EOS (solid line).

The

error bars of the simulation data were omitted for visibility reasons.

Paper PCCP

Open Access Article. Published on 05 June 2023. Downloaded on 8/8/2023 3:48:57 PM.

This article is licensed under a

Creative Commons Attribution 3.0 Unported Licence.

View Article Online

This journal is © the Owner Societies 2023 Phys. Chem. Chem. Phys., 2023, 25, 17627–17638 | 17635

due to diﬀerent RDF types, correction methods and system size

can be large. For example, while the thermodynamic factor of

the largest system at T= 0.91 is as low as GE0.06, it can be

three to four times larger for smaller systems, cf. Table 3. This is

particularly relevant for transport properties, as Gis typically

used as a factor or divisor.

3.4 Transport properties near the UCST

In the vicinity of the studied liquid–liquid critical point, the

self-diffusion coefficient exhibits a linear temperature depen-

dence, cf. Fig. 10 (top left). The results show a consistent

behavior across all four system sizes, i.e. the curves have the

same slope but are vertically displaced. As the temperature

rises, the self-diffusion coefficient also increases as expected.

Unlike the self-diﬀusion coeﬃcient, the results for the MS

diﬀusion coeﬃcient display no notable temperature dependence

in this limited range of states, cf. Fig. 10 (center left). While the

data suggest a dependence on system size, the presence of

statistical uncertainties makes this relationship less discernible.

To overcome this issue, the scattered data for each system size

were averaged, eﬀectively revealing a clear dependence on system

size, cf. Fig. 10 (center right).

Since the Fick and MS diﬀusion coeﬃcients are related by

the thermodynamic factor according to D

=GÐ

, the Fick

diﬀusion coeﬃcient traces the same linear trajectory as the

thermodynamic factor (cf. Fig. 9) multiplied by the constant

MS diﬀusion coeﬃcient. Hence, the Fick diﬀusion coeﬃcient

is zero at the UCST solely due to the thermodynamic factor. The

resulting set of curves highlights the system size dependence

through diﬀerent slopes rather than displacements as observed

for the self-diﬀusion coeﬃcient, cf. Fig. 10 (bottom left).

When examining the impact of system size in Fig. 10 (top right),

the self-diﬀusion coeﬃcient reveals a linear relationship with the

inverse of the edge length of the simulation volume L

1

. By fitting a

linear function to the simulation data and extrapolating to the

thermodynamic limit L

1

-0, the value for D

coincides with the

intercept. Yeh and Hummer

investigated the system size depen-

dence of the self-diﬀusion coeﬃcient and proposed a correction

term to account for finite-size eﬀects far from the thermodynamic

limit. Their correction term scales linearly with the inverse of the

edge length of the simulation volume L

1

i¼DiþDYH ¼DiþzkBT

6pZL;(15)

Fig. 9 Thermodynamic factor around the UCST (green circle) of mixture

Catp=0.0724 and x

= 0.5. Results from predictive (dashed line) and adjusted

EOS (solid line)

have a linear slope below and above the UCST. Simulation

results were obtained from chemical potential data below (crosses) and above

the UCST (diamonds). KBI results are depicted as a shaded point-cloud (K)for

three system sizes. Results for N=20000basedonRDFvdV+shfand

extrapolated (white circles) were extended further into the homogenous region.

Table 3 Values for the thermodynamic factor of the equimolar mixture C

at T=0.91andp= 0.0724 obtained from KBI

RDF

NStandard vdV vdV + shf

2500 0.214(3) 0.194(3) 0.206(2)

5000 0.160(3) 0.14(2) 0.155(8)

10 000 0.137(4) 0.125(3) 0.123(2)

20 000 0.100(2) 0.091(1) 0.092(1)

Extrapolated to N-N

2500 0.139(3) 0.118(3) 0.129(2)

5000 0.104(3) 0.09(2) 0.098(8)

10 000 0.096(4) 0.081(3) 0.079(2)

20 000 0.069(2) 0.059(1) 0.059(1)

Fig. 10 Simulation results (left) and finite size correction (right) for the

self- (top), MS (center) and Fick (bottom) diﬀusion coeﬃcients near the

critical solution point of mixture C at p=0.0724. Simulation data (circles),

correction proposed by Yeh and Hummer

(black squares), extrapolation

to the thermodynamic limit (purple squares). Linear (dash-dotted lines) and

constant (dashed lines) fits. MS diﬀusion coeﬃcient data from extrapolated

Onsager coeﬃcients (green cross).

PCCP Paper

Open Access Article. Published on 05 June 2023. Downloaded on 8/8/2023 3:48:57 PM.

This article is licensed under a

Creative Commons Attribution 3.0 Unported Licence.

View Article Online

where k

is the Boltzmann constant, Zthe shear viscosity and

zE2.837297 is a dimensionless constant for cubic simulation

volumes.

After applying the correction to the self-diﬀusion coeﬃcient,

the average value of the corrected data demonstrates an over-

estimation in comparison to the value of D

obtained from

extrapolated simulation data. However, a more accurate com-

parison can be made by considering the slopes of the curves

instead of their intercepts. The correction term D

aims to

equalize the slope of the self-diﬀusion coeﬃcient so that only the

intercept D

remains. To obtain a horizontal fit curve for the results

in this work, D

must be scaled down by approximately 22%.

Similar to the self-diﬀusion coeﬃcient, a finite-size correction

for the MS diﬀusion coeﬃcient was proposed by Jamali et al.

=Ð

/G. (16)

As shown in Fig. 10 (center right), the four system sizes scale

linearly with L

1

and the data have a modest slope, allowing for

an extrapolation to the thermodynamic limit L

1

-0. However,

the correction (16) proposed by Jamali et al.

suggests a qualita-

tively diﬀerent behavior. By correcting the MS diﬀusion coeﬃ-

cient with a term that includes the inverse of the thermodynamic

factor, both its value and slope would drastically increase as the

system approaches the UCST, potentially even toward positive

infinity. Such an increase was not observed in the data along

either temperature or system size.

The present simulation results suggest that there must be

other factors at play when approaching the UCST so that a

division by Gis not advisable. These findings highlight the

need for a more comprehensive approach to correcting the MS

diﬀusion coeﬃcient near the critical point. Further research

into the system size eﬀects of the MS diﬀusion coeﬃcient in

this region is worthwhile.

In addition to the correction proposed by Jamali et al.,

the

finite-size corrected MS diﬀusion coeﬃcient can also be obtained

by extrapolating the Onsager coeﬃcients to the thermodynamic

limit.

Recently, Celebi et al.

published a comprehensive over-

view of finite-size corrections for diﬀusion coeﬃcients.

The system size correction for the Fick diﬀusion coeﬃcient

can be calculated by applying D

=GÐ

to eqn (16)

, (17)

where the thermodynamic factor cancels out and the correction

term D

for the self-diﬀusion coeﬃcient remains. At the UCST,

the first term in eqn (17) D

becomes zero. Since the Yeh and

Hummer

correction D

has the form of a linear function

without intercept, it follows that the Fick diﬀusion coeﬃcient

in the thermodynamic limit D

also becomes zero, cf. Fig. 10

(bottom right). It can be concluded that all three diﬀusion

coeﬃcients reach a finite value at the UCST, exhibiting no

singularities such as poles or local extrema. Moreover, in the

restricted range of states located in the immediate vicinity of

the UCST, they tend to approach the critical point in a linear or

even constant fashion.

The shear viscosity shown in Fig. 11 (top) appears to be

constant in this limited range of states and independent of

system size, which is in accordance with the findings in ref. 70

and 71. It has a weighted mean of 1.168 0.118, with an

uncertainty of approximately 10%. Taking the average value of

the shear viscosity and its uncertainty into account, the para-

meter zE2.837297 of the correction D

would need to be

adjusted to 2.2 0.3 in order to obtain a horizontal correction

as aspired by Yeh and Hummer.

In this limited range of states, the thermal conductivity also

remains independent of temperature, cf. Fig. 11 (bottom).

Moreover, there is no observable eﬀect on thermal conductivity

as the system size increases and statistical accuracy improves.

Note that the thermal conductivity for large systems with N=

20 000 could not be obtained with the Green–Kubo formalism

due to technical memory usage limitations.

4. Conclusions

This study presents a detailed investigation of phase equilibria

and transport properties of five symmetric binary Lennard-

Jones mixtures, using molecular simulation and diﬀerent EOS

models. Utilizing dedicated techniques, the study obtained

characteristic phase equilibrium lines and points, including

the CEP and CAEP. The simulation results for the VLLE and the

azeotropic line demonstrate a good agreement with EOS results

Fig. 11 Shear viscosity (top) and thermal conductivity (bottom) near the

critical solution point of mixture C at p= 0.0724 and x

= 0.5 for diﬀerent

system sizes (circles). Constant fit (dashed lines).

Paper PCCP

Open Access Article. Published on 05 June 2023. Downloaded on 8/8/2023 3:48:57 PM.

This article is licensed under a

Creative Commons Attribution 3.0 Unported Licence.

View Article Online

in terms of temperature and pressure. However, the critical LLE

line is significantly overestimated by the EOS due to a strong

sensitivity of the CEP and UCST on the binary interaction

parameter x.To address this issue, an empirical correlation

between the binary interaction parameter of the EOS x

EOS

and

the true xvalue was introduced. The CEP and CAEP were found

to have opposing trends, which exacerbates the challenge of

representing them concurrently with EOS, compared to the

already demanding task of characterizing several phase equili-

brium types simultaneously.

Unlike the vapor–liquid critical points, the transport and

thermophysical properties near the liquid–liquid critical point

were found to exhibit no significant anomalies or singularities,

such as local extrema or poles. Instead, all properties approach

the critical point either linearly with temperature or are indepen-

dent of it. The Fick diffusion coefficient is zero there only because

the thermodynamic factor vanishes. To address system-size effects

on diffusion coefficients, the simulation data were extrapolated to

the thermodynamic limit and analytical finite-size corrections

were applied. However, these analytical corrections showed an

overestimation, particularly when they incorporate the thermo-

dynamic factor near the liquid–liquid critical point.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors gratefully acknowledge the financial support by the

BMBF under the grant WindHPC. The authors I. A. and J. V.

gratefully acknowledge financial support by the Deutsche For-

schungsgemeinschaft (DFG) under grant no. VR 6/16. Equili-

brium molecular dynamics simulations were performed on the

Atos BullSequana XH2000 system (Noctua 2) at the Paderborn

Center for Parallel Computing (PC2), contributing to the project

MMHBF2. The authors S. St. and J. S. gratefully acknowledge

financial support by the DFG within IRTG 2057 ‘‘Physical

Modelling for Virtual Manufacturing Systems and Processes’’.

The EOS calculations were carried out on the ELWE computer

at Regional University Computing Center Kaiserslautern

(RHRK) under the grant TUKL-TLMV.

Notes and references

1 T. W. Leland, J. S. Rowlinson and G. A. Sather, Trans.

Faraday Soc., 1968, 64, 1447–1460.

2 S. Stephan and H. Hasse, Phys. Rev. E, 2020, 101, 012802.

3 J.-P. Hansen and I. McDonald, Theory of Simple Liquids,

Academic Press, Cambridge, 4th edn, 2013.

4 G. A. Mansoori, Fluid Phase Equilib., 1993, 87, 1–22.

5 J. K. Lee, J. A. Barker and G. M. Pound, J. Chem. Phys., 1974,

60, 1976–1980.

6 B. S. Carey, L. E. Scriven and H. T. Davis, AIChE J., 1980, 26,

705–711.

7 T. Wadewitz and J. Winkelmann, Ber. Bunsenges. Phys.

Chem., 1996, 100, 1825–1832.

8 D. Fertig, H. Hasse and S. Stephan, J. Mol. Liq., 2022,

367, 120401.

9 O. Loetgering-Lin, M. Fischer, M. Hopp and J. Gross, Ind.

Eng. Chem. Res., 2018, 57, 4095–4114.

10 K. C. Mo and K. E. Gubbins, Mol. Phys., 1976, 31, 825–847.

11 K. C. Mo, K. E. Gubbins, G. Jacucci and I. R. McDonald, Mol.

Phys., 1974, 27, 1173–1183.

12 R. Fingerhut, G. Herres and J. Vrabec, Mol. Phys., 2020,

118, e1643046.

13 D. Fertig and S. Stephan, Mol. Phys., 2023, e2162993.

14 F. Blas and I. Fujihara, Mol. Phys., 2002, 100, 2823–2838.

15 J. Vrabec, A. Lotfi and J. Fischer, Fluid Phase Equilib., 1995,

112, 173–197.

16 S. Stephan and H. Hasse, Mol. Phys., 2020, 118, e1699185.

17 J. J. Potoﬀ and A. Z. Panagiotopoulos, J. Chem. Phys., 1998,

109, 10914–10920.

18 A. Z. Panagiotopoulos, N. Quirke, M. Stapleton and

D. J. Tildesley, Mol. Phys., 1988, 63, 527–545.

19 S. P. Protsenko and V. G. Baidakov, Fluid Phase Equilib.,

2016, 429, 242–253.

20 S. P. Protsenko, V. G. Baidakov and V. M. Bryukhanov, Fluid

Phase Equilib.,2016, 430, 67–74.

21 J. B. Buhn, P. A. Bopp and M. J. Hampe, Fluid Phase Equilib.,

2004, 224, 221–230.

22 P. Geysermans, N. Elyeznasni and V. Russier, J. Chem. Phys.,

2005, 123, 204711.

23 M. Mecke, J. Winkelmann and J. Fischer, J. Chem. Phys.,

1999, 110, 1188–1194.

24 S. Stephan and H. Hasse, Phys. Chem. Chem. Phys., 2020, 22,

12544–12564.

25 R. Fingerhut and J. Vrabec, Fluid Phase Equilib., 2019, 485,

270–281.

26 J. Vrabec, A. Kumar and H. Hasse, Fluid Phase Equilib., 2007,

258, 34–40.

27 R. S. Chatwell, M. Heinen and J. Vrabec, Int. J. Heat Mass

Transfer, 2019, 132, 1296–1305.

28 S. Stephan, D. Schaefer, K. Langenbach and H. Hasse, Mol.

Phys., 2021, 119, e1810798.

29 V. G. Baidakov, S. P. Protsenko and V. M. Bryukhanov, Fluid

Phase Equilib., 2019, 481, 1–14.

30 J. Staubach and S. Stephan, J. Chem. Phys., 2022,

157, 124702.

31 J. M. Garrido, M. M. Pin

˜eiro, A. Mejı

´a and F. J. Blas, Phys.

Chem. Chem. Phys., 2016, 18, 1114–1124.

32 E. L. Granados-Baza

´n, S. E. Quin

˜ones-Cisneros and

U. K. Deiters, J. Chem. Phys., 2021, 154, 084704.

33 I. Nezbeda, J. Kolafa and W. R. Smith, J. Chem. Soc., Faraday

Trans., 1997, 93, 3073–3080.

34 F. J. Martinez-Ruiz, A. I. Moreno-Ventas Bravo and F. J. Blas,

J. Chem. Phys., 2015, 143, 104706.

35 M. Heier, S. Stephan, F. Diewald, R. Mu

¨ller, K. Langenbach

and H. Hasse, Langmuir, 2021, 37, 7405–7419.

36 M. P. Allen and D. J. Tildesley, Computer Simulation of

Liquids, Oxford University Press, Oxford, 1989.

PCCP Paper

Open Access Article. Published on 05 June 2023. Downloaded on 8/8/2023 3:48:57 PM.

This article is licensed under a

Creative Commons Attribution 3.0 Unported Licence.

View Article Online

37 D. A. Jonah and U. K. Deiters, Mol. Phys., 1992, 77, 1071–1083.

38 U. K. Deiters and T. Kraska, High-Pressure Fluid Phase

Equilibria - Phenomenology and Computation, Elsevier,

Amsterdam, 2012.

39 V. Harismiadis, N. Koutras, D. Tassios and

A. Panagiotopoulos, Fluid Phase Equilib., 1991, 65, 1–18.

40 A. M. Georgoulaki, L. V. Ntouros, D. P. Tassios and

A. Z. Panagiotopoulos, Fluid Phase Equilib., 1994, 100, 153–170.

41 J. Kolafa and I. Nezbeda, Fluid Phase Equilib., 1994, 100,

1–34.

42 J. Binney, The Theory of Critical Phenomena: An Introduction to

the Renormalization Group, Clarendon Press, Oxford, 1992.

43 D. I. Uzunov, Introduction to the theory of critical phenomena:

mean field, fluctuations and renormalization, World Scientific,

Singapore, 1993.

44 M. Barmatz, I. Hahn, J. A. Lipa and R. V. Duncan, Rev. Mod.

Phys., 2007, 79, 1–52.

45 M. M. Telo da Gama and R. Evans, Mol. Phys., 1983, 48,

229–250.

46 H. A. Lorentz, Ann. Phys., 1881, 248, 127–136.

47 D. Berthelot, C. R. Hebd. Seances Acad. Sci., 1898, 126,

1703–1706.

48 V. A. Mazur, L. Z. Boshkov and V. G. Murakhovsky, Phys.

Lett. A, 1984, 104, 415–418.

49 U. K. Deiters, High-pressure fluid phase equilibria: Phenom-

enology and computation, Elsevier, Amsterdam, 1st edn,

2012, vol. 2.

50 P. H. van Konynenburg and R. L. Scott, Philos. Trans. R. Soc.,

A, 1980, 298, 495–540.

51 J. G. Kirkwood and F. P. Buﬀ, J. Chem. Phys., 1951, 19,

774–777.

52 I. Antolovic

´and J. Vrabec, Ind. Eng. Chem. Res., 2022, 61,

3104–3112.

53 B. Smit, J. Chem. Phys., 1992, 96, 8639–8640.

54 J. Algaba, B. Mendiboure, P. Go

´mez-A

´lvarez and F. J. Blas,

RSC Adv.,2022, 12, 18821–18833.

55 D. M. Heyes, Comput. Methods Sci. Eng., 2015, 21, 169–179.

56 A. Lotfi, J. Vrabec and J. Fischer, Mol. Phys., 1992, 76,

1319–1333.

57 S. Stephan, M. Thol, J. Vrabec and H. Hasse, J. Chem. Inf.

Model., 2019, 59, 4248–4265.

58 S. Stephan, J. Staubach and H. Hasse, Fluid Phase Equilib.,

2020, 523, 112772.

59 S. Stephan and U. K. Deiters, Int. J. Thermophys., 2020,

41, 147.

60 J. K. Johnson, J. A. Zollweg and K. E. Gubbins, Mol. Phys.,

1993, 78, 591–618.

61 T. Lafitte, A. Apostolakou, C. Avendan

˜o, A. Galindo,

C. S. Adjiman, E. A. Mu

¨ller and G. Jackson, J. Chem. Phys.,

2013, 139, 154504.

62 M. Benedict, G. B. Webb and L. C. Rubin, J. Chem. Phys.,

1940, 8, 334–345.

63 R. T. Jacobsen and R. B. Stewart, J. Phys. Chem. Ref. Data,

1973, 2, 757–922.

64 J. A. Barker and D. Henderson, Annu. Rev. Phys. Chem., 1972,

23, 439–484.

65 J. A. Barker and D. Henderson, J. Chem. Phys., 1967, 47,

2856–2861.

66 T. W. Leland, P. S. Chappelear and B. W. Gamson, Trans.

Faraday Soc., 1962, 8, 482–489.

67 U. Deiters and G. M. Schneider, Ber. Bunsenges. Phys. Chem.,

1976, 80, 1316–1321.

68 R. Lustig, Mol. Phys., 2012, 110, 3041–3052.

69 P. Ganguly and N. F. A. van der Vegt, J. Chem. Theory

Comput., 2013, 9, 1347–1355.

70 I.-C. Yeh and G. Hummer, J. Phys. Chem. B, 2004, 108,

15873–15879.

71 S. H. Jamali, L. Wolﬀ, T. M. Becker, A. Bardow, T. J. H. Vlugt

and O. A. Moultos, J. Chem. Theory Comput., 2018, 14,

2667–2677.

72 G. Guevara-Carrion, R. Fingerhut and J. Vrabec, J. Phys.

Chem. B, 2020, 124, 4527–4535.

73 A. T. Celebi, S. H. Jamali, A. Bardow, T. J. H. Vlugt and

O. A. Moultos, Mol. Simul., 2021, 47, 831–845.

Paper PCCP

Open Access Article. Published on 05 June 2023. Downloaded on 8/8/2023 3:48:57 PM.

This article is licensed under a

Creative Commons Attribution 3.0 Unported Licence.

View Article Online