Update CMakeFiles, add initial pair potential implementation and tests

.gitignore

@@ -1,5 +1,7 @@
 # Builds
 build/
+Debug/
+Testing/
 
 # Google Tests
 tests/lib/
@@ -7,3 +9,6 @@ tests/lib/
 # Jet Brains
 .idea/
 cmake-build-debug/
+
+# Cache dir
+.cache

CMakeLists.txt

@@ -1,5 +1,6 @@
 cmake_minimum_required(VERSION 3.9)
-project(cudaCAC LANGUAGES CUDA CXX)
+set(NAME "cudaCAC")
+project(${NAME} LANGUAGES CUDA CXX)
 
 enable_testing()
 set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
@@ -12,7 +13,6 @@ set(CMAKE_CUDA_ARCHITECTURES 61)
 set(CUDA_SEPARABLE_COMPILATION ON)
 
 # Add Vec3  as a dependency
-add_subdirectory(tests)
 include(FetchContent)
 FetchContent_Declare(Vec3
     GIT_REPOSITORY https://www.alexselimov.com/git/aselimov/Vec3.git
@@ -30,16 +30,17 @@ include_directories(/usr/local/cuda-12.8/include)
 
 add_subdirectory(src)
 add_subdirectory(kernels)
+add_subdirectory(tests)
 
-add_executable(${CMAKE_PROJECT_NAME} main.cpp)
+add_executable(${NAME} main.cpp)
 install(DIRECTORY src/ DESTINATION src/)
 
 
 target_link_libraries(
-    ${CMAKE_PROJECT_NAME} 
+    ${NAME} 
     PRIVATE
-    ${CMAKE_PROJECT_NAME}_lib 
+    ${NAME}_lib 
-    ${CMAKE_PROJECT_NAME}_cuda_lib 
+    ${NAME}_cuda_lib 
     
     ${CUDA_LIBRARIES}
 )

kernels/CMakeLists.txt

@@ -1,4 +1,4 @@
-project(${CMAKE_PROJECT_NAME}_cuda_lib CUDA CXX)
+project(${NAME}_cuda_lib CUDA CXX)
 
 set(HEADER_FILES
     hello_world.h
@@ -8,11 +8,9 @@ set(SOURCE_FILES
 )
 
 # The library contains header and source files.
-add_library(${CMAKE_PROJECT_NAME}_cuda_lib STATIC
+add_library(${NAME}_cuda_lib STATIC
     ${SOURCE_FILES}
     ${HEADER_FILES}
 )
 
-if(CMAKE_COMPILER_IS_GNUCXX)
+target_compile_options(${CMAKE_PROJECT_NAME}_cuda_lib PRIVATE -Wno-gnu-line-marker -Wno-pedantic)
-  target_compile_options(${CMAKE_PROJECT_NAME}_cuda_lib PRIVATE -Wno-gnu-line-marker)
-endif()

main.cpp

@@ -1,10 +1,10 @@
-#include "hello_world.h"
+#include "particle.hpp"
+#include "vec3.h"
 #include <iostream>
 
 int main() {
-  std::cout << "Starting CUDA example..." << std::endl; // Using endl to flush
+  Particle<float> test = {
-  check_cuda();
+      {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, 10};
-  launch_hello_cuda();
+  std::cout << test.pos.x << " " << test.pos.y << " " << test.pos.z;
-  std::cout << "Ending CUDA example" << std::endl; // Using endl to flush
   return 0;
 }

src/CMakeLists.txt

@@ -1,4 +1,4 @@
-project(${CMAKE_PROJECT_NAME}_lib CUDA CXX)
+project(${NAME}_lib CUDA CXX)
 
 set(HEADER_FILES
     particle.hpp
@@ -11,7 +11,7 @@ set(SOURCE_FILES
 )
 
 # The library contains header and source files.
-add_library(${CMAKE_PROJECT_NAME}_lib 
+add_library(${NAME}_lib 
     ${HEADER_FILES} 
     ${SOURCE_FILES}
 )

src/pair_potentials.cpp

@@ -1,6 +1,7 @@
-#include "potentials.hpp"
+#include "pair_potentials.hpp"
 #include <cmath>
 
+PairPotential::~PairPotential() {};
 /**
  * Calculate the Lennard-Jones energy and force for the current particle pair
  * described by displacement vector r
@@ -12,17 +13,16 @@ ForceAndEnergy LennardJones::calc_force_and_energy(Vec3<real> r) {
 
     // Pre-Compute the terms (doing this saves on multiple devisions/pow
     // function call)
-    real sigma_r = m_sigma / inv_rmag;
+    real sigma_r = m_sigma * inv_rmag;
-    real sigma_r5 = sigma_r * sigma_r * sigma_r * sigma_r * sigma_r;
+    real sigma_r6 = sigma_r * sigma_r * sigma_r * sigma_r * sigma_r * sigma_r;
-    real sigma_r6 = sigma_r5 * sigma_r;
-    real sigma_r11 = sigma_r5 * sigma_r5 * 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 * (6.0 * sigma_r5 - 12.0 * sigma_r11);
+    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};

src/pair_potentials.hpp

@@ -24,6 +24,9 @@ struct ForceAndEnergy {
 struct PairPotential {
   real m_rcutoffsq;
 
+  PairPotential(real rcutoff) : m_rcutoffsq(rcutoff * rcutoff) {};
+  virtual ~PairPotential() = 0;
+
   /**
    * Calculate the force and energy for a specific atom pair based on a
    * displacement vector r.
@@ -35,7 +38,12 @@ struct LennardJones : PairPotential {
   real m_epsilon;
   real m_sigma;
 
+  LennardJones(real sigma, real epsilon, real rcutoff)
+      : PairPotential(rcutoff), m_epsilon(epsilon), m_sigma(sigma) {};
+
   ForceAndEnergy calc_force_and_energy(Vec3<real> r);
+
+  ~LennardJones() {};
 };
 
 #endif

tests/unit_tests/CMakeLists.txt

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

tests/unit_tests/test_potential.cpp

@@ -1,4 +1,4 @@
-#include "potentials.hpp"
+#include "pair_potentials.hpp"
 #include "precision.hpp"
 #include "gtest/gtest.h"
 #include <cmath>
@@ -9,22 +9,22 @@ protected:
     // Default parameters
     sigma = 1.0;
     epsilon = 1.0;
-    rCutoff = 2.5;
+    r_cutoff = 2.5;
 
     // Create default LennardJones object
-    lj = new LennardJones(sigma, epsilon, rCutoff);
+    lj = new LennardJones(sigma, epsilon, r_cutoff);
   }
 
   void TearDown() override { delete lj; }
 
   real sigma;
   real epsilon;
-  real rCutoff;
+  real r_cutoff;
   LennardJones *lj;
 
   // Helper function to compare Vec3 values with tolerance
-  void expectVec3Near(const Vec3<real> &expected, const Vec3<real> &actual,
+  void expect_vec3_near(const Vec3<real> &expected, const Vec3<real> &actual,
-                      real tolerance) {
+                        real tolerance) {
     EXPECT_NEAR(expected.x, actual.x, tolerance);
     EXPECT_NEAR(expected.y, actual.y, tolerance);
     EXPECT_NEAR(expected.z, actual.z, tolerance);
@@ -33,36 +33,36 @@ protected:
 
 TEST_F(LennardJonesTest, ZeroDistance) {
   // At zero distance, the calculation should return zero force and energy
-  Vec3<real> r(0.0, 0.0, 0.0);
+  Vec3<real> r = {0.0, 0.0, 0.0};
   auto result = lj->calc_force_and_energy(r);
 
   EXPECT_EQ(0.0, result.energy);
-  expectVec3Near(Vec3<real>(0.0, 0.0, 0.0), result.force, 1e-10);
+  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 > rCutoff (2.5)
+  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);
-  expectVec3Near(Vec3<real>(0.0, 0.0, 0.0), result.force, 1e-10);
+  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);
+  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);
-  expectVec3Near(Vec3<real>(0.0, 0.0, 0.0), result.force, 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);
+  Vec3<real> r = {sigma, 0.0, 0.0};
   auto result = lj->calc_force_and_energy(r);
 
   EXPECT_NEAR(0.0, result.energy, 1e-10);
@@ -74,7 +74,7 @@ TEST_F(LennardJonesTest, AtEquilibrium) {
 
 TEST_F(LennardJonesTest, RepulsiveRegion) {
   // Test in the repulsive region (r < sigma)
-  Vec3<real> r(0.8 * sigma, 0.0, 0.0);
+  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
@@ -84,7 +84,7 @@ TEST_F(LennardJonesTest, RepulsiveRegion) {
 
 TEST_F(LennardJonesTest, AttractiveRegion) {
   // Test in the attractive region (sigma < r < r_min)
-  Vec3<real> r(1.5 * sigma, 0.0, 0.0);
+  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
@@ -95,15 +95,15 @@ TEST_F(LennardJonesTest, AttractiveRegion) {
 
 TEST_F(LennardJonesTest, ArbitraryDirection) {
   // Test with a vector in an arbitrary direction
-  Vec3<real> r(1.0, 1.0, 1.0);
+  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 rmag = std::sqrt(r.squared_norm2());
+  real r_mag = std::sqrt(r.squared_norm2());
 
   // Calculate expected force direction (should be along -r)
-  Vec3<real> normalized_r = r.scale(1.0 / rmag);
+  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;
@@ -120,11 +120,11 @@ TEST_F(LennardJonesTest, ParameterVariation) {
   // Test with different parameter values
   real new_sigma = 2.0;
   real new_epsilon = 0.5;
-  real new_rCutoff = 5.0;
+  real new_r_cutoff = 5.0;
 
-  LennardJones lj2(new_sigma, new_epsilon, new_rCutoff);
+  LennardJones lj2(new_sigma, new_epsilon, new_r_cutoff);
 
-  Vec3<real> r(2.0, 0.0, 0.0);
+  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);
 
@@ -136,7 +136,7 @@ TEST_F(LennardJonesTest, ParameterVariation) {
 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);
+  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
@@ -155,11 +155,11 @@ TEST_F(LennardJonesTest, ExactValueCheck) {
 
 TEST_F(LennardJonesTest, NearCutoff) {
   // Test behavior just inside and just outside the cutoff
-  real inside_cutoff = rCutoff - 0.01;
+  real inside_cutoff = r_cutoff - 0.01;
-  real outside_cutoff = rCutoff + 0.01;
+  real outside_cutoff = r_cutoff + 0.01;
 
-  Vec3<real> r_inside(inside_cutoff, 0.0, 0.0);
+  Vec3<real> r_inside = {inside_cutoff, 0.0, 0.0};
-  Vec3<real> r_outside(outside_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);
@@ -170,5 +170,5 @@ TEST_F(LennardJonesTest, NearCutoff) {
 
   // Outside should be zero
   EXPECT_EQ(0.0, result_outside.energy);
-  expectVec3Near(Vec3<real>(0.0, 0.0, 0.0), result_outside.force, 1e-10);
+  expect_vec3_near({0.0, 0.0, 0.0}, result_outside.force, 1e-10);
 }