Unable to compute RDF for .gsd files

[Lynn3512138]

Hi All,

I'm using MDAnalysis 2.10.0 and try to get RDFs from a .gsd trajectories. The code is as below:

u=mda.Universe("2_B50_try7.gsd","2_B50_try7.gsd")

c=u.select_atoms("name C")
a=u.select_atoms("name A")

print(u.trajectory)
print(len(c))
print(len(a))

rdf = InterRDF(c,a,range=(0,5))
rdf.run()

plt.plot(rdf.bins, rdf.rdf)

However it gets the problems:

<GSDReader 2_B50_try7.gsd with 50 frames of 100 atoms>
50
50

TypeError Traceback (most recent call last)
Cell In[45], line 11
8 print(len(c))
10 rdf = InterRDF(c,a,range=(0,5))
---> 11 rdf.run()
13 plt.plot(rdf.bins, rdf.rdf)

File /opt/homebrew/Caskroom/miniconda/base/envs/hoomdenv/lib/python3.13/site-packages/MDAnalysis/analysis/base.py:898, in AnalysisBase.run(self, start, stop, step, frames, verbose, n_workers, n_parts, backend, unsupported_backend, progressbar_kwargs)
891 computation_groups = self._setup_computation_groups(
892 start=start, stop=stop, step=step, frames=frames, n_parts=n_parts
893 )
895 # get all results from workers in other processes.
896 # we need AnalysisBase classes
897 # since they hold frames, times and results attributes
--> 898 remote_objects: list["AnalysisBase"] = executor.apply(
899 worker_func, computation_groups
900 )
901 self.frames = np.hstack([obj.frames for obj in remote_objects])
902 self.times = np.hstack([obj.times for obj in remote_objects])

File /opt/homebrew/Caskroom/miniconda/base/envs/hoomdenv/lib/python3.13/site-packages/MDAnalysis/analysis/backends.py:208, in BackendSerial.apply(self, func, computations)
192 def apply(self, func: Callable, computations: list) -> list:
193 """
194 Serially applies func to each task object in computations.
195
(...)
206 list of results of the function
207 """
--> 208 return [func(task) for task in computations]

File /opt/homebrew/Caskroom/miniconda/base/envs/hoomdenv/lib/python3.13/site-packages/MDAnalysis/analysis/base.py:554, in AnalysisBase._compute(self, indexed_frames, verbose, progressbar_kwargs)
551 if len(frames) == 0: # if frames were empty in run or stop=0
552 return self
--> 554 for idx, ts in enumerate(
555 ProgressBar(
556 self._sliced_trajectory, verbose=verbose, **progressbar_kwargs
557 )
558 ):
559 self._frame_index = idx # accessed later by subclasses
560 self._ts = ts

File /opt/homebrew/Caskroom/miniconda/base/envs/hoomdenv/lib/python3.13/site-packages/tqdm/notebook.py:250, in tqdm_notebook.iter(self)
248 try:
249 it = super().iter()
--> 250 for obj in it:
251 # return super(tqdm...) will not catch exception
252 yield obj
253 # NB: except ... [ as ...] breaks IPython async KeyboardInterrupt

File /opt/homebrew/Caskroom/miniconda/base/envs/hoomdenv/lib/python3.13/site-packages/tqdm/std.py:1169, in tqdm.iter(self)
1166 # If the bar is disabled, then just walk the iterable
1167 # (note: keep this check outside the loop for performance)
1168 if self.disable:
-> 1169 for obj in iterable:
1170 yield obj
1171 return

File /opt/homebrew/Caskroom/miniconda/base/envs/hoomdenv/lib/python3.13/site-packages/MDAnalysis/coordinates/base.py:333, in FrameIteratorIndices.iter(self)
331 def iter(self):
332 for frame in self.frames:
--> 333 yield self.trajectory._read_frame_with_aux(frame)
334 self.trajectory.rewind()

File /opt/homebrew/Caskroom/miniconda/base/envs/hoomdenv/lib/python3.13/site-packages/MDAnalysis/coordinates/base.py:870, in ProtoReader._read_frame_with_aux(self, frame)
868 def _read_frame_with_aux(self, frame):
869 """Move to frame, updating ts with trajectory, transformations and auxiliary data."""
--> 870 ts = self._read_frame(frame) # pylint: disable=assignment-from-no-return
871 for aux in self.aux_list:
872 ts = self._auxs[aux].update_ts(ts)

File /opt/homebrew/Caskroom/miniconda/base/envs/hoomdenv/lib/python3.13/site-packages/MDAnalysis/coordinates/GSD.py:132, in GSDReader._read_frame(self, frame)
130 def _read_frame(self, frame):
131 try:
--> 132 myframe = self._file[frame]
133 except IndexError:
134 raise IOError from None

File /opt/homebrew/Caskroom/miniconda/base/envs/hoomdenv/lib/python3.13/site-packages/gsd/hoomd.py:1058, in HOOMDTrajectory.getitem(self, key)
1055 raise IndexError()
1056 return self._read_frame(key)
-> 1058 raise TypeError

TypeError:

I think the software can successfully read .gsd because when I directly analyze u.trajectory in other occasions everything is fine. But I don't know if my code is wrong or there's some internal problems with rdf. I also attached my trajectory in this post. 2_B50_try7.gsd.zipI would really appreciate your advice. Thank you in advance!

[orbeckst]

I am not sure if you're loading your universe correctly

u = mda.Universe("2_B50_try7.gsd","2_B50_try7.gsd")

GSD files contain both topology and trajectory, so it makes some sense that you're loading the universe this way. Could you try if just u = mda.Universe("2_B50_try7.gsd") works?

Can you iterate over your trajectory, e.g.,

c = u.select_atoms("name C")
for ts in u.trajectory[:10]:
   print(ts.frame, ts.time, c.center_of_geometry())

Can you randomly access a frame

print(u.trajectory)
ts = u.trajectory[0]
print(ts.frame, ts.time)

ts = u.trajectory[-3]
print(ts.frame, ts.time)

[Lynn3512138]

Dear Oliver,

Thank you for your reply. I tested your suggestions as following:

u=mda.Universe("2_B50_try7.gsd")

c=u.select_atoms("name C")
a=u.select_atoms("name A")


for ts in u.trajectory[:10]:
   print(ts.frame, ts.time, c.center_of_geometry())

print(u.trajectory)
ts = u.trajectory[0]
print(ts.frame, ts.time)

ts = u.trajectory[-3]
print(ts.frame, ts.time)


#rdf = InterRDF(c,a,range=(0,5),start=25,end=28)
#rdf.run()
#
#plt.plot(rdf.bins, rdf.rdf)

The results are given:

0 0.0 [ 7.23059524 -1.29476592 -6.99141042]
1 1.0 [ 7.27865256 -1.106382 -6.84938921]
2 2.0 [ 7.47069027 -1.18301824 -6.81299074]
3 3.0 [ 7.26618802 -1.29809961 -6.68284354]
4 4.0 [ 7.32785551 -1.13341289 -6.77474081]
5 5.0 [ 7.28747968 -0.90813285 -6.67149149]
6 6.0 [ 7.18111347 -1.04081353 -6.72767919]
7 7.0 [ 7.14965885 -0.87958148 -6.65188987]
8 8.0 [ 7.19026589 -0.68062474 -6.71825575]
9 9.0 [ 7.09201965 -0.87805006 -6.60678759]
<GSDReader 2_B50_try7.gsd with 50 frames of 100 atoms>
0 0.0
47 47.0
/var/folders/jj/mz6sd8bj3_976yy7bd95mp9c0000gr/T/ipykernel_6436/2506976544.py:8: UserWarning: Reader has no dt information, set to 1.0 ps
print(ts.frame, ts.time, c.center_of_geometry())
/var/folders/jj/mz6sd8bj3_976yy7bd95mp9c0000gr/T/ipykernel_6436/2506976544.py:12: UserWarning: Reader has no dt information, set to 1.0 ps
print(ts.frame, ts.time)
/var/folders/jj/mz6sd8bj3_976yy7bd95mp9c0000gr/T/ipykernel_6436/2506976544.py:15: UserWarning: Reader has no dt information, set to 1.0 ps
print(ts.frame, ts.time)

However the rdf still didn't work.

[orbeckst]

I tried it with another analysis function (RMSD) and get the same error:

from MDAnalysis.analysis import rms
rr = rms.RMSD(c)
rr.run()   # fails with TypeError as in original report

The issue appears to be that we are trying to access frames in the GSD file return self._read_frame(key) with a key = np.int64(0). The analysis classes create numpy arrays with ints and they have a special data type that may not be compatible with the GSD reader.

I could verify this hypothesis with

import numpy as np
import gsd
ht = gsd.hoomd.HOOMDTrajectory(gsd.fl.open(name="./2_B50_try7.gsd", mode="r"))

# the following succeeds
ht[0]

# indexing with np.int64(0) fails with TypeError
ht[np.int64(0)]

This demonstrates that this is bug in MDAnalysis.

[orbeckst]

I raised issue #5224 .

[orbeckst]

@Lynn3512138 I am sorry, until #5224 is fixed you might have to use an older version of MDAnalysis (try 2.9.0 – I am guessing that the issue was introduced with adding parallelization to RDF).

Or convert your gsd to another topology and trajectory file, e.g., TRR (should always work) or H5MD (which is fancier and faster but may require additional installed packages):

u = mda.Universe("2_B50_try7.gsd")
u.residues.resnames = "UNK"   # PDB needs string as a resname
u.atoms.write("2_B50_try7.pdb")
u.atoms.write("2_B50_try7.trr", frames="all")  # TRR should always work
u.atoms.write("2_B50_try7.h5md", frames="all", timeunit="ps", lengthunit="A") # may require hdf5/h5py/h5md packages

and then use

u2 = mda.Universe("2_B50_try7.pdb", "2_B50_try7.h5md")

to do your analysis with u2:

from MDAnalysis.analysis.rdf import InterRDF
import matplotlib.pyplot as plt

c = u2.select_atoms("name C")
a = u2.select_atoms("name A")

rdf = InterRDF(c,a,range=(0,5))
rdf.run(verbose=True)
plt.plot(rdf.results.bins, rdf.results.rdf)

(PS: Use rdf.results!)

[Lynn3512138]

Dear Oliver,

Thank you very much! I'll try converting the .gsd into other formats. Yet I first used v2.8.0 and encountered this problem. Then I updated to v2.10.0 and found that the problem still existed, after which I posted this question. Thank you for your suggestion and development of the software!

[orbeckst]

Good to know that the bug existed earlier. Thanks for letting us know.

[orbeckst]

Did any of the above answer your question @Lynn3512138 ? (If yes, please "Mark as answer" as it helps others using the discussion as a resource. If no, let me know what remains unclear. Thanks!)

[orbeckst]

The bug #5224 was fixed in PR #5260 . If you install the latest development version then it should work for you.