Maximum Entropy

The maximum entropy method (MEM) is needed to continue imaginary-time spectral data (e.g. from quantum monte carlo simulations) to the real axis. Warning: be careful with all spectra obtained from MEM! Never publish without full understanding and careful checks.

Sandvik Code

The program descending from code written by Sandvik is a relatively simple MEM program which does not diagonalize the covariance matrix and uses a full search algorithm. The advantage is its simplicity and the fact that few measurements of G(tau) suffice to use it.

This (local) directory contains several versions of the program itself (multi..., cont...) and helpers to prepare the input (dat...). I am using datband (source: datband.f) for input preparation and multi1 (source: multi1.f, needs cont1.h) for actual MEM computation. In the latter I have included the possibility to use a gaussian (rather than flat) default model, characterized by its variance sigma.

Example: Hubbard model, 1-band, n=1, beta=7, U=4

The usage is demonstrated in an example. From QMC we get the imaginary-time Green function B7dt01U40.gtau, which looks like this: (10 measurements)

The measurements are averaged for each time slice (while neglecting off-diagonal elements of the covariance-matrix) using the program datband (input is expected in the form frequency/G_up/G_down/G_up+G_down):

~nilsb/Depot/MaxEnt-Sandvik/datband < B7dt01U40.gtau && mv g12.dat B7dt01U40.gtau.av

Result is a number of files: gsum.dat, g2M.dat, g2.dat, g1M.dat, g1.dat. The important one (originally g12.dat) has been renamed to B7dt01U40.gtau.av. It is computed by averaging over G_up and G_down seperately (nevertheless symmetric). Therefore even the exactly known values at tau=0,beta obtain error bars. Theoretically one should average over the sum of G_up and G_down since some error of the calculation exactly cancels out in this sum (and this knowlegde should go into MaxEnt); in practice the method doesn't converge in the latter way.

Now the stage is set for the real MaxEnt-part. We use this input file. The parameters are explained in another sample input file, which we print out here:

   1    # momenta
  71    # times
  71    # times in MC sampling
  7. 0  beta, sigma
  200   # frequencies
  0.1   Delta frequency
 6000   # annealing steps
 1200   starting alpha
    1   # smoothed runs
 12345  random number seed

Here, # momenta refers to the number of independent data sets (different momenta in finite-d calculations, different bands, etc.), # times and # times in MC sampling is the number of imaginary time slices taken for MEM and present in the input from QMC, respectively (the first number can be smaller for taking only every nth time slice; note: includes tau=0 and tau=beta, therefore times=L+1), beta is the inverse temperature, sigma (if larger than 0) is the width of the gaussian default model, # frequencies and Delta frequency are number of and distance between frequency points in the result (on each side of the axis), respectively. Intrinsic parameters of the MEM procedure are # annealing steps, starting alpha, # smoothed runs, and random number seed; they should only be changed by experts.

Here the actual run (takes a few minutes on qmc):

ln -s B7dt01U40.gtau.av gav.dat
multi1 < B7dt01.multi.i
mv aw.dat B7dt01U40.aw.dat && rm sums.dat

The spectrum B7dt01U40.aw.dat looks like this (metallic with well-developed Kondo peak, but no pinning to the non-interacting peak value of 1/Pi=0.318 due to the relatively high temperature and, possibly, MEM errors):

Nils Blümer