Merge branch 'develop' into hbond_update

.github/workflows/gh-ci.yaml

@@ -122,8 +122,9 @@ jobs:
 
     - name: codecov
       if: matrix.codecov
-      uses: codecov/codecov-action@v2
+      uses: codecov/codecov-action@v3
       with:
+        token: ${{ secrets.CODECOV_TOKEN }}
         file: coverage.xml
         fail_ci_if_error: True
         verbose: True

azure-pipelines.yml

@@ -26,7 +26,6 @@ jobs:
           PYTHON_ARCH: 'x86'
           BUILD_TYPE: 'normal'
           imageName: 'windows-2019'
-          PIP_PRE_EXTRAS: 'scipy==1.9.1'
         Win-Python38-64bit-full:
           PYTHON_VERSION: '3.8'
           PYTHON_ARCH: 'x64'
@@ -81,11 +80,8 @@ jobs:
     condition: and(succeeded(), eq(variables['BUILD_TYPE'], 'wheel'))
   - script: python -m pip install --upgrade pip setuptools wheel
     displayName: 'Install tools'
-  - script: python -m pip install $(PIP_PRE_EXTRAS)
-    displayName: 'Pre-deps extras'
-    condition: ne(variables['PIP_PRE_EXTRAS'], '')
   - script: >-
-      python -m pip install --only-binary=h5py
+      python -m pip install --only-binary=h5py,scipy
       cython
       hypothesis
       matplotlib
@@ -95,7 +91,6 @@ jobs:
       pytest-cov
       pytest-xdist
       scikit-learn
-      scipy
       h5py>=2.10
       tqdm
       threadpoolctl

package/AUTHORS

@@ -198,6 +198,7 @@ Chronological list of authors
   - Raymond Zhao
   - Haleema Khan
   - Jennifer A Clark
+  - Jake Fennick
 
 
 External code

package/CHANGELOG

@@ -14,11 +14,12 @@ The rules for this file:
 
 ------------------------------------------------------------------------------
 ??/??/?? IAlibay, Luthaf, hmacdope, rafaelpap, jbarnoud, BFedder, aya9aladdin,
-         jaclark5
+         jaclark5, jfennick
 
  * 2.4.0
 
 Fixes
+  * LAMMPSDump Reader translates the box to the origin (#3741)
   * hbond analysis: access hbonds results through new results member in count_by_ids() and count_by_type()
   * Added isolayer selection method (Issue #3845)
   * Auxiliary; determination of representative frames: Removed undefined
@@ -33,6 +34,9 @@ Enhancements
   * Update hbond analysis doc string to use exclusive bound language (#3847)
   * Added ability for hbond anaylsis to use types when resnames aren't there
   * Added explanatory warnings when hbond analysis doesn't find any hbonds
+  * LAMMPSDump Reader optionally unwraps trajectories with image flags upon 
+    loading (#3843
+  * LAMMPSDump Reader now imports velocities and forces (#3843)
   * Minimum NumPy version for py3.11 is now set to 1.23.2.
   * Added decorator to check for an empty atom group, applied in topological 
     attributes(Issue #3837)
@@ -61,6 +65,8 @@ Changes
   * adding element attribute to TXYZParser if all atom names are valid element symbols (PR #3826)
 
 Deprecations
+  * Add deprecation warning for inclusive stop indexing in
+    MemoryReader.timeseries (#PR 3894)
 
 
 08/29/22 IAlibay, PicoCentauri, orbeckst, hmacdope, rmeli, miss77jun, rzhao271,

package/MDAnalysis/analysis/encore/cutils.pyx

@@ -20,7 +20,7 @@
 # MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations.
 # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787
 #
-#cython embedsignature=True
+# cython embedsignature=True
 
 """
 Mixed Cython utils for ENCORE
@@ -38,19 +38,25 @@ from libc.math cimport sqrt
 
 np.import_array()
 
+
 @cython.boundscheck(False)
 @cython.wraparound(False)
-def PureRMSD(np.ndarray[np.float64_t,ndim=2] coordsi,
-             np.ndarray[np.float64_t,ndim=2] coordsj,
+def PureRMSD(np.ndarray[np.float64_t, ndim=2] coordsi,
+             np.ndarray[np.float64_t, ndim=2] coordsj,
              int atomsn,
-             np.ndarray[np.float64_t,ndim=1] masses,
+             np.ndarray[np.float64_t, ndim=1] masses,
              double summasses):
 
     cdef  int k
-    cdef double normsum, totmasses
+    cdef double normsum
 
     normsum = 0.0
 
     for k in xrange(atomsn):
-        normsum += masses[k]*((coordsi[k,0]-coordsj[k,0])**2 + (coordsi[k,1]-coordsj[k,1])**2 + (coordsi[k,2]-coordsj[k,2])**2)
+        normsum += masses[k] * (
+            (coordsi[k, 0] - coordsj[k, 0])**2 +
+            (coordsi[k, 1] - coordsj[k, 1])**2 +
+            (coordsi[k, 2] - coordsj[k, 2])**2
+        )
+
     return sqrt(normsum/summasses)

package/MDAnalysis/coordinates/LAMMPS.py

@@ -452,25 +452,51 @@ def write(self, selection, frame=None):
 
 
 class DumpReader(base.ReaderBase):
-    """Reads the default `LAMMPS dump format`_
-
-    Supports coordinates in the LAMMPS "unscaled" (x,y,z), "scaled" (xs,ys,zs),
-    "unwrapped" (xu,yu,zu) and "scaled_unwrapped" (xsu,ysu,zsu) coordinate
-    conventions (see https://docs.lammps.org/dump.html for more details).
-    If `lammps_coordinate_convention='auto'` (default),
-    one will be guessed. Guessing checks whether the coordinates fit each convention in the order "unscaled",
-    "scaled", "unwrapped", "scaled_unwrapped" and whichever set of coordinates
-    is detected first will be used. If coordinates are given in the scaled
-    coordinate convention (xs,ys,zs) or scaled unwrapped coordinate convention
-    (xsu,ysu,zsu) they will automatically be converted from their
-    scaled/fractional representation to their real values.
-
-    Supports both orthogonal and triclinic simulation box dimensions (for more
-    details see https://docs.lammps.org/Howto_triclinic.html). In either case,
-    MDAnalysis will always use ``(*A*, *B*, *C*, *alpha*, *beta*, *gamma*)``
-    to represent the unit cell. Lengths *A*, *B*, *C* are in the MDAnalysis
-    length unit (Å), and angles are in degrees.
-
+    """Reads the default `LAMMPS dump format 
+    <https://docs.lammps.org/dump.html>`__
+
+    Supports coordinates in various LAMMPS coordinate conventions and both 
+    orthogonal and triclinic simulation box dimensions (for more details see 
+    `documentation <https://docs.lammps.org/Howto_triclinic.html>`__). In 
+    either case, MDAnalysis will always use ``(*A*, *B*, *C*, *alpha*, *beta*,
+    *gamma*)`` to represent the unit cell. Lengths *A*, *B*, *C* are in the 
+    MDAnalysis length unit (Å), and angles are in degrees.
+
+    Parameters
+    ----------
+    filename : str
+        Filename of LAMMPS dump file
+    lammps_coordinate_convention : str (optional) default="auto"
+        Convention used in coordinates, can be one of the following according
+        to the `LAMMPS documentation <https://docs.lammps.org/dump.html>`__:
+
+         - "auto" - Detect coordinate type from file column header. If auto 
+           detection is used, the guessing checks whether the coordinates
+           fit each convention in the order "unscaled", "scaled", "unwrapped", 
+           "scaled_unwrapped" and whichever set of coordinates is detected 
+           first will be used.
+         - "scaled" - Coordinates wrapped in box and scaled by box length (see
+            note below), i.e., xs, ys, zs
+         - "scaled_unwrapped" - Coordinates unwrapped and scaled by box length,
+           (see note below) i.e., xsu, ysu, zsu
+         - "unscaled" - Coordinates wrapped in box, i.e., x, y, z
+         - "unwrapped" - Coordinates unwrapped, i.e., xu, yu, zu
+
+        If coordinates are given in the scaled coordinate convention (xs,ys,zs)
+        or scaled unwrapped coordinate convention (xsu,ysu,zsu) they will 
+        automatically be converted from their scaled/fractional representation
+        to their real values.
+    unwrap_images : bool (optional) default=False
+        If `True` and the dump file contains image flags, the coordinates 
+        will be unwrapped. See `read_data 
+        <https://docs.lammps.org/read_data.html>`__  in the lammps 
+        documentation for more information.
+    **kwargs
+       Other keyword arguments used in :class:`~MDAnalysis.coordinates.base.ReaderBase`
+
+    .. versionchanged:: 2.4.0
+       Now imports velocities and forces, translates the box to the origin,
+       and optionally unwraps trajectories with image flags upon loading.
     .. versionchanged:: 2.2.0
        Triclinic simulation boxes are supported.
        (Issue `#3383 <https://github.com/MDAnalysis/mdanalysis/issues/3383>`__)
@@ -489,7 +515,9 @@ class DumpReader(base.ReaderBase):
     }
 
     @store_init_arguments
-    def __init__(self, filename, lammps_coordinate_convention="auto",
+    def __init__(self, filename, 
+                 lammps_coordinate_convention="auto",
+                 unwrap_images=False,
                  **kwargs):
         super(DumpReader, self).__init__(filename, **kwargs)
 
@@ -503,6 +531,8 @@ def __init__(self, filename, lammps_coordinate_convention="auto",
                              " is not a valid option. "
                              f"Please choose one of {option_string}")
 
+        self._unwrap = unwrap_images
+
         self._cache = {}
 
         self._reopen()
@@ -606,9 +636,27 @@ def _read_next_timestep(self):
         convention_to_col_ix = {}
         for cv_name, cv_col_names in self._coordtype_column_names.items():
             try:
-                convention_to_col_ix[cv_name] = [attr_to_col_ix[x] for x in cv_col_names]
+                convention_to_col_ix[cv_name] = [attr_to_col_ix[x] 
+                    for x in cv_col_names]
             except KeyError:
                 pass
+
+        if self._unwrap:
+            try:
+                image_cols = [attr_to_col_ix[x] for x in ["ix", "iy", "iz"]]
+            except:
+                raise ValueError("Trajectory must have image flag in order "
+                                 "to unwrap.")
+
+        self._has_vels = all(x in attr_to_col_ix for x in ["vx", "vy", "vz"])
+        if self._has_vels:
+            ts.has_velocities = True
+            vel_cols = [attr_to_col_ix[x] for x in ["vx", "vy", "vz"]]
+        self._has_forces = all(x in attr_to_col_ix for x in ["fx", "fy", "fz"])
+        if self._has_forces:
+            ts.has_forces = True
+            force_cols = [attr_to_col_ix[x] for x in ["fx", "fy", "fz"]]
+
         # this should only trigger on first read of "ATOM" card, after which it
         # is fixed to the guessed value. Auto proceeds unscaled -> scaled
         # -> unwrapped -> scaled_unwrapped
@@ -620,22 +668,46 @@ def _read_next_timestep(self):
             except IndexError:
                 raise ValueError("No coordinate information detected")
         elif not self.lammps_coordinate_convention in convention_to_col_ix:
-            raise ValueError(f"No coordinates following convention {self.lammps_coordinate_convention} found in timestep")
+            raise ValueError(f"No coordinates following convention "
+                             "{self.lammps_coordinate_convention} found in "
+                             "timestep")
 
         coord_cols = convention_to_col_ix[self.lammps_coordinate_convention]
+        if self._unwrap:
+            coord_cols.extend(image_cols)
 
         ids = "id" in attr_to_col_ix
         for i in range(self.n_atoms):
             fields = f.readline().split()
             if ids:
                 indices[i] = fields[attr_to_col_ix["id"]]
-            ts.positions[i] = [fields[dim] for dim in coord_cols]
+            coords = np.array([fields[dim] for dim in coord_cols], 
+                              dtype=np.float32)
+
+            if self._unwrap:
+                images = coords[3:]
+                coords = coords[:3]
+                coords += images * ts.dimensions[:3]
+            else:
+                coords = coords[:3]
+            ts.positions[i] = coords
+
+            if self._has_vels:
+                ts.velocities[i] = [fields[dim] for dim in vel_cols]
+            if self._has_forces:
+                ts.forces[i] = [fields[dim] for dim in force_cols]
 
         order = np.argsort(indices)
         ts.positions = ts.positions[order]
+        if self._has_vels:
+            ts.velocities = ts.velocities[order]
+        if self._has_forces:
+            ts.forces = ts.forces[order]
         if (self.lammps_coordinate_convention.startswith("scaled")):
             # if coordinates are given in scaled format, undo that
             ts.positions = distances.transform_StoR(ts.positions,
                                                     ts.dimensions)
+        # Transform to origin after transformation of scaled variables
+        ts.positions -= np.array([xlo, ylo, zlo])[None,:]
 
         return ts

package/MDAnalysis/coordinates/XDR.py

@@ -197,11 +197,14 @@ def _load_offsets(self):
         try:
             with fasteners.InterProcessLock(lock_name) as filelock:
                 pass
-        except PermissionError:
-            warnings.warn(f"Cannot write lock/offset file in same location as "
-                           "{self.filename}. Using slow offset calculation.")
-            self._read_offsets(store=True)
-            return
+        except OSError as e:
+            if isinstance(e, PermissionError) or e.errno == errno.EROFS:
+                warnings.warn(f"Cannot write lock/offset file in same location as "
+                               "{self.filename}. Using slow offset calculation.")
+                self._read_offsets(store=True)
+                return
+            else:
+                raise
 
         with fasteners.InterProcessLock(lock_name) as filelock:
             if not isfile(fname):

package/MDAnalysis/coordinates/memory.py

@@ -501,9 +501,16 @@ def timeseries(self, asel=None, start=0, stop=-1, step=1, order='afc'):
             of the underlying numpy array is returned, while a copy of the
             data is returned whenever `asel` is different from ``None``.
         start : int (optional)
+            the start trajectory frame
         stop : int (optional)
+            the end trajectory frame
+
+            .. deprecated:: 2.4.0
+               Note that `stop` is currently *inclusive* but will be
+               changed in favour of being *exclusive* in version 3.0.  
+
         step : int (optional)
-            range of trajectory to access, `start` and `stop` are *inclusive*
+            the number of trajectory frames to skip
         order : {"afc", "acf", "caf", "fac", "fca", "cfa"} (optional)
             the order/shape of the return data array, corresponding
             to (a)tom, (f)rame, (c)oordinates all six combinations
@@ -515,6 +522,11 @@ def timeseries(self, asel=None, start=0, stop=-1, step=1, order='afc'):
         .. versionchanged:: 1.0.0
            Deprecated `format` keyword has been removed. Use `order` instead.
         """
+        if stop != -1:
+            warnings.warn("MemoryReader.timeseries inclusive `stop` "
+                      "indexing will be removed in 3.0 in favour of exclusive "
+                      "indexing", category=DeprecationWarning)
+
         array = self.get_array()
         if order == self.stored_order:
             pass

package/MDAnalysis/coordinates/timestep.pyx

@@ -107,7 +107,7 @@ MDAnalysis.
 
       .. versionadded:: 0.11.0
          Added as optional to :class:`Timestep`
-  
+
    .. autoattribute:: dtype
    .. autoattribute:: dimensions
    .. autoattribute:: triclinic_dimensions
@@ -127,21 +127,10 @@ MDAnalysis.
 """
 
 
-from ..lib.util import asiterable, Namespace
-from ..auxiliary.core import auxreader
-from ..auxiliary.base import AuxReader
-from .. import units
-from .. import (
-    _READERS, _READER_HINTS,
-    _SINGLEFRAME_WRITERS,
-    _MULTIFRAME_WRITERS,
-    _CONVERTERS
-)
+from ..lib.util import Namespace
 from ..exceptions import NoDataError
 from . import core
-from libcpp cimport bool
-from libc.stdint cimport int64_t, uint64_t
-import cython
+from libc.stdint cimport uint64_t
 import weakref
 import warnings
 import copy
@@ -535,8 +524,6 @@ cdef class Timestep:
         """Set the unitcell for this Timestep as defined by triclinic vectors
         .. versionadded:: 0.11.0
         """
-        cdef int flag = 0
-
         if new_dimensions is None:
             self.dimensions = None
         else:

package/MDAnalysis/lib/_augment.pyx

@@ -27,7 +27,7 @@ import numpy as np
 from .mdamath import triclinic_vectors
 cimport numpy as np
 cimport MDAnalysis.lib._cutil
-from MDAnalysis.lib._cutil cimport _dot ,_norm, _cross
+from MDAnalysis.lib._cutil cimport _dot, _norm, _cross
 
 from libcpp.vector cimport vector