Direct Molecular Dynamics Simulations Research Articles

Full Temperature-Dependent Potential and Anharmonicity in Metallic Hydrogen: Colossal NQE and the Consequences

The temperature-dependent effective potential (TDEP) method for anharmonic phonon dispersion is generalized to the full potential case by combining with path integral formalism. This extension naturally resolves the intrinsic difficulty in the original TDEP at low temperature. The new method is applied to solid metallic hydrogen at high pressure. A colossal nuclear quantum effect (NQE) and subsequent anharmonicity are discovered, which not only leads to unexpectedly large drift of protons, but also slows down the convergence rate substantially when computing the phonon dispersions. By employing direct ab initio path integral molecular dynamics simulations as the benchmark, a possible breakdown of phonon picture in metallic hydrogen due to colossal NQE is indicated, implying novel lattice dynamical phenomena might exist. Inspired by this observation, a general theoretical formalism for quantum lattice dynamics beyond phonon is sketched, with the main features being discussed.

Molecular dynamics analysis on bending mechanical behavior of alumina nanowires at different loading rates

The molecular dynamics (MD) model of α-Al2O3 nanowires in bending is established by using LAMMPS to calculate the atomic stress and strain at different loading rates in order to study the effect of loading rate on the bending mechanical behaviors of the α-Al2O3 nanowires. Research results show that the maximum surface stress—rotation angle curves of α-Al2O3 nanowires at different loading rates are all divided into three stages of elastic deformation, plastic deformation and failure, where the elastic limit point can be determined by the curve symmetry during loading and unloading cycle. The loading rate has great influence on the plastic deformation but little on the elastic modulus of α-Al2O3 nanowires. When the loading rate is increased, the plastic deformation stage is shortened and the material is easier to fail in brittle fracture. Therefore, the elastic limit and the strength limit (determined by the direct and indirect MD simulation methods) are closer to each other. The MD simulation result of α-Al2O3 nanowires is verified to be valid by the good agreement with the improved loop test results. The direct MD method becomes an effective way to determine the elastic limit and the strength limit of nanoscale whiskers failed in brittle or ductile fracture at arbitrary loading rate.

On the Use of a Non-Constant Non-Affine or Slip Parameter in Polymer Rheology Constitutive Modeling

Since its introduction in the late 1970s, the non-affine or slip parameter, ξ, has been routinely employed by numerous constitutive models as a constant parameter. However, the evidence seems to imply that it should be a function of polymer deformation. In the present work, we phenomenologically modify a constitutive model for the rheology of unentangled polymer melts [P. S. Stephanou et al. J. Rheol. 53, 309 (2009)] to account for a non-constant slip parameter. The revised model predictions are compared against newly accumulated rheological data for a C48 polyethylene melt obtained via direct non-equilibrium molecular dynamics simulations in shear. We find that the conformation tensor data are very well predicted; however, the predictions of the material functions are noted to deviate from the NEMD data, especially at large shear rates.

Instability dynamics of nonlinear normal modes in the Fermi–Pasta–Ulam–Tsingou chains

Nonlinear normal modes are periodic orbits that survive in nonlinear chains, whose instability plays a crucial role in the dynamics of many-body Hamiltonian systems toward thermalization. Here we focus on how the stability of nonlinear modes depends on the perturbation strength and the system size to observe whether they have the same behavior in different models. To this end, as illustrating examples, the instability dynamics of the N/2 mode in both the Fermi–Pasta–Ulam–Tsingou-α and -β chains under fixed boundary conditions are studied systematically. Applying the Floquet theory, we show that for both models the stability time T as a function of the perturbation strength λ follows the same behavior; i.e., , where λ c is the instability threshold. The dependence of λ c on N is also obtained. The results of T and λ c agree well with those obtained by the direct molecular dynamics simulations. Finally, the effect of instability dynamics on the thermalization properties of a system is briefly discussed.

Misfolding-Associated Exposure of Natively Buried Residues in Mutant SOD1 Facilitates Binding to TRAF6

Amyotrophic lateral sclerosis (ALS) is a fatal neurodegenerative disease primarily impacting motor neurons. Mutations in superoxide dismutase 1 (SOD1) are the second most common cause of familial ALS. Several of these mutations lead to misfolding or toxic gain of function in the SOD1 protein. Recently, we reported that misfolded SOD1 interacts with TNF receptor-associated factor 6 (TRAF6) in the SOD1G93A rat model of ALS. Further, we showed in cultured cells that several mutant SOD1 proteins, but not wildtype SOD1 protein, interact with TRAF6 via the MATH domain. Here, we sought to uncover the structural details of this interaction through molecular dynamics (MD) simulations of a dimeric model system, coarse grained using the AWSEM force field. We used direct MD simulations to identify buried residues, and predict binding poses by clustering frames from the trajectories. Metadynamics simulations were also used to deduce preferred binding regions on the protein surfaces from the potential of the mean force in orientation space. Well-folded SOD1 was found to bind TRAF6 via co-option of its native homodimer interface. However, if loops IV and VII of SOD1 were disordered, as typically occurs in the absence of stabilizing Zn2+ ion binding, these disordered loops now participated in novel interactions with TRAF6. On TRAF6, multiple interaction hot-spots were distributed around the equatorial region of the MATH domain beta barrel. Expression of TRAF6 variants with mutations in this region in cultured cells demonstrated that TRAF6T475 facilitates interaction with different SOD1 mutants. These findings contribute to our understanding of the disease mechanism and uncover potential targets for the development of therapeutics.

Dynamic compaction of aluminum with nanopores of varied shape: MD simulations and machine-learning-based approximation of deformation behavior

We compare two machine-learning-based approaches, artificial neural network (ANN) and micromechanical model with automatic Bayesian identification of the model parameters, in application to mimicking the deformation behavior of nanoporous aluminum extracted from molecular dynamics (MD) simulations. Reference data are generated by means of MD simulation of both hydrostatic and uniaxial deformation with compression of representative volume elements of aluminum single crystal with nanopores of spherical, cubic and cylindrical shapes at the temperatures of 300, 500, 700 and 900 K. Several typical sizes of the nanopores are considered: The smallest one corresponds to the initial porosity less than 1%, while the largest one gives the initial porosity in the range of 30–50%. Plastic collapse of nanopores in all cases occurred by the mechanism of emission of partial Shockley dislocation half-loops from the pore surface. The emission of initial loop occurred earlier in the case of cylindrical pore; however, this did not lead to an explosive increase in the number of dislocations in the system at this stage. In general, flat free surfaces of pores are less subjected to the dislocation nucleation than the rounded ones. A new physically-based micromechanical model of the plastic compaction of nanoporous metal is formulated with accounting of different pore shapes and anisotropy of the compaction process. Both tested machine-learning approaches show an adequate approximation of MD data. The developed ANN and parameterized micromechanical model are applied to simulate the propagation of a shock wave in nanoporous aluminum in comparison with direct MD simulations of this process; this comparison shows an adequate description of the shock wave structure by means of both approaches incorporated into continuum mechanics modeling. Thus, the developed machine-learning-based approaches can be applied as constitutive equations of nanoporous aluminum in macroscopic simulations of the dynamic compaction and shock-wave processes in this material.

Freezing point depression of salt aqueous solutions using the Madrid-2019 model.

Salt aqueous solutions are relevant in many fields, ranging from biological systems to seawater. Thus, the availability of a force-field that is able to reproduce the thermodynamic and dynamic behavior of salt aqueous solutions would be of great interest. Unfortunately, this has been proven challenging, and most of the existing force-fields fail to reproduce much of their behavior. In particular, the diffusion of water or the salt solubility are often not well reproduced by most of the existing force-fields. Recently, the Madrid-2019 model was proposed, and it was shown that this force-field, which uses the TIP4P/2005 model for water and non-integer charges for the ions, provides a good description of a large number of properties, including the solution densities, viscosities, and the diffusion of water. In this work, we assess the performance of this force-field on the evaluation of the freezing point depression. Although the freezing point depression is a colligative property that at low salt concentrations depends solely on properties of pure water, a good model for the electrolytes is needed to accurately predict the freezing point depression at moderate and high salt concentrations. The coexistence line between ice and several salt aqueous solutions (NaCl, KCl, LiCl, MgCl2, and Li2SO4) up to the eutectic point is estimated from direct coexistence molecular dynamics simulations. Our results show that this force-field reproduces fairly well the experimentally measured freezing point depression with respect to pure water freezing for all the salts and at all the compositions considered.

Understanding Hydrogenation Chemistry at MgB2 Reactive Edges from Ab Initio Molecular Dynamics.

Solid-state hydrogen storage materials often operate via transient, multistep chemical reactions at complex interfaces that are difficult to capture. Here, we use direct ab initio molecular dynamics simulations at accelerated temperatures and hydrogen pressures to probe the hydrogenation chemistry of the candidate material MgB2 without a priori assumption of reaction pathways. Focusing on highly reactive (101̅0) edge planes where initial hydrogen attack is likely to occur, we track mechanistic steps toward the formation of hydrogen-saturated BH4- units and key chemical intermediates, involving H2 dissociation, generation of functionalities and molecular complexes containing BH2 and BH3 motifs, and B-B bond breaking. The genesis of higher-order boron clustering is also observed. Different charge states and chemical environments at the B-rich and Mg-rich edge planes are found to produce different chemical pathways and preferred speciation, with implications for overall hydrogenation kinetics. The reaction processes rely on B-H bond polarization and fluctuations between ionic and covalent character, which are critically enabled by the presence of Mg2+ cations in the nearby interphase region. Our results provide guidance for devising kinetic improvement strategies for MgB2-based hydrogen storage materials, while also providing a template for exploring chemical pathways in other solid-state energy storage reactions.

Understanding the role of anharmonic phonons in diffusion of bcc metals

Diffusion in the high-temperature bcc phase of IIIB-IVB metals such as Zr, Ti, and their alloys is observed to be orders of magnitude higher than bcc metals of group VB-VIB, including Cr, Mo, and W. The underlying reason for this higher diffusion is still poorly understood. To explain this observation, we compare the first-principles-calculated parameters of monovacancy-mediated diffusion between bcc Ti, Zr, and dilute Zr- Sn alloys and bcc Cr, Mo, and W. Our results indicate that strongly anharmonic vibrations promote both the vacancy concentration and the diffusive jump rate in bcc IVB metals and can explain their markedly faster diffusion compared to bcc VIB metals. Additionally, we provide an efficient approach to calculate diffusive jump rates according to the transition state theory (TST). The use of standard harmonic TST is impractical in bcc IIIB/IVB metals due to the existence of ill-defined harmonic phonons, and most studies use classical or ab initio molecular dynamics for direct simulation of diffusive jumps. Here, instead, we use a stochastically sampled temperature-dependent phonon analysis within the transition state theory to study diffusive jumps without the need of direct molecular dynamics simulations. We validate our first-principles diffusion coefficient predictions with available experimental measurements and explain the underlying reasons for the promotion of diffusion in bcc IVB metals/alloys compared to bcc VIB metals.

Assessing the Accuracy of Machine Learning Thermodynamic Perturbation Theory: Density Functional Theory and Beyond.

Machine learning thermodynamic perturbation theory (MLPT) is a promising approach to compute finite temperature properties when the goal is to compare several different levels of ab initio theory and/or to apply highly expensive computational methods. Indeed, starting from a production molecular dynamics trajectory, this method can estimate properties at one or more target levels of theory from only a small number of additional fixed-geometry calculations, which are used to train a machine learning model. However, as MLPT is based on thermodynamic perturbation theory (TPT), inaccuracies might arise when the starting point trajectory samples a configurational space which has a small overlap with that of the target approximations of interest. By considering case studies of molecules adsorbed in zeolites and several different density functional theory approximations, in this work we assess the accuracy of MLPT for ensemble total energies and enthalpies of adsorption. It is shown that problematic cases can be detected even without knowing reference results and that even in these situations it is possible to recover target level results within chemical accuracy by applying a machine-learning-based Monte Carlo (MLMC) resampling. Finally, on the basis of the ideas developed in this work, we assess and confirm the accuracy of recently published MLPT-based enthalpies of adsorption at the random phase approximation level, whose high computational cost would completely hinder a direct molecular dynamics simulation.

CARNOT: a Fragment-Based Direct Molecular Dynamics and Virtual-Reality Simulation Package for Reactive Systems.

Traditionally, the study of reaction mechanisms of complex reaction systems such as combustion has been performed on an individual basis by optimizations of transition structure and minimum energy path or by reaction dynamics trajectory calculations for one elementary reaction at a time. It is effective, but time-consuming, whereas important and unexpected processes could have been missed. In this article, we present a direct molecular dynamics (DMD) approach and a virtual-reality simulation program, CARNOT, in which plausible chemical reactions are simulated simultaneously at finite temperature and pressure conditions. A key concept of the present ab initio molecular dynamics method is to partition a large, chemically reactive system into molecular fragments that can be adjusted on the fly of a DMD simulation. The theory represents an extension of the explicit polarization method to reactive events, called ReX-Pol. We propose a highest-and-lowest adapted-spin approximation to define the local spins of individual fragments, rather than treating the entire system by a delocalized wave function. Consequently, the present ab initio DMD can be applied to reactive systems consisting of an arbitrarily varying number of closed and open-shell fragments such as free radicals, zwitterions, and separate ions found in combustion and other reactions. A graph-data structure algorithm was incorporated in CARNOT for the analysis of reaction networks, suitable for reaction mechanism reduction. Employing the PW91 density functional theory and the 6-31+G(d) basis set, the capabilities of the CARNOT program were illustrated by a combustion reaction, consisting of 28 650 atoms, and by reaction network analysis that revealed a range of mechanistic and dynamical events. The method may be useful for applications to other types of complex reactions.

Sweep-tracing algorithm: in silico slip crystallography and tension-compression asymmetry in BCC metals

Direct Molecular Dynamics (MD) simulations are being increasingly employed to model dislocation-mediated crystal plasticity with atomic resolution. Thanks to the dislocation extraction algorithm (DXA), dislocation lines can be now accurately detected and positioned in space and their Burgers vector unambiguously identified in silico, while the simulation is being performed. However, DXA extracts static snapshots of dislocation configurations that by themselves present no information on dislocation motion. Referred to as a sweep-tracing algorithm (STA), here we introduce a practical computational method to observe dislocation motion and to accurately quantify its important characteristics such as preferential slip planes (slip crystallography). STA reconnects pairs of successive snapshots extracted by DXA and computes elementary slip facets thus precisely tracing the motion of dislocation segments from one snapshot to the next. As a testbed for our new method, we apply STA to the analysis of dislocation motion in large-scale MD simulations of single crystal plasticity in BCC metals. We observe that, when the crystal is subjected to uniaxial deformation along its [001] axis, dislocation slip predominantly occurs on the {112} maximum resolved shear stress plane under tension, while in compression slip is non-crystallographic (pencil) resulting in asymmetric mechanical response. The marked contrast in the observed slip crystallography is attributed to the twinning/anti-twinning asymmetry of shears in the {112} planes relatively favoring dislocation motion in the twinning sense while hindering dislocations from moving in the anti-twinning directions.

Barrierless HNO3 formation from the hydrolysis reaction of NO2 with Cl atom in the atmosphere

Nitric acid (HNO3) is one of the main pollutants in the atmosphere and has an important effect on the air pollution. However, the formation pathways of HNO3 are still poorly understood in the marine boundary layer (MBL). Here, we show evidence of the formation of HNO3 from chlorine (Cl) radical and nitrogen dioxide (NO2), with 1–3 water molecules, as well as on the air-water interface, based on the direct Born-Oppenheimer molecular dynamics (BOMD) simulation. The relative free energy change along the reaction coordinates for the (Cl)(NO2)(H2O)1 (monohydrate) system was calculated using thermodynamic integration, revealing that the monohydrate system exhibited no free energy barrier and was 26.5 kcal mol−1 exergonic. This suggests that the formation of HNO3 from the monohydrate system is thermodynamically favorable. For all the systems studied, the formation of HNO3 was directly observed using BOMD simulations, and only one water molecule participates in the reaction. However, loop structures, which were expected to be ubiquitous on the air-water interface, were not observed. This study identifies the importance of the reactions of Cl atom, NO2, and water molecules in the formation of HNO3 in the MBL region.

Electro-osmotic diode based on colloidal nano-valves between double membranes

The rectification of electro-osmotic flows is important in micro/nano fluidics applications such as micro-pumps and energy conversion devices. Here, we propose a simple electro-osmotic diode in which colloidal particles are contained between two parallel membranes with different pore densities. While the flow in the forward direction just pushes the colloidal particles toward the high-pore-density membrane, the backward flow is blocked by the particles near the low-pore-density membrane, which clog the pores. Nonequilibrium molecular dynamics simulations show a strong nonlinear dependence on the electric field for both the electric current and electro-osmotic flow, indicating diode characteristics. A mathematical model to reproduce the electro-osmotic diode behavior is constructed, introducing an effective pore diameter as a model for pores clogged by the colloidal particles. Good agreement is obtained between the proposed model with estimated parameter values and the results of direct molecular dynamics simulations. The proposed electro-osmotic diode has potential application in downsized microfluidic pumps, e.g., the pump induced under AC electric fields.

Coarse‐Grained Molecular Simulation of Polymers Supported by the Use of the SAFT‐γ$\gamma$ Mie Equation of State

AbstractA framework to self‐consistently combine a classical equation of state (EoS) and a molecular force field to model polymers and polymer mixtures is presented. The statistical associating fluid theory (SAFT‐ Mie) model is used to correlate the thermophysical properties of oligomers and generate robust and transferrable coarse‐grained (CG) molecular parameters which can be used both in particle based molecular simulations and in equations of state (EoS) calculations. Examples are provided for polyethylene, polypropylene, polyisobutylene atactic polystyrene, 1,4‐cis‐butadiene, polyisoprene, their blends and mixtures with low molecular weight solvents. Different types of liquid–liquid phase behavior are quantitatively captured both by the EoS and by direct molecular dynamics simulations. The use of CG models following this top‐down approach extends the time and length scales accessible to molecular simulation while retaining quantitative accuracy as compared to experimental results.

Nuclear Quantum Effect and Its Temperature Dependence in Liquid Water from Random Phase Approximation via Artificial Neural Network.

We report structural and dynamical properties of liquid water described by the random phase approximation (RPA) correlation together with the exact exchange energy (EXX) within density functional theory. By utilizing thermostated ring polymer molecular dynamics, we examine the nuclear quantum effects and their temperature dependence. We circumvent the computational limitation of performing direct first-principles molecular dynamics simulation at this high level of electronic structure theory by adapting an artificial neural network model. We show that the EXX+RPA level of theory accurately describes liquid water in terms of both dynamical and structural properties.

Vapor-liquid equilibrium of water with the MB-pol many-body potential.

Among the many existing molecular models of water, the MB-pol many-body potential has emerged as a remarkably accurate model, capable of reproducing thermodynamic, structural, and dynamic properties across water's solid, liquid, and vapor phases. In this work, we assessed the performance of MB-pol with respect to an important set of properties related to vapor-liquid coexistence and interfacial behavior. Through direct coexistence classical molecular dynamics simulations at temperatures of 400K < T < 600K, we calculated properties such as equilibrium coexistence densities, vapor-liquid interfacial tension, vapor pressure, and enthalpy of vaporization and compared the MB-pol results to experimental data. We also compared rigid vs fully flexible variants of the MB-pol model and evaluated system size effects for the properties studied. We found that the MB-pol model predictions are in good agreement with experimental data, even for temperatures approaching the vapor-liquid critical point; this agreement was largely insensitive to system sizes or the rigid vs flexible treatment of the intramolecular degrees of freedom. These results attest to the chemical accuracy of MB-pol and its high degree of transferability, thus enabling MB-pol's application across a large swath of water's phase diagram.

Grain-Size-Governed Shear Failure Mechanism of Polycrystalline Methane Hydrates

The shear failure mechanism of polycrystalline gas hydrates is critical for understanding marine geohazards related to gas hydrates under a changing climate and for safe gas recovery from gas hydrate reservoirs. Since current experimental techniques cannot resolve the mechanism on a spatial and temporal nanoscale, molecular simulations can assist with proposing and substantiating nanoscale failure mechanisms. Here, we report the shear failure of polycrystalline methane hydrates using direct molecular dynamics simulations. Based on these simulations, we suggest two modes of shear behavior, depending on the grain sizes, d, in the polycrystal: grain-size-strengthening behavior with a d1/3 grain size dependence for small grain sizes and grain-size-weakening behavior for large grain sizes. Through the crossover from strengthening to weakening behavior, the failure mode changes from shear failure with a failure plane parallel to the applied shear to tensile failure with a failure plane lying at an angle with the applied shear, spanning a network of grain boundaries. The existence of such a change in mechanism suggests that the Hall–Petch breakdown in methane hydrates is due to a change from grain boundary sliding to tensile opening being the most important failure mechanism when the grain size increases.

Stochastic Modelling of 13C NMR Spin Relaxation Experiments in Oligosaccharides.

A framework for the stochastic description of relaxation processes in flexible macromolecules including dissipative effects has been recently introduced, starting from an atomistic view, describing the joint relaxation of internal coordinates and global degrees of freedom, and depending on parameters recoverable from classic force fields (energetics) and medium modelling at the continuum level (friction tensors). The new approach provides a rational context for the interpretation of magnetic resonance relaxation experiments. In its simplest formulation, the semi-flexible Brownian (SFB) model has been until now shown to reproduce correctly correlation functions and spectral densities related to orientational properties obtained by direct molecular dynamics simulations of peptides. Here, for the first time, we applied directly the SFB approach to the practical evaluation of high-quality 13C nuclear magnetic resonance relaxation parameters, and , and the heteronuclear NOE of several oligosaccharides, which were previously interpreted on the basis of refined ad hoc modelling. The calculated NMR relaxation parameters were in agreement with the experimental data, showing that this general approach can be applied to diverse classes of molecular systems, with the minimal usage of adjustable parameters.

Barrierless HONO and HOS(O)2-NO2 Formation via NH3-Promoted Oxidation of SO2 by NO2.

In the troposphere, the knowledge about nitrous acid (HONO) sources is incomplete. The missing source of sulfate and fine particles cannot be explained during haze events. Air quality models cannot predict high levels of secondary fine-particle pollution. Despite extensive studies, one challenging issue in atmospheric chemistry is identifying the source of HONO. Here, we present direct ab initio molecular dynamics simulation evidence and typical air pollution events of the formation of gaseous HONO, nitrogen dioxide/hydrogen sulfite (HOS(O)2-NO2 or NO2-HSO3) from nitrogen dioxide (NO2), sulfur dioxide (SO2), water (H2O), and ammonia (NH3) molecules in a proportion of 2:1:3:3. The reactions show a new mechanism for the formation of HONO and NO2-HSO3 in the troposphere, especially when the concentration of NO2, SO2, H2O, and NH3 is high (e.g., 2:1:3:3 or higher) in the air. Contrary to the proportion NO2, SO2, H2O, and NH3 equaling to 1:1:3:1 and 1:1:3:2, the proportion (2:1:3:3) enables barrierless reactions and weak interactions between molecules via the formation of HONO, NO2-HSO3, and NH3/H2O. In addition, field observations are carried out, and the measured data are summarized. Correlation analysis supported the conversion of NO2 to HONO during observational studies. The weak interactions promote proton transfer, resulting in the generation of HONO, NO2-HSO3, and NH3/H2O pairs.