crystod.md

crystod.md#

Public MD-trajectory API of CrystOD (the crystod-md domain).

This module mirrors crystod-md --adp: the building blocks that turn a molecular-dynamics XDATCAR trajectory into a time-averaged structure with symmetry-constrained anisotropic displacement parameters (ADPs, the U_ij tensors of the _atom_site_aniso_U_* loop of a CIF file). The command itself (crystod.xdatcar_adp.main) chains these functions with the folding of the MD supercell onto one unit cell, the spglib symmetry search of the averaged structure and the CIF writer; the API exposes the reusable pieces, so the same site-symmetry constraints can be applied to displacement statistics obtained elsewhere.

Trajectory input:

  • read_xdatcar – read an XDATCAR (fixed-cell, or NpT with repeated headers) into the element list, the per-frame lattices and the fractional coordinates of every frame.

Site-symmetry constraints on U_ij (crystod-md --adp):

  • get_site_symmetry_operations – the space-group rotations that leave a site fixed;

  • build_symmetry_projector – the 6x6 projector onto the U tensors invariant under those rotations;

  • apply_symmetry_constraints – symmetrize a 3x3 U tensor with such a projector;

  • get_constraint_description – spell a projector out as U11=U22, U12=0, ... for reports.

The functions chain in that order:

import spglib
from crystod import md

symbols, lattices, frames = md.read_xdatcar("XDATCAR")
symmetry = spglib.get_symmetry(cell)        # the time-averaged unit cell
site_ops = md.get_site_symmetry_operations(
    site, symmetry["rotations"], symmetry["translations"])
projector = md.build_symmetry_projector(site_ops)
u_ij = md.apply_symmetry_constraints(u_raw, projector)

Attributes resolve lazily (PEP 562): importing this module is instant, and the implementation module crystod.xdatcar_adp (numpy and spglib) is loaded on first use. Functions reached through this namespace report bad input as ValueError instead of the SystemExit the implementation raises for the command line.

crystod.md.apply_symmetry_constraints(u_cryst, projector)[source]#

Symmetrize a U tensor with a site-symmetry projector.

This is the step of crystod-md --adp that turns the raw displacement covariance of a Wyckoff position into the U_ij written to the CIF file: the six independent components of u_cryst (upper triangle, order [U11, U22, U33, U12, U13, U23]) are multiplied by projector and reassembled into a symmetric matrix.

Parameters:

  • u_cryst (ndarray[tuple[Any, ...], dtype[float64]]) – Symmetric (3, 3) tensor on the fractional axes, in the basis the projector’s rotations act on. Only the upper triangle is read.

  • projector (ndarray[tuple[Any, ...], dtype[float64]]) – The (6, 6) matrix from build_symmetry_projector.

Returns:

The constrained (3, 3) tensor, symmetric by construction.

Return type:

ndarray[tuple[Any, …], dtype[float64]]

Example

>>> import numpy as np
>>> from crystod import md
>>> projector = md.build_symmetry_projector(          # two-fold along c
...     [np.eye(3, dtype=int), np.diag([-1, -1, 1])])
>>> u = np.array([[0.010, 0.001, 0.002],
...               [0.001, 0.012, 0.003],
...               [0.002, 0.003, 0.015]])
>>> md.apply_symmetry_constraints(u, projector)
array([[0.01 , 0.001, 0.   ],
       [0.001, 0.012, 0.   ],
       [0.   , 0.   , 0.015]])
crystod.md.build_symmetry_projector(rotations)[source]#

Build the 6x6 projector onto U tensors with a given site symmetry.

crystod-md --adp calls this once per Wyckoff position with the rotations found by get_site_symmetry_operations. A symmetric tensor is handled as the component vector [U11, U22, U33, U12, U13, U23] on the fractional axes; every rotation R maps it through R @ U @ R.T, and the projector is the average of those 6x6 maps over the group. Multiplying a component vector by it leaves the closest tensor that has the full site symmetry, which is what apply_symmetry_constraints does.

Parameters:

rotations (Sequence[numpy.ndarray]) – The rotation matrices of the site-symmetry group, each of shape (3, 3) and acting on fractional coordinates (the spglib convention). The sequence must be a complete group; [identity] gives the identity projector.

Returns:

The (6, 6) projector P (idempotent, P @ P == P).

Return type:

ndarray[tuple[Any, …], dtype[float64]]

Example

A two-fold axis along c forbids U13 and U23:

>>> import numpy as np
>>> from crystod import md
>>> projector = md.build_symmetry_projector(
...     [np.eye(3, dtype=int), np.diag([-1, -1, 1])])
>>> md.get_constraint_description(projector)
'U13=0, U23=0'
crystod.md.get_constraint_description(projector, tol=1e-06)[source]#

Spell out the relations a site-symmetry projector imposes on U.

crystod-md --adp prints the result next to every Wyckoff position and in its Site / Ueq / Constraint table, for example U11=U33, U12=0, U13=0, U23=0 for the F site of cubic ScF3. Each unit component is sent through the projector: a component that projects to zero is reported as Uij=0, and one that projects onto a later component with the same (opposite) coefficient as Uij=Ukl (Uij=-Ukl). Other linear relations, such as U12 = U11/2 on a three-fold axis of a hexagonal cell, are enforced by the projector but not spelled out.

Parameters:

  • projector (ndarray[tuple[Any, ...], dtype[float64]]) – The (6, 6) matrix from build_symmetry_projector.

  • tol (float) – Absolute tolerance below which a coefficient counts as zero and within which two coefficients count as equal.

Returns:

The relations joined by ", ", or "no constraint" for a site of symmetry 1.

Return type:

str

Example

>>> import numpy as np
>>> from crystod import md
>>> three_fold = np.array([[0, -1, 0], [1, -1, 0], [0, 0, 1]])
>>> group = [np.eye(3, dtype=int), three_fold, three_fold @ three_fold]
>>> md.get_constraint_description(md.build_symmetry_projector(group))
'U11=U22, U13=0, U23=0'
crystod.md.get_site_symmetry_operations(coords, rotations, translations, symprec=0.1)[source]#

Select the space-group operations that leave a site fixed.

This is where the ADP constraints of crystod-md --adp start: the site-symmetry group of a Wyckoff position is the subset of the space-group operations (R, t) for which R @ x + t equals x up to a lattice translation. Only the rotation parts are returned, which is all build_symmetry_projector needs.

Parameters:

  • coords (array-like) – Fractional coordinates (x, y, z) of the site.

  • rotations (array-like) – Rotation parts of the space-group operations, shape (n_ops, 3, 3), as in spglib.get_symmetry(cell).

  • translations (array-like) – The matching translation parts, shape (n_ops, 3).

  • symprec (float) – Tolerance on the Euclidean norm of the wrapped fractional difference R @ x + t - x.

Returns:

The rotation matrices of the site-symmetry group, a list of (3, 3) arrays in the order of rotations.

Example

The F site of cubic ScF3 (Pm-3m) has site symmetry 4/mmm, order 16:

>>> import numpy as np, spglib
>>> from crystod import md
>>> cell = (4.0 * np.eye(3),
...         [[0, 0, 0], [0.5, 0, 0], [0, 0.5, 0], [0, 0, 0.5]],
...         [21, 9, 9, 9])
>>> symmetry = spglib.get_symmetry(cell)
>>> site_ops = md.get_site_symmetry_operations(
...     [0.5, 0, 0], symmetry["rotations"], symmetry["translations"])
>>> len(site_ops)
16
>>> md.get_constraint_description(md.build_symmetry_projector(site_ops))
'U22=U33, U12=0, U13=0, U23=0'
crystod.md.read_xdatcar(path)[source]#

Read a VASP XDATCAR trajectory into arrays.

This is the input step of crystod-md --adp. Both layouts written by VASP are handled: the fixed-cell one, a single header followed by Direct configuration= blocks, and the variable-cell (NpT) one, in which every configuration repeats the header. The scaling factor of the header is folded into the lattice. A truncated trailing frame (a run that is still writing) is dropped silently.

Parameters:

path (str) – Path of the XDATCAR file.

Returns:

  • symbols: the element symbol of every atom, in file order (one entry per atom, so len(symbols) == n_atoms);

  • lattices: one (3, 3) lattice matrix per frame, rows being the lattice vectors in Angstrom; all identical for a fixed-cell run;

  • frames: the fractional coordinates as written by VASP (wrapped into the cell, not unwrapped), shape (n_frames, n_atoms, 3).

Return type:

A tuple (symbols, lattices, frames)

Raises:

  • FileNotFoundError – path does not exist.

  • ValueError – The file holds no complete configuration.

Example

>>> from crystod import md
>>> symbols, lattices, frames = md.read_xdatcar("XDATCAR")
>>> frames.shape          # example/30_xdatcar2adp: ScF3, 4x4x4 cell, NpT
(4001, 256, 3)
>>> symbols[0], symbols[-1], lattices[0].shape
('Sc', 'F', (3, 3))