diff --git a/src/main.cpp b/src/main.cpp index 1ff4753..a608bd7 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -20,11 +20,8 @@ int main() { System system(0.001); - std::vector atoms = generateFCCAtoms(5.26, 3); - system.addAtom(std::move(atoms)); - - atoms = generateFCCAtoms(5.26, 3); - offsetAtoms(atoms, 15, 15, 15); + std::vector atoms = generateFCCAtoms(5.26, 3, "Ar", 1); + initMaxwellBoltzmannVelocities(atoms, 300); system.addAtom(std::move(atoms)); system.printAtoms(); diff --git a/src/md/Atom.cpp b/src/md/Atom.cpp index 219a856..4c6e2f7 100644 --- a/src/md/Atom.cpp +++ b/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 atoms) { diff --git a/src/md/Atom.hpp b/src/md/Atom.hpp index f6aa6ad..b4dc214 100644 --- a/src/md/Atom.hpp +++ b/src/md/Atom.hpp @@ -11,6 +11,7 @@ #pragma once +#include #include /** @@ -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) { } }; diff --git a/src/md/AtomsGenerator.cpp b/src/md/AtomsGenerator.cpp index 07dd3d1..1d31a23 100644 --- a/src/md/AtomsGenerator.cpp +++ b/src/md/AtomsGenerator.cpp @@ -9,11 +9,18 @@ * */ +#include #include #include "AtomsGenerator.hpp" -std::vector generateFCCAtoms(double latticeConstant, std::size_t numberOfCells) { +// Private declarations +double gaussrand(); + +std::vector 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 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 &atoms, double x, double y, double z) { atom.y += y; atom.z += z; } +} + +void initMaxwellBoltzmannVelocities(std::vector &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; } \ No newline at end of file diff --git a/src/md/AtomsGenerator.hpp b/src/md/AtomsGenerator.hpp index ddc4ff1..92d21a5 100644 --- a/src/md/AtomsGenerator.hpp +++ b/src/md/AtomsGenerator.hpp @@ -22,7 +22,18 @@ * @param numberOfCells * @return std::vector */ -std::vector generateFCCAtoms(double latticeConstant, std::size_t numberOfCells); +std::vector 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 &atoms, double temperature); /** * Utilitary functions diff --git a/src/md/System.cpp b/src/md/System.cpp index 06c5c86..24059fb 100644 --- a/src/md/System.cpp +++ b/src/md/System.cpp @@ -9,14 +9,35 @@ * */ +#include #include #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) { diff --git a/src/md/System.hpp b/src/md/System.hpp index f03812e..442a41d 100644 --- a/src/md/System.hpp +++ b/src/md/System.hpp @@ -11,33 +11,58 @@ #pragma once +#include #include #include "Atom.hpp" +typedef std::function PotentialFunction; +typedef std::function ForceFunction; + class System { private: std::vector atoms; ///< Atoms in the system std::vector atomsOld; ///< Atoms in the system in the previous step + std::vector forceTable; ///< Force table + std::vector forceTableOld; ///< Force table in the previous step + std::vector 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(); /**