Move pair potentials and change to traditional Functor

kernels/CMakeLists.txt

@@ -1,7 +1,7 @@
 project(${NAME}_cuda_lib CUDA CXX)
 
 set(HEADER_FILES
-    pair_potentials.cuh
+    potentials/pair_potentials.cuh
     forces.cuh
 )
 set(SOURCE_FILES

kernels/pair_potentials.cuh

@@ -1,91 +0,0 @@
-#ifndef POTENTIALS_CUH
-#define POTENTIALS_CUH
-
-#include "precision.hpp"
-#include "vec3.h"
-
-#ifdef __CUDACC__
-#define CUDA_CALLABLE __host__ __device__
-#else
-#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}};
-  };
-};
-
-/**
- * Abstract implementation of a Pair Potential.
- * Pair potentials are potentials which depend solely on the distance
- * between two particles. These do not include multi-body potentials such as
- * EAM
- *
- */
-struct PairPotential {
-  real m_rcutoffsq;
-
-CUDA_CALLABLE  PairPotential(real rcutoff) : m_rcutoffsq(rcutoff * rcutoff) {};
-#ifdef __CUDACC__
-  CUDA_CALLABLE ~PairPotential();
-#else
-  virtual ~PairPotential() = 0;
-#endif
-
-  /**
-   * Calculate the force and energy for a specific atom pair based on a
-   * displacement vector r.
-   */
-  CUDA_CALLABLE virtual ForceAndEnergy calc_force_and_energy(Vec3<real> r) = 0;
-};
-
-/**
- * Calculate the Lennard-Jones energy and force for the current particle pair
- * described by displacement vector r
- */
-struct LennardJones : PairPotential {
-  real m_epsilon;
-  real m_sigma;
-
-  CUDA_CALLABLE LennardJones(real sigma, real epsilon, real rcutoff)
-      : PairPotential(rcutoff), m_epsilon(epsilon), m_sigma(sigma) {};
-
-  CUDA_CALLABLE ForceAndEnergy calc_force_and_energy(Vec3<real> r) {
-    real rmagsq = r.squared_norm2();
-    if (rmagsq < this->m_rcutoffsq && rmagsq > 0.0) {
-      real inv_rmag = 1 / std::sqrt(rmagsq);
-
-      // Pre-Compute the terms (doing this saves on multiple devisions/pow
-      // function call)
-      real sigma_r = m_sigma * inv_rmag;
-      real sigma_r6 = sigma_r * sigma_r * sigma_r * sigma_r * sigma_r * sigma_r;
-      real sigma_r12 = sigma_r6 * sigma_r6;
-
-      // Get the energy
-      real energy = 4.0 * m_epsilon * (sigma_r12 - sigma_r6);
-
-      // Get the force vector
-      real force_mag =
-          4.0 * m_epsilon *
-          (12.0 * sigma_r12 * inv_rmag - 6.0 * sigma_r6 * inv_rmag);
-      Vec3<real> force = r.scale(force_mag * inv_rmag);
-
-      return {energy, force};
-
-    } else {
-      return ForceAndEnergy::zero();
-    }
-  };
-
-  CUDA_CALLABLE inline ~LennardJones(){};
-};
-
-inline PairPotential::~PairPotential() {};
-#endif

kernels/potentials/pair_potentials.cuh

@@ -0,0 +1,118 @@
+#ifndef POTENTIALS_CUH
+#define POTENTIALS_CUH
+
+#include "precision.hpp"
+#include "vec3.h"
+#include <cmath>
+#include <cstdio>
+#include <variant>
+
+#ifdef __CUDACC__
+#define CUDA_CALLABLE __host__ __device__
+#else
+#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
+ */
+struct LennardJones {
+  real m_sigma;
+  real m_epsilon;
+  real m_rcutoffsq;
+
+  CUDA_CALLABLE LennardJones(real sigma, real epsilon, real rcutoff) {
+    m_sigma = sigma;
+    m_epsilon = epsilon;
+    m_rcutoffsq = rcutoff * rcutoff;
+  };
+
+  CUDA_CALLABLE ForceAndEnergy 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);
+
+      // Pre-Compute the terms (doing this saves on multiple devisions/pow
+      // function call)
+      real sigma_r = m_sigma * inv_rmag;
+      real sigma_r6 = sigma_r * sigma_r * sigma_r * sigma_r * sigma_r * sigma_r;
+      real sigma_r12 = sigma_r6 * sigma_r6;
+
+      // Get the energy
+      real energy = 4.0 * m_epsilon * (sigma_r12 - sigma_r6);
+
+      // Get the force vector
+      real force_mag =
+          4.0 * m_epsilon *
+          (12.0 * sigma_r12 * inv_rmag - 6.0 * sigma_r6 * inv_rmag);
+      Vec3<real> force = r.scale(force_mag * inv_rmag);
+
+      return {energy, force};
+
+    } else {
+      return ForceAndEnergy::zero();
+    }
+  };
+};
+
+/**
+ * Calculate the Morse potential energy and force for the current particle pair
+ * described by displacement vector r
+ */
+struct Morse {
+  real m_D;         // Depth of the potential well
+  real m_a;         // Width of the potential
+  real m_r0;        // Equilibrium bond distance
+  real m_rcutoffsq; // Cutoff distance squared
+
+  CUDA_CALLABLE Morse(real D, real a, real r0, real rcutoff) {
+    m_D = D;
+    m_a = a;
+    m_r0 = r0;
+    m_rcutoffsq = rcutoff * rcutoff;
+  };
+
+  CUDA_CALLABLE ForceAndEnergy calc_force_and_energy(Vec3<real> r) {
+    real rmagsq = r.squared_norm2();
+    if (rmagsq < m_rcutoffsq && rmagsq > 0.0) {
+      real rmag = sqrt(rmagsq);
+      real dr = rmag - m_r0;
+
+      // Compute exponentials
+      real exp_a_dr = exp(-m_a * dr);
+      real exp_2a_dr = exp_a_dr * exp_a_dr;
+
+      // Energy: V(r) = D * (exp(-2a(r - r0)) - 2*exp(-a(r - r0)))
+      real energy = m_D * (exp_2a_dr - 2.0 * exp_a_dr);
+
+      // Force magnitude: F(r) = 2aD * (exp(-2a(r - r0)) - exp(-a(r - r0)))
+      real force_mag = 2.0 * m_a * m_D * (exp_2a_dr - exp_a_dr);
+
+      // Direction: normalized vector
+      Vec3<real> force = r.scale(force_mag / rmag);
+
+      return {energy, force};
+
+    } else {
+      return ForceAndEnergy::zero();
+    }
+  };
+};
+
+// Variant type for storing pair potential types
+using PairPotentials = std::variant<LennardJones, Morse>;
+
+#endif

tests/cuda_unit_tests/test_potential.cu

@@ -1,4 +1,4 @@
-#include "pair_potentials.cuh"
+#include "potentials/pair_potentials.cuh"
 #include "precision.hpp"
 #include "gtest/gtest.h"
 #include <cmath>
@@ -69,6 +69,7 @@ __global__ void lennard_jones_test_kernel(TestResults *results) {
     auto result = lj.calc_force_and_energy(r);
     results->energy_values[2] = result.energy;
     results->force_values[2] = result.force;
+
     results->at_minimum_pass =
         (fabs(result.energy + epsilon) < tolerance) &&
         vec3_near(Vec3<real>{0.0, 0.0, 0.0}, result.force, tolerance);

tests/unit_tests/test_potential.cpp

@@ -1,4 +1,4 @@
-#include "pair_potentials.cuh"
+#include "potentials/pair_potentials.cuh"
 #include "precision.hpp"
 #include "gtest/gtest.h"
 #include <cmath>