NPT Molecular Dynamics Simulations Research Articles

Harnessing machine learning for efficient large-scale interatomic potential for sildenafil and pharmaceuticals containing H, C, N, O, and S.

In this study a cutting-edge approach to producing accurate and computationally efficient interatomic potentials using machine learning algorithms is presented. Specifically, the study focuses on the application of Allegro, a novel machine learning algorithm, running on high-performance GPUs for training potentials. The choice of training parameters plays a pivotal role in the quality of the potential functions. To enable this methodology, the "Solvated Protein Fragments" dataset, containing nearly 2.7 million Density Functional Theory (DFT) calculations for many-body intermolecular interactions involving protein fragments and water molecules, encompassing H, C, N, O, and S elements, is considered as the training dataset. The project optimizes computational efficiency by reducing the initial dataset size according to the intended application. To assess the efficacy of the approach, the sildenafil citrate, iso-sildenafil, aspirin, ibuprofen, mebendazole and urea, representing all five relevant elements, serve as the test bed. The results of the Allegro-trained potentials demonstrate outstanding performance, benefiting from the combination of an appropriate training dataset and parameter selection. This notably enhanced computational efficiency when compared to the computationally intensive DFT method aided by GPU acceleration. Validation of the produced interatomic potentials is achieved through Allegro's own evaluation mechanism, yielding exceptional accuracy. Further verification is carried out through LAMMPS molecular dynamics simulations. Structural optimization by energy minimization and NPT Molecular Dynamics simulations are performed for each potential, assessing relaxation processes and energy reduction. Additional structures, including urea, ammonia, uracil, oxalic acid, and acetic acid, are tested, highlighting the potential's versatility in describing systems containing the aforementioned elements. Visualization of the results confirms the scientific accuracy of each structure's relaxation. The findings of this study demonstrate strong scaling and the potential for applications in pharmaceutical research, allowing the exploration of larger molecular structures not previously amenable to computational analysis at this level of accuracy The success of the machine learning approach underscores its potential to revolutionize computational solid-state physics.

Evidence of caged to normal diffusion transition in benzene along supercritical isobars: Insights from molecular dynamics simulations

Understanding the behavior of supercritical fluids holds both fundamental and practical significance. From a practical standpoint, these fluids find numerous technological applications across environmental, mechanical, chemical, biological, and geothermal industries. However, from a fundamental scientific perspective, unraveling the intricate structure, thermodynamics, and dynamics of supercritical fluids remains an ongoing challenge. In this work, we address this task by employing NPT molecular dynamics simulations on a rigid benzene model under high pressures. Our observations reveal a striking dynamical transition, with a caged regime at low temperatures and a normal diffusion regime at high temperatures. Notably, the subdiffusive region aligns with the compressible region from the experimental phase diagram, while the normal diffusion occurs in the supercritical region. Furthermore, we present compelling evidence linking the subdiffusion to a caged trajectory and illustrate how uncaging and clustering play a key role in the dynamical transition. The identification of the subdiffusive and normal diffusion regimes provides valuable information for optimizing processes involving supercritical fluids in various industries. Moreover, the link between caged trajectories and dynamical transitions adds a new dimension to our understanding of molecular behavior in supercritical systems.

Effects of thermostats/barostats on physical properties of liquids by molecular dynamics simulations

Thermostats/barostats that maintain temperature/pressure as constant (on average) in molecular dynamics (MD) simulations are essential for simulatingisothermal-isochoric (NVT) or isothermal-isobaric (NPT) systems. In order to simulate accurate physical properties, thermostats/barostats need not only sample a correct NVT/NPT ensemble but also minimally disturb the particles’ (Newtonian) dynamics. In other words, the NVT/NPT simulations need to yield accurate physical properties (dynamics properties included, e.g. diffusivity and viscosity) and also their fluctuations. However, few studies have studied in-depth the effects of thermostats/barostats on a comprehensively wide range of simulated properties. In this work, commonly used thermostats, e.g. Andersen, stochastic dynamics (SD), Berendsen, V-rescale and Nosé-Hoover thermostats, and barostats, e.g. Berendsen and Parrinello-Rahman barostats, with different coupling strengths were studied to see whether they can yield accurate simulated properties and fluctuations. The theoretical values of physical properties’ fluctuations were calculated via statistical mechanics to provide a comparison benchmark for MD simulations. Particularly, the accuracy of thermostats was studied in the context of NPT or non-equilibrium molecular dynamics (NEMD) simulations, which has long been overlooked. Berendsen thermostat/barostat suppresses the fluctuations of energy/volume (due to the exponential decay of deviation of the simulated temperature/pressure from the target value) and yields inaccurate simulated properties. In addition, Andersen and SD thermostats perturb the particles’ dynamics so violently (due to the random component in the algorithm) that they fail to accurately simulate dynamics properties. Overall, Nosé-Hoover/V-rescale thermostat, and Parrinello-Rahman barostat with moderate coupling strength are recommended for common NVT/NPT production simulations. NPT or NEMD simulations would require more efficient (or stronger coupling) thermostats than NVT equilibrium molecular dynamics (EMD) to maintain the simulated temperature close to the target.

Patterns in 2D core-softened systems: From sphere to dumbbell colloids

In this paper we analyze the structural behavior of two species of core-softened colloids: spherical and dumbbell shaped molecules with distinct separations from almost complete overlap to no-overlap. Using NPT Molecular Dynamics simulations, we identified that spherical colloids can show ordered structures as a Low Density Triangular phase, Stripes, a Kagome lattice and a High Density Triangular phase. For dumbbells colloids the structures are strongly affected by the anisotropy. For dimers with low intramolecular distance we identified the Low Density Triangular phase, a cluster phase and side-to-side Stripes, while for the cases with intermediate distance we identified the Low Density Triangular, and stripes with distinct orientational patterns. When there is no overlap between the monomers in a dumbbell, a behavior similar to the spherical symmetric case is recovered. We show that the myriad of stripes patterns that arises by introducing anisotropy into the system adds extra competition beyond the existing competition given by the core-softened interaction potential. Our results shed some light in the complex behavior of colloids assembly and show how we can use the dumbbell anisotropy to control the assembly patterns.

Nucleation parameters of SPC/E and TIP4P/2005 water vapor measured in NPT molecular dynamics simulations.

Nucleation rates for droplet formation in water vapor are measured in molecular dynamics (MD) simulations of SPC/E and TIP4P/2005 water by monitoring individual nucleation events. The nucleation process is simulated in the NPT ensemble to evaluate the steady-state nucleation rate in accordance with the assumptions of classical nucleation theory (CNT). Nucleation rates measured between 300 and 425 K for the SPC/E model, and between 325 and 475 K for the TIP4P/2005 model, agree with the CNT predictions roughly within the standard deviation of the MD measurements of the nucleation rates.

Insight into the structural, thermal and ion transport properties of solid and liquid Mg3N2: a model potential and NPT molecular dynamics simulation

ABSTRACT A model potential has been developed for magnesium nitride (Mg3N2) by fitting the parameters to an experimental lattice constant and bulk modulus. The mechanical properties have been calculated at 0 K. The potential was tested by molecular dynamics simulation at constant number, volume and energy ensemble at 300 K. The structural, ion transport and thermodynamic properties were investigated at a constant number and pressure–temperature ensemble between 300 and 2500 K. The simulated XRD diffractogram at 300 K shows that the crystalline peaks with corresponding (hkl) planes are correctly reproduced. A solid–liquid transition occurs at about 1800 K with a marked variation in all calculated physical properties. The potential successfully models the interactions in Mg3N2 at the solid and liquid phases. The results presented in this paper fill the gap of theoretical study in literature and extend our knowledge further, which could be useful to pave the way for the search of functional materials.

Molecular Dynamics Simulations of 2d Size‐Polydisperse Fluid Close to the Freezing Transition

AbstractThe effect of particle size‐polydispersity on the static equilibrium properties of a model 2d‐fluid is studied considering two types of size‐distributions, namely, Gaussian (G) and uniform (U) using computer simulations. In this study, constant NPT molecular dynamics simulations is performed and the results are compared with that of the corresponding size‐monodisperse fluid. In particular, the effect of size‐disparity on the liquid‐solid transition, local and global particle ordering around the freezing/melting transition is investigated. It is found that the transition temperature (T*) is elevated by the introduction of size polydispersity in the system. The value of global hexatic order parameter for polydisperse systems is far below 50%, while it is around 80% for size‐monodisperse fluid at the respective T*.

Molecular Dynamics Simulations of Polyethylene Bilayers

By conducting molecular dynamics (MD) simulations of polyethylene (PE) melts consisting of two different films, each comprised of chains of different molar mass, in molecular contact with each other, we compare the structural and dynamical properties between these subsystems. Joining layers of the same chemical constitution but different molecular weights is explored as a route towards packaging materials that combine good mechanical and barrier properties with recyclability and therefore provide more sustainable solutions for contemporary industrial needs. Initially, we construct two independent PE thin films characterized by periodic boundary conditions in two directions, but of finite thickness in the third direction. An “amorphous builder” is used for this purpose, which constructs the chains bead by bead (united atom), using the TraPPE force field. A slab is made by joining the two films in the thickness direction and periodic boundaries are introduced in this direction. The slab is energy minimized and then subjected to NPT molecular dynamics (MD) simulation at 350 K and 1 bar, over times longer than the longest relaxation times of both films. Both structural and dynamical properties of the films are calculated, including the self-diffusion coefficient for the chain centers of mass.

Phase Behaviors of Ionic Liquids Heating from Different Crystal Polymorphs toward the Same Smectic-A Ionic Liquid Crystal by Molecular Dynamics Simulation

Five distinct crystal structures, based on experimental data or constructed manually, of ionic liquid [C14Mim][NO3] were heated in NPT molecular dynamics simulations under the same pressure such that they melted into the liquid crystal (LC) phase and then into the liquid phase. It was found that the more entropy-favored structure had a higher solid-LC transition temperature: Before the transition into the LC, all systems had to go through a metastable state with the side chains almost perpendicular to the polar layers. All those crystals finally melted into the same smectic-A LC structure irrelevant of the initial crystal structure.

Two-Dimensional WS2@Nitrogen-Doped Graphite for High-Performance Lithium Ion Batteries: Experiments and Molecular Dynamics Simulations

As promising candidates for anode materials in lithium ion batteries (LIB), two-dimensional tungsten disulfide (WS2) and WS2@(N-doped) graphite composites were synthesized, and their electrochemical properties were comprehensibly studied in conjunction with calculations. The WS2 nanosheets, WS2@graphite, and WS2@N-doped graphite (N-graphite) exhibit outstanding cycling performance with capacities of 633, 780, and 963 mA h g-1, respectively. To understand their lithium storage mechanism, first-principles calculations involving a series of ab initio NVT- NPT molecular dynamics simulations were conducted. The calculated discharge curves for amorphous phase are well matched with the experimental ones, and the capacities reach 620, 743, and 915 mA h g-1 for WS2, WS2@graphite, and WS2@N-graphite, respectively. The large capacities of the two composites can be attributed to the tendency of W and Li atoms to interact with graphite, suppressing the formation of W metal clusters. In the case of WS2@N-graphite, vigorous amorphization of the N-graphite enhances the interaction of W and Li atoms with the fragmented N-graphite in such a way that unfavorable Li-W repulsion is avoided at very early stage of lithiation. As a result, the volume expansion in WS2@graphite and WS2@N-graphite is calculated to be remarkably small (only 6 and 44%, respectively, versus 150% for WS2). Therefore WS2@(N-)graphite composites are expected to be almost free of mechanical pulverization after repeated cycles, which makes them promising and excellent candidates for high-performance LIBs.

Local volume as a robust structural measure and its connection to icosahedral content in a model binary amorphous system

Molecular dynamics simulations are performed for a model binary glass to investigate the temperature dependence of the thermally driven localized structural excitations mediating glass relaxation. This is done under fixed zero pressure by increasing the temperature from zero to values close to and beyond the glass transition regime. It is found that for a given ramp rate, fixed pressure simulations admit a similar amount of structural relaxation as fixed volume simulations (Acta Materialia 143 (2018) 205–213), but at a reduced temperature scale corresponding to a reduced glass transition temperature. Furthermore, analysis of the volume evolution under fixed pressure, and stress evolution under fixed volume, demonstrates an equivalence between NPT and NVT ensemble molecular dynamics simulations in achieving a particular structural state. The observed relaxation induces a local densification which may be quantitatively connected to the creation of icosahedral structure that minimizes both bond energy frustration and maximizes local atomic packing. Together these results establish a robust link between the structural state of an amorphous solid and its free volume content, both of which now directly relate to the degree of relaxation.

A multiscale approach to predict the mixed gas separation performance of glassy polymeric membranes for CO2 capture: the case of CO2/CH4 mixture in Matrimid®

The gas solubility in polymeric membranes affects the separation performance, particularly in the case of CO2 capture processes. Solubility and solubility-selectivity in membranes of multicomponent mixtures can deviate rather markedly from the corresponding pure gas values, due to swelling and competition phenomena, and require dedicated time-consuming measurements. Many experiments can be avoided by using a suitable thermodynamic tool, such as an Equation of State (EoS) model, to represent the gases sorption in the membrane. Such models require, for parameterization, knowledge of the polymer behavior above the glass transition Tg, which is a limit for membrane modeling, because the most attractive polymers for gas separation are rigid matrices characterized by very high Tg values, difficult to reach experimentally. In this work, we study the sorption of CO2/CH4 mixtures in a high-Tg polyimide membrane (Matrimid®) using a bottom-up approach. Pressure-volume-temperature data for Matrimid® above Tg are generated using NPT Molecular Dynamics simulations: the results are regressed to find Matrimid® parameters for the PC-SAFT Equation of State. Finally, the Non Equilibrium PC-SAFT macroscopic model (NE-PC-SAFT) is used to calculate CO2 and CH4 solubility and solubility-selectivity as a function of gas mixture pressure, composition and temperature. The approach is tested successfully over many experimental pure gas and vapor sorption data in Matrimid®. Mixed gas calculations predict a marked competition, which affects more methane than CO2 sorption, and results in a higher-than-ideal value of solubility-selectivity. Combined with the fact that experimental mixed gas permeability-selectivity is lower than the ideal value, such results indicate that the diffusivity of CH4 in Matrimid® is significantly enhanced in presence of CO2, causing a decrease of diffusivity-selectivity.

The competition of densification and structure ordering during crystallization of HCP-Mg in the framework of layering

In this paper, we performed an NPT molecular dynamics simulation of crystallization process of HCP-Mg to probe the competition between densification and structural ordering. Two opposite layering patterns, i.e. outward and inward, were designed for analysis. From the perspective of solid-like cluster (SLC) itself, structural ordering always precedes densification; but from the perspective of SLC’s precursor, structural ordering always lags behind densification; the reversion occurs at the closest two liquid layers around SLC. We call it dip-rebound phenomenon. This phenomenon is a completely new finding. It resolves, to some extent, recent debate about whether densification or structural ordering triggers crystallization.

A priori predictions of molecular density by EFP2-MD

A priori density predictions of 32 molecules were attempted by combining the effective fragment potential version 2 (EFP2) and NPT molecular dynamics (MD) simulations. Our EFP2-MD procedure accurately predicted the density maximum of water as a function of temperature, showing its promising performance in the prediction of molecular density. With the help of a uniform scale factor of 0.9099, the mean absolute deviation and root-mean-square deviation of 32 molecular density predictions as compared to experiments are 0.037 and 0.002 g/cm3, respectively. They exhibit remarkable accuracy considering the fact that EFP2 is a purely theoretical parameter.

On the composition dependence of thermodynamic, dynamic and dielectric properties of water-dimethyl sulfoxide model mixtures. NPT molecular dynamics simulation results

We have systematically examined available molecular dynamics (MD) simulation results for thermodynamic, dynamic and dielectric properties of water-dimethyl sulfoxide (DMSO) mixtures with the purpose of analyzing their discrepancies with experimental data. To get a more profound insight into the behavior of these mixtures, we have performed constant pressure-constant temperature simulations in the entire composition range at 1bar and 298.15K and compared the obtained results with experimental data and previous findings. The SPC-E water model was used in combination with the P1, P2 and OPLS united atom models for dimethyl sulfoxide. The dependencies of density, of excess mixing volume and potential energy, of the self-diffusion coefficients of species, and of the dielectric constant on the mixture composition are discussed. The evolution of the hydrogen bonding network was analyzed in terms of the average number of hydrogen bonded molecules and of the distribution of bonding states. In conclusions we discuss the capability of the united atom models for DMSO to adequately reproduce an ample set of properties of these mixtures and outline future work.

Anion Exchange Membrane for Fuel Cell: Multi-Scale Modeling Approach

It is required to form nanophase-segregated structures in anion exchange polymer electrolyte membrane to develop and enhance the ion transport. In this study, we employ multi-scale modeling method to investigate the correlation of the ion transport with the nanophase-segregation in anion exchange polymer electrolyte membrane. The atomistic force field is developed using density functional theory (DFT) method, and used for molecular dynamics (MD) simulations. The annealing procedure is applied to model the nanophase-segregation. The equilibration is achieved through NPT MD simulation, and the extent of nanophase-segregation is evaluated by the structure factor analysis. The ion transport is evaluated using mean-square displacement. The nanophase-segregated structures and the transport properties are compared to the proton exchange membrane consisting of the same polymer backbone except for the acidic functional group. The coarse-grained model is also constructed based on atomistic model. Using dissipative particle dynamics (DPD) simulation, the meso-scale structures are simulated to elucidate large-scale correlation in structure.

Melting point trends and solid phase behaviors of model salts with ion size asymmetry and distributed cation charge.

The melting point trends of model salts composed of coarse grain ions are examined using NPT molecular dynamics simulations. The model salts incorporate ion size asymmetry and distributed cation charge, which are two common features in ionic liquids. A series of single-phase and two-phase simulations are done at set temperatures with 50 K intervals for each salt, and the normal melting point is estimated within 50 K. The melting point trends are then established relative to a charge-centered, size symmetric salt with a normal melting point between 1250 K and 1300 K. We consider two sets of size asymmetric salts with size ratios up to 3:1; the melting point trends are different in each set. The lowest melting point we find is between 450 K and 500 K, which is a reduction of over 60% from the charge-centered, size symmetric case. In both sets, we find diversity in the solid phase structures. For all size ratios with small cation charge displacements, the salts crystallize with orientationally disordered cations. When the partial cation charge is far enough off-center in salts with ion size ratios near 1:1, the salts can become trapped in glassy states and have underlying crystal structures that are orientationally ordered. At ion size ratios near 3:1, the salts with large cation charge displacements show premelting transitions at temperatures as low as 300 K. After the premelting transition, these salts exist either as fast ion conductors, where the smaller anions move through a face centered cubic (fcc) cation lattice, or as plastic crystals, where ion pairs rotate on a fcc lattice.

Theoretical investigation on the restoring step of the carbonic anhydrase catalytic cycle for natural and promiscuous substrates

In the present study steered molecular dynamics simulations were applied to investigate the unbinding process of the complex of human carbonic anhydrase with the natural HCO3− and promiscuous H2NCOHN− products. This process is crucial for restoring the catalytic cycle of the enzyme. This investigation set out to give further insights on the release mechanism involved in the case of the promiscuous product believed suicide inhibitor for the hCAII against the natural final product. In particular, on the basis of the NPT molecular dynamics simulations performed on the bicarbonate, the penta-coordinated complex with the water is observed, while in the case of the ureate the same event does not take place. At this purpose the calculated potential of mean force based on the steered molecular dynamics (SMD) simulations shed light on an optimal pathway for the releasing of the products.

Effects of pH on the association between the inhibitor cystatin and the proteinase chymopapain.

Cysteine proteinases are involved in many aspects of physiological regulation. In humans, some cathepsins have shown another function in addition to their role as lysosomal proteases in intracellular protein degradation; they have been implicated in the pathogenesis of several heart and blood vessel diseases and in cancer development. In this work, we present a fluorometric and computational study of the binding of one representative plant cysteine proteinase, chymopapain, to one of the most studied inhibitors of these proteinases: chicken cystatin. The binding equilibrium constant, Kb, was determined in the pH range between 3.5 and 10.0, revealing a maximum in the affinity at pH 9.0. We constructed an atomic model for the chymopapain-cystatin dimer by docking the individual 3D protein structures; subsequently, the model was refined using a 100 ns NPT molecular dynamics simulation in explicit water. Upon scrutiny of this model, we identified 14 ionizing residues at the interface of the complex using a cutoff distance of 5.0 Å. Using the pKa values predicted with PROPKA and a modified proton-linkage model, we performed a regression analysis on our data to obtain the composite pKavalues for three isoacidic residues. We also calculated the electrostatic component of the binding energy (ΔGb,elec) at different pH values using an implicit solvent model and APBS software. The pH profile of this calculated energy compares well with the experimentally obtained binding energy, ΔGb. We propose that the residues that form an interchain ionic pair, Lys139A from chymopapain and Glu19B from cystatin, as well as Tyr61A and Tyr67A from chymopapain are the main residues responsible for the observed pH dependence in the chymopapain- cystatin affinity.

Molecular mechanics and equation of state modeling of compressible polyolefin solutions: Impact of pressure and cut-off radius of intermolecular potentials

This paper reports on the density of solutions of polyethylene (PE) in hexane using molecular dynamics (MD) simulations based upon the accurate OPLS-AA force field at high pressures and 425K. The NPT–MD simulations are carried out at polymer concentration of 20wt% in hexane in the pressure range from 100 to 3000bar. For PE solutions of 10 and 20wt% the pressure is varied from 100 to 1000bar. Additionally, the PE solution densities are calculated based on the modified Sanchez–Lacombe (MSL) equation of state (EOS) model to examine the accuracy of the MD computations. The simulated densities increase monotonically with increasing external pressure and compare quite favorably with the experimental and EOS data. It is also revealed that the MSL EOS model produces identical mixture densities regardless of the type of the b parameter. The effect of cut-off radius to density is investigated and it is shown that the solution density increases as cut-off radius increases. A minimum cut-off radius of 1.1nm is suggested for the intermolecular forces for accurate densities at pressures below 100bar. For higher pressures density and non-bonded interactions display less sensitivity to cut-off distance. Analysis of the pair distribution function versus pressure is carried out where the height of the first peak increases and the radial distribution function shifts to shorter separations reflecting structural change of the condensed phase. The molecular modeling approach employed in this research provides a good insight into the polymer–polymer, polymer–solvent, and solvent–solvent interactions. The implemented methodology using the OPLS-AA force field and constant pressure/temperature algorithms compare well with the literature data, suggesting the validity of the proposed method.