Full example: Monte-Carlo simulation of the 2D Ising model

The following example is the full implementation of a Monte-Carlo sampling of the magnetization of a 2D Ising model with statistical analysis.

For more information about the workings of the Monte-Carlo worker, we refer the reader to the documentation/manual/triqs.

#include <triqs/mc_tools/random_generator.hpp>
#include <triqs/mc_tools/mc_generic.hpp>
#include <triqs/utility/callbacks.hpp>
#include <triqs/arrays.hpp>
#include <triqs/statistics.hpp>
#include <vector>
#include <iostream>
#include <fstream>
//#define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK
// H = -J \sum_<ij> s_i s_j - h \sum_i s_i
// theoretical T_c = 2/log(1+sqrt(2)) for J = 1.0
using namespace triqs::statistics;
/**************
 * config
 **************/

struct configuration {
  // N is the linear size of spin matrix, M the total magnetization,
  // beta the inverse temperature, J the coupling,
  // field the magnetic field and energy the energy of the configuration
  int N, M;
  double beta, J, field, energy;
  // the chain of spins: true means "up", false means "down"
  triqs::arrays::array<bool, 2> chain;
  observable<double> M_stack;

  // constructor
  configuration(int N_, double beta_, double J_, double field_)
     : N(N_), M(-N * N), beta(beta_), J(J_), field(field_), energy(-J * 4 * N / 2 + N * field), chain(N, N), M_stack() {
    chain() = false;
  }
};

/**************
 * move
 **************/

// A move flipping a random spin
struct flip {
  configuration *config;
  triqs::mc_tools::random_generator &RNG;

  struct site {
    int i, j;
  }; //small struct storing indices of a given site
  site s;
  double delta_energy;

  // constructor
  flip(configuration &config_, triqs::mc_tools::random_generator &RNG_) : config(&config_), RNG(RNG_) {}

  // find the neighbours with periodicity
  std::vector<site> neighbors(site s, int N) {
    std::vector<site> nns(4);
    int counter = 0;
    for (int i = -1; i <= 1; i++) {
      for (int j = -1; j <= 1; j++) {
        if ((i == 0) != (j == 0)) //xor
          nns[counter++] = site{(s.i + i + N) % N, (s.j + j + N) % N};
      }
    }
    return nns;
  }
  double attempt() {
    // pick a random site
    int index = RNG(config->N * config->N);
    s         = {index % config->N, index / config->N};

    // compute energy difference from field
    delta_energy         = (config->chain(s.i, s.j) ? 2 : -2) * config->field;
    auto nns             = neighbors(s, config->N); //nearest-neighbors
    double sum_neighbors = 0.0;
    for (auto &x : nns) sum_neighbors += ((config->chain(x.i, x.j)) ? 1 : -1);
    // compute energy difference from J
    delta_energy += -sum_neighbors * config->J * (config->chain(s.i, s.j) ? -2 : 2);

    // return Metroplis ratio
    return std::exp(-config->beta * delta_energy);
  }

  // if move accepted just flip site and update energy and magnetization
  double accept() {
    config->M += (config->chain(s.i, s.j) ? -2 : 2);
    config->chain(s.i, s.j) = !config->chain(s.i, s.j);
    config->energy += delta_energy;
    return 1.0;
  }

  // nothing to do if the move is rejected
  void reject() {}
};

/**************
 * measure
 **************/
struct compute_m {

  configuration *config;
  double Z, M;

  compute_m(configuration &config_) : config(&config_), Z(0), M(0) {}

  // accumulate Z and magnetization
  void accumulate(int sign) {

    Z += sign;
    M += config->M;
    //config->M_stack << double(config->M/(config->N*config->N));
    config->M_stack << config->M;
  }

  // get final answer M / (Z*N)
  void collect_results(mpi::communicator const &c) {

    double sum_Z = mpi_reduce(Z, c);
    double sum_M = mpi_reduce(M, c);

    if (c.rank() == 0) {
      std::cout << "#Beta:\t" << config->beta << "\tMagnetization:\t" << sum_M / (sum_Z * (config->N * config->N)) << std::endl;
      std::cout << "average_and_error(M) = " << average_and_error(config->M_stack) << std::endl;
      std::cout << "#Beta:\t" << config->beta << "\tAutocorr_time:\t" << autocorrelation_time_from_binning(config->M_stack) << std::endl;
      std::ofstream outfile("magnetization_series.dat");
      for (int i = 0; i < config->M_stack.size(); i++) outfile << config->M_stack[i] << std::endl;
      outfile.close();
    }
  }
};

int main(int argc, char *argv[]) {

  // initialize mpi
  mpi::environment env(argc, argv);
  mpi::communicator world;

  double H = 0.0, B = 0.5;
  int N  = 20;
  int nc = 100000;
  if (argc == 4) {
    H  = atof(argv[1]); //field
    B  = atof(argv[2]); //inverse temp
    N  = atoi(argv[3]); //size along one dimension
    nc = 1000000;
  }
  if (world.rank() == 0) std::cout << "2D Ising with field = " << H << ", beta = " << B << ", N = " << N << std::endl;

  // Prepare the MC parameters
  int n_cycles            = nc;
  int length_cycle        = 100;
  int n_warmup_cycles     = 100000;
  std::string random_name = "";
  int random_seed         = 374982 + world.rank() * 273894;
  int verbosity           = 0;

  // Construct a Monte Carlo loop
  triqs::mc_tools::mc_generic<double> IsingMC(random_name, random_seed, verbosity);

  // parameters of the model
  int length   = N;
  double J     = 1.0;
  double field = H;
  double beta  = B;

  // construct configuration
  configuration config(length, beta, J, field);

  // add moves and measures
  IsingMC.add_move(flip(config, IsingMC.get_rng()), "spin flip");
  IsingMC.add_measure(compute_m(config), "measure magnetization");

  // Run and collect results

  IsingMC.warmup_and_accumulate(n_warmup_cycles, n_cycles, length_cycle, triqs::utility::clock_callback(-1));
  IsingMC.collect_results(world);

  return 0;
}