Why race condition is hapenning in this MPI example?

When I execute and print the output of this MPI version of the NBody problem, this output differs from 1 thread execution and OpenMP execution. Can anyone help with the race condition in this example? I guess the problem is when I compute the forces or move the bodies, because the MPI_Scatter and MPI_Gather takes correctly the data from the bodies.

#include <stdlib.h>
#include <algorithm>
#include <iostream>
#include <fstream>
#include <sstream>
#include <cstring>
#include <vector>
#include <cstdlib>
#include <chrono>
#include <stdio.h>
#include <string.h>
#include <mpi.h>
#include "vector2.h"

/*
 * Compute forces of particles exerted on one another
 */
void ComputeForces(std::vector<Particle> &p_bodies, float p_gravitationalTerm, float p_deltaT)
{
    Vector2 direction,
        force, acceleration;

  Vector2 position1, position2,
    velocity1, velocity2;

    float distance;

    for (size_t j = 0; j < p_bodies.size(); ++j)
    {
    position1 = 0.f, position2 = 0.f;
        
    Particle &p1 = p_bodies[j];
      position1 = Vector2(p1.PositionX, p1.PositionY);
    
        force = 0.f, acceleration = 0.f;
        for (size_t k = 0; k < p_bodies.size(); ++k)
        {
            if (k == j) continue;
        
            Particle &p2 = p_bodies[k];
            position2 = Vector2(p2.PositionX, p2.PositionY);
      
            // Compute direction vector
            direction = position2 - position1;
            
            // Limit distance term to avoid singularities
            distance = std::max<float>( 0.5f * (p2.Mass + p1.Mass), fabs(direction.Length()) );
            
            // Accumulate force
            force += direction / (distance * distance * distance) * p2.Mass; 
        }
        // Compute acceleration for body 
        acceleration = force * p_gravitationalTerm;
        
        // Integrate velocity (m/s)
        p1.VelocityX += acceleration[0] * p_deltaT;
    p1.VelocityY += acceleration[1] * p_deltaT;
    }
}

/*
 * Update particle positions
 */
void MoveBodies(std::vector<Particle> &p_bodies, float p_deltaT)
{
    for (size_t j = 0; j < p_bodies.size(); ++j)
    {
        p_bodies[j].PositionX += p_bodies[j].VelocityX * p_deltaT;
     p_bodies[j].PositionY += p_bodies[j].VelocityY * p_deltaT;
    }
}


int main(int argc, char** argv)
{
    int maxIteration = 1000;
    float deltaT = 0.05f;
    float gTerm = 1.f;

    std::vector<Particle> bodies;
  
  //Variables for the struct.
  MPI_Datatype particletype, oldtypes[1]; 
  MPI_Aint offsets[1];
  int blockcounts[1];
  
  //Start MPI.
  int rc = MPI_Init(&argc, &argv); 
  
  //Variables for MPI.
  int rank, numtasks;
  
  MPI_Comm_rank(MPI_COMM_WORLD, &rank); 
  MPI_Comm_size(MPI_COMM_WORLD, &numtasks);
  
  //Initialize variables for Type_struct.
  offsets[0] = 0;
  oldtypes[0] = MPI_FLOAT;
  blockcounts[0] = 5;
  
  //Define the struct.
  MPI_Type_create_struct(1, blockcounts, offsets, oldtypes, &particletype);
  MPI_Type_commit(&particletype);
   
  std::vector<Particle> particles;
  std::vector<Particle> final;
  
  if (rank == 0)
    final.resize(bodies.size());
    
  particles.resize(bodies.size() / numtasks);
  
  //Start the timer.
    double start = MPI_Wtime();
    for (int iteration = 0; iteration < maxIteration; ++iteration)
    {
    MPI_Scatter(bodies.data(), particles.size(), particletype, particles.data(), particles.size(), particletype, 0, MPI_COMM_WORLD);

      ComputeForces(particles, gTerm, deltaT);
      MoveBodies(particles, deltaT);
      
    MPI_Gather(particles.data(), particles.size(), particletype, final.data(), particles.size(), particletype, 0, MPI_COMM_WORLD);

      if(enableOutput && rank == 0){
            fileOutput.str(std::string());
            fileOutput << "nbody_" << iteration << ".txt";
            PersistPositions(fileOutput.str(), final);
        }
    }
  //Finish the timer.
    double finish = MPI_Wtime();

  //Calculate and show the time.
  if(rank == 0){
      double elapsed = finish - start;
      std::cout << "Elapsed time: " << elapsed << " s\n";
  }
    
  //End MPI.
  MPI_Finalize();
 
    return 0;
}