Tuning the Ising model to arbitrary magnetization

This page contains a demo for our new grand canonical ensemble Monte Carlo tuning algorithm, using the two-dimensional Ising model as a test case. The Ising model is a simple theoretical model for ferromagnetism, the phenomena where magnetic dipoles of a material will tend to align with each other, forming a permanent magnet.

At each site \(i\) of a lattice, we have a binary variable \(\sigma_i = \pm 1\) called a spin. The Hamiltonian for the Ising model describing how energetically favorable different configurations \(\{\sigma_i\}\) of all of the spins are: $$H = -J \sum_{\langle ij \rangle} \sigma_i \sigma_j - B \sum_{i} \sigma_i$$ where \(J, B\) are both real parameters giving the strength of the interaction between spins and the external magnetic field, respectively. Statistical mechanics tells us that if the system is coupled to a thermal bath at temperature \(T\), the equilibrium distribution of the spin configurations is given by: $$P(\{\sigma_i\}) = \exp\Big(H(\{\sigma_i\}) / k_B T\Big) / Z $$ with \(k_B\) being Boltzmann's constant and \(Z = \sum_{\{\sigma_i\}} \exp\Big(H\{\sigma_i\})/k_B T\Big)\) being the partition function (essentially a normalization constant).

The goal is then as follows: Given a fixed value of \(J\), find the value of \(B = B^\ast\) such that the average magnetization $$ \langle m \rangle = \sum_{\{\sigma_i\}} \left(\sum_{i} \sigma_i\right) \exp\Big(-H(\{\sigma_i\}) / k_B T\Big) / Z $$ is equal to some target value \(m^\ast\). This is difficult, because we don't have an analytic solution for the Ising model probability distribution in nonzero field. An alternative/complement to analytic tools is Monte Carlo simulation, which produces information about the distribution by producing a sequence of samples. In physics, the dominant form of this approach is Markov-chain Monte Carlo (MCMC). An unfortunate drawback of MCMC is that sequential samples tend to be highly correlated with each other -- it takes a considerable amount of time for the Markov chain to "explore" the full distribution. If this is our only tool we can use to find \(B^\ast\), how can we do it?

The straightforward (but slow) type of approach would be to treat the simulation as a black-box, and use an off-the-shelf black-box optimization algorithm. For any value of \(B\), we run the simulation for some long enough period of time to estimate \(m(B)\), use the optimization algorithm to propose a new trial \(B\), and repeat until we have converged. However, this is (1) slow, (2) throws away information from previous simulations and (3) leaves precious information on the table. Taking the derivative of \(\langle m\rangle\) with respect to \(B\) gives a fundamental result of statistical mechanics that the gradient of an observable with respect to its coupling field (called the susceptibility, \(\chi\)) is exactly proportional to the fluctuations of that observable at the current field: $$ \chi \equiv \frac{\text{d}\langle m\rangle}{\text{d}B} = \frac{1}{k_B T} \left( \langle m^2 \rangle - \langle m \rangle^2 \right) $$ where the right hand side should be recognizable as a constant times the variance of \(m\).

With this in hand, we propose an algorithm to dynamically tune the field during the course of a single Monte Carlo simulation such that at sufficiently long times \(\langle m\rangle\) converges to \(m^\ast\). The idea behind the dynamic tuning algorithm is very simple:

1. Maintain running averages \(\overline{B}_t, \overline{m}_t, \overline{\chi}_t^{\mathrm{fluc}}\) of all observables, averaging over the previous half of the full simulation data. In particular, \(\overline{\chi}_t^{\mathrm{fluc}}\) is calculated from the variance of \(m\) over that data. The idea here is that we would like to use an increasingly large fraction of the simulation data to estimate all quantites, but samples from early-time dynamics should eventually "fall off" our estimates due to being nonrepresentative of the current parameters.
2. To handle poor early-time statistics and bad estimates of \(\chi\) during nonequilibrium dynamics, we place reasonable upper and lower bounds \(\overline{\chi}_t = \max[\chi_t^{\mathrm{min}}, \min(\chi_t^{\mathrm{max}}, \overline{\chi}_t^{\mathrm{fluc}})]\). We choose these bounds to be \(\chi_t^{\mathrm{min}} = \alpha / \sqrt{t+1}\) for some constant \(\alpha\) proportional to the system size, and \(\chi_t^{\mathrm{max}} = \mathrm{Std}[m]/\mathrm{Std}[B]\). For motivation behind these choices, see our paper.
3. For the next timestep, set $$B_{t+1} = \overline{B}_t + (m^\ast - \overline{m}_t) / \overline{\chi}_t .$$ Repeat.

For details on how this is actually implemented (in particular, how to achieve (1) in constant-time at each step), see the Rust source on Github, and for further results (including those on quantum Monte Carlo simulations) see our upcoming paper, Dynamical tuning of the chemical potential to achieve a target particle number in grand canonical Monte Carlo simulations, to be posted soon!

Below is a tool for playing around with our algorithm. The canvas shows the current spin configuration \(\{\sigma_i\}\), with red/blue pixels being \(\sigma_i = \pm 1\). Below the canvas are various sliders and menus to control both the system parameters and the tuning parameters. If the simulation is slow for you, either (1) pick a smaller system size, (2) hit the "Pause Plots" button to stop the plots from updating, or (3) refresh the page to clear the graph data.

Toggle Anim Randomize Tune!

50 100 250 500 Sys Size

Target \(m^\ast\) Running-average fraction \(\alpha\)