
% BB Nov 19  2007. Adapted to material already covered.
% BB Nov 21  2006.
% BB Mar 01  2004.

\documentclass[20pt,landscape,footrule]{foils}
\usepackage[dvips]{graphicx}
\usepackage[pdftex]{geometry}
\geometry{headsep=3ex,hscale=0.9}
\usepackage[pdftex]{color} % for the colored block
\date{\color{blue} Singapore MCMC Meeting, March 2004}
\usepackage[ignore]{pause}
%\usepackage{pause}
\usepackage{background}
%\usepackage{pp4slide}
%\MyLogo{AV HEP,CSIT}
\LogoOff
\begin{document}
\def\r{\color{red}}
\def\b{\color{blue}}
\def\s{\color{black}}

\foilhead{\r Lecturenotes MCMC IV -- Contents }
\b \vskip -10.0pt

\begin{enumerate}
\item Multicanonical Ensemble
\item How to get the Weights?
\item Example Runs ($2d$ Ising and Potts models)
\item Re-Weighting to the Canonical Ensemble
\item Energy and Specific Heat Calculation
\item Free Energy and Entropy Calculation
\item Summary
\end{enumerate}

% \noindent
 
\foilhead{\r Multicanonical Ensemble} 
\vskip -30pt

One of the questions which ought to be addressed before performing
a large scale computer simulation is {\r ``What are suitable weight
factors for the problem at hand?''} So far we used the Boltzmann 
weights as this appears natural for simulating the
Gibbs ensemble. However, a broader view of the issue
is appropriate.  

Conventional canonical simulations can by re-weighting techniques only 
be extrapolated to a vicinity of this temperature. For multicanonical 
simulations this is different: A single simulation allows to obtain 
equilibrium properties of the Gibbs ensemble over a range of 
temperatures. Of particular interest are two situations for which 
canonical simulations do not provide the appropriate implementation 
of importance sampling: 
\begin{enumerate}
\item The physically important configurations are rare in the 
      canonical ensemble.
\item A rugged free energy landscape makes the physically important 
      configurations difficult to reach. 
\end{enumerate}

\newpage
Multicanonical simulations approach these problem by trying to sample, 
in an appropriate energy range, with an {\bf approximation}
\begin{equation}\s \label{wmuca}
  \widehat{w}_{mu}(k) = w_{mu}(E^{(k)}) 
  = e^{-b(E^{(k)})\, E^{(k)} + a (E^{(k)})}
\end{equation}
to the weights 
\begin{equation}\s \label{wspectral}
 \index{multicanonical!weights}\index{weights!multicanonical}
 \widehat{w}_{1/n} (k) = w_{1/n} (E^{(k)}) = {1\over n(E^{(k)})}
\end{equation}
where $\s n(E)$ is the {spectral density}. The function $\s b(E)$ 
defines the inverse \textbf{microcanonical temperature} and $\s a(E)$ 
the {\bf dimensionless, microcanonical free energy}. The function 
$\s b (E)$ has a relatively smooth dependence on its arguments, which 
makes it a useful quantity when dealing with the weight factors.  

Instead of the canonical energy distribution $\s P(E)$,
one samples a new multicanonical distribution, which does not longer 
depend strongly  on the energy, 
\begin{equation}\s \label{Pmuca}
 P_{mu} (E) = c_{mu}\, n(E)\, w_{mu} (E) \approx c_{mu}\ .
\end{equation}
The desired canonical probability density is obtained by
re-weighting
\begin{equation}\s \label{reweight}
 P (E) = {c_{\beta}\over c_{mu}}\, {P_{mu}(E) \over w_{mu} (E)}\,
         e^{-\beta E} .
\end{equation}
This relation is rigorous, as the weights $\s w_{mu}(E)$ used 
in the simulation are exactly known. Statistical errors are those 
of $\s P_{mu}(E)$ times an exactly known factor. 

The multicanonical method requires {\r two steps:}

\begin{enumerate}

\item Obtain a {\r working estimate} $\s \widehat{w}_{mu}(k)$ of 
      the weights $\s \widehat{w}_{1/n}(k)$. Working estimate means
      that the approximation has to be good enough to ensure movement 
      in the desired energy range, while deviations from a flat 
      $\s P_{mu} (E)$ are tolerable.

\item Perform a Markov chain MC simulation with the {\r fixed} weights 
      $\s \widehat{w}_{mu}(k)$. The thus generated configurations
      constitute the {multicanonical ensemble}. Canonical expectation 
      values are found by re-weighting to the Gibbs ensemble and 
      jackknife methods allow reliable error estimates.

\end{enumerate}

\vfil
\newpage

\foilhead{\r How to get the Weights?}\vskip -30pt

To get the weights is at the heart of the method (and a stumbling block 
for beginners). Some approaches used are listed in the following:

\begin{enumerate}
\item Overlapping constrained (microcanonical) MC simulations.
      Potential problem: ergodicity.
\item {\r Finite Size Scaling (FSS) Estimates.} Best when it works!
      Problem: There may 
      be no FSS theory or the prediction may be so inaccurate that the 
      initial simulation will not cover the target region.
\item General Purpose Recursions. Problem: Tend to deteriorate with
      increasing lattice size (large lattices). 
\end{enumerate}

\vfil
\newpage

\noindent
{\large\r Multicanonical Recursion} {\b\small (A variant of Berg, 
J. Stat. Phys. 82 (1996) 323.)}

\noindent
{\b Multicanonical parameterization of the weights:}
\begin{eqnarray}\s \nonumber 
  w(a) = e^{-S(E_a)} &\s =&\s e^{-b(E_a)\, E_a + a(E_a)}\ , \\ \s 
  \nonumber {\b\rm where}~~b(E) &\s =&\s 
  \left[S(E+\epsilon)-S(E)\right]/\epsilon
  ~~(\epsilon\ {\rm smallest\ stepsize}) \\ \s {\rm and}~~\nonumber 
  a(E-\epsilon) &\s =&\s a(E)+\left[ b(E-\epsilon)-b(E) \right]\, E\ . 
\end{eqnarray}
{\r Recursion:}
\begin{eqnarray}\s \nonumber 
   \r b^{n+1}(E) &\r =&\r b^n (E)+\hat{g}^n_0(E)\,
 [ \ln H^n(E+\epsilon)-\ln H^n(E)] / \epsilon \\ \s \nonumber 
 \hat{g}^n_0 (E) &\s =&\s g^n_0(E)\, /\, [g^n(E) + \hat{g}^n_0 (E)]\, ,\\
 \s\nonumber g^n_0 (E) &\s =&\s H^n (E+\epsilon)\, H^n (E)\, /\,
   [H^n (E+\epsilon) + H^n (E)]\, ,\\ \nonumber
  \s g^{n+1} (E) &\s =&\s g^n(E) + g^n_0(E),\ g^0(E)=0\, .
\end{eqnarray}
{\r ANIMATION:} {\b 2d 10-state $80\times 80$ Potts model
(Iterations of 100 sweeps each).}

\foilhead{\r Example Runs ($2d$ Ising and Potts models)} \vskip -30pt
\noindent Ising model on a $\s 20\times 20$ lattice: The multicanonical 
recursion is run in the range
\begin{equation}\s \label{2dI_namin_max}
  {\tt namin} = 400 \le {\tt iact} \le 800 = {\tt namax}\ . 
\end{equation}
The recursion is terminated after a number of so called tunneling 
events. A {\r tunneling event} is defined as an updating process 
which finds its way from
$$\s {\tt iact}={\tt namin}~~~{\rm to}~~~{\tt iact}={\tt namax}
     ~~~{\rm and}~~{\rm back}\ . $$
This notation comes from applications to first order phase transitions. 
An alternative notation for tunneling event is {\r random walk cycle}.  
For most applications 10 tunneling events lead to acceptable weights.

For an example run of the Ising model we find the requested {\s 10
tunneling} events after {\s 787 recursions} and {\s 64,138 sweeps}
(assignment {\tt a0501\_01}, {\tt a0501\_02} for the $2d$ 10-state 
Potts model).

% \newpage
\foilhead{\r Performance}\vskip -30pt

If the multicanonical weighting would remove all relevant {free energy 
barriers}, the updating process would become a free {\r random walk}.  
Therefore, the theoretically optimal performance for the second part 
of the multicanonical simulation is
$$\s \tau_{\rm tun}\ \sim\ V^2\ .$$
Recent work about first order transitions by Neuhaus and Hager shows 
that that the multicanonical procedure removes only the leading free 
energy barrier, while at least one subleading barrier causes still a 
residual supercritical slowing done. Up to certain medium sized lattices 
the behavior {\s $V^{2+\epsilon}$} gives a rather good effective 
description. For large lattices supercritical slowing down of the form
$$\s \tau_{\rm tun}\ \sim\ \exp\left(+{\rm const}\,L^{d-1}\right)$$
dominates again. The slowing down of the weight recursion with the 
volume size is expected to be even (slightly) worse than that of the 
second part of the simulation. 

\foilhead{\r Re-Weighting to the Canonical Ensemble}\vskip -30pt

Let us assume that we have performed a multicanonical simulation
which covers the action histograms needed for a temperature range
$$\s \beta_{\min}\ \le\ \beta = {1\over T}\ \le\ \beta_{\max}\ .$$
Given the {multicanonical time series}, where $i=1,\dots,n$ labels 
the generated configurations. The formula 
$$\s \overline{\mathcal{O}} = { \sum_{i=1}^n \mathcal{O}^{(i)}\, \exp 
  \left[ -\beta\,E^{(i)} + b(E^{(i)})\,E^{(i)}-a(E^{(i)})\right] \over 
  \sum_{i=1}^n \exp \left[ -\beta\,E^{(i)} 
  + b(E^{(i)})\,E^{(i)}-a(E^{(i)}) \right] }\ .$$
replaces the multicanonical weighting of the simulation by the Boltzmann 
factor. The denominator differs from the partition function $\s Z$ by a 
constant factor which drops out.

\newpage For discrete systems it is sufficient to keep histograms
when only functions of the energy are calculated. For an operator 
$\s \mathcal{O}^{(i)} = f(E^{(i)})$ and reweighting simplifies to
$$\s \overline{f} = { \sum_E f(E)\, h_{mu}(E)\, \exp 
  \left[-\beta\,E+b(E)\,E-a(E) \right]\over \sum_E h_{mu}(E)\, 
  \exp \left[-\beta\,E+b(E)\,E-a(E)\right]} $$
where $\s h_{mu}(E)$ is the histogram sampled during the multicanonical
production run and the sums are over all energy values for which
$\s h_{mu}(E)$ has entries.

\newpage
The computer implementation of these equations requires care. The 
differences between the largest and the smallest numbers encountered 
in the exponents can be really large. We can {\r avoid large numbers 
by dealing only with logarithms} of sums and partial sums: 

For $\s C=A+B$ with $\s A>0$ and $\s B>0$ we can calculate 
$$\s \ln C = \ln (A+B)~~~{\rm from}~~~\ln A~~{\rm and}~~ \ln B$$ 
without ever storing either $\s A$ or $\s B$ or $\s C$:
\begin{eqnarray} \s \nonumber
  \ln C &\s =&\s \ln \left[ \max (A,B)\,\left( 1+{\min (A,B)\over 
   \max (A,B)} \right) \right] \\ \nonumber
  &\s =&\s \max \left( \ln A,\ln B \right) + \ln \{ 1 +
  \exp \left[ \min (\ln A,\ln B)-\max(\ln A,\ln B)\right] \} 
  \\ \nonumber &\s =&\s \max \left( \ln A,\ln B \right) + 
  \ln \{ 1 + \exp \left[ - |\ln A - \ln B| \right] \}\ .
\end{eqnarray}

\foilhead{\r Energy and Specific Heat Calculation}\vskip -30pt

We are now ready to produce multicanonical data for the {\bf energy per 
spin} (with jackknife errors) of the $2d$ Ising model on a {\s $20\times 
20$ lattice} and to compare them with the exact results of Ferdinand 
and Fisher (assignment {\tt a0501\_03}):

\begin{figure} \begin{center}
\includegraphics[height=0.64\textheight]{figures/2dI_es.pdf}
\end{center} \end{figure}

\newpage\noindent
The same numerical data allow to calculate the {\bf specific heat} 
defined by
$$\s  C\ =\ {d\,\widehat{E}\over d\,T}\ =\ \beta^2\, \left( \,
      \langle E^2 \rangle - \langle E \rangle^2 \right)\ .$$
\begin{figure} \begin{center}
\includegraphics[height=0.64\textheight]{figures/2dI_C.pdf}
\end{center} \end{figure}

\newpage The energy histogram of the multicanonical simulation 
together its canonically re-weighted descendants at {\s $\beta=0$, 
$\beta = 0.2$ and $\beta =0.4$} (assignment {\tt a0501\_04}): 
\begin{figure} \begin{center}
\includegraphics[height=0.64\textheight]{figures/2dI_muh.pdf}
\end{center} \end{figure}
The normalization of the multicanonical histogram is adjusted so that
it fits into the same figure with the three re-weighted histograms.

\newpage\noindent Similar data for the $2d$ 10-state Potts model, 
action variable {\tt actm} data (assignment {\tt a0501\_05}):
\begin{figure} \begin{center}
\includegraphics[height=0.74\textheight]{figures/2d10q_act.pdf}
\end{center} \end{figure}

\newpage
The canonically re-weighted histogram at $\s \beta=0.71$ together 
with the multicanonical histogram using suitable normalizations
(the ordinate is on a logarithmic scale):
\begin{figure} \begin{center}
\includegraphics[height=0.74\textheight]{figures/2d10q_muh.pdf}
\end{center} \end{figure}

\newpage 
The multicanonical method allows to estimate the interface 
tension of the transition by following the minimum to maximum 
ratio over many orders of magnitude (Berg and Neuhaus 1992):
\vskip 0pt
\begin{figure} \begin{center}
\includegraphics[height=0.80\textheight]{figures/potts2d10q_ln.pdf}
\end{center} \end{figure}

\newpage
\begin{figure} \begin{center}
\includegraphics[height=0.84\textheight]{figures/92bn04.pdf}
\caption{Interface tension fit (Berg and Neuhaus 1992).  }
\end{center} \end{figure}

\foilhead{\r Free Energy and Entropy Calculation}\vskip -30pt

At $\s\beta=0$ the Potts partition function is $Z_0=q^N$. MUCA
simulations allow for proper normalization of the partition function, 
if $\s\beta=0$ is included in the temperature range. The properly 
normalized partition function allows to calculate the properly
normalized
$$\s {\rm\r Helmholtz\ free\ energy}~~~F\ =\ -\beta^{-1}\,\ln(Z)$$
and the
$$\s {\rm\r Entropy}~~~S\ =\ {F-E\over T}\ =\ \beta\,(F - {E})$$
of the canonical ensemble. Here $\s E$ is the expectation value of the 
internal energy and the last equal sign holds because of our choice of 
units for the temperature.  For the $2d$ Ising model as well as for the 
$2d$ 10-state Potts model, we show in the next figure multicanonical 
estimates of the free energy density per site
$$\s f\ =\ F/N\ .$$ \vskip -30pt

\begin{figure} \begin{center}
\includegraphics[height=0.7\textheight]{figures/free_energy.pdf}
\caption{Free energies from multicanonical simulations of the 
Ising and 10-state Potts models on an $20\times 20$ lattice
(assignments {\tt a0501\_03} and {\tt a0501\_05}). The full lines 
are the exact results of Ferdinand and Fischer for the Ising model 
and the asymptotic behavior 
$f_{\rm as} = {2\,d\over q} - 2\,d - \beta^{-1}\, {\ln(q)\over N}$
for the 10-state Potts model.}
\end{center} \end{figure}

\newpage
$$ {\rm Entropy\ Density}~~~\s s\ =\ S/N\ :$$ \vskip -30pt
\begin{figure} \begin{center}
\includegraphics[height=0.67\textheight]{figures/entropy.pdf}
\caption{Entropies from multicanonical simulations of the Ising 
and 10-state Potts models on an $20\times 20$ lattice 
(assignments {\tt a0501\_03} and {\tt a0501\_05}). The full 
line is the exact result of Ferdinand and Fischer for the 
Ising model.}  
\end{center} \end{figure}

\newpage
For the $2d$ Ising model one may also compare directly with the number
of states. Up to medium sized lattices this integer can be calculated 
to all digits by analytical methods (Beale). MC results are only 
sensitive to the first few (not more than six) digits and one finds 
no real advantages over using other physical quantities.

\begin{center}
{\Large\bf\r Time series analysis}
\end{center}

Typically, one prefers in continuous systems time series data over 
keeping histograms, because one avoids then dealing with the 
discretization errors.  Even in discrete systems time series data 
are of importance, as one often wants to measure more physical 
quantities than just the energy. Then the RAM storage requirements 
may require to use a time series instead of histograms. We illustrate 
this point, using the {\r Potts magnetization:}

In assignments {\tt a0501\_08} and {\tt a0501\_09} we create the same 
statistics on $\s 20\times 20$ lattices as before, including time series 
measurements for the energy and also for the Potts magnetization. For 
energy based observables the analysis of the histogram and the time 
series data give consistent results. 

\newpage \b
For zero magnetic field, $\s H=0$, the expectation value of the Potts 
magnetization on a finite lattice is 
$$\s M_{q0}\ =\ \langle\,\delta_{q_i,q_0}\,\rangle\ =\ \frac{1}{q}\ ,$$
independently of the temperature. For the multicanonical simulation 
each Potts state is visited with probability $\s 1/q$ also at low 
temperatures.  In contrast to this, the expectation 
value of the {\r magnetization squared}
$$\s M_{q0}^2\ =\ q\, \left\langle \left({1\over N}\,\sum_{i=1}^N 
     \delta_{q_i,q_0} \right)^2 \right\rangle $$
is a non-trivial quantity. At $\s \beta=0$ its $\s N\to\infty$ value 
is $\s M^2_{q0}=q\,(1/q)^2=1/q$, whereas it approaches $\s 1$ for 
$\s \beta\to\infty$. For $\s q=2$ and $\s q=10$ the next figure shows 
our numerical results and we see that the crossover of $\s M^2_{q0}$ 
from $\s 1/q$ to $\s 1$ happens in the neighborhood of the critical 
temperature. A FSS analysis would reveal that, a singularity 
develops at $\s \beta_c$, which is in the derivative of $\s M^2_{q0}$ 
for the second order phase transitions {\s ($\s q\le 4$)} and in 
$\s M^2_{q0}$ itself for the {first order transitions} \s ($q\ge 5$).

\begin{figure} \begin{center}
\includegraphics[height=0.7\textheight]{figures/Potts_M2.pdf}
\caption{The Potts magnetization per lattice site squared for the $q=2$ 
and $q=10$ Potts models on an $20\times 20$ lattice (assignments
{\tt a0501\_08} and {\tt a0501\_09}). }
\end{center} \end{figure}

\foilhead{\r\Large\bf Summary} \b

\begin{itemize}
\item We considered Statistics, Markov Chain Monte Carlo simulations, 
      the Statistical Analysis of Markov chain data and, finally, 
      Multicanonical Sampling.

\item It is a strength of computer simulations that one can generate
      artificial (not realized by nature) ensembles, which enhance the 
      probabilities of rare events on may be interested in, or speed up 
      the dynamics.

\item Nowadays {\r Generalized Ensembles} (umbrella, multicanonical, 1/k, 
      ...) have found many applications. Besides for first order phase 
      transitions in particular for {\r Complex Systems} such as 
      biomolecules, where they accelerate the dynamics.
\end{itemize}

\end{document}
