Angle-distribution function

The angle-distribution function describes the probability density of finding a triple of atoms with a specific angle. Note that this requires restricting the distance between atoms to a specific bond-length, or some other means to find bonds. Here we describe how to compute the angle distribution if bonds can be determined from a simple distance-based criterion, i.e. two atoms are considered bonded if their distance is smaller than $r_{\text{b}}$. The angle distribution is then

$$f_{\text{a}}(\varphi )=\frac{1}{N}\sum _{i=1}^N\frac{1}{N_i^2}\sum _{j=1}^N\sum _{k=1}^N\theta (r_{ij}-r_{\text{b}})\theta (r_{ik}-r_{\text{b}})\delta (\text{arccos}({\overrightarrow{r}}_{ij}\cdot {\overrightarrow{r}}_{ik}/(r_{ij}{r}_{ik}))-\varphi )$$

where $\theta (x)$ is the Heaviside step function, $N$ is the total number of atoms and $N_i=\sum _j\theta (r_{ij}-r_{\text{b}})$ is the number of neighbors (coordination number) of atom $i$. This angle-distribution function has the property $\int \text{d}\varphi f_{\text{a}}(\varphi )=1$, i.e. it is a probability density. It can be calculated from using the utility angle_distribution function.

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

from ase.io import read
from matscipy.angle_distribution import angle_distribution
from matscipy.neighbours import neighbour_list

# Load a disordered configuration (amorphous carbon)
a = read('../../tests/aC.cfg')

# Two atoms are bonded if their distance is smaller than this value
max_bond_distance = 2.0  # Angstroms

# Get distance between atoms up to a distance of 2 Angstroms which includes only nearest neighbors.
# The capital 'D' returns distance vectors, which are needed for angles.
i, j, d = neighbour_list('ijD', a, cutoff=max_bond_distance)

# Get compute angle distribution
nbins = 50
bin_spacing = np.pi/nbins
angle = (np.arange(nbins) + 0.5) * bin_spacing
nb_angles = angle_distribution(i, j, d, nbins)  # `angle_distribution` returns the number of triangles
probability = nb_angles / (np.sum(nb_angles) * bin_spacing)  # Normalize

# Plot the angle distribution
plt.plot(angle, probability)
plt.xlabel(r'Bond angle $\phi$ (radians)')
plt.ylabel(r'Probability density $f_\text{a}(\phi)$ (radians$^{-1}$)')
plt.xlim(0, np.pi)

(0.0, 3.141592653589793)