docs: def merge_systems(sysa, sysb, bounding=0.2):
      '''Generate a system by merging *sysa* and *sysb*.
  
      Overlapping molecules are removed by cutting the molecules of
      *sysa* that have atoms near the atoms of *sysb*. The cutoff distance
      is defined by the *bounding* parameter.
  
      **Parameters**
  
      sysa: System
         First system
      sysb: System
         Second system
      bounding: float or False
         Extra space used when cutting molecules in *sysa* to make space
         for *sysb*. If it is False, no overlap handling will be performed.
  
      '''
  
      if bounding is not False:
          # Delete overlaps.
          if sysa.box_vectors is not None:
              periodicity = sysa.box_vectors.diagonal()
          else:
              periodicity = False
  
          p = overlapping_points(sysb.r_array, sysa.r_array,
                                 cutoff=bounding, periodic=periodicity)
  
          sel = np.ones(len(sysa.r_array), dtype=np.bool)
          sel[p] = False
  
          # Rebuild sysa without water molecules
          sysa = subsystem_from_atoms(sysa, sel)
      
      sysres = System.empty(sysa.n_mol + sysb.n_mol, sysa.n_atoms + sysb.n_atoms)
      
      # Assign the attributes
      for attr in type(sysa).attributes:
          attr.assign(sysres,
                      attr.concatenate(sysa, sysb))
      
      # edit the mol_indices and n_mol
      offset = sysa.mol_indices[-1] + sysa.mol_n_atoms[-1]
      sysres.mol_indices[0:sysa.n_mol] = sysa.mol_indices.copy()
      sysres.mol_indices[sysa.n_mol:] = sysb.mol_indices.copy() + offset
      sysres.mol_n_atoms = np.concatenate([sysa.mol_n_atoms, sysb.mol_n_atoms])
      
      sysres.box_vectors = sysa.box_vectors
      
      return sysres