# 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 202.

### 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 ','
@@


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.

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.

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 $\frac{X^\mathrm{(fit)} (t_\mathrm{final}) - X_0 (t_\mathrm{initial})}{\bar{X}}$ simplified False includeLog False