Autocorrelation time - ALF Tutorial

In this example we use the pyALF interface to run ALF on a $4\times 4$ Hubbard square lattice, with the Mz choice of Hubbard-Stratonovich transformation, at $U/t=4$ , half-band filling, and inverse temperature $\beta t = 5$ .

Consecutive Monte Carlo bins are, in general, correlated: the update scheme performs local moves, so neighbouring configurations (and thus neighbouring bins) tend to be similar. If this autocorrelation is ignored, the naive bin-to-bin error underestimates the true statistical error.

The error estimate can be made robust against this by rebinning: grouping N_rebin consecutive raw bins into a single effective bin before performing the analysis. As N_rebin increases past the autocorrelation time (measured in units of bins), consecutive rebinned blocks become independent, and the estimated error saturates to its true value.

Bellow we go through the steps of running a simulation and estimating the autocorrelation time of the internal energy and spin-spin correlations from the saturation of the error as a function of N_rebin by the help of the pyalf widget check_rebin.

1. Import Simulation class from the py_alf python module, which provides the interface with ALF, as well as mathematics and plotting packages:

from py_alf import Simulation, ALF_source  # Interface with ALF
from py_alf.utils import find_sim_dirs    # Function for finding QMC bins
from py_alf import check_rebin          # Function for checking warmup time
#
%matplotlib widget
#
import numpy as np                        # Numerical library
import matplotlib.pyplot as plt           # Plotting library
alf_src = ALF_source(branch='master')     # Obtain ALF source code

2. Create an instance of Simulation, specifying the necessary parameters. We deliberately use a small NSweep (number of sweeps per bin), so that consecutive bins remain significantly autocorrelated, making the rebinning analysis below meaningful:

sim = Simulation(
    alf_src,
    'Hubbard',                    # Hamiltonian
    {                             # Model and simulation parameters
        'Lattice_type': 'Square', #    Lattice type
        'L1': 4,                  #    Lattice length in the first unit vector direction
        'L2': 4,                  #    Lattice length in the second unit vector direction
        'Checkerboard': False,    #    Whether checkerboard decomposition is used or not
        'Symm': True,             #    Whether symmetrization takes place
        'ham_T': 1.0,             #    Hopping parameter
        'ham_U': 4.0,             #    Hubbard interaction
        'beta': 5.0,              #    Inverse temperature
        'Dtau': 0.1,              #    Trotter time step, Ltrot=beta/Dtau
        'Mz': True,               #    If true, sets the M_z-Hubbard model: Nf=2, N_sum=1,
                                  #             HS field couples to z-component of magnetization
        'NSweep': 5,              #    Number of sweeps per bin (kept small on purpose)
        'NBin': 2000,             #    Number of bins
    },
)

3. Compile ALF, downloading it first if not found locally. This may take a few minutes:

sim.compile()                            # Compilation needs to be performed only once

Checking out branch master
Your branch is up to date with 'origin/master'.
Compiling ALF... 
Cleaning up Prog/
Cleaning up Libraries/
Cleaning up Analysis/
Compiling Libraries

Already on 'master'
entanglement_mod.F90:35:2:

   35 | #warning "You are compiling entanglement without MPI. No Renyi entropy results possible, all other observables still work!"
      |  1~~~~~~
Warning: #warning "You are compiling entanglement without MPI. No Renyi entropy results possible, all other observables still work!" [-Wcpp]
ar: creating modules_90.a
ar: creating libqrref.a

Compiling Analysis
Compiling Program
Compiling program modules
Parsing Hamiltonian parameters
filenames: Hamiltonians/Hamiltonian_Kondo_smod.F90 Hamiltonians/Hamiltonian_Kondo_read_write_parameters.F90
filenames: Hamiltonians/Hamiltonian_Hubbard_smod.F90 Hamiltonians/Hamiltonian_Hubbard_read_write_parameters.F90
filenames: Hamiltonians/Hamiltonian_Hubbard_Plain_Vanilla_smod.F90 Hamiltonians/Hamiltonian_Hubbard_Plain_Vanilla_read_write_parameters.F90
filenames: Hamiltonians/Hamiltonian_tV_smod.F90 Hamiltonians/Hamiltonian_tV_read_write_parameters.F90
filenames: Hamiltonians/Hamiltonian_LRC_smod.F90 Hamiltonians/Hamiltonian_LRC_read_write_parameters.F90
filenames: Hamiltonians/Hamiltonian_Z2_Matter_smod.F90 Hamiltonians/Hamiltonian_Z2_Matter_read_write_parameters.F90
filenames: Hamiltonians/Hamiltonian_Spin_Peierls_smod.F90 Hamiltonians/Hamiltonian_Spin_Peierls_read_write_parameters.F90
Link program
Done.

4. Perform the simulation and the default analysis:

sim.run()                                # Perform the actual simulation in ALF
sim.analysis()                           # Perform default (N_rebin=1) analysis
obs = sim.get_obs()
obs.iloc[0][['Ener_scal0', 'Ener_scal0_err']]

Prepare directory "/home/hparussel/ALF_Tutorial/jupyterbook/source/05-caveats/ALF_data/Hubbard_Square_L1=4_L2=4_Checkerboard=False_Symm=True_T=1.0_U=4.0_beta=5.0_Dtau=0.1_Mz=True" for Monte Carlo run.
Resuming previous run.
Run /home/hparussel/ALF_Tutorial/jupyterbook/source/05-caveats/ALF/Prog/ALF.out
 ALF Copyright (C) 2016 - 2022 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.

 Setting N_FL = 2 since this is mandatory for Mz-Hubbard model

### Analyzing /home/hparussel/ALF_Tutorial/jupyterbook/source/05-caveats/ALF_data/Hubbard_Square_L1=4_L2=4_Checkerboard=False_Symm=True_T=1.0_U=4.0_beta=5.0_Dtau=0.1_Mz=True ###
/home/hparussel/ALF_Tutorial/jupyterbook/source/05-caveats
Scalar observables:
Ener_scal
Kin_scal
Part_scal
Pot_scal
Histogram observables:
Equal time observables:
Den_eq
Green_eq
SpinT_eq
SpinXY_eq
SpinZ_eq
Time displaced observables:
Den_tau
Green_tau
SpinT_tau
SpinXY_tau
SpinZ_tau
/home/hparussel/ALF_Tutorial/jupyterbook/source/05-caveats/ALF_data/Hubbard_Square_L1=4_L2=4_Checkerboard=False_Symm=True_T=1.0_U=4.0_beta=5.0_Dtau=0.1_Mz=True
No orbital locations saved.

Ener_scal0       -13.185757
Ener_scal0_err     0.007366
Name: /home/hparussel/ALF_Tutorial/jupyterbook/source/05-caveats/ALF_data/Hubbard_Square_L1=4_L2=4_Checkerboard=False_Symm=True_T=1.0_U=4.0_beta=5.0_Dtau=0.1_Mz=True, dtype: object

5. Estimate the autocorrelation time of the internal energy.

We read the raw (unrebinned) Monte Carlo bins of Ener_scal, Kin_scal, Pot_scal, $S_{Z} (\pi, \pi)$ , $S_{XY} (\pi, \pi)$ and $S_{XYZ} (\pi, \pi)$ directly and repeatedly define the custom observables for the spin. Then, we apply the pyalf widget check_rebin which depicts the quantity’s error as a function of N_rebin. N_rebin is the number of raw bins combined into one effective bin. The resulting error should grow while N_rebin increases and then saturate once N_rebin exceeds the autocorrelation time (in units of bins):

custom_obs = {}  # defining custom_obs


def obs_k(obs, back, sign, N_orb, N_tau, dtau, latt,
          ks=((0., 0.),), mat=None):
    """Mean of correlation function at one, or multiple k-points.

    Calculates integral over tau (=susceptibility) if time-displaced
    correlation is supplied.

    Parameters
    ----------
    obs : array of shape (N_orb, N_orb, N_tau, latt.N)
        Correlation function, the background is already subtracted.
    back : array of shape (N_orb,)
        Background of Correlation function.
    sign : float
        Monte Carlo sign.
    N_orb : int
        Number of orbitals per unit cell.
    N_tau : int
        Number of imaginary time slices. 1 for equal-time correlations.
    dtau : float
        Imaginary time step.
    latt : py_alf.Lattice
        Bravais lattice object.
    ks : list of k-points, default=[(0., 0.)]
    mat : array of shape (N_orb, N_orb), default=None
        Orbital structure. Default: Trace over orbitals.

    """
    if mat is None:
        mat = np.identity(N_orb)
    out = 0
    for k in ks:
        n = latt.k_to_n(k)

        if N_tau == 1:
            out += (obs[:, :, 0, n] * mat).sum()
        else:
            out += (obs[..., n].sum(axis=-1) * mat).sum()

    return [out / len(ks)]


# Correlation of Spin z-component at k=(pi, pi)
custom_obs['SpinZ_pipi']= {
'needs': ['SpinZ_eq'],
'function': obs_k,
'kwargs': {'ks': [(np.pi, np.pi)]}
}
# Correlation of Spin x+y-component at k=(pi, pi)
custom_obs['SpinXY_pipi']= {
'needs': ['SpinXY_eq'],
'function': obs_k,
'kwargs': {'ks': [(np.pi, np.pi)]}
}

# Correlation of total Spin at k=(pi, pi)
custom_obs['SpinXYZ_pipi']= {
'needs': ['SpinT_eq'],
'function': obs_k,
'kwargs': {'ks': [(np.pi, np.pi)]}
}

dirs = find_sim_dirs()  # Setting the directories to be considered

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

rebin_widget

The error starts small (the naive, un-rebinned estimate underestimates the true error), rises as N_rebin grows, and then flattens out into a plateau. The value of N_rebin at which the curve saturates is an estimate, in units of bins, of the integrated autocorrelation time of the energy. Once N_rebin is chosen at or beyond this plateau (e.g., by setting N_rebin in the simulation parameters or within this widget), the error becomes reliable.

The widget check_rebin automatically updates the corresponding N_rebin value (here e.g. N_rebin=20) for the analysis in the parameters file.

5.4.1Exercises¶

Reduce NSweep further (e.g., to 1) and observe how the plateau shifts to larger N_rebin, indicating a longer autocorrelation time.
Increase NSweep (e.g., to 50 or 100) and check that the error is now roughly flat already at N_rebin=1, i.e., that individual bins are already effectively independent.
Repeat the analysis for a different observable (e.g., Part_scal) to see whether its autocorrelation time differs from that of the energy.

5.4 Autocorrelation time

5.4.1Exercises¶