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..4806cd9 100644 --- a/orgui/datautils/xrayutils/CTRfilm.py +++ b/orgui/datautils/xrayutils/CTRfilm.py @@ -125,6 +125,161 @@ def layer_behaviour(self): def layer_behaviour(self, behavior): self.layer_behavior = behavior + def _wyckoff_target_unitcells(self, kwargs): + if hasattr(self, "uc_top") and hasattr(self, "uc_bottom"): + if "unitcell" not in kwargs: + raise ValueError( + "Missing unit cell name. Provide unit cell name as kwarg " + "'unitcell'" + ) + 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, + **kwargs, + ): + """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 kwargs: + For ``EpitaxyInterface``, provide ``unitcell="top"``, + ``unitcell="bottom"``, or a list of these names. + :returns: + Created parameter or list of parameters. + """ + return self._forward_wyckoff_call( + "addWyckoffParameter", + site_id, + variable, + limits=limits, + absolute_limits=absolute_limits, + **kwargs, + ) + + def addWyckoffParameters( + self, + site_id, + variables=None, + limits=(-np.inf, np.inf), + absolute_limits=None, + **kwargs, + ): + """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 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. + """ + return self._forward_wyckoff_call( + "addWyckoffParameters", + site_id, + variables=variables, + limits=limits, + absolute_limits=absolute_limits, + **kwargs, + ) + + def addWyckoffShift( + self, + site_id, + axis, + limits=(-np.inf, np.inf), + absolute_limits=None, + **kwargs, + ): + """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 kwargs: + For ``EpitaxyInterface``, provide ``unitcell="top"``, + ``unitcell="bottom"``, or a list of these names. + :returns: + Created parameter or list of parameters. + """ + return self._forward_wyckoff_call( + "addWyckoffShift", + site_id, + axis, + limits=limits, + absolute_limits=absolute_limits, + **kwargs, + ) + + def addWyckoffShifts( + self, + site_id, + axes=("x", "y", "z"), + limits=(-np.inf, np.inf), + absolute_limits=None, + **kwargs, + ): + """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 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. + """ + return self._forward_wyckoff_call( + "addWyckoffShifts", + site_id, + axes=axes, + limits=limits, + absolute_limits=absolute_limits, + **kwargs, + ) + @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 new file mode 100644 index 0000000..e11b570 --- /dev/null +++ b/orgui/datautils/xrayutils/CTRsymmetry.py @@ -0,0 +1,1434 @@ +# /*########################################################################## +# +# 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, replace +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:", + "wyckoff_site_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], ...] + 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 + 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 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. + + :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, ...] = () + site_couplings: tuple[WyckoffSiteCoupling, ...] = () + + +@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, + "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, + "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 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. + + :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 + displaced = 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 + ): + 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" + + 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) + 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( + 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, + atom.site_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 = _with_unique_site_ids( + _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_index, coordinate in enumerate(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, + ) + ) + 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, + 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), + site_couplings=tuple(site_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 _reindex_generated_atoms(generated) + + +def _reindex_generated_atoms(atoms): + """Return generated atom metadata with indices matching list order.""" + 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 + ), + 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(atoms) + ] + + +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("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)) + lines.append("surface_origin: " + _format_values(model.surface_spec.origin)) + lines.append("wyckoff_sites:") + 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( + [ + site.site_id, + site.element, + site.wyckoff_label, + variables, + *representative_values, + _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), + ] + ) + ) + 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 + + +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 = {} + site_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:]) + + 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( + [ + [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]) + 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=occ, + iDW=iDW, + oDW=oDW, + ) + ) + + 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_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]) + 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, ())), + site_couplings=tuple(site_couplings_by_atom.get(atom_index, ())), + ) + ) + + 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 + 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( + 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 + ) + 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, + oDW=oDW, + ) + + +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 + 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), ( + {}, + {}, + {}, + ) + 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 _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 + 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("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 + 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..964e941 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 @@ -872,6 +874,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() @@ -1139,6 +1142,774 @@ 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 _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, + 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 = [] + site_parent_shifts = {} + 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( + dataclasses.replace( + coupling, + atom_index=int(old_to_new[coupling.atom_index]), + constant=( + coupling.constant + + coordinate_shift[coordinate_index[coupling.coordinate]] + ), + ) + ) + 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 + site_parent_shifts.setdefault(atom.site_id, parent_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 + ), + layer=layer_map[atom.layer], + couplings=tuple(couplings), + 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 + + 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 _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. + + 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 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. + :rtype: + UnitCell + :raises ValueError: + If the affine linear part is not identity or any requested + translation is not an integer. + """ + xy_translation, z_steps = self._validated_layered_translation(translation) + 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: + 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: + where = self.basis[:, 7] == layer + transformed.layerpos[layer_key] = float( + np.mean(self.basis[where, 3]) + ) + transformed.layerpos[layer_key] += z_translation + old_to_new = np.arange(self.basis.shape[0], dtype=np.intp) + layer_map = {layer: layer for layer in cycle} + return self._finish_layered_transform( + transformed, + old_to_new, + layer_map, + xy_translation, + z_shift_by_layer, + ) + + 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: z_translation for layer in cycle} + transformed = self._apply_reordered_layered_transform( + transformed, + cycle, + xy_translation, + layer_map, + z_shift_by_layer, + ) + + 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) + 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. + """ + xy_translation, z_steps = self._validated_layered_translation(translation) + 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: + 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} + return self._finish_layered_transform( + transformed, + old_to_new, + layer_map, + xy_translation, + z_shift_by_layer, + ) + + 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 + } + transformed = self._apply_reordered_layered_transform( + transformed, + cycle, + xy_translation, + layer_map, + z_shift_by_layer, + ) + + 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) + 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. + """ + 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'.") + + 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 = [] + 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 + ) + 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.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), + 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 + 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 = [] + 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]] + ), + ) + ) + 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, + 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), + site_couplings=tuple(site_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.""" @@ -1377,6 +2148,375 @@ 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 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. + + :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) + + @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, + variable, + limits=(-np.inf, np.inf), + absolute_limits=None, + **keyargs, + ): + """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``. 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`. + :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. + """ + site = self._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}.") + + limits = self._delta_limits( + limits, + absolute_limits, + site["variables"][variable], + ) + couplings = [ + coupling + for coupling in self.wyckoff_couplings(site_id) + if coupling.variable == variable + ] + 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}." + ), + 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. + """ + 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." + ) + if variables is None: + variables = site_variables + variables = tuple(variables) + + parameters = [] + for variable in variables: + parameter_keyargs = copy.deepcopy(keyargs) + parameter_limits = self._limit_for(limits, variable) + parameter_absolute_limits = ( + None + if absolute_limits is None + else self._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 axis not in {"x", "y", "z"}: + raise ValueError("axis must be one of 'x', 'y', or 'z'.") + 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." + ) + 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 + ] + 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}." + ), + 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 + """ + + parameters = [] + for axis in axes: + parameter_keyargs = copy.deepcopy(keyargs) + parameter_limits = self._limit_for(limits, axis) + parameter_absolute_limits = ( + None + if absolute_limits is None + else self._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()) @@ -2279,6 +3419,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 +3681,27 @@ 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:", + "parent_a:", + "parent_alpha:", + "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 +3867,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 +3895,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_CTRcalc.py b/orgui/datautils/xrayutils/test/test_CTRcalc.py index a3e4154..d3ea957 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,681 @@ 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_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], + [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_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() + + 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_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=(), + representative_parent_fractional=(0.25, 0.0, 0.0), + 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", + ), + ), + site_couplings=( + CTRsymmetry.WyckoffSiteCoupling( + atom_index=0, + coordinate="x", + axis="x", + 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_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], + [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_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() + + 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]) + + 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() + + 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]) + + 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 new file mode 100644 index 0000000..de863b1 --- /dev/null +++ b/orgui/datautils/xrayutils/test/test_CTRsymmetry.py @@ -0,0 +1,520 @@ +# /*########################################################################## +# +# 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 + + +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", +) +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 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(), + 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_coordinate_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"]["kind"], "coordinate") + 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_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 + + 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) + 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") + + 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, + ) + 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) + + +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"]