pyextal.superspace ================== .. py:module:: pyextal.superspace .. autoapi-nested-parse:: Incommensurate crystallography in (3+d)-dimensional superspace. A modulated crystal has no three-dimensional lattice. Its diffraction pattern is still perfectly sharp, but the reflections do not sit on a 3D lattice, so three integers cannot label them. The superspace description of de Wolff, Janner, Janssen and van Smaalen restores a lattice by adding dimensions: every reflection is indexed by ``(h, k, l, m1, ..., md)`` and sits at S(H) = h a* + k b* + l c* + m1 q1 + ... + md qd This subpackage implements that indexing, the modulated structure factor, and a kinematic zone-axis pattern, with 1T-TaS2 as the worked example. Quick start ----------- :: from pyextal.superspace import tas2, saed_pattern, plot_pattern lat, atoms = tas2.structure(d=2) # (3+2)D, incommensurate pattern = saed_pattern(lat, atoms, zone=(0, 0, 1)) plot_pattern(pattern, annotate=12) What is and is not here ----------------------- * ``d = 1, 2, 3`` are all supported. The q-vectors have to be independent modulo the main lattice, and :class:`~pyextal.superspace.lattice. SuperspaceLattice` warns when they are not -- 1T-TaS2's three CDW arms satisfy ``q1 + q2 + q3 = c*``, so they span a rank-2 module and the correct description is (3+2)D, not (3+3)D. * Displacive and occupational modulation, with arbitrarily many harmonics. Modulated ADPs are not implemented. * Structure factors are computed by quadrature over internal space, checked against the Bessel-function closed form to machine precision. * Symmetry: superspace operations, the ``epsilon`` solver, and reflection transformation. Expanding *modulation functions* through an operation is implemented only for the well-behaved case; see :mod:`~pyextal.superspace.symmetry`. * The dynamical side goes through the commensurate approximant (:mod:`~pyextal.superspace.approximant`), which writes an ordinary ``.dat`` file for the **unmodified** Bloch-wave engine. Nothing in this subpackage changes the Fortran, ``pyextal.acom``, or ``jaxEngine``. The 1T-TaS2 modulation amplitudes in :mod:`~pyextal.superspace.tas2` are a labelled stand-in, not the published refinement. See that module. Submodules ---------- .. toctree:: :maxdepth: 1 /autoapi/pyextal/superspace/approximant/index /autoapi/pyextal/superspace/engine/index /autoapi/pyextal/superspace/indexing/index /autoapi/pyextal/superspace/lattice/index /autoapi/pyextal/superspace/modulation/index /autoapi/pyextal/superspace/saed/index /autoapi/pyextal/superspace/scattering/index /autoapi/pyextal/superspace/structurefactor/index /autoapi/pyextal/superspace/symmetry/index /autoapi/pyextal/superspace/tas2/index Classes ------- .. autoapisummary:: pyextal.superspace.ReflectionSet pyextal.superspace.SuperspaceLattice pyextal.superspace.DisplaciveModulation pyextal.superspace.Harmonic pyextal.superspace.ModulatedAtom pyextal.superspace.OccupationalModulation pyextal.superspace.DiffractionPattern pyextal.superspace.SuperspaceOperation Functions --------- .. autoapisummary:: pyextal.superspace.assign_superspace_index pyextal.superspace.dynamical_pattern pyextal.superspace.supercell_atoms pyextal.superspace.supercell_cell pyextal.superspace.supercell_saed pyextal.superspace.write_dat pyextal.superspace.enumerate_reflections pyextal.superspace.format_index pyextal.superspace.coincidence_coset pyextal.superspace.distinct_beams pyextal.superspace.merge_coincident pyextal.superspace.satellite_indices pyextal.superspace.satellite_length pyextal.superspace.satellite_metric pyextal.superspace.cell_to_metric pyextal.superspace.marker_sizes pyextal.superspace.plot_image pyextal.superspace.plot_pattern pyextal.superspace.project_pattern pyextal.superspace.pseudo_voigt_kernel pyextal.superspace.render_image pyextal.superspace.saed_pattern pyextal.superspace.saturation_level pyextal.superspace.electron_form_factor pyextal.superspace.xray_form_factor pyextal.superspace.auto_n_internal pyextal.superspace.bessel_reference pyextal.superspace.internal_grid pyextal.superspace.structure_factor pyextal.superspace.epsilon_for pyextal.superspace.expand_atoms pyextal.superspace.transform_reflection Package Contents ---------------- .. py:function:: assign_superspace_index(lat, hkl_super, M=None, max_order=6, tol=1e-06) Label supercell reflections with their (3+d) superspace indices. This is what makes the dynamical calculation comparable with the kinematic one. The engine works in the supercell and knows nothing about modulation, so it reports a reflection as ``(1, 0, 0)``; in the parent's terms that same reflection is the first-order satellite ``(0 0 0 | 1 0)``. Going back is a matter of undoing ``h' = M h`` and then splitting the result into an integral main index plus a satellite part. The split is not unique for a commensurate structure -- that is the same degeneracy :meth:`SuperspaceLattice.supercell_matrix` exists because of -- so the assignment of lowest satellite order is chosen, which is the physically meaningful one: a reflection is called a first-order satellite rather than a fifth-order one that happens to land in the same place. :param lat: The commensurate parent lattice. :type lat: SuperspaceLattice :param hkl_super: (n, 3) integer supercell reflection indices. :type hkl_super: array-like :param M: Supercell matrix; derived when None. :type M: array-like, optional :param max_order: Largest satellite order to consider. :type max_order: int :param tol: Integrality tolerance on the main part. :type tol: float :returns: ``(H, ok)`` with `H` of shape (n, 3+d) and `ok` a boolean mask. Where `ok` is False no assignment within `max_order` was found and that row of `H` is zero. :rtype: tuple[np.ndarray, np.ndarray] .. py:function:: dynamical_pattern(lat, atoms, thickness=200.0, zone=(0, 0, 1), xaxis=(1, 0, 0), hv=200.0, gmax=0.8, path=None, max_order=6, keep_engine=False, backend='numpy', **dat_kw) A dynamical zone-axis pattern of the approximant, via the Bloch engine. Writes the supercell as a ``.dat``, runs one incident orientation through the **unmodified** engine, and labels the output beams with their superspace indices so the result lines up with :func:`~pyextal.superspace.saed.saed_pattern`. The engine keeps one crystal in Fortran module-global state, so this reparses and resets it. It calls ``terminate()`` at the end unless `keep_engine` is set -- and note that ``terminate()`` ends the engine for the life of the process, so anything else needing it must run first. Leaving it out instead orphans the MPI workers at 100% CPU. :param lat: Must be commensurate. :type lat: SuperspaceLattice :param atoms: The basic structure, listed in P1. :type atoms: list[ModulatedAtom] :param thickness: Specimen thickness in Angstrom. :type thickness: float :param zone: Zone axis in supercell indices. :type zone: tuple :param xaxis: In-plane reference reflection, supercell indices. :type xaxis: tuple :param hv: Accelerating voltage in kV. :type hv: float :param gmax: Beam-selection radius, inverse Angstrom. The supercell reciprocal lattice is dense; raising this is expensive. :type gmax: float :param path: Where to write the ``.dat``; a temporary file when None. :type path: path-like, optional :param max_order: Largest satellite order when labelling reflections. :type max_order: int :param keep_engine: Skip ``terminate()``, leaving the engine live. :type keep_engine: bool :param backend: Eigensolver, ``'numpy'`` (default; any build) or ``'fortran'`` (MPI build only). :type backend: str :param \*\*dat_kw: Passed through to :func:`write_dat`. :returns: With keys `H` (n, 3+d superspace indices), `hkl_super` (n, 3), `intensity` (n,), `order` (n,), `labelled` (n,) bool, `M`, `cell`, `nbeams` and `path`. :rtype: dict .. py:function:: supercell_atoms(lat, atoms, M=None, wrap=True) Freeze the modulation into an ordinary list of 3D atomic positions. Each basic-cell atom is replicated onto every lattice coset of the supercell. Its internal coordinate follows from where it sits, ``xbar_j = q_j . (t + r)``, so each copy gets a different displacement -- which is exactly what makes the supercell a superstructure rather than ``|det M|`` identical copies. :param lat: Must be commensurate. :type lat: SuperspaceLattice :param atoms: The basic structure, listed in P1. :type atoms: list[ModulatedAtom] :param M: Supercell matrix; taken from :meth:`SuperspaceLattice.supercell_matrix` when None. :type M: array-like, optional :param wrap: Wrap fractional coordinates into [0, 1). :type wrap: bool :returns: ``(M, cell, sites)`` with `sites` a list of ``(label, element, position, B)``; `position` is (3,) in supercell fractional coordinates. :rtype: tuple :raises ValueError: If the structure is incommensurate, so no supercell exists. .. py:function:: supercell_cell(lat, M) Lattice parameters of the supercell. :param lat: Supplies the basic metric. :type lat: SuperspaceLattice :param M: (3, 3) integer supercell matrix; rows are the supercell direct lattice vectors in terms of the basic ones. :type M: array-like :returns: ``(a, b, c, alpha, beta, gamma)``, Angstrom and degrees. :rtype: tuple .. py:function:: supercell_saed(lat, atoms, H, thickness=200.0, zone=(0, 0, 1), xaxis=(1, 0, 0), hv=200.0, M=None, path=None, keep_engine=False, backend='numpy', **dat_kw) Dynamical pattern of the **exact supercell**, on the same beams as `H`. The commensurate structure is an ordinary 3D crystal, so this needs nothing from the superspace machinery: the modulation is frozen into ``|det M|`` times as many atomic coordinates and the unmodified engine solves it. What it does need is the *same beam set* as the superspace calculation, or the two dynamical results would not be comparable -- different beams means a different scattering problem. Each superspace reflection maps to the supercell index ``M @ S(H)``, which is integral exactly because the structure is commensurate, and that list is written out explicitly. Agreement between this and :func:`~pyextal.superspace.engine.dynamical_saed` on the same structure is a strong check on the superspace engine: one route freezes the modulation into 507 atoms and never mentions a satellite, the other keeps 3 atoms and integrates over internal space. :param lat: Must be commensurate. :type lat: SuperspaceLattice :param atoms: The basic structure in P1. :type atoms: list[ModulatedAtom] :param H: (n, 3+d) superspace beam list. :type H: array-like :param thickness: Specimen thickness, Angstrom. :type thickness: float :param zone: Zone axis in supercell indices. :type zone: tuple :param xaxis: In-plane reference reflection, supercell indices. :type xaxis: tuple :param hv: Accelerating voltage, kV. :type hv: float :param M: Supercell matrix; derived when None. :type M: array-like, optional :param path: Where to write the ``.dat``. :type path: path-like, optional :param keep_engine: Skip ``terminate()``. :type keep_engine: bool :param backend: Eigensolver, ``'numpy'`` (default; any build) or ``'fortran'`` (MPI build only). :type backend: str :param \*\*dat_kw: Passed to :func:`write_dat`. :returns: ``H`` (the beams, unchanged and in order), ``hkl_super``, ``intensity`` normalised to the strongest diffracted beam, ``pattern`` (a :class:`~pyextal.superspace.saed.DiffractionPattern`, so the result renders through :func:`~pyextal.superspace.saed.plot_image` like any other), ``M``, ``cell``, ``natoms`` and ``nbeams``. :rtype: dict :raises ValueError: If the structure is incommensurate, or a beam does not map to an integral supercell index. .. py:function:: write_dat(path, lat, atoms, zone=(0, 0, 1), xaxis=(1, 0, 0), hv=200.0, conv=0.0, samp=20, M=None, gmax=0.9, sgmax=3.0, sgmin=0.3, omeg=30.0, bmin=0.1, aper=None, abs_pot=0.0, name='approximant', norm=None, kt=(0, 0, 0), header=None, hkl=None) Write the approximant as a ``.dat`` file the Bloch-wave engine can parse. The file uses ``spg 1 1`` -- P1 with every atom listed -- because the supercell's symmetry is lower than the basic structure's, and letting the engine expand a higher group would put atoms in the wrong places. The ``.dat`` format is parsed sequentially and is order-sensitive, so the keywords below are emitted in the order the parser expects; see ``docs/source/configuration.md``. ``sele`` chooses the beams geometrically, which is convenient but means the beam count moves with the voltage and the cell -- and the supercell reciprocal lattice is ``|det M|`` times denser than the basic one, so a `gmax` that is unremarkable for the parent cell can select thousands of beams here. Start small and check what you got. :param path: Destination file. :type path: path-like :param lat: The commensurate lattice. :type lat: SuperspaceLattice :param atoms: Basic-structure atoms, listed in P1. :type atoms: list[ModulatedAtom] :param zone: Zone axis in *supercell* direct indices. With the basis :meth:`SuperspaceLattice.supercell_matrix` returns, the parent's [001] is still (0, 0, 1). :type zone: tuple :param xaxis: In-plane reference reflection, supercell indices. :type xaxis: tuple :param hv: Accelerating voltage in kV. :type hv: float :param conv: Convergence semi-angle in units of ``|g_xaxis|/|K|``. Zero is a parallel beam, which is what a SAED pattern is. :type conv: float :param samp: Sampling across the disc. :type samp: int :param M: Supercell matrix; derived when None. :type M: array-like, optional :param gmax: Beam-selection radius, inverse Angstrom. :type gmax: float :param sgmax: Largest excitation error for beam selection. :type sgmax: float :param sgmin: Excitation error below which a beam is fully diagonalised. :type sgmin: float :param omeg: Bethe-potential parameter; larger is more accurate. :type omeg: float :param bmin: Bethe potential cutoff. :type bmin: float :param aper: Output aperture radius; `gmax` when None. :type aper: float, optional :param abs_pot: Mean absorption potential in volts. :type abs_pot: float :param name: Crystal name written into the file. :type name: str :param norm: Surface normal; `zone` when None. :type norm: tuple, optional :param kt: Tangential incident wavevector. :type kt: tuple :param header: Extra comment lines at the top. :type header: str, optional :param hkl: Explicit beam list in *supercell* indices, written as an ``hkl`` block instead of a ``sele`` line. Use it to make a supercell calculation solve the same scattering problem as a superspace one: pass the images of the superspace beams under `M`, and the two share a beam set so their intensities can be compared reflection by reflection. :type hkl: array-like, optional :returns: ``(path, M, cell, sites)`` -- the file written plus its contents. :rtype: tuple .. py:class:: ReflectionSet Superspace reflections with their derived 3D quantities. :ivar H: (n, 3+d) integer superspace indices. :vartype H: np.ndarray :ivar hkl: (n, 3) effective, generally non-integral, 3D reciprocal indices. :vartype hkl: np.ndarray :ivar gmag: (n,) ``|S(H)|`` in inverse Angstrom. :vartype gmag: np.ndarray :ivar order: (n,) satellite order -- the rank of this reflection's modulation length among the distinct lengths present, so 0 for a main reflection, 1 for the shortest satellites, and so on. Symmetry-equivalent satellites always share an order; see the module docstring for why ``sum |m|`` will not do. :vartype order: np.ndarray :ivar satlen: (n,) ``|sum_j m_j q_j|`` in inverse Angstrom, the quantity `order` ranks. :vartype satlen: np.ndarray :ivar d: Number of modulation dimensions. :vartype d: int .. py:attribute:: H :type: numpy.ndarray .. py:attribute:: hkl :type: numpy.ndarray .. py:attribute:: gmag :type: numpy.ndarray .. py:attribute:: order :type: numpy.ndarray .. py:attribute:: satlen :type: numpy.ndarray .. py:attribute:: d :type: int .. py:property:: index_sum (n,) ``sum |m_j|``. The conventional label, kept because it is what the literature usually quotes. It is not symmetry-invariant -- do not filter on it. :type: np.ndarray .. py:method:: __len__() .. py:method:: __repr__() .. py:method:: select(mask) A new set holding only the reflections `mask` keeps. :param mask: Boolean or integer index array. :type mask: np.ndarray :returns: The subset. :rtype: ReflectionSet .. py:method:: main() ReflectionSet: Just the main reflections. .. py:method:: satellites(order=None) Satellites, optionally of one exact order. :param order: Keep only this order; all satellites if None. :type order: int, optional :returns: The subset. :rtype: ReflectionSet .. py:method:: labels() list[str]: One formatted index per reflection. .. py:function:: enumerate_reflections(lat, gmax, max_index=6, max_order=4, include_origin=True, condition=None) All superspace reflections inside a resolution sphere. :param lat: Supplies the metric and the q-vectors. :type lat: SuperspaceLattice :param gmax: Largest ``|S(H)|`` to keep, inverse Angstrom. :type gmax: float :param max_index: Half-width of the ``(h, k, l)`` search box. :type max_index: int :param max_order: Satellite cutoff, as a multiple of the shortest modulation wave vector: ``|sum_j m_j q_j| <= max_order * min|q_j|``. This is symmetry-invariant, unlike ``sum |m|`` -- see the module docstring. The published 1T-TaS2 data reached order 6. :type max_order: int :param include_origin: Keep ``(0, 0, 0 | 0 ...)``. Defaults to True, since a diffraction pattern has a direct beam. :type include_origin: bool :param condition: Reflection condition -- takes the ``(n, 3+d)`` index array and returns a boolean mask of the allowed reflections. Systematically absent reflections come out with a zero structure factor anyway, since the centred copies cancel, so this saves work and keeps a listing readable rather than correcting anything. See :func:`pyextal.superspace.tas2.reflection_condition`. :type condition: callable, optional :returns: Sorted by ``|S(H)|``, ascending. :rtype: ReflectionSet :raises ValueError: If nothing falls inside `gmax`. .. py:function:: format_index(H) Render one superspace index as ``(h k l | m1 m2)``. The bar keeps the main indices visually separate from the satellite ones, which matters when reading a list where most entries differ only after it. :param H: (3+d,) integer index. :type H: array-like :returns: Formatted label. :rtype: str .. py:function:: coincidence_coset(lat, dh_max=1.05, max_order=24, tol=1e-06) The satellite offsets that leave a reflection where it is. ``c`` belongs to the coset when ``c . q`` is an integer reciprocal-lattice vector, so that ``(h, m)`` and ``(h - c.q, m + c)`` have the same effective index ``S(H)`` and are therefore the same physical reflection. The set is a lattice -- infinite -- so it has to be truncated. **The truncation is on the length of the shift** ``|c . q|``, not on the order ``sum |c_j|``, and that is not a detail. The coset is closed under the crystal's point group, but ``sum |c_j|`` is not invariant under it: cutting there admits some members of a symmetry orbit and not others, which breaks the pattern's symmetry rather than restoring it. Measured on 1T-TaS2's commensurate phase against the exact supercell, an order cut made the agreement *worse* than no coset at all (174% at order 4, against 41% uncorrected) while shell cuts improve monotonically. ``|c . q|`` is a reciprocal-space length and so is invariant, and truncating there keeps whole shells. :func:`merge_coincident` handles the same degeneracy at the level of a finished reflection list; this is the form the *dynamical* calculation needs, because there the coincident amplitudes have to be summed inside every U(g,h) rather than afterwards. For an incommensurate lattice only the identity qualifies, and the returned arrays are the single zero row -- the no-op the engine starts from. :param lat: The lattice. :type lat: SuperspaceLattice :param dh_max: Largest ``|c . q|`` to include, in inverse Angstrom. The default takes the first shell of shifts. :type dh_max: float :param max_order: How far to search in ``c``; raise it only if a shell needs offsets larger than this, which the returned members make obvious. :type max_order: int :param tol: How close ``c . q`` must be to integral. :type tol: float :returns: ``(c, dh)``, each (n, d) and (n, 3) of ints, ordered by increasing ``|c . q|`` with the identity at row 0. :rtype: tuple[np.ndarray, np.ndarray] .. py:function:: distinct_beams(lat, H, who='this calculation', warn=True) Reduce a beam list to one entry per physical reflection. A dynamical calculation must be given distinct reflections. On a commensurate lattice an enumerated list is not: several superspace indices name one reciprocal-lattice point, and handing them over as separate beams breaks the engine twice over. They become identical rows in the structure matrix, which is then singular -- the eigensolver does not report that, it just returns an arbitrary basis in the degenerate subspace. And the `.dat` written from such a list repeats indices in its ``hkl`` block, so ``bloch_parse`` fails to resolve some ``out`` entries against it and *appends* them past the array allocated from the declared ``nbeams``. That one is memory corruption, and it presents from a long way off: an MKL complaint about a `CGEBAL` parameter, a SIGSEGV inside `CGEEV`, or the child simply being killed. Reducing costs nothing physically. `U(g,h)` already sums the coincident amplitudes internally -- see :func:`coincidence_coset` -- so the distinct reflections carry the whole structure. :param lat: The lattice. :type lat: SuperspaceLattice :param H: (n, 3+d) superspace indices, as enumerated. :type H: array-like :param who: Name used in the warning. :type who: str :param warn: Warn when the list actually shrinks. :type warn: bool :returns: (m, 3+d) representatives, m <= n, in the input order. :rtype: np.ndarray .. py:function:: merge_coincident(lat, H, tol=1e-08) Group superspace indices that are the same physical reflection. For an incommensurate structure this never does anything: the map from ``(h k l m1 .. md)`` to ``S(H)`` is injective, which is what makes the indexing unique. For a **commensurate** one it is not -- distinct indices land on the same reciprocal-lattice point, because the Fourier module has collapsed onto a 3D lattice. In 1T-TaS2's commensurate phase, 36 of 85 positions in a modest beam list carry two indices. That matters because the intensity at such a point is ``|F_1 + F_2 + ...|^2``, a *coherent* sum of amplitudes. Treating the indices as independent reflections splits the intensity between them arbitrarily and destroys the pattern's symmetry -- symmetry-equivalent spots end up with different shares. It is the same degeneracy :meth:`SuperspaceLattice.is_commensurate` warns about, seen from the diffraction side. :param lat: The lattice. :type lat: SuperspaceLattice :param H: (n, 3+d) superspace indices. :type H: array-like :param tol: Distance in reciprocal lattice units below which two effective indices count as the same reflection. :type tol: float :returns: ``(representatives, group, ngroups)`` -- `representatives` is (ngroups, 3+d), taking the lowest-order index of each group; `group` is (n,) mapping every input row to its group; and `ngroups` is how many distinct reflections there are. :rtype: tuple[np.ndarray, np.ndarray, int] .. py:function:: satellite_indices(lat, max_order, tol=1e-06) All ``(m_1..m_d)`` whose modulation vector is within `max_order` q-lengths. The cutoff is ``|sum_j m_j q_j| <= max_order * min_j |q_j|`` -- a length, measured in units of the shortest modulation wave vector. Because it is built from the invariant metric it keeps whole symmetry orbits, never part of one. For ``d = 1`` it reduces to ``|m| <= max_order``. :param lat: The lattice. :type lat: SuperspaceLattice :param max_order: Cutoff in units of the shortest q-vector. :type max_order: int :param tol: Relative slack, so a satellite exactly on the boundary is kept rather than lost to rounding. :type tol: float :returns: (n, d) integer satellite indices, including all-zero. :rtype: np.ndarray .. py:function:: satellite_length(lat, H) ``|sum_j m_j q_j|`` for each reflection, in inverse Angstrom. How far the satellite sits from its main reflection -- zero for a main reflection, and equal for every member of a symmetry orbit. :param lat: The lattice. :type lat: SuperspaceLattice :param H: (n, 3+d) superspace indices. :type H: array-like :returns: (n,) lengths. :rtype: np.ndarray .. py:function:: satellite_metric(lat) Metric on satellite indices: ``Q_jk = q_j . q_k`` in inverse Angstrom squared. A symmetry operation permutes the q-vectors, so it preserves this form. Any cutoff or grouping built from it is therefore symmetry-invariant, which ``sum |m|`` is not. :param lat: Supplies the q-vectors and reciprocal metric. :type lat: SuperspaceLattice :returns: (d, d) positive-definite matrix. :rtype: np.ndarray .. py:class:: SuperspaceLattice(cell, q, check_rank=True, tol=1e-06) A 3D lattice plus ``d`` modulation wave vectors. :ivar cell: The six lattice parameters of the *basic* (average) structure. :vartype cell: tuple :ivar q: (d, 3) modulation wave vectors, in reciprocal lattice units of the basic structure. :vartype q: np.ndarray :ivar d: Number of modulation dimensions, 1 to 3. :vartype d: int :ivar gmx: (3, 3) direct metric tensor. :vartype gmx: np.ndarray :ivar gmxr: (3, 3) reciprocal metric tensor. :vartype gmxr: np.ndarray .. py:attribute:: MAX_D :value: 3 .. py:attribute:: cell .. py:attribute:: q .. py:attribute:: gmx .. py:attribute:: gmxr .. py:attribute:: tol .. py:method:: _warn_on_rank() Warn about a rank-deficient q set, distinguishing the two causes. Commensurate q-vectors are rank-deficient too, but that is a legitimate thing to build -- it is the approximant -- so it gets its own, quieter message pointing at the supercell rather than telling the user they made a mistake. .. py:property:: d Number of modulation dimensions. :type: int .. py:property:: dim Total superspace dimension, ``3 + d``. :type: int .. py:method:: __repr__() .. py:method:: effective_hkl(H) The 3D reciprocal-space index a superspace reflection actually sits at. This is the whole point of the superspace description: a satellite is not at an integral position, so it cannot be named by three integers, but it is at a perfectly well defined place. :param H: (n, 3+d) integer superspace indices, or a single (3+d,) index. :type H: array-like :returns: (n, 3) *non-integral* reciprocal indices ``(h, k, l) + sum_j m_j q_j``, in units of the basic reciprocal cell. :rtype: np.ndarray :raises ValueError: If the last axis is not 3+d long. .. py:method:: g_magnitude(H) Length of the scattering vector, inverse Angstrom. :param H: (n, 3+d) superspace indices. :type H: array-like :returns: (n,) values of ``|S(H)|``. :rtype: np.ndarray .. py:method:: d_spacing(H) Interplanar spacing ``1/|S(H)|`` in Angstrom. :param H: (n, 3+d) superspace indices. :type H: array-like :returns: (n,) spacings; ``inf`` for the direct beam. :rtype: np.ndarray .. py:method:: satellite_order(H) :staticmethod: Satellite order ``sum |m_j|``; zero for a main reflection. :param H: (n, 3+d) superspace indices. :type H: array-like :returns: (n,) integer orders. :rtype: np.ndarray .. py:method:: is_main(H) Whether each reflection is a main reflection. :param H: (n, 3+d) superspace indices. :type H: array-like :returns: (n,) booleans. :rtype: np.ndarray .. py:method:: integer_relations(max_n=6, tol=None) Integer combinations of the q-vectors that land on the main lattice. Each one found means two nominally distinct satellites fall at the same place, so the description carries a dimension the data does not. Searches ``n`` in a box of half-width `max_n`, taking the first of each +/- pair and skipping multiples of a relation already found. :param max_n: Half-width of the integer search box. :type max_n: int :param tol: Distance from a reciprocal lattice point, in reciprocal lattice units. The constructor's tolerance when None. :type tol: float, optional :returns: ``(n, g)`` pairs, `n` the (d,) integer coefficients and `g` the (3,) integer reciprocal lattice vector hit. Empty when the q-vectors are properly independent. :rtype: list[tuple] .. py:method:: rationalise(max_denominator=256, tol=None) Exact integer form of the q-vectors, when they are all rational. Returns ``q_j = p_j / N`` with a single common denominator. This is what makes commensurability decidable rather than guessed at: a brute-force search for integer relations would have to reach ``n = 13`` to notice that ``q = (3/13, 1/13, 1/3)`` is rational, and to ``n = 39`` in three dimensions that is half a million candidates. :param max_denominator: Largest denominator to accept. :type max_denominator: int :param tol: How close the rational has to be; the constructor's tolerance when None. :type tol: float, optional :returns: ``(P, N)`` with `P` the (d, 3) integer numerators and `N` the common denominator, or None if any component is irrational to within `tol`. :rtype: tuple or None .. py:method:: rank(tol=None) Rank of the Fourier module beyond the three main dimensions. :param tol: Relation tolerance; the constructor's when None. :type tol: float, optional :returns: The true number of independent modulation dimensions -- ``d`` when the q-vectors are independent, 0 when they are all rational (a commensurate structure, which is a 3D supercell in disguise), and something between when only some directions are redundant. :rtype: int .. py:method:: is_commensurate(tol=None) Whether every q-vector is rational, so an exact 3D supercell exists. :param tol: Rationality tolerance; the constructor's when None. :type tol: float, optional :returns: True when the modulation adds no incommensurate dimension. :rtype: bool .. py:method:: supercell_matrix(max_denominator=64, tol=None) The smallest 3D supercell in which every satellite is a main reflection. For rational q-vectors the modulated structure *is* an ordinary periodic crystal with a larger cell, and that cell is the commensurate approximant -- the one thing about a modulated structure the existing Bloch-wave engine can already handle, since it is plain 3D. The rows of the returned matrix are the supercell's direct lattice vectors in terms of the basic ones, so a reflection's supercell index is ``M @ hkl`` and the volume ratio is ``|det M|``. The rows are the three shortest independent integer vectors ``v`` with every ``v . q_j`` integral, which in a rank-3 lattice is guaranteed to be a basis. :param max_denominator: Largest denominator to consider when rationalising the q-vectors, and the half-width of the search. :type max_denominator: int :param tol: Rationality tolerance; the constructor's when None. :type tol: float, optional :returns: (3, 3) integer matrix, or None if the structure is genuinely incommensurate and no finite supercell exists. :rtype: np.ndarray or None .. py:method:: _relation_text(n, g) Human-readable form of one integer relation, for the warning. .. py:function:: cell_to_metric(cell) Direct-space metric tensor from the six lattice parameters. :param cell: ``(a, b, c, alpha, beta, gamma)`` with lengths in Angstrom and angles in degrees. :type cell: array-like :returns: (3, 3) direct metric tensor ``G``, so that a vector with fractional coordinates ``x`` has length ``sqrt(x @ G @ x)``. :rtype: np.ndarray .. py:class:: DisplaciveModulation(harmonics=(), d=1) Bases: :py:obj:`_FourierSeries` Positional modulation: a displacement field over internal space. :ivar harmonics: Terms with (3,) amplitudes, in fractional coordinates of the basic cell. :vartype harmonics: tuple[Harmonic] :ivar d: Number of modulation dimensions. :vartype d: int .. py:method:: evaluate(xbar) Displacement at each internal coordinate. :param xbar: (npts, d) internal coordinates. :type xbar: np.ndarray :returns: (npts, 3) displacement in fractional coordinates; all zeros when there are no harmonics. :rtype: np.ndarray .. py:method:: max_amplitude() float: Largest single-harmonic amplitude, for grid sizing. .. py:class:: Harmonic One term of a Fourier series in the internal coordinate. :ivar n: (d,) integer harmonic index. :vartype n: tuple :ivar sin: Amplitude of the sine term -- (3,) for a displacive modulation, () for an occupational one. :vartype sin: np.ndarray :ivar cos: Amplitude of the cosine term, same shape as `sin`. :vartype cos: np.ndarray .. py:attribute:: n :type: tuple .. py:attribute:: sin :type: numpy.ndarray .. py:attribute:: cos :type: numpy.ndarray .. py:method:: __post_init__() .. py:property:: order ``sum |n_j|``, the harmonic order. :type: int .. py:class:: ModulatedAtom One atom of the basic structure, with its modulation functions. :ivar label: Free-form name, e.g. ``'Ta1'``. :vartype label: str :ivar element: Symbol used to look up the scattering factor. :vartype element: str :ivar position: (3,) fractional coordinates in the basic cell. :vartype position: np.ndarray :ivar occupancy: Site occupancy of the basic structure; ignored when `occupational` is set, which carries its own average. :vartype occupancy: float :ivar B: Isotropic displacement parameter in Angstrom^2. Modulated ADPs are not implemented -- see the module note in :mod:`~pyextal.superspace.structurefactor`. :vartype B: float :ivar displacive: Positional modulation. :vartype displacive: DisplaciveModulation, optional :ivar occupational: Occupancy modulation. :vartype occupational: OccupationalModulation, optional .. py:attribute:: label :type: str .. py:attribute:: element :type: str .. py:attribute:: position :type: numpy.ndarray .. py:attribute:: occupancy :type: float :value: 1.0 .. py:attribute:: B :type: float :value: 0.0 .. py:attribute:: displacive :type: DisplaciveModulation :value: None .. py:attribute:: occupational :type: OccupationalModulation :value: None .. py:method:: __post_init__() .. py:method:: displacement(xbar) Displacement field, zeros if the atom is not positionally modulated. :param xbar: (npts, d) internal coordinates. :type xbar: np.ndarray :returns: (npts, 3) fractional displacements. :rtype: np.ndarray .. py:method:: occupation(xbar) Occupancy field, constant if the atom is not occupationally modulated. :param xbar: (npts, d) internal coordinates. :type xbar: np.ndarray :returns: (npts,) occupancies. :rtype: np.ndarray .. py:class:: OccupationalModulation(average=1.0, harmonics=(), d=1) Bases: :py:obj:`_FourierSeries` Occupancy as a function of internal coordinate. :ivar average: The ``n = 0`` term, i.e. the mean occupancy. :vartype average: float :ivar harmonics: Terms with scalar amplitudes. :vartype harmonics: tuple[Harmonic] :ivar d: Number of modulation dimensions. :vartype d: int .. py:attribute:: average .. py:method:: evaluate(xbar) Occupancy at each internal coordinate. :param xbar: (npts, d) internal coordinates. :type xbar: np.ndarray :returns: (npts,) occupancies. :rtype: np.ndarray .. py:class:: DiffractionPattern(H, xy, intensity, order, sg, F, zone, xaxis, satlen=None) A computed zone-axis pattern. :ivar H: (n, 3+d) superspace indices. :vartype H: np.ndarray :ivar xy: (n, 2) in-plane coordinates, inverse Angstrom. :vartype xy: np.ndarray :ivar intensity: (n,) relative intensities, normalised so the strongest *diffracted* beam is 1. The direct beam is far above that and is meant to be clipped when plotted. :vartype intensity: np.ndarray :ivar order: (n,) satellite order, ranked by modulation length so symmetry-equivalent satellites share it. See :mod:`~pyextal.superspace.indexing`. :vartype order: np.ndarray :ivar satlen: (n,) ``|sum_j m_j q_j|`` in inverse Angstrom. :vartype satlen: np.ndarray :ivar sg: (n,) excitation errors, inverse Angstrom. :vartype sg: np.ndarray :ivar F: (n,) complex structure factors. :vartype F: np.ndarray :ivar zone: The zone axis. :vartype zone: tuple :ivar xaxis: The in-plane reference reflection. :vartype xaxis: tuple .. py:attribute:: satlen .. py:property:: index_sum (n,) ``sum |m_j|``; not symmetry-invariant, see `order`. :type: np.ndarray .. py:method:: __len__() .. py:method:: __repr__() .. py:method:: visible(rel=0.0001) Reflections above a relative intensity threshold. The threshold is measured against the strongest *diffracted* beam, not the direct one -- intensities are already normalised that way, and the direct beam is orders of magnitude above it, so scaling to that would cut everything. :param rel: Fraction of the strongest diffracted intensity. :type rel: float :returns: Boolean mask. :rtype: np.ndarray .. py:method:: brightest(n=20, exclude_origin=True) Indices of the strongest reflections, for a quick look at a list. :param n: How many to return. :type n: int :param exclude_origin: Skip the direct beam, which dominates. :type exclude_origin: bool :returns: Indices into the pattern arrays. :rtype: np.ndarray .. py:function:: marker_sizes(intensity, size=260, gamma=0.35, floor=1.5) Marker areas for a diffraction pattern, saturating the direct beam. Two things make a linear area-to-intensity map useless here. The direct beam is normalised against the strongest *diffracted* beam and the Mott-Bethe factor diverges as ``s -> 0``, so its "intensity" is of order 1e7 -- scaled linearly it becomes a single marker larger than the axes, and the whole plot comes out a solid block of colour. And satellites are ~1e-3 of a main reflection, so anything that renders the mains sensibly renders the satellites as nothing. Intensity is therefore clipped at 1 (the direct beam saturates, exactly as it does on a recorded plate) and raised to `gamma`, the power stretch normally used to display a diffraction pattern. :param intensity: Relative intensities, mains near 1. :type intensity: array-like :param size: Marker area for a saturated reflection. :type size: float :param gamma: Power stretch. 1.0 is linear; smaller lifts weak reflections into view. :type gamma: float :param floor: Smallest marker area, so a very weak reflection is still a visible dot rather than nothing. :type floor: float :returns: Marker areas in points squared. :rtype: np.ndarray .. py:function:: plot_image(pattern, ax=None, shape=(512, 512), extent=None, sigma=2.0, gamma=1.0, radius=None, stretch=0.4, cmap='inferno', include_direct=True, saturate=True, title=None, rel_min=0.0, colorbar=False, vmax=None) Draw a rendered pattern as an image. Two things have to be handled or the picture is a black square with one white dot. The direct beam is some seven orders of magnitude above the strongest diffracted beam, so the display is scaled to the *diffracted* maximum and the centre is allowed to saturate -- which is what a recorded plate does anyway. And satellites are ~1e-3 of a main reflection, so a power stretch lifts them into view. :param pattern: What to draw. :type pattern: DiffractionPattern :param ax: Target axes. :type ax: matplotlib.axes.Axes, optional :param shape: Image size in pixels. :type shape: tuple :param extent: Half-width of the field, inverse Angstrom. :type extent: float, optional :param sigma: Gaussian width of the point spread, pixels. :type sigma: float :param gamma: Lorentzian width, pixels. :type gamma: float :param radius: Kernel half-width. :type radius: int, optional :param stretch: Display power. 1.0 is linear; smaller lifts the weak reflections. 0.4 shows satellites without washing out the mains. :type stretch: float :param cmap: Colormap. :type cmap: str :param include_direct: Render the direct beam at all. :type include_direct: bool :param saturate: Scale the display to the strongest diffracted beam and let the direct beam clip. Turning this off scales to the true maximum, which is the direct beam, and shows little else. :type saturate: bool :param title: Axes title. :type title: str, optional :param rel_min: Intensity cutoff before rendering. :type rel_min: float :param colorbar: Add a colorbar keyed to the saturation level. :type colorbar: bool :param vmax: Intensity that maps to full scale, in the units of the returned image; everything above it clips. Chosen automatically when None -- the strongest diffracted beam, which is what `saturate` describes. Pass it to brighten (smaller) or dim (larger) a plot, and above all to put several patterns on **one** scale: the automatic choice is made per pattern, so side-by-side panels are otherwise each normalised to their own maximum and cannot be compared. :func:`saturation_level` computes the automatic value without drawing anything. :type vmax: float, optional :returns: ``(ax, image)`` with `image` the unscaled intensity array. :rtype: tuple :raises ValueError: If `vmax` is not positive. .. py:function:: plot_pattern(pattern, ax=None, size=260, rel_min=0.0001, label_orders=True, colors=None, annotate=0, title=None, gamma=0.35) Draw a pattern, distinguishing main reflections from satellite orders. Marker area is a power-stretched, saturating function of intensity; see :func:`marker_sizes` for why a linear map does not work. :param pattern: What to draw. :type pattern: DiffractionPattern :param ax: Target axes; a new figure when None. :type ax: matplotlib.axes.Axes, optional :param size: Marker area for a saturated reflection. :type size: float :param rel_min: Relative intensity below which a reflection is dropped. :type rel_min: float :param label_orders: Show a legend keyed by satellite order. :type label_orders: bool :param colors: Satellite order -> colour. :type colors: dict, optional :param annotate: Label this many of the strongest reflections with their superspace index. Zero draws none. :type annotate: int :param title: Axes title. :type title: str, optional :param gamma: Power stretch on intensity; smaller lifts the satellites. :type gamma: float :returns: The axes drawn on. :rtype: matplotlib.axes.Axes .. py:function:: project_pattern(lat, H, zone=(0, 0, 1), xaxis=None) In-plane coordinates of superspace reflections at a zone axis. Shared by every route that produces a pattern -- kinematic, the superspace engine, the exact supercell -- so their images can be laid side by side and compared pixel for pixel. Two routes computing the same reflection at slightly different coordinates would be a subtle and very annoying bug. :param lat: The lattice. :type lat: SuperspaceLattice :param H: (n, 3+d) superspace indices. :type H: array-like :param zone: Zone axis, direct-space indices. :type zone: tuple :param xaxis: In-plane reference reflection. The shortest in-plane *main* reflection when None, which keeps the axis in the basic structure's terms so it survives a change of q. :type xaxis: tuple, optional :returns: ``(xy, xaxis)`` with `xy` (n, 2) in inverse Angstrom. :rtype: tuple[np.ndarray, tuple] :raises ValueError: If the zone admits no in-plane main reflection. .. py:function:: pseudo_voigt_kernel(radius, sigma, gamma) The pseudo-Voigt point spread the JAX refinement fits, on a pixel grid. Identical in form to ``jaxEngine.postProcess.create_pseudo_voigt_kernel``, including the Thompson-Cox-Hastings width combination and the eta approximation, so a `sigma` and `gamma` refined against experimental data can be handed straight to :func:`render_image` and mean the same thing. Widths are in pixels. The kernel is tapered to zero on a disc rather than left square. Its Lorentzian part decays only as 1/r^2, so at any affordable radius it is still a per cent or so of the peak at the edge -- and truncating that squarely puts a visible square halo around every strong reflection in the convolved image. A raised-cosine taper over the outer quarter removes it and costs nothing, since the discarded intensity is renormalised back. :param radius: Kernel half-width; the kernel is (2r+1) square, with the corners tapered away. :type radius: int :param sigma: Gaussian width, pixels. :type sigma: float :param gamma: Lorentzian width, pixels. :type gamma: float :returns: (2r+1, 2r+1) kernel summing to 1. :rtype: np.ndarray .. py:function:: render_image(pattern, shape=(512, 512), extent=None, sigma=2.0, gamma=1.0, radius=None, rel_min=0.0, include_direct=True) Render a pattern as a 2D intensity image, the way a plate records it. A scatter plot draws markers whose area someone chose; this puts intensity on a pixel grid, so the result is an image in the same sense as the engine's own diffraction output -- something you can threshold, integrate, line-scan or compare against a recorded pattern. Each reflection's intensity is spread bilinearly over the four pixels around its position, so a peak that lands between pixels keeps its centroid and its total counts rather than snapping to a grid point. The whole image is then convolved with the pseudo-Voigt of :func:`pseudo_voigt_kernel`. The direct beam is orders of magnitude above everything else, so plot the result on a log or power scale -- :func:`plot_image` does. :param pattern: What to render. :type pattern: DiffractionPattern :param shape: ``(rows, cols)`` in pixels. :type shape: tuple :param extent: Half-width of the field in inverse Angstrom. Taken from the pattern, with a margin, when None. :type extent: float, optional :param sigma: Gaussian width of the point spread, in pixels. :type sigma: float :param gamma: Lorentzian width, in pixels. :type gamma: float :param radius: Kernel half-width. ``ceil(6*max(sigma, gamma))`` when None -- generous because the Lorentzian component has long tails; :func:`pseudo_voigt_kernel` tapers whatever is left. :type radius: int, optional :param rel_min: Drop reflections below this fraction of the strongest diffracted beam before rendering. :type rel_min: float :param include_direct: Render the direct beam. It saturates everything; turning it off is often the more readable picture. :type include_direct: bool :returns: The ``(rows, cols)`` image, and the ``(left, right, bottom, top)`` extent in inverse Angstrom, ready for ``imshow(..., extent=...)``. :rtype: tuple[np.ndarray, tuple] .. py:function:: saed_pattern(lat, atoms, zone=(0, 0, 1), gmax=1.6, max_index=8, max_order=4, thickness=200.0, hv=200.0, xaxis=None, kind='electron', sg_max=None, n_internal=None, condition=None, merge=None) Compute a kinematic zone-axis pattern. :param lat: The superspace lattice. :type lat: SuperspaceLattice :param atoms: Structure, listed in P1. :type atoms: list[ModulatedAtom] :param zone: Zone axis, direct-space indices. :type zone: tuple :param gmax: Resolution limit, inverse Angstrom. :type gmax: float :param max_index: Half-width of the ``(h, k, l)`` search box. :type max_index: int :param max_order: Largest satellite order. :type max_order: int :param thickness: Specimen thickness in Angstrom, setting the relrod length ``1/t``. Larger means a thinner slice of reciprocal space contributes, so out-of-zone satellites fade. :type thickness: float :param hv: Accelerating voltage in kV, for the Ewald sphere radius. :type hv: float :param xaxis: In-plane reference reflection; the shortest in-plane main reflection when None. :type xaxis: tuple, optional :param kind: ``'electron'`` or ``'xray'`` scattering factors. :type kind: str :param sg_max: Discard reflections with ``|s_g|`` above this. ``5 / thickness`` when None, past which the relrod contributes essentially nothing. :type sg_max: float, optional :param n_internal: Internal quadrature grid; auto when None. :type n_internal: int, optional :param condition: Reflection condition, passed to :func:`~pyextal.superspace.indexing.enumerate_reflections`. :type condition: callable, optional :param merge: Add the structure factors of superspace indices that are the same physical reflection before squaring. Only commensurate lattices have any, and for them it is not optional -- the intensity there is ``|F_1 + F_2|^2``, and treating the indices separately splits it arbitrarily and destroys the pattern's symmetry. Decided from :meth:`SuperspaceLattice.is_commensurate` when None, which is the right answer either way. :type merge: bool, optional :returns: The pattern. :rtype: DiffractionPattern :raises ValueError: If the zone admits no in-plane reflection within `gmax`. .. py:function:: saturation_level(pattern, shape=(512, 512), extent=None, sigma=2.0, gamma=1.0, radius=None, rel_min=0.0) The intensity :func:`plot_image` clips at when it chooses for itself. The peak of the strongest *diffracted* beam once rendered and convolved. Useful for putting two patterns on one scale: `plot_image` picks this separately for each, so two panels drawn side by side are each normalised to their own maximum and their brightnesses cannot be compared. Compute it once and pass it to both as `vmax`:: level = saturation_level(kinematic, shape=(650, 650), sigma=3.0) plot_image(kinematic, vmax=level, shape=(650, 650), sigma=3.0) plot_image(dynamical, vmax=level, shape=(650, 650), sigma=3.0) The render arguments must match the ones `plot_image` is given, since the convolution spreads a fixed number of counts and the peak height therefore depends on the grid and the kernel. :param pattern: The pattern. :type pattern: DiffractionPattern :param shape: Image size in pixels. :type shape: tuple :param extent: Half-width of the field, inverse Angstrom. :type extent: float, optional :param sigma: Gaussian width of the point spread, pixels. :type sigma: float :param gamma: Lorentzian width, pixels. :type gamma: float :param radius: Kernel half-width. :type radius: int, optional :param rel_min: Intensity cutoff before rendering. :type rel_min: float :returns: Peak image intensity excluding the direct beam. :rtype: float .. py:function:: electron_form_factor(element, s, s_floor=1e-05) Electron scattering factor by the Mott-Bethe relation. ``f_e = MOTT * (Z_eff - f_x) / s^2``, which diverges as ``s -> 0``; the engine clamps ``s`` to 1e-5 there ([sfsub.f90:27]) and so does this, so the direct beam gets a large finite value rather than an infinity. ``Z_eff`` is the Cromer-Mann parameterisation's own zero-angle value, ``sum(a) + c``, not the integer nuclear charge -- which is what the Fortran engine uses (``zt`` comes out as 72.95366 for Ta, against Z = 73). The difference is the fit's residual at ``s = 0``, and dividing it by ``s^2`` turns it into a spurious tail that grows without bound as ``s -> 0``: 0.8% of ``f_e`` for Ta already at ``s = 0.125``. Since satellites sit close to their main reflections, that is exactly where they live. Using the integer Z here made the Python disagree with the engine by ~1e-3 across the board. :param element: Element symbol. :type element: str :param s: ``|g| / 2`` in inverse Angstrom. :type s: array-like :param s_floor: Lower clamp on `s`. :type s_floor: float :returns: ``f_e(s)`` in Angstrom, same shape as `s`. :rtype: np.ndarray .. py:function:: xray_form_factor(element, s) X-ray atomic form factor. :param element: Element symbol. :type element: str :param s: ``|g| / 2`` in inverse Angstrom, i.e. sin(theta)/lambda. :type s: array-like :returns: ``f_x(s)`` in electrons, same shape as `s`. :rtype: np.ndarray .. py:function:: auto_n_internal(atoms, hkl, max_order, safety=4, lo=16, hi=256) Choose a grid fine enough that the quadrature does not alias. The integrand carries ``exp(2 pi i m.xbar)``, contributing bandwidth ``max|m|``, and ``exp(2 pi i S.u(xbar))``, whose Jacobi-Anger expansion has significant terms out to roughly ``|z| + 4`` in each harmonic present, where ``z = 2 pi max|S.u|``. The rectangle rule is exact below ``n/2``, so the grid is sized at `safety` times the estimated bandwidth. :param atoms: Supplies the modulation amplitudes. :type atoms: list[ModulatedAtom] :param hkl: (n, 3) effective reciprocal indices. :type hkl: np.ndarray :param max_order: Largest satellite order in the reflection set. :type max_order: int :param safety: Multiplier on the estimated bandwidth. :type safety: int :param lo: Floor on the returned value. :type lo: int :param hi: Ceiling, to keep ``n**d`` affordable. :type hi: int :returns: Points per internal dimension, rounded up to even. :rtype: int .. py:function:: bessel_reference(lat, H, amplitude, element='TA', position=(0.0, 0.0, 0.0), kind='electron', s_floor=1e-05) Closed-form structure factor for one atom with a single harmonic. For ``u(xbar) = A sin(2 pi xbar)`` in (3+1)D, Jacobi-Anger gives Int_0^1 exp{2 pi i [m xbar + S.A sin(2 pi xbar)]} dxbar = J_-m(2 pi S.A) exactly. This is the analytic check on the quadrature in :func:`structure_factor`; see ``test_superspace.py``. :param lat: Must be (3+1)-dimensional. :type lat: SuperspaceLattice :param H: (n, 4) superspace indices. :type H: array-like :param amplitude: (3,) sine amplitude in fractional coordinates. :type amplitude: array-like :param element: Element symbol for the scattering factor. :type element: str :param position: (3,) fractional position of the single atom. :type position: array-like :param kind: ``'electron'`` or ``'xray'``. :type kind: str :param s_floor: Lower clamp on ``s``. :type s_floor: float :returns: (n,) complex structure factors. :rtype: np.ndarray :raises ValueError: If `lat` is not (3+1)-dimensional. .. py:function:: internal_grid(d, n) Uniform sampling of the internal unit cell ``[0, 1)^d``. :param d: Number of modulation dimensions. :type d: int :param n: Points per dimension. :type n: int :returns: (n**d, d) internal coordinates. :rtype: np.ndarray .. py:function:: structure_factor(lat, atoms, H, kind='electron', n_internal=None, max_points=4000000, s_floor=1e-05) Structure factors of a modulated structure. :param lat: Supplies the metric and the q-vectors. :type lat: SuperspaceLattice :param atoms: The structure, listed in P1 -- every atom of the basic cell, not just the asymmetric unit. :type atoms: list[ModulatedAtom] :param H: (n, 3+d) integer superspace indices. :type H: array-like :param kind: ``'electron'`` for Mott-Bethe scattering factors, matching the Bloch-wave engine, or ``'xray'`` for the Cromer-Mann values. :type kind: str :param n_internal: Points per internal dimension. Sized by :func:`auto_n_internal` when None. :type n_internal: int, optional :param max_points: Reflections are processed in chunks so that ``chunk * n_internal**d`` stays under this, bounding memory. :type max_points: int :param s_floor: Lower clamp on ``s = |g|/2``, as the engine uses. :type s_floor: float :returns: (n,) complex structure factors, in the units of the chosen scattering factor. :rtype: np.ndarray :raises ValueError: If `kind` is not recognised, or `H` has the wrong width. .. py:class:: SuperspaceOperation One symmetry operation of a (3+d) superspace group. :ivar R: (3, 3) integer rotation on direct fractional coordinates. :vartype R: np.ndarray :ivar tau: (3,) translation part, fractional. :vartype tau: np.ndarray :ivar epsilon: (d, d) integer action on internal coordinates. :vartype epsilon: np.ndarray :ivar tau_i: (d,) internal translation. :vartype tau_i: np.ndarray :ivar G: (d, 3) integer reciprocal lattice vectors from the q-transformation; all zero for the well-behaved settings. :vartype G: np.ndarray :ivar name: Optional label. :vartype name: str .. py:attribute:: R :type: numpy.ndarray .. py:attribute:: tau :type: numpy.ndarray :value: None .. py:attribute:: epsilon :type: numpy.ndarray :value: None .. py:attribute:: tau_i :type: numpy.ndarray :value: None .. py:attribute:: G :type: numpy.ndarray :value: None .. py:attribute:: name :type: str :value: '' .. py:method:: __post_init__() .. py:property:: d Number of modulation dimensions. :type: int .. py:method:: resolve(lat, max_epsilon=3, tol=1e-06) Fill in `epsilon` and `G` from the lattice, if not already given. :param lat: Supplies the q-vectors. :type lat: SuperspaceLattice :param max_epsilon: Half-width of the integer search for epsilon. :type max_epsilon: int :param tol: Tolerance on the residual being a lattice vector. :type tol: float :returns: A new operation with the internal part filled in. :rtype: SuperspaceOperation :raises ValueError: If `R` does not map the q-set onto itself. .. py:function:: epsilon_for(R, q, max_epsilon=3, tol=1e-06) Solve ``R^T q_j = sum_k epsilon_jk q_k + G_j`` for integers. A 3D operation is only a symmetry of a modulated structure if it permutes the Fourier module, which is what this checks. The search is a small exhaustive one over integer `epsilon` rows -- ``d <= 3`` and the entries are tiny, so there are at most a few hundred candidates. The decomposition is **not unique** when the q-vectors are dependent modulo the main lattice, which is exactly the 1T-TaS2 case: with three CDW arms satisfying ``q1 + q2 + q3 = c*``, any solution can be shifted by a multiple of that relation. The smallest one is returned -- minimum ``|epsilon|`` first, then minimum ``|G|`` -- which picks out the plain permutation with ``G = 0`` rather than some equivalent monster. Taking whichever candidate an unordered search happened to reach first would make the answer depend on iteration order. :param R: (3, 3) integer rotation on direct fractional coordinates. :type R: array-like :param q: (d, 3) modulation wave vectors. :type q: array-like :param max_epsilon: Half-width of the integer search per component. :type max_epsilon: int :param tol: How close the residual has to be to an integer vector. :type tol: float :returns: ``(epsilon, G)``, shapes (d, d) and (d, 3). :rtype: tuple[np.ndarray, np.ndarray] :raises ValueError: If some ``R^T q_j`` cannot be written this way, meaning `R` is not a symmetry of the modulated structure. .. py:function:: expand_atoms(atoms, ops, lat, tol=1e-06) Generate the symmetry-equivalent atoms and their modulation functions. Under ``(R | tau)`` an atom moves to ``R r + tau`` and its displacement field becomes ``u'(xbar) = R u(epsilon^-1 (xbar - tau_i))``. Written out in harmonics, the index maps as ``n -> epsilon^-T n`` and the amplitudes mix by the phase ``phi = 2 pi n' . tau_i``:: A' = R (A cos phi + B sin phi) B' = R (-A sin phi + B cos phi) :param atoms: The asymmetric unit. :type atoms: list[ModulatedAtom] :param ops: Operations, resolved or not. :type ops: list[SuperspaceOperation] :param lat: Supplies the q-vectors for resolving. :type lat: SuperspaceLattice :param tol: Tolerance for recognising a duplicate position. :type tol: float :returns: The expanded set, duplicates removed. :rtype: list[ModulatedAtom] :raises NotImplementedError: If any operation has a nonzero `G`. See the module docstring -- list the structure in P1 instead. .. py:function:: transform_reflection(op, H) Apply a superspace operation to reflection indices. From ``R^T S(H) = (R^T h + G^T m) + sum_k (epsilon^T m)_k q_k``, so the main part picks up the ``G`` contribution and the satellite part transforms by ``epsilon^T``. Valid for any `G`. :param op: Operation with its internal part resolved. :type op: SuperspaceOperation :param H: (n, 3+d) integer superspace indices. :type H: array-like :returns: (n, 3+d) transformed indices. :rtype: np.ndarray :raises ValueError: If `op` has no internal part, or the widths disagree.