Investigating the solid-liquid phase transition of water nanofilms using the generalized replica exchange method

Effect of confinement on water properties in super-hydrophilic pores using MD simulations with the mW model

CONTEXT

We explore the influence of strongly hydrophilic confinement on various properties of water, such as density, enthalpy, potential energy, radial distribution function, entropy, specific heat capacity, structural dynamics, and transition temperatures (freezing and melting temperatures), using monatomic water (mW) model. The properties of water are found to be dependent on confinement and the wall-fluid surface interaction. Hysteresis loops are observed for density, enthalpy, potential energy, and entropy around the transition temperatures, while the size of hysteresis loops varies with confinement and surface interaction. In smaller pore sizes (H ≤ 20), the solid phase displays a higher density compared to the liquid phase, which is unconventional behavior compared to bulk water systems due to the pronounced hydrophilic properties of the confinement surface. Specific heat capacity exhibits more oscillations in the confined system compared to bulk water, stemming from uneven enthalpy differences across equal temperature intervals. During phase transformation in both heating and quenching processes, there is an abrupt change observed in specific heat capacity. Confinement exerts a notable impact on entropy in the solid phase, but its influence is negligible in the liquid phase. At lower pore sizes (H < 25 Å), there is more fluctuation in freezing temperature for all wall-fluid interactions, which diminishes beyond pore sizes of H > 25 Å. Similarly, more oscillatory behavior is observed in melting temperatures at lower pore sizes (H < 40 Å), which diminishes at higher pore sizes (H > 40 Å). During the quenching process, a sudden jump in the in-plane orientational and tetrahedral order parameters indicates the formation of an ordered phase, specifically a diamond crystalline structure. The percentages of different crystalline structures (cubic diamond, hexagonal diamond, and 2D hexagonal) vary with both the confinement size and the wall-fluid interaction strength.

METHODS

Cooling and heating simulations are conducted with the mW water model using LAMMPS for different nanoscale confinement separation sizes ranging from 8.5 to 70 Å within the temperature range of 100-350 K. The water is modeled using two-body and three-body interaction potential (Stillinger-Weber potential) and the confinement is introduced using LJ 9-3 water-wall interaction potential. Entropy is calculated using RDF data obtained from the simulation experiments for each temperature point with increments or decrements of 2.5 K. The transition temperatures are estimated using the specific heat capacity analysis.

Simulating the nematic-isotropic phase transition of liquid crystal model via generalized replica-exchange method

The nematic-isotropic (NI) phase transition of 4-cyano-4'-pentylbiphenyl was simulated using the generalized replica-exchange method (gREM) based on molecular dynamics simulations. The effective temperature is introduced in the gREM, allowing for the enhanced sampling of configurations in the unstable region, which is intrinsic to the first-order phase transition. The sampling performance was analyzed with different system sizes and compared with that of the temperature replica-exchange method (tREM). It was observed that gREM is capable of sampling configurations at sufficient replica-exchange acceptance ratios even around the NI transition temperature. A bimodal distribution of the order parameter at the transition region was found, which is in agreement with the mean-field theory. In contrast, tREM is ineffective around the transition temperature owing to the potential energy gap between the nematic and isotropic phases.

Alchemical Free-Energy Calculations by Multiple-Replica λ-Dynamics: The Conveyor Belt Thermodynamic Integration Scheme

A new method is proposed to calculate alchemical free-energy differences based on molecular dynamics (MD) simulations, called the conveyor belt thermodynamic integration (CBTI) scheme. As in thermodynamic integration (TI), K replicas of the system are simulated at different values of the alchemical coupling parameter λ. The number K is taken to be even, and the replicas are equally spaced on a forward-turn-backward-turn path, akin to a conveyor belt (CB) between the two physical end-states; and as in λ-dynamics (λD), the λ-values associated with the individual systems evolve in time along the simulation. However, they do so in a concerted fashion, determined by the evolution of a single dynamical variable Λ of period 2π controlling the advance of the entire CB. Thus, a change of Λ is always associated with K/2 equispaced replicas moving forward and K/2 equispaced replicas moving backward along λ. As a result, the effective free-energy profile of the replica system along Λ is periodic of period 2 πK^-1, and the magnitude of its variations decreases rapidly upon increasing K, at least as K^-1 in the limit of large K. When a sufficient number of replicas is used, these variations become small, which enables a complete and quasi-homogeneous coverage of the λ-range by the replica system, without application of any biasing potential. If desired, a memory-based biasing potential can still be added to further homogenize the sampling, the preoptimization of which is computationally inexpensive. The final free-energy profile along λ is calculated similarly to TI, by binning of the Hamiltonian λ-derivative as a function of λ considering all replicas simultaneously, followed by quadrature integration. The associated quadrature error can be kept very low owing to the continuous and quasi-homogeneous λ-sampling. The CBTI scheme can be viewed as a continuous/deterministic/dynamical analog of the Hamiltonian replica-exchange/permutation (HRE/HRP) schemes or as a correlated multiple-replica analog of the λD or λ-local elevation umbrella sampling (λ-LEUS) schemes. Compared to TI, it shares the advantage of the latter schemes in terms of enhanced orthogonal sampling, i.e. the availability of variable-λ paths to circumvent conformational barriers present at specific λ-values. Compared to HRE/HRP, it permits a deterministic and continuous sampling of the λ-range, is expected to be less sensitive to possible artifacts of the thermo- and barostating schemes, and bypasses the need to carefully preselect a λ-ladder and a swapping-attempt frequency. Compared to λ-LEUS, it eliminates (or drastically reduces) the dead time associated with the preoptimization of a biasing potential. The goal of this article is to provide the mathematical/physical formulation of the proposed CBTI scheme, along with an initial application of the method to the calculation of the hydration free energy of methanol.

Conformational analysis of replica exchange MD: Temperature-dependent Markov networks for FF amyloid peptides

Recent molecular modeling methods using Markovian descriptions of conformational states of biomolecular systems have led to powerful analysis frameworks that can accurately describe their complex dynamical behavior. In conjunction with enhanced sampling methods, such as replica exchange molecular dynamics (REMD), these frameworks allow the systematic and accurate extraction of transition probabilities between the corresponding states, in the case of Markov state models, and of statistically-optimized transition rates, in the case of the corresponding coarse master equations. However, applying automatically such methods to large molecular dynamics (MD) simulations, with explicit water molecules, remains limited both by the initial ability to identify good candidates for the underlying Markovian states and by the necessity to do so using good collective variables as reaction coordinates that allow the correct counting of inter-state transitions at various lag times. Here, we show that, in cases when representative molecular conformations can be identified for the corresponding Markovian states, and thus their corresponding collective evolution of atomic positions can be calculated along MD trajectories, one can use them to build a new type of simple collective variable, which can be particularly useful in both the correct state assignment and in the subsequent accurate counting of inter-state transition probabilities. In the case of the ubiquitously used root-mean-square deviation (RMSD) of atomic positions, we introduce the relative RMSD (RelRMSD) measure as a good reaction coordinate candidate. We apply this method to the analysis of REMD trajectories of amyloid-forming diphenylalanine (FF) peptides-a system with important nanotechnology and biomedical applications due to its self-assembling and piezoelectric properties-illustrating the use of RelRMSD in extracting its temperature-dependent intrinsic kinetics, without a priori assumptions on the functional form (e.g., Arrhenius or not) of the underlying conformational transition rates. The RelRMSD analysis enables as well a more objective assessment of the convergence of the REMD simulations. This type of collective variable may be generalized to other observables that could accurately capture conformational differences between the underlying Markov states (e.g., distance RMSD, the fraction of native contacts, etc.).

Peptide dimerization-dissociation rates from replica exchange molecular dynamics

We show how accurate rates of formation and dissociation of peptide dimers can be calculated using direct transition counting (DTC) from replica-exchange molecular dynamics (REMD) simulations. First, continuous trajectories corresponding to system replicas evolving at different temperatures are used to assign conformational states. Second, we analyze the entire REMD data to calculate the corresponding rates at each temperature directly from the number of transition counts. Finally, we compare the kinetics extracted directly, using the DTC method, with indirect estimations based on trajectory likelihood maximization using short-time propagators and on decay rates of state autocorrelation functions. For systems with relatively low-dimensional intrinsic conformational dynamics, the DTC method is simple to implement and leads to accurate temperature-dependent rates. We apply the DTC rate-extraction method to all-atom REMD simulations of dimerization of amyloid-forming NNQQ tetrapetides in explicit water. In an assessment of the REMD sampling efficiency with respect to standard MD, we find a gain of more than a factor of two at the lowest temperature.

Two Dimensional Ice from First Principles: Structures and Phase Transitions

Despite relevance to disparate areas such as cloud microphysics and tribology, major gaps in the understanding of the structures and phase transitions of low-dimensional water ice remain. Here, we report a first principles study of confined 2D ice as a function of pressure. We find that at ambient pressure hexagonal and pentagonal monolayer structures are the two lowest enthalpy phases identified. Upon mild compression, the pentagonal structure becomes the most stable and persists up to ∼2  GPa, at which point the square and rhombic phases are stable. The square phase agrees with recent experimental observations of square ice confined within graphene sheets. This work provides a fresh perspective on 2D confined ice, highlighting the sensitivity of the structures observed to both the confining pressure and the width.

Molecular simulations of heterogeneous ice nucleation. II. Peeling back the layers

Coarse grained molecular dynamics simulations are presented in which the sensitivity of the ice nucleation rate to the hydrophilicity of a graphene nanoflake is investigated. We find that an optimal interaction strength for promoting ice nucleation exists, which coincides with that found previously for a face centered cubic (111) surface. We further investigate the role that the layering of interfacial water plays in heterogeneous ice nucleation and demonstrate that the extent of layering is not a good indicator of ice nucleating ability for all surfaces. Our results suggest that to be an efficient ice nucleating agent, a surface should not bind water too strongly if it is able to accommodate high coverages of water.