Coordination number — and its smarter cousins#
The plain coordination number (CN) — how many neighbours does this atom have? — is a useful first descriptor but a blunt one. It is an integer, it depends on a cutoff radius, and it doesn’t know anything about the quality of those neighbours.
pyscal exposes three smarter variants on top of the plain CN:
descriptor |
what it does |
reference |
|---|---|---|
|
integer count |
— |
|
distance-weighted, continuous |
Hoppe (1979) |
|
weights neighbours by their CN |
Calle-Vallejo et al. (2015) |
|
atoms per unit volume from neighbour distances |
— |
Below we walk through each one with intuition-building examples.
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. Bulk crystals — the four numbers#
Compute all four descriptors for FCC, BCC and HCP. Watch how BCC’s effective coordination is below 14 — that’s the algorithm spotting the two distinct neighbour shells.
rows = []
for name, lc in [('FCC', 4.05), ('BCC', 2.87), ('HCP', 3.21)]:
atoms = make_crystal(name.lower(), lattice_constant=lc, repetitions=(4,4,4))
pyscal.find_neighbors(atoms, method='cutoff', cutoff=0)
cn = pyscal.coordination_number(atoms)
econ = pyscal.effective_coordination_number(atoms)
gcn = pyscal.generalized_coordination_number(atoms)
rho = pyscal.local_density(atoms)
rows.append((name, int(cn[0]), float(econ[0]), float(gcn[0]), float(rho[0])))
hdr_density = 'density(1/A^3)'
print(f"{'':5s} {'CN':>5} {'ECoN':>8} {'GCN':>8} {hdr_density:>16}")
for r in rows:
print(f'{r[0]:5s} {r[1]:>5d} {r[2]:>8.3f} {r[3]:>8.3f} {r[4]:>16.6f}')
CN ECoN GCN density(1/A^3)
FCC 12 12.000 12.000 0.121975
BCC 14 11.276 14.000 0.179541
HCP 12 12.000 12.000 0.086612
Observations
FCC and HCP both have CN = 12, but slightly different densities.
BCC has 8 first-shell + 6 second-shell neighbours = CN 14, but the weaker second shell drops ECoN below 14.
GCN equals CN for any atom in a bulk crystal (everyone is fully coordinated).
2. Effective coordination number (ECoN)#
Hoppe’s ECoN replaces the binary in/out cutoff by a smooth weight that decays for far neighbours. For BCC the first shell sits at \(r_1 \approx 0.866 a\), the second at \(r_2 = a\) — close together. Look at the per-atom weight contributions to see how ECoN “feels” both shells without double-counting.
bcc = make_crystal('bcc', lattice_constant=2.87, repetitions=(4, 4, 4))
pyscal.find_neighbors(bcc, method='cutoff', cutoff=0)
cn = int(pyscal.coordination_number(bcc)[0])
econ = float(pyscal.effective_coordination_number(bcc)[0])
print(f'BCC: CN (binary) = {cn}')
print(f'BCC: ECoN (Hoppe) = {econ:.3f}')
print('-> 2nd-shell neighbours contribute, but with reduced weight.')
BCC: CN (binary) = 14
BCC: ECoN (Hoppe) = 11.276
-> 2nd-shell neighbours contribute, but with reduced weight.
3. Generalized coordination number (GCN)#
GCN looks at how well coordinated each of an atom’s neighbours is. In a perfect bulk this matches CN, but near a defect the values diverge — which is exactly where GCN earns its name.
Remove a single atom from FCC (a vacancy). The 12 atoms that used to surround it now have CN = 11. Their neighbours, in turn, have lost a piece of their environment too — and GCN captures this second-order effect.
atoms = make_crystal('fcc', lattice_constant=4.05, repetitions=(5, 5, 5))
# place vacancy at an interior site
centre = atoms.get_positions().mean(axis=0)
d2 = ((atoms.get_positions() - centre) ** 2).sum(axis=1)
vac_idx = int(np.argmin(d2))
vac_pos = atoms.get_positions()[vac_idx].copy()
del atoms[vac_idx]
pyscal.find_neighbors(atoms, method='cutoff', cutoff=3.5)
cn = pyscal.coordination_number(atoms)
gcn = pyscal.generalized_coordination_number(atoms)
# distance of each remaining atom to the vacancy site
r_to_vac = np.linalg.norm(atoms.get_positions() - vac_pos, axis=1)
fig = go.Figure()
fig.add_trace(go.Scatter(x=r_to_vac, y=cn, mode='markers', name='CN',
marker=dict(size=6, opacity=0.6)))
fig.add_trace(go.Scatter(x=r_to_vac, y=gcn, mode='markers', name='GCN',
marker=dict(size=6, opacity=0.6)))
fig.update_layout(xaxis_title='distance from vacancy (\u00c5)',
yaxis_title='coordination',
title='CN and GCN around an FCC vacancy',
height=420, template='simple_white')
fig
Atoms far from the vacancy sit at CN = GCN = 12 (perfect bulk). The 12 first-shell neighbours drop to CN = 11. GCN drops further than CN, because second-shell atoms also lose a coordinated neighbour. That compounding is exactly what catalysis studies exploit at step edges and surface defects.
4. Local density#
local_density estimates atoms per unit volume from neighbour distances
around each atom. Useful for spotting density variations in
non-uniform systems (interfaces, voids, amorphous regions).
perfect = make_crystal('fcc', lattice_constant=4.05, repetitions=(4, 4, 4))
pyscal.find_neighbors(perfect, method='cutoff', cutoff=0)
rho_p = pyscal.local_density(perfect)
rng = np.random.default_rng(0)
noisy = make_crystal('fcc', lattice_constant=4.05, repetitions=(4, 4, 4))
noisy.positions += rng.normal(scale=0.20, size=noisy.positions.shape)
pyscal.find_neighbors(noisy, method='cutoff', cutoff=0)
rho_n = pyscal.local_density(noisy)
fig = go.Figure()
fig.add_trace(go.Histogram(x=rho_p, nbinsx=40, name='perfect FCC',
opacity=0.7, histnorm='probability density'))
fig.add_trace(go.Histogram(x=rho_n, nbinsx=40, name='noisy FCC (\u03c3=0.2 \u00c5)',
opacity=0.7, histnorm='probability density'))
fig.update_layout(barmode='overlay', xaxis_title='local density (1/\u00c5\u00b3)',
yaxis_title='density', height=380, template='simple_white',
title='Per-atom local density distribution')
fig
5. Where everything is stored#
Each call writes into atoms.arrays under a pyscal_ prefix. Mix and
match descriptors freely; they do not overwrite each other.
atoms = make_crystal('fcc', lattice_constant=4.05, repetitions=(3, 3, 3))
pyscal.find_neighbors(atoms, method='cutoff', cutoff=0)
for fn in [pyscal.coordination_number, pyscal.effective_coordination_number,
pyscal.generalized_coordination_number, pyscal.local_density]:
fn(atoms)
for k in sorted(atoms.arrays):
if k.startswith('pyscal_'):
v = atoms.arrays[k]
print(f' {k:40s} shape={v.shape} first={v[0]}')
pyscal_cn shape=(108,) first=12
pyscal_cutoff shape=(108,) first=3.43653895656662
pyscal_diff shape=(108, 12, 3) first=[[ 0. 2.025 2.025]
[ 2.025 0. 2.025]
[ 2.025 2.025 0. ]
[ 2.025 0. -2.025]
[ 2.025 -2.025 0. ]
[-2.025 2.025 0. ]
[ 0. 2.025 -2.025]
[ 0. -2.025 2.025]
[-2.025 0. 2.025]
[ 0. -2.025 -2.025]
[-2.025 0. -2.025]
[-2.025 -2.025 0. ]]
pyscal_econ shape=(108,) first=11.99999999999998
pyscal_gcn shape=(108,) first=12.0
pyscal_local_density shape=(108,) first=0.1219754869558647
pyscal_neighbordist shape=(108, 12) first=[2.86378246 2.86378246 2.86378246 2.86378246 2.86378246 2.86378246
2.86378246 2.86378246 2.86378246 2.86378246 2.86378246 2.86378246]
pyscal_neighbors shape=(108, 12) first=[98 81 35 9 11 27 26 74 73 2 1 3]
pyscal_neighborweight shape=(108, 12) first=[1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
pyscal_phi shape=(108, 12) first=[ 1.57079633 0. 0.78539816 0. -0.78539816 2.35619449
1.57079633 -1.57079633 3.14159265 -1.57079633 3.14159265 -2.35619449]
pyscal_r shape=(108, 12) first=[2.86378246 2.86378246 2.86378246 2.86378246 2.86378246 2.86378246
2.86378246 2.86378246 2.86378246 2.86378246 2.86378246 2.86378246]
pyscal_temp_neighbordist shape=(108, 42) first=[2.86378246 2.86378246 2.86378246 2.86378246 2.86378246 2.86378246
2.86378246 2.86378246 2.86378246 2.86378246 2.86378246 2.86378246
4.05 4.05 4.05 4.05 4.05 4.05
4.96021673 4.96021673 4.96021673 4.96021673 4.96021673 4.96021673
4.96021673 4.96021673 4.96021673 4.96021673 4.96021673 4.96021673
4.96021673 4.96021673 4.96021673 4.96021673 4.96021673 4.96021673
4.96021673 4.96021673 4.96021673 4.96021673 4.96021673 4.96021673]
pyscal_temp_neighbors shape=(108, 42) first=[ 98 81 35 9 11 27 26 74 73 2 1 3 72 24 8 36 12 4
105 106 34 33 83 82 107 99 75 97 25 10 21 47 63 30 102 93
85 78 71 39 13 6]
pyscal_theta shape=(108, 12) first=[0.78539816 0.78539816 1.57079633 2.35619449 1.57079633 1.57079633
2.35619449 0.78539816 0.78539816 2.35619449 2.35619449 1.57079633]
Take-aways#
Use CN for a quick integer count.
Use ECoN when neighbour shells overlap or when continuity matters (defect characterisation, MD).
Use GCN for surfaces and catalysis — it’s a known descriptor of activity for many transition-metal facets.
Use local_density for non-uniform systems.
References#
R. Hoppe, Effective coordination numbers (ECoN), Z. Kristallogr. 150, 23 (1979).
F. Calle-Vallejo et al., Finding optimal surface sites on heterogeneous catalysts by counting nearest neighbors, Science 350, 185 (2015).