So I am trying a molecular dynamics simulation, and trying to populate a 3D rectangular volume with random particles which are uniformly distributed throughout the volume. Each molecule has a fixed radius r(sphere), where r can be different for each molecule. So I would want the sphere to be generated in the volume, and then no other sphere should exist within ( 2(r+tolerance) ) distance from that point(the tolerance will be really small, like 10^-6).

Also, the channel will have a specified length, breadth and width, so the randomization should be done within those limits in a single direction with a similar condition above i.e. the molecule should not be generated if the molecule is closer to the wall than its radius+tolerance.

I was initially trying a structured lattice, but that would mean I the number of molecules would be correlated to the dimensions, which is does not work for me. So I wrote the following algorithm(I wanted to paste my code, but it doesn't work and its a big, exaggerated mess right now, and has almost crashed my computer once.

So the logic is

1.) For each particle, make max_X=X-(radius+tol). Do it same for Y and Z.

2.) For each particle, generate random number between 0 and max_X.

3.) calculate distance between each particle and create a list of particles which violate any of the above conditions.

4.) Go over the list and re-generate those particles.

5.) Create another list of violating particle pairs.

6.) Go over. Rinse and repeat until the size of the list is zero.

So, it doesn't work. And I need this code to perform eventually for a large number of particles, so I need it as efficient and OpenMP paralleizable as possible. I am using the standard cliched #pragma omp parallel for method in front of the for loop for parallelizing the calculations, but it doesn't work with or without the pragma, at runtime, after like 5-10 minutes of a superpowered fan noise.

I am someone from a non-coding background who recently started learning writing complex code in C++, so I cannot do something fancy like structures or classes or pointers right now. I am working with std::vectors though. Would love it if you guys could show me a way out, and if you are in and around Raleigh, I will personally drive over and hand you an ice cold beer. Have been trying to do this for days, and is the only major non-functional thing in my code.


PS: Just to be clear, this isn't homework, else I would have asked the professor by now. Its part of an independent project which I will be freely distributing after completion.

  • $\begingroup$ How many particles are you expecting to have this run for? Have you considered trying to parallel the number crunching using your GPU with CUDA or OpenCL? Additionally, you could try to generate the positions of particles using a Latin Hypercube sampling that satisfies the various constraints but scales independently of the number of spatial dimensions of your problem. Also, are you able to push your cod to Github or something so we can view it there? $\endgroup$
    – spektr
    May 8 '17 at 20:40
  • $\begingroup$ 1.) Right now, about 10^4 would be fine, later I would like to scale it to about 10^5 or 10^6. 2.)Well at the moment I am pretty new to parallel programming, so I am relying on OpenMP while I learn. But I wanted to make a parallelizable logic first, and then start using CUDA or MPI. $\endgroup$ May 8 '17 at 21:46
  • $\begingroup$ 3.) LHS seems it might work, but I am really unsure of what LHS is and how I would achieve the minimum distance constraint with it. Is there a good resource that takes that constraint into account too? Or if you could help me out with a pseudocode which might work within my constraints...I would be really grateful. 4.) I have never used GitHub yet and don't know much about it. I don't want to be weird, but my code is too embarrassingly messed up and distributed among multiple functions. $\endgroup$ May 8 '17 at 21:46
  • $\begingroup$ Also, is there a similar algorithm which would do well for cylindrical and irregular geometries as well? $\endgroup$ May 8 '17 at 22:29

So one thing you can do is implicitly represent a structured region, like a 3D rectangular volume, as a set of chopped up sub-volumes where each sub-volume automatically enforces the spatial constraints you require. You will do this implicitly by just storing how many pieces to slice up each dimension into individually and referencing each sub-volume using a unique tuple of indices.

You can then compute how many slices you should have for each dimension, such that your spatial constraints are met, using your values for a particle's radius and the tolerance.

Then from here you can iteratively Monte Carlo the next sub-volume to stick a particle in. To efficiently check if you Monte Carlo a sub-volume you have used already, you can compute a unique hash, given a sub-volume's tuple of indices, and store it inside an efficient data structure, like std::map or std::unordered_map in C++, to be able to check for existence against quickly.

A sample code I wrote in C++ for this is below:


#ifndef _particle_cloud_
#define _particle_cloud_

#include <vector>
#include <cstdint>

namespace particle {

    class cloud {

        // useful typedefs
        typedef std::vector<double>     position;
        typedef std::vector<position>   positions;

        // ctors / initializers
        cloud(int num_particles, double xrng[2], double yrng[2], double zrng[2], double radius, double tolerance);
        void init(int num_particles, double xrng[2], double yrng[2], double zrng[2], double radius, double tolerance);

        // getter for particle position
        const position & getParticlePositionAt(int idx) const;

        // get number of particles in particle cloud
        int numParticles() const;

        typedef uint64_t h_int;
        enum Bounds { Low = 0, High = 1};
        positions ps;

        // unique hash function
        static h_int hash(h_int ix, h_int iy, h_int iz, h_int Nx, h_int Ny);




#include "particle_cloud.hpp"
#include <map>
#include <random>
#include <cmath>

namespace particle {


    cloud::cloud(int num_particles, double xrng[2], double yrng[2], double zrng[2], double radius, double tolerance)
        init(num_particles, xrng, yrng, zrng, radius, tolerance);

    void cloud::init(int num_particles, double xrng[2], double yrng[2], double zrng[2], double radius, double tolerance)
        // init constant quantities
        const double r = radius;
        const double e = tolerance;
        const double lb[3]      = { xrng[Low], yrng[Low], zrng[Low] };
        const double ub[3]      = { xrng[High], yrng[High], zrng[High] };
        const double span[3]    = { ub[0]- lb[0], ub[1]- lb[1], ub[2]- lb[2] };
        const h_int N[3]        = { std::floor(0.5*span[0]/(r+e)), std::floor(0.5*span[1] / (r + e)), std::floor(0.5*span[2] / (r + e)) };
        const double del[3]     = { span[0] / static_cast<double>(N[0]), span[1] / static_cast<double>(N[1]) , span[2] / static_cast<double>(N[2]) };

        // init array of temp map indices
        h_int idx[3]            = { 0 };

        // init temp position vector
        position tmp(3, 0.0);

        // size up the positions container correctly
        ps.resize(num_particles, tmp);

        // init the uniform distribution random number generator
        std::default_random_engine generator;
        std::uniform_real_distribution<double> U(0, 1);

        // create a map container to check if a given hash
        // has already been created
        std::map<h_int, bool> exist_map;

        // init variable that will say if the new position generated
        // is valid or not
        bool noGoodPos = false;

        // loop and generate the particle positions
        for (unsigned int i = 0; i < num_particles; ++i) {
            noGoodPos = true;

            // try to generate a valid position via Monte Carlo sampling until one has been made
            while (noGoodPos) {
                for (int d = 0; d < 3; ++d) { idx[d] = U(generator)*N[d]; }
                h_int h_ = hash(idx[0], idx[1], idx[2], N[0], N[1]);

                if (exist_map.count(h_) != 0) { }
                else { 
                    exist_map[h_] = true;
                    noGoodPos = false; 

            // compute each component of position
            // based on indices and dimensions of container
            for (int d = 0; d < 3; ++d) {
                tmp[d] = lb[d] + (del[d] * idx[d]) + (del[d] * 0.5);

            // copy temp position to ith position in positions container
            ps[i] = tmp;


    const typename cloud::position & cloud::getParticlePositionAt(int idx) const
        return ps[idx];

    int cloud::numParticles() const
        return (int)ps.size();

    typename cloud::h_int cloud::hash(h_int ix, h_int iy, h_int iz, h_int Nx, h_int Ny)
        return ix + Nx*(iy + Ny*iz); // unique hash based on spatial indices



#include <stdio.h>
#include "particle_cloud.hpp"

int main(int argc, char** argv) {

    particle::cloud pc;
    double v[2] = { 0, 1 };
    pc.init(100000, v, v, v, 1e-4, 1e-8);

    FILE * file_ = fopen("points.csv", "w");
    if (file_) {
        for (int i = 0; i < pc.numParticles(); ++i) {
            const particle::cloud::position & p = pc.getParticlePositionAt(i);
            fprintf(file_, "%0.8e, %0.8e, %0.8e\n", p[0], p[1], p[2]);

        fclose(file_); file_ = nullptr;

    return 0;

A sample graphic based on generating 1000 particle positions using the above code is shown below as well:

Sample Particle Distribution using Algorithm

  • $\begingroup$ @AyushAgrawal Let me know if you come across questions as you go along trying to understand the code! And nah, I'm not around Raleigh haha. Good luck with getting your code together! $\endgroup$
    – spektr
    May 9 '17 at 1:25
  • $\begingroup$ Not sure it would be uniformly per dV distributed $\endgroup$ May 11 '17 at 17:45
  • $\begingroup$ @SeverinPappadeux What makes you uncertain? $\endgroup$
    – spektr
    May 11 '17 at 18:05
  • $\begingroup$ I suspect the density in the subvolume would be different from density in the same volume but including the subdivision boundary. At least it requires proof/simulation that density is the same $\endgroup$ May 12 '17 at 17:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.