Putting N hard spheres randomly in given volume

Question

I need to put $N$ spheres with given radius $R$ randomly in a Volume $[-0.5,0.5]^3$, without any overlap of spheres.

If I choose values so that all the spheres will occupy ~57% of the total volume, I find it difficult to get results. The basic scheme to place the $i$-th particle after $i-1$ particles have been placed:

• choose random position
• check for overlap with any other sphere
• if(overlap): check the new randomly chosen position
• repeat until $i$ has been placed
• place next particle

However, the more particles have been placed, the less likely my random generator will find a position with enough space AND the more time the overlap check needs. For $N=500$, $R=0.5/7.7$, the algorithm is barely able to place 340 particles (overlap is called billions of times at this point).

C++ code:

using namespace std;

//this function places N hard spheres of radius R on random coordinates in a [-0.5,0.5]^3 box

vector <particle> random_positions_hard_spheres (int N, double R){
vector <particle> positions;

static seed_seq seed_sequence { 100, 78639, 193, 55555, 3, 2348089, 6474367, 264441578 };
static mt19937 gen (seed_sequence);
uniform_real_distribution<double> random_double(-0.5, 0.5);
double x,y,z;
double d=2*R;

for(size_t i=0; i<N; ++i){

    x=random_double(gen);
    y=random_double(gen);
    z=random_double(gen);

    positions.push_back({x,y,z});

    //find position with no overlap
    int j=0;
    while(overlap_index(positions, i, d)){

        positions[i]={random_double(gen), random_double(gen), random_double(gen)};
        ++j;
    }
    cout <<"overlap was called "+to_string(j)+" times"<<endl;
    cout <<"particle "+to_string(i)+" was placed"<<endl;

}

return positions;
}

the overlap function checks whether particle $i$ overlaps with any previously placed particles.

Answer 1

Welcome to Computational Science StackExchange!

What you are trying to do is fairly difficult, and will need a more sophisticated method. That volume fraction of 0.57 lies just within the crystalline branch of the hard sphere phase diagram. The generation of dense random packings is an active area of current research.

In principle, a Monte Carlo simulation, or a collision-by-collision molecular dynamics program, can "equilibrate" an initial configuration prepared by randomly placing spheres which are significantly smaller than your desired size. Then, the idea is to steadily increase the sphere diameter, during the course of the simulation, until the desired volume fraction is reached.

A Monte Carlo code would not be too difficult to write, based on what you have already:

1. Insert the desired number of smaller spheres into the volume, at positions chosen at random, checking for overlap as you do at present.
2. Conduct Monte Carlo: loop over all $N$ spheres sequentially. Attempt to move each sphere by a small amount in a random direction. Test for overlap with the other $N-1$ spheres. If there is any overlap, reject the move and leave the sphere in its original place.
3. At intervals, attempt to increase the diameter of all the spheres by a very small amount, checking all $\frac{1}{2}N(N-1)$ overlaps. If there are any overlaps, keep the diameter unchanged and continue.

The alternative, of steadily increasing the diameter in the course of a molecular dynamics program, is one of the standard computational approaches to "random packing". Due to Lubachevsky and Stillinger, an example of its use can be found on arXiv (also published in Phys Rev Lett, 84 2064 (2000)).

The problem, of course, is that in principle, at that volume fraction, the system should be stable in the crystalline state, not in a randomly packed state. In practice, this may not be a problem: the system is likely to get "jammed" in a randomly packed state, at a volume fraction that depends on the compression rate, as discussed in the Phys Rev Lett.

MD programs for hard spheres can be found online, for example DynamO, or you can find guidance on how to write one in various books and papers. There are some standard tricks to improve the program speed, basically by making lists of near neighbours of the spheres: these would be worthwhile, but would not give huge speedups, for a system size of order 500.