Моделирование модели Изинга смещает критическую температуру

Я пишу симуляцию модели Изинга в 2D. Модель ведет себя так, как и предполагалось, за исключением одного: критическая температура составляет примерно 3,5, в то время как она должна быть около 2/ln(2 + sqrt (2)).

Проект представляет собой программу на C++, которая генерирует данные, и сценарий оболочки, который выполняет эту программу. Полный код можно найти здесь. Также вот lattice.cpp

#include <iostream>
#include "include/lattice.h"

using namespace std;

/*
Copy assignment operator, too long to include in the header.
*/
lattice &lattice::operator=(const lattice &other) {
  size_ = other.size_;
  spins_ = other.spins_;
  J_ = other.J_;
  H_ = other.H_;
  delete spins_;
  return *this;
}

void lattice::print() {
  unsigned int area = size_ * size_;
  for (unsigned int i = 0; i < area; i++) {
    cout << to_symbol(spins_->at(i));
    if (i % size_ == size_ - 1)
      cout << endl;
  }
  cout << endl;
}

/*
Computes the energy associated with a spin at the given point.

It is explicitly float as that would allow the compiler to make use of multiple
registers instead of keeping track of unneeded precision.  (typically J, H ~ 1).
*/
float lattice::compute_point_energy(int row, int col) {
  int accumulator = get(row + 1, col) + get(row - 1, col) + get(row, col - 1) +
                    get(row, col + 1);
  return -get(row, col) * (accumulator * J_ + H_);
}

/*
Computes total magnetisation in O(n^2). Thread safe
*/
int lattice::total_magnetisation() {
  int sum = 0;
  #pragma omp parallel for reduction(+ : sum)
  for (unsigned int i = 0; i < size_ * size_; i++) {
    sum += spins_->at(i);
  }
  return sum;
}

int inline to_periodic(int row, int col, int size) {
  if (row < 0 || row >= size)
    row = abs(size - abs(row));
  if (col < 0 || col >= size)
    col = abs(size - abs(col));
  return row * size + col;
}

с решетка.h

#ifndef lattice_h
#define lattice_h

#include <cmath>
#include <vector>

/* Converts spin up/down to easily printable symbols. */
char inline to_symbol(int in) { return in == -1 ? '-' : '+'; }

/* Converts given pair of indices to those with periodic boundary conditions. */
int inline to_periodic(int row, int col, int size) {
  if (row < 0 || row >= size)
    row = abs(size - abs(row));
  if (col < 0 || col >= size)
    col = abs(size - abs(col));
  return row * size + col;
}

class lattice {
private:
  unsigned int size_;
  // vector<bool> would be more space efficient, but it would not allow
  // multithreading
  std::vector<short> *spins_;
  float J_;
  float H_;

public:
  lattice() noexcept : size_(0), spins_(NULL), J_(1.0), H_(0.0) {}
  lattice(int new_size, double new_J, double new_H) noexcept
      : size_(new_size), spins_(new std::vector<short>(size_ * size_, 1)),
        J_(new_J), H_(new_H) {}
  lattice(const lattice &other) noexcept
      : lattice(other.size_, other.J_, other.H_) {
#pragma omp parallel for
    for (unsigned int i = 0; i < size_ * size_; i++)
      spins_->at(i) = other.spins_->at(i);
  }
  lattice &operator=(const lattice &);

  ~lattice() { delete spins_; }
  void print();
  short get(int row, int col) {
    return spins_->at(to_periodic(row, col, size_));
  }
  unsigned int get_size() { return size_; }
  void flip(int row, int col) { spins_->at(to_periodic(row, col, size_)) *= -1; }
  int total_magnetisation();
  float compute_point_energy(int row, int col);
};

#endif

и simulation.cpp

#include <iostream>
#include <math.h>
#include "include/simulation.h"

using namespace std;

/*
Advances the simulation a given number of steps, and updates/prints the statistics
into the given file pointer.

Defaults to stdout.

The number of time_steps is explcitly unsigned, so that linters/IDEs remind
the end user of the file that extra care needs to be taken, as well as to allow
advancing the simulation a larger number of times.
*/
void simulation::advance(unsigned int time_steps, FILE *output) {
  unsigned int area = spin_lattice_.get_size() * spin_lattice_.get_size();
  for (unsigned int i = 0; i < time_steps; i++) {
    // If we don't update mean_energy_ every time, we might get incorrect
    // thermodynamic behaviour.
    total_energy_ = compute_energy(spin_lattice_);
    double temperature_delta = total_energy_/area - mean_energy_;
    if (abs(temperature_delta) < 1/area){
      cerr<<temperature_delta<<"! Reached equilibrium "<<endl;
    }
    temperature_ += temperature_delta;
    mean_energy_ = total_energy_ / area;
    if (time_ % print_interval_ == 0) {
      total_magnetisation_ = spin_lattice_.total_magnetisation();
      mean_magnetisation_ = total_magnetisation_ / area;
      print_status(output);
    }
    advance();
  }
}

/*
Advances the simulation a single step.

DOES NOT KEEP TRACK OF STATISTICS. Hence private.
*/
void simulation::advance() {
  #pragma omp parallel for collapse(2)
  for (unsigned int row = 0; row < spin_lattice_.get_size(); row++) {
    for (unsigned int col = 0; col < spin_lattice_.get_size(); col++) {
      double dE = compute_dE(row, col);
      double p = r_.random_uniform();
      float rnd = rand() / (RAND_MAX + 1.);
      if (exp(-dE / temperature_) > rnd) {
        spin_lattice_.flip(row, col);
      }
    }
  }
  time_++;
}

/*
Computes change in energy due to flipping one single spin.

The function returns a single-precision floating-point number, as data cannot under
most circumstances make use of greater precision than that (save J is set to a
non-machine-representable value).

The code modifies the spin lattice, as an alternative (copying the neighborhood
of a given point), would make the code run slower by a factor of 2.25
*/
float simulation::compute_dE(int row, int col) {
  float e_0 =  spin_lattice_.compute_point_energy(row, col);
  return -4*e_0;


}
/*
Computes the total energy associated with spins in the spin_lattice_.

I originally used this function to test the code that tracked energy as the lattice
itself was modified, but that code turned out to be only marginally faster, and
not thread-safe. This is due to a race condition: when one thread uses a neighborhood
of a point, while another thread was computing the energy of one such point in
the neighborhood of (row, col).
*/
double simulation::compute_energy(lattice &other) {
  double energy_sum = 0;
  unsigned int max = other.get_size();
  #pragma omp parallel for reduction(+ : energy_sum)
  for (unsigned int i = 0; i < max; i++) {
    for (unsigned int j = 0; j < max; j++) {
      energy_sum += other.compute_point_energy(i, j);
    }
  }
  return energy_sum;
}

void simulation::set_to_chequerboard(int step){
  if (time_ !=0){
    return;
  }else{
    for (unsigned int i=0; i< spin_lattice_.get_size(); ++i){
      for (unsigned int j=0; j<spin_lattice_.get_size(); ++j){
        if ((i/step)%2-(j/step)%2==0){
          spin_lattice_.flip(i, j);
        }
      }
    }
  }
}

с simulation.h

#ifndef simulation_h
#define simulation_h

#include "lattice.h"
#include "rng.h"
#include <gsl/gsl_rng.h>

/*
The logic of the entire simulation of the Ising model of magnetism.

This simulation will run and print statistics at a given time interval.
A simulation can be advanced a single time step, or many at a time,
*/
class simulation {
private:
  unsigned int time_ = 0;  // Current time of the simulation.
  rng r_ = rng();
  lattice spin_lattice_;
  double temperature_;
  double mean_magnetisation_ = 1;
  double mean_energy_;
  double total_magnetisation_;
  double total_energy_;
  unsigned int print_interval_ = 1;
  void advance();

public:
  void set_print_interval(unsigned int new_print_interval) { print_interval_ = new_print_interval; }

  simulation(int new_size, double new_temp, double new_J, double new_H)
      : time_(0), spin_lattice_(lattice(new_size, new_J, new_H)), temperature_(new_temp),
        mean_energy_(new_J * (-4)), total_magnetisation_(new_size * new_size),
        total_energy_(compute_energy(spin_lattice_)) {}

  void print_status(FILE *f) {
    f = f==NULL? stdout : f;
    fprintf(f, "%4d\t%e \t%e\t%e\n", time_, mean_magnetisation_,
            mean_energy_, temperature_);
  }
  void advance(unsigned int time_steps, FILE *output);
  double compute_energy(lattice &other);
  float compute_dE(int row, int col);
  void set_to_chequerboard(int step);
  void print_lattice(){
    spin_lattice_.print();
  };
  // void load_custom(const lattice& custom);
};

#endif

Результат сейчас выглядит примерно так: в то время как он должен быть ниже 2,26