Phase equilibria of symmetric Lennard-Jones mixtures and a look at the transport properties near the upper critical solution temperature

Abstract

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 different 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 different 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 significant 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 development 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 different types of phase behavior. These mixtures are throughout symmetric, i.e. both components have identical size, energy and mass parameters. The five mixtures only differ in their cross-interaction parameter. Symmetric systems are a particularly 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), azeotropic 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 substantial 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 conditions hold at these points (see Section 2.2 for details). A binary mixture may exhibit a single CEP and a single CAEP, which 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 different 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 origin of this issue is not fully understood, 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 ξ_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 thermophysical 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.⁴⁴

Fig. 1 T–x diagram at constant pressure p = 0.0724 (a) and p–x diagram at constant temperature T = 0.85 (b) of mixture C (ξ = 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).

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 different 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 (σ₁ = σ₂), energy (ε₁ = ε₂) and mass (m₁ = m₂) parameters of their pure components, while variations are introduced by different cross-interactions. Interactions between unlike LJ particles were described with the modified Lorentz–Berthelot combining rules^46,47where ξ is the binary interaction parameter. Although the pure components are identical, non-ideal thermodynamic mixture behavior can be introduced through ξ ≠ 1. 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,49Table 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_B/ε, pressure p* = pσ³/ε or self-diffusion coefficient . 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 ξ, the LJ mixtures belong to varying types according to the classification of van Konynenburg and Scott⁵⁰ and have different 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 simulative means. The steep slope of the critical LLE line required horizontal cuts in the T–x diagram, while the VLLE and azeotropic lines were analyzed vertically in p–x diagrams as presented in Fig. 1.

Table 1 Mixtures and properties that were investigated in this work

Mixture	A	B	C	D	E
Binary parameter ξ	0.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
[scr A, script letter A]mn, cv, cp, …, Γ
D i , Đij, Dij
η, λ

---

Fig. 2 Characteristic lines and points of mixture C in the p–T projection. Horizontal (a) and vertical (b) cuts correspond to T–x and p–x diagrams 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).

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 mechanical equilibrium, the equality of chemical potentials is required. For LLE and VLE in a symmetric binary mixture, this requirement is expressed as follows:

μ^L1₁ = μ^L2₁ = μ^L1₂ = μ^L2₂, (3)

μ^L₁ = μ^V₁, μ^L₂ = μ^V₂. (4)

Note that the conditions for LLE (3) and VLE (4) differ. 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 μ^LLE. Similarly, under VLE (4), the chemical potentials can be summarized, but it has to be differentiated between the components μ^VLE₁ and μ^VLE₂. Moving on to the VLLE, the necessary conditions are

μ^L1₁ = μ^L2₁ = μ^V₁, (5)

μ^L1₂ = μ^L2₂ = μ^V₂. (6)

Since all chemical potentials in the two coexisting liquid phases are identical, eqn (5) and (6) can be combined to

μ^LLE = μ^V₁ = μ^V₂. (7)

Eqn (7) offers 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 μ^V₁ = μ^V₂ = μ^V(x₁ = 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 μ^V₁ = μ^VLE₁ and μ^V₂ = μ^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₁ = μ^VLE₂ = μ^V(x₁ = 0.5). (8)

In summary, due to the symmetry of the mixtures, four otherwise different 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 μ–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^σ of the VLLE. The numerical results are listed in the ESI.†

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 (×).

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 spinodal, which are characterized by a vanishing thermodynamic factor Γ = 0, meet. For a binary mixture, the thermodynamic factor is given by

In this work, the values for Γ were 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 thermodynamic factor was additionally obtained by utilizing the derivative 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_i = 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

Γ(T_UCST, x₁ = 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 x₁ = 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 thermodynamic 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 Γ–T diagram 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 Γ as 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 as a reference point. Once the VLLE line and the critical line of the LLE have been determined as discussed above, their stark slope difference 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 (ξ = 0.8) once its approximate location was determined. However, this trial and error approach was found to be ineffective 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 converge. In order to determine the critical temperature T_c, pressure p_c and density ρ_c at the CAEP, the scaling law for density and the rectilinear diameter were employed^53,54

ρ_L − ρ_V = A·(T_c − T)^β, (11)

where β = 0.326 is the Ising critical exponent⁵⁵ and ρ_L, ρ_V are the saturated liquid and vapor densities at a given temperature T. The parameters A, B and critical values T_c, ρ_c 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.

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
ξ	0.75	0.80	0.85	1.20
T CEP	1.000	0.904	0.776	—
p CEP	0.049	0.022	0.006	—
T CAEP	1.153	1.181	1.211	1.439
p CAEP	0.111	0.112	0.114	0.135

---

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 reference 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 + a_pert, (13)

which depend on temperature T and density ρ. BH-type EOS uses the ideal gas contribution a_id and a hard sphere contribution a_hs as the reference

a_ref = a_id + a_hs. (14)

The dispersive interactions are captured within the perturbation contribution a_pert and are calculated by BH theory^64,65 in 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–T diagrams 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 simulation with a focus on transport 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 comprised 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 consists 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 thermodynamic factor route. The azeotropic line consists of individual azeotropic points as well. The critical points CEP and CAEP require all of the above. Therefore, the presentation of these characteristic lines and points in a p–T projection 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 differentiation 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 interaction parameter ξ causes the characteristic lines and points to shift along the pressure and temperature axes, either above or below the pure component vapor pressure curve (for ξ < 1 or ξ > 1, respectively). For ξ = 1, all lines and points coincide with the pure component vapor pressure curve. The CEP and CAEP exhibit opposing trends, approaching each other as ξ decreases. 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 ξ were 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.

Fig. 4 p–T projection 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 (+).

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 ξ_EOS was introduced into the EOS that differs from the force field parameter ξ. 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 × 4 grid is sufficient 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 ξ. This is demonstrated by the significant difference in the rate with which it propagates through the p–T plane compared to the VLLE and azeotropic lines.

The results presented in Fig. 5 demonstrate that the CEP and CAEP exhibit opposing trends. As ξ increases, the temperature 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 ξ < 0.70. Extrapolating the simulation data leads to an estimate of ξ ≈ 0.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 investigated mixtures are listed in Table 2.

Fig. 5 Temperature of the critical points as a function of the binary interaction parameter ξ. 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).

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 computer 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 and extrapolates well to extreme temperature and pressure conditions 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 ξ. Fig. 6 shows the adjusted binary interaction parameter ξ_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 T ≤ 0.92. For higher temperatures, only VLE data were used for the fit. For each temperature studied by molecular simulation, the ξ_EOS value was fitted individually. Thereby, the temperature dependence of ξ_EOS(T) was obtained, which can be taken as a measure for the deficiencies of the mixture model. The values of ξ_EOS close to the true value ξ = 0.8 indicate a good predictive capability of the mixture model. All studied EOS show the same qualitative behavior of ξ_EOS: in the low-temperature regime, it increases linearly with rising temperature and the slope is approximately the same. The EOS of Johnson et al.⁶⁰ and Kolafa and Nezbeda⁴¹ require very similar ξ_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 EOS of Stephan et al.⁵⁸ exhibits the largest ξ_EOS values.

Fig. 6 Binary interaction parameter of the EOS ξ_EOS for mixture C (ξ = 0.8) adjusted to simulation data as a function of temperature. VLE and LLE data were included for T ≤ 0.92 and VLE data only for T > 0.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.

In the temperature range where only VLE data were considered in the fitting procedure (T > 0.92), only a minor temperature dependence of ξ_EOS was observed. For the EOS of Stephan et al.,⁵⁸ there is hardly any temperature dependence of ξ_EOS. The EOS of Lafitte et al.⁶¹ exhibits the largest ξ_EOS values in the high temperature regime, whereas the smallest ξ_EOS values were obtained for the Kolafa and Nezbeda⁴¹ EOS.

Important differences were observed for the ξ_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, ξ_EOS values obtained from the fit to the simulation data exhibit a discontinuity in the vicinity of the CEP. This is the case for all four 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⁴¹ was used as its binary interaction parameter ξ_EOS required the smallest adjustments.

Fig. 7 shows the results for the ξ_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 ξ_EOS is observed. The binary interaction parameter ξ_EOS for the type II mixtures, on the other hand, was fitted to VLE and LLE data, but shows no temperature dependence. The average offset of ξ_EOS from ξ is 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.

Fig. 7 Binary interaction parameter ξ_EOS for mixtures A to E adjusted to simulation data as a function of temperature. Results were obtained by fitting ξ_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.

3.3 Equilibrium properties near the UCST

The attention is now turned to the Helmholtz energy derivatives, which have been obtained with the formalism of Lustig,⁶⁸ and the thermodynamic factor, which serves for conversion between transport diffusion coefficients. 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 (Γ < 0) and is now approached from the homogenous region (Γ > 0), yielding a value of T_UCST = 0.9084 in both cases.

The partial derivatives of the reduced Helmholtz energy, denoted as , can be used to determine all thermodynamic equilibrium properties, such as the isochoric heat capacity c_v, speed of sound w or the Joule–Thomson coefficient μ_JT. For additional information on the subject, the reader is referred to the ESI.†Fig. 8 shows the results for the residual contributions to the reduced Helmholtz energy derivatives. The respective ideal gas contributions for the same states are given in the ESI.† EOS and simulation results indicate that the derivatives 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 , than for the ones associated with thermal properties , . The application of the adjusted binary interaction parameter ξ_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. 8 Residual Helmholtz energy derivatives (top) and thermophysical properties α, β, γ, c_v, c_p, w, Γ_G, μ_JT (bottom) near the UCST of mixture C at p = 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.

At the critical point, the thermodynamic factor vanishes (Γ = 0), which has specific implications for the Maxwell-Stefan (MS) and Fick diffusion coefficients and their finite-size corrections. Fig. 9 shows the results for the thermodynamic factor obtained from simulation and EOS. It is apparent that both approaches predict a 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 different manner. 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 calculated 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 on the unique advantages offered by each method.

Fig. 9 Thermodynamic factor around the UCST (green circle) of mixture C at p = 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 (●) for three system sizes. Results for N = 20 000 based on RDF vdV + shf and extrapolated (white circles) were extended further into the homogenous region.

Fig. 1 and 9 show the deviation between simulation data and EOS in predictive mode, which is off by approximately ΔT = 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 Γ and 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 Γ, 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 Γ values are relatively close to zero, but the difference among them can be significant. Variations due to different 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 Γ ≈ 0.06, it can be three to four times larger for smaller systems, cf.Table 3. This is particularly relevant for transport properties, as Γ is typically used as a factor or divisor.

Table 3 Values for the thermodynamic factor of the equimolar mixture C at T = 0.91 and p = 0.0724 obtained from KBI

	RDF
N	Standard	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 → ∞
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)

---

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 dependence, 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.

Fig. 10 Simulation results (left) and finite size correction (right) for the self- (top), MS (center) and Fick (bottom) diffusion coefficients 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 diffusion coefficient data from extrapolated Onsager coefficients (green cross).

Unlike the self-diffusion coefficient, the results for the MS diffusion coefficient 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, effectively revealing a clear dependence on system size, cf.Fig. 10 (center right).

Since the Fick and MS diffusion coefficients are related by the thermodynamic factor according to D_ij = Γ· Đ_ij, the Fick diffusion coefficient traces the same linear trajectory as the thermodynamic factor (cf.Fig. 9) multiplied by the constant MS diffusion coefficient. Hence, the Fick diffusion coefficient is zero at the UCST solely due to the thermodynamic factor. The resulting set of curves highlights the system size dependence through different slopes rather than displacements as observed for the self-diffusion coefficient, cf.Fig. 10 (bottom left).

When examining the impact of system size in Fig. 10 (top right), the self-diffusion coefficient reveals a linear relationship with the inverse of the edge length of the simulation volume L⁻¹. By fitting a linear function to the simulation data and extrapolating to the thermodynamic limit L⁻¹ → 0, the value for D^∞_i coincides with the intercept. Yeh and Hummer⁷⁰ investigated the system size dependence of the self-diffusion coefficient and proposed a correction term to account for finite-size effects far from the thermodynamic limit. Their correction term scales linearly with the inverse of the edge length of the simulation volume L⁻¹

where k_B is the Boltzmann constant, η the shear viscosity and ζ ≈ 2.837297 is a dimensionless constant for cubic simulation volumes.⁷⁰

After applying the correction to the self-diffusion coefficient, the average value of the corrected data demonstrates an overestimation in comparison to the value of D^∞_i obtained from extrapolated simulation data. However, a more accurate comparison can be made by considering the slopes of the curves instead of their intercepts. The correction term D^YH aims to equalize the slope of the self-diffusion coefficient so that only the intercept D^∞_i remains. To obtain a horizontal fit curve for the results in this work, D^YH must be scaled down by approximately 22%.

Similar to the self-diffusion coefficient, a finite-size correction for the MS diffusion coefficient was proposed by Jamali et al.⁷¹

Đ^∞_ij = Đ_ij + D^YH/Γ. (16)

As shown in Fig. 10 (center right), the four system sizes scale linearly with L⁻¹ and the data have a modest slope, allowing for an extrapolation to the thermodynamic limit L⁻¹ → 0. However, the correction (16) proposed by Jamali et al.⁷¹ suggests a qualitatively different behavior. By correcting the MS diffusion coefficient 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 Γ is not advisable. These findings highlight the need for a more comprehensive approach to correcting the MS diffusion coefficient near the critical point. Further research into the system size effects of the MS diffusion coefficient in this region is worthwhile.

In addition to the correction proposed by Jamali et al.,⁷¹ the finite-size corrected MS diffusion coefficient can also be obtained by extrapolating the Onsager coefficients to the thermodynamic limit.⁷² Recently, Celebi et al.⁷³ published a comprehensive overview of finite-size corrections for diffusion coefficients.

The system size correction for the Fick diffusion coefficient can be calculated by applying D_ij = Γ·Đ_ij to eqn (16)

D^∞_ij = D_ij + D^YH, (17)

where the thermodynamic factor cancels out and the correction term D^YH for the self-diffusion coefficient remains. At the UCST, the first term in eqn (17)D_ij becomes zero. Since the Yeh and Hummer⁷⁰ correction D^YH has the form of a linear function without intercept, it follows that the Fick diffusion coefficient in the thermodynamic limit D^∞_ij also becomes zero, cf.Fig. 10 (bottom right). It can be concluded that all three diffusion coefficients 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 parameter ζ ≈ 2.837297 of the correction D^YH would need to be adjusted to 2.2 ± 0.3 in order to obtain a horizontal correction as aspired by Yeh and Hummer.⁷⁰

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 different system sizes (circles). Constant fit (dashed lines).

In this limited range of states, the thermal conductivity also remains independent of temperature, cf.Fig. 11 (bottom). Moreover, there is no observable effect 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 different 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 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 ξ. To address this issue, an empirical correlation between the binary interaction parameter of the EOS ξ_EOS and the true ξ value 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 equilibrium 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 independent 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 thermodynamic 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 Forschungsgemeinschaft (DFG) under grant no. VR 6/16. Equilibrium 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. Potoff 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. Piñeiro, A. Mejía and F. J. Blas, Phys. Chem. Chem. Phys., 2016, 18, 1114–1124.
32. E. L. Granados-Bazán, S. E. Quiñ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. Mü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.
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: Phenomenology 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. Buff, J. Chem. Phys., 1951, 19, 774–777.
52. I. Antolović 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. Gómez-Á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. Avendaño, A. Galindo, C. S. Adjiman, E. A. Mü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. Wolff, 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.

---

Footnote

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

---