Scalable molecular dynamics with NAMD

Abstract

NAMD is a parallel molecular dynamics code designed for high-performance simulation of large biomolecular systems. NAMD scales to hundreds of processors on high-end parallel platforms, as well as tens of processors on low-cost commodity clusters, and also runs on individual desktop and laptop computers. NAMD works with AMBER and CHARMM potential functions, parameters, and file formats. This article, directed to novices as well as experts, first introduces concepts and methods used in the NAMD program, describing the classical molecular dynamics force field, equations of motion, and integration methods along with the efficient electrostatics evaluation algorithms employed and temperature and pressure controls used. Features for steering the simulation across barriers and for calculating both alchemical and conformational free energy differences are presented. The motivations for and a roadmap to the internal design of NAMD, implemented in C++ and based on Charm++ parallel objects, are outlined. The factors affecting the serial and parallel performance of a simulation are discussed. Finally, typical NAMD use is illustrated with representative applications to a small, a medium, and a large biomolecular system, highlighting particular features of NAMD, for example, the Tcl scripting language. The article also provides a list of the key features of NAMD and discusses the benefits of combining NAMD with the molecular graphics/sequence analysis software VMD and the grid computing/collaboratory software BioCoRE. NAMD is distributed free of charge with source code at www.ks.uiuc.edu.

Figure 1

Growth in practical simulation size illustrated by comparison of estrogen receptor DNA-binding domain simulation of [5], left, 36,000 atoms simulated over 100 ps, published 1997, with multiscale LacI-DNA complex simulation of [6], right, 314,000 atoms simulated over 10 ns, published 2005 and described in Sec. 4.3.

Figure 2

Force field bonded interactions: r governs bond stretching; θ represents the bond angle term; φ gives the dihedral angle; the out-of-plane angle α may be controlled by an “improper” dihedral ϕ.

Figure 3

In PME, a charge (denoted by an empty circle with label “q”) in the figure is distributed over grid (here a mesh in two dimensions) points with weighting functions chosen according to the distance of the respective grid points to the location of the charge. Positioning all charges on a grid enables the application of the FFT method and significantly reduces the computation time. In real applications, the grid is three-dimensional.

Figure 4

Smoothed electrostatic potential of decalanine in vacuum as calculated with the PME plugin of VMD. Atoms are colored by charge (blue is positive, red is negative). The helix dipole is clearly visible from the two potential isosurfaces +20k_BT/e (blue, left lobe) and −20k_BT/e (red, right lobe).

Figure 5

A symplectic method demonstrates superior long-time stability: the symplectic method used here is the sympletic Euler method [, Theorem 3.3], whose local error is proportional to Δt²; the non-symplectic method used is the Runge-Kutta method [, §8.3.3] whose local error is proportional to Δt⁵. The system described is a one-dimension harmonic oscillator. The equation of motion is mẍ = −mω²x, with m = ω = 1, and initial conditions x = 1, υ = 0. The exact trajectory is the unit circle. The chosen timestep is Δt = 0.5 and 10000 steps were integrated. The trajectory of the Runge-Kutta method collapses toward the center while the symplectic Euler method maintains a stable orbit.

Figure 6

Constant velocity pulling in a one-dimensional case. The dummy atom is colored red, and the SMD atom blue. As the dummy atom moves at constant velocity, the SMD atom experiences a force that depends linearly on the distance between both atoms.

Figure 7

In IMD, the user applies forces to atoms in the simulation via a force-feedback haptic device.

Figure 8

Homodimerization of the transmembrane (TM) domain of glycophorin A (GpA): (a): Contact surface of the two α-helices forming the TM domain of GpA. The strongest contacts are observed in the heptad of amino acids, Leu⁷⁵, Ile⁷⁶, Gly⁷⁹, Val⁸⁰, Gly⁸³, Val⁸⁴ and Thr⁸⁷. Residue Leu⁷⁵, which participates in the association of the two α-helices through dispersive interactions, is featured as transparent van der Waals spheres. (b): Free energy profile delineating the reversible dissociation of the two α-helices, obtained using an adaptive biasing force. The ordering parameter, ξ, corresponds to the distance separating the center of mass of the TM segments. The entire pathway was broken down into ten windows, in which up to 15 ns of MD sampling was performed, corresponding to a total simulation time of 125 ns. (c): Thermodynamic cycle utilized to perform the “alchemical transformation” of residues Leu⁷⁵ and Ile⁷⁶ into alanine, demonstrating the participation of these amino acids in the homodimerization of the two α-helices. The left vertical leg of the cycle represents the transformation in a single α-helix from the wild-type (WT) to the mutant form. The right vertical leg denotes a simultaneous point mutation in the two interacting α-helices. Using the notation of the figure, the free energy difference arising from this perturbation is equal to $\mathrm{\Delta }{G}_2^{\text{mutation}}-2\mathrm{\Delta }{G}_1^{\text{mutation}}$. Each leg of the thermodynamic cycle consists of 50 intermediate states and a total MD sampling of 6 ns. For clarity, the environment of the α-helical dimer, formed by a dodecane slab in equilibrium between two lamellae of water, is not shown.

Figure 9

NAMD performance on modern platforms. The cost in CPU seconds of a single timestep for a representative system of 92,000 atoms with a 12 Å short-range cutoff and full electrostatics evaluated via PME every four steps is plotted on the vertical axis, while the wait for a 1 ns simulation (one million timesteps of 1 fs each) may be read relative to the broken diagonal lines. Perfect parallel scaling is a horizontal line. The older PSC Lemieux Alpha cluster has lower serial performance but scales very well. The HPCx IBM POWER 4 cluster at Edinburgh scales similarly, with comparable serial performance to the Xeon, Opteron, and Apple Xserve G5 clusters at Urbana. The two NCSA Itanium platforms have excellent serial performance but lose efficiency beyond 128 processors.

Figure 10

Pulling ubiquitin and tetra-ubiquitin. (A)–(D) mono-ubiquitin is shown at various stages of a constant velocity SMD simulation with the N and C termini highlighted in green. The different structures are seen to unfold at different moments in the simulation. (A′)–(C′) tetra-ubiquitin is shown being pulled with constant force. Each color represents a different monomer and the transparent portion is the surface of the protein with the first and fourth segments not completely shown. Again, the different subunits can be seen to come apart at different times.

Figure 11

Side view of a simulated aquaporin tetramer in the cell membrane. The E. coli water/glycerol channel GlpF is embedded in a patch of POPE lipid bilayer and fully hydrated by water on both sides. Lipid head groups are shown in CPK and the hydrophobic tail region is drawn using licorice representation. The four AQP monomers, each forming an independent water pore, are shown in different colors. The single file of water formed inside the pores is shown in one of the monomers. The characteristic conformational inversion of water at the center of the channel that contributes to the barrier against proton transfer is discernible.

Figure 12

Snapshots taken over the course of the LacI-DNA complex multiscale simulation: (a) The evolution of the structure of the DNA loop; (b) The structure of LacI remains unchanged, with the exception of the rotation of the head groups, which allows the DNA loop to adopt a more relaxed configuration.