Fix for NetCDF files with no time (#4074)

* attempt to deal with netcdf without time

* NetCDF update

* edited to my real name

* time is now setted based on frame index

* added the tests for cpptraj netcdf trajs

* shrinked the cpptraj traj test files

* if time is not set, use default time assignment

* rely just on one traj instead of two

* made the linter happyer

* removed unnecessary commented code

Co-authored-by: Irfan Alibay <[email protected]>

* renamed misnamed method

Co-authored-by: Irfan Alibay <[email protected]>

* Update testsuite/MDAnalysisTests/coordinates/test_netcdf.py

Co-authored-by: Irfan Alibay <[email protected]>

* Update package/CHANGELOG

Co-authored-by: Irfan Alibay <[email protected]>

* Update testsuite/MDAnalysisTests/datafiles.py

Co-authored-by: Irfan Alibay <[email protected]>

* Update test_netcdf.py

* fixed and improved dt checking

* warn the user about missing time

---------

Co-authored-by: Irfan Alibay <[email protected]>

package/AUTHORS

@@ -209,6 +209,7 @@ Chronological list of authors
   - Vishal Parmar
   - Moritz Schaeffler
   - Xu Hong Chen
+  - Domenico Marson
 
 External code
 -------------

package/CHANGELOG

@@ -16,10 +16,11 @@ The rules for this file:
 
 ??/??/?? IAlibay, pgbarletta, mglagolev, hmacdope, manuel.nuno.melo, chrispfae, 
          ooprathamm, MeetB7, BFedder, v-parmar, MoSchaeffler, jbarnoud, jandom,
-         xhgchen, jaclark5
+         xhgchen, jaclark5, DrDomenicoMarson
  * 2.5.0
 
 Fixes
+  * Fix for NetCDF trajectories without `time` variable (Issue #4073)
   * Allows shape_parameter and asphericity to yield per residue quantities
     (Issue #3002, PR #3905)
   * Add tests for "No path data" exception raise in test_psa.py (Issue #4036)

package/MDAnalysis/coordinates/TRJ.py

@@ -528,7 +528,7 @@ def __init__(self, filename, n_atoms=None, mmap=None, **kwargs):
             # len(dimensions['frame']) ==  10: in any case, we need to get
             # the number of frames from somewhere such as the time variable:
             if self.n_frames is None:
-                self.n_frames = self.trjfile.variables['time'].shape[0]
+                self.n_frames = self.trjfile.variables['coordinates'].shape[0]
         except KeyError:
             errmsg = (f"NCDF trajectory {self.filename} does not contain "
                       f"frame information")
@@ -544,7 +544,17 @@ def __init__(self, filename, n_atoms=None, mmap=None, **kwargs):
 
         # checks for not-implemented features (other units would need to be
         # hacked into MDAnalysis.units)
-        self._verify_units(self.trjfile.variables['time'].units, 'picosecond')
+        try:
+            self._verify_units(self.trjfile.variables['time'].units, 'picosecond')
+            self.has_time = True
+        except KeyError:
+            self.has_time = False
+            wmsg = ("NCDF trajectory does not contain `time` information;"
+                    " `time` will be set as an increasing index")  
+            warnings.warn(wmsg)
+            logger.warning(wmsg)
+
+
         self._verify_units(self.trjfile.variables['coordinates'].units,
                            'angstrom')
 
@@ -640,7 +650,8 @@ def _read_frame(self, frame):
                 self.n_frames))
         # note: self.trjfile.variables['coordinates'].shape == (frames, n_atoms, 3)
         ts._pos[:] = self._get_var_and_scale('coordinates', frame)
-        ts.time = self._get_var_and_scale('time', frame)
+        if self.has_time:
+            ts.time = self._get_var_and_scale('time', frame)
         if self.has_velocities:
             ts._velocities[:] = self._get_var_and_scale('velocities', frame)
         if self.has_forces:
@@ -685,7 +696,7 @@ def _get_dt(self):
         try:
             t1 = self.trjfile.variables['time'][1]
             t0 = self.trjfile.variables['time'][0]
-        except IndexError:
+        except (IndexError, KeyError):
             raise AttributeError
         return t1 - t0
 

testsuite/MDAnalysisTests/coordinates/test_netcdf.py

@@ -35,7 +35,8 @@
 
 from MDAnalysisTests.datafiles import (PFncdf_Top, PFncdf_Trj,
                                        GRO, TRR, XYZ_mini,
-                                       PRM_NCBOX, TRJ_NCBOX, DLP_CONFIG)
+                                       PRM_NCBOX, TRJ_NCBOX, DLP_CONFIG,
+                                       CPPTRAJ_TRAJ_TOP, CPPTRAJ_TRAJ)
 from MDAnalysisTests.coordinates.test_trj import _TRJReaderTest
 from MDAnalysisTests.coordinates.reference import (RefVGV, RefTZ2)
 from MDAnalysisTests import make_Universe
@@ -299,6 +300,34 @@ def test_box(self, universe, index, expected):
         assert_almost_equal(self.box_refs[expected], universe.dimensions)
 
 
+class TestNCDFReader4(object):
+    """NCDF Trajectory exported by cpptaj, without `time` variable."""
+    prec = 3
+
+    @pytest.fixture(scope='class')
+    def u(self):
+        return mda.Universe(CPPTRAJ_TRAJ_TOP,
+                            [CPPTRAJ_TRAJ, CPPTRAJ_TRAJ])
+
+    def test_chain_times(self, u):
+        """Check times entries for a chain of trajectories without
+        a defined time variable"""
+        ref_times = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0]
+        time_list = [ts.time for ts in u.trajectory]
+        assert ref_times == time_list
+
+    def test_dt(self, u):
+        ref = 1.0
+        assert u.trajectory.dt == pytest.approx(ref)
+        assert u.trajectory.ts.dt == pytest.approx(ref)
+
+    def test_warn_user_no_time_information(self, u):
+        wmsg = ("NCDF trajectory does not contain `time` information;"
+                " `time` will be set as an increasing index")  
+        with pytest.warns(UserWarning, match=wmsg):
+            u2 = mda.Universe(CPPTRAJ_TRAJ_TOP, CPPTRAJ_TRAJ)
+
+
 class _NCDFGenerator(object):
     """A class for generating abitrary ncdf files and exhaustively test
     edge cases which might not be found in the wild"""