# How to calculate a structure factor from a radial distribution function?

I have been trying to take an FFT of the radial distribution function to calculate a structure factor. My understanding is this is the same as the regular FFT from the time to frequency domain. The calculation should be:

$$S(k) = 1 + \rho \int \exp(-ikr)g(r)dr$$

which often in physics is taken as a radial Fourier transform, i.e.,

$$S(k) = 1 + 4\pi \rho \int r^{2} g(r)\frac{\sin(kr)}{kr}dr$$

where $$\rho$$ is the density of the fluid. This is what it should look like for the simplest such system, argon: I imagine I can just follow Eq. 1. I have the data $$g(r)$$ for that model system, so I do the following:

nfft = 2^nextpow2(length(gr));
dt = r(2) - r(1);
df = 1/dt;
Freq = (df/2)*linspace(0,1,nfft/2+1);
sk = fft(gr,nfft)/length(gr);


And then I think I should take the one-sided FFT and the amplitude following Eq. 1, so I plot the following:

plot(Freq,1+rho*2*abs(sk(1:numel(Freq))))


But I get something radically different from the known case, and it doesn't oscillate about 1. Note that $$q$$ is interchangeable with $$k$$. I'm doing something wrong, but any ideas what? In posting here, I feel that my question is equal parts Matlab and computational physics.

• Looks like you are doing a 1D FT - the first equation above should be a 3D transform with a complex factor of Exp(-ik.r) - i.e. k and r are vectors. This reduces to the second form for spherically symmetric g(r). Mar 23, 2022 at 22:45
• math.stackexchange.com/questions/3527070/… might be of use Mar 24, 2022 at 9:33

$$S({\bf k}) = 1 + \rho \int \int \int \exp(-i{\bf k}\cdot{\bf r})g({\bf r })dr_x dr_y dr_z$$
The second equation is correct for spherically symmetric $$g({\bf r})$$, which is often assumed to be the case. Thus for such functions the triple integral in equation one can be reduced to the single integral in equation two - possibly the simplest way to show this is to expand the plane wave in spherical harmonics and then note as $$g({\bf r})$$ is spherically symmetric only the $$l=0$$ terms survive due to the orthogonality properties of spherical harmonics.
To use equation 2 to evaluate your structure factor, see this answer on the sister maths site. You will need to perform a Sine transform of $$r g({\bf r})$$ and halve the result for this - I don't speak much python, but according to this site it looks like dst() is what you are looking for.