Direct access to thermodynamics from equilibrium compositional heterogeneities

Abstract

Equilibrium compositional heterogeneities contain thermodynamic information that is typically accessed through multiple bulk states. Here, we show that this information can be extracted from a single system, providing access to thermodynamic properties over a broad range of local compositions. In practice, a single equilibrium simulation is sufficient to recover this behavior. In the system considered here, the excess chemical potential of water follows a simple and thermodynamically consistent dependence on local composition, enabling the recovery of the concentration-dependent extraction behavior across extractant concentrations from a single liquid–liquid system. Within the present framework, this approach provides access to equilibrium distributions across compositions from a single system. The resulting relationship quantitatively reproduces experimental data and retains predictive power well beyond the range over which it is established. Remarkably, the same functional dependence is recovered from local compositional variations in homogeneous systems, indicating that it suggests an intrinsic thermodynamic property rather than a feature specific to interfacial environments. The corresponding response coefficient suggests a possible connection with correlation-based statistical mechanical descriptions such as Kirkwood–Buff-type approaches. Altogether, these results show that compositional heterogeneities can be used as a thermodynamic sampling approach, enabling access to thermodynamic properties without sampling multiple bulk states or explicitly computing intermolecular correlation functions. Because this approach relies solely on equilibrium local composition variations, it is not restricted to molecular simulations and can, in principle, be extended to experimental measurements of local composition.

---

Introduction

Understanding how thermodynamic properties emerge from molecular organization remains a central challenge in multicomponent liquids.^1,2 In particular, establishing quantitative connection between local composition and chemical potentials is essential to describe processes such as molecular partitioning, transfer, and extraction, and to predict how species distribute between phases under equilibrium conditions.^3–5 In homogeneous systems, thermodynamic properties can be related to density fluctuations and intermolecular correlations, notably within the framework of Kirkwood–Buff theory.^6,7 However, extracting such information typically requires extensive sampling or multiple simulations at different bulk compositions, limiting both efficiency and practical applicability in molecular simulations.

Multicomponent systems at equilibrium often exhibit spatial compositional heterogeneities, such as those found at liquid–liquid interfaces.^8–10 In such systems, strong interfacial enrichment of specific species, such as extractant molecules, is commonly observed, leading to pronounced gradients in local composition across the interface. These heterogeneities span a continuous range of local environments and therefore contain, within a single system, information that is traditionally accessed by sampling multiple bulk states at different compositions. Importantly, such variations arise spontaneously under equilibrium conditions and are inherently sampled in standard molecular simulations. Yet, the thermodynamic information accessible from these routinely generated configurations remains largely unexploited.

In particular, the concentration-dependent extraction behavior of a molecular species can be recovered from a single heterogeneous system, as illustrated here for water in a liquid–liquid extraction system. By relating the excess chemical potential of a molecular species to its local environment, a simple and thermodynamically consistent relationship emerges between chemical potential and composition. This relationship is quantitatively predictive across a wide range of extractant concentrations and retains predictive power well beyond the range over which it is established, suggesting its underlying thermodynamic origin. Moreover, it is recovered from both interfacial gradients and local equilibrium composition variations, and is validated against independent experimental data.

These results suggest that compositional heterogeneities can be viewed as a thermodynamic sampling framework, in which equilibrium compositional variations contain the underlying thermodynamic response of the system. In this perspective, information governing equilibrium distributions between phases, including partitioning behavior and extraction selectivity, can be inferred from a single heterogeneous system, without the need to compute intermolecular correlation functions or sample multiple bulk states.

Results

To probe how compositional gradients provide access to thermodynamic information, we consider a model system consisting of two phases in contact: a pure aqueous phase and an organic phase containing an aliphatic solvent (here n-dodecane) and extractant molecules representative of amphiphilic extractants commonly employed in nuclear fuel reprocessing (TBP molecules).^11,12 As we show below, this framework shows that concentration-dependent water extraction behavior as a function of extractant concentration is accessible from the composition variations of a single system. The chosen extractant molecules strongly solubilize water in the organic phase, providing a stringent test case that has been extensively investigated experimentally and computationally.^13–17 For clarity, water molecules are hereafter denoted as species A, and extractant molecules as species B when using a simplified notation.

Fig. 1 depicts the molecular system at thermodynamic equilibrium, highlighting the liquid–liquid interface that mediates water extraction. At this interface, extractant molecules align with their polar headgroups oriented toward the aqueous phase. Water transfer into the organic phase occurs primarily via protrusion events, i.e., interfacial fluctuations that lead to the formation and insertion of water clusters of varying sizes.^16,17 These clusters are subsequently stabilized in the organic phase through hydrogen bonding with the extractant molecules. The system therefore samples a range of molecular environments, from simple water-extractant interactions to larger, dynamically evolving aggregates, whose detailed characterization is beyond the scope of the present work (representative examples of such structures can be found in ref. 16 and 17). The interface thus acts as the key region governing exchange between the aqueous reservoir and the organic phase. During the simulations, water molecules continuously cross the interface in both directions, with no net flux between the two phases, indicating that equilibrium is reached (see SI, Section S2). For the analysis below, the z axis is defined as perpendicular to the liquid–liquid interface (see Fig. 1).

Fig. 1 Representative snapshots of the simulated system illustrating the liquid–liquid interface and the spatial organization of the different species. For clarity and visualization purposes, only half of the periodic simulation box is shown here. The full simulation box, including both liquid–liquid interfaces, is provided in the supporting information (SI). (a) Snapshot containing all components (water, n-dodecane and TBP molecules). (b) The system where only water molecules (species A) are shown for sake of clarity. Water and n-dodecane molecules are shown in blue and orange color respectively. For TBP molecules, oxygen, phosphorous and carbon atoms are shown in red, tan and black. The aliphatic hydrogen atoms are not shown for more clarity.

The excess chemical potential of water (µ^exc_A) and the concentration of extractant molecules (c_B) along the z axis are depicted in Fig. 2 and reveal how spatial variations in composition are coupled to systematic changes in thermodynamic stability. The excess chemical potential is obtained from the equilibrium distribution along z, corresponding to an unbiased potential of mean force (see SI, Section S3 for details). In the bulk aqueous phase, µ^exc_A remains constant and equal to 0 kcal mol⁻¹, consistent with the choice of pure water reference state. Across the interfacial region, µ^exc_A increases progressively, indicating a reduced thermodynamic stability of water in this region. Upon entering the organic phase, µ^exc_A reveals a local stabilization associated with the transfer of water clusters across the interface. Deeper into the organic phase, µ^exc_A increases gradually until reaching a plateau at 3.36 kcal mol⁻¹ corresponding to the bulk organic region at its equilibrium extractant concentration, and thus defining the free energy cost of water extraction. The concentration profile c_B exhibits a complementary behavior. While extractant molecules are essentially absent from the aqueous phase, c_B increases sharply across the interface, reaching a maximum of 3.40 mol L⁻¹. This interfacial enrichment reflects the strong affinity of the extractant for the interface and its role in stabilizing transferred water. Upon entering the organic phase, c_B displays a local minimum (2.23 mol L⁻¹), consistent with the region where water clusters are exchanged between phases. Further into the organic phase, c_B gradually decreases toward a plateau at 0.61 mol L⁻¹, corresponding to the bulk organic composition.

Fig. 2 Excess chemical potential of A (µ^exc_A) and concentration of B (c_B) species along the axis perpendicular to the liquid–liquid interface (z axis). The dotted black vertical line indicates the approximate location of the liquid–liquid interface.

In the organic phase, both profiles relax over several tens of angstroms along the z-axis before reaching the bulk region, indicating that water extraction is not localized at the interface but extends deep into the organic phase. This behavior reflects the progressive stabilization of water by extractant molecules. As extractants accumulate near the interface, they stabilize water species, leading to a local enrichment. These stabilized species then diffuse into the organic phase, in dynamic exchange with the aqueous phase, resulting in the gradual relaxation of both concentration and chemical potential profiles toward their bulk values. From a conventional perspective, this system would simply correspond to water extraction at a fixed extractant concentration (∼0.61 mol L⁻¹ in the bulk organic phase), requiring independently equilibrated systems at different compositions to reconstruct the concentration dependence of extraction behavior. However, as shown below, the presence of compositional variations provides access to a broader thermodynamic picture.

It is worth noting that resolving such compositional variations requires extensive sampling. In practice, several microseconds of simulation were necessary to reach equilibrium and obtain well-converged profiles (see Fig. S3, SI), which may explain why such behavior has rarely been reported. Importantly, this extended relaxation could be attributed to electrostatic effects arising from interfacial charge accumulation, as commonly observed in systems driven by an applied potential.^23–25 However, the present system is at equilibrium, and no significant charge density is detected in the region where the relaxation occurs (see Fig. S5, SI) confirming that this behavior originates from equilibrium compositional variations.

Although the profiles along the z-axis describe the spatial organization of the system, they do not directly reveal the thermodynamic dependence of µ^exc_A. To uncover this relationship, we express µ^exc_A as a function of the local concentration of extractant molecules. This mapping is particularly well-defined in the present system due to the negligible solubility of water in the aliphatic solvent in the absence of extractant molecules. When expressed as a function of local composition, the data collapse onto a single, well-defined thermodynamic relationship (see Fig. 3). Remarkably, the same behavior is recovered from experimental data reported in the literature for TBP/water/n-dodecane systems.^13,18–20 The TBP/water extraction system has been extensively investigated experimentally in the literature using a variety of organic diluents.^13,15,21,22 In the present work, the quantitative comparison primarily focuses on n-dodecane systems to maintain consistency with the molecular environment used in the simulations. For completeness, additional literature data obtained in kerosene or closely related aliphatic hydrocarbon diluents are also included in Fig. 3 and exhibit a broadly consistent trend. These results suggest that local equilibrium composition variations contain sufficient thermodynamic information to reconstruct the concentration-dependent response observed in the present system. The collapse onto a single curve indicates that the thermodynamic response of the system is governed by a simple functional dependence on the local extractant concentration c_B:

Fig. 3 Excess chemical potential of A (µ^exc_A) as a function of the concentration of B (c_B) showing the collapse of all data onto a single master curve. Experimental data are included for comparison and exhibit quantitative agreement with the simulation results.^{13,15,18–22} Experimental values were converted into excess chemical potentials using µ^exc_A = −RT ln(c_A/c°), where c° = 55.5 mol L⁻¹ corresponds to pure water. Data marked with an asterisk (*) were obtained in kerosene or related aliphatic hydrocarbon diluents rather than in n-dodecane.

The expression involves a small number of physically meaningful parameters, providing a compact description of the thermodynamic response. The quantity µ^exc_A denotes the excess chemical potential of species A, defined relative to the chosen reference state (pure water in the present work). The variable c_B corresponds to the local concentration of extractant molecules, which governs the stabilization of A in the organic phase. The parameter µ^exc_A,0 represents the excess chemical potential of A at a reference extractant concentration (c°), while the dimensionless coefficient α quantifies the sensitivity of µ^exc_A to variations in c_B. The parameter α thus emerges as an effective thermodynamic susceptibility, as it quantifies the response of µ^exc_A to variations in local composition within the present system. The thermal energy scale is set by RT. Water extraction behavior in other extractant systems has also been rationalized using explicit thermodynamic modeling approaches based on equilibrium solution theories.²⁶ These experimental values are typically obtained from bulk systems at different compositions. In contrast, the present approach samples a continuous range of concentrations within a single simulation through the spatial gradient. Each position along the z-axis thus corresponds to a local environment analogous to a bulk system at a given composition, effectively probing a continuous range of local equilibrium compositions. More generally, Kirkwood–Buff-type fluctuation theories are known to relate chemical potential variations to intermolecular correlation integrals in multicomponent systems, suggesting that the present behavior may reflect an underlying correlation-driven thermodynamic response. The present logarithmic form can therefore be viewed as a physically motivated effective description corresponding to the leading-order concentration dependence of µ^exc_A on local composition. Using the thermodynamic definition of the chemical potential, this relationship can be directly expressed in terms of the partition coefficient D leading to a power-law dependence on extractant concentration: D ∝ (c_B/c°)^α, up to a constant prefactor (see Section S3.4, SI).

The proposed expression accurately describes the data shown in Fig. 3, with a root mean square deviation of 9.86 × 10⁻³ kcal mol⁻¹. This value is more than an order of magnitude smaller than the thermal energy RT (∼0.596 kcal mol⁻¹ at 300 K), indicating that the functional form accurately captures the thermodynamic behavior of the system. Remarkably, this single curve, obtained from a single equilibrium simulation, reproduces the concentration-dependent thermodynamic response observed over a broad range of extractant concentrations. In this sense, it effectively contains the equivalent of an entire set of extraction data across compositions, which would traditionally require multiple independent bulk simulations or experiments. The parameter µ^exc_A,0 is found to be 2.944 kcal mol⁻¹, corresponding to the excess chemical potential of water at the reference concentration c° (1 mol L⁻¹). The thermodynamic susceptibility α is found to be 1.400. While µ^exc_A,0 depends on the chosen reference state and concentration c°, α remains invariant, highlighting its intrinsic thermodynamic character. Importantly, α is not a local structural quantity and cannot be reduced to a coordination number. Instead, it reflects a collective response of the system, extending beyond short-range interactions. These results suggest that heterogeneous equilibrium simulations can provide access to concentration-dependent thermodynamic response properties without explicitly computing intermolecular correlation functions. More broadly, it highlights the ability of the present framework to extract thermodynamic information that is typically accessible only through more involved theoretical or computational approaches.²⁷

Fig. 4 shows µ^exc_A as a function of c_B using the fitted relationship. To assess the predictive capability, the model was evaluated well beyond the calibration range by extending the curve up to c_B ∼ 3.6 mol L⁻¹, corresponding to a pure extractant phase (i.e., in the absence of aliphatic solvent). An independent simulation of pure water in contact with a pure extractant solution was then performed to validate this extrapolation (see Fig. 5). This simulation was also carried out on the microsecond timescale to ensure thorough sampling of equilibrium and a well-converged value of µ^exc_A (see SI, Section S1 and S2). The resulting value of 1.885 kcal mol⁻¹ for c_B ∼ 3.53 mol L⁻¹ is also shown in Fig. 4. Remarkably, the independently obtained value lies directly on the extrapolated curve (deviation < 7 × 10⁻³ kcal mol⁻¹), demonstrating that the relationship retains predictive power far beyond the range of compositions directly sampled.

Fig. 4 Excess chemical potential of A (µ^exc_A) as a function of the extractant concentration (c_B). The shaded orange region, delimited by dashed grey lines, indicates the composition range used to establish the relationship. The value of µ^exc_A in the pure extractant phase, obtained from an independent molecular simulation, is also shown for comparison, highlighting the predictive capability of the approach beyond this range.

Fig. 5 Excess chemical potential of A (µ^exc_A) and concentration of B (c_B) species along the axis perpendicular to the liquid–liquid interface (z axis) for a system consisting of a pure water phase in contact with a pure extractant phase. The dotted black vertical line indicates the approximate location of the liquid–liquid interface.

We now compare these results with available experimental data. As demonstrated above, the relationship presented in eqn (1) accurately captures the thermodynamic behavior of the system. Consistent with Fig. 3, excellent agreement with experimental data is obtained over a broad range of intermediate extractant concentrations. Deviations emerge only near the pure extractant limit. In this regime, experimental values differ from the present prediction (experimentally ∼3.40 mol L⁻¹ of water into pure extractant phase,^{13,15,19,21,22} here 2.35 mol L⁻¹), which may reflect limitations of the force field used for the molecular simulations. Importantly, this discrepancy does not affect the validity of the proposed thermodynamic relationship, which remains internally consistent and quantitatively predictive within the simulation framework. Furthermore, good agreement is also obtained for significantly lower extractant concentrations than those directly sampled in the simulations and therefore outside the fitting range used to establish the thermodynamic relationship (e.g. 0.06 mol L⁻¹ experimentally¹⁹ vs. 0.09 mol L⁻¹ at c_B = 0.35 mol L⁻¹, and, 0.13 mol L⁻¹ experimentally in n-octane²⁸ vs. 0.15 mol L⁻¹ at c_B = 0.50 mol L⁻¹). This further supports the robustness of the proposed concentration-dependent behavior in the extrapolated dilute regime.

While the above analysis relies on spatial concentration gradients, the observed relationship is not specific to strongly inhomogeneous conditions. Instead, it reflects the existence of local composition-dependent thermodynamic environments accessible within equilibrium configurations of the system. To illustrate this point, we analyze local composition variations in a configuration where pure water and pure extractant phases are in contact. Fig. 5 shows the profiles of µ^exc_A and c_B in this system. In this case, the concentration gradients in the organic phase are much weaker than in the diluted system. Nevertheless, finite local variations in c_B and µ^exc_A remain present within the bulk organic region (from z ∼ 75 to 0 Å) and can be analyzed using the same thermodynamic framework. The z-averaged concentration profiles shown here provide a coarse-grained description of the local compositional environments sampled during the simulations. The present analysis relies on finite-time averaged concentration profiles. In the limit of infinite sampling, these profiles would converge toward the corresponding bulk equilibrium composition. The present work does not attempt to quantify the underlying equilibrium concentration fluctuations directly. Nevertheless, such fluctuations and the associated intermolecular correlations remain present in homogeneous systems, suggesting that the observed behavior may ultimately be connected to a more fundamental fluctuation-based thermodynamic description. Despite the limited amplitude of the accessible composition variations in this nearly homogeneous system (see Fig. 5), the same functional dependence is recovered over the explored concentration range, from pure extractant conditions to dilute regimes. The fitted parameters remain consistent, with µ^exc_A,0 = 2.965 kcal mol⁻¹ (vs. 2.944 kcal mol⁻¹ with the gradient) and α = 1.419 (vs. 1.400), and a low RMS deviation of 6.79 × 10⁻³ kcal mol⁻¹. This result suggests that the observed relationship is not generated by the macroscopic concentration gradient itself. Rather, the gradient primarily acts as an efficient mechanism to broaden the range of locally accessible compositions and thereby improve the statistical sampling of the concentration-dependent thermodynamic response.

Importantly, the relationship presented in eqn (1) remains robust despite the underlying microscopic complexity of the system, including transient aggregates of varying size and composition. This robustness suggests that the excess chemical potential is governed by collective thermodynamic effects extending beyond specific local structures. In this context, the thermodynamic susceptibility α can be written as:

which quantifies the response of the excess chemical potential to variations in local composition. More generally, statistical mechanical theories of solutions, including Kirkwood–Buff-type approaches, relate chemical potential variations to intermolecular correlation functions and concentration-dependent collective organization.^1,29 In this perspective, the present logarithmic thermodynamic response may reflect an underlying correlation-driven behavior, with α acting as an effective thermodynamic response parameter associated with the local molecular environment. However, establishing a rigorous formal connection between the present framework and fluctuation-based statistical mechanical theories remains beyond the scope of the present work.

Conclusion

In this work, we demonstrate that the concentration-dependent extraction behavior of a molecular species across compositions can be recovered from a single equilibrium system by exploiting compositional heterogeneities. By leveraging spatial composition variations within a heterogeneous system, we show that partitioning, solvation, and extraction properties can be inferred without the need to sample multiple bulk states or explicitly compute intermolecular correlation functions. We further show that the excess chemical potential can be expressed as a simple function of local composition, leading to a collapse of the data onto a single master curve and enabling quantitative predictions across a wide range of extractant concentrations. This relationship is validated against independent experimental data and retains predictive power well beyond the range over which it is established, including in extreme compositional regimes (pure water phase in contact with the pure extractant phase).

More fundamentally, this behavior is not specific to concentration gradients, but reflects an intrinsic property of local equilibrium composition variations. The introduction of a thermodynamic susceptibility provides a compact and physically meaningful description of the response of the system, which may be viewed as a coarse-grained representation of the underlying intermolecular correlations. Altogether, these results suggest that local equilibrium composition variations may contain useful thermodynamic information related to macroscopic concentration-dependent behavior in complex molecular systems. In the context of liquid–liquid extraction, this approach provides a practical route to predict extraction behavior across compositions from a single equilibrium simulation.

This framework can be extended to more complex multicomponent systems, including liquid–liquid extraction and multicomponent partitioning, where competing interactions and speciation effects can be captured from local equilibrium composition variations or concentration gradients. Importantly, while molecular simulations provide a convenient route to access well-resolved compositional gradients, the underlying principle is not restricted to computational approaches. Liquid–liquid interfaces inherently sample a continuum of local environments that can, in principle, be probed experimentally. Recent advances in techniques such as X-ray reflectivity, which provide access to interfacial composition profiles,^10,30 therefore open the possibility of extracting similar thermodynamic information from experiments.

Methods

Molecular dynamics simulations were performed using the AMBER20 software package.³¹ The systems consist of a liquid–liquid interfaces between an aqueous phase and an organic phase composed of n-dodecane (in the diluted case only) and TBP extractant molecules. Large-scale simulations comprising ∼200 000 atoms were carried out over several microseconds (up to ∼8 µs for the diluted extractant and ∼5 µs for the pure extractant phase) to ensure adequate sampling of equilibrium compositional fluctuations. Classical molecular dynamics methodologies were employed, including validated force fields and established algorithms for temperature and pressure control, long-range electrostatics, and bond constraints (see SI, Section S1 and S2, additional methodological discussion is provided in the SI).^32–41 Simulations were performed at 300 K and 1 atm. The excess chemical potential of species A (water) was obtained from equilibrium concentration profiles, without the use of biasing techniques. In practice, this quantity, µ^exc_A, is obtained from the local probability of occurrence of species A, P_A(z), along the interface (see SI, Section S3), such that µ^exc_A(z) ∝ −RT ln(P_A(z)).⁴²

Author contributions

The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting the findings of this study are available within the article and its supporting information (SI), and from the corresponding author upon reasonable request. See DOI: https://doi.org/10.1039/d6cp01301e.

Acknowledgements

The author thanks Dr Philippe Guilbaud for comments and careful reading of the manuscript.

References

1. J. G. Kirkwood and F. P. Buff, The Statistical Mechanical Theory of Solutions. I, J. Chem. Phys., 1951, 19(6), 774–777,  DOI:10.1063/1.1748352.
2. K. E. Newman, Kirkwood–Buff Solution Theory: Derivation and Applications, Chem. Soc. Rev., 1994, 23(1), 31–40,  10.1039/CS9942300031.
3. A. Ben-Naim, Molecular Theory of Solutions, OUP Oxford, 2006.
4. B. Widom, Some Topics in the Theory of Fluids, J. Chem. Phys., 1963, 39(11), 2808–2812,  DOI:10.1063/1.1734110.
5. A. Grossfield, P. N. Patrone, D. R. Roe, A. J. Schultz, D. W. Siderius and D. M. Zuckerman, Best Practices for Quantification of Uncertainty and Sampling Quality in Molecular Simulations [Article v1.0], Living J. Comput. Mol. Sci., 2018, 1(1), 5067,  DOI:10.33011/livecoms.1.1.5067.
6. M. Sevilla and R. Cortes-Huerto, Connecting Density Fluctuations and Kirkwood–Buff Integrals for Finite-Size Systems, J. Chem. Phys., 2022, 156(4), 044502,  DOI:10.1063/5.0076744.
7. M. Sevilla, L. A. Baptista, K. Kremer and R. Cortes-Huerto, Density Fluctuations, Solvation Thermodynamics, and Coexistence Curves in Grand Canonical Molecular Dynamics Simulations, J. Chem. Phys., 2025, 162(8), 080901,  DOI:10.1063/5.0243895.
8. I. Benjamin, Chemical Reactions and Solvation at Liquid Interfaces:  A Microscopic Perspective, Chem. Rev., 1996, 96(4), 1449–1476,  DOI:10.1021/cr950230+.
9. P. Jungwirth and D. J. Tobias, Specific Ion Effects at the Air/Water Interface, Chem. Rev., 2006, 106(4), 1259–1281,  DOI:10.1021/cr0403741.
10. E. Scoppola, E. B. Watkins, R. A. Campbell, O. Konovalov, L. Girard, J.-F. Dufrêche, G. Ferru, G. Fragneto and O. Diat, Solvent Extraction: Structure of the Liquid–Liquid Interface Containing a Diamide Ligand, Angew. Chem., 2016, 128(32), 9472–9476,  DOI:10.1002/ange.201603395.
11. A. Wright and P. Paviet-Hartmann, Review of Physical and Chemical Properties of Tributyl Phosphate/Diluent/Nitric Acid Systems, Sep. Sci. Technol., 2010, 45(12–13), 1753–1762,  DOI:10.1080/01496395.2010.494087.
12. L. L. Burger, The Chemistry of Tributyl Phosphate: A Review, Hanford Atomic Products Operation, Chemical Separations Unit, Chemical Research Sub-Section, 1955.
13. T. H. Russell, R. M. Counce, J. S. Watson, B. B. Spencer and G. D. Del Cul, Water Content of Organic Solvents and Their Relationship to Extraction of Nitric and Acetic Acid from UREX+ Streams, Sep. Sci. Technol., 2013, 48(14), 2164–2169,  DOI:10.1080/01496395.2013.794361.
14. H. Naganawa and S. Tachimori, Hydration and Dimerization of Tributyl Phosphate in Dodecane, Anal. Sci., 1994, 10(2), 309–314,  DOI:10.2116/analsci.10.309.
15. L. Burger and R. Forsman, The Solubility of Tributyl Phosphate in Aqueous Solutions, HW-20936, 4349304; 1951; p. HW-20936, 4349304 DOI:10.2172/4349304.
16. M. J. Servis and A. E. Clark, Surfactant-Enhanced Heterogeneity of the Aqueous Interface Drives Water Extraction into Organic Solvents, Phys. Chem. Chem. Phys., 2019, 21(6), 2866–2874,  10.1039/C8CP06450D.
17. X. Ye, S. Cui, V. F. de Almeida and B. Khomami, Molecular Simulation of Water Extraction into a Tri-n-Butylphosphate/n-Dodecane Solution, J. Phys. Chem. B, 2013, 117(47), 14835–14841,  DOI:10.1021/jp409332b.
18. R. J. Ellis, T. L. Anderson, M. R. Antonio, A. Braatz and M. Nilsson, A SAXS Study of Aggregation in the Synergistic TBP–HDBP Solvent Extraction System, J. Phys. Chem. B, 2013, 117(19), 5916–5924,  DOI:10.1021/jp401025e.
19. P. Guilbaud, L. Berthon, W. Louisfrema, O. Diat and N. Zorz, Determination of the Structures of Uranyl–Tri-n-Butyl-Phosphate Aggregates by Coupling Experimental Results with Molecular Dynamic Simulations, Chem. – Eur. J., 2017, 23(65), 16660–16670,  DOI:10.1002/chem.201703967.
20. A. Suresh, T. G. Srinivasan and P. R. Vasudeva Rao, The Effect of the Structure of Trialkyl Phosphates on Their Physicochemical Properties and Extraction Behavior, Solvent Extr. Ion Exch., 2009, 27(2), 258–294,  DOI:10.1080/07366290802674481.
21. C. J. Hardy, D. Fairhurst, H. A. C. McKay and A. M. Willson, Extraction of Water by Tri-n-Butyl Phosphate, Trans. Faraday Soc., 1964, 60, 1626,  10.1039/tf9646001626.
22. K. Alcock, S. S. Grimley, T. V. Healy, J. Kennedy and H. A. C. McKay, The Extraction of Nitrates by Tri-n-Butyl Phosphate (TBP). Part 1.—The System TBP + Diluent + H₂O + HNO₃, Trans. Faraday Soc., 1956, 52, 39–47.
23. S. K. Reed, O. J. Lanning and P. A. Madden, Electrochemical Interface between an Ionic Liquid and a Model Metallic Electrode, J. Chem. Phys., 2007, 126(8), 084704,  DOI:10.1063/1.2464084.
24. C. Merlet, C. Péan, B. Rotenberg, P. A. Madden, B. Daffos, P.-L. Taberna, P. Simon and M. Salanne, Highly Confined Ions Store Charge More Efficiently in Supercapacitors, Nat. Commun., 2013, 4(1), 2701,  DOI:10.1038/ncomms3701.
25. A. A. Kornyshev, Double-Layer in Ionic Liquids:  Paradigm Change?, J. Phys. Chem. B, 2007, 111(20), 5545–5557,  DOI:10.1021/jp067857o.
26. B. L. 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(26), 18089–18101,  10.1039/D4CP01369G.
27. N. Dawass, P. Krüger, S. K. Schnell, J.-M. Simon and T. J. H. Vlugt, Kirkwood-Buff Integrals from Molecular Simulation, Fluid Phase Equilib., 2019, 486, 21–36,  DOI:10.1016/j.fluid.2018.12.027.
28. Y. Ueda, K. Kikuchi, T. Sugita and R. Motokawa, Extraction Performance of a Fluorous Phosphate for Zr(IV) from HNO3 Solution: Comparison with Tri-n-Butyl Phosphate, Solvent Extr. Ion Exch., 2019, 37(5), 347–359,  DOI:10.1080/07366299.2019.1638015.
29. J. G. Kirkwood, Statistical Mechanics of Fluid Mixtures, J. Chem. Phys., 1935, 3(5), 300–313,  DOI:10.1063/1.1749657.
30. M. L. Schlossman, Liquid–Liquid Interfaces: Studied by X-Ray and Neutron Scattering, Curr. Opin. Colloid Interface Sci., 2002, 7(3), 235–243,  DOI:10.1016/S1359-0294(02)00053-5.
31. D. A. Case; H. M. Aktulga; K. Belfon; I. Y. Ben-Shalom; S. R. Brozell; D. S. Cerutti; T. E. Cheatham III; V. W. D. Cruzeiro; T. A. Darden; R. E. Duke; G. Giambasu; M. K. Gilson; H. Gohlke; A. W. Goetz; R. Harris; S. Izadi; S. A. Izmailov; C. Jin; K. Kasavajhala; M. C. Kaymak; E. King; A. Kovalenko; T. Kurtzman; T. S. Lee; S. LeGrand; P. Li; C. Lin; J. Liu; T. Luchko; R. Luo; M. Machado; V. Man; M. Manathunga; K. M. Merz; Y. Miao; O. Mikhailovskii; G. Monard; H. Nguyen; K. A. O’Hearn; A. Onufriev; F. Pan; S. Pantano; R. Qi; A. Rahnamoun; D. R. Roe; A. Roitberg; C. Sagui; S. Schott-Verdugo; J. Shen; C. L. Simmerling; N. R. Skrynnikov; J. Smith; J. Swails; R. C. Walker; J. Wang; H. Wei; R. M. Wolf; X. Wu; Y. Xue; D. M. York; S. Zhao and P. A. Kollman, Amber 20, 2021.
32. 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(9), 1157–1174,  DOI:10.1002/jcc.20035.
33. J. Wang and P. A. Kollman, Automatic parameterization of force field by systematic search and genetic algorithms, J. Comput. Chem., 2001, 22(12), 1219–1228,  DOI:10.1002/jcc.1079.
34. S. W. I. Siu, K. Pluhackova and R. A. Böckmann, Optimization of the OPLS-AA Force Field for Long Hydrocarbons, J. Chem. Theory Comput., 2012, 8(4), 1459–1470,  DOI:10.1021/ct200908r.
35. S. Izadi, R. Anandakrishnan and A. V. Onufriev, Building Water Models: A Different Approach, J. Phys. Chem. Lett., 2014, 5(21), 3863–3871,  DOI:10.1021/jz501780a.
36. C. I. Bayly, P. Cieplak, W. Cornell and P. A. Kollman, A Well-Behaved Electrostatic Potential Based Method Using Charge Restraints for Deriving Atomic Charges: The RESP Model, J. Phys. Chem., 1993, 97(40), 10269–10280,  DOI:10.1021/j100142a004.
37. M. J. Frisch; G. W. Trucks; H. B. Schlegel; G. E. Scuseria; M. A. Robb; J. R. Cheeseman; G. Scalmani; V. Barone; B. Mennucci and G. A. Petersson, et al., Gaussian 09, Revision D.01, Gaussian. Inc., Wallingford CT, 2009.
38. 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(13), 2157–2164,  DOI:10.1002/jcc.21224.
39. 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(12), 10089–10092,  DOI:10.1063/1.464397.
40. J.-P. Ryckaert, G. Ciccotti and H. J. C. Berendsen, Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of n-Alkanes, J. Comput. Phys., 1977, 23(3), 327–341,  DOI:10.1016/0021-9991(77)90098-5.
41. 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(8), 3684–3690,  DOI:10.1063/1.448118.
42. B. Roux, The Calculation of the Potential of Mean Force Using Computer Simulations, Comput. Phys. Commun., 1995, 91(1), 275–282,  DOI:10.1016/0010-4655(95)00053-I.

---