docs: def random_box(molecules, total=None, proportions=None, size=[1.,1.,1.], maxtries=100):
      '''Create a System made of a series of random molecules.
      
      Parameters:
      
      total:
      molecules:
      proportions:
      '''
      
      # Setup proportions to be right
      if proportions is None:
          proportions = np.ones(len(molecules)) / len(molecules)
      else:
          proportions = np.array(proportions)
      
      size = np.array(size)
      
      tree = CoverTree(metric="periodic", metric_args={'cell_lengths': size})
      
      type_array = []
      result = []
      vdw_radii = []
      max_vdw = max(vdw_radius(np.concatenate([m.type_array for m in molecules])))
      
      first = True
      for l, n in enumerate((proportions * total).astype(int)):
          
          # We try to insert each molecule    
          for i in range(n):
              
              # Attempts
              for k in range(maxtries):
                  template = molecules[l].copy()
                  reference = np.random.uniform(0, 1, 3) * size
                  r_array = template.r_array + reference
  
                  # Find all collision candidates
                  pts_list, distances_list = tree.query_ball_many(r_array, vdw_radius(template.type_array) + max_vdw)
                  # print pts_list, distances_list
                  # Check if there is any collision
                  ok = True
                  for i, (dist, pts) in enumerate(zip(distances_list, pts_list)):
                      if len(dist) == 0:
                          break
  
                      found_vdw = np.array([vdw_radii[p] for p in pts])
                      ok &= all(dist > found_vdw + vdw_radius(template.type_array[i]))
  
                  if ok:
                      tree.insert_many(r_array)
                      template.r_array = r_array
                      result.append(template)
                      vdw_radii.extend(vdw_radius(template.type_array))
                      break
              if not ok:
                  raise Exception("Trials exceeded")
      
      system = System(result)
      system.box_vectors[0, 0] = size[0]
      system.box_vectors[1, 1] = size[1]
      system.box_vectors[2, 2] = size[2]
      return system