changing DL_Poly units

package/CHANGELOG

@@ -22,6 +22,7 @@ The rules for this file:
  * 2.10.0
 
 Fixes
+ * Fixed force unit conversion in DL_Poly ConfigReader and HistoryReader to account for conversion to 10 J/mol.A (Issue #2393, PR #4983)
  * Fixed an integer overflow in large DCD file seeks on Windows
    (Issue #4879, PR #5086)
  * Fix compile failure due to numpy issues in transformations.c (Issue #5061, PR #5068)

package/MDAnalysis/coordinates/DLPoly.py

@@ -35,7 +35,8 @@
 from ..lib import util
 from ..lib.util import cached, store_init_arguments
 
-_DLPOLY_UNITS = {'length': 'Angstrom', 'velocity': 'Angstrom/ps', 'time': 'ps'}
+_DLPOLY_UNITS = {"length": "Angstrom", "velocity": "Da*Angstrom/ps", "time": "ps",
+                 'force': "Da*Angstrom/ps"}
 
 
 class ConfigReader(base.SingleFrameReaderBase):
@@ -46,14 +47,17 @@ class ConfigReader(base.SingleFrameReaderBase):
     .. versionchanged:: 2.0.0
        coordinates, velocities, and forces are no longer stored in 'F' memory
        layout, instead now using the numpy default of 'C'.
+    .. versionchanged:: 2.10.0
+        Forces are now converted to 10 J/(mol.Å) as per DL_Poly
     """
-    format = 'CONFIG'
+
+    format = "CONFIG"
     units = _DLPOLY_UNITS
 
     def _read_first_frame(self):
         unitcell = np.zeros((3, 3), dtype=np.float32)
 
-        with open(self.filename, 'r') as inf:
+        with open(self.filename, "r") as inf:
             self.title = inf.readline().strip()
             levcfg, imcon, megatm = np.int64(inf.readline().split()[:3])
             if not imcon == 0:
@@ -99,7 +103,10 @@ def _read_first_frame(self):
         if has_vels:
             velocities = np.array(velocities, dtype=np.float32)
         if has_forces:
+            # DL_POLY forces are in units of 10 J/(mol.Å)
             forces = np.array(forces, dtype=np.float32)
+            if self.convert_units:
+                forces = self.convert_forces_from_native(forces)
         self.n_atoms = len(coords)
 
         if ids:
@@ -113,10 +120,9 @@ def _read_first_frame(self):
             if has_forces:
                 forces = forces[order]
 
-        ts = self.ts = self._Timestep(self.n_atoms,
-                                      velocities=has_vels,
-                                      forces=has_forces,
-                                      **self._ts_kwargs)
+        ts = self.ts = self._Timestep(
+            self.n_atoms, velocities=has_vels, forces=has_forces, **self._ts_kwargs
+        )
         ts._pos = coords
         if has_vels:
             ts._velocities = velocities
@@ -132,8 +138,11 @@ class HistoryReader(base.ReaderBase):
     """Reads DLPoly format HISTORY files
 
     .. versionadded:: 0.11.0
+    .. versionchanged:: 2.10.0
+        Forces are now converted to 10 J/(mol.Å) as per DL_Poly
     """
-    format = 'HISTORY'
+
+    format = "HISTORY"
     units = _DLPOLY_UNITS
 
     @store_init_arguments
@@ -142,7 +151,7 @@ def __init__(self, filename, **kwargs):
         self._cache = {}
 
         # "private" file handle
-        self._file = util.anyopen(self.filename, 'r')
+        self._file = util.anyopen(self.filename, "r")
         self.title = self._file.readline().strip()
         header = np.int64(self._file.readline().split()[:3])
         self._levcfg, self._imcon, self.n_atoms = header
@@ -157,26 +166,28 @@ def __init__(self, filename, **kwargs):
             self._has_cell = False
         self._file.seek(rwnd)
 
-        self.ts = self._Timestep(self.n_atoms,
-                                 velocities=self._has_vels,
-                                 forces=self._has_forces,
-                                 **self._ts_kwargs)
+        self.ts = self._Timestep(
+            self.n_atoms,
+            velocities=self._has_vels,
+            forces=self._has_forces,
+            **self._ts_kwargs
+        )
         self._read_next_timestep()
 
     def _read_next_timestep(self, ts=None):
         if ts is None:
             ts = self.ts
 
         line = self._file.readline()  # timestep line
-        if not line.startswith('timestep'):
+        if not line.startswith("timestep"):
             raise IOError
 
         if self._has_cell:
             unitcell = np.zeros((3, 3))
             unitcell[0] = self._file.readline().split()
             unitcell[1] = self._file.readline().split()
             unitcell[2] = self._file.readline().split()
-            ts.dimensions = core.triclinic_box(*unitcell)            
+            ts.dimensions = core.triclinic_box(*unitcell)
 
         # If ids are given, put them in here
         # and later sort by them
@@ -209,6 +220,9 @@ def _read_next_timestep(self, ts=None):
                     ts._velocities[:] = ts._velocities[order]
                 if self._has_forces:
                     ts._forces[:] = ts._forces[order]
+
+        if self._has_forces and self.convert_units:
+            ts._forces = self.convert_forces_from_native(ts._forces)
 
         ts.frame += 1
         return ts
@@ -220,17 +234,17 @@ def _read_frame(self, frame):
         return self._read_next_timestep()
 
     @property
-    @cached('n_frames')
+    @cached("n_frames")
     def n_frames(self):
         # Second line is traj_key, imcom, n_atoms, n_frames, n_records
         offsets = []
 
-        with open(self.filename, 'r') as f:
+        with open(self.filename, "r") as f:
             f.readline()
             f.readline()
             position = f.tell()
             line = f.readline()
-            while line.startswith('timestep'):
+            while line.startswith("timestep"):
                 offsets.append(position)
                 if self._has_cell:
                     f.readline()
@@ -251,7 +265,7 @@ def n_frames(self):
 
     def _reopen(self):
         self.close()
-        self._file = open(self.filename, 'r')
+        self._file = open(self.filename, "r")
         self._file.readline()  # header is 2 lines
         self._file.readline()
         self.ts.frame = -1

package/MDAnalysis/coordinates/TRZ.py

@@ -605,4 +605,4 @@ def close(self):
         if self.trzfile is None:
             return
         self.trzfile.close()
-        self.trzfile = None
+        self.trzfile = None

package/MDAnalysis/units.py

@@ -347,6 +347,7 @@ def __getitem__(self, key):
     "N": 1e13 / constants["N_Avogadro"],
     "J/m": 1e13 / constants["N_Avogadro"],
     "kcal/(mol*Angstrom)": 1 / constants["calorie"],
+    "Da*Angstrom/ps": 0.01,
 }
 # (TODO: build this combinatorically from lengthUnit and energyUnit)
 

testsuite/MDAnalysisTests/coordinates/test_dlpoly.py

@@ -71,7 +71,7 @@ def test_velocities(self, ts):
         assert_allclose(ts._velocities[0], ref)
 
     def test_forces(self, ts):
-        ref = np.array([-1979.558687, 739.7961625, 1027.996603])
+        ref = np.array([-197955.8687, 73979.61625, 102799.6603])
         assert_allclose(ts._forces[0], ref)
 
 
@@ -119,7 +119,7 @@ def test_vel(self, u):
         assert_allclose(u.atoms[2].velocity, ref)
 
     def test_for(self, u):
-        ref = np.array([150.3309776, -812.6932914, 1429.413120])
+        ref = np.array([15033.09776, -81269.32914, 142941.3120])
         assert_allclose(u.atoms[2].force, ref)
 
     def test_number(self, u):
@@ -189,9 +189,9 @@ def test_velocity(self, u):
     def test_force(self, u):
         ref = np.array(
             [
-                [-2621.386432, 1579.334443, 1041.103241],
-                [-1472.262341, 2450.379615, -8149.916193],
-                [2471.802059, -3828.467296, 3596.679326],
+                [-262138.6432, 157933.4443, 104110.3241],
+                [-147226.2341, 245037.9615, -814991.6193],
+                [247180.2059, -382846.7296, 359667.9326],
             ]
         )
         for ts, r in zip(u.trajectory, ref):