How to calculate structure factor correctly using solely radial distribution function statistics?

Below is my code trying to calculate the structure factor (sq as written in this code) from the radial distribution function statistics (Gr).

The first column of data in Gr is scalar distance data with a step size of 2. Since it is not refined enough, I use linear interpolation to refine the data here. Then I tried to calculate the structure factor using the formula in the loop, but I didn't get a good result.

As you can see, there are several problems with my current structure factor calculations:

First, there will be negative values in the loglog graph, which will be automatically ignored by matlab, resulting in discontinuity of the graph.

Second, under theoretical circumstances, the result of my system should be that as the wave vector decreases, when it approaches 0, in the loglog plot, the structure factor will also approach 0 linearly. But you can see from my results that something is obviously wrong here.

Could anybody tell me how can I solve this problem?

Also, I need to do some explanation for the code I'm showing so far.

You can see that the values of my interpolated x start from 24. This is because my rdf result will be 0 in a certain range. I think this will have an impact on the result, so I remove it before calculation. And when r is very large, large vibrations will occur, so I only got 400.

interp_x = linspace(24, 400, 4000);

interp_y = interp1(Gr(:, 1), Gr(:, 2), interp_x, 'linear');

interp_x = interp_x'; interp_y = interp_y';

newGr(:, 1) = interp_x; newGr(:, 2) = interp_y;

q(:, 1) = 0:0.0001:3;

for m = 1:length(q)

J0 = [];

J0 = 2*pi*besselj(0, (q(m, 1) * newGr(:, 1)));

sq(m, 1) = 1 + trapz(newGr(:, 1), rho0 * newGr(:, 1) .* newGr(:, 2) .* J0 .* ...

(sin(pi * newGr(:, 1) / newGr(end, 1)) ./ (pi * newGr(:, 1) / newGr(end, 1))));

end

figure;plot(newGr(:, 1), newGr(:, 2), '.-');

figure; loglog(q, sq, '.-');

[Xuemao Zhou]

Hi, I am also trying to calculate structure factor from the radial distribution function.But, the equation I found is related to sin(qr)/qr, instead of the Bessel function of 0 order and the first kind as you have used in your code "besselj(0,qr). The second thing is, I also have negative values in the structure factor. I have no clue at all. Have you solve this problem? Can you share your experience? Thank you in advance.

Answer 1

[SOUMNATH PAUL]

Hi,

To my understanding, you need help in computing the structure factor (sq) from radial distribution function statistics (Gr) using linear interpolation , while avoiding the negative values and unexpected behaviours in the plot.

Here are some points to improve your code and resolve the issues you mentioned:

1. After interpolation, check for negative values in “interp_y”. If there are negative values, consider setting them to zero or apply a threshold to avoid issues during further calculations.

Example:

interp_y(interp_y < 0) = 0; 

2. If you want to handle the discontinuity in the log-log plot, you can consider using “semilogy” instead of “loglog”. This will plot the y-axis on a logarithmic scale while keeping the x-axis linear.

Example:

semilogy(q, sq, '.-'); 

3. When interpolating, handle edge cases where interp_x or interp_y might contain NaN or Inf values. You can use “isnan” and “isinf” functions to identify and handle these cases appropriately.

Example:

interp_x = interp_x(~isnan(interp_y) & ~isinf(interp_y)); 

interp_y = interp_y(~isnan(interp_y) & ~isinf(interp_y)); 

4. When performing integration, ensure that the integration limits are appropriate for your data. If there are specific ranges where you expect the integral to be meaningful, adjust the integration limits accordingly.

Example:

integration_limits = [min(interp_x), max(interp_x)]; % Adjust limits based on your data 

sq(m, 1) = 1 + trapz(interp_x, ... , 'abstol', 1e-10, 'reltol', 1e-10); 

5. For very small values in sq, you can consider using the “eps” function to set them to a small positive value to prevent issues related to numerical precision.

6. Add additional plots to visualize intermediate results, such as the interpolated data and the integrand, to help identify potential issues in your calculations.

Example:

figure; 

plot(interp_x, interp_y, '.-'); 

title('Interpolated Data'); 

 
figure; 

plot(interp_x, integrand, '.-'); 

title('Integrand for Integration');

By implementing these suggestions, you can improve the accuracy and reliability of your calculations and address the issues you mentioned.

Hope it helps!

Regards,

Soumnath