Molecular dynamics on a distributed-memory multiprocessor

Dynamics simulations of molecular systems are notoriously computationally intensive. Using parallel computers for these simulations is important for reducing their turnaround time. In this article we describe a parallelization of the simulation program CHARMM for the Intel iPSC/860, a distributed memory multiprocessor. In the parallelization, the computational work is partitioned among the processors for core calculations including the calculation of forces, the integration of equations of motion, the correction of atomic coordinates by constraint, and the generation and update of data structures used to compute nonbonded interactions. Processors coordinate their activity using synchronous communication to exchange data values. Key data structures used are partitioned among the processors in nearly equal pieces, reducing the memory requirement per node and making it possible to simulate larger molecular systems. We examine the effectiveness of the parallelization in the context of a case study of a realistic molecular system. While effective speedup was achieved for many of the dynamics calculations, other calculations fared less well due to growing communication costs for exchanging data among processors. The strategies we used are applicable to parallelization of similar molecular mechanics and dynamics programs for distributed memory multiprocessors. © 1992 by John Wiley & Sons, Inc.     

Reactive Bond-Order Simulations Using Both Spatial and Temporal Approaches to Parallelism

We describe two different approaches to exploiting parallel computing architecture that have been used successfully for reactive molecular simulation using bond-order potentials. These potentials are based on the Tersoff bond-order formalism, and allow accurate treatement of covalent bonding reactions in the framework of a classical potential. They include both Brenner's reactive empirical bond order (REBO) potential and our adaptive intermolecular version of this potential (AIREBO). Traditional spatial and atom-based parallel decompositioon techniques have been employed in the RMD-CE program developed for parallel molecular dynamics simulations with a variety of reactive potentials. Key features of this implementation, including the object-oriented approach and novel algorithms for the integrator and neighbor lists, are discussed. The resulting code provides efficient scaling down to system sizes of 400 atoms per processor, and has been applied successfully to systems of as many as half a million atoms. For smaller systems, the parallel replica dynamics algorithm has been successfully applied to take advantage of parallelism in the time domain for rare-event systems. This approach takes advantage of the independence of different parts of a dynamics trajectory, and provides excellent parallel efficiencies for systems as small as tens of atoms, where other parallel simulation techniques are not applicable. This technique has been used to model the pyrolysis of hexadecane on the microsecond timescale, at more realistic temperatures than are achievable with other simulation methods.     

Hierarchical parallelization of divide-and-conquer density functional tight-binding molecular dynamics and metadynamics simulations

Massively parallel divide-and-conquer density functional tight-binding (DC-DFTB) molecular dynamics and metadynamics simulations are efficient approaches for describing various chemical reactions and dynamic processes of large complex systems via quantum mechanics. In this study, DC-DFTB simulations were combined with multi-replica techniques. Specifically, multiple walkers metadynamics, replica exchange molecular dynamics, and parallel tempering metadynamics methods were implemented hierarchically into the in-house Dcdftbmd program. Test simulations in an aqueous phase of the internal rotation of formamide and conformational changes of dialanine showed that the newly developed extensions increase the sampling efficiency and the exploration capabilities in DC-DFTB configuration space.     

An efficient parallelization scheme for molecular dynamics simulations with many-body, flexible, polarizable empirical potentials: application to water

An efficient parallelization scheme for classical molecular dynamics simulations with flexible, polarizable empirical potentials is presented. It is based on the standard Ewald summation technique to handle the long-range electrostatic and induction interactions. The algorithm for this parallelization scheme is designed for systems containing several thousands of polarizable sites in the simulation box. Its performance is evaluated during molecular dynamics simulations under periodic boundary conditions with unit cell sizes ranging from 128 to 512 molecules employing two flexible polarizable water models [DC(F) and TTM2.1-F] containing 1 and 3 polarizable sites, respectively. The time-to-solution for these two polarizable models is compared with the one for a flexible, pairwise-additive water model (TIP4F). The benchmarks were performed on both shared and distributed memory platforms. As a result of the efficient calculation of the induced dipole moments, a superlinear scaling as a function of the number of the processors is observed. To the best of our knowledge, this is the first reported results of parallel scaling and performance for simulations of liquid water with a polarizable potential under periodic boundary conditions.     

Development of massive multilevel molecular dynamics simulation program,platypus (PLATform for dYnamic protein unified simulation), for the elucidation of protein functions

A massively parallel program for quantum mechanical‐molecular mechanical (QM/MM) molecular dynamics simulation, called Platypus (PLATform for dYnamic Protein Unified Simulation), was developed to elucidate protein functions. The speedup and the parallelization ratio of Platypus in the QM and QM/MM calculations were assessed for a bacteriochlorophyll dimer in the photosynthetic reaction center (DIMER) on the K computer, a massively parallel computer achieving 10 PetaFLOPs with 705,024 cores. Platypus exhibited the increase in speedup up to 20,000 core processors at the HF/cc‐pVDZ and B3LYP/cc‐pVDZ, and up to 10,000 core processors by the CASCI(16,16)/6‐31G** calculations. We also performed excited QM/MM‐MD simulations on the chromophore of Sirius (SIRIUS) in water. Sirius is a pH‐insensitive and photo‐stable ultramarine fluorescent protein. Platypus accelerated on‐the‐fly excited‐state QM/MM‐MD simulations for SIRIUS in water, using over 4000 core processors. In addition, it also succeeded in 50‐ps (200,000‐step) on‐the‐fly excited‐state QM/MM‐MD simulations for the SIRIUS in water. © 2016 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

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.     

Computation of Absolute Hydration and Binding Free Energy with Free Energy Perturbation Distributed Replica-Exchange Molecular Dynamics (FEP/REMD)

Distributed Replica (REPDSTR) is a powerful parallelization technique enabling simulations of a group of replicas in a parallel/parallel fashion, where each replica is distributed to different nodes of a large cluster [Theor. Chem. Acc. 109: 140 (2003)]. Here, we use the framework provided by REPDSTR to combine a staged free energy perturbation protocol with different values of the thermodynamic coupling parameters with replica-exchange molecular dynamics (FEP/REMD. The structure of REPDSTR, which allows multiple parallel input/output (I/O), facilitates the treatment of replica-exchange to couple the N window simulations. As a result, each of the N synchronous window simulations benefit from the sampling carried out by the N-1 others. As illustrative examples of the FEP/REMD strategy, calculations of the absolute hydration and binding free energy of small molecules were performed using the biomolecular simulation program CHARMM adapted for the IBM Blue Gene/P platform. The computations show that a FEP/REMD strategy significantly improves the sampling and accelerate the convergence of absolute free energy computations.     

Replica exchanging self-guided Langevin dynamics for efficient and accurate conformational sampling

This work presents a replica exchanging self-guided Langevin dynamics (RXSGLD) simulation method for efficient conformational searching and sampling. Unlike temperature-based replica exchanging simulations, which use high temperatures to accelerate conformational motion, this method uses self-guided Langevin dynamics (SGLD) to enhance conformational searching without the need to elevate temperatures. A RXSGLD simulation includes a series of SGLD simulations, with simulation conditions differing in the guiding effect and∕or temperature. These simulation conditions are called stages and the base stage is one with no guiding effect. Replicas of a simulation system are simulated at the stages and are exchanged according to the replica exchanging probability derived from the SGLD partition function. Because SGLD causes less perturbation on conformational distribution than high temperatures, exchanges between SGLD stages have much higher probabilities than those between different temperatures. Therefore, RXSGLD simulations have higher conformational searching ability than temperature based replica exchange simulations. Through three example systems, we demonstrate that RXSGLD can generate target canonical ensemble distribution at the base stage and achieve accelerated conformational searching. Especially for large systems, RXSGLD has remarkable advantages in terms of replica exchange efficiency, conformational searching ability, and system size extensiveness.     

High-temperature decomposition of the cellulose molecule: a stochastic molecular dynamics study

The kinetics and products of cellulose pyrolysis can be studied using large-scale molecular dynamics simulations at high temperatures, where the reaction rates are high enough to make the simulation times practical. We carried out molecular dynamics simulations employing the ReaxFF reactive force field to study the initial step of the thermal decomposition process. We gathered statistics of simulated reactive events at temperatures ranging from 1400 to 2200 K, considering cellulose molecules with different molecular weights and initial conformations. Our simulations suggest that, in gas-phase conditions at these high temperatures, the decomposition occurs primarily through random cleavage of the β(1 → 4)-glycosidic bonds, for which we obtained an activation energy of (171 ± 2) kJ mol^?1 and a frequency factor of \(\left( {1.07 \pm 0.12} \right) \times 10^{15}\) s^?1. We did not observe dependency of the kinetic parameters on the molecular weight or initial conformation. Some of the decomposition reactions involved the release of low-molecular-weight products. Excluding radicals, the most commonly observed species were glycolaldehyde, water, formaldehyde and formic acid. Many of our observations are supported by the existing experimental and theoretical knowledge. We did not, however, observe the formation of levoglucosan, which is the dominant product in conventional pyrolysis experiments at much lower temperatures. This is understandable, since the high temperatures can force the dominance of radical reactions over pericyclic reactions. Nevertheless, our results support further use of ReaxFF-based molecular dynamics simulations in the study of cellulose pyrolysis.     

Computing Alchemical Free Energy Differences with Hamiltonian Replica Exchange Molecular Dynamics (H-REMD) Simulations

Alchemical free energy calculations play a very important role in the field of molecular modeling. Efforts have been made to improve the accuracy and precision of those calculations. One of the efforts is to employ a Hamiltonian replica exchange molecular dynamics (H-REMD) method to enhance conformational sampling. In this paper, we demonstrated that HREMD method not only improves convergence in alchemical free energy calculations but also can be used to compute free energy differences directly via the Free Energy Perturbation (FEP)algorithm. We show a direct mapping between the H-REMD and the usual FEP equations, which are then used directly to compute free energies. The H-REMD alchemical free energy calculation (Replica exchange Free Energy Perturbation, REFEP) was tested on predicting the pK(a) value of the buried Asp26 in thioredoxin. We compare the results of REFEP with TI and regular FEP simulations. REFEP calculations converged faster than those from TI and regular FEP simulations. The final predicted pK(a) value from the H-REMD simulation was also very accurate, only 0.4 pK(a) unit above the experimental value. Utilizing the REFEP algorithm significantly improves conformational sampling, and this in turn improves the convergence of alchemical free energy simulations.     

A new parallel method for molecular dynamics simulation of macromolecular systems

Short-range molecular dynamics simulations of molecular systems are commonly parallelized by replicated-data methods, in which each processor stores a copy of all atom positions. This enables computation of bonded 2-, 3-, and 4-body forces within the molecular topology to be partitioned among processors straightforwardly. A drawback to such methods is that the interprocessor communication scales as N (the number of atoms) independent of P (the number of processors). Thus, their parallel efficiency falls off rapidly when large numbers of processors are used. In this article a new parallel method for simulating macromolecular or small-molecule systems is presented, called force-decomposition. Its memory and communication costs scale as N/√P, allowing larger problems to be run faster on greater numbers of processors. Like replicated-data techniques, and in contrast to spatial-decomposition approaches, the new method can be simply load balanced and performs well even for irregular simulation geometries. The implementation of the algorithm in a prototypical macromolecular simulation code ParBond is also discussed. On a 1024-processor Intel Paragon, ParBond runs a standard benchmark simulation of solvated myoglobin with a parallel efficiency of 61% and at 40 times the speed of a vectorized version of CHARMM running on a single Cray Y-MP processor. © 1996 by John Wiley & Sons, Inc.     

Midpoint cell method for hybrid (MPI+OpenMP) parallelization of molecular dynamics simulations

We have developed a new hybrid (MPI+OpenMP) parallelization scheme for molecular dynamics (MD) simulations by combining a cell‐wise version of the midpoint method with pair‐wise Verlet lists. In this scheme, which we call the midpoint cell method, simulation space is divided into subdomains, each of which is assigned to a MPI processor. Each subdomain is further divided into small cells. The interaction between two particles existing in different cells is computed in the subdomain containing the midpoint cell of the two cells where the particles reside. In each MPI processor, cell pairs are distributed over OpenMP threads for shared memory parallelization. The midpoint cell method keeps the advantages of the original midpoint method, while filtering out unnecessary calculations of midpoint checking for all the particle pairs by single midpoint cell determination prior to MD simulations. Distributing cell pairs over OpenMP threads allows for more efficient shared memory parallelization compared with distributing atom indices over threads. Furthermore, cell grouping of particle data makes better memory access, reducing the number of cache misses. The parallel performance of the midpoint cell method on the K computer showed scalability up to 512 and 32,768 cores for systems of 20,000 and 1 million atoms, respectively. One MD time step for long‐range interactions could be calculated within 4.5 ms even for a 1 million atoms system with particle‐mesh Ewald electrostatics. © 2014 Wiley Periodicals, Inc.     

The flexible, polarizable, thole-type interaction potential for water (TTM2-F) revisited

We present a revision of the flexible, polarizable, Thole-type interaction potential for water [J. Chem. Phys.2002, 116, 5115], which allows for condensed-phase simulations. The revised version (TTM2.1-F) of the potential correctly describes the individual water molecular dipole moment and alleviates problems arising at short intermolecular separations that can be sampled in the course of molecular dynamics and Monte Carlo simulations of condensed environments. Furthermore, its parallel implementation under periodic boundary conditions enables the efficient calculation of the macroscopic structural and thermodynamic properties of liquid water, as its performance scales superlinearly with up to a number of 64 processors for a simulation box of 512 molecules. We report the radial distribution functions, average energy, internal geometry, and dipole moment in the liquid as well as the density, dielectric constant, and self-diffusion coefficient at T = 300 K from (NVT) and (NPT) classical molecular dynamics simulations by using the revised version of the potential.     

Insights into the dynamics of HIV-1 protease: a kinetic network model constructed from atomistic simulations

The conformational dynamics in the flaps of HIV-1 protease plays a crucial role in the mechanism of substrate binding. We develop a kinetic network model, constructed from detailed atomistic simulations, to determine the kinetic mechanisms of the conformational transitions in HIV-1 PR. To overcome the time scale limitation of conventional molecular dynamics (MD) simulations, our method combines replica exchange MD with transition path theory (TPT) to study the diversity and temperature dependence of the pathways connecting functionally important states of the protease. At low temperatures the large-scale flap opening is dominated by a small number of paths; at elevated temperatures the transition occurs through many structurally heterogeneous routes. The expanded conformation in the crystal structure 1TW7 is found to closely mimic a key intermediate in the flap-opening pathways at low temperature. We investigated the different transition mechanisms between the semi-open and closed forms. The calculated relaxation times reveal fast semi-open ? closed transitions, and infrequently the flaps fully open. The ligand binding rate predicted from this kinetic model increases by 38-fold from 285 to 309 K, which is in general agreement with experiments. To our knowledge, this is the first application of a network model constructed from atomistic simulations together with TPT to analyze conformational changes between different functional states of a natively folded protein.     

Hybrid particle‐field molecular dynamics simulations: Parallelization and benchmarks

The parallel implementation of a recently developed hybrid scheme for molecular dynamics (MD) simulations (Milano and Kawakatsu, J Chem Phys 2009, 130, 214106) where self‐consistent field theory (SCF) and particle models are combined is described. Because of the peculiar formulation of the hybrid method, considering single particles interacting with density fields, the most computationally expensive part of the hybrid particle‐field MD simulation can be efficiently parallelized using a straightforward particle decomposition algorithm. Benchmarks of simulations, including comparisons of serial MD and MD‐SCF program profiles, serial MD‐SCF and parallel MD‐SCF program profiles, and parallel benchmarks compared with efficient MD program GROMACS 4.5.4 are tested and reported. The results of benchmarks indicate that the proposed parallelization scheme is very efficient and opens the way to molecular simulations of large scale systems with reasonable computational costs. © 2012 Wiley Periodicals, Inc.     

A scalable divide-and-conquer algorithm combining coarse and fine-grain parallelization

We describe an efficient algorithm for carrying out a “divide-and-conquer” fit of a molecule's electronic density on massively parallel computers. Near linear speedups are achieved with up to 48 processors on a Cray T3E, and our results indicate that similar efficiencies could be attained on an even greater number of processors. To achieve optimum efficiency, the algorithm combines coarse and fine-grain parallelization and adapts itself to the existing ratio of processors to subsystems. The subsystems employed in our divide-and-conquer approach can also be made smaller or bigger, depending on the number of processors available. This allows us to further reduce the wallclock time and improve the method's overall efficiency. The strategies implemented in this paper can be extended to any other divide-and-conquer method used within an ab initio, density functional, or semi-empirical quantum mechanical program. Received: 15 September 1997 / Accepted: 21 January 1998     

Reversible folding simulation by hybrid Hamiltonian replica exchange

Reversible foldings of a beta-hairpin peptide, chignolin, by recently invented hybrid Hamiltonian replica exchange molecular dynamics simulations based on Poisson-Boltzmann model in explicit water are demonstrated. Initiated from extended structures the peptide folded and unfolded a couple of times in seven out of eight replica trajectories during 100 nanoseconds simulation. The folded states have the lowest all-atom root mean squared deviation of 1.3 A with respect to the NMR structures. At T=300 K the occurrence of folded states was converged to 62% during 80 ns simulation which agrees well with experimental data. Especially, a detailed structural evolution map was constructed based on 800,000 structural snapshots and from where a unique folding doorway emerges. Compared with 130 ns standard replica exchange simulation using 24 replicas on the same system, the hybrid Hamiltonian replica exchange molecular dynamics simulation presents consistent results.     

Efficient evaluation of sampling quality of molecular dynamics simulations by clustering of dihedral torsion angles and Sammon mapping

A direct conformational clustering and mapping approach for peptide conformations based on backbone dihedral angles has been developed and applied to compare conformational sampling of Met-enkephalin using two molecular dynamics (MD) methods. Efficient clustering in dihedrals has been achieved by evaluating all combinations resulting from independent clustering of each dihedral angle distribution, thus resolving all conformational substates. In contrast, Cartesian clustering was unable to accurately distinguish between all substates. Projection of clusters on dihedral principal component (PCA) subspaces did not result in efficient separation of highly populated clusters. However, representation in a nonlinear metric by Sammon mapping was able to separate well the 48 highest populated clusters in just two dimensions. In addition, this approach also allowed us to visualize the transition frequencies between clusters efficiently. Significantly, higher transition frequencies between more distinct conformational substates were found for a recently developed biasing-potential replica exchange MD simulation method allowing faster sampling of possible substates compared to conventional MD simulations. Although the number of theoretically possible clusters grows exponentially with peptide length, in practice, the number of clusters is only limited by the sampling size (typically much smaller), and therefore the method is well suited also for large systems. The approach could be useful to rapidly and accurately evaluate conformational sampling during MD simulations, to compare different sampling strategies and eventually to detect kinetic bottlenecks in folding pathways.     

SIMONA 1.0: An efficient and versatile framework for stochastic simulations of molecular and nanoscale systems

Molecular simulation methods have increasingly contributed to our understanding of molecular and nanoscale systems. However, the family of Monte Carlo techniques has taken a backseat to molecular dynamics based methods, which is also reflected in the number of available simulation packages. Here, we report the development of a generic, versatile simulation package for stochastic simulations and demonstrate its application to protein conformational change, protein–protein association, small-molecule protein docking, and simulation of the growth of nanoscale clusters of organic molecules. Simulation of molecular and nanoscale systems (SIMONA) is easy to use for standard simulations via a graphical user interface and highly parallel both via MPI and the use of graphical processors. It is also extendable to many additional simulations types. Being freely available to academic users, we hope it will enable a large community of researchers in the life- and materials-sciences to use and extend SIMONA in the future. SIMONA is available for download under http://int.kit.edu/nanosim/simona . © 2012 Wiley Periodicals, Inc.