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

Question

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.

Answer 1

Without data to test exactly what you are doing it's difficult to say for absolute certain, but it looks to me that your main problem is that your first equation, the equation you are using for your program, is wrong. It should be a 3D transform. Writing it out in full to make the point it should be

$$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.