Structure factor in a homogenous system I want to calculate the structure factor for a homogenous system. The system that I am dealing with is the results of a Vicsek type model simulation.
The structure factor is defined as :

$$S(q) = \frac{1}{N}\left|\sum_{j=1}^{N}{e^{i \vec{q}\cdot\vec{r_j}}}\right|^2.$$

where $q$ is the wave vector and $r$ is the position of the particles. You may have seen other forms of the structure factor as well. Now let's go to my question.
I read somewhere that for a system like what I mentioned above, where all particles are the same as each other and there is not any specific lattice spacing, there is no need to behave $\vec{q}$ like what people do in crystallography (like doing $\vec{q} \propto (n_x,n_y,n_z)$. All that is important is the magnitude of $\vec{q}$. Now my confusion is that for calculating $S(q)$, should I replace $\vec{q}$ with $|q|$ and take it as a scalar number.
You can take a look at the following papers as references :
https://pdfs.semanticscholar.org/b3ce/f60a57b11bc07bc53187a29e395de9a71dfe.pdf?_ga=2.253183808.812496979.1554644890-228537717.1509991838
https://iopscience.iop.org/article/10.1088/1361-6463/aa7f8e/pdf
 A: You need to include the ensemble average (denoted usually by angle brackets) in your definition of $S(q)$.
If you are performing a Vicsek model simulation, then most likely you are using periodic boundary conditions, as they do in the first paper that you cited. In that paper they explain clearly what to do: choose a set of $\vec{q}$-vectors compatible with the boundary conditions (5000 of them, in their case), evaluate $S(\vec{q})$ according to their eqn (2), which is equivalent to the formula in your question if we include the ensemble averaging, for every wave vector; and then bin up the final results according to $|\vec{q}|$. In this last step you will be averaging over the direction of the wave vector.
To be clear, in a cubic periodic box of side $L$, the wave-vector will have components
$$
\vec{q}= (n_x,n_y,n_z) \frac{2\pi}{L}
$$
where $n_x,n_y,n_z$ are integers. It would be usual to calculate $S(\vec{q})$ for all positive and negative integer values of $n_x,n_y,n_z$, up to some maximum magnitude $q_{\text{max}}$, i.e. all values for which 
$$
n_x^2+n_y^2+n_z^2 \leq  \frac{q_\text{max}^2 L^2}{4\pi^2} .
$$
Actually symmetry allows you to omit nearly half of these, since $(-n_x,-n_y,-n_z)$ is equivalent to $(n_x,n_y,n_z)$. I presume that the 5000 wave-vectors used by Löwen, Smallenburg and their colleagues will have been chosen in this way, i.e. all the distinct ones that are allowed by the periodic boundaries, counting upwards in magnitude, up to a maximum $q_{\text{max}}$. If we imagine a sphere of radius $n_\text{max}$ in $(n_x,n_y,n_z)$-space sufficient to contain 10000 points (half of which will be redundant) on the cubic reciprocal lattice, then
$$
n_x^2+n_y^2+n_z^2 \leq n_\text{max}^2 \approx \left(\frac{3}{4\pi}\times 10000\right)^{2/3} \approx 180 .
$$
After calculating all these, the results need to be expressed in terms of the magnitude of $\vec{q}$, averaging over the direction of $\vec{q}$. Each vector $\vec{q}$ will have a magnitude $q=|\vec{q}|$. Because of the cubic symmetry, many of them will have the same magnitude. Moreover, the magnitudes of these vectors will not be evenly spaced. So the averaging should also be over vectors $\vec{q}$ that lie within a certain (narrow) range of magnitudes. The idea is to construct a histogram $\bar{S}(q)$ with regularly spaced bins of width $\Delta q$:  $[0\ldots \Delta q]$, $[\Delta q\ldots 2\Delta q]$, $[2\Delta q\ldots 3\Delta q]$,  etc. You would post-process your $S(\vec{q})$ results such that each histogram bin $\bar{S}(q)$ contains the average of all your results $S(\vec{q})$ which satisfy the condition that $|\vec{q}|$ lies within the corresponding range of $q$. So the result $S(\vec{q})$ for every vector $\vec{q}$ will end up as a contribution to the average in the appropriate bin of your final result $\bar{S}(q)$.
It is possible to rewrite the double sum over pairs of particles in terms of the magnitude of the wave vector, and the distances between the pairs of particles. This involves conducting an average over the angle between $\vec{q}$ and $\vec{r}_j-\vec{r}_k$, relying on the isotropy of the system. The resulting expression is a bit more complicated, but can be computed as a double sum over $j$ and $k$. But then you lose the nice form involving a single sum over particles, so it is less common to do it that way.
