Questions about MDAnalysis H5MD implementation vs. H5MD specification

[ljwoods2]

After discussing with @orbeckst, it seems like there might be some undocumented limitations with the H5MD implementation, but I'm not 100% sure. I would love some feedback from someone more familiar.

First, reader limitations:

1. Reader and sampling rate

The H5MD spec (https://www.nongnu.org/h5md/h5md.html#particles-group) defines step like this:

step
A dataset with dimensions [variable] that contains the time steps at which the corresponding data were sampled. It is of Integer type to allow exact temporal matching of data from one H5MD element to another. The values of the dataset are in monotonically increasing order.

From my understanding, this means it should be possible to have an h5md file in which velocities and positions are sampled at different rates for example. However, when I try to test this, using this code:

import MDAnalysis as mda
import h5py
from MDAnalysis.tests.datafiles import COORDINATES_TOPOLOGY, COORDINATES_H5MD

with h5py.File(COORDINATES_H5MD, "a") as f:
    with h5py.File("testh5md.h5md", 'w') as new_file:  # Open or create the new file in write mode
        # Copy each item from the original file to the new file
        for item in f.keys():
            f.copy(item, new_file)

        del new_file["particles/trajectory/velocity/step"]
        del new_file["particles/trajectory/velocity/time"]
        # step was previously [0, 1, 2, 3 ,4]
        new_file.create_dataset("particles/trajectory/velocity/step", data=[1,3], dtype=int)
        # time was previously [0, 1, 2, 3, 4]
        new_file.create_dataset("particles/trajectory/velocity/time", data=[1,3], dtype=float)

        vel = new_file["particles/trajectory/velocity/value"][[1, 3]]
        del new_file["particles/trajectory/velocity/value"]
        new_file.create_dataset("particles/trajectory/velocity/value", data=vel, dtype=float)
        new_file["particles/trajectory/velocity/value"].attrs["unit"] = "Angstrom ps-1"

u = mda.Universe(COORDINATES_TOPOLOGY, "testh5md.h5md")
ts = u.trajectory[0]
print(ts.data["step"])
print(ts.positions)
print(ts.velocities)
ts = u.trajectory[1]
print(ts.data["step"])
print(ts.positions)
print(ts.velocities)
ts = u.trajectory[2]
print(ts.data["step"])
print(ts.positions)
print(ts.velocities)

I don't get what I would intuitively expect. It seems like the resulting trajectory should have 5 frames, and the second and fourth frame should contain velocities. However, this code prints velocity data for the first two frames and then fails on the third frame with

IndexError: Index (2) out of range for (0-1)

I think the H5MD reader was implemented to only allow files where positions, velocities, and forces were all sampled at the same rate and therefore share the same step and time arrays. Maybe this is documented somewhere already that I'm missing, not sure.

2. Writer and sampling rate

The writer has a similar issue where the self._has dictionary which is used in each write call to decide whether positions, velocities, and forces should be written is set during the first write call and never updated. This means that if the trajectory being written does not have position/velocity/force in the first time step, it won't be written at all. Again, only an issue if the trajectory has different sampling rates.

3. Multiple trajectories per file

The H5MD spec allows for multiple subgroups of the 'particles' group:

Particles group
Information on each particle, i.e., particle trajectories, is stored in the particles group. The particles group is a container for subgroups that represent different subsets of the system under consideration, and it may hold one or several subgroups, as needed. These subsets may overlap and their union may be incomplete, i.e., not represent all particles of the simulation volume. The subgroups contain the trajectory data for each particle as time-dependent or time-independent data, depending on the situation. Each subgroup contains a specification of the simulation box, see below. For each dataset, the particle index is accommodated by the second (first, in the case of time-independence) array dimension.

The reader simply pulls out the first key, however, and I'm not sure this is documented anywhere:

 # pulls first key out of 'particles'
        # allows for arbitrary name of group1 in 'particles'
        self._particle_group = self._file['particles'][
            list(self._file['particles'])[0]]

Let me know if I'm missing anything or if I should just make a documentation PR!

[orbeckst]

@edisj would be great if you could have a look here!

[edisj]

Hi @ljwoods2 , I can take a look at this over the weekend. Sorry it's been a little while since I've looked at H5MDReader stuff

[edisj]

Hey @ljwoods2 , thanks for looking into this. After reviewing the H5MD documentation, I think I agree with all three of your points. Admittedly, when I first wrote H5MDReader and H5MDWriter, I made some decisions so that I could produce something that worked and didn't violate the H5MD standard, but you're totally right that the current implementation is restricting the total number of possible H5MD files that MDAnalysis can read. I think we should fix this! H5MD in MDAnalysis has always been an unfinished project lurking in the very back of my mind, so I guess there's no time like the present to start chipping away at it...

I'm going to open two new discussions to focus on each of your points individually, namely:
1) How should H5MDReader/Writer handle an input trajectory with data sampled at different rates?
2) How should H5MDReader handle an input file that contains multiple trajectories?

this way, we can have a discussion -> issue -> pull request history for each problem separately, which I think is much cleaner on the development side and easier for code reviewers

[orbeckst]

I’d initially just focus on documenting the known limitations and possibly fail explicitly if an unsupported file layout is encountered.

…

Am 4/7/24 um 14:05 schrieb edisj ***@***.***>: I'm going to open two new discussions to focus on each of your points individually

[edisj]

is there anything I should do for this in the immediate future, or is @ljwoods2 on the case?
edit: never mind, I opened an issue

[orbeckst]

Thanks for raising the issue!I assume that Lawson is working on streaming but please talk to each other. Am 4/7/24 um 14:23 schrieb edisj ***@***.***>: is there anything I should do for this in the immediate future, or is @ljwoods2 on the case? —Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you were mentioned.Message ID: ***@***.***>