From 6162b27a89b4f938c6a6848870e9b4852a70c352 Mon Sep 17 00:00:00 2001 From: Alex Selimov Date: Wed, 16 Apr 2025 14:01:40 -0400 Subject: [PATCH] Add pair potential and tests --- src/CMakeLists.txt | 13 ++- src/box.hpp | 22 ++++ src/pair_potentials.cpp | 33 ++++++ src/pair_potentials.hpp | 41 +++++++ src/particle.hpp | 18 +++ src/precision.hpp | 15 +++ src/simulation.hpp | 17 +++ src/test.cpp | 4 - src/test.h | 2 - tests/unit_tests/CMakeLists.txt | 5 +- tests/unit_tests/test_potential.cpp | 174 ++++++++++++++++++++++++++++ 11 files changed, 330 insertions(+), 14 deletions(-) create mode 100644 src/box.hpp create mode 100644 src/pair_potentials.cpp create mode 100644 src/pair_potentials.hpp create mode 100644 src/particle.hpp create mode 100644 src/precision.hpp create mode 100644 src/simulation.hpp delete mode 100644 src/test.cpp delete mode 100644 src/test.h create mode 100644 tests/unit_tests/test_potential.cpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 6c6cd8f..a0cc382 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,16 +1,17 @@ project(${CMAKE_PROJECT_NAME}_lib CUDA CXX) set(HEADER_FILES - ./test.h + particle.hpp + simulation.hpp + box.hpp + pair_potentials.hpp ) set(SOURCE_FILES - ./test.cpp + pair_potentials.cpp ) # The library contains header and source files. -add_library(${CMAKE_PROJECT_NAME}_lib +add_library(${CMAKE_PROJECT_NAME}_lib + ${HEADER_FILES} ${SOURCE_FILES} - ${HEADER_FILES} ) - -target_include_directories(${CMAKE_PROJECT_NAME}_lib PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) diff --git a/src/box.hpp b/src/box.hpp new file mode 100644 index 0000000..816b53e --- /dev/null +++ b/src/box.hpp @@ -0,0 +1,22 @@ +#ifndef BOX_H +#define BOX_H + +/** + * Struct representing the simulation box. + * Currently the simulation box is always assumed to be perfectly rectangular. + * This code does not support shearing the box. This functionality may be added + * in later. + */ +template struct Box { + T xlo; + T xhi; + T ylo; + T yhi; + T zlo; + T zhi; + bool x_is_periodic; + bool y_is_periodic; + bool z_is_periodic; +}; + +#endif diff --git a/src/pair_potentials.cpp b/src/pair_potentials.cpp new file mode 100644 index 0000000..fc5f6ba --- /dev/null +++ b/src/pair_potentials.cpp @@ -0,0 +1,33 @@ +#include "potentials.hpp" +#include + +/** + * Calculate the Lennard-Jones energy and force for the current particle pair + * described by displacement vector r + */ +ForceAndEnergy LennardJones::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_r5 = 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); + Vec3 force = r.scale(force_mag * inv_rmag); + + return {energy, force}; + + } else { + return ForceAndEnergy::zero(); + } +}; diff --git a/src/pair_potentials.hpp b/src/pair_potentials.hpp new file mode 100644 index 0000000..13d98c7 --- /dev/null +++ b/src/pair_potentials.hpp @@ -0,0 +1,41 @@ +#ifndef POTENTIALS_H +#define POTENTIALS_H + +#include "precision.hpp" +#include "vec3.h" + +/** + * Result struct for the Pair Potential + */ +struct ForceAndEnergy { + real energy; + Vec3 force; + + 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; + + /** + * Calculate the force and energy for a specific atom pair based on a + * displacement vector r. + */ + virtual ForceAndEnergy calc_force_and_energy(Vec3 r) = 0; +}; + +struct LennardJones : PairPotential { + real m_epsilon; + real m_sigma; + + ForceAndEnergy calc_force_and_energy(Vec3 r); +}; + +#endif diff --git a/src/particle.hpp b/src/particle.hpp new file mode 100644 index 0000000..84fd9b8 --- /dev/null +++ b/src/particle.hpp @@ -0,0 +1,18 @@ +#ifndef PARTICLE_H +#define PARTICLE_H + +#include "vec3.h" + +/** + * Class representing a single molecular dynamics particle. + * This class is only used on the host side of the code and is converted + * to the device arrays. + */ +template struct Particle { + Vec3 pos; + Vec3 vel; + Vec3 force; + T mass; +}; + +#endif diff --git a/src/precision.hpp b/src/precision.hpp new file mode 100644 index 0000000..c132c09 --- /dev/null +++ b/src/precision.hpp @@ -0,0 +1,15 @@ +#ifndef PRECISION_H +#define PRECISION_H + +#ifdef USE_FLOATS + +/* + * If macro USE_FLOATS is set then the default type will be floating point + * precision. Otherwise we use double precision by default + */ +typedef float real; +#else +typedef double real; +#endif + +#endif diff --git a/src/simulation.hpp b/src/simulation.hpp new file mode 100644 index 0000000..5b468b9 --- /dev/null +++ b/src/simulation.hpp @@ -0,0 +1,17 @@ +#ifndef SIMULATION_H +#define SIMULATION_H + +#include "box.hpp" +#include "particle.hpp" +#include + +template class Simulation { + // Simulation State variables + T timestep; + Box box; + + // Host Data + std::vector> particles; +}; + +#endif diff --git a/src/test.cpp b/src/test.cpp deleted file mode 100644 index eec7934..0000000 --- a/src/test.cpp +++ /dev/null @@ -1,4 +0,0 @@ -#include "test.h" -#include - -void test_hello() { std::cout << "Hello!"; } diff --git a/src/test.h b/src/test.h deleted file mode 100644 index b9c7ab1..0000000 --- a/src/test.h +++ /dev/null @@ -1,2 +0,0 @@ -#include -void test_hello(); diff --git a/tests/unit_tests/CMakeLists.txt b/tests/unit_tests/CMakeLists.txt index e7273a3..f13ed79 100644 --- a/tests/unit_tests/CMakeLists.txt +++ b/tests/unit_tests/CMakeLists.txt @@ -1,8 +1,9 @@ include_directories(${gtest_SOURCE_DIR}/include ${gtest_SOURCE_DIR}) add_executable(Unit_Tests_run - test_example.cpp + test_potential.cpp ) target_link_libraries(Unit_Tests_run gtest gtest_main) -target_link_libraries(Unit_Tests_run ${CMAKE_PROJECT_NAME}_lib) \ No newline at end of file +target_link_libraries(Unit_Tests_run ${CMAKE_PROJECT_NAME}_lib) +add_test(Name Tests COMMAND Unit_Tests_run) diff --git a/tests/unit_tests/test_potential.cpp b/tests/unit_tests/test_potential.cpp new file mode 100644 index 0000000..429e903 --- /dev/null +++ b/tests/unit_tests/test_potential.cpp @@ -0,0 +1,174 @@ +#include "potentials.hpp" +#include "precision.hpp" +#include "gtest/gtest.h" +#include + +class LennardJonesTest : public ::testing::Test { +protected: + void SetUp() override { + // Default parameters + sigma = 1.0; + epsilon = 1.0; + rCutoff = 2.5; + + // Create default LennardJones object + lj = new LennardJones(sigma, epsilon, rCutoff); + } + + void TearDown() override { delete lj; } + + real sigma; + real epsilon; + real rCutoff; + LennardJones *lj; + + // Helper function to compare Vec3 values with tolerance + void expectVec3Near(const Vec3 &expected, const Vec3 &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 r(0.0, 0.0, 0.0); + auto result = lj->calc_force_and_energy(r); + + EXPECT_EQ(0.0, result.energy); + expectVec3Near(Vec3(0.0, 0.0, 0.0), result.force, 1e-10); +} + +TEST_F(LennardJonesTest, BeyondCutoff) { + // Distance beyond cutoff should return zero force and energy + Vec3 r(3.0, 0.0, 0.0); // 3.0 > rCutoff (2.5) + auto result = lj->calc_force_and_energy(r); + + EXPECT_EQ(0.0, result.energy); + expectVec3Near(Vec3(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 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(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 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 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 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 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()); + + // Calculate expected force direction (should be along -r) + Vec3 normalized_r = r.scale(1.0 / rmag); + 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_rCutoff = 5.0; + + LennardJones lj2(new_sigma, new_epsilon, new_rCutoff); + + Vec3 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 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 = rCutoff - 0.01; + real outside_cutoff = rCutoff + 0.01; + + Vec3 r_inside(inside_cutoff, 0.0, 0.0); + Vec3 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); + expectVec3Near(Vec3(0.0, 0.0, 0.0), result_outside.force, 1e-10); +}