
\frame{\frametitle{\r Metropolis Algorithm and Hastings' Extension}\b
Markov Chain Monte Carlo (MCMC) simulations started in earnest
with the following famous paper: \vskip -12pt
\begin{figure} \begin{center}
 \includegraphics[width=1.1\textwidth]{figures/1953MetropolisAbstract.pdf}
\end{center} \end{figure} \vskip -12pt  
\begin{figure} \begin{center}
 \includegraphics[width=0.6\textwidth]{figures/Metropolis1.pdf}
% \includegraphics[width=0.8\textwidth]{figures/1953MetropolisText.pdf}
\end{center} \end{figure}   
\vskip 20pt }


\frame{\frametitle{\r }\b
\begin{figure} \begin{center}
 \includegraphics[width=0.8\textwidth]{figures/Metropolis2.pdf}
\end{center} \end{figure} \b 
Famous {\r accept/reject} criterion. \bigskip

When a move is rejected: Count present configuration again!
\vskip 30pt }


\frame{\frametitle{Notable early developments:}\b

{\r 1959} {\s Salsburg, Jacobson, Fickett, and Wood, J. Chem. Phys. 30, 
65--72:} An attempt to calculate the partition function by MCMC using a 
method called in the modern language {\r reweighting:} 
$$\s \exp(-\beta\,E)\ \to\ \exp(-\beta\,E-\triangle\beta\,E) =
   \exp(-\beta'\,E)\ . $$
As noticed by the authors their method is restricted to very small 
lattices. 
\bigskip

{\r 1963} {\s Glauber, J. Math. Phys. 4, 294-307:} discussion of 
heatbath and other dynamics for the Ising model.
\bigskip

{\r 1967} {\s McDonald and Singer, Discussions Faraday Soc. 43, 40--49:}
Used reweighting to evaluate physical quantities over a small range 

of temperatures.  \bigskip

{\r 1972} {\s Valleau and Card, J. Chem. Phys. 37, 5457--5462:}
{\r Multistage sampling} using overlapping bridging distributions.
\vskip 20pt }


\frame{\frametitle{\r }\b
Extension:
\begin{figure} \begin{center}
 \includegraphics[width=1.0\textwidth]{figures/1970HastingsAbstract.pdf}
\end{center} \end{figure} 
\vskip 30pt }


\frame{\frametitle{\r }\b
Arbitrary distributions $\pi$: 
\begin{figure} \begin{center} 
 \includegraphics[width=1.0\textwidth]{figures/1970HastingsTransition1.pdf}
\end{center} \end{figure}   
\vskip 30pt }

\frame{\frametitle{\r }\b
\begin{figure} \begin{center}
 \includegraphics[width=1.0\textwidth]{figures/1970HastingsTransition2.pdf}
\end{center} \end{figure}   
{\s Reversibility conditions:} Detailed balance.

{\s Irreducibility:} Ergodicity.

Normalization of the distribution $\s \pi_i$ drops out. \vskip 10pt

Metropolis updating is a special case: \s
\begin{equation}\s\nonumber
  s^M_{ij} = \begin{cases} 1+\frac{\pi_i q_{ij}}{\pi_j q_{ji}}~~
             {\rm for}~~\frac{\pi_j q_{ji}}{\pi_i p_{ij}}\ge 1~~
             \Rightarrow~~\alpha^M_{ij}=1\,,
             \cr           1+\frac{\pi_j q_{ji}}{\pi_i q_{ij}}~~
             {\rm for}~~\frac{\pi_j q_{ji}}{\pi_i p_{ij}}< 1~~
 \Rightarrow~~\alpha^M_{ij}=\frac{\pi_j q_{ji}}{\pi_i q_{ij}}\,,
 \end{cases}
\end{equation}\b
which goes for $\s q_{ij}\ne q_{ji}$ sometimes under the name {\r 
biased Metropolis} updating.
\vskip 40pt }


\frame{\frametitle{\r }\b

Symmetry: 
\begin{equation}\s\nonumber
  q_{ij}=q_{ji}~~\Rightarrow~~\alpha^M_{ij} = \begin{cases} 
  1~~{\rm for}~~\frac{\pi_j}{\pi_i}\ge 1\,,\cr
  \frac{\pi_j}{\pi_i}~~{\rm for}~~\frac{\pi_j}{\pi_i}<1\,. 
                                              \end{cases}
\end{equation}
Boltzmann distribution:
$$\s \pi_i = \exp(-\beta E_i)~~\Rightarrow~~{\rm\b Metropolis\ paper.}$$
\vskip 20pt }


\frame{\frametitle{\r Umbrella Sampling}\s
Torrie and Valleau, Chem. Phys. Lett. 28 (1974) 578: estimation of 
{\b free energies by sampling on distributions designed for this 
purpose (no reference to Hastings). Further Elaboration:}  % \vskip -6pt
\begin{figure} \begin{center}
\includegraphics[width=1.0\textwidth]{figures/PateyValleau1.pdf}
\end{center} \end{figure}  \vskip -6pt
Since we want to explore a wider region of ionic separations, we used
a modified importance sampling scheme. The configuration were chosen
with a frequency proportional to $w(r_{12})\,\exp(-U/kT)$, where the
{\r weighting function} $w(r_{12})$ was chosen to spread the sampling
over the desired range of $r_{12}$. ...
Accepted or rejected depending on the ratio
$$\s \frac{\pi'}{\pi} = \frac{w(r'_{12})}{w(r_{12})}\,
     \exp\left(-\frac{U'-U}{kT}\right)\, .$$
\vskip 20pt }


\frame{\frametitle{}\s 
J. Comp. Phys 23, 187--199 (1977): \vskip -6pt
\begin{figure} \begin{center}
\includegraphics[width=1.0\textwidth]{figures/TorrieValleau2.pdf}
\end{center} \end{figure}  \vskip -6pt
... statistical properties such as the entropy and free energy, because
they cannot be expressed as ensemble averages, have not been so easily
accessible. The conventional technique has been numerical integration,
...~. This somewhat cumbersome method is least efficient or altogether
unworkable when the system undergoes a phase transition, because of
the difficulty of defining a path of integration on which the necessary
ensemble averages can be reliable measured, ...
\vskip 20pt }


\frame{\frametitle{\r }\b
\begin{figure} \begin{center}
\includegraphics[width=1.0\textwidth]{figures/TorrieValleau3a.pdf}
\end{center} \end{figure} 
\begin{figure} \begin{center}
\includegraphics[width=1.0\textwidth]{figures/TorrieValleau3b.pdf}
\end{center} \end{figure} 
\vskip 20pt }


\frame{\frametitle{\r }\b
\begin{figure} \begin{center}
\includegraphics[width=1.0\textwidth]{figures/TorrieValleau4.pdf}
\end{center} \end{figure}
Summary:
\begin{enumerate}
\item Weighting functions $\r w(q_1,\dots,q_N)$ depend on
      (generalized) coordinates.
\item Determination via trial and error.
\end{enumerate}

\b Conclusion eleven years later \s (Li and Scheraga, J. Mol. Struct. 
(Theochem) 179, 333-352 (1988)): \medskip

\r The difficulty of finding such 
weighting factors has prevented wide applications of umbrella sampling.
\vskip 20pt }


\frame{\frametitle{\b Umbrella Sampling in Practice } \s % \vskip -5pt
J. Am. Chem. Soc. 107, 154--163 (1985):
\begin{figure} \begin{center}
\includegraphics[width=1.0\textwidth]{figures/JACS_UMBRELLA1.pdf}
\end{center} \end{figure} \vskip -13pt
\begin{figure} \begin{center}
\includegraphics[width=0.9\textwidth]{figures/JACS_UMBRELLA2.pdf}
\end{center} \end{figure} 
\vskip 20pt }


\frame{\frametitle{\r }\b \vskip -20pt
\begin{figure} \begin{center}
\includegraphics[width=0.7\textwidth]{figures/JACS_UMBRELLA3a.pdf}
\end{center} \end{figure} \vskip -12pt
\begin{figure} \begin{center}
\includegraphics[width=0.9\textwidth]{figures/JACS_UMBRELLA3b.pdf}
\end{center} \end{figure} 
\vskip 20pt }

\frame{\frametitle{\r }

{\b \bf \Large

Summary: \vskip 20pt

Tedious patching of umbrella potentials.} \vskip 20pt

Let's now backup a bit in time to follow another  line of
developments:

\vskip 20pt }


\frame{\frametitle{\r Binder's Method for Estimating Interface Tensions}
Binder, Phys. Rev. A 25, 1699--1709 (1982). \vskip 15pt
\b \centerline{Simulations with periodic boundary conditions.} \vskip -5pt
$${\rm From\ equal\ heights\ double\ peaks:}~~~~\s F^s_L\ =\ -
\frac{1}{L^{(D-1)}}\,\ln\left(\frac{P_L^{\min}}{P_L^{\max}}\right)\,.$$
\begin{figure} \begin{center}
 \includegraphics[width=0.7\textwidth]{figures/BinderPRA1982.jpg} 
\end{center} \end{figure}   
First results with canonical simulations remained pitiful. Reason: 
$\s P_L^{\min}$ is exponentially suppressed in the canonical ensemble.
\vskip 20pt }


\frame{\frametitle{\r } \r 
Chandler: Introduction to Modern Statistical Mechanics, Oxford 
University Press 1987. \b Chapter 6.3: Non-Boltzmann Sampling:
\bigskip \s

Non-Boltzmann sampling is a powerful tool when the reference or
unperturbed energy, $E^{(0)}_{\nu}$, creates a  trajectory that 
is close to that of $E_{\nu}$. 
\r Comment: \b There is no need for the trajectories to stay close.  
\bigskip \s

Non-Boltzmann sampling can also be useful in removing the bottlenecks
that create quasi-ergodic problems and in focusing attention on rare
events. 
\r Comment: \b This is on the right track.
\bigskip

\r The practice: \s To illustrate this methodology, we consider ...
the computation of the free energy function, $\tilde{A}(M)$. This
function is defined for the Ising magnet by
\begin{equation} \nonumber
  \exp\left[-\beta\tilde{A}(M)\right] = \sum_v \Delta\left(M-
  \mu\sum_{i=1}^Ns_i\right)\,\exp\left(-\beta E_{\nu}\right)\,,
\end{equation}
where $\Delta(x)$ is the Kronecker delta ...~. 
\vskip 20pt

Patching of of weighting factors from umbrella or window potentials 
is then done (e.g., Chandler 1987).  But applications remained rather 
limited, although MCMC simulations flourished in the 1980s.

Furthermore, there appeared to be no focus on identifying challenging
problems, which could be overcome by using umbrella sampling.
\vskip 20pt }


\frame{\frametitle{ }  
Clearly,
\begin{equation} \nonumber
  \exp\left[-\beta\tilde{A}(M)\right] \propto P(M) =
  \left\langle \Delta\left(M-\mu\sum_{i=1}^Ns_i\right)\right\rangle\,,
\end{equation}
where $P(M)$ is the probability of observing the Ising magnet with
magnetization $M$. ... if we consider the situation of broken symmetry
(i.e., when $T<T_c$), and plan to compute $\tilde{A}(M)$ for a wide
range of $M$ values, we immediately encounter a serious problem ...
$\tilde{A}(M)$ is a bistable function of $M$, ... visitation of
states with $M=0$ is an infrequent event ...

The method of umbrella sampling, however, avoids this difficulty. 
We chose a set of umbrella or window potentials
\begin{eqnarray} \nonumber
  W_{\nu} &=& 0,~~{\rm for}\ M_i-w/2\le\mu\sum_{j=1}^Ns_j\le 
                   M_i\le M-i+w/2 \\ \nonumber
          &=& \infty,~~{\rm otherwise.}
\end{eqnarray}
... $\tilde{A}(M)$ is determined in each window to within an additive 
constant ... , which must be adjusted from one window to the next, ...
\vskip 20pt }


\frame{\frametitle{ } 
Let $\tau$ denote the required computer time to acquire such (accurate)
statistics in each window. ... Note that
\begin{equation} \nonumber
  \tau \propto w^2\,.
\end{equation}
Therefore, the total computation time required to determine $\tilde{A}
(M)$ by the method of umbrella sampling is
\begin{equation} \nonumber
  t_{\rm CPU} \propto n w^2\,.
\end{equation}
Now, how much time would it have taken if we did not use this method?
As a lower bound, let us assume that $\tilde{A}(M)$ does not vary
more than a few $k_BT$ over the entire range of $M$. The size of 
this range is $nw$. Therefore, the time to sample this range is
proportional to $(nw)^2=n\,t_{\rm CPU}$. Hence, without the windows,
the computation time would be $n$ times longer than that with the
$n$ windows. \bigskip

{\r Exercise BB1:} \b What is wrong with Chandler's argument? 
\smallskip \s

(Replacing Exercise 6.8 of the book.)
\vskip 20pt }


\frame{\frametitle{ } 
The advantage (i.e., lower computation time) of umbrella sampling is,
of course, even greater than this when regions of $M$ have relatively
high values of $\tilde{A}(M)$ and thus relatively low probabilities.
\smallskip \r

Comment: \b This is the (only) major advantage! \bigskip \s

400 spin Ising magnet, ... 10 windows were used between $M=0$ and
$M=400\mu$. From the figure (next transparency): 
$$ \frac{P^{\min}}{P^{\max}}\approx e^{-14} = 8.3\times 10^{-7}\,.$$ 

\r Comment: \b Configurations sampled, which are in Boltzmann
simulations suppressed by \r 6 to 7 orders of magnitude. \b
\bigskip

\r No mention of Binder \b or any attempt to calculate the interface
tension (known exactly since Onsager's 1944 solution of the

2D Ising model). Umbrella potentials by Tedious patching of 
windows remains (in Chandler's example chosen to be constant).
\vskip 20pt }


\frame{\frametitle{\r Multicanonical and Multimagnetical Simulation} 

\b Also in 1987: Umbrella sampling in Chandler's version was 
discovered independently by Bhanot and collaborators, and used 
in a series of papers: \s Phys. Rev. Lett. 59 (1987) 803 and 
references given therein.
\bigskip \b

Due to the tedious patching of windows their simulations remained 
confined to rather small lattices.
\bigskip

This changed with
\begin{figure} \begin{center}
 \includegraphics[width=1.1\textwidth]{figures/BergNeuhausMUCA1.pdf} 
\end{center} \end{figure}   
\vskip 20pt }

\frame{\frametitle{\r } 
\begin{figure} \begin{center}
 \includegraphics[width=1.0\textwidth]{figures/BergNeuhausMUCA2.pdf} 
\end{center} \end{figure}   
\vskip 20pt }

\frame{\frametitle{\r } 
... In practice the appropriate choice of the parameters in Eq.~(7) 
is obtained by making from the given systems an FSS prediction of the
density distribution $P_L(S)$ for the next larger system. A second
run may then be performed with optimized parameters [Citation End].
\bigskip \r

Changes compared to umbrella sampling:
\bigskip \b

No windows anymore, all in one ensemble.
\bigskip \b

Weights are not umbrella potentials $\s w(q_1,\dots,q_N)$, but 
functions of the energy: \s $w(S)$, here $S$ energy.
\bigskip \b

Weights are obtained in one (plus one) step(s).
\vskip 20pt }

\frame{\frametitle{\r } \b
\b \centerline{2D 10-state Potts model:} \s
\begin{figure} \begin{center}
\includegraphics[width=0.6\textwidth]{figures/92bn01.pdf} 
\end{center} \end{figure}   
\vskip 20pt }

\frame{ \vskip -30pt
\begin{figure} \begin{center}
\includegraphics[width=0.9\textwidth]{figures/color_potts2.jpg} 
% \caption{Data/C00.plt}
\end{center} \end{figure}   
}

\frame{ \vskip -30pt
\begin{figure} \begin{center}
\includegraphics[width=0.9\textwidth]{figures/color_potts3.jpg} 
% \caption{Data/C00.plt}
\end{center} \end{figure}   
}

\frame{\frametitle{\r Energy per spin histogram}  \b
\centerline{allows then for accurate estimates of interface tension}
\vskip -75pt \s
\begin{figure} \begin{center}
\includegraphics[width=0.8\textwidth]{figures/potts2d10q_ln.pdf} 
\end{center} \end{figure} \b \vskip -20pt
\vskip 20pt }


\frame{\frametitle{\b Back to Chandler's Case} \b
\begin{figure} \begin{center}
 \includegraphics[width=1.1\textwidth]{figures/BHN_MUMA1.pdf} 
\end{center} \end{figure}   
\vskip 20pt }

\frame{\frametitle{ } \s
\begin{figure} \begin{center}
 \includegraphics[width=1.0\textwidth]{figures/BHN_MUMA2.pdf} 
\end{center} \end{figure} \vskip -10pt
\centerline{ $L= 10,\, 20,\, 30,\, 40,\, 50,\, 60,\, 74,\, 80,\, 100.$}
\vskip 20pt }

\frame{\frametitle{ } \s
\begin{figure} \begin{center}
 \includegraphics[width=1.1\textwidth]{figures/BHN_MUMA3.pdf} 
\end{center} \end{figure}   
\vskip 20pt }

\frame{\frametitle{\b 3D Ising model} \s
Berg, Hansmann, and Neuhaus, Z. Phys. B 90 (1993) 229-239 \vskip -60pt
\begin{figure} \begin{center}
\includegraphics[width=0.8\textwidth]{figures/I3d_fs.pdf} 
\end{center} \end{figure}  \vskip -85pt 
\rightline{\r $10^{60}$ improvement over simulation with 
Boltzmann weights!$\qquad~$}
\vskip 20pt }


\frame{\frametitle{\r New Horizons}
\begin{figure} \begin{center}
 \includegraphics[width=1.1\textwidth]{figures/BergCelik1.pdf} 
\end{center} \end{figure}   
We present a {\b recursive procedure} to calculate the parameters of 

the recently introduced multicanonical ensemble and explore the 
approach for spin glasses. \bigskip

We suggest that in a large {\b class of situations,} in particular 
those where canonical simulations face severe ergodicity problems, it 
is more efficient to reconstruct the Gibbs ensemble form a simulation 

of a multicanonical ensemble [7] than to simulate is directly.
\bigskip

\centerline{[Citation End].}
\vskip 20pt }

\frame{\frametitle{\b Illustration of Rugged Free Energy Landscapes} \b
\begin{figure} \begin{center}
\includegraphics[width=0.8\textwidth]{figures/color_ising_mag.jpg} 
\end{center} \end{figure}  \s \vskip -2pt
Distribution of the 2D Ising Model Magnetization versus Temperature.
\vskip 20pt }

\frame{\frametitle{\b One Edwards-Anderson Ising (EAI) Spin Glass 
Realization:}
\s Parisi order parameter distribution versus temperature. 
\begin{figure} \begin{center}
\includegraphics[width=0.9\textwidth]{figures/color_spg.jpg} 
\end{center} \end{figure}   
\vskip 20pt }

\frame{\frametitle{}

\r Exercise BB2: \b Explain why patching windows cannot work well for
exploring these free energy landscapes. \bigskip

Berg-Celik concluding remarks:
\begin{figure} \begin{center}
\includegraphics[width=1.0\textwidth]{figures/BergCelik2.pdf} 
\end{center} \end{figure}   
J\"ulich Conference, \s Berg, Int. J. Mod. Phys. C 3 (1992) 1083-1098:
\b Review up to then.
% , weights in the form \s
% \begin{equation} \nonumber
%   w_i^{mu} = N^{-1} n_i^{-1} ~~{\rm with}~~ n_i=n(E_i)
% \end{equation}
\vskip 20pt }


\frame{\frametitle{\r }
\begin{figure} \begin{center}
\includegraphics[width=1.0\textwidth]{figures/1993HansmannOkamoto.pdf} 
\end{center} \end{figure}   
\bigskip \b

Their review {\s The Generalized-Ensemble Approach for Protein Folding
Simulations, Ann. Rev. Comp. Phys. 6, 129-157 (1999)} coined the name
\r Generalized Ensembles.  \bigskip
\vskip 20pt }


\frame{\frametitle{\r Poly-alanine:} \b % (also animation)} \b 
\vskip -10pt
\begin{figure} \begin{center}
\includegraphics[width=0.8\textwidth]{figures/ala_hel.jpg} 
\end{center} \end{figure}  \s
Hansmann and Okamoto, J. Chem. Phys. 110, 1267-1276 (1999);
erratum 111, 1339.
\vskip 20pt }


\frame{\frametitle{MUCA Performance} \b
{\Large Slowing down in \r units of updates:} \bigskip
\begin{enumerate} \b
  \item Optimum: $\s\propto N^2$ perfect random walk in the energy.
  \item 2D 10-state Potts model (transition range): Effective $\s 
        \propto N^{2.325\,(10)}$ observed (Berg and Neuhaus 1992). 
  \item 2D EAI spin glass: $\s \propto N^{3.2\,(2)}$ observed (Berg
        and Celik 1992). 
\end{enumerate}
\b Still exponential slowing down expected due to \r hidden energy 
barriers. \b Even for simple first order phase transitions: \s Hager 
and Neuhaus, J. Stat. Phys. 113 (2003) 47--83: $~~\Leftrightarrow$
\begin{figure} \centering
\includegraphics[height=1.5in,width=1.5in]{figures/NHager1.pdf}\hfill
\includegraphics[height=1.5in,width=1.5in]{figures/NHager.pdf}\hskip -0.07in
\includegraphics[height=1.5in,width=1.5in,rotate=0]{figures/NHager2.pdf}
\caption{ lattice (right).} 
\end{figure}
% Transitions only through rare or suppressed configurations.
\vskip 20pt }
