Direct Molecular Dynamics Simulations Research Articles

Intramolecular Vibrational Energy Redistribution in the Reaction H3+ + CO → H2 + HCO.

An ab initio direct molecular dynamics (MD) calculation at the RS2/aug-cc-pVQZ level, followed by vibration mapping, has been applied to the H3+ + CO → H2 + HCO+ reaction to elucidate the intramolecular vibrational energy redistribution (IVR) processes during the reaction. Direct MD calculations were carried out for 20 K (a typical temperature for interstellar dark clouds) and 330 K (a typical translational temperature for ions in a glow discharge). Under the Cs symmetry constraint, the approach of H3+ turned out to be the H-C stretching mode of the [H···CO]+ part, which invoked the C-O stretching and then the H-C-O bending modes. Under no symmetry constraint, a strong bending mode was first invoked, and the intensities of the subsequent H-C and C-O stretching modes were kept relatively small. The detailed analyses of the IVR during the reaction, in terms of vibration mixing, gave a clue to understanding experimentally observed anomalies in the bending modes, such as population inversion at some bending states. In the MD simulation at 20 K, less than two-thirds of the reaction energy was converted to the vibrational energy of the resultant HCO+ part and one-third to the translational and rotational energies of the leaving H2 molecule. These direct MD simulations, when combined with the experimental spectroscopy data, shed light on a clear understanding of the reaction mechanism, including the IVR during the reaction.

Simulation of lithium hydroxide decomposition using deep potential molecular dynamics.

Chemical reactions and vapor-liquid equilibria for molten lithium hydroxide (LiOH) were studied using molecular dynamics simulations and a deep potential (DP) model. The neural network for the model was trained on quantum density functional theory data for a range of conditions. The DP model allows simulations over timescales of hundreds of ns, which provide equilibrium compositions for the systems of interest. Single-phase NPT simulations of the liquid show the decomposition of LiOH into lithium oxide (Li2O) and dissolved water (H2O). These DP results were validated by direct abinitio molecular dynamics simulations that confirmed the accuracy of the model with respect to reaction kinetics and equilibrium properties of the melt. The reactive vapor-liquid behavior of this system was subsequently studied using direct coexistence interfacial DP simulations. Partial pressures of H2O in the vapor are found to be in close agreement with available experimental measurements. By fitting temperature-dependent expressions for the reaction equilibrium and Henry's law constants, the equilibrium composition for any given initial composition and temperature can be quantitatively modeled. For high initial concentrations of Li2O or H2O, mixtures of LiOH + Li2O/H2O are found to undergo phase separation. The present study illustrates how DP-based molecular dynamics simulations can be used for quantitative modeling of multiphase reactive behavior with the accuracy of the underlying abinitio quantum chemical methods.

Contact-Map-Driven Exploration of Heterogeneous Protein-Folding Paths.

We have recently shown how physically realizable protein-folding pathways can be generated using directed walks in the space of inter-residue contact-maps; combined with a back-transformation to move from protein contact-maps to Cartesian coordinates, we have demonstrated how this approach can generate protein-folding trajectory ensembles without recourse to molecular dynamics. In this article, we demonstrate that this framework can be used to study a challenging protein-folding problem that is known to exhibit two different folding paths which were previously identified through molecular dynamics simulation at several different temperatures. From the viewpoint of protein-folding mechanism prediction, this particular problem is extremely challenging to address, specifically involving folding to an identical nontrivial compact native structure along distinct pathways defined by heterogeneous secondary structural elements. Here, we show how our previously reported contact-map-based protein-folding strategy can be significantly enhanced to enable accurate and robust prediction of heterogeneous folding paths by (i) introducing a novel topologically informed metric for comparing two protein contact maps, (ii) reformulating our graph-represented folding path generation, and (iii) introducing a new and more reliable structural back-mapping algorithm. These changes improve the reliability of generating structurally sound folding intermediates and dramatically decrease the number of physically irrelevant folding intermediates generated by our previous simulation strategy. Most importantly, we demonstrate how our enhanced folding algorithm can successfully identify the alternative folding mechanisms of a multifolding-pathway protein, in line with direct molecular dynamics simulations.

Nonlinear Vlasov and Fokker-Planck Dynamics in Confined Systems with Drag: A Numerical Study

Abstract We explore numerically the behavior of a one-dimensional many-body system consisting of particles that interact through short range repulsive forces, and are also under the effects of drag forces, and of an external confining potential. The statistical dynamics of systems of this kind exhibits interesting links with the thermostatistical formalism based on the Sq non-additive entropies. In the regime of overdamped motion, these systems admit an effective description in terms of a non-linear Fokker-Planck equation. When the overdamped condition is relaxed, and inertial effects are explicitly taken into account, the system can be described by a Vlasov-like effective mean field dynamics. The Vlasov-like description of this type of systems has been recently investigated in the literature from an analytical point of view. In the present contribution we explore the behaviour of these system numerically, through direct molecular dynamics simulations. We consider examples of systems with four different short range repulsive forces, with a dependence on distance given by exponential, Heaviside, Lorentzian, and Bessel functions. The results of our numerical simulations are fully consistent with the predictions derived from the Vlasov-like mean field description. In particular, we verify that, in the asymptotic limit of large times, the system evolves towards a state exhibiting a spatial distribution of particles that coincides with the stationary solution of an appropriate nonlinear Fokker-Planck equation. This limit spatial distribution has the form of a q-Gaussian.

Learning dislocation dynamics mobility laws from large-scale MD simulations

By dispensing with all the atoms and only focusing on dislocation lines, the computational method of Discrete Dislocation Dynamics (DDD) gains greatly over Molecular Dynamics (MD) in simulation efficiency of metal plasticity. But whereas in MD dislocations follow natural dynamics of atomic motion, DDD must rely on a dislocation mobility function to prescribe how a dislocation line should respond to the driving force exerted on it. However, reflecting our still incomplete understanding of ways in which dislocations move, mobility functions presently employed in DDD simulations entail simplifications and approximations of limited or, worse still, unknown accuracy and applicability. Here we introduce a data-driven approach in which the dislocation mobility function is modeled as a graph neural network (GNN) trained on large-scale MD simulations of crystal plasticity. We apply our proposed approach to predicting plastic strength of body-centered-cubic (BCC) metal tungsten and show that, once implemented in a DDD model, our GNN dislocation mobility function accurately reproduces the challenging tension/compression asymmetry of plastic flow observed both in ground-truth MD simulations and in experiment. Furthermore, subsequently validated by MD simulations, the same function accurately predicts plastic response of tungsten under conditions not previously seen in training. By demonstrating its ability to learn relevant physics of dislocation motion, our DDD+ML approach opens a promising avenue to bringing fidelity of DDD models closer in line with direct MD simulations at a much reduced computational cost.

First principles validation of energy barriers in Ni75Al25

Precipitates in nickel-based superalloys form during heat treatment on a time scale inaccessible to direct molecular dynamics simulation, but can be studied using kinetic Monte Carlo (KMC) modelling. This requires reliable values for the barrier energies separating distinct configurations over the trajectory of the system. In this study, we validate vacancy migration barriers found with the Activation-Relaxation Technique nouveau (ARTn) method in partially ordered Ni75Al25 with a monovacancy using published potentials for the atomic interactions against first-principles methods. In a first step, we confirm that the ARTn barrier energies agree with those determined with the nudged elastic band (NEB) method. As the number of atoms used in those calculations is too great for direct ab initio calculations, we cut the cell size to 255 atoms, thus controlling finite size effects. We then use the plane-wave density functional theory code CASTEP and its inbuilt NEB method in the smaller cells. This provides us with a continuous validation chain from first principles to KMC simulations with interatomic potentials (IPs). We evaluate the barrier energies of five further IPs with NEB, demonstrating that none yields values with sufficient reliability for KMC simulations, with some of them failing completely. This is a first step towards quantifying the errors incurred in KMC simulations of precipitate formation and evolution.

Atomistic Details of Methyl Linoleate Pyrolysis: Direct Molecular Dynamics Simulation of Converting Biodiesel to Petroleum Products

Dependence on petroleum and petrochemical products is unsustainable; it is both a finite resource and an environmental hazard. Biodiesel has many attractive qualities, including a sustainable feedstock; however, it has its complications. The pyrolysis (a process already in common use in the petroleum industry) of biodiesel has demonstrated the formation of smaller hydrocarbons comprising many petrochemical products but experiments suffer from difficulty quantifying the myriad reaction pathways followed and products formed. A computational simulation of pyrolysis using “ab initio molecular dynamics” offers atomic-level detail of the reaction pathways and products formed. Herein, the most prevalent fatty-acid ester (methyl linoleate) from the most prevalent feedstock for biodiesel in the United States (soybean oil) is studied. Temperature acceleration within the atom-centered density matrix propagation formalism (Car–Parrinello) utilizing the D3-M06-2X/6-31+G(d,p) model chemistry is used to compose an ensemble of trajectories. The results are grounded in comparison to experimental studies through agreement in the following: (1) the extent of reactivity (40% in the experimental and 36.1% in this work), (2) the homology of hydrocarbon products formed (wt % of C6–C10 products), and (3) the CO/CO2 product ratio. Deoxygenation pathways are critically analyzed (as the presence of oxygen in biodiesel represents a disadvantage in its current use). Within this ensemble, deoxygenation was found to proceed through two subclasses: (1) spontaneous deoxygenation, following one of four possible pathways; or (2) induced deoxygenation, following one of three possible pathways.

High throughput screening of thermal interface materials by machine learning

Till now, it remains a challenge for effective prediction and screening of novel materials with high thermal conductivity, as well as further optimization of the interface thermal resistance. Normally, people have to spend long time on tedious calculations when predicting and screening these materials. In this paper, I combined machine learning with molecular dynamics simulations to investigate the thermal conductive properties of materials with the aim of significantly reducing computational consumption. I first applied molecular dynamics simulations to obtain the relevant properties of materials, then generated models for predicting physical properties by machine learning, and finally made predictions of thermophysical properties of materials. The use of machine learning significantly reduces the prediction time compared to direct molecular dynamics simulations. Especially when the XGBoost and the neural network models are employed, the prediction efficiency is significantly improved. This work guides a new way for the future screening of high-performance thermal interface materials.

Revisiting Textbook Azide-Clock Reactions: A "Propeller-Crawling" Mechanism Explains Differences in Rates.

An ongoing challenge to chemists is the analysis of pathways and kinetics for chemical reactions in solution, including transient structures between the reactants and products that are difficult to resolve using laboratory experiments. Here, we enabled direct molecular dynamics simulations of a textbook series of chemical reactions on the hundreds of ns to μs time scale using the weighted ensemble (WE) path sampling strategy with hybrid quantum mechanical/molecular mechanical (QM/MM) models. We focused on azide-clock reactions involving addition of an azide anion to each of three long-lived trityl cations in an acetonitrile-water solvent mixture. Results reveal a two-step mechanism: (1) diffusional collision of reactants to form an ion-pair intermediate; (2) "activation" or rearrangement of the intermediate to the product. Our simulations yield not only reaction rates that are within error of experiment but also rates for individual steps, indicating the activation step as rate-limiting for all three cations. Further, the trend in reaction rates is due to dynamical effects, i.e., differing extents of the azide anion "crawling" along the cation's phenyl-ring "propellers" during the activation step. Our study demonstrates the power of analyzing pathways and kinetics to gain insights on reaction mechanisms, underscoring the value of including WE and other related path sampling strategies in the modern toolbox for chemists.

Mechanism of thermally-activated prismatic slip in Mg

Prismatic slip of the screw <a> dislocation in magnesium at temperatures ≳ 150 K is understood to be governed by double-cross-slip of the stable basal screw through the unstable prism screw and back to the basal screw, with the activation energy controlled by the formation energy of two basal-basal kinks. However, atomistic studies of the double-kink process predict activation energies roughly twice those derived experimentally. Here, a new mechanism of prism glide is proposed, analyzed theoretically, and demonstrated qualitatively and quantitatively via direct molecular dynamics (MD) simulations. The new mechanism is intrinsically 3d, and involves the nucleation of a single kink at the junction where a 3d prismatic dislocation loop transitions from the basal screw segment to non-screw prismatic character. The relevant kink energies are calculated using recently-developed Neural-Network Potentials (NNPs) for Mg that show good agreement versus DFT for basal and prism <a> dislocations, enabling a parameter-free analytic model for the activation barrier. Direct MD simulations show both operation of the precise proposed mechanism and a stress-dependent activation barrier that agrees reasonably with the analytic model. Predictions of dislocation velocity compare very well with in-situ TEM data, and macroscopic strength versus temperature can also be understood. Overall, the new intrinsically 3d mechanism for dislocation glide due to single kink nucleation rather than double-kink nucleation explains key features of the prismatic slip in Mg and may have broader applicability in other metals where kink nucleation processes control thermally-activated flow.

Scalable simulation of coupled adsorption and transport of methane in confined complex porous media with density preconditioning

The growing significance of shales and tight formations in the transition to less carbon-intensive and clean energy drives the research endeavor to understand the physics of gas flow within these systems. However, shales are composed of massively heterogeneous physical and chemical features. Most nano-sized pores connect to millimeter-scale fractures, leading to multiscale transport. These nano-scale pore throats demonstrate non-classical flow behavior, such as non-negligible slip velocities and adsorbed gas layers at the boundary. As a result, classical computational fluid dynamics models do not capture the physics. In this work, we develop a coupling scheme for the multiple-relaxation-time (MRT) lattice Boltzmann (LB) method that integrates the Peng-Robinson equation of state into a pseudo-potential interaction model to capture the physics of methane flow in irregular networks of channels that represent nano-scale porous media. We use atomistic simulations to calibrate and validate our model in slit nano-channels. We propose a preconditioning scheme to initialize the coupled transport and adsorption simulation of methane in complex porous media. The results of this implementation of LB agree with Direct Simulation Monte Carlo (DSMC) and Molecular Dynamics (MD) simulations. We then scale up the LB implementation through vectorization and indirect addressing. We parallelize it using Message Passing Interface (MPI) and OpenMP frameworks to simulate transport and adsorption in complex media with a million lattices. We analyze the differences between coupled and transport-only simulations in two case studies and show that considering phase behavior, i.e., adsorption, can significantly change the flow behavior. This work constitutes an important step towards bridging the gap between molecular flow and system-scale behavior of complex disordered porous media.

Collisional-energy-cascade model for nonthermal velocity distributions of neutral atoms in plasmas.

Nonthermal velocity distributions with much greater tails than a Maxwellian have been observed for radical atoms in plasmas for a long time. Historically, such velocity distributions have been modeled by a two-temperature Maxwell distribution. In this paper, I propose a model based on collisional energy cascade, which has been studied in the field of granular materials. In the collisional energy cascade, a particle ensemble undergoes energy input at the high-energy region, entropy production by elastic collisions among particles, and energy dissipation. For radical atoms, energy input may be caused by the Franck-Condon energy of molecular dissociation or charge-exchange collision with hot ions, and the input energy is eventually dissipated by collisions with the walls. I show that the steady-state velocity distribution in the collisional energy cascade is approximated by the generalized Mittag-Leffler distribution, which is a one-parameter extension of the Maxwell distribution. This parameter indicates the degree of the nonthermality and is related to the relative importance of energy dissipation over entropy production. This model is compared with a direct molecular dynamics simulation for simplified gaseous systems with energy input, as well as some experimentally observed velocity distributions of light radicals in plasmas.

Calculation of Self, Corrected, and Transport Diffusivities of Isopropyl Alcohol in UiO-66.

The UiO-6x family of metal-organic frameworks has been extensively studied for applications in chemical warfare agent (CWA) capture and destruction. An understanding of intrinsic transport phenomena, such as diffusion, is key to understanding experimental results and designing effective materials for CWA capture. However, the relatively large size of CWAs and their simulants makes diffusion in the small-pored pristine UiO-66 very slow and hence impractical to study directly with direct molecular simulations because of the time scales required. We used isopropanol (IPA) as a surrogate for CWAs to investigate the fundamental diffusion mechanisms of a polar molecule within pristine UiO-66. IPA can form hydrogen bonds with the μ3-OH groups bound to the metal oxide clusters in UiO-66, similar to some CWAs, and can be studied by direct molecular dynamics simulations. We report self, corrected, and transport diffusivities of IPA in pristine UiO-66 as a function of loading. Our calculations highlight the importance of the accurate modeling of the hydrogen bonding interactions on diffusivities, with about an order of magnitude decrease in diffusion coefficients when the hydrogen bonding between IPA and the μ3-OH groups is included. We found that a fraction of the IPA molecules have very low mobility during the course of a simulation, while a small fraction are highly mobile, exhibiting mean square displacements far greater than the ensemble average.

Molecular dynamics analysis of particle number fluctuations in the mixed phase of a first-order phase transition

Molecular dynamics simulations are performed for a finite non-relativistic system of particles with Lennard-Jones potential. We study the effect of liquid-gas mixed phase on particle number fluctuations in coordinate subspace. A metastable region of the mixed phase, the so-called nucleation region, is analyzed in terms of a non-interacting cluster model. Large fluctuations due to spinodal decomposition are observed. They arise due to the interplay between the size of the acceptance region and that of the liquid phase. These effects are studied with a simple geometric model. The model results for the scaled variance of particle number distribution are compared with those obtained from the direct molecular dynamic simulations.

Assessing phonon coherence using spectroscopy

As a fundamental physical quantity of thermal phonons, temporal coherence participates in a broad range of thermal and phononic processes, while a clear methodology for the measurement of phonon coherence is still lacking. In this paper, we derive a theoretical model for the experimental exploration of phonon coherence based on spectroscopy, which is then validated by comparing Brillouin light scattering data and direct molecular dynamic simulations of confined modes in nanostructures. The proposed model highlights that confined modes exhibit a pronounced wavelike behavior characterized by a higher ratio of coherence time to lifetime. The dependence of phonon coherence on system size is also demonstrated from the spectroscopy data. The proposed theory allows for reassessing the data of conventional spectroscopy to yield coherence times, which are essential for the understanding and the estimation of phonon characteristics and heat transport in solids in general.

Unpolarized laser method for infrared spectrum calculation of amide I C[dbnd]O bonds in proteins using molecular dynamics simulation

The investigation of the strong infrared (IR)-active amide I modes of peptides and proteins has received considerable attention because a wealth of detailed information on hydrogen bonding, dipole-dipole interactions, and the conformations of the peptide backbone can be derived from the amide I bands. The interpretation of experimental spectra typically requires substantial theoretical support, such as direct ab-initio molecular dynamics simulation or mixed quantum-classical description. However, considering the difficulties associated with these theoretical methods and their applications are limited in small peptides, it is highly desirable to develop a simple yet efficient approach for simulating the amide I modes of any large proteins in solution. In this work, we proposed a comprehensive computational method that extends the well-established molecular dynamics (MD) simulation method to include an unpolarized IR laser for exciting the CO bonds of proteins. We showed the amide I frequency corresponding to the frequency of the laser pulse which resonated with the CO bond vibration. At this frequency, the protein energy and the CO bond length fluctuation were maximized. Overall, the amide I bands of various single proteins and amyloids agreed well with experimental data. The method has been implemented into the AMBER simulation package, making it widely available to the scientific community. Additionally, the application of the method to simulate the transient amide I bands of amyloid fibrils during the IR laser-induced disassembly process was discussed in details.

Generating Protein Folding Trajectories Using Contact-Map-Driven Directed Walks.

Recent advances in machine learning methods have had a significant impact on protein structure prediction, but accurate generation and characterization of protein-folding pathways remains intractable. Here, we demonstrate how protein folding trajectories can be generated using a directed walk strategy operating in the space defined by the residue-level contact-map. This double-ended strategy views protein folding as a series of discrete transitions between connected minima on the potential energy surface. Subsequent reaction-path analysis for each transition enables thermodynamic and kinetic characterization of each protein-folding path. We validate the protein-folding paths generated by our discretized-walk strategy against direct molecular dynamics simulations for a series of model coarse-grained proteins constructed from hydrophobic and polar residues. This comparison demonstrates that ranking discretized paths based on the intermediate energy barriers provides a convenient route to identifying physically sensible folding ensembles. Importantly, by using directed walks in the protein contact-map space, we circumvent several of the traditional challenges associated with protein-folding studies, namely, long time scales required and the choice of a specific order parameter to drive the folding process. As such, our approach offers a useful new route for studying the protein-folding problem.

A universal code for phase separation of prion-like low complexity domains: Insights from quantitative computer simulations.

Biomolecular phase separation is a critical mechanism that contributes to intracellular spatiotemporal organization via the formation of biomolecular condensates. In addition to their widespread implications in aggregation-related disorders, prion-like low complexity domains (PLCDs) have been shown to exhibit demixing behavior and their aberrant condensation has been linked to disease. Here, we investigate the relationship between amino acid sequence and phase behaviour of the PLCDs exploring an unprecedented set of 140 different protein sequences. To this end, we perform direct coexistence molecular dynamics simulations of our residue-resolution coarse-grained protein model, which we have shown accurately predicts quantitative phase diagrams. For variants where experimental phase behaviors have been characterized, we show that our predictions are in excellent agreement with the experimental results. Analysis of the full simulation set of 140 variants reveals that the identity and patterning of sticker residues (e.g., Y, F) are crucial to define the phase behaviour of PLCDs. Notably, aromatic and charged residues have context-dependent effects on the critical temperature of PLCD condensation. Regarding spacer residues (e.g., S, G, N, T), PLCDs condensates are remarkably stable in front of mutations of glycine and serine but are more sensitive to mutations of asparagine and glutamine. Remarkably, these results extend to the entire family of PLCDs tested and, therefore, reveal universal codes for PLCDs phase behavior. Furthermore, the findings here help shed light on the functions of PLCDs and aid in the design of modifications to counter their aberrant phase separation.

Enabling molecular dynamics simulations of helium bubble formation in tritium-containing austenitic stainless steels: An Fe-Ni-Cr-H-He potential

Formation of helium bubbles impacts mechanical properties of materials used in nuclear applications. An Fe-Ni-Cr-H-He quinary potential capable of direct molecular dynamics simulations of nucleation and growth of helium bubbles from randomly-born helium interstitial atoms has been developed. This is accomplished by incorporating helium into an existing Fe-Ni-Cr-H quaternary potential while addressing three challenging paradoxes characterizing helium in austenitic stainless steels: (a) interstitial He atoms form tightly bound dimers and larger clusters in the lattice but He atoms are only bound by weak van de Waals forces in the pure gas phase, (b) He atoms diffuse readily in host metals yet significantly distort the lattice causing large volume expansions, and (c) He atoms prefer tetrahedral interstitial sites as opposed to the larger octahedral sites despite large repulsive interactions with metal atoms within the lattice. We demonstrate that our potential reproduces density functional theory results on important properties relevant to helium bubble nucleation and growth. In addition to molecular statics validation of static properties, molecular dynamics simulation tests establish that our potential leads to the nucleation of helium bubbles from an initial random distribution of He interstitial atoms while at the same time capturing the equation of state in the pure He phase.

Extending and validating bubble nucleation rate predictions in a Lennard-Jones fluid with enhanced sampling methods and transition state theory.

We calculate bubble nucleation rates in a Lennard-Jones fluid through explicit molecular dynamics simulations. Our approach-based on a recent free energy method (dubbed reweighted Jarzynski sampling), transition state theory, and a simple recrossing correction-allows us to probe a fairly wide range of rates in several superheated and cavitation regimes in a consistent manner. Rate predictions from this approach bridge disparate independent literature studies on the same model system. As such, we find that rate predictions based on classical nucleation theory, direct brute force molecular dynamics simulations, and seeding are consistent with our approach and one another. Published rates derived from forward flux sampling simulations are, however, found to be outliers. This study serves two purposes: First, we validate the reliability of common modeling techniques and extrapolation approaches on a paradigmatic problem in materials science and chemical physics. Second, we further test our highly generic recipe for rate calculations, and establish its applicability to nucleation processes.