Header guards

[lastchance]

Would you mind explaining why you think it matters, assuming the particle coordinates are indeed scaled correctly? I will say the mix of scaled coordinates but unscaled L & R is very weird, though.

It makes little sense to scale your coordinates by L (presumably the length of the box) and not scale every other length (including radius R) by the same. If you don't work in consistent units then you are asking for trouble ( https://en.wikipedia.org/wiki/Mars_Climate_Orbiter ). On another thread this was apparently because "L will change because the pressure changes", which is physically bonkers: L is here only to prescribe periodicity; it has no relation to global size. If you keep your coordinate fractions fixed whilst changing L then your absolute coordinates will move instantaneously in space ... otherwise known as teleportation.

That was a relatively minor point. The major one was also suggested on another thread - to do a quick check that only depended on simple comparisons first, since they are cheaper than multiplies. That check quickly excludes anything where the displacement vector doesn't fall within the bounding cube, which is overwhelmingly the majority case.

Sure, I stuck it in a loop, because that's what I thought @PhysicsIsFun wanted to do. If you want to compare one-by-one then simply remove the loop, pass the routine two particles rather than a vector and an index, and use return false; rather than continue; to signify that that particular comparison has concluded.

If you absolutely need to express radius R in units/scaling different from all your coordinates (which I strongly disagree with) then you could simply change the definition of diameter below to
double diameter = ( R + R ) / L;. Then at least all your overlap comparisons could be done in consistent units using diameter alone.

To be honest, if you intend to scale/non-dimensionalise then you should (a) do it consistently; (b) do it with respect to something fixed (for example, the radius).

#include <cmath>     
#include <algorithm> 
#include <vector>
using namespace std;

struct particle
{
   double x, y, z;
   // ... other properties of particle
};


// ...


// Somewhere in the calling routine ...
double diameter = R + R;
double diameterSquared = diameter * diameter;


//...


bool overlap( particle &p1, particle &p2, double diameter, double diameterSquared )
{
   double rx, ry, rz;

   // Calculate shortest distance between spheres (periodic boundaries!)
   rx = abs( p1.x - p2.x );
   ry = abs( p1.y - p2.y );
   rz = abs( p1.z - p2.z );

   // Periodic boundaries (assumes that one box, in NON-DIMENSIONAL COORDINATES, goes from 0 to 1)
   if ( rx > 0.5 ) rx = 1 - rx;
   if ( ry > 0.5 ) ry = 1 - ry;
   if ( rz > 0.5 ) rz = 1 - rz;

   // FAST CHECK for by far the most likely outcome (no overlap) - only uses comparisons
   if ( max( { rx, ry, rz } ) > diameter ) return false;     // Nowhere near, so doesn't need the next check

   // More expensive check (involves multiplies); will be true about pi/6 of the time if it actually gets here
   if ( rx * rx + ry * ry + rz * rz < diameterSquared ) return true;

   return false;
}

[PhysicsIsFun]

Ah well, @PhysicsIsFun, you aren't going to believe a word I say, anyway. I'll stay away from your posts in future.

I already explained the usage of coordinates (in the other thread) and even asked you where it's wrong according to you. You never explained, you never asked, you just kept assuming... don't know why you are so mad now.

Ok, so imagine a box of length L.
My particles of radius R are positioned in this box, thus their position vectors r lie in [-L/2, L/2]³, right?
The overlap condition for two spheres is given as |r₁-r₂|<2R.
Now, divide the positions by L.
Now we have coordinate vectors s in [-1/2, 1/2].
We must write the overlap condition in these new coordinates:
|r₁-r₂|= |Ls₁-Ls₂|<2R
This is equivalent to |s₁-s₂|< 2R / L.

So my pos-vector is the collection of s-vectors. My overlap condition is simply the squared version of that equation.

Why doing all this?
Because the box of my system, its volume, is fluctuating in time! It does not stay constant. L is changing (--> the volume changes in a certain way so that the pressure can stay constant). If L gets larger, i.e. if the volume gets larger, all the particles are drawn away from each other (their realtive position inside of the volume must stay constant). If L gets smaller, the particles get pushed towards each other.
The advantage of my scaled coordinates s in [-1/2, 1/2]³ is that they do not change, if L changes. Only the r vectors change. And L. But the fraction, namely the s vectors do not change. That's also the reason for the boundary conditions looking as they do.

So instead of calculating new r vectors any time my L changes, I can work with constant s vectors and only have to pass the new L to the overlap equation.

Does it make sense?

[helios]

Why not keep coordinates constant and normalize L to 1, and instead scale R? Then you wouldn't need to keep dividing R by L, and your collision check wouldn't need to know you're implicitly scaling the coordinates.

bool collide(const Circle &a, const Circle &b){
    return (a.center - b.center)^2 < (a.radius + b.radius)^2;
}

is a common collision check. No surprises.
This:

bool collide(const Circle &a, const Circle &b, double L){
    return (a.center - b.center)^2 < ((a.radius + b.radius) * L)^2;
}

is not.

[PhysicsIsFun]

Well, the L is basically the scaling factor of the radius, right?
And someone already hinted that I could bring the L to the other side so that I would not have the division, only a multiplication.
Other than that, I also need L at other parts of the simulation. In order to figure out what I would have to put in there, if I were to work with fluctuating R and L=1, I would need to calculate the fraction of current R and original R and multiply it with the "original L". So I would get divisions elsewhere.

The neighbor list (that I am going to implement) and @lastchance's tipp should do the trick.

[helios]

I say just don't pass L to the collision check function. It doesn't need to know about it.

[PhysicsIsFun]

Ah, you mean just scale R before calling overlap and only give R to the function?
This is a good idea, then I can easier reuse the overlap function in other programs. Thanks!