#ifndef FORCES_CUH
#define FORCES_CUH
#include "potentials/pair_potentials.cuh"
#include "precision.hpp"
#include <cstdio>
#include <type_traits>
#include <variant>
#include <vector>

namespace CAC {

inline void reset_forces_and_energies(int n_particles, real *forces,
                                      real *energies) {
  cudaMemset(forces, 0, n_particles * sizeof(real) * 3);
  cudaMemset(energies, 0, n_particles * sizeof(real));
}

template <typename PotentialType>
__global__ void calc_forces_and_energies(real *xs, real *forces, real *energies,
                                         int n_particles, real *box_len,
                                         PotentialType 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;
      }
    }
  }
}

inline void launch_force_kernels(real *xs, real *forces, real *energies,
                                 int n_particles, real *box_len,
                                 std::vector<PairPotentials> potentials,
                                 int grid_size, int block_size) {

  reset_forces_and_energies(n_particles, forces, energies);

  for (const auto &potential : potentials) {
    std::visit(
        [&](const auto &potential) {
          using PotentialType = std::decay_t<decltype(potential)>;
          calc_forces_and_energies<PotentialType><<<grid_size, block_size>>>(
              xs, forces, energies, n_particles, box_len, potential);
        },
        potential);
    cudaDeviceSynchronize();
  }
}
} // namespace CAC

#endif