-
Notifications
You must be signed in to change notification settings - Fork 27
Description
Describe the bug
The new simulations with a Phase does not (necessarily) take symmetry into account when calculating intensities. Only the explicitly included atoms in the underlying structure are part of the structure factor calculations.
To Reproduce
Steps to reproduce the behavior:
from numpy import allclose
from matplotlib import pyplot as plt
from diffpy.structure import Structure, Atom, Lattice
from orix.crystal_map import Phase
from orix.quaternion import Rotation
from diffsims.generators.simulation_generator import SimulationGenerator
l = Lattice(10, 10, 10, 90, 90, 90)
a = [Atom("Au", (0, 0, 0))]
s = Structure(a, l)
P1 = Phase(structure=s, space_group=1)
Fd3m = Phase(structure=s, space_group=227)
# Rotate slightly off-zone to see extinct reflections
rot = Rotation.from_euler((0, 0, 1), degrees=True)
gen = SimulationGenerator()
P1_sim = gen.calculate_diffraction2d(P1, rotation=rot)
Fd3m_sim = gen.calculate_diffraction2d(Fd3m, rotation=rot)
# Use ReciprocalLatticeVector.sanitize_phase
Fd3m_sim.coordinates.sanitise_phase()
Fd3m_expanded = Fd3m_sim.coordinates.phase
Fd3m_sim_expanded = gen.calculate_diffraction2d(Fd3m_expanded, rotation=rot)
# Intensities should be different due to the added symmetry,
# even when the symmetrically equivalent atoms are missing from the structure
assert not allclose(P1_sim.coordinates.intensity, Fd3m_sim_expanded.coordinates.intensity) # Passes
# assert not allclose(P1_sim.coordinates.intensity, Fd3m_sim.coordinates.intensity) # Fails
# Plot patterns
plt.figure()
plt.subplot(1, 3, 1)
plt.title("P1")
P1_sim.plot(ax=plt.gca())
plt.subplot(1, 3, 2)
plt.title("Fd3m default")
Fd3m_sim.plot(ax=plt.gca())
plt.subplot(1, 3, 3)
plt.title("Fd3m expanded")
Fd3m_sim_expanded.plot(ax=plt.gca())Currently, intensities are calculated as follows:
- Find intersecting vectors for a given reflection
- Calculate the structure factor of each intersected reflection
- Square the structure factor, and scale it according to the excitation error
I suggest instead to calculate the structure factors with ReciprocalLatticeVector.calculate_structure_factor(), after calling ReciprocalLatticeVector.sanitize_phase(). This reduces code duplication, and lets us avoid a lot of re-computing. The loop over rotations then only needs to find the excitation error, and scale the squared structure factor using that.
Expected behavior
Expand the unit cell to ensure structure factors are calculated correctly, and then scale using the excitation error
