diff --git a/.gitignore b/.gitignore
index fd57df9..7d0b8c3 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,5 +1,8 @@
 # Builds
 build/
+Debug/
+Testing/
+compile_commands.json
 
 # Google Tests
 tests/lib/
@@ -7,3 +10,7 @@ tests/lib/
 # Jet Brains
 .idea/
 cmake-build-debug/
+
+# Cache dir
+.cache
+
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 36cd7d8..fb27a81 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,23 +1,63 @@
 cmake_minimum_required(VERSION 3.9)
-project(MyProject)
+set(NAME "cudaCAC")
+project(${NAME} LANGUAGES CUDA CXX)
+
+enable_testing()
+set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
+
+# Default settings 
+add_compile_options(-Wall -Wextra -Wpedantic)
+
+add_compile_options($<$<COMPILE_LANGUAGE:CUDA>:-Wno-pedantic>)
+
+# Add pedantic just for 
 
 set(CMAKE_CXX_STANDARD 17)
 
-set(SOURCE_FILES main.cpp)
-add_executable(${CMAKE_PROJECT_NAME}_run ${SOURCE_FILES})
+# Cuda Settings
+set(CMAKE_CUDA_ARCHITECTURES 61)
+set(CUDA_SEPARABLE_COMPILATION ON)
+# Cuda settings to get correct compile_commands.json
+set(CMAKE_CUDA_USE_RESPONSE_FILE_FOR_INCLUDES 0)
+set(CMAKE_CUDA_USE_RESPONSE_FILE_FOR_LIBRARIES 0)
+set(CMAKE_CUDA_USE_RESPONSE_FILE_FOR_OBJECTS 0)
 
+# Add Vec3  as a dependency
+include(FetchContent)
+FetchContent_Declare(Vec3
+    GIT_REPOSITORY https://www.alexselimov.com/git/aselimov/Vec3.git
+)
+
+FetchContent_GetProperties(Vec3)
+if(NOT Vec3_POPULATED)
+    FetchContent_MakeAvailable(Vec3)
+    include_directories(${Vec3_SOURCE_DIR})
+endif()
+
+include_directories(/usr/local/cuda-12.8/include)
 include_directories(src)
+include_directories(kernels)
 
 add_subdirectory(src)
+add_subdirectory(kernels)
 add_subdirectory(tests)
 
-target_link_libraries(${CMAKE_PROJECT_NAME}_run ${CMAKE_PROJECT_NAME}_lib)
+add_executable(${NAME} main.cpp)
+install(DIRECTORY src/ DESTINATION src/)
+
+target_link_libraries(
+    ${NAME} 
+    PRIVATE
+    ${NAME}_lib 
+    ${NAME}_cuda_lib 
+    ${CUDA_LIBRARIES}
+)
 
 # Doxygen Build
-option(BUILD_DOC "Build Documentation" ON)
+option(BUILD_DOC "Build Documentation" OFF)
 
 find_package(Doxygen)
-if(DOXYGEN_FOUND)
+if(DOXYGEN_FOUND AND BUILD_DOC)
     set(BUILD_DOC_DIR ${CMAKE_SOURCE_DIR}/build/docs)
     if(NOT EXISTS ${BUILD_DOC_DIR})
         file(MAKE_DIRECTORY ${BUILD_DOC_DIR})
@@ -33,6 +73,6 @@ if(DOXYGEN_FOUND)
             WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
             COMMENT "Generating API documentation with Doxygen"
             VERBATIM)
-else(DOXYGEN_FOUND)
+else(DOXYGEN_FOUND AND BUILD_DOC)
     message("Doxygen needs to be installed to generate the documentation.")
-endif(DOXYGEN_FOUND)
\ No newline at end of file
+endif(DOXYGEN_FOUND AND BUILD_DOC)
diff --git a/README.md b/README.md
index 816511c..9236f1b 100644
--- a/README.md
+++ b/README.md
@@ -6,75 +6,7 @@ following components:
 
 - Directory Structure
 - Make Build (CMake)
+- CUDA integration
 - Unit Test Framework (Google Test)
 - API Documentation (Doxygen)
 
-Feel free to fork this repository and tailor it to suit you.
-
-## Procedure
-1. Download Bash script to create new C++ projects 
-    ```bash
-    curl -O https://raw.githubusercontent.com/TimothyHelton/cpp_project_template/master/new_cpp_project.sh
-    chmod u+x new_cpp_project.sh
-    ```
-1. Create new C++ project
-    ```bash
-    ./new_cpp_project.sh NewProjectName
-    ```
-1. In the project top level **CMakeLists.txt**:
-    1. Line 2: Change the variable **MyProject** to the name of your project.
-        ```cmake
-        project(NewProject)
-        ```
-        - This variable will be used in a couple of different places.
-            - MyProject_run: will be the main executable name
-            - MyProject_lib: will be the project library name
-    1. Line 4: Set the version of C++ to use.  For example, let's set up the
-    NewProject to use C++ 11.
-        ```cmake
-        set(CMAKE_CXX_STANDARD 11)
-        ```
-1. Update project name and description in the `Doxyfile` located in the `docs`
-directory.
-    1. Update line `PROJECT_NAME`
-        1. This name will appear on each documentation page.
-    1. Update line `PROJECT_NUMBER`
-        1. This is the version number of your project.
-    1. Update line `PROJECT_BRIEF`
-        1. Any text entered here will also appear on each documentation page.
-        Try not to make this one too long.
-1. Reload the top CMake file.
-
-## CLION IDE Specific Instructions
-I started using an IDE from [JET Brains](https://www.jetbrains.com/) tailored
-for Python called [PyCharm](https://www.jetbrains.com/pycharm/) and thought
-it helped me write better code.
-I'd been wanting to learn C++ and decided to give JET Brains C/C++ IDE called
-[CLion](https://www.jetbrains.com/clion/) a try.
-The code completion, interactive suggestions, debugger, introspection tools,
-and built-in test execution are very handy.
-There are a couple extra details to set when using this IDE.
-
-1. The IDE allows you to mark directories with their desired purpose.
-To mark a directory right click on the directory name in the `Project` window
-and select `Mark Directory as` from the drop-down menu.
-    1. Mark the `src` directory as `Project Sources and Headers`
-    1. Mark the `tests/lib/googletest` directory as  `Library Files`
-1. Setup the `Run/Debug Configuration` by selecting `Edit Configurations...`
-from the pull-down menu from the run button (green triangle) in the upper right
-corner.
-    1. Update Doxygen Build to execute the unit test suite.
-        1. Select Doxygen from the Application menu on the left.
-        1. Choose the **executable** for Doxygen to be `Unit_Tests_run`.
-    1. Create a `Google Test` configuration
-        1. In the upper left corner select the plus symbol.
-        1. Chose `Google Test` from the drop-down menu.
-        1. Set **Name** to `Unit Tests`.
-        1. Set **Target** to `Unit_Tests_run`.
-
-## Wrap Up
-That should be all it takes to start writing code.
-If you find any issues or bugs with this repository please file an issue on
-[GitHub](https://github.com/TimothyHelton/cpp_project_template/issues).
-
-Hope you find this template useful and enjoy learning C++!
\ No newline at end of file
diff --git a/create_project.sh b/create_project.sh
deleted file mode 100644
index b92cf96..0000000
--- a/create_project.sh
+++ /dev/null
@@ -1,33 +0,0 @@
-#!/usr/bin/env bash
-
-# Exit if name argument is not given
-if [ -z "$*" ]; then
-    echo "A project name argument must be provided."
-    exit 0
-fi
-
-NAME=$1
-
-
-################################################################################
-
-
-# Clone template repository
-git clone https://github.com/TimothyHelton/cpp_project_template
-
-# Create bare repository
-git --bare init ${NAME}
-
-# Push template master branch to bare repository
-cd cpp_project_template
-git push ../${NAME} +master:master
-
-# Convert bare repository into a normal repository
-cd ../${NAME}
-mkdir .git
-mv * .git
-git config --local --bool core.bare false
-git reset --hard
-
-# Clean Up
-rm -rf ../cpp_project_template ../create_project.sh
\ No newline at end of file
diff --git a/kernels/CMakeLists.txt b/kernels/CMakeLists.txt
new file mode 100644
index 0000000..baa8a60
--- /dev/null
+++ b/kernels/CMakeLists.txt
@@ -0,0 +1,14 @@
+project(${NAME}_cuda_lib CUDA CXX)
+
+set(HEADER_FILES
+    pair_potentials.cuh
+)
+set(SOURCE_FILES
+)
+
+# The library contains header and source files.
+add_library(${NAME}_cuda_lib INTERFACE
+    ${SOURCE_FILES}
+    ${HEADER_FILES}
+)
+
diff --git a/kernels/pair_potentials.cuh b/kernels/pair_potentials.cuh
new file mode 100644
index 0000000..052a079
--- /dev/null
+++ b/kernels/pair_potentials.cuh
@@ -0,0 +1,91 @@
+#ifndef POTENTIALS_H
+#define POTENTIALS_H
+
+#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;
+
+  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 ~LennardJones(){};
+};
+
+PairPotential::~PairPotential() {};
+#endif
diff --git a/main.cpp b/main.cpp
index d269c89..1fcf9ba 100644
--- a/main.cpp
+++ b/main.cpp
@@ -1,3 +1,9 @@
+#include "particle.hpp"
+#include "vec3.h"
+#include <iostream>
+
 int main() {
+  Particle test = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, 10};
+  std::cout << test.pos.x << " " << test.pos.y << " " << test.pos.z;
   return 0;
-}
\ No newline at end of file
+}
diff --git a/new_cpp_project.sh b/new_cpp_project.sh
deleted file mode 100644
index 31924f5..0000000
--- a/new_cpp_project.sh
+++ /dev/null
@@ -1,20 +0,0 @@
-#!/usr/bin/env bash
-
-# Exit if name argument is not given
-if [ -z "$*" ]; then
-    echo "A project name argument must be provided."
-    exit 0
-fi
-
-NAME=$1
-
-
-################################################################################
-
-
-# Download latest version of the build file
-curl -O https://raw.githubusercontent.com/TimothyHelton/cpp_project_template/master/create_project.sh
-chmod u+x create_project.sh
-
-# Create Project
-./create_project.sh ${NAME}
\ No newline at end of file
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 94717aa..48c2307 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,17 +1,16 @@
-project(${CMAKE_PROJECT_NAME}_lib)
+project(${NAME}_lib CUDA CXX)
 
 set(HEADER_FILES
+    particle.hpp
+    simulation.hpp
+    box.hpp
+
 )
 set(SOURCE_FILES
 )
 
-if (EXISTS ${SOURCE_FILES})
-    # The library contains header and source files.
-    add_library(${CMAKE_PROJECT_NAME}_lib STATIC
-        ${SOURCE_FILES}
-        ${HEADER_FILES}
-    )
-else()
-    # The library only contains header files.
-    add_library(${CMAKE_PROJECT_NAME}_lib INTERFACE)
-endif()
\ No newline at end of file
+# The library contains header and source files.
+add_library(${NAME}_lib INTERFACE
+    ${HEADER_FILES} 
+    ${SOURCE_FILES}
+)
diff --git a/src/box.hpp b/src/box.hpp
new file mode 100644
index 0000000..b588c49
--- /dev/null
+++ b/src/box.hpp
@@ -0,0 +1,23 @@
+#ifndef BOX_H
+#define BOX_H
+
+#include "precision.hpp"
+/**
+ * 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.
+ */
+struct Box {
+  real xlo;
+  real xhi;
+  real ylo;
+  real yhi;
+  real zlo;
+  real zhi;
+  bool x_is_periodic;
+  bool y_is_periodic;
+  bool z_is_periodic;
+};
+
+#endif
diff --git a/src/particle.hpp b/src/particle.hpp
new file mode 100644
index 0000000..b9e3464
--- /dev/null
+++ b/src/particle.hpp
@@ -0,0 +1,19 @@
+#ifndef PARTICLE_H
+#define PARTICLE_H
+
+#include "precision.hpp"
+#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.
+ */
+struct Particle {
+  Vec3<real> pos;
+  Vec3<real> vel;
+  Vec3<real> force;
+  real 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..0d69d34
--- /dev/null
+++ b/src/simulation.hpp
@@ -0,0 +1,18 @@
+#ifndef SIMULATION_H
+#define SIMULATION_H
+
+#include "box.hpp"
+#include "particle.hpp"
+#include "precision.hpp"
+#include <vector>
+
+class Simulation {
+  // Simulation State variables
+  real timestep;
+  Box box;
+
+  // Host Data
+  std::vector<Particle> particles;
+};
+
+#endif
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 85a8157..7f994a6 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -10,4 +10,5 @@ if(NOT EXISTS ${GOOGLETEST_DIR})
 endif()
 
 add_subdirectory(lib/googletest)
-add_subdirectory(unit_tests)
\ No newline at end of file
+add_subdirectory(unit_tests)
+add_subdirectory(cuda_unit_tests)
diff --git a/tests/cuda_unit_tests/CMakeLists.txt b/tests/cuda_unit_tests/CMakeLists.txt
new file mode 100644
index 0000000..27490a0
--- /dev/null
+++ b/tests/cuda_unit_tests/CMakeLists.txt
@@ -0,0 +1,9 @@
+include_directories(${gtest_SOURCE_DIR}/include ${gtest_SOURCE_DIR})
+
+add_executable(${NAME}_cuda_tests
+    test_potential.cu
+)
+
+target_link_libraries(${NAME}_cuda_tests gtest gtest_main)
+target_link_libraries(${NAME}_cuda_tests ${CMAKE_PROJECT_NAME}_cuda_lib)
+add_test(NAME ${NAME}CudaTests COMMAND ${CMAKE_BINARY_DIR}/tests/unit_tests/${NAME}_tests)
diff --git a/tests/cuda_unit_tests/test_potential.cu b/tests/cuda_unit_tests/test_potential.cu
new file mode 100644
index 0000000..cab3216
--- /dev/null
+++ b/tests/cuda_unit_tests/test_potential.cu
@@ -0,0 +1,316 @@
+#include "pair_potentials.cuh"
+#include "precision.hpp"
+#include "gtest/gtest.h"
+#include <cmath>
+#include <cuda_runtime.h>
+
+// Structure to hold test results from device
+struct TestResults {
+  bool zero_distance_pass;
+  bool beyond_cutoff_pass;
+  bool at_minimum_pass;
+  bool at_equilibrium_pass;
+  bool repulsive_region_pass;
+  bool attractive_region_pass;
+  bool arbitrary_direction_pass;
+  bool parameter_variation_pass;
+  bool exact_value_check_pass;
+  bool near_cutoff_pass;
+
+  // Additional result data for exact checks
+  real energy_values[10];
+  Vec3<real> force_values[10];
+};
+
+// Check if two Vec3 values are close within tolerance
+__device__ bool vec3_near(const Vec3<real> &a, const Vec3<real> &b,
+                          real tolerance) {
+  return (fabs(a.x - b.x) < tolerance) && (fabs(a.y - b.y) < tolerance) &&
+         (fabs(a.z - b.z) < tolerance);
+}
+
+// Device kernel to run all tests
+__global__ void lennard_jones_test_kernel(TestResults *results) {
+  // Default parameters
+  real sigma = 1.0;
+  real epsilon = 1.0;
+  real r_cutoff = 2.5;
+  real tolerance = 1e-10;
+
+  // Create LennardJones object on device
+  LennardJones lj(sigma, epsilon, r_cutoff);
+
+  // 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;
+    results->zero_distance_pass =
+        (result.energy == 0.0) &&
+        vec3_near(Vec3<real>{0.0, 0.0, 0.0}, result.force, 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;
+    results->beyond_cutoff_pass =
+        (result.energy == 0.0) &&
+        vec3_near(Vec3<real>{0.0, 0.0, 0.0}, result.force, 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;
+    results->at_minimum_pass =
+        (fabs(result.energy + epsilon) < tolerance) &&
+        vec3_near(Vec3<real>{0.0, 0.0, 0.0}, result.force, 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);
+  }
+
+  // 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);
+  }
+
+  // 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);
+  }
+
+  // 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;
+
+    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;
+
+    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);
+  }
+
+  // Parameter Variation Test
+  {
+    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->energy_values[7] = result2.energy;
+    results->force_values[7] = result2.force;
+
+    results->parameter_variation_pass = (result1.energy != result2.energy) &&
+                                        (result1.force.x != result2.force.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);
+
+    results->energy_values[8] = result.energy;
+    results->force_values[8] = result.force;
+
+    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);
+  }
+
+  // Near Cutoff Test
+  {
+    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);
+
+    results->energy_values[9] = result_inside.energy;
+    results->force_values[9] = result_inside.force;
+
+    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);
+  }
+}
+
+// Helper class for CUDA error checking
+class CudaErrorCheck {
+public:
+  static void checkAndThrow(cudaError_t err, const char *msg) {
+    if (err != cudaSuccess) {
+      std::string error_message =
+          std::string(msg) + ": " + cudaGetErrorString(err);
+      throw std::runtime_error(error_message);
+    }
+  }
+};
+
+// Google Test wrapper that runs the device tests
+class LennardJonesCudaTest : public ::testing::Test {
+protected:
+  void SetUp() override {
+    // Allocate device memory for results
+    CudaErrorCheck::checkAndThrow(
+        cudaMalloc(&d_results, sizeof(TestResults)),
+        "Failed to allocate device memory for test results");
+  }
+
+  void TearDown() override {
+    if (d_results) {
+      cudaFree(d_results);
+      d_results = nullptr;
+    }
+  }
+
+  // Helper function to run tests on device and get results
+  TestResults runDeviceTests() {
+    TestResults h_results;
+
+    // Clear device memory
+    CudaErrorCheck::checkAndThrow(cudaMemset(d_results, 0, sizeof(TestResults)),
+                                  "Failed to clear device memory");
+
+    // Run kernel with a single thread
+    lennard_jones_test_kernel<<<1, 1>>>(d_results);
+
+    // Check for kernel launch errors
+    CudaErrorCheck::checkAndThrow(cudaGetLastError(), "Kernel launch failed");
+
+    // Wait for kernel to complete
+    CudaErrorCheck::checkAndThrow(cudaDeviceSynchronize(),
+                                  "Kernel execution failed");
+
+    // Copy results back to host
+    CudaErrorCheck::checkAndThrow(cudaMemcpy(&h_results, d_results,
+                                             sizeof(TestResults),
+                                             cudaMemcpyDeviceToHost),
+                                  "Failed to copy results from device");
+
+    return h_results;
+  }
+
+  TestResults *d_results = nullptr;
+};
+
+// Define the actual test cases
+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
+      << ")";
+}
+
+TEST_F(LennardJonesCudaTest, DeviceBeyondCutoff) {
+  auto results = runDeviceTests();
+  EXPECT_TRUE(results.beyond_cutoff_pass)
+      << "Beyond cutoff test failed on device. Energy: "
+      << results.energy_values[1];
+}
+
+TEST_F(LennardJonesCudaTest, DeviceAtMinimum) {
+  auto results = runDeviceTests();
+  EXPECT_TRUE(results.at_minimum_pass)
+      << "At minimum test failed on device. Energy: "
+      << results.energy_values[2];
+}
+
+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;
+}
+
+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;
+}
+
+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;
+}
+
+TEST_F(LennardJonesCudaTest, DeviceArbitraryDirection) {
+  auto results = runDeviceTests();
+  EXPECT_TRUE(results.arbitrary_direction_pass)
+      << "Arbitrary direction test failed on device.";
+}
+
+TEST_F(LennardJonesCudaTest, DeviceParameterVariation) {
+  auto results = runDeviceTests();
+  EXPECT_TRUE(results.parameter_variation_pass)
+      << "Parameter variation test failed on device.";
+}
+
+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;
+}
+
+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];
+}
diff --git a/tests/unit_tests/CMakeLists.txt b/tests/unit_tests/CMakeLists.txt
index e7273a3..c396ab7 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
+add_executable(${NAME}_tests
+    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(${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)
diff --git a/tests/unit_tests/test_potential.cpp b/tests/unit_tests/test_potential.cpp
new file mode 100644
index 0000000..a94f022
--- /dev/null
+++ b/tests/unit_tests/test_potential.cpp
@@ -0,0 +1,174 @@
+#include "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);
+}