diff --git a/package/CHANGELOG b/package/CHANGELOG index 960ec3a82c2..3ececd834b5 100644 --- a/package/CHANGELOG +++ b/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) diff --git a/package/MDAnalysis/coordinates/DLPoly.py b/package/MDAnalysis/coordinates/DLPoly.py index 71297f978fa..b4b8ba1a823 100644 --- a/package/MDAnalysis/coordinates/DLPoly.py +++ b/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,10 +166,12 @@ 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): @@ -168,7 +179,7 @@ def _read_next_timestep(self, ts=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: @@ -176,7 +187,7 @@ def _read_next_timestep(self, ts=None): 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 diff --git a/package/MDAnalysis/coordinates/TRZ.py b/package/MDAnalysis/coordinates/TRZ.py index 82ada24e6a3..b7c8086103c 100644 --- a/package/MDAnalysis/coordinates/TRZ.py +++ b/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 \ No newline at end of file diff --git a/package/MDAnalysis/units.py b/package/MDAnalysis/units.py index 4f186521190..8cb09c187c9 100644 --- a/package/MDAnalysis/units.py +++ b/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) diff --git a/testsuite/MDAnalysisTests/coordinates/test_dlpoly.py b/testsuite/MDAnalysisTests/coordinates/test_dlpoly.py index 63d0a09464c..f0b25420d54 100644 --- a/testsuite/MDAnalysisTests/coordinates/test_dlpoly.py +++ b/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):