This webpage contains executable versions of Fortran programs to generate and analyze two- and three-dimensional hard-particle packings using event-driven molecular dynamics, as explained in this paper:

"Neighbor List Collision-Driven Molecular Dynamics Simulation for Nonspherical Particles. I. Algorithmic Details II. Applications to Ellipses and Ellipsoids", by A. Donev, F. H. Stillinger, and S. Torquato, J. Comp. Phys, 202(2): 737-764 (part I) and 765-793 (part II),  2005

We have also released C++ source codes (this is a locally-maintained copy of the Torquato group site) for a much stripped-down and less-efficient version that only handles packings of spheres, however, in an arbitrary dimensional Euclidian space.

Disclaimer: 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. I would be glad to assist in preparing the input configuration files to perform whatever specific task you have, just e-mail me at aleks.donev@gmail.com.

The tar file of all Linux executables and input files needed to play around with event-driven MD and analyze the resulting packings is PackLSD.tar.gz, but you can also look at the list of files directly since I sometimes add new exectuables there.

For now I have only documented the sphere packing options and hidden the rest, in order to avoid additional complexities for ellipsoids. I have also hidden options related to neighbour lists and other optimizations that require close familiarity with the algorithm and are hard to explain in detail. Contact me if you need more, or see the last section in this webpage (Advanced Usage). If you use this to produce published results, please reference the above paper.

Please read the following instructions carefully before asking questions, and also take a look at the examples I provided!

Producing Packings using MD

The binaries for generating jammed packing using hard-particle MD are PackLSD.static.x (or get the 64-bit version for x86_64 systems) for packings which consist of a small number of particle species (monodisperse, bidisperse, tridisperse, etc.), or PackLSD.poly.static.x for truly polydisperse packings where each particle has a different shape/size. The executables only link in standard (OpenGL, libtiff, libjpeg, etc.) Linux libraries, and link openglut statically, which may cause failures on some systems. The programs are compiled at medium optimization, for a generic IA32 processor, on an x86 Linux OS. They should run on most x86 Linux distros support just fine, and hopefully other similar distributions as well. To run these on x86_64 systems you will need to have some compatibility support for 32-bit executables; see this script for the packages you need to install for ubuntu 22.04 (jammy).

The runtime input file PackLSD.input or PackLSD.poly.input is to be piped into the executable. The main input file is a collection of Fortran 90 namelists as in the examples Input/PackLSD.nml and Input/PackLSD.poly.nml. The name of the packing can be changed inside PackLSD.input and is used as the base of all input/output filenames (as in the filename PackLSD.nml). The program usually prints when files are opened or closed. The executable needs to be run in a directory where there will be sub-directories called Input, Output, Stats, Data, Scenes, Plots, and VRML, and the files are placed or read from these subdirectories. The nml file must be in the Input directory, and the output packing files will be saved to Data. I have documented the options you will need in the nml file and removed most others (some need to be left in since the defaults are not appropriate for generic uses). Most options are the same for both programs PackLSD.static.x and PackLSD.poly.static.x, and the minor differences can be seen in the example Input/PackLSD.poly.nml. Please ask me if some option is not clear.
Important: The option lsdpg_options%lsdpg_n_events determines the number of events per stage. This must be larger than the number of particles, typically, you want it to be a multiple of the size of the system, for example, you want 10 events per particle.

The simulation processes a certain number of MD events (collisions), called a stage. At the end of every stage certain statistics are output, possibly a temporary file saved, the OpenGL view is refreshed, particle velocities rescaled, etc. Many stages can be run until jamming is reached, as defined by the options in the namelist file. The executable can use OpenGL to render to screen the packing process so you can see what is happening. You can use the mouse to zoom in, rotate, etc. There are also options to generate snapshot images that can be used to make animations (I highly recommend the MNG file format!)

Analyzing Packings

The produced packings can later be viewed/analyzed by using the programs ProcessPacking.static.x or ProcessPacking.poly.static.x. In particular, you can render using OpenGL, make plots, or compute the pair correlation function g2(r) or the structure factor S(k), as well as construct and render NNLs (near-neighbor lists) based on proximity or Delaunay/Voronoi neighbors. The input file is the namelist Input/ProcessPacking.nml in the Input subdirectory. It contains comments describing the various input fields. In 2D one can generate nice EPS figures to include in papers, and in 3D I provide some VRML models. These are not documented in detail (learn by looking at examples, although knowledge of VRML is required). The (anisotropic, vector k) structure factor in 3D is saved as a VTK file that can be looked at using the VisIt rendering program.

File Formats

The input/output packing files use a format which should be simple to understand, see the example dat file Data/LSD.HS.mono.3D.N=1000.1.dat. For packings that contain a small number of particle species (mono, bidisperse, etc.), this format is an ASCII text file as follows:

______________________
n_dims ! Euclidian space dimensionality (2 or 3)
N dispersity ! Total number of particles and number of species
N_1 N_2 ... ! Number of particles of each specie
D_1 D_2... ! Sphere diameters for each specie
Lambda ! Lattice vectors for unit cell
T T T ! Periodic or hard wall boundaries
x1 y1 z1 ! Coordinates of first sphere center
x2 y2 z2 ! second, etc...
...
______________________

The format is slightly different for truly polydisperse packings, the particle-specific diameter is after each particle's coordinates:
______________________
n_dims
N dispersity
N_1 N_2 ... ! One for each specie (each specie is a distribution for polydisperse packings, rather than a fixed size)
Lambda ! Unit cell
T T T ! Periodic or hard wall
x1 y1 z1 D1
x2 y2 z2 D2
...
______________________

Contact/Neighbour networks are written to files in ending with ".nnl.dat" and have the following format:
______________________
M2 N M
1 N_1
j1_1  k1_1  l1_1
j1_2  k1_2  l1_2
...
2 N_2
j2_1 k2_1 l2_1
....
______________________
Here M2 is the total number of contacts listed in the file. This includes both directions for contacts between particles, and M gives the total number of contacts without double-counting (M=M2/2 with periodic BCs, but it may change if hard-wall neighbors are included). N is the total number of particles. Then comes a list of all N_1 contacts of particle 1, listed one after the other, where j1_1 is the neighbour particle, k1_1 is an image unit cell identifier (explained shortly) and l1_1 is the strength of the interaction between the two particles (such as force---it can be ignored if not needed, but it has to be there as a floating-point number), etc. The image identifiers  is used for periodic boundary conditions to indicate contacts crossing the torus and is explained in the paper cited at the top of this page. It is zero if the contact is between particles in the base unit cell. After all the contacts for particle 1 have been listed, the same is done for particle 2, etc. Note that if a particle neighbour is negative that indicates a hard wall, with
neighbour=-(2*N+k)
where k indicates which of the 2d walls it is.

Advanced Usage: Full list of namelist options

Once familiar with the basic operation of the above programs, and one has read the necessary references, one can try to use also more advanced features such as near-neighbor lists (NNLs). For the packing-generation codes, the input file PackLSD.full.input contains all the input options. The namelist input file Input/PackLSD.2D.nml contains all of the input options (it is set up for 2D, with obvious generalization to 3D). For analyzing packings, the input file ProcessLSD.full.input contains the full list of input options, and the namelist Input/ProcessPacking.2D.nml contains a listing of all the available options. I recommend adding options gradually, becoming familiar with them first, and also making sure that the version of the executable program you are using is up to date and accepts the given namelist option (since I add new options all the time---if an option is not recognized, just comment it out and ask for an updated executable). It is not trivial to decide which option is compatible with other options or which option is required when some flag is set. Play around and ask me questions if needed.

There is also a version of the packing codes that generates packing of superellipses and superellipsoids. This is implemented by hacking the ellipsoid code to incorporate superellipsoid exponents. Note that this code is numerically sensitive for high exponents, and also that it is much slower than the ellipsoid codes. Note also that near neighbor lists and thus anything that requires them does not work for superellipsoids because the so-called inner overlap potential has not yet been implemented (or theoretically generalized) to superellipsoids. Executable programs that work with superellipsoids are PackLSD.super.static.x and ProcessPacking.super.static.x. As an example usage, try:
PackLSD.super.static.x < PackLSD.super.input
ProcessPacking.super.static.x < ProcessPacking.full.input