Class methods for Phase based on Lattice/Crystal system

[argerlt]

Description of the change

With Regard to the Docstrings:

The docstrings in Phase have been updated to:

1. Better match numpydoc's guidelines.
2. Clarify some confusing terms and use more consistent language
3. Create a standard copy-paste descriptor for the crystal-to-Cartesian conversion.

The last point was done to clarify a common frustration I was hearing at TMS this year when discussing convention standards, as apparently a lot of synchrotron and neutron data uses a different convention, making it frustrating to align with EBSD.

With Regard to Class Methods:

Adds 7 new methods for creating Phase classes, none of which require importing diffpy:

import orix.crystal_map as ocm

ocm.Phase.triclinic()
ocm.Phase.monoclinic()
ocm.Phase.monoclinic()
ocm.Phase.orthorhombic()
ocm.Phase.tetragonal()
ocm.Phase.cubic()
ocm.Phase.trigonal()
ocm.Phase.hexagonal()

This is a cleaned up version some commits I had in #623. In all seven cases, all or some of the lattice and symmetry information will be filled in, and the remainder will be filled in with defaults. This also means Miller examples can be done without having to import diffpy and then define a Lattice and Structure, which was an early on frustration for me when learning ORIX, and a complaint I've anecdotally heard from time to time from new users.

With Regard to new defaults:

Previously, there was some oddity in how Miller objects worked wrt. a lazily defined Phase

import orix.crystal_map as ocm
import orix.vector as ove

p_cubic = ocm.Phase(point_group='m3m')
p_hex = ocm.Phase(point_group='622')
m1 = ove.Miller.random(phase=p_cubic)
m2 = ove.Miller.random(phase=p_hex)

attempting to print m1 throws an error, as does performing any rotation or plotting, but m2 works just fine. this is because ocm._phase.default_lattice sets a default lattice for trigonal and hexagonal crystals, but not for any other systems. Furthermore, the defaults don't make sense for hexagonal, as a system where a=c would be trigonal.

To fix this, I added defaults that make crystallographic sense for all values. This was also necessary to make class creation methods that don't require importing diffpy to operate properly.

With Regard to new examples:

I would prefer not to make examples for this functionality yet, as I think the better plan is to review/approve this PR, then some additional parts of #623 as a second PR, then handle the examples as part of the rewrite of the Phase tutorials. Writing some examples now just to redo them when we convert the tutorial section seems redundant.

Progress of the PR

• Docstrings for all functions
• Unit tests with pytest for all lines
• Clean code style by running black via pre-commit

For reviewers

• The PR title is short, concise, and will make sense 1 year later.
• New functions are imported in corresponding __init__.py.
• New features, API changes, and deprecations are mentioned in the unreleased
  section in CHANGELOG.rst.
• Contributor(s) are listed correctly in __credits__ in orix/__init__.py and in
  .zenodo.json.

[argerlt]

pre-commit.ci autofix

[CSSFrancis]

@argerlt did you want to link the guidance on f strings in docs? I kind of agree that while things like f strings are nice for some developers it does make the code harder to read sometimes.

[hakonanes]

Thanks for adding these methods, @argerlt! They make many things easier, as you say.

I suggest we stick to the default lattices of

• (1, 1, 1, 90, 90, 90) for all but
• (1, 1, 1, 90, 90, 120) for trigonal and hexagonal

They're valid. They're also the ones used in MTEX.

attempting to print m1 throws an error, as does performing any rotation or plotting, but m2 works just fine

I can print the string representation of the crystal vectors without issue. I'm using v0.15.dev2 (current develop branch). Which error do you get?

[hakonanes]

I suggest to use regular language as much as we can. So I'd rephrase to something like:

The name to give the phase. If not given, a name will be inhereted from structure, if possible. If a phase is given, a copy of that phase is returned instead, and all further arguments are ignored.

So I'd not capitalize, and also not talk about "objects", and rather talk about a phase (Phase), an orientation (Orientation), a crystal vector (Miller) and so on. As long as it's not imprecise, of course, which I don't think it is here.

[argerlt]

Point of clarification: I still think docstrings like this one from Phase.structure are hugely helpful:

        The unit cell is defined using the diffpy
        :class:`~diffpy.structure.Structure` class. For reciprocal
        space calculations, the structure must include a
        :class:`~diffpy.structure.Lattice`. For atomic distance
        calculations, the structure must include one or more
        :class:`~diffpy.structure.atom.Atom`.

As long as you are still on board with this type of cross-linking when specific class types need to be referenced, yes, I agree, regular language is much better.

Actually, maybe a good rule is: if it's important enough to call out a specific class, it's also important enough to link to it's documentation.

[hakonanes]

Yes, intersphinx references are very useful in the built docs and in at least VS Code's docstring preview.

[hakonanes]

This should be under an Examples section.

[hakonanes]

This docstring should be kept. Otherwise, we don't show in the documentation that the name can also be set to a string (the setter isn't included by Sphinx):

Before:

After:

[hakonanes]

To better understand your changes here, could you say why the new docstring is better?

We should use "Return", not "Returns". Numpydoc will complain about "Returns". See e.g. the PEP 257 on one-line docstrings.

[hakonanes]

The "... or not" should be removed, so I think the docstring should be:

Return whether the crystal structure is hexagonal/trigonal.

[argerlt]

My understanding from the numpydoc style guide was that every function should try to have a single line descriptor. This comment and the others like it were an attempt at sticking to that convention.

If you don't like it, we can revert it. If sticking to a convention dilutes the value of the docstring, it's not worth it.

[hakonanes]

I think there are mainly two reasons to keep it short:

• When reading it, you get the gist right away
• It displays nicely in the Sphinx autosummary table: https://orix.readthedocs.io/en/stable/reference/generated/orix.quaternion.html#module-orix.quaternion

I interpret "single line" in that it must be a single sentence and it should be short, but it doesn't have to be a single line. A sentence over two-three lines can still be short, I think.

[argerlt]

That's reasonable. I recently switched to a Ruff linter that attempts to enforce it, but I think the linter is bad.

[hakonanes]

I'm using Ruff in other projects (with a version controlled config file, (#568)), and there I don't experience this. Must be some setting, perhaps related to pydoclint settings?

[hakonanes]

After we agree on how to proceed based on my comments, is it OK if I push some commits?

[argerlt]

I can print the string representation of the crystal vectors without issue. I'm using v0.15.dev2 (current develop branch). Which error do you get?

This sent me on a long journey. Quick summary though: the issue hasn't existed since this commit. Since then, the new if/then/else logic ensures self.structure is always a Structure, and thus there is always a self.structure.name that can be set and accessed, so the error I used to get no longer happens.

A completely unrelated error though I found along the way: You can get a similar error message when you define neither a space group nor point group Edit: Nevermind, this issue is deeper than I thought and shouldn't be solved here. opening a seperate issue.

I suggest we stick to the default lattices of

* (1, 1, 1, 90, 90, 90) for all but

* (1, 1, 1, 90, 90, 120) for trigonal and hexagonal

They're valid. They're also the ones used in MTEX.

Fair. It bothers me that using a default cubic lattice for triclinic systems will produce impossible structure factors and symmetry relations, but an arbitrary dictionary of ad-hoc constants is also bad in retrospect. I'll revert this.

[argerlt]

@argerlt did you want to link the guidance on f strings in docs? I kind of agree that while things like f strings are nice for some developers it does make the code harder to read sometimes.

@CSSFrancis, I don't think I do here, I'm not sure anyone is going to read and gain value from a comment on a style choice in a docstring. That said, it might be worth adding to the style guide page in a later PR.

[argerlt]

pre-commit.ci autofix

[argerlt]

After we agree on how to proceed based on my comments, is it OK if I push some commits?

Yes, push whatever you want. My only concern is making sure we are on the same page w.r.t. this comment, but everything else is good.

I'll also add, I'm not super attached to the updates I made to the docstrings. I felt some of the old ones read roughly and I wanted to try enforcing the one-line style suggestion from numpydoc, but I wont be bothered if you alter or revert several of them.

[hakonanes]

It bothers me that using a default cubic lattice for triclinic systems will produce impossible structure factors and symmetry relations

Hm, this doesn't sound good. Can you provide an example snippet or two?

[argerlt]

It bothers me that using a default cubic lattice for triclinic systems will produce impossible structure factors and symmetry relations

Hm, this doesn't sound good. Can you provide an example snippet or two?

Turns out more accurate concern would be "... impossible structure factors diffraction results and symmetry relations ... ".

Also, I'm adding examples here since you asked, but I think the takeaway is if we DO care about this, we should open a separate issue, and I don't think we should for the sole fact that a default dictionary creates more ad-hoc defaults and assumptions, and will cause ORIX to disagree with other software (MTEX for example) when doing identical lazy operations.

High level though: A cubic lattice cannot have, for example, triclinic symmetry. Any motifs placed on the lattice will result in SG#195 (P23) at minimum. You can show this with a toy example creating a superlattice.

import matplotlib.pyplot as plt
import numpy as np

from diffsims.crystallography import ReciprocalLatticeVector

import orix.crystal_map as ocm
import diffpy.structure as dst
from scipy.stats import norm

# %% create two triclinic phases, but one with an actual triclinic lattice.

# a phase with a default lattice added
default_ph = ocm.Phase('default',
                    space_group = 1,                    
                    )
# a phase with an explicitly triclinic lattice added
triclinic_ph = ocm.Phase('triclinic',
                         space_group = 1,
                    structure = dst.Structure(
                        lattice=dst.Lattice(0.85,1,1.2,87,84,85)
                        )
                    )
# add atoms to both
for phase in [default_ph,triclinic_ph]:
    phase.structure.addNewAtom('Ni',[0,0,0])

ijk = []
for ii in np.arange(-3,3):
    for jj in np.arange(-3,3):
        for kk in np.arange(-3,3):
            ijk.append([ii,jj,kk])
ijk = np.array(ijk)

default_positions = np.einsum('ij,jk',ijk,default_ph.structure.lattice.base)
triclinic_positions = np.einsum('ij,jk',ijk,triclinic_ph.structure.lattice.base)

fig = plt.figure()
ax1 = fig.add_subplot(1,2,1,projection='3d')
ax2 = fig.add_subplot(1,2,2,projection='3d')

ax1.scatter(*default_positions.T)
ax2.scatter(*triclinic_positions.T)

The screen capture doesn't show up well, but the plot on the left clearly has rotational symmetry, while the plot on
the right does not no matter how it is spun. Variations of this could show up any time someone were to define a crystal vector then apply a real space rotation.

A diffraction based example: if I try to simulate diffraction patterns on a 2D flatbed detector, it might look something like this:

import matplotlib.pyplot as plt
import numpy as np

from diffsims.crystallography import ReciprocalLatticeVector

import orix.crystal_map as ocm
import diffpy.structure as dst
from scipy.stats import norm

# %% create two triclinic phases, but one with an actual triclinic lattice.

# a phase with a default lattice added
default_ph = ocm.Phase('default',
                    space_group = 1,                    
                    )
# a phase with an explicitly triclinic lattice added
triclinic_ph = ocm.Phase('triclinic',
                         space_group = 1,
                    structure = dst.Structure(
                        lattice=dst.Lattice(0.85,1,1.2,87,84,85)
                        )
                    )
# add atoms to both
for phase in [default_ph,triclinic_ph]:
    phase.structure.addNewAtom('Ni',[0,0,0])


# %% calculate structure factors using method from Kikuchipy tutorial.
all_ref_a = ReciprocalLatticeVector.from_min_dspacing(default_ph, 0.3)
all_ref_b = ReciprocalLatticeVector.from_min_dspacing(triclinic_ph, 0.3)
all_ref_a.calculate_structure_factor()
all_ref_b.calculate_structure_factor()
F_a = np.abs(all_ref_a.structure_factor)
F_b = np.abs(all_ref_b.structure_factor)
ref_a = all_ref_a[F_a>0.2*F_a]
ref_b = all_ref_b[F_b>0.2*F_b]
# numbers for CHESS 1a3 beamline, semi-arbitrary.
ref_a.calculate_theta(78e3)
ref_b.calculate_theta(78e3)

# create a fake detector (arbitrary units, but these vaguely match a GE flatpanel X-ray detector in cm)
x_length = np.arange(-100,100.01,0.1)+0.05
z_distance = 600

x_arr, y_arr = np.meshgrid(x_length,x_length)
pixel_dist = (x_arr**2 + y_arr**2)**0.5
detector_theta = np.atan(pixel_dist/z_distance)/2

# assume a fwhm of 0.001 degrees and a gaussian distribution (not true, but good approximation)
det_a = np.zeros_like(detector_theta)
for i,ref in enumerate(ref_a):
    st_dev = np.abs(detector_theta-ref.theta)/0.0005
    mask = st_dev<5
    det_a[mask] = norm.pdf(st_dev[mask])*np.abs(ref.structure_factor)

det_b = np.zeros_like(detector_theta)
for i,ref in enumerate(ref_b):
    st_dev = np.abs(detector_theta-ref.theta)/0.0005
    mask = st_dev<5
    det_b[mask] = norm.pdf(st_dev[mask])*np.abs(ref.structure_factor)

# plot the two detectors
fig,ax = plt.subplots(1,2)

ax[0].imshow(det_a**0.5)
ax[0].set_title("triclinic reflections on a cubic lattice")
ax[1].imshow(det_b**0.5)
ax[1].set_title("triclinic reflections on a triclinic lattice")

So in simulations, you lose the difference between [100] vs [010] vs [001] for example, and instead get a single [100] ring with 3x intensity.

There are some other examples but they would be a lot of work to code up. The takeaway though is that lattices have symmetry.

[argerlt]

To maybe summarize my point, I think

They're valid safe and done add unnecessary assumptions. They're also the ones used in MTEX.

Is a solid reason to leave as-is. Maybe not crystallographically valid, but safer than the alternative.

[hakonanes]

Thank you, I appreciate the effort to stress the fact that orix is useful in diffraction work and that we should consider this.

You've convinced me that we can use different lattices for the default ones. The question is which. I'd prefer them to not be random. The options I see are:

1. Use valid but simple ones, like (a, b, c, alpha, beta, gamma) = (1, 2, 3, 90, 91, 92) for a triclinic lattice
2. Use fairly common ones as default. I think they should be with the highest symmetry possible for the system. For triclinic (-1), we could use Wollastonite (https://en.wikipedia.org/wiki/Triclinic_crystal_system). For cubic we could use Nickel. And so on.

Thoughts?

[argerlt]

I think option 2 makes more sense.

2. Use fairly common ones as default. I think they should be with the highest symmetry possible for the system. For triclinic (-1), we could use Wollastonite (https://en.wikipedia.org/wiki/Triclinic_crystal_system). For cubic we could use Nickel. And so on.

This side-steps.my worry of "making more arbitrary standards" as well. I like this option a lot.

Wollastonite and Nickel seem like safe bets. I will look for others, but if you have good suggestions, drop them here as well.

[hakonanes]

Of course, if we're returning a real phase, a valid question is then whether we should instead supply a small database of example phases via orix.data.phase("nickel") etc. We can do both:

1. orix.crystal_map.Phase.triclinic() returns a phase with a (1, 1, 1, 90, 90, 90) lattice and not atoms
2. orix.data.phase("wollastonite") returns a real triclinic phase

The idea is that the former can be used to quickly create an arbitrary lattice. The latter provides a databse of real phases we can use for diffraction examples in diffsims, pyxem, kikuchupy, and so on.

Sorry for the de-tour, here.

[argerlt]

Sorry for the de-tour, here.

Eh, it's a frustrating detour but an important one. Better to get it right now than fix it later.

Your last suggestion is a good plan. As you said:

The idea is that the former can be used to quickly create an arbitrary lattice. The latter provides a databse of real phases we can use for diffraction examples in diffsims, pyxem, kikuchupy, and so on.

Which seems like good compromise of dealing with my concerns while also not deviating from the MTEX defaults or forcing assumptions about phase.

Also, I didn't think about it until days later, but anecdotally, I've had considerable problems with HeXRD assuming Nickel as the default phase. The solution then was "force users to explicitly define everything", and this seems like a kinder but similar approach.