# Krauth 2010

### From Werner KRAUTH

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**

- Sampling
- Direct sampling, sample transformation
- Markov-chain sampling
- Perfect sampling

- Classical hard-sphere systems
- Molecular dynamics
- Direct sampling, Markov chain sampling
- Cluster algorithms, birth-and-death processes

- Quantum Monte Carlo simulations
- Density matrices, naive quantum simulations
- Direct sampling of a quantum particle: Lévy construction
- Ideal bosons: Landsberg recursion and direct sampling

- Spin systems: samples and exact solutions
- Ising Markov-chains: local moves, cluster moves
- Perfect sampling: semi-order and patches
- Direct sampling, and the Onsager solution

- Appendix (Python program for exact sampling)
- References

Electronic version (from arXiv, original version)

# Illustration

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()

Categories: Publication | 2010 | Algorithm | Hard spheres | Bosons | Perfect sampling