N-Body Simulations on GPUs
Erich Elsen
Vijay Pande

V. Vishal
Mike Houston
Pat Hanrahan
Eric Darve
Stanford University

arXiv:0706.3060v1 [] 20 Jun 2007

Abstract
Commercial graphics processors (GPUs) have high compute capacity at very low cost, which makes them attractive for general purpose scientific computing. In this
paper we show how graphics processors can be used for
N -body simulations to obtain improvements in performance over current generation CPUs. We have developed a highly optimized algorithm for performing the
O(N 2 ) force calculations that constitute the major part
of stellar and molecular dynamics simulations. In some
of the calculations, we achieve sustained performance of
nearly 100 GFlops on an ATI X1900XTX. The performance on GPUs is comparable to specialized processors
such as GRAPE-6A and MDGRAPE-3, but at a fraction of the cost. Furthermore, the wide availability of
GPUs has significant implications for cluster computing
and distributed computing efforts like Folding@Home.
Keywords:
Celestial Mechanics, N-Body Simulations, Stellar Dynamics, Molecular Dynamics, Molecular Simulation, Data Parallel Computing, Stream Computing, Programmable Graphics Hardware, GPU Computing, Brook

1

Introduction

The classical N -body problem consists of obtaining the
time evolution of a system of N mass particles interacting according to a given force law. The problem arises
in several contexts, ranging from molecular scale calculations in structural biology to stellar scale research in
astrophysics.
Molecular dynamics (MD) has been successfully
used to understand how certain proteins fold and
function, which have been outstanding questions
in biology for over three decades [Snow et al. 2005;
Gomez et al. 2004].
Exciting new developments in
MD methods offer hope that such calculations

will play a significant role in future drug research [Fujitani et al. 2005]. In stellar dynamics where
experimental observations are hard, if not impossible,
theoretical calculations may often be the only way to
understand the formation and evolution of galaxies.
Analytic solutions to the equations of motion for more
than 2 particles or complicated force functions are intractable which forces one to resort to computer simulations. A typical simulation consists of a force evaluation
step, where the force law and the current configuration
of the system are used to the compute the forces on
each particle, and an update step, where the dynamical equations (usually Newton’s laws) are numerically
stepped forward in time using the computed forces. The
updated configuration is then reused to calculate forces
for the next time step and the cycle is repeated as many
times as desired.
The simplest force models are pairwise additive, that
is the force of interaction between two particles is independent of all the other particles, and the individual
forces on a particle add linearly. The force calculation
for such models is of complexity O(N 2 ). Since typical
studies involve a large number of particles (103 to 106 )
and the desired number of integration steps is usually
very large (106 to 1015 ), the computational requirements
often limit both the problem size as well as the simulation time and consequently, the useful information
that may be obtained from such simulations. Numerous
methods have been developed to deal with these issues.
For molecular simulations, it is common to reduce the
number of particles by treating the solvent molecules
as a continuum. In stellar simulations, one uses individual time stepping or tree algorithms to minimize the
number of force calculations. Despite such algorithmic
approximations and optimizations, the computational
capabilities of current hardware remain a limiting factor.
Typically N -body simulations utilize neighborlists,
tree methods or other algorithms to reduce the order of
the force calculations. Previous work [Elsen et al. 2005]
demonstrated a GPU implementation of a neighbor list
based method to compute non-bonded forces. However,
since the GPU so far outperformed the CPU, the neighborlist creation quickly became a limiting factor. Building the neighborlist on the GPU is extremely difficult
due to the lack of specific abilities (namely indirected

output) and research on computing the neighborlist on
the GPU is still in progress. Other simplistic simulations that do not need neighborlist updates have been
implemented by others [Juekuan Yang 2006]. However,
for small N, we find we can do an O(N 2 ) calculation
significantly faster on the GPU than an O(N ) method
using the CPU (or even with a combination of the GPU
and CPU). This has direct applicability to biological
simulations that use continuum models for the solvent.
We note also that in many of the reduced order methods such as tree based schemes, at some stage an O(N 2 )
calculation is performed on a subsystem of the particles,
so our method can be used to improve the performance
of such methods as well. When using GRAPE accelerator cards for tree based algorithms, the host processor
takes care of building the tree and the accelerator cards
are used to speed up the force calculation step; GPUs
could be used in a similar way in place of the GRAPE
accelerator boards.
Using the methods described below, we are able to
accelerate the force calculation on GPUs over 25 times
compared to highly optimized SSE code running on an
Intel Pentium 4. This performance is in the range of the
specially designed GRAPE-6A [Fukushige et al. 2005]
and MDGRAPE-3 [Taiji et al. 2003] processors, but
uses a commodity processor at a much better performance/cost ratio.

2

Algorithm

General purpose CPUs are designed for a wide variety of applications and take limited advantage of the
inherent parallelism in many calculations. Improving
performance in the past has relied on increasing clock
speeds and the size of high speed cache memories. Programming a CPU for high performance scientific applications involves careful data layout to utilize the cache
optimally and careful scheduling of instructions.
In contrast, graphics processors are designed for intrinsically parallel operations, such as shading pixels,
where the computations on one pixel are completely independent of another. GPUs are an example of streaming processors, which use explicit data parallelism to
provide high compute performance and hide memory
latency. Data is expressed as streams and data parallel operations are expressed as kernels. Kernels can
be thought of as functions that transform each element
of an input stream into a corresponding element of an
output stream. When expressed this way, the kernel
function can be applied to multiple elements of the input stream in parallel. Instead of blocking data to fit
caches, the data is streamed into the compute units.
Since streaming fetches are predetermined, data can be
fetched in parallel with computation. We describe below how the N -body force calculation can be mapped

to streaming architectures.
In its simplest form the N -body force calculation can
be described by the following pseudo-code:
for i = 1 to N
force[i] = 0
ri = coordinates[i]
for j = 1 to N
rj = coordinates[j]
force[i] = force[i] + force_function( ri, rj )
end
end

Since all coordinates are fixed during the force calculation, the force computation can be parallelized for the
different values of i. In terms of streams and kernels,
this can be expressed as follows:
stream coordinates;
stream forces;
kernel kforce( ri )
force = 0
for j = 1 to N
rj = coordinates[j]
force = force + force_function( ri, rj )
end
return force
end kernel
forces = kforce( coordinates )

The kernel kforce is applied to each element of the
stream coordinates to produce an element of the
forces stream. Note that the kernel can perform an
indexed fetch from the coordinates stream inside the
j-loop. An out-of-order indexed fetch can be slow, since
in general, there is no way to prefetch the data. However
in this case the indexed accesses are sequential. Moreover, the j-loop is executed simultaneously for many ielements; even with minimal caching, rj can be reused
for many N i-elements without fetching from memory
thus the performance of this algorithm would be expected to be high. The implementation of this algorithm
on GPUs and GPU-specific performance optimizations
are described in the following section.
There is however one caveat in using a streaming
model. Newton’s Third law states that the force on
particle i due to particle j is the negative of the force
on particle j due to particle i. CPU implementations
use this fact to halve the number of force calculations.
However, in the streaming model, the kernel has no
ability to write an out-of-sequence element (scatter), so
forces[j] can not be updated while summing over the
j-loop to calculate forces[i]. This effectively doubles
the number of computations that must be done on the
GPU compared to a CPU.
Several commonly used force functions were implemented to measure and compare performance. For stellar dynamics, depending on the integration scheme being used, one may need to compute just the forces, or
the forces as well as the time derivative of the forces
(jerk). We have designated the corresponding kernels
GA (Gravitational Acceleration) and GAJ (Gravitational Acceleration and Jerk). In molecular dynamics,

it is not practical to use O(N 2 ) approaches when the
solvent is treated explicitly, so we restrict ourselves to
continuum solvent models. In such models, the quantum interaction of non-bonded atoms is given by a
Lennard-Jones function and the electrostatic interaction is given by Coulomb’s Law suitably modified to
account for the solvent. The LJC(constant) kernel calculates the Coulomb force with a constant dielectric,
while the LJC(linear) and LJC(sigmoidal) kernels use
distance dependent dielectrics. The equations used for
each kernel as well as the arithmetic complexity of the
calculation are shown in Table 1.

3
3.1

Implementation and Optimization on GPUs
Brook

BrookGPU [Buck et al. 2004b] is a C-like high-level language that can be used to program GPUs as streaming processors. Streams are stored as textures and
kernels are implemented as fragment programs. The
BrookGPU run-time library can utilize a number of
graphics interfaces; for this work we used the Microsoft
DirectX 9.0c API and the Pixel Shader 3.0 specification [Microsoft 2006]. DirectX [Microsoft 2003] provides a vendor-independent abstraction of hardware
features. In the Pixel Shader 3.0 specification, the
shader has access to 32 general purpose, 4-component,
single precision floating point (float4) registers, 16
float4 input textures, 4 float4 render targets (output
streams) and 32 float4 constant registers. A shader
consists of a number of assembly-like instructions. Current GPUs have a maximum static program length of
512 (ATI) or 1024 (NVIDIA) instructions. There is a
loop limit of 255 iterations of a loop body, but loops can
be nested to increase the total numbers of iterations.
NVIDIA is limited to 65,535 dynamic instructions and
ATI can support an unlimited number. The BrookGPU
compiler translates kernels into a high level shader language like CG or HLSL, which is then compiled into
pixel shader assembly by an appropriate shader compiler like Microsoft’s fxc or NVIDIA’s cgc. The graphics driver finally maps the Pixel Shader assembly code
into hardware instructions as appropriate to the architecture.

3.2

Precision

Recent graphics boards have 32-bit floating point arithmetic. Consequently we have done all the calculations
in single precision. Whether or not this is sufficiently accurate for the answers being sought from the simulation
is often a subject of much debate and the authors do not
intend to settle it here. We are of the opinion that in

many cases, though certainly not all, single precision is
enough to obtain useful results. Furthermore, if double
precision is necessary, it is usually not required throughout the calculation, but rather only in a select few instances. For reference, GRAPE-6 [Makino et al. 2003a]
performs the accumulation of accelerations, subtraction
of position vectors and update of positions in 64-bit fixed
point arithmetic with everything else in either 36, 32 or
29 bit floating point precision. It is quite common to do
the entire force calculation in single precision for molecular simulations while using double precision for some
operations in the update step. If and where necessary,
the appropriate precision could be emulated on graphics boards [Göddeke et al. 2005]. The impact on performance would depend on where and how often it would
be necessary to do calculations in double precision.

3.3

General Optimization

The algorithm was implemented for several force models. For simplicity, in the following discussion, we
only talk about the GA kernel, which corresponds to
the gravitational attraction between two mass particles,
given by
X
mj
(1)
ai = −G
2 + ǫ2 )3/2 rij
(r
ij
i6=j
where ai is the acceleration on particle i, G is a constant
(often normalized to one), mj is the mass of particle j,
and rij is the vector displacement between particles i
and j. The performance of the kernel for various input
sizes are shown in Figure 1.
The algorithm outlined in Section 2 was implemented in BrookGPU and targeted for the ATI
X1900XTX. Even this naive implementation performs very well, achieving over 40 GFlops, but its
performance can be improved.
This kernel executes 48 Giga-instructions/sec and has a memory bandwidth of 33 GB/sec.
Using information
from GPUBench [Buck et al. 2004a], we expect the
X1900XTX to be able to execute approximately 30-50
Giga-instruction/sec (it depends heavily on the pipelining of commands) and have a cache memory bandwidth
of 41GB/sec. The nature of the algorithm is such that
almost all the memory reads will be from the cache since
all the pixels being rendered at a given time will be accessing the same j-particle. Thus this kernel is limited
by the rate at which the GPU can issue instructions
(compute bound).
To achieve higher performance, we used the standard
technique of loop unrolling. This naive implementation
is designated as a 1×1 kernel because it is not unrolled
in either i or j. The convention followed hereafter when
designating the amount of unrolling will be that A×B
means i unrolled A times and j unrolled B times. The
second GA kernel (1×4) which was written unrolled the

Kernel

Formula

Gravity
(accel)

mj
(r2 +ǫ2 )3/2
ij

rij

mj

Gravity
(accel & jerk)

mj

vij
(r2 +ǫ2 )3/2
ij

−3

Unroll

Input
(bytes)

Inner
Loop
Insns.

BW
(GB/s)

Useful
GFLOPS

Giga
Intrxns
per sec.

System
Size

19

4×4

64

125

19.9

94.3

4.97

65,536

42

1×4

128

104

40.6

53.5

1.27

65,536

rij

(r2 +ǫ2 )3/2
ij

h

Flops
per
Intrxn.

(rij ·vij )rij
(r2 +ǫ2 )5/2
ij

i

LJC
(constant)

qi qj
ǫr3
ij

rij + ǫij

h

σij
rij

6

−



σij
rij

12 i

30

2×4

104

109

33.6

77.6

2.59

4096

LJC
(linear)

qi qj
r4
ij

rij + ǫij

h

σij
rij

6

−



σij
rij

12 i

30

2×4

104

107

34.5

79.5

2.65

4096

43

2×4

104

138

27.3

90.3

2.10

4096

LJC
(sigmoidal)

ǫij

h

qi qj
rij +
ζ(rij )r3
ij
6
12
σ
σij
− r ij
rij
ij
(αr3 +βr2 +γ+δ)





 i

ζ(r) = e

Giga-Interactions per second

Table 1: Values for the maximum performance of each kernel on the X1900XTX. The instructions are counted as the
number of pixel shader assembly arithmetic instructions in the inner loop. Intrxn = Interaction; Insns = Instructions;
BW = Bandwidth.
5
4
3

1x1
1x4
4x4

2

stream coordinates;
stream index = range( 1 to N skip 4 );
stream forces1, forces2, forces4, forces4;

1
0
100

requirements to one quarter of original bandwidth because each j-particle is only read by one-quarter as
many fragments. To make this more clear, we show
the pseudo-code for this kernel below. This kernel is
designated as a 4×4 kernel.

1000

10000

Output Stream Size

Figure 1: GA Kernel with varying amounts of unrolling
j-loop four times, enabling the use of the 4-way SIMD
instructions on the GPU. This reduces instructions that
must be issued by around a factor of 3. (We cannot
reduce instructions by a factor of 4 because some Pixel
Shader instructions are scalar). The performance for
this kernel is shown in Figure 1. It achieves a modest
speedup compared to the previous one, and we have
now switched from being compute bound to bandwidth
bound (35 Giga-Instructions/sec and ≈40GB/sec).
Further reducing bandwidth usage is somewhat more
difficult. It involves using the multiple render targets
(MRT) capability of recent GPUs which is abstracted as
multiple output streams by BrookGPU. By reading in
4 i-particles into each kernel invocation and outputting
the force on each into a separate output stream, we reduce by a factor of four the size of each output stream
compared with original. This reduces input bandwidth

kernel kforce4x4( i )
force1 = 0
force2 = 0
force3 = 0
force4 = 0
ri1 = coordinates[i]
ri2 = coordinates[i+1]
ri3 = coordinates[i+2]
ri4 = coordinates[i+3]
for j = 1 to N skip 4
rj1 = coordinates[j]
rj2 = coordinates[j+1]
rj3 = coordinates[j+2]
rj4 = coordinates[j+3]
force1
force2
force3
force4

+=
+=
+=
+=

force_function4(
force_function4(
force_function4(
force_function4(

ri1,
ri2,
ri3,
ri4,

rj1,
rj1,
rj1,
rj1,

rj2,
rj2,
rj2,
rj2,

rj3,
rj3,
rj3,
rj3,

rj4
rj4
rj4
rj4

)
)
)
)

end
return force1, force2, force3, force4
end kernel
forces1, forces2, forces3, forces4 = kforce4x4( indices )

In the above code, the input is the sequence of integers 1, 5, 9, ...N and the output is 4 force streams.
force function4 uses the 4-way SIMD math available
on the GPU to compute 4 forces at a time. The four
output streams can be trivially merged into a single one
if needed. Results for this kernel can be seen in Figure 1.

3.4

Giga-Interactions per second

Once more the kernel has become instruction-rate limited and its bandwidth is half that of the maximum
bandwidth of the ATI board, but the overall performance has increased significantly.

Optimization for small systems

In all cases, performance is severely limited when the
number of particles is less than about 4000. This is due
to a combination of fixed overhead in executing kernels and the lack of sufficiently many parallel threads
of execution. It is sometimes necessary to process small
systems or subsystems of particles (N ≈ 100 − 1000).
For example, in molecular dynamics where forces
tend to be short-range in nature, it is more common
to use O(N ) methods by neglecting or approximating
the interactions beyond a certain cutoff distance. However, when using continuum solvent models, the number of particles is small enough (N ≈ 1000) that the
O(N 2 ) method is comparable in complexity while giving greater accuracy than O(N ) methods.
It is common in stellar dynamics to parallelize the
individual time step scheme by using the block time
step method [McMillan 1986]. In this method forces
are calculated on only a subset of the particles at any
one time. In some simulations a small core can form
such that the smallest subset might have less than 1000
particles in it. To take maximal advantage of GPUs it is
therefore important to get good performance for small
output stream sizes.
To do this, we can increase the number of parallel
threads by decreasing the j-loop length. For example,
the input stream can be replicated twice, with the j-loop
looping over the first N/2 particles for the first half of
the replicated stream and looping over the second N/2
particles for the second half of the stream. Consider the
following pseudocode that replicates the stream size by
a factor of 2:
stream coordinates;
stream indices = range( 1 to 2N );
stream partial_forces;
kernel kforce( i )
force = 0
if i <= N:
ri = coordinates[i]
for j = 1 to N/2
rj = coordinates[j]
force = force + force_function( ri, rj )
end
else
ri = coordinates[i-N+1]
for j = N/2+1 to N
rj = coordinates[j]
force = force + force_function( ri, rj )
end
endif
return force
end kernel
partial_forces = kforce( indices )

2.5

2.0
4096
2048
1024
768

1.5

1.0

0.5
2

4

6

8

Replication of i-particles

Figure 2: Performance improvement for LJC(sigmoidal)
kernel with i-particle replication for several values of N
In this example, the stream indices is twice as long
as the coordinates stream and contains integers in sequence from 1 to 2N . After applying the kernel kforce
on indices to get partial forces, the force on particle
i can be obtained with by adding partial forces[i]
and partial forces[i+N], which can be expressed as
a trivial kernel. The performance of the LJC(sigmoidal)
kernel for different number of replications of the iparticles is shown in Figure 2 for several system sizes.

4

Results

All kernels were run on an ATI X1900XTX PCIe graphics card on Dell Dimension 8400 with ATI Catalyst version 7.2 drivers and the latest DirectX SDK (December
2006). A number of different force models were implemented with varying compute-to-bandwidth ratios
(see Table 1). A sample code listing is provided in
the appendix (A.1) to show the details of how flops are
counted.
To compare against the CPU, a specially optimized
version of the GA and GAJ kernels were written since
no software suitable for a direct comparison to the GPU
existed. The work of [Nitadori et al. 2005] uses SSE for
the GAJ kernel but does some parts of the calculation in
double precision which makes it unsuitable for a direct
comparison. The performance they achieved is comparable to the performance achieved here. Using SSE
intrinsics and Intel’s C++ Compiler v9.0, we were able
to obtain sustained performance of 3.8 GFlops on a 3.0
GHz Pentium 4.
GROMACS [Lindahl et al. 2001] is currently the
fastest performing molecular dynamics software with
hand-written SSE assembly loops. As mentioned in Section 2 the CPU can do out-of-order writes without a significant penalty. GROMACS uses this fact to halve the
number of calculations needed in each force calculation
step. In the comparison against the GPU in Table 2 the
interactions per second as reported by GROMACS have

been doubled to reflect this. In MD it is common to use
neighborlists to reduce the order of the force computation to O(N ). The performance of GROMACS doing
an O(N 2 ) calculation as well as an O(N ) calculation
for a 80 residue protein (lambda repressor, 1280 atoms)
is shown in Table2. Despite using a fairly modest cutoff
length of 1.2 nm for the O(N ) calculation, the O(N 2 )
GPU calculation represents an order-of-magnitude performance improvement over existing methods on CPUs.

Comparison to other Architectures

In Figure 3 is a comparison of interactions/sec between the ATI X1900XTX, GRAPE-6A and a Pentium
4 3.0GHz. The numbers for the GPU and CPU are
observed values, those for GRAPE-6A are for its theoretical peak. Compared to GRAPE-6A, the GPU can
calculate over twice as many interactions when only the
acceleration is computed, and a little over half as many
when both the acceleration and jerk are computed. The
GPU bests the CPU by 35x, 39x and 15x for the GA,
LJC(constant) and GAJ kernels respectively.
Another important metric is performance per unit
of power dissipated. These results can be seen in
Figure 5. Here the custom design and much smaller
on-board memory allows GRAPE-6A to better the
GPU by a factor of 4 for the GAJ kernel, although
they are still about equal for the GA kernel. The
power dissipation of the Intel Pentium 4 3.0 GHz
is 82W [Intel 2006], the X1900XTX is measured to
be 85W, we estimate GRAPE-6A’s dissipation to be
48W since each of the 4 processing chips on the board
dissipates approximately 12W [Makino et al. 2003b]
and MDGRAPE-3’s (MD3-PCIX) dissipation is
40W [Peta Computing Institute 2006].
The advantages of the GPU become readily apparent when the metric of performance per dollar is examined (Figure 4). The current price of
an Intel Pentium 4 630 3.0GHz is $100, an ATI
X1900XTX is $350, and an MDGRAPE-3 board costs
$16000 [Peta Computing Institute 2006]. The GPU
outperforms GRAPE-6A by a factor of 22 for the GA
kernel and 6 for the GAJ kernel.

5.2

Hardware Constraints

The 4×4 unrolling that is possible with the GA kernel
does not work for the other, more complicated kernels.
For example, the GAJ kernel requires two outputs per
particle (jerk in addition to acceleration). This reduces
the maximum unrolling possibility to 2×4 because the
GPU is limited to a maximum of 4 outputs per kernel.
However, even this amount of unrolling doesn’t work

Billion Interactions/s

CPU (observed)
GPU (observed)

4

GRAPE-6A (theoretical peak)
MDGRAPE-3 (observed)

3

2

1

0

GA

GAJ

LJC(constant)

Figure 3: Speed comparison of CPU, GPU, GRAPE-6A
and MDGRAPE-3

Million Interactions/sec/USD

5.1

Discussion

CPU (observed)
GPU (observed)
GRAPE-6A (theoretical peak)
MDGRAPE-3 (observed)

10

5

0

GA

GAJ

LJC(constant)

Figure 4: Useful MFlops per second per U.S. Dollar of
CPU, GPU, GRAPE-6A and MDGRAPE-3

Million Interactions/Watt

5

5

CPU (observed)
GPU (observed)
GRAPE-6A (theoretical peak)
MDGRAPE-3 (observed)

60

40

20

0

GA

GAJ

LJC(constant)

Figure 5: Millions of Interactions per Watt of CPU,
GPU, GRAPE-6A and MDGRAPE-3

Kernel
LJC(constant)
LJC(linear)*
LJC(sigmoidal)

GMX Million
Interactions/sec

GMX O(N 2 )
ns/day

GMX O(N )
ns/day

GPU Million
Interactions/sec

GPU
ns/day

66
33
40

5.6
2.06
2.5

18.2
9.08
11

1327
1327
1203

140
140
127

Table 2: Comparison of GROMACS(GMX) running on a 3.2 GHz Pentium 4 vs. the GPU showing the simulation time
per day for an 80 residue protein (lambda repressor) *GROMACS does not have an SSE inner loop for LJC(linear)

5.3

On-board Memory vs. Cache Usage

As mentioned in Section 3.3 we expect the kernels to
make very efficient use of the cache on the boards. There
are a maximum of 512 threads in flight on the ATI
X1900XTX at any one time [ATI 2006], and in the ideal
situation, each of these threads will try and access the
same j-particle at approximately the same time. The
first thread to request a j-particle will miss the cache and
cause the particle to be fetched from on-board memory,
however once it is in the cache, all the threads should
be able to read it without it having to be fetched from
on-board memory again.
For example, in the case of the GA kernel with 65,536
particles, there would be 16,384 fragments to be processed, and if fragments were processed in perfectly sep-

100
80

GFlops

because the compiler cannot fit the kernel within the 32
available registers. The number of registers is also what
prevents the LJC kernels from being unrolled by 4×4
instead of 2×4.
This apparent limitation due to the number of registers appears to result from compiler inefficiencies; the
authors are currently hand coding a 2×4 GAJ kernel directly in pixel shader assembly which should cause the
kernel to become compute bound and greatly increase
its performance. The performance gain of unrolling the
LJC kernels to 4×4 by rewriting them in assembly would
most likely be small since these kernels are already compute bound.
While the maximum texture size of 4096×4096 and
512 MB would make it possible to store up to 16 million
particles on the board at a time, this really isn’t necessary. In fact, GRAPE-6A only has storage for 131,000
particles on the board at any one time. This is small
enough to occasionally seem restrictive - a good balance
is around 1 million particles which could easily be accommodated by 64MB. If board manufacturers wanted
to produce cheaper boards specifically for use in these
kinds of computations they could significantly reduce
the cost without affecting the functionality by reducing
the amount of onboard RAM.
The current limits on the number of instructions also
impacts the efficiency of large GPGPU programs. On
ATI hardware, the maximum shader length of 512 instructions limits the amount of loop unrolling and the
complexity of the force functions we can handle.

GA
GAJ
LJC (sig)
LJC (linear)
LJC (const)

60
40
20
0
0

200

400

600

Memory Speed (MHz)

Figure 6: GFlops achieved as a function of memory
speed

arate groups of 512, then 32 groups would need to be
processed. Each group would need to bring in 65,536
particles from main memory to the cache resulting in an
extremely low memory bandwidth requirement of 38.2
MB/sec.
Of course, the reality is that particles are not processed in perfectly separate groups of 512 particles that
all request the same particle at the same time, but by using ATITool [ATITool 2006] to adjust the memory clock
of the board we can determine how much bandwidth
each kernel actually needs to main memory. The results of this testing can be seen in Figure 6.
The performance degradation occurs at approximately 11.3, 5.2, and 2.1 GB/sec for the LJC, GAJ
and GA kernels respectively. The LJC kernels must also
read in an exclusion list for each particle which does not
cache as well as the other reads, and is the reason why
their bandwidth to main memory is higher than that
of the gravity kernels. The number for the GA kernel
suggests that approximately 10 particles are accessing
the same j-particle at once.
At memory speeds above 500MHz all the kernels run
very near their peak speed, thus board manufacturers could not only use less RAM, they could also use
cheaper RAM if they were to produce a number of
boards that would only be used for these calculations.
This would reduce the cost and power requirements over
the standard high end versions used for gaming.

5.4

Distributed Computation

Most biological phenomena of interest occur on
timescales currently beyond the reach of MD simulations. For example, the simplest proteins fold on a
timescale of 5 to 20 microseconds, while more complex
proteins may take milliseconds to seconds. MD simulations on current generation CPUs are usually limited
to simulating about 10 nanoseconds per day - it would
take several years to obtain a 10 microsecond simulation. However, with the speed increases afforded by the
algorithms and hardware discussed here, we are now be
able to simulate protein dynamics with individual trajectories on the 10 microsecond timescale in under three
months. This will allow the direct simulation of the folding of fast-folding proteins. Moreover, by incorporating
this methodology into a distributed computing framework, we are now situated to build Markovian State
Models to simulate even longer timescales, likely approaching seconds [Jaychandran et al. 2006]. Thus with
the combined effort of GPUs and distributed computing, one would be able to reach timescales for folding of
essentially all single-domain, two-state folding proteins.
Compared to the donations of CPUs from over 150,000
Windows computers currently producing 145 TFlops,
we have 550 GPUs donated to the project producing
over 34 TFlops. Thus each GPU is providing roughly
60 times the performance of the average donated x86
CPU.

A

Appendix

A.1

Flops Accounting

To detail how we count flops we present a snippet of
the actual Brook code for the GA kernel. The calculation of the acceleration on the first i-particle has been
commented with our flop counts for each instruction.
In total, the calculation of the acceleration on the first
i-particle performs 76 flops. Since four interactions are
computed, this amounts to 19 flops per interaction.
float3 d1, d2, d3, d4, outaccel1;
float4 jmass, r, rinv, rinvcubed, scalar;
d1
d2
d3
d4

=
=
=
=

r.x
r.y
r.z
r.w

jpos1
jpos2
jpos3
jpos4
=
=
=
=

dot(
dot(
dot(
dot(

-

ipos1;
ipos1;
ipos1;
ipos1;

d1,
d2,
d3,
d4,

d1
d2
d3
d4

)
)
)
)

//3
//3
//3
//3
+
+
+
+

eps;
eps;
eps;
eps;

//6
//6
//6
//6

rinv = rsqrt( r );
rinvcubed = rinv*rinv*rinv;
scalar = jmass * rinvcubed;
outaccel1 += scalar.y * d2 +

//4
//8
//4
scalar.z * d3 + scalar.w * d4; //18

if ( Ilist.x != Jlist1.x ){
//don’t add force due to ourself
outaccel1 += scalar.x * d1; //6
}

References
ATI, 2006. Radeon X1900 product site.
http://www.ati.com/products/radeonx1900/index.html.

6

Conclusion

We have successfully taken advantage of the processing
power available on GPUs to accelerate pairwise force
calculations for several commonly used force models in
stellar and molecular dynamics simulations. In some
cases the GPU is more than 25 times as fast as a highly
optimized SSE-based CPU implementation and exceeds
the performance of custom processors specifically designed for these tasks such as GRAPE-6A. Furthermore, our performance is compute bound, so we are
well poised to take advantage of further increases in the
number of ALUs on GPUs, even if memory subsystem
speeds do not increase significantly. Because GPUs are
mass produced, they are relatively inexpensive and their
performance to cost ratio is an order of magnitude better than the alternatives. The wide availability of GPUs
will allow distributed computing initiatives like Folding@Home to utilize the combined processing power of
tens of thousands of GPUs to address problems in structural biology that were hitherto computationally infeasible. We believe that the future will see some truly
exciting applications of GPUs to scientific computing.

ATITool, 2006. techPowerUp.com.
http://www.techpowerup.com/atitool.
Buck, I., Fatahalian, K., and Hanrahan, P. 2004.
Gpubench: Evaluating gpu performance for numerical
and scientific applications. In Poster Session at GP2
Workshop on General Purpose Computing on Graphics
Processors. http://gpubench.sourceforge.net.
Buck, I., Foley, T., Horn, D., Sugerman, J., Fatahalian, K., Houston, M., and Hanrahan, P.
2004. Brook for GPUs: Stream Computing on Graphics Hardware. In Proceedings of the ACM SIGGRAPH
2004.
Elsen, E., Vishal, V., Darve, E., Hanrahan,
P.,
Pande,
V.,
and
Buck,
I.,
2005.
GROMACS on the
GPU.
http://bcats.stanford.edu/pdf/BCATS 2005 abstract book.pdf.
Fujitani, H., Tanida, Y., Ito, M., Jayachandran,
G., Snow, C. D., Shirts, M. R., Sorin, E. J.,
and Pande, V. S. 2005. Direct Calculation of the
Binding Free Energies of FKBP Ligands. J. Chem.
Phys. 123(8), 84108.
Fukushige, T., Makino, J., and Kawai, A. 2005.
GRAPE-6A: A Single-Card GRAPE-6 for Parallel PCGRAPE Cluster Systems. Publications of the Astronomical Society of Japan 57 (dec), 1009–1021.

Göddeke, D., Strzodka, R., and Turek, S. 2005.
Accelerating double precision FEM simulations with
GPUs. In Proceedings of ASIM 2005 - 18th Symposium
on Simulation Technique.
Gomez, D. R., Darve, E., and Pohorille, A. 2004.
Assessing the efficiency of free energy calculation methods. Journal of Chemical Physics 120, 8 (Feb), 3563–
78.

Taiji, M., Namuri, T., Ohno, Y., Futatsugi, N.,
Suenaga, A., Takada, N., and Konagaya, A.
2003. Protein explorer: A petaflops special-purpose
computer system for molecular dynamics simulations.
In SC ’03: Proceedings of the 2003 ACM/IEEE conference on Supercomputing.

Intel, 2006. Intel Pentium 4 Thermal Management.
http://www.intel.com/support/processors/pentium4/sb/CS007999.htm.
Jaychandran, G., Vishal, V., and Pande, V. S.
2006. Using Massively Parallel Simulations and Markovian Models to Study Protein Folding: Examining the
Villin Head-Piece. J. Chem. Phys. 124(6), 164903–
164914.
Juekuan Yang, Yujuan Wang, Y. C. 2006. Accelerated molecular dynamics simulation of thermal
conductivities. Journal of Computational Physics.
doi:10.1016/j.jcp.2006.06.039.
Lindahl, E., Hess, B., and van der Spoel, D.
2001. GROMACS 3.0: A package for molecular simulation and trajectory analysis. J. Mol. Mod. 7 , 306–317.
Makino, J., Fukushige, T., Koga, M., and Namura, K. 2003. GRAPE-6: Massively-Parallel
Special-Purpose Computer for Astrophysical Particle
Simulations. Publications of the Astronomical Society
of Japan 55 (dec), 1163–1187.
Makino, J., Kokubo, E., and Fukushige, T. 2003.
Performance evaulation and tuning of grape-6 - towards 40 ”real” tflops. In SC ’03: Proceedings of
the 2003 ACM/IEEE conference on Supercomputing,
IEEE Computer Society, Washington, DC, USA, 2.
McMillan, S. 1986. The Vectorization of Small-N
Integrators. In The Use of Supercomputers in Stellar
Dynamics, P. Hut and S. McMillan, Eds., 156–161.
Microsoft, 2003. DirectX home page.
http://www.microsoft.com/windows/directx/default.asp.
Microsoft, 2006. Pixel Shader 3.0 specification on
MSDN.
http://msdn.microsoft.com/library/default.asp?url=/library/enus/directx9 c/directx sdk.asp.
Nitadori, K., Makino, J., and Hut, P., 2005. Performance tuning of n-body codes on modern microprocessors: I. direct integration with a hermite scheme on
x86 64 architectures, Nov. http://arxiv.org/abs/astroph/0511062.
Peta Computing Institute, 2006. MDGRAPE-3
PCI-X.
Snow, C. D., Sorin, E. J., Rhee, Y. M., and
Pande, V. S. 2005. How Well Can Simulation Predict
Protein Folding Kinetics and Thermodynamics? Ann.
Rev. Biophys. Biomol. Struc. 34 , 43–69.

