IsingModel.cpp

Go to the documentation of this file.

/** @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"
 
IsingModel::IsingModel()
{
    this->initialize_engine();
}
 
IsingModel::IsingModel(int L, double T)
{
    this->L = L;
    this->T = T;
    this->initialize_engine();
    this->initialize_lattice();
    this->initialize_neighbors();
    this->initialize_energy_diff();
    this->initialize_magnetization();
    this->initialize_energy();
}
 
IsingModel::IsingModel(int L, double T, int val)
{
    this->L = L;
    this->T = T;
    this->initialize_engine();
    this->initialize_lattice(val);
    this->initialize_neighbors();
    this->initialize_energy_diff();
    this->initialize_magnetization();
    this->initialize_energy();
}
 
void IsingModel::initialize_engine()
{
    std::random_device rd{};
    this->engine = std::mt19937{rd()};
}
 
void IsingModel::initialize_lattice()
{
    this->lattice.set_size(this->L, this->L);
 
    std::uniform_int_distribution<> coin_flip(0, 1);
 
    for (size_t i = 0; i < this->lattice.n_elem; i++)
        this->lattice(i) = 2 * coin_flip(this->engine) - 1;
}
 
void IsingModel::initialize_lattice(int val)
{
    // If val is neither 1 or -1, then initialize random values.
    if (val != 1 && val != -1) {
        this->initialize_lattice();
        return;
    }
    this->lattice.set_size(this->L, this->L);
    this->lattice.fill(val);
}
 
void IsingModel::initialize_neighbors()
{
    this->neighbors.set_size(this->L, 2);
 
    // Having i as a signed integer is necessary in this case.
    for (int i = 0; i < (int)this->L; i++) {
        this->neighbors(i, UP) = INDEX(i - 1, this->L);
        this->neighbors(i, DOWN) = INDEX(i + 1, this->L);
    }
}
 
void IsingModel::initialize_energy_diff()
{
    for (int i = -8; i <= 8; i += 4) {
        this->energy_diff[i+8] = std::exp(-(double)i / this->T);
    }
}
 
void IsingModel::initialize_magnetization()
{
    this->M = 0.;
    for (size_t i = 0; i < this->lattice.n_elem; i++) {
        this->M += this->lattice(i);
    }
}
 
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));
        }
    }
}
 
data_t IsingModel::Metropolis()
{
    int ri, rj;
    int dE;
 
    // Create random distribution for indeces
    std::uniform_int_distribution<> random_index(0, this->L - 1);
    // Create random distribution for acceptance
    std::uniform_real_distribution<> random_number(0., 1.);
 
    // Loop over the number of spins
    for (size_t 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+8]) {
            // Update if the configuration is accepted
            this->lattice(ri, rj) *= -1;
            this->M += 2 * this->lattice(ri, rj);
            this->E += dE;
        }
    }
 
    return data_t((double)this->E, (double)(this->E * this->E), (double)this->M,
                  (double)(this->M * this->M), std::fabs((double)this->M));
}