Hi Rob,

To calculate the structure factor O(N), where N is the number

of atoms, you can utilize S(\vec q)= | F(\vec q) |^2 where

the form factor amplitude is defined as

F(\vec q)=sum_i f_i exp(i \vec q\dot \vec r_i )

and f_i denote the scattering length of the i'th atom which

acounts for the details of the radiation-matter interaction.

Hence each process for each \vec q have to sum

f_i cos(q r_i) and f_i sin(q r_i) over all locally owned

atoms and store them in two arrays. Afterwards do

an MPI_ALLreduce(.., MPI_SUM , .. ) on the array. Then

the amplitude is known globally and cheap to square,

afterwards you can bin the structure factor after

\vec q| and average.

\vec q can only be points on the reciprocal lattice

corresponding to the box due to periodic boundary

conditions. For a fixed rectangular box (Lx,Ly,Lz):

\vec q (nx,ny,nz) = (nx 2Pi/Lx, ny 2Pi/Ly, nz 2Pi/Lz),

similar can be defined for fixed triclinic boxes.

Here FFTs is one option, but a fast and transparent

approach is to realise for a given atom and q vector \vec r_i \dot \vec q

is of the form nx a + ny b + nz c. Since cartesian directions factorise

you are left with the task of calculating cos(nx a) and sin(nx a) for

nx=0, 1, .., M. You can do nx=0,1 explicitly, and then express the case

=2 recursively with products and sums of previous values. In particlar,

for even nx cos(nx a) can be expressed with cos(a nx/2) and sin(a nx/2)

which are already known. For odd nx cos(nx a) can be expressed though

cos(a),sin(a), cos((nx-1)a) and sin((nx-1)a) which have already been

computed. That can be used to derive all the cos, sin terms up to

q><qmax evaluating only 6 trigonometric functions for each atom.

That should produce an algorithm which is O( N M^3) where N is number

of atoms and M^3 is the number of q vectors you want to evaluate.