Skip to content

Commit

Permalink
shuffled parameters, implemented simplified lorentzian
Browse files Browse the repository at this point in the history
  • Loading branch information
@din14970
din14970 committed Nov 17, 2020
1 parent 9735db8 commit c957383
Show file tree
Hide file tree
Showing 3 changed files with 180 additions and 63 deletions.
155 changes: 99 additions & 56 deletions diffsims/generators/diffraction_generator.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,23 @@
get_intensities_params,
)
from diffsims.utils.fourier_transform import from_recip
from diffsims.utils.shape_factor_models import linear
from diffsims.utils.shape_factor_models import (
linear,
atanc,
lorentzian,
sinc,
sin2c,
lorentzian_precession,
)


_shape_factor_model_mapping = {
"linear": linear,
"atanc": atanc,
"sinc": sinc,
"sin2c": sin2c,
"lorentzian": lorentzian,
}


def _z_sphere_precession(phi, r_spot, wavelength, theta):
Expand Down Expand Up @@ -148,27 +164,69 @@ class DiffractionGenerator(object):
----------
accelerating_voltage : float
The accelerating voltage of the microscope in kV.
max_excitation_error : float
Removed in this version, defaults to None
debye_waller_factors : dict of str:value pairs
Maps element names to their temperature-dependent Debye-Waller factors.
scattering_params : str
"lobato" or "xtables"
minimum_intensity : float
Minimum intensity for a peak to be considered visible in the pattern
precession_angle : float
Angle about which the beam is precessed. Default is no precession.
approximate_precession : boolean
When using precession, whether to precisely calculate average
excitation errors and intensities or use an approximation.
shape_factor_model : function or string
A function that takes excitation_error and
`max_excitation_error` (and potentially kwargs) and returns
an intensity scaling factor. If None defaults to
`shape_factor_models.lorentzian`. A number of pre-programmed functions
are available via strings.
kwargs
Keyword arguments passed to `shape_factor_model`.
Notes
-----
* A full calculation is much slower and is not recommended for calculating
a diffraction library for precession diffraction patterns.
* When using precession and approximate_precession=True, the shape factor
model defaults to lorentzian; shape_factor_model is ignored. Only with
approximate_precession=False the custom shape_factor_model is used.
"""

def __init__(
self,
accelerating_voltage,
max_excitation_error=None,
debye_waller_factors={},
scattering_params="lobato",
precession_angle=None,
shape_factor_model=None,
approximate_precession=True,
minimum_intensity=1e-20,
**kwargs,
):
if max_excitation_error is not None:
print(
"This class changed in v0.3 and no longer takes a maximum_excitation_error"
)
self.wavelength = get_electron_wavelength(accelerating_voltage)
self.debye_waller_factors = debye_waller_factors
self.precession_angle = precession_angle
if self.precession_angle is None:
self.precession_angle = 0
if self.precession_angle < 0:
raise ValueError("The precession angle must be >= 0")
self.approximate_precession = approximate_precession
if isinstance(shape_factor_model, str):
if shape_factor_model in _shape_factor_model_mapping.keys():
self.shape_factor_model = _shape_factor_model_mapping[shape_factor_model]
else:
raise NotImplementedError(
f"{shape_factor_model} is not a recognized shape factor "
f"model, choose from: {_shape_factor_model_mapping.keys()} "
f"or provide your own function."
)
elif callable(shape_factor_model):
self.shape_factor_model = shape_factor_model
elif shape_factor_model is None:
self.shape_factor_model = lorentzian
else:
raise TypeError(
"shape_factor_model must be a recognized string "
"or a callable object.")
self.minimum_intensity = minimum_intensity
self.shape_factor_kwargs = kwargs
if scattering_params in ["lobato", "xtables"]:
self.scattering_params = scattering_params
else:
Expand All @@ -183,13 +241,9 @@ def calculate_ed_data(
structure,
reciprocal_radius,
rotation=(0, 0, 0),
shape_factor_model=None,
max_excitation_error=1e-2,
with_direct_beam=True,
minimum_intensity=1e-20,
precession_angle=0,
full_calculation=False,
**kwargs
max_excitation_error=1e-2,
debye_waller_factors={},
):
"""Calculates the Electron Diffraction data for a structure.
Expand All @@ -206,26 +260,14 @@ def calculate_ed_data(
rotation : tuple
Euler angles, in degrees, in the rzxz convention. Default is
(0, 0, 0) which aligns 'z' with the electron beam.
shape_factor_model : function or None
A function that takes excitation_error and
`max_excitation_error` (and potentially kwargs) and returns
an intensity scaling factor. If None defaults to
`shape_factor_models.atanc`.
max_excitation_error : float
The extinction distance for reflections, in reciprocal
Angstroms.
with_direct_beam : bool
If True, the direct beam is included in the simulated
diffraction pattern. If False, it is not.
minimum_intensity : float
Minimum intensity for a peak to be considered in the pattern
precession_angle : float
Angle about which the beam is precessed
full_calculation : boolean
When using precession, whether to precisely calculate average
excitation errors for filtering diffraction spots.
kwargs
Keyword arguments passed to `shape_factor_model`.
max_excitation_error : float
The extinction distance for reflections, in reciprocal
Angstroms. Roughly equal to 1/thickness.
debye_waller_factors : dict of str:value pairs
Maps element names to their temperature-dependent Debye-Waller factors.
Returns
-------
Expand Down Expand Up @@ -257,20 +299,14 @@ def calculate_ed_data(
r_spot = np.sqrt(np.sum(np.square(cartesian_coordinates[:, :2]), axis=1))
z_spot = cartesian_coordinates[:, 2]

if precession_angle is None:
precession_angle = 0

if shape_factor_model is None:
shape_factor_model = linear

if precession_angle > 0 and full_calculation:
if self.precession_angle > 0 and not self.approximate_precession:
# We find the average excitation error - this step can be
# quite expensive
excitation_error = _average_excitation_error_precession(
z_spot,
r_spot,
wavelength,
precession_angle,
self.precession_angle,
)
else:
z_sphere = -np.sqrt(r_sphere ** 2 - r_spot ** 2) + r_sphere
Expand All @@ -279,25 +315,32 @@ def calculate_ed_data(
intersection = np.abs(excitation_error) < max_excitation_error
intersection_coordinates = cartesian_coordinates[intersection]
excitation_error = excitation_error[intersection]
r_spot = r_spot[intersection]
g_indices = spot_indices[intersection]
g_hkls = spot_distances[intersection]
# take into consideration rel-rods
if precession_angle > 0:
if self.precession_angle > 0 and not self.approximate_precession:
shape_factor = _shape_factor_precession(
intersection_coordinates[:, 2],
r_spot[intersection],
r_spot,
wavelength,
precession_angle,
shape_factor_model,
max_excitation_error,
**kwargs,
self.precession_angle,
self.shape_factor_model,
self.max_excitation_error,
**self.shape_factor_kwargs,
)
elif self.precession_angle > 0 and self.approximate_precession:
shape_factor = lorentzian_precession(
excitation_error,
max_excitation_error,
r_spot,
self.precession_angle,
)
elif precession_angle == 0:
shape_factor = shape_factor_model(
excitation_error, max_excitation_error, **kwargs
)
else:
raise ValueError("The precession angle must be >= 0")
shape_factor = self.shape_factor_model(
excitation_error, max_excitation_error,
**self.shape_factor_kwargs
)

# Calculate diffracted intensities based on a kinematical model.
intensities = get_kinematical_intensities(
Expand All @@ -306,11 +349,11 @@ def calculate_ed_data(
g_hkls,
prefactor=shape_factor,
scattering_params=self.scattering_params,
debye_waller_factors=self.debye_waller_factors,
debye_waller_factors=debye_waller_factors,
)

# Threshold peaks included in simulation based on minimum intensity.
peak_mask = intensities > minimum_intensity
peak_mask = intensities > self.minimum_intensity
intensities = intensities[peak_mask]
intersection_coordinates = intersection_coordinates[peak_mask]
g_indices = g_indices[peak_mask]
Expand Down
20 changes: 14 additions & 6 deletions diffsims/generators/library_generator.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,8 @@ def get_diffraction_library(
reciprocal_radius,
half_shape,
with_direct_beam=True,
max_excitation_error=1e-2,
debye_waller_factors={},
):
"""Calculates a dictionary of diffraction data for a library of crystal
structures and orientations.
Expand All @@ -74,6 +76,11 @@ def get_diffraction_library(
The half shape of the target patterns, for 144x144 use (72,72) etc
with_direct_beam : bool
Include the direct beam in the library.
max_excitation_error : float
The extinction distance for reflections, in reciprocal
Angstroms.
debye_waller_factors : dict of str:value pairs
Maps element names to their temperature-dependent Debye-Waller factors.
Returns
-------
Expand Down Expand Up @@ -103,6 +110,8 @@ def get_diffraction_library(
reciprocal_radius,
rotation=orientation,
with_direct_beam=with_direct_beam,
max_excitation_error=max_excitation_error,
debye_waller_factors=debye_waller_factors
)

# Calibrate simulation
Expand All @@ -117,12 +126,11 @@ def get_diffraction_library(
intensities[i] = simulation.intensities

diffraction_library[phase_name] = {
"simulations": simulations,
"orientations": orientations,
"pixel_coords": pixel_coords,
"intensities": intensities,
}

"simulations": simulations,
"orientations": orientations,
"pixel_coords": pixel_coords,
"intensities": intensities,
}
# Pass attributes to diffraction library from structure library.
diffraction_library.identifiers = structure_library.identifiers
diffraction_library.structures = structure_library.structures
Expand Down
68 changes: 67 additions & 1 deletion diffsims/utils/shape_factor_models.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,7 @@ def linear(excitation_error, max_excitation_error):
intensities : array-like or float
"""

return 1 - excitation_error / max_excitation_error
return 1 - np.abs(excitation_error) / max_excitation_error


def sinc(excitation_error, max_excitation_error, minima_number=5):
Expand Down Expand Up @@ -130,3 +130,69 @@ def atanc(excitation_error, max_excitation_error):
np.arctan(fac*excitation_error)/(fac*excitation_error),
nan=1,
)


def lorentzian(excitation_error, max_excitation_error):
"""
Lorentzian intensity profile that should approximate
the two-beam rocking curve. After [1].
Parameters
----------
excitation_error : array-like or float
The distance (reciprocal) from a reflection to the Ewald sphere
max_excitation_error : float
The distance at which a reflection becomes extinct
Returns
-------
intensity_factor : array-like or float
Vector representing the rel-rod factor for each reflection
Notes
-----
[1] L. Palatinus, P. Brázda, M. Jelínek, J. Hrdá, G. Steciuk, M. Klementová, Specifics of the data processing of precession electron diffraction tomography data and their implementation in the program PETS2.0, Acta Crystallogr. Sect. B Struct. Sci. Cryst. Eng. Mater. 75 (2019) 512–522. doi:10.1107/S2052520619007534.
"""
# in the paper, sigma = pi*thickness.
# We assume thickness = 1/max_exitation_error
sigma = np.pi/max_excitation_error
fac = sigma/(np.pi*(sigma**2*excitation_error**2+1))
return fac


def lorentzian_precession(excitation_error, max_excitation_error,
r_spot, precession_angle):
"""
Intensity profile factor for a precessed beam assuming a Lorentzian
intensity profile for the un-precessed beam. After [1].
Parameters
----------
excitation_error : array-like or float
The distance (reciprocal) from a reflection to the Ewald sphere
max_excitation_error : float
The distance at which a reflection becomes extinct
r_spot : array-like or float
The distance (reciprocal) from each reflection to the origin
precession_angle : float
The beam precession angle in degrees; the angle the beam makes
with the optical axis.
Returns
-------
intensity_factor : array-like or float
Vector representing the rel-rod factor for each reflection
Notes
-----
[1] L. Palatinus, P. Brázda, M. Jelínek, J. Hrdá, G. Steciuk, M. Klementová, Specifics of the data processing of precession electron diffraction tomography data and their implementation in the program PETS2.0, Acta Crystallogr. Sect. B Struct. Sci. Cryst. Eng. Mater. 75 (2019) 512–522. doi:10.1107/S2052520619007534.
"""
sigma = np.pi/max_excitation_error
u = sigma**2*(r_spot**2*precession_angle**2 - excitation_error**2)+1
z = np.sqrt(u**2 + 4*sigma**2 + excitation_error**2)
fac = (sigma/np.pi)*np.sqrt(2*(u+z)/z**2)
return fac

0 comments on commit c957383

Please sign in to comment.