Add Space class

The Space class handles the simulation space

src/main.cpp

@@ -13,18 +13,35 @@
 
 #include "md/Atom.hpp"
 #include "md/AtomsGenerator.hpp"
+#include "md/Space.hpp"
 #include "md/System.hpp"
 
 int main() {
   std::cout << "kk eae men" << std::endl;
 
-  System system(0.001);
+  // TODO: get these from command line arguments
+  double lattice       = 5.26;
+  int    numberOfCells = 5;
+  double temperature   = 300;
 
-  std::vector<Atom> atoms = generateFCCAtoms(5.26, 3, "Ar", 1);
-  initMaxwellBoltzmannVelocities(atoms, 300);
-  system.addAtom(std::move(atoms));
+  Space space(3 * 1.5);
 
-  system.printAtoms();
+  std::vector<Atom> atoms = generateFCCAtoms(lattice, numberOfCells, "Ar", 1);
+  initMaxwellBoltzmannVelocities(atoms, temperature);
+
+  std::cout << "Adding atoms to space..." << std::endl;
+  space.addAtom(atoms);
+
+  std::cout << "Setting box..." << std::endl;
+  space.setBox(lattice * numberOfCells, lattice * numberOfCells, lattice * numberOfCells);
+
+  std::cout << "Preparing space..." << std::endl;
+  space.prepareSpace();
+
+  std::cout << "Building cells..." << std::endl;
+  space.buildCells();
+
+  space.printCells();
 
   return 0;
 }

src/md/CMakeLists.txt

@@ -6,4 +6,5 @@ add_library(md
   MolecularDynamics.cpp
   Atom.cpp
   AtomsGenerator.cpp
-  System.cpp)
+  System.cpp
+  Space.cpp)

src/md/Potentials.hpp

@@ -11,8 +11,21 @@
 
 #pragma once
 
+#include <cmath>
+
 class Potential {
 public:
   virtual double potential(double r) = 0;
   virtual double force(double r)     = 0;
 };
+
+class LennardJones : public Potential {
+public:
+  double potential(double r) override {
+    return 4 * (std::pow(r, -12) - std::pow(r, -6));
+  }
+
+  double force(double r) override {
+    return 24 * (2 * std::pow(r, -13) - std::pow(r, -7));
+  }
+};

src/md/Space.cpp

@@ -0,0 +1,119 @@
+/**
+ * @file Space.cpp
+ * @author jun <jun@firmwarejun.net>
+ * @brief Space definitions
+ * @version 0.1
+ * @date 2023-08-06
+ *
+ * Copyright (c) 2023 jun <https://git.firmwarejun.net/jun/MolecularDynamics2>
+ *
+ */
+
+#include <iostream>
+#include <stdexcept>
+
+#include "Space.hpp"
+
+Space::Space(double cutoff) {
+  if (cutoff <= 0.0)
+    throw std::invalid_argument("cutoff must be positive");
+
+  this->m_cutoff = cutoff;
+}
+
+void Space::addAtom(Atom atom) {
+  if (this->m_locked)
+    throw std::runtime_error("cannot add atom to locked space");
+
+  this->m_atoms.push_back(atom);
+}
+
+void Space::addAtom(std::vector<Atom> atoms) {
+  if (this->m_locked)
+    throw std::runtime_error("cannot add atoms to locked space");
+
+  this->m_atoms.insert(this->m_atoms.end(), atoms.begin(), atoms.end());
+}
+
+void Space::setBox(double x, double y, double z) {
+  if (this->m_locked)
+    throw std::runtime_error("cannot set box of locked space");
+
+  if (x <= 0.0 || y <= 0.0 || z <= 0.0)
+    throw std::invalid_argument("box dimensions must be positive");
+
+  this->m_box.x = x;
+  this->m_box.y = y;
+  this->m_box.z = z;
+}
+
+void Space::prepareSpace() {
+  if (this->m_locked)
+    throw std::runtime_error("cannot prepare locked space");
+
+  if (this->m_atoms.empty())
+    throw std::runtime_error("cannot prepare space with no atoms");
+
+  if (this->m_box.x <= 0.0 || this->m_box.y <= 0.0 || this->m_box.z <= 0.0)
+    throw std::runtime_error("cannot prepare space with invalid box");
+
+  if (this->m_cutoff <= 0.0)
+    throw std::runtime_error("cannot prepare space with invalid cutoff");
+
+  // Calculate the number of cells in each dimension
+  this->m_cellsDim.x = static_cast<int>(this->m_box.x / this->m_cutoff);
+  this->m_cellsDim.y = static_cast<int>(this->m_box.y / this->m_cutoff);
+  this->m_cellsDim.z = static_cast<int>(this->m_box.z / this->m_cutoff);
+
+  this->m_locked = true;
+}
+
+void Space::buildCells() {
+  if (!this->m_locked)
+    throw std::runtime_error("cannot build cells of unlocked space");
+
+  this->m_cells.clear();
+
+  for (int i = 0; i < this->m_cellsDim.x; i++)
+    for (int j = 0; j < this->m_cellsDim.y; j++)
+      for (int k = 0; k < this->m_cellsDim.z; k++) {
+        Cell cell;
+        cell.m_x   = i * this->m_cutoff;
+        cell.m_y   = j * this->m_cutoff;
+        cell.m_z   = k * this->m_cutoff;
+        cell.m_idx = i * this->m_cellsDim.y * this->m_cellsDim.z + j * this->m_cellsDim.z + k;
+        this->m_cells.push_back(cell);
+      }
+
+  std::cout << this->m_cells.size() << " cells built" << std::endl;
+
+  for (std::size_t a = 0; a < this->m_atoms.size(); a++) {
+    Atom &atom = this->m_atoms[a];
+
+    int i = static_cast<int>(atom.x / this->m_cutoff);
+    int j = static_cast<int>(atom.y / this->m_cutoff);
+    int k = static_cast<int>(atom.z / this->m_cutoff);
+
+    Cell &cell =
+        this->m_cells[i * this->m_cellsDim.y * this->m_cellsDim.z + j * this->m_cellsDim.z + k];
+
+    cell.x.push_back(atom.x);
+    cell.y.push_back(atom.y);
+    cell.z.push_back(atom.z);
+    cell.idx.push_back(a);
+  }
+}
+
+void Space::printCells() const {
+  for (const Cell &cell : this->m_cells)
+    cell.printCell();
+}
+
+void Cell::printCell() const {
+  std::cout << "Cell " << this->m_idx << " (" << this->m_x << ", " << this->m_y << ", " << this->m_z
+            << "):" << std::endl;
+
+  for (std::size_t i = 0; i < this->idx.size(); i++)
+    std::cout << "  " << this->idx[i] << ": (" << this->x[i] << ", " << this->y[i] << ", "
+              << this->z[i] << ")" << std::endl;
+}

src/md/Space.hpp

@@ -0,0 +1,111 @@
+/**
+ * @file Grid.hpp
+ * @author jun <jun@firmwarejun.net>
+ * @brief Grid declarations
+ * @version 0.1
+ * @date 2023-08-05
+ *
+ * Copyright (c) 2023 jun <https://git.firmwarejun.net/jun/MolecularDynamics2>
+ *
+ */
+
+#pragma once
+
+#include <vector>
+
+#include "Atom.hpp"
+
+class Space;
+
+/**
+ * @brief Cell structure
+ *
+ * The simulation space is partitioned into cells. Each cell contains a list of
+ * atoms that are inside the cell in a SoA format.
+ */
+class Cell {
+private:
+  friend Space;
+
+  double m_x, m_y, m_z; ///< The position of the cell
+  int    m_idx;         ///< The index of the cell
+
+public:
+  std::vector<double>      x, y, z; ///< Positions of atoms in the cell
+  std::vector<std::size_t> idx;     ///< Original indices of the atoms in AoS
+
+  /**
+   * @brief Print cell information along with the atoms in it
+   *
+   */
+  void printCell() const;
+};
+
+/**
+ * @brief The space of the simulation
+ *
+ * This class contains functions for manipulating the simulation space in various ways. It maintains
+ * a copy of the atoms in an Array of Structures (AoS).
+ */
+class Space {
+private:
+  /**
+   * @brief Box structure
+   *
+   * @tparam T
+   */
+  template <typename T> struct Box { T x, y, z; };
+
+  Box<double> m_box      = {}; ///< The box of the simulation
+  Box<int>    m_cellsDim = {}; ///< The number of cells in each dimension
+
+  double m_cutoff = 0;     ///< The cutoff for building the cells
+  bool   m_locked = false; ///< Whether the space is locked
+
+  std::vector<Atom> m_atoms; ///< The atoms in the simulation in AoS format
+  std::vector<Cell> m_cells; ///< The cells of the simulation
+
+public:
+  explicit Space(double cutoff = 0.0);
+
+  /**
+   * @brief Add an atom to the simulation space
+   *
+   * @param atom
+   */
+  void addAtom(Atom atom);
+
+  /**
+   * @brief Add atoms to the simulation space
+   *
+   * @param atoms
+   */
+  void addAtom(std::vector<Atom> atoms);
+
+  /**
+   * @brief Build the cells of the simulation based on the cutoff and simulation box
+   *
+   */
+  void buildCells();
+
+  /**
+   * @brief Set the size of the simulation box
+   *
+   * @param x
+   * @param y
+   * @param z
+   */
+  void setBox(double x, double y, double z);
+
+  /**
+   * @brief Locks the space and prepares it for simulation
+   *
+   */
+  void prepareSpace();
+
+  /**
+   * @brief Print the cells of the simulation
+   *
+   */
+  void printCells() const;
+};

src/md/System.cpp

@@ -14,13 +14,10 @@
 
 #include "System.hpp"
 
-System::System(double            timeDelta,
-               double            rcut,
-               PotentialFunction potentialFunction,
-               ForceFunction     forceFunction,
-               std::size_t       tableResolution)
-    : timeDelta(timeDelta), rcut(rcut), rcut2(rcut * rcut), potentialGenerator(potentialFunction),
-      forceGenerator(forceFunction), tableResolution(tableResolution) {
+System::System(double timeDelta, double rcut, Potential *potential, std::size_t tableResolution)
+    : timeDelta(timeDelta), rcut(rcut), rcut2(rcut * rcut), potential(potential),
+      tableResolution(tableResolution) {
+  this->initTables();
 }
 
 void System::initTables() {
@@ -29,12 +26,20 @@ void System::initTables() {
 
   for (std::size_t i = 0; i < this->tableResolution; i++) {
     double r2 = i * dr;
-    this->forceTable.push_back(this->forceGenerator(sqrt(r2)));
-    this->potentialTable.push_back(this->potentialGenerator(sqrt(r2)));
+    this->forceTable.push_back(this->potential->force(sqrt(r2)));
+    this->potentialTable.push_back(this->potential->potential(sqrt(r2)));
   }
 }
 
 void System::stepFirst() {
+  this->buildGrid();
+  this->evaulateForces();
+
+  // for (auto &atom : this->atoms) {
+  //   atom.vx += 0.5 * atom.fx * this->timeDelta;
+  //   atom.vy += 0.5 * atom.fy * this->timeDelta;
+  //   atom.vz += 0.5 * atom.fz * this->timeDelta;
+  // }
 }
 
 void System::evaulateForces() {
@@ -52,3 +57,12 @@ void System::printAtoms() {
   for (auto atom : this->atoms)
     printAtom(atom);
 }
+
+void System::setBoxSize(double lx, double ly, double lz) {
+  if (lx < 0 || ly < 0 || lz < 0)
+    throw std::runtime_error("Box dimensions cannot be negative");
+
+  this->lx = lx;
+  this->ly = ly;
+  this->lz = lz;
+}

src/md/System.hpp

@@ -15,9 +15,7 @@
 #include <vector>
 
 #include "Atom.hpp"
-
-typedef std::function<double(double)> PotentialFunction;
-typedef std::function<double(double)> ForceFunction;
+#include "Potentials.hpp"
 
 class System {
 private:
@@ -28,8 +26,7 @@ private:
   std::vector<double> forceTableOld;  ///< Force table in the previous step
   std::vector<double> potentialTable; ///< Potential table
 
-  ForceFunction     forceGenerator;     ///< Force generator
-  PotentialFunction potentialGenerator; ///< Potential generator
+  Potential *potential; ///< Potential function
 
   double rcut;  ///< Cutoff radius
   double rcut2; ///< Cutoff radius squared
@@ -38,6 +35,8 @@ private:
 
   double timeDelta = 0; ///< Time delta between steps
 
+  double lx, ly, lz; ///< Box size, if periodic boundary conditions are used
+
   /**
    * @brief Initialize the force and potential tables
    *
@@ -50,11 +49,17 @@ private:
    */
   void evaulateForces();
 
+  /**
+   * @brief Build the grid based on the cutoff radius.
+   *
+   * This function partitions the space into cells of size rcut x rcut x rcut.
+   */
+  void buildGrid();
+
 public:
   explicit System(double      timeDelta,
                   double      rcut,
-                  PotentialFunction potentialFunction,
-                  ForceFunction     forceFunction,
+                  Potential  *potential,
                   std::size_t tableResolution = 1000);
 
   void step();
@@ -84,4 +89,13 @@ public:
    *
    */
   void printAtoms();
+
+  /**
+   * @brief Set the simulation box size
+   *
+   * @param lx
+   * @param ly
+   * @param lz
+   */
+  void setBoxSize(double lx, double ly, double lz);
 };