# 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

## 1 Answer

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.

• +1 Many thanks. I'll dig into this some more. So I guess the standard way to take this transform is to do the radial transform rather than the full 3D transform in Eq. 1? Mar 24, 2022 at 14:44
• Equation 2 will be a lot cheaper both in terms of memory and time Mar 24, 2022 at 15:24
• Nice answer Ian! +1. Mar 24, 2022 at 15:26