Computer Physics Communications 182 (2011) 2307–2313 Contents lists available at ScienceDirect Computer Physics Communications www.elsevier.com/locate/cpc Rigid body constraints realized in massively-parallel molecular dynamics on graphics processing units Trung Dac Nguyen a , Carolyn L. Phillips b , Joshua A. Anderson a , Sharon C. Glotzer a,c,∗ a Department of Chemical Engineering, University of Michigan, Ann Arbor, MI 48109, USA b Applied Physics Program, University of Michigan, Ann Arbor, MI 48109, USA c Department of Materials Science and Engineering, University of Michigan, Ann Arbor, MI 48109, USA article info abstract Article history: Received 25 February 2011 Accepted 11 June 2011 Available online 14 June 2011 Keywords: GPU GPGPU CUDA Molecular dynamics Rigid body Molecular dynamics (MD) methods compute the trajectory of a system of point particles in response to a potential function by numerically integrating Newton’s equations of motion. Extending these basic methods with rigid body constraints enables composite particles with complex shapes such as anisotropic nanoparticles, grains, molecules, and rigid proteins to be modeled. Rigid body constraints are added to the GPU-accelerated MD package, HOOMD-blue, version 0.10.0. The software can now simulate systems of particles, rigid bodies, or mixed systems in microcanonical (NVE), canonical (NVT), and isothermal- isobaric (NPT) ensembles. It can also apply the FIRE energy minimization technique to these systems. In this paper, we detail the massively parallel scheme that implements these algorithms and discuss how our design is tuned for the maximum possible performance. Two different case studies are included to demonstrate the performance attained, patchy spheres and tethered nanorods. In typical cases, HOOMD- blue on a single GTX 480 executes 2.5–3.6 times faster than LAMMPS executing the same simulation on any number of CPU cores in parallel. Simulations with rigid bodies may now be run with larger systems and for longer time scales on a single workstation than was previously even possible on large clusters.  2011 Elsevier B.V. All rights reserved. 1. Introduction Molecular dynamics (MD) simulations and related methods are powerful tools for modeling systems of particles [1]. The basic MD technique computes the trajectory of n particles under the influ- ence of a potential V ( r 1 ,  r 2 ,...,  r n ), the negative gradient of which gives a conservative force  F =−  ∇ V , by integrating Newton’s equa- tions of motion over discrete time steps that each advance the state of the system from [ r i (t ),  p i (t )] to [ r i (t + t ),  p i (t + t )]. The quantities  r i and  p i are the position and momentum of the i -th particle, respectively, t is the current simulation time, and t is the step size. Many applications of MD, such as soft matter self- assembly [2–4] and protein folding [5–7], necessitate running hun- dreds of millions of time steps per run and thousands of individual runs. Accelerating the rate at which time steps are performed re- duces the time to discovery and enables better predictions through the use of higher fidelity models. Classical MD breaks the potential into pair-wise and bond terms V = ∑ pairs i, j V p (r ij ) + ∑ bonds i, j V b (r ij ). Smooth, “soft” potentials V p (r ) and V b (r ) can be used in conjunction with a large step size. On the other hand, steep, “hard” potentials, such as bonds with * Corresponding author at: Department of Chemical Engineering, University of Michigan, Ann Arbor, MI 48109, USA. E-mail address: sglotzer@umich.edu (S.C. Glotzer). stiff spring constants, require using a prohibitively small step size to maintain accuracy and stability. Potentials with infinitely steep interaction terms can only be achieved with extensions to the basic MD framework. One such extension is SHAKE [8]. The SHAKE algorithm en- forces fixed bond distances between two particles. Via an iterative method, any number of bonds in the system can be constrained. A set of particles may be combined into a single rigid body with an appropriate choice of bond constraints while taking special care not to over-constrain the system. However, certain rigid shapes, such as planar and linear molecules, cannot be created in three di- mensions by setting bond distances alone because the constraint matrix is singular. Although the SHAKE algorithm has been ex- tended to handle arbitrary shapes, for example, via angle and dihedral angle constraints, [9,10] the computational cost of these algorithms often becomes prohibitive for parallel simulation codes as the number of constraints per cluster increases. Modeling large or generic rigid arrangements of particles can also be achieved by treating each defined set of particles as a sin- gle rigid body with only three translational and three orientational degrees of freedom (or two and one, respectively, for 2D simula- tions) [11]. Such a method can be added to an MD package with minimal modifications by taking advantage of the existing code that computes particle–particle interactions. Rigid body constraints are available in MD software packages such as DLPOLY [12] and 0010-4655/$ – see front matter  2011 Elsevier B.V. All rights reserved. doi:10.1016/j.cpc.2011.06.005