Berk Hess

GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation

Molecular simulation is an extremely useful, but computationally very expensive tool for studies of chemical and biomolecular systems. Here, we present a new implementation of our molecular simulation toolkit GROMACS which now both achieves extremely high performance on single processors from algorithmic optimizations and hand-coded routines and simultaneously scales very well on parallel machines. The code encompasses a minimal-communication domain decomposition algorithm, full dynamic load balancing, a state-of-the-art parallel constraint solver, and efficient virtual site algorithms that allow removal of hydrogen atom degrees of freedom to enable integration time steps up to 5 fs for atomistic simulations also in parallel. To improve the scaling properties of the common particle mesh Ewald electrostatics algorithms, we have in addition used a Multiple-Program, Multiple-Data approach, with separate node domains responsible for direct and reciprocal space interactions. Not only does this combination of algorithms enable extremely long simulations of large systems but also it provides that simulation performance on quite modest numbers of standard cluster nodes.

GROMACS: Fast, flexible, and free

This article describes the software suite GROMACS (Groningen MAchine for Chemical Simulation) that was developed at the University of Groningen, The Netherlands, in the early 1990s. The software, written in ANSI C, originates from a parallel hardware project, and is well suited for parallelization on processor clusters. By careful optimization of neighbor searching and of inner loop performance, GROMACS is a very fast program for molecular dynamics simulation. It does not have a force field of its own, but is compatible with GROMOS, OPLS, AMBER, and ENCAD force fields. In addition, it can handle polarizable shell models and flexible constraints. The program is versatile, as force routines can be added by the user, tabulated functions can be specified, and analyses can be easily customized. Nonequilibrium dynamics and free energy determinations are incorporated. Interfaces with popular quantum‐chemical packages (MOPAC, GAMES‐UK, GAUSSIAN) are provided to perform mixed MM/QM simulations. The package includes about 100 utility and analysis programs. GROMACS is in the public domain and distributed (with source code and documentation) under the GNU General Public License. It is maintained by a group of developers from the Universities of Groningen, Uppsala, and Stockholm, and the Max Planck Institute for Polymer Research in Mainz. Its Web site is http://www.gromacs.org.

LINCS : A linear constraint solver for molecular simulations

In this article, we present a new LINear Constraint Solver (LINCS) for molecular simulations with bond constraints. The algorithm is inherently stable, as the constraints themselves are reset instead of derivatives of the constraints, thereby eliminating drift. Although the derivation of the algorithm is presented in terms of matrices, no matrix matrix multiplications are needed and only the nonzero matrix elements have to be stored, making the method useful for very large molecules. At the same accuracy, the LINCS algorithm is three to four times faster than the SHAKE algorithm. Parallelization of the algorithm is straightforward. © 1997 John Wiley & Sons, Inc. J Comput Chem 18: 1463–1472, 1997

GROMACS 4.5

MOTIVATION Molecular simulation has historically been a low-throughput technique, but faster computers and increasing amounts of genomic and structural data are changing this by enabling large-scale automated simulation of, for instance, many conformers or mutants of biomolecules with or without a range of ligands. At the same time, advances in performance and scaling now make it possible to model complex biomolecular interaction and function in a manner directly testable by experiment. These applications share a need for fast and efficient software that can be deployed on massive scale in clusters, web servers, distributed computing or cloud resources. RESULTS Here, we present a range of new simulation algorithms and features developed during the past 4 years, leading up to the GROMACS 4.5 software package. The software now automatically handles wide classes of biomolecules, such as proteins, nucleic acids and lipids, and comes with all commonly used force fields for these molecules built-in. GROMACS supports several implicit solvent models, as well as new free-energy algorithms, and the software now uses multithreading for efficient parallelization even on low-end systems, including windows-based workstations. Together with hand-tuned assembly kernels and state-of-the-art parallelization, this provides extremely high performance and cost efficiency for high-throughput as well as massively parallel simulations. AVAILABILITY GROMACS is an open source and free software available from http://www.gromacs.org. SUPPLEMENTARY INFORMATION Supplementary data are available at Bioinformatics online.

P-LINCS: A Parallel Linear Constraint Solver for Molecular Simulation

By removing the fastest degrees of freedom, constraints allow for an increase of the time step in molecular simulations. In the last decade parallel simulations have become commonplace. However, up till now efficient parallel constraint algorithms have not been used with domain decomposition. In this paper the parallel linear constraint solver (P-LINCS) is presented, which allows the constraining of all bonds in macromolecules. Additionally the energy conservation properties of (P-)LINCS are assessed in view of improvements in the accuracy of uncoupled angle constraints and integration in single precision.

Improving efficiency of large time‐scale molecular dynamics simulations of hydrogen‐rich systems

A systematic analysis is performed on the effectiveness of removing degrees of freedom from hydrogen atoms and/or increasing hydrogen masses to increase the efficiency of molecular dynamics simulations of hydrogen‐rich systems such as proteins in water. In proteins, high‐frequency bond‐angle vibrations involving hydrogen atoms limit the time step to 3 fs, which is already a factor of 1.5 beyond the commonly used time step of 2 fs. Removing these degrees of freedom from the system by constructing hydrogen atoms as dummy atoms, allows the time step to be increased to 7 fs, a factor of 3.5 compared with 2 fs. Additionally, a gain in simulation stability can be achieved by increasing the masses of hydrogen atoms with remaining degrees of freedom from 1 to 4 u. Increasing hydrogen mass without removing the high‐frequency degrees of freedom allows the time step to be increased only to 4 fs, a factor of two, compared with 2 fs. The net gain in efficiency of sampling configurational space may be up to 15% lower than expected from the increase in time step due to the increase in viscosity and decrease in diffusion constant. In principle, introducing dummy atoms and increasing hydrogen mass do not influence thermodynamical properties of the system and dynamical properties are shown to be influenced only to a moderate degree. Comparing the maximum time step attainable with these methods (7 fs) to the time step of 2 fs that is routinely used in simulation, and taking into account the increase in viscosity and decrease in diffusion constant, we can say that a net gain in simulation efficiency of a factor of 3 to 3.5 can be achieved. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 786–798, 1999

Implementation of the CHARMM Force Field in GROMACS : Analysis of Protein Stability Effects from Correction Maps, Virtual Interaction Sites, and Water Models

CHARMM27 is a widespread and popular force field for biomolecular simulation, and several recent algorithms such as implicit solvent models have been developed specifically for it. We have here implemented the CHARMM force field and all necessary extended functional forms in the GROMACS molecular simulation package, to make CHARMM-specific features available and to test them in combination with techniques for extended time steps, to make all major force fields available for comparison studies in GROMACS, and to test various solvent model optimizations, in particular the effect of Lennard-Jones interactions on hydrogens. The implementation has full support both for CHARMM-specific features such as multiple potentials over the same dihedral angle and the grid-based energy correction map on the ϕ, ψ protein backbone dihedrals, as well as all GROMACS features such as virtual hydrogen interaction sites that enable 5 fs time steps. The medium-to-long time effects of both the correction maps and virtual sites have been tested by performing a series of 100 ns simulations using different models for water representation, including comparisons between CHARMM and traditional TIP3P. Including the correction maps improves sampling of near native-state conformations in our systems, and to some extent it is even able to refine distorted protein conformations. Finally, we show that this accuracy is largely maintained with a new implicit solvent implementation that works with virtual interaction sites, which enables performance in excess of 250 ns/day for a 900-atom protein on a quad-core desktop computer.

Determining the shear viscosity of model liquids from molecular dynamics simulations

Several methods are available for calculating shear viscosities of liquids from molecular dynamics simulations. There are equilibrium methods based on pressure or momentum fluctuations and several nonequilibrium methods. For the nonequilibrium method using a periodic shear flow, all relevant quantities, including the accuracy, can be estimated before performing the simulation. We compared the applicability, accuracy and efficiency of this method with two equilibrium methods and another nonequilibrium method, using simulations of a Lennard-Jones fluid and the SPC and SPC/E [(extended) simple point charge] water models.

Cation specific binding with protein surface charges

Biological organization depends on a sensitive balance of noncovalent interactions, in particular also those involving interactions between ions. Ion-pairing is qualitatively described by the law of “matching water affinities.” This law predicts that cations and anions (with equal valence) form stable contact ion pairs if their sizes match. We show that this simple physical model fails to describe the interaction of cations with (molecular) anions of weak carboxylic acids, which are present on the surfaces of many intra- and extracellular proteins. We performed molecular simulations with quantitatively accurate models and observed that the order K+ < Na+ < Li+ of increasing binding affinity with carboxylate ions is caused by a stronger preference for forming weak solvent-shared ion pairs. The relative insignificance of contact pair interactions with protein surfaces indicates that thermodynamic stability and interactions between proteins in alkali salt solutions is governed by interactions mediated through hydration water molecules.

A flexible algorithm for calculating pair interactions on SIMD architectures

Calculating interactions or correlations between pairs of particles is typically the most time-consuming task in particle simulation or correlation analysis. Straightforward implementations using a double loop over particle pairs have traditionally worked well, especially since compilers usually do a good job of unrolling the inner loop. In order to reach high performance on modern CPU and accelerator architectures, single-instruction multiple-data (SIMD) parallelization has become essential. Avoiding memory bottlenecks is also increasingly important and requires reducing the ratio of memory to arithmetic operations. Moreover, when pairs only interact within a certain cut-off distance, good SIMD utilization can only be achieved by reordering input and output data, which quickly becomes a limiting factor. Here we present an algorithm for SIMD parallelization based on grouping a fixed number of particles, e.g. 2, 4, or 8, into spatial clusters. Calculating all interactions between particles in a pair of such clusters improves data reuse compared to the traditional scheme and results in a more efficient SIMD parallelization. Adjusting the cluster size allows the algorithm to map to SIMD units of various widths. This flexibility not only enables fast and efficient implementation on current CPUs and accelerator architectures like GPUs or Intel MIC, but it also makes the algorithm future-proof. We present the algorithm with an application to molecular dynamics simulations, where we can also make use of the effective buffering the method introduces.