{'fraction_non_alphanumeric': 0.03196718532831622, 'fraction_numerical': 0.014068273500796638, 'mean_word_length': 4.5765595463138, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Phase behaviour of attractive and repulsive ramp fluids: integral equation and computer simulation studies

Using computer simulations and a thermodynamically self consistent integral equation we investigate the phase behaviour and thermodynamic anomalies of a fluid composed of spherical particles interacting via a two-scale ramp potential (a hard core plus a repulsive and an attractive ramp) and the corresponding purely repulsive model. Both simulation and integral equation results predict a liquid-liquid de-mixing when attractive forces are present, in addition to a gas-liquid transition. Furthermore, a fluid-solid transition emerges in the neighbourhood of the liquid-liquid transition region, leading to a phase diagram with a somewhat complicated topology. This solidification at moderate densities is also present in the repulsive ramp fluid, thus preventing fluid-fluid separation.


I. INTRODUCTION
The existence of a liquid-liquid (LL) equilibrium and density, diffusivity and structural anomalies in single component fluids has attracted considerable interest in the last decade.
The density and diffusion anomalies present in liquid water (i.e. the existence of a maximum in the density within the liquid phase, and an increase in diffusion upon compression, up to a certain point) have long since been known, and have been reproduced by several of the existing interaction models 1,2,3,4,5 . Experimental evidence for LL coexistence has been found for phosphorous 6 , triphenyl phosphite 7 , and n-butanol 8 , and it has been suggested that this LL equilibrium might be the source of the anomalies encountered in water. Computer simulations predict the existence of LL equilibria not only in water 1,2,4 but also in other loosely coordinated fluids, such as silicon 9 , carbon 10 and silica 11 . Whether these transitions physically exist or not is still open to debate, since in most cases they correspond to supercooled states which are rendered experimentally inaccessible by crystallisation. However, closely related first order transitions between low and high density amorphous phases have indeed been found for water 12 , silica 13 and germanium oxide 14 .
While the use of realistic models can provide reasonable explanations for the experimental behaviour of a physical system, a more thorough description of the mechanisms underlying both LL phase transitions and density, diffusivity and related anomalies can be acquired via the study of simplified models. In the case of LL equilibrium and polyamorphism in molecular fluids 15 the simple model of Roberts and Debenedetti 16 has successfully accounted for the behaviour of network forming fluids 17 . On the other hand, since the pioneering work of Hemmer and Stell 18,19 it is known that a simple spherically symmetric potential in which the repulsive interaction has been softened (in this case by the addition of a repulsive ramp to the hard core) can lead to the existence of a second (LL) critical point, as long as a first (liquid-vapour, LV) critical point existed due to the presence of a long range attractive component in the interaction potential. As well as the ramp potential, other simple models with two distinct ranges of interaction, such as the hard-sphere square shoulder-square well potential studied by Skibinsky et al. 20 , have also been shown to exhibit LL equilibria. Indeed the presence of two interaction ranges explains the competition between two locally preferred structures (LPS) -a LPS being defined as an arrangement of particles which, for a given state point, minimises some local Helmholtz energy 21 . This competition between two LPS helps to rationalise the existence of polyamorphism and LL equilibria in single component glassy systems and fluids 21 .
More recently, the original model proposed by Hemmer and Stell has regained attention, especially since Jagla 22 stressed the similarities between the behaviour of the Hemmer-Stell ramp potential and the anomalous properties of liquid water. This has been further explored by Xu et al. 3,23 who analysed the relationship between the LL transition and changes in the dynamic behaviour of fluids interacting via a soft core ramp potential with attractive dispersive interactions added. Gibson and Wilding 24 have recently presented an exhaustive study of a series of ramp potentials exhibiting LL transitions and density anomalies, whose relative position and stability with respect to freezing might be tuned by judicious changes in the interaction parameter. Moreover, Caballero and Puertas 25 have also focused on the relation between the density anomaly and the LL transition for this model system by means of a Monte Carlo based perturbation approach. The aforementioned authors find that, in this case, the density anomaly is absent when the range of the attractive interaction is sufficiently small.
In this work we shall refer to this interaction potential as the 'attractive two-scale ramp potential' (A2SRP). Additionally, when one considers the system without any attractive contribution, i.e. a hard-sphere core plus a repulsive ramp, as the 'two-scale ramp potential' (2SRP). Although it has been found that the system exhibits density anomalies, it is a possibility that LL transition is preempted by crystallisation 26 . In this purely repulsive model, the relation between static and dynamic anomalies has been explored in detail 26,27,28 (similar results have been found for a dumbbell fluid with repulsive ramp site-site interactions 29 ) and the connection between these anomalies and structural order has been investigated by Yan et al. 30 with the aid of the Errington-Debenedetti order map 31,32 . These authors found that there is a region of structural anomalies (in which both translational and orientational order decrease as density is increased) that encapsulates the diffusivity and density anomaly region. The intimate relationship between transport coefficients, and the excess entropy was made clear first by Rosenfeld 33,34 , and again later, specifically for the atomic diffusion, by Dzugutov 35 . Needless to say, a corollary of this relation is that any anomalous diffusion will be accompanied by an anomalous excess entropy. Recently this link has been shown to hold for both liquid silica and for the 2SRP model by Sharma et al. 27 , and for the discontinuous core-softened model by Errington et al. 28 .
The principal objective of this work is to extend our knowledge of the phase behaviour of systems interacting via either attractive or purely repulsive two scale ramp potentials (A2SRP and 2SRP). To that purpose exhaustive Monte Carlo calculations have been performed in order to determine the phase boundaries of the gas, liquid and solid phases for both model systems up to moderate densities -slightly beyond the high density branch of the LL equilibrium. Self-consistent integral equation calculations performed on the A2SRP model complement the Monte Carlo results and agree qualitatively as to the location of the LL equilibria and quantitatively for the LV equilibria. The location of Widom's line (the locus of heat capacity maxima) is obtained for both models. The temperature of maximum density (TMD) curve is obtained for the repulsive 2SRP model and correlated with the location of Widom's line. The rich variety of phases present in these simple models will be illustrated in the calculated phase diagrams.
The structure of the paper is as follows; the model and computational procedures used herein are introduced in Section II. In Section III the most significant results are presented and discussed. Finally, the main conclusions and future prospects can be found in Section IV.

II. MODEL AND COMPUTATIONAL METHODS
The first model system consists of hard spheres of diameter σ, with a repulsive soft core and an attractive region. The interaction potential reads In Figure 1 both the attractive and purely repulsive interactions are illustrated.

A. The self consistent integral equation approach
The integral equation calculations are based on the Ornstein-Zernike relation, which for simple fluids reads 37 where h(r) is the total correlation function (related to the pair distribution function, g, by h = g−1) and c(r) is the direct correlation function. This equation requires a supplementary relation, whose general form is 38 with β = (k B T ) −1 . In Eq.(4) the bridge function, B(r), is a diagrammatic sum of convolutions of h(r), and must be approximated. The simplest instance is the hyper-netted chain with the interpolating function f and Consistency could be implemented on a local level by equating Alternatively, one might use a global consistency. In this case, in the absence of spinodal lines, one can simply enforce the condition In the presence of spinodal lines one might use a mixed thermodynamic integration path, as suggested by Caccamo et al. 44 . In this paper extensive use is made of the local consistency (LC) approach. We have also explored the use of global consistently (GC) as an alternative to overcome the deficiencies exhibited by the LC approach. From a practical point of view, in order to calculate the derivative of the virial pressure in Eq.(8), we shall assume that the parameter α is locally independent of the density, as is customary. For density independent potentials this leads to minor discrepancies when comparing the integrated pressures in Eq.
(9), so both LC and GC approaches yield similar results 45 . Here the situation is somewhat different, as will be seen in the proceeding Section. Finally, the right hand side of Eq. (8) is where the density derivative of the pair distribution function is obtained by iterative solution of the integral equation that results from differentiation of (3) and (5) with respect to density. This procedure, first proposed by Belloni 46 is substantially more efficient than the classical finite differences approach when it comes to evaluating the derivative of βP v . In addition to the density derivatives of the pair correlation function we have also evaluated the corresponding temperature derivatives in order to calculate the heat capacity 47 where ω = h − c + B. Note that the integral equations in terms of the derivatives of the correlation functions (either with respect to ρ or β) are extremely well conditioned and are rapidly solved, even when starting from ideal gas initial estimates.

B. Simulation algorithms
Initially the principal aim was to establish whether the A2SRP model does indeed exhibit a fluid-fluid equilibrium. Therefore, in order to compute both the liquid-vapour and the LL equilibria, a procedure 48 based on the algorithm of Wang and Landau 49,50 was implemented.
Subsequently, it was necessary to compute the equilibrium of the fluid phases with the low density solid phase. To that end use was made of thermodynamic integration 51 and the socalled Gibbs-Duhem integration techniques 51,52,53 , which took advantage of the behaviour of the repulsive model in the low temperature limit. The use of two different methods in independent calculations provides a check on the quality of the results.

Wang-Landau-type algorithm
A more detailed explanation of the Wang-Landau procedure can be found elsewhere 48 .
Here, we shall restrict ourselves to outlining its most salient features. The algorithm fixes the volume and temperature and then samples over a range of densities (N) of the system. For T * c = 1.487 ± 0.003, ρ c σ 3 = 0.103 ± 0.001, and p c σ 3 /|W a | ≡ p * c ≃ 0.042. At low temperatures the LL equilibrium was found, but as the temperature was further reduced we were unable to obtain convergence in the estimates of the Helmholtz energy function, and the systems eventually formed a low density solid phase. Given this situation it was surmised that a triple point, if it exists, between the vapour and the two liquid phases, could not be stable, and the presence of a stable low density solid phase is most likely. After performing simulations at several temperatures, with different system sizes (0 ≤ ρσ 3 ≤ 0.55, L * = 6, 8, 10), it was found that the LL equilibrium appeared for temperatures below T * ≃ 0.38. Using the results for L * = 10, we can estimate the critical point for the LL equilibrium

Thermodynamic Integration: 2SRP model
In order to compute the phase equilibria between the low density solid and the fluid phases, it is useful to firstly consider the repulsive model. Jagla 22 computed the stability of different crystalline structures at T = 0 for the repulsive model supplemented by a mean field attraction term, for a ramp model with d a /σ = 1.75, and found that the stable crystalline phases at low pressure correspond to the face-centred cubic (fcc) and hexagonal close packed (hcp) structures. On the other hand, in the limit T → 0, and not too high a pressure the 2SRP model approaches an effective hard sphere system with diameter d a . We can, then, expect the fcc structure to be the stable phase in the range 51 : 3 , with the transition from the low density fluid to the solid at a pressure 51 of βpσ 3 ≃ 11.6 (σ/d a ) 3 . Assuming that the fcc is indeed the structure of the low density solid phase, and taking into account the low temperature limit for the 2SRP model one can use standard methods to study the equilibrium of this solid with the fluid phase(s).
In order to relate the two ramp models considered in this work, and to explain the procedures used to compute the phase diagram of the systems using thermodynamic and Gibbs-Duhem integration, it is useful to write down a generalised interaction potential, u 1 (r) using an additional 'perturbation' parameter, λ: When λ = 0, there are two interesting limits; when T → ∞ the model approaches that of a system of hard spheres of diameter σ, whereas in the limit T → 0 one once again encounters a hard sphere system, but now having a diameter of d a . In order to compute the Helmholtz energy for condensed phases of the ramp potential systems, one can perform thermodynamic integration over the variable β and the parameter λ. Thus one is able to compute the difference between the Helmholtz energy of a particular state, and the aforementioned hard sphere systems.
In the canonical ensemble the partition function can be written as: where Λ, the de Broglie thermal wavelength, depends on the temperature; R * represents the reduced coordinates of the N particles: ( dR * = 1), and U 1 is the potential energy (given by the sum of pair interactions u 1 (r|λ)). The Helmholtz energy function is given by: One can then write the Helmholtz energy per particle, f = F/N, as: where the ideal contribution f id can be written as and the excess part as In what follows we shall consider the Helmholtz energy as a function of the variables N, V, T , and λ (or, equivalently, as N, ρ, β and λ). We can then compute the derivatives of βf ex with respect to β and λ. I.e. and where the averaged quantity, U a , is given by: On the other hand, using the well known thermodynamic relation: one obtains For the 2SRP model (λ = 0), one can compute the Helmholtz energy of the low density solid using thermodynamic integration. As mentioned earlier, in the limit T → 0 one has an effective hard sphere system of diameter d a . The phase diagram of the hard sphere system 51 , the equation of state for the fluid 55,56 and the equation of state for the solid 57 are all well known. Thus where βf ex hs (ρd 3 ) is the excess Helmholtz energy per particle in the fcc-solid phase and U/N is the excess internal energy per particle. The integrand in (23) We can then compute the Helmholtz energy of the low density solid, as a function of the density on each isotherm, using where ρ 0 = ρ (s) 0 , and βp(T, ρ) = βp s (T, ρ). The Helmholtz energy of the low-density fluid phase is computed using the results of several simulations, having N = 500, at several temperatures and densities; typically ρσ 3 = i×0.025, with i = 1, 2 . . . , 10. In order to guarantee that the samples were indeed within the fluid phase the simulation runs were initiated from equilibrated configurations of temperature T * ≃ 2. The pressure was fitted to a virial equation of state: which is then used in the calculation of the excess Helmholtz energy of the gas: Using (15,16,(26)(27) one can compute the Helmholtz energy of the gas phase as a function of the density. The Helmholtz energies of the high density fluid can be computed using thermodynamic integration from the high temperature limit. Performing a number of simulations, with N = 500, for different values of β at a fixed reference density, ρ l , we have For selected temperatures several simulation runs are performed for various densities, and the results then are used to fit the equation of state of the liquid branch, Similarly, we can compute the Helmholtz energy of the high density fluid at a given temperature by again performing simulations at several densities, and then fitting the results of the pressure as a function of the density and using (25)

Thermodynamic Integration: A2SRP Model
The Helmholtz energy of the fluid phases can be computed following the same procedures outlined in the previous subsection. In order to compute the Helmholtz energy of the low density solid phase one can perform thermodynamic integration starting from the repulsive model and integrating (19) at constant β, ρ, and N: where the indices a and r indicate the A2SRP and the 2SRP models respectively. Such an integration was carried out using λ = k × 0.10 (with k = 0, 1, 2, · · · , 10) with N = 500 at ρσ 3 = 0.24 at the temperatures T * = 0.30 and T * = 0.35. The integrand U a is well behaved in both cases. The Helmholtz energy of the A2SRP model low density solid phase at the reference density ρ

Gibbs-Duhem Integration: 2SRP Model
The partition function, Q(N, p, T, λ) in the isothermal-isobaric (NpT ) ensemble can be written as 51 : The chemical potential, µ is given by Considering constant N, and µ = µ(p, T, λ) one has: where E is the internal energy (including kinetic and potential contributions). Let us consider two phases, α and β, in conditions of thermodynamic equilibrium (i.e. equal values of T , p and µ); if one makes a differential change in some of the variables β, (βp) and λ, thermodynamic equilibrium will be preserved if: where ∆X expresses the difference between the values of the property X in both phases, andX ≡ X/N. In (34) we have considered equal values of the kinetic energy per particle in both phases (classical statistics).
In order to compute the fluid-low density solid equilibrium of the 2SRP model we make The calculation ofŪ andV for both phases at the estimated coexistence conditions was For higher pressures the fluid density becomes higher than the low density solid density, and the temperature decreases with p/T . The third stage is initiated upon reaching T * = 0.40; which happens at (p * /T * ) ≃ 6.33. This involves reverting to the integration scheme given by (35), this time using δT * = −0.01.

Gibbs-Duhem Integration: A2SRP Model
In order to compute the fluid-low density solid equilibria of the A2SRP system, this i.e.: using δλ = 0.025. The final point was found to be (λ = 1, p * /T * = −0.324), corresponding to a equilibrium between a low density solid with density ρ * = 0.248 and a low density liquid with ρ * = 0.239. Given the negative sign of the pressure the equilibrium found is metastable.
From this point we can retake the Gibbs-Duhem integration (now for the attractive model, i.e. setting λ = 1); using (36), with an integration step of δ(p * /T * ) = 0.01. As with the repulsive model a temperature maximum was found at T * M ≃ 0.345, ρσ 3 ≃ 0.255. However, in this case such a maximum seems to be metastable: (p * /T * ) ≃ −0.01. The solid-fluid equilibrium becomes thermodynamically stable at p * /T * ≈ 10 −4 (where this coexistence line crosses the liquid-vapour equilibrium line). At higher pressures and lower temperatures this low density liquid-low density solid coexistence line will eventually coincide with the LL equilibrium line at a triple point. Beyond this a stable low density solid-high density liquid equilibrium will appear.
In order to obtain a precise estimate of the conditions for which the two liquid phases and the low density solid are in equilibria (i.e. the location of the triple point), Gibbs-Duhem integration was performed for the LL equilibrium. Both thermodynamic integration and the Wang-Landau procedure were used to calculate the LL equilibrium, at T * = 0.35. Both of the results agreed to within the error bars: (βpσ 3 ≃ 0.465, with reduced densities of the liquid phase: ρ 1 σ 3 ≃ 0.30, and ρ 2 σ 3 ≃ 0.46). From this initial point the Gibbs-Duhem integrations of the phase equilibrium were carried out using (35) with δT * = 0.0025. The LL equilibrium line met the low density liquid-low density solid equilibrium at T * ≃ 0.331, (p * /T * ) ≃ 0.437, with the three phases having reduced densities: ρ (low density solid) σ 3 ≃ 0.263, ρ (low density liq.) σ 3 ≃ 0.293 and ρ (high density liq.) σ 3 ≃ 0.48 . Gibbs-Duhem integration was then used to calculate the low density solid-high density liquid equilibrium, starting from this triple point, using (35) with δT * = 0.01

Vapour-solid equilibrium
The computation of the density of the low density solid at equilibrium with the vapour phase for the A2SRP model is straightforward. In practice, due to the low pressure at which this equilibrium occurs, the NpT simulations of the low density solid are undertaken at βp ≃ 0, again with N = 500. This immediately yields an accurate estimate of the density of the solid in equilibrium with the vapour phase.

Details of the Gibbs-Duhem integration scheme
A simple predictor-corrector scheme was used to build up the coexistence lines. Generically one seeks the solution of: from some initial condition x 0 , y 0 . In the present case the function f (x, y) has to be computed via a pair of computer simulations. This result will be affected by statistical error. Using a set of discrete values x i = x 0 + i · h, the following scheme was used: with: where (39) and (40) are respectively the predictor and corrector steps. This simple algorithm is both accurate and robust.

III. RESULTS
In order to illustrate the two competing LPS that the A2SRP model system exhibits at low and high densities (or low and high pressures respectively) in Figure 2 we have presented the pair distribution functions for two representative states of the A2SRP model, obtained from the HMSA integral equation. The high density state is close to a hard sphere fluid of diameter σ, whereas the low density g(r) corresponds to a fluid of soft particles of diameter d a = 1.72σ. The purely repulsive system leads to similar results, with only minor differences, due to the lack of dispersive forces (the sharp minima at 3σ are absent). It is clear that the LL transition will result from the 'chemical equilibrium' between two essentially different 'fluids' that stem from the same interaction.

A. Liquid-vapour and liquid-liquid equilibria
In order to calculate the phase diagram using the HMSA integral equation, one has to resort to thermodynamic integration. Although direct closed formulae for the chemical potential are available in the work of Belloni 46 , it has been the authors experience that better results are obtained when using thermodynamic integration. For a given state point, at density ρ (on the right hand side of the phase diagram) and a subcritical inverse temperature β, one can obtain the free energy per particle along a mixed path of the form where β 0 is a supercritical inverse temperatureWith this expression one is able to evaluate the chemical potential by using the relation The phase equilibria can be determined by equating pressures and chemical potentials for both the gas and liquid phases, and for both the low density liquid and high density liquid phases.
The phase diagram thus obtained is plotted in Figure 3  This is best illustrated in Figure 4 where the virial and fluctuation theorem compressibilities are plotted. Here it is observed that the local consistency criterion of Eq.(8) leads to a good agreement for densities ρσ 3 > 0.4, but a considerable inconsistency shows up at lower densities, corresponding to the LL transition (see Fig.3). The assumption of a density independent α parameter in the HMSA closure (Eq. 5) lies at the root of this failure. Thus, whereas the virial pressure predicts a LL transition, no diverging correlations appear near the thermodynamic spinodal, and the equation has solutions throughout the two phase region. This is in marked contrast to the results for the LV transition. On the other hand, in Fig.4 the large values of (∂βP/∂ρ) T in the region 0.2 < ρσ 3 < 0.3 show that the fluid is almost incompressible, which is somewhat surprising, especially at these relatively low densities. This indicates that a solid phase (or a glassy state) may well be about to appear.
We have attempted to explore alternatives to this breakdown in the LC approach. After an unsuccessful attempt to incorporate a linear density dependence on the α parameter in the HMSA interpolating function, in conjunction with a two parameter optimisation, we focused on the GC alternative expressed in Eq. (9). This can be easily implemented using the strategy proposed by Caccamo et al. 44 . To make matters simpler, the GC HMSA calculations were restricted to the purely repulsive 2SRP model. At relatively high temperatures (T * ≈ 2) the procedure worked well and lead to results slightly better that those of the LC HMSA. However, as the temperature is lowered, one finds that the optimisation loop diverges even for relatively low densities. The reason behind this divergence is illustrated in Fig.5. It was observed that for two neighbouring states, a slight decrease in temperature upwardly displaces the residual inconsistency curve. As a consequence of this, when looking for consistency, the minima are now found to be greater than zero, thus the numerical iterations will lead to either a divergence or result in an oscillating behaviour. This also explains why the alternative procedure based on the implementation of a density dependent α failed as well. A completely parallel situation is to be found for the attractive potential model.
It is now clear that a more accurate treatment of this type of soft core models requires not only the implementation of GC conditions on the pressure, but also the use of more flexible closure relations.
Finally, a few words regarding the comparison of the integral equation estimates with the simulation results. As far as the LV transition is concerned, the HMSA predictions are satisfactory. The re-entrant behaviour that is clearly seen on the high density side of HMSA LV curve is also found in the MC estimates, and can be more clearly appreciated in the zero pressure isobar depicted in Fig. 6. Here we see that the errors in the density do not exceed as the curves of compressibility and heat capacity maxima are hardly affected, the phase diagram simplifies substantially. Obviously, the LV transition disappears, and one is left with a rather peculiar fluid-solid transition (see Fig. 10). First, one observes that when density is increased along an isotherm starting from the low density fluid, the systems crystallises into an fcc phase, which melts upon further compression. At higher density values one could Legend as in Fig.8.
find the rhombohedric, cubic and hexagonal phases predicted by Jagla 22 . This melting upon compression is similar to the behaviour of water. In the 2SRP model this is a purely energy driven transition: for certain densities the interparticle distance is necessarily r < d a and the repulsive spheres overlap, thus an ordered configuration no longer corresponds to an energy minimum. At sufficiently large densities (ρσ 3 > 0.9), the hard cores will regain their controlling role and an entropy driven crystallisation will take place (into either fcc or hcp structures). In the intermediate region (ρσ 3 > 0.35) the interplay between entropy and energy will give rise to a much more complex scenario with different solid phases in equilibrium with the fluid phase.

IV. CONCLUSIONS
In summary, we have presented a detailed study of both the attractive and the purely repulsive two scale ramp potential models, using both a self consistent integral equation Despite these limitations, the theoretical approach yields quantitatively correct estimates of the LV transition, the density anomaly and predicts the existence of a LL transition (although with a substantial underestimation of the LL critical density).
By means of extensive Monte Carlo simulations, we have unveiled the rather complicated shape of a significant part of the fluid-solid diagram for the A2SRP models, in which three triple points have been detected. In the purely repulsive model, one finds that the crystallisation preempts the liquid-liquid phase separation. This system has, in common with water, a solid phase that melts upon compression. In the intermediate region of the phase diagram, in which the hard core the the repulsive ramp are competing (i.e. when neither entropy nor energy but a subtle combination of both magnitudes leads the system behaviour) one should expect even more complex phase diagrams for both attractive and repulsive models. This will be the subject of future work.