Radial distribution function for an ideal gas¶
This week we will compute the radial distribution function for an ideal gas, which will get us ready to do it for a real system next time
First we have to generate a "configuration" of an ideal gas. We'll do it in 2D. Fortunately, this is easy because it just means the x and y for each particle are random between 0 and L.
#setup the notebook
%pylab inline
import numpy as np
Intro, random number generator reminder¶
# last time we saw that if we want N random numbers we can do
# np.random.random(size=N)
test_array = np.random.random(size=6)
print("raw array:",test_array)
# if we want an array of dimensions 3 by 2 instead, we can use the reshape function
print( "reshaped array:\n",test_array.reshape(3,2) )
#we can also generate a random array with this shape right away
test_array_2 = np.random.random(size=(3,2))
print("raw array 2:\n",test_array_2)
# we can get the minimum and maximum of a numpy array using the min and max functions
long_array = np.random.random(size=1000)
print( "long array min:", long_array.min() )
print( "long array max:", long_array.max() )
#these should be very close to 0 and 1, because np.random.random returns numbers in the range 0 to 1.
#Let's say we want to convert long_array into a new array, long_array_2, with numbers in the range 0 to 5.
# **Figure out how to do that** and make sure that you get min and max close to 0 and 5
long_array_2 = 5*long_array
print( "long array 2 min:", long_array_2.min() )
print( "long array 2 max:", long_array_2.max() )
Part 1, Generating configurations¶
#this will be a function that returns a configuration
def generate_configuration(N_particles, box_L):
#use np.random.random and other things from the previous part to generate a configuration
#which is an array with N_particles rows and 2 columns for x and y
# the particle positions in x and y should be between 0 and box_L
configuration = box_L * np.random.random(size=(N_particles,2))
return configuration
#this will draw a single configuration
def plot_configuration(configuration,box_L):
plt.figure(figsize=(4,4))
plt.axis('equal')
plt.scatter(configuration[:,0],configuration[:,1])
plt.xlim(0,box_L)
plt.ylim(0,box_L)
#this will test if your function is working correctly
test_box_L = 20
test_num_particles = 100
configuration1 = generate_configuration(test_num_particles,test_box_L)
plot_configuration(configuration1,test_box_L)
Part 2, generating a list of distances¶
#let's write a function that gives a list of distances between particles, and test it on configuration 1
# there are fancier ways to do this, but let's do it with 2 for loops
def compute_distances( configuration ):
distance_list = []
num_particles = configuration.shape[0]
for i in range(num_particles):
for j in range(num_particles):
# we don't want to count distances of particles to themselves
if i == j: continue
posi = configuration[i]
posj = configuration[j]
# compute the euclidian distance between pos1 and pos2 and call it 'dist'
# there are many ways to do this
# you can certainly look up how to do this online if you can't figure it out right away
#dr is a vector (dx,dy)
dr = posj-posi
#dr2 is a vector (dx*dx,dy*dy)
dr2 = dr*dr
#dist = sqrt( dx^2 + dy^2)
dist = np.sqrt( dr2.sum() )
distance_list.append(dist)
return np.array(distance_list)
#the length of the returned distance list should be N*(N-1), let's check
distance_list_1 = compute_distances(configuration1)
print("Does the number of distances equal the number expected?")
print(len(distance_list_1),"=", test_num_particles*(test_num_particles-1) )
Part 3, histogramming the distances¶
Now we will histogram the distance list that we just made. That means count how often a distance occurs between r and r+dr. We can use the numpy histogram function. We can tell it which bins to count inside of.
We will generate a histogram only from 0 to box_L/2 for reasons that will make more sense in part 4
First, look at the documentation for the [numpy histogram] (https://docs.scipy.org/doc/numpy/reference/generated/numpy.histogram.html) function, and the examples!
#first, let's check how close and far apart the particles are
print("Min distance:",distance_list_1.min())
print("Max distance:",distance_list_1.max())
def histogram_distances_default(distance_list):
hist, bin_edges = np.histogram( distance_list )
return hist, bin_edges
def plot_histogram(hist,bin_edges):
#for N bins, there are N+1 bin edges. The centers can be found by averaging the positions of
# bin edge0 and 1, 1 and 2, ..., N-1 and N
bin_centers = (bin_edges[:-1]+bin_edges[1:])/2.0
plot(bin_centers,hist,marker='o')
ylabel("N(r)")
xlabel("$r$")
dist_hist_1, bin_edges_1 = histogram_distances_default(distance_list_1)
plot_histogram(dist_hist_1,bin_edges_1)
As we talked about in class, the number of particles in distance away between r and r+dr grows with r. but then it goes back down. This is because of boundary effects. You get to the edge pretty fast.
Morevoer, the maximum distance 2 paricles can be apart is L*sqrt(2) (do you see why?). We will only histogram out to L/2 for reasons that will make more sense in the next section
def histogram_distances(distance_list, max_dist, bin_size):
# this is the list of bins in which to calculate
bins = np.arange(0, max_dist+bin_size, bin_size)
hist, bin_edges = np.histogram( distance_list, bins=bins )
return hist, bin_edges
#choose a bin size so that this becomes relatively smooth without throwing away too much information
#set max dist to the box size of this configuration *sqrt(2)/2.
test_bin_size = 0.5
max_dist = test_box_L/2.0*np.sqrt(2)
dist_hist_1, bin_edges_1 = histogram_distances(distance_list_1,max_dist=max_dist, \
bin_size=test_bin_size)
plot_histogram(dist_hist_1,bin_edges_1)
Part 4, Periodic Boundary Conditions¶
In simulations, we always have a limit of a finite system size. To get rid of some of the boundary effects, we can use periodic boundary conditions. That means we compute distances as if the boundaries wrap around like a donut. The distance calculation we use is called minimum image convention. That means the distance between to particles is the smallest distance between those two particles in infiniate copies of the system. This make more sense if you make the plot in the next cell:
def plot_configuration_pbc(configuration,box_L):
plt.figure(figsize=(4,4))
plt.axis('equal')
for shift_x in range(-1,2):
for shift_y in range(-1,2):
if shift_x == 0 and shift_y == 0:
#this is for plotting transparency
alpha = 1
else:
alpha = .3
plt.scatter(shift_x*box_L+configuration[:,0],shift_y*box_L+configuration[:,1],alpha=alpha)
plt.xlim(-box_L,2*box_L)
plt.ylim(-box_L,2*box_L)
plot_configuration_pbc(configuration1,test_box_L)
Write a minimum image distance calculator¶
Replace your old dist= line with one that takes into account that there might be a closer copy of particle j to particle i than the one in the main box. There are several ways to do this, and you can easily google how to do it
def compute_distances_minimum_image( configuration, box_size ):
distance_list = []
num_particles = configuration.shape[0]
for i in range(num_particles):
for j in range(num_particles):
if i == j: continue
posi = configuration[i]
posj = configuration[j]
# compute the euclidian distance between pos1 and pos2 and call it 'dist'
# there are many ways to do this
# you can certainly look up how to do this online if you can't figure it out right away
#dr is a vector (dx,dy)
dr = posj-posi
#minimum image dr - can you figure out why this works?
dr = dr - box_size*np.floor(dr/box_size+0.5)
#dr2 is a vector (dx*dx,dy*dy)
dr2 = dr*dr
#dist = sqrt( dx^2 + dy^2)
dist = np.sqrt( dr2.sum() )
distance_list.append(dist)
return np.array(distance_list)
distance_list_1_pbc = compute_distances_minimum_image(configuration1, box_size = test_box_L)
#now, no two particles should be farther than L*sqrt(2)/2, do yo use why?
print("Min distance:",distance_list_1_pbc.min())
print("Max distance:",distance_list_1_pbc.max())
test_bin_size = .25
#set max dist to the box size of this out to L*sqrt(2)/2
max_dist = test_box_L/2.0*np.sqrt(2)
dist_hist_1_pbc, bin_edges_1_pbc = histogram_distances(distance_list_1_pbc,max_dist=max_dist, bin_size=test_bin_size)
plot_histogram(dist_hist_1_pbc,bin_edges_1_pbc)
plt.axvline(test_box_L/2.,color='black',linestyle='--')
You should see something interesting happen in this plot at L/2. Try imagining the 3x3 copies of our system shown above, and see which particles are within a circle of radius r.
What you see in the plot above has to do with the difference in area between a circle of radius r and a square of radius r. For this reason, we will only look up until distance L/2.
In the next cell, remake and plot this histogram up to distance L/2 and choose a bin size so that it is relatively smooth. We will use this data in the next part
test_bin_size = 0.5
max_dist = test_box_L/2.0
dist_hist_1_pbc, bin_edges_1_pbc = histogram_distances(distance_list_1_pbc,max_dist=max_dist, bin_size=test_bin_size)
plot_histogram(dist_hist_1_pbc,bin_edges_1_pbc)
Part 5, computing the radial distribution function¶
the radial distribution function is defined as this histogram divided by the number of things in a shell of radius r around the point. In 2d, this is 2 pi r dr * density.
You also did the distance calculation for every particle, which means you essentially have computed 1 histogram for each particle and combined them. Hence we have to divide by N also.
In the get_gofr function, compute the bin_centers array and the number density. Finally compute the radial distribution function. Store it as an array called gofr, and it will be plotted.
Hint: the dr should be the spacing in your bin_edges, and the r's should be the bin_centers (see above for how to compute those).
Hint: What do you expect for an ideal gas? Is or working correctly?
def plot_rdf(gofr,bin_centers):
plot(bin_centers,gofr,marker='o')
ylabel("g(r)")
xlabel("$r$")
def get_gofr(hist,bin_edges,num_particles, box_size):
rho = num_particles/box_size/box_size
bin_centers = (bin_edges[1:]+bin_edges[:-1])/2.0
dr = bin_edges[1]-bin_edges[0]
denominator = 2.*np.pi*bin_centers*dr*rho*( num_particles )
gofr = hist/denominator
return gofr, bin_centers
gofr, bin_centers = get_gofr( dist_hist_1_pbc, bin_edges_1_pbc, test_num_particles, test_box_L )
plot_rdf(gofr, bin_centers)
Part 6, averaging the radial distribution function¶
When we do a simulation, we will generate many configurations, and average their radial distribution functions. Let's do that for our ideal gas. In the cell below, generate N random configurations, compute their g(r) using a combination of the functions above. Then average the g(r) together. Increase num_configurations until the error in g(r) is below 5% at every point
num_configurations = 50
gofr_list = []
for i in range(num_configurations):
#generate a configuration
configuration_i = generate_configuration(test_num_particles,test_box_L)
#compute pbc distances
distance_list_i_pbc = compute_distances_minimum_image(configuration_i, box_size = test_box_L)
#histogram the distances
dist_hist_i_pbc, bin_edges_i_pbc = histogram_distances(distance_list_i_pbc,max_dist=test_box_L/2., bin_size=test_bin_size)
# get_gofr (fill in the parentheses with your variable names)
gofr, bin_centers = get_gofr( dist_hist_i_pbc, bin_edges_i_pbc, test_num_particles, test_box_L )
gofr_list.append(gofr)
avg_gofr = np.mean(gofr_list, axis=0)
plot_rdf(avg_gofr, bin_centers)
axhline(1.0,linestyle='--',color='black')
ylim(0.8,1.2)
