.. currentmodule:: crystals
The simplest indexing routine, :func:`index_dirax`, is based on the DirAx algorithm presented in:
For this example, we will need a list of reflections oriented in reciprocal space. Let's start with a very simple example using the simplest crystal structure, cubic polonium:
>>> from crystals import Crystal >>> import numpy as np >>> >>> indices = [ (0,0,0), (1,0,0), (0,1,0), (0,0,1) ] >>> >>> cryst = Crystal.from_database('Pu-epsilon') # cubic polonium >>> peaks = [cryst.scattering_vector(hkl) for hkl in indices]
The list of peaks
are peak positions in three-dimensional reciprocal space. To index:
.. testsetup:: from crystals import Crystal import numpy as np indices = [ (0,0,0), (1,0,0), (0,1,0), (0,0,1) ] cryst = Crystal.from_database('Pu-epsilon') # cubic polonium peaks = [cryst.scattering_vector(hkl) for hkl in indices]
>>> from crystals import index_dirax
>>> lattice, hkls = index_dirax(peaks)
>>> lattice
< Lattice object with parameters 3.638Å, 3.638Å, 3.638Å, 90.00°, 90.00°, 90.00° >
>>> hkls.astype(int)
array([[0, 0, 0],
[1, 0, 0],
[0, 1, 0],
[0, 0, 1]])
Here, hkls
are the indices of the reflections in peaks
according to the indexing lattice, lattice
.
As expected, the Miller indices in hkls
match the ones we created in indices
.
Indexing a cubic lattice with perfect peak placement is pretty easy: we only needed four reflections. The DirAx algorithm is robust against missing reflections, but also alien reflections that do not belong there. In the next example, we will introduce 10% of non-fitting reflections:
.. testsetup:: import numpy as np np.random.seed(0) # ensure reproducible example
>>> from crystals import Crystal, index_dirax
>>> import numpy as np
>>>
>>> cryst = Crystal.from_database('Pu-epsilon') # cubic polonium
>>> indices = list(cryst.bounded_reflections(3))
>>> num_aliens = len(indices) // 10 # 10% alien reflections
>>> aliens = [ np.random.random(size=(3,)) for _ in range(num_aliens) ]
>>>
>>> peaks = [cryst.scattering_vector(hkl) for hkl in indices + aliens]
>>> lattice, hkls = index_dirax(peaks)
>>> lattice
< Lattice object with parameters 3.638Å, 3.638Å, 3.638Å, 90.00°, 90.00°, 90.00° >
>>> hkls.round(decimals=3)
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. , 0. ],
[-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. ],
[ 0.549, 0.715, 0.603],
[ 0.545, 0.424, 0.646]])
As you can see, the last indexed reflections in hkls
don't have integer
Miller indices; those are the non-fitting (alien) reflections!
Of course, the DirAx algorithm is also resistant to noise. Here is an example where we add noise to our perfect peaks:
.. testsetup:: import numpy as np np.random.seed(0)
>>> from crystals import Crystal, Lattice, index_dirax
>>> import numpy as np
>>>
>>> cryst = Crystal.from_database("Pu-epsilon")
>>> indices = list(cryst.bounded_reflections(2))
>>>
>>> peaks = [
... cryst.scattering_vector(hkl) + np.random.normal(0, scale=0.01, size=(3,))
... for hkl in indices
... ]
>>> lattice, hkls = index_dirax(peaks)
>>> hkls.round(decimals=2)
array([[-0.98, -0. , -0. ],
[ 0. , -1. , 0.01],
[ 0. , -0. , -1. ],
[ 0. , 0.01, -0. ],
[ 0. , -0. , 0.99],
[ 0. , 1.01, -0. ],
[ 1. , -0.01, 0. ]])
:func:`index_dirax` can also deal with a restricted set of reflections. For example, in electron diffraction measurements, it might occur that the only reflections measured are of the form \{ (hk0) \}. This happens, for example, when measuring diffraction from layered compounds like graphite.
Let's try to index without prior knowledge first:
>>> from crystals import Crystal, index_dirax >>> import numpy as np >>> from itertools import product >>> >>> cryst = Crystal.from_database('C') # graphite >>> indices = product( ... [-2,-1,0,1,2], #-3 < h < 3 ... [-2,-1,0,1,2], #-3 < k < 3 ... [0] # l = 0 ... ) >>> peaks = [cryst.scattering_vector(hkl) for hkl in indices] >>> lattice, hkls = index_dirax(peaks) Traceback (most recent call last): crystals.indexing.common.IndexingError: No candidate lattice vectors could be determined.
As you can see, there not enough information to index. Let's now use an initial guess:
.. testsetup:: from crystals import Crystal, index_dirax from itertools import product cryst = Crystal.from_database('C') # graphite indices = product([-2,-1,0,1,2],[-2,-1,0,1,2],[0]) peaks = [cryst.scattering_vector(hkl) for hkl in indices]
>>> lattice, hkls = index_dirax(peaks, initial=cryst) # Lattice or Crystal works here
>>> lattice
< Lattice object with parameters 2.464Å, 2.464Å, 6.711Å, 90.00°, 90.00°, 120.00° >
Great!
More specialized and performant indexing methods are planned to be included in :mod:`crystals`. If there's a particular method you would like to see included, please do not hesitate to raise an issue!