Add a kernel_config to calculate blocks and threads for launching kernels

kernels/CMakeLists.txt

@@ -3,12 +3,14 @@ project(${NAME}_cuda_lib CUDA CXX)
 set(HEADER_FILES
     potentials/pair_potentials.cuh
     forces.cuh
+    kernel_config.cuh
 )
 set(SOURCE_FILES
+    kernel_config.cu
 )
 
 # The library contains header and source files.
-add_library(${NAME}_cuda_lib INTERFACE
+add_library(${NAME}_cuda_lib STATIC
     ${SOURCE_FILES}
     ${HEADER_FILES}
 )

kernels/forces.cuh

@@ -1,74 +1,75 @@
 #ifndef FORCES_CUH
 #define FORCES_CUH
+#include "kernel_config.cuh"
 #include "potentials/pair_potentials.cuh"
 #include "precision.hpp"
 #include <cstdio>
-#include <type_traits>
-#include <variant>
+#include <cuda_runtime.h>
 #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));
+inline void reset_forces_and_energies(int n_particles,
+                                      float4 *forces_energies) {
+  cudaMemset(forces_energies, 0, n_particles * sizeof(float4));
 }
 
 template <typename PotentialType>
-__global__ void calc_forces_and_energies(real *xs, real *forces, real *energies,
+__global__ void calc_forces_and_energies(float4 *pos, float4 *force_energies,
                                          int n_particles, real *box_len,
                                          PotentialType potential) {
-  int i = blockIdx.x * blockDim.x + threadIdx.x;
+
+  int i = get_thread_id();
 
   if (i < n_particles) {
-    real xi = xs[3 * i];
-    real yi = xs[3 * i + 1];
-    real zi = xs[3 * i + 2];
+    float4 my_pos = pos[i]; // Loads 16 bytes in one transaction
+    real xi = my_pos.x;
+    real yi = my_pos.y;
+    real zi = my_pos.z;
+
+    real total_fx = 0, total_fy = 0, total_fz = 0, total_energy = 0;
 
     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;
+        float4 other_pos = pos[j];
+        real dx = xi - other_pos.x;
+        real dy = yi - other_pos.y;
+        real dz = zi - other_pos.z;
 
         // 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;
+        float4 sol = potential.calc_force_and_energy({dx, dy, dz});
+        total_fx += sol.x;
+        total_fy += sol.y;
+        total_fz += sol.z;
+        total_energy += sol.w;
       }
     }
+
+    force_energies[i] = make_float4(total_fx, total_fy, total_fz, total_energy);
   }
 }
-
-inline void launch_force_kernels(real *xs, real *forces, real *energies,
+inline void launch_force_kernels(float4 *xs, float4 *force_energies,
                                  int n_particles, real *box_len,
                                  std::vector<PairPotentials> potentials,
-                                 int grid_size, int block_size) {
+                                 dim3 blocks, dim3 threads_per_block) {
 
-  reset_forces_and_energies(n_particles, forces, energies);
+  reset_forces_and_energies(n_particles, force_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);
+          calc_forces_and_energies<PotentialType>
+              <<<blocks, threads_per_block>>>(xs, force_energies, n_particles,
+                                              box_len, potential);
         },
         potential);
     cudaDeviceSynchronize();
   }
 }
 } // namespace CAC
-
 #endif

kernels/kernel_config.cu

@@ -0,0 +1,73 @@
+#include "kernel_config.cuh"
+#include <algorithm>
+#include <cstdio>
+
+size_t KernelConfig::total_threads() const {
+  return (size_t)blocks.x * blocks.y * blocks.z * threads.x * threads.y *
+         threads.z;
+}
+
+void KernelConfig::print() const {
+  printf("Grid: (%u, %u, %u), Block: (%u, %u, %u), Total threads: %zu\n",
+         blocks.x, blocks.y, blocks.z, threads.x, threads.y, threads.z,
+         total_threads());
+}
+
+KernelConfig get_launch_config(size_t n_elements, int threads_per_block,
+                               int max_blocks_per_dim) {
+
+  // Ensure threads_per_block is valid
+  threads_per_block = std::min(threads_per_block, 1024);
+  threads_per_block = std::max(threads_per_block, 32);
+
+  // Calculate total blocks needed
+  size_t total_blocks =
+      (n_elements + threads_per_block - 1) / threads_per_block;
+
+  dim3 threads(threads_per_block);
+  dim3 blocks;
+
+  if (total_blocks <= max_blocks_per_dim) {
+    // Simple 1D grid
+    blocks = dim3(total_blocks);
+  } else {
+    // Use 2D grid
+    blocks.x = max_blocks_per_dim;
+    blocks.y = (total_blocks + max_blocks_per_dim - 1) / max_blocks_per_dim;
+
+    // If still too big, use 3D grid
+    if (blocks.y > max_blocks_per_dim) {
+      blocks.y = max_blocks_per_dim;
+      blocks.z =
+          (total_blocks + (size_t)max_blocks_per_dim * max_blocks_per_dim - 1) /
+          ((size_t)max_blocks_per_dim * max_blocks_per_dim);
+    }
+  }
+
+  return KernelConfig(blocks, threads);
+}
+
+int get_optimal_block_size(int device_id) {
+  cudaDeviceProp prop;
+  cudaGetDeviceProperties(&prop, device_id);
+
+  // Use a fraction of max threads per block for better occupancy
+  // Typically 256 or 512 work well for most kernels
+  if (prop.maxThreadsPerBlock >= 1024) {
+    return 256; // Good balance of occupancy and register usage
+  } else if (prop.maxThreadsPerBlock >= 512) {
+    return 256;
+  } else {
+    return prop.maxThreadsPerBlock / 2;
+  }
+}
+
+KernelConfig get_launch_config_advanced(size_t n_elements, int device_id) {
+  cudaDeviceProp prop;
+  cudaGetDeviceProperties(&prop, device_id);
+
+  int threads_per_block = get_optimal_block_size(device_id);
+  int max_blocks_per_dim = prop.maxGridSize[0];
+
+  return get_launch_config(n_elements, threads_per_block, max_blocks_per_dim);
+}

kernels/kernel_config.cuh

@@ -0,0 +1,83 @@
+#ifndef KERNEL_CONFIG_CUH
+#define KERNEL_CONFIG_CUH
+#include <cstdio>
+#include <cuda_runtime.h>
+
+/**
+ * Structure to hold grid launch configuration
+ */
+struct KernelConfig {
+  dim3 blocks;
+  dim3 threads;
+
+  // Convenience constructor
+  KernelConfig(dim3 b, dim3 t) : blocks(b), threads(t) {}
+
+  // Total number of threads launched
+  size_t total_threads() const;
+
+  // Print configuration for debugging
+  void print() const;
+};
+
+/**
+ * Calculate optimal CUDA launch configuration for 1D problem
+ *
+ * @param n_elements Number of elements to process
+ * @param threads_per_block Desired threads per block (default: 256)
+ * @param max_blocks_per_dim Maximum blocks per grid dimension (default: 65535)
+ * @return LaunchConfig with optimal grid and block dimensions
+ */
+KernelConfig get_launch_config(size_t n_elements, int threads_per_block = 256,
+                               int max_blocks_per_dim = 65535);
+
+/**
+ * Calculate 1D thread index for kernels launched with get_launch_config()
+ * Use this inside your CUDA kernels
+ */
+__device__ inline size_t get_thread_id() {
+  return (size_t)blockIdx.z * gridDim.x * gridDim.y * blockDim.x +
+         (size_t)blockIdx.y * gridDim.x * blockDim.x +
+         (size_t)blockIdx.x * blockDim.x + threadIdx.x;
+}
+
+/**
+ * Alternative version that takes grid dimensions as parameters
+ * Useful if you need the index calculation in multiple places
+ */
+__device__ inline size_t get_thread_id(dim3 gridDim, dim3 blockDim,
+                                       dim3 blockIdx, dim3 threadIdx) {
+  return (size_t)blockIdx.z * gridDim.x * gridDim.y * blockDim.x +
+         (size_t)blockIdx.y * gridDim.x * blockDim.x +
+         (size_t)blockIdx.x * blockDim.x + threadIdx.x;
+}
+
+/**
+ * GPU device properties helper - gets optimal block size for current device
+ */
+int get_optimal_block_size(int device_id = 0);
+
+/**
+ * Advanced version that considers device properties
+ */
+KernelConfig get_launch_config_advanced(size_t n_elements, int device_id = 0);
+
+// Example usage in your kernel:
+/*
+template <typename PotentialType>
+__global__ void calc_forces_and_energies(float4 *pos, float4 *force_energies,
+                                         size_t n_particles, real *box_len,
+                                         PotentialType potential) {
+
+    size_t i = get_thread_id();
+
+    if (i < n_particles) {
+        // Your existing force calculation code here...
+        float4 my_pos = pos[i];
+        // ... rest of kernel unchanged
+    }
+}
+
+*/
+
+#endif

kernels/potentials/pair_potentials.cuh

@@ -5,6 +5,7 @@
 #include "vec3.h"
 #include <cmath>
 #include <cstdio>
+#include <cuda_runtime.h>
 #include <variant>
 
 #ifdef __CUDACC__
@@ -13,18 +14,6 @@
 #define CUDA_CALLABLE
 #endif
 
-/**
- * Result struct for the Pair Potential
- */
-struct ForceAndEnergy {
-  real energy;
-  Vec3<real> force;
-
-  CUDA_CALLABLE inline static ForceAndEnergy zero() {
-    return {0.0, {0.0, 0.0, 0.0}};
-  };
-};
-
 /**
  * Calculate the Lennard-Jones energy and force for the current particle
  * pair described by displacement vector r
@@ -40,7 +29,7 @@ struct LennardJones {
     m_rcutoffsq = rcutoff * rcutoff;
   };
 
-  CUDA_CALLABLE ForceAndEnergy calc_force_and_energy(Vec3<real> r) {
+  CUDA_CALLABLE float4 calc_force_and_energy(Vec3<real> r) {
     real rmagsq = r.squared_norm2();
     if (rmagsq < m_rcutoffsq && rmagsq > 0.0) {
       real inv_rmag = 1 / sqrt(rmagsq);
@@ -60,10 +49,10 @@ struct LennardJones {
           (12.0 * sigma_r12 * inv_rmag - 6.0 * sigma_r6 * inv_rmag);
       Vec3<real> force = r.scale(force_mag * inv_rmag);
 
-      return {energy, force};
+      return make_float4(force.x, force.y, force.z, energy);
 
     } else {
-      return ForceAndEnergy::zero();
+      return make_float4(0.0f, 0.0f, 0.0f, 0.0f);
     }
   };
 };
@@ -85,7 +74,7 @@ struct Morse {
     m_rcutoffsq = rcutoff * rcutoff;
   };
 
-  CUDA_CALLABLE ForceAndEnergy calc_force_and_energy(Vec3<real> r) {
+  CUDA_CALLABLE float4 calc_force_and_energy(Vec3<real> r) {
     real rmagsq = r.squared_norm2();
     if (rmagsq < m_rcutoffsq && rmagsq > 0.0) {
       real rmag = sqrt(rmagsq);
@@ -104,10 +93,10 @@ struct Morse {
       // Direction: normalized vector
       Vec3<real> force = r.scale(force_mag / rmag);
 
-      return {energy, force};
+      return make_float4(force.x, force.y, force.z, energy);
 
     } else {
-      return ForceAndEnergy::zero();
+      return make_float4(0.0f, 0.0f, 0.0f, 0.0f);
     }
   };
 };

src/precision.hpp

@@ -1,15 +1,15 @@
 #ifndef PRECISION_H
 #define PRECISION_H
 
-#ifdef USE_FLOATS
+#ifdef USE_DOUBLE
 
 /*
- * If macro USE_FLOATS is set then the default type will be floating point
- * precision. Otherwise we use double precision by default
+ * If macro USE_DOUBLE is set then the default type will be double
+ * precision. Otherwise we use floats by default
  */
-typedef float real;
-#else
 typedef double real;
+#else
+typedef float real;
 #endif
 
 #endif

tests/CMakeLists.txt

@@ -10,5 +10,4 @@ if(NOT EXISTS ${GOOGLETEST_DIR})
 endif()
 
 add_subdirectory(lib/googletest)
-add_subdirectory(unit_tests)
 add_subdirectory(cuda_unit_tests)

tests/cuda_unit_tests/test_forces.cu

@@ -5,14 +5,12 @@
 
 // Include your header files
 #include "forces.cuh"
+#include "kernel_config.cuh"
 #include "potentials/pair_potentials.cuh"
 #include "precision.hpp"
 
 class CudaForceKernelTest : public ::testing::Test {
 protected:
-  const int GRID_SIZE = 1;
-  const int BLOCK_SIZE = 4;
-
   void SetUp() override {
     // Set up CUDA device
     cudaError_t err = cudaSetDevice(0);
@@ -55,33 +53,32 @@ protected:
   }
 
   // Helper function to run the force calculation kernel
-  std::pair<std::vector<real>, std::vector<real>>
-  run_force_calculation(int n_particles, const std::vector<real> &positions,
+  std::vector<float4>
+  run_force_calculation(int n_particles, const std::vector<float4> &positions,
                         const std::vector<real> &box_dimensions) {
-    std::vector<real> forces(3 * n_particles, 0.0);
-    std::vector<real> energies(n_particles, 0.0);
+    std::vector<float4> force_energies(n_particles,
+                                       make_float4(0.0, 0.0, 0.0, 0.0));
 
-    real *d_positions = allocateAndCopyToGPU(positions);
-    real *d_forces = allocateAndCopyToGPU(forces);
-    real *d_energies = allocateAndCopyToGPU(energies);
+    KernelConfig kernel_config = get_launch_config(n_particles);
+    float4 *d_positions = allocateAndCopyToGPU(positions);
+    float4 *d_force_energies = allocateAndCopyToGPU(force_energies);
     real *d_box_len = allocateAndCopyToGPU(box_dimensions);
 
     std::vector<PairPotentials> potentials = {LennardJones(1.0, 1.0, 3.0)};
-    CAC::launch_force_kernels(d_positions, d_forces, d_energies, n_particles,
-                              d_box_len, potentials, GRID_SIZE, BLOCK_SIZE);
+    CAC::launch_force_kernels(d_positions, d_force_energies, n_particles,
+                              d_box_len, potentials, kernel_config.blocks,
+                              kernel_config.threads);
 
     checkCudaError(cudaGetLastError(), "kernel launch");
     checkCudaError(cudaDeviceSynchronize(), "kernel execution");
 
-    std::vector<real> result_forces =
-        copyFromGPUAndFree(d_forces, 3 * n_particles);
-    std::vector<real> result_energies =
-        copyFromGPUAndFree(d_energies, n_particles);
+    std::vector<float4> result_force_energies =
+        copyFromGPUAndFree(d_force_energies, n_particles);
 
     checkCudaError(cudaFree(d_positions), "cudaFree positions");
     checkCudaError(cudaFree(d_box_len), "cudaFree box_len");
 
-    return {result_forces, result_energies};
+    return result_force_energies;
   }
 };
 
@@ -90,14 +87,14 @@ TEST_F(CudaForceKernelTest, BasicFunctionalityTest) {
   const real tolerance = 1e-5;
 
   // Set up test data - simple 2x2 grid of particles
-  std::vector<real> positions = {
-      0.0, 0.0, 0.0, // particle 0
-      0.5, 0.0, 0.0, // particle 1
+  std::vector<float4> positions = {
+      make_float4(0.0, 0.0, 0.0, 0.0), // particle 0
+      make_float4(0.5, 0.0, 0.0, 0.0), // particle 1
   };
 
   std::vector<real> box_dimensions = {10.0, 10.0, 10.0};
 
-  auto [result_forces, result_energies] =
+  auto result_force_energies =
       run_force_calculation(n_particles, positions, box_dimensions);
 
   // Verify results - forces should be non-zero and energies should be
@@ -105,17 +102,14 @@ TEST_F(CudaForceKernelTest, BasicFunctionalityTest) {
   bool has_nonzero_force = false;
   bool has_nonzero_energy = false;
 
-  for (int i = 0; i < 3 * n_particles; i++) {
-    if (std::abs(result_forces[i]) > tolerance) {
-      has_nonzero_force = true;
-      break;
-    }
-  }
-
   for (int i = 0; i < n_particles; i++) {
-    if (std::abs(result_energies[i]) > tolerance) {
+    if (std::abs(result_force_energies[i].x) > tolerance ||
+        std::abs(result_force_energies[i].y) > tolerance ||
+        std::abs(result_force_energies[i].z) > tolerance) {
+      has_nonzero_force = true;
+    }
+    if (std::abs(result_force_energies[i].w) > tolerance) {
       has_nonzero_energy = true;
-      break;
     }
   }
 
@@ -130,60 +124,61 @@ TEST_F(CudaForceKernelTest, PeriodicBoundaryConditionsTest) {
   const real tolerance = 1e-5;
 
   // Place particles near opposite edges of a small box
-  std::vector<real> positions = {
-      0.1, 0.0, 0.0, // particle 0 near left edge
-      4.9, 0.0, 0.0  // particle 1 near right edge
+  std::vector<float4> positions = {
+      make_float4(0.1, 0.0, 0.0, 0.0), // particle 0 near left edge
+      make_float4(4.9, 0.0, 0.0, 0.0)  // particle 1 near right edge
   };
   std::vector<real> box_dimensions = {5.0, 5.0, 5.0}; // Small box to test PBC
 
-  auto [result_forces, result_energies] =
+  auto result_force_energies =
       run_force_calculation(n_particles, positions, box_dimensions);
 
   // With PBC, particles should interact as if they're close (distance ~0.2)
   // rather than far apart (distance ~4.8)
-  EXPECT_GT(std::abs(result_forces[0]), tolerance)
+  EXPECT_GT(std::abs(result_force_energies[0].x), tolerance)
       << "Expected significant force due to PBC";
-  EXPECT_GT(std::abs(result_energies[0]), tolerance)
-      << "Expected significant energy due to PBC";
 }
 
 TEST_F(CudaForceKernelTest, SingleParticleTest) {
   const int n_particles = 1;
 
-  std::vector<real> positions = {0.0, 0.0, 0.0};
+  std::vector<float4> positions = {make_float4(0.0, 0.0, 0.0, 0.0)};
   std::vector<real> box_dimensions = {10.0, 10.0, 10.0};
 
-  auto [result_forces, result_energies] =
+  auto result_force_energies =
       run_force_calculation(n_particles, positions, box_dimensions);
   // Single particle should have zero force and energy
-  EXPECT_NEAR(result_forces[0], 0.0, 1e-10);
-  EXPECT_NEAR(result_forces[1], 0.0, 1e-10);
-  EXPECT_NEAR(result_forces[2], 0.0, 1e-10);
-  EXPECT_NEAR(result_energies[0], 0.0, 1e-10);
+  EXPECT_NEAR(result_force_energies[0].x, 0.0, 1e-10);
+  EXPECT_NEAR(result_force_energies[0].y, 0.0, 1e-10);
+  EXPECT_NEAR(result_force_energies[0].z, 0.0, 1e-10);
+  EXPECT_NEAR(result_force_energies[0].w, 0.0, 1e-10);
 }
 
 TEST_F(CudaForceKernelTest, ForceSymmetryTest) {
   const int n_particles = 2;
   const real tolerance = 1e-5;
 
-  std::vector<real> positions = {
-      0.0, 0.0, 0.0, // particle 0
-      1.5, 0.0, 0.0  // particle 1
+  std::vector<float4> positions = {
+      make_float4(0.0, 0.0, 0.0, 0.0), // particle 0
+      make_float4(1.5, 0.0, 0.0, 0.0)  // particle 1
   };
   std::vector<real> box_dimensions = {10.0, 10.0, 10.0};
 
-  auto [result_forces, result_energies] =
+  auto result_force_energies =
       run_force_calculation(n_particles, positions, box_dimensions);
 
   // Newton's third law: forces should be equal and opposite
-  EXPECT_NEAR(result_forces[0], -result_forces[3], tolerance)
+  EXPECT_NEAR(result_force_energies[0].x, -result_force_energies[1].x,
+              tolerance)
       << "Force x-components should be opposite";
-  EXPECT_NEAR(result_forces[1], -result_forces[4], tolerance)
+  EXPECT_NEAR(result_force_energies[0].y, -result_force_energies[1].y,
+              tolerance)
       << "Force y-components should be opposite";
-  EXPECT_NEAR(result_forces[2], -result_forces[5], tolerance)
+  EXPECT_NEAR(result_force_energies[0].z, -result_force_energies[1].z,
+              tolerance)
       << "Force z-components should be opposite";
 
   // Energies should be equal for symmetric particles
-  EXPECT_NEAR(result_energies[0], result_energies[1], tolerance)
+  EXPECT_NEAR(result_force_energies[0].w, result_force_energies[1].w, tolerance)
       << "Energies should be equal";
 }

tests/cuda_unit_tests/test_potential.cu

@@ -2,6 +2,7 @@
 #include "precision.hpp"
 #include "gtest/gtest.h"
 #include <cmath>
+#include <cstdio>
 #include <cuda_runtime.h>
 
 // Structure to hold test results from device
@@ -18,8 +19,7 @@ struct TestResults {
   bool near_cutoff_pass;
 
   // Additional result data for exact checks
-  real energy_values[10];
-  Vec3<real> force_values[10];
+  float4 force_energy_values[10];
 };
 
 // Check if two Vec3 values are close within tolerance
@@ -35,7 +35,7 @@ __global__ void lennard_jones_test_kernel(TestResults *results) {
   real sigma = 1.0;
   real epsilon = 1.0;
   real r_cutoff = 2.5;
-  real tolerance = 1e-10;
+  real tolerance = 1e-5;
 
   // Create LennardJones object on device
   LennardJones lj(sigma, epsilon, r_cutoff);
@@ -43,87 +43,78 @@ __global__ void lennard_jones_test_kernel(TestResults *results) {
   // Zero Distance Test
   {
     Vec3<real> r = {0.0, 0.0, 0.0};
-    auto result = lj.calc_force_and_energy(r);
-    results->energy_values[0] = result.energy;
-    results->force_values[0] = result.force;
+    float4 result = lj.calc_force_and_energy(r);
+    results->force_energy_values[0] = result;
     results->zero_distance_pass =
-        (result.energy == 0.0) &&
-        vec3_near(Vec3<real>{0.0, 0.0, 0.0}, result.force, tolerance);
+        (result.w == 0.0) &&
+        vec3_near(Vec3<real>{0.0, 0.0, 0.0},
+                  Vec3<real>{result.x, result.y, result.z}, tolerance);
   }
 
   // Beyond Cutoff Test
   {
     Vec3<real> r = {3.0, 0.0, 0.0};
-    auto result = lj.calc_force_and_energy(r);
-    results->energy_values[1] = result.energy;
-    results->force_values[1] = result.force;
+    float4 result = lj.calc_force_and_energy(r);
+    results->force_energy_values[1] = result;
     results->beyond_cutoff_pass =
-        (result.energy == 0.0) &&
-        vec3_near(Vec3<real>{0.0, 0.0, 0.0}, result.force, tolerance);
+        (result.w == 0.0) &&
+        vec3_near(Vec3<real>{0.0, 0.0, 0.0},
+                  Vec3<real>{result.x, result.y, result.z}, tolerance);
   }
 
   // At Minimum Test
   {
     real min_dist = pow(2.0, 1.0 / 6.0) * sigma;
     Vec3<real> r = {min_dist, 0.0, 0.0};
-    auto result = lj.calc_force_and_energy(r);
-    results->energy_values[2] = result.energy;
-    results->force_values[2] = result.force;
+    float4 result = lj.calc_force_and_energy(r);
+    results->force_energy_values[2] = result;
 
     results->at_minimum_pass =
-        (fabs(result.energy + epsilon) < tolerance) &&
-        vec3_near(Vec3<real>{0.0, 0.0, 0.0}, result.force, tolerance);
+        (fabs(result.w + epsilon) < tolerance) &&
+        vec3_near(Vec3<real>{0.0, 0.0, 0.0},
+                  Vec3<real>{result.x, result.y, result.z}, tolerance);
   }
 
   // At Equilibrium Test
   {
     Vec3<real> r = {sigma, 0.0, 0.0};
-    auto result = lj.calc_force_and_energy(r);
-    results->energy_values[3] = result.energy;
-    results->force_values[3] = result.force;
-    results->at_equilibrium_pass = (fabs(result.energy) < tolerance) &&
-                                   (result.force.x > 0.0) &&
-                                   (fabs(result.force.y) < tolerance) &&
-                                   (fabs(result.force.z) < tolerance);
+    float4 result = lj.calc_force_and_energy(r);
+    results->force_energy_values[3] = result;
+    results->at_equilibrium_pass =
+        (fabs(result.w) < tolerance) && (result.x > 0.0) &&
+        (fabs(result.y) < tolerance) && (fabs(result.z) < tolerance);
   }
 
   // Repulsive Region Test
   {
-    Vec3<real> r = {0.8 * sigma, 0.0, 0.0};
-    auto result = lj.calc_force_and_energy(r);
-    results->energy_values[4] = result.energy;
-    results->force_values[4] = result.force;
-    results->repulsive_region_pass =
-        (result.energy > 0.0) && (result.force.x > 0.0);
+    Vec3<real> r = {0.8f * sigma, 0.0, 0.0};
+    float4 result = lj.calc_force_and_energy(r);
+    results->force_energy_values[4] = result;
+    results->repulsive_region_pass = (result.w > 0.0) && (result.x > 0.0);
   }
 
   // Attractive Region Test
   {
-    Vec3<real> r = {1.5 * sigma, 0.0, 0.0};
-    auto result = lj.calc_force_and_energy(r);
-    results->energy_values[5] = result.energy;
-    results->force_values[5] = result.force;
-    results->attractive_region_pass =
-        (result.energy < 0.0) && (result.force.x < 0.0);
+    Vec3<real> r = {1.5f * sigma, 0.0, 0.0};
+    float4 result = lj.calc_force_and_energy(r);
+    results->force_energy_values[5] = result;
+    results->attractive_region_pass = (result.w < 0.0) && (result.x < 0.0);
   }
 
   // Arbitrary Direction Test
   {
     Vec3<real> r = {1.0, 1.0, 1.0};
-    auto result = lj.calc_force_and_energy(r);
-    results->energy_values[6] = result.energy;
-    results->force_values[6] = result.force;
+    float4 result = lj.calc_force_and_energy(r);
+    results->force_energy_values[6] = result;
 
     real r_mag = sqrt(r.squared_norm2());
     Vec3<real> normalized_r = r.scale(1.0 / r_mag);
-    real force_dot_r = result.force.x * normalized_r.x +
-                       result.force.y * normalized_r.y +
-                       result.force.z * normalized_r.z;
+    real force_dot_r = result.x * normalized_r.x + result.y * normalized_r.y +
+                       result.z * normalized_r.z;
 
     results->arbitrary_direction_pass =
-        (force_dot_r < 0.0) &&
-        (fabs(result.force.x - result.force.y) < tolerance) &&
-        (fabs(result.force.y - result.force.z) < tolerance);
+        (force_dot_r < 0.0) && (fabs(result.x - result.y) < tolerance) &&
+        (fabs(result.y - result.z) < tolerance);
   }
 
   // Parameter Variation Test
@@ -135,34 +126,31 @@ __global__ void lennard_jones_test_kernel(TestResults *results) {
     LennardJones lj2(new_sigma, new_epsilon, new_r_cutoff);
 
     Vec3<real> r = {2.0, 0.0, 0.0};
-    auto result1 = lj.calc_force_and_energy(r);
-    auto result2 = lj2.calc_force_and_energy(r);
+    float4 result1 = lj.calc_force_and_energy(r);
+    float4 result2 = lj2.calc_force_and_energy(r);
 
-    results->energy_values[7] = result2.energy;
-    results->force_values[7] = result2.force;
+    results->force_energy_values[7] = result2;
 
-    results->parameter_variation_pass = (result1.energy != result2.energy) &&
-                                        (result1.force.x != result2.force.x);
+    results->parameter_variation_pass =
+        (result1.w != result2.w) && (result1.x != result2.x);
   }
 
   // Exact Value Check Test
   {
     LennardJones lj_exact(1.0, 1.0, 3.0);
     Vec3<real> r = {1.5, 0.0, 0.0};
-    auto result = lj_exact.calc_force_and_energy(r);
+    float4 result = lj_exact.calc_force_and_energy(r);
 
-    results->energy_values[8] = result.energy;
-    results->force_values[8] = result.force;
+    results->force_energy_values[8] = result;
 
     real expected_energy = 4.0 * (pow(1.0 / 1.5, 12) - pow(1.0 / 1.5, 6));
     real expected_force =
         24.0 * (pow(1.0 / 1.5, 6) - 2.0 * pow(1.0 / 1.5, 12)) / 1.5;
 
     results->exact_value_check_pass =
-        (fabs(result.energy - expected_energy) < tolerance) &&
-        (fabs(result.force.x + expected_force) < tolerance) &&
-        (fabs(result.force.y) < tolerance) &&
-        (fabs(result.force.z) < tolerance);
+        (fabs(result.w - expected_energy) < tolerance) &&
+        (fabs(result.x + expected_force) < tolerance) &&
+        (fabs(result.y) < tolerance) && (fabs(result.z) < tolerance);
   }
 
   // Near Cutoff Test
@@ -173,16 +161,18 @@ __global__ void lennard_jones_test_kernel(TestResults *results) {
     Vec3<real> r_inside = {inside_cutoff, 0.0, 0.0};
     Vec3<real> r_outside = {outside_cutoff, 0.0, 0.0};
 
-    auto result_inside = lj.calc_force_and_energy(r_inside);
-    auto result_outside = lj.calc_force_and_energy(r_outside);
+    float4 result_inside = lj.calc_force_and_energy(r_inside);
+    float4 result_outside = lj.calc_force_and_energy(r_outside);
 
-    results->energy_values[9] = result_inside.energy;
-    results->force_values[9] = result_inside.force;
+    results->force_energy_values[9] = result_inside;
 
     results->near_cutoff_pass =
-        (result_inside.energy != 0.0) && (result_inside.force.x != 0.0) &&
-        (result_outside.energy == 0.0) &&
-        vec3_near(Vec3<real>{0.0, 0.0, 0.0}, result_outside.force, tolerance);
+        (result_inside.w != 0.0) && (result_inside.x != 0.0) &&
+        (result_outside.w == 0.0) &&
+        vec3_near(
+            Vec3<real>{0.0, 0.0, 0.0},
+            Vec3<real>{result_outside.x, result_outside.y, result_outside.z},
+            tolerance);
   }
 }
 
@@ -250,44 +240,48 @@ TEST_F(LennardJonesCudaTest, DeviceZeroDistance) {
   auto results = runDeviceTests();
   EXPECT_TRUE(results.zero_distance_pass)
       << "Zero distance test failed on device. Energy: "
-      << results.energy_values[0] << ", Force: (" << results.force_values[0].x
-      << ", " << results.force_values[0].y << ", " << results.force_values[0].z
-      << ")";
+      << results.force_energy_values[0].w << ", Force: ("
+      << results.force_energy_values[0].x << ", "
+      << results.force_energy_values[0].y << ", "
+      << results.force_energy_values[0].z << ")";
 }
 
 TEST_F(LennardJonesCudaTest, DeviceBeyondCutoff) {
   auto results = runDeviceTests();
   EXPECT_TRUE(results.beyond_cutoff_pass)
       << "Beyond cutoff test failed on device. Energy: "
-      << results.energy_values[1];
+      << results.force_energy_values[1].w;
 }
 
 TEST_F(LennardJonesCudaTest, DeviceAtMinimum) {
   auto results = runDeviceTests();
   EXPECT_TRUE(results.at_minimum_pass)
       << "At minimum test failed on device. Energy: "
-      << results.energy_values[2];
+      << results.force_energy_values[2].w;
 }
 
 TEST_F(LennardJonesCudaTest, DeviceAtEquilibrium) {
   auto results = runDeviceTests();
   EXPECT_TRUE(results.at_equilibrium_pass)
       << "At equilibrium test failed on device. Energy: "
-      << results.energy_values[3] << ", Force x: " << results.force_values[3].x;
+      << results.force_energy_values[3].w
+      << ", Force x: " << results.force_energy_values[3].x;
 }
 
 TEST_F(LennardJonesCudaTest, DeviceRepulsiveRegion) {
   auto results = runDeviceTests();
   EXPECT_TRUE(results.repulsive_region_pass)
       << "Repulsive region test failed on device. Energy: "
-      << results.energy_values[4] << ", Force x: " << results.force_values[4].x;
+      << results.force_energy_values[4].w
+      << ", Force x: " << results.force_energy_values[4].x;
 }
 
 TEST_F(LennardJonesCudaTest, DeviceAttractiveRegion) {
   auto results = runDeviceTests();
   EXPECT_TRUE(results.attractive_region_pass)
       << "Attractive region test failed on device. Energy: "
-      << results.energy_values[5] << ", Force x: " << results.force_values[5].x;
+      << results.force_energy_values[5].w
+      << ", Force x: " << results.force_energy_values[5].x;
 }
 
 TEST_F(LennardJonesCudaTest, DeviceArbitraryDirection) {
@@ -306,12 +300,13 @@ TEST_F(LennardJonesCudaTest, DeviceExactValueCheck) {
   auto results = runDeviceTests();
   EXPECT_TRUE(results.exact_value_check_pass)
       << "Exact value check test failed on device. Energy: "
-      << results.energy_values[8] << ", Force x: " << results.force_values[8].x;
+      << results.force_energy_values[8].w
+      << ", Force x: " << results.force_energy_values[8].x;
 }
 
 TEST_F(LennardJonesCudaTest, DeviceNearCutoff) {
   auto results = runDeviceTests();
   EXPECT_TRUE(results.near_cutoff_pass)
       << "Near cutoff test failed on device. Inside energy: "
-      << results.energy_values[9];
+      << results.force_energy_values[9].w;
 }

tests/unit_tests/CMakeLists.txt

@@ -1,9 +0,0 @@
-include_directories(${gtest_SOURCE_DIR}/include ${gtest_SOURCE_DIR})
-
-add_executable(${NAME}_tests
-    test_potential.cpp
-)
-
-target_link_libraries(${NAME}_tests gtest gtest_main)
-target_link_libraries(${NAME}_tests ${CMAKE_PROJECT_NAME}_cuda_lib)
-add_test(NAME ${NAME}Tests COMMAND ${CMAKE_BINARY_DIR}/tests/unit_tests/${NAME}_tests)

tests/unit_tests/test_example.cpp

@@ -1,5 +0,0 @@
-#include "gtest/gtest.h"
-
-TEST(Example, Equals) {
-  EXPECT_EQ(1, 1);
-}

tests/unit_tests/test_potential.cpp

@@ -1,174 +0,0 @@
-#include "potentials/pair_potentials.cuh"
-#include "precision.hpp"
-#include "gtest/gtest.h"
-#include <cmath>
-
-class LennardJonesTest : public ::testing::Test {
-protected:
-  void SetUp() override {
-    // Default parameters
-    sigma = 1.0;
-    epsilon = 1.0;
-    r_cutoff = 2.5;
-
-    // Create default LennardJones object
-    lj = new LennardJones(sigma, epsilon, r_cutoff);
-  }
-
-  void TearDown() override { delete lj; }
-
-  real sigma;
-  real epsilon;
-  real r_cutoff;
-  LennardJones *lj;
-
-  // Helper function to compare Vec3 values with tolerance
-  void expect_vec3_near(const Vec3<real> &expected, const Vec3<real> &actual,
-                        real tolerance) {
-    EXPECT_NEAR(expected.x, actual.x, tolerance);
-    EXPECT_NEAR(expected.y, actual.y, tolerance);
-    EXPECT_NEAR(expected.z, actual.z, tolerance);
-  }
-};
-
-TEST_F(LennardJonesTest, ZeroDistance) {
-  // At zero distance, the calculation should return zero force and energy
-  Vec3<real> r = {0.0, 0.0, 0.0};
-  auto result = lj->calc_force_and_energy(r);
-
-  EXPECT_EQ(0.0, result.energy);
-  expect_vec3_near({0.0, 0.0, 0.0}, result.force, 1e-10);
-}
-
-TEST_F(LennardJonesTest, BeyondCutoff) {
-  // Distance beyond cutoff should return zero force and energy
-  Vec3<real> r = {3.0, 0.0, 0.0}; // 3.0 > r_cutoff (2.5)
-  auto result = lj->calc_force_and_energy(r);
-
-  EXPECT_EQ(0.0, result.energy);
-  expect_vec3_near({0.0, 0.0, 0.0}, result.force, 1e-10);
-}
-
-TEST_F(LennardJonesTest, AtMinimum) {
-  // The LJ potential has a minimum at r = 2^(1/6) * sigma
-  real min_dist = std::pow(2.0, 1.0 / 6.0) * sigma;
-  Vec3<real> r = {min_dist, 0.0, 0.0};
-  auto result = lj->calc_force_and_energy(r);
-
-  // At minimum, force should be close to zero
-  EXPECT_NEAR(-epsilon, result.energy, 1e-10);
-  expect_vec3_near({0.0, 0.0, 0.0}, result.force, 1e-10);
-}
-
-TEST_F(LennardJonesTest, AtEquilibrium) {
-  // At r = sigma, the energy should be zero and force should be repulsive
-  Vec3<real> r = {sigma, 0.0, 0.0};
-  auto result = lj->calc_force_and_energy(r);
-
-  EXPECT_NEAR(0.0, result.energy, 1e-10);
-  EXPECT_GT(result.force.x,
-            0.0); // Force should be repulsive (positive x-direction)
-  EXPECT_NEAR(0.0, result.force.y, 1e-10);
-  EXPECT_NEAR(0.0, result.force.z, 1e-10);
-}
-
-TEST_F(LennardJonesTest, RepulsiveRegion) {
-  // Test in the repulsive region (r < sigma)
-  Vec3<real> r = {0.8 * sigma, 0.0, 0.0};
-  auto result = lj->calc_force_and_energy(r);
-
-  // Energy should be positive and force should be repulsive
-  EXPECT_GT(result.energy, 0.0);
-  EXPECT_GT(result.force.x, 0.0); // Force should be repulsive
-}
-
-TEST_F(LennardJonesTest, AttractiveRegion) {
-  // Test in the attractive region (sigma < r < r_min)
-  Vec3<real> r = {1.5 * sigma, 0.0, 0.0};
-  auto result = lj->calc_force_and_energy(r);
-
-  // Energy should be negative and force should be attractive
-  EXPECT_LT(result.energy, 0.0);
-  EXPECT_LT(result.force.x,
-            0.0); // Force should be attractive (negative x-direction)
-}
-
-TEST_F(LennardJonesTest, ArbitraryDirection) {
-  // Test with a vector in an arbitrary direction
-  Vec3<real> r = {1.0, 1.0, 1.0};
-  auto result = lj->calc_force_and_energy(r);
-
-  // The force should be in the same direction as r but opposite sign
-  // (attractive region)
-  real r_mag = std::sqrt(r.squared_norm2());
-
-  // Calculate expected force direction (should be along -r)
-  Vec3<real> normalized_r = r.scale(1.0 / r_mag);
-  real force_dot_r = result.force.x * normalized_r.x +
-                     result.force.y * normalized_r.y +
-                     result.force.z * normalized_r.z;
-
-  // In this case, we're at r = sqrt(3) * sigma which is in attractive region
-  EXPECT_LT(force_dot_r, 0.0); // Force should be attractive
-
-  // Force should be symmetric in all dimensions for this vector
-  EXPECT_NEAR(result.force.x, result.force.y, 1e-10);
-  EXPECT_NEAR(result.force.y, result.force.z, 1e-10);
-}
-
-TEST_F(LennardJonesTest, ParameterVariation) {
-  // Test with different parameter values
-  real new_sigma = 2.0;
-  real new_epsilon = 0.5;
-  real new_r_cutoff = 5.0;
-
-  LennardJones lj2(new_sigma, new_epsilon, new_r_cutoff);
-
-  Vec3<real> r = {2.0, 0.0, 0.0};
-  auto result1 = lj->calc_force_and_energy(r);
-  auto result2 = lj2.calc_force_and_energy(r);
-
-  // Results should be different with different parameters
-  EXPECT_NE(result1.energy, result2.energy);
-  EXPECT_NE(result1.force.x, result2.force.x);
-}
-
-TEST_F(LennardJonesTest, ExactValueCheck) {
-  // Test with pre-calculated values for a specific case
-  LennardJones lj_exact(1.0, 1.0, 3.0);
-  Vec3<real> r = {1.5, 0.0, 0.0};
-  auto result = lj_exact.calc_force_and_energy(r);
-
-  // Pre-calculated values (you may need to adjust these based on your specific
-  // implementation)
-  real expected_energy =
-      4.0 * (std::pow(1.0 / 1.5, 12) - std::pow(1.0 / 1.5, 6));
-  real expected_force =
-      24.0 * (std::pow(1.0 / 1.5, 6) - 2.0 * std::pow(1.0 / 1.5, 12)) / 1.5;
-
-  EXPECT_NEAR(expected_energy, result.energy, 1e-10);
-  EXPECT_NEAR(-expected_force, result.force.x,
-              1e-10); // Negative because force is attractive
-  EXPECT_NEAR(0.0, result.force.y, 1e-10);
-  EXPECT_NEAR(0.0, result.force.z, 1e-10);
-}
-
-TEST_F(LennardJonesTest, NearCutoff) {
-  // Test behavior just inside and just outside the cutoff
-  real inside_cutoff = r_cutoff - 0.01;
-  real outside_cutoff = r_cutoff + 0.01;
-
-  Vec3<real> r_inside = {inside_cutoff, 0.0, 0.0};
-  Vec3<real> r_outside = {outside_cutoff, 0.0, 0.0};
-
-  auto result_inside = lj->calc_force_and_energy(r_inside);
-  auto result_outside = lj->calc_force_and_energy(r_outside);
-
-  // Inside should have non-zero values
-  EXPECT_NE(0.0, result_inside.energy);
-  EXPECT_NE(0.0, result_inside.force.x);
-
-  // Outside should be zero
-  EXPECT_EQ(0.0, result_outside.energy);
-  expect_vec3_near({0.0, 0.0, 0.0}, result_outside.force, 1e-10);
-}