Hard spheres - check for overlap

[PhysicsIsFun]

Greetings,

in a Monte Carlo simulation of N hard spheres, radius R>0, I randomly choose a sphere and propose a new position for this particle in the direct vicinity of its original position.
Then I need to check wether the particle, at its proposed new position, overlaps with any other sphere.

Usually, one would write a simple loop that iterates through all the other spheres and checks wether the 'chosen' particle overlaps with any of them.

However, since I know that the new position of my particle is in the direct vicinity of its original position, I thought I could somehow use this to my advantage by only considering other spheres that are close to my particle's original position.
But then I would have to iterate through all spheres as well in order to figure out what spheres fit in that category, so the problem was just shifted.

Here is my code:

#include <cmath>
#include <vector>

using namespace std;

struct particle{
    double x;
    double y;
    double z;
};

//this function checks wether a chosen sphere u overlaps with any other sphere
bool overlap(vector<particle> pos, size_t u, const double R, double L){

    double rx, ry, rz;

    for(size_t j=0; j<pos.size(); j++){
        if(u==j) continue;
        //calculate shortest distance between spheres (periodic boundaries!)
        rx=abs(pos.at(u).x-pos.at(j).x);
        ry=abs(pos.at(u).y-pos.at(j).y);
        rz=abs(pos.at(u).z-pos.at(j).z);
        if(rx>0.5) rx=1-rx; //periodic boundaries
        if(ry>0.5) ry=1-ry;
        if(rz>0.5) rz=1-rz;

        if(rx*rx+ry*ry+rz*rz < 4*(R/L)*(R/L)){
            return true;
        }
    }
    return false;
}

pos-vector stores the positions of all the particles (the chosen particle is already on its proposed new position), size_t u is the vector index of the chosen particle. I am in relative coordinates here, which is why the overlap condition if(rx*rx+ry*ry+rz*rz < 4*(R/L)*(R/L))
looks like this (L=system length).

Under the condition that the new position of particle u, that I check for overlaps, is chosen randomly from within a small box of length a around particle u's original position:
Are there any ways to reduce computational effort here?

Regards!

[dhayden]

First, you should pass pos by reference. Otherwise the program makes a copy of the whole vector.

I think the way this is normally handled is to divide the space into cubic regions and keep track of which objects are in which regions. Note that an object can be in multiple regions at once. Then after you move the object and update the regions it's in, you check to see if it intersects the other objects in that (those) regions.

[keskiverto]

Websearch "collision tree algorithm".

[PhysicsIsFun]

dhayden wrote:

First, you should pass pos by reference. Otherwise the program makes a copy of the whole vector.

Good point, thanks! Up until now, I always passed by value, if the function was not supposed to change the passed argument, just to make sure that I don't accidentally change the passed argument in the function. But from a computational point of view, it seems to be better to pass by reference, at least when it comes to large objects.

How is it with single variables, like double or int? Is copying more expensive than passing by reference?

dhayden wrote:

I think the way this is normally handled is to divide the space into cubic regions and keep track of which objects are in which regions. Note that an object can be in multiple regions at once. Then after you move the object and update the regions it's in, you check to see if it intersects the other objects in that (those) regions.

keskiverto wrote:

Websearch "collision tree algorithm".

Thanks, I'll look into it.

[jonnin]

you don't have to iterate all the spheres. you can store a list of spheres eg central point and radius sorted by distance from center of your scene (0,0,0 or something). Then figure out where the new sphere would go and you only have to check spheres that could potentially overlap (some computation of new sphere's dist from center and its radius, existing sphere(s) values, ... if your new one is at 5,5,5 and its radius is 2, then it can only overlap spheres that are within those boundaries... the others are known to be too far away to bother checking against).

Whether such a home-cooked solution is better than the above is questionable, I am not sure. Spheres are a special, simple case of a harder problem and you may be able to get away with a simple and efficient algorithm.

If you enhanced that to know the quadrant each was in, it would reduce the number of candidate intersections farther, at the cost of a little more work/complexity -- a crude version of the cube tracking approach.

[tpb]

Up until now, I always passed by value, if the function was not supposed to change the passed argument, just to make sure that I don't accidentally change the passed argument in the function. But from a computational point of view, it seems to be better to pass by reference, at least when it comes to large objects.

If you don't want to change it then you make it a const reference:

bool overlap(const vector<particle>& pos, etc...)

[Jaybob66]

I always collision detect spheres by checking whether the distance between their centre's is more or less than their combined radiuses.

i tend to go for an octree when hit testing.

if there are more than a handful of objects in the volume
divide the volume into 8 octants (3d quadrants) and recurse into each octant that has more than a handful of items.
else
check if any 2 items in this octant are closer than their combined radiuses


i found it quite efficient.


EDIT

sorted by distance from center of your scene (0,0,0 or something). Then figure out where the new sphere would go and you only have to check spheres that could potentially overlap

i'm stealing that thought Jonnin :)

[PhysicsIsFun]

Jonnin wrote:

you can store a list of spheres eg central point and radius sorted by distance from center of your scene (0,0,0 or something). Then figure out where the new sphere would go and you only have to check spheres that could potentially overlap

But wouldn't I need to iterate over all spheres again to see which ones could "potentially overlap"?

tpb wrote:

If you don't want to change it then you make it a const reference:

I do want to change the vector, but not within the function. So it musn't be a const.

Jaybob66 wrote:

if there are more than a handful of objects in the volume
divide the volume into 8 octants (3d quadrants) and recurse into each octant that has more than a handful of items.
else
check if any 2 items in this octant are closer than their combined radiuses

So suppose my particle was in octant i. Then I would still need to "search" all particles that also are in octant i, so I still would need to go through all of them, or am I wrong.

And how do you take into account that a sphere can reach over into another box due to its R>0? You would have to do lots of extra checks, so is this really more efficient?

[jonnin]

If you do it right, you only need to check a subset of spheres, not all of them. That is really all I was saying... there are many, many ways to do that... but you can find a way to organize your sphere collection such that spheres that cannot intersect at all are not checked, and only potential collisions are checked, reducing the work to near zero assuming its not a degenerate problem of shoving billions of spheres in a small space. If you are tying to do something very crowded, you need to organize the spheres more cleanly (my approach was rough, and I forgot to say you need to check / manage a combination of radius and dist to center, not just the 1 item). The more crowded and higher the number of objects, the more precise you need your algorithm to be to do the least work possible. A crude method that discards 75% of a big open space with few spheres may work fine for some problems but when you have 10 billion spheres you may need to discard 99% of them with better boundary conditions, see?

[PhysicsIsFun]

Well, I am simulating a fixed number of spheres in a fluctuating volume. The volume can potentially grow so small that you reach around 80% of occupied space.

Now how about this?
-distribute space into cubes, each denoted by an index
-save not only coordinates of particles, but also cube index of its position in the vector
-propose a new position for particle u
-calculate current box index of u
-go through the particle vector and for each index, test wether a particle is in the same or in a neighboring box
-only for these particles perform the overlap check

But if the total volume becomes small, the spheres might become bigger than the (also fluctuating) size of a cube ~~

[jonnin]

if the radius of the spheres is bigger than the "radius" of the cubes, you end up with complexity. Ideally you should only have to check the cube where the new one would go, and all the neighboring cubes. If you have to check neighbors of neighbors, the cubes are too small. If a sphere is larger than a cube, you have to check them.

what you say will work, apart from that.

[PhysicsIsFun]

The problem is that my total V fluctuates. If I wanted to prevent the cubes becoming too small, I would have to implement them at a fixed size (so that they don't fluctuate with the total volume). Then this method would work well at large total V, when the spheres are distributed over several cubes, and it would work less well when the total V gets small, because then most of the spheres will be in the same cube and the method would reduce to the "brute force" method I used before.

It might still be a significant improvement. I am going to test it.

edit: Just realize that the total V can also become greater than the system of cubes of fixed size.

If there are any more ideas, feel free to let me now.

[jonnin]

the only other idea I have is to collapse it to a pair of 2-d spaces and check those in parallel via 2 threads. Then you can split to squares instead of cubes, or use quadrants, or try various other approaches to weed out things you don't need to check OR to do the checks quickly, one or the other... totally different approach though, you would need a fair bit of work.

[PhysicsIsFun]

I think that goes beyond my programming skill right now :D

But back to the original plan. I realized that my total V can not become smaller than the total volume of the spheres. If I do implement the cubes in relative size, so they fluctuate with the volume, and if I take a number of cubes that will have a length a>2R at the smallest possible total V (given by N times the number of sphere volume), a single sphere can never be bigger than a cube.
Then the approach in checking the spheres inside of the cube my target particle is in, as well as the spheres of the nearest neighbor cubes, would work, right?

[jonnin]

Its easier, you just would need to work to set it up because you are not looking at it this way.

Basically, solve the same type of problem in 2d, then reuse that solution a second time for a different look (eg if you solved xy plane, do it again for xz or zy) at your space. Its actually a simpler problem. The only complexity was threading it, which we can help you with if you want to take a crack at this approach. Then you just check both solutions and intersections in both are 'real' and its a collision.

I don't know that this is any better than anything else suggested.

[PhysicsIsFun]

The thing is, I am in a bit of a hurry and I can't afford to learn completely new stuff now.
I checked the simulation with my good computer at home (not on laptop) and it takes 'only' 6 minutes. I think with only 400 spheres, the speed gained through implementing a better overlap check isn't worth the time needed for me to set it up right now.

But thank you for your input!

[lastchance]

Have a look at the lines

        if(rx*rx+ry*ry+rz*rz < 4*(R/L)*(R/L)){
            return true;

How many multiplies (and, worse, divides) does that do?


It may do to have a quick and INEXPENSIVE check first to eliminate what ... in the majority of cases ... is going to return false.

Something like (in different notation)

if ( x < xmid - R || x > xmid + R || y < ymid - R || ... ) return false;
// Followed by more detailed comparison only for spheres that are definitely close. 

Shortcutting will allow this to return false (by far the most the likely result) as soon as any conditions are not satisfied, and you have done nothing more than add and subtract.


If you want your code to run fast you should be using normal [] indexing, not the .at() operation.

[PhysicsIsFun]

Thank you for the tipps @lastchance!

I had to drastically increase my simulation steps, so I definitely have to save time in this overlap stuff.


My Prof. told me about the neighbor-list concept:

I imagine another sphere with radius r_nl around the interaction sphere of my particle. The interaction sphere is the sphere of radius 2R around my particle (with radius R), because other particles only interact with my particle, if they come as close as 2R where they will touch each other.
For each particle, I have to calculate the spheres that are within r_nl and store these "neighbors" in a list.

Now, if I want to check wether particles overlap, I only need to check each particle and its corresponding list.
This will save effort, if I don't have to update the lists every simulation step.
I keep reusing the lists until the accumulated movement of one particle since the last list update is greater then 0,5*(r_nl - 2R). This is the first point where a particle might interact with another particle without that pair being considered in the list.

In C one would do this with pointers and linked lists.
Or using a 2D array with (i,j) being either 0 or 1, depending on wether the particles are close to each other.
How would I do this in C++ most efficiently?
Using such a classical array? Or storing a pointer in the structure of my particle, that leads to its list?

[lastchance]

The problem is that particles move. Unless you are in a crystalline solid the neighbour-list will be continually changing. Moreover, it will use more memory than the position, size and velocity you already store for each particle.

Since the overwhelming majority of particle pairs are NOT touching you need to be able to return false very early from an inexpensive preliminary check (only involving comparisons < and >, possibly simple additions/subtractions) rather than anything which will require multiplies to find the squared distance.

I have a similar problem when I am doing interpolation in a finite-volume mesh: the first thing you need to know is which cell a point is in. I rule out most cells by quick checks first before doing expensive precise checks.

The only aspect of neighbour lists that might come in handy is the fact that what was close on one timestep is likely still to be close on the next. It can also help to break domains down into smaller boxes, especially if you are using parallel processing. One processor then only needs to check particles in one box (and possibly its immediate neighbours).

[PhysicsIsFun]

Yeah, but it will not change with every step. And I am in a MC simulation where I always move a single particle at a time, there are no velocities.

I have no clue of parallel processing and the neighbor list seems to be more straightforward than the box method. My total V is fluctuating, that makes it harder.

Shoudn't the list in combination with your approach be good?

edit:
I tried to insert your tipps into my routine, but the crude test you suggested leads to slightly more computing time:

//this function checks wether a chosen sphere u overlaps with any other sphere
bool overlap(vector<particle> &pos, size_t u, const double R, double L){

    double rx, ry, rz;

    for(size_t j=0; j<pos.size(); j++){
        if(u==j) continue;
        //calculate shortest distance between spheres (periodic boundaries!)
        rx=abs(pos[u].x-pos[j].x);
        ry=abs(pos[u].y-pos[j].y);
        rz=abs(pos[u].z-pos[j].z);
        if(rx>0.5) rx=1-rx; //periodic boundaries
        if(ry>0.5) ry=1-ry;
        if(rz>0.5) rz=1-rz;
        if(rx*L>2*R || ry*L>2*R || rz*L>2*R) continue;
        if(rx*rx+ry*ry+rz*rz < 4*(R/L)*(R/L)){
            return true;
        }
    }
    return false;
}

I don't know what you meant with "return false". I only want to return false if there is not a single pair of interacting particles. I can only skip the more precise if-question with "continue".