From 2d948a7e76cc434ad7bfecd9d802aebec04df2a4 Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Wed, 10 Sep 2025 22:03:11 -0400 Subject: [PATCH] Move pair potentials and change to traditional Functor --- kernels/CMakeLists.txt | 2 +- kernels/pair_potentials.cuh | 91 ------------------ kernels/potentials/pair_potentials.cuh | 118 ++++++++++++++++++++++++ tests/cuda_unit_tests/test_potential.cu | 3 +- tests/unit_tests/test_potential.cpp | 2 +- 5 files changed, 122 insertions(+), 94 deletions(-) delete mode 100644 kernels/pair_potentials.cuh create mode 100644 kernels/potentials/pair_potentials.cuh diff --git a/kernels/CMakeLists.txt b/kernels/CMakeLists.txt index fac4474..bef60a3 100644 --- a/kernels/CMakeLists.txt +++ b/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 diff --git a/kernels/pair_potentials.cuh b/kernels/pair_potentials.cuh deleted file mode 100644 index d6c97ac..0000000 --- a/kernels/pair_potentials.cuh +++ /dev/null @@ -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 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 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 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 force = r.scale(force_mag * inv_rmag); - - return {energy, force}; - - } else { - return ForceAndEnergy::zero(); - } - }; - - CUDA_CALLABLE inline ~LennardJones(){}; -}; - -inline PairPotential::~PairPotential() {}; -#endif diff --git a/kernels/potentials/pair_potentials.cuh b/kernels/potentials/pair_potentials.cuh new file mode 100644 index 0000000..537b03c --- /dev/null +++ b/kernels/potentials/pair_potentials.cuh @@ -0,0 +1,118 @@ +#ifndef POTENTIALS_CUH +#define POTENTIALS_CUH + +#include "precision.hpp" +#include "vec3.h" +#include +#include +#include + +#ifdef __CUDACC__ +#define CUDA_CALLABLE __host__ __device__ +#else +#define CUDA_CALLABLE +#endif + +/** + * Result struct for the Pair Potential + */ +struct ForceAndEnergy { + real energy; + Vec3 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 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 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 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 force = r.scale(force_mag / rmag); + + return {energy, force}; + + } else { + return ForceAndEnergy::zero(); + } + }; +}; + +// Variant type for storing pair potential types +using PairPotentials = std::variant; + +#endif diff --git a/tests/cuda_unit_tests/test_potential.cu b/tests/cuda_unit_tests/test_potential.cu index cab3216..9511ea5 100644 --- a/tests/cuda_unit_tests/test_potential.cu +++ b/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 @@ -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{0.0, 0.0, 0.0}, result.force, tolerance); diff --git a/tests/unit_tests/test_potential.cpp b/tests/unit_tests/test_potential.cpp index a94f022..d6bf23b 100644 --- a/tests/unit_tests/test_potential.cpp +++ b/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