[BUG] Trajectory reloading from Gro87 files leads to indexing errors

[lohedges]

I am hoping to use the Sire trajectory interface to parallelise the setup of ABFE replica exchange simulations in BioSimSpace. This would give a highly parallelised way of reading and writing the input and output coordinates for each replica. This seems to work fine for AMBER format coordinate files, e.g.:

In [1]: import sire as sr

In [2]: mols = sr.load("input/system.prm7", "input/traj.dcd")

In [3]: mols.num_frames()
Out[3]: 100

In [4]: sr.save(mols.trajectory(), "frame.rst7")
Out[4]: ['/home/lester/Code/openbiosim/biosimspace/abfe_setup/frame.rst7']

In [5]: new_mols = sr.load("input/system.prm7", "frame.rst7/*")

In [6]: new_mols.num_frames()
Out[6]: 100

In [7]: for frame in new_mols.trajectory()[:10]:
   ...:     print(frame)
   ...:
System( name=ACE num_molecules=631 num_residues=633 num_atoms=1912 )
System( name=ACE num_molecules=631 num_residues=633 num_atoms=1912 )
System( name=ACE num_molecules=631 num_residues=633 num_atoms=1912 )
System( name=ACE num_molecules=631 num_residues=633 num_atoms=1912 )
System( name=ACE num_molecules=631 num_residues=633 num_atoms=1912 )
System( name=ACE num_molecules=631 num_residues=633 num_atoms=1912 )
System( name=ACE num_molecules=631 num_residues=633 num_atoms=1912 )
System( name=ACE num_molecules=631 num_residues=633 num_atoms=1912 )
System( name=ACE num_molecules=631 num_residues=633 num_atoms=1912 )
System( name=ACE num_molecules=631 num_residues=633 num_atoms=1912 )

but fails when using Gro87 format:

In [1]: import sire as sr

In [2]: mols = sr.load("input/system.top", "input/traj.xtc")

In [3]: mols.num_frames()
Out[3]: 100

In [4]: sr.save(mols.trajectory(), "frame.gro")
Out[4]: ['/home/lester/Code/openbiosim/biosimspace/abfe_setup/frame.gro']

In [5]: new_mols = sr.load("input/system.top", "frame.gro/*")

In [6]: new_mols.num_frames()
Out[6]: 100

In [7]: for frame in new_mols.trajectory()[:10]:
   ...:     print(frame)
   ...:
╭─────────────────────────────── Traceback (most recent call last) ────────────────────────────────╮
│ in <module>:1                                                                                    │
│                                                                                                  │
│ /home/lester/.conda/envs/openbiosim/lib/python3.12/site-packages/sire/mol/_trajectory.py:159 in  │
│ __next__                                                                                         │
│                                                                                                  │
│    156 │   │                                                                                     │
│    157 │   │   self._frame = self._iter.__next__()                                               │
│    158 │   │                                                                                     │
│ ❱  159 │   │   return self.current()                                                             │
│    160 │                                                                                         │
│    161 │   def __len__(self):                                                                    │
│    162 │   │   return len(self._values)                                                          │
│                                                                                                  │
│ /home/lester/.conda/envs/openbiosim/lib/python3.12/site-packages/sire/mol/_trajectory.py:283 in  │
│ current                                                                                          │
│                                                                                                  │
│    280 │   │   │   map = self._map.clone()                                                       │
│    281 │   │   │   map.set("transform", self._aligner[self._frame])                              │
│    282 │   │                                                                                     │
│ ❱  283 │   │   ret.load_frame(self._frame, map=map)                                              │
│    284 │   │                                                                                     │
│    285 │   │   try:                                                                              │
│    286 │   │   │   mol = ret.molecule()                                                          │
│                                                                                                  │
│ /home/lester/.conda/envs/openbiosim/lib/python3.12/site-packages/sire/system/_system.py:185 in   │
│ load_frame                                                                                       │
│                                                                                                  │
│    182 │   │   """Load the ith frame into this System"""                                         │
│    183 │   │   from ..base import create_map                                                     │
│    184 │   │                                                                                     │
│ ❱  185 │   │   self._system.load_frame(i, map=create_map(map))                                   │
│    186 │   │   self._molecules = None                                                            │
│    187 │                                                                                         │
│    188 │   def save_frame(self, i=None, map=None):                                               │
╰──────────────────────────────────────────────────────────────────────────────────────────────────╯
ValueError: SireError::incompatible_error: Cannot subset a Frame of length 0 using start 259, count 3. (call
sire.error.get_last_error_details() for more info)

The Gro87 files were definitely written correctly. I'm wondering whether this is because the files don't contain a time stamp, which is somehow being used to index the frames?

Inputs can be found here.

[lohedges]

Weirdly, it does work if you get the trajectory iterator and access the frames directly. I'm not sure what's different about doing it in a loop:

In [8]: traj = new_mols.trajectory()

In [9]: frame1 = traj[10].current()

In [10]: frame2 = traj[20].current()

In [11]: frame_1.coordinates()
Out[11]: ( 13.3248 Å, 14.5172 Å, 12.6031 Å )

In [12: frame_2.coordinates()
Out[12]: ( 13.3263 Å, 14.7003 Å, 12.5236 Å )

The main issue for my purposes is that you can't write the trajectory back to GROMACS format, since you get the indexing error above, i.e.:

In [13]: sr.save(new_mols.trajectory(), "tmp.gro")
╭─────────────────────────────── Traceback (most recent call last) ────────────────────────────────╮
│ in <module>:1                                                                                    │
│                                                                                                  │
│ /home/lester/.conda/envs/openbiosim/lib/python3.12/site-packages/sire/_load.py:646 in save       │
│                                                                                                  │
│   643 │   │   if molecules._is_trajectory_iterator():                                            │
│   644 │   │   │   # we are saving a trajectory - not just the molecules                          │
│   645 │   │   │   map = molecules._populate_map(map)                                             │
│ ❱ 646 │   │   │   molecules = molecules.current()                                                │
│   647 │                                                                                          │
│   648 │   molecules = _to_legacy_system(molecules)                                               │
│   649                                                                                            │
│                                                                                                  │
│ /home/lester/.conda/envs/openbiosim/lib/python3.12/site-packages/sire/mol/_trajectory.py:283 in  │
│ current                                                                                          │
│                                                                                                  │
│    280 │   │   │   map = self._map.clone()                                                       │
│    281 │   │   │   map.set("transform", self._aligner[self._frame])                              │
│    282 │   │                                                                                     │
│ ❱  283 │   │   ret.load_frame(self._frame, map=map)                                              │
│    284 │   │                                                                                     │
│    285 │   │   try:                                                                              │
│    286 │   │   │   mol = ret.molecule()                                                          │
│                                                                                                  │
│ /home/lester/.conda/envs/openbiosim/lib/python3.12/site-packages/sire/system/_system.py:185 in   │
│ load_frame                                                                                       │
│                                                                                                  │
│    182 │   │   """Load the ith frame into this System"""                                         │
│    183 │   │   from ..base import create_map                                                     │
│    184 │   │                                                                                     │
│ ❱  185 │   │   self._system.load_frame(i, map=create_map(map))                                   │
│    186 │   │   self._molecules = None                                                            │
│    187 │                                                                                         │
│    188 │   def save_frame(self, i=None, map=None):                                               │
╰──────────────────────────────────────────────────────────────────────────────────────────────────╯
ValueError: SireError::incompatible_error: Cannot subset a Frame of length 0 using start 307, count 3. (call
sire.error.get_last_error_details() for more info)

[lohedges]

Specifically, the error comes from calling load_frames on the underlying Sire system, e.g.:

In [14]: new_mols._system.load_frame(0)
╭─────────────────────────────── Traceback (most recent call last) ────────────────────────────────╮
│ in <module>:1                                                                                    │
╰──────────────────────────────────────────────────────────────────────────────────────────────────╯
ValueError: SireError::incompatible_error: Cannot subset a Frame of length 0 using start 307, count 3. (call
sire.error.get_last_error_details() for more info)

I'm really confused why this exception doesn't trigger in the cases shown above, i.e. getting a frame from the trajectory iterator directly and calling .current(), since this code path is also used.

[lohedges]

Okay, it turns out that I was using the wrong trajectory above. The main issue here is that the Gro87::getFrame() isn't implemented, i.e.

sire/corelib/src/libs/SireIO/gro87.cpp

Lines 692 to 697 in e77ddd1

	Frame Gro87::getFrame(int i) const
	{
	i = SireID::Index(i).map(this->nFrames());
	return SireMol::Frame();
	}

This just returns a null Frame, with no atoms, hence the indexing errors seen above. It took me about an hour of working through debugging statements in the Sire::Mol::Trajectory code to work out that this was actually where the frame was being returned from. I'm not sure if this bit of code wasn't added for a reason, or whether it was just forgotten. I'll work through SireIO::AmberRst7::getFrame() to see if I can copy the logic across.

[lohedges]

Now fixed in e2b9f62. Will test further and begin to think about and develop ways to sort trajectory frames, i.e. tackle #191.

[lohedges]

Note that the linked commit also raises an unsupported error for Gro87::operator[](int i) since the existing implementation was broken.

[lohedges]

Okay, I'll certainly need to add sorting of some kind since the GROMACS frames aren't returned in the same order as the AMBER ones, which have a time stamp. The order is reproducible, but not correct.

In [1]: import sire as sr

In [2]: amb_mols = sr.load("system.prm7", "frame.rst7/*")

In [3]: gro_mols = sr.load("system.top", "frame.gro/*")

In [4]: for amb_frame, amb_frame in zip(amb_mols.trajectory()[:10], gro_mols.trajectory()[:10]):
   ...:     print(amb_frame.coordinates(), gro_frame.coordinates())
   ...:
( 13.4532 Å, 14.3844 Å, 12.404 Å ) ( 13.2826 Å, 14.3805 Å, 12.4475 Å )
( 13.2404 Å, 14.6071 Å, 12.6041 Å ) ( 13.4123 Å, 14.3929 Å, 12.4422 Å )
( 13.4123 Å, 14.3929 Å, 12.4422 Å ) ( 13.2445 Å, 14.3434 Å, 12.4812 Å )
( 13.3309 Å, 14.4438 Å, 12.5948 Å ) ( 13.329 Å, 14.3945 Å, 12.4025 Å )
( 13.2864 Å, 14.6141 Å, 12.7545 Å ) ( 13.2406 Å, 14.4258 Å, 12.6037 Å )
( 13.0802 Å, 14.6261 Å, 12.4385 Å ) ( 13.2405 Å, 14.3347 Å, 12.604 Å )
( 13.4157 Å, 14.5373 Å, 12.2817 Å ) ( 13.3309 Å, 14.4438 Å, 12.5948 Å )
( 13.3726 Å, 14.489 Å, 12.439 Å ) ( 13.5411 Å, 14.6278 Å, 12.4376 Å )
( 13.0774 Å, 14.391 Å, 12.3642 Å ) ( 13.4955 Å, 14.3439 Å, 12.4422 Å )
( 13.4974 Å, 14.5781 Å, 12.4012 Å ) ( 13.4142 Å, 14.352 Å, 12.4005 Å )

[lohedges]

Okay, it works if you just point sr.load at the directory containing the files, i.e.:

In [1]: import sire as sr

In [2]: amb_mols = sr.load("system.prm7", "frame.rst7")

In [3]: gro_mols = sr.load("system.top", "frame.gro")

In [4]: for amb_frame, amb_frame in zip(amb_mols.trajectory()[:10], gro_mols.trajectory()[:10]):
   ...:     print(amb_frame.coordinates(), gro_frame.coordinates())
   ...:
( 13.4532 Å, 14.3844 Å, 12.404 Å ) ( 13.4532 Å, 14.3844 Å, 12.404 Å )
( 13.2404 Å, 14.6071 Å, 12.6041 Å ) ( 13.2404 Å, 14.6071 Å, 12.6041 Å ) 
( 13.4123 Å, 14.3929 Å, 12.4422 Å ) ( 13.4123 Å, 14.3929 Å, 12.4422 Å )
( 13.3309 Å, 14.4438 Å, 12.5948 Å ) ( 13.3309 Å, 14.4438 Å, 12.5948 Å )
( 13.2864 Å, 14.6141 Å, 12.7545 Å ) ( 13.2864 Å, 14.6141 Å, 12.7545 Å ) 
( 13.0802 Å, 14.6261 Å, 12.4385 Å ) ( 13.0802 Å, 14.6261 Å, 12.4385 Å )
( 13.4157 Å, 14.5373 Å, 12.2817 Å ) ( 13.4157 Å, 14.5373 Å, 12.2817 Å ) 
( 13.3726 Å, 14.489 Å, 12.439 Å )  ( 13.3726 Å, 14.489 Å, 12.439 Å ) 
( 13.0774 Å, 14.391 Å, 12.3642 Å ) ( 13.0774 Å, 14.391 Å, 12.3642 Å ) 
( 13.4974 Å, 14.5781 Å, 12.4012 Å ) ( 13.4974 Å, 14.5781 Å, 12.4012 Å ) 

The issue here is that we don't just want to point to a directory, since that would contain other files that we wouldn't want to parse. I'll figure out how I can reproduce the same sorting logic when you just glob files, or give it a list of file names. Ideally we could just give it a pre-sorted list, then it loads the frames (or sorts them internally) based on the specified ordering.

[lohedges]

I was wrong, sorting the file names does work, so that's half the problem solved, i.e. if there are no time-stamps, then sorting the input files is all that's needed to get a trajectory with arbitrary order. Now to figure out how to write to a specific path.

In [1]: import sire as sr

In [2]: import glob

In [3]: amb_mols = sr.load("system.prm7", "frame.rst7")

In [4]: gro_frames = sorted(glob.glob("frame.gro/*"))

In [5]: gro_mols = sr.load("system.top", *gro_frames)

In [6]: for amb_frame, gro_frame in zip(amb_mols.trajectory()[:10], gro_mols.trajectory()[:10]):
   ...:     print(amb_frame.coordinates(), gro_frame.coordinates())
   ...:
( 13.1149 Å, 14.5162 Å, 12.6433 Å ) ( 13.1149 Å, 14.5162 Å, 12.6433 Å )
( 13.4494 Å, 14.2891 Å, 12.6433 Å ) ( 13.4494 Å, 14.2891 Å, 12.6433 Å )
( 13.1985 Å, 14.4708 Å, 12.8382 Å ) ( 13.1985 Å, 14.4708 Å, 12.8382 Å )
( 13.2404 Å, 14.6071 Å, 12.6041 Å ) ( 13.2404 Å, 14.6071 Å, 12.6041 Å )
( 13.2405 Å, 14.3347 Å, 12.604 Å ) ( 13.2405 Å, 14.3347 Å, 12.604 Å )
( 13.0733 Å, 14.6073 Å, 12.5649 Å ) ( 13.0733 Å, 14.6073 Å, 12.5649 Å )
( 13.2406 Å, 14.4258 Å, 12.6037 Å ) ( 13.2406 Å, 14.4258 Å, 12.6037 Å )
( 13.2826 Å, 14.3805 Å, 12.4475 Å ) ( 13.2826 Å, 14.3805 Å, 12.4475 Å )
( 13.3245 Å, 14.2898 Å, 12.5644 Å ) ( 13.3245 Å, 14.2898 Å, 12.5644 Å )
( 13.2828 Å, 14.6079 Å, 12.5642 Å ) ( 13.2828 Å, 14.6079 Å, 12.5642 Å )

[lohedges]

Okay, I've now hacked the MoleculeParser.write() code to allow writing to user specified trajectory frame names via the map option. This now allows:

In [1]: import sire as sr

In [2]: import glob

In [3]: gro_frames = sorted(glob.glob("frame.gro/*"))

In [4]: gro_mols = sr.load("system.top", *gro_frames)

In [5]: sr.save(gro_mols.trajectory()[:2], "tmp0.gro")
Out[5]: ['/home/lester/Code/openbiosim/sire/trajectory_debug/tmp0.gro']

In [6]: sr.save(gro_mols.trajectory()[:2], "tmp0.gro", map={"frame_names":"cat/cat.gro,dog/dog.gro"})
Out[6]:
['/home/lester/Code/openbiosim/sire/trajectory_debug/cat/cat.gro',
 '/home/lester/Code/openbiosim/sire/trajectory_debug/dog/dog.gro']

In the first instance we write two frames to a name called tmp.gro. In this case Sire does its usual thing and creates a frame directory called tmp.gro and writes two frames in there. In the second case, we explicitly specify the names via the map kwarg, so Sire writes to those instead. There is internal validation to make sure that the number of frame names matches the number of frames.

@chryswoods: What do you think about this approach? I could obviously make this more general and robust, but it pretty much does everything that I want and I think everything else can be handled from the Python API side, e.g. we could allow a user to pass a list of filenames to sr.save when the object they are saving is a trajectory, etc. This would then be wrapped by BioSimSpace to create a MultiSystem object to allow fast setup and tear down of ABFE repex simulations.

[lohedges]

I would still like a way for the user to control how frames are read, e.g. in the case of time-stamps being present but the actual file name being the meaningful think to sort on.

[lohedges]

The main complication in doing something more fancy is that MoleculeParser:: write is static. Hence using the map and hacking the filename passed through to MoleculeParser::writeToFile.

[chryswoods]

It's a hacky solution, but you are right that the design is a bit rigid here...