From 3e290fa912010f9d9ecaeeb097bc43e3a5e34d5e Mon Sep 17 00:00:00 2001 From: "T. Fuchs" Date: Wed, 8 Jul 2026 10:10:35 -0400 Subject: [PATCH 01/10] feat(CTR): add symmmetry and Wyckoff detection using pyxtal --- meson_options.txt | 2 +- orgui/datautils/xrayutils/CTRsymmetry.py | 1203 +++++++++++++++++ orgui/datautils/xrayutils/CTRuc.py | 167 +++ orgui/datautils/xrayutils/__init__.py | 1 + .../xrayutils/test/test_CTRsymmetry.py | 244 ++++ pyproject.toml | 5 +- 6 files changed, 1619 insertions(+), 3 deletions(-) create mode 100644 orgui/datautils/xrayutils/CTRsymmetry.py create mode 100644 orgui/datautils/xrayutils/test/test_CTRsymmetry.py diff --git a/meson_options.txt b/meson_options.txt index 7538b5f..b59468b 100644 --- a/meson_options.txt +++ b/meson_options.txt @@ -1,6 +1,6 @@ option( 'native_optimization', type: 'boolean', - value: false, + value: true, description: 'Build C++ extensions with native CPU flags for local benchmarking' ) diff --git a/orgui/datautils/xrayutils/CTRsymmetry.py b/orgui/datautils/xrayutils/CTRsymmetry.py new file mode 100644 index 0000000..9da4e6c --- /dev/null +++ b/orgui/datautils/xrayutils/CTRsymmetry.py @@ -0,0 +1,1203 @@ +# /*########################################################################## +# +# Copyright (c) 2020-2026 Timo Fuchs +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +# THE SOFTWARE. +# +# ###########################################################################*/ +"""Symmetry metadata helpers for CTR surface unit cells. + +The construction helpers in this module can use PyXtal to assign parent-cell +atoms to Wyckoff positions. The generated surface atoms are ordinary +``CTRuc.UnitCell`` atoms; symmetry information is kept as metadata so the +structure-factor calculation remains unchanged. +""" + +from __future__ import annotations + +from dataclasses import dataclass, field +import itertools +import math + +import numpy as np + +from .HKLVlieg import Lattice + + +COORDINATE_NAMES = ("x", "y", "z") +DEFAULT_VARIABLE_NAMES = ("u", "v", "w") +SYMMETRY_SECTION_HEADERS = { + "surface_transform:", + "wyckoff_sites:", + "wyckoff_atoms:", + "wyckoff_couplings:", +} + + +@dataclass(frozen=True) +class AffineExpression: + """Represent one fractional coordinate as an affine expression. + + :param float constant: + Constant coordinate offset in fractional coordinates. + :param dict coefficients: + Mapping from variable name to fractional-coordinate coefficient. + """ + + constant: float + coefficients: dict[str, float] = field(default_factory=dict) + + def evaluate(self, variables): + """Evaluate the expression for a variable dictionary. + + :param dict variables: + Mapping from variable name to current fractional-coordinate value. + :returns: + Coordinate value in fractional units. + :rtype: + float + """ + value = self.constant + for variable, factor in self.coefficients.items(): + value += factor * variables[variable] + return value + + def shifted(self, shift): + """Return a copy with an added constant offset. + + :param float shift: + Offset in fractional coordinates. + :returns: + Shifted expression. + :rtype: + AffineExpression + """ + return AffineExpression(self.constant + shift, dict(self.coefficients)) + + +@dataclass(frozen=True) +class SurfaceCellSpec: + """Describe an oriented surface cell derived from a parent cell. + + :param parent_a: + Parent conventional-cell lengths in Angstrom. + :param parent_alpha: + Parent conventional-cell angles in degrees. + :param transform: + 3-by-3 matrix whose columns are the surface-cell vectors expressed in + parent conventional fractional coordinates. + :param origin: + Surface-cell origin in parent conventional fractional coordinates. + :param z_bounds: + Inclusive lower and exclusive upper bounds in surface fractional z. + :param layer_origins: + Optional fractional z positions used as layer starts. If provided, + atoms are assigned to intervals between successive origins. + :param translation_range: + Parent-cell integer translations tested in each direction before the + surface cell is wrapped and deduplicated. + """ + + parent_a: tuple[float, float, float] + parent_alpha: tuple[float, float, float] + transform: np.ndarray + origin: tuple[float, float, float] = (0.0, 0.0, 0.0) + z_bounds: tuple[float, float] = (0.0, 1.0) + layer_origins: tuple[float, ...] | None = None + translation_range: int | tuple[int, int] = 1 + + def __post_init__(self): + transform = np.asarray(self.transform, dtype=np.float64) + if transform.shape != (3, 3): + raise ValueError("Surface-cell transform must be a 3-by-3 matrix.") + if abs(np.linalg.det(transform)) < 1e-12: + raise ValueError("Surface-cell transform must be invertible.") + object.__setattr__(self, "transform", transform) + object.__setattr__( + self, + "parent_a", + tuple(float(v) for v in self.parent_a), + ) + object.__setattr__( + self, + "parent_alpha", + tuple(float(v) for v in self.parent_alpha), + ) + object.__setattr__(self, "origin", tuple(float(v) for v in self.origin)) + object.__setattr__(self, "z_bounds", tuple(float(v) for v in self.z_bounds)) + if self.layer_origins is not None: + object.__setattr__( + self, + "layer_origins", + tuple(float(v) % 1.0 for v in self.layer_origins), + ) + + @property + def inverse_transform(self): + """Return the parent-to-surface fractional-coordinate transform. + + :returns: + Inverse of ``transform``. + :rtype: + numpy.ndarray + """ + return np.linalg.inv(self.transform) + + def parent_translations(self): + """Generate parent conventional-cell translations to search. + + :returns: + Iterator over integer translation vectors in parent fractional + coordinates. + :rtype: + iterator + """ + if isinstance(self.translation_range, int): + lo = -self.translation_range + hi = self.translation_range + else: + lo, hi = self.translation_range + for translation in itertools.product(range(lo, hi + 1), repeat=3): + yield np.asarray(translation, dtype=np.float64) + + def surface_lattice_parameters(self): + """Return surface-cell lengths and angles. + + :returns: + Tuple ``(a, alpha)`` where ``a`` contains lengths in Angstrom and + ``alpha`` contains angles in degrees. + :rtype: + tuple + """ + parent_lattice = Lattice(self.parent_a, self.parent_alpha) + surface_matrix = parent_lattice.R_mat @ self.transform + lengths = np.linalg.norm(surface_matrix, axis=0) + alpha = _angle_between(surface_matrix[:, 1], surface_matrix[:, 2]) + beta = _angle_between(surface_matrix[:, 2], surface_matrix[:, 0]) + gamma = _angle_between(surface_matrix[:, 0], surface_matrix[:, 1]) + return lengths, np.asarray([alpha, beta, gamma], dtype=np.float64) + + def parent_to_surface(self, parent_fractional): + """Transform parent conventional fractional coordinates to surface coordinates. + + :param parent_fractional: + Fractional coordinate in the parent conventional cell. + :returns: + Fractional coordinate in the surface cell. + :rtype: + numpy.ndarray + """ + parent_fractional = np.asarray(parent_fractional, dtype=np.float64) + return self.inverse_transform @ (parent_fractional - np.asarray(self.origin)) + + +@dataclass(frozen=True) +class WyckoffSiteSpec: + """Describe one affine Wyckoff site in the parent conventional cell. + + :param str site_id: + Stable site identifier used by query and parameter APIs. + :param str element: + Element symbol for generated atoms. + :param str wyckoff_label: + Crystallographic Wyckoff label, for example ``"4f"``. + :param tuple coordinates: + Parent-cell coordinate expressions for all symmetry-generated atoms. + :param dict variables: + Current values of the independent Wyckoff variables. + :param float occ: + Site occupancy. + :param float iDW: + In-plane Debye-Waller parameter. + :param float oDW: + Out-of-plane Debye-Waller parameter. + """ + + site_id: str + element: str + wyckoff_label: str + coordinates: tuple[tuple[AffineExpression, AffineExpression, AffineExpression], ...] + variables: dict[str, float] = field(default_factory=dict) + occ: float = 1.0 + iDW: float = 0.5 + oDW: float = 0.5 + + def parent_positions(self): + """Evaluate the parent conventional-cell positions. + + :returns: + List of fractional coordinates in the parent conventional cell. + :rtype: + list + """ + positions = [] + for coordinate in self.coordinates: + positions.append( + np.asarray( + [expression.evaluate(self.variables) for expression in coordinate], + dtype=np.float64, + ) + ) + return positions + + +@dataclass(frozen=True) +class WyckoffCoupling: + """Describe one atom-coordinate dependency on a Wyckoff variable. + + :param int atom_index: + Index of the generated atom in ``UnitCell.basis``. + :param str coordinate: + Coordinate name, one of ``"x"``, ``"y"``, or ``"z"``. + :param str variable: + Wyckoff variable name. + :param float constant: + Constant part in surface fractional coordinates. + :param float factor: + Variable coefficient in surface fractional coordinates. + :param str site_id: + Site identifier owning the coupling. + """ + + atom_index: int + coordinate: str + variable: str + constant: float + factor: float + site_id: str + + def value(self, variable_value): + """Evaluate this coordinate for a variable value. + + :param float variable_value: + Variable value in fractional coordinates. + :returns: + Surface fractional coordinate. + :rtype: + float + """ + return self.constant + self.factor * variable_value + + +@dataclass(frozen=True) +class GeneratedWyckoffAtom: + """Metadata for one generated surface-cell atom. + + :param int atom_index: + Index in ``UnitCell.basis``. + :param str element: + Element symbol. + :param str site_id: + Source Wyckoff site identifier. + :param str wyckoff_label: + Source Wyckoff label. + :param numpy.ndarray parent_fractional: + Parent conventional fractional coordinate including the selected + parent-cell translation. + :param numpy.ndarray surface_fractional: + Wrapped surface-cell fractional coordinate. + :param int layer: + Assigned layer index stored in ``UnitCell.basis[:, 7]``. + :param tuple couplings: + Coordinate couplings for this atom. + """ + + atom_index: int + element: str + site_id: str + wyckoff_label: str + parent_fractional: np.ndarray + surface_fractional: np.ndarray + layer: int + couplings: tuple[WyckoffCoupling, ...] = () + + +@dataclass +class SurfaceSymmetryModel: + """Hold generated surface-cell atoms and their symmetry metadata.""" + + surface_spec: SurfaceCellSpec + sites: tuple[WyckoffSiteSpec, ...] + atoms: list[GeneratedWyckoffAtom] = field(default_factory=list) + spacegroup_number: int | None = None + spacegroup_symbol: str | None = None + + def wyckoff_sites(self, parameters=None): + """Return summary dictionaries for generated Wyckoff sites. + + :param dict parameters: + Optional ``UnitCell.parameters`` dictionary used to report whether + individual atom-coordinate parameters break a site's symmetry. + :returns: + List of site summary dictionaries. + :rtype: + list + """ + output = [] + for site in self.sites: + atoms = [ + atom.atom_index + for atom in self.atoms + if atom.site_id == site.site_id + ] + variables = dict(site.variables) + output.append( + { + "site_id": site.site_id, + "element": site.element, + "wyckoff_label": site.wyckoff_label, + "variables": variables, + "atom_indices": atoms, + "spacegroup_number": self.spacegroup_number, + "spacegroup_symbol": self.spacegroup_symbol, + "status": self.symmetry_status(site.site_id, parameters), + } + ) + return output + + def wyckoff_couplings(self, site_id=None): + """Return Wyckoff variable couplings for generated atom coordinates. + + :param str site_id: + Optional site identifier. If omitted, all couplings are returned. + :returns: + List of :class:`WyckoffCoupling` objects. + :rtype: + list + """ + couplings = [] + for atom in self.atoms: + if site_id is None or atom.site_id == site_id: + couplings.extend(atom.couplings) + return couplings + + def atom_wyckoff_metadata(self, atom_index): + """Return metadata for one generated atom. + + :param int atom_index: + Index in ``UnitCell.basis``. + :returns: + Atom metadata or ``None`` when the atom has no symmetry metadata. + :rtype: + GeneratedWyckoffAtom or None + """ + for atom in self.atoms: + if atom.atom_index == atom_index: + return atom + return None + + def symmetry_status(self, site_id, parameters=None): + """Return whether a site is still fully symmetry-preserving. + + :param str site_id: + Site identifier. + :param dict parameters: + Optional ``UnitCell.parameters`` dictionary. + :returns: + ``"metadata_only"``, ``"symmetry_preserving"``, or + ``"partially_overridden"``. + :rtype: + str + """ + if not parameters: + return "metadata_only" + + site_pairs = { + (atom.atom_index, parameter_index) + for atom in self.atoms + if atom.site_id == site_id + for parameter_index in (1, 2, 3) + } + if not site_pairs: + return "metadata_only" + + preserving = False + for kind in ("absolute", "relative"): + for parameter in parameters.get(kind, []): + pairs = set(zip(parameter.indices[0], parameter.indices[1])) + wyckoff_settings = parameter.settings.get("wyckoff", {}) + if ( + kind == "relative" + and wyckoff_settings.get("site_id") == site_id + and pairs <= site_pairs + ): + preserving = True + elif pairs & site_pairs: + return "partially_overridden" + if preserving: + return "symmetry_preserving" + return "metadata_only" + + def build_unitcell(self, name, layer_behavior="ignore"): + """Create a ``CTRuc.UnitCell`` and attach this model as metadata. + + :param str name: + Unit-cell name. + :param str layer_behavior: + Layer behavior passed to ``CTRuc.UnitCell``. + :returns: + Generated unit cell containing ordinary CTR atoms. + :rtype: + CTRuc.UnitCell + """ + from .CTRuc import UnitCell + + a, alpha = self.surface_spec.surface_lattice_parameters() + unitcell = UnitCell(a, alpha, name=name, layer_behavior=layer_behavior) + generated = generate_surface_atoms(self.surface_spec, self.sites) + for atom in generated: + site = self._site(atom.site_id) + unitcell.addAtom( + atom.element, + atom.surface_fractional, + site.iDW, + site.oDW, + site.occ, + layer=atom.layer, + ) + if self.surface_spec.layer_origins is not None: + unitcell.layerpos = { + float(index + 1): origin + for index, origin in enumerate(self.surface_spec.layer_origins) + } + self.atoms = [ + GeneratedWyckoffAtom( + atom.atom_index, + atom.element, + atom.site_id, + atom.wyckoff_label, + atom.parent_fractional, + atom.surface_fractional, + atom.layer, + atom.couplings, + ) + for atom in generated + ] + unitcell.symmetry_metadata = self + return unitcell + + def _site(self, site_id): + for site in self.sites: + if site.site_id == site_id: + return site + raise KeyError(site_id) + + +def possible_wyckoff_positions(spacegroup, style="pyxtal"): + """Return possible Wyckoff positions for a space group using PyXtal. + + PyXtal is imported lazily by this function. + + :param int spacegroup: + International space-group number. + :param str style: + PyXtal Wyckoff setting style. + :returns: + List of dictionaries with ``label``, ``multiplicity``, and ``dof``. + :rtype: + list + :raises ImportError: + If PyXtal is not installed. + """ + Group = _import_pyxtal_group() + group = Group(spacegroup, style=style) + positions = [] + for wp in group: + label = wp.get_label() + positions.append( + { + "label": label, + "multiplicity": int(label[:-1]), + "letter": label[-1], + "dof": int(wp.get_dof()), + } + ) + return positions + + +def sites_from_seed( + seed, + tol=1e-4, + a_tol=5.0, + backend="pymatgen", + style="pyxtal", + iDW=0.5, + oDW=0.5, + occ=1.0, + variable_names=DEFAULT_VARIABLE_NAMES, +): + """Assign parent atoms to Wyckoff sites using PyXtal. + + PyXtal is imported lazily by this function. ``seed`` may be any input + accepted by ``pyxtal.pyxtal().from_seed()``, including a CIF path, + pymatgen ``Structure``, or ASE ``Atoms``. + + :param seed: + Structure seed accepted by PyXtal. + :param float tol: + Coordinate tolerance passed to PyXtal. + :param float a_tol: + Angle tolerance in degrees passed to PyXtal. + :param str backend: + PyXtal seed backend. + :param str style: + PyXtal Wyckoff setting style. + :param float iDW: + In-plane Debye-Waller parameter assigned to generated sites. + :param float oDW: + Out-of-plane Debye-Waller parameter assigned to generated sites. + :param float occ: + Occupancy assigned to generated sites. + :param tuple variable_names: + Names for Wyckoff free variables in order of PyXtal free coordinates. + :returns: + Tuple ``(sites, spacegroup_number, spacegroup_symbol)``. + :rtype: + tuple + :raises ImportError: + If PyXtal is not installed. + """ + pyxtal = _import_pyxtal_class() + crystal = pyxtal() + crystal.from_seed(seed, tol=tol, a_tol=a_tol, backend=backend, style=style) + sites = tuple( + _site_from_pyxtal_site( + atom_site, + variable_names=variable_names, + iDW=iDW, + oDW=oDW, + occ=occ, + ) + for atom_site in crystal.atom_sites + ) + return sites, int(crystal.group.number), crystal.group.symbol + + +def model_from_seed( + seed, + surface_spec, + tol=1e-4, + a_tol=5.0, + backend="pymatgen", + style="pyxtal", + iDW=0.5, + oDW=0.5, + occ=1.0, + variable_names=DEFAULT_VARIABLE_NAMES, +): + """Build a surface symmetry model from a PyXtal-compatible seed. + + :param seed: + Structure seed accepted by PyXtal. + :param SurfaceCellSpec surface_spec: + Surface-cell transform and parent lattice. + :returns: + Surface symmetry model ready to build a ``UnitCell``. + :rtype: + SurfaceSymmetryModel + """ + sites, spacegroup_number, spacegroup_symbol = sites_from_seed( + seed, + tol=tol, + a_tol=a_tol, + backend=backend, + style=style, + iDW=iDW, + oDW=oDW, + occ=occ, + variable_names=variable_names, + ) + return SurfaceSymmetryModel( + surface_spec, + sites, + spacegroup_number=spacegroup_number, + spacegroup_symbol=spacegroup_symbol, + ) + + +def surface_unitcell_from_seed( + seed, + surface_spec, + name, + layer_behavior="ignore", + **kwargs, +): + """Build a generated CTR surface ``UnitCell`` from a PyXtal-compatible seed. + + :param seed: + Structure seed accepted by PyXtal. + :param SurfaceCellSpec surface_spec: + Surface-cell transform and parent lattice. + :param str name: + Unit-cell name. + :param str layer_behavior: + Layer behavior passed to ``CTRuc.UnitCell``. + :param kwargs: + Additional keyword arguments passed to :func:`model_from_seed`. + :returns: + Generated CTR unit cell with attached symmetry metadata. + :rtype: + CTRuc.UnitCell + """ + model = model_from_seed(seed, surface_spec, **kwargs) + return model.build_unitcell(name, layer_behavior=layer_behavior) + + +def rutile_110_surface_spec(parent_a, parent_alpha=(90.0, 90.0, 90.0)): + """Return the standard rutile ``(110)`` commensurate surface-cell spec. + + The surface basis columns are ``[1 -1 0]``, ``[0 0 1]``, and ``[1 1 0]`` + in parent conventional fractional coordinates. + + :param parent_a: + Parent conventional-cell lengths in Angstrom. + :param parent_alpha: + Parent conventional-cell angles in degrees. + :returns: + Surface-cell specification. + :rtype: + SurfaceCellSpec + """ + transform = np.asarray( + [ + [1.0, 0.0, 1.0], + [-1.0, 0.0, 1.0], + [0.0, 1.0, 0.0], + ], + dtype=np.float64, + ) + return SurfaceCellSpec( + tuple(parent_a), + tuple(parent_alpha), + transform, + layer_origins=(0.0, 0.5), + translation_range=1, + ) + + +def generate_surface_atoms(surface_spec, sites, tolerance=1e-8): + """Generate surface-cell atom metadata from parent Wyckoff sites. + + :param SurfaceCellSpec surface_spec: + Surface-cell transform and bounds. + :param tuple sites: + Parent-cell Wyckoff site specifications. + :param float tolerance: + Fractional-coordinate tolerance for bounds and deduplication. + :returns: + Generated atom metadata. ``atom_index`` values match the returned list + order and therefore the order used by ``SurfaceSymmetryModel``. + :rtype: + list + """ + generated = [] + seen = set() + inverse = surface_spec.inverse_transform + origin = np.asarray(surface_spec.origin, dtype=np.float64) + z_min, z_max = surface_spec.z_bounds + + for site in sites: + for coordinate in site.coordinates: + for translation in surface_spec.parent_translations(): + parent_exprs = tuple( + expression.shifted(translation[index]) + for index, expression in enumerate(coordinate) + ) + transform_exprs = tuple( + expression.shifted(translation[index] - origin[index]) + for index, expression in enumerate(coordinate) + ) + surface_exprs = _transform_expressions(inverse, transform_exprs) + surface_values = np.asarray( + [ + expression.evaluate(site.variables) + for expression in surface_exprs + ], + dtype=np.float64, + ) + if surface_values[2] < z_min - tolerance: + continue + if surface_values[2] >= z_max - tolerance: + continue + + wrapped_values, wrapped_exprs = _wrap_surface_expressions( + surface_values, + surface_exprs, + ) + key = ( + site.site_id, + tuple(np.round(wrapped_values / tolerance).astype(np.int64)), + ) + if key in seen: + continue + seen.add(key) + + layer = _assign_layer(wrapped_values[2], surface_spec.layer_origins) + atom_index = len(generated) + couplings = [] + for axis, expression in enumerate(wrapped_exprs): + for variable, factor in expression.coefficients.items(): + if not math.isclose(factor, 0.0, abs_tol=tolerance): + couplings.append( + WyckoffCoupling( + atom_index=atom_index, + coordinate=COORDINATE_NAMES[axis], + variable=variable, + constant=expression.constant, + factor=factor, + site_id=site.site_id, + ) + ) + generated.append( + GeneratedWyckoffAtom( + atom_index=atom_index, + element=site.element, + site_id=site.site_id, + wyckoff_label=site.wyckoff_label, + parent_fractional=np.asarray( + [ + expression.evaluate(site.variables) + for expression in parent_exprs + ], + dtype=np.float64, + ), + surface_fractional=wrapped_values, + layer=layer, + couplings=tuple(couplings), + ) + ) + + generated.sort( + key=lambda atom: ( + -atom.layer, + atom.site_id, + round(atom.surface_fractional[2], 10), + round(atom.surface_fractional[1], 10), + round(atom.surface_fractional[0], 10), + ) + ) + return [ + GeneratedWyckoffAtom( + atom_index=index, + element=atom.element, + site_id=atom.site_id, + wyckoff_label=atom.wyckoff_label, + parent_fractional=atom.parent_fractional, + surface_fractional=atom.surface_fractional, + layer=atom.layer, + couplings=tuple( + WyckoffCoupling( + atom_index=index, + coordinate=coupling.coordinate, + variable=coupling.variable, + constant=coupling.constant, + factor=coupling.factor, + site_id=coupling.site_id, + ) + for coupling in atom.couplings + ), + ) + for index, atom in enumerate(generated) + ] + + +def symmetry_metadata_to_lines(model): + """Serialize symmetry metadata as editable plain-text table lines. + + :param SurfaceSymmetryModel model: + Symmetry model to serialize. + :returns: + Lines without trailing newline characters. + :rtype: + list + """ + lines = [] + if model.spacegroup_number is not None: + symbol = model.spacegroup_symbol or "" + lines.append(f"spacegroup: {model.spacegroup_number} {symbol}".rstrip()) + lines.append("surface_transform:") + for row in model.surface_spec.transform: + lines.append(" " + _format_values(row)) + lines.append("surface_origin: " + _format_values(model.surface_spec.origin)) + lines.append("wyckoff_sites:") + lines.append(" site_id element wyckoff_label variables occ iDW oDW") + for site in model.sites: + variables = _format_variables(site.variables) + lines.append( + " " + + " ".join( + [ + site.site_id, + site.element, + site.wyckoff_label, + variables, + _format_float(site.occ), + _format_float(site.iDW), + _format_float(site.oDW), + ] + ) + ) + lines.append("wyckoff_atoms:") + lines.append( + " atom_index element site_id wyckoff_label " + "parent_x parent_y parent_z surface_x surface_y surface_z layer" + ) + for atom in model.atoms: + values = [ + str(atom.atom_index), + atom.element, + atom.site_id, + atom.wyckoff_label, + *_format_values(atom.parent_fractional).split(), + *_format_values(atom.surface_fractional).split(), + _format_float(atom.layer), + ] + lines.append(" " + " ".join(values)) + lines.append("wyckoff_couplings:") + lines.append(" atom_index coordinate site_id variable constant factor") + for coupling in model.wyckoff_couplings(): + lines.append( + " " + + " ".join( + [ + str(coupling.atom_index), + coupling.coordinate, + coupling.site_id, + coupling.variable, + _format_float(coupling.constant), + _format_float(coupling.factor), + ] + ) + ) + return lines + + +def symmetry_metadata_from_lines(lines, unitcell=None): + """Deserialize plain-text table lines into a symmetry metadata model. + + This parser has no PyXtal dependency and is used when loading `.xtal` and + `.xpr` files containing resolved symmetry metadata. + + :param list lines: + Symmetry metadata lines. + :param CTRuc.UnitCell unitcell: + Optional unit cell providing fallback lattice and atom values. + :returns: + Parsed symmetry model, or ``None`` when no symmetry data is present. + :rtype: + SurfaceSymmetryModel or None + """ + cleaned = [line.strip() for line in lines if line.strip()] + if not cleaned: + return None + + spacegroup_number = None + spacegroup_symbol = None + transform = np.identity(3, dtype=np.float64) + origin = (0.0, 0.0, 0.0) + site_specs = [] + atoms = [] + couplings_by_atom = {} + sections = _collect_sections(cleaned) + + if "spacegroup" in sections: + parts = sections["spacegroup"][0].split() + if parts: + spacegroup_number = int(parts[0]) + if len(parts) > 1: + spacegroup_symbol = " ".join(parts[1:]) + + if "surface_transform" in sections: + transform = np.asarray( + [ + [float(value) for value in row.split()] + for row in sections["surface_transform"] + ], + dtype=np.float64, + ) + if "surface_origin" in sections: + origin = tuple(float(value) for value in sections["surface_origin"][0].split()) + + for row in _data_rows(sections.get("wyckoff_sites", [])): + parts = row.split() + variables = _parse_variables(parts[3]) + site_specs.append( + WyckoffSiteSpec( + site_id=parts[0], + element=parts[1], + wyckoff_label=parts[2], + coordinates=(), + variables=variables, + occ=float(parts[4]), + iDW=float(parts[5]), + oDW=float(parts[6]), + ) + ) + + for row in _data_rows(sections.get("wyckoff_couplings", [])): + parts = row.split() + coupling = WyckoffCoupling( + atom_index=int(parts[0]), + coordinate=parts[1], + site_id=parts[2], + variable=parts[3], + constant=float(parts[4]), + factor=float(parts[5]), + ) + couplings_by_atom.setdefault(coupling.atom_index, []).append(coupling) + + for row in _data_rows(sections.get("wyckoff_atoms", [])): + parts = row.split() + atom_index = int(parts[0]) + atoms.append( + GeneratedWyckoffAtom( + atom_index=atom_index, + element=parts[1], + site_id=parts[2], + wyckoff_label=parts[3], + parent_fractional=np.asarray(parts[4:7], dtype=np.float64), + surface_fractional=np.asarray(parts[7:10], dtype=np.float64), + layer=int(float(parts[10])), + couplings=tuple(couplings_by_atom.get(atom_index, ())), + ) + ) + + parent_a = (1.0, 1.0, 1.0) + parent_alpha = (90.0, 90.0, 90.0) + if unitcell is not None: + parent_a = tuple(float(value) for value in unitcell.a) + parent_alpha = tuple(float(value) for value in np.rad2deg(unitcell.alpha)) + surface_spec = SurfaceCellSpec( + parent_a, + parent_alpha, + transform, + origin=origin, + layer_origins=_layer_origins_from_unitcell(unitcell), + ) + return SurfaceSymmetryModel( + surface_spec, + tuple(site_specs), + atoms=atoms, + spacegroup_number=spacegroup_number, + spacegroup_symbol=spacegroup_symbol, + ) + + +def _import_pyxtal_class(): + try: + from pyxtal import pyxtal + except ImportError as exc: + raise ImportError( + "PyXtal is required for symmetry construction. Install orGUI with " + "the optional 'symmetry' dependencies to use this feature." + ) from exc + return pyxtal + + +def _import_pyxtal_group(): + try: + from pyxtal.symmetry import Group + except ImportError as exc: + raise ImportError( + "PyXtal is required for Wyckoff lookup. Install orGUI with the " + "optional 'symmetry' dependencies to use this feature." + ) from exc + return Group + + +def _site_from_pyxtal_site(atom_site, variable_names, iDW, oDW, occ): + wp = atom_site.wp + free_values = list(wp.get_free_xyzs(atom_site.position)) + if len(free_values) > len(variable_names): + raise ValueError( + f"Not enough Wyckoff variable names for {wp.get_label()}: " + f"{len(free_values)} required." + ) + variables = { + variable_names[index]: float(value) + for index, value in enumerate(free_values) + } + base_position, coefficients = _primary_position_affine( + wp, + free_values, + tuple(variables), + ) + primary = tuple( + AffineExpression(base_position[axis], dict(coefficients[axis])) + for axis in range(3) + ) + coordinates = tuple( + _apply_operation_to_expressions(operation.affine_matrix, primary) + for operation in wp.ops + ) + element = str(atom_site.specie) + return WyckoffSiteSpec( + site_id=f"{element}_{wp.get_label()}", + element=element, + wyckoff_label=wp.get_label(), + coordinates=coordinates, + variables=variables, + occ=occ, + iDW=iDW, + oDW=oDW, + ) + + +def _primary_position_affine(wp, free_values, variable_names): + if not free_values: + return np.asarray(wp.get_position_from_free_xyzs([]), dtype=np.float64), ( + {}, + {}, + {}, + ) + zeros = np.zeros(len(free_values), dtype=np.float64) + base_position = np.asarray(wp.get_position_from_free_xyzs(zeros), dtype=np.float64) + coefficients = [] + for axis in range(3): + coefficients.append({}) + step = 1e-5 + for index, variable in enumerate(variable_names): + probe = np.zeros(len(free_values), dtype=np.float64) + probe[index] = step + position = np.asarray(wp.get_position_from_free_xyzs(probe), dtype=np.float64) + delta = (position - base_position) / step + for axis, factor in enumerate(delta): + if not math.isclose(factor, 0.0, abs_tol=1e-12): + coefficients[axis][variable] = float(factor) + return base_position, tuple(coefficients) + + +def _apply_operation_to_expressions(matrix, expressions): + matrix = np.asarray(matrix, dtype=np.float64) + transformed = [] + for row in matrix[:3]: + constant = row[3] + coefficients = {} + for factor, expression in zip(row[:3], expressions): + constant += factor * expression.constant + for variable, coefficient in expression.coefficients.items(): + coefficients[variable] = coefficients.get(variable, 0.0) + ( + factor * coefficient + ) + transformed.append(AffineExpression(float(constant), coefficients)) + return tuple(transformed) + + +def _angle_between(v1, v2): + cosine = np.dot(v1, v2) / (np.linalg.norm(v1) * np.linalg.norm(v2)) + return np.rad2deg(np.arccos(np.clip(cosine, -1.0, 1.0))) + + +def _transform_expressions(matrix, expressions): + transformed = [] + for row in matrix: + constant = 0.0 + coefficients = {} + for factor, expression in zip(row, expressions): + constant += factor * expression.constant + for variable, coefficient in expression.coefficients.items(): + coefficients[variable] = coefficients.get(variable, 0.0) + ( + factor * coefficient + ) + transformed.append(AffineExpression(constant, coefficients)) + return tuple(transformed) + + +def _wrap_surface_expressions(values, expressions): + wrapped_values = np.mod(values, 1.0) + wrapped_values[np.isclose(wrapped_values, 1.0, atol=1e-12)] = 0.0 + wrapped_expressions = [] + for value, wrapped_value, expression in zip(values, wrapped_values, expressions): + shift = wrapped_value - value + wrapped_expressions.append(expression.shifted(round(shift))) + return wrapped_values, tuple(wrapped_expressions) + + +def _assign_layer(z_value, layer_origins): + if not layer_origins: + return 0 + origins = np.asarray(sorted(layer_origins), dtype=np.float64) + z_value = z_value % 1.0 + matches = np.flatnonzero(z_value + 1e-12 >= origins) + if matches.size == 0: + return len(origins) + return int(matches[-1] + 1) + + +def _format_float(value): + return f"{float(value):.12g}" + + +def _format_values(values): + return " ".join(_format_float(value) for value in values) + + +def _format_variables(variables): + if not variables: + return "-" + return ",".join( + f"{name}={_format_float(value)}" for name, value in variables.items() + ) + + +def _parse_variables(text): + if text == "-": + return {} + variables = {} + for item in text.split(","): + name, value = item.split("=", 1) + variables[name] = float(value) + return variables + + +def _collect_sections(lines): + sections = {} + current = None + for line in lines: + if line.startswith("spacegroup:"): + sections["spacegroup"] = [line.split(":", 1)[1].strip()] + current = None + elif line.startswith("surface_origin:"): + sections["surface_origin"] = [line.split(":", 1)[1].strip()] + current = None + elif line in SYMMETRY_SECTION_HEADERS: + current = line[:-1] + sections[current] = [] + elif current is not None: + sections[current].append(line) + return sections + + +def _data_rows(rows): + for row in rows: + stripped = row.strip() + if not stripped: + continue + if stripped.split()[0] in { + "site_id", + "atom_index", + }: + continue + yield stripped + + +def _layer_origins_from_unitcell(unitcell): + if unitcell is None or not getattr(unitcell, "layerpos", None): + return None + return tuple(unitcell.layerpos[layer] for layer in sorted(unitcell.layerpos)) diff --git a/orgui/datautils/xrayutils/CTRuc.py b/orgui/datautils/xrayutils/CTRuc.py index 9d6cc64..421ef5b 100644 --- a/orgui/datautils/xrayutils/CTRuc.py +++ b/orgui/datautils/xrayutils/CTRuc.py @@ -872,6 +872,7 @@ def __init__(self, a, alpha, **keyargs): self.coherentDomainOccupancy = [1.0] self.dw_increase_constraint = np.array([], dtype=np.bool_) self._special_formfactors_present = False + self.symmetry_metadata = keyargs.get("symmetry_metadata", None) def parametersToDict(self): d = dict() @@ -1377,6 +1378,137 @@ def addRelParameter(self, indexarray, factors, limits=(-np.inf, np.inf), **keyar self.parameters["relative"].append(par) return par + def wyckoff_sites(self): + """Return Wyckoff site metadata for this unit cell. + + :returns: + List of site dictionaries. The list is empty when no symmetry + metadata is attached. + :rtype: + list + """ + if self.symmetry_metadata is None: + return [] + return self.symmetry_metadata.wyckoff_sites(self.parameters) + + def wyckoff_couplings(self, site_id=None): + """Return symmetry couplings for generated Wyckoff atom coordinates. + + :param str site_id: + Optional site identifier. If omitted, all couplings are returned. + :returns: + List of coordinate coupling metadata objects. + :rtype: + list + """ + if self.symmetry_metadata is None: + return [] + return self.symmetry_metadata.wyckoff_couplings(site_id) + + def atom_wyckoff_metadata(self, atom_index): + """Return symmetry metadata for one atom. + + :param int atom_index: + Index in ``basis``. + :returns: + Atom metadata or ``None`` when no metadata is available. + :rtype: + object or None + """ + if self.symmetry_metadata is None: + return None + return self.symmetry_metadata.atom_wyckoff_metadata(atom_index) + + def addWyckoffParameter( + self, + site_id, + variable, + limits=(-np.inf, np.inf), + absolute_limits=None, + **keyargs, + ): + """Add a symmetry-preserving relative fit parameter for a Wyckoff variable. + + The stored fit value is the change in the Wyckoff variable from the + generated coordinates. For example, fitting rutile oxygen ``u`` adds + ``factor * delta_u`` to every generated coordinate that depends on + ``u``. + + :param str site_id: + Site identifier returned by :meth:`wyckoff_sites`. + :param str variable: + Wyckoff variable name, for example ``"u"``. + :param tuple limits: + Fit limits for the variable change in fractional units. + :param tuple absolute_limits: + Optional absolute variable limits. These are converted to change + limits around the metadata variable value. + :returns: + Created relative fit parameter. + :rtype: + CTRutil.Parameter + :raises ValueError: + If no matching affine couplings exist. + """ + if self.symmetry_metadata is None: + raise ValueError(f"UnitCell {self.name} has no symmetry metadata.") + if absolute_limits is not None: + if limits != (-np.inf, np.inf): + raise ValueError("Use either limits or absolute_limits, not both.") + site = next( + ( + site + for site in self.wyckoff_sites() + if site["site_id"] == site_id + ), + None, + ) + if site is None or variable not in site["variables"]: + raise ValueError( + f"Wyckoff site {site_id} has no variable {variable}." + ) + variable_value = site["variables"][variable] + limits = ( + absolute_limits[0] - variable_value, + absolute_limits[1] - variable_value, + ) + + couplings = [ + coupling + for coupling in self.wyckoff_couplings(site_id) + if coupling.variable == variable + ] + if not couplings: + raise ValueError( + f"No affine couplings found for Wyckoff site {site_id} " + f"and variable {variable}." + ) + atoms = np.asarray( + [coupling.atom_index for coupling in couplings], + dtype=np.intp, + ) + coordinates = tuple(coupling.coordinate for coupling in couplings) + factors = np.asarray( + [coupling.factor for coupling in couplings], + dtype=np.float64, + ) + + keyargs.setdefault("name", f"{self.name} {site_id}_{variable}_wyckoff") + settings = keyargs.setdefault("wyckoff", {}) + settings.update( + { + "site_id": site_id, + "variable": variable, + "value_kind": "delta", + } + ) + return self.addRelParameter( + (atoms, coordinates), + factors, + limits=limits, + **keyargs, + ) + def showFitparameters(self): print(self.fitparameterList()) @@ -2279,6 +2411,10 @@ def parameterStr(self, showErrors=True): st += "layer_cycle: {}\n".format( ", ".join(str(layer) for layer in self._explicit_layer_cycle) ) + if self.symmetry_metadata is not None: + from .CTRsymmetry import symmetry_metadata_to_lines + + st += "\n".join(symmetry_metadata_to_lines(self.symmetry_metadata)) + "\n" return st def parameterStrRod(self): @@ -2537,8 +2673,25 @@ def fromStr(cls, string): layerpos = dict() layer_behaviour = "ignore" layer_cycle = None + symmetry_lines = [] + reading_symmetry = False for l in f: # noqa: E741 line = l.rsplit("//")[0] + stripped = line.strip() + if reading_symmetry or stripped.startswith( + ( + "spacegroup:", + "surface_transform:", + "surface_origin:", + "wyckoff_sites:", + "wyckoff_atoms:", + "wyckoff_couplings:", + ) + ): + reading_symmetry = True + if stripped: + symmetry_lines.append(stripped) + continue if line.startswith("layerpos:"): if "=" in line: try: @@ -2704,6 +2857,13 @@ def fromStr(cls, string): uc.basis_0 = np.copy(uc.basis) uc.dw_increase_constraint = np.ones(uc.basis.shape[0], dtype=np.bool_) uc._test_special_formfactors() + if symmetry_lines: + from .CTRsymmetry import symmetry_metadata_from_lines + + uc.symmetry_metadata = symmetry_metadata_from_lines( + symmetry_lines, + uc, + ) return uc elif len(basis) == 1: if ( @@ -2725,6 +2885,13 @@ def fromStr(cls, string): uc.basis_0 = np.copy(uc.basis) uc.dw_increase_constraint = np.ones(uc.basis.shape[0], dtype=np.bool_) uc._test_special_formfactors() + if symmetry_lines: + from .CTRsymmetry import symmetry_metadata_from_lines + + uc.symmetry_metadata = symmetry_metadata_from_lines( + symmetry_lines, + uc, + ) return uc else: raise ValueError( diff --git a/orgui/datautils/xrayutils/__init__.py b/orgui/datautils/xrayutils/__init__.py index 1b35c71..5a5560f 100644 --- a/orgui/datautils/xrayutils/__init__.py +++ b/orgui/datautils/xrayutils/__init__.py @@ -31,6 +31,7 @@ __all__ = [ "CTRcalc", "CTRplotutil", + "CTRsymmetry", "DetectorCalibration", "HKLVlieg", "id31_tools", diff --git a/orgui/datautils/xrayutils/test/test_CTRsymmetry.py b/orgui/datautils/xrayutils/test/test_CTRsymmetry.py new file mode 100644 index 0000000..5be561b --- /dev/null +++ b/orgui/datautils/xrayutils/test/test_CTRsymmetry.py @@ -0,0 +1,244 @@ +# /*########################################################################## +# +# Copyright (c) 2020-2026 Timo Fuchs +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +# THE SOFTWARE. +# +# ###########################################################################*/ +import importlib.util +import subprocess +import sys +import unittest +from unittest import mock + +import numpy as np + +from .. import CTRsymmetry +from ..CTRuc import UnitCell + + +HAS_PYXTAL = importlib.util.find_spec("pyxtal") is not None +HAS_PYMATGEN = importlib.util.find_spec("pymatgen") is not None + + +@unittest.skipUnless( + HAS_PYXTAL and HAS_PYMATGEN, + "PyXtal symmetry tests require PyXtal", +) +class TestPyxtalRutileSurfaceSymmetry(unittest.TestCase): + def make_rutile_seed(self, oxygen_u=0.30569): + """Return a conventional rutile parent structure.""" + from pymatgen.core import Lattice, Structure + + u = oxygen_u + return Structure( + Lattice.tetragonal(4.653255, 2.9692), + ["Ru", "Ru", "O", "O", "O", "O"], + [ + [0.0, 0.0, 0.0], + [0.5, 0.5, 0.5], + [u, u, 0.0], + [1.0 - u, 1.0 - u, 0.0], + [0.5 - u, 0.5 + u, 0.5], + [0.5 + u, 0.5 - u, 0.5], + ], + ) + + def make_rutile_110_unitcell(self): + surface_spec = CTRsymmetry.rutile_110_surface_spec( + (4.653255, 4.653255, 2.9692), + ) + model = CTRsymmetry.model_from_seed( + self.make_rutile_seed(), + surface_spec, + tol=1e-3, + iDW=0.5033, + oDW=0.5033, + ) + return model.build_unitcell("RuO2") + + def test_pyxtal_assigns_rutile_wyckoff_sites(self): + sites, number, symbol = CTRsymmetry.sites_from_seed( + self.make_rutile_seed(), + tol=1e-3, + ) + + self.assertEqual(number, 136) + self.assertEqual(symbol, "P42/mnm") + self.assertEqual([site.site_id for site in sites], ["Ru_2a", "O_4f"]) + self.assertEqual(sites[1].variables, {"u": 0.30569}) + + def test_possible_wyckoff_positions_exposes_group_table(self): + positions = CTRsymmetry.possible_wyckoff_positions(136) + + labels = [position["label"] for position in positions] + self.assertIn("2a", labels) + self.assertIn("4f", labels) + oxygen_position = next( + position for position in positions if position["label"] == "4f" + ) + self.assertEqual(oxygen_position["dof"], 1) + + def test_rutile_110_generation_reproduces_surface_basis(self): + unitcell = self.make_rutile_110_unitcell() + + self.assertEqual(unitcell.basis.shape, (12, 8)) + self.assertEqual(unitcell.parameters["absolute"], []) + self.assertEqual(unitcell.parameters["relative"], []) + np.testing.assert_allclose(unitcell.a, [6.5807, 2.9692, 6.5807], atol=5e-5) + self.assertEqual(unitcell.layerpos, {1.0: 0.0, 2.0: 0.5}) + + expected = { + ("O", 0.5, 0.0, 0.80569, 2.0), + ("O", 0.0, 0.0, 0.69431, 2.0), + ("Ru", 0.5, 0.0, 0.5, 2.0), + ("Ru", 0.0, 0.5, 0.5, 2.0), + ("O", 0.30569, 0.5, 0.5, 2.0), + ("O", 0.69431, 0.5, 0.5, 2.0), + ("O", 0.0, 0.0, 0.30569, 1.0), + ("O", 0.5, 0.0, 0.19431, 1.0), + ("Ru", 0.0, 0.0, 0.0, 1.0), + ("Ru", 0.5, 0.5, 0.0, 1.0), + ("O", 0.19431, 0.5, 0.0, 1.0), + ("O", 0.80569, 0.5, 0.0, 1.0), + } + actual = { + ( + name, + round(row[1], 5), + round(row[2], 5), + round(row[3], 5), + row[7], + ) + for name, row in zip(unitcell.names, unitcell.basis) + } + self.assertEqual(actual, expected) + + def test_wyckoff_query_exposes_oxygen_u_couplings(self): + unitcell = self.make_rutile_110_unitcell() + + sites = unitcell.wyckoff_sites() + self.assertEqual([site["site_id"] for site in sites], ["Ru_2a", "O_4f"]) + self.assertEqual(sites[0]["spacegroup_number"], 136) + self.assertEqual(sites[1]["variables"], {"u": 0.30569}) + self.assertEqual(sites[1]["status"], "metadata_only") + + couplings = unitcell.wyckoff_couplings("O_4f") + self.assertEqual(len(couplings), 8) + self.assertTrue(all(coupling.variable == "u" for coupling in couplings)) + expressions = { + (round(coupling.constant, 5), round(coupling.factor, 5)) + for coupling in couplings + } + self.assertEqual( + expressions, + {(0.0, 1.0), (1.0, -1.0), (0.5, 1.0), (0.5, -1.0)}, + ) + + def test_wyckoff_parameter_preserves_rutile_u_symmetry(self): + unitcell = self.make_rutile_110_unitcell() + couplings = unitcell.wyckoff_couplings("O_4f") + original = unitcell.basis.copy() + + parameter = unitcell.addWyckoffParameter( + "O_4f", + "u", + absolute_limits=(0.2, 0.4), + ) + + self.assertEqual(parameter.settings["wyckoff"]["value_kind"], "delta") + np.testing.assert_allclose(unitcell.getInitialParameters(), [0.0]) + np.testing.assert_allclose( + unitcell.getStartParamAndLimits()[1:], + ([-0.10569], [0.09431]), + ) + self.assertEqual(unitcell.wyckoff_sites()[1]["status"], "symmetry_preserving") + + unitcell.setFitParameters([0.01]) + + for coupling in couplings: + column = unitcell.parameterLookup[coupling.coordinate] + expected = original[coupling.atom_index, column] + 0.01 * coupling.factor + self.assertAlmostEqual( + unitcell.basis[coupling.atom_index, column], + expected, + ) + + def test_manual_atom_parameter_marks_wyckoff_site_as_partially_overridden(self): + unitcell = self.make_rutile_110_unitcell() + atom_index = unitcell.wyckoff_couplings("O_4f")[0].atom_index + + unitcell.addFitParameter( + ([atom_index], ["z"]), + limits=(-1.0, 2.0), + name="manual_oxygen_coordinate", + ) + + self.assertEqual(unitcell.wyckoff_sites()[1]["status"], "partially_overridden") + + def test_symmetry_metadata_round_trips_without_pyxtal_construction(self): + unitcell = self.make_rutile_110_unitcell() + text = unitcell.toStr() + + restored = UnitCell.fromStr(text) + + self.assertEqual(restored.wyckoff_sites()[1]["variables"], {"u": 0.30569}) + self.assertEqual(len(restored.wyckoff_couplings("O_4f")), 8) + restored.addWyckoffParameter("O_4f", "u", absolute_limits=(0.2, 0.4)) + self.assertEqual(restored.wyckoff_sites()[1]["status"], "symmetry_preserving") + + +class TestOptionalSymmetryImports(unittest.TestCase): + def test_missing_pyxtal_reports_optional_dependency(self): + real_import = __import__ + + def import_without_pyxtal(name, *args, **kwargs): + if name.startswith("pyxtal"): + raise ImportError("simulated missing pyxtal") + return real_import(name, *args, **kwargs) + + with mock.patch("builtins.__import__", side_effect=import_without_pyxtal): + with self.assertRaisesRegex(ImportError, "optional 'symmetry'"): + CTRsymmetry.possible_wyckoff_positions(136) + + def test_importing_ctr_modules_does_not_import_heavy_symmetry_dependencies(self): + command = [ + sys.executable, + "-c", + ( + "import sys; " + "import orgui.datautils.xrayutils.CTRcalc; " + "import orgui.datautils.xrayutils.CTRuc; " + "import orgui.datautils.xrayutils.CTRsymmetry; " + "print(any(name in sys.modules for name in " + "('pyxtal', 'spglib', 'pymatgen', 'ase')))" + ), + ] + completed = subprocess.run( + command, + check=True, + capture_output=True, + text=True, + ) + + self.assertEqual(completed.stdout.strip(), "False") + + +if __name__ == "__main__": + unittest.main() diff --git a/pyproject.toml b/pyproject.toml index b7f83cf..19ad146 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -51,8 +51,9 @@ Citation = "https://doi.org/10.5281/zenodo.12592485" console = ['qtconsole'] speedup = ['numba'] extendedfilesupport = ['ase', 'hdf5plugin'] -full = ['qtconsole', 'ase', 'numba' , 'PyOpenGL', 'hdf5plugin'] -full-devel = ['qtconsole', 'ase', 'numba' , 'PyOpenGL', 'hdf5plugin', 'ruff', 'pyupgrade'] +symmetry = ['pyxtal'] +full = ['qtconsole', 'ase', 'numba' , 'PyOpenGL', 'hdf5plugin', 'pyxtal'] +full-devel = ['qtconsole', 'ase', 'numba' , 'PyOpenGL', 'hdf5plugin', 'pyxtal', 'ruff', 'pyupgrade'] [build-system] requires = ["meson-python", "pybind11", "setuptools-scm>=8", "numpy"] From 7b516056bba30e8ab18fc64060655f587ac5cfa1 Mon Sep 17 00:00:00 2001 From: "T. Fuchs" Date: Wed, 8 Jul 2026 11:50:15 -0400 Subject: [PATCH 02/10] feat(CTR): add layer-aware and Wyckoff-aware unit cell translate --- orgui/datautils/xrayutils/CTRuc.py | 238 ++++++++++++++++++ .../datautils/xrayutils/test/test_CTRcalc.py | 227 ++++++++++++++++- 2 files changed, 464 insertions(+), 1 deletion(-) diff --git a/orgui/datautils/xrayutils/CTRuc.py b/orgui/datautils/xrayutils/CTRuc.py index 421ef5b..cb5bcab 100644 --- a/orgui/datautils/xrayutils/CTRuc.py +++ b/orgui/datautils/xrayutils/CTRuc.py @@ -37,6 +37,8 @@ import warnings import random import math +import copy +import dataclasses from scipy.ndimage import gaussian_filter1d import os import re @@ -1140,6 +1142,242 @@ def layer_cycle(self): def layer_cycle(self, layers): self._explicit_layer_cycle = tuple(layers) + @staticmethod + def _translation_from_affine_input(translation): + translation = np.asarray(translation, dtype=np.float64) + if translation.shape == (3,): + return translation + if translation.shape == (3, 4): + linear = translation[:, :3] + vector = translation[:, 3] + elif translation.shape == (4, 4): + if not np.allclose(translation[3], [0.0, 0.0, 0.0, 1.0]): + raise ValueError( + "Homogeneous affine transforms must end with " + "[0, 0, 0, 1]." + ) + linear = translation[:3, :3] + vector = translation[:3, 3] + else: + raise ValueError( + "translation must be a vector with shape (3,), " + "a 3-by-4 matrix, or a homogeneous 4-by-4 matrix." + ) + if not np.allclose(linear, np.identity(3)): + raise ValueError( + "UnitCell.translate_layered currently supports " + "translations only; the affine linear part must be identity." + ) + return vector + + @staticmethod + def _remap_parameter_atom_indices(parameter, old_to_new): + if not isinstance(parameter.indices, tuple): + return + atoms, parindexes = parameter.indices + atoms = np.asarray(atoms, dtype=np.intp) + parameter.indices = (old_to_new[atoms], parindexes) + + @staticmethod + def _translate_symmetry_metadata(metadata, old_to_new, layer_map, coordinate_shift): + if metadata is None: + return None + metadata = copy.deepcopy(metadata) + coordinate_index = {"x": 0, "y": 1, "z": 2} + atoms = [] + for atom in metadata.atoms: + couplings = [] + for coupling in atom.couplings: + couplings.append( + dataclasses.replace( + coupling, + atom_index=int(old_to_new[coupling.atom_index]), + constant=( + coupling.constant + + coordinate_shift[coordinate_index[coupling.coordinate]] + ), + ) + ) + atoms.append( + dataclasses.replace( + atom, + atom_index=int(old_to_new[atom.atom_index]), + surface_fractional=( + np.asarray(atom.surface_fractional, dtype=np.float64) + + coordinate_shift + ), + layer=layer_map[atom.layer], + couplings=tuple(couplings), + ) + ) + metadata.atoms = sorted(atoms, key=lambda atom: atom.atom_index) + return metadata + + def translate_layered(self, translation, name=None): + """Return a translated copy with cyclic z-layer wrapping. + + The translation is expressed in unit-cell fractional coordinates. The + ``x`` and ``y`` translations must be integer unit-cell offsets, and + ``z`` is an integer number of structural-layer steps. Atom fractional + z coordinates and layer origins are shifted by ``z / len(layer_cycle)`` + while layer identifiers are reassigned through the ordered + :class:`CTRstacking.LayerCycle`. A positive z step moves lower layers + upward and wraps the former top layer to the bottom. + + Symmetry metadata is copied only when atom ordering is unchanged. When + layer cycling regroups atom rows, generated Wyckoff metadata is dropped + because it describes the pre-transform atom indices. + + :param translation: + Translation vector with shape ``(3,)``, affine matrix with shape + ``(3, 4)``, or homogeneous affine matrix with shape ``(4, 4)``. + :param str name: + Optional name for the returned unit cell. Defaults to this unit + cell's name. + :returns: + Transformed unit-cell copy. + :rtype: + UnitCell + :raises ValueError: + If the affine linear part is not identity or any requested + translation is not an integer. + """ + translation = self._translation_from_affine_input(translation) + xy_translation = translation[:2] + if not np.allclose(xy_translation, np.rint(xy_translation)): + raise ValueError("x and y affine translations must be integers.") + if not np.isclose(translation[2], np.rint(translation[2])): + raise ValueError("z affine translation must be an integer layer step.") + + xy_translation = np.rint(xy_translation).astype(np.float64) + z_steps = int(np.rint(translation[2])) + + transformed = copy.deepcopy(self) + if name is not None: + transformed.name = name + + if self.basis.size == 0: + if z_steps: + raise ValueError("Cannot cyclically shift layers in an empty UnitCell.") + return transformed + + cycle = self.layer_cycle.layers + layer_count = len(cycle) + z_steps_effective = z_steps % layer_count + z_translation = z_steps / layer_count + + if z_steps_effective == 0 and np.all(xy_translation == 0): + if z_steps == 0: + return transformed + if z_steps_effective == 0: + transformed.basis[:, 1:3] += xy_translation + transformed.basis[:, 3] += z_translation + transformed.basis_0[:, 1:3] += xy_translation + transformed.basis_0[:, 3] += z_translation + if transformed._basis_parvalues is not None: + transformed._basis_parvalues[:, 1:3] += xy_translation + transformed._basis_parvalues[:, 3] += z_translation + for layer in cycle: + layer_key = float(layer) + if layer_key not in transformed.layerpos: + where = self.basis[:, 7] == layer + transformed.layerpos[layer_key] = float( + np.mean(self.basis[where, 3]) + ) + transformed.layerpos[layer_key] += z_translation + coordinate_shift = np.array( + [xy_translation[0], xy_translation[1], z_translation], + dtype=np.float64, + ) + old_to_new = np.arange(self.basis.shape[0], dtype=np.intp) + layer_map = {layer: layer for layer in cycle} + transformed.symmetry_metadata = self._translate_symmetry_metadata( + self.symmetry_metadata, + old_to_new, + layer_map, + coordinate_shift, + ) + return transformed + + layer_positions = { + layer: float(self.layerpos.get(float(layer), np.nan)) for layer in cycle + } + for layer, position in layer_positions.items(): + if np.isnan(position): + where = self.basis[:, 7] == layer + layer_positions[layer] = float(np.mean(self.basis[where, 3])) + + layer_map = { + layer: cycle[(index + z_steps_effective) % layer_count] + for index, layer in enumerate(cycle) + } + + def transform_basis_values(values): + if values is None or values.size == 0: + return values + values_new = np.array(values, copy=True) + values_new[:, 1:3] += xy_translation + values_new[:, 3] += z_translation + for layer in cycle: + where = values[:, 7] == layer + if not np.any(where): + continue + new_layer = layer_map[layer] + values_new[where, 7] = new_layer + return values_new + + basis = transform_basis_values(self.basis) + basis_0 = transform_basis_values(self.basis_0) + basis_parvalues = transform_basis_values(self._basis_parvalues) + + order_lookup = {layer: index for index, layer in enumerate(cycle)} + row_order = np.array( + sorted(range(basis.shape[0]), key=lambda i: (order_lookup[basis[i, 7]], i)), + dtype=np.intp, + ) + old_to_new = np.empty_like(row_order) + old_to_new[row_order] = np.arange(row_order.size) + + transformed.basis = basis[row_order] + transformed.basis_0 = basis_0[row_order] + if basis_parvalues is not None: + transformed._basis_parvalues = basis_parvalues[row_order] + if self.errors is not None: + transformed.errors = np.array(self.errors, copy=True)[row_order] + if self._errors_parvalues is not None: + transformed._errors_parvalues = np.array( + self._errors_parvalues, copy=True + )[row_order] + transformed.names = [self.names[i] for i in row_order] + transformed.dw_increase_constraint = self.dw_increase_constraint[row_order] + if hasattr(self, "f"): + transformed.f = self.f[row_order] + + for parameters in transformed.parameters.values(): + for parameter in parameters: + self._remap_parameter_atom_indices(parameter, old_to_new) + + transformed.layerpos = copy.deepcopy(self.layerpos) + for layer in cycle: + new_layer = layer_map[layer] + transformed.layerpos[float(new_layer)] = ( + layer_positions[layer] + z_translation + ) + if transformed._explicit_layer_cycle is not None: + transformed._explicit_layer_cycle = tuple(cycle) + coordinate_shift = np.array( + [xy_translation[0], xy_translation[1], z_translation], + dtype=np.float64, + ) + transformed.symmetry_metadata = self._translate_symmetry_metadata( + self.symmetry_metadata, + old_to_new, + layer_map, + coordinate_shift, + ) + transformed._test_special_formfactors() + return transformed + @property def layer_behavior(self): """Return how layer selection is handled during crystal stacking.""" diff --git a/orgui/datautils/xrayutils/test/test_CTRcalc.py b/orgui/datautils/xrayutils/test/test_CTRcalc.py index a3e4154..2d98a60 100644 --- a/orgui/datautils/xrayutils/test/test_CTRcalc.py +++ b/orgui/datautils/xrayutils/test/test_CTRcalc.py @@ -38,7 +38,7 @@ import numpy as np from ... import util -from .. import CTRcalc, CTRfilm, CTRplotutil, CTRuc +from .. import CTRcalc, CTRfilm, CTRplotutil, CTRsymmetry, CTRuc from ..CTRdistributions import PoissonProfile, SkellamProfile @@ -344,6 +344,231 @@ def test_unitcell_selects_layer_after_below_layer(self): self.assertEqual(unitcell.end_layer_number, 1) self.assertEqual(unitcell.layer_behaviour, "select") + def test_unitcell_translate_layered_cycles_top_to_bottom(self): + unitcell = self.make_layered_unitcell() + unitcell.basis[:, 1] = [0.1, 0.2, 0.3] + unitcell.basis[:, 3] += 0.05 + unitcell.basis_0[:] = unitcell.basis + + transformed = unitcell.translate_layered([0, 0, 1]) + + np.testing.assert_allclose(transformed.basis[:, 1], [0.3, 0.1, 0.2]) + np.testing.assert_allclose( + transformed.basis[:, 3], + [1.0 + 0.05, 1 / 3 + 0.05, 2 / 3 + 0.05], + ) + np.testing.assert_array_equal(transformed.basis[:, 7], [1, 2, 3]) + self.assertEqual( + transformed.layerpos, + {0.0: 0.0, 1.0: 1.0, 2.0: 1 / 3, 3.0: 2 / 3}, + ) + self.assertTrue(transformed.layer_cycle.is_rotation_of(unitcell.layer_cycle)) + + def test_unitcell_translate_layered_cycles_bottom_to_top(self): + unitcell = self.make_layered_unitcell() + unitcell.basis[:, 1] = [0.1, 0.2, 0.3] + unitcell.basis[:, 3] += 0.05 + unitcell.basis_0[:] = unitcell.basis + + transformed = unitcell.translate_layered([0, 0, -1]) + + np.testing.assert_allclose(transformed.basis[:, 1], [0.2, 0.3, 0.1]) + np.testing.assert_allclose( + transformed.basis[:, 3], + [0.05, 1 / 3 + 0.05, -1 / 3 + 0.05], + ) + np.testing.assert_array_equal(transformed.basis[:, 7], [1, 2, 3]) + + def test_unitcell_translate_layered_preserves_noop_metadata(self): + unitcell = self.make_layered_unitcell() + + transformed = unitcell.translate_layered([0, 0, 0]) + + self.assertIsNot(transformed, unitcell) + np.testing.assert_allclose(transformed.basis, unitcell.basis) + np.testing.assert_allclose(transformed.basis_0, unitcell.basis_0) + self.assertEqual(transformed.layerpos, unitcell.layerpos) + self.assertEqual(transformed.toStr(), unitcell.toStr()) + + def test_unitcell_translate_layered_translates_full_two_layer_cycle(self): + unitcell = CTRcalc.UnitCell( + [3.0, 3.0, 6.0], + [90.0, 90.0, 90.0], + name="two_layer", + ) + for layer, z in zip((1, 2), (0.0, 0.5)): + unitcell.addAtom("C", [0.0, 0.0, z], 0.1, 0.1, 1.0, layer=layer) + unitcell.layerpos[float(layer)] = z + + transformed = unitcell.translate_layered([0, 0, -2]) + + np.testing.assert_allclose(transformed.basis[:, 3], [-1.0, -0.5]) + np.testing.assert_allclose(transformed.basis_0[:, 3], [-1.0, -0.5]) + np.testing.assert_array_equal(transformed.basis[:, 7], unitcell.basis[:, 7]) + self.assertEqual( + transformed.layerpos, + {0.0: 0.0, 1.0: -1.0, 2.0: -0.5}, + ) + + def test_unitcell_translate_layered_preserves_full_cycle_symmetry_metadata(self): + unitcell = CTRcalc.UnitCell( + [3.0, 3.0, 6.0], + [90.0, 90.0, 90.0], + name="two_layer_reversed", + ) + unitcell.addAtom("C", [0.0, 0.0, 0.5], 0.1, 0.1, 1.0, layer=2) + unitcell.layerpos[2.0] = 0.5 + unitcell.addAtom("C", [0.0, 0.0, 0.0], 0.1, 0.1, 1.0, layer=1) + unitcell.layerpos[1.0] = 0.0 + unitcell.symmetry_metadata = CTRsymmetry.SurfaceSymmetryModel( + CTRsymmetry.SurfaceCellSpec( + (3.0, 3.0, 6.0), + (90.0, 90.0, 90.0), + np.identity(3), + ), + (), + atoms=[ + CTRsymmetry.GeneratedWyckoffAtom( + atom_index=0, + element="C", + site_id="C_0", + wyckoff_label="1a", + parent_fractional=np.array([0.0, 0.0, 0.5]), + surface_fractional=np.array([0.0, 0.0, 0.5]), + layer=2, + couplings=( + CTRsymmetry.WyckoffCoupling( + atom_index=0, + coordinate="z", + variable="u", + constant=0.5, + factor=1.0, + site_id="C_0", + ), + ), + ), + CTRsymmetry.GeneratedWyckoffAtom( + atom_index=1, + element="C", + site_id="C_1", + wyckoff_label="1a", + parent_fractional=np.array([0.0, 0.0, 0.0]), + surface_fractional=np.array([0.0, 0.0, 0.0]), + layer=1, + ), + ], + ) + + transformed = unitcell.translate_layered([0, 0, -2]) + + np.testing.assert_array_equal(transformed.basis[:, 7], [2, 1]) + np.testing.assert_allclose(transformed.basis[:, 3], [-0.5, -1.0]) + self.assertEqual(transformed.layerpos[2.0], -0.5) + self.assertEqual(transformed.layerpos[1.0], -1.0) + self.assertIsNotNone(transformed.symmetry_metadata) + atom = transformed.atom_wyckoff_metadata(0) + self.assertEqual(atom.layer, 2) + np.testing.assert_allclose(atom.surface_fractional, [0.0, 0.0, -0.5]) + self.assertEqual(atom.couplings[0].atom_index, 0) + self.assertAlmostEqual(atom.couplings[0].constant, -0.5) + + def test_unitcell_translate_layered_translates_in_plane(self): + unitcell = self.make_layered_unitcell() + original_basis = np.copy(unitcell.basis) + affine = np.array( + [ + [1.0, 0.0, 0.0, 2.0], + [0.0, 1.0, 0.0, -1.0], + [0.0, 0.0, 1.0, 0.0], + [0.0, 0.0, 0.0, 1.0], + ] + ) + transformed = unitcell.translate_layered(affine, name="shifted") + + np.testing.assert_allclose(transformed.basis[:, 1], original_basis[:, 1] + 2) + np.testing.assert_allclose(transformed.basis[:, 2], original_basis[:, 2] - 1) + np.testing.assert_array_equal(transformed.basis[:, 7], original_basis[:, 7]) + np.testing.assert_allclose(unitcell.basis, original_basis) + self.assertEqual(transformed.name, "shifted") + + def test_unitcell_translate_layered_remaps_fit_parameter_atoms(self): + unitcell = self.make_layered_unitcell() + unitcell.addFitParameter(([2], "x"), name="top_x") + + transformed = unitcell.translate_layered([0, 0, 1]) + + atoms, parameters = transformed.parameters["absolute"][0].indices + np.testing.assert_array_equal(atoms, [0]) + np.testing.assert_array_equal(parameters, [1]) + + def test_unitcell_translate_layered_remaps_symmetry_metadata_atoms(self): + unitcell = self.make_layered_unitcell() + unitcell.symmetry_metadata = CTRsymmetry.SurfaceSymmetryModel( + CTRsymmetry.SurfaceCellSpec( + (3.0, 3.0, 6.0), + (90.0, 90.0, 90.0), + np.identity(3), + ), + (), + atoms=[ + CTRsymmetry.GeneratedWyckoffAtom( + atom_index=2, + element="C", + site_id="C_2", + wyckoff_label="1a", + parent_fractional=np.array([0.0, 0.0, 2 / 3]), + surface_fractional=np.array([0.0, 0.0, 2 / 3]), + layer=3, + couplings=( + CTRsymmetry.WyckoffCoupling( + atom_index=2, + coordinate="z", + variable="u", + constant=2 / 3, + factor=1.0, + site_id="C_2", + ), + ), + ), + ], + ) + + transformed = unitcell.translate_layered([0, 0, 1]) + + atom = transformed.atom_wyckoff_metadata(0) + self.assertIsNotNone(atom) + self.assertEqual(atom.layer, 1) + np.testing.assert_allclose(atom.surface_fractional, [0.0, 0.0, 1.0]) + self.assertEqual(atom.couplings[0].atom_index, 0) + self.assertAlmostEqual(atom.couplings[0].constant, 1.0) + + def test_unitcell_translate_layered_rejects_unsupported_affines(self): + unitcell = self.make_layered_unitcell() + + with self.assertRaisesRegex(ValueError, "x and y"): + unitcell.translate_layered([0.5, 0, 0]) + with self.assertRaisesRegex(ValueError, "linear part"): + unitcell.translate_layered( + np.array( + [ + [0.0, -1.0, 0.0, 0.0], + [1.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 1.0, 0.0], + ] + ) + ) + + def test_unitcell_translate_layered_stacks_with_existing_layers(self): + transformed = self.make_layered_unitcell().translate_layered([0, 0, 1]) + film = CTRfilm.Film(transformed) + film.basis[0] = 2 + film.start_layer_number = 1 + film.createLayers() + + self.assertEqual(film.start_layer_number, 2) + self.assertEqual(film.end_layer_number, 3) + np.testing.assert_array_equal(film.layer_order, [2, 3, 1]) + def test_film_and_poisson_rotate_cyclic_layer_order(self): film = CTRfilm.Film(self.make_layered_unitcell("film")) film.basis[0] = 2 From 880a7b75000730388931bae7b62fe2a128607f32 Mon Sep 17 00:00:00 2001 From: "T. Fuchs" Date: Wed, 8 Jul 2026 14:12:11 -0400 Subject: [PATCH 03/10] feat(CTR): symmetry affine transforms and supercell creation --- orgui/datautils/xrayutils/CTRuc.py | 453 +++++++++++++++++- .../datautils/xrayutils/test/test_CTRcalc.py | 195 ++++++++ 2 files changed, 626 insertions(+), 22 deletions(-) diff --git a/orgui/datautils/xrayutils/CTRuc.py b/orgui/datautils/xrayutils/CTRuc.py index cb5bcab..d609a1e 100644 --- a/orgui/datautils/xrayutils/CTRuc.py +++ b/orgui/datautils/xrayutils/CTRuc.py @@ -1179,13 +1179,24 @@ def _remap_parameter_atom_indices(parameter, old_to_new): parameter.indices = (old_to_new[atoms], parindexes) @staticmethod - def _translate_symmetry_metadata(metadata, old_to_new, layer_map, coordinate_shift): + def _transform_symmetry_metadata( + metadata, + old_to_new, + layer_map, + xy_translation, + z_shift_by_layer, + ): if metadata is None: return None metadata = copy.deepcopy(metadata) coordinate_index = {"x": 0, "y": 1, "z": 2} atoms = [] for atom in metadata.atoms: + z_shift = z_shift_by_layer[atom.layer] + coordinate_shift = np.asarray( + [xy_translation[0], xy_translation[1], z_shift], + dtype=np.float64, + ) couplings = [] for coupling in atom.couplings: couplings.append( @@ -1198,10 +1209,15 @@ def _translate_symmetry_metadata(metadata, old_to_new, layer_map, coordinate_shi ), ) ) + parent_shift = metadata.surface_spec.transform @ coordinate_shift atoms.append( dataclasses.replace( atom, atom_index=int(old_to_new[atom.atom_index]), + parent_fractional=( + np.asarray(atom.parent_fractional, dtype=np.float64) + + parent_shift + ), surface_fractional=( np.asarray(atom.surface_fractional, dtype=np.float64) + coordinate_shift @@ -1213,6 +1229,16 @@ def _translate_symmetry_metadata(metadata, old_to_new, layer_map, coordinate_shi metadata.atoms = sorted(atoms, key=lambda atom: atom.atom_index) return metadata + def _layer_positions_for_cycle(self, cycle): + layer_positions = { + layer: float(self.layerpos.get(float(layer), np.nan)) for layer in cycle + } + for layer, position in layer_positions.items(): + if np.isnan(position): + where = self.basis[:, 7] == layer + layer_positions[layer] = float(np.mean(self.basis[where, 3])) + return layer_positions + def translate_layered(self, translation, name=None): """Return a translated copy with cyclic z-layer wrapping. @@ -1224,9 +1250,9 @@ def translate_layered(self, translation, name=None): :class:`CTRstacking.LayerCycle`. A positive z step moves lower layers upward and wraps the former top layer to the bottom. - Symmetry metadata is copied only when atom ordering is unchanged. When - layer cycling regroups atom rows, generated Wyckoff metadata is dropped - because it describes the pre-transform atom indices. + Symmetry metadata is updated by remapping atom indices, atom layers, + surface coordinates, parent coordinates, and Wyckoff coupling + constants. :param translation: Translation vector with shape ``(3,)``, affine matrix with shape @@ -1285,27 +1311,19 @@ def translate_layered(self, translation, name=None): np.mean(self.basis[where, 3]) ) transformed.layerpos[layer_key] += z_translation - coordinate_shift = np.array( - [xy_translation[0], xy_translation[1], z_translation], - dtype=np.float64, - ) old_to_new = np.arange(self.basis.shape[0], dtype=np.intp) layer_map = {layer: layer for layer in cycle} - transformed.symmetry_metadata = self._translate_symmetry_metadata( + z_shift_by_layer = {layer: z_translation for layer in cycle} + transformed.symmetry_metadata = self._transform_symmetry_metadata( self.symmetry_metadata, old_to_new, layer_map, - coordinate_shift, + xy_translation, + z_shift_by_layer, ) return transformed - layer_positions = { - layer: float(self.layerpos.get(float(layer), np.nan)) for layer in cycle - } - for layer, position in layer_positions.items(): - if np.isnan(position): - where = self.basis[:, 7] == layer - layer_positions[layer] = float(np.mean(self.basis[where, 3])) + layer_positions = self._layer_positions_for_cycle(cycle) layer_map = { layer: cycle[(index + z_steps_effective) % layer_count] @@ -1365,19 +1383,410 @@ def transform_basis_values(values): ) if transformed._explicit_layer_cycle is not None: transformed._explicit_layer_cycle = tuple(cycle) - coordinate_shift = np.array( - [xy_translation[0], xy_translation[1], z_translation], - dtype=np.float64, + z_shift_by_layer = {layer: z_translation for layer in cycle} + transformed.symmetry_metadata = self._transform_symmetry_metadata( + self.symmetry_metadata, + old_to_new, + layer_map, + xy_translation, + z_shift_by_layer, ) - transformed.symmetry_metadata = self._translate_symmetry_metadata( + transformed._test_special_formfactors() + return transformed + + def affine_layer_transform(self, translation, name=None): + """Return a canonical layered translation copy. + + The translation input follows :meth:`translate_layered`, but z layer + steps are wrapped back into this unit cell's structural-layer + positions. This preserves a ``0 <= z < 1`` unit-cell representation + for cells whose layer origins are inside that interval. + + Symmetry metadata is updated by remapping atom indices, atom layers, + surface coordinates, parent coordinates, and Wyckoff coupling + constants. + + :param translation: + Translation vector with shape ``(3,)``, affine matrix with shape + ``(3, 4)``, or homogeneous affine matrix with shape ``(4, 4)``. + :param str name: + Optional name for the returned unit cell. Defaults to this unit + cell's name. + :returns: + Transformed unit-cell copy with wrapped layer coordinates. + :rtype: + UnitCell + :raises ValueError: + If the affine linear part is not identity or any requested + translation is not an integer. + """ + translation = self._translation_from_affine_input(translation) + xy_translation = translation[:2] + if not np.allclose(xy_translation, np.rint(xy_translation)): + raise ValueError("x and y affine translations must be integers.") + if not np.isclose(translation[2], np.rint(translation[2])): + raise ValueError("z affine translation must be an integer layer step.") + + xy_translation = np.rint(xy_translation).astype(np.float64) + z_steps = int(np.rint(translation[2])) + + transformed = copy.deepcopy(self) + if name is not None: + transformed.name = name + + if self.basis.size == 0: + if z_steps: + raise ValueError("Cannot cyclically shift layers in an empty UnitCell.") + return transformed + + cycle = self.layer_cycle.layers + layer_count = len(cycle) + z_steps_effective = z_steps % layer_count + if z_steps_effective == 0: + transformed.basis[:, 1:3] += xy_translation + transformed.basis_0[:, 1:3] += xy_translation + if transformed._basis_parvalues is not None: + transformed._basis_parvalues[:, 1:3] += xy_translation + old_to_new = np.arange(self.basis.shape[0], dtype=np.intp) + layer_map = {layer: layer for layer in cycle} + z_shift_by_layer = {layer: 0.0 for layer in cycle} + transformed.symmetry_metadata = self._transform_symmetry_metadata( + self.symmetry_metadata, + old_to_new, + layer_map, + xy_translation, + z_shift_by_layer, + ) + return transformed + + layer_positions = self._layer_positions_for_cycle(cycle) + layer_map = { + layer: cycle[(index + z_steps_effective) % layer_count] + for index, layer in enumerate(cycle) + } + z_shift_by_layer = { + layer: layer_positions[layer_map[layer]] - layer_positions[layer] + for layer in cycle + } + + def transform_basis_values(values): + if values is None or values.size == 0: + return values + values_new = np.array(values, copy=True) + values_new[:, 1:3] += xy_translation + for layer in cycle: + where = values[:, 7] == layer + if not np.any(where): + continue + values_new[where, 3] += z_shift_by_layer[layer] + values_new[where, 7] = layer_map[layer] + return values_new + + basis = transform_basis_values(self.basis) + basis_0 = transform_basis_values(self.basis_0) + basis_parvalues = transform_basis_values(self._basis_parvalues) + + order_lookup = {layer: index for index, layer in enumerate(cycle)} + row_order = np.array( + sorted(range(basis.shape[0]), key=lambda i: (order_lookup[basis[i, 7]], i)), + dtype=np.intp, + ) + old_to_new = np.empty_like(row_order) + old_to_new[row_order] = np.arange(row_order.size) + + transformed.basis = basis[row_order] + transformed.basis_0 = basis_0[row_order] + if basis_parvalues is not None: + transformed._basis_parvalues = basis_parvalues[row_order] + if self.errors is not None: + transformed.errors = np.array(self.errors, copy=True)[row_order] + if self._errors_parvalues is not None: + transformed._errors_parvalues = np.array( + self._errors_parvalues, copy=True + )[row_order] + transformed.names = [self.names[i] for i in row_order] + transformed.dw_increase_constraint = self.dw_increase_constraint[row_order] + if hasattr(self, "f"): + transformed.f = self.f[row_order] + + for parameters in transformed.parameters.values(): + for parameter in parameters: + self._remap_parameter_atom_indices(parameter, old_to_new) + + transformed.layerpos = copy.deepcopy(self.layerpos) + for layer in cycle: + transformed.layerpos[float(layer)] = layer_positions[layer] + if transformed._explicit_layer_cycle is not None: + transformed._explicit_layer_cycle = tuple(cycle) + transformed.symmetry_metadata = self._transform_symmetry_metadata( self.symmetry_metadata, old_to_new, layer_map, - coordinate_shift, + xy_translation, + z_shift_by_layer, + ) + transformed._test_special_formfactors() + return transformed + + def supercell(self, repeats, symmetry="preserve", name=None): + """Return a repeated unit-cell copy. + + :param repeats: + Integer repeat counts along the ``a``, ``b``, and ``c`` lattice + directions. + :param str symmetry: + ``"preserve"`` keeps generated atoms on the original Wyckoff + sites, so a later Wyckoff fit parameter is shared across all + repeated copies. ``"independent"`` creates one copied Wyckoff site + per generated unit-cell repeat, so each copy can be fitted + independently. + :param str name: + Optional name for the returned unit cell. Defaults to this unit + cell's name. + :returns: + Repeated unit-cell copy. + :rtype: + UnitCell + :raises ValueError: + If repeat counts are not positive integers or ``symmetry`` is + unknown. + """ + repeats = np.asarray(repeats, dtype=np.intp) + if repeats.shape != (3,) or np.any(repeats <= 0): + raise ValueError("repeats must contain three positive integers.") + if symmetry not in {"preserve", "independent"}: + raise ValueError("symmetry must be 'preserve' or 'independent'.") + + new_a = self.a * repeats + transformed = UnitCell( + new_a, + np.rad2deg(self.alpha), + name=self.name if name is None else name, + layer_behavior=self.layer_behavior, + ) + transformed.layer_transition = copy.deepcopy(self.layer_transition) + transformed.refHKLTransform = np.copy(self.refHKLTransform) + transformed.refRealTransform = np.copy(self.refRealTransform) + transformed.coherentDomainMatrix = copy.deepcopy(self.coherentDomainMatrix) + transformed.coherentDomainOccupancy = copy.deepcopy( + self.coherentDomainOccupancy + ) + transformed._special_formfactors_present = self._special_formfactors_present + if hasattr(self, "_E"): + transformed._E = self._E + + if self.basis.size == 0: + return transformed + + cycle = self.layer_cycle.layers + layer_positions = self._layer_positions_for_cycle(cycle) + layer_ids = self._supercell_layer_ids(cycle, int(repeats[2])) + + rows = [] + for iz in range(int(repeats[2])): + for iy in range(int(repeats[1])): + for ix in range(int(repeats[0])): + cell_offset = np.asarray([ix, iy, iz], dtype=np.float64) + copy_index = ( + iz * int(repeats[0]) * int(repeats[1]) + + iy * int(repeats[0]) + + ix + ) + for atom_index, row in enumerate(self.basis): + layer = row[7] + layer_index = cycle.index(layer) + new_row = np.array(row, copy=True) + new_row[1:4] = (row[1:4] + cell_offset) / repeats + new_row[7] = layer_ids[(iz, layer)] + rows.append( + { + "basis": new_row, + "basis_0": np.array( + self.basis_0[atom_index], copy=True + ), + "name": self.names[atom_index], + "constraint": self.dw_increase_constraint[atom_index], + "old_atom": atom_index, + "copy_index": copy_index, + "cell_offset": cell_offset, + "layer_order": iz * len(cycle) + layer_index, + } + ) + rows[-1]["basis_0"][1:4] = ( + self.basis_0[atom_index, 1:4] + cell_offset + ) / repeats + rows[-1]["basis_0"][7] = layer_ids[(iz, layer)] + if self.errors is not None: + rows[-1]["errors"] = np.array( + self.errors[atom_index], copy=True + ) + if self._errors_parvalues is not None: + rows[-1]["errors_parvalues"] = np.array( + self._errors_parvalues[atom_index], copy=True + ) + if self._basis_parvalues is not None: + rows[-1]["basis_parvalues"] = np.array( + self._basis_parvalues[atom_index], copy=True + ) + rows[-1]["basis_parvalues"][1:4] = ( + self._basis_parvalues[atom_index, 1:4] + cell_offset + ) / repeats + rows[-1]["basis_parvalues"][7] = layer_ids[(iz, layer)] + if hasattr(self, "f"): + rows[-1]["f"] = np.array(self.f[atom_index], copy=True) + + rows = sorted( + enumerate(rows), + key=lambda item: (item[1]["layer_order"], item[0]), + ) + rows = [row for _, row in rows] + + transformed.basis = np.vstack([row["basis"] for row in rows]) + transformed.basis_0 = np.vstack([row["basis_0"] for row in rows]) + transformed.names = [row["name"] for row in rows] + transformed.dw_increase_constraint = np.asarray( + [row["constraint"] for row in rows], + dtype=np.bool_, + ) + if self.errors is not None: + transformed.errors = np.vstack([row["errors"] for row in rows]) + if self._errors_parvalues is not None: + transformed._errors_parvalues = np.vstack( + [row["errors_parvalues"] for row in rows] + ) + if self._basis_parvalues is not None: + transformed._basis_parvalues = np.vstack( + [row["basis_parvalues"] for row in rows] + ) + if hasattr(self, "f"): + transformed.f = np.vstack([row["f"] for row in rows]) + + transformed.layerpos = {} + for iz in range(int(repeats[2])): + for layer in cycle: + transformed.layerpos[float(layer_ids[(iz, layer)])] = ( + layer_positions[layer] + iz + ) / repeats[2] + transformed._explicit_layer_cycle = tuple( + layer_ids[(iz, layer)] + for iz in range(int(repeats[2])) + for layer in cycle + ) + transformed.parameters = {"absolute": [], "relative": []} + transformed.symmetry_metadata = self._supercell_symmetry_metadata( + rows, + repeats.astype(np.float64), + layer_ids, + symmetry, ) transformed._test_special_formfactors() return transformed + def _supercell_layer_ids(self, cycle, z_repeats): + if z_repeats == 1: + return {(0, layer): layer for layer in cycle} + numeric = np.asarray(cycle, dtype=np.float64) + expected = np.arange(1, len(cycle) + 1, dtype=np.float64) + if np.allclose(numeric, expected): + return { + (iz, layer): float(layer + iz * len(cycle)) + for iz in range(z_repeats) + for layer in cycle + } + return { + (iz, layer): float(iz * len(cycle) + index) + for iz in range(z_repeats) + for index, layer in enumerate(cycle) + } + + def _supercell_symmetry_metadata(self, rows, repeats, layer_ids, symmetry): + if self.symmetry_metadata is None: + return None + metadata = copy.deepcopy(self.symmetry_metadata) + site_by_copy = {} + if symmetry == "preserve": + sites = tuple(copy.deepcopy(metadata.sites)) + else: + sites = [] + for row in rows: + atom = metadata.atom_wyckoff_metadata(row["old_atom"]) + if atom is None: + continue + key = (row["copy_index"], atom.site_id) + if key in site_by_copy: + continue + site = next( + site for site in metadata.sites if site.site_id == atom.site_id + ) + site_id = f"{site.site_id}_copy{row['copy_index']}" + site_by_copy[key] = site_id + sites.append(dataclasses.replace(site, site_id=site_id)) + sites = tuple(sites) + + atoms = [] + coordinate_index = {"x": 0, "y": 1, "z": 2} + for new_index, row in enumerate(rows): + atom = metadata.atom_wyckoff_metadata(row["old_atom"]) + if atom is None: + continue + site_id = atom.site_id + if symmetry == "independent": + site_id = site_by_copy[(row["copy_index"], atom.site_id)] + surface_fractional = ( + np.asarray(atom.surface_fractional, dtype=np.float64) + + row["cell_offset"] + ) / repeats + parent_fractional = ( + np.asarray(atom.parent_fractional, dtype=np.float64) + + metadata.surface_spec.transform @ row["cell_offset"] + ) + couplings = [] + for coupling in atom.couplings: + new_site_id = site_id + couplings.append( + dataclasses.replace( + coupling, + atom_index=new_index, + site_id=new_site_id, + constant=( + coupling.constant + + row["cell_offset"][ + coordinate_index[coupling.coordinate] + ] + ) + / repeats[coordinate_index[coupling.coordinate]], + factor=( + coupling.factor + / repeats[coordinate_index[coupling.coordinate]] + ), + ) + ) + atoms.append( + dataclasses.replace( + atom, + atom_index=new_index, + site_id=site_id, + parent_fractional=parent_fractional, + surface_fractional=surface_fractional, + layer=layer_ids[(int(row["cell_offset"][2]), atom.layer)], + couplings=tuple(couplings), + ) + ) + + metadata.surface_spec = dataclasses.replace( + metadata.surface_spec, + transform=metadata.surface_spec.transform @ np.diag(repeats), + layer_origins=tuple( + (self._layer_positions_for_cycle(self.layer_cycle.layers)[layer] + iz) + / repeats[2] + for iz in range(int(repeats[2])) + for layer in self.layer_cycle.layers + ), + ) + metadata.sites = sites + metadata.atoms = atoms + return metadata + @property def layer_behavior(self): """Return how layer selection is handled during crystal stacking.""" diff --git a/orgui/datautils/xrayutils/test/test_CTRcalc.py b/orgui/datautils/xrayutils/test/test_CTRcalc.py index 2d98a60..260abb7 100644 --- a/orgui/datautils/xrayutils/test/test_CTRcalc.py +++ b/orgui/datautils/xrayutils/test/test_CTRcalc.py @@ -569,6 +569,201 @@ def test_unitcell_translate_layered_stacks_with_existing_layers(self): self.assertEqual(film.end_layer_number, 3) np.testing.assert_array_equal(film.layer_order, [2, 3, 1]) + def test_unitcell_affine_layer_transform_wraps_layers_to_unit_cell(self): + unitcell = self.make_layered_unitcell() + unitcell.basis[:, 1] = [0.1, 0.2, 0.3] + unitcell.basis[:, 3] += 0.05 + unitcell.basis_0[:] = unitcell.basis + + transformed = unitcell.affine_layer_transform([0, 0, 1]) + + np.testing.assert_allclose(transformed.basis[:, 1], [0.3, 0.1, 0.2]) + np.testing.assert_allclose( + transformed.basis[:, 3], + [0.05, 1 / 3 + 0.05, 2 / 3 + 0.05], + ) + np.testing.assert_array_equal(transformed.basis[:, 7], [1, 2, 3]) + self.assertEqual(transformed.layerpos, unitcell.layerpos) + + def test_unitcell_affine_layer_transform_full_cycle_preserves_wrapped_z(self): + unitcell = CTRcalc.UnitCell( + [3.0, 3.0, 6.0], + [90.0, 90.0, 90.0], + name="two_layer", + ) + for layer, z in zip((1, 2), (0.0, 0.5)): + unitcell.addAtom("C", [0.0, 0.0, z], 0.1, 0.1, 1.0, layer=layer) + unitcell.layerpos[float(layer)] = z + + transformed = unitcell.affine_layer_transform([0, 0, -2]) + + np.testing.assert_allclose(transformed.basis, unitcell.basis) + np.testing.assert_allclose(transformed.basis_0, unitcell.basis_0) + self.assertEqual(transformed.layerpos, unitcell.layerpos) + + def test_unitcell_affine_layer_transform_remaps_symmetry_metadata_atoms(self): + unitcell = self.make_layered_unitcell() + unitcell.symmetry_metadata = CTRsymmetry.SurfaceSymmetryModel( + CTRsymmetry.SurfaceCellSpec( + (3.0, 3.0, 6.0), + (90.0, 90.0, 90.0), + np.identity(3), + ), + (), + atoms=[ + CTRsymmetry.GeneratedWyckoffAtom( + atom_index=2, + element="C", + site_id="C_2", + wyckoff_label="1a", + parent_fractional=np.array([0.0, 0.0, 2 / 3]), + surface_fractional=np.array([0.0, 0.0, 2 / 3]), + layer=3, + couplings=( + CTRsymmetry.WyckoffCoupling( + atom_index=2, + coordinate="z", + variable="u", + constant=2 / 3, + factor=1.0, + site_id="C_2", + ), + ), + ), + ], + ) + + transformed = unitcell.affine_layer_transform([0, 0, 1]) + + atom = transformed.atom_wyckoff_metadata(0) + self.assertIsNotNone(atom) + self.assertEqual(atom.layer, 1) + np.testing.assert_allclose(atom.surface_fractional, [0.0, 0.0, 0.0]) + np.testing.assert_allclose(atom.parent_fractional, [0.0, 0.0, 0.0]) + self.assertEqual(atom.couplings[0].atom_index, 0) + self.assertAlmostEqual(atom.couplings[0].constant, 0.0) + + @staticmethod + def make_single_wyckoff_cell(): + unitcell = CTRcalc.UnitCell( + [3.0, 4.0, 5.0], + [90.0, 90.0, 90.0], + name="single", + ) + unitcell.addAtom("C", [0.25, 0.0, 0.0], 0.1, 0.1, 1.0, layer=1) + unitcell.layerpos[1.0] = 0.0 + unitcell.symmetry_metadata = CTRsymmetry.SurfaceSymmetryModel( + CTRsymmetry.SurfaceCellSpec( + (3.0, 4.0, 5.0), + (90.0, 90.0, 90.0), + np.identity(3), + layer_origins=(0.0,), + ), + ( + CTRsymmetry.WyckoffSiteSpec( + site_id="C_1", + element="C", + wyckoff_label="1a", + coordinates=(), + variables={"u": 0.25}, + occ=1.0, + iDW=0.1, + oDW=0.1, + ), + ), + atoms=[ + CTRsymmetry.GeneratedWyckoffAtom( + atom_index=0, + element="C", + site_id="C_1", + wyckoff_label="1a", + parent_fractional=np.array([0.25, 0.0, 0.0]), + surface_fractional=np.array([0.25, 0.0, 0.0]), + layer=1, + couplings=( + CTRsymmetry.WyckoffCoupling( + atom_index=0, + coordinate="x", + variable="u", + constant=0.0, + factor=1.0, + site_id="C_1", + ), + ), + ), + ], + ) + return unitcell + + def test_unitcell_supercell_scales_lattice_and_coordinates(self): + unitcell = self.make_single_wyckoff_cell() + + supercell = unitcell.supercell((2, 3, 1), name="super") + + np.testing.assert_allclose(supercell.a, [6.0, 12.0, 5.0]) + self.assertEqual(supercell.name, "super") + self.assertEqual(supercell.basis.shape[0], 6) + np.testing.assert_allclose( + supercell.basis[:2, 1:4], + [[0.125, 0.0, 0.0], [0.625, 0.0, 0.0]], + ) + np.testing.assert_allclose( + supercell.symmetry_metadata.surface_spec.transform, + np.diag([2.0, 3.0, 1.0]), + ) + self.assertEqual(supercell.parameters, {"absolute": [], "relative": []}) + + def test_unitcell_supercell_repeats_layers_along_z(self): + unitcell = self.make_layered_unitcell() + + supercell = unitcell.supercell((1, 1, 2)) + + np.testing.assert_array_equal(supercell.layers, [1.0, 2.0, 3.0, 4.0, 5.0, 6.0]) + self.assertEqual(set(supercell.layerpos), {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}) + np.testing.assert_allclose( + [supercell.layerpos[layer] for layer in supercell.layers], + [0.0, 1 / 6, 1 / 3, 0.5, 2 / 3, 5 / 6], + ) + + def test_unitcell_supercell_preserves_shared_wyckoff_site(self): + unitcell = self.make_single_wyckoff_cell() + + supercell = unitcell.supercell((2, 1, 1), symmetry="preserve") + + self.assertEqual( + [site["site_id"] for site in supercell.wyckoff_sites()], + ["C_1"], + ) + couplings = supercell.wyckoff_couplings("C_1") + self.assertEqual([coupling.atom_index for coupling in couplings], [0, 1]) + np.testing.assert_allclose( + [coupling.constant for coupling in couplings], + [0.0, 0.5], + ) + np.testing.assert_allclose( + [coupling.factor for coupling in couplings], + [0.5, 0.5], + ) + + supercell.addWyckoffParameter("C_1", "u", limits=(-0.2, 0.2)) + supercell.setFitParameters([0.1]) + + np.testing.assert_allclose(supercell.basis[:, 1], [0.175, 0.675]) + + def test_unitcell_supercell_independent_wyckoff_sites_fit_separately(self): + unitcell = self.make_single_wyckoff_cell() + + supercell = unitcell.supercell((2, 1, 1), symmetry="independent") + + self.assertEqual( + [site["site_id"] for site in supercell.wyckoff_sites()], + ["C_1_copy0", "C_1_copy1"], + ) + supercell.addWyckoffParameter("C_1_copy0", "u", limits=(-0.2, 0.2)) + supercell.setFitParameters([0.1]) + + np.testing.assert_allclose(supercell.basis[:, 1], [0.175, 0.625]) + def test_film_and_poisson_rotate_cyclic_layer_order(self): film = CTRfilm.Film(self.make_layered_unitcell("film")) film.basis[0] = 2 From f50969427f69d229288b2dcf36a411fda8e4f6d8 Mon Sep 17 00:00:00 2001 From: "T. Fuchs" Date: Wed, 8 Jul 2026 23:11:58 -0400 Subject: [PATCH 04/10] feat: generalized Wyckoff parameters, implement fit of parameter --- doc/source/api.rst | 1 + doc/source/wyckoff_fitting.rst | 272 +++++++++++++++ orgui/datautils/xrayutils/CTRcalc.py | 194 +++++++++++ orgui/datautils/xrayutils/CTRfilm.py | 160 +++++++++ orgui/datautils/xrayutils/CTRsymmetry.py | 199 ++++++++++- orgui/datautils/xrayutils/CTRuc.py | 310 +++++++++++++++++- .../datautils/xrayutils/test/test_CTRcalc.py | 190 ++++++++++- .../xrayutils/test/test_CTRsymmetry.py | 227 ++++++++++++- 8 files changed, 1529 insertions(+), 24 deletions(-) create mode 100644 doc/source/wyckoff_fitting.rst diff --git a/doc/source/api.rst b/doc/source/api.rst index fe77ed4..ff93af9 100644 --- a/doc/source/api.rst +++ b/doc/source/api.rst @@ -5,6 +5,7 @@ API :maxdepth: 1 ctr_structure_factors + wyckoff_fitting GUI API ------- diff --git a/doc/source/wyckoff_fitting.rst b/doc/source/wyckoff_fitting.rst new file mode 100644 index 0000000..5abf3bb --- /dev/null +++ b/doc/source/wyckoff_fitting.rst @@ -0,0 +1,272 @@ +Wyckoff Symmetry And Fitting +============================ + +The CTR symmetry helpers make Wyckoff information available as metadata on +ordinary ``UnitCell`` atoms. Structure-factor calculations still use normal +``UnitCell.basis`` coordinates; symmetry is used only for queries, file +persistence, and optional fit-parameter construction. + +The symmetry construction path uses optional PyXtal support. Importing +``CTRcalc``, ``CTRuc``, or ``CTRsymmetry`` does not import PyXtal. PyXtal is +only required when assigning symmetry from a structure seed or querying +Wyckoff tables. Files that already contain saved symmetry metadata can be +loaded and fitted without PyXtal. + +Building A Surface Symmetry Model +--------------------------------- + +Use ``CTRsymmetry.model_from_seed`` to assign parent-crystal Wyckoff sites and +then build a commensurate CTR surface unit cell: + +.. code-block:: python + + import numpy as np + from orgui.datautils.xrayutils import CTRsymmetry + + surface_spec = CTRsymmetry.SurfaceCellSpec( + parent_a=(4.653255, 4.653255, 2.9692), + parent_alpha=(90.0, 90.0, 90.0), + transform=np.asarray( + [ + [1.0, 0.0, 1.0], + [-1.0, 0.0, 1.0], + [0.0, 1.0, 0.0], + ] + ), + layer_origins=(0.0, 0.5), + translation_range=1, + ) + + model = CTRsymmetry.model_from_seed(parent_structure, surface_spec, tol=1e-3) + unitcell = model.build_unitcell("RuO2") + +The ``transform`` columns are surface-cell vectors expressed in parent +conventional fractional coordinates. Parent coordinates are mapped to the +surface cell as: + +.. math:: + + p_{\mathrm{surface}} = S^{-1}(p_{\mathrm{parent}} - p_{\mathrm{origin}}). + +If PyXtal standardizes or reorients the input structure, the surface transform +must be expressed in that PyXtal parent basis. + +Querying Wyckoff Metadata +------------------------- + +``UnitCell`` exposes resolved metadata: + +.. code-block:: python + + unitcell.wyckoff_sites() + unitcell.wyckoff_couplings("O_4f") + unitcell.wyckoff_site_couplings("O_4f") + unitcell.atom_wyckoff_metadata(0) + +``wyckoff_sites()`` returns site dictionaries with the site id, element, +Wyckoff label, free variables, generated atom indices, space group, and status. +The status is one of: + +``metadata_only`` + Symmetry metadata is present, but no Wyckoff-related fit parameter has been + added. + +``symmetry_preserving`` + A true Wyckoff variable fit was added with ``addWyckoffParameter``. + +``site_displaced`` + A representative-site shift was added with ``addWyckoffShift``. This moves + all atoms in the assigned site together, but may lower the crystallographic + site symmetry. + +``partially_overridden`` + A manual atom-coordinate parameter touches atoms in the site. + +Two Fit Modes +------------- + +``addWyckoffParameter`` fits a true Wyckoff variable such as ``u``, ``v``, or +``w``. It preserves the assigned Wyckoff-site symmetry because all equivalent +atom coordinates are updated through the stored affine Wyckoff couplings: + +.. code-block:: python + + unitcell.addWyckoffParameter( + "O_4f", + "u", + absolute_limits=(0.28, 0.34), + ) + +The fitted value is a delta from the stored variable value. For rutile oxygen, +coordinates such as ``u``, ``1-u``, ``0.5+u``, and ``0.5-u`` therefore move +with the correct signs. + +``addWyckoffShift`` fits a parent conventional fractional displacement of the +representative site coordinate: + +.. code-block:: python + + unitcell.addWyckoffShift( + "O_4f", + "x", + limits=(-0.02, 0.02), + ) + +Here ``"x"``, ``"y"``, and ``"z"`` are parent conventional fractional axes, +not generated surface-cell axes. The shift is propagated through the stored +space-group operation for each generated atom and then through the surface-cell +transform. This intentionally keeps the original atom assignment and atom count +while allowing the site to move away from its exact Wyckoff constraint. + +Plural helpers are available for fitting several variables or shifts: + +.. code-block:: python + + unitcell.addWyckoffParameters("Mg_8d") # all free variables + unitcell.addWyckoffShifts("Mg_8d", axes=("x", "z")) + +Container Classes +----------------- + +``Film`` and ``PoissonSurface`` forward Wyckoff helper calls to their single +template unit cell: + +.. code-block:: python + + film.addWyckoffParameter("O_4f", "u", limits=(-0.05, 0.05)) + poisson.addWyckoffShift("O_4f", "x", limits=(-0.02, 0.02)) + +``EpitaxyInterface`` follows the same selector convention as its existing +``addFitParameter`` API. Select the internal unit cell with ``unitcell``: + +.. code-block:: python + + interface.addWyckoffParameter( + "O_4f", + "u", + limits=(-0.05, 0.05), + unitcell="top", + ) + + interface.addWyckoffShift( + "O_4f", + "x", + limits=(-0.02, 0.02), + unitcell=["top", "bottom"], + ) + +When a list is supplied, one parameter is added to each selected internal unit +cell. They remain separate parameters unless a higher-level ``SXRDCrystal`` +coupled parameter links them. + +Linking Across Unit Cells +------------------------- + +``SXRDCrystal`` provides the same two Wyckoff helper modes at the crystal +level. These helpers create one coupled crystal parameter and ordinary +relative fit parameters inside each selected unit cell. + +For a symmetry-preserving Wyckoff variable: + +.. code-block:: python + + crystal.addWyckoffParameter( + { + "RuO2": ("O_4f", "u"), + "TiO2": ("O_4f", "u"), + }, + name="shared_rutile_oxygen_u", + limits=(-0.05, 0.05), + ) + +For a representative-site shift: + +.. code-block:: python + + crystal.addWyckoffShift( + { + "RuO2": ("O_4f", "x"), + "TiO2": ("O_4f", "x"), + }, + name="shared_oxygen_x_shift", + limits=(-0.02, 0.02), + ) + +For ``EpitaxyInterface`` components, select the internal unit cell in the +value dictionary, or use a ``(component, unitcell)`` key: + +.. code-block:: python + + crystal.addWyckoffParameter( + { + "interface": { + "site_id": "O_4f", + "variable": "u", + "unitcell": ("top", "bottom"), + } + }, + name="shared_interface_u", + limits=(-0.05, 0.05), + ) + + crystal.addWyckoffShift( + {("interface", "top"): ("O_4f", "x")}, + name="top_interface_oxygen_shift", + limits=(-0.02, 0.02), + ) + +Persistence +----------- + +Symmetry-bearing ``.xtal`` and ``.xpr`` files store resolved plain-text tables +after the layer metadata. The persisted information includes: + +* space group number and symbol; +* surface transform and origin; +* Wyckoff site ids, labels, variables, and representative parent coordinates; +* generated atom metadata; +* Wyckoff variable couplings; +* representative-site shift couplings. + +Reloaded files can be queried and fitted with ``addWyckoffParameter`` and +``addWyckoffShift`` without PyXtal because the resolved couplings are already +present in the file. + +API Reference +------------- + +.. automodule:: orgui.datautils.xrayutils.CTRsymmetry + :members: SurfaceCellSpec, SurfaceSymmetryModel, possible_wyckoff_positions, + sites_from_seed, model_from_seed, surface_unitcell_from_seed + :member-order: bysource + :no-index: + +.. autoclass:: orgui.datautils.xrayutils.CTRuc.UnitCell + :members: wyckoff_sites, wyckoff_couplings, wyckoff_site_couplings, + atom_wyckoff_metadata, addWyckoffParameter, addWyckoffParameters, + addWyckoffShift, addWyckoffShifts + :member-order: bysource + :no-index: + +.. autoclass:: orgui.datautils.xrayutils.CTRfilm.Film + :members: addWyckoffParameter, addWyckoffParameters, addWyckoffShift, + addWyckoffShifts + :member-order: bysource + :no-index: + +.. autoclass:: orgui.datautils.xrayutils.CTRfilm.EpitaxyInterface + :members: addWyckoffParameter, addWyckoffParameters, addWyckoffShift, + addWyckoffShifts + :member-order: bysource + :no-index: + +.. autoclass:: orgui.datautils.xrayutils.CTRfilm.PoissonSurface + :members: addWyckoffParameter, addWyckoffParameters, addWyckoffShift, + addWyckoffShifts + :member-order: bysource + :no-index: + +.. autoclass:: orgui.datautils.xrayutils.CTRcalc.SXRDCrystal + :members: addWyckoffParameter, addWyckoffShift + :member-order: bysource + :no-index: diff --git a/orgui/datautils/xrayutils/CTRcalc.py b/orgui/datautils/xrayutils/CTRcalc.py index 287dbc4..e453877 100644 --- a/orgui/datautils/xrayutils/CTRcalc.py +++ b/orgui/datautils/xrayutils/CTRcalc.py @@ -502,6 +502,200 @@ def addFitParameter(self, parameter, limits=(-np.inf, np.inf), **keyargs): return par + def addWyckoffParameter(self, parameter, limits=(-np.inf, np.inf), **keyargs): + """Add a coupled symmetry-preserving Wyckoff parameter. + + ``parameter`` maps crystal component names to Wyckoff selections. A + selection can be ``(site_id, variable)`` or a dictionary with + ``site_id`` and ``variable`` entries. For ``EpitaxyInterface`` + components, either use a key ``(component, unitcell)`` or provide + ``unitcell="top"``, ``unitcell="bottom"``, or a list in the selection + dictionary. + + :param dict parameter: + Component-specific Wyckoff variable selections. + :param tuple limits: + Shared delta limits in fractional units. + :param keyargs: + Additional coupled-parameter settings such as ``name`` or + ``prior``. + :returns: + Created coupled parameter. + :rtype: + Parameter + """ + return self._add_wyckoff_coupled_parameter( + parameter, + kind="coordinate", + value_key="variable", + coupling_getter="wyckoff_couplings", + limits=limits, + **keyargs, + ) + + def addWyckoffShift(self, parameter, limits=(-np.inf, np.inf), **keyargs): + """Add a coupled representative Wyckoff-site shift parameter. + + ``parameter`` maps crystal component names to Wyckoff shift + selections. A selection can be ``(site_id, axis)`` or a dictionary with + ``site_id`` and ``axis`` entries. For ``EpitaxyInterface`` components, + either use a key ``(component, unitcell)`` or provide + ``unitcell="top"``, ``unitcell="bottom"``, or a list in the selection + dictionary. + + :param dict parameter: + Component-specific Wyckoff shift selections. + :param tuple limits: + Shared parent-coordinate delta limits in fractional units. + :param keyargs: + Additional coupled-parameter settings such as ``name`` or + ``prior``. + :returns: + Created coupled parameter. + :rtype: + Parameter + """ + return self._add_wyckoff_coupled_parameter( + parameter, + kind="site_displacement", + value_key="axis", + coupling_getter="wyckoff_site_couplings", + limits=limits, + **keyargs, + ) + + def _add_wyckoff_coupled_parameter( + self, + parameter, + kind, + value_key, + coupling_getter, + limits=(-np.inf, np.inf), + **keyargs, + ): + if not parameter: + raise ValueError("At least one Wyckoff target must be provided.") + + prior = keyargs.get("prior", None) + name = keyargs.get("name", f"{self.name}_par_{self._parIdNo}") + keyargs["name"] = name + keyargs["prior"] = prior + + component_indices = [] + for component_key, selection in parameter.items(): + for component_index, unitcell, site_id, value in self._wyckoff_targets( + component_key, + selection, + value_key, + ): + couplings = [ + coupling + for coupling in getattr(unitcell, coupling_getter)(site_id) + if getattr(coupling, value_key) == value + ] + if not couplings: + raise ValueError( + f"No Wyckoff {value_key} couplings found for " + f"{component_key!r}, site {site_id!r}, {value_key} " + f"{value!r}." + ) + atoms = np.asarray( + [coupling.atom_index for coupling in couplings], + dtype=np.intp, + ) + coordinates = tuple(coupling.coordinate for coupling in couplings) + factors = np.asarray( + [coupling.factor for coupling in couplings], + dtype=np.float64, + ) + settings = dict(keyargs) + settings["wyckoff"] = { + "site_id": site_id, + "kind": kind, + value_key: value, + "value_kind": "delta", + } + unitcell.addRelParameter( + (atoms, coordinates), + factors, + limits, + **settings, + ) + component_indices.append(component_index) + + par = Parameter( + name, + np.asarray(component_indices, dtype=np.intp), + ParameterType.RELATIVE, + limits, + prior, + keyargs, + ) + self.parameters["coupled"].append(par) + self._parIdNo += 1 + self.fit_metadata_cache = None + return par + + def _wyckoff_targets(self, component_key, selection, value_key): + internal = None + if isinstance(component_key, tuple): + if len(component_key) != 2: + raise ValueError( + "Tuple Wyckoff component keys must be " + "(component, unitcell)." + ) + component_key, internal = component_key + + if isinstance(selection, dict): + site_id = selection["site_id"] + value = selection[value_key] + internal = selection.get("unitcell", internal) + else: + site_id, value = selection + + component = self[component_key] + component_index = self.getUcIndex(component_key) + if isinstance(internal, list | tuple): + for internal_name in internal: + yield ( + component_index, + self._wyckoff_target_unitcell(component, internal_name), + site_id, + value, + ) + else: + yield ( + component_index, + self._wyckoff_target_unitcell(component, internal), + site_id, + value, + ) + + @staticmethod + def _wyckoff_target_unitcell(component, internal): + if isinstance(component, UnitCell): + if internal is not None: + raise ValueError("UnitCell components do not use 'unitcell'.") + return component + if isinstance(component, Film | PoissonSurface): + if internal is not None: + return component[internal] + return component.unitcell + if isinstance(component, EpitaxyInterface): + if internal is None: + raise ValueError( + "EpitaxyInterface Wyckoff parameters require " + "unitcell='top' or unitcell='bottom'." + ) + return component[internal] + if internal is not None and hasattr(component, "__getitem__"): + return component[internal] + if hasattr(component, "unitcell"): + return component.unitcell + raise TypeError( + f"Component {component!r} does not expose Wyckoff-bearing unit cells." + ) + def getSurfaceBasis(self): return np.concatenate([uc.basis for uc in self if isinstance(uc, UnitCell)]) diff --git a/orgui/datautils/xrayutils/CTRfilm.py b/orgui/datautils/xrayutils/CTRfilm.py index 8063a02..4151087 100644 --- a/orgui/datautils/xrayutils/CTRfilm.py +++ b/orgui/datautils/xrayutils/CTRfilm.py @@ -125,6 +125,166 @@ def layer_behaviour(self): def layer_behaviour(self, behavior): self.layer_behavior = behavior + def _wyckoff_target_unitcells(self, kwarg): + if hasattr(self, "uc_top") and hasattr(self, "uc_bottom"): + if "unitcell" not in kwarg: + raise ValueError( + "Missing unit cell name. Provide unit cell name as kwarg " + "'unitcell'" + ) + unitcell = kwarg["unitcell"] + if isinstance(unitcell, list): + return [self[ucn] for ucn in unitcell] + return [self[unitcell]] + return [self.unitcell] + + def addWyckoffParameter( + self, + site_id, + variable, + limits=(-np.inf, np.inf), + absolute_limits=None, + **kwarg, + ): + """Add a Wyckoff variable parameter on contained unit cells. + + :param str site_id: + Wyckoff site identifier. + :param str variable: + Wyckoff variable name, for example ``"u"``. + :param tuple limits: + Delta fit limits in fractional units. + :param tuple absolute_limits: + Optional absolute variable limits. + :param kwarg: + For ``EpitaxyInterface``, provide ``unitcell="top"``, + ``unitcell="bottom"``, or a list of these names. + :returns: + Created parameter or list of parameters. + """ + parameters = [ + unitcell.addWyckoffParameter( + site_id, + variable, + limits=limits, + absolute_limits=absolute_limits, + **kwarg, + ) + for unitcell in self._wyckoff_target_unitcells(kwarg) + ] + return parameters if len(parameters) > 1 else parameters[0] + + def addWyckoffParameters( + self, + site_id, + variables=None, + limits=(-np.inf, np.inf), + absolute_limits=None, + **kwarg, + ): + """Add several Wyckoff variable parameters on contained unit cells. + + :param str site_id: + Wyckoff site identifier. + :param variables: + Iterable of Wyckoff variable names. If ``None``, fit all variables. + :param tuple limits: + Delta fit limits in fractional units. + :param absolute_limits: + Optional absolute variable limits. + :param kwarg: + For ``EpitaxyInterface``, provide ``unitcell="top"``, + ``unitcell="bottom"``, or a list of these names. + :returns: + Created parameters. A list of lists is returned when multiple + contained unit cells are selected. + """ + parameters = [ + unitcell.addWyckoffParameters( + site_id, + variables=variables, + limits=limits, + absolute_limits=absolute_limits, + **kwarg, + ) + for unitcell in self._wyckoff_target_unitcells(kwarg) + ] + return parameters if len(parameters) > 1 else parameters[0] + + def addWyckoffShift( + self, + site_id, + axis, + limits=(-np.inf, np.inf), + absolute_limits=None, + **kwarg, + ): + """Add a representative Wyckoff-site shift on contained unit cells. + + :param str site_id: + Wyckoff site identifier. + :param str axis: + Parent conventional representative coordinate, one of ``"x"``, + ``"y"``, or ``"z"``. + :param tuple limits: + Delta fit limits in parent fractional units. + :param tuple absolute_limits: + Optional absolute parent-coordinate limits. + :param kwarg: + For ``EpitaxyInterface``, provide ``unitcell="top"``, + ``unitcell="bottom"``, or a list of these names. + :returns: + Created parameter or list of parameters. + """ + parameters = [ + unitcell.addWyckoffShift( + site_id, + axis, + limits=limits, + absolute_limits=absolute_limits, + **kwarg, + ) + for unitcell in self._wyckoff_target_unitcells(kwarg) + ] + return parameters if len(parameters) > 1 else parameters[0] + + def addWyckoffShifts( + self, + site_id, + axes=("x", "y", "z"), + limits=(-np.inf, np.inf), + absolute_limits=None, + **kwarg, + ): + """Add representative Wyckoff-site shifts on contained unit cells. + + :param str site_id: + Wyckoff site identifier. + :param axes: + Parent conventional representative coordinate axes. + :param tuple limits: + Delta fit limits in parent fractional units. + :param absolute_limits: + Optional absolute parent-coordinate limits. + :param kwarg: + For ``EpitaxyInterface``, provide ``unitcell="top"``, + ``unitcell="bottom"``, or a list of these names. + :returns: + Created parameters. A list of lists is returned when multiple + contained unit cells are selected. + """ + parameters = [ + unitcell.addWyckoffShifts( + site_id, + axes=axes, + limits=limits, + absolute_limits=absolute_limits, + **kwarg, + ) + for unitcell in self._wyckoff_target_unitcells(kwarg) + ] + return parameters if len(parameters) > 1 else parameters[0] + @property def start_layer_number(self): """Return the first cyclic layer created by this object.""" diff --git a/orgui/datautils/xrayutils/CTRsymmetry.py b/orgui/datautils/xrayutils/CTRsymmetry.py index 9da4e6c..274612a 100644 --- a/orgui/datautils/xrayutils/CTRsymmetry.py +++ b/orgui/datautils/xrayutils/CTRsymmetry.py @@ -47,6 +47,7 @@ "wyckoff_sites:", "wyckoff_atoms:", "wyckoff_couplings:", + "wyckoff_site_couplings:", } @@ -233,6 +234,8 @@ class WyckoffSiteSpec: element: str wyckoff_label: str coordinates: tuple[tuple[AffineExpression, AffineExpression, AffineExpression], ...] + representative_parent_fractional: tuple[float, float, float] | None = None + operation_matrices: tuple[np.ndarray, ...] = () variables: dict[str, float] = field(default_factory=dict) occ: float = 1.0 iDW: float = 0.5 @@ -295,6 +298,31 @@ def value(self, variable_value): return self.constant + self.factor * variable_value +@dataclass(frozen=True) +class WyckoffSiteCoupling: + """Describe one atom-coordinate dependency on site representative motion. + + :param int atom_index: + Index of the generated atom in ``UnitCell.basis``. + :param str coordinate: + Surface coordinate name, one of ``"x"``, ``"y"``, or ``"z"``. + :param str axis: + Parent conventional representative coordinate, one of ``"x"``, + ``"y"``, or ``"z"``. + :param float factor: + Surface-coordinate coefficient for a parent representative-coordinate + displacement. + :param str site_id: + Site identifier owning the coupling. + """ + + atom_index: int + coordinate: str + axis: str + factor: float + site_id: str + + @dataclass(frozen=True) class GeneratedWyckoffAtom: """Metadata for one generated surface-cell atom. @@ -326,6 +354,7 @@ class GeneratedWyckoffAtom: surface_fractional: np.ndarray layer: int couplings: tuple[WyckoffCoupling, ...] = () + site_couplings: tuple[WyckoffSiteCoupling, ...] = () @dataclass @@ -363,6 +392,11 @@ def wyckoff_sites(self, parameters=None): "element": site.element, "wyckoff_label": site.wyckoff_label, "variables": variables, + "representative_parent_fractional": ( + None + if site.representative_parent_fractional is None + else tuple(site.representative_parent_fractional) + ), "atom_indices": atoms, "spacegroup_number": self.spacegroup_number, "spacegroup_symbol": self.spacegroup_symbol, @@ -387,6 +421,22 @@ def wyckoff_couplings(self, site_id=None): couplings.extend(atom.couplings) return couplings + def wyckoff_site_couplings(self, site_id=None): + """Return site-displacement couplings for generated atom coordinates. + + :param str site_id: + Optional site identifier. If omitted, all couplings are returned. + :returns: + List of :class:`WyckoffSiteCoupling` objects. + :rtype: + list + """ + couplings = [] + for atom in self.atoms: + if site_id is None or atom.site_id == site_id: + couplings.extend(atom.site_couplings) + return couplings + def atom_wyckoff_metadata(self, atom_index): """Return metadata for one generated atom. @@ -428,6 +478,7 @@ def symmetry_status(self, site_id, parameters=None): return "metadata_only" preserving = False + displaced = False for kind in ("absolute", "relative"): for parameter in parameters.get(kind, []): pairs = set(zip(parameter.indices[0], parameter.indices[1])) @@ -437,9 +488,16 @@ def symmetry_status(self, site_id, parameters=None): and wyckoff_settings.get("site_id") == site_id and pairs <= site_pairs ): - preserving = True + if wyckoff_settings.get("kind") == "coordinate": + preserving = True + elif wyckoff_settings.get("kind") == "site_displacement": + displaced = True + else: + return "partially_overridden" elif pairs & site_pairs: return "partially_overridden" + if displaced: + return "site_displaced" if preserving: return "symmetry_preserving" return "metadata_only" @@ -486,6 +544,7 @@ def build_unitcell(self, name, layer_behavior="ignore"): atom.surface_fractional, atom.layer, atom.couplings, + atom.site_couplings, ) for atom in generated ] @@ -713,7 +772,7 @@ def generate_surface_atoms(surface_spec, sites, tolerance=1e-8): z_min, z_max = surface_spec.z_bounds for site in sites: - for coordinate in site.coordinates: + for coordinate_index, coordinate in enumerate(site.coordinates): for translation in surface_spec.parent_translations(): parent_exprs = tuple( expression.shifted(translation[index]) @@ -764,6 +823,28 @@ def generate_surface_atoms(surface_spec, sites, tolerance=1e-8): site_id=site.site_id, ) ) + site_couplings = [] + if coordinate_index < len(site.operation_matrices): + operation = np.asarray( + site.operation_matrices[coordinate_index], + dtype=np.float64, + ) + site_displacements = inverse @ operation[:3, :3] + for parent_axis, axis_name in enumerate(COORDINATE_NAMES): + for surface_axis, coordinate_name in enumerate( + COORDINATE_NAMES + ): + factor = site_displacements[surface_axis, parent_axis] + if not math.isclose(factor, 0.0, abs_tol=tolerance): + site_couplings.append( + WyckoffSiteCoupling( + atom_index=atom_index, + coordinate=coordinate_name, + axis=axis_name, + factor=float(factor), + site_id=site.site_id, + ) + ) generated.append( GeneratedWyckoffAtom( atom_index=atom_index, @@ -780,6 +861,7 @@ def generate_surface_atoms(surface_spec, sites, tolerance=1e-8): surface_fractional=wrapped_values, layer=layer, couplings=tuple(couplings), + site_couplings=tuple(site_couplings), ) ) @@ -812,6 +894,16 @@ def generate_surface_atoms(surface_spec, sites, tolerance=1e-8): ) for coupling in atom.couplings ), + site_couplings=tuple( + WyckoffSiteCoupling( + atom_index=index, + coordinate=coupling.coordinate, + axis=coupling.axis, + factor=coupling.factor, + site_id=coupling.site_id, + ) + for coupling in atom.site_couplings + ), ) for index, atom in enumerate(generated) ] @@ -836,9 +928,18 @@ def symmetry_metadata_to_lines(model): lines.append(" " + _format_values(row)) lines.append("surface_origin: " + _format_values(model.surface_spec.origin)) lines.append("wyckoff_sites:") - lines.append(" site_id element wyckoff_label variables occ iDW oDW") + lines.append( + " site_id element wyckoff_label variables " + "representative_x representative_y representative_z occ iDW oDW" + ) for site in model.sites: variables = _format_variables(site.variables) + representative = site.representative_parent_fractional + representative_values = ( + ["-", "-", "-"] + if representative is None + else _format_values(representative).split() + ) lines.append( " " + " ".join( @@ -847,6 +948,7 @@ def symmetry_metadata_to_lines(model): site.element, site.wyckoff_label, variables, + *representative_values, _format_float(site.occ), _format_float(site.iDW), _format_float(site.oDW), @@ -885,6 +987,21 @@ def symmetry_metadata_to_lines(model): ] ) ) + lines.append("wyckoff_site_couplings:") + lines.append(" atom_index coordinate site_id axis factor") + for coupling in model.wyckoff_site_couplings(): + lines.append( + " " + + " ".join( + [ + str(coupling.atom_index), + coupling.coordinate, + coupling.site_id, + coupling.axis, + _format_float(coupling.factor), + ] + ) + ) return lines @@ -914,6 +1031,7 @@ def symmetry_metadata_from_lines(lines, unitcell=None): site_specs = [] atoms = [] couplings_by_atom = {} + site_couplings_by_atom = {} sections = _collect_sections(cleaned) if "spacegroup" in sections: @@ -937,16 +1055,27 @@ def symmetry_metadata_from_lines(lines, unitcell=None): for row in _data_rows(sections.get("wyckoff_sites", [])): parts = row.split() variables = _parse_variables(parts[3]) + if len(parts) >= 10: + representative = ( + None + if parts[4] == "-" + else tuple(float(value) for value in parts[4:7]) + ) + occ, iDW, oDW = (float(value) for value in parts[7:10]) + else: + representative = None + occ, iDW, oDW = (float(value) for value in parts[4:7]) site_specs.append( WyckoffSiteSpec( site_id=parts[0], element=parts[1], wyckoff_label=parts[2], coordinates=(), + representative_parent_fractional=representative, variables=variables, - occ=float(parts[4]), - iDW=float(parts[5]), - oDW=float(parts[6]), + occ=occ, + iDW=iDW, + oDW=oDW, ) ) @@ -962,6 +1091,17 @@ def symmetry_metadata_from_lines(lines, unitcell=None): ) couplings_by_atom.setdefault(coupling.atom_index, []).append(coupling) + for row in _data_rows(sections.get("wyckoff_site_couplings", [])): + parts = row.split() + coupling = WyckoffSiteCoupling( + atom_index=int(parts[0]), + coordinate=parts[1], + site_id=parts[2], + axis=parts[3], + factor=float(parts[4]), + ) + site_couplings_by_atom.setdefault(coupling.atom_index, []).append(coupling) + for row in _data_rows(sections.get("wyckoff_atoms", [])): parts = row.split() atom_index = int(parts[0]) @@ -975,6 +1115,7 @@ def symmetry_metadata_from_lines(lines, unitcell=None): surface_fractional=np.asarray(parts[7:10], dtype=np.float64), layer=int(float(parts[10])), couplings=tuple(couplings_by_atom.get(atom_index, ())), + site_couplings=tuple(site_couplings_by_atom.get(atom_index, ())), ) ) @@ -1023,6 +1164,7 @@ def _import_pyxtal_group(): def _site_from_pyxtal_site(atom_site, variable_names, iDW, oDW, occ): wp = atom_site.wp + representative = _wrap_fractional(atom_site.position) free_values = list(wp.get_free_xyzs(atom_site.position)) if len(free_values) > len(variable_names): raise ValueError( @@ -1046,12 +1188,22 @@ def _site_from_pyxtal_site(atom_site, variable_names, iDW, oDW, occ): _apply_operation_to_expressions(operation.affine_matrix, primary) for operation in wp.ops ) + operation_matrices = _site_operation_matrices( + wp, + representative, + coordinates, + variables, + ) element = str(atom_site.specie) return WyckoffSiteSpec( site_id=f"{element}_{wp.get_label()}", element=element, wyckoff_label=wp.get_label(), coordinates=coordinates, + representative_parent_fractional=tuple( + float(value) for value in representative + ), + operation_matrices=operation_matrices, variables=variables, occ=occ, iDW=iDW, @@ -1059,6 +1211,30 @@ def _site_from_pyxtal_site(atom_site, variable_names, iDW, oDW, occ): ) +def _site_operation_matrices(wp, representative, coordinates, variables): + group = _import_pyxtal_group()(wp.number) + general_operations = group[0].ops + operation_matrices = [] + for coordinate in coordinates: + target = _wrap_fractional( + [expression.evaluate(variables) for expression in coordinate] + ) + for operation in general_operations: + matrix = np.asarray(operation.affine_matrix, dtype=np.float64) + generated = _wrap_fractional( + matrix[:3, :3] @ representative + matrix[:3, 3] + ) + if _fractional_positions_close(generated, target): + operation_matrices.append(matrix) + break + else: + raise ValueError( + f"Could not match a full symmetry operation for Wyckoff " + f"position {wp.get_label()} at {target}." + ) + return tuple(operation_matrices) + + def _primary_position_affine(wp, free_values, variable_names): if not free_values: return np.asarray(wp.get_position_from_free_xyzs([]), dtype=np.float64), ( @@ -1129,6 +1305,17 @@ def _wrap_surface_expressions(values, expressions): return wrapped_values, tuple(wrapped_expressions) +def _wrap_fractional(values): + wrapped = np.mod(np.asarray(values, dtype=np.float64), 1.0) + wrapped[np.isclose(wrapped, 1.0, atol=1e-12)] = 0.0 + return wrapped + + +def _fractional_positions_close(first, second, tolerance=1e-6): + diff = np.mod(np.asarray(first) - np.asarray(second) + 0.5, 1.0) - 0.5 + return np.allclose(diff, 0.0, atol=tolerance) + + def _assign_layer(z_value, layer_origins): if not layer_origins: return 0 diff --git a/orgui/datautils/xrayutils/CTRuc.py b/orgui/datautils/xrayutils/CTRuc.py index d609a1e..336f5bd 100644 --- a/orgui/datautils/xrayutils/CTRuc.py +++ b/orgui/datautils/xrayutils/CTRuc.py @@ -1209,6 +1209,14 @@ def _transform_symmetry_metadata( ), ) ) + site_couplings = [] + for coupling in atom.site_couplings: + site_couplings.append( + dataclasses.replace( + coupling, + atom_index=int(old_to_new[coupling.atom_index]), + ) + ) parent_shift = metadata.surface_spec.transform @ coordinate_shift atoms.append( dataclasses.replace( @@ -1224,6 +1232,7 @@ def _transform_symmetry_metadata( ), layer=layer_map[atom.layer], couplings=tuple(couplings), + site_couplings=tuple(site_couplings), ) ) metadata.atoms = sorted(atoms, key=lambda atom: atom.atom_index) @@ -1720,7 +1729,19 @@ def _supercell_symmetry_metadata(self, rows, repeats, layer_ids, symmetry): ) site_id = f"{site.site_id}_copy{row['copy_index']}" site_by_copy[key] = site_id - sites.append(dataclasses.replace(site, site_id=site_id)) + representative = site.representative_parent_fractional + if representative is not None: + representative = tuple( + np.asarray(representative, dtype=np.float64) + + metadata.surface_spec.transform @ row["cell_offset"] + ) + sites.append( + dataclasses.replace( + site, + site_id=site_id, + representative_parent_fractional=representative, + ) + ) sites = tuple(sites) atoms = [] @@ -1761,6 +1782,19 @@ def _supercell_symmetry_metadata(self, rows, repeats, layer_ids, symmetry): ), ) ) + site_couplings = [] + for coupling in atom.site_couplings: + site_couplings.append( + dataclasses.replace( + coupling, + atom_index=new_index, + site_id=site_id, + factor=( + coupling.factor + / repeats[coordinate_index[coupling.coordinate]] + ), + ) + ) atoms.append( dataclasses.replace( atom, @@ -1770,6 +1804,7 @@ def _supercell_symmetry_metadata(self, rows, repeats, layer_ids, symmetry): surface_fractional=surface_fractional, layer=layer_ids[(int(row["cell_offset"][2]), atom.layer)], couplings=tuple(couplings), + site_couplings=tuple(site_couplings), ) ) @@ -2052,6 +2087,20 @@ def wyckoff_couplings(self, site_id=None): return [] return self.symmetry_metadata.wyckoff_couplings(site_id) + def wyckoff_site_couplings(self, site_id=None): + """Return couplings for representative Wyckoff-site displacement. + + :param str site_id: + Optional site identifier. If omitted, all couplings are returned. + :returns: + List of site-displacement coupling metadata objects. + :rtype: + list + """ + if self.symmetry_metadata is None: + return [] + return self.symmetry_metadata.wyckoff_site_couplings(site_id) + def atom_wyckoff_metadata(self, atom_index): """Return symmetry metadata for one atom. @@ -2074,12 +2123,13 @@ def addWyckoffParameter( absolute_limits=None, **keyargs, ): - """Add a symmetry-preserving relative fit parameter for a Wyckoff variable. + """Add a symmetry-preserving fit parameter for a Wyckoff variable. The stored fit value is the change in the Wyckoff variable from the generated coordinates. For example, fitting rutile oxygen ``u`` adds ``factor * delta_u`` to every generated coordinate that depends on - ``u``. + ``u``. Multi-variable Wyckoff sites are fitted by adding one parameter + per independent coordinate variable. :param str site_id: Site identifier returned by :meth:`wyckoff_sites`. @@ -2099,21 +2149,26 @@ def addWyckoffParameter( """ if self.symmetry_metadata is None: raise ValueError(f"UnitCell {self.name} has no symmetry metadata.") + site = next( + ( + site + for site in self.wyckoff_sites() + if site["site_id"] == site_id + ), + None, + ) + if site is None: + raise ValueError(f"Unknown Wyckoff site {site_id}.") + if not site["variables"]: + raise ValueError( + f"Wyckoff site {site_id} has no positional coordinate variables." + ) + if variable not in site["variables"]: + raise ValueError(f"Wyckoff site {site_id} has no variable {variable}.") + if absolute_limits is not None: if limits != (-np.inf, np.inf): raise ValueError("Use either limits or absolute_limits, not both.") - site = next( - ( - site - for site in self.wyckoff_sites() - if site["site_id"] == site_id - ), - None, - ) - if site is None or variable not in site["variables"]: - raise ValueError( - f"Wyckoff site {site_id} has no variable {variable}." - ) variable_value = site["variables"][variable] limits = ( absolute_limits[0] - variable_value, @@ -2140,11 +2195,15 @@ def addWyckoffParameter( dtype=np.float64, ) - keyargs.setdefault("name", f"{self.name} {site_id}_{variable}_wyckoff") + keyargs.setdefault( + "name", + f"{self.name} {site_id}_{variable}_wyckoff_coordinate", + ) settings = keyargs.setdefault("wyckoff", {}) settings.update( { "site_id": site_id, + "kind": "coordinate", "variable": variable, "value_kind": "delta", } @@ -2156,6 +2215,225 @@ def addWyckoffParameter( **keyargs, ) + def addWyckoffParameters( + self, + site_id, + variables=None, + limits=(-np.inf, np.inf), + absolute_limits=None, + **keyargs, + ): + """Add symmetry-preserving fit parameters for Wyckoff variables. + + :param str site_id: + Site identifier returned by :meth:`wyckoff_sites`. + :param variables: + Iterable of coordinate variables to fit. If ``None``, all free + variables on the site are fitted. + :type variables: + iterable or None + :param limits: + Either one ``(lower, upper)`` tuple applied to every variable or a + dictionary mapping variable names to delta limits. + :param absolute_limits: + Optional dictionary mapping variable names to absolute limits. + :returns: + Created relative fit parameters in variable order. + :rtype: + list + :raises ValueError: + If the site has no free coordinate variables. + """ + sites = {site["site_id"]: site for site in self.wyckoff_sites()} + if site_id not in sites: + raise ValueError(f"Unknown Wyckoff site {site_id}.") + site_variables = tuple(sites[site_id]["variables"]) + if not site_variables: + raise ValueError( + f"Wyckoff site {site_id} has no positional coordinate variables." + ) + if variables is None: + variables = site_variables + variables = tuple(variables) + + def limit_for(limit_spec, variable_name): + if isinstance(limit_spec, dict): + return limit_spec.get(variable_name, (-np.inf, np.inf)) + return limit_spec + + parameters = [] + for variable in variables: + parameter_keyargs = copy.deepcopy(keyargs) + parameter_limits = limit_for(limits, variable) + parameter_absolute_limits = ( + None + if absolute_limits is None + else limit_for(absolute_limits, variable) + ) + parameters.append( + self.addWyckoffParameter( + site_id, + variable, + limits=parameter_limits, + absolute_limits=parameter_absolute_limits, + **parameter_keyargs, + ) + ) + return parameters + + def addWyckoffShift( + self, + site_id, + axis, + limits=(-np.inf, np.inf), + absolute_limits=None, + **keyargs, + ): + """Add a symmetry-lowering shift for representative site motion. + + ``axis`` is a parent conventional fractional coordinate of the + representative atom. The stored fit value is a delta from the + representative coordinate; generated atoms move through the stored + space-group operation and surface-cell transform factors. + + :param str site_id: + Site identifier returned by :meth:`wyckoff_sites`. + :param str axis: + Parent representative coordinate, one of ``"x"``, ``"y"``, or + ``"z"``. + :param tuple limits: + Fit limits for the coordinate change in parent fractional units. + :param tuple absolute_limits: + Optional absolute parent-coordinate limits, converted to changes + around the stored representative coordinate. + :returns: + Created relative fit parameter. + :rtype: + CTRutil.Parameter + :raises ValueError: + If no matching site-displacement couplings exist. + """ + if self.symmetry_metadata is None: + raise ValueError(f"UnitCell {self.name} has no symmetry metadata.") + if axis not in {"x", "y", "z"}: + raise ValueError("axis must be one of 'x', 'y', or 'z'.") + site = next( + ( + site + for site in self.wyckoff_sites() + if site["site_id"] == site_id + ), + None, + ) + if site is None: + raise ValueError(f"Unknown Wyckoff site {site_id}.") + representative = site.get("representative_parent_fractional") + if representative is None: + raise ValueError( + f"Wyckoff site {site_id} has no representative parent coordinate." + ) + if absolute_limits is not None: + if limits != (-np.inf, np.inf): + raise ValueError("Use either limits or absolute_limits, not both.") + axis_index = {"x": 0, "y": 1, "z": 2}[axis] + axis_value = representative[axis_index] + limits = ( + absolute_limits[0] - axis_value, + absolute_limits[1] - axis_value, + ) + + couplings = [ + coupling + for coupling in self.wyckoff_site_couplings(site_id) + if coupling.axis == axis + ] + if not couplings: + raise ValueError( + f"No site-displacement couplings found for Wyckoff site " + f"{site_id} and parent axis {axis}." + ) + atoms = np.asarray( + [coupling.atom_index for coupling in couplings], + dtype=np.intp, + ) + coordinates = tuple(coupling.coordinate for coupling in couplings) + factors = np.asarray( + [coupling.factor for coupling in couplings], + dtype=np.float64, + ) + + keyargs.setdefault( + "name", + f"{self.name} {site_id}_{axis}_wyckoff_site", + ) + settings = keyargs.setdefault("wyckoff", {}) + settings.update( + { + "site_id": site_id, + "kind": "site_displacement", + "axis": axis, + "value_kind": "delta", + } + ) + return self.addRelParameter( + (atoms, coordinates), + factors, + limits=limits, + **keyargs, + ) + + def addWyckoffShifts( + self, + site_id, + axes=("x", "y", "z"), + limits=(-np.inf, np.inf), + absolute_limits=None, + **keyargs, + ): + """Add representative site-displacement shifts for several axes. + + :param str site_id: + Site identifier returned by :meth:`wyckoff_sites`. + :param axes: + Parent representative coordinate axes to fit. + :type axes: + iterable + :param limits: + Either one ``(lower, upper)`` tuple applied to every axis or a + dictionary mapping axis names to delta limits. + :param absolute_limits: + Optional dictionary mapping axis names to absolute limits. + :returns: + Created relative fit parameters in axis order. + :rtype: + list + """ + + def limit_for(limit_spec, axis_name): + if isinstance(limit_spec, dict): + return limit_spec.get(axis_name, (-np.inf, np.inf)) + return limit_spec + + parameters = [] + for axis in axes: + parameter_keyargs = copy.deepcopy(keyargs) + parameter_limits = limit_for(limits, axis) + parameter_absolute_limits = ( + None + if absolute_limits is None + else limit_for(absolute_limits, axis) + ) + parameters.append( + self.addWyckoffShift( + site_id, + axis, + limits=parameter_limits, + absolute_limits=parameter_absolute_limits, + **parameter_keyargs, + ) + ) + return parameters + def showFitparameters(self): print(self.fitparameterList()) diff --git a/orgui/datautils/xrayutils/test/test_CTRcalc.py b/orgui/datautils/xrayutils/test/test_CTRcalc.py index 260abb7..46f690e 100644 --- a/orgui/datautils/xrayutils/test/test_CTRcalc.py +++ b/orgui/datautils/xrayutils/test/test_CTRcalc.py @@ -665,6 +665,7 @@ def make_single_wyckoff_cell(): element="C", wyckoff_label="1a", coordinates=(), + representative_parent_fractional=(0.25, 0.0, 0.0), variables={"u": 0.25}, occ=1.0, iDW=0.1, @@ -690,6 +691,15 @@ def make_single_wyckoff_cell(): site_id="C_1", ), ), + site_couplings=( + CTRsymmetry.WyckoffSiteCoupling( + atom_index=0, + coordinate="x", + axis="x", + factor=1.0, + site_id="C_1", + ), + ), ), ], ) @@ -750,6 +760,12 @@ def test_unitcell_supercell_preserves_shared_wyckoff_site(self): np.testing.assert_allclose(supercell.basis[:, 1], [0.175, 0.675]) + supercell = unitcell.supercell((2, 1, 1), symmetry="preserve") + supercell.addWyckoffShift("C_1", "x", limits=(-0.2, 0.2)) + supercell.setFitParameters([0.1]) + + np.testing.assert_allclose(supercell.basis[:, 1], [0.175, 0.675]) + def test_unitcell_supercell_independent_wyckoff_sites_fit_separately(self): unitcell = self.make_single_wyckoff_cell() @@ -759,11 +775,183 @@ def test_unitcell_supercell_independent_wyckoff_sites_fit_separately(self): [site["site_id"] for site in supercell.wyckoff_sites()], ["C_1_copy0", "C_1_copy1"], ) - supercell.addWyckoffParameter("C_1_copy0", "u", limits=(-0.2, 0.2)) + supercell.addWyckoffParameter( + "C_1_copy0", + "u", + limits=(-0.2, 0.2), + ) + supercell.setFitParameters([0.1]) + + np.testing.assert_allclose(supercell.basis[:, 1], [0.175, 0.625]) + + supercell = unitcell.supercell((2, 1, 1), symmetry="independent") + supercell.addWyckoffShift("C_1_copy0", "x", limits=(-0.2, 0.2)) supercell.setFitParameters([0.1]) np.testing.assert_allclose(supercell.basis[:, 1], [0.175, 0.625]) + def test_film_forwards_wyckoff_parameters_to_template_unitcell(self): + film = CTRfilm.Film(self.make_single_wyckoff_cell()) + + film.addWyckoffParameter("C_1", "u", limits=(-0.2, 0.2)) + film.setFitParameters([0.1]) + + self.assertEqual(film.unitcell.wyckoff_sites()[0]["status"], "symmetry_preserving") + np.testing.assert_allclose(film.unitcell.basis[:, 1], [0.35]) + + def test_poisson_surface_forwards_wyckoff_shifts_to_template_unitcell(self): + surface = CTRfilm.PoissonSurface(self.make_single_wyckoff_cell()) + + surface.addWyckoffShift("C_1", "x", limits=(-0.2, 0.2)) + surface.setFitParameters([0.1]) + + self.assertEqual(surface.unitcell.wyckoff_sites()[0]["status"], "site_displaced") + np.testing.assert_allclose(surface.unitcell.basis[:, 1], [0.35]) + + def test_epitaxy_interface_uses_unitcell_selector_for_wyckoff_parameters(self): + interface = CTRfilm.EpitaxyInterface( + self.make_single_wyckoff_cell(), + self.make_single_wyckoff_cell(), + ) + + with self.assertRaisesRegex(ValueError, "Missing unit cell name"): + interface.addWyckoffParameter("C_1", "u", limits=(-0.2, 0.2)) + + interface.addWyckoffShift( + "C_1", + "x", + limits=(-0.2, 0.2), + unitcell="top", + ) + interface.setFitParameters([0.1]) + + self.assertEqual(interface.uc_top.wyckoff_sites()[0]["status"], "site_displaced") + self.assertEqual(interface.uc_bottom.wyckoff_sites()[0]["status"], "metadata_only") + np.testing.assert_allclose(interface.uc_top.basis[:, 1], [0.35]) + np.testing.assert_allclose(interface.uc_bottom.basis[:, 1], [0.25]) + + def test_epitaxy_interface_accepts_unitcell_list_for_wyckoff_parameters(self): + interface = CTRfilm.EpitaxyInterface( + self.make_single_wyckoff_cell(), + self.make_single_wyckoff_cell(), + ) + + parameters = interface.addWyckoffParameter( + "C_1", + "u", + limits=(-0.2, 0.2), + unitcell=["top", "bottom"], + ) + interface.setFitParameters([0.1, 0.1]) + + self.assertEqual(len(parameters), 2) + self.assertEqual(interface.uc_top.wyckoff_sites()[0]["status"], "symmetry_preserving") + self.assertEqual(interface.uc_bottom.wyckoff_sites()[0]["status"], "symmetry_preserving") + np.testing.assert_allclose(interface.uc_top.basis[:, 1], [0.35]) + np.testing.assert_allclose(interface.uc_bottom.basis[:, 1], [0.35]) + + def test_sxrdcrystal_links_wyckoff_parameters_across_unitcells(self): + bulk = self.make_single_wyckoff_cell() + film1 = self.make_single_wyckoff_cell() + film2 = self.make_single_wyckoff_cell() + bulk.name = "bulk_template" + film1.name = "film1" + film2.name = "film2" + crystal = CTRcalc.SXRDCrystal(bulk, film1, film2) + + parameter = crystal.addWyckoffParameter( + { + "film1": ("C_1", "u"), + "film2": ("C_1", "u"), + }, + name="shared_u", + limits=(-0.2, 0.2), + ) + crystal.setParameters([0.1]) + + self.assertEqual(parameter.name, "shared_u") + self.assertEqual(crystal.fitparnames, ["shared_u"]) + self.assertEqual(film1.wyckoff_sites()[0]["status"], "symmetry_preserving") + self.assertEqual(film2.wyckoff_sites()[0]["status"], "symmetry_preserving") + np.testing.assert_allclose(film1.basis[:, 1], [0.35]) + np.testing.assert_allclose(film2.basis[:, 1], [0.35]) + + def test_sxrdcrystal_links_wyckoff_shifts_across_unitcells(self): + bulk = self.make_single_wyckoff_cell() + film1 = self.make_single_wyckoff_cell() + film2 = self.make_single_wyckoff_cell() + bulk.name = "bulk_template" + film1.name = "film1" + film2.name = "film2" + crystal = CTRcalc.SXRDCrystal(bulk, film1, film2) + + crystal.addWyckoffShift( + { + "film1": ("C_1", "x"), + "film2": ("C_1", "x"), + }, + name="shared_shift", + limits=(-0.2, 0.2), + ) + crystal.setParameters([0.1]) + + self.assertEqual(crystal.fitparnames, ["shared_shift"]) + self.assertEqual(film1.wyckoff_sites()[0]["status"], "site_displaced") + self.assertEqual(film2.wyckoff_sites()[0]["status"], "site_displaced") + np.testing.assert_allclose(film1.basis[:, 1], [0.35]) + np.testing.assert_allclose(film2.basis[:, 1], [0.35]) + + def test_sxrdcrystal_links_wyckoff_parameters_inside_interface(self): + bulk = self.make_single_wyckoff_cell() + top = self.make_single_wyckoff_cell() + bottom = self.make_single_wyckoff_cell() + bulk.name = "bulk_template" + interface = CTRfilm.EpitaxyInterface(top, bottom, name="interface") + crystal = CTRcalc.SXRDCrystal(bulk, interface) + + with self.assertRaisesRegex(ValueError, "EpitaxyInterface"): + crystal.addWyckoffParameter({"interface": ("C_1", "u")}) + + crystal.addWyckoffParameter( + { + "interface": { + "site_id": "C_1", + "variable": "u", + "unitcell": ("top", "bottom"), + } + }, + name="shared_interface_u", + limits=(-0.2, 0.2), + ) + crystal.setParameters([0.1]) + + self.assertEqual(crystal.fitparnames, ["shared_interface_u"]) + self.assertEqual(interface.uc_top.wyckoff_sites()[0]["status"], "symmetry_preserving") + self.assertEqual(interface.uc_bottom.wyckoff_sites()[0]["status"], "symmetry_preserving") + np.testing.assert_allclose(interface.uc_top.basis[:, 1], [0.35]) + np.testing.assert_allclose(interface.uc_bottom.basis[:, 1], [0.35]) + + def test_sxrdcrystal_links_wyckoff_shift_to_interface_unitcell_key(self): + bulk = self.make_single_wyckoff_cell() + top = self.make_single_wyckoff_cell() + bottom = self.make_single_wyckoff_cell() + bulk.name = "bulk_template" + interface = CTRfilm.EpitaxyInterface(top, bottom, name="interface") + crystal = CTRcalc.SXRDCrystal(bulk, interface) + + crystal.addWyckoffShift( + {("interface", "top"): ("C_1", "x")}, + name="top_shift", + limits=(-0.2, 0.2), + ) + crystal.setParameters([0.1]) + + self.assertEqual(crystal.fitparnames, ["top_shift"]) + self.assertEqual(interface.uc_top.wyckoff_sites()[0]["status"], "site_displaced") + self.assertEqual(interface.uc_bottom.wyckoff_sites()[0]["status"], "metadata_only") + np.testing.assert_allclose(interface.uc_top.basis[:, 1], [0.35]) + np.testing.assert_allclose(interface.uc_bottom.basis[:, 1], [0.25]) + def test_film_and_poisson_rotate_cyclic_layer_order(self): film = CTRfilm.Film(self.make_layered_unitcell("film")) film.basis[0] = 2 diff --git a/orgui/datautils/xrayutils/test/test_CTRsymmetry.py b/orgui/datautils/xrayutils/test/test_CTRsymmetry.py index 5be561b..a6f5031 100644 --- a/orgui/datautils/xrayutils/test/test_CTRsymmetry.py +++ b/orgui/datautils/xrayutils/test/test_CTRsymmetry.py @@ -73,6 +73,54 @@ def make_rutile_110_unitcell(self): ) return model.build_unitcell("RuO2") + def make_pyxtal_generated_unitcell( + self, + spacegroup, + wyckoff_label, + variable_values, + lattice, + transform=None, + element="Si", + ): + """Generate a UnitCell from one PyXtal Wyckoff position.""" + from pymatgen.core import Structure + from pyxtal.symmetry import Group + + wyckoff_position = Group(spacegroup).get_wyckoff_position(wyckoff_label) + coordinates = wyckoff_position.get_all_positions(variable_values) + seed = Structure(lattice, [element] * len(coordinates), coordinates) + surface_spec = CTRsymmetry.SurfaceCellSpec( + lattice.abc, + lattice.angles, + np.identity(3) if transform is None else transform, + layer_origins=(0.0,), + translation_range=0, + ) + model = CTRsymmetry.model_from_seed(seed, surface_spec, tol=1e-3) + return model.build_unitcell(element) + + @staticmethod + def coupling_factor_vectors(unitcell, site_id): + grouped = {} + for coupling in unitcell.wyckoff_couplings(site_id): + key = (coupling.atom_index, coupling.coordinate) + grouped.setdefault(key, {})[coupling.variable] = coupling.factor + return { + tuple(round(factors.get(variable, 0.0), 8) for variable in ("u", "v", "w")) + for factors in grouped.values() + } + + @staticmethod + def site_coupling_factor_vectors(unitcell, site_id): + grouped = {} + for coupling in unitcell.wyckoff_site_couplings(site_id): + key = (coupling.atom_index, coupling.coordinate) + grouped.setdefault(key, {})[coupling.axis] = coupling.factor + return { + tuple(round(factors.get(axis, 0.0), 8) for axis in ("x", "y", "z")) + for factors in grouped.values() + } + def test_pyxtal_assigns_rutile_wyckoff_sites(self): sites, number, symbol = CTRsymmetry.sites_from_seed( self.make_rutile_seed(), @@ -151,7 +199,7 @@ def test_wyckoff_query_exposes_oxygen_u_couplings(self): {(0.0, 1.0), (1.0, -1.0), (0.5, 1.0), (0.5, -1.0)}, ) - def test_wyckoff_parameter_preserves_rutile_u_symmetry(self): + def test_wyckoff_coordinate_parameter_preserves_rutile_u_symmetry(self): unitcell = self.make_rutile_110_unitcell() couplings = unitcell.wyckoff_couplings("O_4f") original = unitcell.basis.copy() @@ -162,6 +210,7 @@ def test_wyckoff_parameter_preserves_rutile_u_symmetry(self): absolute_limits=(0.2, 0.4), ) + self.assertEqual(parameter.settings["wyckoff"]["kind"], "coordinate") self.assertEqual(parameter.settings["wyckoff"]["value_kind"], "delta") np.testing.assert_allclose(unitcell.getInitialParameters(), [0.0]) np.testing.assert_allclose( @@ -180,6 +229,177 @@ def test_wyckoff_parameter_preserves_rutile_u_symmetry(self): expected, ) + def test_wyckoff_site_parameter_displaces_fixed_rutile_site(self): + unitcell = self.make_rutile_110_unitcell() + couplings = [ + coupling + for coupling in unitcell.wyckoff_site_couplings("Ru_2a") + if coupling.axis == "x" + ] + original = unitcell.basis.copy() + + parameter = unitcell.addWyckoffShift( + "Ru_2a", + "x", + absolute_limits=(-0.1, 0.1), + ) + + self.assertEqual(parameter.settings["wyckoff"]["kind"], "site_displacement") + self.assertEqual(unitcell.wyckoff_sites()[0]["status"], "site_displaced") + np.testing.assert_allclose( + unitcell.getStartParamAndLimits()[1:], + ([-0.1], [0.1]), + ) + + unitcell.setFitParameters([0.02]) + + for coupling in couplings: + column = unitcell.parameterLookup[coupling.coordinate] + expected = original[coupling.atom_index, column] + 0.02 * coupling.factor + self.assertAlmostEqual( + unitcell.basis[coupling.atom_index, column], + expected, + ) + + def test_wyckoff_site_parameter_lowers_oxygen_site_symmetry(self): + unitcell = self.make_rutile_110_unitcell() + + unitcell.addWyckoffParameter("O_4f", "u") + self.assertEqual(unitcell.wyckoff_sites()[1]["status"], "symmetry_preserving") + + unitcell = self.make_rutile_110_unitcell() + unitcell.addWyckoffShift("O_4f", "x") + + self.assertEqual(unitcell.wyckoff_sites()[1]["status"], "site_displaced") + + def test_wyckoff_coordinate_parameters_fit_all_site_variables(self): + from pymatgen.core import Lattice + + unitcell = self.make_pyxtal_generated_unitcell( + 62, + "8d", + [0.12, 0.23, 0.34], + Lattice.orthorhombic(4.0, 5.0, 6.0), + element="Mg", + ) + original = unitcell.basis.copy() + + parameters = unitcell.addWyckoffParameters("Mg_8d") + + self.assertEqual( + [par.settings["wyckoff"]["variable"] for par in parameters], + ["u", "v", "w"], + ) + self.assertEqual(unitcell.wyckoff_sites()[0]["status"], "symmetry_preserving") + + unitcell.setFitParameters([0.01, -0.02, 0.03]) + + deltas = {"u": 0.01, "v": -0.02, "w": 0.03} + for coordinate_name in ("x", "y", "z"): + column = unitcell.parameterLookup[coordinate_name] + expected = original[:, column].copy() + for coupling in unitcell.wyckoff_couplings("Mg_8d"): + if coupling.coordinate == coordinate_name: + expected[coupling.atom_index] += ( + coupling.factor * deltas[coupling.variable] + ) + np.testing.assert_allclose(unitcell.basis[:, column], expected) + + def test_fixed_wyckoff_site_rejects_coordinate_parameter(self): + unitcell = self.make_rutile_110_unitcell() + + with self.assertRaisesRegex(ValueError, "no positional coordinate variables"): + unitcell.addWyckoffParameter("Ru_2a", "u") + + def test_trigonal_general_site_exposes_u_minus_v_couplings(self): + from pymatgen.core import Lattice + + unitcell = self.make_pyxtal_generated_unitcell( + 152, + "6c", + [0.12, 0.27, 0.34], + Lattice.hexagonal(4.0, 6.0), + ) + + self.assertEqual( + unitcell.wyckoff_sites()[0]["variables"], + {"u": 0.12, "v": 0.27, "w": 0.34}, + ) + vectors = self.coupling_factor_vectors(unitcell, "Si_6c") + + self.assertIn((1.0, -1.0, 0.0), vectors) + self.assertTrue( + any(abs(vector[0]) == 1.0 and abs(vector[1]) == 1.0 for vector in vectors) + ) + + def test_surface_transform_can_mix_u_v_w_into_one_coordinate(self): + from pymatgen.core import Lattice + + parent_to_surface = np.asarray( + [ + [1.0, 1.0, 1.0], + [0.0, 1.0, 0.0], + [0.0, 0.0, 1.0], + ], + dtype=np.float64, + ) + transform = np.linalg.inv(parent_to_surface) + unitcell = self.make_pyxtal_generated_unitcell( + 62, + "8d", + [0.12, 0.23, 0.34], + Lattice.orthorhombic(4.0, 5.0, 6.0), + transform=transform, + element="Mg", + ) + + vectors = self.coupling_factor_vectors(unitcell, "Mg_8d") + + self.assertIn((1.0, 1.0, 1.0), vectors) + + def test_general_site_displacement_exposes_operation_factor_vectors(self): + from pymatgen.core import Lattice + + unitcell = self.make_pyxtal_generated_unitcell( + 62, + "8d", + [0.12, 0.23, 0.34], + Lattice.orthorhombic(4.0, 5.0, 6.0), + element="Mg", + ) + + vectors = self.site_coupling_factor_vectors(unitcell, "Mg_8d") + + self.assertIn((1.0, 0.0, 0.0), vectors) + self.assertIn((-1.0, 0.0, 0.0), vectors) + self.assertIn((0.0, 0.0, 1.0), vectors) + self.assertIn((0.0, 0.0, -1.0), vectors) + + def test_surface_transform_can_mix_site_displacement_axes(self): + from pymatgen.core import Lattice + + parent_to_surface = np.asarray( + [ + [1.0, 1.0, 1.0], + [0.0, 1.0, 0.0], + [0.0, 0.0, 1.0], + ], + dtype=np.float64, + ) + transform = np.linalg.inv(parent_to_surface) + unitcell = self.make_pyxtal_generated_unitcell( + 62, + "8d", + [0.12, 0.23, 0.34], + Lattice.orthorhombic(4.0, 5.0, 6.0), + transform=transform, + element="Mg", + ) + + vectors = self.site_coupling_factor_vectors(unitcell, "Mg_8d") + + self.assertIn((1.0, 1.0, 1.0), vectors) + def test_manual_atom_parameter_marks_wyckoff_site_as_partially_overridden(self): unitcell = self.make_rutile_110_unitcell() atom_index = unitcell.wyckoff_couplings("O_4f")[0].atom_index @@ -200,9 +420,14 @@ def test_symmetry_metadata_round_trips_without_pyxtal_construction(self): self.assertEqual(restored.wyckoff_sites()[1]["variables"], {"u": 0.30569}) self.assertEqual(len(restored.wyckoff_couplings("O_4f")), 8) + self.assertGreater(len(restored.wyckoff_site_couplings("Ru_2a")), 0) restored.addWyckoffParameter("O_4f", "u", absolute_limits=(0.2, 0.4)) self.assertEqual(restored.wyckoff_sites()[1]["status"], "symmetry_preserving") + restored = UnitCell.fromStr(text) + restored.addWyckoffShift("Ru_2a", "x", absolute_limits=(-0.1, 0.1)) + self.assertEqual(restored.wyckoff_sites()[0]["status"], "site_displaced") + class TestOptionalSymmetryImports(unittest.TestCase): def test_missing_pyxtal_reports_optional_dependency(self): From ea76c7f7e6520314fb03de18719d5db47b7cb041 Mon Sep 17 00:00:00 2001 From: "T. Fuchs" Date: Wed, 8 Jul 2026 23:20:44 -0400 Subject: [PATCH 05/10] fix: disable default native build option --- meson_options.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/meson_options.txt b/meson_options.txt index b59468b..7538b5f 100644 --- a/meson_options.txt +++ b/meson_options.txt @@ -1,6 +1,6 @@ option( 'native_optimization', type: 'boolean', - value: true, + value: false, description: 'Build C++ extensions with native CPU flags for local benchmarking' ) From 3446de66f1cc6954778203d79e89f1fa6c10b7d4 Mon Sep 17 00:00:00 2001 From: "T. Fuchs" Date: Wed, 8 Jul 2026 23:45:16 -0400 Subject: [PATCH 06/10] fix: make Wyckoff ids unique, Apply symmetry also to abs parameters, coherent domains --- orgui/datautils/xrayutils/CTRsymmetry.py | 24 ++++++++- orgui/datautils/xrayutils/CTRuc.py | 49 ++++++++++++++++++- .../datautils/xrayutils/test/test_CTRcalc.py | 31 ++++++++++++ .../xrayutils/test/test_CTRsymmetry.py | 31 ++++++++++++ 4 files changed, 132 insertions(+), 3 deletions(-) diff --git a/orgui/datautils/xrayutils/CTRsymmetry.py b/orgui/datautils/xrayutils/CTRsymmetry.py index 274612a..34cf89c 100644 --- a/orgui/datautils/xrayutils/CTRsymmetry.py +++ b/orgui/datautils/xrayutils/CTRsymmetry.py @@ -31,7 +31,7 @@ from __future__ import annotations -from dataclasses import dataclass, field +from dataclasses import dataclass, field, replace import itertools import math @@ -635,7 +635,7 @@ def sites_from_seed( pyxtal = _import_pyxtal_class() crystal = pyxtal() crystal.from_seed(seed, tol=tol, a_tol=a_tol, backend=backend, style=style) - sites = tuple( + sites = _with_unique_site_ids( _site_from_pyxtal_site( atom_site, variable_names=variable_names, @@ -1211,6 +1211,26 @@ def _site_from_pyxtal_site(atom_site, variable_names, iDW, oDW, occ): ) +def _with_unique_site_ids(sites): + sites = tuple(sites) + site_id_counts = {} + for site in sites: + site_id_counts[site.site_id] = site_id_counts.get(site.site_id, 0) + 1 + if all(count == 1 for count in site_id_counts.values()): + return sites + + occurrences = {} + unique_sites = [] + for site in sites: + if site_id_counts[site.site_id] == 1: + unique_sites.append(site) + continue + occurrence = occurrences.get(site.site_id, 0) + 1 + occurrences[site.site_id] = occurrence + unique_sites.append(replace(site, site_id=f"{site.site_id}_{occurrence}")) + return tuple(unique_sites) + + def _site_operation_matrices(wp, representative, coordinates, variables): group = _import_pyxtal_group()(wp.number) general_operations = group[0].ops diff --git a/orgui/datautils/xrayutils/CTRuc.py b/orgui/datautils/xrayutils/CTRuc.py index 336f5bd..6f29255 100644 --- a/orgui/datautils/xrayutils/CTRuc.py +++ b/orgui/datautils/xrayutils/CTRuc.py @@ -1178,6 +1178,28 @@ def _remap_parameter_atom_indices(parameter, old_to_new): atoms = np.asarray(atoms, dtype=np.intp) parameter.indices = (old_to_new[atoms], parindexes) + @staticmethod + def _update_absolute_parameter_values_after_basis_transform( + transformed, + original_basis, + old_to_new, + ): + new_to_old = np.empty_like(old_to_new) + new_to_old[old_to_new] = np.arange(old_to_new.size, dtype=np.intp) + for parameter in transformed.parameters["absolute"]: + if parameter.value is None or not isinstance(parameter.indices, tuple): + continue + atoms, parindexes = parameter.indices + atoms = np.asarray(atoms, dtype=np.intp) + old_atoms = new_to_old[atoms] + old_values = original_basis[(old_atoms, parindexes)] + new_values = transformed.basis[(atoms, parindexes)] + shifts = np.asarray(new_values - old_values, dtype=np.float64) + if np.allclose(shifts, shifts.flat[0]): + parameter.value += float(shifts.flat[0]) + else: + parameter.value = None + @staticmethod def _transform_symmetry_metadata( metadata, @@ -1321,6 +1343,11 @@ def translate_layered(self, translation, name=None): ) transformed.layerpos[layer_key] += z_translation old_to_new = np.arange(self.basis.shape[0], dtype=np.intp) + self._update_absolute_parameter_values_after_basis_transform( + transformed, + self.basis, + old_to_new, + ) layer_map = {layer: layer for layer in cycle} z_shift_by_layer = {layer: z_translation for layer in cycle} transformed.symmetry_metadata = self._transform_symmetry_metadata( @@ -1383,6 +1410,11 @@ def transform_basis_values(values): for parameters in transformed.parameters.values(): for parameter in parameters: self._remap_parameter_atom_indices(parameter, old_to_new) + self._update_absolute_parameter_values_after_basis_transform( + transformed, + self.basis, + old_to_new, + ) transformed.layerpos = copy.deepcopy(self.layerpos) for layer in cycle: @@ -1457,6 +1489,11 @@ def affine_layer_transform(self, translation, name=None): if transformed._basis_parvalues is not None: transformed._basis_parvalues[:, 1:3] += xy_translation old_to_new = np.arange(self.basis.shape[0], dtype=np.intp) + self._update_absolute_parameter_values_after_basis_transform( + transformed, + self.basis, + old_to_new, + ) layer_map = {layer: layer for layer in cycle} z_shift_by_layer = {layer: 0.0 for layer in cycle} transformed.symmetry_metadata = self._transform_symmetry_metadata( @@ -1521,6 +1558,11 @@ def transform_basis_values(values): for parameters in transformed.parameters.values(): for parameter in parameters: self._remap_parameter_atom_indices(parameter, old_to_new) + self._update_absolute_parameter_values_after_basis_transform( + transformed, + self.basis, + old_to_new, + ) transformed.layerpos = copy.deepcopy(self.layerpos) for layer in cycle: @@ -1576,7 +1618,12 @@ def supercell(self, repeats, symmetry="preserve", name=None): transformed.layer_transition = copy.deepcopy(self.layer_transition) transformed.refHKLTransform = np.copy(self.refHKLTransform) transformed.refRealTransform = np.copy(self.refRealTransform) - transformed.coherentDomainMatrix = copy.deepcopy(self.coherentDomainMatrix) + transformed.coherentDomainMatrix = [] + repeats_float = repeats.astype(np.float64) + for matrix in self.coherentDomainMatrix: + matrix_new = np.array(matrix, copy=True) + matrix_new[:, -1] = matrix_new[:, -1] / repeats_float + transformed.coherentDomainMatrix.append(matrix_new) transformed.coherentDomainOccupancy = copy.deepcopy( self.coherentDomainOccupancy ) diff --git a/orgui/datautils/xrayutils/test/test_CTRcalc.py b/orgui/datautils/xrayutils/test/test_CTRcalc.py index 46f690e..13dee1b 100644 --- a/orgui/datautils/xrayutils/test/test_CTRcalc.py +++ b/orgui/datautils/xrayutils/test/test_CTRcalc.py @@ -410,6 +410,17 @@ def test_unitcell_translate_layered_translates_full_two_layer_cycle(self): {0.0: 0.0, 1.0: -1.0, 2.0: -0.5}, ) + def test_unitcell_translate_layered_shifts_absolute_parameter_values(self): + unitcell = self.make_layered_unitcell() + unitcell.addFitParameter((0, "z"), limits=(-0.1, 2.0), name="z_abs") + unitcell.setFitParameters([0.05]) + + transformed = unitcell.translate_layered([0, 0, 3]) + + np.testing.assert_allclose(transformed.getInitialParameters(), [1.05]) + transformed.setFitParameters([1.05]) + self.assertAlmostEqual(transformed.basis[0, 3], 1.05) + def test_unitcell_translate_layered_preserves_full_cycle_symmetry_metadata(self): unitcell = CTRcalc.UnitCell( [3.0, 3.0, 6.0], @@ -723,6 +734,26 @@ def test_unitcell_supercell_scales_lattice_and_coordinates(self): ) self.assertEqual(supercell.parameters, {"absolute": [], "relative": []}) + def test_unitcell_supercell_rescales_coherent_domain_translations(self): + unitcell = CTRcalc.UnitCell( + [2.0, 3.0, 4.0], + [90.0, 90.0, 90.0], + name="domain", + ) + unitcell.addAtom("C", [0.25, 0.0, 0.0], 0.1, 0.1, 1.0, layer=1) + unitcell.layerpos[1.0] = 0.0 + domain = np.eye(3, 4) + domain[0, 3] = 0.5 + unitcell.coherentDomainMatrix = [domain] + + supercell = unitcell.supercell((2, 1, 1)) + + np.testing.assert_allclose( + supercell.coherentDomainMatrix[0][:, -1], + [0.25, 0.0, 0.0], + ) + np.testing.assert_allclose(supercell.pos_cart(0), unitcell.pos_cart(0)) + def test_unitcell_supercell_repeats_layers_along_z(self): unitcell = self.make_layered_unitcell() diff --git a/orgui/datautils/xrayutils/test/test_CTRsymmetry.py b/orgui/datautils/xrayutils/test/test_CTRsymmetry.py index a6f5031..ef500f3 100644 --- a/orgui/datautils/xrayutils/test/test_CTRsymmetry.py +++ b/orgui/datautils/xrayutils/test/test_CTRsymmetry.py @@ -37,6 +37,37 @@ HAS_PYMATGEN = importlib.util.find_spec("pymatgen") is not None +class TestSymmetryUtilities(unittest.TestCase): + def test_duplicate_wyckoff_site_ids_are_disambiguated(self): + sites = CTRsymmetry._with_unique_site_ids( + ( + CTRsymmetry.WyckoffSiteSpec( + site_id="O_4f", + element="O", + wyckoff_label="4f", + coordinates=(), + ), + CTRsymmetry.WyckoffSiteSpec( + site_id="O_4f", + element="O", + wyckoff_label="4f", + coordinates=(), + ), + CTRsymmetry.WyckoffSiteSpec( + site_id="Ru_2a", + element="Ru", + wyckoff_label="2a", + coordinates=(), + ), + ) + ) + + self.assertEqual( + [site.site_id for site in sites], + ["O_4f_1", "O_4f_2", "Ru_2a"], + ) + + @unittest.skipUnless( HAS_PYXTAL and HAS_PYMATGEN, "PyXtal symmetry tests require PyXtal", From f8938dd4f129f7690e5835a2a6ad214379aa9d1c Mon Sep 17 00:00:00 2001 From: "T. Fuchs" Date: Thu, 9 Jul 2026 20:37:28 -0400 Subject: [PATCH 07/10] feat: handle fit parameter coupling for translate_layered, affine_layer_transform and supercell calls --- orgui/datautils/xrayutils/CTRuc.py | 19 ++++++++++++ .../datautils/xrayutils/test/test_CTRcalc.py | 30 +++++++++++++++++++ 2 files changed, 49 insertions(+) diff --git a/orgui/datautils/xrayutils/CTRuc.py b/orgui/datautils/xrayutils/CTRuc.py index 6f29255..819341f 100644 --- a/orgui/datautils/xrayutils/CTRuc.py +++ b/orgui/datautils/xrayutils/CTRuc.py @@ -1213,6 +1213,7 @@ def _transform_symmetry_metadata( metadata = copy.deepcopy(metadata) coordinate_index = {"x": 0, "y": 1, "z": 2} atoms = [] + site_parent_shifts = {} for atom in metadata.atoms: z_shift = z_shift_by_layer[atom.layer] coordinate_shift = np.asarray( @@ -1240,6 +1241,7 @@ def _transform_symmetry_metadata( ) ) parent_shift = metadata.surface_spec.transform @ coordinate_shift + site_parent_shifts.setdefault(atom.site_id, parent_shift) atoms.append( dataclasses.replace( atom, @@ -1257,6 +1259,21 @@ def _transform_symmetry_metadata( site_couplings=tuple(site_couplings), ) ) + sites = [] + for site in metadata.sites: + representative = site.representative_parent_fractional + if representative is not None and site.site_id in site_parent_shifts: + representative = tuple( + np.asarray(representative, dtype=np.float64) + + site_parent_shifts[site.site_id] + ) + sites.append( + dataclasses.replace( + site, + representative_parent_fractional=representative, + ) + ) + metadata.sites = tuple(sites) metadata.atoms = sorted(atoms, key=lambda atom: atom.atom_index) return metadata @@ -1729,6 +1746,8 @@ def supercell(self, repeats, symmetry="preserve", name=None): for layer in cycle ) transformed.parameters = {"absolute": [], "relative": []} + transformed.basis_0 = np.copy(transformed.basis) + transformed._basis_parvalues = np.copy(transformed.basis) transformed.symmetry_metadata = self._supercell_symmetry_metadata( rows, repeats.astype(np.float64), diff --git a/orgui/datautils/xrayutils/test/test_CTRcalc.py b/orgui/datautils/xrayutils/test/test_CTRcalc.py index 13dee1b..5e201e4 100644 --- a/orgui/datautils/xrayutils/test/test_CTRcalc.py +++ b/orgui/datautils/xrayutils/test/test_CTRcalc.py @@ -553,6 +553,23 @@ def test_unitcell_translate_layered_remaps_symmetry_metadata_atoms(self): self.assertEqual(atom.couplings[0].atom_index, 0) self.assertAlmostEqual(atom.couplings[0].constant, 1.0) + def test_unitcell_translate_layered_updates_wyckoff_representatives(self): + unitcell = self.make_single_wyckoff_cell() + + transformed = unitcell.translate_layered([1, -1, 0]) + + site = transformed.wyckoff_sites()[0] + np.testing.assert_allclose( + site["representative_parent_fractional"], + [1.25, -1.0, 0.0], + ) + parameter = transformed.addWyckoffShift( + "C_1", + "x", + absolute_limits=(1.0, 1.5), + ) + self.assertEqual(parameter.limits, (-0.25, 0.25)) + def test_unitcell_translate_layered_rejects_unsupported_affines(self): unitcell = self.make_layered_unitcell() @@ -754,6 +771,19 @@ def test_unitcell_supercell_rescales_coherent_domain_translations(self): ) np.testing.assert_allclose(supercell.pos_cart(0), unitcell.pos_cart(0)) + def test_unitcell_supercell_uses_current_coordinates_as_new_reference(self): + unitcell = self.make_single_wyckoff_cell() + unitcell.addRelParameter(([0], "x"), [1.0], (-1.0, 1.0), name="dx") + unitcell.setFitParameters([0.05]) + + supercell = unitcell.supercell((2, 1, 1)) + + np.testing.assert_allclose(supercell.basis[:, 1], [0.15, 0.65]) + np.testing.assert_allclose(supercell.basis_0, supercell.basis) + supercell.addRelParameter(([0], "x"), [1.0], (-1.0, 1.0), name="dx2") + supercell.setFitParameters([0.0]) + np.testing.assert_allclose(supercell.basis[:, 1], [0.15, 0.65]) + def test_unitcell_supercell_repeats_layers_along_z(self): unitcell = self.make_layered_unitcell() From adc09f3725ba011577293317f90a552310dbb0f7 Mon Sep 17 00:00:00 2001 From: "T. Fuchs" Date: Fri, 10 Jul 2026 00:29:55 -0400 Subject: [PATCH 08/10] fix(phys): preserve symmetry parent lattice metadata Persist parent lattice parameters in Wyckoff symmetry metadata, round-trip them through UnitCell parsing, and consolidate the new symmetry/layer helper paths without changing fitted behavior. --- orgui/datautils/xrayutils/CTRfilm.py | 107 ++-- orgui/datautils/xrayutils/CTRsymmetry.py | 26 +- orgui/datautils/xrayutils/CTRuc.py | 581 +++++++++--------- .../xrayutils/test/test_CTRsymmetry.py | 17 + 4 files changed, 386 insertions(+), 345 deletions(-) diff --git a/orgui/datautils/xrayutils/CTRfilm.py b/orgui/datautils/xrayutils/CTRfilm.py index 4151087..4806cd9 100644 --- a/orgui/datautils/xrayutils/CTRfilm.py +++ b/orgui/datautils/xrayutils/CTRfilm.py @@ -125,26 +125,33 @@ def layer_behaviour(self): def layer_behaviour(self, behavior): self.layer_behavior = behavior - def _wyckoff_target_unitcells(self, kwarg): + def _wyckoff_target_unitcells(self, kwargs): if hasattr(self, "uc_top") and hasattr(self, "uc_bottom"): - if "unitcell" not in kwarg: + if "unitcell" not in kwargs: raise ValueError( "Missing unit cell name. Provide unit cell name as kwarg " "'unitcell'" ) - unitcell = kwarg["unitcell"] - if isinstance(unitcell, list): + unitcell = kwargs["unitcell"] + if isinstance(unitcell, list | tuple): return [self[ucn] for ucn in unitcell] return [self[unitcell]] return [self.unitcell] + def _forward_wyckoff_call(self, method_name, *args, **kwargs): + parameters = [ + getattr(unitcell, method_name)(*args, **kwargs) + for unitcell in self._wyckoff_target_unitcells(kwargs) + ] + return parameters if len(parameters) > 1 else parameters[0] + def addWyckoffParameter( self, site_id, variable, limits=(-np.inf, np.inf), absolute_limits=None, - **kwarg, + **kwargs, ): """Add a Wyckoff variable parameter on contained unit cells. @@ -156,23 +163,20 @@ def addWyckoffParameter( Delta fit limits in fractional units. :param tuple absolute_limits: Optional absolute variable limits. - :param kwarg: + :param kwargs: For ``EpitaxyInterface``, provide ``unitcell="top"``, ``unitcell="bottom"``, or a list of these names. :returns: Created parameter or list of parameters. """ - parameters = [ - unitcell.addWyckoffParameter( - site_id, - variable, - limits=limits, - absolute_limits=absolute_limits, - **kwarg, - ) - for unitcell in self._wyckoff_target_unitcells(kwarg) - ] - return parameters if len(parameters) > 1 else parameters[0] + return self._forward_wyckoff_call( + "addWyckoffParameter", + site_id, + variable, + limits=limits, + absolute_limits=absolute_limits, + **kwargs, + ) def addWyckoffParameters( self, @@ -180,7 +184,7 @@ def addWyckoffParameters( variables=None, limits=(-np.inf, np.inf), absolute_limits=None, - **kwarg, + **kwargs, ): """Add several Wyckoff variable parameters on contained unit cells. @@ -192,24 +196,21 @@ def addWyckoffParameters( Delta fit limits in fractional units. :param absolute_limits: Optional absolute variable limits. - :param kwarg: + :param kwargs: For ``EpitaxyInterface``, provide ``unitcell="top"``, ``unitcell="bottom"``, or a list of these names. :returns: Created parameters. A list of lists is returned when multiple contained unit cells are selected. """ - parameters = [ - unitcell.addWyckoffParameters( - site_id, - variables=variables, - limits=limits, - absolute_limits=absolute_limits, - **kwarg, - ) - for unitcell in self._wyckoff_target_unitcells(kwarg) - ] - return parameters if len(parameters) > 1 else parameters[0] + return self._forward_wyckoff_call( + "addWyckoffParameters", + site_id, + variables=variables, + limits=limits, + absolute_limits=absolute_limits, + **kwargs, + ) def addWyckoffShift( self, @@ -217,7 +218,7 @@ def addWyckoffShift( axis, limits=(-np.inf, np.inf), absolute_limits=None, - **kwarg, + **kwargs, ): """Add a representative Wyckoff-site shift on contained unit cells. @@ -230,23 +231,20 @@ def addWyckoffShift( Delta fit limits in parent fractional units. :param tuple absolute_limits: Optional absolute parent-coordinate limits. - :param kwarg: + :param kwargs: For ``EpitaxyInterface``, provide ``unitcell="top"``, ``unitcell="bottom"``, or a list of these names. :returns: Created parameter or list of parameters. """ - parameters = [ - unitcell.addWyckoffShift( - site_id, - axis, - limits=limits, - absolute_limits=absolute_limits, - **kwarg, - ) - for unitcell in self._wyckoff_target_unitcells(kwarg) - ] - return parameters if len(parameters) > 1 else parameters[0] + return self._forward_wyckoff_call( + "addWyckoffShift", + site_id, + axis, + limits=limits, + absolute_limits=absolute_limits, + **kwargs, + ) def addWyckoffShifts( self, @@ -254,7 +252,7 @@ def addWyckoffShifts( axes=("x", "y", "z"), limits=(-np.inf, np.inf), absolute_limits=None, - **kwarg, + **kwargs, ): """Add representative Wyckoff-site shifts on contained unit cells. @@ -266,24 +264,21 @@ def addWyckoffShifts( Delta fit limits in parent fractional units. :param absolute_limits: Optional absolute parent-coordinate limits. - :param kwarg: + :param kwargs: For ``EpitaxyInterface``, provide ``unitcell="top"``, ``unitcell="bottom"``, or a list of these names. :returns: Created parameters. A list of lists is returned when multiple contained unit cells are selected. """ - parameters = [ - unitcell.addWyckoffShifts( - site_id, - axes=axes, - limits=limits, - absolute_limits=absolute_limits, - **kwarg, - ) - for unitcell in self._wyckoff_target_unitcells(kwarg) - ] - return parameters if len(parameters) > 1 else parameters[0] + return self._forward_wyckoff_call( + "addWyckoffShifts", + site_id, + axes=axes, + limits=limits, + absolute_limits=absolute_limits, + **kwargs, + ) @property def start_layer_number(self): diff --git a/orgui/datautils/xrayutils/CTRsymmetry.py b/orgui/datautils/xrayutils/CTRsymmetry.py index 34cf89c..472c062 100644 --- a/orgui/datautils/xrayutils/CTRsymmetry.py +++ b/orgui/datautils/xrayutils/CTRsymmetry.py @@ -923,6 +923,8 @@ def symmetry_metadata_to_lines(model): if model.spacegroup_number is not None: symbol = model.spacegroup_symbol or "" lines.append(f"spacegroup: {model.spacegroup_number} {symbol}".rstrip()) + lines.append("parent_a: " + _format_values(model.surface_spec.parent_a)) + lines.append("parent_alpha: " + _format_values(model.surface_spec.parent_alpha)) lines.append("surface_transform:") for row in model.surface_spec.transform: lines.append(" " + _format_values(row)) @@ -1041,6 +1043,19 @@ def symmetry_metadata_from_lines(lines, unitcell=None): if len(parts) > 1: spacegroup_symbol = " ".join(parts[1:]) + parent_a = (1.0, 1.0, 1.0) + parent_alpha = (90.0, 90.0, 90.0) + if "parent_a" in sections: + parent_a = tuple(float(value) for value in sections["parent_a"][0].split()) + elif unitcell is not None: + parent_a = tuple(float(value) for value in unitcell.a) + if "parent_alpha" in sections: + parent_alpha = tuple( + float(value) for value in sections["parent_alpha"][0].split() + ) + elif unitcell is not None: + parent_alpha = tuple(float(value) for value in np.rad2deg(unitcell.alpha)) + if "surface_transform" in sections: transform = np.asarray( [ @@ -1119,11 +1134,6 @@ def symmetry_metadata_from_lines(lines, unitcell=None): ) ) - parent_a = (1.0, 1.0, 1.0) - parent_alpha = (90.0, 90.0, 90.0) - if unitcell is not None: - parent_a = tuple(float(value) for value in unitcell.a) - parent_alpha = tuple(float(value) for value in np.rad2deg(unitcell.alpha)) surface_spec = SurfaceCellSpec( parent_a, parent_alpha, @@ -1380,6 +1390,12 @@ def _collect_sections(lines): if line.startswith("spacegroup:"): sections["spacegroup"] = [line.split(":", 1)[1].strip()] current = None + elif line.startswith("parent_a:"): + sections["parent_a"] = [line.split(":", 1)[1].strip()] + current = None + elif line.startswith("parent_alpha:"): + sections["parent_alpha"] = [line.split(":", 1)[1].strip()] + current = None elif line.startswith("surface_origin:"): sections["surface_origin"] = [line.split(":", 1)[1].strip()] current = None diff --git a/orgui/datautils/xrayutils/CTRuc.py b/orgui/datautils/xrayutils/CTRuc.py index 819341f..7aba1f0 100644 --- a/orgui/datautils/xrayutils/CTRuc.py +++ b/orgui/datautils/xrayutils/CTRuc.py @@ -1287,6 +1287,129 @@ def _layer_positions_for_cycle(self, cycle): layer_positions[layer] = float(np.mean(self.basis[where, 3])) return layer_positions + def _validated_layered_translation(self, translation): + translation = self._translation_from_affine_input(translation) + xy_translation = translation[:2] + if not np.allclose(xy_translation, np.rint(xy_translation)): + raise ValueError("x and y affine translations must be integers.") + if not np.isclose(translation[2], np.rint(translation[2])): + raise ValueError("z affine translation must be an integer layer step.") + return np.rint(xy_translation).astype(np.float64), int( + np.rint(translation[2]) + ) + + @staticmethod + def _shift_layered_basis_values( + values, + cycle, + xy_translation, + z_shift_by_layer, + layer_map=None, + ): + if values is None or values.size == 0: + return values + values_new = np.array(values, copy=True) + values_new[:, 1:3] += xy_translation + for layer in cycle: + where = values[:, 7] == layer + if not np.any(where): + continue + values_new[where, 3] += z_shift_by_layer[layer] + if layer_map is not None: + values_new[where, 7] = layer_map[layer] + return values_new + + def _finish_layered_transform( + self, + transformed, + old_to_new, + layer_map, + xy_translation, + z_shift_by_layer, + test_special_formfactors=False, + ): + for parameters in transformed.parameters.values(): + for parameter in parameters: + self._remap_parameter_atom_indices(parameter, old_to_new) + self._update_absolute_parameter_values_after_basis_transform( + transformed, + self.basis, + old_to_new, + ) + transformed.symmetry_metadata = self._transform_symmetry_metadata( + self.symmetry_metadata, + old_to_new, + layer_map, + xy_translation, + z_shift_by_layer, + ) + if test_special_formfactors: + transformed._test_special_formfactors() + return transformed + + def _apply_reordered_layered_transform( + self, + transformed, + cycle, + xy_translation, + layer_map, + z_shift_by_layer, + ): + basis = self._shift_layered_basis_values( + self.basis, + cycle, + xy_translation, + z_shift_by_layer, + layer_map, + ) + basis_0 = self._shift_layered_basis_values( + self.basis_0, + cycle, + xy_translation, + z_shift_by_layer, + layer_map, + ) + basis_parvalues = self._shift_layered_basis_values( + self._basis_parvalues, + cycle, + xy_translation, + z_shift_by_layer, + layer_map, + ) + + order_lookup = {layer: index for index, layer in enumerate(cycle)} + row_order = np.array( + sorted(range(basis.shape[0]), key=lambda i: (order_lookup[basis[i, 7]], i)), + dtype=np.intp, + ) + old_to_new = np.empty_like(row_order) + old_to_new[row_order] = np.arange(row_order.size) + + transformed.basis = basis[row_order] + transformed.basis_0 = basis_0[row_order] + if basis_parvalues is not None: + transformed._basis_parvalues = basis_parvalues[row_order] + if self.errors is not None: + transformed.errors = np.array(self.errors, copy=True)[row_order] + if self._errors_parvalues is not None: + transformed._errors_parvalues = np.array( + self._errors_parvalues, + copy=True, + )[row_order] + transformed.names = [self.names[i] for i in row_order] + transformed.dw_increase_constraint = self.dw_increase_constraint[row_order] + if hasattr(self, "f"): + transformed.f = self.f[row_order] + + return self._finish_layered_transform( + transformed, + old_to_new, + layer_map, + xy_translation, + z_shift_by_layer, + test_special_formfactors=True, + ) + def translate_layered(self, translation, name=None): """Return a translated copy with cyclic z-layer wrapping. @@ -1316,16 +1439,7 @@ def translate_layered(self, translation, name=None): If the affine linear part is not identity or any requested translation is not an integer. """ - translation = self._translation_from_affine_input(translation) - xy_translation = translation[:2] - if not np.allclose(xy_translation, np.rint(xy_translation)): - raise ValueError("x and y affine translations must be integers.") - if not np.isclose(translation[2], np.rint(translation[2])): - raise ValueError("z affine translation must be an integer layer step.") - - xy_translation = np.rint(xy_translation).astype(np.float64) - z_steps = int(np.rint(translation[2])) - + xy_translation, z_steps = self._validated_layered_translation(translation) transformed = copy.deepcopy(self) if name is not None: transformed.name = name @@ -1344,13 +1458,25 @@ def translate_layered(self, translation, name=None): if z_steps == 0: return transformed if z_steps_effective == 0: - transformed.basis[:, 1:3] += xy_translation - transformed.basis[:, 3] += z_translation - transformed.basis_0[:, 1:3] += xy_translation - transformed.basis_0[:, 3] += z_translation - if transformed._basis_parvalues is not None: - transformed._basis_parvalues[:, 1:3] += xy_translation - transformed._basis_parvalues[:, 3] += z_translation + z_shift_by_layer = {layer: z_translation for layer in cycle} + transformed.basis = self._shift_layered_basis_values( + transformed.basis, + cycle, + xy_translation, + z_shift_by_layer, + ) + transformed.basis_0 = self._shift_layered_basis_values( + transformed.basis_0, + cycle, + xy_translation, + z_shift_by_layer, + ) + transformed._basis_parvalues = self._shift_layered_basis_values( + transformed._basis_parvalues, + cycle, + xy_translation, + z_shift_by_layer, + ) for layer in cycle: layer_key = float(layer) if layer_key not in transformed.layerpos: @@ -1360,77 +1486,27 @@ def translate_layered(self, translation, name=None): ) transformed.layerpos[layer_key] += z_translation old_to_new = np.arange(self.basis.shape[0], dtype=np.intp) - self._update_absolute_parameter_values_after_basis_transform( - transformed, - self.basis, - old_to_new, - ) layer_map = {layer: layer for layer in cycle} - z_shift_by_layer = {layer: z_translation for layer in cycle} - transformed.symmetry_metadata = self._transform_symmetry_metadata( - self.symmetry_metadata, + return self._finish_layered_transform( + transformed, old_to_new, layer_map, xy_translation, z_shift_by_layer, ) - return transformed layer_positions = self._layer_positions_for_cycle(cycle) - layer_map = { layer: cycle[(index + z_steps_effective) % layer_count] for index, layer in enumerate(cycle) } - - def transform_basis_values(values): - if values is None or values.size == 0: - return values - values_new = np.array(values, copy=True) - values_new[:, 1:3] += xy_translation - values_new[:, 3] += z_translation - for layer in cycle: - where = values[:, 7] == layer - if not np.any(where): - continue - new_layer = layer_map[layer] - values_new[where, 7] = new_layer - return values_new - - basis = transform_basis_values(self.basis) - basis_0 = transform_basis_values(self.basis_0) - basis_parvalues = transform_basis_values(self._basis_parvalues) - - order_lookup = {layer: index for index, layer in enumerate(cycle)} - row_order = np.array( - sorted(range(basis.shape[0]), key=lambda i: (order_lookup[basis[i, 7]], i)), - dtype=np.intp, - ) - old_to_new = np.empty_like(row_order) - old_to_new[row_order] = np.arange(row_order.size) - - transformed.basis = basis[row_order] - transformed.basis_0 = basis_0[row_order] - if basis_parvalues is not None: - transformed._basis_parvalues = basis_parvalues[row_order] - if self.errors is not None: - transformed.errors = np.array(self.errors, copy=True)[row_order] - if self._errors_parvalues is not None: - transformed._errors_parvalues = np.array( - self._errors_parvalues, copy=True - )[row_order] - transformed.names = [self.names[i] for i in row_order] - transformed.dw_increase_constraint = self.dw_increase_constraint[row_order] - if hasattr(self, "f"): - transformed.f = self.f[row_order] - - for parameters in transformed.parameters.values(): - for parameter in parameters: - self._remap_parameter_atom_indices(parameter, old_to_new) - self._update_absolute_parameter_values_after_basis_transform( + z_shift_by_layer = {layer: z_translation for layer in cycle} + transformed = self._apply_reordered_layered_transform( transformed, - self.basis, - old_to_new, + cycle, + xy_translation, + layer_map, + z_shift_by_layer, ) transformed.layerpos = copy.deepcopy(self.layerpos) @@ -1441,15 +1517,6 @@ def transform_basis_values(values): ) if transformed._explicit_layer_cycle is not None: transformed._explicit_layer_cycle = tuple(cycle) - z_shift_by_layer = {layer: z_translation for layer in cycle} - transformed.symmetry_metadata = self._transform_symmetry_metadata( - self.symmetry_metadata, - old_to_new, - layer_map, - xy_translation, - z_shift_by_layer, - ) - transformed._test_special_formfactors() return transformed def affine_layer_transform(self, translation, name=None): @@ -1478,16 +1545,7 @@ def affine_layer_transform(self, translation, name=None): If the affine linear part is not identity or any requested translation is not an integer. """ - translation = self._translation_from_affine_input(translation) - xy_translation = translation[:2] - if not np.allclose(xy_translation, np.rint(xy_translation)): - raise ValueError("x and y affine translations must be integers.") - if not np.isclose(translation[2], np.rint(translation[2])): - raise ValueError("z affine translation must be an integer layer step.") - - xy_translation = np.rint(xy_translation).astype(np.float64) - z_steps = int(np.rint(translation[2])) - + xy_translation, z_steps = self._validated_layered_translation(translation) transformed = copy.deepcopy(self) if name is not None: transformed.name = name @@ -1501,26 +1559,34 @@ def affine_layer_transform(self, translation, name=None): layer_count = len(cycle) z_steps_effective = z_steps % layer_count if z_steps_effective == 0: - transformed.basis[:, 1:3] += xy_translation - transformed.basis_0[:, 1:3] += xy_translation - if transformed._basis_parvalues is not None: - transformed._basis_parvalues[:, 1:3] += xy_translation - old_to_new = np.arange(self.basis.shape[0], dtype=np.intp) - self._update_absolute_parameter_values_after_basis_transform( - transformed, - self.basis, - old_to_new, + z_shift_by_layer = {layer: 0.0 for layer in cycle} + transformed.basis = self._shift_layered_basis_values( + transformed.basis, + cycle, + xy_translation, + z_shift_by_layer, + ) + transformed.basis_0 = self._shift_layered_basis_values( + transformed.basis_0, + cycle, + xy_translation, + z_shift_by_layer, ) + transformed._basis_parvalues = self._shift_layered_basis_values( + transformed._basis_parvalues, + cycle, + xy_translation, + z_shift_by_layer, + ) + old_to_new = np.arange(self.basis.shape[0], dtype=np.intp) layer_map = {layer: layer for layer in cycle} - z_shift_by_layer = {layer: 0.0 for layer in cycle} - transformed.symmetry_metadata = self._transform_symmetry_metadata( - self.symmetry_metadata, + return self._finish_layered_transform( + transformed, old_to_new, layer_map, xy_translation, z_shift_by_layer, ) - return transformed layer_positions = self._layer_positions_for_cycle(cycle) layer_map = { @@ -1531,54 +1597,12 @@ def affine_layer_transform(self, translation, name=None): layer: layer_positions[layer_map[layer]] - layer_positions[layer] for layer in cycle } - - def transform_basis_values(values): - if values is None or values.size == 0: - return values - values_new = np.array(values, copy=True) - values_new[:, 1:3] += xy_translation - for layer in cycle: - where = values[:, 7] == layer - if not np.any(where): - continue - values_new[where, 3] += z_shift_by_layer[layer] - values_new[where, 7] = layer_map[layer] - return values_new - - basis = transform_basis_values(self.basis) - basis_0 = transform_basis_values(self.basis_0) - basis_parvalues = transform_basis_values(self._basis_parvalues) - - order_lookup = {layer: index for index, layer in enumerate(cycle)} - row_order = np.array( - sorted(range(basis.shape[0]), key=lambda i: (order_lookup[basis[i, 7]], i)), - dtype=np.intp, - ) - old_to_new = np.empty_like(row_order) - old_to_new[row_order] = np.arange(row_order.size) - - transformed.basis = basis[row_order] - transformed.basis_0 = basis_0[row_order] - if basis_parvalues is not None: - transformed._basis_parvalues = basis_parvalues[row_order] - if self.errors is not None: - transformed.errors = np.array(self.errors, copy=True)[row_order] - if self._errors_parvalues is not None: - transformed._errors_parvalues = np.array( - self._errors_parvalues, copy=True - )[row_order] - transformed.names = [self.names[i] for i in row_order] - transformed.dw_increase_constraint = self.dw_increase_constraint[row_order] - if hasattr(self, "f"): - transformed.f = self.f[row_order] - - for parameters in transformed.parameters.values(): - for parameter in parameters: - self._remap_parameter_atom_indices(parameter, old_to_new) - self._update_absolute_parameter_values_after_basis_transform( + transformed = self._apply_reordered_layered_transform( transformed, - self.basis, - old_to_new, + cycle, + xy_translation, + layer_map, + z_shift_by_layer, ) transformed.layerpos = copy.deepcopy(self.layerpos) @@ -1586,14 +1610,6 @@ def transform_basis_values(values): transformed.layerpos[float(layer)] = layer_positions[layer] if transformed._explicit_layer_cycle is not None: transformed._explicit_layer_cycle = tuple(cycle) - transformed.symmetry_metadata = self._transform_symmetry_metadata( - self.symmetry_metadata, - old_to_new, - layer_map, - xy_translation, - z_shift_by_layer, - ) - transformed._test_special_formfactors() return transformed def supercell(self, repeats, symmetry="preserve", name=None): @@ -2181,6 +2197,83 @@ def atom_wyckoff_metadata(self, atom_index): return None return self.symmetry_metadata.atom_wyckoff_metadata(atom_index) + @staticmethod + def _limit_for(limit_spec, parameter_name): + if isinstance(limit_spec, dict): + return limit_spec.get(parameter_name, (-np.inf, np.inf)) + return limit_spec + + @staticmethod + def _delta_limits(limits, absolute_limits, reference_value): + if absolute_limits is None: + return limits + if limits != (-np.inf, np.inf): + raise ValueError("Use either limits or absolute_limits, not both.") + return ( + absolute_limits[0] - reference_value, + absolute_limits[1] - reference_value, + ) + + def _wyckoff_site(self, site_id): + if self.symmetry_metadata is None: + raise ValueError(f"UnitCell {self.name} has no symmetry metadata.") + site = next( + ( + site + for site in self.wyckoff_sites() + if site["site_id"] == site_id + ), + None, + ) + if site is None: + raise ValueError(f"Unknown Wyckoff site {site_id}.") + return site + + @staticmethod + def _coupling_parameter_arrays(couplings): + atoms = np.asarray( + [coupling.atom_index for coupling in couplings], + dtype=np.intp, + ) + coordinates = tuple(coupling.coordinate for coupling in couplings) + factors = np.asarray( + [coupling.factor for coupling in couplings], + dtype=np.float64, + ) + return atoms, coordinates, factors + + def _add_wyckoff_relative_parameter( + self, + site_id, + selector_name, + selector_value, + kind, + couplings, + limits, + default_name, + empty_message, + keyargs, + ): + if not couplings: + raise ValueError(empty_message) + atoms, coordinates, factors = self._coupling_parameter_arrays(couplings) + keyargs.setdefault("name", default_name) + settings = keyargs.setdefault("wyckoff", {}) + settings.update( + { + "site_id": site_id, + "kind": kind, + selector_name: selector_value, + "value_kind": "delta", + } + ) + return self.addRelParameter( + (atoms, coordinates), + factors, + limits=limits, + **keyargs, + ) + def addWyckoffParameter( self, site_id, @@ -2213,18 +2306,7 @@ def addWyckoffParameter( :raises ValueError: If no matching affine couplings exist. """ - if self.symmetry_metadata is None: - raise ValueError(f"UnitCell {self.name} has no symmetry metadata.") - site = next( - ( - site - for site in self.wyckoff_sites() - if site["site_id"] == site_id - ), - None, - ) - if site is None: - raise ValueError(f"Unknown Wyckoff site {site_id}.") + site = self._wyckoff_site(site_id) if not site["variables"]: raise ValueError( f"Wyckoff site {site_id} has no positional coordinate variables." @@ -2232,53 +2314,29 @@ def addWyckoffParameter( if variable not in site["variables"]: raise ValueError(f"Wyckoff site {site_id} has no variable {variable}.") - if absolute_limits is not None: - if limits != (-np.inf, np.inf): - raise ValueError("Use either limits or absolute_limits, not both.") - variable_value = site["variables"][variable] - limits = ( - absolute_limits[0] - variable_value, - absolute_limits[1] - variable_value, - ) - + limits = self._delta_limits( + limits, + absolute_limits, + site["variables"][variable], + ) couplings = [ coupling for coupling in self.wyckoff_couplings(site_id) if coupling.variable == variable ] - if not couplings: - raise ValueError( + return self._add_wyckoff_relative_parameter( + site_id, + "variable", + variable, + "coordinate", + couplings, + limits, + f"{self.name} {site_id}_{variable}_wyckoff_coordinate", + ( f"No affine couplings found for Wyckoff site {site_id} " f"and variable {variable}." - ) - atoms = np.asarray( - [coupling.atom_index for coupling in couplings], - dtype=np.intp, - ) - coordinates = tuple(coupling.coordinate for coupling in couplings) - factors = np.asarray( - [coupling.factor for coupling in couplings], - dtype=np.float64, - ) - - keyargs.setdefault( - "name", - f"{self.name} {site_id}_{variable}_wyckoff_coordinate", - ) - settings = keyargs.setdefault("wyckoff", {}) - settings.update( - { - "site_id": site_id, - "kind": "coordinate", - "variable": variable, - "value_kind": "delta", - } - ) - return self.addRelParameter( - (atoms, coordinates), - factors, - limits=limits, - **keyargs, + ), + keyargs, ) def addWyckoffParameters( @@ -2310,10 +2368,7 @@ def addWyckoffParameters( :raises ValueError: If the site has no free coordinate variables. """ - sites = {site["site_id"]: site for site in self.wyckoff_sites()} - if site_id not in sites: - raise ValueError(f"Unknown Wyckoff site {site_id}.") - site_variables = tuple(sites[site_id]["variables"]) + site_variables = tuple(self._wyckoff_site(site_id)["variables"]) if not site_variables: raise ValueError( f"Wyckoff site {site_id} has no positional coordinate variables." @@ -2322,19 +2377,14 @@ def addWyckoffParameters( variables = site_variables variables = tuple(variables) - def limit_for(limit_spec, variable_name): - if isinstance(limit_spec, dict): - return limit_spec.get(variable_name, (-np.inf, np.inf)) - return limit_spec - parameters = [] for variable in variables: parameter_keyargs = copy.deepcopy(keyargs) - parameter_limits = limit_for(limits, variable) + parameter_limits = self._limit_for(limits, variable) parameter_absolute_limits = ( None if absolute_limits is None - else limit_for(absolute_limits, variable) + else self._limit_for(absolute_limits, variable) ) parameters.append( self.addWyckoffParameter( @@ -2379,73 +2429,39 @@ def addWyckoffShift( :raises ValueError: If no matching site-displacement couplings exist. """ - if self.symmetry_metadata is None: - raise ValueError(f"UnitCell {self.name} has no symmetry metadata.") if axis not in {"x", "y", "z"}: raise ValueError("axis must be one of 'x', 'y', or 'z'.") - site = next( - ( - site - for site in self.wyckoff_sites() - if site["site_id"] == site_id - ), - None, - ) - if site is None: - raise ValueError(f"Unknown Wyckoff site {site_id}.") + site = self._wyckoff_site(site_id) representative = site.get("representative_parent_fractional") if representative is None: raise ValueError( f"Wyckoff site {site_id} has no representative parent coordinate." ) - if absolute_limits is not None: - if limits != (-np.inf, np.inf): - raise ValueError("Use either limits or absolute_limits, not both.") - axis_index = {"x": 0, "y": 1, "z": 2}[axis] - axis_value = representative[axis_index] - limits = ( - absolute_limits[0] - axis_value, - absolute_limits[1] - axis_value, - ) + axis_index = {"x": 0, "y": 1, "z": 2}[axis] + limits = self._delta_limits( + limits, + absolute_limits, + representative[axis_index], + ) couplings = [ coupling for coupling in self.wyckoff_site_couplings(site_id) if coupling.axis == axis ] - if not couplings: - raise ValueError( + return self._add_wyckoff_relative_parameter( + site_id, + "axis", + axis, + "site_displacement", + couplings, + limits, + f"{self.name} {site_id}_{axis}_wyckoff_site", + ( f"No site-displacement couplings found for Wyckoff site " f"{site_id} and parent axis {axis}." - ) - atoms = np.asarray( - [coupling.atom_index for coupling in couplings], - dtype=np.intp, - ) - coordinates = tuple(coupling.coordinate for coupling in couplings) - factors = np.asarray( - [coupling.factor for coupling in couplings], - dtype=np.float64, - ) - - keyargs.setdefault( - "name", - f"{self.name} {site_id}_{axis}_wyckoff_site", - ) - settings = keyargs.setdefault("wyckoff", {}) - settings.update( - { - "site_id": site_id, - "kind": "site_displacement", - "axis": axis, - "value_kind": "delta", - } - ) - return self.addRelParameter( - (atoms, coordinates), - factors, - limits=limits, - **keyargs, + ), + keyargs, ) def addWyckoffShifts( @@ -2475,19 +2491,14 @@ def addWyckoffShifts( list """ - def limit_for(limit_spec, axis_name): - if isinstance(limit_spec, dict): - return limit_spec.get(axis_name, (-np.inf, np.inf)) - return limit_spec - parameters = [] for axis in axes: parameter_keyargs = copy.deepcopy(keyargs) - parameter_limits = limit_for(limits, axis) + parameter_limits = self._limit_for(limits, axis) parameter_absolute_limits = ( None if absolute_limits is None - else limit_for(absolute_limits, axis) + else self._limit_for(absolute_limits, axis) ) parameters.append( self.addWyckoffShift( @@ -3672,6 +3683,8 @@ def fromStr(cls, string): if reading_symmetry or stripped.startswith( ( "spacegroup:", + "parent_a:", + "parent_alpha:", "surface_transform:", "surface_origin:", "wyckoff_sites:", diff --git a/orgui/datautils/xrayutils/test/test_CTRsymmetry.py b/orgui/datautils/xrayutils/test/test_CTRsymmetry.py index ef500f3..472b5e4 100644 --- a/orgui/datautils/xrayutils/test/test_CTRsymmetry.py +++ b/orgui/datautils/xrayutils/test/test_CTRsymmetry.py @@ -459,6 +459,23 @@ def test_symmetry_metadata_round_trips_without_pyxtal_construction(self): restored.addWyckoffShift("Ru_2a", "x", absolute_limits=(-0.1, 0.1)) self.assertEqual(restored.wyckoff_sites()[0]["status"], "site_displaced") + def test_symmetry_metadata_round_trip_preserves_parent_lattice(self): + unitcell = self.make_rutile_110_unitcell() + + restored = UnitCell.fromStr(unitcell.toStr()) + rebuilt = restored.symmetry_metadata.build_unitcell("rebuilt") + + np.testing.assert_allclose( + restored.symmetry_metadata.surface_spec.parent_a, + unitcell.symmetry_metadata.surface_spec.parent_a, + ) + np.testing.assert_allclose( + restored.symmetry_metadata.surface_spec.parent_alpha, + unitcell.symmetry_metadata.surface_spec.parent_alpha, + ) + np.testing.assert_allclose(rebuilt.a, unitcell.a, atol=1e-8) + np.testing.assert_allclose(rebuilt.alpha, unitcell.alpha, atol=1e-8) + class TestOptionalSymmetryImports(unittest.TestCase): def test_missing_pyxtal_reports_optional_dependency(self): From 549681ad6c5d82b1cca856bfbc7fb8d14727ec43 Mon Sep 17 00:00:00 2001 From: "T. Fuchs" Date: Sat, 11 Jul 2026 20:23:25 -0400 Subject: [PATCH 09/10] fix(phys): preserve deserialized symmetry atoms --- orgui/datautils/xrayutils/CTRsymmetry.py | 12 ++++++++++-- orgui/datautils/xrayutils/test/test_CTRsymmetry.py | 3 +++ 2 files changed, 13 insertions(+), 2 deletions(-) diff --git a/orgui/datautils/xrayutils/CTRsymmetry.py b/orgui/datautils/xrayutils/CTRsymmetry.py index 472c062..e11b570 100644 --- a/orgui/datautils/xrayutils/CTRsymmetry.py +++ b/orgui/datautils/xrayutils/CTRsymmetry.py @@ -518,7 +518,10 @@ def build_unitcell(self, name, layer_behavior="ignore"): a, alpha = self.surface_spec.surface_lattice_parameters() unitcell = UnitCell(a, alpha, name=name, layer_behavior=layer_behavior) - generated = generate_surface_atoms(self.surface_spec, self.sites) + if self.atoms and any(not site.coordinates for site in self.sites): + generated = _reindex_generated_atoms(self.atoms) + else: + generated = generate_surface_atoms(self.surface_spec, self.sites) for atom in generated: site = self._site(atom.site_id) unitcell.addAtom( @@ -874,6 +877,11 @@ def generate_surface_atoms(surface_spec, sites, tolerance=1e-8): round(atom.surface_fractional[0], 10), ) ) + return _reindex_generated_atoms(generated) + + +def _reindex_generated_atoms(atoms): + """Return generated atom metadata with indices matching list order.""" return [ GeneratedWyckoffAtom( atom_index=index, @@ -905,7 +913,7 @@ def generate_surface_atoms(surface_spec, sites, tolerance=1e-8): for coupling in atom.site_couplings ), ) - for index, atom in enumerate(generated) + for index, atom in enumerate(atoms) ] diff --git a/orgui/datautils/xrayutils/test/test_CTRsymmetry.py b/orgui/datautils/xrayutils/test/test_CTRsymmetry.py index 472b5e4..de863b1 100644 --- a/orgui/datautils/xrayutils/test/test_CTRsymmetry.py +++ b/orgui/datautils/xrayutils/test/test_CTRsymmetry.py @@ -473,6 +473,9 @@ def test_symmetry_metadata_round_trip_preserves_parent_lattice(self): restored.symmetry_metadata.surface_spec.parent_alpha, unitcell.symmetry_metadata.surface_spec.parent_alpha, ) + self.assertEqual(len(rebuilt.basis), len(restored.basis)) + self.assertEqual(rebuilt.names, restored.names) + np.testing.assert_allclose(rebuilt.basis[:, 1:4], restored.basis[:, 1:4]) np.testing.assert_allclose(rebuilt.a, unitcell.a, atol=1e-8) np.testing.assert_allclose(rebuilt.alpha, unitcell.alpha, atol=1e-8) From 352719c7a8fe04f94e0ed959237a62cf3e6abb64 Mon Sep 17 00:00:00 2001 From: "T. Fuchs" Date: Sat, 11 Jul 2026 20:51:28 -0400 Subject: [PATCH 10/10] fix(phys): validate integral supercell repeat counts --- orgui/datautils/xrayutils/CTRuc.py | 10 ++++++++-- orgui/datautils/xrayutils/test/test_CTRcalc.py | 6 ++++++ 2 files changed, 14 insertions(+), 2 deletions(-) diff --git a/orgui/datautils/xrayutils/CTRuc.py b/orgui/datautils/xrayutils/CTRuc.py index 7aba1f0..964e941 100644 --- a/orgui/datautils/xrayutils/CTRuc.py +++ b/orgui/datautils/xrayutils/CTRuc.py @@ -1635,9 +1635,15 @@ def supercell(self, repeats, symmetry="preserve", name=None): If repeat counts are not positive integers or ``symmetry`` is unknown. """ - repeats = np.asarray(repeats, dtype=np.intp) - if repeats.shape != (3,) or np.any(repeats <= 0): + repeat_values = np.asarray(repeats, dtype=np.float64) + if ( + repeat_values.shape != (3,) + or not np.all(np.isfinite(repeat_values)) + or np.any(repeat_values <= 0) + or not np.all(repeat_values == np.rint(repeat_values)) + ): raise ValueError("repeats must contain three positive integers.") + repeats = repeat_values.astype(np.intp) if symmetry not in {"preserve", "independent"}: raise ValueError("symmetry must be 'preserve' or 'independent'.") diff --git a/orgui/datautils/xrayutils/test/test_CTRcalc.py b/orgui/datautils/xrayutils/test/test_CTRcalc.py index 5e201e4..d3ea957 100644 --- a/orgui/datautils/xrayutils/test/test_CTRcalc.py +++ b/orgui/datautils/xrayutils/test/test_CTRcalc.py @@ -751,6 +751,12 @@ def test_unitcell_supercell_scales_lattice_and_coordinates(self): ) self.assertEqual(supercell.parameters, {"absolute": [], "relative": []}) + def test_unitcell_supercell_rejects_fractional_repeat_counts(self): + unitcell = self.make_single_wyckoff_cell() + + with self.assertRaisesRegex(ValueError, "positive integers"): + unitcell.supercell((1.5, 1, 1)) + def test_unitcell_supercell_rescales_coherent_domain_translations(self): unitcell = CTRcalc.UnitCell( [2.0, 3.0, 4.0],