Skip to article frontmatterSkip to article content

Chapter 23: Radial Distribution Function

Learning Objectives

By the end of this lecture, you will be able to:

  1. Define the radial distribution function.
  2. Compute the radial distribution function for a given configuration of particles.
  3. Understand the physical significance of the radial distribution function.

Radial Distribution Function

Source
import numpy as np
import matplotlib.pyplot as plt

# Parameters
central_particle = (0, 0)
num_particles = 100
box_size = 10
r = 3.0  # Shell radius
dr = 0.5  # Shell thickness

# Generate random particle positions
np.random.seed(42)  # For reproducibility
particle_positions = np.random.uniform(-box_size/2, box_size/2, size=(num_particles, 2))

# Plot setup
fig, ax = plt.subplots(figsize=(8, 8))
ax.set_xlim(-box_size/2, box_size/2)
ax.set_ylim(-box_size/2, box_size/2)
ax.set_aspect('equal', adjustable='datalim')

# Draw central particle
ax.plot(*central_particle, 'o', color='red', label='Central Particle')

# Draw surrounding particles
for pos in particle_positions:
    ax.plot(*pos, 'o', color='blue', markersize=5, alpha=0.7)

# Draw the shell
circle_inner = plt.Circle(central_particle, r, color='green', fill=False, linestyle='--', label=f'$r = {r}$')
circle_outer = plt.Circle(central_particle, r + dr, color='orange', fill=False, linestyle='--', label=f'$r + \Delta r = {r + dr}$')
ax.add_artist(circle_inner)
ax.add_artist(circle_outer)

# Annotate particles within the shell
for pos in particle_positions:
    distance = np.linalg.norm(np.array(pos) - np.array(central_particle))
    if r <= distance < r + dr:
        ax.plot(*pos, 'o', color='purple', markersize=7)

# Add labels and legend
ax.set_title('Radial Distribution Function Demonstration')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.legend()

# Show plot
plt.show()
Ignoring fixed x limits to fulfill fixed data aspect with adjustable data limits.
<Figure size 800x800 with 1 Axes>

The radial distribution function, denoted by g(r)g(r), is a measure of the “structure” of a fluid or solid. It quantifies the average number of particles at a distance rr from a central particle relative to the ideal gas case. The radial distribution function is defined as

g(r)=N(r)4πr2Δrρg(r) = \frac{\langle N(r) \rangle}{4 \pi r^2 \Delta r \rho}

where N(r)\langle N(r) \rangle is the average number of particles in a shell of radius rr and thickness Δr\Delta r around a central particle, ρ\rho is the number density of particles, and rr is the distance from the central particle. The radial distribution function provides information about the local structure of a system, such as the presence of short-range order, long-range order, or disorder.

Computing the Radial Distribution Function

To compute the radial distribution function for a given configuration of particles, we need to follow these steps:

  1. Compute the distance between all pairs of particles in the system.
  2. Bin the distances into radial bins.
  3. Compute the radial distribution function using the formula given above.

Let’s consider an example to illustrate how to compute the radial distribution function for a simple system of particles.

Example: Lennard-Jones Fluid

Lennard-Jones Fluid

Consider a three-dimensional Lennard-Jones fluid (ϵ=1\epsilon = 1, σ=1\sigma = 1, rc=2.5r_c = 2.5) with periodic boundary conditions and N=1500N = 1500 particles in a cubic box of side length L=20L = 20. The final configuration of the particles, after a molecular dynamics simulation at a temperature of T=0.5T = 0.5, is given in the file lj.xyz. We want to compute the radial distribution function for this configuration of particles.

Let’s start by loading the configuration of particles from the file lj.xyz.

# Load the configuration of particles from the file lj.xyz
import numpy as np

def read_xyz(filename):
    with open(filename, 'r') as f:
        lines = f.readlines()
        n_atoms = int(lines[0])
        data = np.zeros((n_atoms, 3))
        for i in range(2, n_atoms + 2):
            # Skip the first and second columns (id and species) and extract the x, y, z coordinates
            data[i - 2] = np.array([float(x) for x in lines[i].split()[2:5]])
    return data

# Load the configuration of particles from the file lj.xyz
filename = 'lj.xyz'
positions = read_xyz(filename)
n_atoms = len(positions)
print(f'Number of particles: {n_atoms}')
Number of particles: 1500

Now that we have loaded the configuration of particles, we can compute the distance between all pairs of particles in the system.

# Compute the distance between all pairs of particles in the system
def compute_distance(positions, box_length):
    n_atoms = len(positions)
    distances = []
    for i in range(n_atoms):
        for j in range(n_atoms):
            if i >= j:
                continue
            dr = positions[i] - positions[j]
            dr = dr - box_length * np.round(dr / box_length)
            distance = np.linalg.norm(dr)
            distances.append(distance)
    return distances

# Compute the distance between all pairs of particles in the system
box_length = 20
distances = compute_distance(positions, box_length)
print(f'Number of distances: {len(distances)}')

# Compute the statistics of the distances
print(f'Minimum distance: {min(distances)}')
print(f'Maximum distance: {max(distances)}')
print(f'Mean distance: {np.mean(distances)}')
print(f'Standard deviation of distance: {np.std(distances)}')
Number of distances: 1124250
Minimum distance: 0.9396699365202659
Maximum distance: 17.240959235657396
Mean distance: 9.489382004836067
Standard deviation of distance: 3.0900910059195477

Next, we need to bin the distances into radial bins. We will use a bin width of 0.2 and a maximum distance of 10 for this example.

# Bin the distances into radial bins
def bin_distances(distances, bin_width, max_distance):
    bins = np.arange(bin_width, max_distance + bin_width * 2, bin_width)
    hist, _ = np.histogram(distances, bins=bins)
    return hist

# Bin the distances into radial bins
bin_width = 0.01
max_distance = 3.5
hist = bin_distances(distances, bin_width, max_distance)
print(f'Number of bins: {len(hist)}')
Number of bins: 350

Finally, we can compute the radial distribution function using the formula given above.

def radial_distribution_function(hist, n_atoms, box_length, bin_width):
    rho = n_atoms / box_length**3
    r = (np.arange(1, len(hist) + 1) - 0.5) * bin_width
    shell_volumes = 4 * np.pi * r**2 * bin_width
    g = hist / (rho * n_atoms * shell_volumes)
    return r, g

# Compute the radial distribution function
r, g = radial_distribution_function(hist, n_atoms, box_length, bin_width)

Now that we have computed the radial distribution function, we can plot it to visualize the structure of the fluid.

Source
import matplotlib.pyplot as plt

# Plot the radial distribution function
plt.plot(r, g)
plt.xlabel('r')
plt.ylabel('g(r)')
plt.title('Radial Distribution Function')

# Annotate the first peak
first_peak_label = 'First Nearest Neighbor Peak'
plt.annotate(
    first_peak_label,
    xy=(1.2, 3),
    xytext=(1.5, 4),
    arrowprops=dict(facecolor='black', shrink=0.05))

# Annotate the second peak
second_peak_label = 'Second Nearest Neighbor Peak'
plt.annotate(
    second_peak_label,
    xy=(2.2, 1.5),
    xytext=(2.5, 2),
    arrowprops=dict(facecolor='black', shrink=0.05))

# Annotate the third peak
third_peak_label = 'Third Nearest Neighbor Peak'
plt.annotate(
    third_peak_label,
    xy=(3.2, 1),
    xytext=(3.5, 1.5),
    arrowprops=dict(facecolor='black', shrink=0.05))

# Horizontal line at g(r) = 1
plt.axhline(y=1, color='r', linestyle='--')
plt.text(0, 1.1, 'Ideal Gas', color='red')

plt.show()
<Figure size 640x480 with 1 Axes>

Interpretation of the Radial Distribution Function

The radial distribution function provides information about the local structure of a system. Here are some key points to keep in mind when interpreting the radial distribution function:

Summary

In this lecture, we introduced the radial distribution function as a measure of the local structure of a fluid or solid. We discussed how to compute the radial distribution function for a given configuration of particles and interpret the results. The radial distribution function provides valuable insights into the interactions between particles in a system and will helps us understand the structure and properties of bead-spring polymers in the next lecture and project.