# ALPS 2 Tutorials:MC-01 Equilibration

## Contents

# Equilibration

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

## Example: Quantum Monte Carlo (directed worm algorithm) simulations

As an example, we consider a Quantum Monte Carlo simulation implemented in the directed worm algorithm for boson Hubbard model in square lattice geometry of size 20^{2}.

### Using command line

The parameter file parm1a:

LATTICE="square lattice" MODEL="boson Hubbard" L=20 Nmax=20 t=1. U=16. mu=32. THERMALIZATION=10000 SWEEPS=100000 SKIP=400 {T=1.0}

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

parameter2xml parm1a dwa parm1a.in.xml

Detailed information regarding the *Density* measurements, for example, can be extracted from:

h5dump -g /simulation/results/Density parm1a.task1.out.h5

or its (binned) timeseries from:

h5dump -g /simulation/results/Density/timeseries parm1a.task1.out.h5

We can extract the timeseries data into a CSV file:

h5dump -d /simulation/results/Density/timeseries/data -w 1 -y -o parm1a.task1.out.density.timeseries.csv parm1a.task1.out.h5

and plot it using your favourite graphical plotter, say xmgrace:

xmgrace parm1a.task1.out.density.timeseries.csv

or, say gnuplot:

gnuplot <<@@ set datafile separator ',' plot "parm1a.task1.out.density.timeseries.csv" @@

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 tutorial1a.py.

The headers:

import pyalps

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

parms = [{ 'LATTICE' : "square lattice", 'MODEL' : "boson Hubbard", 'L' : 20, 'Nmax' : 20, 't' : 1., 'U' : 16., 'mu' : 32., 'T' : 1., 'THERMALIZATION' : 10000, 'SWEEPS' : 100000, 'SKIP' : 400 }]

Write into XML input file:

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

and run the application **dwa**:

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

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

files = pyalps.getResultFiles(prefix='parm1a', format='hdf5')

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 plt.plot(density_timeseries) plt.show()

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'])

**Description**

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 | 0.01 |

simplified | False | |

includeLog | False |