MSDs are incorrect for wrapped coordinates..

[hmacdope]

Expected behavior

Means squared displacements should respect PBCs, but we currently do not have the right unwrapping transform to do this for "wrapped" trajectories, ie trajectories where all the coordinates are in the primary box such as those produced by GROMACS.

This makes us far too reliant on correct user input, with inadequate docs surrounding the potential sources of errors.

We need to improve warnings and docs and/or possibly throw errors. We should also look at actually implementing the correct transform. (pbc nojump in GROMACS).

Actual behavior

Our current implementation is technically correct, but relies on the coordinates being correctly set out in the box.

Code to reproduce the behavior

This is best demonstrated on a small box with fast diffusion, such as a small box of water. The MSD will plateau rapidly and for a cubic box should be near l * 3**(1/2)) where l is the box length.

import MDAnalysis as mda
import MDAnalysis.analysis.msd as msd
import matplotlib.pyplot as plt
import numpy as np


u = mda.Universe('waterbox.tpr', 'waterbox.trr')
MSD = msd.EinsteinMSD(u, select='all', msd_type='xyz', fft=True)
MSD.run()

nframes = MSD.n_frames
timestep = 1 # this needs to be the actual time between frames
lagtimes = np.arange(nframes)*timestep # make the lag-time axis

# find the min and max of the positions
# if the trajectory is wrapped, these will all be inside the primary box
print(np.max(MSD._position_array, axis=1))
print(np.min(MSD._position_array, axis=1))


#Plot the actual MSD
plt.xlabel('Time (ps)')
plt.ylabel('MSD')
plt.title('MSD')
plt.plot(lagtimes,MSD.timeseries)
plt.show()

....

Current version of MDAnalysis

• Which version are you using? (run python -c "import MDAnalysis as mda; print(mda.__version__)")'
• 2.0.0-dev
• Which version of Python (python -V)?
• 3.7
• Which operating system?
• Ubuntu