pyextal.superspace.engine ========================= .. py:module:: pyextal.superspace.engine .. autoapi-nested-parse:: Driving the Fortran Bloch-wave engine with a modulated structure. This is the Python side of Stage B. The engine carries superspace indices and a modulated structure factor of its own: * ``blochwave.superspace`` holds the modulation dimension ``dmod``, the wave vectors ``qvec``, the satellite indices ``hklm`` parallel to ``difpar.hkl``, and the per-atom harmonics. * ``difpar.gvec`` holds each beam's *effective* reciprocal index, ``hkl + sum_j m_j q_j``, which is generally not integral and so cannot live in ``difpar.hkl``. Every geometry calculation reads ``gvec``, so an ordinary crystal and a modulated one follow one path. * ``sfact_mod`` computes the modulated structure factor by quadrature over internal space, and ``setugh`` calls it instead of ``sfact`` when ``dmod > 0``. ``dmod = 0`` is the default and switches all of it off: an ordinary crystal takes exactly the path it always did, and the result is bit-for-bit unchanged (``test_superspace_engine.py`` pins this). Three limits worth knowing before using it: * **P1 only.** Relating the modulation functions of symmetry-equivalent atoms needs a phase that depends on each atom's own position; that is not implemented, and :func:`load` refuses anything else rather than quietly putting the wrong modulation on every generated atom. * **The beam list must be explicit.** ``sele`` picks beams by enumerating integer indices, which cannot reach a satellite. Write the ``hkl`` block out. * **One calculation per process.** ``terminate()`` ends the engine for the life of the process, and skipping it does not help either. :func:`run_isolated` works around it. Functions --------- .. autoapisummary:: pyextal.superspace.engine.is_loaded pyextal.superspace.engine.load pyextal.superspace.engine.unload pyextal.superspace.engine.structure_factor pyextal.superspace.engine.engine_reflections pyextal.superspace.engine.write_dat pyextal.superspace.engine.dynamical_saed pyextal.superspace.engine.run_isolated Module Contents --------------- .. py:function:: is_loaded() bool: Whether a modulation is currently active in the engine. .. py:function:: load(lat, atoms, H, n_internal=32, coset=None, coset_dh=1.05) Push a modulated structure into the engine and rebuild U(g,h). Call after ``bloch_parse``: the crystal, the atom list and the beam list all have to be in place, because this writes alongside them. `H` supplies the satellite indices for every beam. Its main indices must match the ``hkl`` block of the ``.dat`` row for row -- the engine keeps them in ``difpar.hkl`` and this only adds the satellite part, so a mismatch would silently pair the wrong ``m`` with the wrong ``(h k l)``. It is checked. :param lat: Supplies ``d`` and the wave vectors. :type lat: SuperspaceLattice :param atoms: The structure in P1, in the same order as the ``atom`` lines of the ``.dat``. :type atoms: list[ModulatedAtom] :param H: (nbeams, 3+d) superspace indices, ordered as the ``.dat`` beam list. :type H: array-like :param coset: Sum coincident superspace indices into every U(g,h). Only a commensurate structure has any, and there it is required: several indices name one physical reflection, and the crystal's Fourier coefficient is the coherent sum of their amplitudes. Without it the superspace engine reproduces the *t-averaged* structure rather than the actual crystal, which shows up as satellites off by tens of percent against the exact supercell and as a broken point group. None (the default) turns it on exactly when :meth:`SuperspaceLattice.is_commensurate` says so. :type coset: bool, optional :param coset_dh: Largest coincidence shift ``|c . q|`` to include, in inverse Angstrom. The set is a lattice and so infinite; see :func:`~pyextal.superspace.indexing.coincidence_coset` for why the cut has to be on the shift and not on the satellite order. :type coset_dh: float :param n_internal: Internal-space quadrature points per dimension. The integrand's bandwidth grows with satellite order and modulation amplitude; too few aliases high harmonics onto low ones silently. :type n_internal: int :returns: ``nbeams``, ``natoms``, ``d``, ``n_internal`` and ``harmonics`` (the count loaded per atom). :rtype: dict :raises ValueError: If `H` disagrees with the engine's beam list or atom count, or the modulation dimension is out of range. :raises RuntimeError: If the space group is not P1, or the engine rejects the setup. .. py:function:: unload() Switch modulation off and restore the ordinary U(g,h) matrix. .. py:function:: structure_factor(H, kind=3) Modulated structure factors straight from the Fortran. The engine's own ``sfact_mod``, for checking against :func:`pyextal.superspace.structurefactor.structure_factor` or for reading off individual values. Returns U_g in the engine's units, which is the Python result times ``ratio(3)``. The direct beam is not meaningful here: ``sfact`` divides by ``s^2``, which is exactly zero for ``(0 0 0 | 0 ...)``, and the engine takes U_0 from the mean inner potential instead. :param H: (n, 3+d) superspace indices. :type H: array-like :param kind: 1 x-ray, 2 electron in volts, 3 electron in Angstrom^-2. :type kind: int :returns: ``(elastic, absorptive)``, each (n,) complex. :rtype: tuple[np.ndarray, np.ndarray] :raises RuntimeError: If no modulation is loaded. :raises ValueError: If `H` has the wrong width. .. py:function:: engine_reflections() The engine's current beam list as (3+d) superspace indices. :returns: (nbeams, 3+d) integers, or (nbeams, 3) when no modulation is loaded. :rtype: np.ndarray .. py:function:: write_dat(path, lat, atoms, H, zone=(0, 0, 1), xaxis=(1, 0, 0), hv=200.0, conv=0.0, samp=20, sgmax=3.0, sgmin=0.3, omeg=30.0, out=None, aper=None, abs_pot=0.0, name='modulated', norm=None, kt=(0, 0, 0), header=None) Write a ``.dat`` for a modulated structure: P1, explicit beam list. The file describes the *basic* cell and its atoms. The modulation itself is not in the file -- it goes in through :func:`load` after parsing, which keeps the order-sensitive ``.dat`` grammar untouched. Only the main indices of `H` are written, since that is all ``difpar.hkl`` holds; `load` supplies the satellite part in the same order. :param path: Destination. :type path: path-like :param lat: Supplies the cell. :type lat: SuperspaceLattice :param atoms: The structure in P1. :type atoms: list[ModulatedAtom] :param H: (nbeams, 3+d) superspace indices for the beam list. :type H: array-like :param zone: Zone axis. :type zone: tuple :param xaxis: In-plane reference reflection. :type xaxis: tuple :param hv: Accelerating voltage, kV. :type hv: float :param conv: Convergence semi-angle in units of ``|g_xaxis|/|K|``. :type conv: float :param samp: Sampling across the disc. :type samp: int :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. :type omeg: float :param out: Main indices of the output beams; the direct beam alone when None, and ignored when `aper` is given. :type out: array-like, optional :param aper: Output-aperture radius in inverse Angstrom. Preferred over `out` for a modulated structure, because `out` matches on the three main indices and two beams differing only in their satellite part share those. Note it selects **at parse time, on the main index**: the modulation is loaded afterwards, so ``gvec`` is still ``real(hkl)`` when the aperture is applied. Make it large enough to keep every beam and cut on the effective |g| in Python instead -- :func:`dynamical_saed` does. :type aper: float, optional :param abs_pot: Mean absorption potential, volts. :type abs_pot: float :param name: Crystal name. :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. :type header: str, optional :returns: The file written. :rtype: Path .. py:function:: dynamical_saed(lat, atoms, H, thickness=200.0, zone=(0, 0, 1), xaxis=None, hv=200.0, aper=None, n_internal=32, coset=None, coset_dh=1.05, path=None, keep_engine=False, backend='numpy', **dat_kw) A dynamical zone-axis pattern of a *modulated* structure. The whole of Stage B in one call: write the basic cell out, parse it, push the modulation in, solve one incident orientation, and come back with a :class:`~pyextal.superspace.saed.DiffractionPattern` whose intensities are dynamical rather than ``|F|^2``. It is directly comparable with :func:`~pyextal.superspace.saed.saed_pattern` on the same structure, and renders the same way through :func:`~pyextal.superspace.saed.render_image`. Output beams are chosen by aperture, not by an ``out`` list: ``out`` matches on the three main indices, and two beams differing only in their satellite part share those. The engine holds one crystal in Fortran module-global state, so this reparses and resets it. ``terminate()`` runs at the end unless `keep_engine` is set -- and it ends the engine for the life of the process, so use :func:`run_isolated` when more than one calculation is needed. :param lat: The superspace lattice. :type lat: SuperspaceLattice :param atoms: The structure in P1. :type atoms: list[ModulatedAtom] :param H: (nbeams, 3+d) beam list. :type H: array-like :param thickness: Specimen thickness in Angstrom. :type thickness: float :param zone: Zone axis. :type zone: tuple :param xaxis: In-plane reference reflection; the shortest in-plane main reflection when None. :type xaxis: tuple, optional :param hv: Accelerating voltage, kV. :type hv: float :param aper: Keep only reflections with effective ``|g| <= aper``, applied in Python after the solve. All of `H` when None. :type aper: float, optional :param n_internal: Internal-space quadrature points per dimension. :type n_internal: int :param path: Where to write the ``.dat``; a temporary file when None. :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). Both use the effective beam positions and the modulated structure factors set by :func:`load`. :type backend: str :param \*\*dat_kw: Passed to :func:`write_dat`. :returns: With `intensity` normalised to the strongest diffracted beam, and `F` holding the modulated structure factors. :rtype: DiffractionPattern .. py:function:: run_isolated(func, *args, retries=2, settle=2.0, **kwargs) Run one engine calculation in a fresh interpreter and return its result. The engine keeps a single crystal in Fortran module-global state, and there is no way back from that inside one process: ``terminate()`` ends it for the life of the process, and *skipping* it does not help either -- a second ``bloch_parse`` + ``tilt_run`` hangs, because the parked MPI workers cannot reconcile a re-parsed crystal with a different beam count. One dynamical calculation per process, and that is that. So when a notebook needs to compare two, each goes in its own process:: from pyextal.superspace import engine, supercell_saed a = engine.run_isolated(engine.dynamical_saed, lat, atoms, H, ...) b = engine.run_isolated(supercell_saed, lat, atoms, H, ...) The child is started with :func:`os.posix_spawn`, **not** ``subprocess``. Importing anything from ``pyextal`` runs ``MPI_Init`` in the parent, and ``subprocess`` reaches ``fork()``; forking a process with a live MPI runtime and then spawning MPI again in the descendant corrupts its heap. It shows up as ``Fatal glibc error: malloc.c: assertion failed`` from a child that runs perfectly on its own, which points nowhere near the cause. Two engine runs in quick succession are **flaky** even in separate processes: the second child spawns its MPI workers while the first child's are still tearing down, and it surfaces as ``oneMKL ERROR: Parameter 3 was incorrect on entry to CGEBAL`` or a glibc heap assertion, at roughly one failure in two on back-to-back calls. `retries` and `settle` are a mitigation for that race, not a fix for it -- the calculation itself is deterministic, so a retry either works or the error was real. Arguments and the return value travel by pickle, so `func` has to be importable by name in the child -- a module-level function, not a lambda or a closure -- and the arguments have to be picklable. Everything in this subpackage is. :param func: A module-level function. :type func: callable :param \*args: Positional arguments, pickled. :param retries: Extra attempts after a failure, for the spawn race above. Zero disables it, which is what you want when debugging a genuine failure in `func`. :type retries: int :param settle: Seconds to wait before a retry, letting the previous run's MPI workers exit. :type settle: float :param \*\*kwargs: Keyword arguments, pickled. :returns: Whatever `func` returns. :raises RuntimeError: If every attempt fails; the last child's output is included.