ALPS 2 Tutorials:MC-01 Equilibration

Revision as of 12:55, 5 September 2013 by Tamama (talk | contribs) (Using Python)

Jump to: navigation, search


Rule of thumb: All Monte Carlo simulations have to be equilibrated before taking measurements.

Example: Classical Monte Carlo (local updates) simulations

As an example, we will implement a classical Monte Carlo simulation implemented in the Ising model on a finite square lattice of size 482.

Using command line

The parameter file parm1a:

LATTICE="square lattice"

We first convert the input parameters to XML and then run the application spinmc:

parameter2xml parm1a
spinmc --Tmin 10 --write-xml

Add in timeseries analysis here after Python

Based on the timeseries, the user will then judge for himself/herself whether the simulation has reached equilibration.

Using Python

The following describes what is going on within the script file

The headers:

import pyalps

Set up a python list of parameters (python) dictionaries:

parms = [{
  'LATTICE'         : "square lattice",
  'MODEL'           : "Ising",
  'L'               : 48,
  'J'               : 1.,
  'T'               : 2.269186,
  'THERMALIZATION'  : 10000,
  'SWEEPS'          : 50000,

Write into XML input file:

input_file = pyalps.writeInputFiles('parm1a',parms)

and run the application dwa:

pyalps.runApplication('spinmc', input_file, Tmin=10, writexml=True)

We first get the list of all hdf5 result files via:

files = pyalps.getResultFiles(prefix='parm1a')

and then extract, say the timeseries of the Density measurements:

ar = pyalps.hdf5.h5ar(files[0])
density_timeseries = ar['/simulation/results']['Density']['timeseries']['data']

We can then visualize graphically:

import matplotlib.pyplot as plt

Based on the timeseries, the user will then judge for himself/herself whether the simulation has reached equilibration.

A convenient tool: steady_state_check

ALPS Python provides a convenient tool to check whether a measurement observable(s) has (have) reached steady state equilibrium.

Here is one example:

pyalps.steady_state_check(files[0], 'Density')

and another one:

pyalps.steady_state_check(files[0], ['Density', 'Energy Density'])

1. steady_state_check first performs a linear fit on the timeseries, and decides whether the measurement observable has reached steady state equilibrium based on the gradient/slope of the fitted line.

2. The optional arguments of steady_state_check are:

argument default remark
tolerance 0.01  \mathrm{tolerance} = \frac{X^\mathrm{(fit)} (t_\mathrm{final}) - X^\mathrm{(fit)} (t_\mathrm{initial})}{\bar{X}}
simplified False shall we combine the checks of all observables as 1 final boolean answer?
includeLog False shall we print the detailed log?

3. To see the complete log for instance:

pyalps.steady_state_check(files[0], ['Density', 'Energy Density'], includeLog=True)

Using Vistrails

To run the simulation in Vistrails open the file mc-01b-equilibration.vt.