diff --git a/package/.pylintrc b/package/.pylintrc index ab966147016..da780ce8ac3 100644 --- a/package/.pylintrc +++ b/package/.pylintrc @@ -200,7 +200,6 @@ enable=abstract-class-instantiated, xrange-builtin, zip-builtin-not-iterating, map-builtin-not-iterating, - range-builtin-not-iterating, # Things we'd like to try. # Procedure: diff --git a/package/CHANGELOG b/package/CHANGELOG index 05e9b078d1c..3742ad94e2f 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -14,7 +14,7 @@ The rules for this file: ------------------------------------------------------------------------------ ??/??/18 tylerjereddy, richardjgowers, palnabarun, orbeckst, kain88-de, zemanj, - VOD555, davidercruz, jbarnoud, ayushsuhane, hfmull + VOD555, davidercruz, jbarnoud, ayushsuhane, hfmull, micaela-matta * 0.18.1 @@ -68,6 +68,7 @@ Fixes * PDBWriter now properly sets start value * Zero length TopologyGroup now return proper shape array via to_indices (Issue #1974) + * Added periodic boundary box support to the Tinker file reader (Issue #2001) Changes * TopologyAttrs are now statically typed (Issue #1876) diff --git a/package/MDAnalysis/coordinates/TXYZ.py b/package/MDAnalysis/coordinates/TXYZ.py index 260d11d6ead..aa604b96ad0 100644 --- a/package/MDAnalysis/coordinates/TXYZ.py +++ b/package/MDAnalysis/coordinates/TXYZ.py @@ -74,12 +74,19 @@ def __init__(self, filename, **kwargs): root, ext = os.path.splitext(self.filename) self.xyzfile = util.anyopen(self.filename) self._cache = dict() - + # Check if file has box information saved + with util.openany(self.filename) as inp: + inp.readline() + line = inp.readline() + # If second line has float at second position, we have box info + try: + float(line.split()[1]) + except ValueError: + self.periodic = False + else: + self.periodic = True self.ts = self._Timestep(self.n_atoms, **self._ts_kwargs) - # Haven't quite figured out where to start with all the self._reopen() - # etc. - # (Also cannot just use seek() or reset() because that would break - # with urllib2.urlopen() streams) + self._read_next_timestep() @property @@ -103,6 +110,8 @@ def _read_xyz_n_frames(self): # the number of lines in the XYZ file will be 1 greater than the # number of atoms linesPerFrame = self.n_atoms + 1 + if self.periodic: + linesPerFrame += 1 counter = 0 offsets = [] @@ -134,6 +143,8 @@ def _read_next_timestep(self, ts=None): try: # we assume that there is only one header line per frame f.readline() + if self.periodic: + ts.dimensions = f.readline().split() # convert all entries at the end once for optimal speed tmp_buf = [] for i in range(self.n_atoms): diff --git a/package/MDAnalysis/topology/TXYZParser.py b/package/MDAnalysis/topology/TXYZParser.py index 63698602d06..74bfdce4513 100644 --- a/package/MDAnalysis/topology/TXYZParser.py +++ b/package/MDAnalysis/topology/TXYZParser.py @@ -43,7 +43,10 @@ """ from __future__ import absolute_import + +import itertools import numpy as np +from six.moves import zip from . import guessers from ..lib.util import openany @@ -88,9 +91,22 @@ def parse(self, **kwargs): names = np.zeros(natoms, dtype=object) types = np.zeros(natoms, dtype=np.int) bonds = [] + # Find first atom line, maybe there's box information + fline = inf.readline() + try: + # If a box second value will be a float + # If an atom, second value will be a string + float(fline.split()[1]) + except ValueError: + # If float conversion failed, we have first atom line + pass + else: + # If previous try succeeded it was a box + # so read another line to find the first atom line + fline = inf.readline() # Can't infinitely read as XYZ files can be multiframe - for i in range(natoms): - line = inf.readline().split() + for i, line in zip(range(natoms), itertools.chain([fline], inf)): + line = line.split() atomids[i]= line[0] names[i] = line[1] types[i] = line[5] diff --git a/testsuite/MDAnalysisTests/coordinates/test_txyz.py b/testsuite/MDAnalysisTests/coordinates/test_txyz.py index 6334352645d..8dd49ac8d7a 100644 --- a/testsuite/MDAnalysisTests/coordinates/test_txyz.py +++ b/testsuite/MDAnalysisTests/coordinates/test_txyz.py @@ -23,8 +23,9 @@ import pytest from numpy.testing import assert_almost_equal +from six.moves import zip -from MDAnalysisTests.datafiles import TXYZ, ARC +from MDAnalysisTests.datafiles import TXYZ, ARC, ARC_PBC import MDAnalysis as mda @@ -39,6 +40,11 @@ def ARC_U(): return mda.Universe(ARC) +@pytest.fixture +def ARC_PBC_U(): + return mda.Universe(ARC_PBC) + + def test_txyz_positions(TXYZ_U): assert_almost_equal(TXYZ_U.atoms.positions[0], [-6.553398, -1.854369, 0.000000]) @@ -56,5 +62,19 @@ def test_arc_positions_frame_2(ARC_U): [-0.231579, -0.350841, -0.037475]) -def test_arc_traj_legnth(ARC_U): +def test_arc_traj_length(ARC_U): assert len(ARC_U.trajectory) == 2 + + +def test_arcpbc_traj_length(ARC_PBC_U): + assert len(ARC_PBC_U.trajectory) == 3 + + +def test_pbc_boxsize(ARC_PBC_U): + ref_dimensions=[[ 38.860761, 38.860761, 38.860761, 90.000000, 90.000000, 90.000000], + [ 39.860761, 39.860761, 39.860761, 90.000000, 90.000000, 90.000000], + [ 40.860761, 40.860761, 40.860761, 90.000000, 90.000000, 90.000000]] + + for ref_box, ts in zip(ref_dimensions, ARC_PBC_U.trajectory): + assert_almost_equal(ref_box, ts.dimensions, decimal=5) + diff --git a/testsuite/MDAnalysisTests/data/coordinates/new_hexane.arc b/testsuite/MDAnalysisTests/data/coordinates/new_hexane.arc new file mode 100644 index 00000000000..6705ae636a9 --- /dev/null +++ b/testsuite/MDAnalysisTests/data/coordinates/new_hexane.arc @@ -0,0 +1,24 @@ + 6 Hexane (267 OPLS-UA in 38.8 Ang Box; JPC, 100, 14511 '96) + 38.860761 38.860761 38.860761 90.000000 90.000000 90.000000 + 1 CH3 -0.533809 8.893843 -17.718901 83 2 + 2 CH2 -0.492066 8.734009 -16.201869 86 1 3 + 3 CH2 0.609320 7.723360 -15.894925 86 2 4 + 4 CH2 0.723163 7.301395 -14.432851 86 3 5 + 5 CH2 1.923300 6.401842 -14.151512 86 4 6 + 6 CH3 1.576645 4.928146 -14.343148 83 5 + 6 Hexane (test file for MDA) + 39.860761 39.860761 39.860761 90.000000 90.000000 90.000000 + 1 CH3 0.039519 9.187867 -17.899888 83 2 + 2 CH2 -0.171280 8.486514 -16.561104 86 1 3 + 3 CH2 0.988507 7.632303 -16.057222 86 2 4 + 4 CH2 0.630146 7.086837 -14.677831 86 3 5 + 5 CH2 1.729866 6.240985 -14.042356 86 4 6 + 6 CH3 2.065117 4.899261 -14.687383 83 5 + 6 Hexane + 40.860761 40.860761 40.860761 90.000000 90.000000 90.000000 + 1 CH3 -0.533809 8.893843 -17.718901 83 2 + 2 CH2 -0.492066 8.734009 -16.201869 86 1 3 + 3 CH2 0.609320 7.723360 -15.894925 86 2 4 + 4 CH2 0.723163 7.301395 -14.432851 86 3 5 + 5 CH2 1.923300 6.401842 -14.151512 86 4 6 + 6 CH3 1.576645 4.928146 -14.343148 83 5 diff --git a/testsuite/MDAnalysisTests/datafiles.py b/testsuite/MDAnalysisTests/datafiles.py index 3109f100527..afc51555d0e 100644 --- a/testsuite/MDAnalysisTests/datafiles.py +++ b/testsuite/MDAnalysisTests/datafiles.py @@ -85,7 +85,7 @@ "XTC_sub_sol", "XYZ", "XYZ_psf", "XYZ_bz2", "XYZ_mini", "XYZ_five", # 3 and 5 atoms xyzs for an easy topology - "TXYZ", "ARC", # Tinker files + "TXYZ", "ARC", "ARC_PBC", # Tinker files "PRM", "TRJ", "TRJ_bz2", # Amber (no periodic box) "INPCRD", "PRMpbc", "TRJpbc_bz2", # Amber (periodic box) @@ -294,6 +294,7 @@ XYZ_five = resource_filename(__name__, 'data/five.xyz') TXYZ = resource_filename(__name__, 'data/coordinates/test.txyz') ARC = resource_filename(__name__, 'data/coordinates/test.arc') +ARC_PBC = resource_filename(__name__, 'data/coordinates/new_hexane.arc') PRM = resource_filename(__name__, 'data/Amber/ache.prmtop') TRJ = resource_filename(__name__, 'data/Amber/ache.mdcrd')