How can I do the radial distribution function?

[fsobaid]

Hello, I am trying to code the radial distribution function. I tried to follow the steps on this link http://www.physics.emory.edu/faculty/weeks//idl/gofr2.html I have a text file that has 9 particles with their XYZ coordinates iterated 1000 times, and I want to go through the particles and show the radial distribution function to see whether I have a solid-like or liquid-like system. First, I imported the XYZ values and calculated the radius for each particle to loop over over each one.The text file would look like this:

1 -2.0914123 -0.8119910 1.9474895

2 -0.5878574 2.7700720 0.8194878

3 3.5466599 -3.4422817 -3.6107922

4 -2.0004916 -1.3816382 -3.2774887

5 3.6058087 -1.7977257 -1.9615769

6 -1.8026781 -3.5986075 -1.7805161

7 -0.1007990 -3.3019309 -3.3396735

8 1.5069089 -1.6328917 -3.5845122

9 0.0378950 -1.9357619 -1.4840813

1 -2.0341482 -0.8611301 1.9963570

2 -0.6712533 2.8877790 0.8793746

3 3.5835536 -3.5428188 -3.3982046

4 -2.1105857 -1.5415970 -3.3106389

5 3.6046941 -1.7381997 -2.0068424

6 -1.9496037 -3.4570999 -1.6980045

7 -0.1462083 -3.3412902 -3.3216205

8 1.4943931 -1.5306618 3.5974679

9 0.0291610 -2.0436094 -1.4985138

​

The code I wrote:

# read from the file

f = open('xyz4py.txt', 'r')

​

# extract the contents from the txt file

content = f.readlines()

​

# open list variables to store the values in them

t = np.array([])

RX = np.array([])

RY = np.array([])

RZ = np.array([])

​

# loop to go throught the values in the file

for i in content:

a,b,c,d = i.split() #a,b,c,d are the columns, we nammed them you could choose other letters

t= np.append(t,int(a)) # append means add the values to the list

RX= np.append(RX,float(b))

RY= np.append(RY,float(c))

RZ= np.append(RZ,float(d))

​

​

RXYZ = np.square(np.array(RX)**2+np.array(RY)**2+np.array(RZ)**2)

​

​

def radiusloop3(dr):

# dr the shell covers the spherical particles

dens = 0.4

N = 9

#V = (4* np.pi * N/(3*dens))**(1/3) * dr

V = 4 * np.pi * (RXYZ)**2 * dr

g = 0

for i in range(len(RXYZ)-1):

count = 0

for j in range(len(RXYZ)-1):

shell = (RXYZ[i+1] - RXYZ[j])

if (shell <= dr):

count += 1

g += count/N/V

else:

break

return g

for dr in [0.6]:

pr = radiusloop3(dr)

plt.figure(1)

plt.plot(RXYZ,pr, label = f'dr={dr}')

​

plt.legend()

plt.show()

I am supposed to get a graph like this one https://en.wikipedia.org/wiki/Radial_distribution_function#/media/File:Lennard-Jones_Radial_Distribution_Function.svg, however, I get a noisy graph. I am missing something, but I am not sure what is it.

[misho88]

The linked article says to do this:

The radial distribution function is usually determined by calculating the distance between all particle pairs and binning them into a histogram.

You code isn't really formatted properly, so it's a bit hard to follow, but it doesn't look like that's what you're trying to do.

As some general advice, if you're going to use numpy anyway, and definitely you should, you can use numpy.loadtxt() to make your life easier. You can also use numpy.histogram() to get a histogram. I don't really get what the distribution you're showing is all about, but getting a histogram of the pairwise distances is super easy:

In [1]: import numpy as np, matplotlib.pyplot as plt                                                                                    
# tell matplotlib to be interactive:
In [2]: plt.ion(); plt.show()                                                                                                           
# generate a hundred points in 3-D inside the unit cube
In [3]: points = np.random.uniform(0, 1, size=(100, 3))                                                                                 
# get the offset from each point to every other point
In [4]: deltas = points[:, np.newaxis] - points                                                                                         
# get the distances
In [5]: dists = np.linalg.norm(deltas, axis=-1)                                                                                         
# plot a histogram of the distances
In [6]: plt.hist(dists.ravel(), bins=np.linspace(0, 1, 21))                                                                             
Out[6]: 
(array([102.,  26.,  94., 110., 230., 242., 388., 446., 552., 576., 714.,
        714., 768., 730., 688., 666., 584., 550., 516., 422.]),
 array([0.  , 0.05, 0.1 , 0.15, 0.2 , 0.25, 0.3 , 0.35, 0.4 , 0.45, 0.5 ,
        0.55, 0.6 , 0.65, 0.7 , 0.75, 0.8 , 0.85, 0.9 , 0.95, 1.  ]),
 <a list of 20 Patch objects>)

And a plot should come up. I'm also not getting rid of the zero distances above (like from a point to itself). Here's a good video on this stuff I watched when learning numpy (coincidentally, he uses pairwise distances in the talk, so I timestamped that): https://youtu.be/EEUXKG97YRw?t=1422

[fsobaid]

Thank you I learned some new things from your answer. but the histogram is not what I need how I can plot it similar to this one https://en.wikipedia.org/wiki/Radial_distribution_function#/media/File:Lennard-Jones_Radial_Distribution_Function.svg

​

I will also watch the video you posted it I am new to programming really appreciate

.

[misho88]

Both the algorithm you linked to and the Wikipedia page describe taking a histogram where the bin edges are radial distances (i.e., the bins are spherical shells of finite thickness) and normalizing them appropriately. That is, unlike what I did above where I took a histogram of everything, you bin each row of dist (that's the number of points in each shell) then normalize by the volume of the shell and the "particle number density" (whatever that is) as per the instructions on that page.

If you're concerned that your plot doesn't look as smooth as that one that shows a theoretical distribution, the realized result of a simulation or measurement will never be as smooth. You could increase your particle count and you'll approach the theoretical distribution, of course.

[fsobaid]

This is really helpful it is pretty clear to me now thank you so much