GPU-accelerated Molecular Dynamics Research Articles

Sire: An interoperability engine for prototyping algorithms and exchanging information between molecular simulation programs.

Sire is a Python/C++ library that is used both to prototype new algorithms and as an interoperability engine for exchanging information between molecular simulation programs. It provides a collection of file parsers and information converters that together make it easier to combine and leverage the functionality of many other programs and libraries. This empowers researchers to use sire to write a single script that can, for example, load a molecule from a PDBx/mmCIF file via Gemmi, perform SMARTS searches via RDKit, parameterize molecules using BioSimSpace, run GPU-accelerated molecular dynamics via OpenMM, and then display the resulting dynamics trajectory in a NGLView Jupyter notebook 3D molecular viewer. This functionality is built on by BioSimSpace, which uses sire's molecular information engine to interconvert with programs such as GROMACS, NAMD, Amber, and AmberTools for automated molecular parameterization and the running of molecular dynamics, metadynamics, and alchemical free energy workflows. Sire comes complete with a powerful molecular information search engine, plus trajectory loading and editing, analysis, and energy evaluation engines. This, when combined with an in-built computer algebra system, gives substantial flexibility to researchers to load, search for, edit, and combine molecular information from multiple sources and use that to drive novel algorithms by combining functionality from other programs. Sire is open source (GPL3) and is available via conda and at a free Jupyter notebook server at https://try.openbiosim.org. Sire is supported by the not-for-profit OpenBioSim community interest company.

Protein Interface Regulating the Inserting Process of Imidazole Ionic Liquids into the Cell Membrane.

Ionic liquids (ILs) have shown promising potential in membrane protein extraction; however, the underlying mechanism remains unclear. Herein, we employed GPU-accelerated molecular dynamics (MD) simulations to investigate the dynamic insertion process of ILs into cell membranes containing membrane proteins. Our findings reveal that ILs spontaneously insert into the membrane, and the presence of membrane proteins significantly decelerates the rate of IL insertion into the membrane. Specifically, the relationship between the insertion rate and inserting free energy exhibits non-monotonic changes, which can be attributed to interfacial effects. The protein-water interface acts as trap for free ions and ionic clusters, while free ions preferentially insert into the membrane from the protein-lipid interface, which limits the insertion rate due to its narrowness. Thus, the insertion rate is governed by a combination of the free energy and interfacial effects. These findings provide valuable insights into the interfacial effects of protein-lipid bilayers and have implications for various biochemical-related applications.

Predicting the Activity of Boronate Inhibitors Against Metallo-beta-lactamase Enzymes

Potency of boronate inhibitors against metallo-β-lactamases (MβLs) has been found to be dependent on the electrophilicity of the boron atom. It forms a covalent bond with the oxygen atom of the catalytic OH− ion in the active site of the enzyme. The ability of the boronate inhibitor to influence the protein conformation also affects the binding potency. Molecular dynamics (MD) simulations of cyclic and non-cyclic boronate complexes with NDM-1 MβL show their higher impact on the inhibitor efficiency compared with the electrophilicity of the boron atom. Therefore, we focus on the hardware impact on the computational speedup of the GPU-accelerated MD. Using this data, we propose a comprehensive protocol for in silico prediction of the activity of boronate molecules against MβL enzymes, which includes MD simulations, combined quantum mechanics / molecular mechanics (QM/MM) computations and molecular dynamics simulations with the QM/MM potentials (QM/MM MD).

Thermodynamical Origin of Nonmonotonic Inserting Behavior of Imidazole Ionic Liquids into the Lipid Bilayer.

The GPU-accelerated molecular dynamics simulations are performed to explore the dynamical inserting process of ionic liquids (ILs) into the lipid bilayer. We found that the free ions and clusters coexist in the system, but only the cation can insert into the lipid bilayer. In specific, after a microsecond-scale simulation (up to 1.16 μs), the inserting rate increases first and then decreases nonmonotonic as side chain of cation (nchain) elongates, peaking at nchain = 10. However, the inserting free energy decreases with nchain, indicating the inserting process is easier for the larger nchain. Such contrary originates from the formation of cluster, where the cluster dissociating energy shows that only cluster for nchain ≤ 10 can dissociate spontaneously. Hence, the inserting rate is determined by the balance between nchain and cluster stability. These quantitative competition mechanisms shed light to the rational design of the biocompatible ILs toward their applications in the biochemical-related fields.

GPU-accelerated molecular dynamics: State-of-art software performance and porting from Nvidia CUDA to AMD HIP

Classical molecular dynamics (MD) calculations represent a significant part of the utilization time of high-performance computing systems. As usual, the efficiency of such calculations is based on an interplay of software and hardware that are nowadays moving to hybrid GPU-based technologies. Several well-developed open-source MD codes focused on GPUs differ both in their data management capabilities and in performance. In this work, we analyze the performance of LAMMPS, GROMACS and OpenMM MD packages with different GPU backends on Nvidia Volta and AMD Vega20 GPUs. We consider the efficiency of solving two identical MD models (generic for material science and biomolecular studies) using different software and hardware combinations. We describe our experience in porting the CUDA backend of LAMMPS to ROCm HIP that shows considerable benefits for AMD GPUs comparatively to the OpenCL backend.

DdcMD: A fully GPU-accelerated molecular dynamics program for the Martini force field.

We have implemented the Martini force field within Lawrence Livermore National Laboratory's molecular dynamics program, ddcMD. The program is extended to a heterogeneous programming model so that it can exploit graphics processing unit (GPU) accelerators. In addition to the Martini force field being ported to the GPU, the entire integration step, including thermostat, barostat, and constraint solver, is ported as well, which speeds up the simulations to 278-fold using one GPU vs one central processing unit (CPU) core. A benchmark study is performed with several test cases, comparing ddcMD and GROMACS Martini simulations. The average performance of ddcMD for a protein-lipid simulation system of 136k particles achieves 1.04 µs/day on one NVIDIA V100 GPU and aggregates 6.19 µs/day on one Summit node with six GPUs. The GPU implementation in ddcMD offloads all computations to the GPU and only requires one CPU core per simulation to manage the inputs and outputs, freeing up remaining CPU resources on the compute node for alternative tasks often required in complex simulation campaigns. The ddcMD code has been made open source and is available on GitHub at https://github.com/LLNL/ddcMD.

Computational Estimation of Microsecond to Second Atomistic Folding Times.

Despite the development of massively parallel computing hardware including inexpensive graphics processing units (GPUs), it has remained infeasible to simulate the folding of atomistic proteins at room temperature using conventional molecular dynamics (MD) beyond the microsecond scale. Here, we report the folding of atomistic, implicitly solvated protein systems with folding times τ ranging from ∼10 μs to ∼100 ms using the weighted ensemble (WE) strategy in combination with GPU computing. Starting from an initial structure or set of structures, WE organizes an ensemble of GPU-accelerated MD trajectory segments via intermittent pruning and replication events to generate statistically unbiased estimates of rate constants for rare events such as folding; no biasing forces are used. Although the variance among atomistic WE folding runs is significant, multiple independent runs are used to reduce and quantify statistical uncertainty. Folding times are estimated directly from WE probability flux and from history-augmented Markov analysis of the WE data. Three systems were examined: NTL9 at low solvent viscosity (yielding τf = 0.8-9 μs), NTL9 at water-like viscosity (τf = 0.2-2 ms), and Protein G at low viscosity (τf = 3-200 ms). In all cases, the folding time, uncertainty, and ensemble properties could be estimated from WE simulation; for Protein G, this characterization required significantly less overall computing than would be required to observe a single folding event with conventional MD simulations. Our results suggest that the use and calibration of force fields and solvent models for precise estimation of kinetic quantities is becoming feasible.

Explorations of Drug Transport by P-Glycoprotein using Molecular Dynamics Enabled by High Resolution Crystal Structures

The integral membrane efflux pump P-glycoprotein (P-gp, ABCB1) is a causative agent of multidrug resistance (MDR) to a wide variety of chemotherapeutics in human cancers. In many cancer subtypes, chemotherapy treatment can both induce and select for the overexpression of P-gp. When P-gp is overexpressed, the pump lowers the intracellular concentration of many drugs to sub-therapeutic levels, thereby conferring general MDR to the cancer, and rendering further chemotherapeutic treatment ineffective. Given the clinical importance of the enzyme, P-gp has long been a target for drug discovery and the focus of extensive mechanistic research. Until recently, the lack of known human P-gp structures has hindered efforts to understand the unusually polyspecific substrate profile and complex transport mechanism of the enzyme. Prior research by our laboratory has demonstrated the possibility of multiple transport pathways through the drug binding domains (DBDs) of the transporter, and of variable substrate binding sites using targeted molecular dynamics (TMD) simulations on homology models of human P-gp. Following the recent release of the cryo-EM structure of human P-gp, we performed novel simulations of the solved structure using four known, low energy conformations of P-gp homologues in various states along a putative ABC transporter catalytic efflux cycle. These simulations were done at long timescales using a combination of free and directed GPU-accelerated molecular dynamics with NAMD. These analyses included the new human P-gp structure and demonstrated the pumping of known P-gp substrates across the cell membrane. Parameterizations and simulations of molecules that are either suspected transport substrates or putative transport inhibitors of P-gp were performed. These computational studies support earlier observations, and explore the hypothetical mechanisms of P-gp-mediated substrate efflux and of P-gp transport inhibition. This work was supported by NIH NIGMS [R15GM094771-02].

GPU-Accelerated Molecular Dynamics and Free Energy Methods in Amber18: Performance Enhancements and New Features.

We report progress in graphics processing unit (GPU)-accelerated molecular dynamics and free energy methods in Amber18. Of particular interest is the development of alchemical free energy algorithms, including free energy perturbation and thermodynamic integration methods with support for nonlinear soft-core potential and parameter interpolation transformation pathways. These methods can be used in conjunction with enhanced sampling techniques such as replica exchange, constant-pH molecular dynamics, and new 12-6-4 potentials for metal ions. Additional performance enhancements have been made that enable appreciable speed-up on GPUs relative to the previous software release.

Assembled Structures of Perfluorosulfonic Acid Ionomers Investigated by Anisotropic Modeling and Simulations.

Nafion, a classic of perfluorosulfonic acid ionomers, has broad applications in proton conduction, attributed from the unique structures. However, a satisfactory structure model from theoretical calculation and simulation that can match with the well-known experimental observations is still absent. We performed GPU-accelerated molecular dynamics simulations to investigate the assembled structures of Nafion at different water contents based on an anisotropic coarse-grained model equipped with Gay-Berne potential. Accurate parameters for the coarse-grained model are collected by matching energy profiles based on density functional theory calculations. The results show that the hydrophilic phase in Nafion assemblies undergoes a crossover from isolated spherical clusters to interconnected cluster/channel networks with the increase of water content. We found the crystalline domains in polymer matrix and they are suppressed at elevated water content. These microphase-separated structures achieve quantitative agreement with existing experimental observations, including morphologies from electron microscopy and intensity profiles from scattering experiments. This work suggests that accurate consideration of the anisotropy is a key to reveal the formation of unique assembled structures of perfluorosulfonic acid ionomers at different water contents.

Atomistic Simulations of Lithium Ion Mobility in Battery Electrolytes

Li-ion battery electrolytes are complex solutions composed of solvents, functional additives, and lithium salts. Lithium ion mobility is a central property affecting battery performance. The chemical and physical factors impacting charge transport in complex solutions are not fully understood. Molecular dynamics (MD) simulation is a powerful technique that can provide insight into the critical details impacting ion transport in complex electrolyte solutions. Although several studies have been reported in the literature examining lithium diffusion using MD, limited computational power has restricted simulation time scales to only tens of nanoseconds. Additionally, the simulation results have been hampered due to lack of accurate force-field parameters. In this work, the enhanced predictive reliability for battery electrolytes provided by GPU-accelerated MD using the highly efficient Desmond MD engine [1], along with the quantitative OPLS3 force-field [2] is demonstrated. Using GPU-Desmond/OPLS3, hundreds of simulations can be performed in a production environment with simulation duration in excess of 1 µs. We discuss the impact that simulation duration and force-field has on the accuracy of Li-ion electrolyte simulations. [1] K.J. Bowers, E. Chow, H. Xu, R.O. Dror, M.P. Eastwood, B.A. Gregersen, J.L. Klepeis, I. Kolossvary, M.A. Moraes, F.D. Sacerdoti, J.K.Salmon, Y. Shan and D.E. Shaw, Proceedings of ACM/IEEE Conference on Supercomputing (SC06), Tampa, Florida November 11-17, 2006. [2] E. Harder, W. Damm, J. Maple, C. Wu, M. Reboul, J.Y. Xiang, L. Wang, D. Lupyan, M.K. Dahlgren, J.L. Knight, J.W. Kaus, D.S. Cerutti, G. Krilov, W.L. Jorgensen, R. Abel, R.A. Friesner, J. Chem. Theory Comput., 2016, 12 (1), pp 281; DOI: 10.1021/acs.jctc.5b00864

GPU-Accelerated Molecular Dynamics Simulation to Study Liquid Crystal Phase Transition Using Coarse-Grained Gay-Berne Anisotropic Potential.

Gay-Berne (GB) potential is regarded as an accurate model in the simulation of anisotropic particles, especially for liquid crystal (LC) mesogens. However, its computational complexity leads to an extremely time-consuming process for large systems. Here, we developed a GPU-accelerated molecular dynamics (MD) simulation with coarse-grained GB potential implemented in GALAMOST package to investigate the LC phase transitions for mesogens in small molecules, main-chain or side-chain polymers. For identical mesogens in three different molecules, on cooling from fully isotropic melts, the small molecules form a single-domain smectic-B phase, while the main-chain LC polymers prefer a single-domain nematic phase as a result of connective restraints in neighboring mesogens. The phase transition of side-chain LC polymers undergoes a two-step process: nucleation of nematic islands and formation of multi-domain nematic texture. The particular behavior originates in the fact that the rotational orientation of the mesogenes is hindered by the polymer backbones. Both the global distribution and the local orientation of mesogens are critical for the phase transition of anisotropic particles. Furthermore, compared with the MD simulation in LAMMPS, our GPU-accelerated code is about 4 times faster than the GPU version of LAMMPS and at least 200 times faster than the CPU version of LAMMPS. This study clearly shows that GPU-accelerated MD simulation with GB potential in GALAMOST can efficiently handle systems with anisotropic particles and interactions, and accurately explore phase differences originated from molecular structures.

MyPresto/omegagene: a GPU-accelerated molecular dynamics simulator tailored for enhanced conformational sampling methods with a non-Ewald electrostatic scheme.

Molecular dynamics (MD) is a promising computational approach to investigate dynamical behavior of molecular systems at the atomic level. Here, we present a new MD simulation engine named “myPresto/omegagene” that is tailored for enhanced conformational sampling methods with a non-Ewald electrostatic potential scheme. Our enhanced conformational sampling methods, e.g., the virtual-system-coupled multi-canonical MD (V-McMD) method, replace a multi-process parallelized run with multiple independent runs to avoid inter-node communication overhead. In addition, adopting the non-Ewald-based zero-multipole summation method (ZMM) makes it possible to eliminate the Fourier space calculations altogether. The combination of these state-of-the-art techniques realizes efficient and accurate calculations of the conformational ensemble at an equilibrium state. By taking these advantages, myPresto/omegagene is specialized for the single process execution with Graphics Processing Unit (GPU). We performed benchmark simulations for the 20-mer peptide, Trp-cage, with explicit solvent. One of the most thermodynamically stable conformations generated by the V-McMD simulation is very similar to an experimentally solved native conformation. Furthermore, the computation speed is four-times faster than that of our previous simulation engine, myPresto/psygene-G. The new simulator, myPresto/omegagene, is freely available at the following URLs: http://www.protein.osaka-u.ac.jp/rcsfp/pi/omegagene/ and http://presto.protein.osaka-u.ac.jp/myPresto4/.

Computational grounding of point mutagenesis to enhance the thermostability of human peroxiredoxin 6

In this study an algorithm for increasing the thermostability of the globular human peroxyredoxin 6 antioxidant enzyme is described. A recombinant form of this protein was produced based on the analysis of the amino-acid composition of this enzyme in rats (Rattus rattus). The original concept of increasing the thermostability of small globular proteins using alternative hydrogen bonding of the peripheral acidic and basic amino-acid residues was substantially improved by GPU-accelerated molecular-dynamics simulations.

Molecular Design of Amorphous Porous Organic Cages for Enhanced Gas Storage

Porous molecular solids are garnering increasing attention with examples of high surface areas and applications in molecular separations. Recently, amorphous networks of molecular cages have shown increased porosity with respect to their crystalline counterparts. However, the structures of amorphous materials cannot be precisely elucidated by X-ray diffraction techniques, thus molecular simulations are vital to understanding their pore structures and the origin of porosity. Here, we use GPU-accelerated molecular dynamics simulations as an efficient methodology to construct representative amorphous network structures. We employ Voronoi network analysis of amorphous networks of seven previously reported cage molecules to provide insight into structure–property relationships. Accordingly, we apply this understanding to delineate synthetic design features that give rise to highly porous analogues of chemically robust cages constructed from carbon–carbon bonds.

Big data in drug design

Big data collection in the pharmaceutical research industry has four sources, high-throughput scientific experiments, high-performance computations, automated information acquisition and office automation, and scientific publications and patents. Big data is the product and the promoter of high-performance scientific experiments, therefore the technology for mining big data is the key to future drug discovery. However, big data brings greater challenges, such as, storage, retrieval, curation and quality assurance, sharing/transfer, analysis, visualization, modeling and computing complexities. This review outlines the current progress of processing big data in the drug design field. These problems may be resolved by adopting cloud computing and high-performance computing technologies, and parallelizing existing chemoinformatics and bioinformatics programs. Machine-learning approaches involving Bayesian learning methods and other methods, such as support vector machine and recursive petitioning, can be used for big data mining. Recent progress includes parallelized and GPU-accelerated molecular dynamics simulation technology, enhanced molecular docking technology, new parallelized algorithms for shape-based virtual screening, free-energy landscape calculations, and machine-learning algorithms for big chemical structural data mining. Big data from drug discoveries will increase, so conventional drug design software and methods need to be upgraded. This is a long-term project and we highlight the tasks that need to be accomplished to meet this goal.

π–π Stacking mediated drug–drug interactions in human CYP2E1

Because of having many low molecular mass substrates, CYP2E1 is of particular interests to the pharmaceutical industry. Many evidences showed that this enzyme can adopt multiple substrates to significantly reduce the oxidation rate of the substrates. The detailed mechanism for this observation is still unclear. In the current study, we employed GPU-accelerated molecular dynamics simulations to study the multiple-binding mode of human CYP2E1, with an aim of offering a mechanistic explanation for the unexplained multiple-substrate binding. Our results showed that Thr303 and Phe478 were key factors for the substrate recognition and multiple-substrate binding. The former can form a significant hydrogen bond to recognize and position the substrate in the productive binding orientation in the active site. The latter acted as a mediator for the substrate communications via π-π stacking interactions. In the multiple-binding mode, the aforementioned π-π stacking interactions formed by the aromatic rings of both substrates and Phe478 drove the first substrate far away from the catalytic center, orienting in an additional binding position and going against the substrate metabolism. All these findings could give atomic insights into the detailed mechanism for the multiple-substrate binding in human CYP2E1, providing useful information for the drug metabolism mechanism and personalized use of clinical drugs.

Efficient GPU-accelerated molecular dynamics simulation of solid covalent crystals

Different from previous molecular dynamics (MD) simulation with pair potentials and many-body potentials, an efficient and highly scalable algorithm for GPU-accelerated MD simulation of solid covalent crystals is described in detail in this paper using sophisticated many-body potentials, such as Tersoff potentials for silicon crystals. The algorithm has effectively taken advantage of the reordering and sorting of atoms and the hierarchical memory of a GPU. The final results indicate that, about 30.5% of the peak performance of a single GPU can be achieved with a speedup of about 650 over a contemporary CPU core, and more than 15 million atoms can be processed by a single GPU with a speed of around 2 ns/day. Furthermore, the proposed algorithm is scalable and transferable, which can be applied to other many-body interactions and related large-scale parallel computation.

Physics Based Protein Structure Refinement

Accurate protein structure predictions are important for a number of purposes ranging from computational drug design, understanding experimental data and designing new experiments, to the emerging technique of de novo phasing in crystallography. The process of protein structure refinement occurs at the end of the structure prediction pipeline. The goal is to take an approximately correct starting model and further refine the details to produce a more accurate prediction. We have developed a physics-based approach to refinement that combines Hamiltonian exchange molecular dynamics with bioinformatics-derived restraints. The use of restraints dramatically reduces the volume of phase space that must be sampled and makes the procedure practical on small to medium size proteins. We combine this approach with GPU-accelerated molecular dynamics, new implicit solvent models, and recent improvements in force fields. The initial results of this protocol are encouraging and we are able to successfully refine several difficult CASP9 refinement targets.View Large Image | View Hi-Res Image | Download PowerPoint Slide

Self-assembly of coarse-grained ionic surfactants accelerated by graphics processing units

Due to the relatively long time scales inherent to ionic surfactant self-assembly (>μs), an aggressive computational approach is needed to obtain converged data on micellar solutions. This work presents a study of micellization using a coarse-grained (CG) model of aqueous ionic surfactants in replicated molecular dynamics (MD) simulations run on graphics processing unit hardware. The performance of our implementation of the CG model with electrostatics into the HOOMD-Blue GPU-accelerated MD software package is comparable to that of a modest sized cluster running a highly optimized parallel CPU code. From 0.36 ms of cumulative trajectory data, we are able to predict equilibrium thermodynamic and morphological properties of ionic surfactant micellar solutions. Estimating the critical micelle concentrations (CMC) from the free monomer (ρ1) and premicellar concentrations obtained from simulations of sodium hexyl sulfate (S6S, CMC of 460 ± 6 mM) at high (1 M) concentration, a value in good agreement with experimental results is obtained; however, the same method applied to simulations of sodium nonyl sulfate (S9S, ρ1 of 2.4 ± 0.01 mM) and sodium dodecyl sulfate (SDS, ρ1 of 0.02 ± 0.01 mM) at the same total concentration systematically underestimates the CMCs. An alternative method for calculating the CMC is presented, where the free monomer concentration computed from high concentration CG-MD data is used as the input to a simple theoretical model which can be used to extrapolate to a more accurate prediction of the CMC. Better agreement between the empirical and predicted CMC is obtained from this theory for S9S (28.7 ± 0.3 mM) and SDS (3.32 ± 0.04 mM), though the CMC for S6S is slightly underestimated (304 ± 3 mM). We also present statistically converged morphological data, including aggregation number distributions and the principal components of the gyration tensor. This data suggest a transition from spherical micelles to rod-like at a specific aggregation number, which increases with increasing hydrocarbon length.