MSD and self-diffusivity of water

[horlust]

Hello everyone!

I am trying to calculate the self-diffusitivy of water with MDAnalysis. I followed this tutorial: https://docs.mdanalysis.org/stable/documentation_pages/analysis/msd.html and the tips of Maginn (2018).

My code is:

import MDAnalysis as mda
import MDAnalysis.analysis.msd as msd

MSD = msd.EinsteinMSD(u, select='all', msd_type='xyz', fft=True)
MSD.run(verbose=True)

import matplotlib.pyplot as plt

msd = MSD.results.timeseries

nframes = MSD.n_frames
timestep = 0.001 # 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()
ax.plot(lagtimes, msd, ls="-", label=r'MSD')
plt.legend()
plt.show()

from scipy.stats import linregress

threshold = 1e-6

start_time = 0
start_index = int(start_time/timestep)

for i in range(len(msd) - 1):
if abs(msd[i + 1] - msd[i]) < threshold:
end_index = i
break

linear_model = linregress(lagtimes[start_index:end_index], msd[start_index:end_index])
slope = linear_model.slope
error = linear_model.stderr
dim_fac = 3 # dim_fac is 3 as we computed a 3D msd with 'xyz'

D = slope * 1 / (2 * dim_fac)
print("Diffusion Coefficient:", D)

My result with this code is kinda close to the experimental value: 2.12 nm^2/ns (I am assuming these are the units - could you confirm me that? I haven't seen it anywhere on the tutorial).

I performed my simulation with GROMACS, with a timestep of 0.001 ps. However, if I do u.trajectory.dt, I get 0.1 (I think the units are ps). If I use 0.1 as timestep in the code, I get 0.021 (again, I wanted to be sure of the units here).

So I don't know which calculation is correct, and I don't know exactly which are the units either.

Thank you in advance for your answer, and I wish you a merry christimas and a happy new year!
Best,
Maria

[lilyminium]

Hi @horlust, that's right, the time units in MDAnalysis should be in picoseconds. However, it is quite common to not save every frame due to the size of the trajectory you would end up with -- for example, might you have only saved every 100th frame to get a time difference of 0.1 ps?

Our length unit is typically angstrom so if there's no conversion happening my default expectations would be units of $Å^2/ps$. However, @hmacdope would know best here since he wrote the analysis!

[hmacdope]

@horlust, what @lilyminium says is all correct.

Units are $Å^2/ps$, however you must set the timestep axis to the real time between observations in your trajectory as @lilyminium noted.

Some back of the envelope maths:

$0.021 Å^2/ps = 0.21 nm^2/ns$

so it seems like we are still off by a factor of 10. What is the real time in between frames? If using GROMACS you can calculate it from your your .mdp file with some simple maths or use gmx dump or gmx trjconv to investigate the trajectory.

[horlust]

Thank you! My input on GROMACS was dt = 0.001, but with that I have an incorrect result. It seems that I have a result closer to the experimental value if I put a timestep of 0.01 on the code, but this is not the timestep that I specified on GROMACS nor the MDAnalysis output, so I'm confused. I'll keep looking! Maybe run a simulation containing more water molecules.

Btw, it might not be a MDAnalysis-related question, so feel free to not answer :-) All MSD graphics I've seen are somewhat linear. Mine is linear in the beginning and then it becomes constant (I've attached the image of the graphic). Is it normal?

Thank you very much for your answers, it was very clarifying. Happy new year!

[hmacdope]

@horlust there are two issues here.

1. by timestep between frames, we don't mean dt as specified in your mdp file. We mean the actual time between frames in your saved trajectory. For example, if nstxout-compressed=1000 writing an XTC file every 1000 steps, and your dt is 0.002 (2fs, very common) the time between frames in your trajectory is $1000*0.002 ps = 2 ps$ and your lagtime axis needs to be constructed as

lagtimes = 2*np.arange(n_frames)

2. The MSD profile you provided is typical of a trajectory where the particle positions are wrapped within the primary unit cell. Due to the way MSD is calculated, you need to have the positions present as if they have not been wrapped into the primary box. This can be acheived by using the NoJump() transformation in MDAnalysis https://docs.mdanalysis.org/stable/documentation_pages/transformations/nojump.html or preprocessing your trajectory into a new one with unwrapped positions using gmx trjconv ... -pbc nojump

Happy new year to you also!

[horlust]

Thank you for your answer! So the timestep is indeed 0.1.

I thought the trajectory was unwrapped because I used -pbc whole, -pbc nojump and finally -pbc mol with trjconv (I was using the -pbc mol .trr file), but I got a very good result using the transformation function you mentionned. :-)

[hmacdope]

Fantastic! Glad to help.