
% BB Oct 02 2007

\chapter{Markov Chain Monte Carlo \label{MC}}

\markboth{Markov Chain Monte Carlo Simulations and Their Statistical 
Analysis}{Markov Chain Monte Carlo}  %CHANGE HERE

So far we sampled statistically independent events. This requires
a method to generate the desired probability distribution directly.
But for many problems such a method does not exist. It 
may still be possible to generate statistically independent event, 
but only in such a way that the statistically important contributions 
are missed, {\it i.e.}, they occur only with tiny probabilities.
The way out are in many 
cases\index{Monte Carlo: see MC|textbf}\index{simulations|textbf}
\index{MC|textbf}\index{Markov chain Monte Carlo}{Markov chain Monte 
Carlo} (MC) simulations\footnote{We follow the notation of the physics 
literature, where MC is widely used to denote Markov chain Monte 
Carlo simulations. In mathematics and computer science the more 
accurate notation MCMC\index{MCMC} is preferred, which we shall
use occasionally to emphasize the Markov chain aspect. In the
previous chapters we performed MC simulation, which are obviously 
not of the MCMC type.}.
A \index{Markov process}{\bf Markov process} 
generates event $n+1$ stochastically from event 
$n$, without information about event $n-1$ needed at that point. 
A \index{MC algorithms|textbf}{\bf MC algorithm} arranges the 
statistical weights of the Markov process so that the 
desired probability distribution is asymptotically ({\it i.e.}, 
for large $n$) approached. This enables the
\index{importance sampling}\index{sampling!importance sampling}{\bf 
importance sampling}
discussed in this chapter. For statistical physics the problem to 
generate configurations of the Gibbs canonical ensemble with their 
correct 
\index{weights!Boltzmann}\index{Boltzmann weights}{Boltzmann weights}
(also called 
\index{weights!Gibbs}\index{Gibbs weights}{Gibbs weights}) was 
solved in the well-known paper by Nicholas 
\index{Metropolis}{Metropolis}, 
Arianna \index{Rosenbluth}{Rosenbluth}, Marshall Rosenbluth, 
Augusta \index{Teller}{Teller} and Edward Teller~\cite{Me53}. 
A price paid in MC simulations is that subsequent configurations 
become correlated -- possibly strongly. This has consequences for 
the error analysis, which are discussed in chapter~\ref{Markov_Error}.

In this chapter we focus on practical aspects of performing MC
calculations. Some basic knowledge of 
\index{statistical physics|textbf}{\it statistical physics} 
will help the understanding, but is not an absolutely necessary 
prerequisite. Preliminaries together with the $2d$ 
\index{Ising model}{Ising model} are summarized in the next section. 
Good statistical physics textbooks are those by \cite{Ch87}, 
\cite{LaLi51}, \cite{Hu87} and \cite{Pa72}.

Convergence of importance sampling (including the Metropolis method) 
is proven in section~\ref{Importance_sampling}. Additional details and 
the heatbath algorithm are discussed for the $O(3)$ $\sigma$-model, 
which is also known as the Heisenberg ferromagnet. Section~\ref{Potts} 
introduces generalized Potts models in arbitrary dimensions and 
presents Fortran code for MC simulations of these models. Potts 
models, which include the Ising model as a special case, are 
sufficiently involved to exhibit many relevant features of MC
simulations, while they are simple enough to allow us to focus 
on the essentials. MC simulations of {\bf continuous systems}
face a number of additional difficulties. In 
section~\ref{MC_continuous} we illustrate them using 
$O(n)$ $\sigma$-models. The generalization of the concepts of
this chapter to more involved systems is rather straightforward. 
See, for instance, chapter~\ref{LGT} on lattice gauge theory.

\section{Preliminaries and the Ising Model \label{sec_2dI} }

MC simulations of systems described by the 
\index{Gibbs ensemble}\index{canonical ensemble}{\bf Gibbs
canonical ensemble}\index{ensemble!canonical}
\index{ensemble!Gibbs-Boltzmann}aim at calculating estimators of 
physical observables at a\index{temperature}{temperature} $T$. In the 
following we choose units so that the 
\index{Boltzmann constant}{Boltzmann constant}
becomes one and, consequently,
\begin{equation} \label{Temperature}
\beta = 1/T\ .
\end{equation}
Let us consider the numerical calculation of the 
\index{exceptation value!canonical ensemble}
\index{canonical ensemble!expectation value}{\bf expectation 
value} of an \index{canonical ensemble!observable} 
\index{observable, canonical ensemble}{\bf observable} 
$\mathcal{O}$. Mathematically all systems on a computer are 
\index{systems!discrete}{discrete}, because a finite word length 
has to be used. Hence, the expectation value is given by the sum
\begin{equation} \label{O}
  \widehat{\mathcal{O}} = \widehat{\mathcal{O}} (\beta) = 
  \langle \mathcal{O} \rangle  = Z^{-1} \sum_{k=1}^K
  \mathcal{O}^{(k)}\,e^{-\beta\,E^{(k)} }
\end{equation}
where
\begin{equation} \label{Z}
 Z = Z(\beta) = \sum_{k=1}^K e^{-\beta\,E^{(k)} }
\end{equation}
is called \index{partition function|textbf}{\bf partition function}. 
The index $k=1,\dots , K$ labels {\bf all} 
\index{configuration}{\bf configurations} of the system and $E^{(k)}$ 
is the (internal) \index{energy}{energy} of configuration $k$. The 
configurations are also called 
\index{state!microstate}\index{microstate}{\bf microstates}. To 
distinguish the \index{configuration!index}{configuration index}
from other indices, we put it in parenthesis.

A simple model of a \index{ferromagnet!Ising}{ferromagnet} is the 
\index{Ising}\index{Ising model|textbf}Ising model \cite{Is25},
which was introduced by Lenz at Hamburg University, who left its
investigations to his student Ising. The
\index{energy!Ising model}\index{Ising model!energy}{energy} 
of this spin model is given by\index{discrete systems!Ising model}
\index{systems!discrete!Ising model}
\begin{equation} \label{Ising_E}
  E_0 = - \sum_{<ij>} s_i\,s_j\ .
\end{equation}
The sum is over the \index{site}{sites} of a 
hypercubic\footnote{Other lattice geometries are occasionally used 
in the literature.} $d$-dimensional \index{lattice}{lattice} and the 
\index{Ising model!spin}{Ising spins} take the values
\begin{equation} \label{Ising_spin}
  s_i=\pm 1 ~~{\rm for}~~ i=1,\dots ,N ~~{\rm with}~~ 
 N = \prod_{i=1}^d N_i
\end{equation}
where $N_i$ are the edge lengths of the lattice. For the Ising model
the number of configurations is $K=2^N$. Already for a moderately sized 
lattice, say $100\times 100$, this is a tremendously large number, 
namely $2^{10\,000}$ which by far exceeds the estimated number of 
$10^{80}$ protons in the universe. With the exception of very small 
systems it is therefore impossible to carry out the sums in equations
(\ref{O}) and~(\ref{Z}). Instead, the large number suggests a 
statistical evaluation through generating a representative subsample. 

On an infinite lattice the \index{Ising model!$2d$}{$2d$ Ising model}
was solved by\index{Onsager} \cite{On44}. See, for instance, 
\cite{Ch87} or \cite{Hu87} for textbook treatments. 
A\index{second oder transition:!see phase transition!second order}
\index{phase transition!second order|textbf}{second oder phase 
transition} is found and exact 
solutions are known for many observables. This makes the $2d$ Ising 
model an ideal testing ground for simulation methods as well as for
perturbative expansions. Of course, the availability of an exact 
solution is related to the simplicity of the model in $2d$. In $3d$ 
the solution is not known, what is the typical case for realistic 
models of physical processes. The $1d$ 
\index{Ising model!$1d$}Ising model is still far simpler than the 
$2d$ model, and to derive its exact solution~\cite{Is25} is often 
given as an exercise to graduate students. The $1d$ model does not 
exhibit a phase transition at a finite temperature, although some 
interesting physics is found for $\beta\to\infty$.

Among the prominent features of the $2d$ model are the exact 
equation for\index{Ising model!$2d$!critical temperature} 
its\index{critical temperature|textbf} 
\index{temperature!critical|textbf}{\bf critical temperature}, 
also called\index{Curie temperature|textbf}
\index{temperature!Curie|textbf}{\bf Curie temperature},
\begin{equation} \label{2d_Ising_bc}
  \beta_c\ =\ {1\over T_c}\ =\ {1\over 2}\, \ln (1 + \sqrt{2} )
  = 0.4406867935\dots 
\end{equation}
and exact results 
for\index{critical exponents|textbf}\index{exponents, critical}
\index{exponents, critical!$\alpha,\,\beta,\,\gamma,\,\delta,\,
\nu,\, \eta$}{\bf critical exponents}, 
which govern the singularities at the phase transition:
\index{Ising model!critical exponents}
\begin{equation} \label{2d_Ising_exponents}
  \alpha = 0\,, ~~\beta = {1\over 8}\,, ~~\gamma = {7\over 4}\,,~~
  ~~\nu = 1\,, ~~\delta = 15 ~~{\rm and}~~ \eta = {1\over 4}\,.
\end{equation}
Using the\index{reduced temperature}\index{temperature!reduced} 
\index{reduced temperature}{\bf reduced temperature} $t$ defined by
\begin{equation} \label{reduced_temperature}
 t\ =\ {T-T_c\over T_c}\ =\ \frac{\beta_c-\beta}{\beta}
\end{equation}
the (positive) critical exponents are given by the behavior of 
physical observables near the critical temperature:
\begin{eqnarray} \label{def1_exponents}
  \index{critical exponents!specific heat}\index{specific heat} 
  \index{critical exponents!magnetization}\index{magnetization|textbf} 
  \index{critical exponents!susceptibility}\index{susceptibility}
  {\rm the\ specific\ heat}~~ C &\sim& |t|^{-\alpha}\,,\\
  {\rm the\ magnetization}~~  M &\sim& |t|^{\beta}\,,~~t<0\,,\\
  {\rm the\ susceptibility}~~ \chi &\sim& |t|^{\-\gamma}\,,\\
  {\rm the\ correlation\ length}~~ \xi &\sim& |t|^{-\nu}\,. 
  \label{def3_exponents} \index{correlation length|textbf}
  \index{critical exponents!correlation length}
\end{eqnarray}
The equations can be written this way, because it can be shown that the 
Ising model exponents agree in the $t>0$ high and $t<0$ low temperature 
phase, with exception of $M$ for which on has to have $t<0$. Two more 
exponents are defined for $t=0$: The magnetization as function of an 
external magnetic field
\begin{equation} \label{def2_exponents}
   M \sim H^{-1/\delta}\,,\ (t=0)\,, 
\end{equation}
and $t=0$ power law correlations
\begin{equation} \label{def4_exponents}
  \left< s_i s_j\right> \sim |\vec{x}_i-\vec{x}_j|^{-d+2-\eta}\ 
  {\rm for}\ |\vec{x}_i-\vec{x}_j| \to \infty\,,\ (t=0)\,. 
\end{equation}
See the review by \cite{PeVi02} for a detailed exposition of the 
subject, including a discussion of universal amplitudes of the 
encountered singularities. In MC simulations the evaluation of 
critical exponents is possible by many methods. Particularly useful 
are \index{finite size scaling}\index{FSS}{\bf Finite Size Scaling 
(FSS)} extrapolations, relying on the FSS theory originally introduced 
by\index{Fisher, M.E.} \cite{Fi71}. For the $3d$ Ising model exact
critical exponents are not known. Estimates from p.604 of \cite{PeVi02} 
\index{Ising model!critical exponents}are 
\begin{eqnarray} \label{3d_Ising_exponents}
  \alpha=0.110\,(1)\,,~~ & \beta=0.3265\,(3)\,,&~~ \gamma=1.2372\,(5)\,,\\
  \nonumber \delta=.4.789\,(2)\,,~~ & \nu=0.6301\,(4) & 
            ~~{\rm and}~~~\eta=0.0364\,(4)\,.
\end{eqnarray}
As always, error bars are given in parenthesis.

It appears that the \index{partition function!analyticity}partition 
function~(\ref{Z}), being a sum of analytic terms, must be an analytic 
function of the temperature. Indeed, this is the case for a finite 
volume. However, one is interested in the infinite volume limit, where 
the sum contains an infinite number of terms and analyticity can be 
lost. When analyticity is lost, this must be due to some subtlety in 
taking the infinite volume limit, expected to be associated with 
the development of long-range correlations between spins (or other 
dynamical variables). A powerful method to identify non-analytic points 
is to construct \index{order parameters}{\bf order parameters}. Usually 
an order parameter characterizes a symmetry of the system, so that the 
order parameter is zero (or constant) in one phase and not in the 
other.  For the Ising ferromagnet this is the magnetization, 
which vanishes in a state with reflection 
symmetry $s_i\to-s_i$.  If one has an order parameter of this 
kind one can be sure that there is a singularity in passing from 
one phase to the other, because an analytic function cannot vanish 
in a finite range without vanishing identically. On the other hand,
one finds phase transitions without order parameters.

If the transition is temperature driven and
\index{phase transition!second order}{second order}, 
the (infinite volume) \index{correlation length}{correlation length} 
goes, as function of the temperature, continuously to infinity 
at the critical point. Sufficiently close to the transition point,
the correlation length will greatly exceed any microphysical 
length associated with the problem ({\it e.g.}, the lattice
spacing). These observations lead to the concepts of
\index{scale invariance}{\bf scale invariance} and 
\index{universality|textbf}{\bf universality}, which form the
cornerstone of the modern theory of phase transitions.
Scale invariance follows from the idea that non-trivial correlations
on length scales much larger than the microscopic one can only
depend on the ratio of the length in question to the diverging
correlation length. There should be just one fundamental diverging
length, since there is only one parameter, which must be tuned to
reach the critical point, {\it i.e}, all divergent quantities 
should be given as functions of the 
\index{reduced temperature}{reduced temperature} 
(\ref{reduced_temperature}). Universality
follows from the idea that likely only modes with long-wavelength
fluctuations in equilibrium are important which are required by
the symmetries of the order parameter. Together these principles
tell us that for the investigation of singularities of thermodynamic 
functions near second order phase transitions, it is sufficient
to look for the simplest scale-invariant model in the appropriate
dimension with the required symmetry. In this way symmetries 
define universality classes, which share the same critical 
exponents. 

The order parameter concept works also for predicting that a 
particular \index{phase transition!first order}{phase transition}
is first order. In contrast to the second order case one has no 
powerful predictions for the physical quantities, like the latent
heat, of a first order transition.
 
\begin{figure}[tb]
\vspace{5pc}
\centerline{\hbox{ \psfig{figure=figures/2dI_es.eps,width=10cm} }}
\caption{Expectation value of the energy per spin $\widehat{e}_{0s}$ 
versus $\beta$ for the $2d$ Ising model on an $20\times 20$ lattice. 
The full line is the exact result 
of\index{Ising model!energy!per spin}\index{Ferdinand and Fisher}
\index{Ising model!$2d$}[Ferdinand and Fisher (1969)]. 
The data points are from a multicanonical simulation, see 
assignment~\ref{A_MUCA3} in chapter~\ref{assign_MUCA}. This 
figure was first published in [Berg (2003)]. \label{fig_2dI_es} }
\end{figure}
 
For the purposes of testing our first MC algorithms the paper 
by\index{Ising model!$2d$}\index{Ferdinand and Fisher} \cite{FeFi69}, 
build on work by \cite{Ka49}, is most useful. Ferdinand and Fisher
 derive $2d$ Ising model results for $m\times n$ lattices, which 
allow for direct comparisons with stochastic results. The program
\begin{equation} \label{ferdinand}
\index{ferdinand.f}{\tt ferdinand.f}
\end{equation}
in the folder {\tt ForProg/Ferdinand} is based on their paper. It 
calculates important physical observables and prepares their graphical 
presentation. To give a example, figure~\ref{fig_2dI_es} compares for 
a $20\times 20$ lattice the exact 
\index{expectation value!energy per spin}{\bf expectation value 
of the energy per spin}
\begin{equation} \label{expect_e0s}
 \widehat{e}_{0s} (\beta) = \widehat{E_0} (\beta) / N
\end{equation}
from {\tt ferdinand.f} with its stochastic estimates of a multicanonical 
simulation (compare assignment~\ref{A_MUCA3} of chapter~\ref{assign_MUCA}).

\subsection{Lattice labeling \label{lattice_labeling} }
\index{lattice!labeling|textbf}

\begin{figure}[tb]
\vspace{5pc}
\centerline{\hbox{ \psfig{figure=figures/lat.eps,width=12cm} }}
\caption{An $8\times 8$ lattice in the Cartesian coordinates.
\label{fig_labeling} }
\end{figure}

Let us consider lattices with periodic boundary conditions. It is 
convenient to label all sites of the lattice by a single integer 
variable, {\tt is}, which we call 
\index{site!number}{\bf site number}. The 
dimension of the lattice is given by another integer variable, 
{\tt nd}. The Cartesian coordinates of a site are chosen to be
\begin{equation} \label{site_coordinate}
 x^i = 0, \dots ,n^i-1 ~~~{\rm for}~~~ i=1, \dots , d.
\end{equation}
The lattice size dimensions $n^i$ are stored in the array 
{\tt nla(nd)} (we allow an asymmetric lattice).  Let the 
coordinates of one site be given by an integer array
{\tt ix(nd)}. The Fortran function 
\begin{equation} \label{isfun}
  \index{isfun.f}{\tt isfun(ix,nla,nd)}
\end{equation}
converts the coordinates of the given site to a site number
according to the formula
\begin{equation} \label{is}
  {\tt is} = 1 + \sum_{\tt id=1}^{\tt nd} {\tt ix(id)}\times
  {\tt nsa(id)} ~~{\rm with}~~ {\tt nsa(id)} = 
\begin{cases} 1 ~~{\rm for}~~ {\tt id=1}, \\
  \prod_{\tt i=1}^{\tt id-1} {\tt nla(i)} ~~{\rm for}~~{\tt id} > 1.
\end{cases}
\end{equation}
See figure~\ref{fig_labeling} for {\tt nd=2}. Once we have set
the values for {\tt nla(i)} we obtain the number of sites of the
lattice, {\tt ns}, from the function 
\begin{equation} \label{nsfun}
\index{nsfun.f}
  {\tt nsfun(nla,nd)} = \prod_{\tt i=1}^{\tt nd} {\tt nla(i)}\ .
\end{equation}
To get the coordinates for a given site number {\tt is} is achieved
by the subroutine 
\begin{equation} \label{ixcor}
  \index{ixcor.f}{\tt ixcor(ix,nla,is,nd) }
\end{equation}
which returns in the array {\tt ix(nd)} the Cartesian 
coordinates corresponding to the site number ${\tt is}$. Used is 
an iteration procedure relying on Fortran integer division
({\it i.e.}, $1=[5/3]=[5/5]$, $0=[4/5]$, etc.):
$$ \tt{ix(nd)} = \left[{\tt (is-1) \over [ns\,/\,nla(nd)]}\right] $$
and for ${\tt id} = {\tt nd} - 1,\dots ,1$
$$ {\tt ix} ({\tt id}) = [ ({\tt js}({\tt id}) - 1)\, /\,
   {\tt nspart} ({\tt id}) ]\, ,$$
where
$$ {\tt js} ({\tt id}) = {\tt is} -
  \sum_{\tt i=nd}^{\tt id+1}\, {\tt ix}({\tt i})\,
  {\tt nspart} ({\tt i}) ~~{\rm and}~~  {\tt nspart}({\tt id}) =
  {\tt ns}\prod_{\tt i=nd}^{\tt id+1} (\tt nla(i))^{-1}\ .$$
The subroutine 
\begin{equation} \label{ipointer}
\index{ipointer.f}{\tt ipointer(ipf,ipb,nla,ix,ns,is,nd) } 
\end{equation}
calculates the site numbers of the neighbors of a given site
{\tt is}. Input are {\tt nla, ix, ns, is}\ and\ {\tt nd}. On output 
the \index{pointer arrays|textbf}{\bf pointer arrays}
$$ {\tt ipf(id,is)} ~~{\rm and}~~ {\tt ipb(id,is)} $$
are returned for ${\tt id} = 1,\dots, {\tt nd}$, where
{\tt ipf} gives the site numbers of the {\tt nd} neighbors in
forward direction and ${\tt ipb}$ the site numbers of the {\tt nd}
neighbors in backward direction. To achieve this goal the subroutine
{\tt ixcor}~(\ref{ixcor})  is first used to calculate the Cartesian
coordinates for the site number {\tt is}. From it the coordinates of
the neighbors are calculated, taking into account the periodic boundary
conditions:
\begin{enumerate}
\item For the forward direction ({\tt ipf}) 
 $${\tt ix(id) \to mod (ix(id)+1, nla(id)), \quad id=1,...,nd.}$$
\item For the backward direction ({\tt ipb})
$${\tt ix(id)\to mod (ix(id)-1+nla(id),nla(id)),\quad id=1,...,nd.}$$
\end{enumerate}
Then, the function {\tt isfun}~(\ref{isfun}) is used to get the site
numbers, which are stored in the arrays {\tt ipd} and {\tt ipb}.  For 
a given lattice the pointer arrays {\tt ipf} and {\tt ipb} are 
initialized by a call to the subroutine
\begin{equation} \label{lat_init}
  {\tt lat\_ init(ns,nd,ipf,ipb,nla,ix)}\ .
\end{equation}
The routine needs {\tt ns, nd} and {\tt nla} assigned on input and 
{\tt ix} as workspace. It uses the routines (\ref{ipointer}), 
(\ref{isfun}), (\ref{nsfun}) and (\ref{ixcor}). 

\begin{table}[tbl]  
\caption{Bookkeeping for a $2^24$ lattice: The site number {\tt is}, 
its coordinates {\tt ix}, the forward pointers {\tt ipf} and the 
backward pointers {\tt ipb} are given. These results are obtained 
with the program ${\tt lat3d.f}$ of the directory {\tt ForProg/Lattice}, 
\index{lattice!labeling}see assignment~\ref{assign_pointers}.}
{\small
\centering
\begin{verbatim}
  is   ix(1) (2) (3)   ipf(is,1) (,2) (,3)   ipb(is,1) (,2) (,3)
   1      0   0   0           2    3    5           2    3   13
   2      1   0   0           1    4    6           1    4   14
   3      0   1   0           4    1    7           4    1   15
   4      1   1   0           3    2    8           3    2   16
   5      0   0   1           6    7    9           6    7    1
   6      1   0   1           5    8   10           5    8    2
   7      0   1   1           8    5   11           8    5    3
   8      1   1   1           7    6   12           7    6    4
   9      0   0   2          10   11   13          10   11    5
  10      1   0   2           9   12   14           9   12    6
  11      0   1   2          12    9   15          12    9    7
  12      1   1   2          11   10   16          11   10    8
  13      0   0   3          14   15    1          14   15    9
  14      1   0   3          13   16    2          13   16   10
  15      0   1   3          16   13    3          16   13   11
  16      1   1   3          15   14    4          15   14   12
\end{verbatim}}
\label{tab_pointers}
\end{table}

For $d=3$ the program 
\begin{equation} \label{lat3d} \index{lat3d.f}
  {\tt lat3d.f} ~~{\rm of}~~ {\tt ForProg/Lattice}
\end{equation}
illustrates the use of \index{lat\_init.f}{\tt lat\_init.f}. 
Table~\ref{tab_pointers} collects the thus obtained site numbers, their 
coordinates and their pointers for a $2^24$ lattice. That the forward 
and backward pointers agree in some instances is an artifact of the 
small lattice size, see assignment~\ref{assign_pointers} for more 
cases. For use in MC simulations the results are stored in the common 
block {\tt common\_lat.f} 
\begin{equation} \label{common_lat.f}
   {\tt   common~/lat/~ns,nla(nd),ipf(nd,ms),ipb(nd,ms)}
\end{equation}
which is transferred into the programs through 
$$ {\tt include\ '../../ForLib/common\_lat.f'}$$
statements, {\it i.e.}, the file {\tt common\_lat.f} is kept in the 
directory {\tt ForLib}. In this common block {\tt ns} is the total 
number of sites, {\tt nla} stores the lattice size dimensions, {\tt ipf} 
and {\tt ipb} are the forward and backward pointer arrays. The 
dimensional arguments needed in this common block are assigned in 
the parameter statement {\tt lat.par}\index{\tt lat.par} shown in
table~\ref{tab_latpar}.
%
\begin{table}[th]  
\caption{Parameter file {\tt lat.par}.}
{\small \centering \begin{verbatim}
      character cfln*5,cfile*13 ! Filename (prefix, extension)
      character(cfln="potts")
      parameter(nd=3,ms=20**nd,mlink=nd*ms) 
      parameter(n2d=2*nd,np=ms*(nd*(nd-1))/2) ! Adjust also lat.dat.
\end{verbatim}} \label{tab_latpar}
\end{table}
%
The first row shows character variables, which will be used to name 
data files. Their size should be adjusted to fit the names used. In
the second row {\tt cfln} is set to ``potts''. Other set-ups, like
{\tt cfln*2}, {\tt clnfl="U1"}, will be introduced later. The
variable {\tt cfile} is used to name files in long simulations,
which require backups, see chapter~\ref{LSP}.

In the third row {\tt nd} is the dimension of the 
system, {\tt ms} the number of sites of the lattice, and {\tt mlink} 
the number links. The parameters {\tt nd} and {\tt ms} have to be 
chosen by the user. The choice in table~\ref{tab_latpar} implies 
a $3d$ lattices with a total number of sites $20^3$. The statement 
{\tt lat.par} is transferred into the programs through 
{\tt include 'lat.par'} statements, {\it i.e.}, the file {\tt lat.par} 
is kept in the folder where the program {\tt lat3d.f} is 
compiled. The rationale for keeping {\tt lat.par} and {\tt common\_lat.f} in 
different directories is: We do not intend to change {\tt common\_lat.f}, but 
if such a change is unavoidable it should propagate to all routines 
which include {\tt common\_lat.f}. On the other hand, we expect frequent 
changes of the lattice parameters of {\tt lat.par}, and the actual 
values are chosen for a particular run of the program.

Besides the parameters in {\tt lat.par} 
\index{lattice!size array nla}{\bf the array {\tt nla} needs to be 
initialized}, which contains the {\tt nd} lattice lengths so that
$$\prod_{\tt id=1}^{\tt nd} {\tt nla(id)}\ =\ {\tt ms}$$
holds. This is done by including a file\index{\tt lat.dat} 
{\tt lat.dat}. The purpose of the {\tt lat.dat} complication is 
to allow for asymmetric lattices. An example of {\tt lat.dat}, 
consistent with the choice in table~\ref{tab_latpar} is
\begin{equation} \label{lat.dat}
    {\tt  data~nla/10,20,40/ }
\end{equation}
which is permissible because of $10\times 20\times 40 = 20^3 = 8\,000$. 
The most common choice in our examples are symmetric lattices, 
{\it i.e.}, ${\tt  data~nla/3*20/}$ for {\tt lat.par} of 
table~\ref{tab_latpar}. The statement {\tt lat.dat} is transferred 
into the program in the same way as {\tt lat.par}, {\it i.e.}, 
{\tt lat.dat} is kept in the folder where the program is compiled.

The routines (\ref{isfun}), {\ref{nsfun}), (\ref{ixcor}),
(\ref{ipointer}) and (\ref{lat_init}) of this section are part of 
{\tt ForLib}. We are now ready to turn to simulations of statistical 
physics models.

\subsection{Sampling and Reweighting\label{sec_rew}}

Having established our notation for labeling the lattice points and 
defined the pointers to the nearest neighbors, we are ready to 
calculate the energy of an Ising configuration. Instead of relying 
on the Ising energy function~(\ref{Ising_E}) we use the 
\index{action variable}{\bf action variable}\footnote{The name of the 
variable comes from field theoretic versions of statistical mechanics 
models, which use an action instead of an energy.}
\begin{equation} \label{iact}
  {\tt iact}^{(k)}\ =\ \sum_{<ij>} \delta (q_i^{(k)},q_j^{(k)})\ ,
\end{equation}
which for $q_i=0,\,1,\dots,q-1$ generalizes to the Potts models
introduced in section~\ref{Potts}. The negative {\bf energy per 
link}\index{energy per link} is given by ${\tt 2*iact/mlink}$. For 
$q=2$ the action variable is equivalent to the configurational energy 
$E_0^{(k)}$ (\ref{Ising_E}) of the Ising model due to the equation
\begin{equation} \label{Ising_iact}
  {\tt iact} = -(E_0-d\,N)/2\ .
\end{equation}
We use the Potts notation, because a computer program for the Potts
model is not more involved than a program for the Ising mode, but 
more general due to the variable number of $q$ spin values. The 
variable {\tt iact} is convenient, because it takes on integers in 
the range
\begin{equation} \label{range_iact}
  {\tt iact} \in [0,{\tt mlink}]
\end{equation}
where {\tt mlink} is the number of links of the lattice. The smallest 
possible stepsize is~1, while values larger than 1 are also possible. 
Irregularities are encountered of large ${\tt iact}$. The largest 
{\tt iact} is always {\tt mlink} for $\delta(q_i^{(k)},q_j^{(k)})=1$ 
at all links. The next smaller one is ${\tt mlink}-{\tt nd}$ when the 
spin at one site flips. For $q=2$ the smallest stepsize of {\tt iact} 
is~2: Due to the limited number of spin values a flip of one spin 
changes then the contribution to {\tt iact} (\ref{iact}) always on
an even number of links. For $q\ge 3$ changes of an odd number of 
links become possible. 

\begin{table}[tbl]  
\caption{Subroutine {\tt potts\_act.f} of {\tt ForLib}, 
         which calculates the action variable (\ref{iact}).}  
{\centering\small \begin{verbatim}
      subroutine potts_act(ista,ipf,idel,ns,nqm1,nd,iact)
C Bernd Berg, Sep 5, 2000. Calculates the action variable iact.
      dimension ipf(nd,ns),ista(ns),idel(0:nqm1,0:nqm1)
      iact=0
      do is=1,ns
        ista1=ista(is)
        do id=1,nd
          ista2=ista(ipf(id,is))
          iact=iact+idel(ista2,ista1)
        end do
      end do
      return
      end
\end{verbatim} } \label{tab_potts_act}
\end{table}

For a given configuration the numerical calculation of {\tt iact} 
is achieved by the subroutine\index{potts\_act.f}
\begin{equation} \label{potts_act}
  {\tt potts\_ act.f}~~~{\rm of\ ForLib}
\end{equation}
shown in table~\ref{potts_act}. Note that only forward direction 
pointers are needed. 

\begin{figure}[tb]
\vspace{5pc}
\centerline{\hbox{ \psfig{figure=figures/2dI_h_es.eps,width=10cm} }}
\caption{Energy histograms of 100$\,$000 entries each for the $2d$
\index{Ising model!energy!histograms} 
\index{histograms!$2d$ Ising energy}{Ising model}
on an $20\times 20$ lattice are depicted. By Random
Sampling we obtain the histogram of statistically independent
configurations at $\beta=0$. Histograms at $\beta=0.2$ and $\beta=0.4$ 
are generated by Markov chain MC. Reweighting of the $\beta=0$ random 
configurations to $\beta=0.2$ is shown to fail. See 
assignment~\ref{assign_b0} of the present section and 
assignment~\ref{assign_2dI_H_es} of section~\ref{assign_MC}.  
\label{fig_2dI_H_es} }
\end{figure}

Next, it is straightforward to {\bf sample statistically independent 
configurations}. This holds not only for Potts model but for most 
systems of statistical physics. For the Ising model one simply has 
to generate $N$ spins, each either up or down with 50\% likelihood 
(in the Potts model notation 0 or 1, each with 50\% likelihood). This 
is called \index{random sampling}\index{sampling!random sampling}{\bf 
random sampling}. On its left-hand side figure~\ref{fig_2dI_H_es} 
depicts a thus obtained histogram of 100$\,$000 entries for the $2d$ 
Ising model\index{Ising model!energy!per spin}
\index{energy!per spin}{\bf energy per spin}
\begin{equation} \label{e0s}
  e_{0s} (\beta) = E_0 (\beta) / N\ .
\end{equation}
One has to distinguish the energy measurements on single configurations 
(\ref{e0s}) from their expectation value (\ref{expect_e0s}). The 
expectation value $\widehat{e}_{0s}$ is a single number, while $e_{0s}$ 
fluctuates. From the measurement of many $e_{0s}$ values one finds 
estimators of its moments~(\ref{sample_mom}), in particular the mean 
for which the estimator is denoted by $\overline{e}_{0s}$. In contrast 
to the expectation value $\widehat{e}_{0s}$ itself, its estimator 
$\overline{e}_{0s}$ fluctuates. The histogram from random sampling 
corresponds to $\beta=0$ in the partition function~(\ref{Z}). Also 
shown in figure~\ref{fig_2dI_H_es} are histograms of 100$\,$000 entries 
each at $\beta=0.2$ and $\beta=0.4$, obtained with the Metropolis 
sampling procedure which is a main subject of this chapter. The rather 
small-sized $20\times 20$ system is used to make the results easily 
reproducible.

\begin{table}[tbl]  
\caption{Subroutine {\tt potts\_rwght.f} of {\tt ForLib}. It 
\index{reweighting}{reweights} histograms of the variable 
{\tt iact} (\ref{iact}) from $\beta$ to $\beta_0$.}  
{\centering\small \begin{verbatim}
      subroutine potts_rwght(nlink,ha,hae,beta,beta0,hb,hbe)
C Bernd Berg Dec 12 2000. Reweighting beta -> beta0.
C Reweighting of the histogram ha() and its error bars hae() at beta 
C             to the histogram hb() and its error bars hbe() at beta0.
      include 'implicit.sta'
      include 'constants.par'
      dimension ha(0:nlink),hae(0:nlink),hb(0:nlink),hbe(0:nlink)
c
      hasum=zero
      act0m=zero
      do iact=0,nlink
        hasum=hasum+ha(iact)
        act0m=act0m+iact*ha(iact)
      end do
      act0m=act0m/hasum ! Used to avoid large arguments of exp().
c     
      hbsum=zero
      do iact=0,nlink
        hb(iact)=zero
        hbe(iact)=zero
        if(ha(iact).gt.em12) then
          a=(beta0-beta)*2*(iact-act0m)
          hb(iact)=exp(log(ha(iact))+a)
          hbe(iact)=exp(log(hae(iact))+a)
          hbsum=hbsum+hb(iact)
        end if ! Else hb(iact) is, as ha(iact), zero.
      end do
      factor=hasum/hbsum
      do iact=0,nlink
        hb(iact)=factor*hb(iact)
        hbe(iact)=factor*hbe(iact)
      end do
c     
      RETURN
      END
\end{verbatim} }
\label{tab_potts_rwght}
\end{table}

The histogram entries at $\beta=0$ can be reweighted so that they 
correspond to other $\beta$ values. We simply have to multiply the 
entry corresponding to energy $E$ by $c_{\beta}\exp (-\beta E)$ as 
follows from~(\ref{Z}). The constant $c_{\beta}$ is adjusted to keep 
the normalization unchanged ({\it i.e.}, the sum over the reweighted 
histogram values has to agree with the sum over the original histogram 
values). Similarly histograms corresponding to the Gibbs ensemble at 
some value $\beta$ can be reweighted to a new $\beta_0$ value. This 
is implemented by our subroutine\index{potts\_rwght.f}
\begin{equation} \label{potts_rwght}
  {\tt potts\_ rwght.f}~~~{\rm of}~~~{\tt ForLib}
\end{equation}
which is reproduced in table~\ref{tab_potts_rwght}. 
Care has to be taken to ensure that the involved arguments of the 
exponential function do not become too large. This is done in 
{\tt potts\_rwght.f} by first calculating the mean action, 
{\tt act0m}, and reweighting with respect to the difference 
from the mean. 

The usefulness of the reweighting for FSS investigations has been 
stressed by \cite{FeSw88}.\index{reweighting|textbf} 
\index{Ferrenberg-Swendsen reweighting}{\bf Ferrenberg-Swendsen 
reweighting} is an important technical tool when one likes to 
estimate critical exponents, like those given in equation
(\ref{2d_Ising_exponents}), numerically. The first use of such 
histogram reweighting procedures can be traced back to a paper 
by \cite{Sa59}. \cite{McSi67} were the first to use these technique 
to evaluate physical quantities over a range of temperatures.
Afterwards the method appeared to be dormant and was rediscovered
by \cite{FaMa82} in a context of calculating complex zerors of the
partition function. A time-series reweighting program, 
{\tt rwght\_ts\_zln.f}, which prepares jackknife error analysis, 
can be found in appendix~\ref{More_Fortran}.

In figure~\ref{fig_2dI_H_es} reweighting is done from $\beta_0=0$ to
$\beta=0.2$. By comparison to the histogram from a Metropolis 
MC calculation at $\beta=0.2$, the result is seen to be disastrous. 
The reason is easily identified: In the range where the $\beta=0.2$ 
histogram takes on its maximum, the $\beta=0$ histogram  has not a 
single entry, {\it i.e.}, our naive sampling procedure misses the 
{\bf important configurations} at $\beta =0.2$. As stressed in
\cite{FeSw88}, reweighting to new $\beta$ values works only
in a range $\beta_0\pm \triangle \beta$, where $\triangle\beta\to 0$
in the infinite volume limit (details are given below). 
These are the $\beta$ values for which the energies measured in the
$\beta_0$ simulation still cover the expectation values 
$\widehat{e}_{0s}(\beta)$. For the $\beta_0=0$ histogram of 
100$\,$000 entries a reasonably large statistic is only available 
down to $e_{0s}=-0.2$. According to the $\widehat{e}_{0s}(\beta)$ 
curve depicted in figure~\ref{fig_2dI_es} this would allow to cover 
the range up to approximately $\beta=0.1$.

In particular it is clear that a groundstate configuration, which is 
the physically relevant configuration at a sufficiently low 
temperature, may never be obtained by random sampling. {\it E.g.},
for the Ising model there are two groundstates (all spins up or all 
spins down) in altogether $2^N$ configurations. For large $N$ they 
will never be found when each configuration is sampled with
probability $2^{-N}$.

\subsubsection{Important configurations and canonical reweighting 
               range}

\index{configuration!important}Let us determine the important 
contributions to the partition function~(\ref{Z}). The partition 
function can be re-written as a sum over energies
\begin{equation} \label{Z_sum_E}
  Z = Z(\beta) = \sum_E n(E)\, e^{-\beta\,E}
\end{equation}
where $n(E)$ is defined as the \index{number of states}{\bf  number of 
microstates} $k$ with energy $E$. The \index{spectral density|textbf}{
\bf spectral density} $g(E)$ normalizes this quantity to one:
\begin{equation} \label{spectral_density}
  g(E) = \frac{n(E)}{\sum_E n(E)}\ .
\end{equation}
For a fixed value of $\beta$, away from first and second order phase 
transitions, the canonical energy probability density
\begin{equation} \label{P_E}
  P_{\beta} (E) = c_{\beta}\, n(E)\, e^{-\beta E} 
\end{equation}
is peaked around the average value $\widehat{E}(\beta)$ with a width 
proportional to $\sqrt{V}$, where $V$ is the volume of the system. The 
normalization constant $c_{\beta}$ ensures that $\sum_E P_{\beta}(E)=1$ 
holds. 

That the width of the energy distribution is $\triangle E\sim\sqrt{V}$ 
follows from the fact that, away from phase transition points, the
fluctuations of the $N\sim V$ lattice spins are only locally correlated,
so that the magnitude of a typical fluctuations is $\sim \sqrt{N}$.
This proportionality is the same as the one encountered for the sum 
of uniformly distributed $\pm 1$ random numbers. Therefore, the 
reweighting range is $\triangle\beta\sim 1/\sqrt{V}$, so that 
$\triangle \beta E\sim \sqrt{V}$ stays within the fluctuations 
of the system (the energy $E$ is an extensive quantity $\sim V$).

The reweighting range increases at second order phase transition
\index{phase transition!second order}point, because critical 
fluctuations are larger than non-critical fluctuations. One has 
$\triangle E \sim V^x$ with $1/2<x<1$ and the requirement 
$\triangle\beta E\sim V^x$ yields $\triangle\beta \sim V^{x-1}$. 
For \index{phase transition!first order}first order phase transitions 
one has a latent heat $\triangle V\sim V$, but this does not mean 
that the reweighting range becomes of order one. In essence, the 
fluctuations collapse, because the two phases become separated by
an interfacial tension. One is back to fluctuations within either
of the two phases, {\it i.e.} $\triangle\beta\sim 1/\sqrt{V}$.

The {\bf important configurations} at temperature $T=1/\beta$ are at 
the energy values for which the probability density $P_{\beta}(E)$ 
is large. To sample them efficiently, one needs a procedure which 
generates them with their Boltzmann weights
\begin{equation} \label{w_B}
  w_B^{(k)} = e^{-\beta E^{(k)}}\ . 
\end{equation}
The number of configurations $n(E)$ and the weights combine then 
so that the probability to generate a configuration at energy $E$
becomes precisely $P_{\beta}(E)$ as given by~(\ref{P_E}). In 
section~\ref{Importance_sampling} we introduce a Markov process 
which generates configurations with the desired weights. Before 
we illustrate reweighting and the importance of jackknife error
analysis for a system which is small enough so that reweighting
of random sampling works down to low temperatures.

\subsection{Assignments for section~\ref{sec_2dI} 
\label{assign_sampling}}

\begin{enumerate}

{ \item Use the program {\tt lat3d.f} of {\tt ForProg/Lattice} to 
reproduce the results of table~\ref{tab_pointers}. Next, produce the 
same table for $2^3$ and $3^3$ lattices. Compile and run also the 
program \index{lat4d.f}{\tt lat4d.f} of {\tt ForProg/lat} and interpret 
the output. \label{assign_pointers} }

{ \item Reproduce the $\beta=0$ histogram of figure~\ref{fig_2dI_H_es} 
and its associated error bars.  Subsequently, reweight the histogram
(including error bars) to $\beta_0=0.2$ and $\beta_0=0.05$. Note: The 
histogram relies on 100$\,$000 independent data points, each generated 
by a call of ${\tt potts\_ ran.f}$ (default seeds are used for the 
random number generator).  The error bars are found from the variance 
of the bimodal distribution~(\ref{bimodal_variance}).\label{assign_b0}} 

\end{enumerate}


\section{Importance Sampling and Markov Chain Monte Carlo
         \label{Importance_sampling} }

For the canonical ensemble \index{importance sampling|textbf} 
\index{sampling!importance sampling|textbf}{\bf Importance sampling} 
generates configurations $k$
with probability
\begin{equation} \label{P_B}
   P_B^{(k)} = c_B\, w^{(k)}_B = c_B\, e^{-\beta E^{(k)}}
\end{equation}
where the constant $c_B$ is determined by the normalization condition
$\sum_k P_B^{(k)} = 1$. The \index{state!vector}{\bf state vector} 
$(P_B^{(k)})$, for which the configurations are the vector indices, 
is called 
\index{Boltzmann state}\index{state!Boltzmann}{\bf Boltzmann state}.  
When configurations are {\bf stochastically} generated with probability 
$P_B^{(k)}$, the \index{expectation value!arithmetic average}{\bf 
expectation value}~(\ref{O}) becomes the {\bf arithmetic average:}
\begin{equation} \label{O_B}
  {\widehat {\cal O}} = {\widehat {\cal O}} (\beta) =
  \langle {\cal O} \rangle  = \lim_{N_K\to \infty} 
  {1\over N_K} \sum_{n=1}^{N_K} {\cal O}^{(k_n)}\ .
\end{equation}
When the sum is truncated at some finite value of $N_K$ (as in
laboratory measurements), we obtain an {\bf estimator of the 
expectation value}\index{expectation value!estimator}
\begin{equation} \label{estimate_O_B}
  {\overline {\cal O}} =
  {1\over N_K} \sum_{k=1}^{N_K} {\cal O}^{(k_n)}
\end{equation}
whose errors bars are ultimately determined by the statistical 
analysis methods of chapter~\ref{Markov_Error}.

Normally, we cannot generate configurations $k$ directly with
probability~(\ref{P_B}). But they may be found as members
of the equilibrium distribution of a dynamic process. In practice
\index{Markov chain Monte Carlo|textbf}{Markov chains} are used. A 
\index{Markov process|textbf}{\bf Markov process} is a particularly
simple dynamic process, which generates configuration $k_{n+1}$
stochastically from configuration $k_n$, so that no information
about previous configurations $k_{n-1}, k_{n-2}, \dots$ is needed.
The elements of the Markov chain \index{time series|textbf}{\bf time 
series} are the configurations.  Assume that the configuration $k$ is 
given. Let the\index{Markov process!transition probability}
\index{transition probability}{transition probability}to create the 
configuration $l$ in one step from $k$ be given by $W^{(l)(k)} = 
W[k\to l]$. The\index{Markov process!transition matrix|textbf}
\index{transition matrix|textbf}{\bf transition matrix}
\begin{equation} \label{W}
  W = \left( W^{(l)(k)}\right)
\end{equation}
defines the Markov process. Note, that this is a very big matrix, 
because its labels are the configurations (it is never stored in the 
computer!).  To achieve our goal to generate configurations with the 
desired probabilities~(\ref{P_B}), the matrix $W$ is required 
to satisfy the following properties \cite{Wo57,Sa59}:

\begin{description}

\item{(i)} \index{ergodicity}{\bf Ergodicity}:
\begin{equation} \label{ergodicity}
  e^{-\beta E^{(k)}} > 0 ~~{\rm and}~~ e^{-\beta E^{(l)}} > 0
  ~~~{\rm imply:}
\end{equation}
an integer number $n>0$ exists so that $(W^n)^{(l)(k)}>0$ holds.

\item{(ii)} \index{Markov process!normalization|textbf}{\bf 
Normalization}:
\begin{equation} \label{normalization}
  \sum_l W^{(l)(k)} = 1\ .
\end{equation}

\item{(iii)} \index{balance|textbf}{\bf Balance}: 
\begin{equation} \label{balance}
  \sum_k W^{(l)(k)}\, e^{-\beta E^{(k)}}\ =\ e^{-\beta E^{(l)}}\ .
\end{equation}
Balance means: The Boltzmann state~(\ref{P_B}) is an eigenvector with 
eigenvalue $1$ of the matrix $W = (W^{(l)(k)})$. It is an excercise
in appendix~\ref{Exercises} to find some transition matrices, which
fulfill balance but not ergodicity.

\end{description}

To proceed, we have to elaborate on the ensemble notation. By 
definition, an \index{ensemble}{\bf ensemble} is a collection of
configurations so that to each configuration $k$ a probability
$P^{(k)}$ is assigned, $\sum_k P^{(k)} = 1$. The
\index{ensemble!Gibbs-Boltzmann}{\bf Gibbs or Boltzmann ensemble} 
\index{Gibbs ensemble}\index{Boltzmann ensemble}$E_B$ is defined to 
be the ensemble with probability distribution~(\ref{P_B}). 
An\index{ensemble!equilibrium} 
\index{equilibrium ensemble}{\bf equilibrium ensemble} 
$E_{eq}$ of the Markov process is defined by its probability 
distribution $P_{eq}$ satisfying
\begin{equation} \label{E_eq}
 W\, P_{eq} = P_{eq}\, , ~~{\rm in\ components}~~
 P_{eq}^{(l)}=\sum_k W^{(l)(k)} P_{eq}^{(k)}\ .
\end{equation}
The following statement is central for the MC method.
\smallskip

\noindent{{\bf Statement:} Under} the conditions (i), (ii) and (iii)
the Boltzmann ensemble is the {\bf only} equilibrium ensemble of the
Markov process and an attractive\index{Markov process!fixed-point} 
fixed-point\index{fixed-point}\index{attractive fixed-point} (compare 
equation~(\ref{MC_chain}), (\ref{MC_convergence})).
\smallskip

Proof: Let us first  define a distance between ensembles. Suppose 
we have two ensembles $E$ and $E'$. Denote the probability for 
configuration $k$ in $E$ by $P^{(k)}$ and in $E'$ by ${P'}^{(k)}$. 
We define the \index{ensemble!distance between}{\bf distance} between 
$E$ and $E'$ to be
\begin{equation} \label{distance}
  || E - E' ||\ =\ \sum_k | P^{(k)} - {P'}^{(k)} |\ , 
\end{equation}
where the sum goes over all configurations. Suppose that $E'$
resulted from the application of the transition matrix $W$ to the
ensemble $E$. We can compare the distance of $E'$ from the Boltzmann
ensemble with the distance of $E$ from the Boltzmann ensemble:
\begin{eqnarray} \index{balance}\index{triangle inequality}
  & &|| E' - E_B ||\ =\ \sum_l\, \left|\, \sum_k W^{(l)(k)}\, 
  (P^{(k)}-P_B^{(k)})\, \right|\ \qquad ({\rm using\ balance}) 
  \nonumber\\ &\le & 
  \sum_l \sum_k \left|\, W^{(l)(k)}\, (P^{(k)}-P_B^{(k)})\, \right|\ 
  \qquad ({\rm using\ the\ triangle\ inequality}) \nonumber\\
  &=& \sum_k \left|\, P^{(k)}-P_B^{(k)} \,\right|\ =\ 
  ||E-E_B||\ . \label{contraction}
\end{eqnarray}
The last line is obtained by making use of the
\index{Markov process!normalization}{normalization}
condition $\sum_l W^{(l)(k)}=1$ and of $W^{(l)(k)}\ge 0$. This shows 
that the Markov process can reduce the distance between an ensemble $E$
and the Boltzmann ensemble $E_B$ and will never increase this distance. 
Ergodicity is needed to proof that the Boltzmann ensemble is the only 
equilibrium ensemble and to rule out the equality in the $\le$ sign in 
equation (\ref{contraction}). We now study the approach to the fixed
point and derive and explicit equation for the convergence.

The matrix $W$ by itself is normally
not ergodic. Instead, due to the finiteness of the system, a number
$n$ exists, so that the matrix ${\cal W}=W^{n}$ is ergodic. This 
means, all matrix elements of ${\cal W}$ are larger than zero. In 
particular, excluding states with infinite energy by condition~(i)
from our considerations,
\begin{equation} \label{w_min}
  1 > w_{\min} = \min_{k,l} \left({\cal W}^{(k)(l)} \right) > 0
\end{equation}
holds. Let us assume that the state vectors $P$ and $P'$ are related by
\begin{equation} 
  P'\ =\ {\cal W}\, P\,.
\end{equation}
As shown in the following, this implies the inequality
\begin{equation} 
  || E' - E_B ||\ \le\ (1-w_{\min})\, || E - E_B ||\ .
\end{equation}
Let $\epsilon = ||E-E_B||$. We can decompose the contributions to
$||E-E_B||$ in $P^{(k)}-P_B^{(k)}\ge 0$ and in $P^{(k)}-P_B^{(k)}<0$ 
contributions,
\begin{equation} \label{decompose}
  \epsilon = ||E-E_B|| = \sum_{k\in K^+}\left(P^{(k)}-P_B^{(k)}\right) 
                    + \sum_{k\in K^-}\left(P_B^{(k)}-P^{(k)}\right)\,.
\end{equation}
The normalization $\sum_k P^{(k)}=\sum_k P_B^{(k)}=1$ implies
\begin{equation} \label{Psum_epsilon}
  \sum_{k\in K^+}\left(P^{(k)}-P_B^{(k)}\right) =
  \sum_{k\in K^-}\left(P_B^{(k)}-P^{(k)}\right) = \epsilon/2
\end{equation}
as is seen by adding 
$$ 0= \sum_{k\in K^+}\left(P^{(k)}-P_B^{(k)}\right) 
    + \sum_{k\in K^+}\left(P_B^{(k)}-P^{(k)}\right) $$ 
to equation~(\ref{decompose}) and, similarly, this equation with $K^+$ 
replaced by $K^-$. Using ${\cal W} P_B=P_B$, we have
\begin{eqnarray} \label{W_decomposition}
&~& ||E'-E_B||\ = \\ \nonumber &~& \sum_l \left|\,
  \sum_{k\in K^+} {\cal W}^{(l)(k)}\,\left( P^{(k)}-P_B^{(k)}\right) -
  \sum_{k\in K^-} {\cal W}^{(l)(k)}\,\left( P_B^{(k)}-P^{(k)}\right)\,
  \right| \\ \nonumber
&=& \sum_{l\in L^+} \left( 
    \sum_{k\in K^+}{\cal W}^{(l)(k)}\left(P^{(k)}-P_B^{(k)}\right) 
  - \sum_{k\in K^-}{\cal W}^{(l)(k)}\left(P_B^{(k)}-P^{(k)}\right) 
    \right) \\ \nonumber
&+& \sum_{l\in L^-} \left( 
    \sum_{k\in K^-}{\cal W}^{(l)(k)}\left(P_B^{(k)}-P^{(k)}\right) 
  - \sum_{k\in K^+}{\cal W}^{(l)(k)}\left(P^{(k)}-P_B^{(k)} \right)
    \right)
\end{eqnarray}
Here $L^+$ is defined as the set of configurations $l$ which fulfills
the inequality
$$ \sum_{k\in K^+}{\cal W}^{(l)(k)}\left(P^{(k)}-P_B^{(k)}\right) 
 \ge \sum_{k\in K^-}{\cal W}^{(l)(k)}\left(P_B^{(k)}-P^{(k)}\right)\ ,$$
while $L^-$ is the set of configurations $l$ with
$$ \sum_{k\in K^-}{\cal W}^{(l)(k)}\left(P_B^{(k)}-P^{(k)}\right)
 > \sum_{k\in K^+}{\cal W}^{(l)(k)}\left(P^{(k)}-P_B^{(k)}\right)\ .$$
The normalization (\ref{normalization}) and that the smallest matrix 
element is larger than zero (\ref{w_min}) imply
$$ 1\ \ge\ \sum_{l\in L^+} {\cal W}^{(l)(k)}\ \ge\ w_{\min} $$
and the same equation with $L^+$ replaced by $L^-$. Inserting this
together with equation (\ref{Psum_epsilon}) into (\ref{W_decomposition})
yields the result
\begin{equation} \label{W_inequality}
  ||E'-E_B|| \le
   (1-w_{\min})\,\epsilon/2\ +\ (1-w_{\min})\,\epsilon/2\ 
     =\ (1-w_{\min})\,\epsilon\ .
\end{equation}
Under repeated application of the matrix we obtain states
\begin{equation} \label{MC_chain}
  P^t\ =\ {\cal W}^t\,P,~~{\rm with}~~ t=1,2,\dots 
\end{equation}
and find for the {\bf approach to the} 
\index{equilibration}{\bf equilibrium ensemble}
\begin{equation} \label{MC_convergence}
||E^t - E_B||\ \le\  \exp ( - \lambda\,t)\,||E-E_B||
\end{equation}
where $\lambda=-\ln(1-w_{\min})>0$. Therefore the equilibrium ensemble 
is an \index{Markov process!fixed-point|textbf}{\bf attractive 
fixed-point.}\index{fixed-point|textbf}
\index{attractive fixed-point|textbf} \smallskip

There are many ways to construct a Markov process satisfying (i), (ii)
and (iii). For instance, MC algorithms can be  based on products of
steps which satisfy a stronger condition than balance~(\ref{balance}), 
namely:

\begin{description}

\item{(iii')} \index{detailed balance}\index{balance!detailed}{\bf 
Detailed balance} (microscopic reversibility):
\begin{equation} \label{detailed_balance}
  W^{(l)(k)}\, e^{-\beta E^{(k)}}\ =\ W^{(k)(l)} e^{-\beta E^{(l)}}\ .
\end{equation}
Using the normalization condition $\sum_k W^{(k)(l)} = 1$ detailed
balance implies balance (iii).

\end{description}

At this point we have replaced the canonical ensemble average by a 
time average over an artificial dynamics. Calculating averages over 
large times, like one does in real experiments, is then equivalent to 
calculating averages of the ensemble. One 
distinguishes\index{dynamic universality classes}
\index{universality!dynamic}{\it dynamic universality classes}.  
The Metropolis and heatbath algorithms discussed in the following
fall into the class of so called 
\index{Glauber dynamics}{\it Glauber dynamics}~\cite{Gl63,Kw72}, 
imitating \index{fluctuations!thermal}{thermal fluctuations} of nature 
to some extent. In a frequently used classification~\cite{ChLu97} this 
is called \index{model A}{\it model~A}. Cluster algorithms discussed in 
chapter~\ref{AUS} constitute another universality class~\cite{SwWa87}.
Some recent attention has focused on dynamical universality classes
of non-equilibrium systems, see \cite{Od03} for a review.

\subsection{The Metropolis algorithm} \label{Metropolis_algorithm}

Detailed balance still does not uniquely fix the transition
probabilities $W^{(l)(k)}$. 
The\index{algorithms!Metropolis}\index{MC algorithms!Metropolis} 
\index{Metropolis algorithm}{\bf Metropolis algorithm}~\cite{Me53} 
is a popular choice because of its generality and computational 
simplicity. It can be used whenever one knows how to calculate the 
energy of a configuration. The original formulation uses Boltzmann 
weights~(\ref{w_B}). But, replacing the Boltzmann vection in 
(\ref{balance}) by another vector, the algorithm generalizes 
immediately to those weights, what is of importance for 
generalized ensemble simulations, see chapter~\ref{MUCA}.
Given a configuration $k$, the Metropolis algorithm proposes
new configurations $l$ with some transition probabilities
\begin{equation} \label{Metropolis_f}
 f(l,k) ~~{\rm with\ normalization}~~~ \sum_l f(l,k)=1\ .
\end{equation}
We show that the symmetry condition~(\ref{symmetry}) below is needed 
to ensure detailed balance. 

The new configuration $l$ is accepted with probability
\begin{equation} \label{Metropolis}
  w^{(l)(k)} = \min \left[ 1,\, {P_B^{(l)}\over P_B^{(k)}} \right] 
  = \left\{ \begin{matrix} 1  & {\rm for}& E^{(l)}<E^{(k)}\\
  e^{-\beta (E^{(l)}-E^{(k)})} &{\rm for}& E^{(l)}>E^{(k)}\,.
  \end{matrix}\right.
\end{equation}
If the new configuration is rejected ({\it i.e.}, not accepted), the 
old configuration is kept, and has to be counted again\footnote{
Although this point is already outlined in the original \cite{Me53}
paper, it has remained a frequent mistake of beginners to count 
configurations only after accepted changes.}. 
The \index{Metropolis algorithm!acceptance rate}
\index{acceptance rate|textbf}{\bf acceptance rate}
is defined as the ratio of {\bf accepted changes over proposed 
moves}. With this definition we do not count a move as accepted
when it proposes the at hand configuration (some authors
may count these proposals as accepted).

The Metropolis procedure gives rise to the transition probabilities 
\begin{eqnarray} \label{Metropolis_W_ne}
 W^{(l)(k)} &=& f(l,k)\, w^{(l)(k)}\ \ {\rm for}\ \ l \ne k\\
{\rm and}~~~ W^{(k)(k)} &=& f(k,k) + \sum_{l\ne k} f(l,k)\, 
             (1-w^{(l)(k)})\ .  \label{Metropolis_W_eq}
\end{eqnarray}
Therefore, the ratio $\left( W^{(l)(k)}/W^{(k)(l)} \right)$ satisfies 
detailed balance~(\ref{detailed_balance}) if
\begin{equation} \label{symmetry}
 f(l,k)\ =\ f(k,l) 
\end{equation}
holds. Otherwise the probability density $f(l,k)$ is unconstrained
and can be chosen conveniently, {\it i.e.}, there is an amazing
flexibility in the choice of the transition probabilities. One 
can even use acceptance probabilities distinct from those of
equation~(\ref{Metropolis}) and the proposal probabilities are then
not necessarily symmetric anymore \cite{Ha70}. For reference we 
introduce a {\bf biased Metropolis}\index{biased Metropolis}
\index{Metropolis algorithm!biased}acceptance rate
\begin{equation} \label{Metropolis_biased}
  w^{(l)(k)}_b = \min \left[ 1,\, {P_B^{(l)}\,f(k,l)\over 
  P_B^{(k)}\,f(l,k)} \right]\ .
\end{equation}
Obviously this equation agrees with (\ref{Metropolis}) when $f(l,k)$
is symmetric and it is an easy excersise to proof that the detailed
balance condition (\ref{detailed_balance}) now holds even when
$f(l,k)$ is no longer symmetric. This has occasionally been used
in the statistical physics \cite{Br85,MiHe86} and bio-chemical
\cite{DeBa96,Fa01,Jo03} literature. In chapter~\ref{LGT} on lattice 
gauge theory we use this acceptance probability to immitate heatbath 
probabilities while keeping up with many of the conveniences of a 
Metropolis type algorithm. This leads to a Biased Metropolis-Heatbath 
Algorithm (BMHA), which can be used in situations when an efficient 
heatbath algorithm (see below) is not available \cite{BaBe05}.

At this point, the reader is recommended to go through the simple
exercise of writing a Metropolis code which generates the 
\index{Metropolis algorithm!Gaussian distribution}{Gaussian 
distribution} (assignment~\ref{assign_Metro_Gauss}). Next we 
discuss the Metropolis algorithm for a more involved continuous 
example, which we also use to introduce the heatbath algorithm. 
The reader interested in getting quickly started with simulations 
may prefer to jump to section~\ref{Potts}.

\subsection{The $O(3)$ $\sigma$-model and the heatbath algorithm}
\label{O3_heatbath}

Here we give one example of a model with a continuous\footnote{
Until digitized by the computer.} energy function, the 
\index{O(3) $\sigma$-model}{$O(3)$ $\sigma$-model}. 
The $2d$ version of this model is of interest to field theorists 
because of its analogies with the four-dimensional Yang-Mills 
theory~\cite{Po75,BePo75,BrZi76}.
In statistical physics the $d$-dimensional model is known as the 
\index{ferromagnet!Heisenberg}\index{Heisenberg ferromagnet}{Heisenberg 
ferromagnet}, see \cite{Hu87}. The $3d$, $n=3$ model has a phase 
transition, which is likely third order (although second order is 
not entirely ruled out) according to extensive investigations of the 
critical exponents of this transition, see \cite{HoJa97} for 
references.

Expectation values of the model are calculated with respect to the 
partition function
\begin{equation} \label{O3_Z}
 Z\ =\ \int \prod_i ds_i\ e^{-\beta E(\lbrace s_i \rbrace )}\ .
\end{equation}
The spins
\begin{equation} \label{O3_spin}
 \vec{s}_i = \left( \begin{matrix}s_{i,1}\\ s_{i,2}\\ s_{i,3}
                    \end{matrix} \right) 
           ~~{\rm are\ normalized\ to}~~ (\vec{s}_i)^2=1
\end{equation}
and are defined on the sites of a lattice of size
$n_1\times n_2 \times \dots \times n_{d-1} \times n_d$. 
The measure $ds_i$ is defined by
\begin{equation} \label{O3_measure}
 \int ds_i\ =\ {1\over 4\pi} \int_{-1}^{+1} d \cos (\theta_i) 
                              \int_0^{2\pi}  d \phi_i\ , 
\end{equation}
where $\theta_i$ and $\phi_i$ are polar and azimuth angles defining
the spin $s_i$ on the unit sphere. The energy is
\begin{equation} \label{O3_energy}
\index{continuous systems!O(n) models} 
\index{systems!continuous!O(n) models}
\index{O(3) $\sigma$-model!energy} 
\index{energy!Heisenberg ferromagnet}\index{energy!O(n) models}
 E\ =\ \sum_{<ij>} ( 1 - \vec{s}_i \vec{s}_j )\ ,  
\end{equation}
where the sum goes over the nearest neighbor sites of the lattice. 
We would like to change \index{update}(update) a single spin 
$\vec{s}$ (dropping site indices $i$ for the moment). The sum 
of its $2d$ neighbors is
$$\vec{S} = \vec{s}_1+\vec{s}_2+\dots +\vec{s}_{2d-1}+\vec{s}_{2d}\ .$$
Hence, the contribution of spin $\vec{s}$ to the energy is
$2d-\vec{s} \vec{S}$. We propose a new spin $\vec{s}^{\,'}$ with the
measure~(\ref{O3_measure}). This defines the probability function
$f(\vec{s}^{\,'},\vec{s})$, equation~(\ref{Metropolis_f}), of the 
Metropolis process. In practice, we get the new spin $\vec{s}^{\,'}$
by simply drawing two uniformly distributed random numbers
$$ \phi^r \in [0,\pi) ~~{\rm for\ the\ azimuth\ angle\ and} $$
$$ \cos (\theta^r) = x^r \in [-1,+1) ~~{\rm for\ the\ cosine\ of\
 the\ polar\ angle}.  $$
The Metropolis algorithm accepts the proposed spin
$\vec{s}^{\,'}$ with probability 
$$ w(\vec{s}\to \vec{s}^{\,'})\ =\ \left\{ \begin{matrix} 
   1 & {\rm for} & \vec{S} \vec{s}^{\,'} \ge \vec{S} \vec{s} , \\
 e^{ -\beta (\vec{S} \vec{s} - \vec{S}\vec{s}^{\,'}) } & {\rm for} &
 \vec{S}\vec{s}^{\,'} < \vec{S}\vec{s} . \end{matrix} \right.    $$
The algorithm proceeds by choosing another site and
repeating the above procedure. If sites are chosen with the 
uniform probabilities $1/N$ per site, where $N$ is the total 
number of spins, the procedure fulfills detailed balance: 
Just combine the spin acceptance probabilities with the probability 
function $f(\vec{s}^{\,'},\vec{s})$ and the
thus obtained transition 
probabilities~(\ref{Metropolis_W_ne}),~(\ref{Metropolis_W_eq}) 
fulfill~(\ref{detailed_balance}). It is noteworthy that the 
procedure remains valid when the spins are chosen in the 
systematic order $1,\dots ,N$, then $1,\dots ,N$ again, and 
so on. Balance~(\ref{balance}) still holds, whereas detailed 
balance~(\ref{detailed_balance}) is violated. This is an
excercise in appendix~\ref{Exercises}. 
\medskip

\noindent {\bf Heatbath algorithm} \medskip

Repeating the Metropolis algorithm again and again for the same
spin $\vec{s}$ leads (as has been proven) to the equilibrium 
distribution of this spin with respect to its (fixed) neighbors, 
which reads
$$ P(\vec{s}^{\,'};\vec{S}) = {\rm const}\,
  e^{\beta \vec{S}\vec{s}{\,'}} ~~{\rm with}~~
   \int P(\vec{s}^{\,'};\vec{S})\, ds' = 1\ .  $$
One would prefer to choose $\vec{s}^{\,'}$ directly with the
probability
$$ W(\vec{s}\to \vec{s}^{\,'})\ =\ P(\vec{s}^{\,'};\vec{S})\ =\
  {\rm const}\, e^{\beta\, \vec{s}^{\,'} \vec{S}}\ . $$
The algorithm, which creates this distribution, is called the 
\index{algorithms!heatbath}\index{MC algorithms!heatbath} 
\index{MC algorithms!O(n) models}\index{heatbath algorithm}{\bf 
heatbath} algorithm, because it has the physical interpretation that 
the spin is brought in contact with a heatbath at the particular 
temperature $\beta = 1/T$. It was already known to \cite{Gl63}.
The first computer implementation seems to be due to \cite{Cr80} in
lattice gauge theory. Afterwards it was frequently used by physicists,
{\it e.g.}, \cite{Pi81,BeLu81,CaMa82}.
Applied mathematicians \cite{Ge84} named 
the algorithm the \index{Gibbs Sampler: see heatbath}{\bf Gibbs 
Sampler}, when applying it successfully to Baysian restoration of
images.

Implementation of the heatbath algorithm becomes feasible
when the energy function is sufficiently simple to allow for
a fast calculation of the probability $P(\vec{s}^{\,'};\vec{S})$.
Setting it up is an easy task for the $O(3)$ $\sigma$-model. Let
$$ \alpha = {\rm angle} (\vec{s}^{\,'},\vec{S}),\ \
   x=\cos (\alpha )\ \ {\rm and}\ \  S_{\beta}=\beta |\vec{S}|\ . $$
For $S_{\beta}=0$ a new spin $\vec{s}^{\,'}$ is simply obtained 
by random sampling. We assume in the following $S_{\beta}>0$. The 
Boltzmann weight becomes $\exp (x S_{\beta})$ and the normalization
constant follows from
$$ \int_{-1}^{+1} dx\, e^{x S_{\beta}}\ =\ 
   {2\over S_{\beta}}\, \sinh (S_{\beta})\ .$$
Therefore, the desired probability is
$$ \index{heatbath algorithm!O(3) $\sigma$-model|textbf}
   \index{O(3) $\sigma$-model!heatbath algorithm}
   \index{heatbath algorithm!Heisenberg ferromagnet|textbf}
   \index{Heisenberg ferromagnet!heatbath algorithm|textbf}
  P(\vec{s}^{\,'};\vec{S})\ =\  {S_{\beta}\over 2\sinh(S_{\beta})}\, 
  e^{x S_{\beta}}\ =: f(x)\ . $$
The method of chapter~\ref{Sampling}, equation~(\ref{gd}), can
be used to generate events with the probability density $f(x)$. With
\begin{eqnarray} \label{O3_cdf}
  y &=& 1-F(x) \,=\, 1 - \int_{-1}^x dx'\, f(x') \,=\,
  \int_x^1 dx'\, f(x')\\ \nonumber
  &=& \int_x^1 dx'\, {S_{\beta} \over 2\sinh (S_{\beta})}\, 
  e^{x' S_{\beta}}\,=\,{\exp(+S_{\beta})-\exp(xS_{\beta})
                        \over 2\,\sinh (S_{\beta})}
\end{eqnarray}
a uniformly distributed random number $y^r\in [0,1)$ translates into
\begin{equation} \label{O3_heat_bath-xr}
  x^r = \cos{\alpha^r} = {1\over S_{\beta}} \ln\, [\, \exp (+S_{\beta}) 
      - y^r\, \exp (+S_{\beta}) + y^r\, \exp (-S_{\beta}) ]\ . 
\end{equation}
Finally, one has to give $\vec{s}^{\,'}$ a direction in the plane
orthogonal to $\vec{S}$. This is done by choosing a random angle 
$\beta^r$ uniformly distributed in the range $0\le\beta^r < 2\pi$. 
Then, $x^r=\cos{\alpha^r}$ and $\beta^r$ completely determine
$\vec{s}^{\,'}$ with respect to $\vec{S}$. Before storing
$\vec{s}^{\,'}$ in the computer memory, we have to calculate
coordinates of $\vec{s}^{\,'}$ with respect to a Cartesian
coordinate system $K$, which is globally used for all spins of 
the lattice. This is achieved by a \index{linear transformation}
\index{transformations, linear}{linear transformation}. We define
$$ \cos \theta = {S_3\over S},\
   \sin \theta = \sqrt{1-\cos^2\theta},\
   \cos\phi = {S_1\over S\,\sin\theta} ~~{\rm and}~~
   \sin\phi = {S_2\over S\,\sin\theta}\,, $$
where $S_1$, $S_2$ and $S_3$ are the components of $\vec{S}$ in $K$. 
Unit vectors of a coordinate frame $K'$, with $\hat{z}'$ in the 
direction of $\vec{S}$ and $\hat{y}'$ in the $x-y$ plane, are 
then defined by
$$ \hat{z}'= \left( \begin{matrix} \sin\theta\,\cos\phi \\
             \sin\theta\,\sin\phi \\ \cos\theta \end{matrix} \right),~~
   \hat{x}'= \left( \begin{matrix} \cos\theta\,\cos\phi \\
            \cos\theta\,\sin\phi \\ -\sin\theta \end{matrix} \right)
~~{\rm and}~~ \hat{y}'= \left( \begin{matrix} - \sin\phi \\
                        \cos\phi \\ 0 \end{matrix} \right)\ .$$
Expanding $\vec{s}^{\,'}$ in these units vectors, we have
\begin{eqnarray} \label{stored_vector}
   &\vec{s}^{\, '}& =\ \sin\alpha^r\,\cos\beta^r\,\hat{x}' +
   \sin\alpha^r\sin\beta^r\,\hat{y}' +
   \cos\alpha^r\,\hat{z}' \\ \nonumber
   &=& \left( \begin{matrix}
   \sin\alpha^r\cos\beta^r\,\cos\theta\,\cos\phi
 - \sin\alpha^r\sin\beta^r\,\sin\phi
 + \cos\alpha^r\,\sin\theta\,\cos\phi \\
   \sin\alpha^r\cos\beta^r\,\cos\theta\,\sin\phi
 + \sin\alpha^r\cos\beta^r\,\cos\phi
 + \cos\alpha^r\,\sin\theta\,\sin\phi \\
 - \sin\alpha^r\cos\beta^r\,\sin\theta
 + \cos\alpha^r\,cos\theta \end{matrix} \right)
\end{eqnarray}
and the three numbers of the column vector (\ref{stored_vector}) are 
stored in the computer. The Fortran implementation of this heatbath 
algorithm is given in section~\ref{O3_hb}.

\subsection{Assignments for section~\ref{Importance_sampling} 
\label{assign_Metropolis} }

\begin{enumerate}

{\item Write a Metropolis program to generate the normal Gaussian
distribution through the Markov process
$$ x\, \to\, x' = x + 2\, a\, (x^r-0.5) $$
where $x^r$ is a uniformly distributed random number in the range 
$[0,1)$ and $a$ is a real constant.
Use $a=3.0$ and the initial value $x=0.0$ to test your program.
(i)~Generate a chain of 2$\,$000 events. Plot the empirical 
\index{peaked distribution function!Gaussian}{peaked distribution 
function} for these event in comparison with the exact
peaked normal distribution function. Monitor also the acceptance 
rate. (ii)~Repeat the above for a chain of 20$\,$000 events. 
\label{assign_Metro_Gauss} }

\end{enumerate}
 
\section{Potts Models Metropolis and Heatbath Code\label{Potts}}

We introduce generalized\index{Potts models|textbf}
\index{simulations!Potts models}{Potts models}
in an  external \index{Potts models!magnetic field}{magnetic field}
on $d$-dimensional hypercubic lattices with periodic boundary
conditions. Without being overly complicated, these models are general
enough to illustrate the essential features we are interested in.
Various special cases of these models are by themselves of physical 
interest, see references in \cite{BlNi82,Wu82}. Generalizations of 
the algorithmic concepts to other models are straightforward, 
although technical complications may arise. 

We define the\index{energy!Potts models}\index{energy!Potts models}
\index{Potts models!energy}{energy} function of the system by
\begin{equation} \label{E}
  \index{discrete systems!Potts models|textbf}
  \index{systems!discrete!Potts models|textbf}
  -\beta\,E^{(k)} = -\beta\,E_0^{(k)} + H\,M^{(k)} 
\end{equation}
where
\begin{eqnarray} \label{E0}
  E_0^{(k)} &=& - 2\, \sum_{<ij>} J_{ij}(q_i^{(k)},q_j^{(k)})\
  \delta (q_i^{(k)},q_j^{(k)}) + {2\,d\,N\over q}\\ \nonumber
  &{\rm with}&~\ \delta (q_i,q_j)
  = \left\{ \begin{array}{c} 1\ {\rm for}\ q_i=q_j \\
  0\ {\rm for}\ q_i\ne q_j \end{array}  \right.
\end{eqnarray}
and
\begin{equation} \label{M}
  M^{(k)} = 2 \sum_{i=1}^N \delta (0,q_i^{(k)})\, ,  
\end{equation}
so that the Boltzmann factor is $\exp(-\beta\,E^{(k)}) =
\exp(-\beta\,E_0^{(k)}+H\,M^{(k)})$. In equation (\ref{E0}) the sum 
$<ij>$ is over the nearest neighbor lattice sites and $q_i^{(k)}$ 
is the \index{Potts models!spin (state)}{\bf Potts spin state}
\index{Potts models!state (spin)}\index{state!Potts spin}at site 
$i$. This should not be confused with the microstate which is the 
configuration $k$. Occasionally the notation\index{q-state}
\index{Potts models!q-state}{\bf $q$-state} is used for the Potts 
spin states at some site. The $q$-state $q^{(k)}_i$ takes on the
values\footnote{In our Fortran programs the instruction {\tt 
int(nq*x)}, where {\tt x} is a uniformly distributed random number, 
generates random integers in the range ${\tt iq=0,}\dots {\tt ,nq-1}$.}
$0,\dots ,q-1$. The sum in equation (\ref{M}) goes over the $N$ sites 
of the lattice.  The external magnetic field is chosen to interact 
with the state $q_i=0$ at each site $i$, but not with the other states 
$q_i\ne 0$. The $J_{ij} (q_i,q_j)$ functions define the exchange
coupling constants between the spin states at sites $i$ and~$j$.

The energy function~(\ref{E0}) describes a number of physically
interesting situations. With
\begin{equation} \label{Jij}
  J_{ij} (q_i,q_j) \equiv J > 0~~\ ({\rm conventionally}\ J=1) 
\end{equation}
the original 
\index{Potts}\index{Potts models!original Potts model}model
\cite{Po52} is recovered\footnote{Actually, the article by Potts 
deals with the clock model (which C. Domb proposed to him as his 
Ph.D. thesis) and he introduces the model~(\ref{E0}) only as
a footnote in his thesis.} and $q=2$ becomes equivalent to the Ising 
\index{ferromagnet!Potts}{ferromagnet}, see \cite{Wu82} for a 
detailed review. The Ising case of \cite{EdAn75} spin 
glasses and quadrupolar Potts glasses~\cite{Bi97} are obtained when the 
exchange constants are quenched random variables.  Other choices of the 
$J_{ij}$ include anti-ferromagnets and the fully frustrated Ising 
model~\cite{Vi77}.

For the \index{Potts models!energy!per spin}{\bf energy per spin} 
\index{energy!per spin}we use the notations
\begin{equation} \label{Potts_es}
  e_{0s} = E_0/N ~~{\rm and} ~~ e_s = E/N\,,
\end{equation}
depending on whether the external magnetic field $H$ is zero or not 
(the configuration indices are often often suppressed). In equation
(\ref{E0}) the factor of two in front of the sum for $E_0^{(k)}$ is 
introduced, so that $e_{0s}$ agrees for $q=2$ with our Ising model 
definition~(\ref{e0s}) with energy~(\ref{Ising_E}). For sampling with 
respect to the partition function (\ref{Z}), our definition of $\beta$ 
agrees also with the one of the Ising model as commonly used in the 
literature, {\it e.g.} \cite{Ch87} ($J=1$ there) or \cite{Hu87}. But 
it disagrees by a factor of two with the $\beta$ definition used in 
the Potts model literature \cite{Ba73,BoJa92}:
\begin{equation} \label{beta_Potts}
  \beta\ =\ \beta^{\rm Ising}\ =\ {1\over 2}\, \beta^{\rm Potts}\, .
\end{equation}
Again by reasons of consistency with Ising model notation, a factor 
of two multiplies the magnetic field in equation~(\ref{M}). 

The expectation value 
\begin{equation} \label{Mq}
  M_{q_0}\ =\ \langle\, \delta_{q_i,q_0}\, \rangle 
\end{equation}
to find the Potts spin at site $i$ to be $q_0$ is called {\bf 
Potts magnetization}\index{Potts models!magnetization|textbf} 
\index{magnetization|textbf}with respect to $q_0$. For the 2-state 
model the usual Ising magnetization is then found as $m = 
\langle\delta_{q_i,0}\rangle - \langle\delta_{q_i,1}\rangle$.
For $\beta =0$ it is rather trivial to calculate the dependence of 
the Potts magnetization on the external magnetic field~(\ref{M}). 
One finds (this an exercise in appendix~\ref{Exercises})
\begin{equation} \label{beta0_H}
  \langle \delta_{q_i,0} \rangle = {\exp (2H)\over \exp(2H)+q-1} 
\end{equation}
at any site $i$ independently of the lattice dimension and size. 

For $2d$ Potts models a number of exact results are known in the
infinite volume limit, mainly due to work by\index{Baxter} \cite{Ba73}. 
The critical temperature~(\ref{2d_Ising_bc}) of the $2d$ Ising model 
generalizes to the \index{Potts models!phase transitions}{transition}
temperatures
\begin{equation} \label{2d_Potts_bc}
{1\over 2}\, \beta_c^{\rm Potts}\ =\ \beta_c\ =\ {1\over T_c}\ 
=\ {1\over 2}\, \ln (1 + \sqrt{q} ), ~~q=2,3,\dots\ .
\end{equation}
At $\beta_c$ the average energy per state is given by
\begin{equation} \label{Potts_e0s}
 e_{0s}^c = E_0^c/N = {4\over q} - 2 - 2/\sqrt{q}\ .
% = {4\over q} + e_{0s}[{\rm Baxter}]
\end{equation}
As before, by reasons of consistency with the Ising model notation
(\ref{e0s}), our definition of $e_{0s}=E_0/N$ differs by factor of 
two and an additive constant from the one used in most Potts model 
literature. \index{Potts models!phase transitions}{The phase 
transition} is second order for $q\le 4$ 
and\index{first order transition: see phase transition!first order} 
\index{phase transition!first order|textbf}{first order} for $q\ge 5$. 
Conjectures for the exact critical\index{Potts model!critical 
exponents} exponents of the 3- and 4-state model are (compare p.260 
of \cite{Wu82}):
\begin{equation} \label{2d3q_Potts_exponents}
  \alpha = 1/3,\ \beta = {1\over 9},\ \gamma = {13\over 9},\
  \delta = 14,\ \nu = 5/6~~{\rm and}~~\eta = {4\over 15}
  ~~{\rm for}~~q=3
\end{equation}
and 
\begin{equation} \label{2d4q_Potts_exponents}
  \alpha = 1/3,\ \beta = {1\over 9},\ \gamma = {13\over 9},\
  \delta = 14,\ \nu = 5/6~~{\rm and}~~\eta = {4\over 15}
  ~~{\rm for}~~q=4\,.
\end{equation}
The exact infinite volume\index{latent heat!Potts models}
\index{Potts models!latent heat}{\bf latent heats} $\triangle e_{0s}$ 
and \index{entropy!$2d$ Potts}\index{Potts models!entropy}{\bf entropy 
gaps} $\triangle s$ were also derived by \cite{Ba73}. The knowledge of 
$e_{0s}^c$ and $\triangle e_{0s}$ at the critical point allows to 
calculate the limiting energies for $\beta\to\beta_c$ from the 
disordered and ordered phase, $e_{0d}$ and $e_{0o}$ respectively,
\begin{equation} \label{e_d_o}
  e_{0s}^c=(e_{0d}+e_{0o})/2 ~~~{\rm and}~~~
  \triangle e_{0s}=e_{0d}-e_{0o}\ .
\end{equation}
For the first order transitions the\index{interface tension}
\index{Potts models!interface tensions}{interface tensions} $f_s$ 
are known too \cite{BoJa92}. For easy reference these quantities are 
compiled in table~\ref{tab_potts_exact} for some $q\ge 5$ values. 
The results were obtained with the programs\index{\tt beta\_c.f}
\begin{equation} \label{Potts_exact} \index{\tt lheat.f}
  {\tt beta\_c.f}~~{\rm and}~~{\tt lheat.f}~~{\rm of}~~
  {\tt ForProg/MC\_Potts.}
\end{equation}
A generalization, allowing for a dependence of the interaction 
energy on the direction, leads to the chiral Potts model for which 
exact results were found by \cite{Ba05}.

\begin{table}[tbl] \begin{center}
\tbl{Exact, analytical results for the $2d$ $q$-state Potts model: 
$\beta_c$,\index{Potts models!phase transitions}\index{Potts models!$2d$}
$e_{0s}$,\index{Potts models!latent heat}\index{Potts models!entropy}
\index{Potts models!interface tensions}$\triangle e_{0s}$ 
 and $\triangle s$ from [Baxter (1973)] and $f_s$ 
 from [Borgs and Janke (1992)]. }
% \medskip
{ \begin{tabular}{|c|c|c|c|c|}
\hline
                         & $q=5$      & $q=7$   & $q=8$   & $q=10$  \\
\hline
$\beta_c=\ln(1+\sqrt{q})/2$&0.5871795 &0.6467812&0.6712270&0.7130312\\

$e_{0s}=4/q-2-2/\sqrt{q}$&$-$2.0944272&$-$2.1845004&$-$2.2071068
                                                   &$-$2.2324555\\
$\Delta e_{0s}$          &0.105838    &0.706553&0.972716&1.39210\\
$\Delta s=2 \beta_c \Delta e_{0s}$
                         &0.124292    &0.913971&1.305828&1.98522 \\
$2 f_s$                  &0.000398    &0.020792&0.041879&0.094701\\
\hline
\end{tabular} \label{tab_potts_exact} } \end{center} \end{table}

\index{simulations!Potts models}Let us return to our goal of performing 
MC simulations of the model. Each configuration $k$ defines a particular 
arrangements of Potts spins and, vice versa, each arrangement of the 
Potts spins determines uniquely a configuration:
\begin{equation} \label{k}
 k = \{ q_1^{(k)}, \dots , q_N^{(k)} \} \ . 
\end{equation}
As there are $q$ possible states at each site, the total number of
microstates is
\begin{equation} \label{K}
 Z(0) = q^N\ ,  
\end{equation}
where we have used the definition (\ref{Z}) of $Z$. Already for
moderately sized systems $q^N$ is an enormously large number.

First we like to implement the \index{Potts models!Metropolis algorithm}
\index{Metropolis algorithm!Potts models}{Metropolis algorithm}. 
Given a configuration $k$, a new configurations $l$ has to be proposed 
with an explicitly known transition probability (\ref{Metropolis_f}). 
The symmetry condition (\ref{symmetry}) implies detailed balance 
(\ref{detailed_balance}), which can be replaced by the somewhat weaker 
condition of balance (\ref{balance}). The newly proposed configuration 
$l$ is accepted with the probability (\ref{Metropolis}), otherwise 
rejected. For Potts models it is straightforward to realize these 
steps.  The transition probabilities $f(l,k)$ can be defined by the 
following procedure:
\begin{enumerate}
\item Pick one site $i$ at random, {\it i.e.} with probability $1/N$.
\item Assign one of the numbers $0,\,1,\,\dots,\,q-1$, each with 
      probability $1/q$, to $q_i^{(l)}$, 
\end{enumerate}
These rules amount to
\begin{equation} \label{Potts_p0}
 f(l,k) =  \begin{cases} 1/q ~~{\rm for}~~l=k,\\ 
   1/(qN)\ {\rm for\ each\ configuration}\ l\ {\rm with}\ 
   q_i^{(l)}\ne q_i^{(k)}\ {\rm at\ site}\ i,\\
   0\ {\rm for\ all\ other\ configurations}\ l.
 \end{cases}
\end{equation}
The new configuration $l$ is then accepted or rejected according
to equation~(\ref{Metropolis}). Ergodicity is fulfilled, because
with (\ref{Potts_p0}) every configuration can be reached
in $N$ update steps. The random choice of a site ensures detailed 
balance. A procedure where the Potts spins are updated in the 
systematic order of a permutation $\pi_1,\dots,\pi_N$ of the
site numbers leads still to the desired distribution.
Detailed balance is then violated but the weaker condition
of balance still holds (see appendix~\ref{Exercises}, exercise~1 
to section~\ref{Importance_sampling}). 

\subsection{Metropolis code}

\index{MC algorithms!Potts models}The Fortran code presented here 
is for the $J_{ij}(q_i,q_j)\equiv 1$ original Potts model (\ref{Jij}) 
in arbitrary dimensions $d$ with $H=0$. Later on code allowing for 
$H\ne 0$ will be introduced. Generalization to arbitrary 
$J_{ij}(q_i,q_j)$ arrays is straightforward, but to use it for 
the situation where all elements are equal does slow down the 
computational speed and occupies additional {\bf RAM}\index{RAM} 
(Random Access Memory). On the other hand, our code is not optimized 
for plain Potts models either, where one should for instance consider 
\index{multispin coding}multispin coding \cite{JaRe81,Zo81}, where one 
uses the bits of a computer word to store many lattices in one computer 
word. This is particularly efficinet for the Ising model, where a spin 
can be presented by a single bit, but should in that case be combined
with cluster algorithms, which are considered in chapter~\ref{Cluster}.
The code of the present chapter is intended to serve as prototype for 
general applications.

Instead of the energy function~(\ref{E0}) the action variable
(\ref{iact}) is used in the following, and
\begin{equation} \label{E0_iact}
  E_0^{(k)}\ =\ {2\,d\,N\over q} - 2\, {\tt iact}^{(k)}
\end{equation}
holds. Recall that the variable {\tt iact} is convenient because 
it takes on integer values 
$$ {\tt iact} = 0, 1, 2, \dots, {\tt mlink}\ , $$
increasing in steps of two the 2-state Potts model. The {\bf action per 
link}\index{action variable!per link} is defined by
\begin{equation} \label{actl}
  {\tt actl}\ =\ { {\tt iact}\over {\tt nlink} }
\end{equation}
so that it is normalized in the range $0\le {\tt actl}\le 1$. 
\begin{equation} \label{actlm}
  {\rm The\ estimator\ of\ the}\ {\tt actl}\ {\rm expectation\ value\ 
  is\ denoted\ by}\ {\tt actlm}.
\end{equation}
There are now two steps: First we initialize the system, then the 
\index{simulations!Metropolis}{Metropolis simulations} are performed. 
Programs are provided in
\begin{equation} \label{ForProg_Potts}
  {\tt ForProg/MC\_ Potts}\ .
\end{equation}
There are subroutines in {\tt ForLib} and
\begin{equation} \label{ForLib_par}
  {\tt ForLib\_ par}\ .
\end{equation}
The subroutines in {\tt ForLib\_par} all use statements like for 
instance
$$ {\tt include}\ {\tt 'lat.par'} $$
to transfer parameter files into the code, which are kept in the run 
directories and frequently changed by the user (an example of {\tt 
lat.par} is given in table~\ref{tab_latpar}). The parameters may 
dimension arrays in common blocks, like {\tt common\_lat.f} and 
{\tt common\_potts.f}, which the routines of {\tt ForLib\_par} share 
with other routines and/or the main program. To keep these routines
in a library (instead of moving them around with the main programs)
is a major improvement over the first version of the STMC code. It
comes with the mild penalty that all programs which inlcude routines
from {\tt ForLib\_par} have to be {\bf compiled by running the 
{\tt make77} (\ref{make77}) script},\index{make77} instead of 
simply using a Fortran~77 compiler. The rationale for this design 
is that it appears not wise to avoid common blocks in MC code, see 
our discussion in chapter~\ref{Fortran}. The structure of the common 
blocks is not supposed to change frequently, while the dimensions of 
the involved arrays have to stay flexible. 

In the following we will be content with introducing simple programs,
which are suitable for generating energy histogram and time series data. 
Large scale simulations require more sophisticated programs and a 
serious error analysis discussed in chapter~\ref{Markov_Error}. One may 
also consider to speed up the code. For instance, replace the Marsaglia 
random number generator by a faster one. The CPU time used for the 
random number generation is substantial in the Potts models, but almost 
negligible for systems which require a lot of number crunching to 
calculate the energy.

\subsubsection{Initialization} 

Initialization is done by the subroutine 
\begin{equation} \label{potts_init}\index{potts\_init.f}
  {\tt potts\_ init.f}~~~{\rm of}~~~{\tt ForLib\_par}
\end{equation}
which is (up to write statements) reproduced in 
table~\ref{tab_potts_init}. It relies mainly on calls to other 
subroutines, which are part of {\tt ForLib}. In detail 
${\tt potts\_ init.f}$ does the following:

\begin{table}[th]  
\caption{Initialization of the Potts model for Metropolis simulations: 
The relevant parts of the subroutine ${\tt potts\_ init.f}$ are listed 
(omitting write and a few other statements).}
{\small \centering \begin{verbatim}
      subroutine potts_init(lpr) ! Bernd Berg, Mar 31, 2005.
C Initializes arrays for Potts Model MC Metropolis calculations.
      include 'implicit.sta'
      include 'constants.par'
      include 'lat.par'
      include 'potts.par'
      include 'mc.par'
      parameter(i0=0)
      include 'common_lat.f'
      include 'common_potts.f'
C Initialize:
      if(lpr) call rmaset(iuo,iud1,iseed1,iseed2,'nexiste.pa') ! ran-num
      if(.not.lpr) call rmaset(-iuo,iud1,iseed1,iseed2,'nexiste.pa') 
      call ranmar(xr) ! first random number
      acpt=zero                                        ! acceptance rate
      ia_min=mlink
      ia_max=0
      ns=nsfun(nla,nd)
      if(ns.gt.ms) stop "potts_init.f: ns.gt.ms."
      call lat_init(ms,nd,ipf,ipb,nla,ix)           ! lattice
      call potts_act_tab(idel,nqm1)                 ! potts action table
      if(istart.eq.1) call potts_order(ista,iact,nd,ms,izero) ! ordered.
      if(istart.eq.2) call potts_ran(ista,ms,nq) ! random start.
      if(istart.ne.1.and.istart.ne.2) stop "potts_init: istart false."
      call potts_act(ista,ipf,idel,ms,nqm1,nd,iact)       ! potts action
      call potts_wght(wrat,beta,mlink,nd)                ! potts weights
      call razero(ha,i0,mlink)             ! initialize action histogram
      return
      end
      end
\end{verbatim} } \label{tab_potts_init}
\end{table}

\begin{table}[tbl]  
\caption{Parameter block {\tt mc.par}: User defined parameters
for a MC calculations.}
{\centering\small 
\input mc.par      
} \label{tab_mc.par}
\end{table}


\begin{enumerate}

\item {\tt rmaset} initializes the Marsaglia random number generator. 
Of its arguments the input/output units {\tt iud}, {\tt iud1} and 
the seeds {\tt iseed1}, {\tt iseed2} are defined in the parameter 
block {\tt mc.par}\index{\tt mc.par}  of which an example is reproduced 
in table~\ref{tab_mc.par}. Some of these parameters will only be
activated in later chapters of the book.

\item {\tt ranmar} is called for the first random number after
initialization of the random number generator. This allows
to print {\tt xr} for check purposes.

\item ${\tt acpt}=0$ is the initial value of the variable {\tt acpt},
which is part of {\tt common\_potts.f} (table~\ref{tab_potts.com})
and  used to monitor the acceptance rate of the Metropolis updates
(\ref{Metropolis}).

\item {\tt nsfun} checks consistency between the total number of
spins {\tt ms} and the same number calculates from the {\tt nla}
array.

\item {\tt potts\_act\_tab} initializes {\tt idel(0:nqm1,0:nqm1)} 
of {\tt common\_potts.f}. This is the delta function used to calculate 
the Potts model action~(\ref{iact}). In the computer code {\tt nq} 
denotes the number of Potts states $q$ and ${\tt nqm1}={\tt nq}-1$. 
The user has to define {\tt nq} in the parameter statement {\tt 
potts.par}, for example
\begin{equation} \label{potts.par}
  {\tt parameter(nq=2, nqm1=nq-1,nqm2=nq-2) }
\end{equation}
for an Ising model simulation.

\item {\tt potts\_order} initializes in {\tt common\_potts.f} 
the array of Potts spins {\tt ista} to the value ${\tt izero}=0$.
Alternatively:

\item {\tt potts\_ran} initializes the array of Potts spins
{\tt ista} with uniformly distributed integer random numbers in 
the range $0,\dots,{\tt nqm1}$. The user has to choose a value 
for {\tt istart} in {\tt mc.par}.

\item {\tt potts\_act} calculates, after initialization of the 
spin array {\tt ista}, the initial action value {\tt iact}, which 
is part of the {\tt links} common block of {\tt common\_potts.f}
(its calculations involves the sum over all links). To illustrate 
the lattice structure, we have included an array {\tt ialink}, 
which allows to assign an integer variable to each link, although 
it is not needed for the programs used in this section. If RAM space
matters, {\tt ialink} may be removed, as is done in the common block 
{\tt common\_potts1.f} discussed below.

\item \index{potts\_wght}{\tt potts\_wght} initializes 
{\tt wrat(-2*nd:2*nd,0:mlink+2*nd)}, the weight factor array of 
{\tt common\_potts.f} with Boltzmann (\ref{w_B}) weights.The use 
of an array has speed advantages as array element look ups are 
normally faster than calculations of the exponential function. 
However, this observation is platform dependent and may vary in 
course of time. The reader is advised to perform his own benchmarks 
whenever performance of the code is important. The array {\tt wrat} 
is sized up to {\tt mlink+2*nd} to allow for multicanonical 
simulations discussed in chapter~\ref{MUCA}. For canonical 
Metropolis simulations one can reduce the array size from 
{\tt (-2*nd:2*nd,0:mlink+2*nd)} to {\tt (-2*nd:2*nd)}. This matters
for simulations on large lattices. The user can switch to this version
of the code by including the common statement {\tt common\_potts1.f} 
instead of {\tt common\_potts.f} and using the appropriate ``version~1''
subroutines discussed in section~\ref{potts1_etc}.

\item {\tt razero}(\index{razero.f} initializes the array of its 
arguments to zero for the array range defined by its other two 
arguments. The call in {\tt potts\_init.f} initializes the array 
{\tt ha} of {\tt common\_potts.f} to zero. The array {\tt ha} is 
used to histogram the action variable {\tt iact}.

\end{enumerate}

\begin{table}[tbl]  
\caption{Common block {\tt common\_potts.f}: Arrays and variables used 
for Metropolis simulations of the $q$-state Potts model.}
{\centering\small \begin{verbatim}
      common /sites/  ista(ms)
c iamin,iamax:   action minimum and maximum of the MC sweep.      
      common /links/ nlink, ialink(nd,ms),iact,iamin,iamax
      common /delfct/  idel(0:nqm1,0:nqm1)
c acpt:          counts accepted updates. 
      common /weights/ wrat(-2*nd:2*nd,0:mlink+2*nd),acpt
c ia_min,ia_max: action minimum and maximum for a series of sweeps.      
      common /haction/ ha(0:mlink),ia_min,ia_max
\end{verbatim} } \label{tab_potts.com}
\end{table}

\subsubsection{Updating routines}

\begin{table}[hp]  
\caption{Metropolis MC updating for the Potts models: 
{\tt potts\_met.f} of {\tt ForLib\_par}.}
{\small \centering \begin{verbatim}
      subroutine potts_met ! Copyright Bernd Berg, Nov 16 2000.
C Metropolis updating with sequential spin choice.
      include 'implicit.sta'
      include 'constants.par'
      include 'lat.par'
      include 'potts.par'
      include 'common_lat.f'
      include 'common_potts.f'
      q=nq*ONE
      iamin=iact
      iamax=iact
      do is=1,ms
        istaold=ista(is)
        call ranmar(xr)
        istanew=int(q*xr)
        if(istanew.ne.istaold) then
          idact=0
          do id=1,nd
            ista2=ista(ipf(id,is))
            idact=idact+idel(ista2,istanew)-idel(ista2,istaold)
            ista2=ista(ipb(id,is))
            idact=idact+idel(ista2,istanew)-idel(ista2,istaold)
          end do
          if(idact.ne.0) then
            call ranmar(xr)
            if(xr.lt.wrat(idact,iact)) then
              ista(is)=istanew
              iact=iact+idact
              acpt=acpt+one
              iamin=min(iamin,iact)
              iamax=max(iamax,iact)
            end if
          else
            ista(is)=istanew
            acpt=acpt+one
          end if
        end if
        ha(iact)=ha(iact)+one
      end do
      ia_min=min(ia_min,iamin)
      ia_max=max(ia_max,iamax)
      return
      end
\end{verbatim} }
\label{tab_potts_met} \end{table}

The subroutine 
\begin{equation} \label{potts_met} \index{potts\_met.f}
  {\tt potts\_ met.f}~~~{\rm of}~~~{\tt ForLib\_ par}
\end{equation}
is reproduced in table~\ref{tab_potts_met}. It implements 
the Metropolis procedure~(\ref{Metropolis}) 
using\index{sequential updating|textbf} 
\index{update!sequential|textbf}{\bf sequential updating}, 
{\it i.e.}, sites are updated in their order ${\tt is}=1,2,\dots ,N$. 
\index{random updating|textbf}\index{update!random|textbf}{\bf Random 
updating}, for which sites are chosen with the uniform probability
$1/N$, is performed by the subroutine 
\begin{equation} \label{potts_met_r} \index{potts\_met\_r.f}
  {\tt potts\_ met\_ r.f}~~~{\rm of}~~~{\tt ForLib\_ par}\,,
\end{equation}
which is otherwise identical with ${\tt potts\_ met.f}$. 

A call to ${\tt potts\_ met}$ performs what is called a
\index{sweep}{\bf sweep}: 
\begin{equation} \label{sweep}
  {\rm Each\ spin\ of\ the\ lattice\ is\ updated\ once.} 
\end{equation} 
For a call to {\tt potts\_met\_r} this modifies to:
\begin{equation} \label{sweep_r}
  {\rm Each\ spin\ is\ updated\ once\ in\ the\ average.}
\end{equation} 
The notion sweep is also used for this situation. The Metropolis 
subroutines (\ref{potts_met}) and (\ref{potts_met_r}) increment also 
the {\bf histogram}
\begin{equation} \label{ha}
  {\tt ha(iact)}\ \to\ {\tt ha(iact)+one}~~{\rm at\ each\ update\ 
  attempt\ of}~~{\tt iact}\, .
\end{equation}
Further, the subroutines keep track\index{acceptance rate!Potts models} 
of the {\bf acceptance rate}\index{Potts models!acceptance rate} {\tt 
acpt} and compute the variables
\begin{equation} \label{ia_min_max}
  {\tt iamin} ~~{\rm and}~~ {\tt iamax}
\end{equation}
which are the minimum and the maximum {\tt iact} values visited
during the sweep. The variables ${\tt ia\_min}$ and ${\tt ia\_max}$
track then the {\tt iact} minima and maxima over many sweeps. All 
these variables are part of the common block {\tt common\_potts.f}.

\subsubsection{Start and equilibration}

Under repeated application of$\,$\footnote{The same considerations 
are valid for the other updating procedures.}
${\tt potts\_ met}$ the probability of states will approach the 
Boltzmann distribution~(\ref{P_B}). However, initially we have 
to start with a microstate which may be far off the Boltzmann 
distribution. Far off means, that the Boltzmann probability (at
temperature $T$) for a typical state of the initially generated
distribution can be very, very small. Suppression factors like
$10^{-10000}$ are well possible. As the Markov process proceeds the 
weight of the initial states decreases $\sim 1/n$ and is ultimately 
swallowed by the statistical noise which decreases $\sim 1/\sqrt{n}$. 
Nevertheless we should exclude the initial states from our equilibrium 
statistics as the proportionality factor of the $1/n$ fall-off can
be (much) larger than the one of the statistical noise. In practice 
this means we should allow for a certain number of sweeps
\begin{equation} \label{nequi} \index{equilibration}
  {\tt nequi}~~{\rm to\ equilibrate\ the\ system.}
\end{equation}
Statistical analysis may determine the appropriate value for 
{\tt nequi}, see chapter~\ref{Markov_Error}.
\smallskip

\begin{figure}[tb]
\vspace{5pc}
\centerline{\hbox{ \psfig{figure=figures/2dI_ts.eps,width=10cm} }}
\caption{ \index{time series!$2d$ Ising model}
\index{Ising model!energy!time series}Ising energy per spin
(\ref{e0s}): Two Metropolis time series of 200 sweeps are shown for an 
Ising model on an $80\times 80$ lattice at $\beta=0.4$. Ordered and 
disordered starts are used. The exact mean value $\widehat{e}_{0s}=
-1.10608$, calculated with {\tt ferdinand.f}~(\ref{ferdinand}), is also 
indicated, see assignment~\ref{assign_2dI_ts}.  \label{fig_2dI_ts} }
\end{figure}

Many ways to generate start configurations exist. Two natural and
easy to implement choices are:

\begin{enumerate}

\item Generate a random configuration with 
\index{potts\_ran.f}{\tt potts\_ran.f} of {\tt ForLib}. Such a 
configuration corresponds to $\beta =0$ and is called 
\index{disordered configuration}{\bf disordered}. This defines 
a \index{random start}\index{Potts models!random start}{\bf random} 
or \index{disordered start}\index{Potts models!disordered start}{\bf 
disordered start} of a MC simulation.

\item Generate a completely 
\index{ordered configuration}{\bf ordered} configuration,
for which all Potts spins take on the same $q$-value. The subroutine 
\index{potts\_order.f}{\tt potts\_order.f} of {\tt ForLib} is called 
for an \index{ordered start}\index{Potts models!ordered start}{\bf 
ordered start} of a MC simulation.

\end{enumerate}

Figure~\ref{fig_2dI_ts} shows two Metropolis time series of 200 
sweeps each for an Ising model on an $80\times 80$ lattice at 
$\beta=0.4$, which may be reproduced using the program
\begin{equation} \label{potts_ts}\index{potts\_ts.f}
  {\tt potts\_ ts.f}~~~{\rm of}~~~{\tt ForProg/MC\_Potts}\ .
\end{equation}
The time series starting off at the minimum value $e_{0s}=-2$
of the energy per spin $e_{0s}$ (\ref{e0s}). The other time series
corresponds to a disordered start ($\langle e_{0s}\rangle =0$ is 
the expectation value for the energy per spin at $\beta=0$). After 
a number of sweeps the two time series mix and we may argue that 
configurations may be counted after neglecting the first 100 sweeps 
for reaching equilibrium. A more careful analysis is carried out in 
the next chapter (compare table~\ref{tab_iat}).

To reproduce the Metropolis histograms of figure~\ref{fig_2dI_H_es} 
(the $2d$ Ising model at $\beta=0.2$ and $\beta=0.4$) the program 
{\tt potts\_hist.f}~(\ref{potts_hist}), discussed in 
subsection~\ref{Energy_code} below, can be used, see 
assignment~\ref{assign_2dI_H_es}.

\subsubsection{More updating routines\label{potts1_etc}}

We describe a variant and extensions of the {\tt potts\_met.f}
(\ref{potts_met}) Metropolis updating routine:
\begin{equation} \label{potts_more}\index{potts\_metn.f}
  \index{potts1\_metn.f}\index{potts1\_metn.f} {\tt potts\_ metn.f}, 
  ~~~{\tt potts1\_ met.f}~~~{\rm and}~~~{\tt potts1\_ metn.f}
\end{equation}
of {\tt ForProg/MC\_Potts}. An extension of {\tt potts\_met.f} 
is {\tt potts\_metn.f}: It hits the same spin $n$ times with the 
Metropolis updating procedure. In this way the 
distribution of the spin approaches the Gibbs ensemble created by its 
neighbors. To optimize a simulation, one has to balance the improved 
performance of the spin updates against the additional CPU time needs, 
an issue addressed in chapter~\ref{Markov_Error}.

The routines {\tt potts1\_met.f} and {\tt potts1\_metn.f} agree in 
their function with their counterparts without the ``1'', but they 
include the common block {\tt common\_potts1.f} instead of {\tt 
common\_potts.f}. The common block {\tt common\_potts1.f} reduces 
the size of the weight array {\tt wrat} (table~\ref{tab_potts.com}) 
to its minimum dimension {\tt (-2*nd:2*nd)}. To save further memory 
space, the array {\tt ialink} of {\tt common\_potts.f} is also removed 
in {\tt common\_potts1.f}.  The inclusion of {\tt common\_potts1.f} 
instead of {\tt common\_potts.f} implies that we have to change the 
main program and other routines accordingly. For instance, the properly 
modified version of the program {\tt potts\_hist.f} of {\tt 
ForProg/MC\_Potts} is 
\index{potts1\_hyst.f}\index{potts1\_init.f}{\tt potts1\_hyst.f}, 
which includes {\tt potts1\_init.f} from {\tt ForProg/MC\_Potts} and 
\index{potts1\_wght.f}{\tt potts1\_wght.f} from {\tt ForLib} (instead 
of {\tt potts\_init.f} and {\tt potts\_wght.f}). Similarly,
for the use of {\tt potts1\_metn.f} a properly modified time series
program, {\tt pottsn1\_ts.f}, is provided. 
It is left as assignment~\ref{assign_potts1_ts} to the reader to 
write a similar modification of the updating routine
{\tt potts\_met\_r.f}~(\ref{potts_met_r}).

\subsection{Heatbath code\label{heat_bath_code}}

The\index{Potts models!heatbath algorithm}
\index{heatbath algorithm!Potts models}{heatbath algorithm} 
chooses a state $q_i$ directly with the local Boltzmann distribution 
defined by its nearest neighbors. Compare our discussion for the $O(3)$ 
$\sigma$-model in section~\ref{O3_heatbath}. For the Potts model the 
state $q_i$ takes on one of the values $0,\dots,q-1$. With all other 
states set it determines the value of the energy function~(\ref{E}). 
For this we denote this energy by $E(q_i)$ and the Boltzmann 
probabilities are
\begin{equation} \label{Potts_B_weights}
  P_B(q_i)\ =\ {\rm const}\ e^{-\beta\, E(q_i)}
\end{equation}
where the constant is determined by the normalization condition
\begin{equation} \label{Potts_B_normalization}
  \sum_{q_i=0}^{q-1} P_B(q_i)\ =\ 1\ .
\end{equation}        
In equation~(\ref{Potts_B_weights}) we can define $E(q_i)$ to be just 
the energy of the interaction of $q_i$ with its nearest neighbors, 
as the other contributions to the total energy are the same for all 
$q_i$ and absorbed into the overall constant. So $E(q_i)$ depend only 
on the the nearest neighbors of $q_i$ are partitioned into values 
$0,\dots,q-1$. For low values of $q$ and the dimension $d$ the most 
efficient implementation of the heatbath algorithm is to tabulate all 
possibilities\footnote{Such tables are discussed in chapter~\ref{EDS}.}. 
Here we follow a different path, because we want a generic code, which 
works for arbitrary values of $q$ and $d$.

\begin{table}[hp]  
\caption{Heatbath MC updating for the Potts model in an
external magnetic field $H$.}  
{\centering \footnotesize \begin{verbatim}
      subroutine potts_mchb ! Copyright Bernd Berg, Oct 23 2002.
C MC Heat Bath updating with sequential spin choice.
C Measurements: the action ha and the magnetization hm histograms.
C Faster when the magnetization measurements are eliminated!
      include 'implicit.sta'
      include 'constants.par'
      include 'lat.par'
      include 'potts.par'
      include '../ForLib_par/common_lat.f'
      include '../ForLib_par/common_potts_hb.f'
      dimension prob_hb(-1:nqm1),iact_a(0:nqm1)
      iamin=iact
      iamax=iact
      prob_hb(-1)=zero
      do is=1,ms
        call ranmar(xr)               
        do iq=0,nqm1
          iact0=0 ! Calculates the action contribution for the state iq.
          do id=1,nd
            iact0=iact0+idel(iq,ista(ipf(id,is)))
     &                 +idel(iq,ista(ipb(id,is)))
          end do
          prob_hb(iq)=whb_tab(iact0,iq)+prob_hb(iq-1)
          iact_a(iq)=iact0                      ! Defines the array iact_a.
        end do
        do iq=0,nqm2
          ptest=prob_hb(iq)/prob_hb(nqm1) ! Heat bath probabilities.
        if(xr.lt.ptest) go to 1
        end do
        iq=nqm1
1       iqold=ista(is)
        if(iq.ne.iqold) then
          achg=achg+one
          iqnew=iq
          ista(is)=iqnew
          iact=iact+iact_a(iqnew)-iact_a(iqold)
          iamin=min(iamin,iact)
          iamax=max(iamax,iact)
          nstate(iqold)=nstate(iqold)-1 ! For magnetization measurements.
          nstate(iqnew)=nstate(iqnew)+1 ! For magnetization measurements.
        end if
        ha(iact)=ha(iact)+one
        do iq=0,nqm1          ! Do loop for magnetization measurements.
          hm(nstate(iq),iq)=hm(nstate(iq),iq)+one
        end do
      end do
      ia_min=min(ia_min,iamin)
      ia_max=max(ia_max,iamax)
      return
      end
\end{verbatim} } \label{tab_potts_mchb}
\end{table}

The heatbath updating routine is
\begin{equation} \label{potts_mchb}\index{potts\_mchb.f}
  {\tt potts\_ mchb.f}~~~{\rm of}~~~{\tt ForLib\_ par}
\end{equation}
and listed in table~\ref{tab_potts_mchb}. It has to be used with 
the initialization routine 
\begin{equation} \label{potts_inithb}\index{potts\_inithb.f}
  {\tt potts\_ inithb.f}~~~{\rm of}~~~{\tt ForLib\_ par}\ .
\end{equation}
Both routines include the common block {\tt common\_potts\_hb.f}.
The major change of {\tt potts\_inithb.f}, compared to {\tt 
potts\_init.f} of table~\ref{tab_potts_init}, is that heatbath 
weights {\tt WHB\_TAB(0:2*ND,0:NQM1)} are initialized by
{\tt potts\_wghb.f}, which is shown in table~\ref{tab_potts_wghb}. The 
weight table includes an external field {\tt H0}, whose value is set 
in the {\tt mc.par} parameter file (table~\ref{tab_mc.par}). The use 
of the heatbath and Metropolis weights is quite different as a 
comparison of the updating  subroutines {\tt potts\_met.f} 
(table~\ref{tab_potts_met}) and {\tt potts\_mchb.f} 
(table~\ref{tab_potts_mchb}) reveals. Whereas the Metropolis weights 
determine directly the updating decision, the heatbath weight table 
is used to calculate the cumulative distribution function for the 
heatbath probabilities $P_B(q_i)$ of equation~(\ref{Potts_B_weights}). 
This is the array {\tt prob\_hb(-1:nqm1)} of the {\tt potts\_mchb.f} 
code, whose final values are
\begin{equation} \label{prob_hb}
  {\tt prob\_ hb(iq)} = \sum_{q_i=0}^{\tt iq} P_B(q_i)\,,\
  {\tt iq}=0,\dots,{\tt nqm1}\ .
\end{equation}
The normalization condition~(\ref{Potts_B_normalization}) implies
{\tt prob\_hb(nqm1)=1}. Comparison of these cumulative probabilities
with the uniform random number {\tt xr} yields the heatbath update
$q_i\to q'_i$. Note that in the heatbath procedure the original value
$q_i$ does not influence the selection of $q'_i$. 

The variable {\tt achg} counts the updates, which lead to a change of 
the spin variable at hand. It is similar to {\tt acpt}, which monitors 
the accepted changes of the Metropolis updating.

\begin{table}[tbl]  
\caption{\index{potts\_wghb.f}{\tt potts\_wghb.f} of ForLib, which 
generates the weight table used for our Potts model heatbath updating.}
{\centering \small \begin{verbatim}
        SUBROUTINE POTTS_WGHB(WHB_TAB,BETA,H0,ND,NQM1)
C Copyright, Bernd Berg, Nov 22 2000.
C Table for Potts model (with magnetic field) heatbath update.
        include 'implicit.sta' 
        include 'constants.par' 
        DIMENSION WHB_TAB(0:2*ND,0:NQM1)
        DO IQ=0,NQM1
          H=ZERO
          IF(IQ.EQ.0) H=H0
          DO IACT=0,2*ND
            WHB_TAB(IACT,IQ)=EXP(TWO*(BETA*IACT+H))
          END DO
        END DO
        RETURN
        END
\end{verbatim} }
\label{tab_potts_wghb}
\end{table}

Our heatbath program accumulates also\index{Potts models!magnetization}
\index{histograms!Potts magnetization}{\bf magnetization histograms}.  
For the Potts model each of the $q$ distinct states defines a direction 
of the Potts magnetization. Hence, $q$ histograms are kept in one 
array. See in table~\ref{tab_potts_mchb} the lines commented:
\begin{equation} \label{Potss_magnetization}
 {\tt For\ magnetization\ measurements.} 
\end{equation}

\subsection{Timing and time series comparison of the routines}

We have to balance the updating performance of our subroutines versus 
their CPU time consumptions. The results will depend on the hardware 
and compilers used so that evaluations have to be done by the final 
user. Here an illustrative example is discussed. A Pentium~II 400~MHz 
processor under Linux with the gnu $g77$ Fortran compiler version 
2.91.66 was used. You get information about your CPU by typing
\index{CPUinfoc}{\tt cat /proc/cpuinfo}. The version number of the 
\index{g77!version number|textbf}$g77$ Fortran compiler under 
\index{Linux}{Linux} is identified with the command
\begin{equation} \label{g77_version}
  {\tt g77\ -v}\ .
\end{equation}
Compilations were done with the ${\tt -O}$ optimization.

Initialization and 10$\,$000 sweeps were performed for the 10-state 
Potts model on an $80\times 80$ lattice. The number of sweeps is chosen 
relatively large, to make the initialization time (in which we are not 
really interested) negligible. The different updating routines consumed 
CPU time as follows (see (\ref{Linux_time}) and (\ref{tcsh_CPU_time}) 
for CPU time measurements under Linux):
\begin{eqnarray} \nonumber \index{CPU time!measurements} 
  &39.1~{\rm s}& ({\tt potts\_ met.f}),~~~~
   83.5~{\rm s}\ ({\rm 2-hit}\ {\tt potts\_ metn.f})\\ \label{CPU_time}
  &{\rm and}&~ 170.2~{\rm s}\ ({\tt potts\_ mchb.f})\ .
\end{eqnarray}
\index{simulations!heatbath}\index{simulations!Metropolis}Of
\index{simulations!Potts models}the three routines used the heatbath 
updating has normally the best dynamics, but one heatbath updating 
needs more than four times the CPU time of one 1-hit Metropolis 
updating. Multi-hit Metropolis can usually only beat the performance 
of the heatbath updating when its CPU time consumption is less than 
that of a heatbath routine. For {\tt nhit=2} the multi-hit Metropolis 
CPU time consumption\footnote{With the gnu g77 compiler our $n$-hit 
Metropolis updating is found to scale in CPU time with the surprisingly 
bad factor of $n$ and this did not improve when the neighbor spins were
first stored in a small intermediate array to convince the compiler to
keep these numbers in cache during the multi-hits.}
is approximately half of that of our heatbath updating. It is
assignment~\ref{assign_CPU_time} to repeat the CPU timing on your 
computer.

\begin{figure}[tb]
\vspace{5pc}
\centerline{\hbox{ \psfig{figure=figures/2d10q_ts.eps,width=10cm} }}
\caption{\index{time series!$2d$ 10-state Potts model}
\index{Potts models!energy!time series}$q=10$ Potts model time 
series of 200 sweeps on an $80\times 80$ lattice at $\beta=0.62$. 
Measurements of the action variable (\ref{actl}) are plotted for 
ordered and disordered starts. Updating was performed with the 
following routines: {\tt potts\_met.f}, {\tt nhit=2 potts\_metn.f} 
and {\tt potts\_mchb.f}, see assignment~\ref{assign_2d10q_ts}. 
\label{fig_2d10q_ts} }
\end{figure}

How can we compare the performance of the different updating procedures? 
As we see in chapter~\ref{Markov_Error}, this requires to measure 
autocorrelation times. At the present state we are content with a look 
at the different time series and find that the performance can indeed 
be quite distinct. For the 10-state model on an $80\times 80$ lattice 
at $\beta=0.62$ we compare in figure~\ref{fig_2d10q_ts} times series 
of 200 sweeps each. The simulation temperature is in the disordered 
phase, but not too far away from the phase transition point of
table~\ref{tab_potts_exact}. Ordered and 
disordered starts were performed. While for the disordered 
starts there is relatively little visible difference in the approach 
of the simulations to the equilibrium process, this is not the case
for the ordered starts. The 1-hit Metropolis process has considerable 
difficulties to find its way to the equilibrium chain. The obvious 
interpretation is that this is due to the finite lattice remnant of 
the first order phase transition. The convergence speed is improved 
for the 2-hit Metropolis process, while the heatbath algorithm 
performs clearly best. The question, whether the improvement due 
to the heatbath updating does for the equilibrium process still
overcompensate for the increased CPU time consumption~(\ref{CPU_time}) 
is addressed in the next chapter (see table~\ref{tab_iathb}).

\subsection{Assignments for section~\ref{Potts}\label{assign_Potts}}

In the following we use always Marsaglia random number with their 
default seed (table~\ref{tab_mar} of chapter~\ref{Sampling}). This 
makes the results easily reproducible. As long as the same random 
numbers are used in different settings the actually produced 
correlations are presumably negligible. Nevertheless, to be on 
the statistically save side, the random numbers should differ 
from simulation to simulation in real research investigations.

\begin{enumerate}

{\item Reproduce the energy per spin time series of 
figure~\ref{fig_2dI_ts}. Initialize the parameter files 
to the requested values and run the program ${\tt potts\_ ts.f}$ 
from {\tt ForProg/MC\_Potts}.  Use {\tt ferdinand.f}~(\ref{ferdinand}) 
to calculate the exact mean value. Repeat the simulation using 
{\tt potts\_ts\_r.f}. \label{assign_2dI_ts} }

{\item Use the program {\tt potts\_hist.f} of {\tt ForProg/MC\_Potts}
to reproduce the $\beta = 0.2$ and $\beta = 0.4$ histograms of 
figure~\ref{fig_2dI_H_es}. Start each run with the default seed of 
the Marsaglia random number generator and allow 10$\,$000 sweeps to 
reach equilibrium. The error bars are with respect to 20 bins of 
5$\,$000 sweeps each. Combine the results with those 
of assignment~\ref{assign_b0} of subsection~\ref{assign_sampling} 
to reproduce figure~\ref{fig_2dI_H_es}. \label{assign_2dI_H_es} } 

{\item Modify the updating routine {\tt potts\_met\_r.f} 
(\ref{potts_met_r}) to {\tt potts1\_met\_r.f}, to use {\tt 
common\_potts1.f} instead of {\tt common\_potts.f}. Check that the 
{\tt potts\_ts\_r.f} time series of assignment~\ref{assign_2dI_ts} 
stays unchanged. \label{assign_potts1_ts} }

{\item Run {\tt potts\_hist.f} at $\beta=0.2$ with 1$\,$000 sweeps for 
reaching equilibrium and $10\times 1\,000$ sweeps for data production. 
Repeat the same calculation with {\tt potts1\_hist.f} and compare the 
generated data files (under Linux with {\tt diff fln1.d fln2.d}). 
Perform the same comparisons for the {\tt nhit} Metropolis routines 
{\tt potts\_metn.f} and {\tt potts1\_metn.f} with {\tt nhit=2}. 
\label{assign_potts1_hist} }

{\item Repeat the CPU timing~(\ref{CPU_time}) of {\tt potts\_met.f}, 
2-hit {\tt potts\_metn.f} and {\tt potts\_mchb.f} on your computer 
using your favorite Fortran compiler. Try different optimization 
levels. Further, include the routines {\tt potts1\_met.f} and 
{\tt potts1\_metn.f} in the competition. \label{assign_CPU_time} }

{\item Use {\tt potts\_ts.f} of {\tt ForProg/MC\_Potts} to
generate the Metropolis time series of figure~\ref{fig_2d10q_ts}
for the 10-state Potts model at $\beta =0.62$ on an $80\times 80$
lattice. Repeat the calculation using (a)~$n$-hit Metropolis updating 
with $n=2$ (the program {\tt pottsn\_ts.f} of {\tt ForProg/MC\_Potts}) 
and (b)~heatbath updating (the program {\tt pottshb\_ts.f} of 
{\tt ForProg/MC\_Potts}). \label{assign_2d10q_ts} }

\end{enumerate}


\section{Potts Model Monte Carlo Simulations \label{Potts_MC}}

\index{Potts models!energy!reference values}Using
\index{data analysis}our MC subroutines, we present some high 
precision calculations of the average internal energy per Potts spin.
The purpose is to test the correctness of the routines versus exact 
results, numerical results from the literature, and against one another. 
This is immediately hit by a snag: The \index{MC data}{MC data} of 
a Markov chain are correlated and we do not yet know how to calculate 
error bars reliably in this situation. At our present level of 
expertise we resort on a heuristic way out of the dilemma: We use 
the binning procedure of chapter~\ref{Binning}. The bins become 
(almost) independent when the binsize is chosen large enough. How 
large is large enough? This is a topic of chapter~\ref{Markov_Error}. 
Here we just choose a number which is (very) large compared to the 
scale on which convergence of the time series ala 
figure~\ref{fig_2dI_ts} is found and hope that the correlation 
between bins are then negligible. Due to the central limit
theorem the distribution of the bins will be approximately Gaussian 
and we can use Gaussian error analysis tools. 
\index{Kolmogorov test}Kolmogorov tests allow to check the Gaussian 
assumption for consistency, see assignment~\ref{assign_kolmo_gau}.

\subsection{Data production, analysis code and energy references
\label{Energy_code} }

\begin{table}[hp]
\caption{The data production program {\tt potts\_hist.f} performs 
 {\tt nequi} equilibration sweeps followed by {\tt nrpt} bins of 
 {\tt nmeas} measurements. For each bin the assembled histogram 
 {\tt ha} is written to an unformatted data file, which is named 
 to reflect the parameters of the run.}
{\centering \small \begin{verbatim}
      program potts_hist ! Copyright, Bernd Berg, Dec 20 2000.
C MC Metropolis production run. Energy (action variable) histograms.
      include '../../ForLib/implicit.sta'
      include '../../ForLib/constants.par'
      parameter(n0=0,lpri=.true.)
      character cd*1,cq*2,cl*3
      include 'lat.par'
      include 'mc.par'
      include 'potts.par'
      include '../../ForLib_par/common_lat.f'
      include '../../ForLib_par/common_potts.f'
      include 'lat.dat'
C
      call potts_init(lpri) ! Initialize Potts Metropolis MC:
C
      do iequi=1,nequi      ! Sweeps for reaching equilibrium.
        call potts_met
      end do
      acpt=acpt/(one*nequi)/(one*ms)
C
      write(cd,'(I1.1)') nd
      write(cq,'(I2.2)') nq
      write(cl,'(I3.3)') nla(1)
      open(iud1,file="p"//cd//"d"//cq//"q"//cl//".D",
     &     form="unformatted",status="unknown")
      write(iud1) beta,nd,nq,nla,mlink,nequi,nrpt,nmeas
      do irpt=1,nrpt ! nrpt repititions of nmeas measurement sweeps.
        call write_progress(iud2,"irpt,iact,acpt:",irpt,iact,acpt)
        acpt=zero
        call razero(ha,n0,mlink)
        do imeas=1,nmeas
          call potts_met
        end do
        acpt=acpt/(one*nmeas)/(one*ms)
        write(iud1) ha,acpt,irpt
      end do
      close(iud1)
C
      stop
      end
\end{verbatim} } \label{tab_potts_hist}
\end{table}

In table~\ref{tab_potts_hist} we list a simple
\index{data production}{data production} program
\begin{equation} \label{potts_hist}\index{potts\_hist.f}
  {\tt potts\_ hist.f}~~{\rm of}~~{\tt ForProg/MC\_ Potts}\ .
\end{equation}
It is only suitable for production runs, which are considerably shorter 
than the medium up-time of you computer system (for large scale 
production runs see chapter~\ref{LSP}).  

After performing {\tt nequi} sweeps for reaching equilibrium, the 
program creates a \index{filename creation}{\bf filename} which 
includes the dimension $d$, the number of states $q$, and the lattice 
size $n^1$ in $x^1$ direction~(\ref{site_coordinate}). For the 
frequently encountered symmetric case $n^1$ determines the entire 
lattice. Using the same programming technique, one can extend or alter 
this name list to fit specific purposes. The format of the data file is 
chosen to be {\tt "unformatted"}, which is the efficient way of handling 
large data sets under Fortran. (For small data sets readable {\tt 
"formatted"} files have, besides being readable, the advantage of being 
easily exchangeable between different computer architectures.) For the 
filename extension we try to be consistnet with the following hierarchy:

\begin{description}

\item{\tt *.d} for files which can be (relatively) easily reproduced.
In MC simulations these are typically files of short test runs or
files created by data analysis programs. Under Linux all theses files 
can then be removed with the instruction {\tt rm *.d}.

\item{\tt *.D} for files whose reproduction may take considerable CPU
time. In a MC simulations of spin systems these are typically files
from production runs. If they are now longer needed, they can all be 
deleted using {\tt rm *.D} under Linux. Be sure that you really want 
this.

\item{\tt *.dat} for files which normally should never be deleted,
{\it e.g.}, {\tt lat.dat}, which defines the the lattice size.

\end{description}

After opening the data file, the program proceeds with writing general 
information 
into the first block of the data file. Subsequently, in the {\tt nrpt}
do loop, it writes one action variable histogram {\tt ha} together with
other information per block. Each block corresponds to {\tt nmeas} 
updates with measurements (the parameter {\tt nsw} of {\tt mc.par} is 
only used in later time series production programs and allows then to 
take measurements less often than after every sweep). In the production 
part the program calls also the subroutine 
\begin{equation} \label{write_progress}\index{write\_progess.f}
  {\tt write\_ progess.f}~~{\rm of}~~{\tt ForLib} 
\end{equation}
which is listed in table~\ref{tab_write_progress}. This subroutine opens 
the formatted data file {\tt progress.d}, writes information about the 
status of the calculation on the data file and closes it immediately 
afterwards. The contents of the file {\tt progress.d} can be looked up 
\index{Fortran!runtime supervision}\index{monitor runtime progress}at 
runtime and thus allows to monitor the simulation. This is of importance 
for long runs, which are normally performed in the background. Without 
opening and closing the file before and after each write, the computer 
may store the data in some buffer, instead of writing them immediately 
into the file.

\begin{table}[th]  
\caption{The subroutine {\tt write\_progress.f} of {\tt ForLib} allows
the user to to monitor the progress of a (long) MC calculation at
run time by looking up the content of the file {\tt progress.d}.}
{\centering \small \begin{verbatim}
      subroutine write_progress(iud,ctxt,irpt,iact,acpt)
C Copyright, Bernd Berg, Dec 14, 2000.
C Write subroutine to monitor progress.
      include 'implicit.sta'
      character*(*) ctxt
      open(iud,file="progress.d",form="formatted",status="unknown")
      rewind iud
      write(iud,'(A33,2I12,G15.5)') ctxt,irpt,iact,acpt
      close(iud)
      return
      end 
\end{verbatim} } \label{tab_write_progress}
\end{table}

\subsubsection{Data analysis}

\begin{table}[hp]  
\caption{Using {\tt read} instructions the subroutine {\tt 
read\_steb0.f} of {\tt ForLib} adapts the error bar subroutine 
{\tt steb0.f}~(\ref{STEB0}) to handle sets of possibly large arrays.}
{\centering \small \begin{verbatim}
      SUBROUTINE READ_STEB0(nrpt,iud,norm,nlink,ha,ham,hae,actl)
C Bernd Berg Jan 20 2002.
C Reads nrpt times the array ha() from unit iud and calculates the
C       histogram mean values ham() and their error bars hae().
C Further, the nrpt action variable averages are calculated.
C norm=0:  no normalization of the ha histogram, otherwise normalization.
      include 'implicit.sta'
      include 'constants.par'
      dimension ha(0:nlink),ham(0:nlink),hae(0:nlink),actl(nrpt)
C
      call razero(ham,izero,nlink)
      factor=one
      do irpt=1,nrpt
        read(iud) ha
	if(norm.ne.0) then ! Normalization of ha() to hasum=nlink.
          hasum=zero
          do ilink=0,nlink
            hasum=hasum+ha(ilink)
          end do
          factor=nlink/hasum
          do ilink=0,nlink ! Normalization of ha() to hasum=nlink.
            ha(ilink)=factor*ha(ilink)
          end do
        end if
        actl(irpt)=potts_actl(nlink,ha) ! Mean action variable.
        do ilink=0,nlink
          ham(ilink)=ham(ilink)+ha(ilink)
        end do
      end do
      do ilink=0,nlink
        ham(ilink)=ham(ilink)/nrpt
      end do
C
      do irpt=1,nrpt
        backspace iud
      end do
C
      call razero(hae,izero,nlink)
      do irpt=1,nrpt
        read(iud) ha
        do ilink=0,nlink
          hae(ilink)=hae(ilink)+(factor*ha(ilink)-ham(ilink))**2
        end do
      end do
      do ilink=0,nlink
        hae(ilink)=sqrt(hae(ilink)/(nrpt*(nrpt-1)))
      end do
C
      return
      end
\end{verbatim} } \label{tab_read_steb0}
\end{table}

For the\index{histograms!data analysis} 
\index{Potts models!histogram analysis}{\bf analysis of the histogram 
data} the program 
\begin{equation} \label{ana_hist.f}\index{ana\_hist.f}
  {\tt ana\_ hist.f}~~{\rm of}~~{\tt ForProg/MC\_ Potts} 
\end{equation}
is provided. After initializing the data filename and reading its
first block, the analysis program relies to a large extent on the 
subroutine 
\begin{equation} \label{read_steb0.f}
  {\tt read\_ steb0.f} ~~{\rm of}~~ {\tt ForLib} 
\end{equation}
which is reproduced in table~\ref{tab_read_steb0}. This subroutine adapts 
our basic statistical error bar routine {\tt steb0.f}~(\ref{STEB0}) to 
the handling of arrays, which are possibly too big to keep {\tt nrpt} of 
them into the RAM. Read instructions are used to get away with just three 
equally sized arrays, {\tt ha}, {\tt ham} and {\tt hae}. Optionally, we 
can normalize the incoming histograms {\tt ha}, so that for each one of 
them the sum of its entries becomes {\tt nlink}. When setting the input 
argument ${\tt norm\ne 0}$ the {\tt ha} arrays are normalized, for 
${\tt norm=0}$ they are used as read.

\begin{table}[th]  
\caption{Fortran function \index{potts\_actl.f}{\tt potts\_actl.f} 
of {\tt Forlib}: The average action per link is calculated from a 
histogram of the {\tt iact} variable~(\ref{iact}). }
% A variant of this routine, which is for some applications much 
% faster, is {\tt potts\_act2l.f} of {\tt ForLib}.}
{\centering \small \begin{verbatim}
      FUNCTION potts_actl(nlink,ha) ! Bernd Berg, Dec 13, 2000.
C Calculation of the mean action variable per link from its histogram. 
C Use hist_aver for arbitrary histograms (no link normalization).
      include 'implicit.sta'
      dimension ha(0:nlink)
      potts_actl=0
      hsum=0
      do ilink=0,nlink
        potts_actl=potts_actl+ilink*ha(ilink)
        hsum=hsum+ha(ilink)
      end do
      potts_actl=potts_actl/(hsum*nlink) ! iact per link.
      return
      end
\end{verbatim} }
\label{tab_potts_actl}
\end{table}
 
For each of the {\tt nrpt} arrays {\tt ha} the subroutine 
{\tt read\_steb0.f} calculates also the average action per link
and returns these {\tt nrpt} numbers in the array {\tt act}. This
calculation relies on calls to the Fortran function {\tt potts\_actl.f} 
of {\tt ForLib} listed in table~\ref{tab_potts_actl}. It is better to 
calculate the average action from the histograms than by using the {\tt 
iact} action values after each sweep, because the updating routine 
generates \index{fluctuations!small}{small fluctuations} during a 
sweep which smoothen the histograms. 

Our analysis program returns the average action per link {\tt actlm},
and also converts it to the average energy per spin $\overline{e}_{0s}$, 
which is according to equations (\ref{E0_iact}), (\ref{actl}), (\ref{actlm})
\begin{equation} \label{act_e0s}
  \overline{e}_{0s} = {2\,d\over q} - 2\,d\, {\tt actlm}\ .
\end{equation}
As explained, to achieve consistency with Ising model conventions of
the literature our definition of $\overline{e}_{0s}$ differs from 
the one used in most of the Potts model literature.

\subsubsection{$2d$ Ising model}

We begin with the $2d$ Ising model and test against the exact results 
of our program {\tt ferdinand.f}~(\ref{ferdinand}). As for the time 
series from an $80^2$ lattice in figure~\ref{fig_2dI_ts}, we use 
the value $\beta=0.4$, but simulate on an $20^2$ lattice as for the 
histograms in figure~\ref{fig_2dI_H_es}. The convergence of the Markov 
chain to its equilibrium process is expected to be even faster than in 
figure~\ref{fig_2dI_ts}, where the first 200 sweeps are depicted. Now 
we use a statistics of 10$\,$000 sweeps for reaching equilibrium and 
assume that this is an overkill.  The statistics for measurement is 
chosen 32 times larger: 64 bins of 5$\,$000 sweeps each. The binsize 
of 5$\,$000 ($\gg 200$) is argued to be large enough so that the 
distribution of the 64 data points is indeed practically Gaussian, 
and correlations between the bins are (hopefully) small enough to be 
neglected. The number of 64 bins ensures that according to the student 
distribution (tables~\ref{tab_stud1} and~\ref{tab_stud2}) confidence 
levels of error bars are Gaussian on the two $\sigma$ level. With this 
statistics we find (assignment~\ref{assign_2dI_es}) 
\begin{equation} \label{2dI_es}
  \overline{e}_{0s} = -1.1172\ (14)~~{\rm (Metropolis)} 
  ~~~{\rm versus}~~~\widehat{e}_s = -1.117834~~{\rm (exact)}\ .
\end{equation}
Performing the Gaussian difference test of chapter~\ref{Gauss_Dif} 
gives $Q=0.66$, a perfectly admissible value. 

\subsubsection{$2d$ 4-state Potts models}

Our next example is the $2d$ 4-state Potts model. We compare with 
a high statistics result on a $120\times 120$ lattice: From table~1 
of \cite{CaTa99} we pick $\beta = 1.06722/2$ (their $\beta^{\rm Potts}
 = 1.06722$) and repeat the calculation with our updating routine 
{\tt potts\_met.f} (Caselle {\it et al}. used Swendsen-Wang cluster
updating, which we introduce in chapter~\ref{AUS}). Our statistics is
400$\,$000 sweeps for reaching equilibrium and $64 \times 200\,000$ sweeps 
for measurement. The results are
\begin{eqnarray} \nonumber
  &{\tt actlm}& =\ 0.646057\ (24)~~{\rm (Metropolis)}~~~{\rm versus}\\
  &{\tt actlm}& =\ 0.64604\ (6)~~[{\rm Caselle}\ {\it et\ al}.\ (1999)]
  \label{2d04q_actlm}
\end{eqnarray}
and $Q=0.79$ follows for the Gaussian difference test, see 
assignment~\ref{assign_2d04q}. 

\subsubsection{$2d$ 10-state Potts models}

For the $2d$ 10-state Potts model at $\beta =0.62$ we test our Metropolis
routine {\tt potts\_met.f} versus our heatbath routine {\tt potts\_mchb.f} 
on a small $20\times 20$ lattice. For the heatbath
updating we use the same statistics as previously for the $2d$ Ising
model~(\ref{2dI_es}), 10$\,$000~sweeps for reaching equilibrium and
$64\times 5\,000$ sweeps for measurements. For the Metropolis
updating we increase these numbers by a factor of four to 40$\,$000 sweeps
for reaching equilibrium and $64\times 20\,000$ sweeps for measurements.
This increase is done, because we expect the performance of Metropolis 
updating for the 10-state model to be far worse than for the 2-state 
model. At low temperature in the ordered phase the likelihood 
to propose the most probable (aligned) Potts state is 1/2 for 
the 2-state model, but only 1/10 for the 10-state model. Although 
$\beta=0.62$ is not yet in the ordered phase, it is sufficiently close 
to it, so that this effect is expected to be of relevance. The results 
of our simulations are
\begin{eqnarray} \nonumber
  &{\tt actlm}& =\ 0.321772\ (75)~~{\rm (Metropolis)}~~~{\rm versus}\\
  &{\tt actlm}& =\ 0.321661\ (70)~~{\rm (heatbath)}\label{2d10q_e0s}
\end{eqnarray}
and $Q=0.28$ for the Gaussian difference test, see 
assignment~\ref{assign_2d10q_e0s}. The heatbath 
\index{data production}{data production} program 
is a variant of {\tt potts\_hist.f}, called 
\index{potts\_hihb.f}{\tt potts\_hihb.f}. It stores besides the 
energy\index{magnetization}\index{histograms!Potts magnetization}
\index{Potts models!magnetization}{\bf also magnetization histograms}.
For the energy histograms the analysis program {\tt ana\_hist.f} 
(\ref{ana_hist.f}) still works.

\subsubsection{$3d$ Ising model}

In \cite{AlBe90} a precise numerical estimate of the internal energy is 
given for the Ising model on a $14^3$ lattice at $\beta = 0.22165$. This 
$\beta$ value is close to the infinite volume critical point, which in 
the same reference\index{Ising model!$3d$!critical temperature}
\index{critical temperature}\index{temperature!critical}is estimated 
to be $\beta_c = 0.22157\ (3)$. We repeat their calculation with a 
statistics of 20$\,$000 sweeps for reaching equilibrium and $64\times 
10\,000$ sweeps with measurements. The results compare as 
follows\footnote{\cite{AlBe90} gives $E_A= -715.97\ (6)$ where
$E_A$ is related to our definition (\ref{E}) by $E=4 E_a$. Therefore, 
$\overline{e}_{0s}=4\overline{E}_a/N = -2863.88/14^2 = -1.04369$. }
\begin{eqnarray} \nonumber \index{Ising model!$3d$}
  &\overline{e}_{0s}& =\ -1.0427\ (14)~~ {\rm (Metropolis\ here)} 
  ~~~{\rm versus}\\ \label{3dI_es}
  &\overline{e}_{0s}& =\ -1.04369\ (9)~~[{\rm Alves}\ {\it et\ al}.\ 
  (1990)]
\end{eqnarray}
and the Gaussian difference test gives $Q=0.48$, see 
assignment~\ref{assign_3dI_es}.

\subsubsection{$4d$ Ising model with non-zero magnetic field}

Finally, let us test the heatbath code for a non-zero value of the 
magnetic field. Equation (\ref{beta0_H}) allows for simple tests 
at $\beta =0$. This equation does not depend on the dimension $d$, 
and we use the opportunity to run the code in $d=4$. In 
assignment~\ref{assign_beta0_H} we choose the external magnetic 
field to correspond to $\langle \delta_{q_i,0}\rangle = 0.5$ 
for the 10-state Potts model and run $32\times 10$ sweeps with 
measurements on an $10^4$ lattice. The numerical result is
\begin{equation}\label{4d10q_M1}\index{Potts models!$4d$ magnetization} 
  \langle \delta_{q_i,1} \rangle = 0.49961\, (71)\,,
\end{equation}
which yields $Q=0.58$ when compared to $0.5$ by means of the Gaussian
difference test.

\subsection{Study of a first order phase transition: 
            $3d$ 3-state Potts model}

To illustrate investigations of first order phase transitions, we 
perform simulations of the $3d$ 3-state Potts model using the 1-hit 
Metropolis algorithm. Typically one starts with small lattices for
which it is relatively easy to scan the temperature dependence. In this 
way one may locate a \index{double peak histogram}{\bf double peak} 
structure in the enery histogram on a $8^3$ lattice at $\beta =0.2745$. 
For a statistics of $64\times 16\,000$ sweeps after discarding another 
$16\,000$ sweeps for reaching equilibrium the histogram of the action 
variable is depicted in figure~\ref{fig_3d3q_h}, where we have chosen
\begin{equation} \label{hist_norm} 
   \sum_{{\tt iact}=0}^{\tt nlink} H({\tt iact}) = {\tt nlink}
\end{equation}
as normalization (see assignment~\ref{assign_3d3q_h}).

\begin{figure}[tb]
\vspace{5pc}
\centerline{\hbox{ \psfig{figure=figures/3d3q08l_eh.eps,width=10cm} }}
\caption{Histogram of the action per link variable {\tt actl} for the 
$3d$ 3-state Potts model on a $8^3$ lattice at $\beta=0.2745$ and its
reweighting to equal heights (see assignment~\ref{assign_3d3q_h}). 
\label{fig_3d3q_h} }
\end{figure}

Such a \index{histograms!double peak}{double peak} structure, is 
typically obtained when systems with first order transitions 
are simulated on finite lattices in the neighborhood of so called 
pseudo-transition temperatures. The distance between the peaks is 
a finite volume estimator of the\index{latent heat!estimator} 
\index{Potts models!$3d$ 3-state latent heat}{latent heat} 
of\index{pseudo-transition temperatures!equal heights}
the transition.\index{pseudo-transition temperatures|textbf}
\index{temperature!pseudo-transition|textbf}{\bf Pseudo-transition
temperatures} are finite lattice temperature definitions, which 
converge with increasing system size towards the infinite volume 
transition temperature. For first order transitions a definition of 
a pseudo-transition temperature is the value which gives 
\index{equal heights pseudo-transition temperature}{\bf equal heights}
to both peaks of the histogram. An occasionally used alternative is
\index{pseudo-transition temperatures!equal weights}the {\bf equal 
weights} defintion, for which the integrals on both side agree 
\cite{BoKa92}. In the following we use the equal heights definition
as it is technically easier to handle.

The reweighting to equal heights is done iteratively by the analysis 
program {\tt ana\_ehgts.f} of {\tt ForProg/MC\_Potts}. It relies on
two routines: The subroutine \index{\tt double\_peak.f}
\begin{equation} \label{double_peak} 
   {\tt double\_ peak.f~~~{\rm of}~~~ForLib}
\end{equation}
is used to locate the two maxima and the in-between minimum of 
the double-peak histogram. From this we obtain an estimate of the 
equal heights $\beta$ values, which is used to calculate the
reweighted histogram with our subroutine {\tt potts\_rwght.f}
(\ref{potts_rwght}). For our $8^3$ lattice the reweighted histogram
is also shown in figure~\ref{fig_3d3q_h}. The result is fed into
(\ref{double_peak}) until convergence is achieved. The jackknife 
estimate of the equal-heigts $\beta$ values is $\beta_c(8) = 0.275295
\, (91)$. It is collected in table~\ref{tab_3d3q} together with 
estimates of the mean action-variable value {\tt actlm}, the distance 
between the two maxima $\triangle\,{\tt actm}$ and ratio of the 
equal-heights maxima divided by the in-between minimum 
$H_{\max}/H_{\min}$.

\begin{table}[tbl] \begin{center}
\tbl{Estimates of equal-heights pseudo-transition temperatures for the 
3-state Potts model on $L^3$ lattices together with other observables 
of the equal-heights distributions (assignment~\ref{assign_3d3q_ehgts}). 
The statistics used was $n$ sweeps for reaching equilibrium and $64\,n$ 
sweeps with measurements, where $n$ is given in the second column.}
% \medskip
{ \begin{tabular}{|c|c|c|c|c|c|}
\hline
$L$& $n$     & $\beta_c$   & {\tt actlm} 
                                  & $\triangle\,{\tt actlm}$ &
                                         $H_{\max}/H_{\min}$ \\ \hline
 8 &$16\,000$& 0.275295 (91)& 0.60036 (72)& 0.092 (19) & 1.105 (13)\\
12 &$32\,000$& 0.275276 (42)& 0.58498 (45)& 0.0847 (15)& 1.227 (21)\\
16 &$64\,000$& 0.275286 (21)& 0.57813 (31)& 0.0747 (23)& 1.495 (27)\\
20 &$80\,000$& 0.275278 (15)& 0.57465 (27)& 0.0688 (09)& 1.837 (46)\\
24 &$96\,000$& 0.275292 (23)& 0.57241 (61)& 0.0631 (08)& 2.32 (14) \\
\hline
\end{tabular} \label{tab_3d3q} } \end{center} \end{table}

\begin{figure}[tb]
\vspace{5pc}
\centerline{\hbox{ \psfig{figure=figures/3d3q_ehgts.eps,width=10cm} }}
\caption{Normalized histograms of the action variable {\tt actl} for 
the $3d$ 3-state Potts model on $L^3$ lattices reweighted to 
equal heights at $\beta_c(L)$ of table~\ref{tab_3d3q}, see 
assignment~\ref{assign_3d3q_ehgts}. \label{fig_3d3q_ehgts} }
\end{figure}

The equal-heights value $\beta =0.275295$ of the $8^3$ lattice is used 
as inverse simulation temperature for the next larger lattice of size 
$12^3$. We create a statistics of $32\,000$ sweeps for reaching 
equilibrium and $64\times 32\,000$ sweeps with measurements. Repeating 
the previous analysis procedure yields the $L=12$ equal heights value
$\beta_c(12) = 0.275276\, (42)$. Iterating the procedure towards larger 
lattices, always using the $\beta_c(L)$ estimate of the previous for
the next MC run gives the data collected in table~\ref{tab_3d3q}, 
where the statistics is given in the second column 
(assignment~\ref{assign_3d3q_ehgts}). However, with the methods 
of the present chapter we cannot decide whether the statistics is 
really sufficient. The reader is refered to later chapters, in 
particular to~\ref{MC_consistent} and~\ref{MUCA}. In 
figure~\ref{fig_3d3q_ehgts} we plot of all reweighted histograms 
together. 

We combine all our pseudo-transition temperature estimates to an 
infinite volume extrapolation by performing a two-parameter fit
with a coorectgion of order 1/Volume:
\begin{equation} \label{3d3q_beta_fit}
  \beta_c (L) = \beta_c + a\,L^{-3}\ .
\end{equation}
The data show little lattice size dependence and our infinite
volume estimate is $\beta_c = 0.275283\,(14)$ with a goodness of 
fit $Q=0.95$ (assignment~\ref{assign_3d3q_beta}).
Accurate estimates of the infinite volume transition temperature
exist\index{Potts models!$3d$ 3-state transition temperature} in
the literature, {\it e.g.}, $\beta_c = 0.2752616\,(48)$ 
\cite{AlBe91}\footnote{The $\beta$ value given in table~v of 
\cite{AlBe91} is $\beta_A = 0.3669727$ and related to our notation 
by $\beta = 3\beta_A/4$.} or $\beta_c = 0.2752825\,(50)$ \cite{JaVi97}
with a second order endpoint at a non-zero magnetic field 
$h_c=0.0003875\,(50)$ for magnetization $M_c=h_c/\beta_c(h_c)$
and $\beta_c(h_c)=0.274690\,(10)$ \cite{KaSt00}.

The Gaussian difference test gives $Q=0.15$ when comparing our fit 
result with $\beta_c = 0.2752616\,(48)$ and $Q=0.97$ when comparing 
it with $\beta_c = 0.2752825\,(50)$. For reference we also give the 
mean action value on our $24^3$ lattice at the critical temperature 
estimate $\beta=0.275278$ (slightly off the equal-heights reweighted 
value $\beta=0.275292$ of table~\ref{tab_3d3q}). It is
\begin{equation} \label{3d3q_actm} 
  {\tt actlm} = 0.57233\ (70) \ .
\end{equation}
Due to the double peak structure of the histogram the error bar is 
relatively large (as well as those of {\tt actlm} in 
table~\ref{tab_3d3q}). One may be interested to check whether the 
64 bins on which the error bar of (\ref{3d3q_actm}) relies are really 
Gaussian distributed. In assignment~\ref{assign_kolmo_gau} a Kolmogorov 
test\index{Kolmogorov test}\index{Kolmogorov test!one-sided} is 
performed, which shows that the statistics of the bins is indeed 
large enough so that the central limit works: $Q_1=0.33$ for the 
positive and $Q_2=0.70$ for the negative deviation of the one-sided
Kolmogorov test~(\ref{kolm1}).


\subsection{Assignments for section~\ref{Potts_MC}\label{assign_MC}}

As in the previous assignment section we use always Marsaglia random 
number with their default seed (see the discussion there).

\begin{enumerate}

{\item Use {\tt potts\_met.f} to reproduce the $2d$ Ising model 
$\overline{e}_{0s}$ Metropolis result of equation~(\ref{2dI_es}). 
Suitable production and data analysis programs are {\tt potts\_hist.f} 
and {\tt ana\_hist.f} of {\tt ForProg/MC\_Potts}. Note that the exact 
value from {\tt ferdinand.f} differs from that of 
assignment~\ref{assign_2dI_ts} in section~\ref{Potts} due to finite 
size effects. Is the difference detectable within the accuracy of 
the MC run? \label{assign_2dI_es} }

{\item Use {\tt potts\_met.f} to reproduce the $2d$ 4-state Potts 
model {\tt actlm} Metropolis result of equation~(\ref{2d04q_actlm}).
Make an estimate of the CPU time needed (compare 
assignment~\ref{assign_CPU_time})
before starting this simulation. \label{assign_2d04q} }

{\item Use {\tt potts\_met.f} versus {\tt potts\_mchb.f} to 
reproduce the $2d$ 10-state Potts model {\tt actlm} results of 
equation~(\ref{2d10q_e0s}). \label{assign_2d10q_e0s} }

{\item Use {\tt potts\_met.f} to reproduce our $3d$ Ising model 
$\overline{e}_{0s}$ Metropolis result of equation~(\ref{3dI_es}). 
\label{assign_3dI_es} }

{\item Consider the 10-state Potts model at $\beta =0$ and determine 
the magnetic field $H$ so that $\langle \delta_{q_i,0}\rangle=0.5$ 
holds according to equation~(\ref{beta0_H}). Subsequently, use this
$H$ value and {\tt potts\_hihb.f} of {\tt ForProg/MC\_Potts} to 
calculate $\langle\delta_{q_i,0}\rangle$ numerically on a $10^4$ 
lattice. Use 32 bins, each of 10 sweeps with measurements. 
No equilibrium sweeps are necessary. Why? Estimate 
$\langle\delta_{q_i,0}\rangle$ from the data (the appropriate 
analysis program is {\tt ana\_mag.f} of {\tt ForProg/MC\_Potts}) 
and perform the Gaussian difference test with respect to $0.5$. 
Calculate also the exact mean action per link and perform the 
Gaussian difference test with respect to its estimate from the 
data. \label{assign_beta0_H} }

% ------------ first-order:

{\item Use {\tt potts\_met.f} to reproduce the $3d$ 3-state Potts 
model histogram of figure~\ref{fig_3d3q_h} by running at $\beta =
0.2745$.  Reweight the histogram to equal heights as also plotted 
in figure~\ref{fig_3d3q_h} and find the $\beta$ value reported for
$L=8$ in table~\ref{tab_3d3q}.  \label{assign_3d3q_h} }

{\item Iterate the calculation of the previous assignment to lattices
of size $12^3,\,16^3,\,20^3$ and $24^3$. Simulate always at the 
equal-heights temperature estimate from the previous lattice and 
generate the statistics given in the third row of table~\ref{tab_3d3q}. 
The pseudocritical $\beta_c(L)$ values of this table as well as 
figure~\ref{fig_3d3q_ehgts} should be obtained. To ensure easy 
reproducibility all simulations are started with our default seed
of the random number generator, what should be avoided in serious
applications. \label{assign_3d3q_ehgts} }

{\item Perform the two-parameter fit (\ref{3d3q_beta_fit}) for the 
$\beta_c(L)$ values of table~\ref{tab_3d3q} to obtain the $L\to\infty$ 
estimate given in the text.  \label{assign_3d3q_beta} }

{\item Check for assignments~\ref{assign_2dI_es}, \ref{assign_2d04q}, 
\ref{assign_2d10q_e0s}, \ref{assign_3dI_es} and~\ref{assign_3d3q_ehgts} 
$(L=24)$, whether the assumption that the binned data are 
Gaussian distributed passes the upper and lower one-sided
\index{Kolmogorov test!one-sided}{Kolmogorov tests} (\ref{kolm1}). 
This is best done by running the program\index{\tt ana\_kol.f} 
{\tt ana\_kol.f} of {\tt ForProg/MC\_Potts} in the folders of 
assignments~\ref{assign_2dI_es}, \ref{assign_2d04q}, 
\ref{assign_2d10q_e0s}, \ref{assign_3dI_es} and~\ref{assign_3d3q_ehgts},
where the parameter and data files ought to be located. For the 
interpretation compare assignment~\ref{assign_kolm1a} in 
section~\ref{sec_assign_kolmo}.  \label{assign_kolmo_gau} }

\end{enumerate}

\section{Specific Heat, Reweighting, Error Bars and Jackknife\label{Crwght}}

With $\widehat{E}=\langle E\rangle$ the\index{specific heat|textbf}
{\bf specific heat} is defined by
\begin{equation} \label{specific_heat} 
  \index{canonical estimator!specific heat}
  C\ =\ {d\,\widehat{E}\over d\,T}\ =\ \beta^2\, \left( \,
  \langle E^2 \rangle - \langle E \rangle^2 \right)\ .
\end{equation}
The last equal sign follows by working out the temperature derivative 
in (\ref{O}) and is known as \index{fluctuation-dissipation}{\bf 
fluctuation-dissipation theorem}. It is often used to estimate the 
specific heat from equilibrium simulations without relying on numerical 
derivatives. An equivalent formulation is
\begin{equation} \label{fluc_dis} 
  C\ =\ \beta^2\, \left\langle \left( E - \widehat{E}\right)^2
                 \right\rangle
\end{equation}
as is easily shown by expanding $\left(E-\widehat{E}\right)^2$ and 
working out the expectation values. Defining
\begin{equation} \label{act2lm2}
  {\tt act2lm} = {\overline{\tt iactm^2}\over {\tt mlink}}\,,
\end{equation}
and using the relation with the energy (\ref{act_e0s}), we have instead 
of (\ref{specific_heat}) a notation, which is close to the computer code:
\begin{equation} \label{C_Ising}
  C\ =\ \frac{\beta^2\,N\,d^2}{n}\,\sum_{i=1}^n
  \left({\tt act2l}_i-{\tt actlm}^2\right)\,,
\end{equation}
where the sum is over all measurements in the times series. 
Translating (\ref{fluc_dis}) similarly gives
\begin{equation} \label{act_fluc_dis}
  C\ =\ \frac{\beta^2\,N\,d^2}{n}\,
  \sum_{i=1}^n ({\tt actl}_i-{\tt actlm})^2\,.
\end{equation}
When energy histograms are available, (\ref{act_fluc_dis}) is
calculated by the equations used in\index{\tt potts\_mu2.f}
{\tt potts\_mu2.f} of {\tt ForLib}.

In the limit of an infinite statistics specific heat estimates
from all these equation agree. But with a finite statistics a number 
of problems emerge. In the simple binning approach (\ref{NBIN}), where 
{\tt nrpt} of the production run defines the number of blocks, one 
may want to use for $\widehat{E}$ estimators $\overline{E}_i$ which 
are constructed from the histograms of the blocks. With $N_b$ the 
number of data in each bin and $H_i$ the energy histogram of block $i$:
\begin{equation} \label{EiEstimate}
  \overline{E}_i = \frac{1}{N_b} \sum_{j \in {\rm block}(i)} E_j\
  = \frac{\sum_E E\,H_i(E)}{\sum_E H_i(E)}\ .
\end{equation}
Estimates of $\overline{C}_i$ from either equation (\ref{specific_heat}) 
or (\ref{fluc_dis}) agree then and Gaussian error bars are expected due 
to the binning. However, for a not so good statistics as bias towards 
too small $\overline{C}_i$ values is occurs, because $E_j$ and 
$\overline{E}_i$ come from the same block in the estimate
\begin{eqnarray} \label{Cest1}
  \overline{C}_i = \frac{\beta^2}{N_b} \sum_{j\in{\rm block}(i)} 
  \left( E_j^{\,2} - \overline{E}_i^{\,2} \right) & = & 
  \frac{\beta^2\, \sum_E \left(E^{\,2}-\overline{E}_i^{\,2}\right)
  \,H_i(E)}{\sum_E H_i(E)} \\ \label{Cest2}
   = \frac{\beta^2}{N_b} \sum_{j\in{\rm block}(i)} 
  \left(E_j-\overline{E}_i\right)^2& =& \frac{\beta^2\,
  \sum_E (E-\overline{E}_i)^2\,H_i(E)}{\sum_E H_i(E)}\ .
\end{eqnarray}
The $\overline{E}_i$ estimators are certainly inferior to the estimate 
$\overline{E}$, which relies on the entire statistics. Therefore, one
may consider to use $\overline{E}$ instead of $\overline{E}_i$ in
equations (\ref{Cest1}) and (\ref{Cest2}). However, then one does not 
know anymore how to calculate the error bar of $\overline{C}$ as the 
$\overline{C}_i$ estimates would rely on overlapping data. The 
situation gets even worse when we include reweighting 
(\ref{potts_rwght}), as this non-linear procedure implies
that the estimators (\ref{Cest1}) and (\ref{Cest2}) will in
general differ. These difficulties are overcome by the jackknife
method of chapter~\ref{jackknife}. When histograms can be used
the fast way to create jackknife bins is to sum first the entire
statistics ({\tt HSUM} in the code):
\begin{equation} 
  H(E) = \sum_{i=1}^{N_b} H_i(E)\ .
\end{equation}
Subsequently jackknife histograms (superscript $J$) are defined by
\begin{equation} 
  H_i^J(E) = H(E) - H_i(E)
\end{equation}
and jackknife estimates $\overline{C}_i^J$ are obtained by using 
$H_i^J(E)$ instead of $H_i(E)$ in equations (\ref{EiEstimate}),
(\ref{Cest1}) and (\ref{Cest2}).

In the following we compare results obtained by several variants of 
simple binning to those of the jackknife method. To pin down the
subtleties of the statistical analysis we first consider the rather
extrem case of reweighting of random sampling on a lattice which is 
small enough to allow for the calculation of its partition function 
at all temperatures. Afterwards we perform a canonical simulations 
of the 2D Ising model on larger lattices and reweight the specific 
heat to a small neighbourhood of the simulation temperature. 

\subsection{Reweighting of random sampling}

\begin{figure}[tb]
\vspace{5pc}
\centerline{\hbox{ \psfig{figure=figures/2dI_4x4.eps,width=10cm} }}
\caption{Energy histogram from random sampling of the $2d$ Ising model 
on a $4\times 4$ lattice together with is reweighting to $\beta_0=0.2$
and $\beta_0=0.4$ (see assignment~\ref{assign_4x4b0}).  
\label{fig_2dI_4x4} }
\end{figure}

To illustrate the jackknife method, we consider again reweighting of 
random sampling. On a $4\times 4$ lattice the Ising model has $2^{16}
=65\,536$ states. This is small enough to generate them all by random 
sampling. For the large statistics of assignment~\ref{assign_4x4b0} 
figure~\ref{fig_2dI_4x4} depicts on a log scale the histogram of random 
sampling together with its reweighting to $\beta_0=0.2$ and $\beta_0=
0.4$. As all energies are covered we are able to calculate the specific 
heat at all temperatures. If the statistics is very large, simple 
binning gives reliable results. However, problems are encountered
when the statistics covers some states barely. 

\begin{table}[th]
\tbl{Calculation of the specific heat by reweighting of random sampling
for the Ising model on a $4\times 4$ lattice 
(assignments~\ref{assign_4x4b0short} to~\ref{assign_4x4b0repeat}).} 
% \medskip \centering
{ \begin{tabular}{||c|c|c|c|c|c|c||}               \hline
Statistics: & \multicolumn{2}{|c|}{(a) $32\times 30\,000$}   
   & \multicolumn{2}{|c|}{$32\times 3\,10^8$}   
   & \multicolumn{2}{|c||}{$10\,000$ times (a)} \\ \hline
& Estimate   & Q  & Estimate       & Q  & Estimate        & Q \\ \hline
\multicolumn{2}{||c|}{ Simple binning: }&&&&& \\ \hline
$C$ & $x.xxx\,(xx)$ &x.xx&$x.xxxxx\,(xx)$&x.xx&$-$& $-$\\ \hline
\multicolumn{3}{||c|}{with mean from entire statistics:}&&&& \\ \hline
$C$ via (\ref{C_Ising})&
$ x.xx\,(xx)$ &x.xx&$x.xxxx\,(xx)$ &x.xx&$-$& $-$\\ \hline
$C$ via (\ref{act_fluc_dis})&
$x.xxxx\,(xx)$&x.xx&$x.xxxxx\,(xx)$&x.xx&$x.xxxxx\,(xx)$& 0.00 \\ \hline
\multicolumn{2}{||c|}{Jacknife binning:}&&&&& \\ \hline
$C$           &$x.xxxx\,(xx)$ &x.xx&$x.xxxxx\,(xx)$&x.xx&
                                    $x.xxxxx\,(xx)$&x.xx       \\ \hline
bias corrected&$x.xxxx\,(xx)$ &x.xx&$x.xxxxx\,(xx)$&x.xx&
                                    $x.xxxxx\,(xx)$&x.xx       \\ \hline
\end{tabular} \label{tab_4x4b0} }
\end{table}

In assignment~\ref{assign_4x4b0short} we generate the rather small 
statistics of $32\times 30\,000$ sweeps. Resulting estimates of the 
specific heat are shown in column two of table~\ref{tab_4x4b0}. First
simple binning is used as in equations (\ref{Cest1}) and (\ref{Cest2}),
which are identical as the $E_j$ and $\overline{E}_i$ rely on the same 
measurement for each bin. The estimate is unsatisfactory, because
comparison with the exact result
$$C=0.812515$$ 
from {\tt ferdinand.f} (\ref{ferdinand}) gives $Q=0.00$ for the 
Gaussian difference test (\ref{gau_dif}). We try to improve on this by 
using the average energy $\overline{E}$ of the entire statistics in 
equations (\ref{Cest1}) and (\ref{Cest2}), which are then non longer 
identical. In the table the computer notation of these equation, 
(\ref{C_Ising}) and (\ref{act_fluc_dis}), is quoted. The observation 
is that (\ref{C_Ising}) fails entirely, while the estimate from 
(\ref{act_fluc_dis}) is good. The bad estimate comes because the error 
of the fluctuations of the second moment used in (\ref{C_Ising}) are 
large, whereas the reduced second moment used in (\ref{act_fluc_dis}) 
enforces a positive result with smaller fluctuations. In either case 
$\overline{E}$ is entered without statistical error, because its 
fluctuations are small compared to those of $E_j$. 

To demonstrate that the bad estimator does still converge towards
the exact result we increase in assignment~\ref{assign_4x4b0long}
our statistics by a factor of~$1\,000$. As shown in column~4 of the 
table, the value has become reasonable when compared with the exact 
result, but the error bar is five times larger than for the other 
estimators, which agree very well with one another. Independently 
of any increase of the statistics the fluctuations of the bad estimator 
stay to large, because they are calculated with respect to zero 
instead of $\overline{E}$.

\begin{figure}[tb]
\vspace{5pc}
\centerline{\hbox{ \psfig{figure=figures/2dI_qdfC.eps,width=10cm} }}
\caption{Empirical peaked distribution function of three estimators
of the specific heat of the Ising model as discussed in the text.
\label{fig_2dI_qdfC} }
\end{figure}

Using jackknife bins, the estimates (\ref{C_Ising}) and 
(\ref{act_fluc_dis}) ar identical again, and the goodness of fit is 
satisfactory. However, it remains inconclusive whether there is an 
advantage compared to the simple binning estimates (\ref{act_fluc_dis}) 
in column two. Also, there is no improvement when we use equation 
(\ref{bias_corrected_estimator}) to correct for the bias (we keep the 
error bars of the jackknife estimators as the bias turns out to be 
small compared to them, so that a next level jackknife analysis is not 
necessary). The results are so close to one another that we cannot 
draw conclusions about the quality of the estimators from a single 
low statistics run. To get to the bottom of this we repeat $10\,000$ 
times in assignment~\ref{assign_4x4b0repeat} the low statistics 
simulation of assignment~\ref{assign_4x4b0short}. From these simulation
averages and error bars with respect to the $10\,000$ repetitions 
are given in column~6 of table~\ref{tab_4x4b0}. {\bf Only the bias 
corrected jackknife estimator averages to an acceptable value.}
The uncorrected jackknife estimator comes in second and the good 
simple binning estimate third. Figure~\ref{fig_2dI_qdfC} depicts
the empirical peaked distribution functions (\ref{qdf}) of these
estimates. As there is no disadvantage in applying jackknife 
instead of simple binning: {\bf The jackknife method should in 
essence always be used.}

\subsection{Reweighting of a canonical simulation}

In real applications reweighting is mostly of importance when one
wants to locate maxima of observables such as the specific heat. In 
assignment~\ref{assign_2dI_rwght} canomnical MC simulations are 
performed at for the $2d$ Ising model at $\beta=0.43$, which is
relatively close to $\beta_c$. After reweighting bias corrected
jackknife estimates of the specific heat are calculated and 
plotted in figure~\ref{fig_2dI_rwght}. 

\begin{figure}[tb]
\vspace{5pc}
\centerline{\hbox{ \psfig{figure=figures/2dI_rwght.eps,width=10cm} }}
\caption{Specific heat of the Ising model on $L\times L$ lattices. 
Exact result produced by the program {\tt ferdinand.f} (\ref{ferdinand}) 
versus reweighting of simulations from $\beta = 0.43$. 
\label{fig_2dI_rwght} }
\end{figure}

Maxima of the specific heat are clearly within the reweighting range 
and the peak is more pronounced on the larger lattice. On the other 
hand the reweighted values start to deviate on the larger lattice from 
the exact curve at the boundary of the chosen $\beta$ range. That
reflects the shrinking of the reweighting range with increasing lattice 
size. A fit of the $L$ depence of the location of the peak allows,
however, to determine a good value for the simulation temperature on 
the next larger lattice and so one. FSS investigation can be performed 
by iterating such a process. After introducing more advanced simulation
algorithms FSS techniques are discussed in chapter~\ref{FSS}.

\subsection{Assignments for section~\ref{Crwght}\label{assign_Crwght}}

As in the previous assignment sections we use always Marsaglia random 
number with their default seed. 

% ------------ 4x4 reweighting:

\begin{enumerate}

{\item Repeat the reweighting to $\beta_0=0.2$ of 
assignment~\ref{assign_b0} of section~\ref{sec_2dI} on a $4\times 4$ 
lattice. Reweight also to $\beta_0=0.4$. Increase the statistics by
a factor of ten.  \label{assign_4x4b0} }

{ \item Use random sampling to create a statistics of $32 \times 
30\,000$ sweeps for the $2d$ Ising model on on a $4\times 4$ lattice. 
Reweight to $\beta_0=0.4$ and perform the following estimates:
{\tt actlm} using simple binning, {\tt actlm} from all data, $C$ 
using simple binning, $C$ using (\ref{C_Ising}) with simple binning 
and {\tt actlm} from all data, $C$ using (\ref{act_fluc_dis}) with 
simple binning and {\tt actlm} from all data, $C$ from jackknife
binning, and the bias corrected $C$.  \label{assign_4x4b0short} }

{ \item Increase the statistics of the previous assignment
to $32\times 3\,10^8$ and show that all estimators give now 
reasonable results. But are their error bars consistent
in the sense of the F-test. If not, which of the error bars
should one trust? \label{assign_4x4b0long} }

{ \item Repeat $10\,000$ times the analysis of 
assignment~\ref{assign_4x4b0short} for the good estimator of simple 
binning, the jacknife and the biased improved jackknife estimator. 
Compute mean values and their standard deviations with respect to 
the $10\,000$ events in each case. \label{assign_4x4b0repeat} }

{\item Use the Metropolis algorithm to simulate the Ising model on 
$10\times 10$ and $20\times 20$ lattices at $\beta=0.43$. Reweight
to the $\beta_0$ range $[0.38:0.48]$ and plot bias corrected
jackknife estimates of the specific heat together with the
exact result for $C(\beta)$. For the $20\times 20$ lattice:
Compare the estimates of {\tt actlm} from simple binning and
from jackknife binning (a)~at $\beta=0.43$ (no reweighting) 
and (b)~at $\beta=0.40$ (reweighting).
\label{assign_2dI_rwght} }

\end{enumerate}


\section{Continuous Systems \label{MC_continuous} }

MC simulations of \index{continuous systems|textbf}systems with
\index{systems!continuous|textbf}continuous variables face additional 
difficulties. Histogramming of the observables leads to discretization 
errors and it is better to keep the time series. However, keeping 
too many data can slow down the code and exhaust available storage 
facilities. One may also become concerned that the discretization 
of the Marsaglia random numbers will matter. Any resolution smaller 
than $\epsilon_{\rm Marsaglia} = 2^{-24}=1/16,777,216$ needs special 
attention and one may want to employ {\tt ranmar2.f} (\ref{ranmar2}). 
Another complication is that rounding errors may limit the bit by bit 
reproducibility of our numerical results, see below.

\index{O(n) models}\index{simulations!O(n) models}To give an 
impression about dealing with continuous systems, we consider 
simulations of $O(n),\, n\ge 2$ spin systems in $d\ge 1$ dimensions. 
The systems are defined by equations (\ref{O3_Z}) to (\ref{O3_energy}), 
only that we allow now\index{energy!O(n) models}
for\index{continuous systems!O(n) models|textbf}
\index{systems!continuous!O(n) models|textbf}{\bf $\bf n$-component 
spins}, instead of just for $n=3$ ($n=1$ is the Ising model). 
Of particular physical interest are the cases $n=3$, which 
we introduced in section~\ref{O3_heatbath}, and the 
$n=2$\index{O(2) model: see XY model}
\index{planar rotator: see YX model}{\bf planar rotator}, 
also called \index{XY model}{\bf $XY$ model}. According 
to\index{Kosterlitz-Thouless transition}
\index{phase transition!Kosterlitz-Thouless}\cite{KoTo73,Ko74} 
there exists a low temperature phase of the 
$d=2$ planar rotator which is characterized
by a power law decay in the pair correlation function, modified by the
presence of pairs of vortices. The vortex pairs unbind at the critical
temperature to create a high temperature phase where the correlations
decay exponentially. In this way the transition does not violate the
\index{Mermin-Wagner theorem}{\it theorem} by \cite{MeWa66}, which 
forbids spontaneous breaking of a continuous symmetry in two dimensions. 
\cite{Co73} explains the rationale behind the Mermin-Wagner theorem 
as follows: The breaking of a continuous symmetry leads to massless 
\index{Goldstone bosons}{Goldstone bosons}, but this is not possible 
in two dimensions where the massless scalar field fails to exist.

It is straightforward to write an $O(n)$ Metropolis code for generic 
$n\ge 2$ values, see the next subsection. The code illustrates the 
simplicity of the Metropolis approach, but is too slow for serious 
applications. For the planar rotator and the $n=3$ Heisenberg model
specialized subroutines are provided in sections~\ref{XY_met} 
and~\ref{O3_hb}: For the $XY$ model an improved Metropolis codes 
and for the Heisenberg model the heatbath algorithm of 
section~\ref{O3_heatbath}. It has to be remarked that in both cases 
the cluster code of chapter~\ref{AUS} is highly efficient and the 
reader should consider it before starting any large scale simulations
of $O(n)$ or similar models.

In physics MC simulations of continuous systems span a wide range of 
applications from particle and astrophysics to chemical and biophysics. 
In particle physics MC code for SU(n) \index{lattice gauge theory}{\it 
lattice gauge theories}, in particular for the physically important 
cases SU(2) and SU(3), can be written along similar lines as those 
discussed in the present section \cite{Cr80,CaMa82}. See 
chapter~\ref{LGT}. {\it All atom}\index{all atom} Metropolis simulations 
of \index{peptides}\index{peptides!all atom simulations}{\it peptides}
\index{proteins}\index{proteins!all atom simulations}(small proteins) 
are an example from biophysics. Although molecular dynamics simulations
are the main stream approach in this field (see chapter~\ref{Conclusions}
for some references) the MCMC approach has received some attention in 
recent years, {\it e.g.}~\cite{BKS92,MiSu01}. Several MC simulation 
packages are available. Due to its (relative) modularity, the Fortran 
code \index{SMMP}{SMMP}~\cite{Ei01} ({\it Simple Molecular Mechanics for 
Proteins}) is suitable for beginners, although there are limitation,
for instance due to the not yet implemented atom endgroups.

Over the wide range of applications our 
data analysis methods are rather universally applicable to MC 
simulations. 

{\bf Rounding errors:} Arithmetic\index{rounding errors|textbf} 
\index{error!rounding|textbf}in the floating point representation 
is not exact. See section~1.2 of {\it Numerical Recipes} \cite{PrFl86} 
for a short discussion. Still one will get identical results on 
distinct platforms as long as all operation are covered by the
\index{IEEE}IEEE (\ref{IEEE}) standard. However, there is no IEEE 
standard for the basic functions, like the exponential function or 
logarithms. The rounding errors become compiler dependent and limit 
the reproducibility of our numerical values to staying within a certain 
computing platform. Only agreement in the statistical sense remains 
then between results obtained on distinct platforms. Whenever this 
happens, we shall give the values obtained under the Cygwin (special)
version 3.4.4 of the g77 Fortran compiler. Even the {version number} 
(\ref{g77_version})\index{g77!version number} of the g77 compiler 
matters. 

\subsection{Simple Metropolis code for the $O(n)$ spin models}
\label{On_met}

\begin{table}[ht]  % \begin{table}[tbl]  
\caption{Metropolis updating for $O(n)$ spin models:
{\tt on\_met.f} of {\tt ForLib\_par}. Continued in the
next table.}
{\small \centering \begin{verbatim}
      subroutine On_met ! Bernd Berg, Feb 13 2002.
C O(n) model: Sequential Metropolis updating (inefficient code).
      include 'implicit.sta'
      include 'constants.par'
      include 'lat.par'
      include 'mc.par'
      include 'on.par'
      parameter(halfmeps=half-one/two**25)
      include '../ForLib_par/common_lat.f'
      include '../ForLib_par/common_on.f'
      dimension stanew(n)
      amin=act
      amax=act
\end{verbatim} } \label{tab1_on_met}
\end{table}

\begin{table}[ht]  % \begin{table}[tbl]  
\caption{Continuation: Metropolis updating for $O(n)$ spin models.}
{\small \centering \begin{verbatim}
      do is=1,ms
1       sum=zero ! Propose new spin stanew().
          do i=1,n
            stanew(i)=two*(rmafun()-halfmeps)
            sum=sum+stanew(i)**2
          end do
        if(sum.gt.one) go to 1 ! stanew() inside the unit sphere.
        dact=zero
        fn=one/sqrt(sum)
        do id=1,nd
          do i=1,n
            dact=dact+sta(i,ipf(id,is))*(fn*stanew(i)-sta(i,is))
            dact=dact+sta(i,ipb(id,is))*(fn*stanew(i)-sta(i,is))
          end do
        end do
        if(rmafun().lt.exp(beta*dact)) then
          do i=1,n
            sta(i,is)=fn*stanew(i)
          end do
          act=act+dact
          acpt=acpt+one
          amin=min(amin,act)
          amax=max(amax,act)
        end if
      end do
      a_min=min(a_min,amin)
      a_max=max(a_max,amax)
      return
      end
\end{verbatim} } \label{tab2_on_met}
\end{table}

Tables~\ref{tab1_on_met} and~\ref{tab2_on_met} depicts our Metropolis 
code\index{MC algorithms!O(n) models} 
\index{Metropolis algorithm!O(n) models}\index{on\_met.f}{\tt on\_met.f} 
for simulating the $O(n)$ model for general ${\tt n}\ge 1$ in arbitrary
dimensions ${\tt nd}\ge 1$. The value of {\tt n} is set in the 
parameter file {\tt on.par}. All other parameters are set in the files 
{\tt lat.par} and {\tt mc.par}, just as in the Potts model simulations. 
The $n$ spin components are stored in the array {\tt sta(n,ms)} of the 
common block {\tt common\_on.f}, which is reproduced in 
table~\ref{tab_common_on}. This common block is, as always, kept in 
{\tt ForLib\_par}. It replaces the {\tt common\_potts.f} block used 
by the Potts model routines, while the {\tt common\_lat.f} common block 
stays the same. The array {\tt alink(nd,ms)} of {\tt common\_on.f}
is not used by the present routine and can be eliminated when storage 
space matters. Should one be interested in measuring energy--energy 
correlations, one would store the link energies in this array.

\begin{table}[tbl]  
\caption{The common block {\tt common\_on.f} of {\tt ForLib\_par} for 
simulations of the $O(n)$ models.}
{\centering \small \begin{verbatim}
      common /sites/  sta(n,ms)
c amin,amax:    action  act  minimum and maximum of the MC sweep.      
c acpt:         counts accepted updates. 
      common /links/ alink(nd,ms),act,amin,amax,acpt,nlink
c tsa: time series action array.
c a_min,a_max:  action minimum and maximum for a series of sweeps.      
      common /ts_action/ tsa(nmeas),a_min,a_max
\end{verbatim} } \label{tab_common_on}
\end{table}

Instead of using histograms we rely now on 
the\index{time series!analysis}\index{time series!O(n) models} 
{\bf time series} and introduce appropriate data analysis techniques. 
The time series approach avoids discretization errors, which are 
encountered when the histogramming approach is used for continuous 
observables. During a production run the time series data are 
intermediately kept in the array {\tt tsa(nmeas)} of {\tt common\_on.f}. 
Whenever the array is filled up, it is written out to a file. To 
store time series data intermediately in reasonably sized arrays 
and to write the arrays out is much faster than writing many small 
data sets separately to disk. To allow to jump some sweeps the
parameter {\tt nsw} of {\tt lat.par} is introduced into these
production programs.

The main program
\begin{equation} \label{on_ts}\index{on\_ts.f}
  {\tt on\_ ts.f}~~~{\rm of}~~~{\tt ForProg/MC\_ On}
\end{equation}
performs equilibration as well as \index{data production}{data 
production}. In its simulation run this program writes the array 
{\tt tsa} {\tt nrpt} times out. {\it I.e.}, a time series of length
\begin{equation} \label{On_time_series}
  {\tt nrpt}\times {\tt nmeas} 
\end{equation}
is kept on file. For initialization the program {\tt on\_ts.f} calls 
the subroutine \index{\tt on\_init.f}{\tt on\_init.f}, which invokes 
\index{on\_order.f}{\tt on\_order.f} or \index{on\_ran.f}{\tt on\_ran.f} 
and \index{on\_act.f}{\tt on\_act.f}.

$O(n)$ spins are proposed in \index{on\_met.f}{\tt on\_met.f} 
between the line with the label~{\tt 1} and the {\tt go~to~1} 
statement. They are stored in the array {\tt stanew(n)} and calculated 
as follows: Spin components are 
generated which are uniformly distributed in the $n$-dimensional 
$[-1,+1]$ box of volume $B_n=2^n$. Spins which do not fit into 
the unit sphere are rejected, until a spin is found which fits. This 
spin is subsequently normalized to one. In the code the products
{\tt fn*stanew(i)} are the components of the normalized spin. This 
method of proposing spins works well for small $n$,
but deteriorates quickly for large $n$. The volume $V_n$ of the unit 
sphere is $V_n=\pi^{n/2}/\Gamma (1+n/2)$, see appendix~C of \cite{Pa72}. 
For large $n$ the acceptance ratio $(V_n/B_n)$ of the selection process
approaches zero. In our following considerations the cases $n=2$ and 
$n=3$ are of importance. With $\Gamma (n+1) = n\,\Gamma(n)$, 
$\Gamma (1)=1$ and $\Gamma (1/2) =\pi^{1/2}$ the familiar results 
for the area of the unit circle and the volume of the unit sphere in 
three dimensions are obtained, so that the \index{acceptance rate!O(n) 
models} \index{O(n) models!acceptance rate}{acceptance ratios} become
\begin{equation} \label{SnBn_rat}
  {V_2\over B_2} = {\pi\over 2^2} = 0.785\dots ~~~{\rm and}~~~
  {V_3\over B_3} = {4\,\pi/3\over 2^3} = 0.523\dots\ .
\end{equation}
The Metropolis criterium~(\ref{Metropolis}) is implemented by the
Fortran statement
\begin{equation} \label{Metro_Fortran}
 {\tt if(rmafun().lt.exp(beta*dact))}\ {\tt then} 
\end{equation}
of the {\tt on\_met.f} subroutine. Note that {\tt on\_met.f} is 
written as a short routine at the price of performance. If a
similar Metropolis code has to be used for long simulations,
it is recommended to test possible improvements. For instance, 
the calculation of the exponential function is quite CPU time 
consuming and a code which by an additional {\tt if} statement 
accepts the proposals with ${\tt dact}\ge 0$ before
using (\ref{Metro_Fortran}) may be faster. Such conjectures need
to be tested in the actual computing environment, because they
depend on the CPU unit used as well as on the Fortran compiler. As
a general \index{programming}{programming} rule, it is a good idea 
to assume working with a dumb compiler, who will miss out on simple 
optimization opportunities. With this rule one will gain on computing 
speed with the dumb compilers without loosing on the others. In the 
code of {\tt on\_met.f} a precaution would be to modify the most inner 
{\tt do} loop by introducing an extra Fortran statement which
assigns a variable to the product ${\tt fn*stanew(i)}$. This would
insure that the multiplication is not done twice. Of course, a
reasonable compiler will figure that out anyway. Illustrations and 
comparisons of the actual performance of our code are given next.

\subsection{Metropolis code for the $XY$ model\label{XY_met}}

\index{Metropolis algorithm!XY model}\index{MC algorithms!XY model}The 
subroutines
\index{o2\_met.f}\index{xy\_met.f}\index{xy\_met0.f}\index{xy\_met1.f}
\begin{equation} \label{xy_met}
  {\tt o2\_ met.f},\ {\tt xy\_ met.f},\ {\tt xy\_ met0.f}~~{\rm and}~~ 
  {\tt xy\_ met1.f}~~{\rm of}~~{\tt FoLib\_ par}
\end{equation}
are all implementing the Metropolis algorithm for the $XY$ model. 
The corresponding main programs for initialization, equilibration 
and data (time series) production are, respectively,
\index{o2\_ts.f}\index{xy\_ts.f}\index{xy\_ts0.f}\index{xy\_ts1.f}
\begin{equation} \label{xy_ts}
  {\tt o2\_ ts.f},\ {\tt xy\_ ts.f},\ {\tt xy\_ ts0.f}~~{\rm and}~~ 
  {\tt xy\_ ts1.f}~~{\rm of}~~ {\tt ForProg/MC\_ On}\ .
\end{equation}
For initialization the programs {\tt xy\_ts.f} and {\tt xy\_ts0.f} use
\index{xy\_init.f}{\tt xy\_init.f}, relying on 
\index{xy\_order.f}{\tt xy\_order.f} or 
\index{xy\_ranf.f}{\tt xy\_ranf.f} and 
\index{xy\_act.f}{\tt xy\_act.f}, whereas {\tt xy\_ts1.f} needs its 
own initialization subroutine, \index{xy\_init1.f}{\tt xy\_init1.f}, 
relying on \index{xy\_ran1.f}{\tt xy\_ran1.f} and 
\index{xy\_act1.f}{\tt xy\_act1.f}.
The program {\tt o2\_ts.f} is a specialized version of {\tt on\_ts.f}
(\ref{on_ts}) and uses the same routines as {\tt on\_ts.f} for
initialization. The analysis program 
\begin{equation} \label{ana_ts_on}
  \index{ana\_ts\_on.f}{\tt ana\_ ts\_ on.f}~~~{\rm of}
  ~~~{\tt ForProg/MC\_ On}
\end{equation}
can be used on the output files of all the programs~(\ref{xy_ts}). 
It allows also to analyze the time series for autocorrelations as 
needed in chapter~\ref{Markov_Error} (compare 
table~\ref{tab_xy_iat}).

The main deviation of the Metropolis routine 
\index{xy\_met0.f}{\tt xy\_met0.f} from 
{\tt on\_met.f} is that it specializes to $n=2$ and uses a direct 
calculation of the newly proposed coordinates, 
\begin{equation} \label{xy_phi}
x=\cos (\phi),~~y=\sin (\phi)~~~{\rm and}~~~\phi=2\,\pi\,x^r\,,
\end{equation}
where $x^r\in[0,1)$ is a uniformly distributed random number. As 
the calculation of the sine and cosine functions is notoriously
slow on some computers, it requires testing to find out, whether
the random procedure with the 0.785\index{acceptance rate!XY model} 
\index{XY model!acceptance rate}{acceptance rate} (\ref{SnBn_rat}) 
or the direct calculation (\ref{xy_phi}) is faster. 
If the random procedure turns out to be faster, a code which is 
specialized to $n=2$ should still be used but the lines for the 
direct calculation should be replaced by those between {\tt label 1} 
and {\tt go to 1} in table~\ref{tab2_on_met}.

The \index{acceptance rate!XY model}{acceptance rate}
for the Metropolis update of {\tt xy\_met0.f} goes
to zero for large $\beta$. The code of {\tt xy\_met.f} adds an 
additional twist by adjusting the new spin proposals so that 
the acceptance rate becomes never smaller than a preset value
{\tt acpt\_rate}. Instead of defining the angle $\phi$ of equation 
(\ref{xy_phi}) with respect to the $x$-axis and proposing it 
uniformly in the range $[0,2\pi)$, $\phi$ is now defined with 
respect to the direction of the sum of the neighboring spins and 
proposed in a more sophisticated way: Whenever the acceptance rate 
falls below the value {\tt acpt\_rate}, which is set to 0.4 by
default (justified by figure~\ref{fig_acpt_tauint} of the next
chapter), the range $[-\phi_{\max},\phi_{\max}]$ of the $\phi$ 
proposals is reduced to an interval smaller than $[-\pi,\pi)$ 
until {\tt acpt\_rate} is reached. The other way round, if the 
acceptance rate increases above {\tt acpt\_rate} the range is 
extended until either the {\tt acpt\_rate} acceptance rate value 
or the maximum $\phi$ range $[-\pi,\pi)$ is reached. At small 
$\beta$ acceptance rates with the spin proposals~(\ref{xy_phi})
are, of course, already larger than 40\%.

To test our code, we compare two energy values with data reported
in a paper by \cite{ToCh79}, which work was one of the first MC 
simulations of the $XY$ model and investigated the Kosterlitz-Thouless 
scenario in some details. In their tables~I and~II Tobochnik and 
Chester give energy values for a $60\times 60$ lattice at different 
temperatures, each data point relying on 
1$\,$000 sweeps for reaching equilibrium followed by 2$\,$000 
sweeps with measurements. For their data in table~I they cool the 
initially random system and for table~II they warm the initially 
ordered system. For comparison, we pick the high temperature 
value $1/T=\beta =0.5$, where the system moves quite freely, and 
$1/T=\beta = 1$, a value in the critical region, slightly below 
the peak in the \index{specific heat}{specific heat}. Assuming equal 
weights we combine the data given in the two tables. This gives
(assignment~\ref{assign_xy_ToChe_s}).
\begin{equation} \label{xy_ToChe_s}
  e_s = - 0.54775\,(43)~~{\rm at}~~\beta = 0.5~~~{\rm and}~~~
  e_s = - 1.3213\, (15)~~{\rm at}~~\beta = 1\ .
\end{equation}
Nowadays it is easy to produce a large statistics. We perform 
10$\,$000 sweeps for reaching equilibrium followed by 
$32 \times 2\,000$ sweeps with measurements. Running the 
general $O(n)$ code {\tt on\_ts.f}~(\ref{on_ts}) as well
as {\tt xy\_ts0.f}~(\ref{xy_ts}) on a laptop with Cygwin we
obtained at $\beta=0.5$ (see  assignment~\ref{assign_xy_hTe_s}). 
\begin{equation} \label{xy_hTe_s}
  e_s = - 0.54756\,(13)~~{\rm with}~~{\tt on\_ ts.f}~~~{\rm and}~~~
  e_s = - 0.54762\,(14)~~{\rm with}~~{\tt xy\_ ts0.f}\ .
\end{equation}
We use the Gaussian difference
test (\ref{gau_dif}) to compare these estimates with one another as
well as with $\beta=0.5$ the estimate (\ref{xy_ToChe_s}) and find 
the satisfactory $Q$ values 0.75, 0.67 and 0.77. As the results
of equation (\ref{xy_hTe_s}) may depend on the computer platform
on which the program is run, the reader may get different numbers,
which ought to be statistically consistent with those reported here.

\subsubsection{Timing, discretization and rounding errors}

\index{CPU time!measurements}{Timing} the {\tt on\_ts.f} and
{\tt xy\_ts0.f} simulations against one another gives machine 
dependent results. We report here some results from runs performed 
in the year~2002. On a \index{Linux}Linux PC with the g77 compiler 
the {\tt on\_ts.f} 
simulation took longer than the {\tt xy\_ts0.f} simulation: $514\,$s 
versus $409\,$s. On a \index{Sun workstations}{Sun workstation} the 
results were closer together and in the reverse order: $313\,$s with 
{\tt on\_ts.f} and $336\,$s with {\tt xy\_ts0.f}. The routine 
{\tt o2\_met.f} illustrates that we can speed up the methods of the 
general {\tt on\_met.f} Metropolis code by specializing to $n=2$. 
In particular, {\tt o2\_met.f} avoids the calculation of 
{\tt exp(beta*dact)} and one random number for the acceptance 
case ${\tt dact}\ge 0$. The corresponding Fortran statement
\begin{equation} \label{if_dact}
{\tt if(dact.lt.zero)\ then }
\end{equation}
of {\tt o2\_met.f} is also used in {\tt xy\_met0.f}. The driving 
program for the simulations with {\tt o2\_met.f} is {\tt o2\_ts.f}. 
With {\tt o2\_ts.f} the simulations for equation~(\ref{xy_hTe_s}) 
took $459\,$s on the Linux PC and $260\,$s on the Sun workstation.

Out of this competition is the program \index{xy\_ts1.f}{\tt xy\_ts1.f} 
(using \index{xy\_met1.f}{\tt xy\_met1.f}) for which the 
same run took $960\,$s on the Linux PC and $1\,199\,$s on the Sun 
workstation. Instead of using the $x,\,y$
coordinates, this program stores only the angles $\phi$. Storage
space requirements are reduced by a factor of two, but the 
calculation of the energy differences from the angles is far
slower than from the coordinates. To use such a program makes only
sense when RAM is a serious problem.

For {\tt xy\_ts0.f} the result on the Sun workstation was 
$e_s=-0.54739\,(13)$, agreeing in the statistical sense of $Q=0.23$ 
for the Gauss difference test with the value~(\ref{xy_hTe_s}). At 
the first look the difference is surprising. Due to the portability 
of the 
Marsaglia's random number generator, the proposed $XY$ coordinates 
differ only by the eventually distinct rounding errors of the sine 
and cosine 
functions in equation (\ref{xy_phi}). Thus the energy differences 
encountered by the Metropolis acceptance steps are expected to be 
governed by the machine rounding error, which is for double 
precision about $2.2\times 10^{-16}$, compare paragraph 20.1 of 
{\it Numerical Recipes}~\cite{PrFl86}. In our MC calculations we 
use $2.304\times 10^8$ random numbers per run for accept/reject 
decisions and the likelihood that the two calculations diverge on
different machines is less than $10^{-7}$. So, what is going on? 

What happens is uncovered in assignment~\ref{assign_xy_ts0}. 
Due to the fact that there are only $2^{24}=16,777,216$ distinct 
Marsaglia random numbers (assignment~\ref{assign_ran_range} of 
section~\ref{RanNum}) there is a discretization of the angle $\phi$, 
and, hence, the $XY$ coordinates (\ref{xy_phi}). Therefore,
it will occasionally happen that precisely the present (old) $XY$
spin is proposed at some site. The energy difference is then zero 
up to a $\pm \epsilon_{\rm machine}$ rounding error, which turns 
out to be exactly zero on the Sun, but around $\pm 10^{-18}$ on 
the Linux PC. This leads to different decisions in the if statement 
(\ref{if_dact}) and an additional random number is or is not used. 
Once this happens, the calculations on the two  machines diverge. If 
one is unhappy with the accuracy of this discretization of the angle 
$\phi$ in equation (\ref{xy_phi}), one can greatly improve on it by 
using the {\tt ranmar2.f} (\ref{ranmar2}) generator based on two 
Marsaglia random numbers.

\subsubsection{Acceptance rate \label{sec3_acpt}}

Repeating our simulation (\ref{xy_hTe_s}) at $\beta=1$ gives
(assignment~\ref{assign_xy_T1e_s}) 
\begin{equation} \label{xy_T1e_s}
e_s = -1.3197\, (11) ~~{\rm with}~~{\tt on\_ ts.f} ~~~{\rm and}~~~
e_s = -1.3207\, (11) ~~{\rm with}~~{\tt xy\_ ts0.f}\ .
\end{equation}
The Gaussian difference tests yield $Q= 0.52$, $0.39$ and 0.75. 

Whereas the Metropolis\index{acceptance rate!XY model}
\index{XY model!acceptance rate}{acceptance rate} for the high 
temperature simulations (\ref{xy_hTe_s}) at $\beta = 0.5$ was about 
$0.64$, it is now down to $0.32$, and it will go further down with 
decreasing temperature. A very low acceptance rate is bad. The routine 
{\tt xx\_met.f} allows to increase the acceptance rate by proposing an 
angle $\phi$ in the range
\begin{equation} \label{phi_add}
  0\ \le\ |\phi|\ < \triangle\Phi~~{\rm with}~~0<\triangle\Phi<\pi 
\end{equation}
with respect to the direction of the sum of the neighboring spins.
The proposed new spin at site $i$ is then given by the equations
\begin{equation} \label{xy_new}
  x_i^{\rm new}=\cos(\phi )\,x_i^{\rm old}-\sin(\phi)\,y_i^{\rm old}
  ~~~{\rm and}~~~
  y_i^{\rm new}=\sin(\phi )\,x_i^{\rm old}+\cos(\phi)\,y_i^{\rm old}\,.
\end{equation}
By decreasing the value of $\triangle \Phi$, one increases the 
acceptance rate, because the proposed energy differences become 
smaller. A large 
acceptance rate is not good too: In the limit where the acceptance 
rate goes to one, we find $\triangle\phi =0$, and the system does not 
move anymore. On the other side, the system does not move either when 
the acceptance rate is zero. What is the optimal acceptance rate? In 
essence, this question can only be answered by a careful study of 
autocorrelation times (see the next chapter). In practice one may 
not care to be somewhat off the optimum. The rule of thumb is: An 
acceptance rate in the range 30\% to 50\% works well (compare
figure~\ref{fig_acpt_tauint}).

We use the notation\index{natural acceptance rate|textbf} 
\index{acceptance rate!natural|textbf}{\bf natural acceptance rate} 
for the value obtained with $\triangle\Phi =\pi$ in 
equation~(\ref{phi_add}). The Metropolis routine {\tt xy\_met.f} 
allows for adjustments of its acceptance rate to the value of its 
parameter {\tt acpt\_rate}, unless {\tt acpt\_rate} is smaller than 
the natural acceptance rate, which is its lower bound. The preset 
default value of {\tt actp\_rate} is 0.4 and can be changed 
in the routine. Whenever the argument {\tt lacpt} of {\tt xy\_met.f} 
is set to {\tt .true.}, the routine adjusts $\triangle \Phi$ to 
attain the desired acceptance rate. For this task it relies on its 
so far assembled statistics to estimate $\triangle \Phi$. The driving 
program {\tt xy\_ts.f} introduces a schedule for {\tt lacpt=.true.}$\!$ 
updates of $\triangle\Phi =\pi$ during the equilibration part of the 
run. 
% Alternatively one may adjust the acceptance rate during the entire
% run as is discussed in appendix~\ref{Exercises}, exercise~1 to the 
% present section 

For the simulation of assignment~\ref{assign_xy_T1e_s} the natural 
acceptance rate is about 32\%. In assignment~\ref{assign_xy_acpt} no 
improvement of the error bar is found by producing the same statistics 
with {\tt xy\_ts.f} at an acceptance rate of 40\%: The result 
\begin{equation} \label{xy_T1e_s50}
  e_s=-1.3202\, (10) 
\end{equation}
is well compatible with the estimates~(\ref{xy_T1e_s}). With much higher 
statistics the integrated autocorrelation times are estimated in 
assignment~\ref{assign_xy_iat} of chapter~\ref{sec_auto_algo}. One 
should however also note that the routine with the larger acceptance 
rate, normally {\tt xy\_ts.f}, will often run faster, because there
are typically fewer evaluations of the time consuming exponential 
function necessary.

Really high $\beta$ values are needed to find a significant improvement 
of error bars due to the adjustment of the acceptance rate to 40\%. In 
assignment~\ref{assign_xy_ts} MC simulations are performed on $10\times 
10$ lattice, for which a high statistics is assembled at $\beta=10$ and 
$\beta =20$, divided into $\tt nrpt =100$ blocks. The energy results at 
$\beta =10$ are (as before on a laptop under Cygwin)
\begin{eqnarray} \nonumber 
  e_s&=&-1.9498\,61\,(39)~~{\rm with}~~{\tt xy\_ ts0.f}~~~{\rm and}\\
  e_s&=&-1.9498\,69\,(17)~~{\rm with}~~{\tt xy\_ ts.f}\ .\label{xy_b10}
\end{eqnarray}
Comparing the error bars with the F-test (\ref{F_test}) gives
$Q<10^{-14}$. The {\tt xy\_ts.f} simulation performs significantly 
better than the one with {\tt xy\_ts0.f}. For the latter the acceptance 
rate was 0.081. At $\beta=20$ the results are
\begin{eqnarray} \nonumber
  e_s&=&-1.7840\,87\,(23)~~{\rm with}~~{\tt xy\_ ts0.f}~~~{\rm and}\\
  e_s&=&-1.9750\,96\,(08)~~{\rm with}~~{\tt xy\_ ts.f}\ . \label{xy_b20}
\end{eqnarray}
%
\begin{figure}[tb]
\vspace{5pc}
 \centerline{\hbox{ \psfig{figure=figures/domain.eps,width=10cm} }}
 \caption{\index{domain wall}Domain wall configuration encountered 
in one of the simulations when quenching the $10\times 10$ $XY$ 
system to $\beta =20$, see assignment~\ref{assign_xy_ts}. 
\label{fig_xyd} } 
\end{figure}
%
The acceptance rate of the {\tt xy\_ts0.f} simulation at $\beta=20$
is 0.060 and the error bar of the {\tt xy\_ts.f} simulation is 
again the significantly smaller one. However, most remarkable is 
that the mean values do not agree at all. The value obtained 
with {\tt xy\_ts0.f} is far too high (several thousand standard
deviations). It turns out that this has nothing to do with the 
performance of the {\tt xy\_met0.f} Metropolis algorithm. 
Figure~\ref{fig_xyd} is taken after the 100$\,$000 equilibration 
sweeps of this simulation and plots the $XY$ spins of the $2d$
lattice with sites labeled by coordinates $x$ and $y$. It shows 
that the system got stuck in a \index{domain wall}{domain wall} 
configuration from which there is practically no escape 
by small fluctuations. The spins at $x=3$ point down, whereas the 
spins at at $x=8$ point up, creating two domains in between. 
This can happen due to rapid \index{quenching}{quenching} 
of the system from a disordered start to almost zero temperature. 
By lowering the temperature more slowly, {\it i.e.} 
\index{annealing}{annealing}, one can avoid this. Or one can
simply use an ordered start configuration for which the result 
of the {\tt xy\_ts0.f} simulation becomes $e_s=-1.9750\,92\,(22)$. 

\subsection{Heatbath code for the $O(3)$ model} \label{O3_hb}

The subroutine\index{O(3) $\sigma$-model!heatbath algorithm}
\index{heatbath algorithm!O(3) $\sigma$-model}\index{o3\_mchb.f}
\index{heatbath algorithm!Heisenberg ferromagnet}
\index{Heisenberg ferromagnet!heatbath algorithm}
\begin{equation} \label{o3_mchb}
  {\tt o3\_ mchb.f}~~~{\rm of}~~~{\tt ForLib\_par}
\end{equation}
is an implementation of the $O(3)$ heatbath algorithm of 
section~\ref{O3_heatbath}. It is run by the driving program
\begin{equation} \label{o3_tshb}\index{o3\_tshb.f}
  {\tt o3\_ tshb.f}~~~{\rm of}~~~{\tt ForProg/MC\_ On}\ .
\end{equation}
Following the usual structure it does first initialization, then
equilibration and finally \index{data production}{data production}.
For the initialization the routines \index{o3\_init.f}{\tt o3\_init.f}, 
\index{o3\_order.f}{\tt o3\_order.f} or \index{o3\_ran.f}{\tt o3\_ran.f} 
and \index{o3\_act.f}{\tt o3\_act.f} are invoked. 
The last two routines can be replaced by the general $O(n)$ 
routines {\tt on\_ran.f} and {\tt on\_act.f}, if desired.

To test the $O(3)$ model code, we simulate a $100\times 100$ lattice
at $\beta=1.1$ to compare with an energy values of table~I in
\cite{BeLu81}, where the 
\index{topological charge}{topological charge} 
of the $O(3)$ $\sigma$-model was 
investigated\footnote{In QCD the $\eta$--$\eta'$ mass difference is 
related to fluctuations of the \index{instanton}{instanton}
topological charge~\cite{tH76}. The $O(3)$ $\sigma$-model in $d=2$ 
serves as a toy model for related investigations. To get physical 
results, the approach of \cite{BeLu81} needs to be refined. See 
\cite{Ne99} for a review of lattice gauge theory work on this 
subject.}. 
We perform 10$\,$000 sweeps for reaching equilibrium and 
$20\times 2\,500$ sweeps with measurement. Our results on the 
1.7~GHz laptop are (see assignment~\ref{assign_o3_e_s})
\begin{eqnarray} \nonumber
  e_s&=&-0.852\,06\,(12)~~{\rm with}~~{\tt on\_ ts.f}~~~
  {\rm and}\\ \label{o3_e_s}
  e_s&=&-0.8519\,87\,(80)~~{\rm with}~~{\tt o3\_ tshb.f}\ .
\end{eqnarray}
They are consistent with one another. The Gaussian difference 
test gives $Q=0.61$. The same is true for the error bars. The
F-test (\ref{F_test}) gives $Q=0.09$, a value too large to call
the smaller error bar of the heatbath simulation statistically
relevant. This is clarified by increasing the production
part of the statistics from twenty to one hundred bins
(assignment~\ref{assign_o3_large}). The error bar of the heatbath
simulation remain the smaller one, 33 to 66 is the ratio, and
the $Q$ value of the F-test decreases to $3.3\times 10^{-10}$. The 
performance of the heatbath simulation is significantly better 
than that of our Metropolis simulation and this improvement wins 
over the increased CPU time consumption of the heatbath code: For 
the simulations of equation (\ref{o3_e_s}) the CPU time used
on the laptop was $325\,$s for {\tt on\_ts.f} and $566\,$s 
for {\tt o3\_ts.f}. 

The $\beta=1.1$ energy value of table~I in \cite{BeLu81} translates 
to $e_s=0.85250\,(25)$\footnote{The relation $e_s=(e_{BL}-4)/2$ holds, 
where $e_{BL}$ is the value given by \cite{BeLu81}. Compare their 
equation (29). The error bar used here assumes maximum rounding 
towards their last digit given.}.
Gaussian difference tests show that this value is consistent with 
those of equation~(\ref{o3_e_s}): $Q=0.11$ and $Q=0.05$.

\begin{figure}[tb] \vspace{5pc}
\centerline{\hbox{ \psfig{figure=figures/qdf_o3.eps,width=10cm} }}
\caption{Illustration of self-averaging: 
\index{distribution!function, peaked!O(3) $\sigma$-model}
\index{peaked distribution function!O(3) $\sigma$-model}Peaked 
distribution functions from MC simulations of the $O(3)$ 
$\sigma$-model at $\beta=1.1$, see assignments~\ref{assign_o3_e_s} 
and~\ref{assign_qdf_o3}.  \label{fig_qdf_o3} }
\end{figure}

In figure~\ref{fig_qdf_o3} we compare the peaked distribution 
function (\ref{qdf}) of the mean action per link from the simulations 
of assignment~\ref{assign_o3_e_s} with the corresponding peaked 
distribution functions from simulations on smaller lattices
(assignment~\ref{assign_qdf_o3}). The property of {\bf self-averaging} 
\index{self-averaging}is observed: The larger the lattice, the smaller 
the range of a given confidence level. When self-averaging is at
work, the peaked distribution function will with increasing lattice 
size approach the function which is $0.5$ at the median and zero 
otherwise. The other way round, the peaked distribution function is 
very well suited to investigate observables for which self-averaging
does not work, as for instance encountered in spin glass 
simulations~\cite{BeBiJa00}.

\subsubsection{Rounding errors}

In the year 2002 we compared runs from a Sun workstation with those on 
a Linux PC and we found that the values obtained with {\tt o3\_ts.f} 
differed, but remained consistent in the statistical sense. For the 
$XY$ model we encountered a sudden divergence of results due to a 
distinct use of a random number as consequence of a rounding error. 
This is not the case here. In the heatbath algorithm all random numbers 
are used in an identical way on both platforms and we found a slow 
divergence of results, which can be traced to the rounding errors in 
the line
$$ {\tt snew3=log(Sexp+rmafun()*(exp(-Sbeta)-Sexp))/Sbeta} $$
of the heatbath routine {\tt o3\_mchb.f}. This line is the Fortran 
implementation of equation~(\ref{O3_heat_bath-xr}).

\subsection{Assignments for section~\ref{MC_continuous}}
\label{assign_MC_continous}

In the following random starts and Marsaglia random number with their
default seed are used unless stated otherwise.

\begin{enumerate}

{\item The energy per spin values in \cite{ToCh79}
at $\beta=0.5$ are $e_s=-0.5470\,(7)$ in table~I and 
$e_s=-0.5485\,(5)$ in table~II.  At $\beta =1$ the values are 
$e_s=-1.3188\,(23)$ in table~I and $e_s=-1.3237\,(19)$ in table~II.
Use the subroutine {\tt steb2.f} (\ref{STEB2}) to reproduce the
estimates of our equation~(\ref{xy_ToChe_s}).
\label{assign_xy_ToChe_s} }

{\item Use the programs {\tt on\_ts.f} (\ref{on_ts}) and {\tt 
xy\_ts0.f} (\ref{xy_ts}) to simulate the $XY$ model at $\beta=0.5$ on 
an $60\times 60$ lattice with 10$\,$000 sweeps for reaching equilibrium 
followed by $32\times 2\,000$ sweeps with measurements. Subsequently, 
run the analysis program {\tt ana\_ts\_on.f} to obtain energy values.
Perform Gaussian difference tests to compare the results with one 
another, with equation~(\ref{xy_hTe_s}) as well as with the $\beta=0.5$ 
estimate (\ref{xy_ToChe_s}). Time the two
simulations against one another, as well as against {\tt o2\_ts.f} 
and {\tt xy\_ts1.f}.  Explain why the energy values obtained with 
{\tt o2\_ts.f} are not to all digits identical with those from 
{\tt on\_ts.f}.  \label{assign_xy_hTe_s} }

{\item Investigate the equilibrium sweep ${\tt iequi} = 3707$
of the {\tt xy\_ts0.f} run from the previous assignment and 
print out the {\tt dact} values of the subroutine {\tt xy\_met0.f}
for $2390\le {\tt is}\le 2395$. Explain the value which you get
for ${\tt is} = 2393$. \label{assign_xy_ts0} }

{\item Repeat assignment~\ref{assign_xy_hTe_s} for $\beta = 1$ 
and compare he energy values with those reported in 
equation~(\ref{xy_T1e_s}).  \label{assign_xy_T1e_s} }

{\item Repeat the previous assignment using {\tt xy\_ts.f}. 
\label{assign_xy_acpt} }

{\item Use {\tt xy\_ts0.f} and {\tt xy\_ts.f} for imulations with 
${\tt nequi}=100\,000$, ${\tt nrpt}=100$ and ${\tt nmeas}=10\,000$ 
on a $10\times 10$ lattice at $\beta=10$ and $\beta=20$. Perform the 
F-test (\ref{F_test}) to compare the error bars of the two simulations. 
Plot the spins after equilibration for the $\beta=20$ run with 
{\tt xy\_ts0.f}. Repeat the $\beta=20$ {\tt xy\_ts0.f} simulation
using an ordered start. \label{assign_xy_ts} }

{\item Use the programs {\tt on\_ts.f} (\ref{on_ts}) and 
{\tt o3\_tshb.f} (\ref{o3_tshb}) to simulate the $O(3)$ model 
at $\beta=1.1$ on a $100\times 100$ lattice with 10$\,$000 sweeps 
for reaching equilibrium followed by $20\times 2,500$ sweeps 
with measurements. Subsequently, run the analysis program {\tt 
ana\_ts\_on.f} to obtain the energy values. Perform Gaussian 
difference tests to compare the results with one another, with
those of equation (\ref{o3_e_s}) as well as with the estimate from 
\cite{BeLu81}. Perform the F-test (\ref{F_test}) to compare the two 
error bars from your simulations.  \label{assign_o3_e_s} }

{\item Increase the production statistics of the previous assignment
to $100\times 2,500$ sweeps and repeat the other parts of that
assignment.  \label{assign_o3_large} }

{\item Repeat the heatbath simulations of assignment~\ref{assign_o3_e_s}
({\it i.e.}, the same statistics and $\beta$ value) on lattices of
size $L=20$, $34$ and $60$. Reproduce figure~\ref{fig_qdf_o3}.
\label{assign_qdf_o3} }

\end{enumerate}

