A GPU solvent–solvent interaction calculation accelerator for biomolecular simulations using the GROMOS software

During the past few years, graphics processing units (GPUs) have become extremely popular in the high performance computing community. In this study, we present an implementation of an acceleration engine for the solvent–solvent interaction evaluation of molecular dynamics simulations. By careful optimization of the algorithm speed‐ups up to a factor of 54 (single‐precision GPU vs. double‐precision CPU) could be achieved. The accuracy of the single‐precision GPU implementation is carefully investigated and does not influence structural, thermodynamic, and dynamic quantities. Therefore, the implementation enables users of the GROMOS software for biomolecular simulation to run the solvent–solvent interaction evaluation on a GPU, and thus, to speed‐up their simulations by a factor 6–9. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010     

Molecular dynamics with the united-residue model of polypeptide chains. II. Langevin and Berendsen-bath dynamics and tests on model alpha-helical systems

The implementation of molecular dynamics (MD) with our physics-based protein united-residue (UNRES) force field, described in the accompanying paper, was extended to Langevin dynamics. The equations of motion are integrated by using a simplified stochastic velocity Verlet algorithm. To compare the results to those with all-atom simulations with implicit solvent in which no explicit stochastic and friction forces are present, we alternatively introduced the Berendsen thermostat. Test simulations on the Ala(10) polypeptide demonstrated that the average kinetic energy is stable with about a 5 fs time step. To determine the correspondence between the UNRES time step and the time step of all-atom molecular dynamics, all-atom simulations with the AMBER 99 force field and explicit solvent and also with implicit solvent taken into account within the framework of the generalized Born/surface area (GBSA) model were carried out on the unblocked Ala(10) polypeptide. We found that the UNRES time scale is 4 times longer than that of all-atom MD simulations because the degrees of freedom corresponding to the fastest motions in UNRES are averaged out. When the reduction of the computational cost for evaluation of the UNRES energy function is also taken into account, UNRES (with hydration included implicitly in the side chain-side chain interaction potential) offers about at least a 4000-fold speed up of computations relative to all-atom simulations with explicit solvent and at least a 65-fold speed up relative to all-atom simulations with implicit solvent. To carry out an initial full-blown test of the UNRES/MD approach, we ran Berendsen-bath and Langevin dynamics simulations of the 46-residue B-domain of staphylococcal protein A. We were able to determine the folding temperature at which all trajectories converged to nativelike structures with both approaches. For comparison, we carried out ab initio folding simulations of this protein at the AMBER 99/GBSA level. The average CPU time for folding protein A by UNRES molecular dynamics was 30 min with a single Alpha processor, compared to about 152 h for all-atom simulations with implicit solvent. It can be concluded that the UNRES/MD approach will enable us to carry out microsecond and, possibly, millisecond simulations of protein folding and, consequently, of the folding process of proteins in real time.     

Molecular dynamics simulations of aqueous ions at the liquid–vapor interface accelerated using graphics processors

Molecular dynamics (MD) simulations are a vital tool in chemical research, as they are able to provide an atomistic view of chemical systems and processes that is not obtainable through experiment. However, large‐scale MD simulations require access to multicore clusters or supercomputers that are not always available to all researchers. Recently, scientists have returned to exploring the power of graphics processing units (GPUs) for various applications, such as MD, enabled by the recent advances in hardware and integrated programming interfaces such as NVIDIA's CUDA platform. One area of particular interest within the context of chemical applications is that of aqueous interfaces, the salt solutions of which have found application as model systems for studying atmospheric process as well as physical behaviors such as the Hoffmeister effect. Here, we present results of GPU‐accelerated simulations of the liquid–vapor interface of aqueous sodium iodide solutions. Analysis of various properties, such as density and surface tension, demonstrates that our model is consistent with previous studies of similar systems. In particular, we find that the current combination of water and ion force fields coupled with the ability to simulate surfaces of differing area enabled by GPU hardware is able to reproduce the experimental trend of increasing salt solution surface tension relative to pure water. In terms of performance, our GPU implementation performs equivalent to CHARMM running on 21 CPUs. Finally, we address possible issues with the accuracy of MD simulaions caused by nonstandard single‐precision arithmetic implemented on current GPUs. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2011     

Accelerating molecular dynamic simulation on graphics processing units

We describe a complete implementation of all‐atom protein molecular dynamics running entirely on a graphics processing unit (GPU), including all standard force field terms, integration, constraints, and implicit solvent. We discuss the design of our algorithms and important optimizations needed to fully take advantage of a GPU. We evaluate its performance, and show that it can be more than 700 times faster than a conventional implementation running on a single CPU core. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2009     

GPU accelerated numerical simulations of viscoelastic phase separation model

We introduce a complete implementation of viscoelastic model for numerical simulations of the phase separation kinetics in dynamic asymmetry systems such as polymer blends and polymer solutions on a graphics processing unit (GPU) by CUDA language and discuss algorithms and optimizations in details. From studies of a polymer solution, we show that the GPU-based implementation can predict correctly the accepted results and provide about 190 times speedup over a single central processing unit (CPU). Further accuracy analysis demonstrates that both the single and the double precision calculations on the GPU are sufficient to produce high-quality results in numerical simulations of viscoelastic model. Therefore, the GPU-based viscoelastic model is very promising for studying many phase separation processes of experimental and theoretical interests that often take place on the large length and time scales and are not easily addressed by a conventional implementation running on a single CPU.     

Potential of mean force of water–proton bath and molecular dynamic simulation of proteins at constant pH

An advanced implicit solvent model of water–proton bath for protein simulations at constant pH is presented. The implicit water–proton bath model approximates the potential of mean force of a protein in water solvent in a presence of hydrogen ions. Accurate and fast computational implementation of the implicit water–proton bath model is developed using the continuum electrostatic Poisson equation model for calculation of ionization equilibrium and the corrected MSR6 generalized Born model for calculation of the electrostatic atom–atom interactions and forces. Molecular dynamics (MD) method for protein simulation in the potential of mean force of water–proton bath is developed and tested on three proteins. The model allows to run MD simulations of proteins at constant pH, to calculate pH‐dependent properties and free energies of protein conformations. The obtained results indicate that the developed implicit model of water–proton bath provides an efficient way to study thermodynamics of biomolecular systems as a function of pH, pH‐dependent ionization‐conformation coupling, and proton transfer events. © 2012 Wiley Periodicals, Inc.     

Role of ions on structure and stability of a synthetic gramicidin ion channel in solution. A molecular dynamics study

We performed a molecular dynamics (MD) simulation to the investigate structure and stability of a synthetic gramicidin-like peptide in solution with and without ions. The starting structures of the MD simulations were taken from two recently solved NMR structures of this peptide in isotropic solution, which forms stable monomers or dimers in the presence or absence of ions, respectively. The monomeric structure is channel-like and is assumed to be stabilized by the presence of two Cs(+) ions bound in the channel, each one close to one channel entrance. In our MD simulations, we observed how the Cs(+) ions bind in the channel formed by the monomeric gramicidin-like peptide using implicit solvent and explicit ions with a concentration of 2 M. MD simulations were performed with and without explicit ions but with an implicit solvent model defined by the generalized Born approximation, which was used to mimic the dielectric properties of the solvent and to speed up the computations.     

Efficient implementation of constant pH molecular dynamics on modern graphics processors

The treatment of pH sensitive ionization states for titratable residues in proteins is often omitted from molecular dynamics (MD) simulations. While static charge models can answer many questions regarding protein conformational equilibrium and protein–ligand interactions, pH‐sensitive phenomena such as acid‐activated chaperones and amyloidogenic protein aggregation are inaccessible to such models. Constant pH molecular dynamics (CPHMD) coupled with the Generalized Born with a Simple sWitching function (GBSW) implicit solvent model provide an accurate framework for simulating pH sensitive processes in biological systems. Although this combination has demonstrated success in predicting pK_a values of protein structures, and in exploring dynamics of ionizable side‐chains, its speed has been an impediment to routine application. The recent availability of low‐cost graphics processing unit (GPU) chipsets with thousands of processing cores, together with the implementation of the accurate GBSW implicit solvent model on those chipsets (Arthur and Brooks, J. Comput. Chem. 2016, 37, 927), provide an opportunity to improve the speed of CPHMD and ionization modeling greatly. Here, we present a first implementation of GPU‐enabled CPHMD within the CHARMM‐OpenMM simulation package interface. Depending on the system size and nonbonded force cutoff parameters, we find speed increases of between one and three orders of magnitude. Additionally, the algorithm scales better with system size than the CPU‐based algorithm, thus allowing for larger systems to be modeled in a cost effective manner. We anticipate that the improved performance of this methodology will open the door for broad‐spread application of CPHMD in its modeling pH‐mediated biological processes. © 2016 Wiley Periodicals, Inc.     

Molecular dynamics simulation of complex multiphase flow on a computer cluster with GPUs

Compute Unified Device Architecture (CUDA) was used to design and implement molecular dynamics (MD) simulations on graphics processing units (GPU). With an NVIDIA Tesla C870, a 20–60 fold speedup over that of one core of the Intel Xeon 5430 CPU was achieved, reaching up to 150 Gflops. MD simulation of cavity flow and particle-bubble interaction in liquid was implemented on multiple GPUs using a message passing interface (MPI). Up to 200 GPUs were tested on a special network topology, which achieves good scalability. The capability of GPU clusters for large-scale molecular dynamics simulation of meso-scale flow behavior was, therefore, uncovered. Supported by the National Natural Science Foundation of China (Grant Nos. 20336040, 20221603 and 20490201), and the Chinese Academy of Sciences (Grant No. Kgcxz-yw-124)     

Fast and flexible gpu accelerated binding free energy calculations within the amber molecular dynamics package

Alchemical free energy (AFE) calculations based on molecular dynamics (MD) simulations are key tools in both improving our understanding of a wide variety of biological processes and accelerating the design and optimization of therapeutics for numerous diseases. Computing power and theory have, however, long been insufficient to enable AFE calculations to be routinely applied in early stage drug discovery. One of the major difficulties in performing AFE calculations is the length of time required for calculations to converge to an ensemble average. CPU implementations of MD‐based free energy algorithms can effectively only reach tens of nanoseconds per day for systems on the order of 50,000 atoms, even running on massively parallel supercomputers. Therefore, converged free energy calculations on large numbers of potential lead compounds are often untenable, preventing researchers from gaining crucial insight into molecular recognition, potential druggability and other crucial areas of interest. Graphics Processing Units (GPUs) can help address this. We present here a seamless GPU implementation, within the PMEMD module of the AMBER molecular dynamics package, of thermodynamic integration (TI) capable of reaching speeds of >140 ns/day for a 44,907‐atom system, with accuracy equivalent to the existing CPU implementation in AMBER. The implementation described here is currently part of the AMBER 18 beta code and will be an integral part of the upcoming version 18 release of AMBER. © 2018 Wiley Periodicals, Inc.     

A high‐performance parallel‐generalized born implementation enabled by tabulated interaction rescaling

Implicit solvent representations, in general, and generalized Born models, in particular, provide an attractive way to reduce the number of interactions and degrees of freedom in a system. The instantaneous relaxation of the dielectric shielding provided by an implicit solvent model can be extremely efficient for high‐throughput and Monte Carlo studies, and a reduced system size can also remove a lot of statistical noise. Despite these advantages, it has been difficult for generalized Born implementations to significantly outperform optimized explicit‐water simulations due to more complex functional forms and the two extra interaction stages necessary to calculate Born radii and the derivative chain rule terms contributing to the force. Here, we present a method that uses a rescaling transformation to make the standard generalized Born expression a function of a single variable, which enables an efficient tabulated implementation on any modern CPU hardware. The total performance is within a factor 2 of simulations in vacuo. The algorithm has been implemented in Gromacs, including single‐instruction multiple‐data acceleration, for three different Born radius models and corresponding chain rule terms. We have also adapted the model to work with the virtual interaction sites commonly used for hydrogens to enable long‐time steps, which makes it possible to achieve a simulation performance of 0.86 μs/day for BBA5 with 1‐nm cutoff on a single quad‐core desktop processor. Finally, we have also implemented a set of streaming kernels without neighborlists to accelerate the non‐cutoff setup occasionally used for implicit solvent simulations of small systems. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010     

Benchmarking implicit solvent folding simulations of the amyloid beta(10-35) fragment

A pathogenetic feature of Alzhemier disease is the aggregation of monomeric beta-amyloid proteins (Abeta) to form oligomers. Usually these oligomers of long peptides aggregate on time scales of microseconds or longer, making computational studies using atomistic molecular dynamics models prohibitively expensive and making it essential to develop computational models that are cheaper and at the same time faithful to physical features of the process. We benchmark the ability of our implicit solvent model to describe equilibrium and dynamic properties of monomeric Abeta(10-35) using all-atom Langevin dynamics (LD) simulations, since Alphabeta(10-35) is the only fragment whose monomeric properties have been measured. The accuracy of the implicit solvent model is tested by comparing its predictions with experiment and with those from a new explicit water MD simulation, (performed using CHARMM and the TIP3P water model) which is approximately 200 times slower than the implicit water simulations. The dependence on force field is investigated by running multiple trajectories for Alphabeta(10-35) using the CHARMM, OPLS-aal, and GS-AMBER94 force fields, whereas the convergence to equilibrium is tested for each force field by beginning separate trajectories from the native NMR structure, a completely stretched structure, and from unfolded initial structures. The NMR order parameter, S2, is computed for each trajectory and is compared with experimental data to assess the best choice for treating aggregates of Alphabeta. The computed order parameters vary significantly with force field. Explicit and implicit solvent simulations using the CHARMM force fields display excellent agreement with each other and once again support the accuracy of the implicit solvent model. Alphabeta(10-35) exhibits great flexibility, consistent with experiment data for the monomer in solution, while maintaining a general strand-loop-strand motif with a solvent-exposed hydrophobic patch that is believed to be important for aggregation. Finally, equilibration of the peptide structure requires an implicit solvent LD simulation as long as 30 ns.     

Convergence and sampling efficiency in replica exchange simulations of peptide folding in explicit solvent

Replica exchange methods (REMs) are increasingly used to improve sampling in molecular dynamics (MD) simulations of biomolecular systems. However, despite having been shown to be very effective on model systems, the application of REM in complex systems such as for the simulation of protein and peptide folding in explicit solvent has not been objectively tested in detail. Here we present a comparison of conventional MD and temperature replica exchange MD (T-REMD) simulations of a beta-heptapeptide in explicit solvent. This system has previously been shown to undergo reversible folding on the time scales accessible to MD simulation and thus allows a direct one-to-one comparison of efficiency. The primary properties compared are the free energy of folding and the relative populations of different conformers as a function of temperature. It is found that to achieve a similar degree of precision T-REMD simulations starting from a random set of initial configurations were approximately an order of magnitude more computationally efficient than a single 800 ns conventional MD simulation for this system at the lowest temperature investigated (275 K). However, whereas it was found that T-REMD simulations are more than four times more efficient than multiple independent MD simulations at one temperature (300 K) the actual increase in conformation sampling was only twofold. The overall gain in efficiency using REMD resulted primarily from the ordering of different conformational states over temperature, as opposed to a large increase of conformational sampling. It is also shown that in this system exchanges are accepted primarily based on (random) fluctuations within the solvent and are not strongly correlated with the instantaneous peptide conformation raising questions in regard to the efficiency of T-REMD in larger systems.     

New faster CHARMM molecular dynamics engine

We introduce a new faster molecular dynamics (MD) engine into the CHARMM software package. The new MD engine is faster both in serial (i.e., single CPU core) and parallel execution. Serial performance is approximately two times higher than in the previous version of CHARMM. The newly programmed parallelization method allows the MD engine to parallelize up to hundreds of CPU cores. © 2013 Wiley Periodicals, Inc.     

Best bang for your buck: GPU nodes for GROMACS biomolecular simulations

The molecular dynamics simulation package GROMACS runs efficiently on a wide variety of hardware from commodity workstations to high performance computing clusters. Hardware features are well‐exploited with a combination of single instruction multiple data, multithreading, and message passing interface (MPI)‐based single program multiple data/multiple program multiple data parallelism while graphics processing units (GPUs) can be used as accelerators to compute interactions off‐loaded from the CPU. Here, we evaluate which hardware produces trajectories with GROMACS 4.6 or 5.0 in the most economical way. We have assembled and benchmarked compute nodes with various CPU/GPU combinations to identify optimal compositions in terms of raw trajectory production rate, performance‐to‐price ratio, energy efficiency, and several other criteria. Although hardware prices are naturally subject to trends and fluctuations, general tendencies are clearly visible. Adding any type of GPU significantly boosts a node's simulation performance. For inexpensive consumer‐class GPUs this improvement equally reflects in the performance‐to‐price ratio. Although memory issues in consumer‐class GPUs could pass unnoticed as these cards do not support error checking and correction memory, unreliable GPUs can be sorted out with memory checking tools. Apart from the obvious determinants for cost‐efficiency like hardware expenses and raw performance, the energy consumption of a node is a major cost factor. Over the typical hardware lifetime until replacement of a few years, the costs for electrical power and cooling can become larger than the costs of the hardware itself. Taking that into account, nodes with a well‐balanced ratio of CPU and consumer‐class GPU resources produce the maximum amount of GROMACS trajectory over their lifetime. © 2015 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

PAPER—Accelerating parallel evaluations of ROCS

Modern graphics processing units (GPUs) are flexibly programmable and have peak computational throughput significantly faster than conventional CPUs. Herein, we describe the design and implementation of PAPER, an open‐source implementation of Gaussian molecular shape overlay for NVIDIA GPUs. We demonstrate one to two order‐of‐magnitude speedups on high‐end commodity GPU hardware relative to a reference CPU implementation of the shape overlay algorithm and speedups of over one order of magnitude relative to the commercial OpenEye ROCS package. In addition, we describe errors incurred by approximations used in common implementations of the algorithm. © 2009 Wiley Periodicals, Inc. J Comput Chem 2010     

Solvent effects on electronic structures and chain conformations of alpha-oligothiophenes in polar and apolar solutions

Solvent effects on electronic structures and chain conformations of alpha-oligothiophenes nTs (n = 1 to 10) are investigated in solvents of n-hexane, 1,4-dioxane, carbon tetrachloride, chloroform, and water by using density functional theory (DFT) and molecular dynamics (MD) simulations. Both implicit and explicit solvent models are employed. The polarized continuum model (PCM) calculations and MD simulations demonstrate the weak solvent effects on the electronic structures of alpha-oligothiophenes. The lowest dipole-allowed vertical excitation energies of nTs, obtained from time-dependent DFT/PCM calculations at the B3LYP/6-31G(d) level, exhibit a red shift as the solvent polarity increases, in agreement with experiments. The studied solvents have little impact on the state order of the low-lying excited states provided that the nTs are kept in C2h or C2v symmetry. The MD simulations demonstrate that the chain conformations are distorted to some extent in polar and nonpolar solvents. A qualitative picture of the distribution of solvent molecules around the solvated nTs is drawn by means of radial and spatial distribution functions. The S...H-O and pi...H-O solute-solvent interactions are insignificant in aqueous solution.     

Stability of ion triplets in ionic liquid/lithium salt solutions: Insights from implicit and explicit solvent models and molecular dynamics simulations

Binding energies of ion triplets formed in ionic liquids by Li⁺ with two anions have been studied using quantum‐chemical calculations with implicit and explicit solvent supplemented by molecular dynamics (MD) simulations. Explicit solvent approach confirms variation of solute‐ionic liquid interactions at distances up to 2 nm, resulting from structure of solvation shells induced by electric field of the solute. Binding energies computed in explicit solvent and from the polarizable continuum model approach differ largely, even in sign, but relative values generally agree between these two models. Stabilities of ion triplets obtained in quantum‐chemical calculations for some systems disagree with MD results; the discrepancy is attributed to the difference between static optimized geometries used in quantum chemical modeling and dynamic structures of triplets in MD simulations. © 2015 Wiley Periodicals, Inc.     

MOIL-opt: Energy-Conserving Molecular Dynamics on a GPU/CPU system

We report an optimized version of the molecular dynamics program MOIL that runs on a shared memory system with OpenMP and exploits the power of a Graphics Processing Unit (GPU). The model is of heterogeneous computing system on a single node with several cores sharing the same memory and a GPU. This is a typical laboratory tool, which provides excellent performance at minimal cost. Besides performance, emphasis is made on accuracy and stability of the algorithm probed by energy conservation for explicit-solvent atomically-detailed-models. Especially for long simulations energy conservation is critical due to the phenomenon known as "energy drift" in which energy errors accumulate linearly as a function of simulation time. To achieve long time dynamics with acceptable accuracy the drift must be particularly small. We identify several means of controlling long-time numerical accuracy while maintaining excellent speedup. To maintain a high level of energy conservation SHAKE and the Ewald reciprocal summation are run in double precision. Double precision summation of real-space non-bonded interactions improves energy conservation. In our best option, the energy drift using 1fs for a time step while constraining the distances of all bonds, is undetectable in 10ns simulation of solvated DHFR (Dihydrofolate reductase). Faster options, shaking only bonds with hydrogen atoms, are also very well behaved and have drifts of less than 1kcal/mol per nanosecond of the same system. CPU/GPU implementations require changes in programming models. We consider the use of a list of neighbors and quadratic versus linear interpolation in lookup tables of different sizes. Quadratic interpolation with a smaller number of grid points is faster than linear lookup tables (with finer representation) without loss of accuracy. Atomic neighbor lists were found most efficient. Typical speedups are about a factor of 10 compared to a single-core single-precision code.