\documentclass{beamer} % BB Nov 13 2006.
% Figures are with the Boston GEnsembles seminar.

\usetheme{fsuheplat}

%AB: color definitions
\definecolor{abtitlecolor}{rgb}{0.0,0.255,0.494}
\definecolor{absecondarycolor}{rgb}{0.0,0.416,0.804}
\definecolor{abprimarycolor}{rgb}{1.0,0.686,0.0}

%AB: redefinition of beamer colors
\setbeamercolor{title}{fg=abtitlecolor}
\setbeamercolor{frametitle}{fg=abtitlecolor}
\setbeamercolor{palette tertiary}{fg=white,bg=abtitlecolor}
\setbeamercolor{palette secondary}{fg=white,bg=absecondarycolor}
\setbeamercolor{palette primary}{fg=black,bg=abprimarycolor}
\setbeamercolor{structure}{fg=abtitlecolor}

\setbeamerfont{section in toc}{series=\bfseries}
\def\b{\color{blue}}
\def\r{\color{red}}
\def\s{\color{black}}

%AB: remove navigation icons
\beamertemplatenavigationsymbolsempty


\title[Generalized Ensembles]{\textbf{Lecgture 10: Computer Simulations 
of Generalized Ensembles}} % From Fall 2006 Boston lecture.

\author{Bernd A. Berg}

\institute[FSU]{Florida State University}

\date{\today}

\begin{document}

\frame{\titlepage}

\frame{   \frametitle{\r Overview} \b
\begin{enumerate} \b 
\item Reweighting
\item Umbrella Sampling
\item Binder's Method for Estimating Interface Tensions
\item Multicanonical Simulations
  \begin{itemize}
  \item First order phase transitions
  \item Groundstates and rugged free energy landscapes
  \item Simulations of peptides
  \item How to Get the Weights (Wang-Landau Recursion)?
  \item MUCA Performance 
  \item Second order phase transitions and cluster algorithms
  \end{itemize}
% \item Systems with Rugged Free Energy Landscapes
\item Replica Exchange Method (Parallel Tempering)
\end{enumerate}
\bigskip

\hfil\break \vfil
}

\frame{\frametitle{\r Reweighting -- Single Histogram Method}\b
\vskip-10pt
$$\s \exp(-\beta\,E)\ \to\ \exp(-\beta\,E-\triangle\beta\,E) =
   \exp(-\beta'\,E) $$ \vskip -8pt

{\r 1959:} A first attempt to calculate the partition function by MCMC 
is made in a paper by Salsburg, Jacobson, Fickett, and Wood, using a 
method called in the modern language {\r reweighting.} As noticed by 
the authors their method is restricted to very small lattices. 
\smallskip

{\r 1967:} McDonald and Singer used reweighting to evaluate physical 
quantities over a small range of temperatures. Thereafter the approach 
appeared to be dormant.
\smallskip

{\r 1982:} Reweighting was rediscovered by Falcioni, Marinari, 
Paciello, {\r Parisi} and Taglienti who focused on calculating 
complex zeros of the partition function.
\smallskip

{\r 1988:} {\r Ferrenberg and Swendsen} formulated crystal-clear for 
what the method is good, and for what it is not: Reweighting allows to 
focus on maxima or minima of appropriate observables, but not to cover 
a finite temperature range in the infinite volume limit.
\smallskip
\vskip 20pt }

\frame{\frametitle{\r Reweighting -- Multi Histogram Method}\b
% \vskip-10pt

{\r 1972:} Valleau and Card introduced the use of overlapping bridging 
distributions and called their method {\r multistage sampling.}
\smallskip

{\r 1989:} Ferrenberg and Swendsen proposed a recursion for 
creating a joint distribution.  \smallskip

{\r 1990:} Alves, Berg and Villanova patched histograms by minimizing 
$\s\chi^2$ for the reweighted overlap: $\s\nu=0.6303\,(14)$ for 3D 
Ising model ($\s\nu=0.6301\,(4)$ nowadays best, Pelissetto~\& Vicari 2002). 
$\s 14^3$ lattice: \vskip -8pt
\begin{figure} \begin{center}
\includegraphics[width=0.6\textwidth]{figures/AlvesPRB1990.jpg} 
% \caption{Data/C00.plt}
\end{center} \end{figure}   
\vskip 20pt }

\frame{\frametitle{\r Umbrella Sampling}\b
Patching histograms of canonical simulations faces a number of limitations:
\begin{enumerate} \b
\item Bottlenecks in the canonical ensemble can create quasi-ergodicity.
\item There is no way to focus on canonically rare events.
\end{enumerate} \medskip

{\r 1976: Torrie and Valleau} introduced umbrella sampling, which 
copes with these difficulties by allowing arbitrary sampling 
distributions.  \medskip

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.
 \medskip

{\r Li and Scheraga, 1988:} {\it\s The difficulty of finding such 
weighting factors has prevented wide applications of umbrella 
sampling.}\medskip

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

\frame{\frametitle{\r Binder's Method for Estimating Interface Tensions}
\b \centerline{Simulations with periodic boundary conditions (Binder 1982).}
$${\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}   
But in practice results remained pitiful: $\s L^{\max}=12$ for the 
3D Ising model with far off estimates like $\s F^s_{\infty}= 0.0050\, 
(25)$ at $\s\beta=0.27227$. Reason: $\s P_L^{\min}$ is exponentially
suppressed in the canonical ensemble.
\vskip 20pt }

\frame{\frametitle{\r Multicanonical Simulations 
(Berg and Neuhaus 1991/1992)} \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{\b $\dots$ achieved by sampling in an appropriate energy range 
$\s E^{\min}\le E\le E^{\max}$ with an {\r approximation of the weights}
$$\s W_{\rm MUCA}(E) = \frac{1}{n(E)} = e^{-b(E)\,E+\alpha(E)} $$
and reweighting to the canonical ensemble. Thus, the multicanonical 
approach requires two steps (to be detailed later):
\begin{enumerate}\b
\item Obtain a {\r working estimate} of the weights $\s W_{\rm MUCA}(E).$
\item Perform a MCMC simulation with {\r fixed} weights.
\end{enumerate}
Working estimate means that the Markov chain {\r cycles (tunnels)}
in the energy range:
$$\s E^{\min} \le E \le E^{\max} $$
which corresponds to a {\r canonical temperature range (multicanonical)}
$$\s \beta^{\min}=b(E^{\max})\le \beta\le b(E^{\min})=\beta^{\max}\,.$$
Rare configurations are now sampled (2D 10-state Potts model):
\vskip 20pt }

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

\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} \vskip -60pt
\begin{figure} \begin{center}
\includegraphics[width=0.8\textwidth]{figures/potts2d10q_ln.pdf} 
% \caption{Data/C00.plt}
\end{center} \end{figure}   
\vskip 20pt }

\frame{
\begin{figure} \begin{center}
\includegraphics[width=0.6\textwidth]{figures/92bn04.pdf} 
\end{center} \end{figure}   
}

\frame{\frametitle{\r Potts models: infinite volume interface tensions 
estimates.} \b \bigskip

MUCA:
\centerline{\s $q=10:~~F^s_{\infty} = 0.09781\ (75)$}
\centerline{(Berg and Neuhaus 1992, submitted July 1991)}
\smallskip

\centerline{\s $q=7:~~~F^s_{\infty} = 0.0241\,(10)$ }
\centerline{(Janke, Berg, and Katoot 1992, submitted March 1992)}
\bigskip

Before:
\centerline{$\s q=7:~~~F^s_{\infty} = 0.1886\,(12)~~$ (Boston 1989)} 
\centerline{$\s q=7:~~~F^s_{\infty} \approx 0.20~~~~~~~~~$ 
            (Helsinki 1989)}
\bigskip

{\r Exact:} (Borgs and Janke 1992, submitted September 1992)\medskip

\centerline{\s $q=10:~~F^s_{\infty} = 0.094701\dots $ }
\centerline{\s $q=7:~~~F^s_{\infty} = 0.020792\dots $ }
\vskip 20pt }

\frame{\frametitle{\r Multimagnetical Simulations -- Back to the
3D Ising Model} \b
\centerline{(Berg, Hansmann, and Neuhaus, 1993)}\s \vskip -60pt
\begin{figure} \begin{center}
\includegraphics[width=0.8\textwidth]{figures/I3d_fs.pdf} 
\end{center} \end{figure}  \vskip -95pt 
\rightline{\r $10^{60}$ improvement in statistics!$\qquad~$}
\vskip 20pt }

\frame{\vskip 5pt $$\s F^s_{\infty} = 0.01293\ (17)
                  ~~~{\rm\b at}~~~ \beta = 0.272 $$ \vskip -10pt
\begin{figure} \begin{center}
\includegraphics[width=0.7\textwidth]{figures/Fig6BHN_ZPB1993.jpg} 
\end{center} \end{figure}  
\vskip 20pt }

\frame{\frametitle{\r Rugged Free Energy Landscapes and Groundstate
Entropy} \b
\centerline{(Berg and Celik 1992)} \s
\begin{figure} \begin{center}
\includegraphics[width=0.8\textwidth]{figures/color_ising_mag.jpg} 
\end{center} \end{figure}  \b \vskip -2pt
Distribution of the 2D Ising Model Magnetization versus Temperature.
\vskip 20pt }

\frame{\frametitle{\r Edward-Anderson Ising Spin Glass:} \b
$$\s E = - \sum_{\langle ij\rangle} J_{ij}\,s_i\,s_j\,,~~ 
     \s_i = \pm 1\,. $$
The exchange coupling constants $\s J_{ij}$ are quenched random 
variables. MUCA  simulations allow to equilibrate in the low 
temperature region. One averages over many $\s J_{ij}$ realizations 
(each requires to similate an Ising model with random couplings).  
Groundstates can be sampled 

on small systems (infinite volume 
extrapolations: Berg, Celik, and Hansmann 1994). \medskip

The magnetization is no longer an order parameter. 
The substitute is the {\r Parisi overlap} parameter:
$$\s q = \sum_i s_i^1\,s_i^2 $$ 
The superscripts $\s 1$ and $\s 2$ label {\r replica} of the same 
realization.
\vskip 20pt }

\frame{\frametitle{\r For one realization:}
\b 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 Residual Entropy of Ordinary Ice:}\b
\centerline{(Berg, Muguruma, and Okamoto 2006)} \s
\begin{figure} \begin{center}
\includegraphics[width=0.8\textwidth]{figures/icepntcl.pdf} 
\end{center} \end{figure}   
\centerline{\b Ice rule: All links between H$_2$O molecules are 
              hydrogen bonds.}
\vskip 20pt }

\frame{\vskip 10pt

{\Large\r Residual entropy per molecule:}\b \bigskip

$\s S^0 = \ln W_1$ where $\s W=(W_1)^N$ is the total number of 
configurations, which fulfill the ice rule. \medskip

{\r Ising-like model
(Linus Pauling):} Hydrogen can take two positions on the bond 
$\s W_0=2^{2N}$. Without correlations: $\s 6/16$ probability 
of correct hydrogen positions at a molecule and (1935)
$$\s W^{\rm Pauling} = \left(\frac{6}{16}\right)^N 4^N\ 
     \Rightarrow\ W_1^{\rm Pauling} = \frac{3}{2}\ .$$
Onsager and Dupuis showed that this is a lower bound. Onsager's
student Nagle devoted his Ph.D. thesis to improve Pauling's 
estimate by a series expansion method and (1966): 
$\s W_1^{\rm Nagle} = 1.50684\, (15).$ \medskip

Experimentalists (1936): $\s W_1^{\rm meas} = 1.507\, (20)$
(reduced error bars). \medskip

{\r Notes:} The residual entropy of ice was experimentally discovered 
by Giauque and Ashley before the hydrogen bond picture was fully 
developed (1933). H-bond clusters are important in water at room 
temperature.  \vskip 20pt }

\frame{\vskip 10pt

{\Large\r MUCA simulation: $\s W_1^{\rm MUCA}=1.50738\,(16).$}
\s \bigskip
\begin{figure} \begin{center}
\includegraphics[width=0.7\textwidth]{figures/fitpl.pdf} 
\end{center} \end{figure}  \b
\centerline{Theoretical prediction with an accurary better than 
one per mille!} \medskip

\centerline{\r Experimental efforts to confirm the prediction!?}
\vskip 20pt }

\frame{\frametitle{\r Multicanonical Simulations of Peptides} \b
\bigskip

{\r 1993 Hansmann and Okamoto:} Predictions of peptide conformation
by multicanonical algorithm (Met-Enkephalin).
\bigskip

{\r 1994 Hao and Scheraga:} Monte Carlo simulations of a first-order 
transition for protein folding using {\r entropic sampling} (lattice
protein). 
\bigskip

Entropic sampling is simply another name for the multicanonical 
method (Berg, Hansmann, and Okamoto 1995). {\r Broad or flat 
histogram, density of states} methods are other names in use. 

Still others names are presumably in the making ...
\vskip 20pt }

\frame{\frametitle{\r Poly-alanine: (also animation)} \b
\begin{figure} \begin{center}
\includegraphics[width=0.8\textwidth]{figures/ala_hel.jpg} 
\end{center} \end{figure}  \b
\centerline{(Hansmann and Okamoto 1999)}
\vskip 20pt }


\frame{\frametitle{\r How to get the Weights?} \b

\begin{enumerate} \b
\item Overlapping constrained (microcanonical) MC simulations.
      Problems: Tedious and ergodicity of constrained simulations 
      likely to break down in a rugged free energy landscapes.
\item {\r Finite Size Scaling (FSS) Estimates.} Best when it works!
      Problem: There may be no FSS theory.
\item {\r Recursions.} In general most convenient. Problem: 
      May deteriorate quickly with increasing system size. 
\end{enumerate}
\medskip

Example: {\r Multicanonical Recursion} (Berg 1996). Iterates
{\r microcanonical temperature} $\s b(E)$. Animation for the {\s 2D 
10-state} Potts model on an $\s 80\times 80$ lattice (100 iterations 
of 100 sweeps each).  \vskip 20pt }


\frame{\frametitle{\r Wang-Landau Recursion (Wang and Landau 2001)}\b
\medskip

Updates are performed with estimators $\s g(E)$ of the density of states:
\begin{equation}\nonumber \s 
  p(E_1\to E_2) = \min \left[ {g(E_1)\over g(E_2)}, 1\right]\ .
\end{equation}
Each time an energy level is visited, the estimator is updated
{\r multiplicatively:}
\begin{equation}\nonumber \s g(E) \to g(E)\,f\,. \end{equation}
Initially $\s g(E)=1$ and $\s f=f_0=e^1$. Once the desired energy range 
is adequately covered, the factor $\s f$ is refined
\begin{equation}\nonumber \s 
  f_1=\sqrt{f},\ f_{n+1}=\sqrt{f_{n+1}}
\end{equation}
until a value sufficiently close to  $\s 1.$ is reached. Means: The 
system keeps on {\r cycling} with frozen weights. \medskip

\centerline{Switch to the usual \r MUCA production run as soon as 
possible!} \medskip

\centerline{(To the contrary, Wang-Landau advertise to keep on iterating
forever.)} \vskip 20pt }

\frame{\frametitle{\r MUCA Performance} \b
{\Large Slowing down in units of sweeps} ($\s L^D$ updates): \bigskip
\begin{enumerate} \b
  \item Optimum: $\s\propto L^D$ perfect random walk in the energy.
        Compare to $\s\propto \exp(+f^s\,L^{D-1})$ canonical slowing
        down for first order phase transitions.
  \item 10-state Potts model (transition range): Effective $\s 
        \propto L^{2.65\,(2)}$ observed (Berg and Neuhaus 1992). 
        Expectation for large {\s $L$:} Subleading exponential 
        (Neuhaus and Hager 2003).
  \item EAI spin glass: $\s \propto L^{4.4\,(3)}$ observed (Berg
        and Celik 1992). Exponential expected for large $\s L$.
  \item 2D Ising model (entire energy range): $\s\propto L^{2.8}$ 
        improved to $\s L^2\,\ln(L)$ by optimizing cycling times 
        (Trebst, Huse, and \break Troyer 2004).
\end{enumerate}

\vskip 20pt }

\frame{\frametitle{\r Second order phase transitions and cluster 
algorithms} \b 
\centerline{(Berg and Janke 2006)} \s \vskip -8pt
\begin{figure} \begin{center}
\includegraphics[width=0.75\textwidth]{figures/rwght80.pdf} 
\end{center} \end{figure}  \b
$\s L^3,\ L=80^3$ Ising model at $\s\beta_c:~$ Desired reweighting range 

(entire axis) versus actual canonical reweighting range (rwght).
\vskip 20pt }

\frame{\b\vskip 10pt
Conventional Wang-Landau/MUCA simulations lack the efficiency 

of {\r cluster algorithms (Swendsen and Wang 1987, Wolff 1989).} 
\medskip 

Combination is possible using bonds instead of energy (MUBO). \s
\begin{figure} \begin{center}
\includegraphics[width=0.70\textwidth]{figures/tau3DI.pdf} 
\end{center} \end{figure}  \b
For 3D Ising: $~~\s \tau_{\rm int}^{\rm MUCA}\propto L^{2.22\,(11)}
~~{\rm and}~~\tau_{\rm int}^{\rm MUBO}=L^{1.05\,(5)}\,. $
\vskip 20pt }


\frame{\frametitle{\r Replica Exchange Method (Parallel Tempering)}\b 
\medskip

{\r 1988} Swendsen and Wang: (Too) General replica exchange method.
\medskip

{\r 1990} Frantz, Freemann, and Doll: {\r J}ump{\r -walking} feeds 
replica from a high into a low temperature simulation. Does not fulfill 
balance.
\medskip

{\r 1991 Geyer:} Multiple Markov chains (later called replica exchange).
\medskip

{\r 1992} Lyubartsev, Martsinovski, Shevkanov, and Vorontsov- Velyaminov:
Expanded ensembles (no exchange of replica).
\medskip

{\r 1992} Marinari and Parisi: Simulated tempering, which is
a special case of the Russian method (coined the tempering
notation).
\medskip

{\r 1996} Hukusima and Nemoto: Replica Exchange, often called 
parallel tempering (when temperatures are exchanged).
\medskip 

{\r 1997} Hansmann: First application to biomolecules (Met-Enkephalin).
\vskip 20pt }


\frame{\frametitle{\r Parallel Tempering Algorithm}\b 
For $\s \beta_1 < \beta_2 < ... < \beta_n$ consider the {\r joint 
partition function} 
\begin{equation}\s\nonumber 
  Z = \prod_{k=1}^n \left\lbrace \prod_{i=1}^N 
      \frac{(2\pi\,m_i\,k_B\,T_k)^{3/2} }{N!} 
      \int d^3 x_i^k\, \exp\left[ -\beta_k\,
      U(\vec{r}_1^{\,k},\dots,\vec{r}_N^{\,k})\right]
      \right\rbrace 
\end{equation}
of a molecular system, where the momenta are integrated out. {\r Swaps}
between (normally neighbouring) temperatures are proposed with uniform 
probability and accepted with the Metropolis probability
\begin{equation}\s\nonumber 
  P_{\rm acpt} = \min \left\lbrace 1,\,\exp \left[+(\beta_k-\beta_l)\,
         \left( U(\vec{r}_1^{\,k},\dots,\vec{r}_N^{\,k}) -
                U(\vec{r}_1^{\,l},\dots,\vec{r}_N^{\,l})\right)
                \right]\right\rbrace
\end{equation}
so that {\r balance holds.} Refreshes low temperature simulations
through high temperature excursions. {\r Refinements:} \medskip

\centerline{Other variables than temperatures can be exchanged as well.}
\medskip

\centerline{Generalized ensembles can be exchanged, ... }
\medskip

\medskip \vskip 20pt }


\frame{\frametitle{\r Temperature Exchange}\b 
\begin{figure} \begin{center}
\includegraphics[width=0.75\textwidth]{figures/ParallelTempering1.jpg} 
\end{center} \end{figure}  \b
\centerline{Typical time series.}
\medskip \vskip 20pt }


\frame{\frametitle{\r Energy Histograms}\b 
\begin{figure} \begin{center}
\includegraphics[width=0.75\textwidth]{figures/ParallelTempering2.jpg} 
\end{center} \end{figure}  \b
Overlap of energy histograms between adjacent replicas allows for 
acceptance of temperature swaps.
\medskip \vskip 20pt }

\end{document}

