Molecular Dynamics Software Research Articles

NE-RDFE: A protocol and toolkit for computing relative dissociation free energies with GROMACS between dissimilar molecules using bidirectional nonequilibrium dual topology schemes.

We describe a step-by-step protocol and toolkit for the computation of the relative dissociation free energy (RDFE) with the GROMACS molecular dynamics package, based on a novel bidirectional nonequilibrium alchemical approach. The proposed methodology does not require any intervention on the code and allows computing with good accuracy the RDFE between small molecules with arbitrary differences in volume, charge, and chemical topology. The procedure is illustrated for the challenging SAMPL9 batch of host-guest pairs. The article is supplemented by a detailed online tutorial, available at https://procacci.github.io/vdssb_gromacs/NE-RDFE and by a public Zenodo repository available at https://zenodo.org/record/6982932.

Enhanced path sampling using subtrajectory Monte Carlo moves.

Path sampling allows the study of rare events, such as chemical reactions, nucleation, and protein folding, via a Monte Carlo (MC) exploration in path space. Instead of configuration points, this method samples short molecular dynamics (MD) trajectories with specific start- and end-conditions. As in configuration MC, its efficiency highly depends on the types of MC moves. Since the last two decades, the central MC move for path sampling has been the so-called shooting move in which a perturbed phase point of the old path is propagated backward and forward in time to generate a new path. Recently, we proposed the subtrajectory moves, stone-skipping (SS) and web-throwing, that are demonstrably more efficient. However, the one-step crossing requirement makes them somewhat more difficult to implement in combination with external MD programs or when the order parameter determination is expensive. In this article, we present strategies to address the issue. The most generic solution is a new member of subtrajectory moves, wire fencing (WF), that is less thrifty than the SS but more versatile. This makes it easier to link path sampling codes with external MD packages and provides a practical solution for cases where the calculation of the order parameter is expensive or not a simple function of geometry. We demonstrate the WF move in a double-well Langevin model, a thin film breaking transition based on classical force fields, and a smaller ruthenium redox reaction at the ab initio level in which the order parameter explicitly depends on the electron density.

Scalable hybrid deep neural networks/polarizable potentials biomolecular simulations including long-range effects.

Deep-HP is a scalable extension of the Tinker-HP multi-GPU molecular dynamics (MD) package enabling the use of Pytorch/TensorFlow Deep Neural Network (DNN) models. Deep-HP increases DNNs' MD capabilities by orders of magnitude offering access to ns simulations for 100k-atom biosystems while offering the possibility of coupling DNNs to any classical (FFs) and many-body polarizable (PFFs) force fields. It allows therefore the introduction of the ANI-2X/AMOEBA hybrid polarizable potential designed for ligand binding studies where solvent-solvent and solvent-solute interactions are computed with the AMOEBA PFF while solute-solute ones are computed by the ANI-2X DNN. ANI-2X/AMOEBA explicitly includes AMOEBA's physical long-range interactions via an efficient Particle Mesh Ewald implementation while preserving ANI-2X's solute short-range quantum mechanical accuracy. The DNN/PFF partition can be user-defined allowing for hybrid simulations to include key ingredients of biosimulation such as polarizable solvents, polarizable counter ions, etc.… ANI-2X/AMOEBA is accelerated using a multiple-timestep strategy focusing on the model's contributions to low-frequency modes of nuclear forces. It primarily evaluates AMOEBA forces while including ANI-2X ones only via correction-steps resulting in an order of magnitude acceleration over standard Velocity Verlet integration. Simulating more than 10 μs, we compute charged/uncharged ligand solvation free energies in 4 solvents, and absolute binding free energies of host-guest complexes from SAMPL challenges. ANI-2X/AMOEBA average errors are discussed in terms of statistical uncertainty and appear in the range of chemical accuracy compared to experiment. The availability of the Deep-HP computational platform opens the path towards large-scale hybrid DNN simulations, at force-field cost, in biophysics and drug discovery.

Research Based on the Selection of Regenerative Agents Based on Molecular Simulation New-Old Asphalt Diffusion Behavior

In order to understand the new-old asphalt diffusion behavior, use the Materials Studio (MS) molecular dynamic simulation software to build a matrix asphalt and aging asphalt molecular model based on the average molecular model of the asphalt, and select the appropriate force field and potential energy function. As well as parameters such as boundary conditions, set up as-phalt-regenerative dot-shaped contact models and laminar contact models. From the perspective of molecular perspectives, the diffusion behavior between the old and new asphalt and the re-generative agent is analyzed. The promotion effect of the asphalt diffusion behavior, the simulation results show that the asphalt's internal diffusion process is mutual, and the promotion effect of aromatherapy regeneratives on the spread of new and old asphalt under the same conditions is better than chain regenerative agents and waste oil regeneration agents.

Peridynamics method for problems solve of solids destruction

Abstract. The article investigates the method of peridynamics, which is an alternative approach to solving destruction problems based on integral equations. It is assumed that particles in a continuum interact with each other at a finite distance, as in molecular dynamics. Damage is part of the theory at the level of two-particle interactions, so damage finding and destruction occurs when solving the equation of motion. During this work, bondbased and state-based peridynamics models of destruction used in the Sandia Laboratory were described and implemented within the framework of the MoDyS molecular dynamics software package. In the bond-based model, the defining relationship is the bond stiffness function, which corrects the force of particle-particle interaction and imposes a restriction on the use of the Poisson’s ratio. The state-based model generalizes the bond-based approach and may be applied to materials with any Poisson’s ratio. The relationship of both models is ascertained. Calculation convergence is demonstrated on the example of a one-dimensional elasticity problem. The possibility of using the implemented models for fracture problems is also shown.

Interfacial free energy at the liquid-solid phase boundary calculated by the method of capillary fluctuations

The interfacial free energy at the liquid-solid phase boundary was calculated for some close packed metals (Au, Cu, Al) and silicon. For calculations, the method of capillary fluctuations at temperatures close to the melting point of the studied materials was used. The formation of the two-phase system at the atomic level was carried out using the method of molecular dynamicssimulations by means of the LAMMS classical molecular dynamics package. The interfacial stiffness for materials with different crystallographic orientation of the phase boundary was determined from the analysis of the fluctuation spectrum, which in turn was obtained by Fourier transformation of the oscillations of the interphase boundary. It was shown that the value of the free energy increases in the series Al→Au→Cu→Si

Molecular dynamics based kinetic Monte Carlo simulation for accelerated diffusion

Molecular dynamics (MD) simulations have long been unable to provide significant insights into large-scale processes involving diffusion. Here, we implement a kinetic Monte Carlo (kMC) algorithm to circumvent such issues. We have implemented this approach in the widely used MD simulation package, LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator). Validation of the MD/kMC algorithm is achieved by replicating results predicted by diffusion equations. Comparisons are made between a traditional MC approach and the kMC algorithm for the case of L12 phase growth in Ni-Al alloys. These examples highlight the unique advantages provided by the kMC approach to unlock new capabilities for MD simulations.

DeepCV: A Deep Learning Framework for Blind Search of Collective Variables in Expanded Configurational Space.

We present Deep learning for Collective Variables (DeepCV), a computer code that provides an efficient and customizable implementation of the deep autoencoder neural network (DAENN) algorithm that has been developed in our group for computing collective variables (CVs) and can be used with enhanced sampling methods to reconstruct free energy surfaces of chemical reactions. DeepCV can be used to conveniently calculate molecular features, train models, generate CVs, validate rare events from sampling, and analyze a trajectory for chemical reactions of interest. We use DeepCV in an example study of the conformational transition of cyclohexene, where metadynamics simulations are performed using DAENN-generated CVs. The results show that the adopted CVs give free energies in line with those obtained by previously developed CVs and experimental results. DeepCV is open-source software written in Python/C++ object-oriented languages, based on the TensorFlow framework and distributed free of charge for noncommercial purposes, which can be incorporated into general molecular dynamics software. DeepCV also comes with several additional tools, i.e., an application program interface (API), documentation, and tutorials.

MetalWalls: Simulating electrochemical interfaces between polarizable electrolytes and metallic electrodes.

Electrochemistry is central to many applications, ranging from biology to energy science. Studies now involve a wide range of techniques, both experimental and theoretical. Modeling and simulations methods, such as density functional theory or molecular dynamics, provide key information on the structural and dynamic properties of the systems. Of particular importance are polarization effects of the electrode/electrolyte interface, which are difficult to simulate accurately. Here, we show how these electrostatic interactions are taken into account in the framework of the Ewald summation method. We discuss, in particular, the formal setup for calculations that enforce periodic boundary conditions in two directions, a geometry that more closely reflects the characteristics of typical electrolyte/electrode systems and presents some differences with respect to the more common case of periodic boundary conditions in three dimensions. These formal developments are implemented and tested in MetalWalls, a molecular dynamics software that captures the polarization of the electrolyte and allows the simulation of electrodes maintained at a constant potential. We also discuss the technical aspects involved in the calculation of two sets of coupled degrees of freedom, namely the induced dipoles and the electrode charges. We validate the implementation, first on simple systems, then on the well-known interface between graphite electrodes and a room-temperature ionic liquid. We finally illustrate the capabilities of MetalWalls by studying the adsorption of a complex functionalized electrolyte on a graphite electrode.

Blockchain technology in quantum chemistry: A tutorial review for running simulations on a blockchain

AbstractSimulations of molecules have recently been performed directly on a blockchain virtual computer at atomic resolution. This tutorial review covers the current applications of blockchain technology for molecular modeling in physics, chemistry, and biology, and provides a step‐by‐step tutorial for computational scientists looking to use blockchain computers to simulate physical and scientific processes in general. Simulations of carbon monoxide have been carried out using molecular dynamics software on the Ethereum blockchain in order to facilitate the tutorial.

CHARMM-GUI Enhanced Sampler for various collective variables and enhanced sampling methods.

Enhanced sampling methodologies modifying underlying Hamiltonians can be used for the systems with a rugged potential energy surface that makes it hard to observe convergence using conventional unbiased molecular dynamics (MD) simulations. We present CHARMM-GUI Enhanced Sampler, a web-based tool to prepare various enhanced sampling simulations inputs with user-selected collective variables (CVs). Enhanced Sampler provides inputs for the following nine methods: accelerated MD, Gaussian accelerated MD, conformational flooding, metadynamics, adaptive biasing force, steered MD, temperature replica exchange MD, replica exchange solute tempering 2, and replica exchange umbrella sampling for the method-implemented MD packages including AMBER, CHARMM, GENESIS, GROMACS, NAMD, and OpenMM. Users only need to select a group of atoms via intuitive web-implementation in order to define commonly used nine CVs of interest: center of mass based distance, angle, dihedral, root-mean-square-distance, radius of gyration, distance projected on axis, two types of angles projected on axis, and coordination numbers. The enhanced sampling methods are tested with several biological systems to illustrate their efficiency over conventional MD. Enhanced Sampler with carefully optimized system-dependent parameters will help users to get meaningful results from their enhanced sampling simulations.

Traumatic Brain Injury Biomarkers, Simulations and Kinetics.

This paper reviews the predictive capabilities of blood-based biomarkers to quantify traumatic brain injury (TBI). Biomarkers for concussive conditions also known as mild, to moderate and severe TBI identified along with post-traumatic stress disorder (PTSD) and chronic traumatic encephalopathy (CTE) that occur due to repeated blows to the head during one's lifetime. Since the pathways of these biomarkers into the blood are not fully understood whether there is disruption in the blood-brain barrier (BBB) and the time it takes after injury for the expression of the biomarkers to be able to predict the injury effectively, there is a need to understand the protein biomarker structure and other physical properties. The injury events in terms of brain and mechanics are a result of external force with or without the shrapnel, in the wake of a wave result in local tissue damage. Thus, these mechanisms express specific biomarkers kinetics of which reaches half-life within a few hours after injury to few days. Therefore, there is a need to determine the concentration levels that follow injury. Even though current diagnostics linking biomarkers with TBI severity are not fully developed, there is a need to quantify protein structures and their viability after injury. This research was conducted to fully understand the structures of 12 biomarkers by performing molecular dynamics simulations involving atomic movement and energies of forming hydrogen bonds. Molecular dynamics software, NAMD and VMD were used to determine and compare the approximate thermodynamic stabilities of the biomarkers and their bonding energies. Five biomarkers used clinically were S100B, GFAP, UCHL1, NF-L and tau, the kinetics obtained from literature show that the concentration values abruptly change with time after injury. For a given protein length, associated number of hydrogen bonds and bond energy describe a lower bound region where proteins self-dissolve and do not have long enough half-life to be detected in the fluids. However, above this lower bound, involving higher number of bonds and energy, we hypothesize that biomarkers will be viable to disrupt the BBB and stay longer to be modeled for kinetics for diagnosis and therefore may help in the discoveries of new biomarkers.

Coarse-Graining ddCOSMO through an Interface between Tinker and the ddX Library.

The domain decomposition conductor-like screening model is an efficient way to compute the solvation energy of solutes within a polarizable continuum medium in a linear scaling computational time. Despite its efficiency, the application to very large systems is still challenging. A possibility to further accelerate the algorithm is resorting to coarse-graining strategies. In this paper we present a preliminary interface between the molecular dynamics package Tinker and the ddX library. The interface was used to test a united atom coarse-graining strategy that allowed us to push ddCOSMO to its limits by computing solvation energies on systems with up to 7 million atoms. We first present benchmarks to find an optimal discretization, and then, we discuss the performance and results obtained with fine- and coarse-grained solvation energy calculations.

Review of Electrostatic Force Calculation Methods and Their Acceleration in Molecular Dynamics Packages Using Graphics Processors.

Molecular dynamics (MD) simulations probe the conformational repertoire of macromolecular systems using Newtonian dynamic equations. The time scales of MD simulations allow the exploration of biologically relevant phenomena and can elucidate spatial and temporal properties of the building blocks of life, such as deoxyribonucleic acid (DNA) and protein, across microsecond (μs) time scales using femtosecond (fs) time steps. A principal bottleneck toward extending MD calculations to larger time scales is the long-range electrostatic force measuring component of the naive nonbonded force computation algorithm, which scales with a complexity of (N, number of atoms). In this review, we present various methods to determine electrostatic interactions in often-used open-source MD packages as well as the implementation details that facilitate acceleration of the electrostatic interaction calculation.

Adiabatic Excited-State Molecular Dynamics with an Explicit Solvent: NEXMD-SANDER Implementation.

We present a method to link the Nonadiabatic EXcited-state Molecular Dynamics (NEXMD) package to the SANDER package supplied by AMBERTOOLS to provide excited-state adiabatic quantum mechanics/molecular mechanics (QM/MM) simulations. NEXMD is a computational package particularly developed to perform simulations of the photoexcitation and subsequent nonadiabatic electronic and vibrational energy relaxation in large multichromophoric conjugated molecules involving several coupled electronic excited states. The NEXMD-SANDER exchange has been optimized in order to achieve excited-state adiabatic dynamics simulations of large conjugated materials in a QM/MM environment, such as an explicit solvent. Dynamics of a substituted polyphenylene vinylene oligomer (PPV3-NO2) in vacuum and different explicit solvents has been used as a test case by performing comparative analysis of changes in its optical spectrum, state-dependent conformational changes, and quantum bond orderings. The method has been tested and compared with respect to previous implicit solvent implementations. Also, the impact on the expansion of the QM region by including a variable number of solvent molecules has been analyzed. Altogether, these results encourage future implementations of NEXMD simulations using the same combination of methods.

Molecular mechanisms from reaction coordinate graph enabled multidimensional free energies illustrated on water dimer hydrogen bonding.

Computing the free energies of molecular mechanisms in multidimensional space relies on combinations of geometrically complex reaction coordinates. We show how a graph theory implementation reduces complexity, and illustrate this on the arrangements of hydrogen bonding of a water dimer. The reaction coordinates and forces are computed using graphs that define the dependencies on the atoms in the Free Energy from Adaptive Reaction Coordinate Forces (FEARCF) library. The library can be interfaced with classical molecular dynamics as well as quantum molecular dynamics packages. Multidimensional interdependent reaction coordinates are constructed to produce complex free energy hypersurfaces. The reaction coordinates are graphed from atomic and molecular components to define points, distances, vectors, angles, planes and combinations thereof. The resultant free energy surfaces that are a function of the distance, angles, planes, and so on, can represent molecular mechanisms in reduced dimensions from the component atomic Cartesian coordinate degrees of freedom. The FEARCF library can be interfaced with any molecular package. Here, we demonstrate the link to NWChem to compute a hyperdimensional DFT (aug‐cc‐pVDZ basis set and X3LYP exchange correlation functionals) free energy space of a water dimer. Analysis of the water dimer free energy hypervolume reveals that while the chain and cyclic hydrogen bonding configurations are located in stable minimum energy wells, the bifurcated hydrogen bond configuration is a gateway to instability and dimer dissociation.

Study on the Effect of Processing Parameters on Residual Stresses of Injection Molded Micro-Pillar Array.

As one of the main methods for fabricating microstructured surfaces, micro-injection molding has the advantages of short cycle time, high production efficiency, and the potential for batch manufacturing. However, non-negligible residual stresses inside the molded part could affect the replication quality, dimensions, and physical properties of the microstructure. Therefore, studying the effects of processing parameters on residual stresses is a necessary prerequisite to ensure the successful fabrication of microstructured parts. In this paper, an injection molding simulation model of micro-pillar arrays was developed using molecular dynamics software, and a series of injection molding experiments were conducted. It was found that increasing the mold temperature and melt temperature can reduce the thermal residual stresses and molecular orientation stresses, and effectively improve the uniformity of residual stress distribution. The increase in the packing pressure can make the shear field of flow more intense and increase the molecular orientation stresses, thus making the residual stresses more severe.

Molecular dynamics analysis of the influence of ion implantation parameters on ultra-precision machining of silicon carbide

This paper aims to study the mechanism of ion implantation. Molecular dynamics software was used to simulate the process of 3C-SiC workpiece being ion implanted with different implant parameters. The influence of ion implantation parameters on workpiece damage, as well as on workpiece temperature, potential energy, dislocation defects, surface morphology, cutting force, stress and other processing parameters, is analyzed. The results show that with the increase of dose, the damage range inside the workpiece shows a periodic change, the potential energy and temperature generally show a rising trend, and the internal normal stress exhibits an increasing trend. In addition, the temperature, potential energy, cutting force and stress of the workpiece increase with the increasing inclination angle. Furthermore, the channel effect is obvious when the inclination angle of ion implantation is between 30° and 50°. Alternatively, for inclination angle greater than 70°, a scattering effect is distinctly observed. Moreover, the average cutting force decreases with the increase of injection distance. Lastly, the average temperature is lowest when the ion implantation region is in the middle of the workpiece, and the farther the implantation region is, the faster is the rise in temperature and potential energy.

Computing the Work of Solid-Liquid Adhesion in Systems with Damped Coulomb Interactions via Molecular Dynamics: Approaches and Insights.

Recently, the dry-surface method [Langmuir2016, 31, 8335−8345] has been developed to compute the work of adhesion of solid–liquid and other interfaces using molecular dynamics via thermodynamic integration. Unfortunately, when long-range Coulombic interactions are present in the interface, a special treatment is required, such as solving additional Poisson equations, which is usually not implemented in generic molecular dynamics software, or as fixing some groups of atoms in place, which is undesirable most of the time. In this work, we replace the long-range Coulombic interactions with damped Coulomb interactions, and explore several thermal integration paths. We demonstrate that regardless of the integration path, the same work of adhesion values are obtained as long as the path is reversible, but the numerical efficiency differs vastly. Simple scaling of the interactions is most efficient, requiring as little as 8 sampling points, followed by changing the Coulomb damping parameter, while modifying the Coulomb interaction cutoff length performs worst. We also demonstrate that switching long-range Coulombic interactions to damped ones results in a higher work of adhesion by about 10 mJ/m2 because of slightly different liquid molecule orientation at the solid–liquid interface, and this value is mostly unchanged for surfaces with substantially different Coulombic interactions at the solid–liquid interface. Finally, even though it is possible to split the work of adhesion into van der Waals and Coulomb components, it is known that the specific per-component values are highly dependent on the integration path. We obtain an extreme case, which demonstrates that caution should be taken even when restricting to qualitative comparison.

A thermalized electrokinetics model including stochastic reactions suitable for multiscale simulations of reaction–advection–diffusion systems

We introduce a scheme to simulate the spatial and temporal evolution of the densities of charged species, taking into account diffusion, thermal fluctuations, coupling to a carrier fluid, and chemical reactions. To this end, the diffusive fluxes in the electrokinetic model by Capuani et al. (2004) are supplemented with thermal fluctuations. Chemical reactions are included via an additional source term in the mass balance equation. The diffusion–reaction model is then coupled to a solver for fluctuating hydrodynamics based on the lattice Boltzmann method. This combination is particularly useful for soft matter simulations, due to the ability to couple particles to the lattice-Boltzmann fluid. These could, e.g., be charged colloids or polymers, which then interact with an ion distribution. We describe one implementations based on the automatic code generation tools pystencils and lbmpy, and another one that is contained in the molecular dynamics package ESPResSo and that allows for an easy coupling of particles to the density fields. We validate our implementations by comparing to several known analytic results. Our method can be applied to coarse-grained catalysis problems as well as to many other multi-scale problems that require the coupling of explicit-particle simulations to flow fields, diffusion, and reaction problems in arbitrary geometries.