A comparison between parallelization approaches in molecular dynamics simulations on GPUs

We test the relative performances of two different approaches to the computation of forces for molecular dynamics simulations on graphics processing units. A “vertex‐based” approach, where a computing thread is started per particle, is compared to an “edge‐based” approach, where a thread is started per each potentially non‐zero interaction. We find that the former is more efficient for systems with many simple interactions per particle while the latter is more efficient if the system has more complicated interactions or fewer of them. By comparing computation times on more and less recent graphics processing unit technology, we predict that, if the current trend of increasing the number of processing cores—as opposed to their computing power—remains, the “edge‐based” approach will gradually become the most efficient choice in an increasing number of cases. © 2014 Wiley Periodicals, Inc.     

Difficulties with multiple time stepping and fast multipole algorithm in molecular dynamics

Numerical experiments are performed on a 36,000-atom protein–DNA–water simulation to ascertain the effectiveness of two devices for reducing the time spent computing long-range electrostatics interactions. It is shown for Verlet-I/r-RESPA multiple time stepping, which is based on approximating long-range forces as widely separated impulses, that a long time step of 5 fs results in a dramatic energy drift and that this is reduced by using an even larger long time step. It is also shown that the use of as many as six terms in a fast multipole algorithm approximation to long-range electrostatics still fails to prevent significant energy drift even though four digits of accuracy is obtained. © 1997 John Wiley & Sons, Inc. J Comput Chem 18 : 1785–1791, 1997     

Reaching multi‐nanosecond timescales in combined QM/MM molecular dynamics simulations through parallel horsetail sampling

We report an enhanced sampling technique that allows to reach the multi‐nanosecond timescale in quantum mechanics/molecular mechanics molecular dynamics simulations. The proposed technique, called horsetail sampling, is a specific type of multiple molecular dynamics approach exhibiting high parallel efficiency. It couples a main simulation with a large number of shorter trajectories launched on independent processors at periodic time intervals. The technique is applied to study hydrogen peroxide at the water liquid–vapor interface, a system of considerable atmospheric relevance. A total simulation time of a little more than 6 ns has been attained for a total CPU time of 5.1 years representing only about 20 days of wall‐clock time. The discussion of the results highlights the strong influence of the solvation effects at the interface on the structure and the electronic properties of the solute. © 2017 Wiley Periodicals, Inc.     

Algorithms for GPU‐based molecular dynamics simulations of complex fluids: Applications to water,mixtures, and liquid crystals

A custom code for molecular dynamics simulations has been designed to run on CUDA‐enabled NVIDIA graphics processing units (GPUs). The double‐precision code simulates multicomponent fluids, with intramolecular and intermolecular forces, coarse‐grained and atomistic models, holonomic constraints, Nosé–Hoover thermostats, and the generation of distribution functions. Algorithms to compute Lennard‐Jones and Gay‐Berne interactions, and the electrostatic force using Ewald summations, are discussed. A neighbor list is introduced to improve scaling with respect to system size. Three test systems are examined: SPC/E water; an n‐hexane/2‐propanol mixture; and a liquid crystal mesogen, 2‐(4‐butyloxyphenyl)‐5‐octyloxypyrimidine. Code performance is analyzed for each system. With one GPU, a 33–119 fold increase in performance is achieved compared with the serial code while the use of two GPUs leads to a 69–287 fold improvement and three GPUs yield a 101–377 fold speedup. © 2015 Wiley Periodicals, Inc.     

Increasing the time step with mass scaling in Born‐Oppenheimer ab initio QM/MM molecular dynamics simulations

Born‐Oppenheimer ab initio QM/MM molecular dynamics simulation with umbrella sampling is a state‐of‐the‐art approach to calculate free energy profiles of chemical reactions in complex systems. To further improve its computational efficiency, a mass‐scaling method with the increased time step in MD simulations has been explored and tested. It is found that by increasing the hydrogen mass to 10 amu, a time step of 3 fs can be employed in ab initio QM/MM MD simulations. In all our three test cases, including two solution reactions and one enzyme reaction, the resulted reaction free energy profiles with 3 fs time step and mass scaling are found to be in excellent agreement with the corresponding simulation results using 1 fs time step and the normal mass. These results indicate that for Born‐Oppenheimer ab initio QM/MM molecular dynamics simulations with umbrella sampling, the mass‐scaling method can significantly reduce its computational cost while has little effect on the calculated free energy profiles. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2009     

Evaluation of atomic pressure in the multiple time‐step integration algorithm

In molecular dynamics (MD) calculations, reduction in calculation time per MD loop is essential. A multiple time‐step (MTS) integration algorithm, the RESPA (Tuckerman and Berne, J. Chem. Phys. 1992, 97, 1990–2001), enables reductions in calculation time by decreasing the frequency of time‐consuming long‐range interaction calculations. However, the RESPA MTS algorithm involves uncertainties in evaluating the atomic interaction‐based pressure (i.e., atomic pressure) of systems with and without holonomic constraints. It is not clear which intermediate forces and constraint forces in the MTS integration procedure should be used to calculate the atomic pressure. In this article, we propose a series of equations to evaluate the atomic pressure in the RESPA MTS integration procedure on the basis of its equivalence to the Velocity‐Verlet integration procedure with a single time step (STS). The equations guarantee time‐reversibility even for the system with holonomic constrants. Furthermore, we generalize the equations to both (i) arbitrary number of inner time steps and (ii) arbitrary number of force components (RESPA levels). The atomic pressure calculated by our equations with the MTS integration shows excellent agreement with the reference value with the STS, whereas pressures calculated using the conventional ad hoc equations deviated from it. Our equations can be extended straightforwardly to the MTS integration algorithm for the isothermal NVT and isothermal–isobaric NPT ensembles. © 2017 Wiley Periodicals, Inc.     

Efficient nonbonded interactions for molecular dynamics on a graphics processing unit

We describe an algorithm for computing nonbonded interactions with cutoffs on a graphics processing unit. We have incorporated it into OpenMM, a library for performing molecular simulations on high‐performance computer architectures. We benchmark it on a variety of systems including boxes of water molecules, proteins in explicit solvent, a lipid bilayer, and proteins with implicit solvent. The results demonstrate that its performance scales linearly with the number of atoms over a wide range of system sizes, while being significantly faster than other published algorithms. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

Ergodicity and model quality in template‐restrained canonical and temperature/Hamiltonian replica exchange coarse‐grained molecular dynamics simulations of proteins

Molecular simulations restrained to single or multiple templates are commonly used in protein‐structure modeling. However, the restraints introduce additional barriers, thus impairing the ergodicity of simulations, which can affect the quality of the resulting models. In this work, the effect of restraint types and simulation schemes on ergodicity and model quality was investigated by performing template‐restrained canonical molecular dynamics (MD), multiplexed replica‐exchange molecular dynamics, and Hamiltonian replica exchange molecular dynamics (HREMD) simulations with the coarse‐grained UNRES force field on nine selected proteins, with pseudo‐harmonic log‐Gaussian (unbounded) or Lorentzian (bounded) restraint functions. The best ergodicity was exhibited by HREMD. It has been found that non‐ergodicity does not affect model quality if good templates are used to generate restraints. However, when poor‐quality restraints not covering the entire protein are used, the improved ergodicity of HREMD can lead to significantly improved protein models. © 2017 Wiley Periodicals, Inc.     

Coulomb replica‐exchange method: Handling electrostatic attractive and repulsive forces for biomolecules

We propose a new type of the Hamiltonian replica‐exchange method (REM) for molecular dynamics (MD) and Monte Carlo simulations, which we refer to as the Coulomb REM (CREM). In this method, electrostatic charge parameters in the Coulomb interactions are exchanged among replicas while temperatures are exchanged in the usual REM. By varying the atom charges, the CREM overcomes free‐energy barriers and realizes more efficient sampling in the conformational space than the REM. Furthermore, this method requires only a smaller number of replicas because only the atom charges of solute molecules are used as exchanged parameters. We performed Coulomb replica‐exchange MD simulations of an alanine dipeptide in explicit water solvent and compared the results with those of the conventional canonical, replica exchange, and van der Waals REMs. Two force fields of AMBER parm99 and AMBER parm99SB were used. As a result, the CREM sampled all local‐minimum free‐energy states more frequently than the other methods for both force fields. Moreover, the Coulomb, van der Waals, and usual REMs were applied to a fragment of an amyloid‐β peptide (Aβ) in explicit water solvent to compare the sampling efficiency of these methods for a larger system. The CREM sampled structures of the Aβ fragment more efficiently than the other methods. We obtained β‐helix, α‐helix, 3₁₀‐helix, β‐hairpin, and β‐sheet structures as stable structures and deduced pathways of conformational transitions among these structures from a free‐energy landscape. © 2012 Wiley Periodicals, Inc.     

A single‐molecule reaction cascade: First‐principles molecular dynamics simulation

The success of mechanochemistry is continued with the targeted organic synthesis of functional nano‐scale devices. In the present theoretical study, first principles molecular dynamics simulations are performed for a recently synthesized three‐ring system that functions as a molecular reaction cascade. Mechanochemical and photochemical reaction conditions are investigated. The system was designed in a way that three bonds would break consecutively in mechanochemical and sonochemical setups. We succeeded to simulate the reaction mechanisms with first‐principles molecular dynamics simulations and discuss the stereochemistry.     

A multiple time step algorithm compatible with a large number of distance classes and an arbitrary distance dependence of the time step size for the fast evaluation of nonbonded interactions in molecular simulations

A new algorithm is introduced to perform the multiple time step integration of the equations of motion for a molecular system, based on the splitting of the nonbonded interactions into a series of distance classes. The interactions between particle pairs in successive classes are updated at a progressively decreasing frequency. Unlike previous multiple time-stepping schemes relying on distance classes, the present algorithm sorts interacting particle pairs by their next update times rather than by their update frequencies. For this reason, the proposed scheme is extremely flexible with respect to the number of classes that can be employed (up to hundred or more) and the distance dependence of the relative time step size (arbitrary integer function of the distance). It can also easily be adapted to classes defined based on a criterion other than the interparticle distance (e.g., interaction magnitude). Different variants of the algorithm are tested in terms of accuracy and efficiency for simulations of a pure water system (6167 molecules) under truncated-octahedral periodic boundary conditions, and compared to the twin-range method standardly used with GROMOS96 (short- and long-range cutoff distances of 0.8 and 1.4 nm, pair list and intermediate-range interactions updated every five steps). In particular, multiple time-stepping schemes with an accuracy comparable to that of the twin-range method can be designed, that permit to increase the effective (long-range) cutoff distance from 1.4 to 3.0 nm with a performance loss of only about a factor 2. This result is quite encouraging, considering the benefits of doubling the cutoff radius in the context of (bio-)molecular simulations.     

GridMAT‐MD: A grid‐based membrane analysis tool for use with molecular dynamics

GridMAT‐MD is a new program developed to aid in the analysis of lipid bilayers from molecular dynamics simulations. It reads a GROMACS coordinate file and generates two types of data: a two‐dimensional contour plot depicting membrane thickness, and a polygon‐based tessellation of the individual lipid headgroups. GridMAT‐MD can also account for proteins or small molecules within the headgroups of the lipids, closely approximating their occupied lateral area. The program requires no installation, is fast, and is freely available. © 2008 Wiley Periodicals, Inc. J Comput Chem, 2009     

Pyrite: A blender plugin for visualizing molecular dynamics simulations using industry‐standard rendering techniques

Molecular dynamics (MD) simulations provide critical insights into many biological mechanisms. Programs such as VMD, Chimera, and PyMOL can produce impressive simulation visualizations, but they lack many advanced rendering algorithms common in the film and video‐game industries. In contrast, the modeling program Blender includes such algorithms but cannot import MD‐simulation data. MD trajectories often require many gigabytes of memory/disk space, complicating Blender import. We present Pyrite, a Blender plugin that overcomes these limitations. Pyrite allows researchers to visualize MD simulations within Blender, with full access to Blender's cutting‐edge rendering techniques. We expect Pyrite‐generated images to appeal to students and non‐specialists alike. A copy of the plugin is available at http://durrantlab.com/pyrite/ , released under the terms of the GNU General Public License Version 3. © 2017 Wiley Periodicals, Inc.     

Chemical reaction dynamics of PeCB and TCDD decomposition: A tight‐binding quantum chemical molecular dynamics study with first‐principles parameterization

The decomposition reaction dynamics of 2,3,4,4′,5‐penta‐chlorinated biphenyl (2,3,4,4′,5‐PeCB), 3,3′,4,4′,5‐penta‐chlorinated biphenyl (3,3′,4,4′,5‐PeCB), and 2,3,7,8‐tetra‐chlorinated dibenzo‐p‐dioxin (2,3,7,8‐TCDD) was clarified for the first time at atomic and electronic levels, using our novel tight‐binding quantum chemical molecular dynamics method with first‐principles parameterization. The calculation speed of our new method is over 5000 times faster than that of the conventional first‐principles molecular dynamics method. We confirmed that the structure, energy, and electronic states of the above molecules calculated by our new method are quantitatively consistent with those by first‐principles calculations. After the confirmation of our methodology, we investigated the decomposition reaction dynamics of the above molecules and the calculated dynamic behaviors indicate that the oxidation of the 2,3,4,4′,5‐PeCB, 3,3′,4,4′,5‐PeCB, and 2,3,7,8‐TCDD proceeds through an epoxide intermediate, which is in good agreement with the previous experimental reports and consistent with our static density functional theory calculations. These results proved that our new tight‐binding quantum chemical molecular dynamics method with first‐principles parameterization is an effective tool to clarify the chemical reaction dynamics at reaction temperatures. © 2004 Wiley Periodicals, Inc. Int J Quantum Chem, 2005     

ORAC: A Molecular dynamics program to simulate complex molecular systems with realistic electrostatic interactions

In this study, we present a new molecular dynamics program for simulation of complex molecular systems. The program, named ORAC, combines state-of-the-art molecular dynamics (MD) algorithms with flexibility in handling different types and sizes of molecules. ORAC is intended for simulations of molecular systems and is specifically designed to treat biomolecules efficiently and effectively in solution or in a crystalline environment. Among its unique features are: (i) implementation of reversible and symplectic multiple time step algorithms (or r-RESPA, reversible reference system propagation algorithm) specifically designed and tuned for biological systems with periodic boundary conditions; (ii) availability for simulations with multiple or single time steps of standard Ewald or smooth particle mesh Ewald (SPME) for computation of electrostatic interactions; and (iii) possibility of simulating molecular systems in a variety of thermodynamic ensembles. We believe that the combination of these algorithms makes ORAC more advanced than other MD programs using standard simulation algorithms. © 1997 John Wiley & Sons, Inc. J Comput Chem 18 : 1848–1862, 1997     

Three pillars for achieving quantum mechanical molecular dynamics simulations of huge systems: Divide‐and‐conquer,density‐functional tight‐binding,and massively parallel computation

The linear‐scaling divide‐and‐conquer (DC) quantum chemical methodology is applied to the density‐functional tight‐binding (DFTB) theory to develop a massively parallel program that achieves on‐the‐fly molecular reaction dynamics simulations of huge systems from scratch. The functions to perform large scale geometry optimization and molecular dynamics with DC‐DFTB potential energy surface are implemented to the program called DC‐DFTB‐K. A novel interpolation‐based algorithm is developed for parallelizing the determination of the Fermi level in the DC method. The performance of the DC‐DFTB‐K program is assessed using a laboratory computer and the K computer. Numerical tests show the high efficiency of the DC‐DFTB‐K program, a single‐point energy gradient calculation of a one‐million‐atom system is completed within 60 s using 7290 nodes of the K computer. © 2016 Wiley Periodicals, Inc.