Issue with Mean Square Displacement

[Farshad Saberi]

Hello,

I have been using "MeanSquareDisplacement" from "MDAnalysis.analysis.waterdynamics" package. For some of my trajectories, I get the following error. What could be the problem? 

  File "/home/venv/lib/python3.6/site-packages/MDAnalysis/analysis/waterdynamics.py", line 1144, in _sameMolecTandDT

    b = set(selection[tf])

IndexError: list index out of range

I have attached the full error message as well as part of my script to run MSD analysis. The script is able to read topology (as a PSF file) and trajectories (as a DCD file) and create a universe object. I have attached a summary report of the number of atoms and trajectories read into the universe object.

I appreciate your help.

Best,

Farshad

[Farshad Saberi]

Hi,

I just wanted to give you an update about my issue. I noticed that there is a problem with my own trajectory. It must have had 2001 frames. But, it actually contains 1271 frames, as I verified in VMD. That's why MSD function stops at frame 1271 (which you can find in the error message).

However, one thing that misled me in the first place is the result of following commands. In particular, "u.trajectory.n_frames" command told me that my trajectory has 2001 frames, while VMD shows that it has 1271 frames. I don't know why the DCD.py module in MDAnalysis (or any other related module) couldn't capture this fact. My trajectory is a DCD file.

    with open(os.path.join(output_path, outfile), 'w+') as f:

        f.write("\n######################\n")

        f.write("Trajectory info")

        f.write("\n######################\n")

        f.write("Number of atoms in system: {}\n".format(u.trajectory.n_atoms))

        f.write("Number of frames in trajectory: {}\n".format(u.trajectory.n_frames))

        f.write("Trajectory timestep (The time difference between frames in picosecond): {}\n".format(u.trajectory.dt))

        f.write('Total time of trajectory: {0:8.3f} ps\n'.format(u.trajectory.totaltime))

[Oliver Beckstein]

Hi Farshad,

The DCD format stores the number of frames in the header of the trajectory. The header is written in the beginning so if the trajectory ended prematurely, the wrong n_frames will be reported. It’s a limitation of the DCD format.

I would rewrite your DCD with the correct n_frames. A simple way is

u = mda.Universe(PSF, oldDCD)

with mda.Writer(“new.dcd”, n_atoms=u.atoms.n_atoms) as W:

   for ts in u.trajectory:

       W.write(u.atoms)

This *should* create a trajectory with the correct header information. If not (I didn’t check the code) then one has to manually set “n_frames” but that might be a bit more involved (I don’t think u.trajectory.n_frames = XXXX will just work… but you can try).

Oliver

--
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 post to this group, send email to mdnalysis-...@googlegroups.com.
Visit this group at https://groups.google.com/group/mdnalysis-discussion.
For more options, visit https://groups.google.com/d/optout.

--

Oliver Beckstein

orbe...@gmail.com * orbe...@gmx.net

[Farshad Saberi]

Hi Oliver,

Thanks for your explanations. That makes sense to me. The problem was with our HPC server that failed to write properly one single frame into DCD trajectories. I've come up with a procedure to do appropriate sanity check on my DCD trajectories before running any MSDanalysis to prevent any similar issue in the future.

Thanks for your support.

Farshad