diff --git a/src/main.cpp b/src/main.cpp index a608bd7..e81df38 100644 --- a/src/main.cpp +++ b/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 atoms = generateFCCAtoms(5.26, 3, "Ar", 1); - initMaxwellBoltzmannVelocities(atoms, 300); - system.addAtom(std::move(atoms)); + Space space(3 * 1.5); - system.printAtoms(); + std::vector 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; } \ No newline at end of file diff --git a/src/md/CMakeLists.txt b/src/md/CMakeLists.txt index 1f8220e..4190bdf 100644 --- a/src/md/CMakeLists.txt +++ b/src/md/CMakeLists.txt @@ -6,4 +6,5 @@ add_library(md MolecularDynamics.cpp Atom.cpp AtomsGenerator.cpp - System.cpp) \ No newline at end of file + System.cpp + Space.cpp) \ No newline at end of file diff --git a/src/md/Potentials.hpp b/src/md/Potentials.hpp index 9141fdf..564a8d1 100644 --- a/src/md/Potentials.hpp +++ b/src/md/Potentials.hpp @@ -11,8 +11,21 @@ #pragma once +#include + 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)); + } }; \ No newline at end of file diff --git a/src/md/Space.cpp b/src/md/Space.cpp new file mode 100644 index 0000000..e51ed10 --- /dev/null +++ b/src/md/Space.cpp @@ -0,0 +1,119 @@ +/** + * @file Space.cpp + * @author jun + * @brief Space definitions + * @version 0.1 + * @date 2023-08-06 + * + * Copyright (c) 2023 jun + * + */ + +#include +#include + +#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 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(this->m_box.x / this->m_cutoff); + this->m_cellsDim.y = static_cast(this->m_box.y / this->m_cutoff); + this->m_cellsDim.z = static_cast(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(atom.x / this->m_cutoff); + int j = static_cast(atom.y / this->m_cutoff); + int k = static_cast(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; +} \ No newline at end of file diff --git a/src/md/Space.hpp b/src/md/Space.hpp new file mode 100644 index 0000000..39eea13 --- /dev/null +++ b/src/md/Space.hpp @@ -0,0 +1,111 @@ +/** + * @file Grid.hpp + * @author jun + * @brief Grid declarations + * @version 0.1 + * @date 2023-08-05 + * + * Copyright (c) 2023 jun + * + */ + +#pragma once + +#include + +#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 x, y, z; ///< Positions of atoms in the cell + std::vector 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 struct Box { T x, y, z; }; + + Box m_box = {}; ///< The box of the simulation + Box 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 m_atoms; ///< The atoms in the simulation in AoS format + std::vector 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 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; +}; \ No newline at end of file diff --git a/src/md/System.cpp b/src/md/System.cpp index 24059fb..3736377 100644 --- a/src/md/System.cpp +++ b/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() { @@ -51,4 +56,13 @@ void System::addAtom(std::vector atoms) { 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; } \ No newline at end of file diff --git a/src/md/System.hpp b/src/md/System.hpp index 442a41d..ece8904 100644 --- a/src/md/System.hpp +++ b/src/md/System.hpp @@ -15,9 +15,7 @@ #include #include "Atom.hpp" - -typedef std::function PotentialFunction; -typedef std::function ForceFunction; +#include "Potentials.hpp" class System { private: @@ -28,8 +26,7 @@ private: std::vector forceTableOld; ///< Force table in the previous step std::vector 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,12 +49,18 @@ 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, - std::size_t tableResolution = 1000); + explicit System(double timeDelta, + double rcut, + 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); }; \ No newline at end of file