diff --git a/mdevaluate/distribution.py b/mdevaluate/distribution.py
index 5f7a27c..5f54a96 100644
--- a/mdevaluate/distribution.py
+++ b/mdevaluate/distribution.py
@@ -360,10 +360,29 @@ def next_neighbor_distribution(atoms, reference=None, number_of_neighbors=4, bin
     
     
 def hbonds(D, H, A, box, DA_lim=0.35, HA_lim=0.35, min_cos=np.cos(30*np.pi/180), full_output=False):
+    """
+    Compute h-bond pairs
+    
+    Args:
+        D: Set of coordinates for donators.
+        H: Set of coordinates for hydrogen atoms. Should have the same 
+            length as D.
+        A:  Set of coordinates for acceptors.
+        DA_lim (opt.): Minimum distance beteen donator and acceptor.
+        HA_lim (opt.): Minimum distance beteen hydrogen and acceptor.
+        min_cos (opt.): Minimum cosine for the HDA angle. Default is 
+            equivalent to a maximum angle of 30 degree.
+        full_output (opt.): Returns additionally the cosine of the 
+        angles and the DA distances
+        
+    Return:
+        List of (D,A)-pairs in hbonds. 
+    """ 
+            
     def dist_DltA(D, H, A, box, max_dist=0.35):
         ppoints, pind = pbc_points(D, box, thickness=max_dist+0.1, index=True)
-        Dtree = cKDTree(ppoints)
-        Atree = cKDTree(A)
+        Dtree = spatial.cKDTree(ppoints)
+        Atree = spatial.cKDTree(A)
         pairs = Dtree.sparse_distance_matrix(Atree, max_dist, output_type='ndarray')
         pairs = np.asarray(pairs.tolist())
         pairs = np.int_(pairs[pairs[:,2] > 0][:,:2])
@@ -372,8 +391,8 @@ def hbonds(D, H, A, box, DA_lim=0.35, HA_lim=0.35, min_cos=np.cos(30*np.pi/180),
 
     def dist_AltD(D, H, A, box, max_dist=0.35):
         ppoints, pind = pbc_points(A, box, thickness=max_dist+0.1, index=True)
-        Atree = cKDTree(ppoints)
-        Dtree = cKDTree(D)
+        Atree = spatial.cKDTree(ppoints)
+        Dtree = spatial.cKDTree(D)
         pairs = Atree.sparse_distance_matrix(Dtree, max_dist, output_type='ndarray')   
         pairs = np.asarray(pairs.tolist())
         pairs = np.int_(pairs[pairs[:,2] > 0][:,:2])
@@ -386,9 +405,9 @@ def hbonds(D, H, A, box, DA_lim=0.35, HA_lim=0.35, min_cos=np.cos(30*np.pi/180),
     else:
         pairs = dist_AltD(D,H,A,box,DA_lim)
     
-    vDH = md.pbc.pbc_diff(D[pairs[:,0]], H[pairs[:,0]], box)
-    vDA = md.pbc.pbc_diff(D[pairs[:,0]], A[pairs[:,1]], box)
-    vHA = md.pbc.pbc_diff(H[pairs[:,0]], A[pairs[:,1]], box)
+    vDH = pbc_diff(D[pairs[:,0]], H[pairs[:,0]], box)
+    vDA = pbc_diff(D[pairs[:,0]], A[pairs[:,1]], box)
+    vHA = pbc_diff(H[pairs[:,0]], A[pairs[:,1]], box)
     angles_cos = np.clip(np.einsum('ij,ij->i', vDH, vDA)/
           np.linalg.norm(vDH,axis=1)/
           np.linalg.norm(vDA,axis=1), -1, 1)