Add Maxwell-Boltzmann method for velocities

src/md/Atom.cpp

@@ -14,7 +14,9 @@
 #include "Atom.hpp"
 
 void printAtom(Atom atom) {
-  std::cout << "Atom(" << atom.x << ", " << atom.y << ", " << atom.z << ")" << std::endl;
+  // Print an atom in XYZ format with extra information
+  std::cout << atom.name << " " << atom.x << " " << atom.y << " " << atom.z << " " << atom.vx << " "
+            << atom.vy << " " << atom.vz << " " << atom.m << std::endl;
 }
 
 void printAtoms(std::vector<Atom> atoms) {

src/md/Atom.hpp

@@ -11,6 +11,7 @@
 
 #pragma once
 
+#include <string>
 #include <vector>
 
 /**
@@ -19,11 +20,13 @@
  */
 class Atom {
 public:
-  double x, y, z;    ///< Positions
-  double vx, vy, vz; ///< Velocities
-  double m;          ///< Mass
+  double      x, y, z;    ///< Positions
+  double      vx, vy, vz; ///< Velocities
+  double      m    = 1;   ///< Mass
+  std::string name = "";  ///< Name
 
-  Atom(double x, double y, double z) : x(x), y(y), z(z), vx(0), vy(0), vz(0), m(1) {
+  Atom(double x, double y, double z, std::string name, double mass)
+      : x(x), y(y), z(z), vx(0), vy(0), vz(0), m(mass), name(name) {
   }
 };
 

src/md/AtomsGenerator.cpp

@@ -9,11 +9,18 @@
  *
  */
 
+#include <cmath>
 #include <stdexcept>
 
 #include "AtomsGenerator.hpp"
 
-std::vector<Atom> generateFCCAtoms(double latticeConstant, std::size_t numberOfCells) {
+// Private declarations
+double gaussrand();
+
+std::vector<Atom> generateFCCAtoms(double      latticeConstant,
+                                   std::size_t numberOfCells,
+                                   std::string element,
+                                   double      mass) {
   // Check for invalid lattice constants
   if (latticeConstant <= 0)
     throw std::invalid_argument("Lattice constant must be positive");
@@ -24,13 +31,23 @@ std::vector<Atom> generateFCCAtoms(double latticeConstant, std::size_t numberOfC
   for (std::size_t i = 0; i < numberOfCells; i++)
     for (std::size_t j = 0; j < numberOfCells; j++)
       for (std::size_t k = 0; k < numberOfCells; k++) {
-        atoms.push_back(Atom(i * latticeConstant, j * latticeConstant, k * latticeConstant));
         atoms.push_back(
-            Atom((i + 0.5) * latticeConstant, (j + 0.5) * latticeConstant, k * latticeConstant));
-        atoms.push_back(
-            Atom((i + 0.5) * latticeConstant, j * latticeConstant, (k + 0.5) * latticeConstant));
-        atoms.push_back(
-            Atom(i * latticeConstant, (j + 0.5) * latticeConstant, (k + 0.5) * latticeConstant));
+            Atom(i * latticeConstant, j * latticeConstant, k * latticeConstant, element, mass));
+        atoms.push_back(Atom((i + 0.5) * latticeConstant,
+                             (j + 0.5) * latticeConstant,
+                             k * latticeConstant,
+                             element,
+                             mass));
+        atoms.push_back(Atom((i + 0.5) * latticeConstant,
+                             j * latticeConstant,
+                             (k + 0.5) * latticeConstant,
+                             element,
+                             mass));
+        atoms.push_back(Atom(i * latticeConstant,
+                             (j + 0.5) * latticeConstant,
+                             (k + 0.5) * latticeConstant,
+                             element,
+                             mass));
       }
 
   // Use std::move to move the vector to the caller
@@ -43,4 +60,52 @@ void offsetAtoms(std::vector<Atom> &atoms, double x, double y, double z) {
     atom.y += y;
     atom.z += z;
   }
+}
+
+void initMaxwellBoltzmannVelocities(std::vector<Atom> &atoms, double temperature) {
+  // Check for invalid temperature
+  if (temperature <= 0)
+    throw std::invalid_argument("Temperature must be positive");
+
+  double vx_sum = 0, vy_sum = 0, vz_sum = 0;
+
+  // Initialize random velocities according to the Maxwell-Boltzmann distribution
+  for (auto &atom : atoms) {
+    atom.vx = gaussrand() * sqrt(temperature / atom.m);
+    atom.vy = gaussrand() * sqrt(temperature / atom.m);
+    atom.vz = gaussrand() * sqrt(temperature / atom.m);
+    vx_sum += atom.vx;
+    vy_sum += atom.vy;
+    vz_sum += atom.vz;
+  }
+
+  // Calculate the average velocity
+  vx_sum /= atoms.size();
+  vy_sum /= atoms.size();
+  vz_sum /= atoms.size();
+
+  // Set the average velocity to zero
+  for (auto &atom : atoms) {
+    atom.vx -= vx_sum;
+    atom.vy -= vy_sum;
+    atom.vz -= vz_sum;
+  }
+}
+
+double gaussrand() {
+  static double U, V;
+  static int    phase = 0;
+  double        Z;
+
+  if (phase == 0) {
+    double u = (rand() + 1.0) / (RAND_MAX + 2.0);
+    double v = rand() / (RAND_MAX + 1.0);
+    Z        = sqrt(-2.0 * log(u)) * sin(2.0 * M_PI * v);
+  } else {
+    Z = sqrt(-2.0 * log(U)) * cos(2.0 * M_PI * V);
+  }
+
+  phase = 1 - phase;
+
+  return Z;
 }

src/md/AtomsGenerator.hpp

@@ -22,7 +22,18 @@
  * @param numberOfCells
  * @return std::vector<Atom>
  */
-std::vector<Atom> generateFCCAtoms(double latticeConstant, std::size_t numberOfCells);
+std::vector<Atom> generateFCCAtoms(double      latticeConstant,
+                                   std::size_t numberOfCells,
+                                   std::string element,
+                                   double      mass);
+
+/**
+ * @brief Initialize random velocities according to the Maxwell-Boltzmann distribution
+ *
+ * @param atoms
+ * @param temperature
+ */
+void initMaxwellBoltzmannVelocities(std::vector<Atom> &atoms, double temperature);
 
 /**
  * Utilitary functions

src/md/System.cpp

@@ -9,14 +9,35 @@
  *
  */
 
+#include <cmath>
 #include <iostream>
 
 #include "System.hpp"
 
-System::System(double timeDelta) : timeDelta(timeDelta) {
+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) {
 }
 
-void System::initForceTable() {
+void System::initTables() {
+  // Generate the force table based on squared distance
+  double dr = this->rcut2 / this->tableResolution;
+
+  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)));
+  }
+}
+
+void System::stepFirst() {
+}
+
+void System::evaulateForces() {
 }
 
 void System::addAtom(Atom atom) {

src/md/System.hpp

@@ -11,33 +11,58 @@
 
 #pragma once
 
+#include <functional>
 #include <vector>
 
 #include "Atom.hpp"
 
+typedef std::function<double(double)> PotentialFunction;
+typedef std::function<double(double)> ForceFunction;
+
 class System {
 private:
   std::vector<Atom> atoms;    ///< Atoms in the system
   std::vector<Atom> atomsOld; ///< Atoms in the system in the previous step
 
+  std::vector<double> forceTable;     ///< Force table
+  std::vector<double> forceTableOld;  ///< Force table in the previous step
+  std::vector<double> potentialTable; ///< Potential table
+
+  ForceFunction     forceGenerator;     ///< Force generator
+  PotentialFunction potentialGenerator; ///< Potential generator
+
+  double rcut;  ///< Cutoff radius
+  double rcut2; ///< Cutoff radius squared
+
+  std::size_t tableResolution = 1000; ///< Resolution of the force table
+
   double timeDelta = 0; ///< Time delta between steps
 
   /**
-   * @brief Initialize the force table
+   * @brief Initialize the force and potential tables
    *
    */
-  void initForceTable();
+  void initTables();
 
   /**
-   * @brief Initialize the potential table
+   * @brief Evaluate interaction forces between atoms
    *
    */
-  void initPotentialTable();
+  void evaulateForces();
 
 public:
-  explicit System(double timeDelta);
+  explicit System(double            timeDelta,
+                  double            rcut,
+                  PotentialFunction potentialFunction,
+                  ForceFunction     forceFunction,
+                  std::size_t       tableResolution = 1000);
 
   void step();
+
+  /**
+   * @brief Step the system to the first step using euler method
+   *
+   */
   void stepFirst();
 
   /**