The generalized Replica Exchange Method (gREM) was applied to study a solid-liquid phase transition in a nanoconfined bilayer water system using the monatomic water (mW) model. Exploiting optimally designed non-Boltzmann sampling weights with replica exchanges, gREM enables an effective sampling of configurations that are metastable or unstable in the canonical ensemble via successive unimodal energy distributions across phase transition regions, often characterized by S-loop or backbending in the statistical temperature. Extensive gREM simulations combined with Statistical Temperature Weighted Histogram Analysis Method (ST-WHAM) for nanoconfined mW water at various densities provide a comprehensive characterization of diverse thermodynamic and structural properties intrinsic to phase transitions. Graph representation of minimized structures of bilayer water systems determined by the basin-hopping global optimization revealed heterogeneous ice structures composed of pentagons, hexagons, and heptagons, consistent with an increasingly ordered solid phase with decreasing density. Apparent crossover from a first-order solid-liquid transition to a continuous one in nanoconfined mW water with increasing density of the system was observed in terms of a diminishing S-loop in the statistical temperature, smooth variation of internal energies and heat capacities, and a characteristic variation of lateral radial distribution functions, and transverse density profiles across transition regions.

## I. INTRODUCTION

The phase behavior of bulk water has received extensive interest due to the rich complexity of structures characterizing liquid, solid and clusters.^{1–4} Nanoscale confined water adds a new dimension of phase behavior and has generated intense interest^{5–10} due to its relevance in biology and materials science. In a nano-confined water system, the separation between two infinite parallel plates serves as a natural variable to tune the degree of confinement. As the interplate separation increases, monolayer, bilayer, three-layer, and four-layer confined water systems can be produced.^{8,11–14} Molinero *et*. *al*. found that, of the series of confined ices ranging up to four layers, only the buckled bilayer has an anomalously high melting point.^{8} Bilayer water can form various crystals, quasicrystals, and amorphous structures, including hexagonal ice, pure pentagonal ice, mixed hexagonal and pentagonal ice, and dodecagonal quasicrystals.^{9} The transitions from liquid to various crystal and quasicrystal states in confined bilayer water were shown to be first-order, based on the sharp drop in the potential energy and discontinuity in the diffusion coefficient.^{9}

First-order phase transitions in finite size systems have a unique feature in the statistical temperature, *T*_{S}(*E*) = [∂ln Ω(*E*)/∂*E*]^{−1}, Ω(*E*) being the density of states, referred to as an S-loop or backbending, through which the temperature decreases upon absorbing energy in the region of metastable and unstable states. The depletion of phase-coexistent configurations stemming from the free energy penalty for forming interfaces^{15} is manifested in a bimodal structure in the energy distribution and the energy states between the two peaks are intrinsically unstable in the canonical ensemble.^{15–19} One effective way to enhance sampling in phase-coexistence regions is to replace a canonical ensemble sampling by a generalized ensemble sampling with non-Boltzmann sampling weights. The recently developed generalized Replica Exchange Method (gREM)^{20} optimally combines multiple non-Boltzmann ensemble samplings with a popular replica exchange scheme.^{21–33} The generalized ensemble sampling weights are determined from tailored effective temperatures through an inverse mapping strategy. The mapping is equivalent to umbrella sampling for a number of energy windows, with a “thermometer” in each window. Since its development, gREM has been successfully applied to study phase transitions in Potts spin systems, an adapted Dzutugov model, Lennard-Jones fluid, and bulk water.^{20,34–36} Here we study the solid-liquid phase transition in bilayer confined water as a demonstration of the utility of the generalized Replica Exchange Method.

## II. METHODS AND MATERIALS

### A. Generalized exchange method

In the gREM simulation, each replica α, (α = 1,⋯, *M*), is associated with an effective temperature *T*_{α}(*E*; λ_{α}), λ_{α} being a set of relevant parameters, and samples a configurational space with the generalized ensemble weight *W*_{α}(*E*; λ_{α}), which is determined from the effective temperature through the inverse mapping

where

where the control parameter γ_{0} is the constant slope, *E*_{0} is a constant in the relevant energy range, and λ_{α} is the *T*-intercept at a chosen *E*_{0}. The optimal parameter values of λ_{α} and γ_{0} can be determined by performing short canonical Monte Carlo (MC) simulations at the lowest and highest temperatures of interest (see Ref. 20 for details).

The linear effective temperature of Eq. (2) produces a generalized ensemble sampling weight

which governs the trial moves within one replica and replica exchanges between neighboring replicas. The acceptance probability of MC trial move in configuration space within replica α is

The acceptance ratio for replica exchange between neighboring replicas α and α + 1 is

where

*N*single particle trial moves).

### B. Statistical temperature weighted histogram analysis method (ST-WHAM)

Since each replica in gREM samples a configurational space with the non-canonical sampling weight *W*_{α}(*E*) all simulation results should be combined to estimate correct canonical averages via a reweighting. The Statistical Temperature Weighted Histogram Analysis Method (ST-WHAM)^{37–39} is utilized to combine multiple generalized ensemble runs to determine the entropy estimate *S*(*E*) = *k*_{B}lnΩ(*E*). Once the density of states is determined all canonical thermodynamic averages are determined for any temperature by the reweighting technique. Contrary to the conventional Weighted Histogram Analysis Method (WHAM),^{40} ST-WHAM does not require an iterative process to determine relevant partition functions, and directly determines the inverse statistical temperature estimate, β_{S}(*E*) = 1/*T*_{S}(*E*), given *W*_{α}(*E*) and associated histograms *H*_{α}(*E*).

The entropy estimate *S*(*E*) or corresponding density of states estimate is computed from a numerical integration of β_{S}(*E*) = ∂*S*(*E*)/∂*E*, enabling a substantial acceleration of the data analysis without loss in accuracy. The ST-WHAM estimate for the inverse statistical temperature is obtained as

where *H*_{α}(*E*) is the energy histogram in replica α, *f*_{α}(*E*) = *H*_{α}(*E*)/∑_{α}*H*_{α}(*E*) is the simulated histogram fraction. Due to the rapid variation of β_{S} for small *E* a direct integration of β_{S} is not desirable so we first approximate the statistical temperature *T*_{S}(*E*) on equally spaced energy grids, allowing a piecewise analytical integration.^{37} All canonical thermodynamic properties are completely determined with the determined *S*(*E*).

The Helmholtz free energy density at a given temperature *T* is calculated by *F*_{T}(*E*) = *E* − *TS*(*E*) and the reweighted probability density function is given by

and the canonical heat capacity is estimated through the energy fluctuation as

### C. Monatomic water (mW) model

We employed the monatomic water model (mW), which is a coarse-grained model that represents a water molecule as an intermediate element between carbon and silicon.^{41} This potential reproduces the structural, thermodynamic, and dynamic properties of liquid water with comparable or better accuracy than the most popular atomistic water models at much less computational cost. mW has been applied in the study of pure bulk water^{3,42,43} and nanoconfined water,^{44,45} as well as biological water^{46} and clathrate hydrates.^{47}

The mW model consists of a sum of pairwise two-body and three-body interactions described by a functional form of the Stillinger-Weber potential.^{48} The Stillinger-Weber potential was originally developed to model silicon, so most parameters in the silicon potential can be used for the mW model without change. Only three parameters are changed to adapted to the monatomic water model, which are the tetrahedrality λ = 23.15, the diameter σ_{mW} = 2.3925 Å, and energy scale ε = 25.895 kJ/mol. In this study, the water-like molecules were confined between two parallel featureless hydrophobic plates separated by a distance *D*. The water-wall interaction was modeled by a Lennard-Jones 9-3 potential as

where δ*z* is the distance in *z* between the water molecule and the plate, σ_{p} = 3.56 Å, and

^{8}

The separation between plates, *D*, was chosen to be 8.5 Å because previous study^{8} showed that 8.5 Å is the optimal distance for which the melting temperature of bilayer ice reaches a maximum in the mW model. Periodic boundary conditions were applied in the parallel (*x* and *y*) directions to mimic infinite parallel plates. Total ten simulations were performed with varying *N* =256, 576, and 800 water molecules at six different densities ranging from 1.08 g cm^{−3} up to 1.23 g cm^{−3}. We have chosen to consider the system at constant volume, following the design of a previous study of phase transitions in nanoconfined water films.^{5} The simulation parameters including numbers of particles, lengths of plates, effective densities, effective temperature parameters in Eq. (2) for *g*REM, and total simulation times were summarized in Table I. The effective densities are estimated using the method in the previous study^{11} as

*N*is the number of water molecules,

*m*is the mass of a water molecule, and

*L*

_{x}and

*L*

_{y}are the lengths of the plates. Here

System . | N . | L_{x} (Å)
. | L_{y} (Å)
. | ρ (g cm^{−3})
. | E_{0} (kJ/mol)
. | λ_{1} (K)
. | λ_{32} (K)
. | γ_{0} (mol K $\rm {kJ^{-1}}$ $ kJ \u22121$)
. |
---|---|---|---|---|---|---|---|---|

System 1 | 256 | 36.0 | 31.2 | 1.2328 | −40.95 | 220.00 | 420.00 | −36.5679 |

System 2 | 256 | 34.5 | 34.5 | 1.1643 | −41.61 | 200.13 | 380.25 | −18.2764 |

System 3 | 256 | 34.8 | 34.8 | 1.1443 | −42.18 | 200.13 | 380.25 | −14.0019 |

System 4 | 256 | 35.0 | 35.0 | 1.1312 | −41.91 | 200.13 | 380.25 | −15.5259 |

System 5 | 256 | 35.5 | 35.5 | 1.0996 | −42.76 | 200.13 | 440.29 | −22.6551 |

System 6 | 256 | 35.8 | 35.8 | 1.0812 | −42.23 | 200.13 | 440.29 | −24.1468 |

System 7 | 576 | 56.1 | 48.6 | 1.1643 | −42.01 | 220.00 | 340.00 | −4.9615 |

System 8 | 576 | 57.7 | 50.0 | 1.0812 | −42.69 | 240.00 | 440.00 | −10.2134 |

System 9 | 800 | 61.5 | 61.5 | 1.1643 | −42.07 | 220.00 | 370.00 | −5.0381 |

System 10 | 800 | 63.3 | 63.3 | 1.0812 | −42.68 | 220.00 | 280.00 | −5.5218 |

System . | N . | L_{x} (Å)
. | L_{y} (Å)
. | ρ (g cm^{−3})
. | E_{0} (kJ/mol)
. | λ_{1} (K)
. | λ_{32} (K)
. | γ_{0} (mol K $\rm {kJ^{-1}}$ $ kJ \u22121$)
. |
---|---|---|---|---|---|---|---|---|

System 1 | 256 | 36.0 | 31.2 | 1.2328 | −40.95 | 220.00 | 420.00 | −36.5679 |

System 2 | 256 | 34.5 | 34.5 | 1.1643 | −41.61 | 200.13 | 380.25 | −18.2764 |

System 3 | 256 | 34.8 | 34.8 | 1.1443 | −42.18 | 200.13 | 380.25 | −14.0019 |

System 4 | 256 | 35.0 | 35.0 | 1.1312 | −41.91 | 200.13 | 380.25 | −15.5259 |

System 5 | 256 | 35.5 | 35.5 | 1.0996 | −42.76 | 200.13 | 440.29 | −22.6551 |

System 6 | 256 | 35.8 | 35.8 | 1.0812 | −42.23 | 200.13 | 440.29 | −24.1468 |

System 7 | 576 | 56.1 | 48.6 | 1.1643 | −42.01 | 220.00 | 340.00 | −4.9615 |

System 8 | 576 | 57.7 | 50.0 | 1.0812 | −42.69 | 240.00 | 440.00 | −10.2134 |

System 9 | 800 | 61.5 | 61.5 | 1.1643 | −42.07 | 220.00 | 370.00 | −5.0381 |

System 10 | 800 | 63.3 | 63.3 | 1.0812 | −42.68 | 220.00 | 280.00 | −5.5218 |

### D. Water structure analysis

To characterize structural properties of water molecules the lateral radial distribution function (RDF) *g*_{xy}(*r*)^{11} is computed as a function of the lateral position *r*_{xy} parallel to the confining plates as

where *V* is the volume, *r*_{ij} is the lateral distance between coarse-grained molecules *i* and *j*, *z* is the *z* coordinate, δ*z* = 2 Å, and δ(*x*) is the Dirac δ function. Note that the Heaviside function θ(*x*) restricts the sum to pairs within the same layer.

The lateral static structure factor *S*(*q*) is the Fourier transform of the lateral RDF^{49,50} determined as

where *q* = 2π*k*/*L*, *k* being an integer number from 1 to *N*, and *L* is the length of the simulation box.

### E. Basin-hopping global optimization

Basin-hopping (BH) global optimization,^{51,52} as implemented in the GMIN^{53} package, was employed to explore a simulated potential energy landscape. The BH scheme used in this work is as follows:

Apply a random Cartesian displacement to the initial coordinates

$\bm {r}_{i}$$ri$;Find the local minimum

$\bm {r}_{n}$$rn$ from the perturbed coordinates$\bm {r}_{i}^\prime$$ri\u2032$;Accept a trial move to the new configuration

$\bm {r}_{n}$$rn$ with the probability of$p(i \rightarrow n) = \mathrm{min}[1,e^{-\beta (E_{n}-E_{i})} ]$$p(i\u2192n)= min [1,e\u2212\beta (En\u2212Ei)]$ where*E*_{i}and*E*_{n}are the energies at the initial and new local minima, and β = 1/*kT*.

6 × 10^{3} BH steps were run for each starting structure. At each step, random Cartesian displacements up to 0.8 Å were applied to each particle. The temperature parameter *T* was fixed at

^{54,55}The root-mean-square gradient of the local minima was converged to

## III. RESULTS AND DISCUSSION

With an optimal choice of effective temperatures *T*_{α}(*E*; λ_{α}) unstable or metastable energy states in the canonical ensemble corresponding to the S-loop region in *T*_{S}(*E*) are transformed into stable ones in the generalized ensemble, resulting in a unimodal probability distribution function (PDF). The necessary and sufficient condition on *T*_{α}(*E*; λ_{α}) to achieve a unimodal PDF is derived by examining a local stability of the generalized free energy density, β*F*_{α}(*E*) = *w*_{α}(*E*) − *S*(*E*), at extrema

where

*T*

_{α}(

*E*) and

*T*

_{S}(

*E*). The stability condition is

where

*E*. Exploiting the linear effective temperatures in Eq. (2) a unimodal PDF is ensured by

*P*

_{α}(

*E*) up to second order at

where

Fig. 1(a) demonstrates that the linear effective temperatures *T*_{α}(*E*) form unique crossing points

*T*

_{S}(

*E*) across the transition region, where

*T*

_{S}(

*E*) displays the S-loop. Setting

_{0}= −24.15 mol K kJ

^{−1}, replica-dependent effective temperatures,

*T*

_{α}(

*E*; λ

_{α}) = λ

_{α}+ γ

_{0}(

*E*−

*E*

_{0}), were assigned using evenly spaced λ

_{α}between λ

_{1}= 200.13 K and λ

_{32}= 440.29 K. As illustrated in Fig. 1(b) the linear effective temperatures with negative slope,

*P*

_{α}(

*E*) centered at

In order to characterize structural differences between liquid and solid phases near the transition temperature, we compared the lateral radial distribution functions *g*_{xy}(*r*) and structure factors *S*(*q*) of two representative replicas 18 and 25 (see Fig. 2). Note that the effective temperatures, *T*_{18}(*E*) and *T*_{25}(*E*), form crossing points with *T*_{S}(*E*) at

*q*≈ 2, a sharp first peak, and a split second peak, while

*S*(

*q*) in replica 25 lacks the prepeak and the split in the second peak.

To examine the layering effects from the confinement we computed the transverse density profile ρ_{z} (TDP) along the *z* direction in Fig. 2(c). The TDPs for replicas 18 and 25 have two pronounced symmetric peaks with respect to the slit center (*z* = 0), confirming that two layers of water molecules are confined in the nanoslit. However, the TDP of replica 25 has a much wider distribution with smaller peaks at each layer center and with non-zero density values across the slit center, while the TDP of replica 18 showed a vanishing density at the slit center. This implies that water molecules in replica 18 can easily move between two layers across the slit center, while transverse movements of water molecules in replica 25 are highly restricted, resulting in a strong layering effect. The differences in RDFs, structure factors, and transverse density profiles reveal the solid-like and liquid-like characteristics of the configurations in replica 18 and 25, respectively, manifesting the coexistence of two structurally distinct states in the canonical ensemble.

It is instructive to examine the free energy density as a function of energy (see Fig. 3). After a sufficiently long simulation with gREM, multiple replica simulations are optimally combined to produce the entropy estimate, *S*(*E*), via ST-WHAM.^{37} Once the entropy is determined, canonical thermodynamic properties, including internal energy *E*(*T*) and heat capacity *C*_{v}(*T*), can be calculated as in Eqs. (7) and (8). In contrast to the statistical temperature, *T*_{S}(*E*), characteristic of the microcanonical ensemble in Fig. 1(a), the canonical ensemble average, *E*(*T*), monotonically decreases as *T* decreases and shows a sharp drop across the phase transition region in Fig. 3(a), which corresponds to the peak in the heat capacity *C*_{v} around the transition temperature *T*_{m} = 288.7 K. The free energy density in Fig. 3(b) at *T*_{m}, *F*(*E*, *T*_{m}) = *E* − *T*_{m}*S*(*E*), exhibits two local minima at *E*_{1} = −38.1 and *E*_{3} = −40.5 kJ mol^{−1}, and one local maximum at *E*_{2} = −39.5 kJ mol^{−1}, resulting in the bimodal structure in *P*_{T}(*E*) in contrast to the unimodal PDFs in gREM across the phase transition region.^{15,17–19}

We also performed gREM simulations at several densities, ρ, varying the system size *N* = 256, 576, and 800 (see Table I for the detailed description of systems), in order to investigate the effect of density on the phase transition in nanoconfined mW water. The reweighted internal energies and heat capacities in the canonical ensemble for System 1–6 are illustrated in Figs. 4(a) and 4(b), respectively. Interestingly, the decrease in internal energy (latent heat) across the transition region from liquid to solid is monotonically decreasing with increasing ρ. This means that the character of the first-order phase transition between liquid and solid in nanoconfined mW water becomes much weaker in the high density regime, in which freezing of water may occur through a continuous transformation.

To investigate structural changes associated with the phase transition at various ρ both profiles of RDF and TDP are examined across replicas near the transition region for System 2 and 6 corresponding to ρ = 1.0812 g cm^{−3} and 1.1643 g cm^{−3}, respectively. As illustrated in Fig. 5(a) the variation of *g*_{xy}(*r*) at high density is almost continuous across replicas, while at low density *g*_{xy}(*r*) is well segregated as clearly seen in first and second peaks in Fig. 5(b). The continuous variation of *g*_{xy}(*r*) at higher density is most apparent in the long-range order. For the lower density system, there are clearly two sets of replicas — those that show substantial long-range order and those that do not — with an abrupt transition between the two. The substantial differences in the behavior of *g*_{xy}(*r*) for the high and low density systems are also apparent in the transverse density profiles in Figs. 6(a) and 6(b). At higher density (System 2), the peaks at each layer center are continuously rising across replicas from liquid to solid phase with a marginal density at the slit center, implying that transverse movements of waters are highly restricted in both solid and liquid phase. At lower density (System 6), the profiles show higher non-zero density at the slit center in liquid phase replicas, implying that water molecules can freely move between the two layers in the liquid phase. These observations are consistent with the Stanley group's conjecture regarding the possibility of a “crossover” from first-order to continuous transition between the liquid and solid state with increasing density.^{5} We note that this crossover behavior was not observed in previous work by Molinero and co-workers.^{9} However, that work was based on an NP_{xy}T ensemble while we have used an NVT ensemble, as was done in the study of the Stanley group.^{5} The Stanley group also did note crossover behavior in an NP_{xy}T ensemble simulation as well.^{5}

The crossover behavior from a first-order transition to a continuous one with increasing density is also seen in the profile of the approximate statistical temperatures, *T*_{S}(*E*), in Fig. 7, in which most probable energies,

*T*

_{α}(

*E*) and

*T*

_{S}(

*E*), were plotted with

*N*= 256, 576, and 800, while maintaining the density of systems at ρ

_{1}= 1.1643 g cm

^{−3}in Fig. 7(a) and ρ

_{2}= 1.0812 g cm

^{−3}in Fig. 7(b). It is apparent that the statistical temperature estimates from several independent gREM simulations at different

*N*are in overall good agreement, implying that finite size effects are marginal in our simulations. Interestingly,

*T*

_{S}(

*E*) for low density systems shows a clear S-loop characteristic of first-order phase transitions, while it varies monotonically without noticeable structure across the transition region in high density systems.

To characterize the global minima structures for System 1–10, basin-hopping global optimization^{51,52,56} was applied starting from equilibrium configurations of the lowest energy replica in each gREM simulation using the GMIN package, yielding representative snap shots in Figs. 8 and 9. The putative global minimum of System 8 at ρ = 1.0812 g cm^{−3} for *N* = 576 corresponds to a regular hexagonal ice form, while in all other systems the ground states were imperfect hexagonal ices with some defects. While this leads to a less regular global minimum structure, our results show that it does not impact the character of the thermodynamics of the solid-to-liquid phase change. The minimized structures for these systems show varying composition of polygonal rings in each layer. We made use of graph representation of the structures, in which vertices represent particles, and an edge between two vertices represents a nearest neighbor contact between two particles with a cut-off of 3.35 Å. A ring is a connected graph in which each vertex shares exactly two edges with other vertices in the ring. For the bilayer structures in this work, we counted rings within individual layers. We calculated a composition of *n*-membered rings on the basis of 100 minimized structures, only one of which was the global minimum.

Table II lists the percentage distributions of *n*-membered rings (*n* = 4, 5, 6, 7) for System 1–10, showing an overall trend that the composition of pentagonal rings drops with decreasing density and the combined composition of hexagonal and heptagonal rings increases with decreasing density. This is consistent with a view of a transition from a more disordered solid at higher density to more ordered solid at lower density. Interestingly, the ring counts for the 100 structures are quite similar to those of the global minimum structure in all cases.

Ring size . | n = 4 (%) . | n = 5 (%) . | n = 6 (%) . | n = 7 (%) . | ρ (g cm^{−3})
. |
---|---|---|---|---|---|

System 1 | 0.5 (0) | 93.3 (92.8) | 6.2 (7.2) | 0 (0) | 1.2328 |

System 2 | 3.2 (8.3) | 25.9 (27.7) | 68.9 (63.9) | 2.0 (0) | 1.1643 |

System 3 | 2.4 (1.6) | 15.8 (9.7) | 77.1 (82.3) | 4.7 (6.4) | 1.1443 |

System 4 | 2.6 (3.0) | 14.6 (14.9) | 78.6 (79.1) | 4.2 (3.0) | 1.1312 |

System 5 | 0.1 (0) | 13.4 (10.9) | 73.5 (78.1) | 13.0 (10.9) | 1.0996 |

System 6 | 1.2 (0) | 17.8 (16.1) | 65.8 (74.2) | 15.2 (9.7) | 1.0812 |

System 7 | 3.4 (2.6) | 25.7 (27.0) | 63.9 (63.2) | 7.0 (7.2) | 1.1643 |

System 8 | 0.03 (0) | 2.8 (0) | 95.2 (99.3) | 2.0 (0.7) | 1.0812 |

System 9 | 2.8 (2.1) | 32.9 (31.5) | 57.7 (59.4) | 6.6 (7.0) | 1.1643 |

System 10 | 0.09 (0) | 12.7 (12.6) | 76.0 (76.3) | 11.2 (11.1) | 1.0812 |

Ring size . | n = 4 (%) . | n = 5 (%) . | n = 6 (%) . | n = 7 (%) . | ρ (g cm^{−3})
. |
---|---|---|---|---|---|

System 1 | 0.5 (0) | 93.3 (92.8) | 6.2 (7.2) | 0 (0) | 1.2328 |

System 2 | 3.2 (8.3) | 25.9 (27.7) | 68.9 (63.9) | 2.0 (0) | 1.1643 |

System 3 | 2.4 (1.6) | 15.8 (9.7) | 77.1 (82.3) | 4.7 (6.4) | 1.1443 |

System 4 | 2.6 (3.0) | 14.6 (14.9) | 78.6 (79.1) | 4.2 (3.0) | 1.1312 |

System 5 | 0.1 (0) | 13.4 (10.9) | 73.5 (78.1) | 13.0 (10.9) | 1.0996 |

System 6 | 1.2 (0) | 17.8 (16.1) | 65.8 (74.2) | 15.2 (9.7) | 1.0812 |

System 7 | 3.4 (2.6) | 25.7 (27.0) | 63.9 (63.2) | 7.0 (7.2) | 1.1643 |

System 8 | 0.03 (0) | 2.8 (0) | 95.2 (99.3) | 2.0 (0.7) | 1.0812 |

System 9 | 2.8 (2.1) | 32.9 (31.5) | 57.7 (59.4) | 6.6 (7.0) | 1.1643 |

System 10 | 0.09 (0) | 12.7 (12.6) | 76.0 (76.3) | 11.2 (11.1) | 1.0812 |

## IV. CONCLUSION

In summary, the effectiveness of the gREM for the study of first-order phase transitions was demonstrated using nanoconfined mW systems at diverse simulation conditions. By employing efficient replica exchanges across unimodal energy distributions spanning phase transition regions, gREM attains a comprehensive sampling for metastable and unstable configurations intrinsic to the phase transitions, regions only rarely accessed in canonical ensemble simulations due to the presence of high free energy barriers.

The graph representation analysis for low-lying solid-phase structures, determined by basin-hopping global optimization using the GMIN package, shows heterogeneous crystalline structures composed of different compositions of pentagonal, hexagonal, and heptagonal rings depending on the simulation conditions. Interestingly, minimized structures of low density systems exhibiting the first-order transition character predominantly contained hexagonal rings, while the percentage of pentagonal rings gradually increases in high density systems showing a continuous transitions.

With an optimal integration of multiple replica simulations of gREM using the ST-WHAM we characterized detailed thermodynamic and structural properties of solid-liquid phase transitions in nanoconfined mW water at various densities. The first-order phase transitions between liquid and hexagonal ice in the low density regime was explicitly illustrated in terms of the existence of an S-loop or backbending in the statistical temperature, *T*_{S}(*E*), sharp drops in internal energies, pronounced peaks in heat capacities, and quantitative differences in profiles of RDF and TDP between solid and liquid phases. Interestingly, characteristics of a first-order transition in nanoconfined mW water were observed to be much weaker at high density, accompanied by an absence of backbending character in *T*_{S}(*E*), gradual changes of both internal energies and heat capacities, and continuous variations of both RDFs and TDPs across replicas sampling the phase transition regions. This suggests the possibility of crossover from a first-order transition to a continuous one in quasi-two-dimensional confined water systems depending on the density, which was recently hypothesized based on a full atomistic simulation of TIP5P water.^{5} Additional studies of nanoconfined water using the gREM approach in alternative ensembles, including the *NP*_{xy}*T* ensemble, may provide additional insight into the possible validity of this intriguing conjecture.

## ACKNOWLEDGMENTS

We are grateful to the National Science Foundation (Grant Nos. CHE-1114676 and CHE-1362524) for the generous support of our research.