Scalable molecular dynamics on CPU and GPU architectures with NAMD

Abstract

NAMDis a molecular dynamics program designed for high-performance simulations of very large biological objects on CPU- and GPU-based architectures. NAMD offers scalable performance on petascale parallel supercomputers consisting of hundreds of thousands of cores, as well as on inexpensive commodity clusters commonly found in academic environments. It is written in C++ and leans on Charm++ parallel objects for optimal performance on low-latency architectures. NAMD is a versatile, multipurpose code that gathers state-of-the-art algorithms to carry out simulations in apt thermodynamic ensembles, using the widely popular CHARMM, AMBER, OPLS, and GROMOS biomolecular force fields. Here, we review the main features of NAMD that allow both equilibrium and enhanced-sampling molecular dynamics simulations with numerical efficiency. We describe the underlying concepts utilized by NAMD and their implementation, most notably for handling long-range electrostatics; controlling the temperature, pressure, and pH; applying external potentials on tailored grids; leveraging massively parallel resources in multiple-copy simulations; and hybrid quantum-mechanical/molecular-mechanical descriptions. We detail the variety of options offered by NAMD for enhanced-sampling simulations aimed at determining free-energy differences of either alchemical or geometrical transformations and outline their applicability to specific problems. Last, we discuss the roadmap for the development of NAMD and our current efforts toward achieving optimal performance on GPU-based architectures, for pushing back the limitations that have prevented biologically realistic billion-atom objects to be fruitfully simulated, and for making large-scale simulations less expensive and easier to set up, run, and analyze. NAMD is distributed free of charge with its source code at www.ks.uiuc.edu.

FIG. 1.

Workflow architecture of a NAMD simulation. The blue boxes denote the core physical details of the simulation, the orange boxes denote the features that are embedded in the code (written in C++), and the green boxes denote the features that are fully implemented through user scripts (written in Tcl or Python).

FIG. 2.

Overview of a Charm++ application on four processor elements (PEs). The solid line indicates an object on PE 3 sending a message to an object on PE 1, and the dashed line shows how the message is sent via per-PE schedulers.

FIG. 3.

NAMD software architecture. “Computes” are objects responsible for force calculations, operating on data received from “patches,” small boxes generated by the spatial decomposition of the computational assay. Long-range electrostatic forces are handled by the PME algorithm.

FIG. 4.

NAMD distributed memory parallel scaling benchmarks performed on satellite tobacco mosaic virus (STMV) in water on leadership supercomputers, Summit, using up to 1024 nodes containing 6144 GPUs (upper panel), and Frontera, using up to 512 CPU-based nodes (lower panel). Two benchmark systems were constructed as a tiled array of a periodic 1.067M-atom system, a 5 × 2 × 2 tiling for 21M atoms (▴), and a 7 × 6 × 5 tiling for 224M atoms (▾).

FIG. 5.

Standard GPU offload approach compared against new GPU-resident execution scheme for a single-node NAMD simulation of apolipoprotein 1 (ApoA1) in water, consisting of 92 224 atoms. The light blue line tracks GPU activity, while the black strip tracks CPU activity. GPU force calculations are labeled “force,” and GPU integration calculations are labeled “int.”

FIG. 6.

Comparison of NAMD GPU computing schemes for simulations with increasing atom counts, namely, (i) reductase enzyme (DHFR) with 23 558 atoms, (ii) apolipoprotein 1 (ApoA1) with 92 224 atoms, (iii) F1-ATPase with 327 506 atoms, and (iv) STMV with 1 066 628 atoms for the GPU-resident scheme in NAMD 3.0 (▴) and for standard GPU offloading scheme in NAMD 2.14 (▾). The simulations were performed in the microcanonical ensemble, using the CHARMM (turquoise line) and AMBER-like (magenta line) cutoff schemes: 12 Å and 8 Å, respectively.

FIG. 7.

NAMD performance compared on two platforms, NVIDIA DGX-2 (upper panel) and Amazon Web Services (AWS) P3.16xlarge (lower panel), using the standard GPU offloading scheme in NAMD 2.14 and the GPU-resident scheme in NAMD 3.0. The DGX-2 and AWS P3.16xlarge platforms consisted, respectively, of 16× NVIDIA Tesla V100 GPUs and 8× NVIDIA Tesla V100 GPUs associated with a 48-core CPU. The benchmarks were conducted on replica simulations of apolipoprotein 1 (ApoA1) in water, representing 92 224 atoms, in the microcanonical ensemble, with one independent replica running on each GPU. The benchmarks were carried out with the CHARMM (turquoise bars) and AMBER-like (magenta bars) cutoff schemes: 12 Å and 8 Å, respectively.

FIG. 8.

NAMD scaling on the Summit supercomputer shows the performance advantage of using stochastic velocity rescaling (blue line) over the Langevin damping thermostat (red line), up to 20% faster for the same number of nodes.

FIG. 9.

Graphical representation of a Colvars configuration.

FIG. 10.

Gradients of the Distance to Bound Configuration (DBC) coordinate for a phenol molecule bound to T4 lysozyme. The reference pose of phenol heavy atoms is shown in light cyan. The arrows are proportional to the derivatives of the CV with respect to atomic Cartesian coordinates. The small gradient contributions on protein C_α carbons illustrate the purely roto-translational counter-forces exerted on the receptor when biasing forces are applied to a DBC coordinate.

FIG. 11.

CVs for solvent reorganization based on volumetric maps. Snapshot of a hydrophobic cavity containing a varying number of water molecules (a). A volumetric map (red transparent contour) is used to define a continuous variable $N_w$ measuring the number of molecules in the cavity. Trajectory of a short SMD simulation on $N_w$ (red), compared to the number of molecules counted with VMD (blue) (b). PMF of $N_w$, computed over ∼400 ns. The use of gridForces maps (see Sec. XIII) in Colvars allows the development of methods for enhanced sampling of fully dynamic aggregates, such as water densities and lipid membranes (c).

FIG. 12.

Sampling of a rugged free-energy landscape. Boltzmann sampling favors low-energy regions (a). Free-energy barriers can be overcome by depositing harmonic potentials, as is done in US (b); by applying a time-dependent bias that yields a Hamiltonian bereft of an average force along ξ(x), as is done in ABF, or its extended-Lagrangian variant, eABF (c); or by flooding valleys by means of Gaussian potentials, or hills, as is done in MtD (d). Combination of MtD and eABF, meta-eABF, concurrently shaves free-energy barriers and floods valley (e). To enhance ergodic sampling, a multiple-walker extension of MtD (f), ABF or eABF (g), and meta-eABF (h) has been implemented.

FIG. 13.

Potential of mean force characterizing the permeation of a 1–palmitoyl–2–oleoyl-phosphatidylcholine:cholesterol model membrane by 2′, 3′–dideoxyadenosine, obtained using ABF (red line) and meta-eABF (blue line) (a). The reaction-coordinate model is the projection onto the z-axis of the distance separating the center of mass of the permeant from that of the lipid bilayer (b).

FIG. 14.

Single- and dual-topology frameworks for alchemical transformations in NAMD. When the reference and the target states are chemically different and their respective charge distributions are distinct, a minimum common substructure is sought in the dual-topology paradigm [(a) and (c)]. The two topologies coexist in a hybrid compound, albeit do not interact, either directly (through nonbonded contributions) or indirectly (through bonded contributions). Unchanged (black), incoming (blue), and outgoing (magenta) atoms are defined in three partitions. In some extreme cases, when the end states are chemically distinct and no common substructure can be found, the two topologies are introduced independently, and geometric restraints are enforced to prevent them from drifting apart from each other. In the single-topology paradigm [(b) and (d)], holonomic constraints are applied to the common substructure (black), and dummy particles are introduced to describe the alternate state. As the transformation proceeds, these dummy particles become real, interacting atoms. To avoid the creation of an unphysical net force, the alternate state is only partially chemically bonded to the rest of the molecule.

FIG. 15.

Hamiltonian-exchange FEP calculations or FEP/($\lambda ,\mathcal{H}$)–REMD (see Sec. X). The computational assay is cloned into replicas corresponding to the different strata of the alchemical transformation whereby the interaction of a small molecule with its environment is progressively turned off. The configurations of adjacent λ-states, i and j, are swapped periodically, and the proposed move is accepted following a Metropolis–Hastings criterion, $p\left(\right.{\lambda }_i\to {\lambda }_j=min\left\{1;{\mathrm{e}}^{-\beta \left\{\left[U_j\left(\mathrm{x},{\lambda }_j\right)-U_i\left(\mathrm{x},{\lambda }_j\right)\right]+\left[U_i\left(\mathrm{x},{\lambda }_i\right)-U_j\left(\mathrm{x},{\lambda }_i\right)\right]\right\}}\right\}$.

FIG. 16.

Generic implementation of MCA in the Charm++ RTS. Multiple NAMD instances run concurrently, and each of them is executed by an independent Charm++ RTS (gray vertical lines). Replica exchange/inter-copy communications (red double arrows) are implemented through a Tcl scripting interface. comm_local denotes a local internal Charm++ partition, dedicated to be the communication layer of a single trajectory run/NAMD instance. comm_cross denotes a global Charm++ communication layer across all local partitions. A scripting interface of replica exchange is implemented with communication-enabled Tcl functions built on top of comm_cross.

FIG. 17.

Sample titration curves computed from a constant-pH MD trajectory of the BBL miniprotein reveal a range of protonation states over a modest pH range. Solid lines represent fits to a Hill curve, from which pK_a values and Hill coefficients, n, are derived. The present results were obtained using the CHARMM36 force field and standard settings (e.g., rigid bonds with hydrogen, PME). An MC step was performed every ps, and a uniform switch time of 15 ps was used for each move (particular residues requiring longer/shorter times can also be specified). About 20 000 moves were attempted at seven pH values, leading to an aggregate of about 2 µs of MD.

FIG. 18.

Schematic representation of a hybrid QM/MM MD step in NAMD. Two trialanines, each containing one QM region, are simulated through simultaneous and independent executions of the chosen QM software (ORCA). In a QM/MM simulation, positions and elements of QM atoms are sent by NAMD to the QM software, along with positions and partial charges representing the local MM environment within a given cutoff (dashed magenta or purple lines, top-left panel). NAMD expects from the QM software forces, total energy, and partial charges for all QM atoms, and forces acting on MM partial charges due to electrostatic interactions. All atoms are moved based on the calculated force gradients. Different atoms will have their resulting force gradient computed either by the QM software, by NAMD, or a combination thereof (top-right panel). If selected, long-range electrostatics can be calculated by NAMD for all atoms using PME, utilizing the updated partial charges calculated by the QM software. For additional details, the reader is referred to Ref. .

FIG. 19.

Using grid-based potentials to steer an MD simulation. (a) A cut-away view of an all-atom simulation system featuring a single DNA strand (orange) passing through an α-hemolysin nanopore (blue) embedded in a lipid membrane (green). Water and ions are not shown for clarity. (b) Electrostatic potential contour map in and around an α-hemolysin nanopore determined using all-atom MD methods. Coupling this potential to atoms of ssDNA only through the gridForces module of NAMD dramatically accelerates the simulation of the DNA translocation process. (c) The cryo-EM density (white wireframe surface) is used as a steering potential for docking the all-atom structure of an α/β (red/blue) tubulin dimer to build an all-atom model of a microtubule. (d) All-atom MD simulation of plasmonic trapping. Electric field-driven translocation of DNA (blue and orange) through a nanopore in a solid-state membrane (gray) is halted by the plasmonic field (semi-transparent red surface) generated by the gold bowtie nanostructure (yellow). The arrows illustrate the magnitude and direction of plasmonic forces acting on the DNA molecules. Water and ions are not shown for clarity.

FIG. 20.

A simple MDFF example of fitting a protein structure into a density map. (a) The x-ray structure of an adenylate kinase protein (PDB:1AKE). (b) The protein structure rigid body docked to the density map. (c) The protein structure inside the density map after running an MDFF simulation. (d) The overall quality of fit (cross correlation) during an MDFF simulation.

FIG. 21.

Continuous development effort over the years toward simulating with NAMD realistic biological objects of increasing complexity from a small, solvated protein, on the thousand-atom size scale, in the early 1990s, to a full protocell, on the billion-atom size scale, nowadays.