Introduction to structuretoolkit

Introduction to structuretoolkit#

structuretoolkit is a Python library for building, analysing and visualising atomistic structures for materials science. It does not introduce a new structure class of its own – instead it operates directly on the ase.atoms.Atoms object from the Atomic Simulation Environment (ASE). You can therefore think of structuretoolkit as a large collection of additional, ASE-compatible functions: anything you can build with ase.build, you can analyse and post-process with structuretoolkit, and the result is again a plain ase.atoms.Atoms object that you can keep using with ASE (writing files, running an ASE calculator, etc.).

How the package is organised#

The functionality is grouped into four submodules:

  • structuretoolkit.build – construct new structures (compounds, surfaces, grain boundaries, SQS, meshes, …)

  • structuretoolkit.analyse – analyse existing structures (neighbors, symmetry, strain, descriptors, …)

  • structuretoolkit.common – small helper functions and converters to other representations (pymatgen, phonopy, pyscal)

  • structuretoolkit.visualize – 3d visualisation of structures

All of these functions are also re-exported from the root of the package (e.g. structuretoolkit.get_neighbors) for backwards compatibility, but it is recommended to import functions through their submodule (structuretoolkit.analyse.get_neighbors) because that is where they are documented and grouped by purpose. Throughout this notebook we therefore import structuretoolkit as stk and always call functions as stk.analyse.xxx(), stk.build.xxx(), stk.common.xxx() or stk.visualize.xxx().

This notebook walks through (almost) every public function of structuretoolkit, grouped by submodule, and shows how it is meant to be combined with ASE.

import warnings

warnings.filterwarnings(
    "ignore"
)  # silence third-party deprecation warnings for this demo

import numpy as np
from ase.build import bulk

import structuretoolkit as stk

print("structuretoolkit version:", stk.__version__)

structuretoolkit version: 0.0.46.dev25+g04a02d8ac.d20260626

structuretoolkit.common – helpers and converters#

The common submodule contains small utility functions that are used throughout the rest of the package, plus converters that translate an ase.atoms.Atoms object into the structure representations used by other materials-science packages (pymatgen, phonopy, pyscal). They are handy on their own whenever you need a quick property of a structure without writing the boilerplate yourself.

Basic structure properties#

get_number_species_atoms() counts how many atoms of each chemical element are present, and select_index() returns the indices of all atoms of a given element – both operate directly on the Atoms object returned by ase.build.bulk().

structure = bulk("Fe", cubic=True).repeat(2) + bulk("Al", cubic=True).repeat((2, 1, 1))

print("Number of atoms per species:", stk.common.get_number_species_atoms(structure))
print(
    "Indices of the Al atoms:",
    stk.common.select_index(structure=structure, element="Al"),
)

Number of atoms per species: {'Al': 8, 'Fe': 16}
Indices of the Al atoms: [16 17 18 19 20 21 22 23]

Cell geometry#

  • get_cell() turns a float, a 3-vector, a 3x3 matrix or an Atoms object into a normalised (3, 3) cell array - this is what most other functions use internally to accept flexible cell input.

  • get_vertical_length() returns, for every cell vector, the height of the box measured perpendicular to the opposite face (the relevant length scale e.g. for slab thickness or periodic cutoffs).

  • get_extended_positions() repeats the atoms across the periodic boundary so that a search for neighbors close to the cell boundary does not miss any periodic image.

  • get_wrapped_coordinates() wraps arbitrary Cartesian coordinates back into the periodic cell.

  • center_coordinates_in_unit_cell() re-wraps all atoms of a structure in place, e.g. after coordinates were shifted outside of [0, 1) in fractional units.

fe_bulk = bulk("Fe", cubic=True)

print("Normalised (3, 3) cell from a lattice constant:\n", stk.common.get_cell(2.87))
print(
    "\nVertical heights of the Fe cell:",
    stk.common.get_vertical_length(structure=fe_bulk),
)

extended_positions, original_indices = stk.common.get_extended_positions(
    structure=fe_bulk, width=3.0, return_indices=True
)
print(
    "\nAtoms after adding a 3 Angstrom buffer layer around the cell:",
    len(extended_positions),
)

shifted = fe_bulk.copy()
shifted.positions += 5 * np.diag(shifted.cell)  # move atoms several cells away
stk.common.center_coordinates_in_unit_cell(structure=shifted)
print("\nPositions after wrapping back into the cell:\n", shifted.positions)

Normalised (3, 3) cell from a lattice constant:
 [[2.87 0.   0.  ]
 [0.   2.87 0.  ]
 [0.   0.   2.87]]

Vertical heights of the Fe cell: [2.87 2.87 2.87]

Atoms after adding a 3 Angstrom buffer layer around the cell: 91

Positions after wrapping back into the cell:
 [[-6.68890263e-16 -6.68890263e-16 -6.68890263e-16]
 [ 1.43500000e+00  1.43500000e+00  1.43500000e+00]]

Applying a homogeneous strain#

apply_strain() deforms the cell (and, by default, the atoms with it) according to a strain matrix epsilon. With mode="linear" the deformation gradient is F = epsilon + 1; with mode="lagrangian" epsilon is interpreted as the Lagrangian strain tensor (F^T F - 1) / 2, which must be symmetric. Passing return_box=True leaves the original structure untouched and returns a strained copy instead.

structure = bulk("Fe", cubic=True)
strained = stk.common.apply_strain(structure=structure, epsilon=0.01, return_box=True)

print("Original volume:  ", structure.get_volume())
print("Strained volume:  ", strained.get_volume())

Original volume:   23.63990300000001
Strained volume:   24.356215700803006

Converting to other structure representations#

Several materials-science packages have their own structure classes. structuretoolkit provides light-weight, round-trip converters so you never have to leave the ASE Atoms representation for long:

  • ase_to_pymatgen() / pymatgen_to_ase() / pymatgen_read_from_file() – convert to/from pymatgen.core.Structure, e.g. to use pymatgen’s symmetry or surface generation tools.

  • ase_to_pyscal() – convert to a pyscal3.core.System, which several of the structuretoolkit.analyse functions use internally for structural fingerprinting.

  • atoms_to_phonopy() / phonopy_to_atoms() – convert to/from phonopy.structure.atoms.PhonopyAtoms, e.g. as input for a phonon calculation with phonopy.

fe_bulk = bulk("Fe", cubic=True)

pmg_structure = stk.common.ase_to_pymatgen(structure=fe_bulk)
print(pmg_structure)

back_to_ase = stk.common.pymatgen_to_ase(structure=pmg_structure)
print("\nRound-tripped through pymatgen:", back_to_ase)

pyscal_system = stk.common.ase_to_pyscal(structure=fe_bulk)
print("\npyscal system:", type(pyscal_system))

phonopy_atoms = stk.common.atoms_to_phonopy(atom=fe_bulk)
print("\nphonopy atoms:", type(phonopy_atoms))
print(
    "Round-tripped through phonopy:",
    stk.common.phonopy_to_atoms(ph_atoms=phonopy_atoms),
)

Full Formula (Fe2)
Reduced Formula: Fe
abc   :   2.870000   2.870000   2.870000
angles:  90.000000  90.000000  90.000000
pbc   :       True       True       True
Sites (2)
  #  SP      a    b    c    magmom
---  ----  ---  ---  ---  --------
  0  Fe    0    0    0         2.3
  1  Fe    0.5  0.5  0.5       2.3

Round-tripped through pymatgen: Atoms(symbols='Fe2', pbc=True, cell=[2.87, 2.87, 2.87], initial_magmoms=...)

pyscal system: <class 'pyscal3.core.System'>

phonopy atoms: <class 'phonopy.structure.atoms.PhonopyAtoms'>
Round-tripped through phonopy: Atoms(symbols='Fe2', pbc=True, cell=[2.87, 2.87, 2.87])

Finally, structuretoolkit.common.SymmetryError is the exception raised by the symmetry-related functions in structuretoolkit.analyse (see below) whenever spglib cannot determine the symmetry of a structure, for example because the cell or the atomic positions are degenerate.

structuretoolkit.analyse – analysing existing structures#

This is the largest submodule. It contains everything needed to characterise the local and global geometry of a structure: neighbor lists, symmetry operations, strain, and a number of structural fingerprints that are commonly used to classify crystal structures or identify defects.

Neighbor lists: get_neighbors() and get_neighborhood()#

get_neighbors() is one of the most used functions in the package. It returns a Neighbors object which lazily exposes distances, vecs (vectors), indices and shells (the indices of coordination shells, found by clustering distances) for every atom of the structure, fully periodic-boundary-aware. You can ask for a fixed number of neighbors (num_neighbors) and/or a cutoff_radius.

structure = bulk("Fe", cubic=True)  # bcc iron
neigh = stk.analyse.get_neighbors(structure=structure, num_neighbors=8)

print("Distance to the 8 nearest neighbors of atom 0:\n", neigh.distances[0])
print("\nCoordination shell of each of those neighbors:\n", neigh.shells[0])
print("\nChemical symbols of the neighbors of atom 0:\n", neigh.chemical_symbols[0])

Distance to the 8 nearest neighbors of atom 0:
 [2.48549291 2.48549291 2.48549291 2.48549291 2.48549291 2.48549291
 2.48549291 2.48549291]

Coordination shell of each of those neighbors:
 [1 1 1 1 1 1 1 1]

Chemical symbols of the neighbors of atom 0:
 ['Fe' 'Fe' 'Fe' 'Fe' 'Fe' 'Fe' 'Fe' 'Fe']

get_neighborhood() answers the same question for arbitrary points in space (not necessarily atomic positions), which is useful to probe e.g. interstitial sites or grid points.

random_point_in_cell = np.dot(np.random.random(3), structure.cell)
neigh_point = stk.analyse.get_neighborhood(
    structure=structure, positions=random_point_in_cell, num_neighbors=4
)
print(
    "Distances from a random point in the cell to its 4 nearest atoms:\n",
    neigh_point.distances,
)

Distances from a random point in the cell to its 4 nearest atoms:
 [1.40751277 1.54226709 1.70976867 1.95412966]

Minimum-image distances: get_distances_array() and find_mic()#

get_distances_array() returns the full pairwise distance (or vector) matrix between two sets of positions, automatically applying the minimum image convention for periodic directions. find_mic() is the lower-level function that maps arbitrary displacement vectors onto their periodic minimum image – it is what get_distances_array() and get_neighbors() use internally.

structure = bulk("Al", cubic=True)
distance_matrix = stk.analyse.get_distances_array(structure=structure)
print("Pairwise distance matrix:\n", distance_matrix)

displacement = structure.positions[0] - structure.positions[1] + 10 * structure.cell[0]
mic_distance = stk.analyse.find_mic(structure=structure, v=displacement, vectors=False)
print("\nDisplacement folded back via the minimum image convention:", mic_distance)

Pairwise distance matrix:
 [[0.         2.86378246 2.86378246 2.86378246]
 [2.86378246 0.         2.86378246 2.86378246]
 [2.86378246 2.86378246 0.         2.86378246]
 [2.86378246 2.86378246 2.86378246 0.        ]]

Displacement folded back via the minimum image convention: 2.8637824638055176

Box symmetry#

get_symmetry() wraps spglib and returns a Symmetry object with all rotation/translation pairs that leave the structure invariant. It is the basis for several convenience functions:

  • get_spacegroup() – the international space-group symbol/number.

  • get_symmetry_dataset() – the full spglib symmetry dataset (Wyckoff positions, point group, …).

  • get_primitive_cell() – reduce a (super)cell to its primitive cell.

  • symmetrize_vectors() / group_points_by_symmetry() / get_equivalent_points() – use the box symmetry to symmetrize per-atom vectors, or to find groups of geometrically equivalent points.

  • get_ir_reciprocal_mesh() – the irreducible set of points of a reciprocal-space mesh (e.g. for a DFT k-point grid).

structure = bulk("Al", cubic=True).repeat(2)

sym = stk.analyse.get_symmetry(structure=structure)
print("Number of symmetry operations:", len(sym.rotations))
print("Space group:", stk.analyse.get_spacegroup(structure=structure))

primitive = stk.analyse.get_primitive_cell(structure=structure)
print(
    "\nAtoms in the supercell vs. its primitive cell:",
    len(structure),
    "->",
    len(primitive),
)

# Forces/displacements on a perfect crystal must symmetrize to (numerically) zero
random_vectors = np.random.random(structure.positions.shape)
print("\nNorm of a random vector field before/after symmetrization:")
print(
    np.linalg.norm(random_vectors),
    "->",
    np.linalg.norm(stk.analyse.symmetrize_vectors(structure, random_vectors)),
)

mapping, grid_points = stk.analyse.get_ir_reciprocal_mesh(
    structure=structure, mesh=[4, 4, 4]
)
print(
    "\nIrreducible points of a 4x4x4 k-point mesh:",
    len(np.unique(mapping)),
    "out of",
    len(grid_points),
)

Number of symmetry operations: 1536
Space group: {'InternationalTableSymbol': 'Fm-3m', 'Number': 225}

Atoms in the supercell vs. its primitive cell: 32 -> 1

Norm of a random vector field before/after symmetrization:
5.660616829549077 -> 8.161356608876601e-17

Irreducible points of a 4x4x4 k-point mesh: 10 out of 64

get_equivalent_atoms() answers a related, but independent, question: which atoms of the structure map onto each other under the box symmetry. It is implemented via phonopy’s PhonopyAtoms and spglib, rather than through the Symmetry class above.

structure = bulk("Al", cubic=True).repeat(2)
print("Equivalent-atom labels:", stk.analyse.get_equivalent_atoms(structure=structure))

Equivalent-atom labels: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]

Layers, tessellations and clustering#

  • get_layers() groups atoms into layers along the cell directions (or arbitrary planes), which is useful for analysing slabs, stacking sequences or counting the number of atomic layers.

  • get_voronoi_vertices(), get_voronoi_neighbors() and get_delaunay_neighbors() expose the Voronoi/Delaunay tessellation of the (periodic) atom positions.

  • get_cluster_positions() merges positions that lie within a given distance of each other (DBSCAN clustering), e.g. to average out thermal noise or symmetric duplicates.

structure = bulk("Al", cubic=True).repeat(3)

layers = stk.analyse.get_layers(structure=structure)
print("Number of (111)-independent layers along x, y, z:", layers.max(axis=0) + 1)

voronoi_vertices = stk.analyse.get_voronoi_vertices(structure=bulk("Al", cubic=True))
print(
    "\nVoronoi vertices of an fcc unit cell (octahedral + tetrahedral sites):",
    len(voronoi_vertices),
)

voronoi_pairs = stk.analyse.get_voronoi_neighbors(structure=structure)
delaunay_pairs = stk.analyse.get_delaunay_neighbors(structure=structure)
print(
    "Voronoi neighbor pairs:",
    len(voronoi_pairs),
    " Delaunay neighbor pairs:",
    len(delaunay_pairs),
)

doubled_positions = np.append(structure.positions, structure.positions, axis=0)
clustered = stk.analyse.get_cluster_positions(
    structure=structure, positions=doubled_positions
)
print(
    "\nDuplicated positions clustered back down to:",
    len(clustered),
    "(original count:",
    len(structure),
    ")",
)

Number of (111)-independent layers along x, y, z: [6 6 6]

Voronoi vertices of an fcc unit cell (octahedral + tetrahedral sites): 12
Voronoi neighbor pairs: 846  Delaunay neighbor pairs: 1051

Duplicated positions clustered back down to: 108 (original count: 108 )

get_average_of_unique_labels() is the small helper used internally by get_layers() and get_cluster_positions(): given an array of (possibly repeated) integer labels, it returns the average value of the elements that share a label.

labels = [0, 1, 0, 2]
values = [0, 1, 2, 3]
print(stk.analyse.get_average_of_unique_labels(labels=labels, values=values))

[1. 1. 3.]

Interstitial sites: get_interstitials()#

get_interstitials() returns an Interstitials object which automatically searches for interstitial sites of a given coordination number, e.g. tetrahedral (num_neighbors=4) or octahedral (num_neighbors=6) sites in a bcc lattice. Besides the positions, it exposes diagnostic quantities like get_distances(), get_volumes() and get_steinhardt_parameters() of the candidate sites – useful to judge how “clean” / symmetric the identified sites are.

bcc_fe = bulk("Fe", cubic=True)

octahedral = stk.analyse.get_interstitials(structure=bcc_fe, num_neighbors=6)
tetrahedral = stk.analyse.get_interstitials(structure=bcc_fe, num_neighbors=4)

print("Octahedral interstitial sites in bcc Fe:", len(octahedral.positions))
print("Tetrahedral interstitial sites in bcc Fe:", len(tetrahedral.positions))
print(
    "\nDistance of the octahedral sites to their nearest neighbor atoms:",
    octahedral.get_distances(),
)

Octahedral interstitial sites in bcc Fe: 6
Tetrahedral interstitial sites in bcc Fe: 12

Distance of the octahedral sites to their nearest neighbor atoms: [1.435 1.435 1.435 1.435 1.435 1.435]

Local strain: get_strain()#

get_strain() compares a (e.g. deformed, defective, or thermally displaced) structure against an ideal reference structure of the same crystal type, and returns the Lagrangian strain tensor of every atom based on its local neighborhood. This is a common way to visualise the strain field around a dislocation, crack tip or grain boundary.

reference = bulk("Fe", cubic=True)
deformed = stk.common.apply_strain(structure=reference, epsilon=0.02, return_box=True)

strain = stk.analyse.get_strain(structure=deformed, ref_structure=reference)
print("Strain tensor shape:", strain.shape)
print("Per-atom isotropic strain (trace / 3):", np.trace(strain, axis1=1, axis2=2) / 3)

Strain tensor shape: (2, 3, 3)
Per-atom isotropic strain (trace / 3): [0.0202 0.0202]

Structural fingerprints (powered by pyscal)#

A group of functions builds on pyscal to classify the local crystal structure around each atom – very useful to identify the crystal structure of an unknown configuration, or to spot defects (grain boundaries, stacking faults, melted regions, …) in an otherwise crystalline structure:

  • get_adaptive_cna_descriptors() – common neighbor analysis (fcc/hcp/bcc/icosahedral/other).

  • get_centro_symmetry_descriptors() – centrosymmetry parameter (large for atoms near a defect).

  • get_diamond_structure_descriptors() – identifies cubic/hexagonal diamond environments.

  • get_steinhardt_parameters() – rotationally invariant bond-orientational order parameters.

  • get_voronoi_volumes() – the Voronoi cell volume of every atom.

  • find_solids() – counts how many atoms are “solid” vs. “liquid”-like, e.g. for melting-point calculations.

bcc_fe = bulk("Fe", cubic=True)
diamond_si = bulk("Si", crystalstructure="diamond", a=5.43, cubic=True)

print(
    "Common neighbor analysis of bcc Fe:",
    stk.analyse.get_adaptive_cna_descriptors(structure=bcc_fe),
)
print(
    "Diamond-structure analysis of Si:  ",
    stk.analyse.get_diamond_structure_descriptors(structure=diamond_si),
)

centro = stk.analyse.get_centro_symmetry_descriptors(structure=bcc_fe, num_neighbors=8)
print("\nCentrosymmetry parameter of perfect bcc Fe (should be ~0):", centro)

q, cluster_id = stk.analyse.get_steinhardt_parameters(structure=bcc_fe.repeat(3))
print("\nSteinhardt parameters q4, q6 shape:", q.shape)

print(
    "\nVoronoi cell volume per atom:", stk.analyse.get_voronoi_volumes(structure=bcc_fe)
)
print(
    "Number of solid-like atoms (out of",
    len(bcc_fe),
    "):",
    stk.analyse.find_solids(structure=bcc_fe),
)

Common neighbor analysis of bcc Fe: {'others': 0, 'fcc': 0, 'hcp': 0, 'bcc': 2, 'ico': 0}
Diamond-structure analysis of Si:   {'others': 0, 'cubic diamond': 8, 'cubic diamond 1NN': 0, 'cubic diamond 2NN': 0, 'hex diamond': 0, 'hex diamond 1NN': 0, 'hex diamond 2NN': 0}

Centrosymmetry parameter of perfect bcc Fe (should be ~0): [2.36658272e-30 0.00000000e+00]

Steinhardt parameters q4, q6 shape: (2, 432)

Voronoi cell volume per atom: [11.8199515 11.8199515]
Number of solid-like atoms (out of 2 ): 2.0

Machine-learning structure descriptors#

Two further analysis tools generate per-atom descriptor vectors that are commonly used as input features for machine-learned interatomic potentials:

  • soap_descriptor_per_atom() computes the Smooth Overlap of Atomic Positions (SOAP) descriptor via the optional dscribe dependency (pip install structuretoolkit[dscribe]).

  • get_snap_descriptors_per_atom(), get_snap_descriptor_derivatives() and get_snap_descriptor_names() compute the bispectrum components used by the SNAP machine-learning potential, via an embedded LAMMPS instance.

cu_fcc = bulk("Cu", "fcc", a=3.6, cubic=True)

soap = stk.analyse.soap_descriptor_per_atom(
    structure=cu_fcc, r_cut=6.0, n_max=8, l_max=6
)
print("SOAP descriptor shape (n_atoms, n_features):", soap.shape)

SOAP descriptor shape (n_atoms, n_features): (4, 252)

ni_fcc = bulk("Ni", cubic=True)

snap_names = stk.analyse.get_snap_descriptor_names(twojmax=6)
print("Number of SNAP bispectrum components:", len(snap_names))

# The exact keyword arguments accepted by lammps.create_atoms() differ between LAMMPS
# python-module releases, so this call is most reliable when structuretoolkit and your
# LAMMPS installation are from a matching, recent release.
try:
    snap_descriptors = stk.analyse.get_snap_descriptors_per_atom(
        structure=ni_fcc, atom_types=["Ni"]
    )
    print("SNAP descriptor shape (n_atoms, n_components):", snap_descriptors.shape)
except TypeError as error:
    print("Skipping execution - incompatible LAMMPS python module installed:", error)

Number of SNAP bispectrum components: 30
SNAP descriptor shape (n_atoms, n_components): (4, 30)

structuretoolkit.build – constructing new structures#

While ase.build already provides bulk crystals, surfaces and molecules, structuretoolkit.build adds generators for several structure types that come up frequently in materials science but are not part of ASE itself: ordered intermetallic compounds, high-index stepped/kinked surfaces, grain boundaries, special quasirandom structures (SQS) for disordered alloys, regular real-space meshes, and a client for the Materials Project database.

Ordered compound prototypes: B2(), C14(), C15(), C36(), D03()#

These functions build the unit cells of common binary intermetallic phases directly from two chemical symbols (optionally a custom lattice constant), instead of having to look up Wyckoff positions and space groups by hand. The result is, as always, a plain ase.atoms.Atoms object.

b2 = stk.build.B2("Fe", "Al")  # CsCl-type FeAl
c15 = stk.build.C15("Mg", "Cu")  # cubic Laves phase MgCu2

print("B2 FeAl: ", b2.get_chemical_formula(), "cell =", b2.cell.lengths())
print("C15 MgCu2:", c15.get_chemical_formula(), "cell =", c15.cell.lengths())

B2 FeAl:  AlFe cell = [2.87 2.87 2.87]
C15 MgCu2: Cu16Mg8 cell = [7.38598547 7.38598547 7.38598547]

High-index stepped & kinked surfaces#

get_high_index_surface_info() derives the Miller index of a high-index surface that exposes a given terrace, step and kink orientation (following the microfacet notation of Van Hove & Somorjai), and high_index_surface() builds the corresponding relaxed slab (using pymatgen’s SpacegroupAnalyzer internally to refine the structure).

miller_index, kink_dir, step_dir = stk.build.get_high_index_surface_info(
    element="Ni",
    crystal_structure="fcc",
    lattice_constant=3.526,
    terrace_orientation=[1, 1, 1],
    step_orientation=[1, 1, 0],
    kink_orientation=[1, 0, 1],
    step_down_vector=[1, 1, 0],
    length_step=2,
    length_terrace=3,
    length_kink=1,
)
print("High-index Miller index for this terrace/step/kink combination:", miller_index)

slab = stk.build.high_index_surface(
    element="Ni",
    crystal_structure="fcc",
    lattice_constant=3.526,
    terrace_orientation=[1, 1, 1],
    step_orientation=[1, 1, 0],
    kink_orientation=[1, 0, 1],
    step_down_vector=[1, 1, 0],
    length_step=2,
    length_terrace=3,
    length_kink=1,
    layers=20,
    vacuum=10,
)
print("Resulting stepped/kinked Ni slab:", slab)

High-index Miller index for this terrace/step/kink combination: [-9 -3 -5]
Resulting stepped/kinked Ni slab: Atoms(symbols='Ni20', pbc=True, cell=[[5.575095514876862, 0.0, 0.0], [-2.7875477574384036, 5.978623880124941, 0.0], [1.6071788640578362e-15, 2.522639151650182e-15, 26.247222712325165]])

Grain boundaries#

get_grainboundary_info() lists the geometrically possible grain boundaries (by CSL sigma value) for a given rotation axis, using the optional aimsgb dependency (pip install structuretoolkit[grainboundary]). Once you picked a sigma value and plane, grainboundary() builds the actual bicrystal as an ase.atoms.Atoms object.

gb_options = stk.build.get_grainboundary_info(axis=[0, 0, 1], max_sigma=5)
print("Available grain boundaries for the [0 0 1] axis up to sigma=5:")
print(list(gb_options.keys()))

al_bulk = bulk("Al", cubic=True)
gb_structure = stk.build.grainboundary(
    axis=[0, 0, 1], sigma=5, plane=[1, 2, 0], initial_struct=al_bulk
)
print(
    "\nNumber of atoms in the sigma=5 [0 0 1](1 2 0) grain boundary:", len(gb_structure)
)

Available grain boundaries for the [0 0 1] axis up to sigma=5:
[np.int64(5)]

Number of atoms in the sigma=5 [0 0 1](1 2 0) grain boundary: 40

Regular real-space meshes: create_mesh()#

create_mesh() builds a regular grid of points spanning a cell (either an Atoms object or a raw cell) – for example to sample a potential energy surface, or as the starting set of candidate positions in a custom interstitial search.

structure = bulk("Al", cubic=True)
mesh = stk.create_mesh(cell=structure, n_mesh=4)
print("Mesh array shape (3 components, nx, ny, nz):", mesh.shape)

Mesh array shape (3 components, nx, ny, nz): (3, 4, 4, 4)

Special quasirandom structures (SQS)#

sqs_structures() generates Special Quasirandom Structures – supercells whose short-range order parameters best approximate those of a fully random alloy with a given composition – via the optional sqsgenerator dependency. This is the standard way to set up a disordered-alloy supercell for DFT calculations. The function returns a SqsResultPack; use .best() to get the single best result, which exposes both the resulting Atoms object (.atoms()) and its short-range order parameters (.sro()).

au_bulk = bulk("Au", cubic=True)

best_sqs = stk.build.sqs_structures(
    structure=au_bulk,
    composition={"Cu": 16, "Au": 16},
    supercell=(2, 2, 2),
).best()

sqs_structure = best_sqs.atoms()
print("SQS composition:", stk.common.get_number_species_atoms(sqs_structure))
print(
    "Short-range order parameter array shape (n_shells, n_species, n_species):",
    best_sqs.sro().shape,
)

SQS composition: {'Cu': 16, 'Au': 16}
Short-range order parameter array shape (n_shells, n_species, n_species): (5, 2, 2)

Materials Project database access#

materialsproject_search() and materialsproject_by_id() wrap the Materials Project REST API (via the optional mp-api dependency, pip install structuretoolkit[mp-api]) and directly return ase.atoms.Atoms objects. Both functions need a (free) Materials Project API key, either passed as api_key=... or exported as the MP_API_KEY environment variable.

import os

if os.environ.get("MP_API_KEY"):
    hits = list(stk.build.materialsproject_search("Fe", is_stable=True))
    print("Stable Fe-containing entries found:", len(hits))
    fe_ground_state = stk.build.materialsproject_by_id("mp-13")
    print(fe_ground_state)
else:
    print(
        "Set the MP_API_KEY environment variable (or pass api_key=...) to query the Materials Project,"
        " e.g.:\n\n"
        "    structures = list(stk.build.materialsproject_search('Fe', is_stable=True))\n"
        "    fe_ground_state = stk.build.materialsproject_by_id('mp-13')"
    )

Set the MP_API_KEY environment variable (or pass api_key=...) to query the Materials Project, e.g.:

    structures = list(stk.build.materialsproject_search('Fe', is_stable=True))
    fe_ground_state = stk.build.materialsproject_by_id('mp-13')

structuretoolkit.visualize – looking at structures#

plot3d() renders any ase.atoms.Atoms object in 3d. The default mode="NGLview" produces an interactive nglview widget – the nicest choice for interactive work in a live Jupyter session, so it is what you will want to use day-to-day. For this notebook we instead use mode="plotly" and render the resulting figure to a static PNG (which needs the optional kaleido package, pip install kaleido), so the plot is also visible in a statically rendered copy of this notebook, e.g. on GitHub.

structure = bulk("Al", cubic=True).repeat(2)
stk.visualize.plot3d(structure)

plot3d() accepts a scalar_field to color atoms by a per-atom property – this is a convenient way to visualise the output of the analysis functions above. As an example, we color a small Al cell by the centrosymmetry parameter after displacing one atom: the displaced atom (and, to a lesser extent, its neighbors) lights up against the otherwise perfectly centrosymmetric lattice.

displaced = bulk("Al", cubic=True)
displaced.positions[0] += [0.4, 0.0, 0.0]

centro = stk.analyse.get_centro_symmetry_descriptors(
    structure=displaced, num_neighbors=12
)
stk.visualize.plot3d(displaced, scalar_field=centro)

Where to go from here#

  • Browse the source code – every function has a docstring with the full list of arguments.

  • Have a look at the test suite, which contains many more focused, runnable usage examples for each function shown here.

  • structuretoolkit is also used as the structure-analysis backend of pyiron_atomistics, where the same functions are available as methods directly on the structure object (e.g. structure.analyse.get_neighbors()).

  • If a function is missing optional dependencies (pyscal3, spglib, aimsgb, dscribe, mp-api, nglview, …), install the corresponding extra, e.g. pip install structuretoolkit[pyscal,symmetry].