Basic analysis

Basic analysis#

As already shown in Minimal example, the basic analysis can be executed through py_alf.Simulation.analysis(), which in turn calls py_alf.analysis(). This section demonstrates how to directly use the latter function and how to access and work with analysis results.

As a first step, some libraries and functions are imported. The command %matplotlib widget enables the Matplotlib Jupyter Widget Backend, which is not necessary in this part, but for the functions used in Checking warmup and autocorrelation times, therefore it makes sense to establish it as a default.

# Enable Matplotlib Jupyter Widget Backend
%matplotlib widget

# Imports
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

The function find_sim_dirs() returns a list of all directories containing a file named data.h5, the file containing all Monte Carlo measurements if ALF has been compiled with HDF5. We use it to conveniently list all simulations run in the previous sections.

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

Looping over this list, we call analysis() for each directory. The function reads Monte Carlo bins from data.h5, or if this file does not exist alternatively from all files ending in _scal, _eq and _tau. Furthermore, n_skip and n_rebin are read from file parameters. The analysis omits the first n_skip bins and combines n_rebin original bins into a new one[1]. On the resulting bins, Jackknife resampling [5] is applied to estimate expectation values and their standard error. For brevity, the resulting printout is truncated.

for directory in dirs:
    analysis(directory)

### Analyzing ./ALF_data/Hubbard ###
/scratch/pyalf-docu/doc/source/usage
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
### Analyzing ./ALF_data/Hubbard_L1=4_L2=4_U=1.0 ###
/scratch/pyalf-docu/doc/source/usage
Scalar observables:
...
...
...

Get analysis results#

The analysis results are saved in each simulation directory, both in plain text in the folder res and as a pickled Python dictionary in the file res.pkl.

The binary data from multiple res.pkl files can be conveniently read with load_res(), which returns a single pandas DataFrame, a tabular data structure. It not only contains analysis results, but also the Hamiltonian-specific parameters. The parameter names are in all lower case.

res = load_res(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_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

The DataFrame has one row per simulation directory, which is also used as the index:

res.index

Index(['./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_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'],
      dtype='object')

Column indices can be accessed through:

res.columns

Index(['continuous', 'ham_chem', 'ham_t', 'ham_t2', 'ham_tperp', 'ham_u',
       'ham_u2', 'mz', 'l1', 'l2',
       ...
       'SpinXY_tau_lattice', 'SpinZ_tauK', 'SpinZ_tauK_err', 'SpinZ_tauR',
       'SpinZ_tauR_err', 'SpinZ_tau_lattice', 'Acc_Temp_scal_sign',
       'Acc_Temp_scal_sign_err', 'Acc_Temp_scal0', 'Acc_Temp_scal0_err'],
      dtype='object', length=114)

In the following, we will only use results from one simulation, corresponding to one row in the DataFrame. It is selected with:

item = res.loc['./ALF_data/Hubbard']

Which is equivalent to

item = res.iloc[0]

Most, but not all of the same data is also stored in plain text form in the folder ALF_data/Hubbard/res:

ls ALF_data/Hubbard/res

Den_eq_K      Green_eq_K      Part_scal       SpinT_tau/       SpinZ_eq_K
Den_eq_K_sum  Green_eq_K_sum  Pot_scal        SpinXY_eq_K      SpinZ_eq_K_sum
Den_eq_R      Green_eq_R      SpinT_eq_K      SpinXY_eq_K_sum  SpinZ_eq_R
Den_eq_R_sum  Green_eq_R_sum  SpinT_eq_K_sum  SpinXY_eq_R      SpinZ_eq_R_sum
Den_tau/      Green_tau/      SpinT_eq_R      SpinXY_eq_R_sum  SpinZ_tau/
Ener_scal     Kin_scal        SpinT_eq_R_sum  SpinXY_tau/

Scalar observables#

Scalar observable results are stored as a number of separate scalars, storing the sign, observable expectation value and their statistical errors. Here are, for example, the results for the internal energy Ener_scal:

for i in item.keys():
    if i.startswith('Ener_scal'):
        print(i, item[i])

Ener_scal_sign 1.0
Ener_scal_sign_err 0.0
Ener_scal0 -29.983503005734427
Ener_scal0_err 0.23268481534861885

Note the 0 in Ener_scal0 and Ener_scal0_err. This is the index in the vector of observables Ener_scal, since a scalar observable can hold a vector of scalars.

The same data is present in this plain text file:

!cat ALF_data/Hubbard/res/Ener_scal

# Sign: 1.0 0.0
-2.998350300573442695e+01 2.326848153486188453e-01

Example#

Here is a simple example that demonstrates the convenience of working with pandas DataFrames. We select out of all simulations the one with \(L_1 = 4\) and plot their internal energy against The value of the Hubbard \(U\).

# Create figure with axis labels
fig, ax = plt.subplots()
ax.set_xlabel(r'Hubbard interaction $U$')
ax.set_ylabel(r'Internal energy $E$')

# Select only rows with l1==4 and sort by ham_u
df = res[res.l1 == 4].sort_values(by='ham_u')

# Plot data
ax.errorbar(df.ham_u, df.Ener_scal0, df.Ener_scal0_err);

Equal-time correlation functions#

ALF and pyALF offer support for correlation functions of the form

(1)#\[\begin{split}C(\boldsymbol{r}, n_1, n_2) &= \frac{1}{N_r} \sum_{\boldsymbol{r}_0} \left\langle O(\boldsymbol{r}_0, n_1) O(\boldsymbol{r}_0 + \boldsymbol{r}, n_2) \right\rangle - \left\langle O(n_1) \right\rangle \left\langle O(n_2) \right\rangle \\ C(\boldsymbol{k}, n_1, n_2) &= \frac{1}{N_r} \sum_{\boldsymbol{r}} e^{i \boldsymbol{k} \boldsymbol{r}} C(\boldsymbol{r}, n_1, n_2)\end{split}\]

Where the sums go over the unit cells of the finite size Bravais lattice, \(N_r\) is the number of unit cells and \(n_1\), \(n_2\) denominate the orbitals within a unit cell.

Each observable produces a set of members in the results, these are for example the ones for the equal-time Green’s function:

for i in item.keys():
    if i.startswith('Green_eq'):
        print(i, np.shape(item[i]))

Green_eqK (1, 1, 36)
Green_eqK_err (1, 1, 36)
Green_eqK_sum (36,)
Green_eqK_sum_err (36,)
Green_eqR (1, 1, 36)
Green_eqR_err (1, 1, 36)
Green_eqR_sum (36,)
Green_eqR_sum_err (36,)
Green_eq_lattice ()

Members ending in K, K_err, R and R_err correspond to Eq. (1) and their errors. They have the shape \((N_\text{orb}, N_\text{orb}, N_r)\), where \(N_\text{orb}\) is the number of orbitals per unit cell. The objects ending in _sum have been traced over the orbital degrees of freedom. To correctly interpret the index over the unit cells, the member ending in _lattice is a dictionary containing the parameters for creating a Bravais lattice object py_alf.Lattice:

item['Green_eq_lattice']

{'L1': array([6., 0.]),
 'L2': array([0., 6.]),
 'a1': array([1., 0.]),
 'a2': array([0., 1.])}

from py_alf import Lattice
latt = Lattice(item['Green_eq_lattice'])

Here is, for example, the equal-time Greens function at \(\boldsymbol{k}= (\pi,\pi)\) with its error:

n = latt.k_to_n([np.pi, np.pi])
print(item.Green_eqK_sum[n], item.Green_eqK_sum_err[n])

0.0655750236373261 0.002778042249941216

The lattice object offers functions for conveniently plotting correlation functions in real and momentum space. Below, we plot the Spin-Spin correlations in real and momentum space, showing signs for antiferromagnetic order.

latt.plot_r(item.SpinZ_eqR_sum)

latt.plot_k(item.SpinZ_eqK_sum)

The plain text result files ending in _K and _R contain momentum and real-space resolved correlations, respectively. Here is an excerpt from the Greens function in momentum space:

!head -n 3 ALF_data/Hubbard/res/Green_eq_K

#    kx       ky   	              (0, 0)              	         trace over n_orb         
-2.09440 -2.09440	 1.4876423115e-01  6.9686819699e-03	 1.4876423115e-01  6.9686819699e-03
-2.09440 -1.04720	 1.0237360120e+00  1.6136697811e-02	 1.0237360120e+00  1.6136697811e-02

Where (0, 0) refers to the orbital indices. Since there is only one orbital per unit cell, this is the only combination and identical to the trace over all orbitals. The first to columns are the coordinates and the come alternatingly expectation value and standard error.

Time-displaced correlation functions#

The structure for time-displaced correlation functions is very similar to equal-time correlations, but by default only the trace over the orbital degrees of freedom is stored. These are the results for the time-displaced Green function:

for i in item.keys():
    if i.startswith('Green_tau'):
        print(i, np.shape(item[i]))

Green_tauK (51, 36)
Green_tauK_err (51, 36)
Green_tauR (51, 36)
Green_tauR_err (51, 36)
Green_tau_lattice ()

Here we plot the time-displaced Greens function at \(\boldsymbol{r}= 0\):

# Create figure with axis labels and logscale on y-axis
fig, ax = plt.subplots()
ax.set_xlabel(r'$\tau$')
ax.set_ylabel(r'$G(r=0, \tau)$')
ax.set_yscale('log')

# Create lattice object
latt = Lattice(item['Green_tau_lattice'])

# Get index corresponding to r=0
n = latt.r_to_n([0, 0])

# Plot data
ax.errorbar(
    item.dtau*range(len(item.Green_tauR[:, n])), 
    item.Green_tauR[:, n],
    item.Green_tauK_err[:, n]
)

<ErrorbarContainer object of 3 artists>

Again, plain text results data are available in the folder res. There is a separate folder for each \(\boldsymbol{k}\)-point and the data for \(\boldsymbol{r}=0\):

ls ALF_data/Hubbard/res/Green_tau

00_0.00/   1.05_0.00/    1.05_-2.09/   -2.09_1.05/   -2.09_3.14/  3.14_3.14/
00_-1.05/  -1.05_-1.05/  1.05_2.09/    2.09_-1.05/   2.09_3.14/   R0
00_1.05/   -1.05_1.05/   -1.05_3.14/   2.09_1.05/    3.14_0.00/
00_-2.09/  1.05_-1.05/   1.05_3.14/    -2.09_-2.09/  3.14_-1.05/
00_2.09/   1.05_1.05/    -2.09_0.00/   -2.09_2.09/   3.14_1.05/
00_3.14/   -1.05_-2.09/  2.09_0.00/    2.09_-2.09/   3.14_-2.09/
-1.05_0.00/  -1.05_2.09/   -2.09_-1.05/  2.09_2.09/    3.14_2.09/

The data is in the following format with tree columns: \(\tau\), expectation value and error:

!head ALF_data/Hubbard/res/Green_tau/0.00_0.00/dat

0000000       0.03455197       0.00477901
1000000       0.01719004       0.00965419
2000000       0.01193452       0.00737570
3000000       0.01320968       0.00500187
4000000       0.00127847       0.00630168
5000000       0.00038241       0.00930724
6000000       0.00744082       0.00568311
7000000       0.00335035       0.00378054
8000000      -0.00131313       0.00768490
9000000      -0.00099027       0.00230007