# Evaluating the MSD of my simulation

I am running a molecular dynamics simulation of water in TIP3P, and I am storing the positions of my particles in a 2D array called relevant_positions. The number of particles in my simulation is numPart. I am running a simulation from t=0 to t=n_time_points-1. In effect, I have snapshots of positions of particles at n_time_points time points.

I am trying to evaluate the MSD of my simulation, and this is the code I am using:

for d in range(1, n_time_points):
for i in range(0, n_time_points-d):
msd[d] += np.sum(np.square(\
relevant_positions[numPart*(d+i):numPart*(d+i+1),:] -\
relevant_positions[numPart*i:numPart*(i+1),:]))
msd[d] = msd[d]/(n_time_points-d)

msd = msd/numPart


The result I am getting with this is:

I expect this to be a straight line, but it clearly is not. What am I doing incorrectly here?

• As someone reading this with no prior context, can you define MSD? Is it the mean square diffusion or something along those lines? Aug 26 '20 at 21:20

I have no prior experience with molecular dynamics really, but to get you started I have found a resource with a pretty detailed method of calculating the MSD (Mean square displacement).

Looking at your plot, I suspect that periodic boundary conditions are somehow causing you a problem. It looks like it converges to a value that happens when everything scatters as far as it can on average from its starting position. When an atom wraps around the cell, its position will jump from a large number (such as 10.5 angstroms) to a small number (0.1 angstroms). This will be seen as displacement by your code, I believe.

If you want a really bad solution to fix your code though, you could throw away any displacements greater than a specific value corresponding to something like half the unit cell (an impossible distance to move in your MD simulation). This would remove the problematic data points.

Using ASE for distance calculations between atoms for each timestep would allow you to account for periodic boundary conditions as it supports taking distances with regards to it. This would be a much better solution.

• +1. Thanks for taking care of this one! We reached 64 unanswered questions, with only ~600 total questions, which means we dipped below 90% answered: area51.stackexchange.com/proposals/122958/matter-modeling. I appreciate your heroic efforts to keep our site in good shape! Aug 26 '20 at 21:48

Alright, so it turns out @TristanMaxson was right - it was the periodic boundary conditions that were messing with my computation. The solution to it was to find the unwrapped coordinates of my system.

This is how I unwrapped my coordinates:

unwrapped_positions = relevant_positions.copy()
for ts in range(0,n_time_points-1):
periodic_displacement = \
relevant_positions[numPart*(ts+1):numPart*(ts+2),:]-\
relevant_positions[numPart*(ts):numPart*(ts+1),:]

boundary_crossing = (np.abs(periodic_displacement) > (L/2))*1
boundary_crossing_sign = np.sign(periodic_displacement)

absolute_displacement = periodic_displacement\
- boundary_crossing*boundary_crossing_sign*L

unwrapped_positions[numPart*(ts+1):numPart*(ts+2),:] = \
unwrapped_positions[numPart*(ts):numPart*(ts+1),:] + \
absolute_displacement


followed by:

msd = np.zeros(np.shape(t))
msd[0] = 0

for d in range(1, n_time_points):
for i in range(0, n_time_points-d):
msd[d] += np.sum(np.square(\
unwrapped_positions[numPart*(d+i):numPart*(d+i+1),:] -\
unwrapped_positions[numPart*i:numPart*(i+1),:]))
msd[d] = msd[d]/(n_time_points-d)

msd = msd/numPart


As you can see, there was an inherent flaw in my code - I wasn't averaging my MSD right, apart from the fact that I had problems with wrapped coordinates.

• +1. Thanks for coming back and answering your own question, it will be helpful for future people if anyone else has the same problem that you had. I'm also glad to know that @Tristan was right :) Aug 27 '20 at 20:03