Crystallizer: the immobile framework¶

The companion Crystallizer notebook reconstructs the mobile sodium sublattice of Na₃(Sb/W)S₄ from its diffusion density. Here we look at the other half of the crystallized structure: the immobile Sb/W–S host framework.

[Crystallizer.to_cif][gemdat.crystallizer.Crystallizer.to_cif] crystallizes the trajectory (mobile sites + time-averaged framework, fitted to the highest space group) and writes the result straight to a CIF. We then read that CIF back and pull out the immobile sites exactly the way we do for the two workshop reference CIFs — drop the sodium and keep the framework — so all three are handled through the same path.

In [1]:

import plotly.io as pio

# Render plotly figures as static HTML so they display on the docs website
pio.renderers.default = 'plotly_mimetype+notebook_connected'

from gemdat import Trajectory

# Same Na3(Sb/W)S4 AIMD run as the shape-analysis tutorial. The data lives in the
# `tests/data` submodule (run `git submodule update --init`); the paths below are
# relative to this notebook.
# `from_vasprun` loads the parse cache directly when it exists, so the (large)
# raw vasprun.xml is not actually needed here.
DATA = '../../tests/data/na3sbs4'
VASPRUN = f'{DATA}/vasprun.xml'

trajectory = Trajectory.from_vasprun(VASPRUN, cache=f'{DATA}/vasprun.xml.a45e320a.cache')
trajectory

import plotly.io as pio # Render plotly figures as static HTML so they display on the docs website pio.renderers.default = 'plotly_mimetype+notebook_connected' from gemdat import Trajectory # Same Na3(Sb/W)S4 AIMD run as the shape-analysis tutorial. The data lives in the # `tests/data` submodule (run `git submodule update --init`); the paths below are # relative to this notebook. # `from_vasprun` loads the parse cache directly when it exists, so the (large) # raw vasprun.xml is not actually needed here. DATA = '../../tests/data/na3sbs4' VASPRUN = f'{DATA}/vasprun.xml' trajectory = Trajectory.from_vasprun(VASPRUN, cache=f'{DATA}/vasprun.xml.a45e320a.cache') trajectory

Out[1]:

Full Formula (Na46 Sb14 W2 S64)
Reduced Formula: Na23Sb7WS32
abc   :  14.437000  14.437028  14.475216
angles:  90.000427  89.999751  90.237694
pbc   :       True       True       True
Constant lattice (True)
Sites (126)
Time steps (25000)

Crystallize to a CIF¶

We build the [Crystallizer][gemdat.crystallizer.Crystallizer] from the full trajectory (sodium is the mobile species; everything else — Sb, W, S — is the static framework) and call [to_cif][gemdat.crystallizer.Crystallizer.to_cif], which crystallizes and writes the symmetry-fitted structure to na3sbs4_crystallized.cif. The resolution and background_level are the same site-segmentation knobs used in the companion notebook.

In [2]:

from gemdat.crystallizer import Crystallizer

background_level = 0.25
CIF = 'na3sbs4_crystallized.cif'

cr = Crystallizer.from_trajectory(trajectory, floating_specie='Na', resolution=0.3)
result = cr.to_cif(CIF, background_level=background_level, snap_to_lower=True)

print(f'space group: {result.spacegroup_symbol} (#{result.spacegroup_number})')
print(f'symprec:     {result.symprec} Å')
print(f'wrote:       {CIF} ({len(result.structure)} sites)')

from gemdat.crystallizer import Crystallizer background_level = 0.25 CIF = 'na3sbs4_crystallized.cif' cr = Crystallizer.from_trajectory(trajectory, floating_specie='Na', resolution=0.3) result = cr.to_cif(CIF, background_level=background_level, snap_to_lower=True) print(f'space group: {result.spacegroup_symbol} (#{result.spacegroup_number})') print(f'symprec: {result.symprec} Å') print(f'wrote: {CIF} ({len(result.structure)} sites)')

space group: I-43m (#217)
symprec:     0.5 Å
wrote:       na3sbs4_crystallized.cif (128 sites)

Retrieving the immobile sites¶

To get the immobile sites, we read a CIF back, drop the sodium, and place the remaining framework in a shared lattice. The helper below does exactly that for any CIF — the crystallized output and the references alike. The reference cells are conventional unit cells, so they pass supercell=(2, 2, 2) to expand to the MD supercell first; the crystallized CIF is already the full supercell and needs no expansion.

We plot the sites with [plot_3d][gemdat.plots.plot_3d]: passing only a structure (and no volume) draws the sites and unit cell without any density isosurface. Sites are coloured by their CIF label, so the symmetry-distinct Wyckoff environments come out in different colours. Each label gets one entry in the interactive legend — click an entry to hide that species, or double-click to isolate it.

In [3]:

import re

from pymatgen.core import Structure

from gemdat.io import read_cif
from gemdat.plots import plot_3d

# Share one coordinate frame across all three structures.
lattice = result.structure.lattice


def immobile_sites(filename: str, *, supercell: tuple | None = None) -> Structure:
    """Read a CIF, keep only the immobile framework (drop the mobile Na), and
    put it in the shared lattice."""
    structure = read_cif(filename)
    if supercell:
        structure.make_supercell(supercell)
    structure.remove_species(['Na'])
    # `make_supercell` appends a per-copy index ('Sb' -> 'Sb_1', 'Sb_2', ...) to
    # every label, which would make each copy a distinct entry/colour in plot_3d.
    # Strip that suffix so we keep the original CIF labels (e.g. 'Sb', 'W', 'S')
    # and plot_3d colours the distinct Wyckoff environments together.
    labels = [re.sub(r'_\d+$', '', site.label) for site in structure]
    return Structure(
        lattice=lattice,
        species=[site.species.elements[0] for site in structure],
        coords=[site.frac_coords for site in structure],
        labels=labels,
    )


framework = immobile_sites(CIF)

print(f'crystallized framework sites: {len(framework)}')
plot_3d(structure=framework, title='Crystallized framework')

import re from pymatgen.core import Structure from gemdat.io import read_cif from gemdat.plots import plot_3d # Share one coordinate frame across all three structures. lattice = result.structure.lattice def immobile_sites(filename: str, *, supercell: tuple | None = None) -> Structure: """Read a CIF, keep only the immobile framework (drop the mobile Na), and put it in the shared lattice.""" structure = read_cif(filename) if supercell: structure.make_supercell(supercell) structure.remove_species(['Na']) # `make_supercell` appends a per-copy index ('Sb' -> 'Sb_1', 'Sb_2', ...) to # every label, which would make each copy a distinct entry/colour in plot_3d. # Strip that suffix so we keep the original CIF labels (e.g. 'Sb', 'W', 'S') # and plot_3d colours the distinct Wyckoff environments together. labels = [re.sub(r'_\d+$', '', site.label) for site in structure] return Structure( lattice=lattice, species=[site.species.elements[0] for site in structure], coords=[site.frac_coords for site in structure], labels=labels, ) framework = immobile_sites(CIF) print(f'crystallized framework sites: {len(framework)}') plot_3d(structure=framework, title='Crystallized framework')

crystallized framework sites: 80

Reference frameworks¶

The workshop ships two reference descriptions of Na₃(Sb/W)S₄. We feed each through the same immobile_sites helper. structure_1.cif is the conventional cubic I-43m cell, with the framework cations on the Sb/W sublattice.

In [4]:

framework_1 = immobile_sites(f'{DATA}/structure_1.cif', supercell=(2, 2, 2))

print(f'structure_1 framework sites: {len(framework_1)}')
plot_3d(structure=framework_1, title='structure_1 framework (I-43m)')

framework_1 = immobile_sites(f'{DATA}/structure_1.cif', supercell=(2, 2, 2)) print(f'structure_1 framework sites: {len(framework_1)}') plot_3d(structure=framework_1, title='structure_1 framework (I-43m)')

structure_1 framework sites: 80

structure_2.cif is the lower-symmetry P-42₁c description. Its framework cations sit on the Sb sublattice (the W is folded into the disordered Sb site), so only Sb and S remain after removing the sodium.

In [5]:

framework_2 = immobile_sites(f'{DATA}/structure_2.cif', supercell=(2, 2, 2))

print(f'structure_2 framework sites: {len(framework_2)}')
plot_3d(structure=framework_2, title='structure_2 framework (P-42_1c)')

framework_2 = immobile_sites(f'{DATA}/structure_2.cif', supercell=(2, 2, 2)) print(f'structure_2 framework sites: {len(framework_2)}') plot_3d(structure=framework_2, title='structure_2 framework (P-42_1c)')

structure_2 framework sites: 80