Example: the Ising chain in a magnetic field

Here is the a simple Monte-Carlo for a one-dimensional Ising chain. The problem is described in detail in this section about the Ising model.

The configuration

We start by defining a configuration class on which the move and measure classes will act. We write this class in a file configuration.hpp:

#pragma once
// The configuration of the system
struct configuration {

  // N is the length of the chain, 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"
  std::vector<bool> chain;

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

The move

The move class should have three methods: attempt(), accept() and reject():

#pragma once
#include <triqs/mc_tools/random_generator.hpp>
#include <vector>
#include "configuration.hpp"

// A move flipping a random spin
struct flip {

  configuration *config;
  triqs::mc_tools::random_generator &RNG;

  int site;
  double delta_energy;

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

  double attempt() {
    // pick a random site
    site = RNG(config->N);

    // find the neighbours with periodicity
    int left  = (site == 0 ? config->N - 1 : site - 1);
    int right = (site == config->N - 1 ? 0 : site + 1);

    // compute energy difference from field
    delta_energy = (config->chain[site] ? 2 : -2) * config->field;

    // compute energy difference from J
    if (config->chain[left] == config->chain[right]) { delta_energy += (config->chain[left] == config->chain[site] ? 4 : -4) * config->J; }

    // 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[site] ? -2 : 2);
    config->chain[site] = !config->chain[site];
    config->energy += delta_energy;

    return 1.0;
  }

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

Measure

The measure class has two methods, accumulate and collect_results:

#pragma once
#include "configuration.hpp"

// The measure of the magnetization
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;
  }

  // get final answer M / (Z*N)
  void collect_results(mpi::communicator c) {
    double sum_Z = mpi::reduce(Z, c);
    double sum_M = mpi::reduce(M, c);
    if (c.rank() == 0) std::cout << "Magnetization: " << sum_M / sum_Z << std::endl;
  }
};

Main program

The Monte-Carlo itself can now be written:

#include <iostream>
#include <triqs/mc_tools/mc_generic.hpp>
#include <triqs/utility/callbacks.hpp>

#include "moves.hpp"
#include "configuration.hpp"
#include "measures.hpp"

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

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

  // greeting
  if (world.rank() == 0) std::cout << "Ising chain" << std::endl;

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

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

  // parameters of the model
  int length   = 100;
  double J     = -1.0;
  double field = 0.5;
  double beta  = 0.5;

  // 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);
}

Ising chain
Magnetization: 9.27603
Total number of measures: 500000