Skip to content

Commit

Permalink
Merge 66da664 into 8482796
Browse files Browse the repository at this point in the history
  • Loading branch information
@isabelwood100
isabelwood100 committed Aug 27, 2020
2 parents 8482796 + 66da664 commit 9af8eac
Show file tree
Hide file tree
Showing 4 changed files with 192 additions and 257 deletions.
116 changes: 48 additions & 68 deletions diffsims/generators/diffraction_generator.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,7 @@
get_points_in_sphere,
get_vectorized_list_for_atomic_scattering_factors,
is_lattice_hexagonal,
get_intensities_params,
)
from diffsims.utils.atomic_diffraction_generator_support.fourier_transform import (
from_recip,
Expand Down Expand Up @@ -127,7 +128,7 @@ def calculate_ed_data(
# Obtain crystallographic reciprocal lattice points within `max_r` and
# g-vector magnitudes for intensity calculations.
recip_latt = latt.reciprocal()
spot_indices, cartesian_coordinates, spot_distances = get_points_in_sphere(
spot_indices, spot_coords, spot_distances = get_points_in_sphere(
recip_latt, reciprocal_radius
)

Expand All @@ -137,41 +138,43 @@ def calculate_ed_data(
np.deg2rad(rotation[2]),
)
R = euler2mat(ai, aj, ak, axes="rzxz")
cartesian_coordinates = np.matmul(R, cartesian_coordinates.T).T
spot_coords = np.matmul(R, spot_coords.T).T

# Identify points intersecting the Ewald sphere within maximum
# excitation error and store the magnitude of their excitation error.
r_sphere = 1 / wavelength
r_spot = np.sqrt(np.sum(np.square(cartesian_coordinates[:, :2]), axis=1))
r_spot = np.sqrt(np.sum(np.square(spot_coords[:, :2]), axis=1))
z_sphere = -np.sqrt(r_sphere ** 2 - r_spot ** 2) + r_sphere
proximity = np.absolute(z_sphere - cartesian_coordinates[:, 2])
intersection = proximity < max_excitation_error
excitation_error = np.absolute(z_sphere - spot_coords[:, 2])
intersection = excitation_error < max_excitation_error
# Mask parameters corresponding to excited reflections.
intersection_coordinates = cartesian_coordinates[intersection]
intersection_indices = spot_indices[intersection]
proximity = proximity[intersection]
intersection_coordinates = spot_coords[intersection]
g_indices = spot_indices[intersection]
excitation_error = excitation_error[intersection]
g_hkls = spot_distances[intersection]
multiplicites = np.ones_like(g_hkls)

# Calculate diffracted intensities based on a kinematical model.
intensities = get_kinematical_intensities(
structure,
intersection_indices,
g_indices,
g_hkls,
proximity,
max_excitation_error,
debye_waller_factors,
multiplicites,
scattering_params,
excitation_error,
max_excitation_error,
)

# Threshold peaks included in simulation based on minimum intensity.
peak_mask = intensities > 1e-20
intensities = intensities[peak_mask]
intersection_coordinates = intersection_coordinates[peak_mask]
intersection_indices = intersection_indices[peak_mask]
g_indices = g_indices[peak_mask]

return DiffractionSimulation(
coordinates=intersection_coordinates,
indices=intersection_indices,
indices=g_indices,
intensities=intensities,
with_direct_beam=with_direct_beam,
)
Expand All @@ -185,6 +188,7 @@ def calculate_profile_data(
):
"""
Calculates a one dimensional diffraction profile for a structure.
For use with the get_plot_x in DiffractionSimulation.
Parameters
----------
Expand Down Expand Up @@ -213,79 +217,55 @@ def calculate_profile_data(
latt = structure.lattice
is_hex = is_lattice_hexagonal(latt)

(
coeffs,
fcoords,
occus,
dwfactors,
) = get_vectorized_list_for_atomic_scattering_factors(
structure, {}, scattering_params=scattering_params
)

# Obtain crystallographic reciprocal lattice points within range
recip_latt = latt.reciprocal()
spot_indices, _, spot_distances = get_points_in_sphere(
recip_latt, reciprocal_radius
)

peaks = {}
mask = np.logical_not((np.any(spot_indices, axis=1) == 0))

for hkl, g_hkl in zip(spot_indices[mask], spot_distances[mask]):
# Force miller indices to be integers.
hkl = [int(round(i)) for i in hkl]

d_hkl = 1 / g_hkl

# Bragg condition
# theta = asin(wavelength * g_hkl / 2)

# s = sin(theta) / wavelength = 1 / 2d = |ghkl| / 2 (d =
# 1/|ghkl|)
s = g_hkl / 2

# Store s^2 since we are using it a few times.
s2 = s ** 2

# Vectorized computation of g.r for all fractional coords and
# hkl.
g_dot_r = np.dot(fcoords, np.transpose([hkl])).T[0]

# Highly vectorized computation of atomic scattering factors.
fs = np.sum(coeffs[:, :, 0] * np.exp(-coeffs[:, :, 1] * s2), axis=1)

dw_correction = np.exp(-dwfactors * s2)

# Structure factor = sum of atomic scattering factors (with
# position factor exp(2j * pi * g.r and occupancies).
# Vectorized computation.
f_hkl = np.sum(fs * occus * np.exp(2j * np.pi * g_dot_r) * dw_correction)
##spot_indicies is a numpy.array of the hkls allowd in the recip radius
unique_hkls, multiplicites, g_hkls = get_intensities_params(
recip_latt, reciprocal_radius
)
g_indices = unique_hkls
debye_waller_factors = self.debye_waller_factors
excitation_error = None
max_excitation_error = None
g_hkls_array = np.asarray(g_hkls)

# Intensity for hkl is modulus square of structure factor.
i_hkl = (f_hkl * f_hkl.conjugate()).real
i_hkl = get_kinematical_intensities(
structure,
g_indices,
g_hkls_array,
debye_waller_factors,
multiplicites,
scattering_params,
excitation_error=None,
max_excitation_error=None,
)

# two_theta = degrees(2 * theta)
if is_hex:
# Use Miller-Bravais indices for hexagonal lattices.
hkl = (hkl[0], hkl[1], -hkl[0] - hkl[1], hkl[2])

if is_hex:
# Use Miller-Bravais indices for hexagonal lattices.
hkl = (hkl[0], hkl[1], -hkl[0] - hkl[1], hkl[2])
hkls_labels = ["".join([str(int(x)) for x in xs]) for xs in unique_hkls]

peaks[g_hkl] = [i_hkl, [tuple(hkl)], d_hkl]
peaks = {}
for l, i, g in zip(hkls_labels, i_hkl, g_hkls):
peaks[l] = [i, g]

# Scale intensities so that the max intensity is 100.

max_intensity = max([v[0] for v in peaks.values()])
x = []
y = []
hkls = []
d_hkls = []
for k in sorted(peaks.keys()):
for k in peaks.keys():
v = peaks[k]
fam = get_unique_families(v[1])
if v[0] / max_intensity * 100 > minimum_intensity:
x.append(k)
if v[0] / max_intensity * 100 > minimum_intensity and (k != "000"):
x.append(v[1])
y.append(v[0])
hkls.append(fam)
d_hkls.append(v[2])
hkls.append(k)

y = y / max(y) * 100

Expand Down
48 changes: 43 additions & 5 deletions diffsims/sims/diffraction_simulation.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,8 @@
import numpy as np
from scipy import ndimage as ndi

from diffsims.utils.sim_utils import get_unique_families


class DiffractionSimulation:
"""Holds the result of a kinematic diffraction pattern simulation.
Expand Down Expand Up @@ -148,16 +150,17 @@ def get_diffraction_pattern(self, size=512, sigma=10):
the order of 0.5nm and a the default size and sigma are used.
"""
side_length = np.min(np.multiply((size / 2), self.calibration))
mask_for_sides = np.all((np.abs(self.coordinates[:, 0:2]) < side_length), axis=1)

mask_for_sides = np.all(
(np.abs(self.coordinates[:, 0:2]) < side_length), axis=1
)

spot_coords = np.add(
self.calibrated_coordinates[mask_for_sides], size / 2
).astype(int)
spot_intens = self.intensities[mask_for_sides]
pattern = np.zeros([size, size])
#checks that we have some spots
if spot_intens.shape[0]==0:
# checks that we have some spots
if spot_intens.shape[0] == 0:
return pattern
else:
pattern[spot_coords[:, 0], spot_coords[:, 1]] = spot_intens
Expand All @@ -174,6 +177,8 @@ class ProfileSimulation:
Magnitudes of scattering vectors.
intensities : array-like, shape [n_peaks, 1]
The kinematic intensity of the diffraction peaks.
##change this to list qq
hkls: [{(h, k, l): mult}] {(h, k, l): mult} is a dict of Miller
indices for all diffracted lattice facets contributing to each
intensity.
Expand All @@ -200,10 +205,13 @@ def get_plot(self, g_max, annotate_peaks=True, with_labels=True, fontsize=12):
fontsize : integer
Fontsize for peak labels.
"""

label_hkl = get_unique_families(self.hkls)

ax = plt.gca()
for g, i, hkls in zip(self.magnitudes, self.intensities, self.hkls):
if g <= g_max:
label = ", ".join([str(hkl) for hkl in hkls.keys()])
label = ", ".join([str(hkl) for hkl in label_hkl.keys()])
ax.plot([g, g], [0, i], color="k", linewidth=3, label=label)
if annotate_peaks:
ax.annotate(label, xy=[g, i], xytext=[g, i], fontsize=fontsize)
Expand All @@ -213,3 +221,33 @@ def get_plot(self, g_max, annotate_peaks=True, with_labels=True, fontsize=12):
ax.set_ylabel("Intensities (scaled)")

return plt

def get_plot_x(self, g_max, annotate_peaks=True, with_labels=True, fontsize=12):

"""Plots the diffraction profile simulation for the
calculate_profile_data method in DiffractionGenerator.
Parameters
----------
g_max : float
Maximum g-vector magnitude to plot.
annotate_peaks : boolean
If True, peaks are annotaed with hkl information.
with_labels : boolean
If True, xlabels and ylabels are added to the plot.
fontsize : integer
Fontsize for peak labels.
"""

ax = plt.gca()
for g, i, hkls in zip(self.magnitudes, self.intensities, self.hkls):
label = hkls
ax.plot([g, g], [0, i], color="k", linewidth=3, label=label)
if annotate_peaks:
ax.annotate(label, xy=[g, i], xytext=[g, i], fontsize=fontsize)

if with_labels:
ax.set_xlabel("A ($^{-1}$)")
ax.set_ylabel("Intensities (scaled)")

return plt
Loading

0 comments on commit 9af8eac

Please sign in to comment.