Add ReciprocalLatticePoint class storing structure, Miller indices and related quantities

[hakonanes]

Signed-off-by: Håkon Wiik Ånes hwaanes@gmail.com

Release Notes

new feature / improvement
Summary: Add a ReciprocalLatticePoint class storing a crystal structure, associated Miller indices and related quantities for computing structure factors

What does this PR do? Please describe and/or link to an open issue.
This was initially a PR to kikuchipy pyxem/kikuchipy#197, but was moved here because of reasons discussed in #80. My plan is to use this functionality to create a suitable set of Miller indices for EBSD stuff, like getting plane traces on the detector, which then can be used for e.g. visually checking whether an indexing is correct, integrate intensities along certain families of planes, etc.

The ReciprocalLatticePoint class simplifies Miller indices handling, like creating a suitable set of families by considering unique indices under symmetry. When a crystal structure (space group, lattice, atoms) and a set of indices are defined, derived class properties include gspacing, dspacing, structure factors, multiplicity, whether the reflection is allowed. One ReciprocalLatticePoint object stores all indices.

It is based upon EMsoft's gnode/reciprocal lattice point class (see e.g. https://github.com/EMsoft-org/EMsoft/blob/develop/Source/EMsoftLib/typedefs.f90#L1238).

Example usage:

>>> from orix.crystal_map import Phase
>>> from diffpy.structure import Lattice, Atom, Structure
>>> from diffsims.crystallography import ReciprocalLatticePoint
>>> lattice = Lattice(2.8665, 2.8665, 2.8665, 90, 90, 90)
>>> atoms = [Atom("Fe", [0, 0, 0]), Atom("Fe", [0.5, 0.5, 0.5])]
>>> structure = Structure(lattice=lattice, atoms=atoms)
>>> p = Phase("ferrite", space_group=229, structure=structure)
>>> rlp = ReciprocalLatticePoint.from_min_dspacing(p, 0.5)  # Ångstrøm
>>> rlp
ReciprocalLatticePoint (19,)
Phase: ferrite (m-3m)
[[3 3 3]
 [3 3 2]
 [3 3 1]
 [3 3 0]
 [3 2 2]
 [3 2 1]
 [3 2 0]
 [3 1 1]
 [3 1 0]
 [3 0 0]
 [2 2 2]
 [2 2 1]
 [2 2 0]
 [2 1 1]
 [2 1 0]
 [2 0 0]
 [1 1 1]
 [1 1 0]
 [1 0 0]]
>>> rlp.multiplicity
array([ 8, 24, 24, 12, 24, 48, 24, 24, 24,  6,  8, 24, 12, 24, 24,  6,  8,
       12,  6])
>>> rlp.dspacing
array([0.55165818, 0.61113985, 0.6576202 , 0.67564053, 0.69522837,
       0.76610435, 0.79502406, 0.86428227, 0.90646689, 0.9555    ,
       0.82748727, 0.9555    , 1.01346079, 1.17024372, 1.28193777,
       1.43325   , 1.65497455, 2.02692159, 2.8665    ])
>>> rlp.gspacing
array([1.8127167 , 1.63628668, 1.52063455, 1.48007699, 1.43837629,
       1.30530521, 1.25782357, 1.15702941, 1.10318425, 1.04657248,
       1.2084778 , 1.04657248, 0.98671799, 0.85452285, 0.78006907,
       0.69771498, 0.6042389 , 0.493359  , 0.34885749])
>>> rlp[rlp.allowed]
ReciprocalLatticePoint (9,)
Phase: ferrite (m-3m)
[[3 3 2]
 [3 3 0]
 [3 2 1]
 [3 1 0]
 [2 2 2]
 [2 2 0]
 [2 1 1]
 [2 0 0]
 [1 1 0]]
>>> rlp.calculate_structure_factor(method="doyleturner", voltage=20e3)
>>> rlp.structure_factor
array([0.        , 1.32809688, 0.        , 1.6808037 , 0.        ,
       2.18277428, 0.        , 0.        , 2.93909638, 0.        ,
       2.51791255, 0.        , 3.49518383, 4.29238125, 0.        ,
       5.58130236, 0.        , 8.09662875, 0.        ])
>>> rlp[rlp.structure_factor > 0]
ReciprocalLatticePoint (9,)
Phase: ferrite (m-3m)
[[3 3 2]
 [3 3 0]
 [3 2 1]
 [3 1 0]
 [2 2 2]
 [2 2 0]
 [2 1 1]
 [2 0 0]
 [1 1 0]]

Describe alternatives you've considered

• It might be that the symmetrise() method should be moved to a new orix.vector.Miller class inheriting from orix.vector.Vector3d?
• The parts in unique() using "symmetry" should perhaps be moved to Vector3d.unique()?

Are there any known issues? Do you need help?
Don't know the best place for this in the code. I feel that this should be its own module (or part of a new module), available via diffsims.crystallography.ReciprocalLatticePoint. Having a clear import structure is important, and hiding this away in e.g. diffsims.utils.sim_utils.ReciprocalLatticePoint is unclear.

In general, I find the diffsims package structure difficult to navigate. Also, I think the use of an ever growing utils module is a bad idea, because whatever the coder doesn't find a suitable place for is "temporarily" put into that module.

Work in progress?
I've not used diffsims previously, therefore I need help with the following:

• Where does it belong in the code (see my thoughts in the previous section)
• In which parts of the code will this class be used
• What properties and methods should be added? This should preferably be added in later PRs.

Other things:

• Make use of existing atomic scattering parametrizations in diffsims Another PR.
• Generalize get_structure_factor() function to use different atomic scattering parametrizations
• Add a ReciprocalLatticePoint.allowed attribute
• Add a ReciprocalLatticePoint.theta attribute to store the Bragg angle (half the Kikuchi band width) by providing an acceleration voltage
• Docstring examples
• Tests yielding 100% line coverage
• Tests showing desired behaviour for all 7 crystal systems (compared to EMsoft, which has existed for over 20 years)

[dnjohnstone]

Why is this xray? Its using atomic scattering factors for electrons right?

[hakonanes]

Yes, that wording is part of the WIP... I've adopted both the X-ray and Doyle-Turner approach from EMsoft, although only the X-ray Doyle-Turner calculation is included in this PR. Will add it to the TODO list until I find some better way to do it.

[dnjohnstone]

I think in the long run we'd like to support both X-rays and electrons, because it's pretty trivial here.

Maybe a "radiation" attribute for the ReciprocalLatticePoint class and this method can just stay as "get_structure_factor" as that terminology applies to both xrays and electrons so long as the right atomic scattering factors are taken?

[hakonanes]

A radiation attribute sounds like a good idea, I'll add it to the todo list.

Should the get_structure_factor() function take a parametrization parameter to chose which atomic scattering parametrization to produce atomic scattering factors with?

[hakonanes]

Sorry, my comments here are not entirely correct... The get_xray_structure_factor() uses kinematical X-ray scattering factors, while I will push a more involved scattering factor calculation soon, based on EMsoft's approach where they use Doyle-Turner atomic scattering parametrization. That is the reason for the naming. The plan was to include a method parameter in the calculate_structure_factor() either passing xray or doyleturner for different calculations of the factors.

Relevant code in EMsoft https://github.com/EMsoft-org/EMsoft/blob/develop/Source/EMsoftLib/diffraction.f90#L456.

[hakonanes]

I've now added the computation of the atomic scattering factor and structure factor using Doyle-Turner parametrization as done in EMsoft (https://github.com/EMsoft-org/EMsoft/blob/develop/Source/EMsoftLib/diffraction.f90#L496).

[hakonanes]

Considering my comments yesterday were incorrect, I want to postpone a radiation attribute.

[dnjohnstone]

I'm a little confused about this - as far as I can see the only parameterization of the atomic scattering factor that you've put in this PR is the Doyle-Turner one? That 1968 paper considered both X-rays and electrons.

So you mean that get_xray_structure_factor is getting the X-ray version of the Doyle-Turner parameterization and the get_doyleturner_structure_factor is getting the electron version of the Doyle-Turner parameterization using the Mott formula?

Where is the reference to kinematical coming from here? I think the Doyle-Turner parameterization you're using is an evaluation of the atomic scattering factor within the first Born approximation (i.e. kinematical) whether you're taking the X-ray version or the electron version?

To me the atomic scattering factor is just defined as being in the first Born approximation and this makes sense because that quantity is relevant both in kinematical theories of electron diffraction and dynamical theories of electron diffraction from crystals. See e.g. https://iopscience.iop.org/article/10.1088/0034-4885/42/11/002/pdf?casa_token=V2KN7WTLBSoAAAAA:VfkDr5GYF5qm4izwAFbAYgOcTwDVuOEBOjmSRk1QNrnfqEdXEEsrTRUf5HWCXVYtrL07fuQgNyLE

[hakonanes]

I must admit that I've adapted the Doyle-Turner structure factor calculation from EMsoft (https://github.com/pyxem/diffsims/pull/111/files#diff-97bd368539d6d77fc0a03e1bc2c8691bR103). I've made this fact clear by referencing EMsoft in relevant docstrings. The atomic scattering factors are the same, for both methods, yes. The "kinematical" comes from the statement in EMsoft comments here: https://github.com/EMsoft-org/EMsoft/blob/develop/Source/EMsoftLib/diffraction.f90#L384.

[pc494]

Sorry to not have been very forthcoming on some of the diffsims questions here, I think the basic answer is something to the effect of, the architecture will be overhauled (although keeping generators,libraries etc) and thus you kind of have a free reign on this

[hakonanes]

No problem. I've made some design choices in the latest commit:

• Renamed ReciprocalLatticePoint to CrystalPlane
• Renamed reciprocal_lattice_point module to crystallography
• Created a new module diffraction, and moved calculation of atomic scattering factors and structure factors there

[hakonanes]

Thoughts on introducing the diffsims.diffraction and diffsims.crystallography modules:

• Module name diffsims.diffraction might be confusing for users, however, the name accurately accompasses calculation of scattering/structure factors. Other stuff from diffsims.utils.sim_utils should go here.
• Content in diffsims.crystallography should be an extension of functionality in orix. Every PR to this module should be considered for orix.

[dnjohnstone]

• Module name diffsims.diffraction might be confusing for users, however, the name accurately accompasses calculation of scattering/structure factors. Other stuff from diffsims.utils.sim_utils should go here.

I agree with the principle of having a module for code relating to calculating the structure factor. But I'm not keen on the name diffsims.diffraction because the computation of the structure factor is just one way to simulate diffraction. Any objections to either diffsims.scattering_factor or diffsims.structure_factor?

• Content in diffsims.crystallography should be an extension of functionality in orix. Every PR to this module should be considered for orix.

+1 on this.

On some other assorted points:

• Why did you switch from ReciprocalLatticePoint to CrystalPlane? - I preferred the former, but not for any particularly good reason except I like thinking about the 3D reciprocal space when imagining diffraction experiments...

• Don't forget to rip the atomic scattering parameters from here https://github.com/pyxem/diffsims/blob/master/diffsims/utils/sim_utils.py

• In terms of where it will be used - I think it will be in the DiffractionGenerator methods get_ed_data and calculate_profile_data initially - as you can see from even just the naming inconsistency revealed by me writing that they need an overhaul anyway. However, we will merge bugfix for calculate_profile_data #102 before we change those methods again, so it might be best to leave that for a separate PR or to wait a week or two?

• One thing I haven't looked in to is how your calculation of the structure factor compares in terms of performance to the vectorized one in get_ed_data, yours may well be better, but we just ought to check.

Overall this is looking really good and will certainly improve the architecture ultimately.

[hakonanes]

Thanks for your suggestions, @dnjohnstone! Yes, let's figure out where this can be used in another PR.

• I agree that the module name is confusing, I'll change it to diffsims.structure_factor.
• I changed the class name to CrystalPlane because EBSD patterns show Kikuchi bands and not spot patterns. But, assuming the majority if not all diffsims users at the moment do TEM/XRD diffraction and not EBSD, I guess ReciprocalLatticePoint is best? Anyway, I'm okay with either, so I'll revert to ReciprocalLatticePoint.
• Haven't looked into the vectorized structure factor calculation (could someone point to specific code lines?), but vectorized computation should be at least faster. Yes, we should compare the values.
• I'll get the other scattering parameters.

[dnjohnstone]

• I changed the class name to CrystalPlane because EBSD patterns show Kikuchi bands and not spot patterns. But, assuming the majority if not all diffsims users at the moment do TEM/XRD diffraction and not EBSD, I guess ReciprocalLatticePoint is best? Anyway, I'm okay with either, so I'll revert to ReciprocalLatticePoint.

Ok, I would argue that the ReciprocalLattice is just an abstract entity associated with the family of direct lattice planes and therefore just a helpful way to represent the crystal, rather than being linked to spot patterns. Kikuchi bands being lines rather than spots is, at least in my mind, more about the source of electrons that are being scattered to produce them rather than the representation of crystal. But for the same reason, it's not really a big deal.

• Haven't looked into the vectorized structure factor calculation (could someone point to specific code lines?), but vectorized computation should be at least faster. Yes, we should compare the values.

diffsims/diffsims/utils/sim_utils.py

Line 276 in 8482796

f_hkls = np.sum(

[pc494]

Just a thought on naming, is there a reason we've stuck to singulars: CrystalPlane over CrystalPlanes and ReciprocalLatticePoint over ReciprocalLatticePoints (or ReciprocalLattice)?

[hakonanes]

Just a thought on naming, is there a reason we've stuck to singulars: CrystalPlane over CrystalPlanes and ReciprocalLatticePoint over ReciprocalLatticePoints (or ReciprocalLattice)?

I find the case where I have one point in ReciprocalLatticePoint better than one point in ReciprocalLatticePoints. Also, it inherits from the orix.vector.Vector3d, a library with singular class names.

[dnjohnstone]

Yeah, I think it's approximately consistent with orix at least to stay singular?

[hakonanes]

I assume that os, warnings and numpy don't need to be imported in diffsims/__init__.py, when they are not used in that script.

[hakonanes]

I'm not sure how to handle the different atomic scattering parameter parametrizations in diffsims. I find three, namely from Kirkland, Lobato and International Table? And the Doyle-Turner ones will make it four?

We should make a unified API, and perhaps place a module atomic_scattering_parameters in the structure_factor module added in this PR?

• I think the module should have one public function taking in a diffpy.structure.Atom with an element ID (like "H" or "Ni") and which parametrization to use, and return the parameters.
• It can be used by another submodule atomic_scattering_factor, which has one public function taking in a diffpy.structure.Atom, scattering parameter and other relevant parameters and return f for those parameters.

But I think this should be done in a separate PR.

[pc494]

I assume that os, warnings and numpy don't need to be imported in diffsims/__init__.py, when they are not used in that script.

Probably not, I was hoping to automate that sort of clean up in the next batch of coding (0.4 for diffsims)

We should make a unified API, and perhaps place a module atomic_scattering_parameters in the structure_factor module added in this PR?

Yes, this section is due a refactor, again probably in the 0.4 cycle. Do you have a route to make progress here without those changes?

[hakonanes]

Yes, this section is due a refactor, again probably in the 0.4 cycle. Do you have a route to make progress here without those changes?

Well, my route forward is basically not to touch existing atomic scattering parameters, and just write tests for the added functionality as is. So the ReciprocalLatticePoint.calculate_structure_factor() method will in this PR only take "kinematical" and "doyleturner" for different calculations, based on the same atomic scattering factors from Doyle and Turner. The latter calculation takes the accelerating voltage into account.

Apart from that I guess the function naming @dnjohnstone referenced in the review above (#111 (comment)), although I'm not entirely sure which naming is best, as commented below his comment.

[dnjohnstone]

Well, my route forward is basically not to touch existing atomic scattering parameters, and just write tests for the added functionality as is. So the ReciprocalLatticePoint.calculate_structure_factor() method will in this PR only take "kinematical" and "doyleturner" for different calculations, based on the same atomic scattering factors from Doyle and Turner. The latter calculation takes the accelerating voltage into account.

Apart from that I guess the function naming @dnjohnstone referenced in the review above (#111 (comment)), although I'm not entirely sure which naming is best, as commented below his comment.

OK - I think we should just say that these things stay as they are for this PR and we revisit.

So does that get us to just finishing the testing and we merge?

[hakonanes]

OK - I think we should just say that these things stay as they are for this PR and we revisit.

Great. I'll open an issue to discuss the module structure and API after this is merged.

So does that get us to just finishing the testing and we merge?

Yup, writing them now.

[hakonanes]

I'm now satisfied with the class, and coverage is 100%. There are lots of potential improvements for the class. E.g. more tests of symmetrically equivalent planes etc. for different crystal systems should be added...