Evaluating the MSD of my simulation

Question

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?

Answer 1

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.

Answer 2

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.