Generate a set of random x,y,z numbers, with a minimum difference between them, between defined limits

Question

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.

Help?!!!

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.

Answer 1

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:

particle_cloud.hpp

#ifndef _particle_cloud_
#define _particle_cloud_

#include <vector>
#include <cstdint>

namespace particle {

    class cloud {
    public:

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

        // ctors / initializers
        cloud();
        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;

    private:
        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);
    };

}


#endif

particle_cloud.cpp

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

namespace particle {


    cloud::cloud()
    {
    }

    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
    }


}

main.cpp

#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: