Molecular Dynamics Algorithm Research Articles

Modeling method of random collision of ellipsoidal particles in mechanical simulation for high energy solid propellant

Abstract It is based on the molecular dynamics algorithm of elliptical optimization to establish a mesostructure model of oval particles of high-energy solid propellant considering gradation and controllable filling rate according to the mesostructure characteristics of high-energy solid propellant. The bilinear cohesion model was applied to describe the interface damage behavior and carry out the high-energy propellant tensile test at three rates. The rationality of the elliptical particle model has been verified between the comparison of the calculated results with the experimental results established to describe the micromechanical behavior of high-energy propellants, which is to reveal the mesoscopic damage law of high-energy solid propellant under different tensile rates in this paper.

SWAP algorithm for lattice spin models.

We adapted the SWAP molecular dynamics algorithm for use in lattice Ising spin models. We dressed the spins with a randomly distributed length and we alternated long-range spin exchanges with conventional single spin flip Monte Carlo updates, both accepted with a stochastic rule which respects detailed balance. We show that this algorithm, when applied to the bidimensional Edwards-Anderson model, speeds up significantly the relaxation at low temperatures and manages to find ground states with high efficiency and little computational cost. The exploration of spin models should help in understanding why SWAP accelerates the evolution of particle systems and sheds light on relations between dynamics and free-energy landscapes.

Integrating Molecular Dynamics and Machine Learning Algorithms to Predict the Functional Profile of Kinase Ligands.

The modulation of protein function via designed small molecules is providing new opportunities in chemical biology and medicinal chemistry. While drugs have traditionally been developed to block enzymatic activities through active site occupation, a growing number of strategies now aim to control protein functions in an allosteric fashion, allowing for the tuning of a target's activation or deactivation via the modulation of the populations of conformational ensembles that underlie its function. In the context of the discovery of new active leads, it would be very useful to generate hypotheses for the functional impact of new ligands. Since the discovery and design of allosteric modulators (inhibitors/activators) is still a challenging and often serendipitous target, the development of a rapid and robust approach to predict the functional profile of a new ligand would significantly speed up candidate selection. Herein, we present different machine learning (ML) classifiers to distinguish between potential orthosteric and allosteric binders. Our approach integrates information on the chemical fingerprints of the ligands with descriptors that recapitulate ligand effects on protein functional motions. The latter are derived from molecular dynamics (MD) simulations of the target protein in complex with orthosteric or allosteric ligands. In this framework, we train and test different ML architectures, which are initially probed on the classification of orthosteric versus allosteric ligands for cyclin-dependent kinases (CDKs). The results demonstrate that different ML methods can successfully partition allosteric versus orthosteric effectors (although to different degrees). Next, we further test the models with FDA-approved CDK drugs, not included in the original dataset, as well as ligands that target other kinases, to test the range of applicability of these models outside of the domain on which they were developed. Overall, the results show that enriching the training dataset with chemical physics-based information on the protein-ligand dynamic cross-talk can significantly expand the reach and applicability of approaches for the prediction and classification of the mode of action of small molecules.

The explicit bonding reaction ensemble Monte Carlo method.

We present the explicit bonding Reaction ensemble Monte Carlo (eb-RxMC) method, designed to sample reversible bonding reactions in macromolecular systems in thermodynamic equilibrium. Our eb-RxMC method is based on the reaction ensemble method; however, its implementation differs from the latter by the representation of the reaction. In the eb-RxMC implementation, we are adding or deleting bonds between existing particles, instead of inserting or deleting particles with different chemical identities. This new implementation makes the eb-RxMC method suitable for simulating the formation of reversible linkages between macromolecules, which would not be feasible with the original implementation. To enable coupling of our eb-RxMC algorithm with molecular dynamics algorithm for the sampling of the configuration space, we biased the sampling of reactions only within a certain inclusion radius. We validated our algorithm using a set of ideally behaving systems undergoing dimerization and polycondensation reactions, for which analytical results are available. For dimerization reactions with various equilibrium constants and initial compositions, the degree of conversion measured in our simulations perfectly matched the reference values given by the analytical equations. We also showed that this agreement is not affected by the arbitrary choice of the inclusion radius or the stiffness of the harmonic bond potential. Next, we showed that our simulations can correctly match the analytical results for the distribution of the degree of polymerization and end-to-end distance of ideal chains in polycondensation reactions. Altogether, we demonstrated that our eb-RxMC simulations correctly sample both reaction and configuration spaces of these reference systems, opening the door to future simulations of more complex interacting macromolecular systems.

Multiple Parameter Replica Exchange Gaussian Accelerated Molecular Dynamics for Enhanced Sampling and Free Energy Calculation of Biomolecular Systems.

This study introduces a novel method named multiple parameter replica exchange Gaussian accelerated molecular dynamics (MP-Rex-GaMD), building on the Gaussian accelerated molecular dynamics (GaMD) algorithm. GaMD enhances sampling and retrieves free energy information for biomolecular systems by adding a harmonic boost potential to smooth the potential energy surface without the need for predefined reaction coordinates. Our innovative approach advances the acceleration power and energetic reweighting accuracy of GaMD by incorporating a replica exchange algorithm that enables the exchange of multiple parameters, including the GaMD boost parameters of force constant and energy threshold, as well as temperature. Applying MP-Rex-GaMD to the three model systems of dialanine, chignolin, and HIV protease, we demonstrate its superior capability over conventional molecular dynamics and GaMD simulations in exploring protein conformations and effectively navigating various biomolecular states across energy barriers. MP-Rex-GaMD allows users to accurately map free energy landscapes through energetic reweighting, capturing the ensemble of biomolecular states from low-energy conformations to rare high-energy transitions within practical computational time scales.

Hamiltonian limited valence model for liquid polyamorphism

Liquid-liquid phase transitions have been found experimentally or by computer simulations in many compounds such as water, hydrogen, sulfur, phosphorus, carbon, silica, and silicon. Limited valence model implemented via event-driven molecular dynamics algorithm provides a simple generic mechanism for the liquid-liquid phase transitions in all these diverse cases. Here, we introduce a variant of the limited valence model with a well defined Hamiltonian, i.e., a unique algorithm by which the potential energy of the system of particles can be computed solely from the coordinates of the particles and is thus equivalent to a complex multi-body potential. We present several examples of the model which can be used to reproduce liquid--liquid phase transition in systems with maximum valence z = 1 (hydrogen), z = 2 (sulfur) and z = 4 (water), where z is the maximum number of bonds an atom is allowed to have. For z = 1, we find a set of parameters for which the system has a liquid-liquid and an isostructural solid-solid critical points. For z = 4, we find a set of parameters for which the phase diagram resembles that of water with a wide region of negative thermal expansion coefficient (density anomaly) extending into the metastable region of negative pressures. The limited valence model can be modified to forbid not only too large valences but also too low valences. In the case of sulfur, we forbid the formation of monomers, thus restricting the valence v of an atom to be within an interval 1 = vmin ≤ v ≤ vmax ≡ z = 2.

Energy-Conserving Surface Hopping for Auger Processes.

Auger-type processes are ubiquitous in nanoscale materials because quantum confinement enhances Coulomb interactions, and there exist large densities of states. Modeling Auger processes requires the modification of nonadiabatic (NA) molecular dynamics algorithms to include transitions caused by both NA and Coulomb couplings. The system is split into quantum and classical subsystems, e.g., electrons and vibrations, and as a result, energy conservation becomes nontrivial. In surface hopping, an electronic transition induced by NA coupling is accompanied by a classical velocity readjustment to ensure conservation of the total quantum-classical energy. A different treatment is needed for Auger transitions driven by Coulomb interactions. We develop a nonadiabatic molecular dynamics methodology that meticulously differentiates the energy redistribution accompanying hops induced by the NA coupling and the Coulomb interaction and correctly conserves the total energy at each transition. If the transition is driven by a Coulomb interaction, the hop energy is redistributed within the quantum electronic subsystem only. If the transition is NA, the energy is redistributed between the quantum and classical subsystems. Properly maintaining energy conservation for both types of transitions is crucial to generate a correct order of events, obtain accurate transition times, maintain a proper statistical distribution of state populations, and reach thermodynamic equilibrium. We test the method with biexciton annihilation and Auger-assisted hot electron relaxation in a CdSe quantum dot. The sequence of Auger and phonon-driven processes and the calculated time scales are in excellent agreement with the experimental results. The developed approach can be coupled with any surface-hopping method and provides a crucial practical advance to study charge-carrier dynamics in the nanoscale and condensed matter systems.

Study on interface characteristics of TC4/7075 explosive welding based on molecular dynamics algorithm

The preparation of TC4/7075 composites for ultra-thin flyer plates is a significant challenge in the field of explosive welding. A weldability window for multi-layer metal explosive welding and a configuration for ultra-thin flyer plate explosive welding were established in this study. TC4/7075 composite materials were successfully prepared with a flyer plate thickness of only 0.3 mm. An analysis was conducted on the material bonding ability, element diffusion, crystal evolution, and microscopic morphology during the explosive welding process of TC4/7075, utilizing the weldability window, molecular dynamics algorithm, and electron backscattered diffraction testing. The results show that various dislocations are present at the interface, predominantly 1/6 ⟨112⟩ dislocations. Element diffusion primarily occurs during the unloading stage at high temperature and zero external pressure; the interface has the best bonding ability when titanium exhibits FCC lattice structure. The prepared composite material demonstrates high intra-grain and grain boundary stresses; the rolling texture is observed on the aluminum side while an equiaxed twin structure forms on the titanium side due to interactions between stacking faults, Stair-rod dislocations, and Hirth immovable dislocations.

A Reverse Nonequilibrium Molecular Dynamics Algorithm for Coupled Mass and Heat Transport in Mixtures.

We present a new method for introducing stable nonequilibrium concentration gradients in molecular dynamics simulations of mixtures. This method extends earlier reverse nonequilibrium molecular dynamics (RNEMD) methods, which use kinetic energy scaling moves to create temperature or velocity gradients. In the new scaled particle flux (SPF-RNEMD) algorithm, energies and forces are computed simultaneously for a molecule existing in two nonadjacent regions of a simulation box, and the system evolves under a linear combination of these interactions. A continuously increasing particle scaling variable is responsible for the migration of the molecule between the regions as the simulation progresses, allowing for simulations under an applied particle flux. To test the method, we investigate diffusivity in mixtures of identical but distinguishable particles and in a simple mixture of multiple Lennard-Jones particles. The resulting concentration gradients provide Fick diffusion constants for mixtures. We also discuss using the new method to obtain coupled transport properties using simultaneous particle and thermal fluxes to compute the temperature dependence of the diffusion coefficient and activation energies for diffusion from a single simulation. Lastly, we demonstrate the use of this new method in interfacial systems by computing the diffusive permeability of a molecular fluid moving through a nanoporous graphene membrane.

Constant-pH Simulations with the Polarizable Atomic Multipole AMOEBA Force Field.

Accurately predicting protein behavior across diverse pH environments remains a significant challenge in biomolecular simulations. Existing constant-pH molecular dynamics (CpHMD) algorithms are limited to fixed-charge force fields, hindering their application to biomolecular systems described by permanent atomic multipoles or induced dipoles. This work overcomes these limitations by introducing the first polarizable CpHMD algorithm in the context of the Atomic Multipole Optimized Energetics for Biomolecular Applications (AMOEBA) force field. Additionally, our implementation in the open-source Force Field X (FFX) software has the unique ability to handle titration state changes for crystalline systems including flexible support for all 230 space groups. The evaluation of constant-pH molecular dynamics (CpHMD) with the AMOEBA force field was performed on 11 crystalline peptide systems that span the titrating amino acids (Asp, Glu, His, Lys, and Cys). Titration states were correctly predicted for 15 out of the 16 amino acids present in the 11 systems, including for the coordination of Zn2+ by cysteines. The lone exception was for a HIS-ALA peptide where CpHMD predicted both neutral histidine tautomers to be equally populated, whereas the experimental model did not consider multiple conformers and diffraction data are unavailable for rerefinement. This work demonstrates the promise polarizable CpHMD simulations for pKa predictions, the study of biochemical mechanisms such as the catalytic triad of proteases, and for improved protein-ligand binding affinity accuracy in the context of pharmaceutical lead optimization.

Electrochemical rewiring through quantum conductance effects in single metallic memristive nanowires.

Memristive devices have been demonstrated to exhibit quantum conductance effects at room temperature. In these devices, a detailed understanding of the relationship between electrochemical processes and ionic dynamic underlying the formation of atomic-sized conductive filaments and corresponding electronic transport properties in the quantum regime still represents a challenge. In this work, we report on quantum conductance effects in single memristive Ag nanowires (NWs) through a combined experimental and simulation approach that combines advanced classical molecular dynamics (MD) algorithms and quantum transport simulations (DFT). This approach provides new insights on quantum conductance effects in memristive devices by unravelling the intrinsic relationship between electronic transport and atomic dynamic reconfiguration of the nanofilment, by shedding light on deviations from integer multiples of the fundamental quantum of conductance depending on peculiar dynamic trajectories of nanofilament reconfiguration and on conductance fluctuations relying on atomic rearrangement due to thermal fluctuations.

Non-equilibrium molecular dynamics of steady-state fluid transport through a 2D membrane driven by a concentration gradient.

We use a novel non-equilibrium algorithm to simulate steady-state fluid transport through a two-dimensional (2D) membrane due to a concentration gradient by molecular dynamics (MD) for the first time. We confirm that, as required by the Onsager reciprocal relations in the linear-response regime, the solution flux obtained using this algorithm agrees with the excess solute flux obtained from an established non-equilibrium MD algorithm for pressure-driven flow. In addition, we show that the concentration-gradient-driven solution flux in this regime is quantified far more efficiently by explicitly applying a transmembrane concentration difference using our algorithm than by applying Onsager reciprocity to pressure-driven flow. The simulated fluid fluxes are captured with reasonable quantitative accuracy by our previously derived continuum theory of concentration-gradient-driven fluid transport through a 2D membrane [D. J. Rankin, L. Bocquet, and D. M. Huang, J. Chem. Phys. 151, 044705 (2019)] for a wide range of solution and membrane parameters, even though the simulated pore sizes are only several times the size of the fluid particles. The simulations deviate from the theory for strong solute-membrane interactions relative to thermal energy, for which the theoretical approximations breakdown. Our findings will be beneficial for a molecular-level understanding of fluid transport driven by concentration gradients through membranes made from 2D materials, which have diverse applications in energy harvesting, molecular separations, and biosensing.

Atomic motion behavior calculation and bonding mechanism analysis of explosive welding of high-strength and high-hardness titanium alloy Ti6Al4V/aluminum alloy 7075

The explosive welding technology of titanium alloy Ti6Al4V/aluminum alloy 7075 is the most valuable research problem in the field of Ti/Al explosive welding. In this paper, the smooth particle hydrodynamics and molecular dynamics algorithms are used to conduct multi-scale numerical calculations on the explosive welding of Ti6Al4V and 7075. The interface morphology, interface atomic movement, lattice changes and crystal defects are comprehensively analyzed. With the calculation as a reference, the Ti6Al4V/7075 composite with a flat interface is prepared. The calculation results indicate that with the transformation of titanium aluminum crystal structure, a large number of various dislocations, mainly 1/6〈110〉dislocations, are generated on the titanium side; During explosive welding, uniform and stable element diffusion occurred. Element diffusion mainly occurred in the unloading stage when the pressure was 0. The melting zone and ablation zone were composed of TiAl and TiAl3, respectively; The fundamental reason for the combination of Ti6Al4V/7075 is the diffusion of elements and the solid-phase binding of atoms.

Over- and Undercoordinated Atoms as a Source of Electron and Hole Traps in Amorphous Silicon Nitride (a-Si3N4).

Silicon nitride films are widely used as the charge storage layer of charge trap flash (CTF) devices due to their high charge trap densities. The nature of the charge trapping sites in these materials responsible for the memory effect in CTF devices is still unclear. Most prominently, the Si dangling bond or K-center has been identified as an amphoteric trap center. Nevertheless, experiments have shown that these dangling bonds only make up a small portion of the total density of electrical active defects, motivating the search for other charge trapping sites. Here, we use a machine-learned force field to create model structures of amorphous Si3N4 by simulating a melt-and-quench procedure with a molecular dynamics algorithm. Subsequently, we employ density functional theory in conjunction with a hybrid functional to investigate the structural properties and electronic states of our model structures. We show that electrons and holes can localize near over- and under-coordinated atoms, thereby introducing defect states in the band gap after structural relaxation. We analyze these trapping sites within a nonradiative multi-phonon model by calculating relaxation energies and thermodynamic charge transition levels. The resulting defect parameters are used to model the potential energy curves of the defect systems in different charge states and to extract the classical energy barrier for charge transfer. The high energy barriers for charge emission compared to the vanishing barriers for charge capture at the defect sites show that intrinsic electron traps can contribute to the memory effect in charge trap flash devices.

Load–displacement relation and gap distribution between rough surfaces: Partial differential equations approach

We develop a theoretical model to predict the load–displacement relation and probability density function for gaps between contacting rough surfaces. We derive partial differential equations from the previous model (Joe et al., 2018), and extend them to the non-adhesive contact problem. The predictions of the present theory are compared with numerical results using the Green’s function molecular dynamics algorithm, and a good agreement is obtained.

Translational and rotational non-Gaussianities in homogeneous freely evolving granular gases.

The importance of roughness in the modeling of granular gases has been increasingly considered in recent years. In this paper, a freely evolving homogeneous granular gas of inelastic and rough hard disks or spheres is studied under the assumptions of the Boltzmann kinetic equation. The homogeneous cooling state is studied from a theoretical point of view using a Sonine approximation, in contrast to a previous Maxwellian approach. A general theoretical description is done in terms of d_{t} translational and d_{r} rotational degrees of freedom, which accounts for the cases of spheres (d_{t}=d_{r}=3) and disks (d_{t}=2,d_{r}=1) within a unified framework. The non-Gaussianities of the velocity distribution function of this state are determined by means of the first nontrivial cumulants and by the derivation of non-Maxwellian high-velocity tails. The results are validated by computer simulations using direct simulation Monte Carlo and event-driven molecular dynamics algorithms.

MD-Bench: A performance-focused prototyping harness for state-of-the-art short-range molecular dynamics algorithms

Molecular dynamics (MD) simulations provide considerable benefits for the investigation and experimentation of systems at atomic level. Their usage is widespread into several research fields, but their system size and timescale are crucially limited by the available computing power. Performance engineering of MD kernels is therefore critical to understand their bottlenecks and investigate possible improvements. For that reason, we developed MD-Bench, a performance-focused prototyping harness for short-range MD kernels that implements state-of-the-art algorithms from multiple production applications such as LAMMPS and GROMACS. The MD-Bench source code is simple, understandable, and extensible, and therefore well suited for benchmarking, teaching, and researching MD algorithms. In this paper we introduce MD-Bench, describe its design, structure, and implemented algorithms. Finally, we show five use-cases of MD-Bench and describe how these are useful to gain a deeper understanding of the performance of MD kernels.

Bio-ESMD: A Data Centric Implementation for Large-Scale Biological System Simulation on Sunway TaihuLight Supercomputer

Molecular dynamics (MD) simulations of biological systems are playing an increasingly important role in the research of pathogens and drugs. Most MD methods for biological simulations rely on the listed bonds which interact among specific groups of atoms identified by atom tags (unique atom tags regardless the storage location). However, efficient mapping of tags to atom locations is often challenging on modern many-core processors because data locality can not always be guaranteed for large-scale systems. In this paper, we present Bio-ESMD, a new MD implementation supporting listed bonds. Bio-ESMD is designed and developed based on our previously designed ESMD framework for many-core processors. In Bio-ESMD, we have introduced a data-centric approach for refactoring MD algorithms by reorganizing the cell list data structure to adopt bond lists with guaranteed data locality. Our implementation achieves speedups of over two compared to SW_GROMACS on Sunway TaihuLight. Furthermore, Bio-ESMD can simulate a system of 308.8 million atoms at 1.33 ns/day or 14.44 million atoms at 17.28 ns/day with linear weak scaling efficiency.

Linear and Nonlinear Viscoelasticities of Dissociative and Associative Covalent Adaptable Networks: Discrepancies and Limits

Covalent adaptable networks (CANs) can be classified into dissociative (Diss-CANs) and associative (Asso-CANs) networks according to the exchange mechanism of covalent bonds. We simulate the exchange reaction by the hybrid Monte Carlo–molecular dynamics (hybrid MC/MD) algorithm, aiming to discover the connection and difference between Diss-CANs and Asso-CANs in viscoelasticity behavior. In the linear regime, a major difference originating from the cross-linking density is reflected in the pre-exponential factor τs0 of the characteristic relaxation time τs. For nonlinear rheology, Diss-CANs show a faster shear thinning behavior under steady shear, while Asso-CANs have a stronger strain hardening under the shear rate start-up. The physics behind the phenomenon results from the different chain conformations and configurations related to the exchange mechanism. Compared with Diss-CANs, the inability for sticker dissociation of Asso-CANs generates a slower relaxation under shear, leading to less chain orientation and tumbling. Meanwhile, we find that multiscale relaxation times obtained from linear viscoelasticity (LVE) can be crucial limits in nonlinear applications for associative polymers (APs). Our work strongly deepens the understanding of APs in terms of both linear and nonlinear viscoelasticities.

Learning Reactive Islands of the Voter97 System

In this paper, we assess the effectiveness of a widely used machine learning technique, support vector machines (SVM) for computing reactive islands in a benchmark system for testing molecular dynamics algorithms, the Voter97 model. Reactive islands are the phase space geometrical structure that mediate chemical reactions dynamics. The Voter97 model contains particular challenges for reaction dynamics methods as the reactant and product potential wells are separated by an intermediate well. We show that SVM can accurately compute the reactive islands in the Voter97 model and we assess the accuracy and the computational effort of the approach by comparing it with brute force methods for computing the reactive islands.