Distance and force calculation for a molecular dynamics simulation

I'm writing an MD-like simulation and I'm having some difficulties in making this code run faster. I profiled the code using callgrind and kcachegrind, and it appears I'm using about 30% of my time on calculating distances, 43% on setting up the linked-cell list and calculating my forces, 20% on malloc and free (I assume these are called when I create new vectors). In my simulation, I'm using two classes: 'particle', which holds two STL vectors (position and velocity); and 'system' which holds a vector of particles and the functions I use for my simulation.

The linked-cell method I'm using consists of dividing the system into blocks of equal size (in my case just '1'), and then only takes into account interactions between particles in neighbouring blocks. Periodic boundary conditions are applied.

The two functions which use pretty much all my CPU time are the following:

bool checkDistance(Particle& p1, Particle& p2){ //checks whether the distance between p1 and p2 is less than 1
double dist = 0, relDist; 
for ( int i = 0; i  halfSize){      //max separation in one dimension is half of system size with PBC
        if(relDist < 0){ //shift the relative distance to the shortest possible one (minimum image convention)
            relDist += systemSize;      
        }
        else {
            relDist -= systemSize;
        }
    }
    dist += relDist * relDist;
 }
 return (dist < 1); 
}

and

```
vector > System::calcForces(vector >& noise){ //calculate the force working on all particles in my system, adds a noise to this force (which is made in other function)

vector > forces (numberOfParticles, vector(dim,0.0)); //contains the force working on every particle

/ this looks up neighbors with linked cell method/

int numberOfCells = systemSize/interactionRange; //number of cells in one dimension. This way, all particles within interactionRange from a certain particle are in the same cell and in the neighboring cells.
vector > header(numbe

Given the additional information mentioned in the comments, I suggest the following changes:
Use Classes with Proper Members and Useful Methods

It appears that your Particle class has an array that represents a coordinate. This is an inefficient way of representing a coordinate, particularly if you know you will never need more than 3 dimensions. Iterating over an array instead of just directly handling values invokes a lot of overhead.

I would make 2 classes: Vector2D and Vector3D each with their coordinates explicitly defined. Something like:

class Vector2D {
    // ...methods, etc.
    
    float x;
    float y;
};

and

class Vector3D {
    // ...methods, etc.

    float x;
    float y;
    float z;
};

If you have proper methods on the class, such as a length() and distanceFrom(otherVector) method, then you don't need to iterate over an array. You can just call the method, which is easier to read, too.

It might make sense to either have an abstract VectorBase base class that has common methods from both Vector2D and Vector3D, or to have your Particle class be a template that takes a Vector2D or Vector3D as the template type.

As mentioned above, there will be common methods between the 2 vector types and this will simplify some of the code in calcForces(). For example, instead of normalizing all those vectors by having 2 for loops inside another for loop, you'd just have 1 loop:

for (int i = 0; i < numberOfParticles; i++) {
    unitVectors[i] = particleList[i].normalize();
}

You'll probably also want to implement things like operator+() and operator-() to simplify other parts of the code, such as at the end when you "calculate social force + friction + noise".
Make Loops More Efficient

You will inevitably need to have some complex loops. But you're repeatedly calculating various positions in those arrays multiple times in each loop. For example, you have this loop:

//calculate distance between particles in central cell + particles in neighboring cell. if smaller than interaction range, add to neighbors
        for (unsigned int k = 0; k < centralCell.size(); k++){
            neighborCount[centralCell[k]] +=1;
            for (int d = 0; d < dim; d++){
                averages[centralCell[k]][d] += unitVectors[centralCell[k]][d];
            }
            for(unsigned int l = 0; l < k; l++){ 
                if(checkDistance(particleList[centralCell[k]], particleList[centralCell[l]])){
                    neighborCount[centralCell[k]] +=1;
                    neighborCount[centralCell[l]] +=1;
                    for (int d = 0; d < dim; d++){
                        averages[centralCell[k]][d] += unitVectors[centralCell[l]][d];
                        averages[centralCell[l]][d] += unitVectors[centralCell[k]][d];
                    }
                }
            }
        }

In it, you calculate centralCell[k] 7 time and centralCell[l] 4 times. And then in the next loop you do something similar. It would be more efficient to do something like this:

for (unsigned int k = 0; k < centralCell.size(); k++)
{
    int index = centralCell[k];
    neighborCount[index] += 1;
    //... etc., replacing centralCell[k] with index everywhere
}

And of course, if you were using a proper Vector* class, you'd have fewer inner loops.
Don't Repeat Yourself

You have 5 nearly identical loops that make up a good portion of calcForces(). The first one is slightly different from the next 4 because you're putting them into the centralCell vector rather than the neighboringCell vector. The last 4 loops, at least, could be combined into 1 by putting all of the indexes into an array and iterating over that array. Something like this:

int neighborIndexes[4] = {
    header[i][j],
    header[(i+1) % numberOfCells][j],
    header[i][(j+1) % numberOfCells],
    header[(i+1) % numberOfCells][((j-1) % numberOfCells + numberOfCells) % numberOfCells]
};

for (int k = 0; k  -1) {
        neighboringCell.emplace_back(tempIndex);
        tempIndex = link[tempIndex];
    }
}

It's shorter and easier to read. There is slightly more overhead because we've added on outer loop, though. You'd need to profile to see whether it makes any difference.
Conclusions

Overall this code is very messy. You have one very long function that should be broken up into smaller, properly named pieces that are easier to understand at a glance.

You are allocating a very large number of vectors at the start of the calcForces function, and your profiling seems to indicate that's a problem. I think that having a Vector* class will reduce that somewhat. For example, instead of a 2D vector at the very beginning of calcForces, you'd only need a 1D vector, like this:

vector forces(numberOfParticles);

or if it was templatized:

vector forces(numberOfParticles);

StackExchange Code Review Q#158843, answer score: 3