Structure factor calculation bug

[viljarjf]

Describe the bug I believe I found the bug making the two simulations different in #201. When a Phase-object is initialized, the structure is transformed to ensure x || a and z || c*. This changes structure.lattice, naturally, but according to the diffpy documentation, the xyz coordinates of the atoms are in the lattice coordinate system, regardless of them being transformed. Still, in https://github.com/CSSFrancis/diffsims/blob/use_phase_class/diffsims/utils/sim_utils.py#L281-L282, the positions of the atoms are transformed according to the lattice basis.

Since these lines are really old, I am probably wrong in thinking they should not be there.

Instead, there might be a bug when setting the transformed lattice in the Phase.structure-setter, where the positions of the atoms are still expressed as if they were in the old basis. This can be mitigated by copying the old cartesian positions and setting them again after transforming the basis.

With either approach (removing the two lines linked, or fixing the coordinates), the simulations are identical, and in both cases they more closely resemble the output of ASTAR diffgen. However, in the case of the transformed basis, the mat in the linked lines is NOT the identity matrix. Therefore, I'm leaning more towards this being a bug in the Phase.structure-setter.

As a sidenote: I believe transposing recbase in https://github.com/CSSFrancis/diffsims/blob/use_phase_class/diffsims/utils/sim_utils.py#L281 is incorrect, and seems to be the cause of the discrepancy between ASTAR and the old simulation. Since transposing is a new change, not present before this PR, and the "old" simulations we have been comparing have been made with the new transposing.

To Reproduce Might not be minimum working example, but at least an example.

import diffpy
import numpy as np
import matplotlib.pyplot as plt
from diffsims.libraries.structure_library import StructureLibrary
from diffsims.generators.diffraction_generator import DiffractionGenerator
from diffsims.generators.library_generator import DiffractionLibraryGenerator
from orix.crystal_map import Phase
from orix.quaternion import Rotation
from diffsims.generators.simulation_generator import SimulationGenerator

# Shared parameters
latt = diffpy.structure.lattice.Lattice(2.464, 2.464, 6.711, 90, 90, 120)
atoms = [diffpy.structure.atom.Atom(atype='C', xyz=[0.0, 0.0, 0.25], lattice=latt),
         diffpy.structure.atom.Atom(atype='C', xyz=[0.0, 0.0, 0.75], lattice=latt),
         diffpy.structure.atom.Atom(atype='C', xyz=[1/3, 2/3, 0.25], lattice=latt),
         diffpy.structure.atom.Atom(atype='C', xyz=[2/3, 1/3, 0.75], lattice=latt),]
structure_matrix = diffpy.structure.Structure(atoms=atoms, lattice=latt)
calibration = 0.0262
reciprocal_radius = 1.6768
with_direct_beam = False
max_excitation_error = 0.1
accelerating_voltage = 200
shape = (128, 128)
sigma = 1.4
generator_kwargs = {
    "accelerating_voltage": accelerating_voltage,
    "scattering_params": "lobato",
    "precession_angle": 0,
    "shape_factor_model": "lorentzian",
    "approximate_precession": True,
    "minimum_intensity": 1e-20,
}

# The euler angles are different, as orix.Phase enforces x||a, z||c*
euler_angles_old = np.array([[0, 90, 120]])
euler_angles_new = np.array([[0, 90, 90]])

# Old
struct_library = StructureLibrary(["Graphite"], [structure_matrix], [euler_angles_old])
diff_gen = DiffractionGenerator(**generator_kwargs)
lib_gen = DiffractionLibraryGenerator(diff_gen)
diff_lib = lib_gen.get_diffraction_library(struct_library,
                                           calibration=calibration,
                                           reciprocal_radius=reciprocal_radius,
                                           with_direct_beam=with_direct_beam,
                                           max_excitation_error=max_excitation_error,
                                           half_shape=shape[0]//2,
                                           )

# New
p = Phase("Graphite", point_group="6/mmm", structure=structure_matrix)
gen = SimulationGenerator(**generator_kwargs)
rot = Rotation.from_euler(euler_angles_new, degrees=True)
sim = gen.calculate_ed_data(phase=p, 
                            rotation=rot, 
                            reciprocal_radius=reciprocal_radius, 
                            max_excitation_error=max_excitation_error, 
                            with_direct_beam=with_direct_beam,
                            )

# Compare
old_data = diff_lib["Graphite"]["simulations"][0].get_diffraction_pattern(shape=shape, sigma=sigma)
new_data = sim.get_diffraction_pattern(shape=shape, sigma=sigma, calibration=calibration)

plt.figure(figsize=(9.5, 3))
plt.subplot(1, 3, 1)
plt.title("Difference |new - old|")
plt.imshow(np.abs(new_data - old_data))
plt.axis("off")
plt.colorbar()
plt.subplot(1, 3, 2)
plt.title("Old simulation, vmax=0.05")
plt.imshow(old_data, vmax=0.05)
plt.axis("off")
plt.colorbar()
plt.subplot(1, 3, 3)
plt.title("New simulation, vmax=0.05")
plt.imshow(new_data, vmax=0.05)
plt.axis("off")
plt.colorbar()
plt.tight_layout()

Intensities differ by around (0.05 - 0.004) $\approx$ 0.05, i.e. 10x

Expected behavior By removing .T in https://github.com/CSSFrancis/diffsims/blob/use_phase_class/diffsims/utils/sim_utils.py#L281, and by adding the following to https://github.com/pyxem/orix/blob/develop/orix/crystal_map/phase_list.py#L137-L139:

...
old_xyz_cartn = value.xyz_cartn
value.lattice.setLatBase(new_matrix)
value.xyz_cartn = old_xyz_cartn
...

the output of the example above becomes:

i.e. both practically identical simulations (intensity differences around 1e-11 can probably be attributed to the imprecise euler angle matrix calculations in the old simulations, which are much better now), and both more closely resemble the output of ASTAR (https://github.com/pyxem/diffsims/pull/201#issuecomment-1926638205):

Additional context Using https://github.com/CSSFrancis/diffsims/tree/use_phase_class

[CSSFrancis]

@viljarjf Wow thank you for digging into this! Let me go through this a little bit more tomorrow just to check. It seems like this is correct, but I just have to run some more tests.

[hakonanes]

I agree, this seems like a bug in orix. Looking forward to review a PR!

I'll have to test that PR change against kikuchipy as well.

[CSSFrancis]

@viljarjf Do you want to submit a PR to orix?

[viljarjf]

Yes, let me just write it up

[hakonanes]

I'd say we can close this as the bug was fixed by @viljarjf in orix in https://github.com/pyxem/orix/pull/469.