Understanding the MSD calculation in MDAnalysis

[Teslim Olayiwola]

Dear All,

I need some clarifications on the MSD code as contained in the v2.0 (available through the GitHub page). Using the approach as provided in refs [1] and [2]. I got the following plots. 

Figure 1 is the plot of MSD (in arbitrary units) against tau while figure 2 is the log-log plot of MSD against time. I am not sure about why my plot (as shown in Figure 1) didn't follow the same trend. Figure 2 is similar to that generated from the gromacs gmx-msd command. Kindly advise accordingly.

PS: my trajectory contains 99991 frames with timestep of 10. The modified is as shown below.

# trajectory information
traj = 'PAM_1000ns.xtc'
grofile = 'em.gro'

# load traj
u = mda.Universe(grofile, traj)

# calc MSD

atom_select = u.select_atoms('resid 25 and name C1')
MSD = msd.EinsteinMSD(atom_select, select='all', msd_type='xyz', fft=True)
MSD.run()

# access msd values

msd =  MSD.timeseries

# plot msd

import matplotlib.pyplot as plt
nframes = MSD.n_frames
print(nframes)
timestep = 10 # this needs to be the actual time between frames
lagtimes = np.arange(nframes)*timestep # make the lag-time axis
fig = plt.figure()
ax = plt.axes()
#plot the actual MSD
ax.plot(lagtimes, msd, ls="-",color='blue')
exact = lagtimes*6
# plot the exact result
ax.plot(lagtimes, exact, ls="--",color='black')

ax.legend(labels = ('3D random walk','$y=2 D\tau$'))
plt.show()

plt.loglog(lagtimes, msd)
plt.show()

[Hugo Macdermott-Opeskin]

Hey Teslim,

Can I see the Fig.1 but where you don't plot the "exact" result.  Just plot the MSD on its own.

As a side not there is no "exact" result for an arbitrary trajectory as only a random walk approaches the `D=2*d` limit, where D is the diffusion coefficient and `d` is the dimensionality of the MSD (3 in your case).

Cheers

Hugo

[Teslim Olayiwola]

See the fig in the attached image.

Your message "As a side not there is no ..." is not clear.

Thanks.

[Hugo Macdermott-Opeskin]

Hey Teslim,

A few points.

• In general the MSD should be a monotonically increasing function. Yours doesn't look that way to me. What is your trajectory of and what did you use for the `select` keyword?
• Are you using concatenated trajectories to measure the MSD? In general this is discouraged.
• Once you set the lagtime axis with `lagtimes = np.arange(nframes)*timestep` your x axis is no longer in arbitrary units. It should have units of time. What is the time between frames in your trajectory?
• My "side note" comment was intended only to explain that there is no "exact" result for most systems that we simulate using molecular dynamics. The "exact" result only holds either for an ensemble of brownian walks, or a single brownian walk in the long time limit.

Let me know if I can clarify anything for you :)

[Hugo Macdermott-Opeskin]

Also did you use a trajectory that has been "wrapped" into the primary unit cell? In general you should use a trajectory that has had no periodic boundary corrections applied.