docs: def voropp_for_non_orthorhombic_cells(_a, q='%v', voropp_path=VOROPP_PATH,
                                        fast=False, dump=None):
      """
      Run voro++ on current configuration and return selected quantities.
      Parameter *q* can be a list of voro++ output quantities.
      Run 'voro++ -hc' to see options. Will take care of Lees-Edwards boundary
      conditions by embedding a sheared cell in its periodic images and then
      throwing away the border (will hence be slower).
      """
      
      # Make a copy because we will modify the Atoms object
      a = _a.copy()
      nat = len(a)
      # Wrap into cell
      a.set_scaled_positions(a.get_scaled_positions()%1.0)
  
      # shear_dx should go into the cell
      if 'shear_dx' in a.info:
          lx, ly, lz = a.get_cell().diagonal()
          cx, cy, cz = a.get_cell()
          assert abs(cz[0]) < 1e-12
          assert abs(cz[1]) < 1e-12
          shear_dx = a.info['shear_dx']
          cz[0] = shear_dx[0]
          cz[1] = shear_dx[1]
          a.set_cell([cx, cy, cz], scale_atoms=False)
          a.set_scaled_positions(a.get_scaled_positions()%1.0)
          
      cx, cy, cz = a.get_cell()
     
      if fast:
          # Make 2 copies of the box in each direction. Could fail!
          a *= ( 2, 2, 2 )
  
          # Translate such that the box with the lowest indices sits in the middle
          a.translate(cx/2+cy/2+cz/2)
      else:
          # Make 3 copies of the box in each direction
          a *= ( 3, 3, 3 )
         
          # Translate such that the box with the lowest indices sits in the middle
          a.translate(cx+cy+cz)
          
      # Wrap back to box
      a.set_scaled_positions(a.get_scaled_positions()%1.0)
  
      # Get enclosing box
      lower, upper = get_enclosing_orthorhombic_box(a.get_cell())
      elx, ely, elz = upper-lower
      
      # Shift and set cell such that the general system is enclosed in the
      # orthorhombic box
      a.translate(-lower)
      a.set_cell([elx,ely,elz], scale_atoms=False)
      
      # Dump snapshot for debugging purposes
      if dump:
          write(dump, a)
      
      # Do Voronoi analysis    
      x, y, z = a.get_positions().T
      f = open('tmp.voronoi', 'w')
      for jn, ( jx, jy, jz ) in enumerate(zip(x, y, z)):
          print >> f, jn, jx, jy, jz
      f.close()
      if isinstance(q, str) or isinstance(q, unicode):
          c = '%s' % format(q)
      else:
          c = reduce(lambda x,y: '{0} {1}'.format(x,y), map(lambda x: '%'+x, q))
      os.system('{0} -o -p -c "%i {1}" 0 {2} 0 {3} 0 {4} tmp.voronoi' \
                .format(voropp_path, c, elx, ely, elz))
      r = np.loadtxt('tmp.voronoi.vol', unpack=True)
      os.remove('tmp.voronoi')
      os.remove('tmp.voronoi.vol')
      
      # Sort particles according to their ids
      r = r[:,np.array(r[0,:], dtype=int)]
      # Use only the lowest indices (i.e. the box in the middle)
      if r.shape[0] == 2:
          return r[1,:nat]
      else:
          return r[1:,:nat]