Molecular Dynamics Simulation Program Research Articles

MolDy: molecular dynamics simulation made easy.

Molecular dynamics (MD) is a computational experiment that is crucial for understanding the structure of biological macro and micro molecules, their folding, and the inter-molecular interactions. Accurate knowledge of these structural features is the cornerstone in drug development and elucidating macromolecules functions. The open-source GROMACS biomolecular MD simulation program is recognized as a reliable and frequently used simulation program for its precision. However, the user requires expertise, and scripting skills to carrying out MD simulations. We have developed an end-to-end interactive MD simulation application, MolDy for Gromacs. This front-end application provides a customizable user interface integrated with the Python and Perl-based logical backend connecting the Linux shell and Gromacs software. The tool performs analysis and provides the user with simulation trajectories and graphical representations of relevant biophysical parameters. The advantages of MolDy are (i) user-friendly, does not requiring the researcher to have prior knowledge of Linux; (ii) easy installation by a single command; (iii) freely available for academic research; (iv) can run with minimum configuration of operating systems; (v) has valid default prefilled parameters for beginners, and at the same time provides scope for modifications for expert users. MolDy is available freely as compressed source code files with user manual for installation and operation on GitHub: https://github.com/AIBResearchMolDy/Moldyv01.git and on https://aibresearch.com/innovations.

Molecular dynamics simulation of 100 keV argon ion radiation effects in metals and comparison with stopping and range of ions in matter code

Simulation is essential to analyse radiation-induced effects because energetic ions have short transit times within materials (Ar ions of 1 keV and 100 keV have transit times of 0.29 picoseconds and 29 femtoseconds respectively, through a 25 nm Al film). The binary collision approximation (BCA) is extensively used in simulations to study damage caused by energetic ions due to its computational efficiency. Stopping and Range of Ions in Matter (SRIM), a widely used Monte Carlo software using BCA, provides accurate ion ranges consistent with experiments. Molecular dynamics simulations are powerful computational method, where classical equations of motion are solved to track atomic movements from radiation damage, but computational resources required are heavier for this method. In our current effort, the energy transfer from 100 keV Ar+ to metal targets is analysed as a function of target properties, including melting temperature and atomic density, using these two methodologies. A detailed investigation of 100 keV Ar+ ion impact on metal targets is done using the classical Molecular Dynamics simulation program ‘MDRange’, which uses the Recoil Interaction Analysis (RIA) approximation, where the interaction between the ion and its closest neighbours is considered via a two-body potential. The correlation of energy transferred to the metal with the target properties is examined considering the physical approximations of each technique.

Reactive molecular dynamics simulation on the disintegration of Kapton and Upilex-S during atomic oxygen impact

Polyimides are polymeric materials that are widely used in spacecraft applications owing to their unique properties. However, exposure to a low-Earth-orbit environment containing atomic oxygen (AO) results in the disintegration of polymeric materials on the surface of spacecraft, thereby affecting the lifespan. Along with the development of theoretical research, the reactive force-field (ReaxFF) interatomic potential has become a robust computational method for exploring, developing and optimizing the material properties. This study employs the ReaxFF reactive-force-field molecular dynamics simulation (ReaxFF MD) program to investigate and compare the performance of two typical polyimide materials, Kapton and Upilex-S, under the impact of AO. Various aspects such as variations in the temperature, mass loss, decomposition products, and damage propagation depth were examined. Although these materials have similar elemental composition (C/H/O/N), they have different structures. Our results indicate that AO is initially adsorbed on the surfaces of both Kapton and Upilex-S. The continuous impact of AO leads to chemical reactions between AO and Kapton/Upilex-S. Erosion proceeds from the surface toward the interior of the materials. Similar to the findings of Experiment 2 conducted by the Materials International Space Station, our results also reveal that Upilex-S exhibits a lower mass loss and erosion yield than Kapton under the same AO conditions. This difference is primarily attributed to the distinct molecular structures of both Kapton and Upilex-S. Our study could provide valuable technical support for the extensive application of Upilex-S in spacecraft.

An exa-scale high-performance molecular dynamics simulation program: MODYLAS.

A new version of the highly parallelized general-purpose molecular dynamics (MD) simulation program MODYLAS with high performance on the Fugaku computer was developed. A benchmark test using Fugaku indicated highly efficient communication, single instruction, multiple data (SIMD) processing, and on-cache arithmetic operations. The system's performance deteriorated only slightly, even under high parallelization. In particular, a newly developed minimum transferred data method, requiring a significantly lower amount of data transfer compared to conventional communications, showed significantly high performance. The coordinates and forces of 101 810 176 atoms and the multipole coefficients of the subcells could be distributed to the 32 768 nodes (1 572 864 cores) in 2.3ms during one MD step calculation. The SIMD effective instruction rates for floating-point arithmetic operations in direct force and fast multipole method (FMM) calculations measured on Fugaku were 78.7% and 31.5%, respectively. The development of a data reuse algorithm enhanced the on-cache processing; the cache miss rate for direct force and FMM calculations was only 2.74% and 1.43%, respectively, on the L1 cache and 0.08% and 0.60%, respectively, on the L2 cache. The modified MODYLAS could complete one MD single time-step calculation within 8.5ms for the aforementioned large system. Additionally, the program contains numerous functions for material research that enable free energy calculations, along with the generation of various ensembles and molecular constraints.

Successes and Failures of Static Aptamer-Target 3D Docking Models.

While Molecular Dynamics simulation programs are probably superior for predicting the binding and affinity of aptamers and their cognate ligands, such molecular dynamics programs require more computing power and analysis time than static docking programs that are more widely accessible to the scientific community on the internet. Static docking programs can be used to investigate the geometric fit of rigid DNA or RNA aptamer 3D structures and their ligands to aid in predicting the relative affinities and cross-reactivity of various potential ligands. Herein, the author describes when such static 3D docking analysis has worked well to produce useful predictions or confirmation of high-affinity aptamer interactions or successful aptamer beacon behavior and when it has not worked well. The analysis of why failures may occur with static 3D computer models is also discussed.

FSATOOL 2.0: An integrated molecular dynamics simulation and trajectory data analysis program.

Molecular dynamics simulation is important in the computational study of the biomolecules. In this paper, we upgrade our previous FSATOOL to version 2.0. It is no longer a plugin as before. Besides the existed enhanced sampling and Markov state model analysis module, FSATOOL 2.0 has three new features now. First, it contains a molecular dynamics simulation engine on both CPU and GPU device. The engine works with an embedded enhanced sampling module. Second, it can do the free energy calculation by various practical methods, including the weighted histogram analysis method and Gaussian mixture model. Third, it has many subroutines to process the trajectory data, such as principal component analysis, time-structure based independent component analysis, contact analysis, and Φ-value analysis. Most importantly, all these calculations are integrated into one package. The trajectory data format is compatible with all the modules. With a proper input parameter file, users can do the molecular dynamics simulation and data analysis work by only a few simplified commands. The capabilities and theoretical backgrounds of FSATOOL 2.0 are introduced in the paper.

Assembly of a Phenylalanine Nanotube by the use of Molecular Dynamics Manipulator

Studies of the processes of self-organization and self-assembly of various complexly organized (including spiral) structures based on amino acids intensively carried out in recent years. Various methods of molecular modeling, including molecular dynamics (MD) methods, are developed. In this paper, we propose a new approach for a relatively simple technique for conducting MD simulation (MDS) of various molecular nanostructures, determining the trajectory of the MD run and forming the final structure: a molecular dynamic manipulator (MD manipulator). It is an imitation of the operation of an existing or imaginary device or structure by applying force to the existing initial structure in order to obtain a new final structure, having the same chemical composition, but with a different geometry (topology). The PUMA-CUDA software package was applied as the main MD modeling program, which uses the physics of the PUMA software package, developed by the laboratory headed by N.K. Balabaev. Using this MDS tool, you can investigate the formation of helical structures from a linear sequence of any amino acids variation. As an example, the applicability of the developed algorithm for assembling nanotubes from linear phenylalanine (Phe) chains of different chirality (left L-Phe and right D-Phe) is considered by including additional force effects in the molecular dynamics simulation program for these structures. During the MD run, the applied actions, which are the same for the left and right helices of the formed nanotubes, lead the system to an α-helix structure. The work was carried out in an interactive mode using a number of additional programs, incl. trajectory analyzer program MD (TAMD). As a result, the modes that are most adequate for the formation of nanotubes of the right chirality D from the initial L-Phe monomer and nanotubes of the left chirality L from the D-Phe amino acid monomer were determined. The results obtained were compared with data from other works on modeling similar nanotubes of different chirality and experimental data. These are fully in line with the law of change in sign of chirality of molecular structures with complication of their hierarchical level of organization. The molecular animation of the assembly of a left-chiral nanotube from D-monomers is freely available at: http://lmd.impb.ru/Supplementary/PHE.avi.

Simulation and experiment of the cooling effect of trapped ion by pulsed laser* *Project supported by the National Key Research and Development Program of China (Grant No. 2017YFA0304401), the National Development Project for Major Scientific Research Facility, China (Grant No. ZDYZ2012-2), the National Natural Science Foundation of China (Grant Nos. 11774388 and 11634013), the Strategic Priority Research Program of the

We investigate the process of pulsed laser cooling using a self-constructed molecular dynamics simulation (MD-Simulation) program. We simulate the Doppler cooling process and pulsed laser Doppler cooling process of a single 40Ca+ ion, and the comparison with the experimental results shows that this self-constructed MD-Simulation program works well in the weak laser intensity situation. Furthermore, we analyze the pulsed laser Doppler cooling process of a single 27Al+ ion. This program can be used to analyze the molecular dynamic process of various situations of Doppler cooling in an ion trap, which could give predictions and experimental guidance.

The interaction and mechanism between threonine-montmorillonite composite and Pb2+ or Cu2+: Experimental study and theory calculation

The maximum adsorption capacities of threonine-montmorillonite (Thr–Mt) composite to Pb2+ and Cu2+ in aqueous solution were higher than that of original sodium-montmorillonite (Na–Mt). Molecular dynamics (MD) simulation and Multiwfn wavefunction program based on density functional theory (DFT) were exploited to probe the interaction and mechanism between threonine-montmorillonite composite and Pb2+ or Cu2+. Moreover, for the first time, Independent Gradient Model (IGM) analysis was utilized to decode the interaction essence among Thr, Mt. and Pb2+ or Cu2+ to reveal the removal mechanism of Thr-Mt to Pb2+ and Cu2+. MD results indicated that the arrangement of Thr+ presented a preferred almost vertical orientation, which was consistent with the X-ray diffraction (XRD) result. There was hydrogen bonding between carboxyl group or hydroxy group of Thr+ and Mt. surface. Multiwfn calculation implied that the minimum points of ESP and ALIE were located near the O atoms of hydroxyl group and carboxyl group, which were most likely providing the lone pair electron to form a complex with the metal ion,which was accordance with the result of basin analysis. IGM analysis vividly revealed that O atoms of carboxyl group and hydroxyl group were binding sites, and the interaction between Pb2+ and Na–Mt or Pb2+ and Thr-Mt was stronger than that of Cu2+, which was well consistent with the experimental conclusion.

Molecular dynamics simulation of single-walled carbon nanotubes based on GPU

A graphics processing unit (GPU)-based molecular dynamics simulation of single-walled carbon nanotube (SWCNT) molecules is proposed. Under the background of GPU technology, the differential unit and stress of an SWCNT molecule are analyzed first. The stress of the SWCNT molecule at a distance from the end of the SWCNT molecule is obtained by integration, and the calculation of the critical length of stress transfer at the interface of the SWCNT molecule is completed. The SWCNT molecule is divided into two steps by temperature control and using a molecular dynamics simulation program. In this paper, GPU technology is used to derive the velocity Verlet algorithm on the basis of the Verlet algorithm, which completes the optimization design of the numerical algorithm of the differential equation of motion of SWCNTs. Finally, through the SWCNTs, the optimization design of the numerical algorithm of the differential equation of motion of SWCNTs is completed. The molecular dynamics simulation of SWCNTs based on GPU is realized by the equilibrium molecular dynamics simulation of nanotubes and the non-equilibrium molecular dynamics simulation of SWCNTs. The simulation results show that the molecular dynamics simulation method based on GPU can reduce the eccentricity of SWCNT molecules.

Computer-aided synthesis of dapsone-phytochemical conjugates against dapsone-resistant Mycobacterium leprae

Leprosy continues to be the belligerent public health hazard for the causation of high disability and eventual morbidity cases with stable prevalence rates, even with treatment by the on-going multidrug therapy (MDT). Today, dapsone (DDS) resistance has led to fear of leprosy in more unfortunate people of certain developing countries. Herein, DDS was chemically conjugated with five phytochemicals independently as dapsone-phytochemical conjugates (DPCs) based on azo-coupling reaction. Possible biological activities were verified with computational chemistry and quantum mechanics by molecular dynamics simulation program before chemical synthesis and spectral characterizations viz., proton-HNMR, FTIR, UV and LC-MS. The in vivo antileprosy activity was monitored using the ‘mouse-foot-pad propagation method’, with WHO recommended concentration 0.01% mg/kg each DPC for 12 weeks, and the host-toxicity testing of the active DPC4 was seen in cultured-human-lymphocytes in vitro. One-log bacilli cells in DDS-resistant infected mice footpads decreased by the DPC4, and no bacilli were found in the DDS-sensitive mice hind pads. Additionally, the in vitro host toxicity study also confirmed that the DCP4 up to 5,000 mg/L level was safety for oral administration, since a minor number of dead cells were found in red color under a fluorescent microscope. Several advanced bioinformatics tools could help locate the potential chemical entity, thereby reducing the time and resources required for in vitro and in vitro tests. DPC4 could be used in place of DDS in MDT, evidenced from in vivo antileprosy activity and in vitro host toxicity study.

Insights into product release dynamics through structural analyses of thymidylate kinase

Several studies on enzyme catalysis have pointed out that the product release event could be a rate limiting step. In this study, we have compared the release event of two products, Adenosine di-phosphate (ADP) and Thymidine di-phosphate (TDP) from the active-site of human and Thermus thermophilus thymidine mono-phosphate kinase (TMPK), referred to as hTMPK and ttTMPK, respectively. TMPK catalyses the conversion of Thymidine mono-phosphate (TMP) to TDP using ATP as phosphoryl donor in the presence of Mg2+ ion. Most of the earlier studies on this enzyme have focused on understanding substrate binding and catalysis, but the critical product release event remains elusive. Competitive binding experiments of the substrates and the products using ttTMPK apo crystals have indicated that the substrate (TMP) can replace the bound product (TDP), even in the presence of an ADP molecule. Further, the existing random accelerated molecular dynamics (RAMD) simulation program was modified to study the release of both the products simultaneously from the active site. The RAMD simulations on product-bound structures of both ttTMPK and hTMPK, revealed that while several exit patterns of the products are permissible, the sequential exit mode is the most preferred pattern for both ttTMPK and hTMPK enzymes. Additionally, the product release from the hTMPK was found to be faster and more directional as compared to ttTMPK. Structural investigation revealed that the critical changes in the residue composition in the LID-region of ttTMPK and hTMPK have an effect on the product release and can be attributed to the observed differences during product release event. Understanding of these dissimilarities is of considerable utility in designing potent inhibitors or prodrugs that can distinguish between eukaryotic and prokaryotic homologues of thymidylate kinase.

Solvation Structure and Thermodynamic Mapping (SSTMap): An Open-Source, Flexible Package for the Analysis of Water in Molecular Dynamics Trajectories.

We have developed SSTMap, a software package for mapping structural and thermodynamic water properties in molecular dynamics trajectories. The package introduces automated analysis and mapping of local measures of frustration and enhancement of water structure. The thermodynamic calculations are based on Inhomogeneous Fluid Solvation Theory (IST), which is implemented using both site-based and grid-based approaches. The package also extends the applicability of solvation analysis calculations to multiple molecular dynamics (MD) simulation programs by using existing cross-platform tools for parsing MD parameter and trajectory files. SSTMap is implemented in Python and contains both command-line tools and a Python module to facilitate flexibility in setting up calculations and for automated generation of large data sets involving analysis of multiple solutes. Output is generated in formats compatible with popular Python data science packages. This tool will be used by the molecular modeling community for computational analysis of water in problems of biophysical interest such as ligand binding and protein function.

Girsanov reweighting for path ensembles and Markov state models.

The sensitivity of molecular dynamics on changes in the potential energy function plays an important role in understanding the dynamics and function of complex molecules. We present a method to obtain path ensemble averages of a perturbed dynamics from a set of paths generated by a reference dynamics. It is based on the concept of path probability measure and the Girsanov theorem, a result from stochastic analysis to estimate a change of measure of a path ensemble. Since Markov state models (MSMs) of the molecular dynamics can be formulated as a combined phase-space and path ensemble average, the method can be extended to reweight MSMs by combining it with a reweighting of the Boltzmann distribution. We demonstrate how to efficiently implement the Girsanov reweighting in a molecular dynamics simulation program by calculating parts of the reweighting factor "on the fly" during the simulation, and we benchmark the method on test systems ranging from a two-dimensional diffusion process and an artificial many-body system to alanine dipeptide and valine dipeptide in implicit and explicit water. The method can be used to study the sensitivity of molecular dynamics on external perturbations as well as to reweight trajectories generated by enhanced sampling schemes to the original dynamics.

4H-SiC Surface Structures and Oxidation Mechanism Revealed by Using First-Principles and Classical Molecular Dynamics Simulations

The 4H-SiC(000-1)C and (0001)Si surface reconstructed structures and the oxidation processes of these surfaces are investigated using a first-principles calculation method. The most stable reconstructed 4H-SiC(000-1)C and (0001)Si surfaces have p-bonded chains. In the topmost SiC bilayer, half of Si and C atoms exchange their positions and C-C or Si-Si bonds formed densely below the surfaces. When we place a SiO2 layer on the p-bonded chain (000-1)C surface, C-C bonds are formed more densely below the interface. We simulate a sequence of O2 molecules arrivals at an interface of tridymite-phase-SiO2 and 4H-SiC(000-1)C. Dissociated O atoms at the interface tended to make bonds with Si atoms. The C-C bonds in the SiC substrate break easily and a local C surface occasionally appears. We have examined how the surface structure changes through an O2 molecule exposure by using a classical molecular dynamics simulation program and confirmed the formation of C-C bonds below the surface and the interface.

CHARMM-GUI Input Generator for NAMD, Gromacs, Amber, Openmm, and CHARMM/OpenMM Simulations using the CHARMM36 Additive Force Field

Proper treatment of nonbonded interactions is essential for the accuracy of molecular dynamics (MD) simulations, especially in studies of lipid mono- and bilayers. The use of the CHARMM36 force field (C36 FF) in different MD simulation programs can result in disagreements with published simulations performed with CHARMM due to differences in the protocols used to treat the long-range and 1-4 nonbonded interactions. In this study, we systematically test the use of the C36 lipid FF in NAMD, GROMACS, AMBER, OpenMM, and CHARMM/OpenMM. A wide range of Lennard-Jones (LJ) cutoff schemes and integrator algorithms were tested to find the optimal simulation protocol to best match lipid properties of six bilayers with varying acyl chain saturation and head groups. MD simulations of 1,2-dipalmitoyl-sn-phosphatidylcholine (DPPC) bilayers were used to obtain optimal protocols for each program. MD simulations with all programs were found to reasonably match the DPPC bilayer membrane properties, such as surface area per lipid, chain order parameters, and area compressibility modulus, obtained using the standard protocol used in CHARMM as well as from experiments. The optimal simulation protocol was then applied to the other five lipid simulations and resulted in excellent agreement between results from most simulation programs as well as with experimental data. AMBER compared least favorably with the expected membrane properties, which appears to due to its use of the hard-truncation of the LJ potential versus a force-based switching function used to smooth the LJ potential as it approaches the cutoff distance. The optimal simulation protocol for each program has been implemented in CHARMM-GUI. We expect this protocol to be applicable to the remainder of the additive C36 FF including the proteins, nucleic acids, carbohydrates, and small molecules.

Early Experiences Porting the NAMD and VMD Molecular Simulation and Analysis Software to GPU-Accelerated OpenPOWER Platforms.

All-atom molecular dynamics simulations of biomolecules provide a powerful tool for exploring the structure and dynamics of large protein complexes within realistic cellular environments. Unfortunately, such simulations are extremely demanding in terms of their computational requirements, and they present many challenges in terms of preparation, simulation methodology, and analysis and visualization of results. We describe our early experiences porting the popular molecular dynamics simulation program NAMD and the simulation preparation, analysis, and visualization tool VMD to GPU-accelerated OpenPOWER hardware platforms. We report our experiences with compiler-provided autovectorization and compare with hand-coded vector intrinsics for the POWER8 CPU. We explore the performance benefits obtained from unique POWER8 architectural features such as 8-way SMT and its value for particular molecular modeling tasks. Finally, we evaluate the performance of several GPU-accelerated molecular modeling kernels and relate them to other hardware platforms.

CHARMM-GUI Input Generator for NAMD, GROMACS, AMBER, OpenMM, and CHARMM/OpenMM Simulations Using the CHARMM36 Additive Force Field.

Proper treatment of nonbonded interactions is essential for the accuracy of molecular dynamics (MD) simulations, especially in studies of lipid bilayers. The use of the CHARMM36 force field (C36 FF) in different MD simulation programs can result in disagreements with published simulations performed with CHARMM due to differences in the protocols used to treat the long-range and 1-4 nonbonded interactions. In this study, we systematically test the use of the C36 lipid FF in NAMD, GROMACS, AMBER, OpenMM, and CHARMM/OpenMM. A wide range of Lennard-Jones (LJ) cutoff schemes and integrator algorithms were tested to find the optimal simulation protocol to best match bilayer properties of six lipids with varying acyl chain saturation and head groups. MD simulations of a 1,2-dipalmitoyl-sn-phosphatidylcholine (DPPC) bilayer were used to obtain the optimal protocol for each program. MD simulations with all programs were found to reasonably match the DPPC bilayer properties (surface area per lipid, chain order parameters, and area compressibility modulus) obtained using the standard protocol used in CHARMM as well as from experiments. The optimal simulation protocol was then applied to the other five lipid simulations and resulted in excellent agreement between results from most simulation programs as well as with experimental data. AMBER compared least favorably with the expected membrane properties, which appears to be due to its use of the hard-truncation in the LJ potential versus a force-based switching function used to smooth the LJ potential as it approaches the cutoff distance. The optimal simulation protocol for each program has been implemented in CHARMM-GUI. This protocol is expected to be applicable to the remainder of the additive C36 FF including the proteins, nucleic acids, carbohydrates, and small molecules.

Free-Energy Landscapes of the Translocation of a Substrate in Four Proteasome - Activator Complexes Analyzed using Molecular Dynamics Simulations

Proteasome is involved in the degradation of proteins. Proteasome activators bind to the proteasome core particle (CP) and facilitate opening a gate of the CP, where Tyr8 and Asp9 in the N-termini tails of the CP form the ordered open gate. In a double mutant (Tyr8Gly/Asp9Gly), the N-termini tails are disordered and the stabilized open-gate conformation cannot be formed. To understand how proteasome is responsible for the efficient proteolysis of the substrate, four different molecular dynamics (MD) simulations were carried out: ordered- and disordered-gate models of the CP complexed with an ATP-independent PA26 and ordered- and disordered-gate models of the CP complexed with an ATP-dependent PAN-like activator. The free-energies of the translocation of a poly-peptide substrate moving through the gate were estimated using an MD simulation program called SCUBA developed by the Japan Atomic Energy Agency (JAEA). In the ordered-gate models, the substrate in the activator was more stable than that in the CP. The conformational entropy of the N-termini tails of the CP was larger when the substrate was in the activator than in the CP. In the disordered-gate models, the substrate in the activator was more destabilized than in the ordered-gate models. The mutated N-termini tails became randomized and their increased conformational entropy could no longer increase further even when the substrate was in the activator, meaning the randomized N-termini tails had lost the ability to stabilize the substrate in the activator. Thus, it was concluded that the dynamics of the N-termini tails entropically play a key role in the translocation of the substrate.

分子動力学プログラムMXDORTO用可視化ソフトウェアMDVISの開発

分子動力学プログラムMXDORTO用可視化ソフトウェアMDVISの開発