docs: def read_atoms(fn, cycfn=None, pos_only=False, conv=1.0):
      """
      Read atom information from an atoms.dat file (i.e., tblmd, MDCORE input file)
      """
      f = paropen(fn, "r")
  
      l = f.readline().lstrip()
      while len(l) > 0 and ( l[0] == '#' or l[0] == '<' ):
          l = f.readline().lstrip()
  
      n_atoms = int(l)
  
      l = f.readline().lstrip()
      while len(l) > 0 and ( l[0] == '#' or l[0] == '<' ):
          l = f.readline().lstrip()
  
      l = f.readline().lstrip()
      while len(l) > 0 and ( l[0] == '#' or l[0] == '<' ):
          l = f.readline().lstrip()
  
      #
      # Read positions
      #
  
      forces = np.zeros( [ n_atoms, 3 ] )
      groups = np.zeros( [ n_atoms ] )
      gamma  = np.zeros( [ n_atoms ] )
      T      = np.zeros( [ n_atoms ] )
  
      ats = [ ]
      for i in range(n_atoms):
          s = l.split()
          #      type   x            y            z
          sym = None
          try:
              Z    = int(s[0])
              sym  = ase.data.chemical_symbols[Z]
          except:
              sym  = s[0]
          a = ase.Atom(sym, ( float(s[2])*conv, float(s[3])*conv, float(s[4])*conv ) )
          groups[i] = int(s[5])
          gamma[i]  = float(s[6])
          T[i]      = float(s[7])
  
          ats += [ a ]
          l = f.readline()
      this = ase.Atoms(ats, pbc=True)
  
      if not pos_only:
          while l and l == "":
              l = f.readline().strip()
  
          while l:
              key = l.strip(" <-#\r\n")
  
              if key.upper() == "VELOCITIES":
                  for i in range(n_atoms):
                      s = f.readline().split()
                      m = this[i].mass
                      if m is None:
                          m = ase.data.atomic_masses[ase.data.chemical_symbols.index(this[i].symbol)]
                      this[i].momentum = ( m*float(s[0]), m*float(s[1]), m*float(s[2]) )
  
                  l = None
  
              elif key.upper() == "FORCES":
                  for i in range(n_atoms):
                      s = f.readline().split()
                      forces[i] =  np.array( [ float(s[0]), float(s[1]), float(s[2]) ] )
  
                  l = None
  
              elif key.upper() == "CHARGES":
                  for i in this:
                      l = f.readline()
                      if l and len(l.split()) == 1:
                          i.charge = float(l)
  
                  l = None
  
              elif key.upper() == "CELL" or key.upper().split()[0:2] == ("BOX", "VECTORS" ):
                  l1 = f.readline()
                  l2 = f.readline()
                  l3 = f.readline()
  
                  this.set_cell( [ [float(x) for x in l1.split()],
                                   [float(x) for x in l2.split()],
                                   [float(x) for x in l3.split()] ] )
  
                  l = None
  
              else:
                  aux = [ ]
                  l = f.readline().strip()
                  while l and l[0] not in [ '<', '#' ]:
                      s = l.split()
  
                      aux += [ [float(x) for x in s] ]
  
                      l = f.readline().strip()
  
                  if len(aux) == n_atoms:
                      this.set_array(key, np.asarray(aux))
                  else:
                      print("Warning: Encountered field '%s' which does not seem to be per-atom data." % key)
  
              if l is None:
                  l = f.readline().strip()
              while l and l == "":
                  l = f.readline().strip()
  
  
      f.close()
  
      this.set_array("forces", forces)
      this.set_array("groups", groups)
      this.set_array("gamma", gamma)
      this.set_array("T", T)
  
      if cycfn:
          read_cyc(this, cycfn, conv=conv)
  
      return this