Simulate diffusion on a grid using Python

by Christoph Gohlke

Updated on June, 28, 2022

This notebook is released under the BSD 3-Clause license.


This Jupyter Notebook uses the Python 3 programming language, the numpy library for n-dimensional array-programming and linear algebra, the numba Python compiler, and the matplotlib library for plotting.

The mean square displacement (MSD) of Brownian particles in n-dimensions Euclidean space is:

$$ MSD = \left<\left(x_1(t) - x_1(0)\right)^2 + \left(x_2(t) - x_2(0)\right)^2 + ... + \left(x_n(t) - x_n(0)\right)^2\right> = 2nDt $$

where $t$ is the time, $D$ the diffusion coefficient, $n$ the dimension, $x_n(t)$ the cartesian coordinate of a particle in dimension $n$ at time $t$, and $\left<\right>$ the ensemble average over all particles.

Some of the functions used in the first part of the simulation code:

Diffusion on a one-dimensional grid

Plot the positions of some particles over time

Plot all particle positions as a color-coded image

Plot a histogram of particle positions at the end of the simulation

Plot MSD and line fit

Remove all names defined above from the global namespace. Global names might interfere with following code.

Import libraries and modules used in this document:

Diffusion on a n-dimensional grid

To make the code more modular, manageable, extensible, and reusable, it is refactored into small functions for simulation, particle counting, analysis, and plotting. The simulation model is extended to handle many dimensions. Hook functions allow customizing the initialization and restriction of particle movements. A simple detector "particle counter box" and methods to analyze and plot particle counts are added.

Diffusion in one, two, and three dimensions are simulated and compared.

Constrained diffusion in a box

To demonstrate diffusion under restrictions, particles are placed in a box. Particles can diffuse unconstrained within the box. Three different cases are explored when particles hit a box boundary:

Hook functions are defined for each case and passed to the simulation function.