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/.

The molecular dynamics (MD) method is a key technology for dissecting dynamical properties of molecular systems at the atomic level. Along with the rapid growing of the high performance computing technologies, the field of MD simulations has been growing and extensively applied for highly complex and large-scale phenomena. A recent milestone in this field is the rise of the special-purpose hardware, named “ANTON” [1,2]. This hardware achieved millisecond-timescale simulations for a variety of molecular systems [3,4]. Yet another MD-specialized hardware, “MD-GRAPE”, was recently developed by Taiji et al. in RIKEN [5]. While specialized hardware pushes the limits of the MD method, it is not widely used due to its very high costs and difficulty in modifying the algorithms set in the hardware.

Consequently, developments in the field of MD simulations have mainly targeted general-purpose cluster machines and a variety of MD software have been developed. At the moment, one of the most efficient implementations of the MD is “GROMACS”, which comes with implementations of low level single-instruction-multiple-data (SIMD) instructions [6,7]. The software “AMBER” also accomplishes a good performance, especially when powered by NVIDIA GPU-based accelerator [8,9]. In parallel, several MD programs based on unique state-of-the-art physics theories have also been developed. The MD program “GENESIS” takes advantage of a variety of replica exchange MD (REMD) methods, e.g., the temperature REMD and the surface tension REMD, and several original algorithms [10-12]. “MODYLAS” applies a non-Ewald electrostatic potential calculation scheme, the fast multipole method [13]. Another MD program “MARBLE” implements a partial rigid-body method under membrane-specific ensembles [14]. “SCUBA” has a unique algorithm for generating the structural ensemble fitted to low resolution electron microscopy data [15]. These programs are extensively tuned for supercomputer systems, such as the “K computer”.

During the past several decades, we developed and maintained “myPresto” package, which consists of a variety of tools for in-silico drug development, including a docking program between a receptor protein and drug candidates (“sievegene”) [16-18], a ligand binding pocket finder (“molsite”) [19], a small-compound structure library (“ligandbox”) [20]. More recently added to myPresto was the MD simulation engine “psygene-G”, which implemented two major original features: a multi-GPU accelerated parallelization scheme [21], and the zero-multipole summation method (ZMM) [22-26]. By utilizing many GPUs to parallelize the space-decomposition routine, myPresto/psygene-G accomplishes a scalable computation for a wide range of system sizes and have been applied for a variety of biomolecular systems [27-29]. As the ZMM estimates the electrostatic potential with a cutoff pair potential, the computational cost is drastically reduced compared with Ewald-based methods, and the scalability of the multi-node parallel computation is improved.

In addition, we developed a series of theories for the non-canonical MD simulations, namely enhanced conformational sampling methods. We have developed the multi-canonical ensemble MD simulation method [30], and here in after we refer this method as multi-canonical MD (McMD). On the basis of the McMD method, several enhanced conformational sampling methods have been developed and applied to analyze a variety of bio-molecular systems [31-34]. These methods are powerful tools to elucidate physicochemical properties of molecular systems at an equilibrium state. We have previously demonstrated that trajectories generated from many independent McMD runs with different initial conditions can be combined for statistical analysis [35]. This theory allows us to execute many relatively short MD runs independently in parallel instead of a single long simulation.

Our enhanced conformational sampling methods and the ZMM have been implemented in myPresto/psygene-G. However, myPresto/psygene-G was tailored for long-term simulations of large molecular systems using multi-GPU parallelization with the space-decomposition algorithm, and thus it is not suitable for running multiple independent simulations of small systems, which are what enhanced conformational sampling methods target.

Here, we developed a new MD simulation program from scratch, named “myPresto/omegagene” that is tailored for our original McMD methods with the ZMM. myPresto/omegagene is distributed under open-source license from http://www.protein.osaka-u.ac.jp/rcsfp/pi/omegagene/, and it is also available from http://presto.protein.osaka-u.ac.jp/myPresto4/. In this report, Section 1 “Software Details” presents requirements and detailed information of the implementation. Section 2 “Simulations” presents results of two types of test simulations, with a 20-mer peptide, Trp-cage. The first simulation with the micro-canonical ensemble demonstrates a base capability: four-fold acceleration compared with myPresto/psygene-G with an energy conservation property. The second simulation with the virtual-system-coupled McMD (V-McMD) method, which is a variant of the McMD method, shows successful application results: an energy landscape of the Trp-cage at 300 K along with near-native conformations sampled as the most stable cluster.

1. Software Details

1.1. Overview

myPresto/omegagene consists of the two parts: i) the core MD engine written in C++ and CUDA languages, and ii) the toolkit for pre- and post-simulation processing written in Python language (Fig. 1). The installation process is semi-automated, by taking advantage of a standard cross-platform building tool “cmake” (version 3.2 or later). myPresto/omegagene is easily built on a variety of environment, e.g., Linux, Windows, and OS X (also known as MacOS). For compiling the programs, C++ compilers compatible to C++11 standard is required. The GPU acceleration requires a NVIDIA GPU board with compute capability ≥3.5.

Figure 1

Overview of the myPresto/omegagene ecosystem.

The input and output files formats of myPresto/omegagene are compatible with those of existing tools in the myPresto package. The input molecular topologies of proteins and nucleic acids are generated by a “tplgene” program, a myPresto tool. A list of atom groups with non-holonomic constraints is prepared, by using the “SHAKEinp” program. In addition, some other tools for pre- and post-simulation analyses and data handling are attached as omega_toolkit. For example, the simulation trajectory written in PRESTO format can be converted into other standard binary formats, e.g., .trr format (GROMACS). The initial velocities of atoms in simulation systems are generated by using a Python script included in omega_toolkit. These input files are then integrated into a single binary file, by using another script. Separation of such pre- and post-simulation lightweight tasks from the core engine makes it easier to develop and maintain the tools.

1.2. Implementation

The MD engine of myPresto/omegagene is implemented in C++11 in an object-oriented manner. SIMD parallelization, powered by NVIDIA GPU with CUDA was applied for the cutoff pair-wise potential calculations, which is the bottle-neck of MD computations. As this program utilizes the ZMM for the electrostatic potential calculation, Lennard-Jones and electrostatic potentials are computed in the same kernel function on GPU. A main difference from myPresto/psygene-G is the routine for pairwise potential calculations on GPU. myPresto/omegagene applies an algorithm based on the method proposed by Páll and Pronk [6] which is utilized in GROMACS. To find atom pairs within the cutoff radius (the neighbor search task), myPresto/omegagene groups the atoms in the system into “sub-cells” with the same number of atoms, in contrast to myPresto/psygene-G, which divides the system into sub-cells with a uniform volume [21].

As the SIMD width of NVIDIA GPU is 32 threads (a bunch of 32 threads is defined as a “warp”), the atoms in the system are grouped into 8-atom sub-cells (Fig. 2A). Accordingly, there are 64 atom-atom interactions between atoms from each sub-cell in a pair, and these pair potentials are calculated in bulk by one warp (32 threads) using two loop iterations. Before the pair potential calculations, a search for all pairs of sub-cells whose inter-cell distance are within a distance threshold (the summation over the cutoff radius and an offset value) is performed at every predefined interval number of steps (test calculations for this interval are described in the next section). The neighbor search process enumerates all sub-cell pairs within this cutoff. For example, in Figure 2A, the neighbor sub-cells of the 1st sub-cell are 0th, 2nd, 3rd, and 4th sub-cells. In one iteration of the pairwise potential calculation algorithm, the m-th warp calculates pairs of the m-th sub-cell and its neighbor sub-cells. In order to balance the number of sub-cell pairs owned by each warp, we developed our original scheduling scheme, in which a sub-cell pair m–n (m and n refer to sub-cell indices) is considered for calculation if either one of following conditions hold: 1) when m is an even number, n is an odd number for mm; and 2) when m is an odd number, n is an even number for nm; or m=n. Using Figure 2B as an example, the 1st warp in the GPU grid calculates pairwise potentials in the sub-cell pairs 1–0, 1–1, and 1–3, while the 2nd and 4th warps calculate pairwise potentials in sub-cell pairs 2–1 and 4–1, respectively. Note that this algorithm does not guarantee a good load balance among warps, because not all sub-cells have the same number of neighbors. While Páll and Pronk’s method limits the number of sub-cell pairs calculated in each warp to eight in order to balance the workload of each warp, we did not apply this scheme here, since dividing the loops requires frequent communications between the registers and the global memory with atomic operations. In addition, as it is expected that the distribution of the particle density is not largely biased in usual systems, we assume the impact of this load imbalance is limited.

Figure 2

Schematic diagram of the algorithm for the pairwise potential calculations. (A) The system consisting of sub-cells, each of which consists of eight atoms (filled circles). The number in the circle indicates the ID of each sub-cell. The red dashed line indicates the cutoff length from the border of the 1st sub-cell. The blue dashed line represents the cutoff added by the offset value. (B) The matrix describing assignments of each sub-cell pairs to each warp. The pairs marked as the orange squares in each row are processed by a warp. Note that the periodic boundary condition is ignored for simplicity in this figure. (C) Pairs of atoms calculated in each thread in warp 1. The numbers in the matrix denote the thread IDs in the warp, processing the corresponding pair of atoms. Each row indicates one of eight atoms in 1st sub-cell. Each column indicates an atom nearby the 1st sub-cell. The filled circles correspond to the atoms in Figure 2A. (D) A pseudo code of the kernel function. num_neighbor_cubcells[m] is the number of neighboring sub-cells of m-th sub-cell. is Masked (i, j) excludes the i–j pairs within four successive covalent bonds and the pairs i ≤ j in the same sub-cell. force[i] is the 3D array keeping the force of the i-th atom along x, y, and z axes. The function cal Force(i, j) calculates the force of the i-th atom acted by the j-th atom. warpShuffle(force[i]) sums up the force values among threads with the same i value in the same warp. atomicAdd adds the force[i] value to the variable globalmem_force[I] in the global memory.

Then, the pairwise potentials are calculated in each warp (Fig. 2C and D). As each sub-cell has eight atoms, there are 64 atom pairs in each pair of sub-cells, and these pair potentials are calculated in two cycles of the loop, by using 32 threads (one warp). The variable l means the iterator variable of the loop. The first cycle (l=0) of the 1st warp, pairs of eight atoms in the 1st sub-cell and four atoms in the 0th sub-cell are computed. In the next cycle (l=1), pairs of the eight atoms in the 1st sub-cell and the remaining four atoms in the sub-cell 0 are computed. The calculated forces for each atom are gathered with a warp shuffle instruction and are written back to the global memory with an atomic operation. While the calculations of force and energy in each term are executed in the single precision, summation of these terms (warpShuffle and atomicAdd in Fig. 2D) is done in the double precision. To exclude the i–j atom pairs connected within four successive covalent bonds and the atom pairs i ≤ j in a sub-cell (m=n; shaded pairs in Fig. 2C), a 64-bit bitmask is prepared for each sub-cell pair to switch the pair potential of each atom pair [36].

2. Simulations

To demonstrate the performance of myPresto/omegagene, we carried out micro-canonical simulations and enhanced conformational sampling with the V-McMD method. The former was performed to evaluate the energy drift and compare the computation speed compared with that of myPresto/psygene-G. The system used in the both simulations includes a 20-mer peptide, Trp-cage with the sequence NLYIQWLKDGGPSSGRPPPS. For the potential parameters, Amber99SB-ILDN force field [37] for the protein, TIP3P water model [38], and the ion model reported by Joung and Cheatham [39] were applied. The electrostatic potential was calculated by ZMM with the zero-dipole condition and dumping parameter α=0 [22]. The integration time step was 2.0 fs and covalent bonds with hydrogen atoms were constrained by the SHAKE algorithm [40].

2.1. Micro-canonical Simulation

We ran micro-canonical simulations with myPresto/omegagene for 1.0 ns, using several settings of intervals for the neighbor search, e.g., 1, 10, 30, and 50 steps. The initial structure of Trp-cage was taken from the solution NMR structure (PDB ID: 1L2Y, model 1) [41], immersed into 150 mM NaCl solution (13,277 atoms; Supplementary Fig. S1A). A 12 Å cutoff radius was applied for the Lennard-Jones and electrostatic interactions and a 1 Å offset value was applied for the neighbor search. The same system was also simulated with myPresto/psygene-G for comparison, whose computation speeds were described for many protein systems [21].

We found that while a 10-step interval for the neighbor search improves the computation time and maintains the accuracy of simulation run that runs the neighbor search routine at every step, further elongation of the interval (30 and 50) causes drifts in the total energy (Fig. 3A and B). Thus, we recommend a 10-step interval as the default setting. Although even the conditions with small intervals (1 and 10 steps) in myPresto/omegagene showed a drift of the energy as well as myPresto/psygene-G does, the drift is so small that the results would be acceptable for usual purposes. Under the 10-step interval setting, the computation time with a single GPU is 4.12-times faster than the execution by myPresto/psygene-G with the same resource, and 1.68-times slower than the execution by myPresto/psygene-G with eight GPUs (the 2×2×2 space decomposition). The bottle-neck of computations was the pair potential functions on GPU (78.5% of the computation time was occupied; Supplementary Fig. S2) and the other potential calculations do not affect the computation speed because they are performed by CPU behind the GPU task. The integration task including the SHAKE algorithm spent 16.6% of the computation time. The remaining 5.9% was for the neighbor search. Although our main targets of myPresto/omegagene are small molecular systems with enhanced conformational sampling methods, we also tested a large molecular system consisting of an 11-subunit assembly of RNA polymerase (PDB ID: 4qiw; Supplementary Figs. S1C and S2B). The computation speed was approximately 0.27 ns/day for this 670,957-atom system by a single CPU/GPU.

Figure 3

Benchmark results of myPresto/omegagene. (A) Time courses of the total energy of the system in various conditions. “NS interval” means the interval steps of the neighbor search. “GPU” is the number of GPU boards used in the simulation. The vertical axis indicates the ratio of the total energy over that at the first step. The value overlaid on each plot is the slope of fitted linear function. In these calculations, the time step was 2.0 fs, and the SHAKE algorithm was applied. (B) The computation speed in each condition.

2.2. McMD simulation

We applied myPresto/omegagene to conformational sampling with the V-McMD method [31]. The initial structure of Trp-cage was generated using “Modeller” without any template, and randomized by a constant temperature simulation at a high temperature (800 K; Supplementary Fig. S1B). The solution contains only one Cl− for neutralization of the system (15,596 atoms in total). The same force field used for the earlier micro-canonical simulations was applied.

We briefly explain the protocols of the V-McMD for the present study below. A detailed elaboration of the methodology can be found in the reference [42]. The initial guess of the density of states was first approximated from a series of constant temperature simulations at various temperatures, ranging from 296 K to 629 K. Based on this result, the first iteration of the V-McMD simulation was performed under a biasing force. The energy distribution obtained from this simulation (Supplementary Fig. S3A) was then used to refine the estimated density of states and the biasing force, and this iterative process was performed 34 times (1.8×109 steps in total), until the energy distribution became nearly flat. After obtaining the biasing force, which produces the flat energy distribution, 50 independent runs with 1.5×107 steps for each were performed as the production runs (Supplementary Fig. S3B). The canonical ensemble at 300 K of this system was generated by reweighting the probability of each snapshot. The conformational ensemble of Trp-cage at 300 K is summarized as the potential of mean force (PMF) landscape, by applying principal component analysis (PCA), where 190 interatomic distances for the all pairs of Cα atoms in the peptide were used for variables to define the variance-covariance matrix [34]. As a result, we obtained the landscape with a major basin and some peripheral shallow basins (Fig. 4A). The most thermodynamically stable basin corresponds to the native structure (Fig. 4B) and the Cα RMSD value between one of the most stable structures and the NMR structure (PDB ID: 1L2Y, model 1) is 0.944 Å. This result demonstrates that myPresto/omegagene reasonably analyzed the conformational ensemble of the 20-mer peptide. In this calculation, the calculation speed was 10.7 ns/day, which is four-times faster than myPresto/psygene-G for the V-McMD with the ZMM.

Figure 4

The results of the V-McMD simulation. (A) The simulated PMF landscape at 300 K. The horizontal and vertical axes are the first and second principal components (PC1 and PC2), respectively. The color graduation represents the PMF values (blue means the most stable). The lowest PMF is set to zero kJ/mol. (B) A superimposed picture of the NMR structure of Trp-cage (gray; PDB ID: 1L2Y, model 1) and the most stable conformation in the V-McMD simulation (blue).

3. Conclusions

We developed a new MD simulation program, myPresto/omegagene, that is tailored for our original enhanced sampling methods [31,32] and the electrostatic potential scheme [22-26]. This software is freely available to users, and the source code is distributed under an open-source license. In contrast to myPresto/psygene-G, which is powered by multi-GPU parallel computations [21], myPresto/omegagene is tailored for a single process execution with a single GPU, in order to optimize the enhanced conformational sampling methods [35]. Elimination of space-decomposition routines from the code greatly simplifies the codebase and the simple object-oriented structure of the code allows for easy development and maintenance of the software. In addition, the compatibility with the myPresto family provides advantages for applying many existing tools accumulated during the past decades and easy to use for myPresto users.

The evaluation study with the micro-canonical ensemble demonstrated the acceptable properties of the energy drift and efficient computation. Furthermore, we demonstrated that the enhanced conformational sampling simulation successfully reproduced experimentally solved structure of Trpcage, as the most thermodynamically stable structure in the simulated conformational ensemble. myPresto/omegagene can be applied for the conformational sampling of such a practical molecular system with computation speed about four-times faster than that of myPresto/psygene-G.

The molecular dynamics (MD) method is a promising approach to investigate a variety of biophysical phenomena at the atomic level. In particular, elucidation of the conformational ensemble of bio-molecules is of paramount importance. We present here a new MD simulation program, “myPresto/omegagene” which is tailored for efficient computation of enhanced conformational sampling powered by GPU acceleration. “myPresto/omegagene” is unique in that it adopts our original multi-canonical ensemble approach and our non-Ewald electrostatic potential scheme named “zero-multipole summation method” to effectively enhance computations. In addition, “myPresto/omegagene” builds upon a wealth of tools and resources provided by the myPresto suite to enable good user experience for myPresto users.