Nearest Image for Particle Simulations and Molecular Dynamics in MATLAB

NearestImage.png

Often when working on particle simulations and Molecular Dynamics you need to find the nearest image of a particle when working with Periodic Boundary Conditions. After working out several implementations in MATLAB, I have come to favor the approach I show here. I will a simple approach with a few assumptions that makes the implementation easier but it can easily be extended. This works great for diffusion or random walk simulations or for your molecular dynamics code and for on lattice or off lattice simulations.

If you don't have MATLAB at home, grab the software off Amazon.

 

 

Nearest Image In Periodic Boundary Conditions

In order to simulate an infinite domain, one can employ Periodic Boundary Conditions (PBCs) in their particle simulations. This, however, requires one to find the nearest image of a particle when computing inter-particle interactions. This is assuming one only needs short range interaction that are on a length scale smaller than the size of the domain. In many cases you can take just the closest particle images, and that is what we will do here. 

Assumptions

In order to make the implementation straight forward, we make the following assumptions:

  • Domain is of equal length in each spatial dimension (square or cube)

  • We only need the nearest image for any particle

  • Domain is from 0 to boxSize (not +-boxSize/2)

I also assume the following variables exist in your code to track particles:

  • particlePosition: MATRIX of SIZE [Spatial Dimension by Number of particles]

  • boxSize: SCALAR  (size of physical domain)

MATLAB Implementation

MATLAB's logical indexing makes this process quick and easy. Once I discovered this form of indexing, life became a lot easier; not to mention logical indexing makes for very readable code. If you haven't already, you should look at the Documented Indexing Options

We have to employ a loop to go through each particle, but the good news is we can use the logical indexing. We will not go through implementing the actual physics, but the code here will leave you with all of the deltas in each spatial direction (3D here, but could be 1D, 2D, or nD). It is important to note that each of the delta vectors, (dx, dy, dz) will have a length of the number of particles in you particlePosition array, and will enclude a 0 for the delta from particle to itself. 

for i = 1 : numParticles
    %Calculate Distance Components
    dx = particlePosition(1,i) - particlePosition(1,:) ;
    dy = particlePosition(2,i) - particlePosition(2,:) ;
    dz = particlePosition(3,i) - particlePosition(3,:) ;
   
    
    %Find nearest image in PBC by reassinging only when image is closer
    % j to the left
    dx(dx > boxSize/2) = dx(dx > boxSize/2) - boxSize ;
    dy(dy > boxSize/2) = dy(dy > boxSize/2) - boxSize ;
    dz(dz > boxSize/2) = dz(dz > boxSize/2) - boxSize ;
    % j to the right
    dx(dx < -boxSize/2) = dx(dx < -boxSize/2) + boxSize ;
    dy(dy < -boxSize/2) = dy(dy < -boxSize/2) + boxSize ;
    dz(dz < -boxSize/2) = dz(dz < -boxSize/2) + boxSize ;
    
    %HERE IS WHERE YOU STORE OR USE
    %Could at least calculate all of the square distances
    %This will be covered in another post

end

In the loop we calculate the deltas for all of the particles in the simulated domain. Once we have these vectors, we can look for the cases where the image to either the left or right (negative or positive directions) is the shorter distance. For each value where this is true, we shift the delta by either plus or minus one half the domain with, giving us the image particle.

A 1D example illustrates the idea here. Say particle i is at position 9 on a domain of length 10, and particle is at position 2. The delta from 9 to 2 is 7, but the closest version of particle j is actually its image at the imaginary position 12. Since 7 > 10/2, we shift the delta by subtracting the domain length, or 7 - 10 = -3 = 9 - 12. The negative here is just showing the image particle is to the right of particle i. This works in both directions (i.e., j to i: 2 - 9 = -7 => -7 + 10 = 3 = 2 - ( -1 ) ).

Questions?

Let me know by leaving a comment or connecting:
YouTube Channel: https://www.youtube.com/user/NicholasMSchneider
Twitter: https://twitter.com/NMSchneider 
Facebook: http://facebook.com/Schneider.NicholasM
LinkedIn: http://www.linkedin.com/in/nmschneider/
Google+: https://plus.google.com/+NicholasSchneiderPennRIT/



1 Comment

Nicholas M Schneider

Nicholas M Schneider is a 2010 graduate from the Kate Gleason College of Engineering who is now a Doctoral Candidate at the University of Pennsylvania. Originally from an obscure town south of Buffalo, New York, he attended the Rochester Institute of Technology where he received concurrent Bachelor of Science and Master of Science degrees. While there he had a number of Co-ops including a six month stay as a Design Engineer at Lockheed Martin and Research positions with Dr. Satish Kandlikar. Nicholas currently works with Dr. Haim H Bau in the field dubbed “in situ electron microscopy of liquid systems” where he studies applications in energy and biological systems. Outside of the lab, Nicholas Schneider is a Graduate Associate in Rodin College House and enjoys running (he ran his second Philly Marathon this past November), cooking, baking, reading, and justifying his coffee addiction by making it a hobby.