Tutorials:Code-06 MonteCarloHOWTO

Revision as of 16:15, 25 May 2009 by Alet (talk | contribs)

Jump to: navigation, search

Work in progress ... The tutorial is neither finished nor corrected ...

The plan for this page is to explain how to write from scratch a Monte Carlo program using the ALPS libraries.

Please send me any comments / suggestions on this page


While several Monte Carlo (MC) programs are already available in the ALPS applications, clearly there is a need for helping developers to program their own MC program in the easiest and quickest possible way. Indeed there are several common tasks performed by all (Markov chain) MC algorithms which do not need to be re-programmed each time. The goal of this article is to convince that the ALPS libraries propose a simple framework to do so, and to demonstrate with examples how to use them. This includes both classical Monte Carlo as well as Quantum Monte Carlo algorithms.

What are the advantages of the ALPS libraries over one's own tools for developing a new MC application ?

  • Automatic calculation of error bars and autocorrelation times
  • Automatic parallelization
  • Automatic output of results (in XML) with the corresponding extraction tools
  • Automatic management of the sign (for Quantum Monte Carlo simulations with a sign problem)
  • Easy checkpointing & Easy restarting of simulations
  • Easy management of input parameters
  • Easy way to add new measurements
  • Easy management of random numbers
  • (Optional) Easy management of lattices
  • (Optional) Easy management of quantum Hamiltonians (for Quantum Monte Carlo)
  • ...

All this comes with only a small amount of programming/modification that we describe now.

Starting your own ALPS MC program

First of all, using the ALPS libraries implies that you use C++ as the programming language. From now on, let us assume that you want to program a new MC application and that you already know what are going to be your internal data structure and the algorithm you will use for a single Monte Carlo step.

To use the ALPS libraries, you will need to create your own C++ class, that will be derived from an internal ALPS class. Let's call this class MyMonteCarlo and create a MyMonteCarlo.hpp header file.

 #ifndef MYMC_HPP
 #define MYMC_HPP
 #include <alps/scheduler.h>
 #include <alps/alea.h>
 #include <alps/scheduler/montecarlo.h>
 #include <alps/osiris.h>
 #include <alps/osiris/dump.h>
 #include <alps/expression.h>
 using namespace std;
 class MyMonteCarlo : public alps::scheduler::MCRun{
 public :
   MyMonteCarlo(const alps::ProcessList &,const alps::Parameters &,int);
   static void print_copyright(std::ostream &);
   void save(alps::ODump &) const;
   void load(alps::IDump &);
   void dostep();
   bool is_thermalized() const;
   double work_done() const;
 private :
 // your own internal data here ...

We now see what all these lines mean.

  • First of all you need to include a few ALPS headers (some more will go here for other ALPS functionalities).
  • Then we will define our MyMonteCarlo class out by deriving it from a ALPS class. The above header file implies that you will use the MCRun ALPS class, which is the simplest one. If you want to enjoy the lattice functionalities of ALPS, replace the class definition line by:
 class MyMonteCarlo : public alps::scheduler::LatticeMCRun<>{

and add in the headers

 #include <alps/lattice.h>

If you moreover want to use the Model library (useful for lattice quantum models), use instead:

 class MyMonteCarlo : public alps::scheduler::LatticeModelMCRun<>{
  • In the public part, MyMonteCarlo(const alps::ProcessList&,const alps::Parameters&,int) is the constructor of your class, and will be needed to initialize all parameters.
  • print_copyright(std::ostream&) is a simple function to output any useful information that you want to be output at the beginning of each simulation.
  • The save and load functions are very useful functions where you will describe what is needed to save (load) on disk to restart the simulations after each checkpoint
  • The dostep function is your main MC function: this is the function that is going to be executed every MC step.
  • The is_thermalized function will tell the ALPS libraries when the thermalization part of your simulation is finished (that is, when can the measurement series start)
  • The work_done function will tell the ALPS libraries what is the percentage of the simulations that is already done.

The rest of the job is now to correctly interface these functions with your program. This will be done in a file called, say, MyMonteCarlo.cpp

Building your own ALPS MC program

The copyright function

We start with the easiest function print_copyright() in the file MyMonteCarlo.cpp

 #include "MyMonteCarlo.hpp";
 /************************************ ALPS functions **********************************************/
 // Copyright statement
 void MyMonteCarlo::print_copyright(std::ostream & out)
   out << "My own ALPS Monte Carlo program v. 0.1\n"
       << "  copyright (c) 2006 by Myself,\n"
       << "  available from the author on request\n\n";

Management of Monte Carlo steps

Let us now concentrate on the management of MC steps. A typical situation is the following: one would like to perform a fixed number of steps for thermalization, then a fixed number of steps for the measurement part. Very often, one would like also to perform a certain amount of Monte Carlo steps between each measurement. Therefore it could be useful to define in your internal data structure the following variables (in MyMonteCarlo.hpp):

 private :
  // your own internal data here ...
  int Nb_Steps; 
  int Nb_Thermalisation_Steps; 
  int Each_Measurement;
  int Steps_Done_Total; 
  int Measurements_Done;
  void do_update();
  void do_measurements();
  // the rest of your own internal data here ...

Here Nb_Thermalisation_Steps will be the number of thermalization steps that you ask for and Nb_Steps the number of steps after thermalization that you ask for. Each_Measurement will be the number of steps between each measurement. These three numbers will be initialized later in the constructor. Steps_Done_Total will store the current number of MC finished steps (including those of thermalization). Finally Measurements_Done will store the intermediate number of steps done between each measurement. The do_update() and do_measurements() function (that you will have to define yourself) will perform respectively one single MC step and one series of measurement.

All these definitions lead to the following simple definitions of your class member functions in MyMonteCarlo.cpp

 bool MyMonteCarlo::is_thermalized() const
   {  return (Steps_Done_Total >= Nb_Thermalisation_Steps); }

This function will indeed return 1 if the current number of MC steps is superior to the total number of thermalization steps asked for , 0 otherwise. The second function

 double MyMonteCarlo::work_done() const
 { return (is_thermalized() ? (Steps_Done_Total-Nb_Thermalisation_Steps)/double(Nb_Steps) :0.); }

will return 0 if the simulation is not thermalized, and a number between 0 and 1 corresponding to the percentage of asked measurement steps already performed. Finally the dostep() function will look like this

 void MyMonteCarlo::dostep()
 { do_update(); // you'll have to define what this function does later
   ++Steps_Done_Total; // increment the number of steps done
   if (is_thermalized()) // do measurements only if simulation thermalized
     { if (++Measurements_Done==Each_Measurement) // do a measurement every Each_Measurement
       { Measurements_Done=0;
       do_measurements(); // you'll have to define what this function does later

The save and load functions

ALPS allows an easy treatment of the internal data to be checkpointed to disk. Let us imagine that in your internal data structure you need to checkpoint a few integers and an array of double to be able to restart your simulations

 private :
  // the rest of your own internal data here ...
  int Number_of_Spins; 
  int MyOwnVariable;
  std::vector<double> SpinArray;