docs: def parse_gro_lines(lines):
      '''Reusable parsing'''
      title = lines.pop(0)
      natoms = int(lines.pop(0))
      atomlist = []
  
      # I need r_array, type_array
      datalist = []
      for l in lines:
          fields = l.split()
          line_length = len(l)
  
          if line_length in (45, 46, 69, 70):
              #Only positions are provided
              molidx = int(l[0:5])
              moltyp = l[5:10].strip()
              attyp = l[10:15].strip()
              atidx = int(l[15:20])
              rx = float(l[20:28])
              ry = float(l[28:36])
              rz = float(l[36:44])
  
              hasvel = False
              if line_length == 69:
                  hasvel = True
                  # Provide velocities
                  vx = float(l[44:52])
                  vy = float(l[52:60])
                  vz = float(l[60:68])
  
              # Do I have to convert back the atom types, probably yes???
              #if attyp.lower() not in symbol_list:
              #    attyp = gro_to_cl[attyp]
              datalist.append((molidx, moltyp, attyp, rx, ry, rz))
          else:
              # This is the box size
              stuff  = [float(f) for f in fields]
              a,b,c = stuff[0], stuff[1], stuff[2]
              box_vectors = np.array([[a, 0, 0], [0, b, 0], [0, 0, c]])
              break
  
      dataarr = np.array(datalist,
                         dtype=np.dtype([('f0', int), ('f1', object),
                                         ('f2', object), ('f3', np.float64),
                                         ('f4', np.float64), ('f5', np.float64)]
                                        )
                         )
  
      # Molecule indices: unique elements in molidx
      mol_id, mol_indices = np.unique(dataarr['f0'], return_index=True)
      
      maps = {('atom', 'molecule') : dataarr['f0'] - 1}
      
      r_array = np.vstack([dataarr['f3'],
                           dataarr['f4'],
                           dataarr['f5']]).transpose()
      grotype_array = dataarr['f2']
      grores_array = dataarr['f1'][mol_indices]
      
      molecule_export = np.array([dict(groname=g)
                             for g in dataarr['f1'][mol_indices]])
      atom_export = np.array([dict(grotype=g) for g in grotype_array])
  
      # Gromacs Defaults to Unknown Atom type
      
      # We need to parse the gromacs type in some way...
      # grotype_array = [re.sub('[0-9+-]+$', '', g) for g in grotype_array]
      type_array = np.array([gro_to_cl.get(re.sub('[0-9+-]+$','', g), "Unknown") for g in grotype_array])
  
      # Molecular Formula Arrays
      mol_formula = []
      end = len(r_array)
      for i, _ in enumerate(mol_indices):
          s = mol_indices[i]
          e = mol_indices[i+1] if i+1 < len(mol_indices) else end
          mol_formula.append(make_formula(type_array[s:e]))
  
      mol_formula = np.array(mol_formula)
  
      # n_mol, n_at
      sys = System.from_arrays(r_array=r_array,
                               maps=maps,
                               type_array=type_array,
                               atom_name=grotype_array,
                               atom_export=atom_export,
                               molecule_export=molecule_export,
                               molecule_name=grores_array,
                               box_vectors=box_vectors)
      return sys