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.

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 interest5–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, TS(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 interfaces15 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.

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

\begin{equation}T_{\alpha }(E;\lambda _{\alpha })=[\partial w_{\alpha }/\partial E]^{-1},\end{equation}
Tα(E;λα)=[wα/E]1,
(1)

where

$w_{\alpha }=-\textrm {ln}W_{\alpha }$
wα= ln Wα is the generalized effective potential. The effective temperatures are conveniently parameterized using linear functions with a uniform slope for all replicas as

\begin{equation}T_{\alpha }(E;\lambda _{\alpha })=\lambda _{\alpha }+\gamma _0(E-E_0),\end{equation}
Tα(E;λα)=λα+γ0(EE0),
(2)

where the control parameter γ0 is the constant slope, E0 is a constant in the relevant energy range, and λα is the T-intercept at a chosen E0. 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

\begin{equation}W_\alpha (E;\lambda _\alpha ) \sim [\lambda _\alpha + \gamma _0 (E-E_0)]^{-1/\gamma _0},\end{equation}
Wα(E;λα)[λα+γ0(EE0)]1/γ0,
(3)

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

\begin{eqnarray}{A}_{\rm gREM}({\bf x} \rightarrow {\bf x}^\prime ) = \textrm {min}[1,e^{{w}_\alpha (E({\bf x}))-{w}_\alpha (E({\bf x^\prime }))}].\end{eqnarray}
A gREM (xx)= min [1,ewα(E(x))wα(E(x))].
(4)

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

\begin{eqnarray}A_{\rm gREM}(\alpha ;{\bf x x}^\prime )= \textrm {min} [1,\exp (\Delta _{\alpha }^x)],\end{eqnarray}
A gREM (α;xx)= min [1,exp(Δαx)],
(5)

where

$\Delta _{\alpha }^{\bf x} \!\!=\!\! w_{\alpha \!+\!1}(E({\bf x^\prime }))\!-\!w_{\alpha \!+\!1}(E({\bf x}))\!+\!w_\alpha (E({\bf x}))\!-\!w_\alpha (E({\bf x^\prime }))$
Δαx=wα+1(E(x))wα+1(E(x))+wα(E(x))wα(E(x))⁠. After all replicas complete one MC sweep, replica exchanges were attempted between neighboring replicas with a random selection of α (a sweep being defined by N single particle trial moves).

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) = kBlnΩ(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/TS(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

\begin{equation}\beta _{S}=\sum _{\alpha }f_{\alpha }(E)\left(\frac{\partial \textrm {ln}H_{\alpha }}{\partial E}-\frac{\partial \textrm {ln}W_{\alpha }}{\partial E}\right),\end{equation}
βS=αfα(E) ln HαE ln WαE,
(6)

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 TS(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 FT(E) = ETS(E) and the reweighted probability density function is given by

$P_T(E)=e^{-{F}_T(E)/T}\break=e^{S(E)-E/T}$
PT(E)=eFT(E)/T=eS(E)E/T⁠. The canonical expectation value for any quantity is computed as

\begin{equation}\langle A(T)\rangle = \frac{\int dE e^{S(E)-\beta E}A(E) }{\int dE e^{S(E)-\beta E}},\end{equation}
A(T)=dEeS(E)βEA(E)dEeS(E)βE,
(7)

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

\begin{equation}C_v(T)=\frac{\langle E(T)^2\rangle -\langle E(T)\rangle ^2}{k_BT^2}.\end{equation}
Cv(T)=E(T)2E(T)2kBT2.
(8)

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 water3,42,43 and nanoconfined water,44,45 as well as biological water46 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

\begin{equation}\phi _{{\rm LJ-93}}=\epsilon \left[\frac{2}{15}\left (\frac{\sigma _p}{\delta z}\right)^9-\left (\frac{\sigma _p}{\delta z}\right)^3\right] ,\end{equation}
ϕ LJ 93=ε215σpδz9σpδz3,
(9)

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

$\epsilon =0.569\,\,\rm {kJ\,mol^{-1}}$
ε=0.569 kJ mol 1⁠.8 

The separation between plates, D, was chosen to be 8.5  Å because previous study8 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 gREM, and total simulation times were summarized in Table I. The effective densities are estimated using the method in the previous study11 as

$\rho ={Nm}/({L_xL_yL_z^{\prime }})$
ρ=Nm/(LxLyLz)⁠, where N is the number of water molecules, m is the mass of a water molecule, and Lx and Ly are the lengths of the plates. Here
$L_z^{\prime }$
Lz
is the effective distance between the plates accessible to the water molecules confined between two hydrophobic plates, defined as
$L_z^{\prime }=D-(\sigma _{\rm {mW}}+\sigma _{p})/2=5.5237\,\text{\AA}$
Lz=D(σ mW +σp)/2=5.5237Å
.

Table I.

Simulation parameters for 10 systems including number of molecules, plate sizes, effective densities and parameters in gREM. Here Lx and Ly are the lengths of the plate in the parallel directions, Ly is fixed to 8.5  Å for all systems, ρ is the effective densities, and E0, λ1, λ32, and γ0 are parameters in Tα as shown in Eq. (2). Thirty two replicas are used in gREM simulations for all systems, and simulation length is 106 MC sweeps.

SystemNLx (Å)Ly (Å)ρ (g cm−3)E0 (kJ/mol)λ1 (K)λ32 (K)γ0 (mol K
$\rm {kJ^{-1}}$
kJ 1
)
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 
SystemNLx (Å)Ly (Å)ρ (g cm−3)E0 (kJ/mol)λ1 (K)λ32 (K)γ0 (mol K
$\rm {kJ^{-1}}$
kJ 1
)
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 

To characterize structural properties of water molecules the lateral radial distribution function (RDF) gxy(r)11 is computed as a function of the lateral position rxy parallel to the confining plates as

\begin{equation}g_{xy}(r)=\frac{1}{\rho ^2V}\sum _{i\ne j} \delta (r-r_{ij}) [1-\theta (|z_i-z_j|-\delta z/2) ],\end{equation}
gxy(r)=1ρ2Vijδ(rrij)[1θ(|zizj|δz/2)],
(10)

where V is the volume, rij 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 RDF49,50 determined as

\begin{equation}S(q)=1+2\pi \rho \int _0^L r\left(\frac{\text{sin}(qr)}{qr} \right)[g_{xy}(r)-1]dr ,\end{equation}
S(q)=1+2πρ0Lrsin(qr)qr[gxy(r)1]dr,
(11)

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

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

  1. Apply a random Cartesian displacement to the initial coordinates

    $\bm {r}_{i}$
    ri⁠;

  2. Find the local minimum

    $\bm {r}_{n}$
    rn from the perturbed coordinates
    $\bm {r}_{i}^\prime$
    ri
    ;

  3. 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(in)= min [1,eβ(EnEi)]
    where Ei and En are the energies at the initial and new local minima, and β = 1/kT.

6 × 103 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

$8.0\,\rm { kJ\,mol^{-1}}$
8.0 kJ mol 1⁠. Local optimization was performed using a modified version of Nocedal's limited memory BFGS (L-BFGS) minimizer.54,55 The root-mean-square gradient of the local minima was converged to
$10^{-4}\,\,\rm {kJ\,mol^{-1}\,{\AA}^{-1}}$
104 kJ mol 1Å1
.

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 TS(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

$E_{\alpha }^*$
Eα*⁠, determined as

\begin{equation}{T}_\alpha (E_\alpha ^*;\lambda _\alpha )=T_S(E_\alpha ^*)=T_\alpha ^*,\end{equation}
Tα(Eα*;λα)=TS(Eα*)=Tα*,
(12)

where

$E_{\alpha }^*$
Eα* corresponds to the crossing point between Tα(E) and TS(E). The stability condition is

\begin{equation}\beta {F}_{\alpha }^{\prime \prime }(E_\alpha ^*)=(\gamma _S -\gamma _\alpha )/{T_\alpha ^*}^2,\end{equation}
βFα(Eα*)=(γSγα)/Tα*2,
(13)

where

$\gamma _S=T_S^\prime (E_\alpha ^*)$
γS=TS(Eα*) and
$\gamma _\alpha =T_\alpha ^\prime (E_\alpha ^*)$
γα=Tα(Eα*)
, the prime being a differentiation with respect to E. Exploiting the linear effective temperatures in Eq. (2) a unimodal PDF is ensured by
$\gamma _\alpha (E_\alpha ^*)= \gamma _0 < \gamma _S(E_\alpha ^*)$
γα(Eα*)=γ0<γS(Eα*)
around
$E_\alpha ^*$
Eα*
. Expanding Pα(E) up to second order at
$E_{\alpha }^*$
Eα*
verifies

\begin{equation}P_\alpha (E;\gamma _0) \approx {\rm exp}\big[-{(E-E_\alpha ^*)^2}/{2\sigma _{\gamma }}\big],\end{equation}
Pα(E;γ0) exp (EEα*)2/2σγ,
(14)

where

$\sigma _\gamma = {T_\alpha ^*}^2/(\gamma _S-\gamma _0)$
σγ=Tα*2/(γSγ0)⁠, yielding a Gaussian PDF centered at
$E_{\alpha }^*$
Eα*
.

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

$E_{\alpha }^*$
Eα* with the statistical temperature TS(E) across the transition region, where TS(E) displays the S-loop. Setting
$E_0=-42.23\,\rm {kJ\,mol^{-1}}$
E0=42.23 kJ mol 1
and γ0 = −24.15 mol K kJ−1, replica-dependent effective temperatures, Tα(E; λα) = λα + γ0(EE0), 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,
$\gamma _0 < \gamma _S(E_\alpha ^*)$
γ0<γS(Eα*)
, result in Pα(E) centered at
$E_\alpha ^*$
Eα*
with a Gaussian shape, naturally bridging between ordered and disordered phases with unimodal energy distributions across the transition region. Figs. 1–3 show simulation results of System 6 in Table I.

FIG. 1.

(a) Effective temperatures Tα(E) (a set of parallel lines with negative slope) form unique crossing points (black open squares) with the statistical temperature TS(E) (black curve). (b) Generalized probability distributions functions Pα(E) of corresponding replicas α = 17, 18, 19,…, 26 of System 6 in Table I.

FIG. 1.

(a) Effective temperatures Tα(E) (a set of parallel lines with negative slope) form unique crossing points (black open squares) with the statistical temperature TS(E) (black curve). (b) Generalized probability distributions functions Pα(E) of corresponding replicas α = 17, 18, 19,…, 26 of System 6 in Table I.

Close modal
FIG. 3.

(a) Energy temperature curve in the canonical ensemble (red line) and molar heat capacity Cv(T) (blue line) of System 6. (b) Probability distribution function PT(E) and free energy FT(E) at the melting temperature Tm = 288.7 K.

FIG. 3.

(a) Energy temperature curve in the canonical ensemble (red line) and molar heat capacity Cv(T) (blue line) of System 6. (b) Probability distribution function PT(E) and free energy FT(E) at the melting temperature Tm = 288.7 K.

Close modal
FIG. 2.

(a) Lateral radial distribution function gxy(r) and (b) structure factor transformed from gxy(r) of replica 18 (blue line) and replica 25 (red line) of the same systems as in Fig. 1. (c) The transverse density profile of water along confinement direction (z direction) for replica 18 (blue) and replica 25 (red).

FIG. 2.

(a) Lateral radial distribution function gxy(r) and (b) structure factor transformed from gxy(r) of replica 18 (blue line) and replica 25 (red line) of the same systems as in Fig. 1. (c) The transverse density profile of water along confinement direction (z direction) for replica 18 (blue) and replica 25 (red).

Close modal

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

$E_{18}^*$
E18* = −40.45 kJ/mol and
$E_{25}^*$
E25*
= −38.20 kJ/mol, respectively, at the same
$T_{18}^*=T_{25}^*=288.7$
T18*=T25*=288.7
K. The RDFs of replicas 18 and 25 in Fig. 2(a) show marked differences in terms of the magnitude of the peaks, and the number of peaks. The RDF of replica 18 corresponding to the solid phase has several pronounced peaks up to 14  Å, while only three peaks are visible in replica 25, indicating an absence of long range order characteristic of the liquid phase. The structure factor of replica 18 displays a prepeak 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 Cv(T), can be calculated as in Eqs. (7) and (8). In contrast to the statistical temperature, TS(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 Cv around the transition temperature Tm = 288.7 K. The free energy density in Fig. 3(b) at Tm, F(E, Tm) = ETmS(E), exhibits two local minima at E1 = −38.1 and E3 = −40.5  kJ mol−1, and one local maximum at E2 = −39.5  kJ mol−1, resulting in the bimodal structure in PT(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.

FIG. 4.

(a) Energy versus temperature curve, E(T), and (b) molar heat capacity, Cv(T), of System 1–6, the parameters of which are given in Table I. The change in internal energy (latent heat) across the transition region from liquid to solid is monotonically decreasing with increasing density.

FIG. 4.

(a) Energy versus temperature curve, E(T), and (b) molar heat capacity, Cv(T), of System 1–6, the parameters of which are given in Table I. The change in internal energy (latent heat) across the transition region from liquid to solid is monotonically decreasing with increasing density.

Close modal

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 gxy(r) at high density is almost continuous across replicas, while at low density gxy(r) is well segregated as clearly seen in first and second peaks in Fig. 5(b). The continuous variation of gxy(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 gxy(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 NPxyT 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 NPxyT ensemble simulation as well.5 

FIG. 5.

Lateral radial distribution function gxy(r) of all replicas in System 2 (a) and System 6 (b). A substantial difference between systems is apparent in the long-range order, which shows a continuous transition at high density and a more abrupt transition at low density.

FIG. 5.

Lateral radial distribution function gxy(r) of all replicas in System 2 (a) and System 6 (b). A substantial difference between systems is apparent in the long-range order, which shows a continuous transition at high density and a more abrupt transition at low density.

Close modal
FIG. 6.

The transverse density profiles (TDPs) of water along the confinement direction (z) of all replicas in System 2 (a) and System 6 (b). A substantial difference between systems is apparent in the density at the slit center, which varies continuously at high density and more abruptly at low density. There is a substantially higher density of water at the slit center for the liquid phase at low density, compared with the liquid phase in the higher density system.

FIG. 6.

The transverse density profiles (TDPs) of water along the confinement direction (z) of all replicas in System 2 (a) and System 6 (b). A substantial difference between systems is apparent in the density at the slit center, which varies continuously at high density and more abruptly at low density. There is a substantially higher density of water at the slit center for the liquid phase at low density, compared with the liquid phase in the higher density system.

Close modal

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, TS(E), in Fig. 7, in which most probable energies,

$E_\alpha ^*$
Eα*⁠, corresponding to the crossing points between Tα(E) and TS(E), were plotted with
$T_\alpha ^*=T_\alpha (E_\alpha ^*;\lambda _\alpha )$
Tα*=Tα(Eα*;λα)
across replicas for 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, TS(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.

FIG. 7.

The energy-temperature curves formed by most probable energy sets

$[E_{\alpha }^*,T_{\alpha }^*]$
[Eα*,Tα*] determined by gREM simulations for systems with densities ρ1 = 1.1643  g cm−3 (a) and ρ2 = 1.0812  g cm−3 (b). The lines and symbols in red, green, and blue show the results for systems with 256, 576, and 800 molecules, respectively.

FIG. 7.

The energy-temperature curves formed by most probable energy sets

$[E_{\alpha }^*,T_{\alpha }^*]$
[Eα*,Tα*] determined by gREM simulations for systems with densities ρ1 = 1.1643  g cm−3 (a) and ρ2 = 1.0812  g cm−3 (b). The lines and symbols in red, green, and blue show the results for systems with 256, 576, and 800 molecules, respectively.

Close modal

To characterize the global minima structures for System 1–10, basin-hopping global optimization51,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.

FIG. 8.

Examples of minimum energy structures derived from basin-hopping global optimization for the smaller (N = 256) water molecule system. (a)–(f) The minimized structures of System 1–6, respectively.

FIG. 8.

Examples of minimum energy structures derived from basin-hopping global optimization for the smaller (N = 256) water molecule system. (a)–(f) The minimized structures of System 1–6, respectively.

Close modal
FIG. 9.

Examples of minimum energy structures derived from basin-hopping global optimization for the intermediate (N = 576) and large (N = 800) water molecule systems. (a)–(d) The minimized structures of System 7–10, respectively.

FIG. 9.

Examples of minimum energy structures derived from basin-hopping global optimization for the intermediate (N = 576) and large (N = 800) water molecule systems. (a)–(d) The minimized structures of System 7–10, respectively.

Close modal

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.

Table II.

Percentage distributions of the n-membered ring in System 1–10, derived from 100 minimized structures. Results for global minimum structure only are shown in parentheses.

Ring sizen = 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 sizen = 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 

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, TS(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 TS(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 NPxyT ensemble, may provide additional insight into the possible validity of this intriguing conjecture.

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

1.
C.
Pérez
,
M. T.
Muckle
,
D. P.
Zaleski
,
N. A.
Seifert
,
B.
Temelso
,
G. C.
Shields
,
Z.
Kisiel
, and
B. H.
Pate
,
Science
336
,
897
(
2012
).
2.
O.
Mishima
and
H. E.
Stanley
,
Nature (London)
396
,
329
(
1998
).
3.
E. B.
Moore
and
V.
Molinero
,
Nature (London)
479
,
506
(
2011
).
4.
J. R.
Errington
,
P. G.
Debenedetti
, and
S.
Torquato
,
Phys. Rev. Lett.
89
,
215503
(
2002
).
5.
S.
Han
,
M. Y.
Choi
,
P.
Kumar
, and
H. E.
Stanley
,
Nat. Phys.
6
,
685
(
2010
).
6.
A. A.
Bakulin
,
D.
Cringus
,
P. A.
Pieniazek
,
J. L.
Skinner
,
T. L. C.
Jansen
, and
M. S.
Pshenichnikov
,
J. Phys. Chem. B
117
,
15545
(
2013
).
7.
J.
Bai
and
X. C.
Zeng
,
Proc. Natl. Acad. Sci. U.S.A
109
,
21240
(
2012
).
8.
N.
Kastelowitz
,
J. C.
Johnston
, and
V.
Molinero
,
J. Chem. Phys.
132
,
124511
(
2010
).
9.
J. C.
Johnston
,
N.
Kastelowitz
, and
V.
Molinero
,
J. Chem. Phys.
133
,
154516
(
2010
).
10.
N.
Giovambattista
,
P. J.
Rossky
, and
P. G.
Debenedetti
,
Phys. Rev. Lett.
102
,
050603
(
2009
).
11.
P.
Kumar
,
S. V.
Buldyrev
,
F. W.
Starr
,
N.
Giovambattista
, and
H. E.
Stanley
,
Phys. Rev. E
72
,
051503
(
2005
).
12.
T.
Kaneko
,
J.
Bai
,
K.
Yasuoka
,
A.
Mitsutake
, and
X. C.
Zeng
,
J. Chem. Theory Comput.
9
,
3299
(
2013
).
13.
K.
Koga
,
X. C.
Zeng
and
H.
Tanaka
,
Phys. Rev. Lett.
79
,
5262
(
1997
).
14.
A. L.
Ferguson
,
N.
Giovambattista
,
P. J.
Rossky
,
A. Z.
Panagiotopoulos
and
P. G.
Debenedetti
,
J. Chem. Phys.
137
,
144501
(
2012
).
15.
D. J.
Wales
and
R. S.
Berry
,
Phys. Rev. Lett.
73
,
2875
(
1994
).
16.
F.
Calvo
,
D. J.
Wales
,
J. P. K.
Doye
,
R. S.
Berry
,
P.
Labastie
, and
M.
Schmidt
,
Europhys. Lett.
82
,
43003
(
2008
).
17.
J. P. K.
Doye
and
D. J.
Wales
,
J. Chem. Phys.
102
,
9673
(
1995
).
18.
D. J.
Wales
and
J. P. K.
Doye
,
J. Chem. Phys.
103
,
3061
(
1995
).
19.
R. M.
Lynden-Bell
and
D. J.
Wales
,
J. Chem. Phys.
101
,
1460
(
1994
).
20.
J.
Kim
,
T.
Keyes
, and
J. E.
Straub
,
J. Chem. Phys.
132
,
224107
(
2010
).
21.
C. J.
Geyer
and
A.
Thompson
,
J. Am. Stat. Assoc.
90
,
909
(
1995
).
22.
Y.
Sugita
and
Y.
Okamoto
,
Chem. Phys. Lett.
314
,
141
(
1999
).
23.
R.
Zhou
and
B. J.
Berne
,
Proc. Natl. Acad. Sci. U.S.A
99
,
12777
(
2002
).
24.
A. E.
Garcia
and
J. N.
Onuchic
,
Proc. Natl. Acad. Sci. U.S.A
100
,
13898
(
2003
).
25.
D.
Paschek
,
S.
Gnanakarnan
, and
A. E.
Garcia
,
Proc. Natl. Acad. Sci. U.S.A
102
,
6765
(
2005
).
26.
R.
Yamamoto
and
W.
Kob
,
Phys. Rev. E
61
,
5473
(
2000
).
27.
E.
Flenner
and
G.
Szamel
,
Phys. Rev. E
73
,
061505
(
2006
).
28.
P. A.
Frantsuzov
and
V. A.
Mandelshtam
,
Phys. Rev. E
72
,
037102
(
2005
).
29.
P.
Poulain
,
F.
Calvo
,
R.
Antoine
,
M.
Broyer
, and
P.
Dugourd
,
Phys. Rev. E
73
,
056704
(
2006
).
30.
K.
Hukushima
and
K.
Nemoto
,
J. Phys. Soc. Jpn.
65
,
1604
(
1996
).
31.
H.
Kamberaj
and
A.
van der Vaart
,
J. Chem. Phys.
127
,
234102
(
2007
).
32.
U. H. E.
Hansmann
,
Chem. Phys. Lett.
281
,
140
(
1997
).
33.
A.
Mitsutake
,
Y.
Sugita
, and
Y.
Okamoto
,
Biopolymers
60
,
96
(
2001
).
34.
Q.
Lu
,
J.
Kim
, and
J.
Straub
,
J. Phys. Chem. B
116
,
8654
(
2012
).
35.
Q.
Lu
,
J.
Kim
, and
J. E.
Straub
,
J. Chem. Phys.
138
,
104119
(
2013
).
36.
W. J.
Cho
,
J.
Kim
,
J.
Lee
,
T.
Keyes
,
J. E.
Straub
, and
K. S.
Kim
,
Phys. Rev. Lett.
112
,
157802
(
2014
).
37.
J.
Kim
,
T.
Keyes
, and
J. E.
Straub
,
J. Chem. Phys.
135
,
061103
(
2011
).
38.
L. G.
Rizzi
and
N. A.
Alves
,
J. Chem. Phys.
135
,
141101
(
2011
).
39.
M.
Church
,
C.
Ferry
, and
A. E.
van Giessen
,
J. Chem. Phys.
136
,
245102
(
2012
).
40.
A. M.
Ferrenberg
and
R. H.
Swenden
,
Phys. Rev. Lett.
63
,
1195
(
1989
).
41.
E. B.
Moore
and
V.
Molinero
,
J. Chem. Phys.
130
,
244505
(
2009
).
42.
L.
Le
and
V.
Molinero
,
J. Phys. Chem. A
115
,
5900
(
2011
).
43.
A.
Reinhardt
and
J. P. K.
Doye
,
J. Chem. Phys.
136
,
054501
(
2012
).
44.
Y.
Lu
,
X.
Zhang
, and
M.
Chen
,
J. Phys. Chem. B
117
,
10241
(
2013
).
45.
J. C.
Johnston
and
V.
Molinero
,
J. Am. Chem. Soc.
134
,
6650
(
2012
).
46.
R. C.
DeMille
,
T. E.
Cheatham
, and
V.
Molinero
,
J. Phys. Chem. B
115
,
132
(
2011
).
47.
A. H.
Nguyen
and
V.
Molinero
,
J. Phys. Chem. B
117
,
6330
(
2013
).
48.
F. H.
Stillinger
and
T. A.
Weber
,
Phys. Rev. B
31
,
5262
(
1985
).
49.
J. J.
Salacuse
,
A. R.
Denton
, and
P. A.
Egelstaff
,
Phys. Rev. E
53
,
2382
(
1996
).
50.
J. N.
Herrera
,
P. T.
Cummings
, and
H.
Ruiz-Estrada
,
Mol. Phys.
96
,
835
(
1999
).
51.
Z.
Li
and
H. A.
Scheraga
,
Proc. Natl. Acad. Sci. U.S.A.
84
,
6611
(
1987
).
52.
D. J.
Wales
and
J. P. K.
Doye
,
J. Phys. Chem. A
101
,
5111
(
1997
).
53.
D. J.
Wales
, GMIN: A program for basin-hopping global optimization, basin-sampling, and parallel tempering.
55.
D.
Liu
and
J.
Nocedal
,
Math. Program.
45
,
503
(
1989
).
56.
D. J.
Wales
and
H. A.
Scheraga
,
Science
285
,
1368
(
1999
).