Motif Normal Form¶

This guide explains how the Motif Normal Form (MNF), the atomic position coordinates in a CNF, is constructed.

Multiple Equivalent Representations¶

Given a crystal's atomic positions, there are many ways to write them down. The origin can be shifted to any atom, the lattice may have symmetries that permute equivalent positions, and atoms of the same element can be listed in any order. The MNF algorithm finds a unique canonical representation by exploring all these choices and selecting the lexicographically smallest.

In [1]:

from pymatgen.core import Structure, Lattice
from cnf import CNFConstructor
from cnf.lattice import Superbasis
from cnf.lattice.selling import SuperbasisSellingReducer
import numpy as np

# Create rutile TiO2
a, c = 4.594, 2.959
lattice = Lattice.tetragonal(a, c)
species = ["Ti", "Ti", "O", "O", "O", "O"]
coords = [[0, 0, 0], [0.5, 0.5, 0.5], [0.305, 0.305, 0], [0.695, 0.695, 0], [0.805, 0.195, 0.5], [0.195, 0.805, 0.5]]
rutile = Structure(lattice, species, coords)

print(f"Rutile TiO2: {len(rutile)} atoms")
for i, site in enumerate(rutile):
    print(f"  {site.specie}: {[f'{x:.3f}' for x in site.frac_coords]}")

from pymatgen.core import Structure, Lattice from cnf import CNFConstructor from cnf.lattice import Superbasis from cnf.lattice.selling import SuperbasisSellingReducer import numpy as np # Create rutile TiO2 a, c = 4.594, 2.959 lattice = Lattice.tetragonal(a, c) species = ["Ti", "Ti", "O", "O", "O", "O"] coords = [[0, 0, 0], [0.5, 0.5, 0.5], [0.305, 0.305, 0], [0.695, 0.695, 0], [0.805, 0.195, 0.5], [0.195, 0.805, 0.5]] rutile = Structure(lattice, species, coords) print(f"Rutile TiO2: {len(rutile)} atoms") for i, site in enumerate(rutile): print(f" {site.specie}: {[f'{x:.3f}' for x in site.frac_coords]}")

Rutile TiO2: 6 atoms
  Ti: ['0.000', '0.000', '0.000']
  Ti: ['0.500', '0.500', '0.500']
  O: ['0.305', '0.305', '0.000']
  O: ['0.695', '0.695', '0.000']
  O: ['0.805', '0.195', '0.500']
  O: ['0.195', '0.805', '0.500']

Stabilizers from the LNF¶

The construction begins with the stabilizers from the LNF. These are lattice permutations that leave the canonical vonorms unchanged. Each stabilizer permutation can correspond to one or more unimodular matrices that transform fractional coordinates while preserving the lattice geometry.

In [2]:

# Get the reduced superbasis and its vonorms
superbasis = Superbasis.from_pymatgen_structure(rutile)
reducer = SuperbasisSellingReducer()
result = reducer.reduce(superbasis)
reduced_superbasis = result.reduced_object
reduced_vonorms = reduced_superbasis.compute_vonorms()

# Find the canonical vonorms
perms = reduced_vonorms.permissible_perms
candidates = []
for perm in perms:
    permuted = reduced_vonorms.apply_permutation(tuple(perm.vonorm_permutation))
    candidates.append((permuted.tuple, perm))
candidates.sort(key=lambda x: x[0])
canonical_vonorms = candidates[0][0]

# Find stabilizers - permutations that give the canonical vonorms
stabilizers = []
for vonorms_tuple, perm in candidates:
    if vonorms_tuple == canonical_vonorms:
        stabilizers.append(perm)

print(f"Number of stabilizer permutations: {len(stabilizers)}")

# Get the reduced superbasis and its vonorms superbasis = Superbasis.from_pymatgen_structure(rutile) reducer = SuperbasisSellingReducer() result = reducer.reduce(superbasis) reduced_superbasis = result.reduced_object reduced_vonorms = reduced_superbasis.compute_vonorms() # Find the canonical vonorms perms = reduced_vonorms.permissible_perms candidates = [] for perm in perms: permuted = reduced_vonorms.apply_permutation(tuple(perm.vonorm_permutation)) candidates.append((permuted.tuple, perm)) candidates.sort(key=lambda x: x[0]) canonical_vonorms = candidates[0][0] # Find stabilizers - permutations that give the canonical vonorms stabilizers = [] for vonorms_tuple, perm in candidates: if vonorms_tuple == canonical_vonorms: stabilizers.append(perm) print(f"Number of stabilizer permutations: {len(stabilizers)}")

Number of stabilizer permutations: 1

In [3]:

# Each stabilizer permutation may correspond to multiple unimodular matrices
print("Unimodular matrices for each stabilizer permutation:")
all_matrices = []
for i, stab in enumerate(stabilizers):
    matrices = stab.all_matrices
    all_matrices.extend(matrices)
    print(f"  Permutation {i+1}: {len(matrices)} matrix(es)")

print(f"\nTotal unimodular matrices: {len(all_matrices)}")

# Each stabilizer permutation may correspond to multiple unimodular matrices print("Unimodular matrices for each stabilizer permutation:") all_matrices = [] for i, stab in enumerate(stabilizers): matrices = stab.all_matrices all_matrices.extend(matrices) print(f" Permutation {i+1}: {len(matrices)} matrix(es)") print(f"\nTotal unimodular matrices: {len(all_matrices)}")

Unimodular matrices for each stabilizer permutation:
  Permutation 1: 4 matrix(es)

Total unimodular matrices: 4

In [4]:

# Show a few example matrices
print("Example unimodular matrices:")
for i, mat in enumerate(all_matrices[:3]):
    print(f"\n  Matrix {i+1}:")
    for row in mat.matrix:
        print(f"    {list(row)}")

# Show a few example matrices print("Example unimodular matrices:") for i, mat in enumerate(all_matrices[:3]): print(f"\n Matrix {i+1}:") for row in mat.matrix: print(f" {list(row)}")

Example unimodular matrices:

  Matrix 1:
    [0, 0, -1]
    [0, -1, 0]
    [-1, 0, 0]

  Matrix 2:
    [0, 0, 1]
    [0, 1, 0]
    [-1, 0, 0]

  Matrix 3:
    [0, 0, 1]
    [0, -1, 0]
    [1, 0, 0]

Applying Stabilizer Transformations¶

Each unimodular matrix transforms the fractional coordinates. We apply all of them to generate candidate motif representations.

In [5]:

from cnf.motif import FractionalMotif

# Get the original motif
motif = FractionalMotif.from_pymatgen_structure(rutile)

print("Original motif:")
for atom, pos in zip(motif.atoms, motif.positions):
    print(f"  {atom}: {[f'{x:.3f}' for x in pos]}")

from cnf.motif import FractionalMotif # Get the original motif motif = FractionalMotif.from_pymatgen_structure(rutile) print("Original motif:") for atom, pos in zip(motif.atoms, motif.positions): print(f" {atom}: {[f'{x:.3f}' for x in pos]}")

Original motif:
  Ti: ['0.000', '0.000', '0.000']
  Ti: ['0.500', '0.500', '0.500']
  O: ['0.305', '0.305', '0.000']
  O: ['0.695', '0.695', '0.000']
  O: ['0.805', '0.195', '0.500']
  O: ['0.195', '0.805', '0.500']

In [6]:

# Apply a stabilizer transformation
if len(all_matrices) > 1:
    example_matrix = all_matrices[1]
else:
    example_matrix = all_matrices[0]
    
transformed_motif = motif.apply_unimodular(example_matrix)

print(f"After applying unimodular matrix:")
print(f"Matrix: {example_matrix.matrix.tolist()}")
print()
for atom, pos in zip(transformed_motif.atoms, transformed_motif.positions):
    wrapped = [p % 1.0 for p in pos]
    print(f"  {atom}: {[f'{x:.3f}' for x in wrapped]}")

# Apply a stabilizer transformation if len(all_matrices) > 1: example_matrix = all_matrices[1] else: example_matrix = all_matrices[0] transformed_motif = motif.apply_unimodular(example_matrix) print(f"After applying unimodular matrix:") print(f"Matrix: {example_matrix.matrix.tolist()}") print() for atom, pos in zip(transformed_motif.atoms, transformed_motif.positions): wrapped = [p % 1.0 for p in pos] print(f" {atom}: {[f'{x:.3f}' for x in wrapped]}")

After applying unimodular matrix:
Matrix: [[0, 0, 1], [0, 1, 0], [-1, 0, 0]]

  Ti: ['0.000', '0.000', '0.000']
  Ti: ['0.500', '0.500', '0.500']
  O: ['0.000', '0.305', '0.305']
  O: ['0.000', '0.695', '0.695']
  O: ['0.500', '0.195', '0.805']
  O: ['0.500', '0.805', '0.195']

Origin Shifts¶

For each transformed motif, we try placing different atoms at the origin. This eliminates translational freedom. By convention, the implementation considers atoms of a specific element type (determined by sorting order) as origin candidates.

In [7]:

# The sorted elements determine which atoms can be origin candidates
print(f"Sorted elements: {motif.sorted_elements}")
print(f"Origin candidates are atoms of element: {motif.sorted_elements[0]}")
print(f"Number of origin candidates: {motif.num_origin_atoms}")

# The sorted elements determine which atoms can be origin candidates print(f"Sorted elements: {motif.sorted_elements}") print(f"Origin candidates are atoms of element: {motif.sorted_elements[0]}") print(f"Number of origin candidates: {motif.num_origin_atoms}")

Sorted elements: [Element Ti, Element O]
Origin candidates are atoms of element: Ti
Number of origin candidates: 2

In [8]:

# Show origin shifts for the first element type
origin_element = motif.sorted_elements[0]
origin_positions = motif.get_element_positions(origin_element)

print(f"Shifting origin to each {origin_element} atom:")
for i, origin_pos in enumerate(origin_positions):
    print(f"\nOrigin at {origin_element} atom {i+1} (shift by {[f'{-x:.3f}' for x in origin_pos]}):")
    for atom, pos in zip(motif.atoms, motif.positions):
        new_pos = [(p - s) % 1.0 for p, s in zip(pos, origin_pos)]
        marker = " ← origin" if np.allclose(pos, origin_pos) else ""
        print(f"  {atom}: {[f'{x:.3f}' for x in new_pos]}{marker}")

# Show origin shifts for the first element type origin_element = motif.sorted_elements[0] origin_positions = motif.get_element_positions(origin_element) print(f"Shifting origin to each {origin_element} atom:") for i, origin_pos in enumerate(origin_positions): print(f"\nOrigin at {origin_element} atom {i+1} (shift by {[f'{-x:.3f}' for x in origin_pos]}):") for atom, pos in zip(motif.atoms, motif.positions): new_pos = [(p - s) % 1.0 for p, s in zip(pos, origin_pos)] marker = " ← origin" if np.allclose(pos, origin_pos) else "" print(f" {atom}: {[f'{x:.3f}' for x in new_pos]}{marker}")

Shifting origin to each Ti atom:

Origin at Ti atom 1 (shift by ['-0.000', '-0.000', '-0.000']):
  Ti: ['0.000', '0.000', '0.000'] ← origin
  Ti: ['0.500', '0.500', '0.500']
  O: ['0.305', '0.305', '0.000']
  O: ['0.695', '0.695', '0.000']
  O: ['0.805', '0.195', '0.500']
  O: ['0.195', '0.805', '0.500']

Origin at Ti atom 2 (shift by ['-0.500', '-0.500', '-0.500']):
  Ti: ['0.500', '0.500', '0.500']
  Ti: ['0.000', '0.000', '0.000'] ← origin
  O: ['0.805', '0.805', '0.500']
  O: ['0.195', '0.195', '0.500']
  O: ['0.305', '0.695', '0.000']
  O: ['0.695', '0.305', '0.000']

Sorting and Lexicographic Selection¶

For each combination of stabilizer matrix and origin shift, atoms are sorted by element and then by position. The coordinates are then flattened to a tuple, excluding the origin atom which is implicitly at (0, 0, 0). Finally, the lexicographically smallest tuple across all combinations is selected as the canonical MNF.

Discretization¶

The fractional coordinates are discretized to integer multiples of 1/δ (delta). For example, with δ=100, a coordinate of 0.305 becomes round(0.305 × 100) = 31.

In [9]:

# Get the actual MNF from CNF construction
constructor = CNFConstructor(xi=1.0, delta=100)
cnf_result = constructor.from_pymatgen_structure(rutile)
mnf = cnf_result.cnf.motif_normal_form

print("Final Motif Normal Form:")
print(f"  Elements: {mnf.elements}")
print(f"  Coordinates: {mnf.coord_list}")
print(f"  Delta: {mnf.delta}")
print()
print(f"The first atom ({mnf.elements[0]}) is at the origin [0, 0, 0]")
print(f"Remaining {len(mnf.elements)-1} atoms × 3 coordinates = {len(mnf.coord_list)} integers")

# Get the actual MNF from CNF construction constructor = CNFConstructor(xi=1.0, delta=100) cnf_result = constructor.from_pymatgen_structure(rutile) mnf = cnf_result.cnf.motif_normal_form print("Final Motif Normal Form:") print(f" Elements: {mnf.elements}") print(f" Coordinates: {mnf.coord_list}") print(f" Delta: {mnf.delta}") print() print(f"The first atom ({mnf.elements[0]}) is at the origin [0, 0, 0]") print(f"Remaining {len(mnf.elements)-1} atoms × 3 coordinates = {len(mnf.coord_list)} integers")

Final Motif Normal Form:
  Elements: ['Ti', 'Ti', 'O', 'O', 'O', 'O']
  Coordinates: (50, 50, 50, 0, 30, 30, 0, 70, 70, 50, 20, 80, 50, 80, 20)
  Delta: 100

The first atom (Ti) is at the origin [0, 0, 0]
Remaining 5 atoms × 3 coordinates = 15 integers

In [10]:

# Show the mapping from integers back to fractional coordinates
delta = mnf.delta
print(f"Discretization with delta = {delta}:")
print()
print("MNF integers → fractional coordinates:")
coords = mnf.coord_list
for i in range(0, len(coords), 3):
    atom_idx = i // 3 + 1  # +1 because first atom is at origin
    int_coords = coords[i:i+3]
    frac_coords = [c / delta for c in int_coords]
    print(f"  Atom {atom_idx}: {list(int_coords)} → {[f'{x:.3f}' for x in frac_coords]}")

# Show the mapping from integers back to fractional coordinates delta = mnf.delta print(f"Discretization with delta = {delta}:") print() print("MNF integers → fractional coordinates:") coords = mnf.coord_list for i in range(0, len(coords), 3): atom_idx = i // 3 + 1 # +1 because first atom is at origin int_coords = coords[i:i+3] frac_coords = [c / delta for c in int_coords] print(f" Atom {atom_idx}: {list(int_coords)} → {[f'{x:.3f}' for x in frac_coords]}")

Discretization with delta = 100:

MNF integers → fractional coordinates:
  Atom 1: [50, 50, 50] → ['0.500', '0.500', '0.500']
  Atom 2: [0, 30, 30] → ['0.000', '0.300', '0.300']
  Atom 3: [0, 70, 70] → ['0.000', '0.700', '0.700']
  Atom 4: [50, 20, 80] → ['0.500', '0.200', '0.800']
  Atom 5: [50, 80, 20] → ['0.500', '0.800', '0.200']

Summary¶

The MNF construction proceeds as follows:

1. Obtain stabilizer matrices from the LNF (permutations that preserve canonical vonorms)
2. Apply each unimodular matrix to the atomic coordinates
3. For each transformed motif, shift the origin to each candidate atom
4. Sort atoms by element and position
5. Discretize coordinates to multiples of 1/δ
6. Select the lexicographically smallest coordinate tuple

The result is 3×(N-1) integers that uniquely characterize the atomic positions.

The Complete CNF¶

Combining LNF (7 integers) and MNF (3×(N-1) integers) gives the complete Crystal Normal Form.

In [11]:

cnf = cnf_result.cnf

print("Complete Crystal Normal Form:")
print(f"\n  CNF = {cnf.coords}")
print(f"\n  Breakdown:")
print(f"    LNF (7 vonorms):        {cnf.coords[:7]}")
print(f"    MNF ({len(cnf.coords)-7} coordinates): {cnf.coords[7:]}")
print(f"\n  Parameters:")
print(f"    xi = {cnf.xi} Ų")
print(f"    delta = {cnf.delta}")
print(f"    elements = {cnf.elements}")

cnf = cnf_result.cnf print("Complete Crystal Normal Form:") print(f"\n CNF = {cnf.coords}") print(f"\n Breakdown:") print(f" LNF (7 vonorms): {cnf.coords[:7]}") print(f" MNF ({len(cnf.coords)-7} coordinates): {cnf.coords[7:]}") print(f"\n Parameters:") print(f" xi = {cnf.xi} Ų") print(f" delta = {cnf.delta}") print(f" elements = {cnf.elements}")

Complete Crystal Normal Form:

  CNF = (9, 21, 21, 51, 30, 30, 42, 50, 50, 50, 0, 30, 30, 0, 70, 70, 50, 20, 80, 50, 80, 20)

  Breakdown:
    LNF (7 vonorms):        (9, 21, 21, 51, 30, 30, 42)
    MNF (15 coordinates): (50, 50, 50, 0, 30, 30, 0, 70, 70, 50, 20, 80, 50, 80, 20)

  Parameters:
    xi = 1.0 Ų
    delta = 100
    elements = ['Ti', 'Ti', 'O', 'O', 'O', 'O']