Hamiltonian limited valence model for liquid polyamorphism

Liquid-liquid phase transitions have been found experimentally or by computer simulations in many compounds such as water, hydrogen, sulfur, phosphorus, carbon, silica, and silicon. Limited valence model implemented via event-driven molecular dynamics algorithm provides a simple generic mechanism for the liquid-liquid phase transitions in all these diverse cases. Here, we introduce a variant of the limited valence model with a well defined Hamiltonian, i.e., a unique algorithm by which the potential energy of the system of particles can be computed solely from the coordinates of the particles and is thus equivalent to a complex multi-body potential. We present several examples of the model which can be used to reproduce liquid-liquid phase transition in systems with maximum valence $z=1$ (hydrogen), $z=2$ (sulfur) and $z=4$ (water), where $z$ is the maximum number of bonds an atom is allowed to have. For $z=1$, we find a set of parameters for which the system has a liquid-liquid and an isostructural solid--solid critical points. For $z=4$, we find a set of parameters for which the phase diagram resembles that of water with a wide region of negative thermal expansion coefficient (density anomaly) extending into the metastable region of negative pressures. The limited valence model can be modified to forbid not only too large valences but also too low valences. In the case of sulfur, we forbid the formation of monomers, thus restricting the valence $v$ of an atom to be within an interval $1=v_{\rm min}\leqslant v\leqslant v_{\rm max}\equiv z=2$.

The limited valence model introduced by Cummings and Stell [34] in 1984 has been used extensively in condensed matter physics [35][36][37][38][39][40][41][42][43][44][45][46][47][48].More recently, it has been shown that a variant of limited valence model can reproduce liquid-liquid phase transitions (LLPT) in some cases listed above [49,50].The main idea of the model is that the system consists of identical atoms interacting via spherically symmetric potentials.Depending on the number of atoms within a certain distance   away from a given atom, this atom is assigned a unique type or valence, :  = 0 with no neighbors within this distance,  = 1 with one neighbor,  = 2 with two neighbors and so on.If an atom has more than  neighbors, where  is the limited valence, the potential energy of the system is assigned to be infinity and hence such a configuration cannot exist.After the types of atoms are determined by counting their neighbors, the potential energy is computed using a table of pair-wise potentials    () defined as the potential energy of the pair of atoms of types  and  separated by distance .Similarly, we can forbid not only too large valences but also too low valences assigning atoms with  <  min an infinite potential energy.The model can be straightforwardly implemented by the Monte-Carlo method.
This approach does not contradict the quantum understanding of inter-atomic interactions which suggests that for a fixed position of ions, the solution of the Schrodinger equation for electrons has a ground state with the energy which depends in a very complex way on the position of the ions and, thus, does not allow simple parametrization of the force field in terms of the ion coordinates.Recently these difficulties were addressed by the machine learning [51] in which the force field is learned from an extensive data base of sample solutions for typical ionic configuration.However, machine learning does not provide us with a simple physical interpretation of the inter-atomic interactions.The maximum valence model, on the other hand, provides a simple cartoon interpretation of this mechanism if we postulate that the changes in the density in a very small vicinity of atoms limited to a length of a covalent bond or a hydrogen bond leads to drastic effects on longer range inter-atomic interactions.In the case of hydrogen, the pioneering quantum calculations of Wigner and Huntington [52] suggested that at high density the energy of the metallic lattice is lower than the energy of the molecular lattice.In other words, at high densities, H 2 , dimers disassociate due to a steric effect; i.e., the electron shells of H 2 dimers are getting larger than the intermolecular distances.This steric repulsion is expected to happen not only for the case of hydrogen but also for other molecular liquids with higher valences.Similar arguments may be applied to the case of water.It was found experimentally that, at high densities, when the fifth neighbor is being pushed into the first coordination shell of a molecule, hydrogen bonds bifurcate and the energy of the bifurcated bonds is roughly half of the energy of the straight bonds [53].Essentially, this implies that the shells of the atoms with coordination number 4 become impenetrable for a fifth intruder at low temperatures.
Historically, the limited valence model was implemented using a reaction scheme in the discrete molecular dynamics (DMD) algorithm [54].DMD which is also known as event driven molecular dynamics was first applied in 1959 to a system of hard spheres [55].The simplest case is the model with  = 1 [50], which mimics the dimerization of metallic hydrogen into molecular hydrogen when the pressure is lowered.The model also reproduces a LLPT in sulfur for  = 2 [49], which is associated with the polymerization at high pressure [12], and a LLPT for  = 4 [50], which mimics the hypothesized LLPT in the four coordinated liquids such as water and silica and also reproduces the region of density anomaly known to exist in these liquids.
In the previous publications [49,50] using the DMD algorithm, the Hamiltonian requirement for the limited valence model was omitted.This omission, obviously contradicts the regular definition of statistical mechanics, since in this case the potential energy of the system depends not only on the configurational space but also on the trajectory of the system in this space, thus making the definition of the entropy ambiguous.This problem also arises in a directional variant of the limited valence model based on the Kern-Frenkel model of patchy colloids [41] in which a pair of particles with overlapping patches may or may not have a bond, but if the patches are so wide that a patch overlaps with multiple patches of other particles, only one bond per patch is allowed.This is solved by introducing random reassignment of bonds among pairs of particles for which their geometrical position allows them to have a bond, so that only at most one bond per patch can exist at each moment of time.The reassignment is produced with a probability proportional to the Boltzmann factor of the bond formation, which ensures correct evaluation of entropy via thermodynamic integration.
In this paper, we make a Hamiltonian formulation of the model by simply forbidding all other particles, except the bonded ones, to be closer to each other than the bond length,   .We reexamine the results obtained previously [49,50] for  = 1,  = 2, and  = 4 and see how the Hamiltonian condition alters the results.
The paper is organized as follows.In section 2 we define the Hamiltonian limited valence model and describe its in detail.In section 3 we compare the results of the Hamiltonian variant of the model with the previously studied non-Hamiltonian variant [50] for  = 1 and study how the LLPT changes with the width of the repulsive shoulder   in the dimerization transition ( = 1).In section 4 we study the same question for  = 4 and show that for the wide shoulder   ≈ 1.34, the phase diagram resembles that of water.In section 5 we study the problem of polymerization for  = 2 as function of the bond potential strength and find the values at which the positions of the LLPT and liquid-gas phase transition (LGPT) resemble those in sulfur.

The model
The original DMD limited valence model for spherically symmetric step potentials is described in reference [49].The potential energy of the system changes when the particles collide, forming bonds or coming to the distance from each other at which the interparticle potential has a step.Moreover, certain pairs of particles when they come to a particular distance from each other can react and change their types.Some of the reactions are bond forming and they are reversible, so that when the bonded pair of particles come to the distance equal to the limited bond length, the bond can break and the particles change their types back to their original types before the bond formation.The model can produce a cascade of reactions in which the type of particles is determined by how many bonds they have, so if a particle  of type  with  bonds react with a particle  with  bonds they can form a new bond and their types become  + 1 and  + 1, respectively, if  + 1 ⩽  and  + 1 ⩽ , where  is the maximum valence.The potential energy of the system changes by Δ, where Δ is the potential energy change of the system associated with the changes of the types of the reacting particles with each other and due to their interactions with all other particles.Δ also includes a term associated to the change of the potential energy of the interacting pair due to the bond formation.Note that if before the reaction a pair potential energy of the reacting pair is   , after the reaction the pair potential energy is   =   +   , where   is the bond energy.In the example of figure 1(b),   = −6.Before the reaction, the pair of reacting particles separated by the reaction distance   has only the energy of van der Waals interaction   = −, while after the reaction the pair has a potential energy   +   = −7.
Conversely, the bond between the two particles of types  > 1 and  > 1 breaks if their distance exceeds the bond length.The particles change their types to  − 1 and  − 1, respectively, and the potential energy of the system changes by the bond energy −Δ.The kinetic energy of the colliding pair changes by the opposite value so that the total energy is conserved.Besides this, the colliding particles change their velocities so that the total linear and angular momenta of the pair are conserved.The laws of energy and momenta conservation give six scalar equations from which the new components of velocities of the colliding pairs can be found by solving a quadratic equation.If this quadratic equation [54] does not have real roots, the particles does not change their types, the bond does not break or form and the particles bounce off like hard spheres by changing the signs of their velocity components in the reference frame of the center of mass of the colliding pair.
This scheme possesses almost all features for the implementation of the Hamiltonian formalism except that we must make sure that the particles of bond forming types cannot exist within the bond forming range without forming a bond, so that the number of bonds of each particle exactly corresponds to its type.This can be achieved by adding an effectively infinite repulsive shoulder with the diameter equal to the bond-forming distance for the interaction potential of the pair of particles with the bond forming types.We chose the value of this quasi-infinite shoulder to be 100, where  is the energy of the van der Waals interactions.
Another feature of the DMD algorithm that was used in references [49,50] is that each atom is split into the "shell" which is a bond forming reactive part and the "core" which is inert and is responsible for the van der Waals long-range interactions.The shell and the core were connected by a short infinite square well bond of length .In this work, we eliminate that unnecessary split and assign both the van der Waals interactions and bond formation ability to a single particle representing an entire atom, which significantly increases the computational efficiency of the model, because it reduces the number of particles by a factor of two.We verify that the phase diagram of the core-shell model converges to the phase diagram of the single-particle model when  → 0.
In principle, the limited valence model allows us to construct the tables of interatomic potentials    () and bond potentials     () for bonds between atoms of different types as complex as we wish, with multiple steps at different distances .This might be necessary in order to accurately reproduce the behavior of various compounds.For example, one may want to make interatomic potentials consisting of many small steps to create a ramp like in the Jagla model [56,57] or make a bond potential resembling a parabola of a harmonic oscillator [58].But in this paper our goal is to keep the model as simple as possible with potentials being a square well potential or a combination of a square well and a square shoulder.Moreover, we choose interatomic potential to be the same for atoms of all types  <  and  < , adding to it a repulsive shoulder if  =  or  =  to model dissociation of bonds under pressure as in hydrogen, or adding an attractive well if both  =  and  =  to model association under pressure as in sulfur.This reduces the number of parameters of the model to six: the maximum valence , the width of van der Waals attraction well, , (the diameter of hard core of this well, , and the depth of this well  serving, respectively, as units of length and energy), the width of the bond   , the energy of the bond,   , the width of the repulsive shoulder or an attractive well,   , and the height of the shoulder (depth of the well),   , with   > 0 for repulsion and   < 0 for attraction.Figure 1(a, b, c) illustrates the potentials of the single-particle model for the case when the atoms with maximum valence have an additional potential shoulder which repel all types of atoms [figure 1(c)].This variant is used to model hydrogen and water.Figure 1(d, e, f) illustrates the potentials for the case when the atoms with maximal valence have an additional attractive well which attracts only atoms of maximum valence [figure 1(f)].This variant is used to model sulfur.
To summarize, the model has three essential potentials, without which it cannot reproduce the phenomenon of liquid polyamorphism for all values of  > 0: the van der Waals potential: the bond potential for 0 <  ⩽  and 0 <  ⩽ : and the modified van der Waals potential: (2.3) Note, that matrix    is symmetric.For the case of dissociation under pressure (  > 0 like in hydrogen and water), we take  <  and  <  in equation ( 2.1) and we take  =  or  =  in equation (2.3).In the case of association under pressure (  < 0 like in sulfur) we take  <  and  ⩽  in equation ( 2.1) and we take  =  =  in equation (2.3).In principle, we can take  = 0 if  ⩾ 3 because the liquid for  ⩾ 3 can form due to multipe bonds between the particles wihout van der Waals attraction [38].
All results below are presented in dimensionless units based on the three basic units [49,50]: , the hard core diameter is the unit of length; , the depth of the van der Waals potential well is the unit of energy; and the mass of a particle, , is the unit of mass.We perform all simulations in the NVT ensemble with a constant number of atoms  = 1000, in a cubic box with periodic boundaries of fixed volume  such that the density  = / has a given value.The temperature is controlled by the Berendsen thermostat [59] modified for the DMD algorithm [54].To estimate the location of the critical points, we first perform preliminary simulations of the single particle model for 10 5 time units by slowly cooling the system from the initial temperature greater than critical temperature to a very low temperature which is below all interesting features of the phase diagram, saving average temperature, pressure, and potential energy every 1000 time units, which is well above the equillibration time, estimated to be below 500 time units in the entire range of temperatures and densities studied.Then, for selected points in the parameter space we perform long equillibration runs of 10 5 time units for each state point.

Liquid-liquid phase transition associated with dimerization (𝒛 = 1)
The results for the single-particle model with the Hamiltonian restriction does not differ much from the previously obtained ones for the core-shell model without the Hamiltonian restriction [50].Here, we tested the case of   = 1.1,   = 1.2,   = −6,   = 12, which corresponds to almost totally impenetrable shoulder.We find the LLCP at   = 2.33,   = 9.07,   = 0.575 for the single-particle Hamiltonian model instead of the previously found [50] for the core-shell model   = 2.12,   = 8.17 and the same density   .The shape of the coexistence curve on the  −  plane in reduced coordinates closely follows the tilted parabola published in figure 3 of reference [50] [figure 2(a)].We also carefully tested the slope of the coexistence line on the  −  phase diagram [figure 2(b)] and found a definitely negative slope d/d = −0.7 ± 0.2, but this value is two times smaller than for the maximum valence core-shell model without the Hamiltonian restriction.
To understand what causes this difference, we perform the simulations of the core-shell model with the Hamiltonian restriction and compare it with the results of the same model without the Hamiltonian restriction as in reference [50].We obtained very similar results:   = 0.575,   = 2.11,   = 8.07, and d/d = −1.6.This is not surprising, because the number of shells which stay within the bond One can see that the results are very close to each other, but the Hamiltonian variant produces a slightly narrower coexisting curve but not significant in reproducing the hypothesized transition for real hydrogen [11].(b) The coexisting lines on the  −  plane.The dashed straight line is the straight line fit of the results of reference [50].One can see that in reduced coordinates, the slope of the coexistence line in the Hamiltonian model is two times smaller than in the non-Hamiltonian core-shell model, However, it is still negative as it is predicted for real hydrogen [11].
distance from another shell without forming a bond and without the potential energy loss associated with the bond is small in the range of temperatures and densities we studied.For example, for  = 2.11,  = 0.58, in the bond range of the shell of type 1, there are on average 0.015 non-bonded intruders, while in the bond range of the shell of type 0, there are on average 0.005 intruders.Such small values are caused by a very high repulsive shoulder   = 12.For smaller values of   , we expect stronger differences between the Hamiltonian and the non-Hamiltonian variants of the model.The larger discrepancy with the Hamiltonian single-particle model is caused by the fact that the core and the shell of the same atom do not exactly coincide with each other but are connected by the permanent bond of distance  = 0.1.In principle, the Hamiltonian core-shell model must converge to the Hamiltonian single-particle model when  → 0. We test this convergence in figure 2. The values of   remain 0.575 for the entire range of , while the values of   ,   and d/d change almost linearly with , smoothly converging to the values of the single-particle model ( = 0).
In addition, we study the effect of the width of the repulsive shoulder   on the location of the critical points and the slope of the coexistence line (figure 4).We observe that the critical density, temperature, and pressure of the LLPT all strongly decrease when the width of the repulsive shoulder increases.By contrast, the slope of the coexistence line behaves non-monotonously with a minimum at   = 1.2, where d/d is definitely negative as hypothesized to be the case for real hydrogen [11].
At the widest repulsive shoulder   = 1.4, the dimeric (molecular) liquid becomes metastable with respect to the monomeric (metallic) liquid.This is because the van der Waals effective attractive range of dimers becomes just  −   = 1.5 − 1.4 = 0.1, compared to practically impenetrable shoulder of   = 1.4.Such a system would have a very low temperature of the liquid-gas critical point (LGCP) which is metastable to crystallization [60].However, the dimeric liquid also becomes metastable with respect to the metallic liquid in which monomers lacking the repulsive shoulder may have 12 or more neighbors in the wide van der Waals well of  = 1.5, which gives the average potential energy per particle −6 = −6.If we assume that an atom in the dimeric liquid has 11 neighbors in a very thin attractive well (one place is taken by the atom forming a bond with energy −6), the potential energy per atom is (−11 − 6)/2 = −8.5, but the entropy cost for forming such a dense configuration is such that the liquid near the critical point     never achieves such a density and its free energy is greater than that of monomeric liquid.
At the narrowest repulsive shoulder   = 1.15, we observe spontaneous crystallization into the face centered cubic (fcc) crystal upon cooling for  ⩾ 0.75, indicated by the abrupt drop of pressure, because the ordered fcc crystal occupies less volume than the disordered liquid.The LLPT is still present at  = 0.735,  = 3.88,  = 0.441.In addition, we observe the isostructural phase transition between two fcc crystals with different lattice constants.This phase transition ends at the isostructural critical point [61] at larger  = 0.78, but smaller  = 3.21 (figure 5).It is worth mentioning that the isostructural critical point was discovered experimentally and was explained theoretically using coresoftened potentials [61,62] long before the LLCP [63].The systems with core-softened potentials with a narrow repulsive shoulder are prone to crystallization at high density, while the maximum valence model based on the same potentials, with the restriction that only  neighbors can enter the shoulder, is not [43].
Here, we observe a novel situation, when both the solid-solid and the liquid-liquid critical points can be observed simultaneously.

Liquid-liquid phase transition in tetrahedral networks (𝒛 = 4)
The hypothetical LLPT in the supercooled water is associated with the reordering of the tetrahedral network of hydrogen bonds [30].The justification of the application of the limited valence model with  = 4 is as follows.In the low density liquid, each oxygen is linked with its four neighbors by 4 hydrogen bonds forming a rigid first coordination shell, penetration into which has a high energy cost modelled by a repulsive shoulder   .These oxygens are assumed to be of type 4. In the high density liquid, some of the bonds are bifurcated.These bonds have much higher potential energy than the normal bonds [53].Thus, in our model we assume that the oxygens which have a bifurcated bond have only 3 bonds and these oxygens belong to type 3.They are not surrounded by a repulsive shoulder and hence they allow other oxygens of type 3 to enter their first coordination shell.The core-shell variant of this model has been already studied in reference [50] and shows the LLPT with a negative slope of the coexistence line, surrounded by a region of density anomaly on the  −  plane.The region of density anomaly existed only at very high pressures and did not reach  = 0 line like in real water.
Here we will study how the change of the repulsive shoulder width   affects the position of the LLCP and the entire phase diagram for  = 4.The motivation for this study is to try to relate the known experimental properties of the water radial oxygen-oxygen distribution function  OO () (see, e.g., figure 6 of reference [30]) to the shoulder width,   .One can see that both high density liquid (HDL) and low density liquid (LDL) have a sharp peak of the radial distribution function at  ≈ 0.27 nm.By contrast, the LDL have a minimum at  = 0.33 nm absent in the HDL which instead has a shoulder at this point and meets the LDL graph at  = 3.9 nm.In terms of the limited valence model, this means that the repulsive shoulder   of the atoms with  = 4 bonds (which prevents the fifth intruder to enter the vicinity of the first coordination shell) extends up to  ≈ 3.9 nm.Thus, we can conclude that in order to achieve good agreement with water we need to set 3.3/2.7 ≈ 1.2 <   < 3.9/2.7 ≈ 1.4.Figure 6 shows the dependence of the location of the critical points on   for  = 4,  = 1.5,   = 1.1,   = −6,   = 6.The behavior of all the variables qualitatively resembles their behavior for  = 1.Except that the dependence is weaker and the critical pressure becomes negative for the shoulder width   = 1.35.The slope of the  One can see that as   increases, the critical density, temperature and pressure decrease, until the pressure becomes negative and the LLCP disappears below the liquid gas spinodal for   ⩾ 1.4.The liquid-gas critical temperature and pressure also decreases.The liquid-liquid critical temperature is always below the liquid-gas critical temperature as in water.
liquid-liquid coexistence line is always negative.For   = 1.4, the LLCP becomes submerged below the liquid-gas spinodal.For   = 1.15, the system crystallizes into a body centered cubic lattice at high densities.The isostructural solid-solid critical point may exist at even higher densities, but this question requires further investigation.Figure 7 shows the  −  phase diagrams which resemble those of water [30], showing the region of density anomaly at negative pressures [64], the compressibility minimum and the heat capacity maximum near the critical point.
The exact location of the hypothetical (LLCP) of water is unknown.The simulations of the ST2 model [65] give the value of   = 195 MPa and   = 245 K.More recent estimates based on the TIP4P/ice model [66] suggest the values   = 125 MPa,   = 195 K. Theoretically, the LLCP can be at negative pressures or can even disappear below the liquid-gas spinodal [30].The LGCP of water is located at   = 22 MPa which is several times smaller than   , and   = 647 K, which is several times larger than   .It is clear that the maximum valence model with  = 4 can reproduce these values after some additional tweaking of parameters.Beside this, the maximum valence model qualitatively correctly describes the water anomalies.The best agreement with water is reached at   = 1.34 which does not contradict our estimates based on  OO ().

Liquid-liquid phase transition associated with polymerization of dimers (𝒛 = 2)
The maximum valence model for  = 2 has been used to mimic the LLPT associated with polymerization [49] as was experimentally observed in sulfur at high pressures [12].One of many unrealistic features of that model is that real sulfur remains molecular S 2 near its liquid-gas critical temperature, while in the limited valence model the dimers at such a temperature dissociate into monomers.Here, we modify the maximum valence model in its Hamiltonian formulation assigning the atoms with zero neighbors in their bond range an infinite potential energy.Within the DMD package, this can be achieved by modifying a bond potential between an atom with one bond and its neighbor which may have either One can see the isochores crossing at the critical points surrounded by the region of the density anomaly as in reference [65] for the ST2 model of water.This region is separated by the TMD (temperature of maximum and minimum density) line connecting the points of minimum and maximum pressure at constant density.Also are shown the line of isothermal compressibility,   maxima and minima (KTM) which crosses the TMD line at its vertical point [67] and the constant pressure specific heat maxima (CPM) line (green) which coincides with the KTM line near the critical point [67].The TMD, KTM and CPM lines are obtained by the polynomial interpolation of the computed points on the isocores, obtained with steps Δ = 0.01 and Δ = 0.01.Some of the inflections and irregularities on these lines can be due to the boundaries of the fitted region.The low density isochores in (b) start to intersect with each other forming a liquid-gas spinodal, which becomes non-monotonious following the isochores with a minimum at the point of maximum density like in the IAPWS95 [68] equation of state of water extrapolated to the metastable region of negative pressures.
one or two bonds a quasi-infinite very narrow square potential of barrier as shown in figure 8.In this new model, the only possible types of atoms are atoms with one bond forming dimers or being at the end of a polymer chain and atoms with two bonds forming a chain.We tested that near the LGCP, the gaseous phase has only 5.5% atoms with two bonds, which are the middle atoms of the trimers and no atoms with zero bonds.The potential of the bond between atoms of type 1 and 1 and atoms of 1 and 2. preventing atoms of type 1 to separate and become monomers.The height of the potential barrier can be taken as high as 1000 and the width  can be taken as small as 0.0001.
Qualitatively, the new Hamiltonian model with  = 2 and forbidden monomers have phase diagrams similar to the cases studied in reference [49].Here, as an example, we present a  −  phase diagram for  = 1.7,   = 0.01,   = 1.3,   = 2 and   = −1 (figure 9).Note that the atoms belonging to polymer chains attract each other due to the additional attraction for atoms of type 2, but the cost of forming a bond is positive.Moreover, the bonds are very narrow.Technically, it is not the bonds which hold the polymer chain together but the additional attractive potential between the atoms with two bonds.These bonds with positive energy are not necessarily very unstable.They can be metastable with as long life time as we wish if we add to the bond potential a large enough potential barrier as in figure 8.In the limit of  → 0, this potential will not change the Hamiltonian, and thus will not affect the phase diagram, but will dramatically slow down the kinetics.Whether this type of bond potential can be justified by quantum chemistry remains an open question.
Nevertheless, the model qualitatively describes the LLPT in sulfur [figure 9(a)].It shows two phase transitions and their spinodals indicated by the loci of isochore crossing which join at critical points: the liquid-liquid one indicated by a circle at   = 0.865 ± 0.005,   = 1.86 ± 0.01, and   = 5.54 ± 0.05 and a liquid-gas one indicated by a square at   = 0.275 ± 0.025,   = 2.30 ± 0.05,   = 0.078.Note that in real sulfur   = 1314 K,   = 20.7 MPa, and   = 563 kg/m 3 [69], while the LLCP is located at   = 1035 K,   = 2.15 GPa, and   ≈ 2000 kg/m 3 [12], making the ratios of critical parameters for real sulfur equal to   /  = 3.55,   /  = 0.79,   /  = 104, while for the model with the given set of parameters they are, respectively, 3.2, 0.81 and 7.1.As one can see from figure 9(b) a further increase of   , can easily increase both   /  and   /  .Note also that the liquid-liquid coexistence line must belong to a narrow region between the spinodals and thus must have a very large positive slope.Moreover, it crosses the liquid-gas coexistence line at a very small temperature and at a small pressure almost equal to zero, where crystallization can be anticipated, making the phase diagram very realistic, because according to reference [12], the liquid-liquid coexistence line in sulfur ends at a triple point with the crystalline phase.In order to improve the agreement with the experimental data on sulfur, we explore how the location of the critical points is affected by the bond strength   ≡ Δ 0 (enthalpy of bond formation), figure 9. We see that for negative bond formation enthalpy (exothermic reaction) leads to a very low critical pressure of the LLCP making it even negative if   < −2.As   grows,   increases and for   ≈ 3 it becomes 100 times larger than   as in real sulfur.When   increases, the critical density,   , also increases from 0.63 at   = −2 to 0.915 at   = 4.The critical temperature   slowly increases with   , becoming greater than   for   > 3, which is not correct for real sulfur.The way to overcome this problem is to increase   from −1.3 to −1.0.This significantly reduces   without any significant effect on other critical parameters.The parameters of the LGCP do not depend on   and   in any significant way.

Conclusion
We develop a Hamiltonian variant of the limited valence model in which each atom is represented by a single particle and implement it via the event driven molecular dynamics.The new model has a well-defined Hamiltonian, which allows it to be treated within the formalism of statistical mechanics.For example, it has a well-defined entropy and can be studied by the Monte Carlo method.We show that this model provides a simple way to reproduce LLPT in various substances.The six independent parameters of the model , ,   ,   ,   , and   strongly change the location of the LLPT and its interplay with the LGPT.We also study the relation of this new model with the previously studied core-shell non-Hamiltonian limited valence model of references [49,50] and show that qualitatively all the phenomena observed there are reproduced by the new model, which is much faster computationally.
For  = 1, we show that the model can reproduce molecular-metallic transition in liquid hydrogen.For  = 2, it can reproduce the LLPT associated with polymerization of dimers in sulfur.And for  = 4, it can reproduce the phase diagram of water with a LLPT at low temperature and high pressure.The limited valence model is thus a powerful, versatile and simple tool to study LLPTs.We hope that in future it can be applied to the LLPT in phosphorus for  = 3 and in many other compounds.
We find that for  = 1 and narrow repulsive shoulder   = 1.15, it is possible to simultaneously observe the LLCP and the solid-solid isostructural critical point at slightly different densities and temperatures.
For the case of  = 4, we find that for the large values of   ≈ 1.34 which can be relatated to the position of the minimum of the oxygen-oxygen radial distribution function in water, the phase diagram of the model resembles the phase diagram of water obtained by all-atom simulations.
For the case of sulfur ( = 2) we modify the limited valence model to also include the minimum valence  min = 1, forbidding the formation of monomers.

Figure 1 .
Figure 1.(Colour online) Panels (a, b, c).Pairwise potentials    which we use to model hydrogen ( = 1) and water ( = 4).Atoms interact with the square well potential of width  >   >  and energy − [black line, panel (a)] if their types are  <  and  < .If they come closer than the bond width   , they change their types to  + 1 an  + 1 and they form a bond [magenta line, panel (b)] and their pair potential becomes   =   +   = −7 (in this example,   = − and   = −6).If the type of at least one atom in the interacting pair is , their interaction potential has an additional square shoulder of width   and height   [blue line, panel (c)].Both blue and black lines go to infinity at the bond distance   , to satisfy the Hamiltonian condition.Panels (d, e, f) are the same as (a, b, c) but for the case of sulfur ( = 2).Atoms interact with the square well potential of width  >   >  and energy − [black line, panel (d)] if their types are not both equal to :  <  and  ⩽ .If they come closer than the bond width   , they change their types to  + 1 an  + 1 and form a bond [magenta line, panel (e)] and their pair potential becomes   =   +   (in this example,   = − and   = 3).If the types of both atoms in the interacting pair are , their interaction potential has an additional square shoulder of width   [blue line, panel (f)].Both blue and black lines go to infinity at the bond distance   , to satisfy the Hamiltonian condition.

Figure 2 .
Figure 2. (Colour online) Comparison of the liquid-liquid coexistence lines of the single-particle Hamiltonian maximum valence model and the core-shell non-Hamiltonian model studied in reference [50] for  = 1,  = 1.5,   = 1.1,   = 1.2,   = −6,   = 12.(a) Coexistence line in the  −  plane in reduced coordinates.Circles are the results obtained by Maxwell construction for the Hamiltonian model.The jumps in the results are due to the discreteness of the values of  which are produced with step Δ = 0.005.The black curve shows the same results after a smoothing procedure.The red curve shows the results for the core-shell model without the Hamiltonian restriction.One can see that the results are very close to each other, but the Hamiltonian variant produces a slightly narrower coexisting curve but not significant in reproducing the hypothesized transition for real hydrogen[11].(b) The coexisting lines on the  −  plane.The dashed straight line is the straight line fit of the results of reference[50].One can see that in reduced coordinates, the slope of the coexistence line in the Hamiltonian model is two times smaller than in the non-Hamiltonian core-shell model, However, it is still negative as it is predicted for real hydrogen[11].

Figure 3 .
Figure 3. Convergence of the phase diagram of the core-shell model to that of the single-particle model as the bond length connecting the core and the shell tends to zero.The parameters of the model are taken for the case attempting to mimic dimerization of metallic hydrogen:  = 1,  = 1.5,   = 1.1,   = 1.2,   = −6,   = 12.Panels (a, b, c) show the parameters of the critical points   ,   and / as functions of .One can see that their dependence on  is almost linear, including the last point for the single-particle model  = 0.

Figure 4 .
Figure 4. (Colour online) (a) Dependence of the location of the critical points on the repulsive shoulder width,   for the dimerization model  = 1,  = 1.5,   = 1.1,   = −6 and   = 13.The critical density decreases from 0.735 at   = 1.15 to 0.285 at   = 1.4.The critical temperature decreases by factor of 3 and pressure decreases even more dramatically, approximately as 0.056(  −   ) −2 .(b) The slope of the coexistence line versus   estimated by the slopes of the two isochores with step  = 0.01 crossing in the  −  plane at the maximum temperature, which is also how we obtain the location of the critical points in panel (a).

Figure 5 .
Figure 5. (Colour online) (a) The  −  phase diagram of the dimerization model ( = 1) with   = 1.15,  = 1.5,   = 1.1,   = −6 and   = 13.The black curves are parts of the isochores corresponding to the liquid state, and the red curves are parts of the isochores corresponding to crystalline state, obtained by slow cooling with equilibrarion time of 1000 time units.The crystallization event is shown at the isochore with  = 0.77.The two critical points are visible: the liquid-liquid one indicated by a circle and the isostructural solid-solid one indicated by a square.(b) The snapshot of the system taken along the the liquid-liquid critical isochore  = 0.74 above the critical point at  = 4.44,  = 49.7,(c) the snapshot of the system taken along the same isochore  = 0.74 below the critical point at  = 3.02,  = 39.5,(d) the snapshot of the system taken along the solid-solid critical isochore  = 0.77 above the critical point at  = 4.8,  = 47.7,(e) the snapshot of the system taken along the the solid-solid critical isochore  = 0.77 below the critical point at  = 2.42,  = 31.4.In (d) and (e) one can clearly see the crystalline plane (1,0,0) of the fcc lattice.Atoms of type 0 (monomers) are shown in red.Atoms of type 1 (belonging to dimers) are shown in white.

Figure 6 .
Figure 6.(Colour online) The dependence of the location of the critical points for the tetrahedral model  = 4,  = 1.5,   = 1.1,   = −6,   = 6 on the width of the repulsive shoulder   .One can see that as   increases, the critical density, temperature and pressure decrease, until the pressure becomes negative and the LLCP disappears below the liquid gas spinodal for   ⩾ 1.4.The liquid-gas critical temperature and pressure also decreases.The liquid-liquid critical temperature is always below the liquid-gas critical temperature as in water.

Figure 7 .
Figure 7. (Colour online)  −  phase diagrams of the  = 4 with  = 1.5,   = 1.1,   = −6 and   = 6 for   = 1.34(a) and   = 1.35 (b).One can see the isochores crossing at the critical points surrounded by the region of the density anomaly as in reference[65] for the ST2 model of water.This region is separated by the TMD (temperature of maximum and minimum density) line connecting the points of minimum and maximum pressure at constant density.Also are shown the line of isothermal compressibility,   maxima and minima (KTM) which crosses the TMD line at its vertical point[67] and the constant pressure specific heat maxima (CPM) line (green) which coincides with the KTM line near the critical point[67].The TMD, KTM and CPM lines are obtained by the polynomial interpolation of the computed points on the isocores, obtained with steps Δ = 0.01 and Δ = 0.01.Some of the inflections and irregularities on these lines can be due to the boundaries of the fitted region.The low density isochores in (b) start to intersect with each other forming a liquid-gas spinodal, which becomes non-monotonious following the isochores with a minimum at the point of maximum density like in the IAPWS95[68] equation of state of water extrapolated to the metastable region of negative pressures.
Figure 8.The potential of the bond between atoms of type 1 and 1 and atoms of 1 and 2. preventing atoms of type 1 to separate and become monomers.The height of the potential barrier can be taken as high as 1000 and the width  can be taken as small as 0.0001.

Figure 9 .
Figure 9. (Colour online) (a) The  −  phase diagram for model of the polymerization of dimers ( = 2,  = 1.7,   = 0.01,   = 1.3,   = 2 and   = −1).The lines are isochores with step Δ = 0.05 from  = 0.05 to  = 0.85 and with step Δ = 0.01 from  = 0.85 to  = 0.95.The LLCP is indicated by a circle while the LGCP is indicated by a square.The results are obtained by slow cooling with equillibration time for each point 1000 time units.(b) The dependence of the critical point location on the repulsion/attraction energy of the bond   for  = 2,  = 1.7,   = 0.01,   = 1.3, and   = −1.3.The values for LLCP (filled symbols) are connected with bold lines, while the values forLGCP (empty symbols) are connected with thin lines.One can see that changing bond from attractive to repulsive significantly increases   , making the ratio   /  ≈ 100 as in real sulfur for   ≈ 3. Two additional points are given for   = −1, which yields a better agreement with the experimental data because reducing   significantly reduces   making it below the   as in real sulfur, without affecting other critical values.These two points for   = 2 and 3 are connected with dashed lines.