diff --git a/mrdna/readers/segmentmodel_from_pdb.py b/mrdna/readers/segmentmodel_from_pdb.py
index 3cbb683de5a95be8f0fa76ea565e449b7cf12b76..53b0c31698a016032fe617c2a353a76a7f6d690b 100644
--- a/mrdna/readers/segmentmodel_from_pdb.py
+++ b/mrdna/readers/segmentmodel_from_pdb.py
@@ -149,16 +149,20 @@ def find_basepairs( u, centers, transforms, selection_text='all' ):
     possible_resI = np.where( ids )[0]
     possible_resJ = possible_basepairs[ ids ].astype(int) 
     resI,resJ = [[],[]]
-    for i,j,R1,c1,c2 in zip(possible_resI,possible_resJ,
+    for i,j,R1,R2,c1,c2 in zip(possible_resI,possible_resJ,
                             transforms[possible_resI],
+                            transforms[possible_resJ],
                             centers[possible_resI],
                             centers[possible_resJ]):
         c2_expected = c1 + ref_bp_position.dot(R1)
         # fh.write("graphics top cylinder {{{}}} {{{}}} radius 0.2 resolution 30\n".format(printv(c1),printv(c2)))
 
-        if np.linalg.norm(c2_expected-c2) < 2:
-            resI.append(i)
-            resJ.append(j)
+        if np.linalg.norm(c2_expected-c2) < 3.5:
+            angle= R1.T.dot(np.array((0,0,1))).dot(R2.T.dot(np.array((0,0,1))))
+            if angle < -0.7:
+                resI.append(i)
+                resJ.append(j)
+
 
     resI = np.array(resI, dtype=np.int)
     resJ = np.array(resJ, dtype=np.int)