Issue reading Gromacs written TNG file containing both positions and velocities

[conor-b-abraham]

Hi All,

I've been playing around with trying to convert some huge TRR trajectory files (with positions and velocities) into TNG files using Gromacs 2021.5 trjconv in order to reduce the file sizes, but am having trouble reading those files into MDAnalysis 2.7.0. For example, when I do: u = mda.Universe("4.2_production.tpr", "4.2_production.tng") I get the following error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[3], line 1
----> 1 u = mda.Universe("4.2_production.tpr", "4.2_unwrapped.tng")

File /projectnb/abaggreg/conor/.conda/envs/py3/lib/python3.9/site-packages/MDAnalysis/core/universe.py:375, in Universe.__init__(self, topology, all_coordinates, format, topology_format, transformations, guess_bonds, vdwradii, fudge_factor, lower_bound, in_memory, in_memory_step, *coordinates, **kwargs)
    370 coordinates = _resolve_coordinates(self.filename, *coordinates,
    371                                    format=format,
    372                                    all_coordinates=all_coordinates)
    374 if coordinates:
--> 375     self.load_new(coordinates, format=format, in_memory=in_memory,
    376                 in_memory_step=in_memory_step, **kwargs)
    378 if transformations:
    379     if callable(transformations):

File /projectnb/abaggreg/conor/.conda/envs/py3/lib/python3.9/site-packages/MDAnalysis/core/universe.py:580, in Universe.load_new(self, filename, format, in_memory, in_memory_step, **kwargs)
    577 # supply number of atoms for readers that cannot do it for themselves
    578 kwargs['n_atoms'] = self.atoms.n_atoms
--> 580 self.trajectory = reader(filename, format=format, **kwargs)
    581 if self.trajectory.n_atoms != len(self.atoms):
    582     raise ValueError("The topology and {form} trajectory files don't"
    583                      " have the same number of atoms!\n"
    584                      "Topology number of atoms {top_n_atoms}\n"
   (...)
    588                          fname=filename,
    589                          trj_n_atoms=self.trajectory.n_atoms))

File /projectnb/abaggreg/conor/.conda/envs/py3/lib/python3.9/site-packages/MDAnalysis/lib/util.py:2553, in store_init_arguments.<locals>.wrapper(self, *args, **kwargs)
   2551             else:
   2552                 self._kwargs[key] = arg
-> 2553 return func(self, *args, **kwargs)

File /projectnb/abaggreg/conor/.conda/envs/py3/lib/python3.9/site-packages/MDAnalysis/coordinates/TNG.py:172, in TNGReader.__init__(self, filename, convert_units, **kwargs)
    169 self.filename = filename
    170 self.convert_units = convert_units
--> 172 self._file_iterator = pytng.TNGFileIterator(self.filename, "r")
    173 self.n_atoms = self._file_iterator.n_atoms
    174 self._n_steps = self._file_iterator.n_steps

File /projectnb/abaggreg/conor/.conda/envs/py3/lib/python3.9/site-packages/pytng/pytng.pyx:639, in pytng.pytng.TNGFileIterator.__cinit__()

File /projectnb/abaggreg/conor/.conda/envs/py3/lib/python3.9/site-packages/pytng/pytng.pyx:708, in pytng.pytng.TNGFileIterator._open()

File /projectnb/abaggreg/conor/.conda/envs/py3/lib/python3.9/site-packages/pytng/pytng.pyx:959, in pytng.pytng.TNGFileIterator._sanitise_block_metadata()

ValueError: Block TNG_TRAJ_POSITIONS has no frames contaning data

but when I do u = mda.Universe("4.2_production.tpr", "4.2_production.tng") everything works fine.

I notice that when I use gmx dump, each frame in the trajectory has a data block named POSITIONS. I haven't looked at the MDA source code yet, but is this what the TNGReader is looking for when searching for TNG_TRAJ_POSITIONS? Does anyone have a solution? If not, can anyone recommend another way for me to reduce the size of these huge TRR files (>200GB), but keep them in an mdanalysis-ready format?

Thanks in advance,

Conor

[orbeckst]

In your example you failed with the "4.2_unwrapped.tng" but you said you succeeded with "4.2_production.tng". What is the difference between the files?

[orbeckst]

Also, just in case, have your tried MDAnalysis 2.9.0?

[conor-b-abraham]

Sorry my mistake, I meant to say that using 4.2_production.trr works. I have not tried MDAnalysis 2.9.0, but will give it a try. Unfortunately, 2.7 seems the highest I can use with some of the other packages I use for analysis.

[orbeckst]

If it works in a separate environment then at least we could start figuring out what changes lead to it working.

[conor-b-abraham]

For the record, using the default format for the .tng file also did not work in a new environment with MDAnalysis 2.9.0. Using chemfiles (the comment I marked as the answer) works in both MDA versions, however.

[orbeckst]

By default, MDAnalysis uses the pytng package https://github.com/MDAnalysis/pytng for handling TNG.

However, you can also use the chemfiles library for many formats, including TNG, so perhaps try

u = mda.Universe("4.2_production.tpr", "4.2_unwrapped.tng", format="CHEMFILES")

Maybe chemfiles can do it.

[conor-b-abraham]

Using chemfiles worked in both MDAnalysis 2.7.0 and 2.9.0. Thank you so much!

[orbeckst]

Glad that this worked.

I opened MDAnalysis/pytng#116 to have a record in the pytng repo.