Wigner–Seitz defect analysis#

Given a current configuration and a reference (perfect-lattice) configuration, the Wigner–Seitz analysis assigns every atom to its nearest reference site. The site-occupancy histogram then tells you the defect content directly:

occupancy = 1 → normal lattice atom
occupancy = 0 → vacancy
occupancy ≥ 2 → interstitial(s)
different type on a single site → antisite

It’s the simplest, most robust defect detector for crystalline systems and works in arbitrary dimension.

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

1. A perfect crystal — nothing to find#

If current == reference, every site has exactly one atom.

ref = make_crystal('fcc', lattice_constant=4.05, repetitions=(4, 4, 4))
result = pyscal.wigner_seitz_analysis(ref, ref)
print(f"vacancies   : {result['vacancy_count']}")
print(f"interstitials: {result['interstitial_count']}")
print(f"unique site occupancies: {set(result['occupancy'].tolist())}")

vacancies   : 0
interstitials: 0
unique site occupancies: {1}

2. Single vacancy#

Delete one atom and re-run. The vacancy count should be 1, and the vacancy index points at the missing reference site.

current = ref.copy()
del current[100]
result = pyscal.wigner_seitz_analysis(current, ref)
print(f"vacancies     : {result['vacancy_count']}")
print(f"interstitials : {result['interstitial_count']}")
print(f"vacancy index : {result['vacancy_indices']}")

vacancies     : 1
interstitials : 0
vacancy index : [100]

3. Frenkel pair: vacancy + interstitial#

Move one atom from its lattice site to a non-lattice position. The analysis sees one empty site (vacancy) and one site with two atoms (the interstitial sits on top of its nearest neighbour).

current = ref.copy()
pos = current.get_positions()
pos[100] = pos[100] + np.array([2.0, 2.0, 0.0])  # kick atom 100 off site
current.set_positions(pos)
result = pyscal.wigner_seitz_analysis(current, ref)
print(f"vacancies     : {result['vacancy_count']}")
print(f"interstitials : {result['interstitial_count']}")
print(f"vacancy idx   : {result['vacancy_indices']}")
print(f"interstitial site idx : {result['interstitial_sites']}")

vacancies     : 1
interstitials : 1
vacancy idx   : [100]
interstitial site idx : [103]

4. Visualising defects in 3D#

Colour every atom by the occupancy of its assigned site. Then overlay vacancy positions in red. The picture makes the defect topology obvious.

rng = np.random.default_rng(0)
current = ref.copy()
# Five vacancies + two interstitials at random positions
vac_indices = rng.choice(len(current), size=5, replace=False)
del current[sorted(vac_indices.tolist(), reverse=True)]
extra_pos = ref.get_positions().mean(axis=0) + rng.normal(scale=2.0, size=(2, 3))
from ase import Atom
for p in extra_pos:
    current.append(Atom('Cu', position=p))

result = pyscal.wigner_seitz_analysis(current, ref)
occ_per_atom = result['occupancy'][result['site_index']]
vac_pos = ref.get_positions()[result['vacancy_indices']]

pos = current.get_positions()
fig = go.Figure()
fig.add_trace(go.Scatter3d(x=pos[:, 0], y=pos[:, 1], z=pos[:, 2],
                           mode='markers',
                           marker=dict(size=3, color=occ_per_atom,
                                       colorscale='Viridis',
                                       colorbar=dict(title='site occ.'),
                                       opacity=0.6),
                           name='atoms'))
fig.add_trace(go.Scatter3d(x=vac_pos[:, 0], y=vac_pos[:, 1], z=vac_pos[:, 2],
                           mode='markers',
                           marker=dict(size=8, color='red', symbol='x'),
                           name='vacancies'))
fig.update_layout(template='simple_white', height=560,
                  scene=dict(xaxis_title='x', yaxis_title='y', zaxis_title='z'),
                  title=f"WS analysis: {result['vacancy_count']} vacancies, "
                        f"{result['interstitial_count']} interstitials")
fig

5. Thermal expansion — use `affine_mapping='to_reference'`#

If the simulation cell has expanded thermally, raw distances no longer make sense — every atom looks slightly displaced. The affine_mapping option rescales the current cell to the reference cell first, removing this trivial offset.

current = ref.copy()
# Uniform thermal expansion of 2 %
current.set_cell(np.array(ref.cell) * 1.02, scale_atoms=True)

no_map  = pyscal.wigner_seitz_analysis(current, ref, affine_mapping='none')
with_map = pyscal.wigner_seitz_analysis(current, ref, affine_mapping='to_reference')
print(f"raw      : vacancies={no_map['vacancy_count']:3d}  "
      f"interstitials={no_map['interstitial_count']:3d}")
print(f"affine map: vacancies={with_map['vacancy_count']:3d}  "
      f"interstitials={with_map['interstitial_count']:3d}")

raw      : vacancies=  0  interstitials=  0
affine map: vacancies=  0  interstitials=  0

Without the affine map, expansion can shift atoms across Wigner–Seitz boundaries and produce spurious defects. With the mapping, the analysis correctly reports zero defects.

6. Antisites in a binary lattice#

Set per_type_occupancies=True and the result includes a per-element occupancy array. Comparing it to the reference types lets identify_defect_atoms flag antisites — right kind of site, wrong type of atom.

from ase.build import bulk
ref = bulk('NaCl', crystalstructure='rocksalt', a=5.64, cubic=True).repeat(3)
current = ref.copy()
syms = current.get_chemical_symbols()
# swap a pair of atoms to create two antisites
syms[0], syms[1] = syms[1], syms[0]
current.set_chemical_symbols(syms)

info = pyscal.identify_defect_atoms(current, ref)
print(info['defect_summary'])

2 antisites

7. Storage#

After the analysis, current.arrays carries the per-atom site assignment and current.info keeps the summary counts.

ref = make_crystal('fcc', lattice_constant=4.05, repetitions=(3, 3, 3))
current = ref.copy(); del current[5]
pyscal.wigner_seitz_analysis(current, ref)
print('arrays:')
for k in sorted(current.arrays):
    if k.startswith('pyscal_ws'):
        print(f'  {k:30s} shape={current.arrays[k].shape}')
print('info:')
for k in sorted(current.info):
    if k.startswith('pyscal_ws'):
        print(f'  {k:30s} = {current.info[k]}')

arrays:
  pyscal_ws_occupancy            shape=(107,)
  pyscal_ws_site_index           shape=(107,)
info:
  pyscal_ws_interstitial_count   = 0
  pyscal_ws_vacancy_count        = 1

Take-aways#

Wigner–Seitz analysis is a one-shot, parameter-free vacancy / interstitial detector — perfect for irradiation cascades, vacancy diffusion, MD post-processing.
Use affine_mapping='to_reference' whenever the cell can change shape (thermal expansion, NPT runs).
For multi-component systems, identify_defect_atoms adds antisite detection on top.