I am working on a pair correlation function

I am working on a pair correlation function, where I have this code so far. I am reading everywhere that in the limit of 'r' to infinity g(r) should go to 1. But I cannot figure out how this works. I have the following code so far. Something is wrong with my normalization factor I think, but I can't see what.. Anyone any idea? (data file is to big to store here..)

% script to calculate the Radial Distribution Function (RDF)
% This is for particles of the same type

% atom type you want to analyse
element = 'Na_pv';

%extract data from MD file
element_numbers = find(strcmp(sim_data.atom_element, element)==1); %rows which correspond to the target elements
all_pos = sim_data.frac_pos(:,element_numbers,:); % Fractional coordinates during the MD, for just the target elements
lattice = sim_data.lattice; % define lattice

% we iterate the following proces over 10 randomly chosen particles to
% create a more accurate RDF

%define bins to store the data in
bins = linspace(1,14,1000);
%generate 1000 random timesteps to limit computation time
time_steps = sort(randi([0 sim_data.nr_steps],1,1000));
% bulk density is total number of particles(48) divided by total volume in
% Angstrom
bulk_density = number_particles/(sim_data.volume*volume_factor);
% resolution is my delta r
resolution = 0.0110;

%Calculating all inter distances for every timestep:
%distance_all_timesteps =  48x48x1000 matrix
for i = 1:length(element_numbers)
    for j = 1:length(element_numbers)
        for k = 1:length(time_steps)
            distance_all_timesteps(i,j,k) = calc_dist_sqrd_frac(all_pos(:,i,time_steps(k)), all_pos(:,j,time_steps(k)), lattice);
        end
    end
end

%bin all the different distances in to the different shells
for i = 1:length(element_numbers)
    for j = 1:length(element_numbers)
        [N, edges] = histcounts(distance_all_timesteps(:,:,:),bins);
    end
N = N/2; % distance between particle 1 and 2 is the same as distance between 2 and 1
end
centers = (edges(1:end-1) + edges(2:end))/2;

%calculate all the normalization factors for a 3D system
for i = 1:length(centers) 
    normalization_factors(i) = bulk_density*4*pi*(centers(i))^2 * resolution;
end

%divide all the binned data from N by the normalization factors
for i = 1:length(centers)
    final(i) = N(i)/normalization_factors(i);
end

plot(centers, final)
title('Radial distribution function')
xlabel('r (Angstrom)')
ylabel('g(r)')
%ylim(0, 1000);

Answer 1

[Shrey Tripathi]

The issue with your code might lie in the formula for calculation of the normalization factors. To compute the normalization factor for a spherical shell, we need to consider the volume of the shell, not just the surface area.

We can calculate the normalization factors like so:

% Calculate all the normalization factors for a 3D system
for i = 1:length(centers) 
    shell_volume = (4/3) * pi * (edges(i+1)^3 - edges(i)^3);
    normalization_factors(i) = bulk_density * shell_volume * resolution;
end

Here, the variable shell_volume represents the volume of the spherical shell between edges(i) and edges(i+1). The formula (4/3) * pi * (edges(i+1)^3 - edges(i)^3) calculates the volume of the shell using the difference between the cubes of the radii.

By multiplying the shell_volume with bulk_density and resolution, we can obtain the correct normalization factor for each shell.

With this correction, the g(r) values should be properly normalized, and in the limit of large r, g(r) should approach 1.

[pepijn kluijtmans]

thanks for you answer! But it still seems to converge to 0.;

I just don't see how the values in final(i) should become 1 if N(i) for large i is 0, as there are no particles pairs detected at such distances in my unitcell. Does it have something today with periodic boundary conditions? Those are implemented in the calc_dist_sqrd_frac() function

% script to calculate the Radial Distribution Function (RDF)
% This is for particles of the same type

% atom type you want to analyse
element = 'Na_pv';

element_numbers = find(strcmp(sim_data.atom_element, element)==1); %rows which correspond to the target elements
all_pos = sim_data.frac_pos(:,element_numbers,:); % Fractional coordinates during the MD, for just the target elements
lattice = sim_data.lattice; % define lattice

% we iterate the following proces over 10 randomly chosen particles to
% create a more accurate RDF
number_particles = length(element_numbers);
volume_factor = 10^30;
reference_particles = sort(randi([min(element_numbers) max(element_numbers)],1,20));
time_steps_1 = sort(randi([0 sim_data.nr_steps],1,1000));
time_steps_2 = 1:1000;
bins = linspace(1,14,1000);
bulk_density = number_particles/(sim_data.volume*volume_factor);
resolution = 0.0110;
volume = sim_data.volume * volume_factor;


for i = 1:length(element_numbers)
    for j = 1:length(element_numbers)
        for k = 1:length(time_steps_1)
            distance_all_timesteps(i,j,k) = calc_dist_sqrd_frac(all_pos(:,i,time_steps_1(k)), all_pos(:,j,time_steps_1(k)), lattice);
        end
    end
end

for i = 1:length(element_numbers)
    for j = 1:length(element_numbers)
        [N, edges] = histcounts(distance_all_timesteps(:,:,:),bins);
    end
N = N/2000;
end


centers = (edges(1:end-1) + edges(2:end))/2;

%calculate all the normalization factors for a 3D system
% for i = 1:length(centers) 
%     normalization_factors(i) = 4*pi*(centers(i))^2 * resolution;
% end

% Calculate all the normalization factors for a 3D system
for i = 1:length(centers) 
    shell_volume = (4/3) * pi * (edges(i+1)^3 - edges(i)^3);
    normalization_factors(i) = bulk_density * shell_volume * resolution;
end

%divide all the binned data from N by the normalization factors
for i = 1:length(centers)
    final(i) = (N(i)/normalization_factors(i))%* (2*volume)/(number_particles * (number_particles-1));
end
%% plot

plot(centers, final)
title('Radial distribution function')
xlabel('r (Angstrom)')
ylabel('g(r)')
%ylim(0, 1000);

[pepijn kluijtmans]

With the following formula attache, I get at least nice values for the y-axis. Result also attached