I am working on a pair correlation function
18 views (last 30 days)
Show older comments
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);
0 Comments
Answers (1)
Shrey Tripathi
on 27 Jun 2023
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.
2 Comments
See Also
Categories
Find more on Matrix Computations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!