From 1e4f6f2894d86bd7b0d16143abad14733e75288a Mon Sep 17 00:00:00 2001 From: Cory Date: Tue, 31 Oct 2023 10:52:00 +0100 Subject: [PATCH] First draft of Ising model --- include/IsingModel.hpp | 103 +++++++++++++++++++++++++++++++++++ src/IsingModel.cpp | 121 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 224 insertions(+) create mode 100644 include/IsingModel.hpp create mode 100644 src/IsingModel.cpp diff --git a/include/IsingModel.hpp b/include/IsingModel.hpp new file mode 100644 index 0000000..74b8323 --- /dev/null +++ b/include/IsingModel.hpp @@ -0,0 +1,103 @@ +/** @file IsingModel.hpp + * + * @author Cory Alexander Balaton (coryab) + * @author Janita Ovidie Sandtrøen Willumsen (janitaws) + * + * @version 0.1 + * + * @brief The definition of the Ising model. + * + * @bug No known bugs + * */ +#ifndef __ISING_MODEL__ +#define __ISING_MODEL__ + +#include "typedefs.hpp" + +#include +#include +#include + +// Faster modulo +#define INDEX(I, N) (I + N) % N + +// Indeces for the neighbor matrix. +#define UP 0 +#define LEFT 0 +#define DOWN 1 +#define RIGHT 1 + +/** @brief The Ising model in 2 dimensions. + * + * @details Here we set \f$ J = 1 \f$, and the Boltzmann constant + * \f$ k_B = 1 \f$. + * */ +class IsingModel { +private: + friend class IsingModelTest; + /** @brief \f$ L \cross L \f$ matrix where element $ x \in {-1, 1}$. + * */ + arma::Mat lattice; + + arma::Mat neighbors; + + /** @brief A hash map containing all possible energy changes. + * */ + std::unordered_map energy_diff; + + /** @brief Temperature + * */ + double T; + + uint L; + + /** @brief The current energy state. unit: \f$ J \f$. + * */ + int E; + + /** @brief The current magnetic strength. unit: Unitless. + * */ + int M; + + /** @brief Energy per spin. unit: \f$ J \f$. + * */ + double eps; + + /** @brief Magnetization per spin. unit: Unitless. + * */ + double m; + + /** @brief Initialize the lattice with a random distribution of 1s and + * -1s. + * */ + void initialize_lattice(); + + void initialize_neighbors(); + + /** @brief Initialize the hashmap with the correct values. + * */ + void initialize_energy_diff(); + + /** @brief Initialize the model. + * */ + void initialize_magnetization(); + + void initialize_energy(); + +public: + /** @brief Constructor for the Ising model. + * + * @param L The size of the lattice. + * */ + IsingModel(uint L, double T); + + /** @brief Constructor for the Ising model. + * + * @param L The size of the lattice. + * */ + IsingModel(uint L, double T, int val); + + void Metropolis(std::mt19937 engine); +}; + +#endif diff --git a/src/IsingModel.cpp b/src/IsingModel.cpp new file mode 100644 index 0000000..ef03f85 --- /dev/null +++ b/src/IsingModel.cpp @@ -0,0 +1,121 @@ +/** @file IsingModel.cpp + * + * @author Cory Alexander Balaton (coryab) + * @author Janita Ovidie Sandtrøen Willumsen (janitaws) + * + * @version 0.1 + * + * @brief The implementation of the Ising model + * + * @bug No known bugs + * */ + +#include "IsingModel.hpp" + +#include "constants.hpp" + +#include + +IsingModel::IsingModel(uint L, double T) +{ + this->L = L; + this->T = T; + this->initialize_lattice(); + this->initialize_neighbors(); + this->initialize_energy_diff(); + this->initialize_magnetization(); + this->initialize_energy(); +} + +IsingModel::IsingModel(uint L, double T, int val) +{ + this->L = L; + this->T = T; + this->lattice.set_size(this->L, this->L); + this->lattice.fill(val); + this->initialize_neighbors(); + this->initialize_energy_diff(); + this->initialize_magnetization(); + this->initialize_energy(); +} + +void IsingModel::initialize_lattice() +{ + this->lattice.set_size(this->L, this->L); + std::random_device rd{}; + std::mt19937 engine(rd()); + + std::uniform_int_distribution coin_flip(0, 1); + + for (size_t i = 0; i < this->lattice.n_elem; i++) + this->lattice(i) = coin_flip(engine) == 1 ? 1 : -1; +} + +void IsingModel::initialize_neighbors() +{ + this->neighbors.set_size(this->L, 2); + + for (int i = 0; i < this->L; i++) { + this->neighbors(i, UP) = INDEX(i - 1, this->L); + this->neighbors(i, DOWN) = INDEX(i + 1, this->L); + } + this->neighbors.print(); +} + +void IsingModel::initialize_energy_diff() +{ + for (size_t i = -8; i <= 8; i += 4) + this->energy_diff.insert({i, std::exp(-((double)i) / this->T)}); +} + +void IsingModel::initialize_magnetization() +{ + this->M = arma::accu(this->lattice); +} + +void IsingModel::initialize_energy() +{ + this->E = 0; + + // Loop through the matrix + for (size_t j = 0; j < this->L; j++) { + for (size_t i = 0; i < this->L; i++) { + this->E -= this->lattice(i, j) + * (this->lattice(i, this->neighbors(j, RIGHT)) + + this->lattice(this->neighbors(i, DOWN), j)); + } + } +} + +void IsingModel::Metropolis(std::mt19937 engine) +{ + uint ri, rj; + int dE; + + // Create random distribution for indeces + std::uniform_int_distribution random_index(0, this->L); + // Create random distribution for acceptance + std::uniform_real_distribution<> random_number(0., 1.); + + // Loop over the number of spins + for (int i = 0; i < this->lattice.n_elem; i++) { + // Get random indeces + ri = random_index(engine); + rj = random_index(engine); + + // Calculate the difference in energy + dE = 2. * this->lattice(ri, rj) + * (this->lattice(ri, this->neighbors(rj, LEFT)) + + this->lattice(ri, this->neighbors(rj, RIGHT)) + + this->lattice(this->neighbors(ri, UP), rj) + + this->lattice(this->neighbors(ri, DOWN), rj)); + + // Choose whether or not to accept the new configuration + if (random_number(engine) <= this->energy_diff[dE]) { + // Update if the configuration is accepted + this->lattice(ri, rj) *= -1; + this->M += 2 * this->lattice(ri, rj); + this->E += dE; + } + } +}