Phase behaviour of coarse-grained fluids

Abstract

Soft condensed matter structures often challenge us with complex many-body phenomena governed by collective modes spanning wide spatial and temporal domains. In order to successfully tackle such problems, mesoscopic coarse-grained (CG) statistical models are being developed, providing a dramatic reduction in computational complexity. CG models provide an intermediate step in the complex statistical framework of linking the thermodynamics of condensed phases with the properties of their constituent atoms and molecules. These allow us to offload part of the problem to the CG model itself and reformulate the remainder in terms of reduced CG phase space. However, such exchange of pawns to chess pieces, or ‘Hamiltonian renormalization’, is a radical step and the thermodynamics of the primary atomic and CG models could be quite distinct. Here, we present a comprehensive study of the phase diagram including binodal and interfacial properties of a dissipative particle dynamics (DPD) model, extended to include finite-range attraction to support the liquid–gas equilibrium. Despite the similarities with the atomic model potentials, its phase envelope is markedly different featuring several anomalies such as an unusually broad liquid range, change in concavity of the liquid coexistence branch with variation of the model parameters, volume contraction on fusion, temperature of maximum density in the liquid phase and negative thermal expansion in the solid phase. These results provide new insight into the connection between simple potential models and complex emergent condensed matter phenomena.

---

Introduction

Coarse-graining (CG) is now well established as an essential ingredient of multiscale modelling frameworks¹ and provides accelerated pathways to characterizing soft condensed matter² including complex fluids,³ biological membranes⁴ and ionic solutions.⁵ Despite the inevitable loss in structural detail at the CG level, its enormous potential in expanding the length and time scales to address large collective phenomena drives current advances in complex condensed phases characterization. Soft matter often contains colloids and polymer solutions⁶ and in order to achieve the desired speedup in computer simulation of such systems both the solutes and solvents need to be coarse-grained congruently. Most forthrightly, CG can be carried out in biological and polymer systems by mapping (macro-)molecular fragments and moieties of chemically connected atomic fragments to ‘beads’ interacting via CG forces.^7,8 Here, CG methods have scored rapid successes from their arrival in the late 1960s, culminating in the 2013 Nobel Prize in Chemistry awarded to Martin Karplus, Michael Levitt and Arieh Warshel ‘for the development of multiscale models for complex chemical systems’.⁹ Their ideas on simplified protein structures¹⁰ are now broadly adopted in biomolecular simulation and software, e.g., in the Martini force field.¹¹

Conversely, CG mapping of solvents, and fluid phases in general, is a less clear-cut problem^12–14 apparently due to the very problem of consistent definition of a fluid element (aka ‘blob’) in statistical mechanics. Unrestricted motion of atoms in the fluid phases stipulates an open system view of CG particles¹⁵ and approaches based on Voronoi tessellation,^14,16,17 iterative k-means clustering of MacQueen,^18,19 ghost ‘blob’ microprobes,²⁰ and Brownian Quasiparticles²¹ have been put forward to tackle this matter, although the problem of an ‘inter-blob’ potential still remains largely open and no closed form solution exists so far.¹³ Crucially, although CG interactions can be cast in the form of an effective pair potential, its functional form is likely to be different from that of atomistic systems.^6,22 Accuracy of CG mapping depends critically on the functional form of the ansatz, which could only be hypothesised. This already poses a problem with bonded interactions,^6,22,23 but presents a fairly non-trivial problem for non-bonded interactions that involves genuine perception.¹³ Existing methods provide the means to optimise the parameters of the model, but not its functional form. Common ansätze validated at the atomic scale, such as the omnipresent Lennard-Jones (LJ) potential, are not necessarily the best choices for generic effective CG potentials,²⁴ where a soft Gaussian would be the likely candidate. In fact, the utility of the LJ potential has been recently questioned even at the atomic scale.²⁵ Furthermore, although it usually goes unquestioned in coarse-graining, various bottom-up approaches^26,27 hold the underlying atomic representation as a ground-truth structure. It is important to realise, however, that the atomic data are themselves based on effective empirical pair potentials, liable to be inaccurate away from their training sets.

For large-scale fluid dynamics problems dissipative particle dynamics (DPD)²⁸ is arguably the most widely used particle-based method describing various aspects of multiphase flow at the mesoscale.¹³ Being generic, fast and simple to implement, DPD expands the range of accessible time and spatial scales by several orders of magnitude.²⁹ DPD introduces a remarkably simple, albeit fundamentally challengeable ansatz: harmonic repulsion between the beads representing fluid elements. The model has been instructive in broad range of applications including colloids, polymers, fluid mixture and amphiphilic systems.¹³ However, employing only soft repulsive forces has a serious drawback – DPD covers both fluid and solid phases, but misses gas–liquid coexistence.³⁰ In other words, it is incapable of capturing such vital liquid surface phenomena as wetting, spreading, and many important capillary effects.^31,32 Attraction is a prerequisite for a stable liquid–gas interface,³³ with different requirements for systems of different dimensionality. The relative range of attraction in CG potentials is likely be reduced,³⁴ though it is believed that for 3D systems any finite-range attraction guarantees a first-order liquid–gas transition.³⁵

A currently adopted way to remedy this serious drawback in the DPD model is to add attraction to its conservative forces by including many-body interactions through an ad hoc density-dependent term, as done in many-body DPD (MB-DPD),³⁶ with a caveat of thermodynamic inconsistency.³⁷ However, explicit many-body terms are not strictly necessary for gas–liquid coexistence, and simple fluid models that include many-body terms only effectively are able to capture a wealth of capillary phenomena within the pairwise approximation.³¹ In a more straightforward approach the attraction can be included by extending the DPD model beyond the harmonic approximation and combining two terms of different powers: a positive term of higher power would give an anharmonic repulsion, and a lower-power term taken with a negative sign would provide the attraction. For example, Groot and Stoyanov³⁰ used a quadratic force in their sticky elastic sphere model to study colloidal dispersions. The standard DPD model specifies the sign of the weight function but not its functional form,²⁸ and lifting the requirement of nonnegativity one can explore this freedom in full. Thus, combining two terms with opposite signs is equivalent to a single term with a weight function that is negative within a certain domain. Furthermore, since the conservative and Langevin weight functions of the model are not generally correlated, the latter can be preserved, retaining the hydrodynamics of the original DPD. Here, we explore this approach to extend the DPD model to all three phases of matter and relate the corresponding thermodynamics at meso- and nanoscales.

In the van der Waals (vdW) picture of simple liquids, the short-range structure and correlations are defined primarily by inter-particle repulsion, while attraction plays only a minor role.³⁸ However, when it comes to thermodynamic properties, both attraction and repulsion actively shape the phase diagram,³⁹ which is liable to differ for the CG system as the result of entropy loss. Even for ‘bottom-up’ derived CG potentials for bonded species,⁴⁰ isomorphism between atomic and CG phase envelopes is not guaranteed. To establish a corresponding connection for systems with coarse-grained fluid particles reliable effective potentials are required. Although their rigorous form is not available, a clue may be found from polymer solutions by taking the view of a fluid blob as a limiting case of a star or ring polymer with the force constant between the monomers k → 0. Both theory²³ and simulation⁶ provide evidence that the effective potential between the polymers is soft, repulsive and short-ranged. Thus, it is well known that the interactions between the polymer chains in solution can be described by a Gaussian.²³ Similarly, ring polymers are coarse-grained to ultrasoft bounded potentials.⁴¹ Star polymers are also shown to be exponentially decaying outside their corona diameter.⁴² Thus one could conclude that the effective potential between the fluid blobs is also bounded and short-ranged. Softness of the CG potential for non-bonded interactions has also been demonstrated from the atomistic perspective,⁴³ and moreover, the DPD method is entirely hinged on this approach.

Simple potentials do not necessarily lead to simple thermodynamics, which could be rich in unusual and anomalous features.⁴⁴ There are several soft potentials including the Stell–Hemmer^33,35 and Jagla models⁴⁵ which are known for several anomalies.^46–48 Their anomalies are ascribed to the presence of regions with negative curvature in the potentials (‘core-softened’ potentials). On the other hand, our extension of the DPD model eschews core-softening and its analytic potential remains convex throughout its entire range, being likely the first of this type to exhibit a full set of anomalies.

The principal objective of the present paper is to examine the phase behaviour of a generic model for the coarse-grained solvent, contrasting it with that of a fully fledged atomistic solvent. Often, the origins of complex thermodynamic phenomena are subtly encoded in the details of interparticle interactions. However, a simple general model could help revealing the underlying physical mechanisms. In the next section we outline the generalised potential, termed nDPD, and provide details of simulations. We then present the results of condensed phase simulations using nDPD potentials, demonstrating stable liquid and solid phases for all three studied orders, n = 2, 3, 4. We further discuss the observed anomalies in the thermodynamic properties and outline our conclusions.

Simulation methods

DPD interactions with attraction

In classical statistical simulations an interaction between spherically-symmetric species is well described by the Mie⁴⁹ potential, which in the current notation is written aswhere ε and σ are the parameters defining the energy and length scales correspondingly, and n and m (n > m > 3) are the exponents defining together the potential well width and the steepness of the potential at short range. The case n = 12, m = 6 is the archetypal Lennard-Jones (LJ) potential,⁷ illustrated in the inset to Fig. 1. Although this generic potential is able to describe quantitatively a wide range of classical atomistic systems including their liquid phases and critical phenomena, its use in CG modelling is hindered by several factors. Firstly, it is too steep at short range⁶ and diverges at the origin. Furthermore, as a potential of infinite range, it requires truncation in numerical schemes,²⁵ and the thermodynamics of a truncated potential are sensitive to the details of truncation.⁵⁰ The last two issues have been recently discussed from the atomic perspective²⁵ but the proposed solution, based on potentials diverging at the origin, is again aimed at atomistic modelling or coarse-graining bonded potentials.

Fig. 1 Standard DPD (orange line) and proposed nDPD interaction model (green, red, and blue lines for n = 2, 3, and 4, respectively) showing the attractive region in r from σ (∼0.5r_c) up to the characteristic length scale (cutoff distance) r_c. In the inset, the standard Lennard-Jones potential is shown in the LJ reduced units, ε = σ = 1.

Unlike CG potentials for chemically bonded species, diverging at the origin, coarse-grained non-bonded interactions in complex fluids are likely to be ‘soft’, i.e., bound everywhere,^6,22,42 which is a consequence of intra-particle motion.^42,43 Perhaps the most general form of the soft ‘overlap’ potential is the Gaussian potential, which has re-appeared in classical statistical studies many times in different guises.^{23,24,42,51,52} However, its usage has one serious shortcoming: a vanishing force at the origin. A bounded potential is also at the crux of the DPD model – there, the force is at its maximum at the origin, and its softness and continuity enables fast local equilibration and exploration of phase space.

A possible way to include both repulsion and attraction in the DPD framework is based on the cubic spline weight functions of smoothed particle hydrodynamics (SPH).⁵³ Such ad hoc solutions would nevertheless provide little generality and insight, and as with the case of Gaussian potential, would lead to a vanishing force at the origin. Instead, we put forward a new potential class for CG simulations building on the DPD model as the lowest term in a series expansion, and construct the CG potential by combining two different-order terms: the higher order, n, describing the repulsion and the lower order, m, describing the attraction and taken with the negative sign. We conjecture that by augmenting the DPD potential with a short-ranged attractive part we retain the DPD advantages of accurately describing the collective dynamics of fluid phases, adding the benefit of including the rich world of capillary phenomena due to integrated positive surface tension. We start from the DPD model, where there are two types of forces between the particles (beads): the conservative (potential) force, which defines the thermodynamic state of the substance (fluid), and the Langevin (dissipative and random) forces, which are responsible for its hydrodynamic behaviour. Here, we focus on the first type, F_c(r_ij), which in the standard DPD model decreases linearly with the interparticle separation r_ij ≡ |r_ij|,

giving rise to a quadratic potentialHere, e_ij ≡ r_ij/r_ij is the unit vector in the interparticle direction and r_c is the cutoff distance beyond which both the force and potential vanish. It defines a characteristic length for the DPD model and our extension, which can be used as a basis for comparison with other interaction models, particularly those that make use of a cutoff distance primarily for computational convenience rather than as a material-defining parameter.

We generalize the conservative force in eqn (2) to arbitrary powers n, m,

where the repulsive exponent, n, represents the force order, which together with the attractive term exponent, 1 ≤ m ≤ n, and two other parameters of the model, A_ij and b_ij, completely define the interactions. A_ij is the repulsive parameter with the unit of force, and b_ij is a scaling factor defining the magnitude of repulsion relative to attractive interactions. Here, we consider only the case of m = 1, and in the following we omit m completely for brevity, denoting the potential as nDPD. The corresponding energy term is, therefore,From this expression it is clear that b_ij > 0.5(n + 1) is required for U⁽ⁿ⁾_c(r_ij) = 0 to have a root σ (0 < σ < r_c), which from now on we identify as the nDPD particle size in a similar manner to the same parameter in the Mie and LJ potentials. Also, the energy at origin, U⁽ⁿ⁾_c(0), depends on both A_ij and b_ij. Finally, setting also n = 1 and b_ij = 2 would recover the standard DPD potential. We re-emphasise that the modification affects only the conservative interactions, preserving the Langevin part (the ‘DPD thermostat’).

We restrict our consideration here to the first three low-order cases only, n = 2, 3, 4, for a single component (A ≡ A_ij, b ≡ b_ij), leaving higher powers for subsequent publication. The cases are sketched in Fig. 1, where they are contrasted with the standard DPD potential and LJ potential (inset) used in atomistic simulations. With the selected parameters given in Table 1, the location of the first zero of the potential, σ_ij, is ∼0.5r_c, decreasing somewhat for higher values of n. By contrast, the standard DPD potential decays monotonically to zero at r_c. A common peculiarity of both DPD and nDPD potentials, distinguishing them from conventional inverse-power ones, is the relationship between the force and the potential exponents: whereas in the former the exponent in the force is lower by one from the potential power, in the latter it is higher by one.

Table 1 Parameters of the nDPD potential in DPD reduced units. Particle size, σ_ij, is given as a fraction of r_c

n	A ij	b ij	σ ij /rc
2	25.0	3.02	0.5033
3	15.0	7.2	0.4730
4	10.0	15.0	0.4497

---

DPD simulations and units

All nDPD calculations were performed using a modified version of the DL_MESO mesoscopic simulation package^54,55 to accommodate the new potential type. We also implemented this potential in the DL_POLY molecular simulation package⁵⁶ and used it to independently verify the obtained pressure in a microcanonical ensemble (leaving Langevin forces out) in certain cases. We used nDPD units throughout, which are the DPD reduced units,⁷ with the following exception: we define the length scale in terms of σ, the shortest distance at which the interparticle potential reaches zero, as is established, e.g., for the LJ potential. We take it as the particle size, as shown in the last column in Table 1. When the properties are expressed in terms of the critical parameters (Table 2), they are denoted by the asterisk superscript (e.g. T* = T/T_c) and referred to as scaled reduced properties.

Table 2 Critical parameters of the model in nDPD units (see text)

n	T c	p c	ρ c
2	1.025	0.2951	0.519
3	1.283	0.3979	0.504
4	1.290	0.4095	0.484

---

Integration of the Langevin equations of motion stipulates isothermal conditions. In setting the system and during the NVT simulations we used the standard velocity Verlet integrator with the DPD thermostat,⁵⁵ and in DL_POLY calculations with the Nosé–Hoover thermostat, while for NpT calculations Langevin pistons were used for both thermostat and barostat.⁵⁷

Since nDPD interactions are smooth, bounded and of finite range, no special numerical treatment is required to deal with their calculation for large inter-particle distances as for Lennard-Jones interactions. Such treatments can involve truncation of the potential within a cutoff distance with tapering of forces and energies to zero at the cutoff to prevent energy fluctuations and estimation of long-range tail corrections for contributions to energies and pressure.⁵⁸ These approaches need special care and are prone to introducing errors in calculations of thermodynamic properties⁵⁹ that can only be eliminated by using more accurate and expensive computational methods for long-range interactions⁶⁰ to account correctly for inhomogeneous phase formation and interfacial phenomena.⁶¹

Without the need to account for long-range effects, the only controllable sources of systematic errors in the results are thus due to the finite timestep and the system size. Both factors were evaluated and in the simulations we used a dimensionless timestep of 0.005. Due to the finite time step, the system temperature estimated from the ensemble-averaged kinetic energy increases quadratically from the thermostat temperature, typically by 0.1% with the selected timestep. However, by extrapolating to an infinitesimal limit the thermostat temperature is recovered. In the fluid part of the nDPD phase diagram the pressure is dominated by the kinetic contribution, and finite time steps therefore lead to a quadratic rise in pressure, which amounted to an increase of one half percent near the critical point. To alleviate this effect we used the same timestep in all calculations, and for low-temperature coexistence simulations extrapolated pressure to zero timestep. With regard to the system size, we used between 8000 and 64 000 particles depending on the estimated property to make the error due to system size of the order of statistical uncertainty in the results or even smaller.

Simulation of orthobaric densities and the solid phase

Coexistence curves were estimated from the series of two-phase calculations performed in slab geometry (see Fig. 2) using rectangular elongated prism-shaped simulation boxes containing between 20 000 and 64 000 particles. Such a setup allows us to simultaneously estimate both the coexisting pressure and densities at a given temperature, and also the surface tension from the difference in the diagonal pressure components.⁷

Fig. 2 A sketch of the system in slab geometry used to estimate the coexisting densities. The graph presents the density profile of the DPD particles (green line) giving the values of two coexisting densities, where the profile is flat, and showing two liquid–vapour interfaces at z ≈ ±10r_c.

Closer to the critical point the coexisting densities can no longer be reliably determined using this setup due to increased fluctuations and spontaneous creation and destruction of interfaces. The coexisting pressure can still, however, be reliably estimated as the normal pressure component along the long axis in a rectangular prism of a high aspect ratio, given that the liquid/gas interfaces have insignificant curvature. This is expected in the proximity of the critical point where the correlation length rapidly grows and the surface tension decreases. We estimated therefore the coexisting densities using a set of isotherms in the subcritical region, as illustrated in Fig. 3. The isotherms were constructed from a set of isochoric calculations of cubic boxes with N = 8000 particles of progressively larger sizes. Below the critical temperature they include the van der Waals (wdW) loop as illustrated in the Figure for the T = 1.27 and 1.28 cases. Although for finite volumes all states along the isotherms are thermodynamically stable,⁶² the regions marked by the dashed lines in two lower curves may or may not phase separate inside the boxes.⁶² These regions, of course, are metastable or completely unstable in the thermodynamic limit, which precludes us from using the Maxwell construction.⁶³ Instead, we used the pressure values obtained in two-phase simulations (the horizontal dashed lines in the Figure) and estimated the coexisting densities from the crossing points with the isotherm where it has a negative slope, as marked by the vertical dashed lines in the Figure. In this way our estimates are based only on thermodynamically safe data from the regions covered by the full lines in the Figure and allowed us to obtain the densities for temperatures as close as T = 0.99T_c.

Fig. 3 Pressure as a function of volume per particle in terms of particle size, σ, for several near-critical isotherms for n = 4 potential. Horizontal dash-dotted lines denote coexisting pressure from the two-phase simulation; vertical dashed lines denote the positions of coexisting volumes; black dot marks the critical point. See text for details of estimating the coexisting densities.

We determined the type of the solid structure by annealing the systems, each containing 4N³ and 2N³ particles (where N is an integer, we used between 12 and 20) and thus promoting formation of BCC and FCC solids correspondingly, and allowing the solid phases to form. We subsequently used common neighbour analysis (CNA), as implemented in the Open Visualization Tool OVITO,⁶⁴ to analyse the particle trajectories and determine the preferred type of the cubic cell. The ratio of FCC to BCC local structures was more than 20 in all cases, from which we concluded that the equilibrium structure is FCC.

In order to determine the ground state (zero temperature) energy and structure we calculated the energy of the FCC lattice as a function of the nearest-neighbour distance r₀, the results are shown in Fig. 4. Thermodynamic stability requires the condition that the total configurational energy of a system of N particles is bounded from below,⁶⁵i.e., that U_N(r^N) ≥ −NB, where B is a fixed number. This condition is controlled by the nDPD parameter b_ij, which for a single-component system with this pairwise potential resolves to⁶⁶

and the values in Table 1 were selected to satisfy it. The limiting values of b in the 2DPD (n = 2), 3DPD (n = 3) and 4DPD (n = 4) cases are thus 3, 7 and 14 respectively.

Fig. 4 Zero-temperature energy per particle as a function of interparticle nearest-neighbour separation. (A) FCC lattice with nDPD potential n = 2 for three values of the repulsion parameter; inset – comparison of FCC and BCC energies for b = 3.02. (B) energies for n = 4 potential for four values of the repulsion parameter.

The oscillations in energy, well resolved for the n = 2 potential (Fig. 4A), result from new members of the lattice entering the range of attractive potential forces as the system compresses. The oscillations indicate an intriguing possibility of different polymorphs in the solid state, something to be looked into in the future. The inset shows the energy per particle density dependence in the FCC and BCC lattices, indicating a series of FCC ↔ BCC transitions as the pressure increases. Choosing a parameter b below the stability limit will lead to a system collapse akin to gravitational collapse in neutron stars, as seen from the cases of b = 2.98 for 2DPD (Fig. 4A) and b = 14.0 for 4DPD (Fig. 4B); it is also supported by our simulation tests.

Critical point and critical constants estimation

We determined the location of the critical point from a set of isotherms at progressively increased temperatures. The first isotherm that does not contain the vdW loop, here the T = 1.29 isotherm, is deemed to be the critical one. To estimate the critical density of the model we fit the Wegner expansion,⁶⁷ρ_l(τ) −ρ_v(τ) = Aτ^β + Bτ^β+Δ, where τ = 1 − T/T_c, to our coexistence data, with β = 0.32653 and Δ = 0.5 assuming the commonly accepted 3D Ising universality class⁶⁸ and using the fact that the lattice gas model and the Ising model are isomorphic.⁶⁹

Solid–liquid phase transition

In order to locate the solid–liquid transition control parameters (transition temperature in the T-scan and transition pressure for pressure-induced melting) we initially established the approximate area of transition by scanning the control parameters and monitoring the density (volume) for a system of N = 32 000 particles. This gave us roughly ±20% margins in the control parameters around the transition. To solve the nucleation problem, we then prepared an initial two-phase state by equilibrating separately bulk liquid and bulk solid samples using the parameters outside of the transition margins and then bringing the two systems in contact, watching for the absence of significant overlaps at the flat boundary. Using a sample of N = 64 000 particles with one half in a typical liquid configuration and the other in a typical solid one, we conducted a series of short NpT simulations at progressively higher values of the control parameter, monitoring the kinetics of the interphase boundary. The direction of the boundary motion (towards freezing or melting of the whole sample) appears to be a sensitive indicator of the ensuing phase, insensible to the initial configurations, and we were able to locate the transition within a ±1% margin for the temperature and within ±10% for the pressure.

Results

Properties along the liquid–gas coexistence curve

Using the nDPD potential with n = 2, 3, 4, we explored both branches of the binodal of their fluid phases from the triple point to the critical point. The binodal curves presented in Fig. 5 in scaled reduced units (with respect to the corresponding critical density and critical temperature, see Table 2) demonstrate several unusual features. First, the liquid phase in the nDPD model has a much wider extent than simple atomic liquids and water. For example, for the n = 4 (4DPD) potential the liquid phase extends down to T* = 0.0709 (T = 0.0915), whereas argon has its triple point at T* = 0.556, and water at T* = 0.422. Second, while the gas branch in all three cases falls on the same curve on the scale of the figure, the liquid branches are markedly different outside of the scaling region. Even the concavity of the liquid branch varies with the parameter b of the potential: with the selected parameters the liquid branch is concave for n = 4, a straight line for n = 3, and a convex curve for n = 2, challenging the law of corresponding states. We have verified that the n = 2 branch becomes concave when b increases from 3.02 to 3.1. The straight line appears to be the limiting case for n ≥ 3. In the scaling region near the critical point all liquid branches collapse onto a single universal curve.

Fig. 5 Coexistence curves for the nDPD potential, calculated for quadratic, cubic, and quartic repulsive forces with parameters given in Table 1. The blue line denotes the gas branch (common for all powers) and the green, orange, and red lines are for the liquid branches. The black dot denotes the critical point (T* = 1). Reduced units have been used with critical values defining the scales. (Insets) Radial distribution functions at the coexistence densities: left – gas branch, right – liquid branch, at T* = 0.4 for the three potential classes.

The insets to the Figure illustrate the radial distribution functions (RDF) for both branches at a representative temperature. For the gas branch all cases display a typical low-density atomic (e.g., LJ fluid) RDF, with its height increasing with the potential well depth: more so for lower powers of potential (n). Conversely, the liquid RDFs show a transition from the typical for simple fluids liquid-density RDFs for the two higher powers to an unusual case for the 2DPD potential with the second minimum more pronounced than the first, where the positions of the first and the third solvation shells are also shifted to shorter distances but the position of the second one remains unchanged. We have found these features of the liquid-density RDF to be sensitive to the value of b – they disappear when b is increased from 3.02 to 3.10 – and are not unique to 2DPD: the second minimum becomes similarly pronounced when b < 7 for 3DPD and b < 14.5 for 4DPD.

We calculated the surface tension γ as a function of temperature and, by fitting the data to a power law γ ∝ (T_c − T)^μ, we obtained an independent estimate of the critical temperature. For the 4DPD potential, illustrated in Fig. 6A, T_c = 1.28(2), in excellent accordance with the value obtained from the binodal data (T_c = 1.29). The second parameter of the fit, the critical exponent μ, has a converged value of 1.38, deviating somewhat from the accepted value for fluids, μ = 1.28(6).⁷⁰ In the Figure we expressed the surface tension in the standard DPD units (length in units of potential range a). This facilitates the comparison with experimental data for key solvents. For example, water modelled as a DPD fluid, with an a = 8.52 Å (7 : 1) mapping would have a surface tension of γ = 12.7 at T* = 0.46 (based on 72.0 mN m⁻¹ at T = 298 K⁷¹).

Fig. 6 (A) Calculated surface tension of 4DPD liquid in DPD units as a function of temperature up to the critical point. Solid line is a power law fit (see text). (B) Temperature dependence of density around the solid–liquid phase transition. Red diamonds – liquid phase, blue diamonds – solid phase, and the temperature of transition is marked by the green dashed line. The large blue dot indicate the position of the temperature of maximum density (TMD).

Contraction on fusion, temperature of maximum density, and negative thermal expansion in the solid

At low temperatures, an nDPD fluid freezes for all three studied classes, which is where it exposes another unusual behaviour. Whereas most substances expand when heated, there is a compelling group of materials that contract upon heating within a certain temperature range. Classic textbook examples include liquid water below 4 °C (277 K) and ice below 63 K.⁷² This phenomenon, termed Negative Thermal Expansion (NTE),⁷³ occurs when materials contract in volume upon heating and is characterised by the negative value of the bulk thermal expansion coefficient, α = (∂ ln V/∂T)_p < 0. NTE in complex materials with directional bonding, where it is more often found,^73,74 received much attention in recent years due to its industrial applications. Molecular mechanisms explaining NTE in these systems include, inter alia, Grüneisen vibrational theory of thermal expansion,^73,75 special phonon properties,⁷⁴ and rigid unit modes.⁷⁴ At the same time, there are a few pure elements where the NTE and contraction of fusion have been found, neither of which can be explained by the above mechanisms. These elements include the semi-metal bismuth, the metals gallium and plutonium, and the metalloids antimony, germanium, and silicon. They crystallise into different crystal structures,⁷⁶ belong to different groups in the periodic table, and exhibit markedly different NTE ranges: it is thus unlikely that their behaviour can be attributed to unique electronic properties alone. There is no clear understanding of the underlying mechanisms in these systems. Furthermore, several classical single-component model systems with isotropic pair potentials, such as the purely repulsive Gaussian core model,^51,77 as well as the ad hoc devised potentials with a softened interior of their basins of attraction,⁷⁸ also feature NTE. However, finding a simple model with an isotropic potential, which reproduces the anomalous behaviour of water—contraction on melting with an immediate negative expansion region in the liquid phase—is known to be a challenging open problem.⁷⁸ Here, we demonstrate that the nDPD model includes all of the essential ingredients required to capture this complex phenomenon.

In order to locate the solid–liquid phase transition we calculated the isobaric density of nDPD at low temperatures. The results for the n = 4 potential at zero external pressure are presented in Fig. 6B. The freezing transition occurs at T = 0.0779 (T* = 0.0760) for the 2DPD potential, at T = 0.0873 (T* = 0.0680) for the 3DPD potential, and at T = 0.0915 (T* = 0.0709) for the 4DPD potential. The system freezes into an FCC solid with an expansion in volume of ca. 5.5% in all three cases. Just before freezing a density maximum is observed at T = 0.11 (T* = 0.0852) for the 4DPD case. On the liquid side, it qualitatively describes water, which has a temperature of maximum density (TMD) of 4 °C (277 K) or, in reduced units, T* ≈ 0.4283. Even taking into account that water expands by 9% upon freezing, nDPD presents alternatives to hydrogen bonding mechanisms that lead to the qualitatively same effect, which could be explored further when developing accurate CG water models.

Pressure-induced melting

If the liquid phase of a substance is more dense than its solid phase, the melting temperature decreases with increased pressure. This general result, which follows from the Clausius-Clapeyron equation, is well known for water, where a regelation occurs for up to 20 K below the normal freezing point. A similar phenomenon is also observed for the chemical elements mentioned in the previous section that feature NTE, as well as in alkali metals. Using the two-phase setup to study the solid–liquid transitions, we observed a monotonous decrease in transition temperature from T* = 0.0709 at effectively zero pressure down to T* ≈ 0.03 at p = 10 for a 4DPD potential. As a consequence of this, there is a maximum temperature for which the solid phase exists for the nDPD model: T* = 0.0709 for 4DPD. For temperatures below this maximum value, the system melts when pressurised. It should be noted that other models also feature pressure-induced melting. For example, in the Gaussian core model⁵¹ reentrant melting has been reported.⁷⁹

Parameterisation of nDPD and thermodynamic consistency

In order to relate thermodynamic properties at CG and atomistic levels a mapping between corresponding units is needed. At equilibrium, this includes at least energy and length scales. In the DPD model the link to a specific chemistry requires a mean-field approximation.²⁸ With only one parameter in the model, A_ij, it is sufficient to fix it by matching the isothermal compressibility of DPD fluid and that of a given material, e.g. water.²⁸ This is justified by the fundamental connection between DPD and fluctuating Landau–Lifshitz Navier–Stokes hydrodynamics.⁸⁰

The elastic responses of a nDPD liquid are richer and also include tensile stresses with the tensile strength controlled by the parameter b_ij. In addition, the liquid–gas coexistence present in the nDPD model with its associated critical point can fix the temperature scale. Thus, A_ij can be fixed from the condition T_c = 1, and b_ij by mapping the surface tension of the nDPD model and the target fluid. The value of b_ij is limited by the need for thermodynamic stability,^65,66 which is satisfied for a given value of n by eqn (6).

The length scale for the nDPD model can be fixed by relating liquid coexisting densities between simulation and physical representations. Although we did not aim for a specific system here, a relation to ambient water can be made as follows. The length scale in DPD is usually fixed by mapping 3 water molecules to one bead.²⁹ There is no compelling reason for this choice, and due to tetrahedral ordering of liquid water at ambient conditions a 5 : 1 mapping, i.e., a bead representing a water molecule together with its first coordination shell, seems equally sensible. Since the liquid density of water is 1/30 Å⁻³,²⁹ and the corresponding 4DPD density is 1.62σ⁻³, both taken at 0.4608T_c (298 K), the length scale using this mapping is σ = 6.24 Å, and the cutoff distance r_c = 13.88 Å.

With these parameters we can link the liquid phases of water and nDPD. Most conveniently, this could be done using a linear transformation T′ = (T − T₀)/(T_c − T₀), where the scaled temperature T′ is set to zero at the zero-pressure melting point (triple point for water), T₀, and increases linearly with T reaching a value of 1 at the critical temperature T_c. We have found that a 4DPD liquid with a b parameter in the range 15.0 ≤ b ≤ 16.5 accurately describes the saturation pressure of water across a wide range of temperatures, as illustrated in Fig. 7 with a comparison between our results and the reference equation of state given by the International Association for the Properties of Water and Steam (IAPWS).⁸¹

Fig. 7 Saturation pressure in scaled reduced units of a 4DPD liquid with three values of the b parameter (symbols), and the water using the reference IAPWS equation of state⁸¹ (full line). The standard error in our data is of the order of symbol size. In the inset, the surface tension for these three cases is compared with the IAPWS reference. The temperature scales from the triple point for IAPWS or zero pressure melting point for 4DPD (T′ = 0), to the critical temperature (T′ = 1).

The surface tension of the nDPD liquid is obtained as a by-product from the difference between normal and lateral pressure components in a coexistence simulation using a slab geometry.⁷ Since the nDPD potential is of finite range, no long-range corrections are required. An inset to Fig. 7 shows the temperature dependence of the surface tension of the 4DPD case for three values of b. The surface tension, given in the Figure in DPD scaled units (T_c = 1), markedly depends on b, which provides a suitable method of tuning it to a specific target. Although the surface tension is not reproduced over the entire temperature range by a single set of parameters, it can be tuned around a particular temperature. Thus, b = 15.0 seems a reasonable choice for ambient water, T′ = 0.0668 (298 K).

Discussion and conclusions

In his milestone paper discussing the hierarchical structure of science, ‘More is different’, Phillip Anderson argued that behaviour of complex systems cannot be entirely ‘understood in terms of a simple extrapolation of the properties of a few particles’, and that ‘at each level of complexity entirely new properties appear’.⁸² As we dial down the magnification moving to larger scales, we would find the new phenomena in focus, and these phenomena might or might not be present at finer resolutions. Liquid structure is defined by correlations between the density fluctuations, which includes both short-range order and long-range hydrodynamic fluctuations. However, short-range and long-range hydrodynamic fluctuations are coupled, and coarse-graining the short-range interactions only will inevitably disrupt the balance. In designing a coarse-grained model there are a couple of specific points that needs to be taken into account.

It is well-known that coarse-graining of any set of degrees of freedom results in changes to the governing equations of motion from the Newtonian form to generalised Langevin dynamics (GLE),⁸³i.e., to the dynamics of interacting Brownian particles. Thus, the structure and dynamics of a CG fluid are coupled: a part of the (conservative) chemically-specific interatomic interactions, not captured by the inter-particle interactions at the CG level, is reinstated in the form of non-conservative generic interactions with the ‘bath’, i.e., in the form of viscous and stochastic forces coupled via the fluctuation–dissipation theorem.⁸⁴ We note in passing that, strictly, even atomic-level dynamics are also of the GLE type since empirical interatomic potentials effectively include averaged-out electronic degrees of freedom. While the neglected forces are likely to be insignificant at the atomic level, this is not so at the mesoscopic scale. Incorporating memory effects is a challenging problem and is the focus of current research,^85,86 although no viable scheme for integration of GLE exists as yet. Existing models are usually based on a Markovian (memoryless) approximation, which enormously simplifies the equations of motion but requires effective (mean-field) friction coefficients though obtainable, as has been demonstrated recently,²¹ by direct ensemble averaging of an atomistic simulation. In the CG process it is the short-range correlations that are suffer most.

Another moot point is related to generating adequate representations of many-body interactions using pairwise CG potentials. In other words, whether it is possible to retain long-range correlations in liquid with depleted short-range ones, since they are related via the Ornstein-Zernike (OZ) equation?⁸⁷ In general, any adjustment in the Hamiltonian alters the thermodynamic phase behaviour.⁵⁰ This is true even if the CG Hamiltonian is obtained using the rigorous Mori-Zwanzig projection operator formalism⁸³—the intrinsic reduction in entropy will affect the thermodynamic functions of the CG system.⁸⁸ With the reduced set of effective degrees of freedom certain correlations are lost and, with them, all related phenomena can be missed. This problem of representability, which has been often skirted in previous CG studies, is now receiving much attention.^40,88,89 The development of rigorous bottom-up approaches with chemical accuracy is further hampered by the limited transferability of CG models and economic factors, which need to be taken into account: honing the CG model to specific chemistry may require more effort than solving the problem in question at the atomistic level. At the same time, there is always great demand for finding models that would provide greatly simplified points of view, and as such there is a continuous quest for simple generic models with extended representability. Simple intermolecular interactions do not necessarily entail simple thermodynamics,⁴⁴ and analytically simple soft potentials may conceal unexpected macroscopic phenomena.^34,45,51,90 In order to shed more light on the impact of CG of a system on its thermodynamics we study the phase diagram of a simple soft pair potential that is a generalization of the standard potential used in dissipative particle dynamics (DPD).

Despite the wide popularity of the DPD model,¹³ little justification has been provided for its harmonic repulsion potential.²⁸ It might therefore be surprising that a simple repulsive force of the DPD model can provide a working solution for many cases.¹³ To qualitatively understand this success, imagine a coarse-graining process with progressively larger scaling factors. The liquid structure is defined in terms of density correlations, which includes—according to the OZ approach⁸⁷—a short-ranged component referred to as a direct correlation function. In the OZ picture the total correlation between two particles is composed of the direct correlation and an indirect part, which is longer-ranged and composed of all possible particle chains of direct correlation in the fluid.⁸⁷ In the CG process the direct correlations are averaged out, which leaves just indirect correlations to be reproduced by a suitable ansatz. At this level of CG the interactions are much softer, and we argue that they are captured by the nDPD model. Since much of the chemical specificity is absorbed within the short range, the longer-range correlations are of a more generic nature. At even higher CG scaling most of the correlation is averaged out, leaving only the soft elastic repulsion between the beads due to the finite compressibility of the fluid. The interactions emerging at this scale are therefore of the DPD type.

As we have demonstrated, the nDPD potential with appropriate choices for its three parameters is able to produce stable liquid and solid phases. At the same time, the phase diagram of nDPD matter is topologically different from that of simple atomic systems, with the solid phase occupying a narrow density region at low temperatures. The thermodynamics of its condensed phases are rich with anomalous features, usually attributed to complex molecular systems. They include expansion on freezing, temperature of maximum density, negative thermal expansion in the solid phase, an unusual liquid–vapour phase envelope and pressure-induced melting. In this work, we have not explored in full the model parameter space, particularly the relative range of the attraction, which is likely to conceal additional anomalies. However, the fact that the anomalies are found in systems with soft coarse-grained isotropic pair potentials indicates that they likely to originate in the medium-range structure, i.e., beyond the first neighbours in the underlying atomistic system. The complex phase behaviours have also been found in models with coarse-grained bonded interactions^34,44 and were observed in colloidal suspensions.⁹¹

Clearly, no CG model can reproduce the whole spectrum of observables of the underlying atomistic model, in the same way as no atomistic model can reproduce the entire set of experimental observables. This is a manifestation of the representability problem,⁸⁹ which arises due to the resolution change at the mean-field level: certain correlations are missing. Thus, a CG model cannot describe the processes of self-assembly or gelation, which involve changes in fine-level molecular topology. However, this is not the purpose of the model. Rather, each model level in the hierarchy defines its own classes of compatible observables that reproduce the experimental observables within the model constraints.⁸⁹ A good model goes beyond this and has predictive power. Paraphrasing the famous designer Yves Saint-Laurent, one can say that a good model can advance the field by ten years.⁹²

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by the STFC's Innovation and Development Programme (IDP) and made use of ARCHER2 computational support by CoSeC, the Computational Science Centre for Research Communities, through the High-End Computing Materials Chemistry Consortium (MCC). Additional computing resources were provided by STFC Scientific Computing Department's SCARF cluster. Contributing funding was also provided via the Virtual Materials Marketplace (VIMMP) project from the European Union's Horizon 2020 Research and Innovation Programme (Grant Agreement no. 760907).

Notes and references

1. W. Noid, G. S. Ayton, S. Izvekov and G. A. Voth, Coarse-graining of condensed phase and biomolecular systems, CRC press, Boca Raton, 2008, pp. 21–40.
2. S. R. Nagel, Rev. Mod. Phys., 2017, 89, 025002.
3. C. Peter and K. Kremer, Soft Matter, 2009, 5, 4357–4366.
4. M. Sadeghi and F. Noé, Nat. Commun., 2020, 11, 2951.
5. L. I. Vazquez-Salazar, M. Selle, A. H. de Vries, S. J. Marrink and P. C. T. Souza, Green Chem., 2020, 22, 7376–7386.
6. C. N. Likos, Phys. Rep., 2001, 348, 267–439.
7. M. P. Allen and D. J. Tildesley, Computer simulation of liquids, Oxford University Press, Oxford, 2017.
8. M. G. Saunders and G. A. Voth, Annu. Rev. Biophys., 2013, 42, 73–93.
9. M. Karplus, Angew. Chem., Int. Ed., 2014, 53, 9992–10005.
10. M. Levitt and A. Warshel, Nature, 1975, 253, 694–698.
11. P. C. T. Souza, R. Alessandri, J. Barnoud, S. Thallmair, I. Faustino, F. Grünewald, I. Patmanidis, H. Abdizadeh, B. M. H. Bruininks, T. A. Wassenaar, P. C. Kroon, J. Melcr, V. Nieto, V. Corradi, H. M. Khan, J. Domanski, M. Javanainen, H. Martinez-Seara, N. Reuter, R. B. Best, I. Vattulainen, L. Monticelli, X. Periole, D. P. Tieleman, A. H. de Vries and S. J. Marrink, Nat. Methods, 2021, 18, 382–388.
12. A. Eriksson, M. Nilsson Jacobi, J. Nyström and K. Tunstrøm, Europhys. Lett., 2009, 86, 44001.
13. P. Español and P. B. Warren, J. Chem. Phys., 2017, 146, 150901.
14. Y. Han, J. F. Dama and G. A. Voth, J. Chem. Phys., 2018, 149, 044104.
15. G. Ciccotti and L. Delle Site, Soft Matter, 2019, 15, 2114–2124.
16. M. Serrano and P. Español, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2001, 64, 046115.
17. E. G. Flekkøy and P. V. Coveney, Phys. Rev. Lett., 1999, 83, 1775–1778.
18. J. MacQueen, Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability, 1967, pp. 281–297.
19. K. R. Hadley and C. McCabe, J. Phys. Chem. B, 2010, 114, 4590–4599.
20. G. S. Ayton, H. L. Tepper, D. T. Mirijanian and G. A. Voth, J. Chem. Phys., 2004, 120, 4074–4088.
21. V. P. Sokhan and I. T. Todorov, Mol. Simul., 2021, 47, 248–256.
22. J.-P. Hansen and H. Löwen, Bridging Time Scales: Molecular Simulations for the Next Decade, Springer-Verlag, Berlin, 2002, pp. 167–196.
23. P. J. Flory and W. R. Krigbaum, J. Chem. Phys., 1950, 18, 1086–1094.
24. B. J. Berne and P. Pechukas, J. Chem. Phys., 1972, 56, 4213–4216.
25. X. Wang, S. Ramrez-Hinestrosa, J. Dobnikar and D. Frenkel, Phys. Chem. Chem. Phys., 2020, 22, 10624–10633.
26. S. Izvekov and G. A. Voth, J. Phys. Chem. B, 2005, 109, 2469–2473.
27. M. S. Shell, J. Chem. Phys., 2008, 129, 144108.
28. R. D. Groot and P. B. Warren, J. Chem. Phys., 1997, 107, 4423–4435.
29. R. Groot and K. Rabone, Biophys. J., 2001, 81, 725–736.
30. R. D. Groot and S. D. Stoyanov, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 78, 051403.
31. J. S. Rowlinson and B. Widom, Molecular theory of capillarity, Clarendon Press, Oxford, 1982.
32. A. Ghoufi, P. Malfreyt and D. J. Tildesley, Chem. Soc. Rev., 2016, 45, 1387–1409.
33. P. C. Hemmer and J. L. Lebowitz, Phase transitions and critical phenomena, Academic Press, London, 1976, vol. 5b, pp. 107–203.
34. C. F. Tejero, A. Daanoun, H. N. W. Lekkerkerker and M. Baus, Phys. Rev. Lett., 1994, 73, 752–755.
35. P. C. Hemmer and G. Stell, Phys. Rev. Lett., 1970, 24, 1284–1287.
36. P. B. Warren, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 68, 066702.
37. A. A. Louis, J. Phys.: Condens. Matter, 2002, 14, 9187–9206.
38. J. D. Weeks, D. Chandler and H. C. Andersen, J. Chem. Phys., 1971, 54, 5237–5247.
39. B. Smit and D. Frenkel, J. Chem. Phys., 1991, 94, 5663–5668.
40. W. G. Noid, J. Chem. Phys., 2013, 139, 090901.
41. A. Narros, A. J. Moreno and C. N. Likos, Soft Matter, 2010, 6, 2435–2441.
42. C. N. Likos, H. Löwen, M. Watzlawek, B. Abbas, O. Jucknischke, J. Allgaier and D. Richter, Phys. Rev. Lett., 1998, 80, 4450–4453.
43. S. H. L. Klapp, D. J. Diestler and M. Schoen, J. Phys.: Condens. Matter, 2004, 16, 7331–7352.
44. G. Malescio, J. Phys.: Condens. Matter, 2007, 19, 073101.
45. E. A. Jagla, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 1998, 58, 1478–1486.
46. M. R. Sadr-Lahijany, A. Scala, S. V. Buldyrev and H. E. Stanley, Phys. Rev. Lett., 1998, 81, 4895–4898.
47. N. B. Wilding and J. E. Magee, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2002, 66, 031509.
48. A. Barros de Oliveira, P. A. Netz, T. Colla and M. C. Barbosa, J. Chem. Phys., 2006, 124, 084505.
49. G. Mie, Ann. Phys., 1903, 316, 657–697.
50. T. Coussaert and M. Baus, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 1995, 52, 862–869.
51. F. H. Stillinger, J. Chem. Phys., 1976, 65, 3968–3974.
52. J. A. Purton and W. Smith, Mol. Simul., 2010, 36, 796–800.
53. M. Liu, P. Meakin and H. Huang, Phys. Fluids, 2006, 18, 017101.
54. M. A. Seaton, R. L. Anderson, S. Metz and W. Smith, Mol. Simul., 2013, 39, 796–821.
55. M. A. Seaton, Mol. Simul., 2021, 47, 228–247.
56. I. T. Todorov, W. Smith, K. Trachenko and M. T. Dove, J. Mater. Chem., 2006, 16, 1911–1918.
57. A. F. Jakobsen, J. Chem. Phys., 2005, 122, 124901.
58. S. V. Lishchuk and J. Fischer, J. Chem. Phys., 2018, 149, 091102.
59. D. Heyes, Comput. Methods Sci. Technol., 2015, 21, 169–179.
60. P. Ewald, Ann. Phys., 1921, 253–287.
61. J. B. Klauda, X. Wu, R. W. Pastor and B. R. Brooks, J. Phys. Chem. B, 2007, 111, 4393–4400.
62. K. Binder, B. J. Block, P. Virnau and A. Tröster, Am. J. Phys., 2012, 80, 1099–1109.
63. J. C. Maxwell, Nature, 1875, 11, 374–377.
64. A. Stukowski, Model. Simul. Mater. Sci. Eng., 2009, 18, 015012.
65. M. E. Fisher and D. Ruelle, J. Math. Phys., 1966, 7, 260–270.
66. A. A. Louis, P. G. Bolhuis and J. P. Hansen, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2000, 62, 7961–7972.
67. F. J. Wegner, Phys. Rev. B: Condens. Matter Mater. Phys., 1972, 5, 4529–4536.
68. F. J. Wegner, Phase transitions and critical phenomena, Academic Press, NY, 1976, vol. 6, pp. 7–124.
69. T. D. Lee and C. N. Yang, Phys. Rev., 1952, 87, 410–419.
70. B. Widom, Phase transitions and critical phenomena, Academic Press, NY, 1976, vol. 2, pp. 79–100.
71. N. R. Pallas and Y. Harrison, Colloids Surf., 1990, 43, 169–194.
72. G. Dantl, Z. Physik, 1962, 166, 115–118.
73. M. T. Dove and H. Fang, Rep. Prog. Phys., 2016, 79, 066503.
74. A. K. A. Pryde, K. D. Hammonds, M. T. Dove, V. Heine, J. D. Gale and M. C. Warren, J. Phys.: Condens. Matter, 1996, 8, 10973–10982.
75. G. D. Barrera, J. A. O. Bruno, T. H. K. Barron and N. L. Allan, J. Phys.: Condens. Matter, 2005, 17, R217–R252.
76. D. A. Young, Phase diagrams of the elements, Univ of California Press, Berkeley, 1991.
77. F. H. Stillinger and D. K. Stillinger, Phys. A, 1997, 244, 358–369.
78. M. C. Rechtsman, F. H. Stillinger and S. Torquato, J. Phys. Chem. A, 2007, 111, 12816–12821.
79. P. V. Giaquinta and F. Saija, ChemPhysChem, 2005, 6, 1768–1771.
80. P. Español and M. Revenga, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 67, 026705.
81. W. Wagner and A. Pruß, J. Phys. Chem. Ref. Data, 2002, 31, 387–535.
82. P. W. Anderson, Science, 1972, 177, 393–396.
83. R. Zwanzig, Nonequilibrium statistical mechanics, Oxford University Press, Oxford, 2001.
84. N. G. Van Kampen, Stochastic processes in physics and chemistry, Elsevier, Amsterdam, The Netherlands, 2007.
85. C. Hijón, P. Español, E. Vanden-Eijnden and R. Delgado-Buscalioni, Faraday Discuss., 2010, 144, 301–322.
86. V. Klippenstein, M. Tripathy, G. Jung, F. Schmid and N. F. A. van der Vegt, J. Phys. Chem. B, 2021, 125, 4931–4954.
87. J. S. Rowlinson and F. L. Swinton, Liquids and Liquid Mixtures, Butterworths, London, 1982.
88. N. J. H. Dunn, T. T. Foley and W. G. Noid, Acc. Chem. Res., 2016, 49, 2832–2840.
89. J. W. Wagner, J. F. Dama, A. E. P. Durumeric and G. A. Voth, J. Chem. Phys., 2016, 145, 044108.
90. S. Prestipino, C. Speranza, G. Malescio and P. V. Giaquinta, J. Chem. Phys., 2014, 140, 084906.
91. P. N. Pusey and W. van Megen, Nature, 1986, 320, 340–342.
92. Ritz Newspaper, 1984.

---