docs: def _reciprocal(coords1, charges1, coords2, charges2, kmax, kappa, box):
      """Calculate ewald reciprocal part. Box has to be a cuboidal box you should
      transform any other box shape to a cuboidal box before using this.
      
      """
      n = coords1.shape[0]
      m = coords2.shape[0]
      result = np.zeros(n, dtype=np.float64)
      need_self = np.zeros(n, dtype=np.uint8)
  
      # Reciprocal unit vectors
      g1, g2, g3 = reciprocal_vectors(box)
      
      V = box_volume(box)
      
      prefac = 1.0 / (np.pi * V)
      for i in range(n):
          q_i = charges1[i]
          r_i = coords1[i]
          
          for j in range(m):
              q_j = charges2[j]
              r_j = coords2[j]
              
              r_ij = _dist(r_i, r_j)
              if r_ij < 1e-10:
                  need_self[i] = 1
          
              for k_i in range(-kmax, kmax + 1):
                  for k_j in range(-kmax, kmax + 1):
                      for k_k in range(-kmax, kmax + 1):
                          if k_i == 0 and k_j == 0 and k_k == 0:
                              continue
                          # Reciprocal vector
                          k = k_i * g1 + k_j * g2 + k_k * g3
                          
                          k_sq = sqsum(k)
                          
                          result[i] += (prefac * q_i * q_j *
                                        4.0 * np.pi ** 2 / k_sq *
                                        math.exp(-k_sq / (4.0 * kappa ** 2)) *
                                        math.cos(np.dot(k, r_i - r_j)))
          
      # Self-energy correction
      # I had to do some FUCKED UP stuff because NUMBA SUCKS BALLS and 
      # breaks compatibility with simple expressions such as that one
      # apparently doing result -= something is too hard in numba
      self_energy = 2 * (need_self * kappa * charges1 ** 2) / (np.pi**0.5)
      return result - self_energy