#include "forces.cuh"

__global__ void CAC::calc_forces_and_energies(real *xs, real *forces,
                                              real *energies, int n_particles,
                                              real *box_len,
                                              PairPotential &potential) {
  int i = blockIdx.x * blockDim.x + threadIdx.x;
  if (i < n_particles) {
    real xi = xs[3 * i];
    real yi = xs[3 * i + 1];
    real zi = xs[3 * i + 2];

    for (int j = 0; j < n_particles; j++) {
      if (i != j) {
        real xj = xs[3 * j];
        real yj = xs[3 * j + 1];
        real zj = xs[3 * j + 2];

        real dx = xi - xj;
        real dy = yi - yj;
        real dz = zi - zj;

        // Apply periodic boundary conditions
        dx -= box_len[0] * round(dx / box_len[0]);
        dy -= box_len[1] * round(dy / box_len[1]);
        dz -= box_len[2] * round(dz / box_len[2]);

        ForceAndEnergy sol = potential.calc_force_and_energy({dx, dy, dz});
        forces[3 * i] += sol.force.x;
        forces[3 * i + 1] += sol.force.y;
        forces[3 * i + 2] += sol.force.z;
        energies[i] = sol.energy;
      }
    }
  }
}