fixed Issue #2001 - support for Tinker periodic boundary box

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:

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)

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):

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]

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)
+

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

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')