Add pbc kwarg to Contacts

package/MDAnalysis/analysis/contacts.py

@@ -376,7 +376,7 @@ class Contacts(AnalysisBase):
 
     """
     def __init__(self, u, select, refgroup, method="hard_cut", radius=4.5,
-                 kwargs=None, **basekwargs):
+                 pbc=False, kwargs=None, **basekwargs):
         """
         Parameters
         ----------
@@ -392,6 +392,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 ``False``.
         kwargs : dict, optional
             dictionary of additional kwargs passed to `method`. Check
             respective functions for reasonable values.
@@ -424,18 +427,29 @@ 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 = []
 
         if isinstance(refgroup[0], AtomGroup):
             refA, refB = refgroup
-            self.r0.append(distance_array(refA.positions, refB.positions))
+            if(self.pbc):
+                self.r0.append(distance_array(refA.positions, refB.positions,
+                                                box=refA.universe.dimensions))
+            else:
+                self.r0.append(distance_array(refA.positions, refB.positions))
+
             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))
+                if(self.pbc):
+                    self.r0.append(distance_array(refA.positions, refB.positions,
+                                                    box=refA.universe.dimensions))
+                else:
+                    self.r0.append(distance_array(refA.positions, refB.positions))
+
                 self.initial_contacts.append(contact_matrix(self.r0[-1],
                                                             radius))
 
@@ -444,9 +458,13 @@ def _prepare(self):
 
     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)
+        if(self.pbc):
+            d = distance_array(self.grA.positions, self.grB.positions,
+                                                    box=self._ts.dimensions)
+        else:
+            d = distance_array(self.grA.positions, self.grB.positions)
 
         for i, (initial_contacts, r0) in enumerate(zip(self.initial_contacts,
                                                        self.r0), 1):
@@ -461,7 +479,8 @@ def _new_selections(u_orig, selections, frame):
     """create stand alone AGs from selections at frame"""
     u = MDAnalysis.Universe(u_orig.filename, u_orig.trajectory.filename)
     u.trajectory[frame]
-    return [u.select_atoms(s) for s in selections]
+    temp_u = u.copy()
+    return [temp_u.select_atoms(s) for s in selections]
 
 
 def q1q2(u, select='all', radius=4.5):

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,26 @@ def test_non_callable_method(self, universe):
         with pytest.raises(ValueError):
             self._run_Contacts(universe, method=2, stop=2)
 
+    def test_distance_box(self):
+        u = mda.Universe(TPR, XTC)
+        sel = "(resname ASP GLU)"
+        acidic = u.select_atoms(sel)
+
+        r = contacts.Contacts(u,
+                              select=(sel, sel),
+                              refgroup=(acidic, acidic), pbc=False)
+
+        q = contacts.Contacts(u,
+                              select=(sel, sel),
+                              refgroup=(acidic, acidic),
+                              pbc=True)
+
+        r.run()
+        q.run()
+        average_contacts_r = np.mean(r.timeseries[:, 1])
+        average_contacts_q = np.mean(q.timeseries[:, 1])
+
+        assert not average_contacts_q == average_contacts_r
 
 def test_q1q2():
     u = mda.Universe(PSF, DCD)