From 9b357657f211fdd4e458a62c9274e09e72464777 Mon Sep 17 00:00:00 2001 From: Colin Date: Wed, 18 Oct 2023 11:54:26 -0700 Subject: [PATCH 1/9] Adding occupancy to crystal class --- py4DSTEM/process/diffraction/crystal.py | 7 +++ py4DSTEM/process/diffraction/crystal_viz.py | 49 ++++++++++++++++----- 2 files changed, 45 insertions(+), 11 deletions(-) diff --git a/py4DSTEM/process/diffraction/crystal.py b/py4DSTEM/process/diffraction/crystal.py index 30bb061f8..af949efbe 100644 --- a/py4DSTEM/process/diffraction/crystal.py +++ b/py4DSTEM/process/diffraction/crystal.py @@ -74,6 +74,7 @@ def __init__( positions, numbers, cell, + occupancy = None, ): """ Args: @@ -139,6 +140,12 @@ def __init__( else: raise Exception("Cell cannot contain " + np.size(cell) + " entries") + # occupancy + if occupancy is not None: + self.occupancy = np.array(occupancy) + else: + self.occupancy = None + # pymatgen flag if "pymatgen" in sys.modules: self.pymatgen_available = True diff --git a/py4DSTEM/process/diffraction/crystal_viz.py b/py4DSTEM/process/diffraction/crystal_viz.py index 9f9336155..3bd050f5a 100644 --- a/py4DSTEM/process/diffraction/crystal_viz.py +++ b/py4DSTEM/process/diffraction/crystal_viz.py @@ -91,18 +91,26 @@ def plot_structure( # Fractional atomic coordinates pos = self.positions + occ = self.occupancy + # x tile sub = pos[:, 0] < tol_distance pos = np.vstack([pos, pos[sub, :] + np.array([1, 0, 0])]) ID = np.hstack([ID, ID[sub]]) + if occ is not None: + occ = np.hstack([occ, occ[sub]]) # y tile sub = pos[:, 1] < tol_distance pos = np.vstack([pos, pos[sub, :] + np.array([0, 1, 0])]) ID = np.hstack([ID, ID[sub]]) + if occ is not None: + occ = np.hstack([occ, occ[sub]]) # z tile sub = pos[:, 2] < tol_distance pos = np.vstack([pos, pos[sub, :] + np.array([0, 0, 1])]) ID = np.hstack([ID, ID[sub]]) + if occ is not None: + occ = np.hstack([occ, occ[sub]]) # Cartesian atomic positions xyz = pos @ self.lat_real @@ -141,17 +149,36 @@ def plot_structure( # atoms ID_all = np.unique(ID) - for ID_plot in ID_all: - sub = ID == ID_plot - ax.scatter( - xs=xyz[sub, 1], # + d[0], - ys=xyz[sub, 0], # + d[1], - zs=xyz[sub, 2], # + d[2], - s=size_marker, - linewidth=2, - facecolors=atomic_colors(ID_plot), - edgecolor=[0, 0, 0], - ) + if occ is None: + for ID_plot in ID_all: + sub = ID == ID_plot + ax.scatter( + xs=xyz[sub, 1], # + d[0], + ys=xyz[sub, 0], # + d[1], + zs=xyz[sub, 2], # + d[2], + s=size_marker, + linewidth=2, + facecolors=atomic_colors(ID_plot), + edgecolor=[0, 0, 0], + ) + else: + for ID_plot in ID_all: + sub = ID == ID_plot + ax.scatter( + xs=xyz[sub, 1], # + d[0], + ys=xyz[sub, 0], # + d[1], + zs=xyz[sub, 2], # + d[2], + s=size_marker, + linewidth=2, + facecolors='none', + edgecolor=[0, 0, 0], + ) + # poly = PolyCollection( + # verts, + # facecolors=['r', 'g', 'b', 'y'], + # alpha = 0.6) + # ax.add_collection3d(poly, zs=zs, zdir='y') + # plot limit if plot_limit is None: From efce6cddb1f930460ff0fcb5e9600d46b383ab30 Mon Sep 17 00:00:00 2001 From: Colin Date: Wed, 18 Oct 2023 13:14:48 -0700 Subject: [PATCH 2/9] updating occupancy plot --- py4DSTEM/process/diffraction/crystal_viz.py | 87 ++++++++++++++++++--- 1 file changed, 76 insertions(+), 11 deletions(-) diff --git a/py4DSTEM/process/diffraction/crystal_viz.py b/py4DSTEM/process/diffraction/crystal_viz.py index 3bd050f5a..ced51ea0e 100644 --- a/py4DSTEM/process/diffraction/crystal_viz.py +++ b/py4DSTEM/process/diffraction/crystal_viz.py @@ -3,6 +3,8 @@ from matplotlib.axes import Axes import matplotlib.tri as mtri from mpl_toolkits.mplot3d import Axes3D, art3d +from mpl_toolkits.mplot3d.art3d import Poly3DCollection + from scipy.signal import medfilt from scipy.ndimage import gaussian_filter from scipy.ndimage.morphology import distance_transform_edt @@ -162,17 +164,80 @@ def plot_structure( edgecolor=[0, 0, 0], ) else: - for ID_plot in ID_all: - sub = ID == ID_plot - ax.scatter( - xs=xyz[sub, 1], # + d[0], - ys=xyz[sub, 0], # + d[1], - zs=xyz[sub, 2], # + d[2], - s=size_marker, - linewidth=2, - facecolors='none', - edgecolor=[0, 0, 0], - ) + # init + tol = 1e-4 + num_seg = 180 + radius = 0.5 + zp = np.zeros(num_seg+1) + + mark = np.ones(xyz.shape[0],dtype='bool') + for a0 in range(xyz.shape[0]): + if mark[a0]: + xyz_plot = xyz[a0,:] + inds = np.argwhere(np.sum((xyz - xyz_plot)**2,axis=1) < tol) + test = np.sum((xyz - xyz_plot)**2,axis=1) + mark[inds] = False + + ID_plot = ID[inds] + occ_plot = occ[inds] + + if np.sum(occ_plot) < 1.0: + occ_plot = np.append(occ_plot, 1 - np.sum(occ_plot)) + ID_plot = np.append(ID_plot, -1) + + # Plot site as series of filled arcs + theta0 = 0 + for a1 in range(occ_plot.size): + theta1 = theta0 + occ_plot[a1] * 2.0 * np.pi + theta = np.linspace(theta0,theta1,num_seg+1) + + xp = np.cos(theta) * radius + yp = np.sin(theta) * radius + xyz_rot = np.vstack((xp,yp,zp)) + xyz_rot = np.append(xyz_rot,np.array((0,0,0))[:,None], axis = 1) + + # Rotate towards camera + xyz_rot = orientation_matrix @ xyz_rot + + # add to plot + verts = [list(zip( + xyz_rot[1,:] + xyz_plot[1], + xyz_rot[0,:] + xyz_plot[0], + xyz_rot[2,:] + xyz_plot[2], + ))] + # ax.add_collection3d( + # Poly3DCollection( + # verts + # ) + # ) + collection = Poly3DCollection( + verts, + linewidths = 2.0, + alpha = 1.0, + edgecolors = 'k', + ) + face_color = [0.5, 0.5, 1] # alternative: matplotlib.colors.rgb2hex([0.5, 0.5, 1]) + if ID_plot[a1] == -1: + collection.set_facecolor((1.0,1.0,1.0)) + else: + collection.set_facecolor(atomic_colors(ID_plot[a1])) + ax.add_collection3d(collection) + + # update start point + if a1 < occ_plot.size: + theta0 = theta1 + + # for ID_plot in ID_all: + # sub = ID == ID_plot + # ax.scatter( + # xs=xyz[sub, 1], # + d[0], + # ys=xyz[sub, 0], # + d[1], + # zs=xyz[sub, 2], # + d[2], + # s=size_marker, + # linewidth=2, + # facecolors='none', + # edgecolor=[0, 0, 0], + # ) # poly = PolyCollection( # verts, # facecolors=['r', 'g', 'b', 'y'], From 6911bd3ddbbec59b6d3edf0d5179bf6b93c2ddfe Mon Sep 17 00:00:00 2001 From: Colin Date: Wed, 18 Oct 2023 14:51:20 -0700 Subject: [PATCH 3/9] Finishing kinematical occupancy --- py4DSTEM/process/diffraction/crystal.py | 18 +++-- py4DSTEM/process/diffraction/crystal_viz.py | 89 +++++++++++---------- 2 files changed, 59 insertions(+), 48 deletions(-) diff --git a/py4DSTEM/process/diffraction/crystal.py b/py4DSTEM/process/diffraction/crystal.py index af949efbe..1e83aea4e 100644 --- a/py4DSTEM/process/diffraction/crystal.py +++ b/py4DSTEM/process/diffraction/crystal.py @@ -590,12 +590,20 @@ def calculate_structure_factors( # Calculate structure factors self.struct_factors = np.zeros(np.size(self.g_vec_leng, 0), dtype="complex64") for a0 in range(self.positions.shape[0]): - self.struct_factors += f_all[:, a0] * np.exp( - (2j * np.pi) - * np.sum( - self.hkl * np.expand_dims(self.positions[a0, :], axis=1), axis=0 + if self.occupancy is None: + self.struct_factors += f_all[:, a0] * np.exp( + (2j * np.pi) + * np.sum( + self.hkl * np.expand_dims(self.positions[a0, :], axis=1), axis=0 + ) + ) + else: + self.struct_factors += f_all[:, a0] * np.exp( + (2j * np.pi * self.occupancy[a0]) + * np.sum( + self.hkl * np.expand_dims(self.positions[a0, :], axis=1), axis=0 + ) ) - ) # Divide by unit cell volume unit_cell_volume = np.abs(np.linalg.det(self.lat_real)) diff --git a/py4DSTEM/process/diffraction/crystal_viz.py b/py4DSTEM/process/diffraction/crystal_viz.py index ced51ea0e..89fd70bc6 100644 --- a/py4DSTEM/process/diffraction/crystal_viz.py +++ b/py4DSTEM/process/diffraction/crystal_viz.py @@ -167,7 +167,7 @@ def plot_structure( # init tol = 1e-4 num_seg = 180 - radius = 0.5 + radius = 0.7 zp = np.zeros(num_seg+1) mark = np.ones(xyz.shape[0],dtype='bool') @@ -175,57 +175,60 @@ def plot_structure( if mark[a0]: xyz_plot = xyz[a0,:] inds = np.argwhere(np.sum((xyz - xyz_plot)**2,axis=1) < tol) - test = np.sum((xyz - xyz_plot)**2,axis=1) + occ_plot = occ[inds] mark[inds] = False - ID_plot = ID[inds] - occ_plot = occ[inds] if np.sum(occ_plot) < 1.0: occ_plot = np.append(occ_plot, 1 - np.sum(occ_plot)) ID_plot = np.append(ID_plot, -1) + else: + occ_plot = occ_plot[0] + ID_plot = ID_plot[0] - # Plot site as series of filled arcs - theta0 = 0 - for a1 in range(occ_plot.size): - theta1 = theta0 + occ_plot[a1] * 2.0 * np.pi - theta = np.linspace(theta0,theta1,num_seg+1) - xp = np.cos(theta) * radius - yp = np.sin(theta) * radius - xyz_rot = np.vstack((xp,yp,zp)) - xyz_rot = np.append(xyz_rot,np.array((0,0,0))[:,None], axis = 1) + # Plot site as series of filled arcs + theta0 = 0 + for a1 in range(occ_plot.shape[0]): + theta1 = theta0 + occ_plot[a1] * 2.0 * np.pi + theta = np.linspace(theta0,theta1,num_seg+1) + xp = np.cos(theta) * radius + yp = np.sin(theta) * radius + - # Rotate towards camera - xyz_rot = orientation_matrix @ xyz_rot - - # add to plot - verts = [list(zip( - xyz_rot[1,:] + xyz_plot[1], - xyz_rot[0,:] + xyz_plot[0], - xyz_rot[2,:] + xyz_plot[2], - ))] - # ax.add_collection3d( - # Poly3DCollection( - # verts - # ) - # ) - collection = Poly3DCollection( - verts, - linewidths = 2.0, - alpha = 1.0, - edgecolors = 'k', - ) - face_color = [0.5, 0.5, 1] # alternative: matplotlib.colors.rgb2hex([0.5, 0.5, 1]) - if ID_plot[a1] == -1: - collection.set_facecolor((1.0,1.0,1.0)) - else: - collection.set_facecolor(atomic_colors(ID_plot[a1])) - ax.add_collection3d(collection) - - # update start point - if a1 < occ_plot.size: - theta0 = theta1 + # Rotate towards camera + xyz_rot = np.vstack((xp.ravel(),yp.ravel(),zp.ravel())) + if occ_plot[a1] < 1.0: + xyz_rot = np.append(xyz_rot,np.array((0,0,0))[:,None], axis = 1) + xyz_rot = orientation_matrix @ xyz_rot + + # add to plot + verts = [list(zip( + xyz_rot[1,:] + xyz_plot[1], + xyz_rot[0,:] + xyz_plot[0], + xyz_rot[2,:] + xyz_plot[2], + ))] + # ax.add_collection3d( + # Poly3DCollection( + # verts + # ) + # ) + collection = Poly3DCollection( + verts, + linewidths = 2.0, + alpha = 1.0, + edgecolors = 'k', + ) + face_color = [0.5, 0.5, 1] # alternative: matplotlib.colors.rgb2hex([0.5, 0.5, 1]) + if ID_plot[a1] == -1: + collection.set_facecolor((1.0,1.0,1.0)) + else: + collection.set_facecolor(atomic_colors(ID_plot[a1])) + ax.add_collection3d(collection) + + # update start point + if a1 < occ_plot.size: + theta0 = theta1 # for ID_plot in ID_all: # sub = ID == ID_plot From ef3a0e3267eb7ac975e8baa21982f3d646cf1f2e Mon Sep 17 00:00:00 2001 From: Colin Date: Wed, 18 Oct 2023 14:52:08 -0700 Subject: [PATCH 4/9] Formatting --- py4DSTEM/process/diffraction/crystal.py | 2 +- py4DSTEM/process/diffraction/crystal_viz.py | 59 ++++++++++++--------- 2 files changed, 34 insertions(+), 27 deletions(-) diff --git a/py4DSTEM/process/diffraction/crystal.py b/py4DSTEM/process/diffraction/crystal.py index 1e83aea4e..a4e45c16d 100644 --- a/py4DSTEM/process/diffraction/crystal.py +++ b/py4DSTEM/process/diffraction/crystal.py @@ -74,7 +74,7 @@ def __init__( positions, numbers, cell, - occupancy = None, + occupancy=None, ): """ Args: diff --git a/py4DSTEM/process/diffraction/crystal_viz.py b/py4DSTEM/process/diffraction/crystal_viz.py index 89fd70bc6..642c3223b 100644 --- a/py4DSTEM/process/diffraction/crystal_viz.py +++ b/py4DSTEM/process/diffraction/crystal_viz.py @@ -168,14 +168,14 @@ def plot_structure( tol = 1e-4 num_seg = 180 radius = 0.7 - zp = np.zeros(num_seg+1) + zp = np.zeros(num_seg + 1) - mark = np.ones(xyz.shape[0],dtype='bool') + mark = np.ones(xyz.shape[0], dtype="bool") for a0 in range(xyz.shape[0]): if mark[a0]: - xyz_plot = xyz[a0,:] - inds = np.argwhere(np.sum((xyz - xyz_plot)**2,axis=1) < tol) - occ_plot = occ[inds] + xyz_plot = xyz[a0, :] + inds = np.argwhere(np.sum((xyz - xyz_plot) ** 2, axis=1) < tol) + occ_plot = occ[inds] mark[inds] = False ID_plot = ID[inds] @@ -186,42 +186,50 @@ def plot_structure( occ_plot = occ_plot[0] ID_plot = ID_plot[0] - # Plot site as series of filled arcs theta0 = 0 for a1 in range(occ_plot.shape[0]): theta1 = theta0 + occ_plot[a1] * 2.0 * np.pi - theta = np.linspace(theta0,theta1,num_seg+1) + theta = np.linspace(theta0, theta1, num_seg + 1) xp = np.cos(theta) * radius yp = np.sin(theta) * radius - # Rotate towards camera - xyz_rot = np.vstack((xp.ravel(),yp.ravel(),zp.ravel())) + xyz_rot = np.vstack((xp.ravel(), yp.ravel(), zp.ravel())) if occ_plot[a1] < 1.0: - xyz_rot = np.append(xyz_rot,np.array((0,0,0))[:,None], axis = 1) - xyz_rot = orientation_matrix @ xyz_rot + xyz_rot = np.append( + xyz_rot, np.array((0, 0, 0))[:, None], axis=1 + ) + xyz_rot = orientation_matrix @ xyz_rot # add to plot - verts = [list(zip( - xyz_rot[1,:] + xyz_plot[1], - xyz_rot[0,:] + xyz_plot[0], - xyz_rot[2,:] + xyz_plot[2], - ))] + verts = [ + list( + zip( + xyz_rot[1, :] + xyz_plot[1], + xyz_rot[0, :] + xyz_plot[0], + xyz_rot[2, :] + xyz_plot[2], + ) + ) + ] # ax.add_collection3d( # Poly3DCollection( # verts # ) # ) collection = Poly3DCollection( - verts, - linewidths = 2.0, - alpha = 1.0, - edgecolors = 'k', - ) - face_color = [0.5, 0.5, 1] # alternative: matplotlib.colors.rgb2hex([0.5, 0.5, 1]) + verts, + linewidths=2.0, + alpha=1.0, + edgecolors="k", + ) + face_color = [ + 0.5, + 0.5, + 1, + ] # alternative: matplotlib.colors.rgb2hex([0.5, 0.5, 1]) if ID_plot[a1] == -1: - collection.set_facecolor((1.0,1.0,1.0)) + collection.set_facecolor((1.0, 1.0, 1.0)) else: collection.set_facecolor(atomic_colors(ID_plot[a1])) ax.add_collection3d(collection) @@ -242,12 +250,11 @@ def plot_structure( # edgecolor=[0, 0, 0], # ) # poly = PolyCollection( - # verts, - # facecolors=['r', 'g', 'b', 'y'], + # verts, + # facecolors=['r', 'g', 'b', 'y'], # alpha = 0.6) # ax.add_collection3d(poly, zs=zs, zdir='y') - # plot limit if plot_limit is None: plot_limit = np.array( From f0979af9dc871bc7c02b173c3e4755f27419a412 Mon Sep 17 00:00:00 2001 From: alex-rakowski Date: Mon, 30 Oct 2023 12:32:51 -0700 Subject: [PATCH 5/9] ASE and prismiatic files; occupancies set as one --- py4DSTEM/process/diffraction/crystal.py | 80 ++++++++++++++++++++----- 1 file changed, 65 insertions(+), 15 deletions(-) diff --git a/py4DSTEM/process/diffraction/crystal.py b/py4DSTEM/process/diffraction/crystal.py index a4e45c16d..6606a0d07 100644 --- a/py4DSTEM/process/diffraction/crystal.py +++ b/py4DSTEM/process/diffraction/crystal.py @@ -85,7 +85,7 @@ def __init__( 3 numbers: the three lattice parameters for an orthorhombic cell 6 numbers: the a,b,c lattice parameters and ɑ,β,ɣ angles for any cell 3x3 array: row vectors containing the (u,v,w) lattice vectors. - + occupancy (np.array): Partial occupancy values for each atomic site. Must match the length of positions """ # Initialize Crystal self.positions = np.asarray(positions) #: fractional atomic coordinates @@ -143,8 +143,13 @@ def __init__( # occupancy if occupancy is not None: self.occupancy = np.array(occupancy) + # check the occupancy shape makes sense + if self.occupancy.shape[0] != self.positions.shape[0]: + raise Warning( + f"Number of occupancies ({self.occupancy.shape[0]}) and atomic positions ({self.positions.shape[0]}) do not match" + ) else: - self.occupancy = None + self.occupancy = np.ones(self.positions.shape[0], dtype="intp") # pymatgen flag if "pymatgen" in sys.modules: @@ -272,6 +277,59 @@ def get_strained_crystal( else: return crystal_strained + def from_ase( + atoms, + ): + """ + Create a py4DSTEM Crystal object from an ASE atoms object + + Args: + atoms (ase.Atoms): an ASE atoms object + + """ + # get the occupancies from the atoms object + occupancies = ( + atoms.arrays["occupancies"] + if "occupancies" in atoms.arrays.keys() + else None + ) + # TODO support getting it from atoms.info['occupancy'] dictionary + xtal = Crystal( + positions=atoms.get_scaled_positions(), # fractional coords + numbers=atoms.numbers, + cell=atoms.cell.array, + occupancy=occupancies, + ) + return xtal + + def from_prismatic(filepath): + """ + Create a py4DSTEM Crystal object from an prismatic style xyz co-ordinate file + + Args: + filepath (str|Pathlib.Path): path to the prismatic format xyz file + + """ + + from ase import io + + # read the atoms using ase + atoms = io.read(filepath, format="prismatic") + + # get the occupancies from the atoms object + occupancies = ( + atoms.arrays["occupancies"] + if "occupancies" in atoms.arrays.keys() + else None + ) + xtal = Crystal( + positions=atoms.get_scaled_positions(), # fractional coords + numbers=atoms.numbers, + cell=atoms.cell.array, + occupancy=occupancies, + ) + return xtal + def from_CIF(CIF, conventional_standard_structure=True): """ Create a Crystal object from a CIF file, using pymatgen to import the CIF @@ -590,20 +648,12 @@ def calculate_structure_factors( # Calculate structure factors self.struct_factors = np.zeros(np.size(self.g_vec_leng, 0), dtype="complex64") for a0 in range(self.positions.shape[0]): - if self.occupancy is None: - self.struct_factors += f_all[:, a0] * np.exp( - (2j * np.pi) - * np.sum( - self.hkl * np.expand_dims(self.positions[a0, :], axis=1), axis=0 - ) - ) - else: - self.struct_factors += f_all[:, a0] * np.exp( - (2j * np.pi * self.occupancy[a0]) - * np.sum( - self.hkl * np.expand_dims(self.positions[a0, :], axis=1), axis=0 - ) + self.struct_factors += f_all[:, a0] * np.exp( + (2j * np.pi * self.occupancy[a0]) + * np.sum( + self.hkl * np.expand_dims(self.positions[a0, :], axis=1), axis=0 ) + ) # Divide by unit cell volume unit_cell_volume = np.abs(np.linalg.det(self.lat_real)) From f50ae93c0e908af20470885ec78847c0d912c139 Mon Sep 17 00:00:00 2001 From: alex-rakowski Date: Mon, 30 Oct 2023 12:34:17 -0700 Subject: [PATCH 6/9] fixing occupancy dtype --- py4DSTEM/process/diffraction/crystal.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/py4DSTEM/process/diffraction/crystal.py b/py4DSTEM/process/diffraction/crystal.py index 6606a0d07..60d4c490d 100644 --- a/py4DSTEM/process/diffraction/crystal.py +++ b/py4DSTEM/process/diffraction/crystal.py @@ -149,7 +149,7 @@ def __init__( f"Number of occupancies ({self.occupancy.shape[0]}) and atomic positions ({self.positions.shape[0]}) do not match" ) else: - self.occupancy = np.ones(self.positions.shape[0], dtype="intp") + self.occupancy = np.ones(self.positions.shape[0], dtype=np.float32) # pymatgen flag if "pymatgen" in sys.modules: From 71651f71b00027412bdfb9a2dfe47ea9192a4e63 Mon Sep 17 00:00:00 2001 From: Steven Zeltmann Date: Tue, 12 Dec 2023 10:43:14 -0500 Subject: [PATCH 7/9] add occupancy to CIF reader --- py4DSTEM/process/diffraction/crystal.py | 36 ++++++++++++++++++++----- 1 file changed, 29 insertions(+), 7 deletions(-) diff --git a/py4DSTEM/process/diffraction/crystal.py b/py4DSTEM/process/diffraction/crystal.py index 60d4c490d..f3e3b7739 100644 --- a/py4DSTEM/process/diffraction/crystal.py +++ b/py4DSTEM/process/diffraction/crystal.py @@ -6,6 +6,7 @@ from typing import Union, Optional from scipy.optimize import curve_fit import sys +import warnings from emdfile import tqdmnd, PointList, PointListArray from py4DSTEM.process.utils import single_atom_scatter, electron_wavelength_angstrom @@ -277,6 +278,7 @@ def get_strained_crystal( else: return crystal_strained + @staticmethod def from_ase( atoms, ): @@ -293,7 +295,12 @@ def from_ase( if "occupancies" in atoms.arrays.keys() else None ) - # TODO support getting it from atoms.info['occupancy'] dictionary + + if "occupancy" in atoms.info.keys(): + warnings.warn( + "This Atoms object contains occupancy information but it will be ignored." + ) + xtal = Crystal( positions=atoms.get_scaled_positions(), # fractional coords numbers=atoms.numbers, @@ -302,6 +309,7 @@ def from_ase( ) return xtal + @staticmethod def from_prismatic(filepath): """ Create a py4DSTEM Crystal object from an prismatic style xyz co-ordinate file @@ -330,7 +338,10 @@ def from_prismatic(filepath): ) return xtal - def from_CIF(CIF, conventional_standard_structure=True): + @staticmethod + def from_CIF( + CIF, primitive: bool = True, conventional_standard_structure: bool = True + ): """ Create a Crystal object from a CIF file, using pymatgen to import the CIF @@ -346,12 +357,13 @@ def from_CIF(CIF, conventional_standard_structure=True): parser = CifParser(CIF) - structure = parser.get_structures()[0] + structure = parser.get_structures(primitive=primitive)[0] return Crystal.from_pymatgen_structure( structure, conventional_standard_structure=conventional_standard_structure ) + @staticmethod def from_pymatgen_structure( structure=None, formula=None, @@ -448,8 +460,6 @@ def from_pymatgen_structure( else selected["structure"] ) - positions = structure.frac_coords #: fractional atomic coordinates - cell = np.array( [ structure.lattice.a, @@ -461,10 +471,22 @@ def from_pymatgen_structure( ] ) - numbers = np.array([s.species.elements[0].Z for s in structure]) + site_data = np.array( + [ + (*site.frac_coords, elem.number, comp) + for site in structure + for elem, comp in site.species.items() + ] + ) + positions = site_data[:, :3] + numbers = site_data[:, 3] + occupancies = site_data[:, 4] - return Crystal(positions, numbers, cell) + return Crystal( + positions=positions, numbers=numbers, cell=cell, occupancy=occupancies + ) + @staticmethod def from_unitcell_parameters( latt_params, elements, From 924f0e157fbe074bb972c69dc634d2084003619a Mon Sep 17 00:00:00 2001 From: Steven Zeltmann Date: Tue, 12 Dec 2023 11:06:42 -0500 Subject: [PATCH 8/9] add occupancy to structure factors --- py4DSTEM/process/diffraction/crystal.py | 12 ++++++++---- py4DSTEM/process/diffraction/crystal_bloch.py | 5 ++--- 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/py4DSTEM/process/diffraction/crystal.py b/py4DSTEM/process/diffraction/crystal.py index a2bf36a02..aa1eb8555 100644 --- a/py4DSTEM/process/diffraction/crystal.py +++ b/py4DSTEM/process/diffraction/crystal.py @@ -662,10 +662,14 @@ def calculate_structure_factors( # Calculate structure factors self.struct_factors = np.zeros(np.size(self.g_vec_leng, 0), dtype="complex64") for a0 in range(self.positions.shape[0]): - self.struct_factors += f_all[:, a0] * np.exp( - (2j * np.pi * self.occupancy[a0]) - * np.sum( - self.hkl * np.expand_dims(self.positions[a0, :], axis=1), axis=0 + self.struct_factors += ( + f_all[:, a0] + * self.occupancy[a0] + * np.exp( + (2j * np.pi) + * np.sum( + self.hkl * np.expand_dims(self.positions[a0, :], axis=1), axis=0 + ) ) ) diff --git a/py4DSTEM/process/diffraction/crystal_bloch.py b/py4DSTEM/process/diffraction/crystal_bloch.py index 6a3c9b1ac..0878feb17 100644 --- a/py4DSTEM/process/diffraction/crystal_bloch.py +++ b/py4DSTEM/process/diffraction/crystal_bloch.py @@ -27,7 +27,6 @@ def calculate_dynamical_structure_factors( tol_structure_factor: float = 0.0, recompute_kinematic_structure_factors=True, g_vec_precision=None, - verbose=True, ): """ Calculate and store the relativistic corrected structure factors used for Bloch computations @@ -92,7 +91,7 @@ def calculate_dynamical_structure_factors( # Calculate the reciprocal lattice points to include based on k_max - k_max = np.asarray(k_max) + k_max: np.ndarray = np.asarray(k_max) if recompute_kinematic_structure_factors: if hasattr(self, "struct_factors"): @@ -215,7 +214,7 @@ def get_f_e(q, Z, thermal_sigma, method): # Calculate structure factors struct_factors = np.sum( - f_e * np.exp(2.0j * np.pi * np.squeeze(self.positions[:, None, :] @ hkl)), + f_e * self.occupancy[:,None] * np.exp(2.0j * np.pi * np.squeeze(self.positions[:, None, :] @ hkl)), axis=0, ) From c70914e7ae7e432d65d451fe45ab329b37595117 Mon Sep 17 00:00:00 2001 From: Steven Zeltmann Date: Tue, 12 Dec 2023 11:06:53 -0500 Subject: [PATCH 9/9] format with black --- py4DSTEM/process/diffraction/crystal_bloch.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/py4DSTEM/process/diffraction/crystal_bloch.py b/py4DSTEM/process/diffraction/crystal_bloch.py index 0878feb17..ce8bb8622 100644 --- a/py4DSTEM/process/diffraction/crystal_bloch.py +++ b/py4DSTEM/process/diffraction/crystal_bloch.py @@ -214,7 +214,9 @@ def get_f_e(q, Z, thermal_sigma, method): # Calculate structure factors struct_factors = np.sum( - f_e * self.occupancy[:,None] * np.exp(2.0j * np.pi * np.squeeze(self.positions[:, None, :] @ hkl)), + f_e + * self.occupancy[:, None] + * np.exp(2.0j * np.pi * np.squeeze(self.positions[:, None, :] @ hkl)), axis=0, )