MSD calculation

[Maria Thereza]

Hello everyone!

I am trying to calculate MSD values for two trajectories. In both trajectories, there is a single molecule of sugar in water, the sugars being different from each other. 

I don't know if I'm doing something wrong, but the MSD curve I'm getting is exactly the same for both molecules.

My code is:

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

msd1 = msd.EinsteinMSD(u1, select='resname UNL', msd_type='xyz', fft=True)
msd1.run()

msd2 = msd.EinsteinMSD(u2, select='resname UNL', msd_type='xyz', fft=True)
msd2.run()

msd_1 = msd1.results.timeseries
msd_2 = msd2.results.timeseries

nframes = msd1.n_frames
timestep = 1 # 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_1, label='1')
ax.plot(lagtimes, msd_2, label='1')
exact = lagtimes*6
# plot the exact result
# ax.plot(lagtimes, exact, label=r'$y=2 D\tau$')
plt.legend()
plt.show()

Is there an explanation for having exactly the same results for both systems? Am I doing something wrong?

Thank you very much!

Best,

Maria

[Hugo Macdermott-Opeskin]

Hi Maria, could you also show where you create u1 and u2? 

--
You received this message because you are subscribed to the Google Groups "MDnalysis discussion" 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/10a4efa2-b67e-4d64-bce0-39411309e35dn%40googlegroups.com.

[Maria Thereza]

Hello Hugo!

Thank you for your reply :)

I do something like:

# path to sorbitol files
%cd /home/maria/Documents/simulations/sorbitol_17102023/output_files/md

# define sorbitol universe with MDAnalysis
u1 = mda.Universe('md.tpr', 'md_mol.trr')

# path to mannitol files
%cd /home/maria/Documents/simulations/mannitol_22102023/output_files/md

# define mannitol universe with MDAnalysis
u2 = mda.Universe('md.tpr', 'md_mol.trr')

It works for the rest of the code (I calculated other properties, like radius of gyration, RMSD values...).

[Hugo Macdermott-Opeskin]

Hi Maria,

would you be able to check that the arrays are binary identical rather than just look similar on a plot. 

Try

```
import numpy as np 

np.assert_allclose(MSD1.timeseries, MSD2.timeseries)

```

Also double check the simulation itself, what is the molecule in question doing? Do the simulations look different.

Something to note that is not technically incorrect but more something to think about is that an MSD is only really well defined when a representative ensemble average is taken. 
This is over both time and molecules, however as I understand you only have one molecule. This is unlikely to be a good estimate for the MSD, except in the long time limit if that makes sense? This paper is an excellent reference: https://livecomsjournal.org/index.php/livecoms/article/view/v1i1e6324

Cheers, 

Hugo