Simulations: Monte Carlo, molecular dynamics and the Ising magnet

Most interacting systems cannot be solved, so we let a computer try them

Intuition

A game of chess has a rule book you can write on one page, but no one can write down the winner of every possible game. Interacting particles are like that. The rules are simple (forces between neighbours) but the number of arrangements is far too big to solve on paper. So you play many games on a computer and watch what tends to happen.

The DNA zipper from the last lecture was solvable. Lecture 20 says most interacting systems are not. We cannot derive exact expressions relating, say, pressure to temperature. Computer simulation fills the gap. The lecture’s reasons:

• Many-body systems of interacting particles need computers.
• The fundamental laws of motion cannot be solved exactly except for a few special cases.
• Simulations give microscopic insight that experiments cannot observe, for macroscopic properties.
• They reach inaccessible regimes: extreme conditions, or materials that have not been made yet.
• They improve data analysis, for example X-ray and neutron diffraction and NMR.
• They give precise results to test a theory against.

Slide 5 draws the loop. A real system gives experimental results. The same system turned into a model gives simulation results and, with approximations, theoretical predictions. Everything is compared to improve the model and the theory.

Importance sampling visits likely configurations more often

The average of a quantity $X$ over the canonical ensemble is

$\langle X\rangle = \frac{\sum_{j=1}^{\tilde N_{\text{trials}}} X_j\, e^{-\beta E_j}}{\sum_{j=1}^{\tilde N_{\text{trials}}} e^{-\beta E_j}}$

There are two ways to compute it. Option 1: generate configurations with equal probability, then weigh each by $e^{-\beta E}$. The lecture marks this as the poor choice, because most random configurations have tiny weights and are wasted. Option 2: generate configurations with probability $e^{-\beta E}$ and weigh them equally. This is importance sampling. Think of polling a country: it is better to ask people in proportion to how many of them there are than to ask equal numbers from each street and fix it afterwards. We want the computer to produce configurations with probability proportional to their equilibrium probability.

Metropolis Monte Carlo accepts every downhill move and some uphill ones

Monte Carlo methods are probabilistic techniques that generate and sample configurations using random trial moves and accept or reject rules. The goal is to sample configurations and evaluate equilibrium properties, not dynamics. The name comes from the casinos of Monte Carlo.

The oldest and simplest version is Metropolis Monte Carlo. Repeat two steps.

1. Choose a random particle, change its state randomly, and calculate the new energy.
2. If the energy went down, accept the move. If it went up, accept the move with probability $\exp(-\beta\Delta E)$.

Why allow uphill moves at all? A ball that only rolled downhill would get stuck in the first dip. Allowing the occasional uphill step lets heat shake the system out, and the Boltzmann factor makes big uphill steps rare. Eventually this samples the canonical ensemble (constant $N$, $V$, $T$), and the configurations give equilibrium properties.

A simulation cell is small, so a large share of its particles sit at the surface and it does not look like bulk. The fix is periodic boundary conditions: surround the cell with virtual copies of itself, so a particle that leaves one side re-enters on the other. The cell must be large enough to avoid self-interactions and multiple interactions.

Exam detail

Learn the two Metropolis steps and the acceptance probability $\exp(-\beta\Delta E)$ for an uphill move. Monte Carlo aims at equilibrium sampling, not dynamics.

The Ising model is a lattice of up and down arrows

Intuition

A field of tiny compass needles that can only point up or down. Each needle prefers to match its neighbours. Heat shakes them randomly. At a low temperature the neighbours win and the needles line up, giving a magnet. At a high temperature the shaking wins and the needles point every which way, so the magnet vanishes.

Heat a magnet above its Curie temperature and it loses its magnetism. That is a transition between the ferromagnetic and paramagnetic phases. The Ising model is the simple model of it. Its assumptions:

• a lattice of spins,
• each spin can only point up or down, $\sigma = \pm 1$,
• interactions only with nearest neighbours.

The total energy is

$E = -B\sum_i \sigma_i - J\sum_{\langle ij\rangle}\sigma_i\sigma_j$

$B$ is the applied magnetic field. $J$ is the interaction between nearest-neighbour spins, and $\langle ij\rangle$ means a pair of neighbours. The minus signs mean that spins lined up with the field, or with each other (for positive $J$), lower the energy. The main quantity we want is the net magnetisation as a function of $T$:

$m = \sum_i \sigma_i$

$m$ is zero if spins point randomly up and down. $m$ is non-zero if one orientation is preferred. There is no exact solution in 3D, so the lecture first looks at computer simulation. Its simulation results are summarised as

$T \gg T_c: \ m \sim 0 \qquad\qquad T \ll T_c: \ |m| \sim 1$

The slide’s snapshots show a random mix at $T \gg T_c$ and almost one colour with a few isolated flips at $T \ll T_c$. Slide 13 lists questions to try on an interactive simulation (is the magnetisation constant at high $T$, why are nearby spins aligned at low $T$, and at what temperature does the magnetisation vary the most). The slides do not give answers.

Magnetic domains. Does every piece of iron below the Curie temperature show a net magnetisation? No. Spins align within domains, and if the domains point in random directions the net magnetisation is zero. Boundaries between domains are essentially grain boundaries. To magnetise iron, put it in a field to align the domains.

Molecular dynamics follows the motion, which Monte Carlo does not

Molecular dynamics (MD) simulates the dynamics of atoms and molecules by solving Newton’s laws numerically for interacting particles.

$\mathbf{F}_i = m_i\mathbf{a}_i = -\frac{\partial U}{\partial \mathbf{r}_i}$

The force on atom $i$ is minus the slope of the potential energy $U$ with respect to its position. Add up the pull from every neighbour, get the acceleration, move the atoms a little, repeat. The lecture’s example is argon-like attractive particles in a fixed-volume cell: starting as a gas, they condense into a liquid droplet as $T$ is lowered, like a cloud. MD is the most widely used method for biomolecule dynamics.

Compare

Monte Carlo

Random trial moves and accept or reject rules. Samples equilibrium configurations. No real time. Gives equilibrium properties.

Molecular dynamics

Solves Newton’s laws step by step. Follows the motion in time. Gives dynamics, such as condensation or biomolecule motion.

Pitfall

Do not say Monte Carlo simulates how a system moves in time. Slide 7 says its goal is equilibrium sampling, not dynamics. Time evolution is for molecular dynamics.

Recall

Why does Metropolis Monte Carlo sometimes accept a move that raises the energy?

So the system can climb out of local dips and sample the whole equilibrium distribution. An uphill move of size $\Delta E$ is accepted with probability $\exp(-\beta\Delta E)$, so big jumps are rare.