Otdfbt

Here is an example code (to calculate the so called structure factor, please see the comment below the question) for a small data set of 2d coordinates. In reality I have a factor of up to 100 times more coordinates.

Code with AbsoluteTiming:

coordinates = Get@"https://pastebin.com/raw/wFxJ1m4U"; 

lX, lZ = 100, 100;

kX = 2 Pi/lX*(Range[lX] - lX/2) // N;
kZ = 2 Pi/lZ*(Range[lZ] - lZ/2) // N;

cx = Flatten[coordinates[[All, 1]]];
cz = Flatten[coordinates[[All, 2]]];

rX = Flatten@Outer[Differences[##] &, cx, cx];
rZ = Flatten@Outer[Differences[##] &, cz, cz];

nP = Length[coordinates];

s = Array[0 &, lZ, lX];

Do[
 Do[
 s[[j, i]] = 1/nP*Total[Cos[kX[[i]]*rX + kZ[[j]]*rZ]]
 , i, 1, lX
 ]
 , j, 1, lZ
 ]; // AbsoluteTiming

 31.3555, Null

How can I speed up this code?

A huge problem occurs for many coordinates: for about 2400 coordinates the calculation takes about 3 days!

Here are 2403 coordinates: https://pastebin.com/raw/0t1MBAgX .
Please use for this test lX, lZ = 1600, 1200

Here is the code for this large data set: https://pastebin.com/raw/MH5Pb4h0

Can the double loop be compiled?

Table works a little bit better than Do:

s = Table[Total[Cos[kX[[i]]*rX + kZ[[j]]*rZ]], j, 1, lZ, i, 1, lX]/nP;

Much more speed can be gained by first rewriting the cosines as imaginary exponentials, and then rewriting exponentials of sums as products of exponentials. As a result, rX and rZ aren't needed, and the calculation takes 0.03 seconds, or just over one second with the large data set:

s = Abs[Exp[I*KroneckerProduct[kZ, cz]].Exp[I*KroneckerProduct[cx, kX]]]^2/nP; //
 AbsoluteTiming // First
(* 1.22367 *)

A little more speed can be gained with this trick of using an internal function to calculate Abs[...]^2:

s = Internal`AbsSquare[
 Exp[I*KroneckerProduct[kZ, cz]].Exp[I*KroneckerProduct[cx, kX]]]/nP; //
 AbsoluteTiming // First
(* 1.09222 *)

A bit more speed could be gained by replacing the Kronecker products by a NestList, as the values in kX and kZ are uniformly spaced and thus the many exponentials can be replaced by repeated multiplications of a single exponential. This is numerically less stable though, and in any case at this point 50% of the time are spent in the Dot product which is unavoidable, so the best remaining speedup you could hope for is another factor of two.

Lesson of the day: don't think about compiling before you've tried rewriting the algorithm.

mathematical derivation

Starting from Eq. (2.6) of https://arxiv.org/pdf/1001.3342.pdf

$$
S_vecq=frac1Nsum_i j cos[vecqcdot(vecr_i-vecr_j)]
= frac1NResum_i j e^texti vecqcdot(vecr_i-vecr_j)
= frac1NResum_i j e^texti vecqcdotvecr_ie^-textivecqcdotvecr_j\
= frac1NReleft[sum_i e^texti vecqcdotvecr_iright]left[sum_j e^-texti vecqcdotvecr_jright]
= frac1Nleft|sum_i e^texti vecqcdotvecr_iright|^2
= frac1Nleft|sum_i e^texti q_xr_i,xe^texti q_zr_i,zright|^2
$$

From there it's a matter of calculating all the possible combinations of $e^texti q_xr_i,x$ and $e^texti q_zr_i,z$ with all desired $vecq$-values (in kX and kZ), and then summing over them. These operations are most efficiently done with vector operations like KroneckerProduct, listable Exp, Dot-products, and finally a listable Abs or AbsSquare.