MSD and self-diffusivity of water

[Maria Thereza]

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

# extract results
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()
# plot the actual MSD
ax.plot(lagtimes, msd, ls="-", label=r'MSD')
plt.legend()
plt.show()

from scipy.stats import linregress

# determine endpoints based on a threshold change in MSD values
threshold = 1e-6

start_time = 0
start_index = int(start_time/timestep)

# find the point where the change in MSD values falls below the threshold
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'

# calculate diffusion coefficient
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

[Rocco Meli]

Hi Maria,

We are retiring the mailing list and moving to GitHub discussions (blog post to come soon). The mailing lists are still active, but only to finish off ongoing discussions. Would you mind posting this new question on GitHub Discussions?

https://github.com/MDAnalysis/mdanalysis/discussions

Many thanks,

Rocco

--
You received this message because you are subscribed to the Google Groups "MDnalysis Discussion (ARCHIVED)" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mdnalysis-discus...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/964e9fdb-3513-4255-afd2-f43fc99ceda6n%40googlegroups.com.