First draft of Ising model

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 <random>
+
+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<int> 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<uint> 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;
+        }
+    }
+}