docs: def _real(coords1, charges1, coords2, charges2, rcut, alpha, box):
      """Calculate ewald real 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]
      # Unit vectors
      a = box[0]
      b = box[1]
      c = box[2]
      
      
      # This is helpful to add the correct number of boxes
      l_max = int(np.ceil(2.0 * rcut / np.min(np.trace(box))))
      result = np.zeros(n)
      
      
      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]
  
              for l_i in range(-l_max, l_max + 1):
                  for l_j in range(-l_max, l_max + 1):
                      for l_k in range(-l_max, l_max + 1):
                          nv = l_i * a + l_j * b + l_k * c
                          r_j_n = r_j + nv
                          
                          r_ij = _dist(r_i, r_j_n)
                          if r_ij < 1e-10 or r_ij > rcut:
                              continue
                          value = q_i * q_j * math.erfc(alpha * r_ij) / r_ij
                          result[i] += value
  
      return result