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.
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'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 i 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 j 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 ) ).
Let me know by leaving a comment or connecting:
YouTube Channel: https://www.youtube.com/user/NicholasMSchneider