# Is the common formulation for (Ising Model) Monte Carlo simulations a bit off?

A common [Ising Model] Monte Carlo simulation repeats the following algorithm:

1. randomly pick a [flip] event
2. compute change in energy $$\Delta E$$
3. if new energy is lower than old energy, accept the change. Otherwise, accept the change with probability $$\exp(-\Delta E/kT)$$

In this formulation, all “downhill” energy moves are treated the same (occur with equal probability). In a physical system, wouldn’t larger “downhill” events occur more frequently? An easy reformulation that better matches with my expectations about a real system would be:

1. if new energy is lower than old energy, accept the change with probability $$1 - \exp(-|\Delta E|/kT)$$. Otherwise, accept the change with probability $$\exp(-|\Delta E|/kT)$$.

The thinking is that if an event has a probability $$p$$ of going to a higher energy state due to thermal energy, then it should have a probability of $$1-p$$ of staying in a higher energy state due to thermal energy.

• You might be interested in Glauber dynamics, which is a way of simulating the Ising model that's maybe more in line with your intuition. But the real answer is in the comments under Norbert Schuch's answer: the goal of Metropolis-Hastings and most related algorithms is not to simulate the physical dynamics at all, but to create an ersatz system that's known to have the same equilibrium distribution, so that we can sample from it. Traditionally we only care about the equilibrium distribution, so we want it to converge as fast as possible. Apr 4 at 23:27

The goal of any sampling algorithm has to be that it is unbiased -- i.e., the random walk it describes is in a configuration $$\vec s$$ with relative probability $$e^{-E(\vec s)/kT}$$.

This is achieved by the standard Metropolis sampling you describe.

It is also achieved by any other random process where the ratio of probabilities to go from $$\vec s$$ to $$\vec r$$ vs. back satisfies $$\frac{p_{\vec s\to \vec r}}{p_{\vec r\to\vec s}} = \frac{e^{-E(\vec r)/kT}}{e^{-E(\vec s)/kT}}\ , \tag{*}$$ the so-called "detailed balance" condition.

In particular, you could choose the transition probabilities equal to $$e^{-E(\vec r)/kT}$$ and $$e^{-E(\vec s)/kT}$$, respectively. This is not the formula you give (which will not work), but maybe more along the lines of your intuition.

Now why don't we choose the probabilities as above? Well, we want to maximize the probability for a move to be accepted (otherwise, we have to take far more samples), keeping the detailed balance condition (*). Thus, it is best to rescale the transition probabilities as much as possible -- that is, until the larger of them becomes equal to 1. This is precisely what the standard Metropolis update, as you describe it at the beginning of your question, does.

• Thank you for the explanation. There is still the conflict with my intuition that a larger drop in energy would occur more frequently than a smaller one. If my intuition holds, than it would only be appropriate to use the Metropolis update if all energy changes were of equal magnitude. Otherwise, only the largest possible energy drop ought to be assigned a probability of 1, and all others, including those that do reach a lower energy state but with less magnitude than a maximal change, should have an appropriately scaled probability. Can you explain why not? Apr 3 at 20:53
• In other words, can we say anything about $\frac{p_{\vec r\to \vec s}}{p_{\vec r\to\vec t}}$ if $\vec s$ is a lower energy state than $\vec t$? The treatment I outlined in the question says it's $1$ as long as both $\vec s$ and $\vec t$ are lower energy than $\vec r$. Apr 3 at 21:25
• @EricStimpson A larger drop in energy does happen (get sampled) more often than a smaller one. To make the discussion less confusing, let's name some states -- we're in state A now, B is a little lower energy, and C is much lower energy. If I go A -> B, then sometimes I go B -> C next step. But if I go A -> C, I almost never go C -> B next step. So I'm much more likely to be in C than B a little bit from now. Apr 4 at 14:09
• @EricStimpson: It's important to remember that the goal of a Monte Carlo "simulation" is not to accurately simulate the dynamics of a system; it's to generate a Boltzmann-distributed sample of states. Schroeder puts this nicely in his Introduction to Thermal Physics: "It is tempting to imagine that you are watching a simulation of what really happens in the magnet, as the dipoles change their alignments back and forth with the passage of time. ... Apr 4 at 15:11
• "But please remember that we have made no attempt to simulate the real time-dependent behavior of a magnet. Instead we have implemented a 'pseudodynamics,' which flips only one dipole at a time and otherwise ignores the true time-dependent dynamics of the system. The only realistic property of our pseudodynamics is that it generates states with probabilities proportional to their Boltzmann factors, just as the real dynamics of a magnet presumably does." Apr 4 at 15:11