Isothermal-isobaric Molecular Dynamics Simulations Research Articles

Molecular dynamics study on relaxation of supercooled liquid water at different cooling rates

The relaxation in supercooled liquids is a phenomena by which the thermodynamic system starts transforming from one state to other after losing its high end thermodynamic properties (i.e. density, potential energy, enthalpy, etc.) with a significant change. Since the relaxation phenomena is quite challenging to observe, specifically near the transition region of liquid water. Hence, the thermodynamic analysis for understanding the relaxation phenomena of liquid water in supercooled region from the high density liquid (HDL) state to low density amorphous (LDA) state are important from both fundamental and as well as applied perspectives. Therefore, we focused our study on performing thermodynamic analysis on the relaxation of super-cooled liquid water by using a monatomic water (mW) potential model with an isothermal-isobaric (NPT) molecular dynamics (MD) simulation technique at zero pressure under different cooling rates. In this work, we examine the relaxation in terms of root mean square (RMS) fluctuation of two and three-body potential energy of the maximum size of four coordinated cluster under different cooling rates. The predicted data over this study shows the independency of the cooling rate, where the relaxation of the liquid state of mW water results in a sharp irreversible drop in potential energy and large fluctuation in two and three body potential energy of the maximum size of four coordinated cluster. Further, we found the relaxation of liquid state of mW water associated with the sharp and continuous increase in the fraction of four coordinated particles as a function of temperature for different cooling rates. The predicted results of this study are consistent with the literature, which helps to understand the relaxation phenomena of mW water and other such thermodynamic systems more significantly.

Active sparse Bayesian committee machine potential for isothermal-isobaric molecular dynamics simulations.

Recent advancements in machine learning potentials (MLPs) have significantly impacted the fields of chemistry, physics, and biology by enabling large-scale first-principles simulations. Among different machine learning approaches, kernel-based MLPs distinguish themselves through their ability to handle small datasets, quantify uncertainties, and minimize over-fitting. Nevertheless, their extensive computational requirements present considerable challenges. To alleviate these, sparsification methods have been developed, aiming to reduce computational scaling without compromising accuracy. In the context of isothermal and isobaric ML molecular dynamics (MD) simulations, achieving precise pressure estimation is crucial for reproducing reliable system behavior under constant pressure. Despite progress, sparse kernel MLPs struggle with precise pressure prediction. Here, we introduce a virial kernel function that significantly enhances the pressure estimation accuracy of MLPs. Additionally, we propose the active sparse Bayesian committee machine (BCM) potential, an on-the-fly MLP architecture that aggregates local sparse Gaussian process regression (SGPR) MLPs. The sparse BCM potential overcomes the steep computational scaling with the kernel size, and a predefined restriction on the size of kernel allows for fast and efficient on-the-fly training. Our advancements facilitate accurate and computationally efficient machine learning-enhanced MD (MLMD) simulations across diverse systems, including ice-liquid coexisting phases, Li10Ge(PS6)2 lithium solid electrolyte, and high-pressure liquid boron nitride.

On the apparent molar volume of methanol in water-methanol mixtures. Composition and temperature effects from molecular dynamics study

We explore the effects of composition and temperature on the apparent molar volumes of species of water-methanol mixtures. Isothermal-isobaric molecular dynamics simulations are used with this purpose. Several combinations of models for water and for methanol are explored. Validity of predictions concerned with a puzzling minimum of apparent molar volume of methanol in water-rich solutions is tested against experimental results.

Second critical point in two realistic models of water.

The hypothesis that water has a second critical point at deeply supercooled conditions was formulated to provide a thermodynamically consistent interpretation of numerous experimental observations. A large body of work has been devoted to verifying or falsifying this hypothesis, but no unambiguous experimental proof has yet been found. Here, we use histogram reweighting and large-system scattering calculations to investigate computationally two molecular models of water, TIP4P/2005 and TIP4P/Ice, widely regarded to be among the most accurate classical force fields for this substance. We show that both models have a metastable liquid-liquid critical point at deeply supercooled conditions and that this critical point is consistent with the three-dimensional Ising universality class.

Liquid-liquid critical point of water

Chemical Physics Despite a broad range of experimental observations that indirectly point out the possible existence of a liquid-liquid phase transition in deeply supercooled water, no unambiguous experiment has shown this yet. The challenges of performing experiments under such conditions are associated with the inevitable rapid crystallization of the metastable liquid state, and therefore computational simulations become a crucial alternative. Using two accurate classical models of water and long (more than 40 microseconds) isobaricisothermal molecular dynamics simulations, Debenedetti et al. provide strong computational evidence of the presence of a second metastable liquid-liquid critical point in the deeply supercooled region that has critical behavior consistent with the three-dimensional Ising universality class. Science this issue p. [289][1] [1]: /lookup/doi/10.1126/science.abb9796

The cocrystal mechanism of HMX and LLM-105 by theoretical simulations

Isothermal-isobaric molecular dynamics (NPT-MD) simulations and density functional theory (DFT) calculations were performed to explore the cocrystal formation mechanism of energetic materials. Based on the morphologically dominant crystal habits of HMX (octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine) and LLM-105 (2,6-diamino-3,5-dinitropyrazine-l-oxide), six interface structures of HMX with LLM-105 have been established to obtain the interaction energy, radial distribution function (RDF), electrostatic potential and mechanical properties. It revealed that the interaction energies in different models are different, owing to the various molecular orientation and stackingmanner of the contacted crystal face of HMX and LLM-105. The (1 1 –1)/(1 1 0) interface is mainly governed by the O···H hydrogen bond and strong van der Waals forces and this model possesses the largest interaction energy. Electrostatic potentials indicate (1 1 –1)/(1 1 0) interface blending the strongest, which reveals that the cocrystal process is stronglyrelevant to the interface behavior of (1 1 –1)/(1 1 0) among the 6 interfaces. Moreover, the mechanical parameters show that the (1 1 –1)/(1 1 0) interface structure (VI) exhibits excellent ductility and fracture strength. Additionally, this work presents a novel insight, which from the crystal face of HMX and LLM-105 as well as the interaction behavior to research the cocrystal mechanism.

Quantitative Prediction of the Structure and Viscosity of Aqueous Micellar Solutions of Ionic Surfactants: A Combined Approach Based on Coarse-Grained MARTINI Simulations Followed by Reverse-Mapped All-Atom Molecular Dynamics Simulations.

We address the problem of the quantitative prediction of micelle formation in dilute aqueous solutions of ionic surfactants using sodium dodecyl sulfate (SDS) as a model system through a computational approach that involves three steps: (a) execution of coarse-grained simulations based on the MARTINI force field (with slightly modified parameters to afford the formation of large micelles); (b) reverse mapping of the final self-assembled coarse-grained configuration into an all-atom configuration; and (c) final relaxation of this all-atom configuration through short-time (on the order of a few tens of nanoseconds), detailed isothermal-isobaric molecular dynamics simulations using the CHARMM36 force field. For a given concentration of the solution in SDS molecules, the modified MARTINI-based coarse-grained simulations lead to the formation of large micelles characterized by mean aggregation numbers above the experimentally observed ones. However, by reintroducing the detailed chemical structure through a strategy that solves a well-defined geometric problem and re-equilibrating, these large micellar aggregates quickly dissolve to smaller ones and equilibrate to sizes that perfectly match the average micelle size measured experimentally at the given surfactant concentration. From the all-atom molecular dynamics simulations, we also deduce the surfactant diffusivity DSDS and the zero-shear rate viscosity, η0, of the solution, which are observed to compare very favorably with the few experimental values that we were able to find in the literature.

Functional uncertainty quantification for isobaric molecular dynamics simulations and defect formation energies

Functional uncertainty quantification (FunUQ) was recently proposed to quantify uncertainties in models and simulations that originate from input functions, as opposed to parameters. This paper extends FunUQ to quantify uncertainties originating from interatomic potentials in isothermal-isobaric molecular dynamics (MD) simulations and to the calculation of defect formation energies. We derive and verify a computationally inexpensive expression to compute functional derivatives in MD based on perturbation theory. We show that this functional derivative of the quantities of interest (average internal energy, volume, and defect energies in our case) with respect to the interatomic potential can be used to predict those quantities for a different interatomic potential, without re-running the simulation. The codes and scripts to perform FunUQ in MD are freely available for download. In addition, to facilitate reproducibility and to enable use of best practices for the approach, we created Jupyter notebooks to perform FunUQ analysis on MD simulations and made them available for online simulation in nanoHUB. The tool uses cloud computing resources and users can view, edit, and run end-to-end workflows from a standard web-browser without the need to download or install any software.

A polymerization scenario of the liquid–liquid transition of GeI4

SnI4 and GeI4 have been shown to exhibit similar polyamorphic nature. We examine the microscopic nature of the liquid–liquid transition in GeI4 by conducting an isothermal–isobaric molecular dynamics simulation for the system composed of rigid tetrahedral molecules. The model allows us to semiquantitatively discuss the structural properties of liquid GeI4 below 1 GPa. We define a physical bond between the nearest intermolecular iodine sites satisfying the conditions of forming the metallic I2 bond. We then focus on the formation of molecular clusters in dynamic networks of the bonds. The clusters are mainly formed by molecules whose nearest pairs are in edge-to-edge, face-to-edge, and vertex-to-edge orientations. The clusters grow as pressure increases, and the onset of percolation is observed below 1 GPa. The finite-size scaling analysis for the percolation probability identifies that the threshold pressure is GPa for the present model, which is near the extension of the boundary between the two liquid phases. We thus speculate that the liquid–liquid transition of GeI4 is attained by polymerization, i.e. percolation of molecular networks. The same percolation scenario with a slight modification such as introducing a bootstrap mechanism is expected to apply to the transition in liquid SnI4.

Molecular dynamics simulations of the properties of water-methanol mixtures. Effects of force fields

Isothermal-isobaric molecular dynamics simulations are used to examine the microscopic structure and some properties of water-methanol liquid mixture. The TIP4P/2005 and SPC/E water models are combined with the united atom TraPPE and the all-atom force field model for methanol. Our principal focus is to evaluate the quality of predictions of different combinations of model force fields concerning the composition dependence of basic properties of this system. Specifically, we explored the composition effects on density, excess molar volume and excess entropy, as well as on the surface tension and static dielectric constant. In addition, the structural properties are described in terms of the coordination numbers and the average number of hydrogen bonds between molecules of constituent species. Finally, the composition dependence of self-diffusion coefficients of the species is evaluated. All theoretical predictions are tested with respect to experimental data.

Successful test of the classical nucleation theory by molecular dynamic simulations of BaS

We used a recently developed two-body interatomic potential for barium sulfide (BaS), and performed molecular dynamics (MD) simulations with 36,000 particles to determine the kinetics of spontaneous, homogeneous nucleation and growth of BaS crystals in the supercooled liquid. Isothermal-isobaric MD simulations were accomplished at three temperatures. The calculated pair correlation function, along with several snapshots, allowed us to quantify the nucleation times, their crystal growth rates, and the time evolution of overall crystallization. Nucleation was spontaneously achieved in the supercooled liquid state; therefore, we computed the average birth times of the critical nuclei (also known as average onset time of the first nucleus) from 15 samples at each temperature, from which we computed the steady-state nucleation rates, MD Jss(T). Then, we independently obtained by MD the diffusion coefficients, D(T), the melting point, Tm, and the enthalpy of melting, ΔHm. Thus, the MD Jss could be compared with the predictions of the Classical Nucleation Theory (CNT) for homogeneous nucleation using only one fitting parameter, the nucleus/liquid interfacial free energy, σ. The calculated critical nucleus is made of only 2 to 3-unit cells, which is consistent with the critical size observed in the simulations, approximately 10–15 atoms. Using a constant (fitted value of) σ, and the D(T) and thermodynamic parameters from the simulations, we found that the MD pre-exponential factor has the same order of magnitude as the theoretical value predicted by the CNT. To the best of our knowledge, this agreement of the predictions of CNT and MD simulations has seldom been reported. Therefore, our results corroborate the validity of the CNT for simple supercooled liquids.

Quantumness and state boundaries hidden in supercritical helium-4: A path integral centroid molecular dynamics study.

Isothermal-isobaric path integral centroid molecular dynamics simulations were conducted for fluid 4He at more than 600 state points in the proximity of the critical point to reveal the detailed states and relevant quantumness underlying the supercritical state. Through intensive analyses of the thermodynamic, dynamic, and quantum properties, we revealed the hidden state boundaries that separate the liquid-like and gas-like states in the supercritical region of this fluid. The Widom line, defined as the locus of the maxima of isobaric heat capacity C P , is also the quantum boundary at which there are changes in the isobaric temperature-dependence of the quantum wavelength, λ quantum, i.e., maximum amplitude of the Feynman imaginary-time paths (necklaces) of individual atoms. The Frenkel line, the famous dynamic state boundary, was observed to start from nearly the same point, 0.73-0.76 T c, on the P-T plane as observed for classical fluids. Several state boundaries based on the new criteria were found to emanate from the critical point or its vicinity on the P-T plane and are discussed in comparison with these boundaries. The quantumness of this fluid was expressed as (a) non-classical significant depression of C P observed in the liquid-like state; (b) the depression of the slopes dP/dT of the Widom line and the liquid-gas coexistence line near the critical point; and (c) the depression of the heat of pseudo-boiling across the Widom line. This is explained in terms of the decreasing kinetic energy with temperature observed in the liquid-like state below the Widom temperature T Widom, or alternatively in terms of the lattice model heat capacity, including the λ quantum.

Effects of ion concentration and solvent composition on the properties of water-methanol solutions of NaCl. NPT molecular dynamics computer simulation results

Isothermal-isobaric molecular dynamics simulations are used to examine the microscopic structure and other properties of a model solution consisting of NaCl salt dissolved in water-methanol mixture. The SPC/E water model and the united atom model for methanol are combined with the force field for ions by Dang [J. Amer. Chem. Soc., 1995, 117, 6954] to describe the entire system. Our principal focus is to study the effects of two variables, namely, the solvent composition and ion concentrations on the solution's density, on the structural properties, self-diffusion coefficients of the species and the dielectric constant. Moreover, we performed a detailed analysis of the first coordination numbers of the species. Trends of the behaviour of the average number of hydrogen bonds between solvent molecules are evaluated.

Thermodynamics of 1-ethyl-3-methylimidazolium methylsulfate co-solvent mixtures using molecular dynamics simulation

In the present work, isothermal-isobaric (NPT) and canonical (NVT) ensemble molecular dynamics (MD) simulations have been used to calculate the binary mixture properties (excess molar heats and volumes) of 1-ethyl-3-methylimidazolium methyl sulfate [EMIM][MeSO4] with six polar covalent molecules (acetic acid, acetone, chloroform, dimethyl sulfoxide, isopropyl alcohol and methanol) over the entire composition range. Excess partial molar properties as a function of mole fraction of co-solvent were calculated by using a Redlich-Kister expansion equation. The excess molar heats and volumes for [EMIM][MeSO4] with all the six co-solvents were all negative over the entire composition range. The size and structure of these molecules, as well as the nature of interaction between the mixture components seems to play a major role in the observed effects. These results have been interpreted in terms of hydrogen bonding and intermolecular interactions between ionic liquid ions and co-solvent molecules via radial distribution functions. The results showed that hydrogen bond interactions between ionic liquid and co-solvents play a major role in the volumetric and energetic properties of the mixtures. Further, radial distribution functions suggest that ionic liquid ions aggregate at higher concentrations of co-solvents.

Molecular mechanism of melting of a helical polymer crystal: Role of conformational order, packing and mobility of polymers

The molecular mechanism of melting of a superheated helical polymer crystal has been investigated using isothermal-isobaric molecular dynamics simulation that allows anisotropic deformation of the crystal lattice. A detailed microscopic analysis of the onset and progression of melting and accompanying changes in the polymer conformational order, translational, and orientation order of the solid along the melting pathway is presented. Upon gradual heating from room temperature to beyond the melting point at ambient pressure, the crystal exhibits signatures of premelting well below the solid-to-liquid melting transition at the melting point. The melting transition is manifested by abrupt changes in the crystal volume, lattice energy, polymer conformation, and dynamical properties. In the premelting stage, the crystal lattice structure and backbone orientation of the polymer chains are retained but with the onset of weakening of long-range helical order and interchain packing of polymers perpendicular to the fibre axis of the crystal. The premelting also marks the onset of conformational defects and anisotropic solid-state diffusion of polymers along the fibre axis. The present study underscores the importance of the interplay between intermolecular packing, interactions, and conformational dynamics at the atomic level in determining the macroscopic melting behavior of polymer crystals.

On the composition dependence of the microscopic structure, thermodynamic, dynamic and dielectric properties of water-dimethyl formamide model mixtures. Molecular dynamics simulation results

Isothermal-isobaric molecular dynamics simulations have been performed to examine an ample set of properties of the model water-N,N-dimethylformamide (DMF) mixture as a function of composition. The SPC-E and TIP4P-Ew water models together with two united atom models for DMF [Chalaris M., Samios J., J. Chem. Phys., 2000, 112, 8581; Cordeiro J., Int. J. Quantum Chem., 1997, 65, 709] were used. Our principal analyses concern the behaviour of structural properties in terms of radial distribution functions, and the number of hydrogen bonds between molecules of different species as well as thermodynamic properties. Namely, we explore the density, excess mixing molar volume and enthalpy, the heat capacity and excess mixing heat capacity. Finally, the self-diffusion coefficients of species and the dielectric constant of the system are discussed. In addition, surface tension of water-DMF mixtures has been calculated and analyzed.

Isobaric-isothermal molecular dynamics computer simulations of the properties of water-1,2-dimethoxyethane model mixtures

Isothermal-isobaric molecular dynamics simulations have been performed to examine a broad set of properties of the model water-1,2-dimethoxyethane (DME) mixture as a function of composition. The SPC-E and TIP4P-Ew water models and the modified TraPPE model for DME were applied. Our principal focus was to explore the trends of behaviour of the structural properties in terms of the radial distribution functions, coordination numbers and number of hydrogen bonds between molecules of different species, and of conformations of DME molecules. Thermodynamic properties, such as density, molar volume, enthalpy of mixing and heat capacity at constant pressure have been examined. Finally, the self-diffusion coefficients of species and the dielectric constant of the system were calculated and analyzed.

How fast fast-folding proteins fold in silico

In reported microcanonical molecular dynamics simulations, fast-folding proteins CLN025 and Trp-cage autonomously folded to experimentally determined native conformations. However, the folding times of these proteins derived from the simulations were more than 4–10 times longer than their experimental values. This article reports autonomous folding of CLN025 and Trp-cage in isobaric–isothermal molecular dynamics simulations with agreements within factors of 0.69–1.75 between simulated and experimental folding times at different temperatures. These results show that CLN025 and Trp-cage can now autonomously fold in silico as fast as in experiments, and suggest that the accuracy of folding simulations for fast-folding proteins begins to overlap with the accuracy of folding experiments. This opens new prospects of developing computer algorithms that can predict both ensembles of conformations and their interconversion rates for a protein from its sequence for artificial intelligence on how and when a protein acts as a receiver, switch, and relay to facilitate various subcellular-to-tissue communications. Then the genetic information that encodes proteins can be better read in the context of intricate biological functions.

Molecular dynamics study of water diffusion in an amphiphilic block copolymer with large difference in the blocks’ glass transition temperatures

Isothermal-isobaric molecular dynamics simulation was used to study the diffusion mechanism of water in polyurethane-block-poly(N-isopropyl acrylamide) (PU-block-PNIPAm) with a hydrophobic PU/hydrophilic PNIPAm mass ratio of 1.4 to 1 at 298 K and 450 K. Here, the experimental glass transition temperature (T g ) of PU is 243 K while that of PNIPAm is 383 K. Different amounts of water up to 15 wt-% were added to PU-block-PNIPAm. We were able to reproduce the specific volumes and glass transition temperatures (250 K and 390 K) of PU-block-PNIPAm. The computed self-diffusion coefficient of water increased exponentially with increasing water concentration at both temperatures (i.e., following the free volume model of Fujita). It suggested that water diffusion in PU-block-PNIPAm depends only on its fractional free volume despite the free volume inhomogeneity. It is noted that at 298 K, PU is rubbery while PNIPAm is glassy. Regardless of temperature, radial distribution functions showed that water formed clusters with sizes in the range of 0.2–0.4 nm in PU-block-PNIPAm. At low water concentrations, more clusters were found in the PU domain but at high water concentrations, more in the PNIPAm domain. It is believed that water molecules diffuse as clusters rather than as individual molecules.

Use of molecular dynamics simulations to estimate the solubility of menadione in supercritical CO2 using Chrastil's model

The binary system of menadione in explicit supercritical carbon dioxide (SC-CO2) was studied using molecular dynamics (MD) simulations, with the objective to understand the nature of interactions between menadione and SC-CO2 at different temperatures and pressures in order to complement experimental solubility measurements. A force field was developed for menadione and tested by comparing computed and experimental monoclinic crystal structures at 283 K and 0.01 MPa. Lattice parameters obtained from anisotropic isothermal-isobaric MD simulations agreed reasonably well with experimental values, with an average absolute relative deviation (AARD%) less than 7%. A previously validated force field for SC-CO2 was used, and simple mixing rules were used to describe cross interactions. Canonical ensemble MD simulations were used to estimate the association number for CO2 about menadione and the enthalpy required to form a SC-CO2 solvate complex with menadione as a function of temperature and CO2 density. Spatial distribution functions were computed to better understand the nature of the molecular-level interactions between menadione and SC-CO2 as well as between associating menadione molecules. This work is the first part of a study that uses MD simulations as the main tool to represent a binary system. The MD methodologies of this work will be applied to represent our future studies of menadione derivatives in SC-CO2.