Perform the reciprocal space summation. Calculates the quantity
        E_recip = 1/(2PiV) sum_{G < Gmax} exp(-(G.G/4/eta))/(G.G) S(G)S(-G)
        where
        S(G) = sum_{k=1,N} q_k exp(-i G.r_k)
        S(G)S(-G) = |S(G)|**2

        This method is heavily vectorized to utilize numpy's C backend for
        speed.

def _calc_recip(self):
        """
        Perform the reciprocal space summation. Calculates the quantity
        E_recip = 1/(2PiV) sum_{G < Gmax} exp(-(G.G/4/eta))/(G.G) S(G)S(-G)
        where
        S(G) = sum_{k=1,N} q_k exp(-i G.r_k)
        S(G)S(-G) = |S(G)|**2

        This method is heavily vectorized to utilize numpy's C backend for
        speed.
        """
        numsites = self._s.num_sites
        prefactor = 2 * pi / self._vol
        erecip = np.zeros((numsites, numsites), dtype=np.float)
        forces = np.zeros((numsites, 3), dtype=np.float)
        coords = self._coords
        rcp_latt = self._s.lattice.reciprocal_lattice
        recip_nn = rcp_latt.get_points_in_sphere([[0, 0, 0]], [0, 0, 0],
                                                 self._gmax)

        frac_coords = [fcoords for (fcoords, dist, i, img) in recip_nn if dist != 0]

        gs = rcp_latt.get_cartesian_coords(frac_coords)
        g2s = np.sum(gs ** 2, 1)
        expvals = np.exp(-g2s / (4 * self._eta))
        grs = np.sum(gs[:, None] * coords[None, :], 2)

        oxistates = np.array(self._oxi_states)

        # create array where q_2[i,j] is qi * qj
        qiqj = oxistates[None, :] * oxistates[:, None]

        # calculate the structure factor
        sreals = np.sum(oxistates[None, :] * np.cos(grs), 1)
        simags = np.sum(oxistates[None, :] * np.sin(grs), 1)

        for g, g2, gr, expval, sreal, simag in zip(gs, g2s, grs, expvals,
                                                   sreals, simags):

            # Uses the identity sin(x)+cos(x) = 2**0.5 sin(x + pi/4)
            m = (gr[None, :] + pi / 4) - gr[:, None]
            np.sin(m, m)
            m *= expval / g2

            erecip += m

            if self._compute_forces:
                pref = 2 * expval / g2 * oxistates
                factor = prefactor * pref * (
                    sreal * np.sin(gr) - simag * np.cos(gr))

                forces += factor[:, None] * g[None, :]

        forces *= EwaldSummation.CONV_FACT
        erecip *= prefactor * EwaldSummation.CONV_FACT * qiqj * 2 ** 0.5
        return erecip, forces