Add pbc kwarg to Contacts

package/CHANGELOG

@@ -73,6 +73,7 @@ Fixes
   * Fixed HydrogenBondAnalysis to return atom indices rather than atom ids (PR #2572).
     Fixed the check for bond information in the _get_dh_pairs method (Issue #2396).
   * Added missing selection module to leaflet.py (Issue #2612)
+  * Contact Analysis class respects PBC (Issue #2368)
 
 Enhancements
   * Added new density.DensityAnalysis (Issue #2502)

package/MDAnalysis/analysis/contacts.py

@@ -373,10 +373,12 @@ class Contacts(AnalysisBase):
     .. versionchanged:: 1.0.0
        ``save()`` method has been removed. Use ``np.savetxt()`` on
        :attr:`Contacts.timeseries` instead.
+    .. versionchanged:: 2.0.0
+        added ``pbc`` attribute to calculate distances using PBC.
 
     """
     def __init__(self, u, select, refgroup, method="hard_cut", radius=4.5,
-                 kwargs=None, **basekwargs):
+                 pbc=True, kwargs=None, **basekwargs):
         """
         Parameters
         ----------
@@ -392,6 +394,9 @@ def __init__(self, u, select, refgroup, method="hard_cut", radius=4.5,
         method : string | callable (optional)
             Can either be one of ``['hard_cut' , 'soft_cut', 'radius_cut']`` or a callable
             with call signature ``func(r, r0, **kwargs)`` (the "Contacts API").
+        pbc : bool (optional)
+            Uses periodic boundary conditions to calculate distances if set to ``True``; the
+            default is ``True``.
         kwargs : dict, optional
             dictionary of additional kwargs passed to `method`. Check
             respective functions for reasonable values.
@@ -424,30 +429,40 @@ def __init__(self, u, select, refgroup, method="hard_cut", radius=4.5,
         self.select = select
         self.grA = u.select_atoms(select[0])
         self.grB = u.select_atoms(select[1])
-
+        self.pbc = pbc
+
         # contacts formed in reference
         self.r0 = []
         self.initial_contacts = []
 
+        #get dimension of box if pbc set to True
+        if self.pbc:
+            self._get_box = lambda ts: ts.dimensions
+        else:
+            self._get_box = lambda ts: None
+
         if isinstance(refgroup[0], AtomGroup):
             refA, refB = refgroup
-            self.r0.append(distance_array(refA.positions, refB.positions))
+            self.r0.append(distance_array(refA.positions, refB.positions,
+                                            box=self._get_box(refA.universe)))
             self.initial_contacts.append(contact_matrix(self.r0[-1], radius))
+
         else:
             for refA, refB in refgroup:
-                self.r0.append(distance_array(refA.positions, refB.positions))
-                self.initial_contacts.append(contact_matrix(self.r0[-1],
-                                                            radius))
+                self.r0.append(distance_array(refA.positions, refB.positions,
+                                                box=self._get_box(refA.universe)))
+                self.initial_contacts.append(contact_matrix(self.r0[-1], radius))
 
     def _prepare(self):
         self.timeseries = np.empty((self.n_frames, len(self.r0)+1))
 
     def _single_frame(self):
         self.timeseries[self._frame_index][0] = self._ts.frame
-
+        
         # compute distance array for a frame
-        d = distance_array(self.grA.positions, self.grB.positions)
-
+        d = distance_array(self.grA.positions, self.grB.positions,
+                            box=self._get_box(self._ts))
+
         for i, (initial_contacts, r0) in enumerate(zip(self.initial_contacts,
                                                        self.r0), 1):
             # select only the contacts that were formed in the reference state

testsuite/MDAnalysisTests/analysis/test_contacts.py

@@ -37,6 +37,8 @@
 from MDAnalysisTests.datafiles import (
     PSF,
     DCD,
+    TPR,
+    XTC,
     contacts_villin_folded,
     contacts_villin_unfolded,
     contacts_file
@@ -217,7 +219,6 @@ def test_villin_folded(self):
         results = soft_cut(f, u, sel, sel)
         assert_almost_equal(q.timeseries[:, 1], results[:, 1])
 
-
     def test_villin_unfolded(self):
         # both folded
         f = mda.Universe(contacts_villin_folded)
@@ -302,6 +303,23 @@ def test_non_callable_method(self, universe):
         with pytest.raises(ValueError):
             self._run_Contacts(universe, method=2, stop=2)
 
+    @pytest.mark.parametrize("pbc,expected", [
+    (True, [1., 0.43138152, 0.3989021, 0.43824337, 0.41948765,
+            0.42223239, 0.41354071, 0.43641354, 0.41216834, 0.38334858]),
+    (False, [1., 0.42327791, 0.39192399, 0.40950119, 0.40902613,
+             0.42470309, 0.41140143, 0.42897862, 0.41472684, 0.38574822])
+    ])
+    def test_distance_box(self, pbc, expected):
+        u = mda.Universe(TPR, XTC)
+        sel_basic = "(resname ARG LYS)"
+        sel_acidic = "(resname ASP GLU)"
+        acidic = u.select_atoms(sel_acidic)
+        basic = u.select_atoms(sel_basic)
+
+        r = contacts.Contacts(u, select=(sel_acidic, sel_basic),
+                        refgroup=(acidic, basic), radius=6.0, pbc=pbc)
+        r.run()
+        assert_array_almost_equal(r.timeseries[:, 1], expected)
 
 def test_q1q2():
     u = mda.Universe(PSF, DCD)