Krauth 2010

From Werner KRAUTH

Jump to: navigation, search

W. Krauth, Four lectures on computational statistical physics (Oxford University Press (2010)) (in: Exact Methods in Low-dimensional Statistical Physics and Quantum Computing Lecture Notes of the Les Houches Summer School: Volume 89, July 2008 Edited by Jesper Jacobsen, Stephane Ouvry, Vincent Pasquier, Didina Serban and Leticia Cugliandolo)


Contents

Paper

In my lectures at the Les Houches Summer School 2008, I discussed central concepts of computational statistical physics, which I felt would be accessible to the very cross-cultural audience at the school. I started with a discussion of sampling, which lies at the heart of the Monte Carlo approach. I specially emphasized the concept of perfect sampling, which offers a synthesis of the traditional direct and Markov-chain sampling approaches. The second lecture concerned classical hard-sphere systems, which illuminate the foundations of statistical mechanics, but also illustrate the curious difficulties that beset even the most recent simulations. I then moved on, in the third lecture, to quantum Monte Carlo methods, that underly much of the modern work in bosonic systems. Quantum Monte Carlo is an intricate subject. Yet one can discuss it in simplified settings (the single-particle free propagator, ideal bosons) and write direct-sampling algorithms for the two cases in two or three dozen lines of code only. These idealized algorithms illustrate many of the crucial ideas in the field. The fourth lecture attempted to illustrate aspects of the unity of physics as realized in the Ising model simulations of recent years. More details on what I discussed in Les Houches, and wrote up (and somewhat rearranged) here, can be found in my book, "Statistical Mechanics: Algorithms and Computations", as well as in recent papers. Computer programs are available for download and perusal at the book's web site.

Table of contents

  1. Sampling
    1. Direct sampling, sample transformation
    2. Markov-chain sampling
    3. Perfect sampling
  2. Classical hard-sphere systems
    1. Molecular dynamics
    2. Direct sampling, Markov chain sampling
    3. Cluster algorithms, birth-and-death processes
  3. Quantum Monte Carlo simulations
    1. Density matrices, naive quantum simulations
    2. Direct sampling of a quantum particle: Lévy construction
    3. Ideal bosons: Landsberg recursion and direct sampling
  4. Spin systems: samples and exact solutions
    1. Ising Markov-chains: local moves, cluster moves
    2. Perfect sampling: semi-order and patches
    3. Direct sampling, and the Onsager solution
  5. Appendix (Python program for exact sampling)
  6. References


Electronic version (from arXiv, original version)

Illustration

50px
This is a picture illustrating the concept of perfect sampling, treated in the first chapter of the lectures. The concept relates to attempts for using Markov chains for sampling the exact stationary distribution, rather than to do that only asymptotically.


Algorithm

Here is a 43-line computer program, written in Python, which implements the direct-sampling algorithm for non-interacting Bosons in a three-dimensional trap. Graphic output is more basic than in the below gif output.

##========+=========+=========+=========+=========+=========+=========+
## PROGRAM : direct_harmonic_bosons.py
## PURPOSE : implement SMAC algorithm 4.6 (direct-harmonic-bosons)
##           with plot of two-dimensional snapshot in pylab (SMAC fig. 4.23)
## LANGUAGE: Python 2.5
##========+=========+=========+=========+=========+=========+=========+
from pylab import *
from math import sqrt, sinh, tanh,exp
from random import uniform as ran,gauss
def z(beta,k):
   sum = 1/(1- exp(-k*beta))**3
   return sum
def canonic_recursion(beta,N):
   Z=[1.]
   for M in range(1,N+1):
      Z.append(sum(Z[k] * z(beta,M-k) for k in range(M))/M)
   return Z
def pi_list_make(Z,M):
   pi_list=[0]+[z(beta,k)*Z[M-k]/Z[M]/M for k  in range(1,M+1)]
   pi_sum=[0]
   for k in range(1,M+1):
      pi_sum.append(pi_sum[k-1]+pi_list[k])
   return pi_sum
def tower_sample(data,Upsilon): #naive tower sampling, cf. SMAC Sect. 1.2.3
   for k in range(len(data)):
      if Upsilon<data[k]: break
   return k
def levy_harmonic_path(Del_tau,N): #
   beta=N*Del_tau
   xN=gauss(0.,1./sqrt(2*tanh(beta/2.)))
   x=[xN]
   for k in range(1,N):
      Upsilon_1 = 1./tanh(Del_tau)+1./tanh((N-k)*Del_tau)
      Upsilon_2 = x[k-1]/sinh(Del_tau)+xN/sinh((N-k)*Del_tau)
      x_mean=Upsilon_2/Upsilon_1
      sigma=1./sqrt(Upsilon_1)
      x.append(gauss(x_mean,sigma))
   return x
#-------------------------------------------------------
# -- main program starts here, don't choose N too large
#-------------------------------------------------------
N=512
T_star=.3
T=T_star*N**(1./3.)
beta=1./T
Z=canonic_recursion(beta,N)
M=N
#Dic_cycles={}
x_config=[]
y_config=[]
while M > 0:
   pi_sum=pi_list_make(Z,M)
   Upsilon=ran(0,1)
   k=tower_sample(pi_sum,Upsilon)
   x_config+=levy_harmonic_path(beta,k)
   y_config+=levy_harmonic_path(beta,k)
   M -= k
#  if Dic_cycles.has_key(k): Dic_cycles[k]+=1
#  else: Dic_cycles[k]=1
#ll = Dic_cycles.keys()
#ll.sort()
#   print k, Dic_cycles[k]
##--begin figure output
figure(1)
axis('scaled')
xlim(-5.,5.)
ylim(-5.,5.)
plot(x_config,y_config,'ro')
xlabel('$x$')
ylabel('$y$')
title('3d bosons in trap (2d projection):'+' $N$= '+str(N)+' T* =' + str(T_star))
show()
Personal tools