Angular and bond-length distribution functions#

Two of the simplest, most informative “global” descriptors of a structure are:

ADF — the histogram of bond angles \(\theta_{jik}\) formed by every pair of neighbours around every atom.
BLDF — the histogram of bond lengths \(r_{ij}\) for the current neighbour list.

They are pre-Steinhardt, pre-machine-learning, but still excellent for spotting phase changes, melting, glass formation, or for sanity-checking a freshly built structure.

import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
import pyscal3 as pyscal
from pyscal3.structures import make_crystal

1. ADF for the three close-packed lattices#

Perfect FCC has bond angles at exactly 60°, 90°, 120° and 180°. BCC peaks at 70.5° and 109.5° (the tetrahedral pair) once the second shell is included; HCP shares FCC’s main angles but with a slightly different weight at 60° (extra in-plane bonds).

structures = {
    'FCC': (make_crystal('fcc', lattice_constant=4.05, repetitions=(4,4,4)), 3.2),
    'BCC': (make_crystal('bcc', lattice_constant=2.87, repetitions=(4,4,4)), 3.0),
    'HCP': (make_crystal('hcp', lattice_constant=3.21, repetitions=(4,4,4)), 3.5),
}

fig = go.Figure()
for name, (atoms, cut) in structures.items():
    pyscal.find_neighbors(atoms, method='cutoff', cutoff=cut)
    adf, angles = pyscal.angular_distribution_function(atoms, bins=180)
    fig.add_trace(go.Scatter(x=angles, y=adf, mode='lines', name=name))

for ang, txt in [(60, '60\u00b0'), (90, '90\u00b0'),
                 (109.5, '109.5\u00b0'), (120, '120\u00b0'), (180, '180\u00b0')]:
    fig.add_vline(x=ang, line_dash='dot', line_color='lightgrey',
                  annotation_text=txt, annotation_position='top')
fig.update_layout(xaxis_title='angle (degrees)', yaxis_title='probability density',
                  title='Angular distribution function for FCC, BCC, HCP',
                  template='simple_white', height=420)
fig

Each lattice has a unique fingerprint of peaks. Use the ADF as a quick visual classifier without needing any spherical-harmonic machinery.

2. Bond-length distribution#

BLDF is just a histogram of distances in the neighbour list. With a small cutoff you isolate the first shell; widen the cutoff to bring in the next shell and watch new peaks appear.

fcc = make_crystal('fcc', lattice_constant=4.05, repetitions=(5,5,5))
fig = go.Figure()
for cut in [3.2, 4.0, 4.6, 5.6]:
    a = fcc.copy()
    pyscal.find_neighbors(a, method='cutoff', cutoff=cut)
    bldf, r = pyscal.bond_length_distribution(a, bins=120, rmin=2.4, rmax=5.7)
    fig.add_trace(go.Scatter(x=r, y=bldf, mode='lines', name=f'cutoff = {cut} \u00c5'))
fig.update_layout(xaxis_title='r (\u00c5)', yaxis_title='probability density',
                  title='BLDF for FCC at increasing cutoff (more shells appear)',
                  template='simple_white', height=400)
fig

Each spike marks a neighbour shell: \(a/\sqrt 2 \approx 2.86\) Å (1st), \(a \approx 4.05\) Å (2nd), \(a\sqrt{3/2} \approx 4.96\) Å (3rd), and so on.

3. Crystal vs. liquid#

Heat the lattice (large random displacements) and the sharp peaks turn into broad humps — the hallmark of a disordered phase. ADF and BLDF make this picture obvious without doing any further analysis.

rng = np.random.default_rng(0)

perfect = make_crystal('fcc', lattice_constant=4.05, repetitions=(5,5,5))
liquid = perfect.copy()
liquid.positions += rng.normal(scale=0.45, size=liquid.positions.shape)

for atoms in (perfect, liquid):
    pyscal.find_neighbors(atoms, method='cutoff', cutoff=3.5)

fig = make_subplots(rows=1, cols=2, subplot_titles=('ADF', 'BLDF'))
for label, atoms in [('crystal', perfect), ('liquid-like', liquid)]:
    adf, ang = pyscal.angular_distribution_function(atoms, bins=180)
    bldf, r  = pyscal.bond_length_distribution(atoms, bins=120, rmin=2.0, rmax=4.0)
    fig.add_trace(go.Scatter(x=ang, y=adf, mode='lines', name=label,
                             legendgroup=label, showlegend=True), row=1, col=1)
    fig.add_trace(go.Scatter(x=r, y=bldf, mode='lines', name=label,
                             legendgroup=label, showlegend=False), row=1, col=2)
fig.update_xaxes(title_text='angle (deg)', row=1, col=1)
fig.update_xaxes(title_text='r (\u00c5)', row=1, col=2)
fig.update_layout(template='simple_white', height=400,
                  title='Crystal vs. liquid-like configuration')
fig

4. What is stored#

Both functions write their last result into atoms.info:

pyscal_adf, pyscal_adf_angles
pyscal_bldf, pyscal_bldf_r

atoms = make_crystal('bcc', lattice_constant=2.87, repetitions=(4,4,4))
pyscal.find_neighbors(atoms, method='cutoff', cutoff=3.0)
pyscal.angular_distribution_function(atoms)
pyscal.bond_length_distribution(atoms)
for k in ('pyscal_adf', 'pyscal_adf_angles', 'pyscal_bldf', 'pyscal_bldf_r'):
    v = atoms.info[k]
    print(f'  {k:20s} shape={np.asarray(v).shape}')

  pyscal_adf           shape=(180,)
  pyscal_adf_angles    shape=(180,)
  pyscal_bldf          shape=(100,)
  pyscal_bldf_r        shape=(100,)

Take-aways#

ADF and BLDF are simple, robust descriptors that often carry more signal than a single scalar order parameter.
They are the natural starting point when inspecting a new structure or trajectory for the first time.
Pair them with Steinhardt or Minkowski \(q_l\) for a more complete characterisation.