Calculate Eq. (7)

        k = K + G + g

        K -> _qpoints_s[_q_index] (in SBZ)
        k -> _qpoints_p[_q_index] (in PBZ)
        G -> _comm_points (in PBZ)
        g -> a reciprocal lattice point in PBZ (unused explicitly)
        j -> Primitive translations in supercell (_trans_p)
        J -> Band indices of supercell phonon modes (axis=1 or eigvecs)

        The phase factor corresponding to K is not included in eigvecs
        with our choice of dynamical matrix.

def _get_unfolding_weights(self):
        """Calculate Eq. (7)

        k = K + G + g

        K -> _qpoints_s[_q_index] (in SBZ)
        k -> _qpoints_p[_q_index] (in PBZ)
        G -> _comm_points (in PBZ)
        g -> a reciprocal lattice point in PBZ (unused explicitly)
        j -> Primitive translations in supercell (_trans_p)
        J -> Band indices of supercell phonon modes (axis=1 or eigvecs)

        The phase factor corresponding to K is not included in eigvecs
        with our choice of dynamical matrix.

        """

        eigvecs = self._eigvecs[self._q_count]
        dtype = "c%d" % (np.dtype('double').itemsize * 2)
        weights = np.zeros(eigvecs.shape[1], dtype=dtype)
        q_p = self._qpoints_p[self._q_count]  # k
        q_s = self._qpoints_s[self._q_count]  # K
        diff = q_p - np.dot(q_s, np.linalg.inv(self._supercell_matrix))

        # Search G points corresponding to k = G + K
        for G in self._comm_points:
            d = diff - G
            d -= np.rint(d)
            if (np.abs(d) < 1e-5).all():
                break

        e = np.zeros(eigvecs.shape[:2], dtype=dtype)
        phases = np.exp(2j * np.pi * np.dot(self._trans_p, G))
        for phase, indices in zip(
                phases, self._atom_mapping[self._index_map_inv]):
            eig_indices = (
                np.c_[indices * 3, indices * 3 + 1, indices * 3 + 2]).ravel()
            e += eigvecs[eig_indices, :] * phase
        e /= self._N
        weights[:] = (e.conj() * e).sum(axis=0)

        # e = np.zeros(eigvecs.shape[:2], dtype=dtype)
        # phases = np.exp(2j * np.pi * np.dot(self._trans_p, G))
        # indices = self._atom_mapping
        # eig_indices_r = (
        #     np.c_[indices * 3, indices * 3 + 1, indices * 3 + 2]).ravel()
        # for phase, indices in zip(
        #         phases, self._atom_mapping[self._index_map_inv]):
        #     eig_indices_l = (
        #         np.c_[indices * 3, indices * 3 + 1, indices * 3 + 2]).ravel()
        #     e += eigvecs[eig_indices_l, :] * eigvecs[eig_indices_r, :].conj() * phase
        # e /= self._N
        # weights[:] = e.sum(axis=0)

        if (weights.imag > 1e-5).any():
            print("Phonopy warning: Encountered imaginary values.")

        # assert (np.abs(weights.real.sum(axis=1) - 1) < 1e-5).all()

        return weights.real