Random walk Monte Carlo seems a bit silly at first glance, since a simulation does not really contain a lot of information. As the name implies, the motion of the particle is completely random, much like Brownian motion of a particle suspended in a fluid. The behavior of a particle in a random walk simulation is very much like the movement of the average particle in a fluid, which randomly and constantly bounces into other particles, every time from a random direction.

A collection on particles will then have average properties much like a liquid or (as we have looked a lot at in the course) the Lennar-Jones fluid. Do note, that there are no interactions between any particles.

A simple one-particle random walk Monte Carlo simulation is shown here.

The equations of motion are shown here:

$$x(t + dt) = x(t) \pm \mathrm{random.random()} \cdot \mathrm{d}l $$

$$y(t + dt) = y(t) \pm \mathrm{random.random()} \cdot \mathrm{d}l $$

where random.random() is the Python default random number generator and $\mathrm{d}l$ is a scaling factor that determines the average step-length. The fact that the equation does not include any interaction terms makes it computationally cheap to simulate a lot of particles.

In the last Molecular Statistics class we used random walk Monte Carlo to obtain the diffusion coefficient via Fick's second law, here in one dimension:

$$\frac{\partial \phi}{\partial t} = D\ \frac{\partial^2 \phi}{\partial x^2}$$

The solution of Fick's second law is a Gaussian function

$$\phi (x, t) = \frac{1}{\sqrt{2\ \pi\ \sigma^2}} \exp \left( - \frac{x^2}{2 \sigma^2}\right)$$

with $\sigma = \sqrt{2\ D\ t}$. So if we can calculate $\sigma$ from the distribution of particles, we can immediately get the diffusion coefficient.

The simple approach is to start a random walk simulation with many particles and then fit the distribution of particles to a Gaussian.

We used this extremely simple Python/Numpy program to simulate random walk and calculate $D$:

**X = numpy.zeros(n_particles)**

**for i in range(n_steps):**

**X += numpy.random.uniform(-dl, dl, n_particles)**

**sigma = X.std()**

**print "D =", sigma**2 / ( 2 * i)**

The simulation looks like this in 2D (with 400 particles):

As you can see, the particles collectively look like they are diffusing through a medium. The diffusion constant is the derived (here only in the $X$-direction by fitting the distribution to a Gaussian function (which was simpy done via Numpy's .std() function). The time evolution of the distribution of particles (blue histogram) and the Gaussian fit (red curve) looks like this:

The time evolution of $D$ should converge pretty quickly. You can see $D$ for a simulation with 10,000 particles below. It took a matter of seconds on Jimmy's old laptop to run this. Imagine doing an long enough MD with 10,000 particles on your laptop!