Implement replicas and dipole field tracking to DipolarDirectSumGpu

src/core/energy.cpp

@@ -158,10 +158,4 @@ std::optional<double> System::particle_bond_energy(int pid, int bond_id,
   }
 }
 
-#ifdef DIPOLE_FIELD_TRACKING
-void System::calculate_long_range_fields() {
-  dipoles.calc_long_range_field(cell_structure->local_particles());
-}
-#endif
-
 } // namespace System

src/core/forces.cpp

@@ -128,6 +128,16 @@ void init_forces_and_thermostat(const CellStructure &cell_structure,
   // Initialize ghost forces (unchanged)
   init_forces_ghosts(cell_structure);
 }
+#ifdef DIPOLE_FIELD_TRACKING
+void invalidate_dip_fld(const CellStructure &cell_structure) {
+#ifdef CALIPER
+  CALI_CXX_MARK_FUNCTION;
+#endif
+
+  cell_structure.for_each_local_particle(
+      [](Particle &p) { p.dip_fld() = {0., 0., 0.}; });
+}
+#endif
 
 void init_forces_ghosts(const CellStructure &cell_structure) {
   cell_structure.for_each_ghost_particle(
@@ -181,6 +191,11 @@ void System::System::calculate_forces() {
     npt_inst_pressure->p_vir = Utils::Vector3d{};
   }
 #endif
+#ifdef DIPOLE_FIELD_TRACKING
+  // reset dipole field
+  invalidate_dip_fld(*cell_structure);
+#endif
+
   // Use combined function instead of two separate calls
   init_forces_and_thermostat(*cell_structure, *this);
 
@@ -258,6 +273,9 @@ void System::System::calculate_forces() {
   CALI_MARK_BEGIN("copy_forces_from_GPU");
 #endif
   gpu.copy_forces_to_host(particles, this_node);
+#ifdef DIPOLE_FIELD_TRACKING
+  gpu.copy_dip_fld_to_host(particles, this_node);
+#endif
 #ifdef CALIPER
   CALI_MARK_END("copy_forces_from_GPU");
 #endif

src/core/magnetostatics/dipolar_direct_sum.cpp

@@ -368,6 +368,9 @@ void DipolarDirectSum::add_long_range_forces(
     (*p)->force() += prefactor * fi.f;
     (*p)->torque() += prefactor * fi.torque;
   }
+#ifdef DIPOLE_FIELD_TRACKING
+  DipolarDirectSum::dipole_field_at_part(particles);
+#endif
 }
 
 /**

src/core/magnetostatics/dipolar_direct_sum_gpu.cpp

@@ -37,8 +37,12 @@ static void get_simulation_box(BoxGeometry const &box_geo, float *box,
   }
 }
 
-DipolarDirectSumGpu::DipolarDirectSumGpu(double prefactor) {
+DipolarDirectSumGpu::DipolarDirectSumGpu(double prefactor, int n_replicas) {
   set_prefactor(prefactor);
+  this->n_replicas = n_replicas;
+  if (n_replicas < 0) {
+    throw std::domain_error("Parameter 'n_replicas' must be >= 0");
+  }
 }
 
 void DipolarDirectSumGpu::on_activation() const {
@@ -47,6 +51,9 @@ void DipolarDirectSumGpu::on_activation() const {
   gpu_particle_data.enable_property(GpuParticleData::prop::torque);
   gpu_particle_data.enable_property(GpuParticleData::prop::pos);
   gpu_particle_data.enable_property(GpuParticleData::prop::dip);
+#ifdef DIPOLE_FIELD_TRACKING
+  gpu_particle_data.enable_property(GpuParticleData::prop::dip_fld);
+#endif
 }
 
 void DipolarDirectSumGpu::add_long_range_forces() const {
@@ -64,9 +71,17 @@ void DipolarDirectSumGpu::add_long_range_forces() const {
   auto const torques_device = gpu.get_particle_torques_device();
   auto const positions_device = gpu.get_particle_positions_device();
   auto const dipoles_device = gpu.get_particle_dipoles_device();
+#ifdef DIPOLE_FIELD_TRACKING
+  auto const dipole_fields_device = gpu.get_particle_dip_fld_device();
+#endif
   DipolarDirectSum_kernel_wrapper_force(
-      static_cast<float>(prefactor), npart, positions_device, dipoles_device,
-      forces_device, torques_device, box, periodicity);
+      static_cast<float>(prefactor), npart, positions_device, dipoles_device
+#ifdef DIPOLE_FIELD_TRACKING
+      ,
+      dipole_fields_device
+#endif
+      ,
+      forces_device, torques_device, box, periodicity, n_replicas);
 }
 
 void DipolarDirectSumGpu::long_range_energy() const {

src/core/magnetostatics/dipolar_direct_sum_gpu.hpp

@@ -26,7 +26,8 @@
 #include "magnetostatics/actor.hpp"
 
 struct DipolarDirectSumGpu : public Dipoles::Actor<DipolarDirectSumGpu> {
-  DipolarDirectSumGpu(double prefactor);
+  int n_replicas;
+  DipolarDirectSumGpu(double prefactor, int n_replicas);
 
   void on_activation() const;
   void on_boxl_change() const {}

src/core/magnetostatics/dipolar_direct_sum_gpu_cuda.cu

@@ -56,14 +56,14 @@ __device__ inline void get_mi_vector_dds(float res[3], float const a[3],
   }
 }
 
-__device__ void dipole_ia_force(float pf, float const *r1, float const *r2,
-                                float const *dip1, float const *dip2, float *f1,
-                                float *torque1, float *torque2, float box_l[3],
-                                int periodic[3]) {
-  // Distance between particles
-  float dr[3];
-  get_mi_vector_dds(dr, r1, r2, box_l, periodic);
-
+__device__ void dipole_ia_force(float pf, float const dr[3], float const *dip1,
+                                float const *dip2, float *f1, float *torque1,
+                                float *torque2
+#ifdef DIPOLE_FIELD_TRACKING
+                                ,
+                                float *dip_fld1, float *dip_fld2
+#endif
+) {
   // Powers of distance
   auto const r_sq = scalar_product(dr, dr);
   auto const r_sq_inv = 1.0f / r_sq;
@@ -78,6 +78,18 @@ __device__ void dipole_ia_force(float pf, float const *r1, float const *r2,
   auto const pe3 = scalar_product(dip2, dr);
   auto const pe4 = 3.0f * r5_inv;
 
+#ifdef DIPOLE_FIELD_TRACKING
+  auto const rep1 = pe3 * pe4;
+  auto const rep2 = pe2 * pe4;
+  dip_fld1[0] = pf * (rep1 * dr[0] - dip2[0] * r3_inv);
+  dip_fld1[1] = pf * (rep1 * dr[1] - dip2[1] * r3_inv);
+  dip_fld1[2] = pf * (rep1 * dr[2] - dip2[2] * r3_inv);
+
+  dip_fld2[0] = pf * (rep2 * dr[0] - dip1[0] * r3_inv);
+  dip_fld2[1] = pf * (rep2 * dr[1] - dip1[1] * r3_inv);
+  dip_fld2[2] = pf * (rep2 * dr[2] - dip1[2] * r3_inv);
+#endif
+
   // Force
   auto const aa = pe4 * pe1;
   auto const bb = -15.0f * pe2 * pe3 * r7_inv;
@@ -135,56 +147,156 @@ __device__ float dipole_ia_energy(float pf, float const *r1, float const *r2,
 // LCOV_EXCL_STOP
 
 __global__ void DipolarDirectSum_kernel_force(float pf, unsigned int n,
-                                              float *pos, float *dip, float *f,
-                                              float *torque, float box_l[3],
-                                              int periodic[3]) {
-
-  auto const i = blockIdx.x * blockDim.x + threadIdx.x;
-
+                                              const float *pos, const float *dip
+#ifdef DIPOLE_FIELD_TRACKING
+                                              ,
+                                              float *dip_fld
+#endif
+                                              ,
+                                              float *f, float *torque,
+                                              const float box_l[3],
+                                              const int periodic[3],
+                                              int n_replicas) {
+  unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;
   if (i >= n)
     return;
 
-  // Kahan summation based on the wikipedia article
-  // Force
-  float fi[3], fsum[3], tj[3];
+  // Load particle i
+  float xi = pos[3 * i + 0], yi = pos[3 * i + 1], zi = pos[3 * i + 2];
 
-  // Torque
-  float ti[3], tsum[3];
+  const float *mi = dip + 3 * i;
 
-  // There is one thread per particle. Each thread computes interactions
-  // with particles whose id is smaller than the thread id.
-  // The force and torque of all the interaction partners of the current thread
-  // is atomically added to global results at once.
-  // The result for the particle id equal to the thread id is atomically added
-  // to global memory at the end.
+  // Per‐thread accumulators
+  float fsum[3] = {0.0f, 0.0f, 0.0f};
+  float tsum[3] = {0.0f, 0.0f, 0.0f};
+#ifdef DIPOLE_FIELD_TRACKING
+  float dfsum[3] = {0.0f, 0.0f, 0.0f};
+#endif
+
+  float fi[3], ti1[3], ti2[3];
+#ifdef DIPOLE_FIELD_TRACKING
+  float dfi[3], dfj[3];
+#endif
 
-  // Clear summation vars
-  for (unsigned int j = 0; j < 3; j++) {
-    // Force
-    fsum[j] = 0;
-    // Torque
-    tsum[j] = 0;
+  // --- Minimal‐image branch (no replicas) ---
+  if (n_replicas == 0) {
+    for (unsigned int j = i + 1; j < n; ++j) {
+      // Compute MIC vector
+      float d0[3];
+      get_mi_vector_dds(d0, (float[]){xi, yi, zi}, pos + 3 * j, box_l,
+                        periodic);
+      // Compute force/torque
+      dipole_ia_force(pf, d0, mi, dip + 3 * j, fi, ti1, ti2
+#ifdef DIPOLE_FIELD_TRACKING
+                      ,
+                      dfi, dfj
+#endif
+      );
+
+      // Deposit on j (atomic)
+      for (int k = 0; k < 3; ++k) {
+        atomicAdd(f + 3 * j + k, -fi[k]);
+        atomicAdd(torque + 3 * j + k, ti2[k]);
+#ifdef DIPOLE_FIELD_TRACKING
+        atomicAdd(dip_fld + 3 * j + k, dfj[k]);
+#endif
+      }
+      // Accumulate on i
+      for (int k = 0; k < 3; ++k) {
+        fsum[k] += fi[k];
+        tsum[k] += ti1[k];
+#ifdef DIPOLE_FIELD_TRACKING
+        dfsum[k] += dfi[k];
+#endif
+      }
+    }
+    // Flush i’s result
+    for (int k = 0; k < 3; ++k) {
+      atomicAdd(f + 3 * i + k, fsum[k]);
+      atomicAdd(torque + 3 * i + k, tsum[k]);
+#ifdef DIPOLE_FIELD_TRACKING
+      atomicAdd(dip_fld + 3 * i + k, dfsum[k]);
+#endif
+    }
+    return;
+  }
+  // --- Replica branch (n_replicas > 0) ---
+  int nx = periodic[0] * n_replicas;
+  int ny = periodic[1] * n_replicas;
+  int nz = periodic[2] * n_replicas;
+
+  // Self‐images of i (exclude primary)
+  for (int dx = -nx; dx <= nx; ++dx) {
+    for (int dy = -ny; dy <= ny; ++dy) {
+      for (int dz = -nz; dz <= nz; ++dz) {
+        if (dx == 0 && dy == 0 && dz == 0)
+          continue;
+        float dr[3] = {static_cast<float>(dx) * box_l[0],
+                       static_cast<float>(dy) * box_l[1],
+                       static_cast<float>(dz) * box_l[2]};
+        dipole_ia_force(pf, dr, mi, mi, fi, ti1, ti2
+#ifdef DIPOLE_FIELD_TRACKING
+                        ,
+                        dfi, dfj
+#endif
+        );
+        for (int k = 0; k < 3; ++k) {
+          fsum[k] += fi[k];
+          tsum[k] += ti1[k];
+#ifdef DIPOLE_FIELD_TRACKING
+          dfsum[k] += dfi[k];
+#endif
+        }
+      }
+    }
   }
 
-  for (unsigned int j = i + 1; j < n; j++) {
-    dipole_ia_force(pf, pos + 3 * i, pos + 3 * j, dip + 3 * i, dip + 3 * j, fi,
-                    ti, tj, box_l, periodic);
-    for (unsigned int k = 0; k < 3; k++) {
-      // Add rhs to global memory
-      atomicAdd(f + 3 * j + k, -fi[k]);
-      atomicAdd((torque + 3 * j + k), tj[k]);
-      tsum[k] += ti[k];
-      fsum[k] += fi[k];
+  // Pairwise i <-> j over all images
+  for (unsigned int j = i + 1; j < n; ++j) {
+    float xj = pos[3 * j + 0], yj = pos[3 * j + 1], zj = pos[3 * j + 2];
+    const float *mj = dip + 3 * j;
+
+    for (int dx = -nx; dx <= nx; ++dx) {
+      for (int dy = -ny; dy <= ny; ++dy) {
+        for (int dz = -nz; dz <= nz; ++dz) {
+          float dr[3] = {(xi - xj) + static_cast<float>(dx) * box_l[0],
+                         (yi - yj) + static_cast<float>(dy) * box_l[1],
+                         (zi - zj) + static_cast<float>(dz) * box_l[2]};
+          dipole_ia_force(pf, dr, mi, mj, fi, ti1, ti2
+#ifdef DIPOLE_FIELD_TRACKING
+                          ,
+                          dfi, dfj
+#endif
+          );
+
+          // Deposit on j (atomic)
+          for (int k = 0; k < 3; ++k) {
+            atomicAdd(f + 3 * j + k, -fi[k]);
+            atomicAdd(torque + 3 * j + k, ti2[k]);
+#ifdef DIPOLE_FIELD_TRACKING
+            atomicAdd(dip_fld + 3 * j + k, dfj[k]);
+#endif
+            fsum[k] += fi[k];
+            tsum[k] += ti1[k];
+#ifdef DIPOLE_FIELD_TRACKING
+            dfsum[k] += dfi[k];
+#endif
+          }
+        }
+      }
     }
   }
 
-  // Add the left hand side result to global memory
-  for (int j = 0; j < 3; j++) {
-    atomicAdd(f + 3 * i + j, fsum[j]);
-    atomicAdd(torque + 3 * i + j, tsum[j]);
+  // Add i’s total to global memory **atomically**,
+  // in case any other asynchronous update might race
+  for (int k = 0; k < 3; ++k) {
+    atomicAdd(f + 3 * i + k, fsum[k]);
+    atomicAdd(torque + 3 * i + k, tsum[k]);
+#ifdef DIPOLE_FIELD_TRACKING
+    atomicAdd(dip_fld + 3 * i + k, dfsum[k]);
+#endif
   }
 }
-
 // LCOV_EXCL_START
 __device__ void dds_sumReduction(float *input, float *sum) {
   auto const tid = static_cast<int>(threadIdx.x);
@@ -245,8 +357,16 @@ inline void copy_box_data(float **box_l_gpu, int **periodic_gpu,
 }
 
 void DipolarDirectSum_kernel_wrapper_force(float k, unsigned int n, float *pos,
-                                           float *dip, float *f, float *torque,
-                                           float box_l[3], int periodic[3]) {
+                                           float *dip
+
+#ifdef DIPOLE_FIELD_TRACKING
+                                           ,
+                                           float *dip_fld
+#endif
+                                           ,
+                                           float *f, float *torque,
+                                           float box_l[3], int periodic[3],
+                                           int n_replicas) {
 
   unsigned int const bs = 64;
   dim3 grid(1, 1, 1);
@@ -267,8 +387,13 @@ void DipolarDirectSum_kernel_wrapper_force(float k, unsigned int n, float *pos,
   int *periodic_gpu;
   copy_box_data(&box_l_gpu, &periodic_gpu, box_l, periodic);
 
-  KERNELCALL(DipolarDirectSum_kernel_force, grid, block, k, n, pos, dip, f,
-             torque, box_l_gpu, periodic_gpu);
+  KERNELCALL(DipolarDirectSum_kernel_force, grid, block, k, n, pos, dip
+#ifdef DIPOLE_FIELD_TRACKING
+             ,
+             dip_fld
+#endif
+             ,
+             f, torque, box_l_gpu, periodic_gpu, n_replicas);
   cudaFree(box_l_gpu);
   cudaFree(periodic_gpu);
 }