# Checking warmup and autocorrelation times#

Two common challenges in Monte Carlo studies are ensuring that the measured bins represent equilibrated configurations and that different bins are statistically independent. In this section, we will briefly explain these issues and present the tools pyALF offers for dealing with them.

## Preparations#

As a first step, we use the same import as in previous sections.

```
# Enable Matplotlib Jupyter Widget Backend
%matplotlib widget
import numpy as np # Numerical libary
import matplotlib.pyplot as plt # Plotting library
from py_alf.analysis import analysis # Analysis function
from py_alf.utils import find_sim_dirs # Function for finding Monte Carlo bins
from py_alf.ana import load_res # Function for loading analysis results
```

We also import the functions `py_alf.check_warmup()`

and `py_alf.check_rebin()`

, which play the main role in this section.

```
from py_alf import check_warmup, check_rebin
```

Finally, from the local file `custom_obs.py`

, we import the same `custom_obs`

defined in the previous section.

```
from custom_obs import custom_obs
```

For demonstration purposes, we run a simulation with very small bins.

```
from py_alf import ALF_source, Simulation
sim = Simulation(
ALF_source(),
'Hubbard',
{
# Model specific parameters
'L1': 4,
'L2': 4,
'Ham_U': 5.0,
# QMC parameters
'Nbin': 5000,
'Nsweep': 5,
'Ltau': 0,
},
)
sim.compile()
sim.run()
```

```
Compiling ALF...
Cleaning up Prog/
Cleaning up Libraries/
Cleaning up Analysis/
Compiling Libraries
```

```
Compiling Analysis
Compiling Program
Parsing Hamiltonian parameters
filename: Hamiltonians/Hamiltonian_Kondo_smod.F90
filename: Hamiltonians/Hamiltonian_Hubbard_smod.F90
filename: Hamiltonians/Hamiltonian_Hubbard_Plain_Vanilla_smod.F90
filename: Hamiltonians/Hamiltonian_tV_smod.F90
filename: Hamiltonians/Hamiltonian_LRC_smod.F90
filename: Hamiltonians/Hamiltonian_Z2_Matter_smod.F90
Compiling program modules
Link program
Done.
Prepare directory "/scratch/pyalf-docu/doc/source/usage/ALF_data/Hubbard_L1=4_L2=4_U=5.0" for Monte Carlo run.
Create new directory.
Run /home/jschwab/Programs/ALF/Prog/ALF.out
ALF Copyright (C) 2016 - 2021 The ALF project contributors
This Program comes with ABSOLUTELY NO WARRANTY; for details see license.GPL
This is free software, and you are welcome to redistribute it under certain conditions.
No initial configuration
```

Set the directories to be considered.

```
dirs = find_sim_dirs()
dirs
```

```
['./ALF_data/Hubbard',
'./ALF_data/Hubbard_L1=4_L2=4_U=1.0',
'./ALF_data/Hubbard_L1=4_L2=4_U=2.0',
'./ALF_data/Hubbard_L1=4_L2=4_U=3.0',
'./ALF_data/Hubbard_L1=4_L2=4_U=4.0',
'./ALF_data/Hubbard_L1=4_L2=4_U=5.0',
'./ALF_data/Hubbard_Square',
'./ALF_data/temper_Hubbard_L1=4_L2=4_U=2.5/Temp_0',
'./ALF_data/temper_Hubbard_L1=4_L2=4_U=2.5/Temp_1']
```

## Check warmup#

A Monte Carlo simulation creates a time series of configurations through stochastic updates. Usually, measurements from a number of updates get combined in one so-called bin. In the case of ALF, `Nsweep`

sweeps create one bin of measurements (for more details on updating procedures we refer to the ALF documentation). Usually, the simulation starts in a non-optimal state and it takes some time to reach equilibrium. Bins from this “warming up” period should be dismissed before calculating results. This is achieved by setting the variable `N_skip`

in the file `parameters`

, which will make the analysis omit the first `N_skip`

bins.

Warning

Different observables can have different warmup and autocorrelation times. For example, charge degrees of freedom may equilibrate much faster than spin degrees of freedom. Or a sum of observables might have much shorter autocorrelation times than an individual observable, e.g. the total spin versus one spin component.

To judge the correct value for `N_skip`

, pyALF offers the function `check_warmup()`

, which plots the time series of bins for a given list of scalar and custom observables. It can be used with the previous simulations as:

```
warmup_widget = check_warmup(
dirs,
['Ener_scal', 'Kin_scal', 'Pot_scal',
'E_pot_kin', 'R_Ferro', 'R_AFM',
'SpinZ_pipi', 'SpinXY_pipi', 'SpinXYZ_pipi'],
custom_obs=custom_obs, gui='ipy'
)
```

The first argument is a list of directories containing simulations, the second argument specifies which observables to plot, the keyword argument `custom_obs`

is needed, when plotting custom observables, e.g. *E_pot_kin* and `gui`

specifies which GUI framework to use. With `gui='ipy'`

, the function returns a Jupyther Widget, which allows to seamlessly work within Jupyter. Another option would be `gui='tk'`

, which open a separate window using `tkinter`

. The latter option might be suitable when working directly from a shell.

The variable `N_skip`

can be directly changed in the GUI, which automatically updates the file `parameters`

.

```
warmup_widget
```

## Check rebin#

When estimating statistical errors, the analysis assumes different bins to be statistically independent. As a result, one bin must span over enough updates to generate statically independent configurations, or in other words, a bin must be larger than the autocorrelation time. Otherwise the statistical errors will be underestimated. To address this issue, the analysis employs so-called rebinning, which combines `N_rebin`

bins into one new bin. The pyALF function `check_rebin()`

helps in determining the correct `N_rebin`

. It plots the errors of the chosen observables against `N_rebin`

. With enough statistics, one should see growing errors with increasing `N_rebin`

until a saturation point is reached, this saturation point marks a suitable value for `N_rebin`

. The usage of `check_rebin()`

is very similar to `check_warmup()`

.

```
rebin_widget = check_rebin(
dirs,
['Ener_scal', 'Kin_scal', 'Pot_scal',
'E_pot_kin', 'R_Ferro', 'R_AFM',
'SpinZ_pipi', 'SpinXY_pipi', 'SpinXYZ_pipi'],
custom_obs=custom_obs, gui='ipy'
)
```

Below, we can see how different observables have different autocorrelation times. While the error of the kinetic energy saturates already at \(N_\text{rebin} = 3\), the correlations of the z component of the spin at \((\pi, \pi)\) (`SpinZ_pipi`

) need \(N_\text{rebin} \sim 40\). For the correlations of the total spin, on the other hand, \(N_\text{rebin} = 1\) is sufficient.

This is a good example for the concept of an improved estimator: The simulated Hubbard model is \(\text{SU}(2)\) symmetric, therefore correlations of the x, y and z components of the spin are equivalent, but with the chosen parameter `Mz=True`

(cf. Compiling and running ALF) the auxiliary field couples to the z component of the spin. As a result, the \(\text{SU}(2)\) symmetry is broken for an individual auxiliary field configuration, but restored by sampling the field configurations. Therefore, even though measuring the spin correlations through the z component, the x-y plane, or the full spin are equivalent, the latter option produces the most precise results and has the shortest autocorrelation times, because it explicitly restores the \(\text{SU}(2)\) symmetry instead of “waiting” for the sampling to do that.

Furthermore, the ferromagnetic correlation ratio \(R_{\text{Ferro}}\) doesn’t seem to converge at all in the considered range of \(N_{\text{rebin}}\). This is connected to the fact that the system is not close the ferromagnetic order and therefore \(R_{\text{Ferro}}\) is a bad observable.

```
rebin_widget
```

The next section will also show options for an improved estimator by employing symmetry operations on the Bravais lattice.