From molecular to macroscopic: predicting liquid–liquid phase equilibria and small-angle scattering of mixtures of organic liquids from atomistic simulation using Kirkwood–Buff theory

Abstract

Macroscopic phase equilibria between solutions define the functionality of many biological and industrial processes, yet they are challenging to predict due to the inherent complexity of liquids containing large molecules. This work introduces an approach for the purely predictive calculation of such phase equilibria in temperature-composition space from molecular dynamics (MD) simulations at one temperature in the single-phase region. We use an approach developed previously to obtain the entropic and enthalpic contributions to the free energy of mixing from the atomic-scale information given by MD simulations via Kirkwood–Buff theory. This allows us to accurately estimate the free energy of mixing as a function of temperature, and thus obtain liquid–liquid phase equilibria, including liquid–liquid critical points, associated binodal and spinodal lines, and composition fluctuations across a region of temperature and composition. Results for binary malonamide–alkane systems are validated by comparison to a direct experimental probe of the fluctuations: the small angle X-ray scattering intensity near zero wavenumber. The MDKB → Phase method demonstrated here provides a significant improvement in predicting liquid–liquid equilibria and free energy as a function of temperature for our systems of interest compared to conventional thermodynamic models. The accurate performance of this purely predictive approach lies in its preservation of atomistic details when determining thermodynamic properties. Furthermore, its inherent extensibility to multi-component systems will likely make the MDKB → Phase approach a valuable general tool for connecting molecular interactions to macroscopic phase equilibria and for the computational screening of materials for targeted thermodynamic behavior.

---

Introduction

Liquid phase equilibria are fundamental to many biological^1,2 and industrial chemical processes.³ Prediction of such phase equilibria is generally performed using models of solution thermodynamics which treat mixtures using simplified representations of their molecular constituents and their interactions.⁴ Meanwhile, molecular simulation is the predominant computational tool for investigating liquids at the molecular level, such as understanding their molecular structure or the atomistic origins of physicochemical properties of liquid mixtures (e.g., density, viscosity).⁵ A key advantage of molecular dynamics (MD) simulations is their generation of entire ensembles of condensed-phase structures with atomistic resolution to all intra- and inter-molecular interactions, thus avoiding the limitations of thermodynamic models that employ simplified interaction parameters that by necessity represent averages over entire molecular species and their configurations in solution.^6,7 However, it is impractical to directly simulate macroscopic phase separation using MD, because of the limited time and length scales possible, and alternative methods are prohibitively expensive for many chemical systems. Here, we seek to bridge the divide between molecular simulation and traditional thermodynamic models of macroscopic phase equilibria, using MD to predict phase equilibria in a manner that retains atomistic details of intramolecular (conformational) and intermolecular (explicit mutual solvation) configurations.

Kirkwood–Buff (KB) theory⁸ provides a means to connect the information encoded in molecular simulation trajectories to the thermodynamic activities of their components.⁹ Combining atomistic simulations, which naturally account for conformational degrees of freedom and all intermolecular interactions explicitly, with this exact statistical mechanical theory provides a powerful framework for deriving macroscopic thermodynamic properties directly from underlying molecular behavior.^6,8,10 Pioneering work by Petris et al.¹¹ showed how thermodynamic quantities such as the free energy of mixing could be obtained from MD simulations using KB theory. Initial studies of solutions of small molecules¹¹ were subsequently extended to large molecules.¹² Crucially, by using the enthalpy of mixing obtained from the MD simulations, the entropy of mixing could be separately determined. We have used this feature to extend the predictions from MD simulations at a single temperature to a region of temperature-composition space, allowing calculation of phase diagrams from single-phase simulations, avoiding computational expenses associated with direct simulation of coexisting phases,¹³ or needing to compute solvation free energies of all species at all compositions with alchemical free energy methods.¹⁴ This approach can be paired with direct experimental measurement of composition fluctuations using small-angle X-ray scattering (SAXS) intensity as wavenumber q → 0, thereby offering a robust means to compare experimental results with theoretical predictions at any point in the single-phase region of the phase diagram. Thus, SAXS provides a more convenient alternative to probe solution thermodynamics than direct observation of phase transitions. In particular, it is useful for systems whose liquid–liquid critical point is inaccessible, e.g. due to freezing. This method therefore enables validation of the results obtained from the KB theory combined with MD simulations, even when direct phase diagram measurements are not feasible.

The integration of KB theory with MD simulations to predict phase equilibria represents a transformative step, offering a direct pathway to predictive capabilities based on detailed atomistic structures. The workflow of this methodology is illustrated in Fig. 1, with the details of each step summarized below. This approach significantly advances the ability to predict phase diagrams and related thermodynamic properties directly from molecular correlations, marking a shift from merely descriptive applications of MD to a robust, predictive role in determining phase behavior.⁵ By grounding thermodynamic predictions in atomistic insights, this approach enhances the accuracy and relevance of simulations in practical applications. Specifically, we demonstrate the combined MDKB → Phase approach to the study of binary organic phase systems, which are central to the solvent extraction of rare earth elements – a process of significant industrial importance.¹⁵ We compare this approach to a traditional thermodynamic method that fits scalar interaction parameters (UNIQUAC), as well as to predictive approaches (UNIFAC and COSMO-RS). The results from these comparisons underscore the power of the MDKB → Phase approach, demonstrating its ability to provide more accurate and reliable predictions than those obtained through conventional empirical methods. This not only validates the theoretical framework but also showcases the practical benefits of integrating detailed atomistic simulations with advanced thermodynamic theories.

Fig. 1 Workflow for MDKB → Phase approach. Six steps consist of: (1) run MD simulations as a function of composition at a particular temperature T_sim; (2) calculate RDFs from each simulation between each pair of species; (3) calculate and correct Kirkwood–Buff integrals and combine to obtain free energy curvature; (4) fit excess free energy curvature and enthalpy to get expressions for excess properties; (5) extrapolate in temperature and generate phase diagram; and (6) compare calculated and measured X-ray scattering at an experimentally convenient (i.e., room) temperature.

Methodology

Systems studied

We selected six binary (extractant/solvent) systems relevant to liquid–liquid extraction of rare earth elements, an important separations process for those technologically critical materials. These systems include two malonamide (MA) extractants, N,N′-dimethyl-N,N′-dibutylpentyl malonamide (DBP) and N,N′-dimethyl-N,N′-dioctylhexylethoxy malonamide (DOHE), as shown in Fig. 1, paired with various n-alkane solvents. Specifically, the systems studied are DBP/n-hexane, DBP/n-octane, DBP/n-decane, DBP/n-dodecane, DOHE/n-heptane, and DOHE/n-dodecane. For each system, MD simulations were conducted at every 0.1 volume fraction of extractant at T_sim = 400 K, resulting in nine MD samples for each system along with simulations for each pure component. Additional details are provided in the SI.

Molecular dynamics simulations were conducted using the GROMACS software package,¹⁶ employing the GAFF2 force field for all molecules,¹⁷ which we have previously found accurately capture behavior in these organic mixtures.^18–21 Random initial configurations were created using Packmol, and the samples were energy minimized using the steepest descent algorithm.²² The simulations utilized a leap-frog Verlet integrator with a 2 fs timestep²³ with hydrogen-containing bonds constrained using the LINCS algorithm. Long-range electrostatic interactions were calculated using Particle-Mesh Ewald summation, with a 1.5 nm cutoff for both short-range electrostatics and Lennard-Jones interactions.²⁴ Each sample was equilibrated in the isobaric–isothermal (NpT) ensemble for 5 ns at 1 atm and 400 K, using the Berendsen barostat and v-rescale thermostat.^25,26 Subsequently, production runs were carried out in the NpT ensemble for 600 or 1100 ns for DBP and DOHE systems, respectively at the same temperature and pressure, employing the c-rescale barostat and Nosé–Hoover thermostat.^27,28

Kirkwood–Buff integrals

KB theory provides a robust statistical mechanical framework that bridges microscopic molecular correlations to macroscopic thermodynamic properties. This connection is facilitated through the Kirkwood–Buff integrals (KBIs).^8,29 KBIs are defined in the open grand canonical ensemble (µVT), evaluating the excess or deficiency of molecule j around molecule i.^7,8,30–32 They can be expressed aswhere g^µVT_ij(r) is the pair correlation function for species i and j.

While KBIs are well-defined, in practice their calculation from MD simulations are typically performed in the NpT ensemble, with constant particle number.⁶ Additionally, for finite samples, the application of KB theory requires that the simulation box size be sufficiently large compared to the solution's correlation length scale.^32,33 Within the NpT ensemble, a version of eqn (1) is often applied, where the integral is truncated at a distance R, up to 1/2 the simulation box linear dimension, to yield the “running KBI,”^32,34,35

At short distances, G^R_ij suffers from convergence issues reflected in the running KBI, specifically as r → R, as a result of the asymptotic behavior in the NpT ensemble, diverging as R³.^32,34–36 Capturing the tail behavior is vital for accurate KBI analysis as long-range fluctuations and seemingly small differences in the tail of the RDF can significantly impact the structural and thermodynamic properties of a sample.³⁷ To correct KBIs from finite systems to the thermodynamic limit, the following two corrections are implemented (details provided in the SI). We apply the KBI correction demonstrated by Milzetti et al. to correct convergence of g_ij(r).³⁵ Then, we use the approach detailed by Simon et al. to extrapolate the KBIs to the thermodynamic limit, where their relationships to macroscopic thermodynamic properties are defined.³⁴

Solution thermodynamics from Kirkwood–Buff integrals

The composition and temperature dependence of the Gibbs free energy of mixing ΔG_mix defines the phase diagram.³⁸ This property is not straightforward to directly calculate, and, experimentally, is typically obtained through indirect measurements, i.e., by its impact on other measurable quantities. KB theory⁸ relates preferential solvation in the grand canonical ensemble to the curvature(s) of ΔG_mix with respect to composition. For a binary mixture, there is only one independent mole fraction x, taken to be that of component 1. The curvature of the mixing free energy is given by^8,11,12,39where RT is the molar Boltzmann constant times the temperature, and ρ is the total molar density of all components. The excess free energy curvature iswhere G^E is the excess free energy and G^id ≡ RT[x log x + (1 − x)log(1 − x)] is the ideal contribution.

Excess thermodynamic quantities have often been modeled using polynomials of the composition. Mukhopadhyay et al. demonstrated the applicability of such polynomials to various excess thermodynamic properties, such as G^E, ΔH_mix and S^E.⁴⁰ While initial research on these polynomial models focused on two or three component mixtures, Tomiska and Neckel showed how a fourth order (quartic) polynomial can describe excess thermodynamic properties in multi-component mixtures. We implement such a quartic polynomial to describe excess thermodynamic quantities in our binary system, namely⁴¹

A^E = a₁₁(x² − x) + a₁₁₁(x³ − x) + a₁₁₁₁(x⁴ − x), (5)

where A^E represents any excess thermodynamic property, including ΔH_mix, S^E, and G^E, which are zero for pure systems (x = 0 and x = 1). The expression for the second derivative is

We can thus fit eqn (6) to the second derivative of the excess free energy curvature obtained from the KBIs at various x using eqn (3) and (4) to obtain the three parameters a₁₁, a₁₁₁, and a₁₁₁₁. Using these parameters in eqn (5) gives an expression for the excess free energy of mixing G^E as a function of composition.

Separating enthalpic and entropic contributions

In order to extrapolate ΔG_mix from the simulation temperature T_sim to other temperatures, as described below, it is necessary to separate the enthalpic and entropic contributions to G^E. The mixing enthalpy can be obtained directly from simulation without KB theory using

ΔH_mix = 〈H〉 − [xH^pc₁ + (1 − x)H^pc₂], (7)

where H^pc_i denotes the pure component enthalpy of species i, and 〈H〉 is the observed enthalpy of the mixture.⁴² This allows us to obtain the excess entropy using

We can fit simulation values of ΔH_mix for various x with the quartic function of eqn (5). Subtracting the quartic functions for ΔH_mix and G^E thus gives S^E as a function of composition according to eqn (8).

Extrapolation in temperature

The methods described above give expressions for G^E, ΔH_mix, and S^E as a function of composition at T = T_sim. To obtain phase equilibria and scattering intensity at other temperatures from ΔG_mix, we can assume that ΔH_mix and S^E are independent of temperature, and usewhere in the last line we have substituted in eqn (8) for S^E. While the total enthalpy H is expected to change with temperature, defined by the constant pressure heat capacity , ΔH_mix and S^E (and by extension ΔS_mix, with ΔS_mix = −R[x ln x + (1 − x)ln(1 − x)] + S^E) do not necessarily vary with temperature.⁴³ The temperature independence of ΔH_mix and ΔS_mix is a standard assumption in thermodynamic models. In the SI, we demonstrate the temperature independence of ΔH_mix and ΔS_mix using MD as a function of composition and temperature for a typical system considered here.

Structure factors and scattering

Equilibrium small-angle scattering is quantitatively related to composition fluctuations and thus to the curvature of the free energy of mixing. The scattering intensity contains terms proportional to composition fluctuations as well as vibrations. A simplification occurs because the term arising from vibrations is typically negligible, which we demonstrate for our systems in the SI. As shown there, the X-ray scattering cross section per unit volume as q → 0 is given by⁴⁴where r_e² = 7.95 × 10⁻²⁶ cm² per electron is the electron scattering cross section, N_A = 6.02 × 10²³ electrons per mole is Avogadro's number, Z_i is the number of electrons per molecule, V_i is the partial molar volume of component i, and V̄ = xV₁ + (1 − x)V₂ = 1/ρ is the molar volume of the solution. Thus I₀ diverges as ∂²ΔG_mix/∂x² → 0 at the phase instability (spinodal). This is readily calculated at any temperature using the quartic equation obtained for ΔG_mix(x, T) as described above.

Results and discussion

This section follows the workflow presented in Fig. 1, taking the results of simulations of mixtures of MA and n-alkanes, starting with the corrected KBIs. Subsequently, expressions for excess properties were determined by fitting the excess curvature from KBIs and the mixing enthalpy from MD simulations. The entropic and enthalpic contributions to ΔG_mix are then obtained. We then compare these MDKB → Phase results to various thermodynamic models. Lastly, we use the MDKB → Phase results to compute liquid–liquid phase diagrams and compare the results to experimental SAXS values of I₀ for each system to validate the approach.

Kirkwood–Buff integrals

The running KBIs for each sample were corrected according to approaches outlined by Milzetti et al. and Simon et al., as described in the SI. Detailed individual fits for each sample, illustrating extrapolation of running KBIs to the thermodynamic limit, are included in the SI. The final extrapolated KBI values, denoted as G^∞_ij, are given in Fig. 2, with error bars calculated as described in the SI. The errors on the individual KBIs are largest at the composition extremes for self–self pairs where there are the small numbers of that molecule in the simulation, suggesting counting statistics of sampling the pairwise correlations dominates the observed errors. In the SI, we validate the corrected KBIs by cross-referencing the molar volumes and compressibilities calculable from the KBIs themselves^8,45 to values obtained directly from the average and variation of the MD simulation box volumes, finding excellent agreement.

Fig. 2 Composition dependence of each G^∞_ij is given, with each system in a different panel (A–F), for extractant–extractant, extractant–solvent, and solvent–solvent correlations, with (A) DBP/n-dodecane, (B) DBP/n-decane, (C) DBP/n-octane, (D) DBP/n-hexane, (E) DOHE/n-dodecane, and (F) DOHE/n-heptane.

Excess thermodynamic properties

As outlined in the Methodology section, the expression for G^E at the simulation T was determined by fitting a second degree polynomial to the weighted excess curvature points determined from the KBIs. Fitting a parabola, rather than a higher-order polynomial, has the advantage of avoiding complex shapes which can lead to artifacts in constructing the phase diagram. The error on the excess curvature can be determined from the fit as described in the SI. These points and fits are shown in Fig. 3, with the shaded region representing the error in ∂²G^E/∂x² determined from the polynomial fit. While the parabola does not always provide a good fit to the full range of data, the errors in the fit are smallest at intermediate compositions which are most important in determining the phase diagram.

Fig. 3 Points show the excess curvatures determined from the KBI for all systems and samples. Curves show the results of fitting to a second-order polynomial in x. Each system is displayed in a different panel (A–F), with (A) DBP/n-dodecane, (B) DBP/n-decane, (C) DBP/n-octane, (D) DBP/n-hexane, (E) DOHE/n-dodecane, and (F) DOHE/n-heptane.

By decoupling the enthalpic and entropic contributions to G^E, it becomes possible to estimate ΔG_mix as a function of temperature, which is required for predicting phase diagrams. To validate this approach, the temperature independence of ΔH_mix and ΔS_mix was first confirmed using the DBP/n-dodecane system as a case study, with simulations conducted at 380, 400, and 410 K (details provided in SI). The lack of temperature dependence in ΔH_mix and ΔS_mix allowed for reliable extrapolation of ΔG_mix to various temperatures.

An example of the ΔH_mix values obtained for a typical system is shown in Fig. 4A, along with the fit of these values to the polynomial given in eqn (5). Calculated curves of the entropy and free energy are also shown. Enthalpy values, fits, and calculated curves for all systems are given in the SI. The resulting analytical expression for ΔG_mix(x, T) greatly facilitates calculation of phase diagrams and scattering.

Fig. 4 Comparison of enthalpic and entropic contributions to the mixing free energy from various thermodynamic models for DBP in n-dodecane. (A) Mixing enthalpy, where the scatter points are results from MD simulation data, (B) mixing entropy, and (C) mixing free energy. Results from polynomial fits from the MDKB → Phase workflow shown in solid black, UNIFAC in dashed green, UNIQUAC in dotted blue, and COSMO-RS in dotted dash pink.

Liquid–liquid phase equilibria

While excess thermodynamic properties reveal the energetics of mixing, their extrapolation in temperature allows us to predict phase transitions directly from the behavior in the single-phase region. To predict the phase diagrams of each liquid mixture, ΔG_mix was estimated as a function of composition and temperature, as depicted for one system (DBP/n-dodecane) in Fig. 5. Fig. 5A displays curves of ΔG_mix(x) at various temperatures, both above and below the critical temperature (T_c). These curves demonstrate the changing curvature of ΔG_mix as the temperature is lowered below T_c, indicative of systems that exhibit an upper critical solution temperature (UCST). Open circles on these curves mark spinodal points, which are the inflection points, i.e. where the curvature equals zero. Closed circles indicate binodal points, which are determined using the common tangent construction.³⁸ At the binodal points, the chemical potentials of each component are equal in each phase.

Fig. 5 Results for phase equilibria using the MDKB → Phase approach for DBP in n-dodecane. (A) Mixing free energy estimated for select temperatures with binodal and spinodal points in black closed and open circles, respectively. (B) Heatmap of mixing free energy. (C) Heatmap of X-ray scattering as q → 0. Panels B and C have the solid binodal and dashed spinodal curves overlayed, with the critical point shown.

Once the spinodals and binodals are identified as functions of temperature, phase diagrams can be constructed. These phase diagrams are illustrated in Fig. 5B, where they are overlaid on a heatmap of ΔG_mix. The heatmap offers a visual representation of the thermodynamic landscape, highlighting regions of stability and instability. Subsequently, T_c was determined as the highest temperature at which spinodal points are observed. Note that only liquid–liquid equilibria are shown, not freezing or vaporization. For the DBP/n-dodecane system, where the critical temperature is experimentally accessible and measured (using the same experimental materials and procedure described in ref. 21), the MDKB → Phase approach exhibits good agreement with the experimental T_c, shown in Table 1.

Table 1 Comparison of T_c (K) of various thermodynamic models to experiment for DBP/n-dodecane, and their corresponding percent error from experiment. Error in the last digit for experimental measurement and the MDKB → Phase model is represented in parenthesis

Model	Tc (K)
Experiment	287(5)
MDKB → Phase	288(1)
UNIQUAC	400
UNIFAC	261
COSMO-RS	235

---

The remaining systems were similarly evaluated with the MDKB → Phase approach, with T_c and critical compositions summarized in SI. These results underscore the significant impact of extractant structure on the critical composition and phase behavior. Notably, the larger extractant, DOHE, forms mixtures that are more resistant to liquid–liquid demixing compared to the smaller DBP extractant when paired with the same diluent, as expected. Furthermore, this study confirms that increasing the length of the n-alkane diluent leads to an increase in both T_c and the critical composition of the extractant, consistent with previous findings reported for phosphate extractants.^19,46 These observations highlight the crucial role of molecular structure and diluent composition in determining phase separation behavior in liquid–liquid extraction systems.

Free energy curvature and small angle scattering

To further validate our thermodynamic modeling of mixing energetics, we compare our temperature-extrapolated mixing free energy curvatures to experimental scattering data. The prediction of X-ray scattering intensity from the MDKB → Phase approach facilitates a direct comparison between theory and experiment, establishing a robust framework for validating the accuracy of the model. This comparison of thermodynamic predictions using the experimentally observable composition fluctuations allows us to investigate all systems, even those whose liquid–liquid phase transitions are not experimentally accessible. Fig. 5C displays a heatmap of I₀, calculated using eqn (10), with the corresponding phase diagram overlaid. As the temperature approaches the spinodal temperature (T → T_sp), the phase instability becomes evident, characterized by I₀ → ∞, signaling the divergence of concentration fluctuations. Experimentally, observations are confined to conditions above the binodal curve, where the system remains stable. There, the largest fluctuations are observed near the critical point, where the binodal and spinodal curves intersect. To validate the MDKB → Phase approach, the values of I₀ calculated for 295 K were compared to experimentally determined I₀ values from SAXS measurements previously reported by us,^19,20 as shown in Fig. 6. Across the six systems studied, the results from the MDKB → Phase approach exhibit overall very good agreement with experimental data, capturing both the shape and the absolute intensity of I₀ particularly well for the DBP systems. This consistency not only confirms the reliability of the MDKB → Phase approach but also underscores its utility in analyzing the thermodynamics of demixing in complex liquid systems.

Fig. 6 Each system is represented in a different subplot (A–F), where the left panel shows the phase diagrams determined from the MDKB → Phase approach with the gray line indicating MD simulation temperature of 400 K and the colored line at 295 K is the temperature used for SAXS comparison as q → 0, which is shown on the right with the scatter points representing experimental data at room temperature. The shaded region on the right panel represents the error in I₀ calculated from MDKB → Phase approach. Subplots correspond to the following systems, (A) DBP/n-dodecane, (B) DBP/n-decane, (C) DBP/n-octane, (D) DBP/n-hexane, (E) DOHE/n-dodecane, and (F) DOHE/n-heptane.

Comparison to thermodynamic models

To compare the results from the MDKB → Phase method against established thermodynamic models, we calculate the predicted values of T_c for UNIQUAC, UNIFAC, and COSMO-RS.^47–49 Models based on binary interaction parameters have often been used to estimate various thermodynamic properties as a function of temperature and composition.⁴ UNIQUAC, a local-composition-based model, incorporates interaction energies into ΔH_mix and calculates S^E based on volume and area fractions.⁴⁷ The method we use to fit UNIQUAC interaction parameters to ΔH_mix from MD is explained in the SI. Similarly, UNIFAC uses a group contribution method to estimate thermodynamic properties.⁴⁸ The detailed equations and formulations for these models, including their derivations, are well documented in the literature and are widely recognized in the field, with relevant equations provided in the SI.^47,48,50,51 Using a different approach than UNIQUAC or UNIFAC, COSMO-RS, implemented using openCOSMO-RS,⁵² provides a procedure to predict excess thermodynamic functions based on description of molecules derived from electronic structure calculations.

For the DBP/n-dodecane system considered here, the MDKB → Phase approach most accurately predicts the experimental T_c, as shown in Table 1. The decomposition of the contributions to the excess free energy also varies significantly between thermodynamic models and our MDKB → Phase results. Fig. 4 compares the MDKB → Phase results for ΔH_mix, ΔS_mix, and ΔG_mix to those obtained from the UNIQUAC, UNIFAC, and COSMO-RS models, for a typical system. Comparisons for all systems are given in the SI. Both predictive models, UNIFAC and COSMO-RS, underpredict the contributions of ΔH_mix and ΔS_mix compared with the MDKB → Phase results. Similarly, UNIQUAC, which uses the same entropic term as UNIFAC, also underpredicted the ΔS_mix term, with its interaction parameters fit to ΔH_mix.

Scalar interaction parameters in established thermodynamic models are often effective but can encounter challenges, particularly in systems with diverse molecular conformations where capturing the full complexity of intermolecular interactions is crucial. The binary interaction parameters are typically determined by fitting expressions from thermodynamic models to experimental data, such as the enthalpy of mixing or activity coefficients.^53,54 However, for a binary system, only a single scalar interaction parameter must account for all inter- and intra-molecular interactions for the entire configurational space of a thermodynamic state. This simplification limits the accuracy of the models, particularly when molecules exists in a variety of conformations resulting in molecular interactions with different orientations and at various distances. For example, COSMO-RS utilizes parameters based on a single molecular conformation.⁴⁹ This simplification limits the accuracy of the models, particularly when molecules exists in a variety of conformations that affect intermolecular interactions and their resulting phase behavior. Meanwhile, the MDKB → Phase approach naturally accounts for the entire thermodynamic ensemble of both intramolecular (conformational) and intermolecular (configurational) degrees of freedom, including any interactions between them. While the quality of the MDKB → Phase prediction naturally depends entirely the choice of MD force field, we might expect this predictive (for a given force field) methodology outperforms standard thermodynamic models, as we have shown for the DBP/n-dodecane system. The atomistic resolution of MD models also provides greater resolution to defining how molecules interact, and the resulting structuring in solution is accounted for explicitly in the RDFs that define the KBIs. By contrasting the MDKB → Phase approach with these established models, this work underscores the potential benefits of preserving molecular-level detail of ensembles of configurations when assessing the accuracy and predictive capabilities of various thermodynamic models.

Conclusions

This work introduces the MDKB → Phase approach to predict macroscopic liquid–liquid phase equilibria from molecular simulations. This purely predictive approach (for a given MD force field) retains atomic-scale information from MD simulations, providing a robust tool for analyzing complex thermodynamic systems. By separately quantifying the entropic and enthalpic contributions to the excess free energy G^E, the MDKB → Phase approach enables accurate estimation of ΔG_mix as a function of temperature from simulations at a single T. Equilibrating MD simulations becomes increasingly problematic at temperatures approaching the spinodal, even in the single-phase region, because of the increasing correlation lengths and times. The ability to extrapolate from MD simulations at temperatures well away from the spinodal solves this problem, and enables calculation of phase equilibria and scattering at any temperature. The temperature extrapolation assumes that the T dependence of the enthalpy and entropy of mixing are negligible, which we verified for the type of systems studied here. As in conventional thermodynamic models, the MDKB → Phase approach gives a linear T dependence of the free energy and scattering intensity, which corresponds to mean-field critical behavior. Thus corrections may be needed close to the spinodal. For systems with a strong temperature dependence of the enthalpy of mixing, extrapolation in temperature will require additional information, e.g. from additional simulations at varying temperature. Because of their direct connection to the free energy of mixing, SAXS measurements can verify the accuracy of the temperature extrapolation.

In this study, we employed polynomial functions to fit the excess thermodynamic properties. Although these simple forms do not uniformly capture the entire composition range, as shown in Fig. 3, they reproduce the behavior near the critical composition, where residual errors are minimized; this region is most consequential for predicting phase behavior. More complex, system-specific functions—particularly those that explicitly incorporate molar volumes—may improve accuracy across the full composition range. In addition to the shorter extrapolation in temperature, this could explain the better agreement with DBP compared to DOHE (Fig. 6). Future work will evaluate such formulations and quantify the trade-offs between added complexity and predictive performance.

The ability to accurately predict liquid–liquid critical points with the MDKB → Phase approach shows improvement for our systems over conventional thermodynamic models such as UNIQUAC, UNIFAC, and COSMO-RS. This is particularly evident in its replication of both the intensity and shape of the SAXS I₀ across different compositions for a number chemical systems with varying critical temperatures, suggesting accuracy is not the result of fortuitous error cancellation in enthalpic and entropic excess free energies, but a reflection of accurate estimation of each contribution to the mixing free energy.

The preservation of atomistic detail makes the MDKB → Phase approach exceptionally valuable for linking molecular interactions to macroscopic thermodynamic properties. While we have demonstrated the method here using binary systems, a notable strength of this approach is its inherent flexibility to extend to systems with multiple components.^55–57 Future work will explore the impact of various intermolecular interactions, such as hydrogen bonding, on the enthalpic and entropic contributions across different temperatures and compositions.

We note that an analogous strategy of combining KB theory with atomistic simulations to predict phase equilibria has also been applied to solid solutions by Miyaji et al.,⁴⁵ further underscoring the generality of this framework across different classes of materials.

Overall, this approach opens new avenues for the computational screening of thermodynamic properties, thereby enhancing the potential for materials discovery and design.

Author contributions

AAP: conceptualization, methodology, formal analysis, investigation, writing – original draft, review & editing; GBS: conceptualization, methodology, writing – review & editing; MJS: conceptualization, writing – original draft, review & editing, supervision, funding aquisition.

Conflicts of interest

There are no conflicts to declare.

Data availability

Data supporting the results of this study are available in the supplementary information (SI). The input parameters for molecular dynamics and code for analysis are available on request from the authors. Kirkwood-Buff integrals (KBIs) and related thermodynamic properties were computed using KBKit, an open-source software package developed by the authors and publicly available on GitHub (https://github.com/anl-sepsci/kbkit).

Supplementary information: system compositions; corrections to Kirkwood–Buff integrals, relations for thermodynamic properties from UNIQUAC and UNIFAC thermodynamic models; validation of MD pure component molar volumes to experiment; excess molar volumes; isothermal compressibility; validation of KB molar volumes and isothermal compressibility; error estimations for free energy derived properties; contribution of isothermal compressibility to scattering; comparison of enthalpy and entropy contributions to mixing free energy from thermodynamic models; temperature independence of mixing enthalpy and excess entropy; liquid–liquid phase equilibria; comparison of critical points for systems studied. See DOI: https://doi.org/10.1039/d6sc03294j.

Acknowledgements

The authors acknowledge support from the U.S. Department of Energy (DOE), Office of Science, Office of Basic Energy Sciences, Chemical Sciences, Geosciences, and Biosciences Division, Separation Science Program, under contract DE-AC02-06CH11357 to UChicago Argonne, LLC, Operator of Argonne National Laboratory. This research used resources of beamline 12-ID-C at the Advanced Photon Source, a U.S. DOE Office of Science User Facility, operated for the DOE Office of Science by Argonne National Laboratory under the same contract number. We gratefully acknowledge the computing resources provided on Improv, a high-performance computing cluster operated by the Laboratory Computing Resource Center at Argonne National Laboratory.

References

1. S. Rekhi, C. G. Garcia, M. Barai, A. Rizuan, B. S. Schuster, K. L. Kiick and J. Mittal, Expanding the molecular language of protein liquid–liquid phase separation, Nat. Chem., 2024, 16, 1113–1124.
2. H. Fu, J. Huang, J. J. van der Tol, L. Su, Y. Wang, S. Dey, P. Zijlstra, G. Fytas, G. Vantomme and P. Y. Dankers, et al., Supramolecular polymers form tactoids through liquid–liquid phase separation, Nature, 2024, 626, 1011–1018.
3. N. Dubouis, C. Park, M. Deschamps, S. Abdelghani-Idrissi, M. Kanduč, A. Colin, M. Salanne, J. Dzubiella, A. Grimaud and B. Rotenberg, Chasing aqueous biphasic systems from simple salts by exploring the LiTFSI/LiCl/H₂O phase diagram, ACS Cent. Sci., 2019, 5, 640–643.
4. J. M. Prausnitz, R. N. Lichtenthaler and E. G. de Azevedo, Molecular Thermodynamics of Fluid-Phase Equilibria, Prentice Hall Inc., Englewood Cliffs, 3rd edn, 1999.
5. D. Abranches and J. Coutinho, Everything You Wanted to Know about Deep Eutectic Solvents but Were Afraid to Be Told, Annu. Rev. Chem. Biomol. Eng., 2023, 14, 141–163.
6. N. Dawass, P. Krüger, S. K. Schnell, J.-M. Simon and T. Vlugt, Kirkwood-Buff integrals from molecular simulation, Fluid Phase Equilib., 2019, 486, 21–36.
7. A. O. Dohn, V. Markmann, A. Nimmrich, K. Haldrup, K. B. Møller and M. M. Nielsen, Eliminating finite-size effects on the calculation of X-ray scattering from molecular dynamics simulations, J. Chem. Phys., 2023, 159, 124115.
8. J. G. Kirkwood and F. P. Buff, The Statistical Mechanics Theory of Solutions. I, J. Chem. Phys., 1951, 19, 774–777.
9. E. A. Ploetz, G. N. Pallewela and P. E. Smith, Fluctuation solution theory of pure fluids, J. Chem. Phys., 2017, 146, 094501.
10. A. A. Chialvo, Linking Solution Microstructure and Solvation Thermodynamics of Mixed-Solvent Systems: Formal Results, Critical Observations, and Modeling Pitfalls, Thermo, 2024, 4, 407–432.
11. P. C. Petris, S. D. Anogiannakis, P.-N. Tzounis and D. N. Theodorou, Thermodynamic Analysis of n-Hexane–Ethanol Binary Mixtures Using the Kirkwood–Buff Theory, J. Phys. Chem. B, 2019, 123, 247–257.
12. F. Venetsanos, S. D. Anogiannakis and D. N. Theodorou, Mixing Thermodynamics and Flory–Huggins Interaction Parameter of Polyethylene Oxide/Polyethylene Oligomeric Blends from Kirkwood–Buff Theory and Molecular Simulations, Macromolecules, 2022, 55, 4852–4862.
13. M. Muniz, T. Gartner, C. Knight, S. Yue, F. Paesani and A. Panagiotopoulos, Vapor-liquid equilibrium of water with the MB-pol many-body potential, J. Chem. Phys., 2021, 154, 211103.
14. D. York, Modern Alchemical Free Energy Methods for Drug Discovery Explained, ACS Phys. Chem. Au, 2023, 3, 478–491.
15. T. Cheisson and E. J. Schelter, Rare earth elements: Mendeleev's bane, modern marvels, Science, 2019, 363, 489–493.
16. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindahl, GROMACS: high performance molecular simulations through multi-level parallelism from laptops to supercomputers, SoftwareX, 2015, 1–2, 19–25.
17. J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman and D. A. Case, Development and testing of a general amber force field, J. Comput. Chem., 2004, 25, 1157–1174.
18. T. Rahman, B. L. Bonnett, D. Poe, P. N. Wimalasiri, S. Seifert, J. Lal, G. B. Stephenson and M. J. Servis, Decomposition of liquid-liquid extraction organic phase structure into critical and pre-peak contributions, J. Mol. Liq., 2024, 393, 123625.
19. B. Bonnett, D. Sheyfer, P. N. Wimalasiri, S. Nayak, J. Lal, Q. Zhang, S. Seifert, G. B. Stephenson and M. J. Servis, Critical fluctuations in liquid-liquid extraction organic phases controlled by extractant and diluent molecular structure, Phys. Chem. Chem. Phys., 2023, 25, 16389–16403.
20. B. Bonnett, T. Rahman, D. Poe, S. Seifert, G. B. Stephenson and M. J. Servis, Insights into water extraction and aggregation mechanisms of malonamide-alkane mixtures, Phys. Chem. Chem. Phys., 2024, 26, 18089–18101.
21. B. L. Bonnett, P. N. Wimalasiri, D. Sheyfer, A. A. Peroutka, J. Lal, Q. Zhang, E. Dufresne, S. Narayanan, S. Seifert, G. B. Stephenson and M. J. Servis, Extracted Water Induces Concentration Fluctuations in Model Ternary Liquid–Liquid Extraction System, J. Phys. Chem. B, 2025, 129, 8177–8191.
22. L. Martínez, R. Andrade, E. G. Birgin and J. M. Martínez, PACKMOL: a package for building initial configurations for molecular dynamics simulations, J. Comput. Chem., 2009, 30, 2157–2164.
23. R. W. Hockney, S. P. Goel and J. W. Eastwood, Quiet high-resolution computer models of a plasma, J. Comput. Phys., 1974, 14, 148–158.
24. T. Darden, D. York and L. Pedersen, Particle mesh Ewald: an N log(N) method for Ewald sums in large systems, J. Chem. Phys., 1993, 98, 10089–10092.
25. H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola and J. R. Haak, Molecular dynamics with coupling to an external bath, J. Chem. Phys., 1984, 81, 3684–3690.
26. G. Bussi, D. Donadio and M. Parrinello, Canonical sampling through velocity rescaling, J. Chem. Phys., 2007, 126, 014101.
27. W. G. Hoover, Canonical dynamics: equilibrium phase-space distributions, Phys. Rev. A: At., Mol., Opt. Phys., 1985, 31, 1695–1697.
28. M. Bernetti and G. Bussi, Pressure control using stochastic cell rescaling, J. Chem. Phys., 2020, 153, 114107.
29. D. Hall, Kirkwood-Buff Theory of Solutions. An Alternative Derivation of Part of It and Some Applications, Trans. Faraday Soc., 1971, 67, 2516–2524.
30. K. Newman, Kirkwood-Buff Solution Theory: Derivation and Applications, Chem. Soc. Rev., 1994, 23, 31–40.
31. K. Nishikawa, The Solution Chemistry of Mixing States Probed via Fluctuations: A Direct Description of Inhomogeneity in Mixing, Bull. Chem. Soc. Jpn., 2021, 94, 2170–2186.
32. P. Ganguly and N. F. A. van der Vegt, Convergence of Sampling Kirkwood–Buff Integrals of Aqueous Solutions with Molecular Dynamics Simulations, J. Chem. Theory Comput., 2013, 9, 1347–1355.
33. M. Sevilla and R. Cortes-Huerto, Connecting density fluctuations and Kirkwood–Buff integrals for finite-size systems, J. Chem. Phys., 2022, 156, 044502.
34. J.-M. Simon, P. Krüger, S. K. Schnell, T. J. H. Vlugt, S. Kjelstrup and D. Bedeaux, Kirkwood–Buff integrals: from fluctuations in finite volumes to the thermodynamic limit, J. Chem. Phys., 2022, 157, 130901.
35. J. Milzetti, D. Nayar and N. F. A. van der Vegt, Convergence of Kirkwood–Buff Integrals of Ideal and Nonideal Aqueous Solutions Using Molecular Dynamics Simulations, J. Phys. Chem. B, 2018, 122, 5515–5526.
36. P. Krüger and T. Vlugt, Size and shape dependence of finite-volume Kirkwood-Buff integrals, Phys. Rev. E, 2018, 97, 051301.
37. A. Banerjee, M. Sevilla, J. Rudzinski and R. Cortes-Huerto, Finite-size scaling and thermodynamics of model supercooled liquids: long-range concentration fluctuations and the role of attractive interactions, Soft Matter, 2022, 18, 2373–2382.
38. N. Shcherbakova, V. Gerbaud and K. Roger, Using the Intrinsic Geometry of Binodal Curves to Simplify the Computation of Ternary Liquid–Liquid Phase Diagrams, Entropy, 2023, 25, 1329.
39. R. Busselez, Robust Kirkwood–Buff inversion in complex mixtures via reciprocal-space methods, J. Chem. Phys., 2025, 163, 134105.
40. B. Mukhopadhyay, B. Sabyasachi and M. J. Holdaway, A discussion of Margules-type formulations for multicomponent solutions with a generalized approach, Geochim. Cosmochim. Acta, 1993, 57, 277–283.
41. J. Tomiska and A. Neckel, The Margules Concept: The Basis of Modern Algebraic Representations of Thermodynamic Excess Properties, J. Phase Equilib., 1996, 17, 11–20.
42. A. Podgorsek, J. Jacquemin, A. Pádua and M. Costa Gomes, Mixing enthalpy for binary mixtures containing ionic liquids, Chem. Rev., 2016, 116, 6075–6106.
43. H. S. Ashbaugh and L. R. Pratt, Colloquium: scaled particle theory and the length scales of hydrophobicity, Rev. Mod. Phys., 2006, 78, 159–178.
44. G. B. Stephenson, W. K. Warburton, W. Haller and A. Bienenstock, Real-time small-angle X-ray scattering study of the early stage of phase separation in the SiO₂-BaO-K₂O system, Phys. Rev. B: Condens. Matter Mater. Phys., 1991, 43, 13417–13437.
45. M. Miyaji, J.-M. Simon and P. Krüger, Thermodynamic Analysis of Ar_xXe_1−x Solid Solutions Based on Kirkwood–Buff Theory, Physchem, 2022, 2, 191–206.
46. M. J. Servis, S. Nayak and S. Seifert, The pervasive impact of critical fluctuations in liquid-liquid extraction organic phases, J. Chem. Phys., 2021, 155, 244506.
47. D. S. Abrams and J. M. Prausnitz, Statistical Thermodynamics of Liquid Mixtures: A New Expression for the Excess Gibbs Energy of Partly or Completely Miscible Systems, AIChE J., 1975, 21, 116–128.
48. A. Fredenslund, R. Jones and J. M. Prausnitz, Group-contribution estimation of activity coefficients in nonideal liquid mixtures, AIChE J., 1975, 21, 1086–1099.
49. A. Klamt and F. Edkert, COSMO-RS: a novel and efficient method for the a priori prediction of thermophysical data of liquids, Fluid Phase Equilib., 2000, 172, 43–72.
50. T. F. Anderson and J. M. Prausnitz, Application of the UNIQUAC Equation to Calculation of Multicomponent Phase Equilibria. 2. Liquid-Liquid Equilibria, Ind. Eng. Chem. Process Des. Dev., 1978, 17, 561–568.
51. J. Li, B. Sundman, J. Winkelman, A. Vakis and F. Picchioni, Implementation of the UNIQUAC model in the OpenCalphad software, Fluid Phase Equilib., 2020, 507, 112398.
52. T. Gerlach, S. Müller, A. G. de Castilla and I. Smirnova, An open source COSMO-RS implementation and parameterization supporting the efficient implementation of multiple segment descriptors, Fluid Phase Equilib., 2022, 560, 113472.
53. H. Renon and J. M. Prausnitz, Local compositions in thermodynamic excess functions for liquid mixtures, AIChE J., 1968, 14, 135–144.
54. J. Gmehling, Present status and potential of group contribution methods for process development, J. Chem. Thermodyn., 2009, 41, 731–747.
55. A. Ben-Naim, Molecular Theory of Solutions, Oxford University Press, 2006.
56. M. Kang and P. E. Smith, Kirkwood–Buff theory of four and higher component mixtures, J. Chem. Phys., 2008, 128, 244511.
57. S. Seo, H.-s. Lee and T. J. Yoon, Kirkwood-Buff Analysis of Binary and Ternary Systems Consisting of Alcohols (Methanol, Ethanol, 1-Propanol, and 2-Propanol), Water, and n-Hexane to Understand the Formation of Surfactant-Free Microemulsions, J. Phys. Chem. B, 2024, 128, 5092–5108.

---