pyextal.adp =========== .. py:module:: pyextal.adp .. autoapi-nested-parse:: Anisotropic displacement parameters and their symmetry constraints. The Bloch engine has always been anisotropic. The elastic scattering factor in ``sfsub.f90`` applies .. math:: T(\mathbf{h}) = \exp(-\mathbf{h}^{T} \beta \mathbf{h}) with a full symmetric :math:`\beta` read from ``blochwave.cryst.atpar`` rows 3..8 (Fortran ``atpar(4:9)``), stored in the order ``(b11, b22, b33, b12, b13, b23)``. What was missing was a Python path to it: :meth:`pyextal.optimize.CoarseOptimize.optimizeDWF` used to write only the three diagonal components, from a single isotropic *B* per atom. ``blochwave.cryst.dw`` is *not* an independent quantity. It is the equivalent isotropic :math:`B_{eq} = \frac{4}{3}\mathrm{Tr}(\beta G)`, which ``xtal0.f90`` derives from :math:`\beta` at parse time, and only the **absorptive** form factor uses it. Both therefore have to be written together, which is what :func:`set_beta` does. Conventions, all taken from the Fortran rather than assumed: =========================== ========================================= :math:`\beta` component order ``b11, b22, b33, b12, b13, b23`` isotropic *B* to :math:`\beta` :math:`\beta = \frac{B}{4} G^{*}` :math:`U_{ij}` to :math:`\beta` :math:`\beta_{ij} = 2\pi^{2} U_{ij} a^{*}_{i} a^{*}_{j}` :math:`\beta` to :math:`B_{eq}` :math:`\frac{4}{3}\mathrm{Tr}(\beta G)` =========================== ========================================= Note that the isotropic mapping uses the **reciprocal** metric tensor :math:`G^{*}` including its off-diagonal elements. ``xatoms`` in ``xtal0.f90`` does exactly this. Using ``1/G_ii`` for the diagonal and leaving the off-diagonals alone -- what the old ``optimizeDWF`` did -- agrees only for orthogonal cells. Symmetry -------- An atom's :math:`\beta` is not free. ``sfsub.f90`` builds the orbit of each asymmetric-unit atom by transforming the *reflection* rather than the atom, so the tensor an orbit-mate effectively carries is :math:`W \beta W^{T}`. For an operation that maps the atom onto itself, that has to reproduce the tensor it started from: .. math:: \beta = W \beta W^{T} \quad \forall (W, \mathbf{w}) \in \mathrm{site\ symmetry} :func:`allowed_beta_basis` turns those conditions into the basis of the subspace they leave free, so a refinement can vary the site's real degrees of freedom and nothing else. For SrTiO3 that is one parameter on Sr (1a) and Ti (1b), and two on O (3c, site symmetry 4/mmm): refining all six on the oxygen would be refining four directions the crystal's symmetry does not have. Attributes ---------- .. autoapisummary:: pyextal.adp.BETA_ROWS pyextal.adp._VOIGT pyextal.adp.TWO_PI_SQ Classes ------- .. autoapisummary:: pyextal.adp.DWFModel Functions --------- .. autoapisummary:: pyextal.adp.beta_matrix pyextal.adp.beta_vector pyextal.adp.isotropic_beta pyextal.adp.uij_to_beta pyextal.adp.beta_to_uij pyextal.adp.equivalent_B pyextal.adp.space_group_ops pyextal.adp.site_symmetry pyextal.adp._transform_matrix pyextal.adp.beta_projector pyextal.adp.allowed_beta_basis pyextal.adp.atom_positions pyextal.adp.atom_labels pyextal.adp.get_beta pyextal.adp.set_beta pyextal.adp.describe pyextal.adp.atom_line Module Contents --------------- .. py:data:: BETA_ROWS :value: (3, 4, 5, 6, 7, 8) .. py:data:: _VOIGT :value: ((0, 0), (1, 1), (2, 2), (0, 1), (0, 2), (1, 2)) .. py:data:: TWO_PI_SQ :value: 19.7392088 .. py:function:: beta_matrix(vec) Expands a 6-component beta into the symmetric 3x3 tensor. :param vec: ``(b11, b22, b33, b12, b13, b23)``. :type vec: array-like :returns: (3, 3) symmetric array. :rtype: np.ndarray .. py:function:: beta_vector(mat) Contracts a symmetric 3x3 tensor to the engine's 6-component order. :param mat: (3, 3) symmetric array. :type mat: np.ndarray :returns: ``(b11, b22, b33, b12, b13, b23)``. :rtype: np.ndarray .. py:function:: isotropic_beta(B, gmxr) Converts an isotropic *B* to the equivalent beta tensor. :param B: Isotropic Debye-Waller factor, Angstrom^2. :type B: float :param gmxr: (3, 3) reciprocal metric tensor, ``param.gmxr``. :type gmxr: np.ndarray :returns: The 6-component beta. :rtype: np.ndarray .. note:: The off-diagonal components are non-zero whenever the cell is not orthogonal, and they matter -- omitting them is not "the isotropic approximation", it is a different tensor. .. py:function:: uij_to_beta(uij, gmxr) Converts crystallographic ``Uij`` to the engine's beta. :param uij: ``(U11, U22, U33, U12, U13, U23)``, Angstrom^2. :type uij: array-like :param gmxr: (3, 3) reciprocal metric tensor. :type gmxr: np.ndarray :returns: The 6-component beta. :rtype: np.ndarray .. py:function:: beta_to_uij(beta, gmxr) Inverse of :func:`uij_to_beta`. :param beta: The 6-component beta. :type beta: array-like :param gmxr: (3, 3) reciprocal metric tensor. :type gmxr: np.ndarray :returns: ``(U11, U22, U33, U12, U13, U23)``, Angstrom^2. :rtype: np.ndarray .. py:function:: equivalent_B(beta, gmx) Equivalent isotropic *B*, matching ``xtal0.f90``. :param beta: The 6-component beta. :type beta: array-like :param gmx: (3, 3) direct metric tensor, ``param.gmx``. :type gmx: np.ndarray :returns: ``B_eq = 4/3 Tr(beta @ gmx)``. :rtype: float .. py:function:: space_group_ops() Reads the space group operations out of the engine. :returns: ``(W, w)`` of shapes ``(n, 3, 3)`` and ``(n, 3)``, acting as ``x' = W x + w`` on fractional coordinates. :rtype: tuple[np.ndarray, np.ndarray] .. note:: Requires a parsed crystal -- call after :func:`pyextal.callBloch.bloch_parse`. ``spgr.gp`` holds only the coset representatives. Two further pieces have to be folded in to get the group ``atmult`` actually generates, and leaving either out silently under-counts site symmetry: * the inversion-related half, when ``spgr.centro`` is set; * the ``spgr.ncsub4`` centring translations in ``spgr.csub4`` -- without them Si in Fd-3m looks like a site of order 6 with two free temperature-factor components instead of order 24 with one, because the operations that return the atom to itself do so only up to a face-centring vector. .. py:function:: site_symmetry(position, ops=None, tol=0.0001) Finds the operations that map an atom onto itself. :param position: (3,) fractional coordinates. :type position: array-like :param ops: ``(W, w)`` from :func:`space_group_ops`. Read from the engine when omitted. :type ops: tuple, optional :param tol: Tolerance on the coordinate match, in fractional units. Defaults to 1e-4. :type tol: float, optional :returns: ``(m, 3, 3)`` array of the rotation parts, always including the identity. :rtype: np.ndarray .. py:function:: _transform_matrix(W) The 6x6 map that ``beta -> W beta W^T`` induces on the 6-vector. .. py:function:: beta_projector(rot) Averages the constraint over a set of operations. :param rot: ``(m, 3, 3)`` rotation parts of a site-symmetry group. :type rot: np.ndarray :returns: (6, 6) projector onto the invariant subspace. :rtype: np.ndarray .. py:function:: allowed_beta_basis(rot, tol=1e-08) Basis of the beta components a site symmetry leaves free. :param rot: ``(m, 3, 3)`` operations of one site, or a list of such arrays to constrain several sites to a shared tensor (the intersection of their invariant subspaces). :type rot: np.ndarray | list :param tol: Singular-value cutoff. Defaults to 1e-8. :type tol: float, optional :returns: ``(6, nfree)`` orthonormal basis. Columns span the allowed tensors; ``nfree`` is 1 for a cubic site and 6 for a site with no symmetry beyond the identity. :rtype: np.ndarray .. py:function:: atom_positions() Fractional coordinates of the asymmetric-unit atoms. :returns: ``(natoms, 3)``. :rtype: np.ndarray .. py:function:: atom_labels() Element label of each asymmetric-unit atom. :returns: ``natoms`` labels, e.g. ``['SR', 'TI', 'O']``. :rtype: list[str] .. py:function:: get_beta(atom) Reads one atom's beta out of the engine. :param atom: Index into the asymmetric unit. :type atom: int :returns: The 6-component beta. :rtype: np.ndarray .. py:function:: set_beta(atom, beta, gmx) Writes one atom's beta into the engine, keeping ``dw`` consistent. :param atom: Index into the asymmetric unit. :type atom: int :param beta: The 6-component beta. :type beta: array-like :param gmx: (3, 3) direct metric tensor, for the ``B_eq`` the absorptive form factor needs. :type gmx: np.ndarray .. note:: ``blochwave.resetugh()`` still has to be called before the change reaches the structure matrix. .. py:function:: describe(gmx=None, gmxr=None, ops=None) Summarises what each atom's symmetry allows. :param gmx: Direct metric tensor, for ``B_eq``. :type gmx: np.ndarray, optional :param gmxr: Reciprocal metric tensor, unused but accepted so callers can pass both. :type gmxr: np.ndarray, optional :param ops: ``(W, w)`` from :func:`space_group_ops`. :type ops: tuple, optional :returns: One entry per atom with its label, position, site-symmetry order, number of free beta components and current values. :rtype: list[dict] .. py:function:: atom_line(label, position, beta, occupancy=1.0) Formats an atom line for a ``.dat`` written with ``dw = bij``. :param label: Element label. :type label: str :param position: (3,) fractional coordinates. :type position: array-like :param beta: The 6-component beta. :type beta: array-like :param occupancy: Site occupancy. Defaults to 1.0. :type occupancy: float, optional :returns: One ``atom`` line. :rtype: str .. note:: A file using these lines must declare ``dw = bij`` on its ``crystal`` line, and then every atom in it needs nine numbers -- the parser reads the whole block one way. .. py:class:: DWFModel(groups, bases, x0, bounds, names, gmx, gmxr, aniso) The free temperature-factor parameters of a refinement. Sits between an optimizer's flat parameter vector and the per-atom tensors the engine reads. Two things it takes care of that a plain "one B per atom" vector cannot: * **Symmetry.** In anisotropic mode each site is refined only along the directions :func:`allowed_beta_basis` says are free, so a cubic site keeps one parameter and a 4/mmm site gets two -- never six. * **Units.** Every coefficient is rescaled to Angstrom^2 of equivalent *B*, so a vector of coefficients sums to :math:`B_{eq}` for that site, the bounds mean the same thing they did for the isotropic refinement, and an optimizer's step length is comparable across parameters. Directions that carry no trace (pure shear) are instead scaled to the magnitude of an isotropic tensor with *B* = 1, and bounded symmetrically about zero. :ivar groups: Atom indices sharing each set of parameters. :vartype groups: list[list[int]] :ivar bases: ``(6, nk)`` rescaled basis per group. :vartype bases: list[np.ndarray] :ivar x0: Starting coefficients, read from the current crystal. :vartype x0: np.ndarray :ivar bounds: Per-coefficient bounds. :vartype bounds: list[tuple] :ivar names: Human-readable name of each coefficient. :vartype names: list[str] .. py:attribute:: groups .. py:attribute:: bases .. py:attribute:: x0 .. py:attribute:: bounds .. py:attribute:: names .. py:attribute:: gmx .. py:attribute:: gmxr .. py:attribute:: aniso .. py:attribute:: slices :value: [] .. py:property:: n Number of free parameters. :type: int .. py:method:: build(param, aniso=False, group='atom', boundary=(0.0, 2.0)) :classmethod: Constructs the model for the crystal currently in the engine. :param param: Supplies ``gmx`` and ``gmxr``. :type param: SimParams :param aniso: Refine the symmetry-allowed anisotropic components rather than one isotropic *B*. Defaults to False. :type aniso: bool, optional :param group: ``'atom'`` for one set of parameters per asymmetric-unit atom, ``'type'`` to tie atoms of the same element together, or an explicit list of index lists. Defaults to ``'atom'``. :type group: str | list, optional :param boundary: Bounds on a *B*-like coefficient. Defaults to (0.0, 2.0). :type boundary: tuple, optional :returns: The model, starting from the crystal's current tensors. :rtype: DWFModel .. py:method:: betas(x) The beta tensor each group takes for a parameter vector. :param x: Coefficients. :type x: array-like :returns: One 6-component beta per group. :rtype: list[np.ndarray] .. py:method:: penalty(x) Rejects unphysical tensors without running a simulation. A beta that is not positive definite means an imaginary mean-square displacement. Catching it here costs microseconds; letting it reach the engine costs a full diagonalization and returns a meaningless number. :param x: Coefficients. :type x: array-like :returns: A penalty above any attainable objective value, or None when every tensor is positive definite. :rtype: float | None .. py:method:: apply(x) Writes a parameter vector into the engine. :param x: Coefficients. :type x: array-like .. note:: ``blochwave.resetugh()`` must follow before the change reaches the structure matrix. .. py:method:: report(x) A one-line-per-group summary of a parameter vector. :param x: Coefficients. :type x: array-like :returns: Formatted table of coefficients, ``B_eq`` and beta. :rtype: str