{'fraction_non_alphanumeric': 0.028067078552515445, 'fraction_numerical': 0.008914386584289497, 'mean_word_length': 4.966824644549763, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Generalized ensemble method applied to study systems with strong first order transitions

At strong first-order phase transitions, the entropy versus energy or, at constant pressure, enthalpy, exhibits convex behavior, and the statistical temperature curve correspondingly exhibits an S-loop or back-bending. In the canonical and isothermal-isobaric ensembles, with temperature as the control variable, the probability density functions become bimodal with peaks localized outside of the S-loop region. Inside, states are unstable, and as a result simulation of equilibrium phase coexistence becomes impossible. To overcome this problem, a method was proposed by Kim, Keyes and Straub [1], where optimally designed generalized ensemble sampling was combined with replica exchange, and denoted generalized replica exchange method (gREM). This new technique uses parametrized effective sampling weights that lead to a unimodal energy distribution, transforming unstable states into stable ones. In the present study, the gREM, originally developed as a Monte Carlo algorithm, was implemented to work with molecular dynamics in an isobaric ensemble and coded into LAMMPS, a highly optimized open source molecular simulation package. The method is illustrated in a study of the very strong solid/liquid transition in water.


Introduction
First order phase transitions (solid↔liquid, liquid↔gas, and solid↔gas) are defined by a discontinuity in the free energy with respect to density. The phase transition starts in one phase and ends up in another, passing through coexisting region with both phases present. It involves heat absorption or release, the so called latent heat, where the temperature does not change with added heat until the phase transition is completed. Due to their ubiquitous presence, phase transitions have been extensively studied. However, computer simulation of phase transitions is a challenging task and here we focus on the source of the problem and applications of a proposed solution.
The problem in simulating first order phase transitions centers on the phase coexistence region and it is caused by unusual behaviour of entropy, S [1,2,3,4,5,6,7]. Entropy is a function that monotonically increases with energy or enthalpy, but in the transition region its second derivative changes sign twice, and so S changes curvature from concave to convex and back to concave, forming "convex intruder". Consequently, the statistical temperature, T S , also does not monotonically grow, but instead exhibits a back-bending, called an S-loop, as shown in black in Fig. 1.  Color dashed lines represents three temperatures chosen for isothermal simulations, with blue one as transition temperature. Below: canonical probability density function (PDF) for equilibrium transition temperature indicated by blue line on the plot above, shows two peaks, each corresponding to a pure phase system, with a region in between that is not sampled.
Isothermal simulations (either as NVT or NPT) equilibrate readily in regions where only one phase is present, and the system is well sampled satisfying the ergodic hypothesis. Such ensembles correspond to red and green dots on Fig. 1, where red and green lines representing chosen temperatures intersect the statistical temperature. For both cases, below and above transition temperature, the constant temperature defining the ensemble intersects the statistical temperature only once. The intersection locates the most probable energy and minimum in free energy.
The situation changes when one would like to study the transition region and sample configurations with coexisting phases. A constant temperature in the vicinity of the transition temperature crosses the statistical temperature three times (see blue line on Fig. 1). Two filled circles correspond to minima in free energy where only one phase is present in the system, and maxima in the energy probability distribution function (PDF). The result is bimodal behavior of the PDF (see lower panel of Fig. 1). The empty circle in the middle indicates a barrier in free energy profile that might lead to broken ergodicity. Depending on the barrier height, sampling in the phase coexisting region is either inefficient or, in fact, impossible, with system in simulations stuck on either side of the barrier in minima indicated by filled blue circles.

The gREM
While constant-temperature ensembles are the most widely used, we would like to suggest a generalized ensemble that is characterized by a temperature function depending on energy, T g (E) (or H in case of constant pressure simulations) [1]. The location of the maxima of the PDF is given by the generalized relation, The gREM utilizes generalized ensembles with a linear T g (E), with negative slope, η, chosen to be large enough so that Eq. 2 has a single solution in the coexistence region, see orange lines on Fig. 2. E 0 is a reference energy. The result is an easilysimulated unimodal PDF in the transition region [1]. Replicas are selected with varying intercepts of the T g (E) function, λ α , such that the coexistence region is spanned. Replicas exchange configurations with an acceptance probability, Acc, calculated to preserve detailed balance of the generalized ensemble, where ω α is the gREM sampling weight of the α-th replica and is calculated as follows: dE .
The original gREM method employed Monte Carlo, while here the approach was combined with the more efficient molecular dynamics method. The method was implemented by our group as a fix type module in LAMMPS [8,9]. The required sampling is obtained from isothermalisobaric MD by fixing the temperature at a constant "kinetic temperature", T 0 , and scaling the forces by the factor, s(H) = T 0 /T g (H).
The new MD gREM algorithm was applied to pure water. A cubic box containing 3337 water molecules in configuration of hexagonal ice was initially equilibrated in the isobaric-isothermal ensemble at 200 K and 1 atm, and then used as the starting configuration in each generalized replica. Periodic boundary conditions were applied. The kinetic temperature was chosen as 330 K and pressure was kept as 1 atm using Nose-Hoover style non-Hamiltonian equations of motion [10,11,12,13] as implemented in LAMMPS.
Water interactions were modeled by Stillinger-Weber potential [14] with parameters adjusted by Molinero [15] and denoted the mW model. In this approach each water molecule is represented as a single bead interacting via two-and three-body potential terms. Such coarse-graining reduces computational time while retaining a good reproduction of the properties of water [15].
Data gathered in all replicas were combined using the statistical temperature weighted histogram analysis method (ST-WHAM) [16] that is efficient and iteration-free method in contrast to the usual WHAM method [17].

Solid↔liquid transition in water
Three of the most important properties of the gREM ensemble are that (1) it can successfully sample the phase coexisting region and overcome barriers, (2) though structures in this region are unstable in canonical simulations, the gREM ensemble stabilizes them and allows for equilibration [1] and (3) by using λ to control the intersection of T g (H) with T S (H), enthalpy, not temperature, becomes the control variable for the simulation.
Characteristic structures observed around the solid↔liquid transition in water are shown in Fig. 3. The lowest enthalpy structure corresponds to hexagonal ice occupying the entire simulation box (structure denoted as A). Next, when ice starts melting, a small liquid ball (B) appears followed by a liquid cylinder (C), both surrounded by hexagonal ice. Once there is a sufficiently large amount of liquid water in the system, topology changes to liquid/ice slabs (D, E), where the interface between these two phases stays flat while the ratio of liquid to ice increases with increasing enthalpy. The next type of structures appears when there is not enough ice to form a slab and we observe instead an ice cylinder in equilibrium with bulk water (F). Finally, ice fully melts and the system consists of only bulk water (G). Topologies of the observed structures are in agreement with predictions from other studies [18]. Every change in the topology is connected to a feature in the calculated statistical temperature presented on Fig. 4. Therefore the expected S-loop in T S (H) for ice↔liquid water exhibits a complex structure in the phase coexistence region. Plateaus (or regions where T S (H) flattens) correspond to structures that belong to the same topology, while every change of topology type is in fact a transition and can be identified by an almost vertical change in T S (H).
There are two special points of interest on the T S (H) curve, namely spinodals, that define the limits of stability of both phases. The hexagonal ice in the mW model in our study can be heated up to 325.8 K before melting process starts. The experimental value of this limit equals 330±10 K [19]. In an earlier study of the mW model [20], a box of hexagonal ice was heated in canonical simulations for 100 ns until it decomposed to liquid water. The estimated limit of stability was found to be between 330 and 340 K. This value is higher compared to our result and is most likely caused the problem of overcoming the melting barrier when a liquid nucleus appears in an ice block.
The other limit of stability describes the situation where ice fully melts and changes into supercooled water at 223.4 K. The homogeneous nucleation temperature of mW water is equal to 202 K under 1 atm [15], below our estimated spinodal temperature, which is nor surprising considering fact that we started the simulation from hexagonal ice, not liquid water, and as such we monitored ice melting instead of ice nucleation. With the statistical temperature available one can apply the Maxwell construction to the inverse of T S (H) to obtain the transition temperature. In this procedure one is looking for a temperature that yields equal areas above and below it. The calculated equilibrium temperature of ice melting in mW model is equal to 269.1 K, compared to temperature of 274.6 K obtained by Molinero [15] using the phase coexistence method. Interestingly, the calculated melting temperature corresponds to a large plateau on the statistical temperature where structures consist of liquid and ice slabs with increasing amount of liquid towards higher enthalpies. This topology is identical to those used in the phase coexistence method where half of the simulation box is filled with one phase and the other half with the second phase. Location of this plateau can in fact be used to estimate the transition temperature.
In the present paper we showed preliminary results from the new MD version of an enhanced sampling algorithm, gREM, applied to the strong first order solid↔liquid transition in water modeled by the mW potential. We successfully calculated the statistical temperature for this transition as well as the transition temperature. We observed a series of structures with coexisting phases in several distinct topologies. Though usage of the gREM simulation method each of the observed structures was in equilibrium, in contrast to results from canonical simulations where the transition region is barely sampled.