code: def _calc_real_and_point(self):
          
          fcoords = self._s.frac_coords
          forcepf = 2.0 * self._sqrt_eta / sqrt(pi)
          coords = self._coords
          numsites = self._s.num_sites
          ereal = np.empty((numsites, numsites), dtype=np.float)
  
          forces = np.zeros((numsites, 3), dtype=np.float)
  
          qs = np.array(self._oxi_states)
  
          epoint = - qs ** 2 * sqrt(self._eta / pi)
  
          for i in range(numsites):
              nfcoords, rij, js, _ = self._s.lattice.get_points_in_sphere(fcoords,
                                      coords[i], self._rmax, zip_results=False)
  
              # remove the rii term
              inds = rij > 1e-8
              js = js[inds]
              rij = rij[inds]
              nfcoords = nfcoords[inds]
  
              qi = qs[i]
              qj = qs[js]
  
              erfcval = erfc(self._sqrt_eta * rij)
              new_ereals = erfcval * qi * qj / rij
  
              # insert new_ereals
              for k in range(numsites):
                  ereal[k, i] = np.sum(new_ereals[js == k])
  
              if self._compute_forces:
                  nccoords = self._s.lattice.get_cartesian_coords(nfcoords)
  
                  fijpf = qj / rij ** 3 * (erfcval + forcepf * rij *
                                           np.exp(-self._eta * rij ** 2))
                  forces[i] += np.sum(np.expand_dims(fijpf, 1) *
                                      (np.array([coords[i]]) - nccoords) *
                                      qi * EwaldSummation.CONV_FACT, axis=0)
  
          ereal *= 0.5 * EwaldSummation.CONV_FACT
          epoint *= EwaldSummation.CONV_FACT
          return ereal, epoint, forces
docs: Determines the self energy -(eta/pi)**(1/2) * sum_{i=1}^{N} q_i**2