Histogram Methods in Monte Carlo Simulations

Belliappa
IMSc

Sep 3, 2010

Preliminaries and Notation

The density of states, the number of microstates available at an energy E, is (E )

If we choose a state at random from the density of states, the probability Pr of picking any one microstate r is Pr = 1 total

Thus the probability that the microstate weve picked has energy E is P(E ) =
r |Er =E

1 total (E ) total

In a canonical ensemble the probability Pr of a microstate r of A with energy Er is no longer equal for all r but depends on the temperature : Pr () e Er e Er = Er r e where the denominator is referred to as the partition function Z ().

Thus, the probability that the microstate weve picked has energy E , PE (), is (E ) e E E E (E ) e

PE () =

Sampling from a canonical ensemble

How do we get the canonical distribution? the energy of a microstate of the Ising model is J <i,j> si sj with J = -1 in all graphs shown below try generating microstates at random and binning them according to their energies to get the density of states (E ) but this doesnt work in practice, because ...

Sampling from a canonical ensemble

20x20 discrete ising; J = 1 3e+119

2.5e+119

2e+119 density of states

1.5e+119

1e+119

5e+118

0 -800

-600

-400

-200

0 E

200

400

600

800

Figure: Density of States

umbrella sampling of the 20x20 Ising model at t = 2; 10 bins, 50000 full steps per bin, edge choice 2 1e+06 900000 800000 700000 Punnormalised(E) 600000 500000 400000 300000 200000 100000 0 -800

-600

-400

-200

0 E

200

400

600

800

Figure: Unnormalized probability distribution

We can see that the peak of the canonical distribution PE () lies on the tail of the density of states (E ) so sampling (E ) gives us terrible statistics in our region of interest, PE () Solution: the Metropolis Algorithm (demonstrated in the next few pages) Samples the canonical probability distribution directly We build up an energy histogram (second gure) to get PE () ...from which, for an observable X , if we know X (E ), we can trivially get its expectation value X () .

Dynamics of the ising model two competing tendencies: energy minimization entropy maximization, which manifests as thermal uctuations . . . U - TS is minimized uctuations are essentially random, so the probability of occurrence of a microstate at each timestep only depends on the microstate that the system was in in the previous timestep (a Markov process) so all we need to carry out a simulation are the transition probabilities Wij = p(j i) that tell you the likelihood of going from state j to i, which when arranged in a matrix form the transition matrix

Consider an ensemble (set of systems) in a steady state steady state the probability of picking a microstate i at tn and at tn+1 is the same, Pi But the ensemble is not static; some systems (in that set) which are in the microstate i will go over to some other microstate j and vice versa For the probabilities to remain the same with respect to time, the probability ow into state i should be the same as that out of it

bl e

se

en

Figure: probability ow

en

se

bl e

Sampling from a canonical ensemble: the Metropolis algorithm

in other words, for an ensemble of size N, the number of systems falling into the microstate i between tn and tn+1 should be equal to the number of systems falling out of it. In our notation, Wij Nj
j(=i)

=
j(=i)

Wji Ni Wji (Pi N)

j(=i)

or ,
j(=i)

Wij (Pj N) =
j(=i)

Wij Pj

=
j(=i)

Wji Pi

Sampling from a canonical ensemble: the Metropolis algorithm

a condition called balance, which can also be written as (Wij Pj Wji Pi ) = 0

j(=i)

Balance is satised if every term in the above sum is equal to zero, a stronger condition called detailed balance: Wij Pi = Wji Pj and its said to be in equilibrium

Sampling from a canonical ensemble: the Metropolis algorithm

ij Were free to tune each ratio Wji separately, within the constraint of detailed balance

since we want to explore phase space as fast as possible, well make the Wij s as large as possible since the upper bound of a probability is one, well make the larger of the Wij s one For the canonical ensemble, letting the index i denote a state with energy Ei : Wij e Ei = E = e (Ei Ej ) j Wji e Assuming Ei > Ej , Wji > Wij so we let Wji = 1 and Wij = e (Ei Ej )

Sampling from a canonical ensemble: the Metropolis algorithm

the algorithm: start with a microstate chosen according to your convenience choose a spin at random calculate the change in energy if you ip it, E if its less than zero, ip it and go on if its more than zero, ip it with a probability e (E ) repeat the process from step 2.

Sampling from a canonical ensemble: the Metropolis algorithm

Shorthand for the algorithm: P(i j) = min e Ej ,1 e Ei

Sampling from a canonical ensemble: the Metropolis algorithm

Some sample C code:

i = ( i n t ) ( ( N) r ) ; //we u s e N b e c a u s e i n t e g e r c a s t i n g i s j = ( i n t ) ( ( N) r ) ; a s s e r t ( ( i <N)&&( j <N ) ) ; a o l d = a [ i ] [ j ] ; // o l d s p i n o l d = e n e r 4 ( i , j , N, a ) ; // o l d e n e r g y a [ i ] [ j ] = a [ i ] [ j ] ; new = e n e r 4 ( i , j , N, a ) ; // new e n e r g y d e l t a = new o l d ; / d e l t a can o n l y t a k e one o f t h r e e v a l u e s , s o t h e e x p o n i f ( ( d e l t a < 0 . 0 ) | | ( c h e c k [ d e l t a ] > r ) ) { // M e t r o p o l i s a l g o r i e += d e l t a ; // e n e r g y m += 2 a [ i ] [ j ] ; // m a g n e t i z a t i o n c o nt i n u e ; // new c o n f i g w o r k s } a [ i ] [ j ] = a o l d ; // new c o n f i g d o e s n t work

Sampling from a canonical ensemble: the Metropolis algorithm

Example: graph of < E > vs kT showing large variance near the critical point
energy vs temperature of 20x20 ising model -0.4

-0.6

-0.8

-1 <E>

-1.2

-1.4

-1.6

-1.8

-2 0 1 2 kT 3 4 5

Figure: < E > vs kT

Sampling from a canonical ensemble: the Metropolis algorithm

binder cumulants of the 4x4, 8x8, and 16x16 Vink model at p = 0.5; 500000 0.7

0.65

0.6 1 - <M4>/(3*<M2><M2>) 0.55 0.5 0.45 0.4 0.35 0.5 0.6 0.7 0.8 0.9 1 kT 1.1 1.2

Histogram reweighting
A problem of resolution (in a discrete set of temperatures) arises if we want to calculate a sharply peaked function of
2 1.8 1.6 1.4 specific heat Cv 1.2 1 0.8 0.6 0.4 0.2 0 1 1.2 1.4 1.6 1.8 2 kT 2.2 2.4 2.6 2.8 3 n=4 8 16 32

Figure: specic heats

Histogram reweighting

From what weve seen so far, one run of the metropolis algorithm at one temperature enables you to calculate quantities of interest only at that temperature. But the histogram PE () contains more more information than that - it encodes the density of states... PE (0 ) = (E ) e 0 E /Z (0 ) so (E ) = Z (0 ) PE (0 ) e 0 E

Histogram reweighting

Substitute for (E ): PE () = = = (E ) e E E E (E ) e Z (0 ) PE (0 ) e 0 E e E 0 E ) e E E (Z (0 ) PE (0 ) e PE (0 ) e (0 )E (0 )E E PE (0 ) e

Histogram reweighting
Note that the above expression for PE () is only accurate for energies E where PE (0 ) has been sampled suciently many times, that is, near 0 . PE (0 ) HE (0 ) so we can replace one with the other In practice, we calculate the reweighted histogram HE () = HE (0 ) e (0 )E and then normalize it to get PE () Using similar logic, if you have a distribution in M and E, PM,E () = PM,E (0 ) e (0 )E (0 )E ( E e M PM,E (0 ))

Histogram reweighting
1

0.8

0.6 n(E) 0.4 0.2 0 100

150

200

250

300

350 E

400

450

500

550

600

Figure: histogram generated at kT = 2.5 and reweighted for kt = 2.25 and 2.75

Umbrella Sampling

Near a critical point, the correlation length between spins diverges (note that this makes the histogram in energies broader, so it is a good thing if you want to reweight the histogram) domains of like (unidirectional) spin grow innitely large ipping a single spin always increases the energy so the system is now in a (local) potential well this is a problem for two reasons : time taken to reach equilibrium is large statistics are bad even after equilibration

Umbrella Sampling
Suppose were well below the critical point. In terms of probabilities, magnetizations MA and MB are very likely in the probability distribution at low tempetature:

Figure: probability vs magnetization

we dont care about MAB ; we want MA and MB to be suciently sampled.

Umbrella Sampling
If we were to choose states truly at random then everything would be ne, wed get good statistics but this is not the case for the single-spin-ip Metropolis algorithm, where successive states are adjacent to each other in magnetization (suppose we) start the simulation with the system in a state with magnetization MA Since we ip one spin at a time, there is no way that the system will be able to jump the gap between the peaks MA and MB ; the only way it can cross over is by passing through a state with magnetization MAB which, as you can see from he gure, is extremely unlikely

Umbrella Sampling

the crucial point: Sampling within the interval A with the correct proportion in relation to B requires many such crossings, which is unlikelier still. few coin tosses, large relative variance many coin tosses, small variance the system should spend equal time around MA and MB solution is to increase the chance of crossing over many ways of achieving this - umbrella sampling, replica exchange / parallel tempering, etc. well look at umbrella sampling

Umbrella Sampling
weighting and reweighting revisited: given an unnormalized histogram ni we calculate the mean of a quantity Xi thus: i ni Xi < X >= i ni multiplying and dividing with another distribution i , < X >= Letting ni = ni i , < X >= i ni (Xi /i ) i ni /i i (ni i ) (Xi /i ) i (ni i )/i

Umbrella Sampling
If we wanted to increase and the transition rate between the peaks corresponding to MA and MB , we could change the probability distribution pM in our canonical ensemble to another probability distribution pM M for which the probability of occurrence of MAB is more (possible choices of M in red below)

Figure: probability vs magnetization

Umbrella Sampling

What to do in practice: Since p(Mi ) is now p(Mi )(Mi ), dene (Mi ) = ln((Mi )) to get P(i j) = min pj pi , 1 = min e (Mj )Ej ,1 e (Mi )Ei

and the average of a quantity XM would now be < X >= M nM (XM /M ) M nM /M

Umbrella Sampling

What happens if we choose M to be a hat function: 1 within a certain range of M and zero everywhere else? We can no longer obtain < X > all we know is that, since ln(1) is zero, within the range of M where M = 1, P(i j) = min e 0Ej ,1 e 0Ei

the simulation is unweighted within that range; in other words all states are reweighted equally within that range.

If we partition the entire range of the order parameter M into windows (not necessarily uniform) with some overlap (the reason for which is explained later) and run the simulation for each window such that each time M is 1 within the window and zero everywhere else, well get a histogram in the magnetization which looks something like this:
umbrella sampling of the 20x20 Ising model at t = 2; 10 bins, 50000 full steps per bin, edge choice 2 15 14 13 12 ln(Pk(M)) 11 10 9 8 7 6 -400

-300

-200

-100

0 M

100

200

300

400

Figure: ln(nM ())

Umbrella Sampling
We know that the histogram in the magnetization is continuous for the model were looking at (the ising model) We could use that information to get the (unnormalized) probability PM () ... and then simply normalize it to get PM () how do we join up the edges? do we add, multiply or something else? we know that, whatever we do, the microstates shared by adjacent bins (since weve left some overlap) at adjacent edges must have the same value in the nal histogram. we also know that if we were to run the simulation for twice as long for some window, each entry would be updated twice as many times... so we multiply

Umbrella Sampling

instead of multiplying each window by a constant, its easier to take the log and then add a constant to each window so we just, in eect, shift the entries of the gure above up or down so that the edges match this can be done from left to right, for instance getting ln(PM ()):

Umbrella Sampling
umbrella sampling of the 20x20 Ising model at t = 2; 10 bins, 50000 full steps per bin, edge choice 2 16

14

12

ln(Pk(M))

10

2 -400

-300

-200

-100

0 M

100

200

300

400

Figure: ln(PM ())

Umbrella Sampling

Now we can nd the expectation values of any quantity that depends upon the magnetization alone But we cant reweight the histogram without knowing the corresponding energies of the microstates we need to record both M and E! use a 2-D histogram nM,E () as before, the histogram is partitioned by Magnetization as before, all histogram entries with M values within the same window are scaled by the same amount the probability distribuiton (PM,E ()) is shown on the next slide.

Umbrella Sampling

Figure: (PM,E ())

Umbrella Sampling

Umbrella Sampling
Successive Umbrella Sampling: For extremely steep probability distributions (eg. at extremely low temperatures) the statistics at one side of the window would be bad If we weight the probability distribution within the window by its inverse, well get an equal number of samples everywhere in the window, improving statistics We must somehow guess what the probability distribution is in the window: the rst window is simulated as usual the probability distribution in the second window is guessed at by quadratically extrapolating the probability distribution from the rst window the p.d. in the third window is guessed at by quadratically extrapolating the p.d. from the second window, and so on

Umbrella Sampling

To conclude: Pathological cases exist for which umbrella sampling could give erroneous results; to check, run the same simulation multiple times with a dierent choice of starting microstate within each window every time. If the (untted) histograms turn out to be signicantly dierent, it means that the system has been trapped in a local potential well.

Finding the Density of States: the Wang-Landau Algorithm

sometimes called the at histogram method for reasons to be shown nds the density of states directly, without recourse to the canonical ensemble was used to calculate the graph for the density of states shown earlier there are no issues of the system being trapped in a local potential well etc.

To see how it works, a) observe that if we were to choose a microstate at random from the density of states, we would pick a state with energy E with a probability proportional to (E ), that is, P1 (E ) = (E ) total

if we now keep the state that weve picked with a probability P2 (E ) = 1 (E )

then the likelihood that were left with a state of energy E is

P(E ) = P1 (E ) P2 (E ) 1 = total Weve got a uniform probability distribution (or at histogram) in energy.

Finding the Density of States: the Wang-Landau Algorithm

b) detailed balance for an equilibrium distribution Pi is Wij Pi = Wji Pj the Metropolis algorithm with the Pi = e Ei satises detailed balance with the steady state (or equilibrium distribution) being the canonical distribution in energy, P(E ) (E ) e E
1 Replace e Ei with (Ei ) and it satises detailed balance with the steady state being a uniform distribution in energy, 1 P(E ) = Pchoose (E ) Paccept (E ) (E ) (E )

Conversely, we know that some function g (E ) is proportional to the density of states if it gives us a uniform distribution P(E ) in energy when 1/g (E ) is substituted in the Metropolis algorithm . . .

Finding the Density of States: the Wang-Landau Algorithm

. . . that is P(E1 E2 ) = min g (E1 ) ,1 g (E2 ) (1)

for comparison, the equivalent expression for the canonical distribution is P(E1 E2 ) = min e E2 ,1 e E1

Well be using (1) as the transition probability for the three strategies below.

Finding the Density of States: the Wang-Landau Algorithm

strategy 0
keep g (E ) = 1 E , i.e. do a random walk in phase space, and record the histogram of energies H(E ), which approaches (E ) when the number of steps is large terrible, as was shown in the beginning, because the wings of (E ) (which is what were interested in) are poorly sampled.

Finding the Density of States: the Wang-Landau Algorithm

strategy 1
start with g (E ) = 1 for all E; record the energy histogram H(E ) for, say, 10000 steps; set g (E ) = H(E ) g (E ), H(E ) to zero, and go to the second step continue as long as possible now g (E ) approaches (E ) when the number of steps is large . . . Samples the wings better, but still too slow.

Finding the Density of States: the Wang-Landau Algorithm

Change from strategy 1: update g (E ) at every step instead of incrementing g (E ), multiply it by a factor f > 1 does not satisfy detailed balance! but tends towards it as f approaches 1.

Finding the Density of States: the Wang-Landau Algorithm

the Wang-Landau algorithm (third and nal)

1. start with g (E ) = 1 for all E, f = 2 2. ip a spin and accept it according to (1) 3. at the (new) value of E, set g (E ) = f g (E ) H(E ) = H(E ) + 1 4. repeat step 3 until H(E ) is at enough 5. set H(E ) to zero, f to f 6. go to step 2, and leave g (E ) alone! 7. continue as long as possible

Finding the Density of States: the Wang-Landau Algorithm

The Finer Points: actually, g (E ) becomes too large to handle so set ln(g (E )) = ln(f ) + ln(g (E )) the ratio of the density of states in (1) becomes exp[ln(g (E1 )) ln(g (E2 ))] at means Hmax /Hmin some number p which depends on f we can calculate the density of states as a function of any parameter, as long as that parameter has a nite value for every microstate . . . for example, if we want to calculate (M), just replace E with M in the entire algorithm above in the same way its straightforward to bin for more than one parameter: (M, E ) for example.

Finding the Density of States: the Wang-Landau Algorithm

20x20 discrete ising; J = 1 70000

60000

50000

E hist (flat)

40000

30000

20000

10000

0 -800

-600

-400

-200

0 E

200

400

600

800

Finding the Density of States: the Wang-Landau Algorithm

Primary justication: it works!
20x20 discrete ising; J = 1 1.6 from DoS T = 2.3 T = 2.4 T = 2.5 from Metropolis

1.4

1.2

1 <E>

0.8

0.6

0.4

0.2

0 0 0.5 1 1.5 2 2.5 kT 3 3.5 4 4.5 5

Finding the Density of States: the Wang-Landau Algorithm

A little more justication . . . its still a Markov process therefore, each individual step can be said to full detailed balance for the g (E ) in that step for a thorough investigation seeUnderstanding and improving the Wang-Landau algorith by Zhou, Chenggang and Bhatt, R. N., Phys. Rev. E 72

Improving the Wang-Landau Algorithm

It has been found 1 that, using the standard Wang-Landau algorithm, estimates of the density of states reach a limiting statistical accuracy which is not improved with further MC steps, because the large number of congurations generated towards the end of the simulation make only a small contribution to the calculated density of states It has been reformulated by Shell, Debenedetti and Panagiotopoulos 2 who record every proposed transition to nd the transition matrix, and thus the density of states.

Q. Yan and J. J. de Pablo, Phys. Rev. Lett. 90, 035701 (2003) S. M. Shell, P. G. Debenedetti, and A. Z. Panagiotopoulos, The Journal of Chemical Physics 119, 9406 (2003)
2

Improving the Wang-Landau Algorithm

The probability that the system goes from a state with energy I to a state with energy J in some monte carlo step is equal to the probability of starting o in a state I times the probability that, given it is in a state I, it jumps to a state J: (I ) T (I J) TOT This is the same as the probability that its in a microstate i with energy I, times the transition probability from state i to some state j with energy J, summed over i and j 1 TOT i j t(i j) s.t. i I, j J

Equating these two, we get T (I J) = (I )1 i j t(i j)iI jJ where iI is 1 when i I and 0 otherwise.

Improving the Wang-Landau Algorithm

So (J) i j t(i j)iI jJ T (I J) = T (J I ) (I ) i j t(j i)iI jJ t(i j) = (i j)Pacc (i j) where (i j) is the probability of proposing the transition i to j and and Pacc (i j) is the probability of accepting it. For a canonical ensemble, the probability Pacc (i j) is e Eij At t = , = 0 and Pacc (i j) becomes 1. For single-spin-ip dynamics, (i j) = (j i) so (J) T (I J) = T (J I ) (I )

But

Improving the Wang-Landau Algorithm

Knowledge of the transition matrix T will allow us to know the density of states But a random walk in microcanonical state space will only sample those energies I which are relatively likely Therefore: once we have built up a rough estimate of the DOS from the standard Wang-Landau algorithm, we can use it to do a random walk in (canonical) energy phase space and thereby better sample all non-zero transition matrix elements T (I J).

Improving the Wang-Landau Algorithm

If we store a matrix CIJ whose (I , J)th element is incremented each time there is a move proposed from a state with energy I to a state with energy J, then T (I J) CIJ within the row for any particular energy I, and adding the normalization condition J T (I J) = 1 we get T (I J) = CIJ /K CIK Because of this normalization, the expression given above for T (I J) is accurate regardless of the weight function (I ).

Improving the Wang-Landau Algorithm

Further improvement: In fact, sampling the TM elements during the time it takes to build up (E ) works quite well, in the sense that the DOS obtained from the TM elements is better than that obtained using the standard WL algorithm. So for faster convergence, we could periodically refresh the DOS obtained during the course of standard WL sampling with the the DOS obtained from the TM elements thus far recorded. For higher accuracy (and because detailed balance is ensured) we could start recording the nal version of T (I J) after weve reached a suciently accurate estimate of (E )

Conclusion

Topics I havent had time to cover: histogram reweighting for the entire range of temperatures: the multiple histogram method the Wol algorithm, and cluster-ip algorithms in general systems with continuous degrees of freedom (XY model, Heisenberg model) approaches to parellelization