Determines the self energy -(eta/pi)**(1/2) * sum_{i=1}^{N} q_i**2

def _calc_real_and_point(self):
        """
        Determines the self energy -(eta/pi)**(1/2) * sum_{i=1}^{N} q_i**2
        """
        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