Example: Wurtzite Ferroelectric (ScAlN)

This example computes the out-of-plane cation displacement (D_z) in a wurtzite-type ferroelectric such as Sc_xAl_1-xN.

In wurtzite, each cation (Ga/Sc/Al) is tetrahedrally coordinated by 4 N atoms: 1 along the c-axis and 3 in the basal plane. The ferroelectric polarization is governed by the cation’s displacement relative to its N neighbors along z. The key challenge is that we only want the 3 in-plane N neighbors (not the axial one) to define the coordination center, so we need a custom neighbor list. For core concepts, see Introduction to ferrodispcalc.

Complete Script

import numpy as np
from ase.io import read
from ferrodispcalc.neighborlist import build_neighbor_list
from ferrodispcalc.compute import calculate_displacement

atoms = read("avg.xsf")

# --- Custom neighbor list for wurtzite ---
# Start with all 4 N neighbors per cation
nl_full = build_neighbor_list(atoms, ["Ga", "Sc"], ["N"], cutoff=4, neighbor_num=4)

center_id = nl_full[:, 0]
neighbor_id = nl_full[:, 1:]

# Filter: keep only the 3 neighbors in the basal plane (small dz)
box = atoms.get_cell().array
box_inv = np.linalg.inv(box)
center_pos = atoms.get_positions()[center_id - 1]
nl = np.full((nl_full.shape[0], 4), -1, dtype=int)  # 1 center + 3 neighbors
nl[:, 0] = center_id

for i, neighbors in enumerate(neighbor_id):
    diff = atoms.get_positions()[neighbors - 1] - center_pos[i]
    diff_frac = diff @ box_inv
    diff_frac[diff_frac > 0.5] -= 1.0
    diff_frac[diff_frac < -0.5] += 1.0
    diff = diff_frac @ box
    in_plane = np.abs(diff[:, 2]) < 1.0  # small z-component = basal plane
    nl[i, 1:] = neighbors[in_plane]

# --- Compute displacement ---
disp = calculate_displacement(atoms, nl)
mean_Dz = np.mean(disp[:, 2])
print(f"Mean Dz: {mean_Dz:.4f} Angstrom")

Step-by-Step Breakdown

Why a Custom Neighbor List?

build_neighbor_list returns the 4 nearest N neighbors for each cation. In wurtzite, these 4 neighbors split into two groups:

• 3 in the basal plane (roughly the same z as the cation)

• 1 along the c-axis (large dz)

For computing D_z, we want the displacement of the cation relative to the 3 in-plane neighbors only. So we filter by |dz| < 1 Angstrom after applying the minimum image convention:

nl_full = build_neighbor_list(atoms, ["Ga", "Sc"], ["N"], cutoff=4, neighbor_num=4)

# For each center, compute MIC-corrected dz to each neighbor
for i, neighbors in enumerate(neighbor_id):
    diff = atoms.get_positions()[neighbors - 1] - center_pos[i]
    diff_frac = diff @ box_inv
    diff_frac[diff_frac > 0.5] -= 1.0
    diff_frac[diff_frac < -0.5] += 1.0
    diff = diff_frac @ box
    in_plane = np.abs(diff[:, 2]) < 1.0
    nl[i, 1:] = neighbors[in_plane]

The resulting nl has shape (n_cations, 4) — column 0 is the center index, columns 1-3 are the 3 in-plane N neighbors.

Compute D_z

With the filtered neighbor list, calculate_displacement gives the cation offset from the center of its 3 basal-plane N neighbors. The z-component is D_z:

disp = calculate_displacement(atoms, nl)  # shape (n_cations, 3)
Dz = disp[:, 2]  # z-component only

See calculate_displacement() and build_neighbor_list() for API details.