docs: def guess_bonds(r_array, type_array, threshold=0.1, maxradius=0.3, radii_dict=None):
      '''Detect bonds given the coordinates (r_array) and types of the 
      atoms involved (type_array), based on their covalent radii.
      
      To fine-tune the detection, it is possible to set a **threshold** and a 
      maximum search radius **maxradius**, and the radii lookup **radii_dict**.
      '''
      if radii_dict is None:
          covalent_radii = cdb.get('data', 'covalentdict')
      else:
          covalent_radii = radii_dict
      
      # Find all the pairs
      ck = cKDTree(r_array)
      pairs = ck.query_pairs(maxradius)
      
      bonds = []
      for i,j in pairs:
          a, b = covalent_radii[type_array[i]], covalent_radii[type_array[j]]
          rval = a + b
          
          thr_a = rval - threshold
          thr_b = rval + threshold 
          
          thr_a2 = thr_a**2
          thr_b2 = thr_b**2
          
          dr2  = ((r_array[i] - r_array[j])**2).sum()
          # print(thr_a, dr2**0.5, thr_b)
          if thr_a2 < dr2 < thr_b2:
              bonds.append((i, j))
      return np.array(bonds)