Redesign of Simulations

[CSSFrancis]

Is your feature request related to a problem? Please describe.
Simulating a Library of Diffraction patterns is (much) slower than it should be. It also could be streamlined to allow for
ASE to be used as well as diffpy.

Describe the solution you'd like
I think there are a couple of things we should do in response to this:

1. Rewrite the Structure class as a diffsims class that extends the functionality of diffpy and ASE?

class Structure:
    @cache
    def reciprocal_space_lattice(max_radius=10):
        pass
    @property
    def lattice():
        pass
    def plot()
        # visualize the atoms using ASE?
        pass

2. Clean up creating orientation lists with new ConstrainedRotation class defined by @din14970

3. Add in some visualization tools to make things a bit more explicit

[hakonanes]

For 1, I suggest to use orix' Phase class. Or extend it. It can also be changed to accommodate new needs. The Phase class was made to combine a Symmetry from orix and a Structure from diffpy.structure.

Note that diffsims also has the ReciprocalLatticeVector class for handling reflectors. This might be of some use?

[CSSFrancis]

For 1, I suggest to use orix' Phase class. Or extend it. It can also be changed to accommodate new needs. The Phase class was made to combine a Symmetry from orix and a Structure from diffpy.structure.

@hakonanes The Phase class looks nice! I am a little curious is there is a way to have the point_group/ space_group initialized from the structure. It seems like spglib might be able to do this. At least this is how ASE handles it.

I like the idea of starting with some list of Phases and returning a list of Phases for some CrystalMap object. It would be nice to abstract out the use of diffpy as I think it is confusing (if only because the name is very similar to diffsims)

Note that diffsims also has the ReciprocalLatticeVector class for handling reflectors. This might be of some use?

Yea I probably need to try and tie some of these things together. I think the currently diffsims has a lot of separate parts but not a great coherent structure so it would be good to start by defining what object we are trying to simulate and then go from there.

[CSSFrancis]

Something like this might work:

The idea would be to combine some of the structure knowledge with some of the functionality to create Zone axis patterns, or create the unique rotations for some crystal structure. It also can handle the reciprocal space points

from orix.crystal_map import Phase as OrixPhase
from diffsims.generators.rotation_list_generators import get_reduced_fundamental_zone_grid
from diffsims.generators.zap_map_generator import generate_zap_map
from diffsims.crystallography.reciprocal_space_sample import ReciprocalSpaceSample


class Phase(OrixPhase):
    def __init__(self, **kwargs):
        super().__init__(**kwargs)
        self._radius = 0
        self._reciprocal_space_lattice = None

    def zap_rotations(self,
                      density: str = "low"):
        """
        Returns a list of rotations that can be used to generate ZAP diffraction
        patterns.

        Parameters
        ----------
        density
            The density of directions to use. Options are 'low' 'high' referring to
            3 for only the corners of the fundamental zone or 7 for the corners,
            midpoints and centroids.
        """
        pass

    def constrained_rotation(self,
                             resolution: float,
                             mesh: str = None,):
        """
        Returns all rotations for some crystal structure reduced to only the unique rotations
        by symmetry.

        Parameters
        ----------
        resolution : float
            The resolution of the grid (degrees).
        mesh : str
            The mesh type to use. Options are 'cuboctahedron', 'icosahedron', 'octahedron',

        """

        return get_reduced_fundamental_zone_grid(resolution,
                                                 mesh,
                                                 self.point_group)

    def reciprocal_space_sample(self,
                                reciprocal_radius: float = 10):
        """
        Returns the reciprocal space lattice vectors for a given radius for the structure.

        This is a cached property, so the first time it is called it will be slow, but
        subsequent calls will be fast.

        Parameters
        ----------
        reciprocal_radius
            The radius of the sphere in reciprocal space (units of reciprocal
            Angstroms) within which reciprocal lattice points are returned
        """
        if reciprocal_radius >= self._radius:
            # recalculate to a higher radius
            # this could be more efficient if we just calculated only the new points
            self._radius = reciprocal_radius
            self._reciprocal_space_sample = ReciprocalSpaceSample.from_radius(self.structure.lattice.reciporical(),
                                                                              reciprocal_radius)
            return self._reciprocal_space_sample
        else:
            # return the existing lattice sliced to the radius
            in_sphere = self._reciprocal_space_sample.spot_distances < reciprocal_radius
            return self._reciprocal_space_sample[in_sphere]

    def plot(self):
        """
        Plots the structure
        """
        pass

[hakonanes]

I am a little curious is there is a way to have the point_group/ space_group initialized from the structure

No, the symmetry of the lattice must be specified. Which use case do you have in mind?

The Phase class is static. As such, it may be considered rather primitive... However, its power lies in convenient storage and access to most quantities related to a crystal structure.

We can initialize a Phase if we have a CIF file using Phase.from_cif() (without importing diffpy.structure). Otherwise we need to specify the structure and lattice "manually". I see a benefit of making it easier to create a Phase instance, though. We could for example allow specifying the lattice and atoms as lists upon initialization. These can be combined into a Structure instance internally.

Regarding your wrapping of Phase:

1. A diffsims class cannot also be named Phase.
2. If plot() is meant to return a plot of the unit cell: This method should be in orix!
3. What will the ReciprocalSpaceSample return?

[CSSFrancis]

I am a little curious is there is a way to have the point_group/ space_group initialized from the structure

No, the symmetry of the lattice must be specified. Which use case do you have in mind?

Ahh never mind. I think I was thinking that the Structure included all of the atoms rather than a reduced representation. I think for most cases just using Phase.from_cif() is the right thing to do.

The Phase class is static. As such, it may be considered rather primitive... However, its power lies in convenient storage and access to most quantities related to a crystal structure.

I agree 100% :)

We can initialize a Phase if we have a CIF file using Phase.from_cif() (without importing diffpy.structure). Otherwise we need to specify the structure and lattice "manually". I see a benefit of making it easier to create a Phase instance, though. We could for example allow specifying the lattice and atoms as lists upon initialization. These can be combined into a Structure instance internally.

Let's keep it as is here and maybe push that off?

Regarding your wrapping of Phase:

1. A diffsims class cannot also be named Phase.

Yea I was struggling with a new name :) I will change that...

2. If plot() is meant to return a plot of the unit cell: This method should be in orix!

Honestly this is the reason that I want ASE. I think that the ASE visualization does a a great job of visualizing some atom. I agree that it should go in orix

3. What will the ReciprocalSpaceSample return?

This is a class the @din14970 defined which simplifies things a bit

from diffpy.structure import Lattice
from functools import cached_property
from typing import Union
import numpy as np


class ReciprocalSpaceSample:

    def __init__(
        self,
        reciprocal_lattice: Lattice,
        spot_indices: np.ndarray,
    ):
        self.__lattice = reciprocal_lattice
        self.__spot_indices = spot_indices

    @property
    def lattice(self) -> Lattice:
        return self.__lattice

    @property
    def spot_indices(self) -> np.ndarray:
        return self.__spot_indices

    @cached_property
    def cartesian_coordinates(self) -> np.ndarray:
        return self.lattice.cartesian(self.spot_indices)

    @cached_property
    def spot_distances(self) -> np.ndarray:
        return self.lattice.dist(self.spot_indices, [0, 0, 0])

    def __getitem__(self, sliced: Union[_T, np.ndarray]) -> ReciprocalSpaceSample:
        return ReciprocalSpaceSample(self.lattice, self.spot_indices[sliced])

    @classmethod
    def from_radius(
        cls,
        reciprocal_lattice: Lattice,
        reciprocal_radius: float,
    ) -> ReciprocalSpaceSample:
        """Finds all reciprocal lattice points inside a given reciprocal sphere.
        Utilised within the DiffractionGenerator.

        Parameters
        ----------
        reciprocal_lattice : diffpy.Structure.Lattice
            The reciprocal crystal lattice for the structure of interest.
        reciprocal_radius  : float
            The radius of the sphere in reciprocal space (units of reciprocal
            Angstroms) within which reciprocal lattice points are returned.

        Returns
        -------
        spot_indices : numpy.array
            Miller indices of reciprocal lattice points in sphere.
        cartesian_coordinates : numpy.array
            Cartesian coordinates of reciprocal lattice points in sphere.
        spot_distances : numpy.array
            Distance of reciprocal lattice points in sphere from the origin.
        """
        n = cls._determine_minimum_repetitions(reciprocal_radius, reciprocal_lattice)
        coordinate_grid = np.mgrid[-n:n+1, -n:n+1, -n:n+1].T.reshape(-1, 3)
        surrounding_grid = cls(reciprocal_radius, coordinate_grid)
        in_sphere = surrounding_grid.spot_distances < reciprocal_radius
        return surrounding_grid[in_sphere]

    @classmethod
    def _determine_minimum_repetitions(
        cls,
        distance: float,
        reciprocal_lattice: Lattice,
    ) -> int:
        """Returns the minimum number of cell repetitions to fill a distance"""
        vertices = np.array([
            [-1, -1, -1],
            [-1, -1,  0],
            [-1, -1,  1],
            [-1,  0, -1],
            [-1,  0,  0],
            [-1,  0,  1],
            [-1,  1, -1],
            [-1,  1,  0],
            [-1,  1,  1],
            [0, -1, -1],
            [0, -1,  0],
            [0, -1,  1],
            [0,  0, -1],
            [0,  0,  1],
            [0,  1, -1],
            [0,  1,  0],
            [0,  1,  1],
            [1, -1, -1],
            [1, -1,  0],
            [1, -1,  1],
            [1,  0, -1],
            [1,  0,  0],
            [1,  0,  1],
            [1,  1, -1],
            [1,  1,  0],
            [1,  1,  1]]
        )
        lengths = np.abs(reciprocal_lattice.dist(vertices, [0, 0, 0]))
        min_length = np.min(lengths)
        return int(np.ceil(distance / min_length))

[hakonanes]

We could change orix' Phase class to use ASE instead of diffpy.structure. It was suggested by @din14970 in pyxem/orix#270.

Is spot_indices a set of (hkl)? If so, this class is exactly like ReciprocalLatticeVector! Here's the docs. Have you considered adding these methods to that class instead?

That class has loads of granular functionality for constructing a suitable reflector list. One can construct a reflector list e.g. from minimum interplanar spacing or from highest (hkl). Everything taking the point group symmetry of your phase into account (which reflections are allowed, etc.), of course! from_radius() could then just be a third class method?

There are also methods for calculating kinematical structure factors and Bragg angles.

ReciprocalLatticeVector wraps Miller in orix. Miller represents Vector3d with respect to a Phase. The class might do more than what is necessary for your ReciprocalSpaceSample. What is the intended use of ReciprocalSpaceSample?

[hakonanes]

Here's an example of constructing a suitable reflector list for nickel: https://kikuchipy.org/en/stable/tutorials/kinematical_ebsd_simulations.html#Nickel.

[CSSFrancis]

We could change orix' Phase class to use ASE instead of diffpy.structure. It was suggested by @din14970 in pyxem/orix#270.

Yea it's something that I might consider trying to do. I think I might try to first change diffsims so that it only uses the orix phase class first and then we should be able to change everything in one place.

Is spot_indices a set of (hkl)? If so, this class is exactly like ReciprocalLatticeVector! Here's the docs. Have you considered adding these methods to that class instead?

You are correct it should be part of ReciprocalLatticeVector. In fact a lot of diffsims should probably be rewritten based on ReciprocalLatticeVector... now that I look at it I think ReciprocalLatticeVector is pretty criminally underused.

That class has loads of granular functionality for constructing a suitable reflector list. One can construct a reflector list e.g. from minimum interplanar spacing or from highest (hkl). Everything taking the point group symmetry of your phase into account (which reflections are allowed, etc.), of course! from_radius() could then just be a third class method?

The from_radius() function would be equivalent to from_min_dspacing (or to the reciprocal of from_min_dspacing)

from_radius()
There are also methods for calculating kinematical structure factors and Bragg angles.

ReciprocalLatticeVector wraps Miller in orix. Miller represents Vector3d with respect to a Phase. The class might do more than what is necessary for your ReciprocalSpaceSample. What is the intended use of ReciprocalSpaceSample?

The ReciprocalSpaceSample is used to get all of the reciprocal space points and then determine if they intersect the Ewald sphere in the simulation. It should be pretty easy to use the ReciprocalLatticeVector for everything. I just need to figure out the caching so that it doesn't get recalculated when doing many simulations.

[din14970]

@CSSFrancis Thanks for picking up where I left off. Simulating the templates faster and making the simulation more generic (e.g. accounting for projection geometry/detector misalignment) is initially why I started rewriting some of these things. @hakonanes is probably right that many concepts in diffsims should just be thrown out, and instead be replaced with something from orix.

I implemented a kinematical diffraction simulator in another project about two years ago and was thinking to add some of those concepts in here, but again never got around to it. In it I use ASE to handle the crystal structures. It's only for simulating and plotting single patterns but with rotation lists it just becomes a loop. I will add you both to the private repo, maybe you can take some inspiration from some things in there.

That said, I think the simulation can be made truly fast by compiling the entire thing with something compiled (numba, C, Cpp, rust, ...). When everything is purely numpy + python a lot of time is wasted doing memory allocations. A very loopy repetitive process like simulating template libraries can be made super fast when fully compiled. It's something I'd like to work on when I have time (which will not be soon).

[CSSFrancis]

@CSSFrancis Thanks for picking up where I left off.

You're welcome! It's a collective effort so any bit helps, your start gave me a good idea of where to start at least.

Simulating the templates faster and making the simulation more generic (e.g. accounting for projection geometry/detector misalignment) is initially why I started rewriting some of these things. @hakonanes is probably right that many concepts in diffsims should just be thrown out, and instead be replaced with something from orix.

That's what I'm starting to think as well :)

I implemented a kinematical diffraction simulator in another project about two years ago and was thinking to add some of those concepts in here, but again never got around to it. In it I use ASE to handle the crystal structures. It's only for simulating and plotting single patterns but with rotation lists it just becomes a loop. I will add you both to the private repo, maybe you can take some inspiration from some things in there.

Thank you! That is helpful.

That said, I think the simulation can be made truly fast by compiling the entire thing with something compiled (numba, C, Cpp, rust, ...). When everything is purely numpy + python a lot of time is wasted doing memory allocations. A very loopy repetitive process like simulating template libraries can be made super fast when fully compiled. It's something I'd like to work on when I have time (which will not be soon).

Yea... It really could be done much much faster. Maybe I'll get the motivation to do that someday.

[hakonanes]

The ReciprocalLatticeVector class is underused because I added it after the spot pattern simulation setup was already in place...

The from_radius() function would be equivalent to from_min_dspacing (or to the reciprocal of from_min_dspacing)

Great! Please let me know if I can help in reviewing additions and changes to this class.

I just need to figure out the caching so that it doesn't get recalculated when doing many simulations.

The ReciprocalLatticeVector class is used in kikuchipy as a thing to create once. It is a high-level class for setup of everything. I would not suggest to initialize many instances in a loop! For computations such as simulation of many patterns, the relevant parts are extracted from the class and used in (fairly) optimized code using Numba and Dask.

@din14970, I think your approach in the mentioned repo is similar to how simulations are handled in kikuchipy. For example, the separation of distinct things such as crystal phase, rotations, and detector. This distinction makes for flexible and iterable workflows. Thanks for sharing.

I've said so before, but I quite like the way simulations are generated in kikuchipy. We do so in three steps:

1. Build a reflector list from a crystal phase and a set of reciprocal lattice vectors.
2. Create a simulator from the reflector list (which carries along the crystal phase). With the simulator, we can visualize the reflectors in the stereographic projection or on the Kikuchi sphere, generate a kinematical master pattern, or plot geometrical simulations on top of a detector.
3. Geometrical simulations are made by giving the simulator a set of rotations and a detector. Kinematical simulations are made by pasing the kinematical master pattern a set of rotations and a detector.

The most pressing current drawback with this setup is the lack of generating a multi-phase map of geometrical or kinematical simulations. diffsims has solved this with the library. Even so, is it possible to use this workflow and separation of work into distinct classes for simulation of spot patterns?

[hakonanes]

With #205 now in, this is done, right? Closing, but please re-open if I'm mistaken.