If you are a student and don't have MATLAB at home, grab the software off Amazon for $99.
Earlier this year I wrote up a simulation for on-lattice diffusion in MATLAB. Since then, several people have asked me about my implementation. It is very hard to pick up someone else's code, so in this post I will walk through the details of the code line-by-line with some examples of how each section works. I will focus on 2D, but the code could easily be extended to 3D.
Often for research we find the need to model systems of diffusing particles. The application might be aggregating nanoparticles, electrodeposition, or even simulations of molecular motors. These simulations are striking for their simplicity and ability produce qualitatively accurate results. All of this is done with a few simple (algorithmic) rules called a "logical geometric description" of the world.
Allowing particles to occupy the same space in a simulation makes implementation significantly easier. This deviation from reality is bothersome in problems where diffusion is important, especially systems like electrodeposition. Here we will walk through the implementation for an on lattice system of particles in two dimensions.
We will cover
- Setting Up the Initial Domain
- Random Distribution
- Algorithmic Distribution
- Arbitrary Imported Distribution
- Random Walk On-Lattice
- Checking for Particle Overlap
- Creating Image Sequences
Random Particle Distribution On-Lattice
First we will look at a uniform distribution of particles. MATLAB makes this very easy for us with the rand( ) function. A call of rand(M,N,P,...) will generate an N x M x P x ... array of random numbers from an uniform distribution on the interval of [0,1], inclusive. Go ahead and make a call in the Command Window and assign it to a variable called domain. As an example, we will work with a small matrix so it is easy to see what is going on.
>> domain = rand(3,2) domain = 0.4039 0.9421 0.0965 0.9561 0.1320 0.5752
This is great because we can exploit this to get any density of randomly distributed particles simply by making a logical comparison. Say we wanted an area density of particles that is, on average, 1/3. We can use the less than (or greater than) logical operator on this matrix to return a new matrix of ones and zeros, where we will make the ones our particles. We want we want to know when domain is less than 1/3:
>> domain < 1/3 ans = 0 0 1 0 1 0
Notice that this returns an Array the same size as the original that is also of the logical type. We could have assigned this result right back to domain, but we can actually do better. We might want to simulate systems with more than one species, and so we may need to track the type of particle. This would most easily be done by storing particles in the array as different numbers. That is, 0 is empty space, 1 is the first species, 2 is the second, and so on.
To do this we will instead first define how large we want our domain to be, and then decides how to randomly place our particles. Now would be a good time to open a new script file and save it with some name. I will call mine NMS_DiffusionWalkthrough.m. If we wanted only one species with a density of 0.03, we would use the following:
% Random Initial Distribution with One Species domainSize = [400 400]; domain = zeros(domainSize); domain(rand(domainSize) < .03) = 1;
If you now open the domain as a figure (figure; imshow(domain)), you will see a nice distribution of particles.
Algorithmic Particle Distribution
We may want to define a region in our domain as particles. Here we will first create a region with several, evenly spaced circles, and then create a square region of particles. Let us first go ahead and define a rectangular domain of zeros.
domainSize = [400 800]; domain = zeros(domainSize);
First, let us define a simple square on the left half of the domain by:
domain(150:250,150:250) = 1;
Now, we will exploit a few built in functions to make it easier to place an arbitrary number of circles in our domain. We will use the morphological operator imdilate to replace points in our domain with a mask of our design. In this case we define an array of points as the circle centers, and then we generate a circle with a 20 pixel radius using the built in function fspecial, a function that is usually used to create image filters.
domain(50:75:350,450:75:750) = 1; % Define centers circle = fspecial('disk', 20) ~= 0; % Create a mask of 1's domain(:,401:end) = imdilate( domain(:,401:end) , circle ); % Replace points with mask
Where we only put the circles on the right side of the domain. When you run this, you should get a domain that looks like this:
Had we run the dilation over the entire domain, we would have a rounded square that is 20 pixels larger.
Arbitrary Imported Distribution
Lastly we will import and arbitrary image to use as our starting domain. We just need to use the import function and make sure we shape our array to fit the form we have set up. Meaning we need to make sure it is an N x M x 1. We accomplish this by taking only a single channel from an RGB image (PNG in this case).
domain = imread('./sample_input.png'); domain = domain(:,:,1); domainSize = size(domain);
On-Lattice Random Walk
We now have a spatial representation of the particles within our domain. Each particle is represented by a white pixel (value of 1) when the omain is displayed as an image. We will now define an array that tracks the location of each particle given its ( x , y ) location. (Remember that the first dimension in MATLAB is down the array and the second is toward the right.) We can accomplish this task in a single line via the find function.
% Generate Position Arrays [particlePosition(:,1), particlePosition(:,2)] =find(domain);
The call of find(), when assigned to two variables as above, returns the row and column of each non-zero entity in the input array. In our case the ones from our particles will each return their location.
Now that we know where each particle is, we can randomly move them around. Since we want to stick with an on-lattice arrangement, the particles will only ever move in one of four directions (+/- x or +/- y). This allows the use a switch structure to make our movement choice for each particle. We loop through each of the particles and move them in one of the defined directions at random:
for i = 1:length(particlePosition)
%Select Direction to Move
switch ceil(4 * rand)
dR = [-1 0];
dR = [+1 0];
dR = [0 -1];
dR = [0 +1];
%New Particle Location
tempPosition = particlePosition + dR;
particlePosition(i,:) = tempPosition;
Where ceil(4 * rand) returns an integer from the set (1, 2, 3, 4) to be used as the selector for the step. If you were wrap this in another for loop to iterate the stepping process, you would notice that the particles quickly leave the originally defined domain and can overlap one another. To solve the first problem, we will implement periodic boundary conditions and we use a temporary variable called tempPosition to hold the place we are moving so we can solve the latter problem.
Periodic Boundary Conditions (Screen Wrapping like Astroids Game)
When we want to simulate an infinite medium, we use periodic boundary conditions. This will assume a mirror image of our domain in every direction/across each of the four boundaries. Since a particle moving can leave the screen to the left, right, top, or bottom, we need to check for each case separately. We can easily do this by adding an if/elseif set of conditionals before we move our particle. Go ahead and add the following code before the line where we move the particle (still in the for loop).
%Periodic Boundary Conditions
if tempPosition(2) < 1 %Left
tempPosition(2) = tempPosition(2) + domainSize(2);
elseif tempPosition(2) > domainSize(2) %Right
tempPosition(2) = tempPosition(2) - domainSize(2);
elseif tempPosition(1) < 1 %Top
tempPosition(1) = tempPosition(1) + domainSize(1);
elseif tempPosition(1) > domainSize(1) %Bottom
tempPosition(1) = tempPosition(1) - domainSize(1);
Where all we do is shift the particle by the coordinating domainSize for the direction we left. Remember that MATLAB starts at an index of one, and so this is why we check for values less than one and greater than the domainSize.
Check for Particle Overlap (Steric Effects)
Lastly, before we move, we want to check to make sure we do not overlap another particle. This provides some steric effects to our simulation. We could loop through all the other particles and check to make sure our move is open, however, MATLAB's way of indexing positions in an array will give us a faster means of performing the check.
Any position in array can be indexed in two different ways. The first is by its coordinates as in array(x1 , x2, x3, ...), also called the subscript notation. Secondly, MATLAB gives each value an increasing integer index by moving through the dimensions (index notation). This means in a 3 x 3 array, each position can be index by its row and column or by an integer from 1 to 9. For example array(2, 1) = array(2) and array(1,2) = array(4).
We shall also track our particle using their indexes within the larger domain. This will help us speed up both check if a location is open to move a particle and in creating an image for a given time step. We use the sub2ind function convert from subscript notation to index notations. Go ahead and add a line after you assign the locations to the particlePosition array and before the start of the for loop to convert the indices:
particleIndex = sub2ind( domainSize, particlePosition(:,1), particlePosition(:,2));
Note that for the first time, you also could have just gone:
particleIndex = find(domain);
Both will return the same results, and find will actually be faster. Going back into our loop, however, we will use sub2ind to convert our new position to an index so we can use a simple logical operator to make sure we will not overlap with an existing particle.
If we make a vector of our new index that is the same length as our particleIndex vector, then we can directly compare them. If we sum along the resulting logical vector and and get a 0, then we know we can move the particle to the new location. If the sum along the logical vector is not 0, then the new location is occupied. Note that this conversion and logical comparison should be much faster than looping through all of the neighbors and seeing if they are filled. Go ahead, comment out the code that you used to move the particle before and add the following check below the periodic boundary conditions and the loop end statement:
% Move Particle
% particlePosition(i,:) = tempPosition;
% Convert new location to indicies for speedy search
tempIndex = sub2ind( domainSize, tempPosition(1), tempPosition(2));
% If where we want to move is empty, move
if ~sum(ones(length(particleIndex),1)*tempIndex == particleIndex)
particlePosition(i,:) = tempPosition;
particleIndex(i) = tempIndex;
Loop for Many Steps
Now you are ready to wrap the whole thing in another for loop and iterate. Start with say 10 steps to see what it looks like. If you use the Algorithmic Distribution as the starting point, you should get a domain that gets something like the following. (This might take a bit as there almost 44k particles)
For fun, if we allow overlap of particles, we see that the structure of our starting configuration deteriorates quickly. Again, after 10 steps:
Creating Image Sequence
The last thing we will want to do is create an image sequence so we can see how the system develops in time. Before your outer for loop, you will want to create a subfolder to store the images in since we will create many of them. Add this code just before your step loop:
% Write image file mkdir('./results') imwrite(domain,['./results/_',num2str(0),'.png'])
This will create a subfolder and write out the initial configuration. Next, after the loop through each particle and before the end of the stepping loop (so between the two end statements), we will add code to write images at a defined interval.
% Write Image File at some rate (set to 1 or every step right now)
if mod(step,1) == 0
% Update Domain Image
domain = zeros(domainSize);
domain(particleIndex) = 1;
% Print to Screen at some other rate
if mod(step,1) == 0
fprintf(1,'Time is %d \n', step )
Here we check if the step just completed is evenly divisible by our rate, and if so, we write the image and print out to the workspace the current step. Right now we have set these to 1, but for long simulations it might make sense to have different values for when these actions are performed. These could be combined inside a single if statement, but have been separated out to allow for different rates for each action.
Creating an AVI
Now that you have the image sequence, you can create an AVI if you would like. One could use MATLAB to create the movie as your program iterates by adding the following code after the imwrite call (or you can replace/comment out the imwrite).
aviobj = addframe(aviobj,domain);
You also need to initialize the object before you loop with:
aviobj = avifile('example.avi','compression','None');
And close the object at the very end with:
aviobj = close(aviobj);
Personally, I prefer to create the image sequence and then use ImageJ to create an AVI when needed since I can change the video properties easily, and you can always import any frame back into the program to start again.
Comparing No Overlap Diffusion to Diffusion with Particle Overlap
Check out the results from what we just created!