Direct Molecular Dynamics Research Articles

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.

C-C Bond Formation Reaction Catalyzed by a Lithium Atom: Benzene-to-Biphenyl Coupling.

Transition-metal-catalyzed carbon-carbon (C-C) bond formation is an important reaction in pharmaceutical and organic chemistry. However, the reaction process is composed of multiple steps and is expensive owing to the presence of transition metals. This study proposes a lithium-catalyzed C-C coupling reaction of two benzene molecules (Bz) to form a biphenyl molecule, which is a transition-metal-free reaction, based on ab initio and direct ab initio molecular dynamics (AIMD) calculations. The static ab initio calculations indicate that the reaction of two Bz molecules with Li- ions (reactant state, RC) can form a stable sandwiched complex (precomplex), where the Li- ion is sandwiched by two Bz molecules. The complex formation reaction can be expressed as 2Bz + Li - → Bz(Li -)Bz, where the C-C distance between the Bz rings is 2.449 Å. This complex moves to the transition state (TS) via the structural deformation of Bz(Li-)Bz, where the C-C distance is shortened to 2.118 Å. The barrier height was calculated to be -9.9 kcal/mol (relative to RC) at the MP2/6-311++G(d,p) level. After TS, the C(sp3)-C(sp3) single bond was completely formed between the Bz rings (the C-C bond distance was 1.635 Å) (late complex). After the dissociation of H2 from the late complex, a biphenyl molecule was formed: the C(sp2)-C(sp2) bond. The calculations suggest that the C-C bond coupling of Bz occurred spontaneously from 2Bz + Li-, and biphenyl molecules were directly formed without an activation barrier. Direct AIMD calculations show that the C-C coupling reaction also takes place under electron attachment to Li(Bz)2: Li(Bz)2 + e- → [Li-(Bz)2]ver → precomplex → TS → late complex, where [Li-(Bz)2]ver is the vertical electron capture species of Li(Bz)2. Namely, the C-C coupling reaction spontaneously occurred in Li(Bz)2 owing to electron attachment. Similar C-C coupling reactions were also observed for halogen-substituted benzene molecules (Bz-X, X = F and Cl). Furthermore, this study discusses the mechanism of C-C bond formation in electron capture based on the theoretical results.

Atomistic model of an oxide film in contact with a liquid metal coolant: Defects concentrations and chemical potentials of dissolved Fe–O

The paper considers a prototype atomistic model of an oxide film with a wüstite structure in contact with a liquid metal coolant containing dissolved iron and oxygen. On the basis of direct molecular dynamics modeling with the Tersoff interatomic potential for Fe–O, we check the correspondence between the formation free energies of defects, determined from the temperature dependencies of their equilibrium concentrations, and the formation free energies determined from the calculation of the chemical potentials of dissolved iron and oxygen, the static calculations for defect energies and thermodynamic integration for defect entropies in the oxide lattice. The influence of liquid metal coolant on the chemical potentials of solvated Fe–O is highlighted. These results pave a way for bridging the gap between multiscale models of oxide film formation and ab initio first principles calculations of atomistic models.

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.

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.

Formation Mechanism of Odd- and Even-Numbered Hydrogen Cluster Cations Using the Direct Ab Initio Molecular Dynamics Approach.

The photoreaction of hydrogen clusters caused by cosmic-ray irradiations plays a crucial role in interstellar molecular clouds because the reaction of molecules takes place in these surrounding clusters. Previous gas-phase experiments introduced the reaction products after the ionization of the neutral hydrogen cluster. In the experiments, the odd-numbered hydrogen cluster cation, H+(H2)m, was the main product, while the even-numbered cluster cation, (H2)n+, was the minor product. However, the formation yield of (H2)n+ was significantly low. Despite many theoretical calculations and experiments, the formation mechanism of odd- and even-numbered cluster ions from ionized hydrogen clusters is still unclear. In this study, the direct ab initio molecular dynamics (AIMD) method was applied to the reaction of the hydrogen cluster cation (H2)n+ to understand the abovementioned formation mechanism. The trajectories of (H2)n+ following the vertical ionization of the neutral cluster were calculated. Direct AIMD calculations indicated that odd-numbered cluster cations, H3+(H2)n-2 + H (dissociation), were formed directly as the main product after the ionization of (H2)n. An even-numbered cluster, (H2)n+, was temporally formed when an H-shaped complex composed of three H2 molecules existed within a neutral cluster. However, it was rapidly dissociated to an odd-numbered cluster. Static ab initio calculations suggested that the odd-numbered cluster complex, H3+(H2)n-2-H, could be transferred to even-numbered ions via a transition state with a low energy barrier. The even-numbered cluster cation was formed stepwise from the odd-numbered cluster cation. The formation mechanism of odd- and even-numbered cluster cations is discussed based on these results.

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.

Towards fully ab initio simulation of atmospheric aerosol nucleation

Atmospheric aerosol nucleation contributes to approximately half of the worldwide cloud condensation nuclei. Despite the importance of climate, detailed nucleation mechanisms are still poorly understood. Understanding aerosol nucleation dynamics is hindered by the nonreactivity of force fields (FFs) and high computational costs due to the rare event nature of aerosol nucleation. Developing reactive FFs for nucleation systems is even more challenging than developing covalently bonded materials because of the wide size range and high dimensional characteristics of noncovalent hydrogen bonding bridging clusters. Here, we propose a general workflow that is also applicable to other systems to train an accurate reactive FF based on a deep neural network (DNN) and further bridge DNN-FF-based molecular dynamics (MD) with a cluster kinetics model based on Poisson distributions of reactive events to overcome the high computational costs of direct MD. We found that previously reported acid-base formation rates tend to be significantly underestimated, especially in polluted environments, emphasizing that acid-base nucleation observed in multiple environments should be revisited.

Mpemba-like effect protocol for granular gases of inelastic and rough hard disks

We study the conditions under which a Mpemba-like effect emerges in granular gases of inelastic and rough hard disks driven by a class of thermostats characterized by the splitting of the noise intensity into translational and rotational counterparts. Thus, granular particles are affected by a stochastic force and a stochastic torque, which inject translational and rotational energy, respectively. We realize that a certain choice of a thermostat of this class can be characterized just by the total intensity and the fraction of noise transferred to the rotational degree of freedom with respect to the translational ones. Firstly, Mpemba effect is characterized by the appearance of a crossing between the temperature curves of the considered samples. Later, an overshoot of the temperature evolution with respect to the steady-state value is observed and the mechanism of Mpemba effect generation is changed. The choice of parameters allows us to design plausible protocols based on these thermostats for generating the initial states to observe the Mpemba-like effect in experiments. In order to obtain explicit results, we use a well-founded Maxwellian approximation for the evolution dynamics and the steady-state quantities. Finally, theoretical results are compared with direct simulation Monte Carlo and molecular dynamics results, and a very good agreement is found.

Intracluster Reaction Dynamics of Ionized Micro-Hydrated Hydrogen Peroxide (H2O2): A Direct Ab Initio Molecular Dynamics Study

Hydrogen peroxide (H2O2) is a unique molecule that is applied in various fields, including energy chemistry, astrophysics, and medicine. H2O2 readily forms clusters with water molecules. In the present study, the reactions of ionized H2O2–water clusters, H2O2+(H2O)n, after vertical ionization of the parent neutral cluster were investigated using the direct ab initio molecular dynamics (AIMD) method to elucidate the reaction mechanism. Clusters with one to five water molecules, H2O2–(H2O)n (n = 1–5), were examined, and the reaction of [H2O2+(H2O)n]ver was tracked from the vertical ionization point to the product state, where [H2O2+(H2O)n]ver is the vertical ionization state (hole is localized on H2O2). After ionization, fast proton transfer (PT) from H2O2+ to the water cluster (H2O)n was observed in all clusters. The HOO radical and H3O+(H2O)n−1 were formed as products. The PT reaction proceeds directly without an activation barrier. The PT times for n = 1–5 were calculated to be 36.0, 9.8, 8.3, 7.7, and 7.1 fs, respectively, at the MP2/6-311++G(d,p) level, indicating that PT in these clusters is a very fast process, and the PT time is not dependent on the cluster size (n), except in the case of n = 1, where the PT time was slightly longer because the bond distance and angle of the hydrogen bond in n = 1 were deformed from the standard structure. The reaction mechanism was discussed based on these results.

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.

Adaptive moving environment for efficient molecular dynamics simulations of high-fluence ion irradiation

Ion-beam processing of materials is widely supported by atomistic simulations by means of molecular dynamics. Although the approach has given several valuable insights, it has a limited operational window in both time and length scales. In particular, for high-fluence ion irradiation with multiple consecutive cascades, the direct molecular dynamics method becomes prohibitively time consuming.In this work, we propose a speed-up algorithm for molecular-dynamics simulations of multiple consecutive collision cascades employing an adaptive moving environment model. In the model, the computational power is primarily focused on calculating the atomic movement in the propagating cascade regions, while thermally equilibrated regions outside the cascades are excluded. Up to a five-fold efficiency increase was seen with the adaptive moving environment compared to classical molecular dynamics, without any significant statistical difference in the results of multiple individual ion-cascade simulations in a heterostructure of alternating Si and SiO2 layers. Simulations of temperature-driven dynamic annealing during high-fluence ion irradiation of Si nanopillars at elevated temperatures using the adaptive moving environment showed similar trends as experiments with respect to temperature dependence. The model is included in the atomistic simulator toolkit, COSIRMA (COmputer Simulator for IRradiation of MAterials), and can easily be enabled through the user-friendly graphical interface.

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.