Sensitivity of solid phase stability to the interparticle potential range: studies of a new Lennard-Jones like model

In a recent article, Wang et al (Phys. Chem. Chem. Phys., 22, 10624 (2020)) introduced a new class of interparticle potential for molecular simulations. The potential is defined by a single range parameter, eliminating the need to decide how to truncate truly long-range interactions like the Lennard-Jones (LJ) potential. The authors explored the phase diagram for a particular value of the range parameter for which their potential is similar in shape to the LJ 12-6 potential. We have reevaluated the solid phase behaviour of this model using both Lattice Switch Monte Carlo and thermodynamic integration. In addition to finding that the boundary between hexagonal close packed (hcp) and face centred cubic (fcc) phases presented by Wang et al was calculated incorrectly, we show that owing to its finite range, the new potential exhibits several 'artifact' reentrant transitions between hcp and fcc phases. The artifact phases, which do not occur in the full (untruncated) LJ system, are also found for typically adopted forms of the truncated and shifted LJ potential. However, whilst in the latter case one can systematically investigate and correct for the effects of the finite range on the calculated phase behaviour, this is not possible for the new potential because the choice of range parameter affects the entire potential shape. Our results highlight that potentials with finite range may fail to represent the crystalline phase behaviour of systems with long-range dispersion interactions, even qualitatively.


Introduction
In a recent paper, Wang et al 1  ( Here σ sets the length scale, and is a coefficient that ensures that the depth of the attractive well is −ε.The potential is formulated such as to vanish smoothly at the specified value of the cutoff parameter r c , thus circumventing the question of which truncation scheme to employ when seeking to render computationally tractable a truly long-ranged interaction such as the well known Lennard-Jones (LJ) potential.However, in contrast to a truncated LJ potential, the choice of r c in (1) sets not just the truncation distance, but also determines the overall shape of the potential.The authors find that for r c = 2σ (for which α = 1), its form is similar to that of the LJ potential as shown in Fig. 1.We shall henceforth refer to Eq. ( 1) with r c = 2σ as the Lennard-Jones like (LJL) potential.Eq. ( 1) with rc = 7σ Fig. 1 Comparison of the standard 12 − 6 Lennard-Jones (LJ) interaction potential with the potential of Wang et al 1 (Eq. 1) with range parameters r c listed in the legend.The case r c = 2.0σ is the Lennard-Jones like (LJL) potential.
Wang et al presented various thermophysical properties of the LJL potential including the phase diagram for the vapor, liquid and crystalline solid phases, all computed from free energy measurements via thermodynamic integration (TI).In the present work we have revisited the crystalline solid region of the LJL model using both TI and Lattice Switch Monte Carlo (LSMC).The latter is a powerful method for determining the free energy between coexisting crystalline phases.Our results reveal an error in the measurement of the phase boundary between hcp and fcc phase presented by Wang et al 1 .Moreover we find that on increasing number density the correct phase diagram for the solid region manifests several reentrant fcc and hcp phases -features that are absent in the full LJ potential.We trace the source of these re-entrant phases to the finite range of the LJL potential ie. the lack of long ranged dispersion interactions.Our results for the phase behaviour are set out in Sec. 2 while Sec. 3 provides a discussion and summarises our conclusions concerning the model.

Phase diagram of the Lennard-Jones like model
We have studied the crystalline solid phase behaviour of the LJL model using both Lattice Switch Monte Carlo and thermodynamic integration in conjunction with Molecular Dynamics simulation.We have also performed ground state energy calculations as a function of number density for both the LJL model and a truncated and shifted LJ potential.Implementation details for these methods and the results of our calculations are reported below and in Supplementary Material (SM) 2 .In all cases, numerical values are expressed in dimensionless units, namely reduced temperature T * = k B T /ε, pressure p * = pσ 3 /ε and particle number density ρ * = ρσ 3 .For notational simplicity we shall henceforth suppress the superscript * on these quantities.

Lattice Switch Monte Carlo calculations
LSMC [3][4][5][6] is a well established and powerful simulation technique for determining solid-solid coexistence parameters.Its high precision stems from the fact that it focuses on the difference in the free energy between two candidate stable phases rather than the absolute free energy of each 7 .In the course of a LSMC simulation the system switches repeatedly between two candidate stable structures allowing the accumulation of statistics on their relative statistical weight.To enable such switching, the sampling must be biased to visit -on a regular basis -certain 'gateway' configurations from which a switch from one phase to the other can be launched.The requisite bias function can be obtained by using the transition matrix method 8,9 .The effects of the bias are unfolded from the statistics a posteriori.
We have deployed the LSMC method to obtain the crystalline solid phase behaviour of the LJL model using the implementation included within the general purpose MC simulation engine DL_MONTE 10 .Details regarding the simulations can be found in the SM 2 .We investigate in particular the relative stability of hcp and fcc structures as a function of the particle number density ρ = N/V and temperature T .To do so, we work within the isobaric-isothermal (constant-NpT) ensemble for which LSMC provides direct access to the Gibbs free energy difference between the hcp and fcc phases ∆G = G hcp − G f cc .This is given by ∆G = −k B T log(P hcp /P f cc ) where P hcp and P f cc are the integrated probabilities for the system to be found in microstates typical of hcp and fcc respectively.
Coexistence state points are defined by ∆G = 0. We located the value of the pressure (p) that corresponds to phase coexistence for a prescribed temperature using a root finding algorithm.Specifically, we applied Newton-Raphson's method to the function ∆G(p) which leads, since dG = V dP, to iterate the pressure according to 6 where ∆V i is the volume difference between the two phases at the i th iteration.At each iteration, we computed the LSMC bias function afresh, though it is possible to avoid doing so by deploying histogram reweighting techniques 9 .
Our LSMC results for the phase diagram are presented in Fig. 2. While we find agreement with the work of Wang et al 1 for the vapor-fluid binodal (which we have determined using separate Gibbs Ensemble Monte Carlo (GEMC) simulations within DL_MONTE), our results for the crystalline region are very different.Specifically, we find at least three separate lines of hcp-fcc transitions, none of which coincide with that of Wang et al.In order to check our findings and throw light on the discrepancy, we have performed TI calculations for the solid region.The results of these calculations are described in the next subsection.1) with r c = 2.0σ .Blue dots are phase boundaries reported in Wang et al 1 , symbols are the present work.Green crosses denote liquid-vapour coexistence densities calculated from GEMC.Our LSMC estimates of the reentrant hcp-fcc coexistence boundaries are denoted by + and were obtained in the reduced temperature range 0 < T ≤ 1 and for the particle numbers N shown in the legend.Note the complete disparity with the single hcp-fcc boundary calculated by Wang et al, as marked by blue dots in the solid region.In all cases, the difference in the coexisting hcp and fcc densities is invisibly small on the scale of the graph.Uncertainties for the coexistence densities were calculated as described in the SM 2 and unless indicated by an error bar, do not exceed the symbol size.

Thermodynamic integration calculations
To elucidate the origin of the discrepancy between our prediction for the hcp-fcc coexistence line(s) based on LSMC calculations and the thermodynamic integration calculations of Ref. 1 , we have computed absolute Helmholtz free energies of both phases using another simulation program, GROMACS 2,11,12 , and the TI method 13 for a system comprising N = 768 particles.We selected a number of (T, ρ) state points where absolute Helmoltz free energies have previously been calculated by Wang et al. and given in the supplementary material of Ref. 1 .
For our TI calculations we adopted the integration path described by Aragones et al 14,15 to connect the crystalline phase of interest and a known reference state, namely the ideal Einstein crystal with one atom fixed at its lattice site (state 0) * .The Helmholtz free energy of state 0 is ) and k E is the spring constant of the (N − 1) harmonic oscillators.The absolute free energy A can be expressed as 14,16 where ∆A 1 is the free energy difference between state 1, which is an Einstein crystal where the spring-bound atoms interact through LJ-like interactions, and state 0. ∆A 2 is the free energy difference between state 2, the crystalline phase of interest but with one atom fixed, and state 1.The last term in Eq. ( 4) compensates for the fixed atom constraint.∆A 1 is computed using Bennett's formula where the thermal average is over equilibrium trajectories of the ideal Einstein crystal (state 0).Notice that E 1 − E 0 = E pair is merely the sum of the LJ-like pair interactions.To avoid overflows in the exponent, the formula is rewritten as is the sum of the pair interactions in the perfect crystal at T = 0.The spring constants were adjusted so that the second term (per particle) − ln⟨exp(−β (E pair − E (T =0) pair ))⟩ 0 /N is close to 0.02, as recommended 14,15 .The integral in the contribution is the reference lattice site of the i th particle and the thermal average is over an interacting Einstein crystal with spring constant k ′ E , was computed by using the change of variable 13,16 x = ln(β 1 2 k ′ E (10 −1 nm) 2 + e 3.5 ) and a Gauss-Legendre quadrature with 15 points.
The results of our TI calculations are compared in Table 1 with LSMC calculations that we have performed in the canonical (constant-NV T ) ensemble.The state points listed are a selection of those also studied by Wang et al 1 (see their SM).At T = 0.1, our LSMC and TI both predict a (first) transition from hcp to fcc around density ρ ≈ 1.18(1) (where ∆A(ρ) changes sign), far removed from the density ρ ≈ 1.29 at which coexistence is predicted in Ref. 1 .At higher temperatures we find consistency between our canonical LSMC and TI estimates of transition densities on this first hcp-fcc phase boundary, and the constant-NpT LSMC calculations of the full phase boundary reported in Fig. 2. While error bars are difficult to assess in the TI method, they can be calculated straightforwardly 2 for LSMC and standard deviations on the LSMC results are reported in the table.Absolute free energies depend on the atomic mass via the thermal de Broglie wavelength Λ T , while the free energy difference ∆A does not.Our TI values * Note that this is a different choice of integration path to that adopted by Wang et al 1 were computed by setting Λ T = σ , separately at each considered temperature.This standard convention was apparently also used in Ref. 1 .At all densities and temperatures reported in Table 1, the free energy difference ∆A = A hcp − A fcc that we obtain with TI deviates slightly and systematically from our LSMC results by a shift of the order −0.02NkB T .This small temperature-independent shift might be due to systematic integration errors in the TI calculations.Much larger discrepancies are observed between our LSMC and TI values for ∆A and those of Ref. 1 (compare columns 3 and 4 with column 5 of Table 1 respectively).Comparing separately the free energies estimates of the fcc and hcp phases, we find generally reasonable agreement between our TI results and those of Wang et al 1 for the fcc phase, but a much larger discrepancy for the hcp phase, as reported in the rightmost two columns of Table 1.This shows that the absolute free energies of the hcp phase at least have very likely been miscalculated in Ref. 1 , leading to a very different and incorrect crystalline phase boundary.Since the missing contribution is proportional to the temperature, the hcpfcc coexistence curve determined in Ref. 1 (for ρ ≤ 1.4) deviates increasingly from our coexistence curve.We note that the error in the calculation of the free energy also raises questions regarding the correctness of the estimates of the liquid-solid phase boundary reported by Wang et al, as shown in Fig. 2, though we have not checked that particular calculation.

Ground state energy calculations.
The corrected ρ − T phase diagram for the LJL model displayed in Fig. 2 shows that the boundaries between hcp and fcc phases are almost linear and independent of temperature over a wide temperature range which extends down to very low T .It is thus reasonable to assume that coexistence properties calculated from ground state (ie.T = 0) energy calculations may aid an understanding of the overall solid state phase behaviour.
In the ground state, particles occupy their lattice sites, there are no entropic contributions to the hcp-fcc free energy difference, and energy calculations suffice to determine the phase behaviour.We have used the potential eq. ( 1) to calculate the ground state energy difference per particle ∆E/N = (E hcp − E f cc )/N as a function of ρ.The criterion ρ(∆E = 0) corresponds to hcp-fcc transition points.Note though that while the stability of each phase is correctly predicted by the sign of ∆E, the hcp and fcc phases have slightly different specific volumes.The coexisting densities are then not given merely by the condition ∆E = 0, but by ∆G = ∆E + p∆V = 0.They can be deduced from the curve E(ρ, T = 0) of each phase by performing the common tangent construction, which imposes a common pressure at coexistence.However as the coexistence densities straddle the density ρ(∆E = 0) and are very close in value to one another -indistinguishably so on the scale of fig.2-it is sufficient in the current context to simply determine ρ(∆E = 0).
Our ground state energy calculations for the LJL model are shown in figure 3, which plots ∆E/N as a function of ρ.One sees that ∆E = 0 occurs for ρ = 1.195(5) (p = 12.48), which is close to the extrapolation to T = 0 of the first hcp-fcc transition line shown Table 1 Free energy difference ∆A, in units of Nk B T , between the hcp and fcc structures of the LJL solid (r c = 2σ ) at a selection of state points.The discrepancy between the LSMC and TI calculations is small (see 6 th column) and might be due to systematic integration errors in the TI method.The difference between the absolute free energy calculated via TI in this work and that in Ref. 1 are shown for both the fcc and the hcp structures in the rightmost two columns.Note the large discrepancy for the hcp phase.
Fig. 3 Dependence on number density ρ of the ground state energy difference per particle ∆E/N = (E hcp − E f cc )/N between hcp and fcc crystalline phase for the LJL potential.Data is shown for a range of particle number N. The results illustrate the reentrant hcp phase at ρ = 1.815 (5), while the accord between the curves for different particle number N demonstrates that the behaviour is not a finite-size effect.
It is instructive to compare the ground state phase behaviour of the LJL model with that of the truncated and shifted LJ potential which has the form where φ LJ (r) = 4ε σ r 12 − σ r 6 is the full LJ potential, a portion of which is shown in fig. 1.We have used eq.( 5) to calculate the ground state energy as a function of the density for various values of the LJ truncation distance r c .The results, which replicate and extend to higher values of r c similar calculations by Jackson et al 5 , are shown in figure 4.This figure also includes the results of calculations for the limiting case of the full LJ potential given by Stillinger 17 for which it is important to note there is only a single hcp-fcc transition, which occurs at ρ = 2.1728.On increasing r c , our results for ∆E/N approach the limiting curve, but only for a large truncation distance of r c ≳ 6σ .For smaller cutoffs r c ≲ 5σ , ∆E fluctuates wildly and changes sign at several densities, indicating that in this regime additional hcp-fcc transitions arises which are artifacts in the sense that they are not characteristic of the full LJ potential.Also shown for comparison is the limiting case for the full (untruncated) LJ potential 17 .Note how the results for the finite cutoff approach the limiting case for r c ≳ 7.0σ .
We have also considered the changes to the T = 0 phase behaviour of the potential of Wang et al that result from varying the range parameter r c in Eq. 1.The results displayed in fig. 5 show that the number and location of phases occurring varies dramatically, similarly to what was found in fig. 4 for the truncated and shifted LJ potential.In order for a single hcp-fcc transition to occur, as in the full LJ potential, one must utilise in Eq. 1 a range parameter r c ≳ 7σ .However, as shown in Fig. 1 and discussed below, for this value of the range parameter, the potential of Eq. 1 is not at all Lennard-Jones like.Fig. 5 Dependence on number density ρ of the ground state energy difference per particle between hcp and fcc crystalline solid phases as calculated for the LJL potential of Eq. 1 at a selection of values of the range parameter r c .For solid like densities (ρ > 0.8), the smallest value of the range parameter for which only a single hcp-fcc transition occurs is r c = 6.89σ .

Discussion and conclusions
In this work we have reconsidered the phase diagram of a recently proposed finite-ranged potential, Eq. ( 1), that has been recommended as an alternative to the use of truncated LJ potentials in molecular simulation 1 .The advantage of this LJL potential is that while it is similar in shape to the LJ potential, it is defined in such a way as to vanish smoothly at the prescribed value of the range parameter.It thus provides a unique truncation scheme, which contrasts with the multiple approaches in common use for the LJ potential.
We have studied the phase diagram of the LJL model both in the liquid-vapor and the crystalline solid regions.While we confirm the original authors' findings for the vapor-liquid binodal, our results in the crystalline phases deviate greatly.Specifically on increasing the number density starting from the solid branch of the vapor-solid coexistence region, we find a very different hcpfcc boundary to that reported.Our calculated hcp-fcc boundary (fig.2) exhibits relatively weak temperature dependence and intersects the freezing line at a fluid-hcp-fcc triple point.We have traced the discrepancy between our results and those of Wang et al to an apparently incorrect calculation of the free energy on their part.
We find, moreover, that on increasing the density further within the crystalline region, the LJL potential exhibits not one but multiple reentrant hcp-fcc transitions that extend over a wide range of temperatures down to T = 0. Indeed a previous study by Choi et al 18 of the LJ potential, truncated and shifted at r c = 2.5σ , reported a qualitatively similar phase diagram, displaying reentrant hcp and fcc transitions whose boundaries depend only weakly on temperature, as well as a fluid-hcp-fcc triple point.Importantly, however, such features are absent in the full LJ potential which except for a very small temperature range exhibits only a single hcp-fcc transition line, which is strongly temperature dependent at low T and does not intersect the melting curve 5,6,19 .Thus the solid phase diagrams of both the LJL potential and the truncated and shifted LJ system exhibit phases that are artifacts in the sense that they do not occur in the full LJ potential.
It is well established that the difference in free energy between hcp and fcc phases for spherically symmetric interparticle potentials can be small across a range of state points 4,5,13 .Previous studies of the ground state (T = 0) energy as a function of number density of the truncated and shifted LJ potential have shown that artifact transitions arise due to the sensitivity of the sign of the hcp-fcc energy difference ∆E/N to the truncation range 5,6 .We note also that the magnitude of ∆E/N for a truncated potential can be much larger than for the full LJ potential.For instance, when truncating the LJ interaction at some short cutoff, say r c = 2σ , ∆E/N is about 10 times that of the full LJ potential.The size of this energy difference would appear to be responsible for the insensitivity to temperature of the measured hcp-fcc coexistence boundaries, as seen for both the truncated LJ potential 18 and the LJL system (fig.2), because it dominates over the hcp-fcc entropy difference.It is therefore not surprising that the much smaller value of ∆E associated with the full LJ potential leads to very different crystalline phase diagram 5,19 .Ground state energy calculations for the truncated and shifted LJ potential performed in the present work (fig.4) demonstrate that, at least at T = 0, one requires r c ≳ 6σ in order to eliminate artifact phases.This would therefore seem to represent a lower bound on the truncation length scale necessary to obtain a solid phase diagram in qualitative agreement with that of the full LJ potential † .
In the case of the LJL potential of Wang et al, Eq. 1, we find (see fig. 5) that the number of stable solid phases occurring at T = 0 can similarly be reduced by increasing the range parameter r c , such that for r c ≳ 7σ only a single hcp-fcc transition remains.However, as fig. 1 shows, for such large values of the range parameter, the potential is no longer similar to the LJ potential.This latter point exposes a significant disadvantage of the LJL potential, Eq. ( 1): whilst for a truncated LJ system one can investigate (and systematically correct for 5,6 ) the effects of the finite potential range when modelling a given physical system, this seems impossible for the LJL potential because the range parameter controls not only the truncation distance but also the entire potential shape.Thus while the LJL potential may be perfectly adequate in situations where one wishes to employ a simple truncated attractive interparticle interaction, it comes with a caveat on its use for crystalline solid phases since it may be unable to represent in a qualitatively correct form the phase behaviour of real atomic sys- † Note, however, the finding of Jackson et al 5,6 that by treating the ground state exactly (ie in untruncated form) so that the effects of the potential truncation are confined to the fluctuation spectrum, one can obtain a finite-temperature phase diagram close to that of the full LJ potential for shorter truncation distances.

1 2
have proposed a new class of interparticle potential for use in molecular simulation.The new potential, which has received considerable attention, takes the form φ (r) ≡ for r ≤ r c , 0 for r > r c .

Fig. 2
Fig.2Phase diagram in the density-temperature plane of the LJL system, eq.(1) with r c = 2.0σ .Blue dots are phase boundaries reported in Wang et al1 , symbols are the present work.Green crosses denote liquid-vapour coexistence densities calculated from GEMC.Our LSMC estimates of the reentrant hcp-fcc coexistence boundaries are denoted by + and were obtained in the reduced temperature range 0 < T ≤ 1 and for the particle numbers N shown in the legend.Note the complete disparity with the single hcp-fcc boundary calculated by Wang et al, as marked by blue dots in the solid region.In all cases, the difference in the coexisting hcp and fcc densities is invisibly small on the scale of the graph.Uncertainties for the coexistence densities were calculated as described in the SM 2 and unless indicated by an error bar, do not exceed the symbol size.

Fig. 4
Fig.4Dependence on number density ρ of the ground state energy difference per particle between hcp and fcc crystalline phases as calculated for the truncated and shifted LJ potential at a selection of cutoff values.Also shown for comparison is the limiting case for the full (untruncated) LJ potential17 .Note how the results for the finite cutoff approach the limiting case for r c ≳ 7.0σ .