docs: def genrecip(a1, a2, a3, encut):
      """
      Args:
          a1, a2, a3: lattice vectors in bohr
          encut: energy cut off in eV
      Returns:
          reciprocal lattice vectors with energy less than encut
      """
      vol = np.dot(a1, np.cross(a2, a3))  # 1/bohr^3
      b1 = (2 * np.pi / vol) * np.cross(a2, a3)  # units 1/bohr
      b2 = (2 * np.pi / vol) * np.cross(a3, a1)
      b3 = (2 * np.pi / vol) * np.cross(a1, a2)
  
      # create list of recip space vectors that satisfy |i*b1+j*b2+k*b3|<=encut
      G_cut = eV_to_k(encut)
      # Figure out max in all recipricol lattice directions
      i_max = int(math.ceil(G_cut / norm(b1)))
      j_max = int(math.ceil(G_cut / norm(b2)))
      k_max = int(math.ceil(G_cut / norm(b3)))
  
      # Build index list
      i = np.arange(-i_max, i_max)
      j = np.arange(-j_max, j_max)
      k = np.arange(-k_max, k_max)
  
      # Convert index to vectors using meshgrid
      indicies = np.array(np.meshgrid(i, j, k)).T.reshape(-1, 3)
      # Multiply integer vectors to get recipricol space vectors
      vecs = np.dot(indicies, [b1, b2, b3])
      # Calculate radii of all vectors
      radii = np.sqrt(np.einsum('ij,ij->i', vecs, vecs))
  
      # Yield based on radii
      for vec, r in zip(vecs, radii):
          if r < G_cut and r != 0:
              yield vec