diff --git a/CMakeLists.txt b/CMakeLists.txt
index 22cecac16f5694c0e7b7610e42a75b611180690b..0512280a636d8ac94b41cb60ed60096c5cc41ec7 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -96,27 +96,33 @@ endif()
 set(CMAKE_MACOSX_RPATH 1)	# Unsure if this works for CMAKE_BUIlD_RPATH, or just CMAKE_INSTALL_RPATH
 set(CMAKE_BUILD_RPATH "${CMAKE_CUDA_IMPLICIT_LINK_DIRECTORIES}")
 
+## Set up doctest
+# option(ENABLE_DOCTESTS "Include tests in the library. Setting this to OFF will remove all doctest related code.
+#                     Tests in tests/*.cpp will still be enabled." ON)
+# set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake/")
+# include(doctest)
+
+## Add subdirectories
+add_subdirectory(src)
+add_subdirectory(src/Tests)
+
 # set(CMAKE_VERBOSE_MAKEFILE True)
 add_executable("${PROJECT_NAME}" src/arbd.cpp
-  src/ParticlePatch.cpp
-  src/Integrator.cpp
-  src/Integrator/CPU.cpp
-  src/Integrator/CUDA.cu
-  src/SignalManager.cpp
-  src/GPUManager.cpp
-  src/SimManager.cpp
-  src/useful.cu
-  )
+src/SimManager.cpp
+src/useful.cu
+)
+
+target_link_libraries("${PROJECT_NAME}" PUBLIC "lib${PROJECT_NAME}")
 
 ## Add optional libraries
 if(USE_CUDA)
-  target_link_libraries("${PROJECT_NAME}" PRIVATE curand)
+  target_link_libraries("${PROJECT_NAME}" PUBLIC curand)
 endif()
 if(USE_NCCL)
-  target_link_libraries("${PROJECT_NAME}" PRIVATE nccl)
+  target_link_libraries("${PROJECT_NAME}" PUBLIC nccl)
 endif()
 if(USE_NVTX)
-  target_link_libraries("${PROJECT_NAME}" PRIVATE nvToolsExt)
+  target_link_libraries("${PROJECT_NAME}" PUBLIC nvToolsExt)
 endif()
 
 install(TARGETS "${PROJECT_NAME}")
diff --git a/src/GPUManager.h b/src/GPUManager.h
index 2086361712b1ddc837091702c16db9a94faf63b9..d4b8d91d654beb68a4041e8e6c063289cb9ecc36 100644
--- a/src/GPUManager.h
+++ b/src/GPUManager.h
@@ -6,7 +6,7 @@
 #include <cuda.h>
 #include <cuda_runtime_api.h>
 
-#include "useful.h"
+// #include "useful.h"
 
 #ifdef USE_NCCL
 #include <nccl.h>
diff --git a/src/Integrator.h b/src/Integrator.h
index 8e71eb99b1ac001617a71325c85bea71a4799978..7ccb0f941b79f8cb9ae3134aea2100702fdc15e1 100644
--- a/src/Integrator.h
+++ b/src/Integrator.h
@@ -1,26 +1,22 @@
+/*********************************************************************
+ * @file  Integrator.h
+ * 
+ * @brief Declaration of Integrator class with factory-like method
+ * GetIntegrator()
+ *
+ * Defines hashable Integrator::Conf struct that allows operator
+ * resuse. Common CPU/GPU kernels implemented in Integrator/kernels.h
+ * with various BD/MD integrators implemented on different backends in
+ * Integrator/CUDA.h and Integrator/CPU.h
+ *********************************************************************/
 #pragma once
 
 #include <cassert>
 #include <iostream>
 #include <map>
-#include "PatchOps.h"
+#include "PatchOp.h"
 
-#ifdef __CUDACC__
-    #define HOST __host__
-    #define DEVICE __device__
-#else
-    #define HOST
-    #define DEVICE
-#endif
-
-namespace IntegratorKernels {
-    HOST DEVICE  void __inline__ BDIntegrate() {
-	// std::cout << "Computes::BDIntegrate_inline" << std::endl;
-	printf("Integrator::BDIntegrate\n");
-    };
-}
-
-class Integrator : public BaseCompute {
+class Integrator : public BasePatchOp {
 public:
     virtual void compute(Patch* patch) = 0;
     int num_patches() const { return 1; };
@@ -45,5 +41,6 @@ protected:
 
 };
 
+#include "Integrator/kernels.h"
 #include "Integrator/CUDA.h"
 #include "Integrator/CPU.h"
diff --git a/src/Interaction.h b/src/Interaction.h
new file mode 100644
index 0000000000000000000000000000000000000000..cad0850688f8f38b27abf8a87bd784f8914ae311
--- /dev/null
+++ b/src/Interaction.h
@@ -0,0 +1,87 @@
+#pragma once
+
+#include <cassert>
+#include <iostream>
+#include <map>
+#include "PatchOp.h"
+
+class LocalInteraction : public BasePatchOp {
+public:
+    virtual void compute(Patch* patch) = 0;
+    int num_patches() const { return 1; };
+
+    // Following relates to lazy initialized factory method
+    struct Conf {
+	enum Object    {Particle, RigidBody };
+	enum DoF {Bond, Angle, Dihedral, Bonded, NeighborhoodPair};
+	enum Form {Harmonic, Tabulated, LennardJones, Coulomb, LJ};
+	enum Backend   { Default, CUDA, CPU };
+	
+	Object object_type;
+	DoF dof;
+	Form form;
+	Backend backend;
+
+	explicit operator int() const {return object_type*64 + dof*16 + form*4 + backend;};
+    };
+
+    static LocalInteraction* GetInteraction(Conf& conf);
+
+private:
+    size_t num_interactions;
+    
+protected:
+    static std::map<Conf, LocalInteraction*> _interactions;
+
+};
+
+// class LocalNonbondedInteraction : public BasePatchOp {
+// public:
+//     virtual void compute(Patch* patch) = 0;
+//     int num_patches() const { return 1; };
+
+//     // Following relates to lazy initialized factory method
+//     struct Conf {
+// 	enum Object    {Particle, RigidBody };
+// 	enum Electrostatics {None, Coulomb, DebyeHuckel}
+// 	enum Form      {Tabulated, LennardJones, Coulomb, LJ, DebyeHueckel};
+// 	enum Backend   { Default, CUDA, CPU };
+	
+// 	Object object_type;
+	
+// 	std::string tabulated_file = "";
+// 	Algorithm algorithm;
+// 	Backend backend;
+
+// 	explicit operator int() const {return object_type*16 + algorithm*4 + backend;};
+//     };
+
+//     static Integrator* GetIntegrator(Conf& conf);
+
+// private:
+//     size_t num_interactions;
+    
+// protected:
+//     static std::map<Conf, Integrator*> _integrators;
+
+// };
+
+// class LocalBonded : public LocalInteraction {
+
+// private:
+//     size_t bondlist;		// Encode bond, angle, dihedral
+    
+// };
+
+// class LocalBonds : public LocalInteraction {
+
+// private:
+//     // static const char* type = "Bond";
+//     size_t bondlist;
+    
+// };
+
+
+#include "Interaction/kernels.h"
+#include "Interaction/CUDA.h"
+#include "Interaction/CPU.h"
diff --git a/src/Interactions.h b/src/Interactions.h
deleted file mode 100644
index d58ae497e9eefdb54463b50689a0f28b12042859..0000000000000000000000000000000000000000
--- a/src/Interactions.h
+++ /dev/null
@@ -1,18 +0,0 @@
-#pragma once
-
-#include "useful.h"
-
-class Interactions {
-    // Object to store all kinds of info about the simulation system, but no particle data
-
-public:
-    size_t num_interactions;
-    
-};
-
-class BondInteractions : public Interactions {
-
-private:
-    static const char* type = "Bond";
-    size_t bondlist;
-};
diff --git a/src/ParticlePatch.cpp b/src/ParticlePatch.cpp
index 1977e1fbcfdffebddc54a55989e69e0f399261df..a4464acaff0c9d677521770f70eb16154012b532 100644
--- a/src/ParticlePatch.cpp
+++ b/src/ParticlePatch.cpp
@@ -1,11 +1,5 @@
 #include "ParticlePatch.h"
 
-// BasePatch::BasePatch(size_t num, short thread_id, short gpu_id) { ;};
-// BasePatch::BasePatch() {};
-// BasePatch::~BasePatch() {};
-
-// Patch::Patch(size_t num, short thread_id, short gpu_id) {};
-
 void Patch::compute() {
     for (auto& c_p: local_computes) {
 	c_p->compute(this);
diff --git a/src/ParticlePatch.cu b/src/ParticlePatch.cu
deleted file mode 100644
index 3628b208e26278975b34866e6cf9272c943524eb..0000000000000000000000000000000000000000
--- a/src/ParticlePatch.cu
+++ /dev/null
@@ -1,14 +0,0 @@
-#include "ParticlePatch.h"
-
-// BasePatch::BasePatch(size_t num, short thread_id, short gpu_id) { ;};
-// BasePatch::BasePatch() {};
-// BasePatch::~BasePatch() {};
-
-Patch::Patch(size_t num, short thread_id, short gpu_id) {};
-
-void Patch::compute() {
-    for (auto& c_p: local_computes) {
-	c_p->compute(this);
-    }
-};
-
diff --git a/src/ParticlePatch.h b/src/ParticlePatch.h
index 3f337868a2f28cdff0344a4b92e4dbd76ade774e..2efcda21b9598d6a17d8bee71c93afeb30e3bd2e 100644
--- a/src/ParticlePatch.h
+++ b/src/ParticlePatch.h
@@ -18,10 +18,10 @@
 #endif
 
 #include "SimSystem.h"
-#include "useful.h"
+// #include "useful.h"
+#include "Types.h"
 
-#include "PatchOps.h"
-//class BaseCompute;
+#include "PatchOp.h"
 
 class BasePatch {
 public:
@@ -56,10 +56,10 @@ public:
     // void addParticles(int n, int typ);
     // template<class T>
     // void add_compute(std::unique_ptr<T>&& p) {
-    // 	std::unique_ptr<BaseCompute> base_p = static_cast<std::unique_ptr<BaseCompute>>(p);
+    // 	std::unique_ptr<BasePatchOp> base_p = static_cast<std::unique_ptr<BasePatchOp>>(p);
     // 	local_computes.emplace_back(p);
     // };
-    void add_compute(std::unique_ptr<BaseCompute>&& p) {
+    void add_compute(std::unique_ptr<BasePatchOp>&& p) {
 	local_computes.emplace_back(std::move(p));
     };
 
@@ -67,8 +67,8 @@ public:
     
 private:
     // std::vector<PatchProxy> neighbors;    
-    std::vector<std::unique_ptr<BaseCompute>> local_computes; // Operations that will be performed on this patch each timestep
-    std::vector<std::unique_ptr<BaseCompute>> nonlocal_computes; // Operations that will be performed on this patch each timestep
+    std::vector<std::unique_ptr<BasePatchOp>> local_computes; // Operations that will be performed on this patch each timestep
+    std::vector<std::unique_ptr<BasePatchOp>> nonlocal_computes; // Operations that will be performed on this patch each timestep
     
     static int patch_idx;		// Unique ID across ranks
 
diff --git a/src/PatchOp.h b/src/PatchOp.h
new file mode 100644
index 0000000000000000000000000000000000000000..b1a4223dde39f3ca6c8acb4e6d4586b83b49c214
--- /dev/null
+++ b/src/PatchOp.h
@@ -0,0 +1,25 @@
+/*********************************************************************
+ * @file  PatchOp.h
+ * 
+ * @brief Declaration of BasePatchOp class. 
+ *
+ * Includes headers of derived classes for convenient access to
+ * factory methods. Virtual method `BasePatchOp::compute(Patch*
+ * patch)` called by `patch->compute()` by any Patch to which the
+ * PatchOp has been added.
+ *********************************************************************/
+#pragma once
+
+class Patch;
+
+/// Abstract base class that operates on Patch data.
+class BasePatchOp {
+public:
+    virtual void compute(Patch* patch) = 0;
+    virtual int num_patches() const = 0;
+private:
+    void* compute_data;
+};
+
+#include "Integrator.h"
+#include "Interaction.h"
diff --git a/src/PatchOps.h b/src/PatchOps.h
deleted file mode 100644
index 13115b01602576f6ad286f11ac8d3b99d67ead1b..0000000000000000000000000000000000000000
--- a/src/PatchOps.h
+++ /dev/null
@@ -1,14 +0,0 @@
-#pragma once
-
-class Patch;
-
-class BaseCompute {
-    // Low level class that operates on Patch data
-public:
-    virtual void compute(Patch* patch) = 0;
-    virtual int num_patches() const = 0;
-private:
-    void* compute_data;
-};
-
-#include "Integrator.h"
diff --git a/src/SimManager.cpp b/src/SimManager.cpp
index 8fee0dc43450a373c057abf40cd39113954d9e6e..defd93cbe810281df01ac503389b26fbee47a4f0 100644
--- a/src/SimManager.cpp
+++ b/src/SimManager.cpp
@@ -1,16 +1,14 @@
 #include "SimManager.h"
 #include <memory>
 
-// class LocalPairForce;
-// class NeighborPairForce;
-// class BDIntegrate;
-// #include "Computes.h"
-
 void SimManager::run() {
     std::cout << "running" << std::endl;
+
     // SimSystem sys = SimSystem();
     // Patch p(10,0,0,sys);
+
     Patch p(10,0,0);
+
     //ProxyPatch p2(10,0,0);
 
     // p.add_compute( std::make_unique<LocalPairForce>() );
@@ -18,11 +16,14 @@ void SimManager::run() {
 
 #ifdef USE_CUDA
     p.add_compute( std::make_unique<BDIntegrateCUDA>() );
+    p.add_compute( std::make_unique<LocalBondedCUDA>() );
 #else
     p.add_compute( std::make_unique<BDIntegrate>() );
+    p.add_compute( std::make_unique<LocalBonded>() );
 #endif
     
     for (size_t step = 0; step < 10; ++step) {
+	printf("Step\n");
 	p.compute();
 #ifdef USE_CUDA
 	cudaDeviceSynchronize();
diff --git a/src/SimManager.h b/src/SimManager.h
index f48f00176af5dddd04ee475bd6a9e91ec14a17ad..a379954a2834ac67b194725e0127399a37327302 100644
--- a/src/SimManager.h
+++ b/src/SimManager.h
@@ -1,7 +1,7 @@
 #pragma once
 #include <iostream>
 #include "ParticlePatch.h"
-#include "PatchOps.h"
+#include "PatchOp.h"
 
 class SimManager {
 public:
diff --git a/version.cmake b/version.cmake
deleted file mode 100644
index c366f35dd4c189dda3c9185b9687cb4b6d93c54d..0000000000000000000000000000000000000000
--- a/version.cmake
+++ /dev/null
@@ -1,43 +0,0 @@
-execute_process(COMMAND git log --pretty=format:'%h' -n 1
-                OUTPUT_VARIABLE GIT_REV
-		                ERROR_QUIET)
-
-# Check whether we got any revision (which isn't
-# always the case, e.g. when someone downloaded a zip
-# file from Github instead of a checkout)
-if ("${GIT_REV}" STREQUAL "")
-    set(GIT_REV "N/A")
-        set(GIT_DIFF "")
-	    set(GIT_TAG "N/A")
-	        set(GIT_BRANCH "N/A")
-		else()
-		    execute_process(
-		            COMMAND bash -c "git diff --quiet --exit-code || echo +"
-			            OUTPUT_VARIABLE GIT_DIFF)
-				        execute_process(
-					        COMMAND git describe --exact-match --tags
-						        OUTPUT_VARIABLE GIT_TAG ERROR_QUIET)
-							    execute_process(
-							            COMMAND git rev-parse --abbrev-ref HEAD
-								            OUTPUT_VARIABLE GIT_BRANCH)
-
-    string(STRIP "${GIT_REV}" GIT_REV)
-        string(SUBSTRING "${GIT_REV}" 1 7 GIT_REV)
-	    string(STRIP "${GIT_DIFF}" GIT_DIFF)
-	        string(STRIP "${GIT_TAG}" GIT_TAG)
-		    string(STRIP "${GIT_BRANCH}" GIT_BRANCH)
-		    endif()
-
-set(VERSION "const char* GIT_REV=\"${GIT_REV}${GIT_DIFF}\";
-const char* GIT_TAG=\"${GIT_TAG}\";
-const char* GIT_BRANCH=\"${GIT_BRANCH}\";")
-
-if(EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/version.cpp)
-    file(READ ${CMAKE_CURRENT_SOURCE_DIR}/version.cpp VERSION_)
-    else()
-        set(VERSION_ "")
-	endif()
-
-if (NOT "${VERSION}" STREQUAL "${VERSION_}")
-    file(WRITE ${CMAKE_CURRENT_SOURCE_DIR}/version.cpp "${VERSION}")
-endif()
\ No newline at end of file