This webpage contains source codes of C++ programs to
generate hard-particle packings in a Euclidean space of arbitrary
dimension. If you use these codes, please reference this paper:
hard spheres in high dimensional Euclidean spaces", by M.
Skoge, A. Donev, F. H. Stillinger, and S. Torquato, Phys. Rev. E, Vol. 74: 041127 (2006)
[ibid 75: 029901 (2007)], [cond-mat/0608362]
The programs given here are general
molecular-dynamics codes. They can be used to prepape disordered jammed
packings of hard particles, if used with the appropriate options. The programs do not automatically generate
packings, yet alone jammed packings. The above paper should
provide background on the algorithm and appropriate parameters to use.
Additional more sophisticated executable codes can be found at http://cims.nyu.edu/~donev/Packing/PackLSD/Instructions.html,
including codes to analyze packings.
Please read the following instructions carefully before asking
questions, and contact the original author Monica Skoge (mskoge
at domain ucsd.edu)
and/or Aleksandar Donev (aleks.donev at
The simplest (but not
maintained!) version of code, suitable if you want to quickly modify it
or learn how the algorithm works, is found in the zip archive spheres.zip,
with periodic boundaries. The
usage is documented below. Reflecting (hard
boundaries are implemented in box_reflecting.C,
which replaces box.C. Packings inside cylindrical
containers in 3D can be generated using some modifications by Todd
Weisgraber (LLNL), which you can find in the archive
The most advanced and latest
code is one that can handle bidisperse packings, found in
the archive spheres_poly.zip
(the input file has the radius and mass ratios and the fraction of each
component). These codes in principle support any distribution of sphere
radii or masses, and support both periodic and reflecting boundaries,
but you will need to customize for your own needs. You can get
monodisperse packings by setting the bidispersity ratio to zero.
Instructions for Monodisperse Codes
The files needed to compile are neighbor.C, spheres.C, box.C,
sphere.C, event.C, heap.C, and read_ input.C
along with corresponding header (.h) files. The main program is
spheres.C and calls the core computational functions in box.C. Simply
g++ -O neighbor.C spheres.C box.C sphere.C event.C heap.C read_input.C
NOTE: With some compilers
you may need to rename our own vector.h into something like my_vector.h
and change the C++ sources accordingly because it may confuse it with
the STL header! Additionally, newer C++ compilers may not include the
standard C library headers by default, in which case you need to
"#include <cstdlib>" in spheres.C.
The codes use
The dimension of the packing is defined through the constant DIM in
vector.h and the input file containing the
parameters and input/output
file names is in input. In this file, the user can
either specify to read in an existing configuration of
hard spheres from a file by specifying the appropriate file name in
readfile or the user can
create a new configuration by setting readfile
to new. In the latter case, a
configuration of N spheres is
created at density initialpf
- A modified Lubachevsky-Stillinger algorithm, in which hard
grown in size and evolved according to Newtonian dynamics.
- Periodic boundary conditions applied to a hypercubic cell of
length one along each dimension.
- The linked-list cell method, in which the computational domain is
cubic cells and only neighboring cells are checked with predicting
collisions for a given sphere.
If the temperature temp is
set to zero, then all spheres are given zero
initial velocities; otherwise velocities are initially chosen from a
Maxwell-Boltzman distribution. The expansion rate growthrate governs
the rate at which the sphere diameter increases in units of the
simulation box size divided by the sphere velocity. Due to expansion,
the kinetic energy increases after collisions. Therefore velocities are
rescaled to 1 after every "cycle". The user specifies the number of
events per sphere in one cycle with the parameter eventspercyle.
Statistics, such as pressure, are calculated after every cycle and
outputted to datafile. The
simulation runs until reaching either the maximum allowed packing
fraction maxpf or the maximum
allowed pressure maxpressure.
terminating, the positions of all spheres are outputted to writefile.
The program produces two output files:
- writefile contains the
final positions of all the spheres. Its format is explained below.
contains the statistics calculated after each cycle. The statistics
include (as columns in this order) the current packing fraction,
pressure, total change
in kinetic energy before rescaling and the total number of MD events
The ASCII file formats used for storing sphere packing configurations
is as follows (in three dimensions):
n_dims ! Euclidian space dimensionality d
N n_species ! Total number of particles and number of species
N1 [N2...] ! Number of particles in each species
D1 [D2...] ! Sphere diameters for each species
Lambda ! Lattice vectors for unit cell, here a unit cube
T T T ! Periodic boundary conditions along all dimensions
x1 y1 z1 ! Coordinates of first sphere center
x2 y2 z2 ! second, etc...
The format may seem weird, but it is chosen to conform to the format
described at http://cims.nyu.edu/~donev/Packing/PackLSD/Instructions.html.