From 72063cdc2d202f3b52aff5b30f5ebe29cd2125b2 Mon Sep 17 00:00:00 2001 From: Chaitanya Kolluru Date: Tue, 29 Jun 2021 09:49:49 -0500 Subject: [PATCH 1/3] updates in separate file --- ingrained/construct.py | 106 +++++------ ingrained/optimize.py | 222 +++++++++++------------ ingrained/structure.py | 394 ++++++++++++++++++++--------------------- ingrained/utilities.py | 21 ++- 4 files changed, 371 insertions(+), 372 deletions(-) diff --git a/ingrained/construct.py b/ingrained/construct.py index 9782447..a57714c 100644 --- a/ingrained/construct.py +++ b/ingrained/construct.py @@ -5,13 +5,13 @@ from scipy.spatial import KDTree, Delaunay, distance_matrix # pymatgen tools -from pymatgen import MPRester +from pymatgen.ext.matproj import MPRester from pymatgen.io.xyz import XYZ from pymatgen.io.ase import AseAtomsAdaptor from pymatgen.symmetry.analyzer import SpacegroupAnalyzer -class Slab(object): - def __init__(self, chemical_formula, space_group, uvw_project, uvw_upward, tilt_angle, max_dimension): +class Slab(object): + def __init__(self, chemical_formula, space_group, uvw_project, uvw_upward, tilt_angle, max_dimension): self.chemical_formula = chemical_formula self.space_group = space_group self.unit_cell = self._query_MP() @@ -30,8 +30,8 @@ def _query_MP(self): """ mpr = MPRester("MAPI_KEY") query = mpr.query(criteria={"pretty_formula": self.chemical_formula}, properties=["structure","icsd_ids","spacegroup"]) - - # First filter by space_group if provided + + # First filter by space_group if provided if self.space_group: query = [query[i] for i in range(len(query)) if SpacegroupAnalyzer(query[i]['structure']).get_space_group_symbol()==self.space_group] @@ -52,7 +52,7 @@ def _construct_slab(self,max_dimension): Returns: A pymatgen supercell (max_dimension preserved when structure cubed after rotation) """ - # + # supercell = self.unit_cell.copy() # Each row should correspond to a lattice vector. @@ -63,17 +63,17 @@ def _construct_slab(self,max_dimension): a_start = 0.5 expansion_vector = np.array([1,1,1]) expansion_matrix = lattice_matrix - expand = np.sqrt(3)# [value_false, value_true][] + expand = np.sqrt(3)# [value_false, value_true][] while a_start < (float(max_dimension+10)*expand)/2: - # Find lattice vector with minimum length + # Find lattice vector with minimum length _ , idx = KDTree(expansion_matrix).query([0,0,0]) - + # Update expansion vector expansion_vector[idx]+= 1 - - # Update expansion matrix + + # Update expansion matrix expansion_matrix = (lattice_matrix.T*expansion_vector).T # Update bounding box @@ -94,7 +94,7 @@ def _construct_slab(self,max_dimension): return supercell def construct_oriented_slab(self,uvw_project,uvw_upward,tilt_angle,max_dimension): - """ + """ Construct a slab, oriented for a specific direction of viewing, and slice into a cube Args: @@ -105,7 +105,7 @@ def construct_oriented_slab(self,uvw_project,uvw_upward,tilt_angle,max_dimension max_dimension: A float representing the max edge length of the supercell Returns: - A pymatgen supercell (cube) oriented for a specific direction of viewing + A pymatgen supercell (cube) oriented for a specific direction of viewing """ slab = self._construct_slab(max_dimension) @@ -117,29 +117,29 @@ def construct_oriented_slab(self,uvw_project,uvw_upward,tilt_angle,max_dimension # Convert u,v,w vector to cartesian along_xyz = np.matmul(np.array(atom.get_cell()).T,uvw_project) - # Rotate coordinates and cell so that 'along_xyz' is coincident with [0,0,1] + # Rotate coordinates and cell so that 'along_xyz' is coincident with [0,0,1] atom.rotate(along_xyz,[0,0,1],rotate_cell=True) # Convert u,v,w vector to cartesian upwrd_xyz = np.matmul(np.array(atom.get_cell()).T,uvw_upward) - # Rotate coordinates and cell so that 'upwrd_xyz' is coincident with [0,1,0] + # Rotate coordinates and cell so that 'upwrd_xyz' is coincident with [0,1,0] atom.rotate(upwrd_xyz,[0,1,0],rotate_cell=True) # Rotate coordinates along z to capture tilt angle atom.rotate(tilt_angle,'z') - + bx_size = ((max_dimension+1)/2) pos = atom.get_positions() - pos -= self._centroid_coordinates(atom.get_positions()) - atom.set_positions(pos) + pos -= self._centroid_coordinates(atom.get_positions()) + atom.set_positions(pos) inidx = np.all(np.logical_and(pos>=[-bx_size,-bx_size,-bx_size],pos<=[bx_size,bx_size,bx_size]),axis=1) if not np.sum(inidx)>0: warnings.warn('Unit cell entirely outside bounding volume') - + del atom[np.logical_not(inidx)] # Enforce (0,0,0) origin and reset cell around rotated/chiseled slab @@ -152,16 +152,16 @@ def construct_oriented_slab(self,uvw_project,uvw_upward,tilt_angle,max_dimension def get_repeat_dist(self,direction="width",mode="tol"): """ - Find the approximate length needed for one full repeat of the structure along width or depth. + Find the approximate length needed for one full repeat of the structure along width or depth. Args: pymatgen_structure: The input pymatgen structure - direction: The direction along which to find repeat length + direction: The direction along which to find repeat length (width = perp to uvw_project and uvw_upward, depth = along uvw_project) mode: The decision used to accept the solution (tol = min tolerance, len = min length) Returns: - A float representing the length needed for structure to repeat. + A float representing the length needed for structure to repeat. """ slab = AseAtomsAdaptor.get_atoms(self.slab) positions_list = slab.get_positions() @@ -170,10 +170,10 @@ def get_repeat_dist(self,direction="width",mode="tol"): dtol = [] for element in chemical_symbols: - # Get all indices of + # Get all indices of indices = [int(i) for i, elem in enumerate(chemical_symbols_list) if element in elem] - position = np.array([positions_list[idx] for idx in indices]) - for select_idx in range(np.shape(position)[0]): + position = np.array([positions_list[idx] for idx in indices]) + for select_idx in range(np.shape(position)[0]): if direction == 'width': # Filter out atoms that have same a position_filt = np.array([a for a in position if a[0] != position[select_idx][0]]) @@ -188,7 +188,7 @@ def get_repeat_dist(self,direction="width",mode="tol"): adis = np.round(np.abs(position[select_idx][0]-position_filt[cb_idx][0]),4) # Compute "b" tolerance to closest proxy atom btol = np.round(np.abs((position_filt[cb_idx][1]-position[select_idx][1])),4) - elif direction == 'depth': + elif direction == 'depth': # Filter out atoms that have same c position_filt = np.array([a for a in position if a[2] != position[select_idx][2]]) # Filter out atoms that are not at same width (has tolerance) @@ -232,7 +232,7 @@ def _in_hull(p, hull): Test if points in `p` are in `hull` `p` should be a `NxK` coordinates of `N` points in `K` dimensions - `hull` is either a scipy.spatial.Delaunay object or the `MxK` array of the + `hull` is either a scipy.spatial.Delaunay object or the `MxK` array of the coordinates of `M` points in `K`dimensions for which Delaunay triangulation will be computed """ @@ -248,7 +248,7 @@ def _centroid_coordinates(crds): class TopGrain(Slab): def __init__(self, chemical_formula, space_group, uvw_project, uvw_upward, tilt_angle, max_dimension, flip_species=False): super().__init__(chemical_formula, space_group, uvw_project, uvw_upward, tilt_angle, max_dimension) - + self.height = max_dimension self.structure = self.set_in_bicrystal(flip_species=flip_species) @@ -257,11 +257,11 @@ def set_in_bicrystal(self, flip_species=False): Given grains, find the proper width/depth dimensions Args: - pymatgen_structure: The input pymatgen grain structure - width: The width of the bicrystal supercell to set grain in - height: The height of the bicrystal supercell to set grain in - depth: The depth of the bicrystal supercell to set grain in - set_in: The position (top/bottom) to set grain in bicrystal supercell + pymatgen_structure: The input pymatgen grain structure + width: The width of the bicrystal supercell to set grain in + height: The height of the bicrystal supercell to set grain in + depth: The depth of the bicrystal supercell to set grain in + set_in: The position (top/bottom) to set grain in bicrystal supercell Returns: A pymatgen grain structure set into a bicrystal supercell @@ -272,10 +272,10 @@ def set_in_bicrystal(self, flip_species=False): chiseled = atomObj.copy() pos = chiseled.get_positions() - pos -= np.min(pos,axis=0) - chiseled.set_positions(pos) + pos -= np.min(pos,axis=0) + chiseled.set_positions(pos) inidx = np.all(np.logical_and(pos>=[0,0,0],pos<[width,height,depth]),axis=1) - + del pos del chiseled[np.logical_not(inidx)] @@ -286,11 +286,11 @@ def set_in_bicrystal(self, flip_species=False): # print(np.max(pos[:,1])) # print(up_shift) pos[:,1] = pos[:,1]+up_shift - chiseled.set_positions(pos) - + chiseled.set_positions(pos) + # Reset the cell so it is big enough to eventually accomodate a second grain - chiseled.set_cell([width,2*height,depth]) - + chiseled.set_cell([width,2*height,depth]) + # Remove structure conflicts at width and depth boundaries as a result of imperfect PBC pmg = AseAtomsAdaptor.get_structure(chiseled) s = pmg.copy() @@ -327,7 +327,7 @@ def to_xyz(self,filename): class BottomGrain(Slab): def __init__(self, chemical_formula, space_group, uvw_project, uvw_upward, tilt_angle, max_dimension,flip_species=False): super().__init__(chemical_formula, space_group, uvw_project, uvw_upward, tilt_angle, max_dimension) - + self.height = max_dimension self.structure = self.set_in_bicrystal(flip_species=flip_species) @@ -336,11 +336,11 @@ def set_in_bicrystal(self,flip_species=False): Given grains, find the proper width/depth dimensions Args: - pymatgen_structure: The input pymatgen grain structure - width: The width of the bicrystal supercell to set grain in - height: The height of the bicrystal supercell to set grain in - depth: The depth of the bicrystal supercell to set grain in - set_in: The position (top/bottom) to set grain in bicrystal supercell + pymatgen_structure: The input pymatgen grain structure + width: The width of the bicrystal supercell to set grain in + height: The height of the bicrystal supercell to set grain in + depth: The depth of the bicrystal supercell to set grain in + set_in: The position (top/bottom) to set grain in bicrystal supercell Returns: A pymatgen grain structure set into a bicrystal supercell @@ -351,23 +351,23 @@ def set_in_bicrystal(self,flip_species=False): chiseled = atomObj.copy() pos = chiseled.get_positions() - pos -= np.min(pos,axis=0) - chiseled.set_positions(pos) + pos -= np.min(pos,axis=0) + chiseled.set_positions(pos) inidx = np.all(np.logical_and(pos>=[0,0,0],pos<[width,height,depth]),axis=1) - + del pos del chiseled[np.logical_not(inidx)] - + # Set grain so that max pos[:,1] is at the interface pos = chiseled.get_positions() pos[:,1] = pos[:,1] - np.min(pos[:,1])+1E-10 up_shift = ((2*height)- 1E-10)/2 - np.max(pos[:,1]) pos[:,1] = pos[:,1]+up_shift - chiseled.set_positions(pos) + chiseled.set_positions(pos) # Reset the cell so it is big enough to eventually accomodate a second grain - chiseled.set_cell([width,2*height,depth]) - + chiseled.set_cell([width,2*height,depth]) + # Remove structure conflicts at width and depth boundaries as a result of imperfect PBC pmg = AseAtomsAdaptor.get_structure(chiseled) s = pmg.copy() diff --git a/ingrained/optimize.py b/ingrained/optimize.py index 3f48027..874e1ae 100644 --- a/ingrained/optimize.py +++ b/ingrained/optimize.py @@ -10,96 +10,96 @@ import ingrained.image_ops as iop -class CongruityBuilder(object): +class CongruityBuilder(object): """ Find optimal correspondence "congruity" between a simulated and experimental image. - Posed as a β€œjigsaw puzzle” problem where the experimental image is fixed and the goal - of optimization is to find a set of simulation parameters, πœƒ, that produce an image - that can be arranged in such a way inside the experimental image that minimizes: - - 𝐽(πœƒ)= 𝛼*d_TC(πœƒ) + 𝛽*d_𝑆𝑆𝐼𝑀(πœƒ) - - - where d_TC(πœƒ) is the taxicab distance required for optimal cross-correlation-based + Posed as a β€œjigsaw puzzle” problem where the experimental image is fixed and the goal + of optimization is to find a set of simulation parameters, πœƒ, that produce an image + that can be arranged in such a way inside the experimental image that minimizes: + + 𝐽(πœƒ)= 𝛼*d_TC(πœƒ) + 𝛽*d_𝑆𝑆𝐼𝑀(πœƒ) + + - where d_TC(πœƒ) is the taxicab distance required for optimal cross-correlation-based registration after upsampling (enforces consistency in the pattern across boundaries) - - where d_SSIM(πœƒ) is the Structural Similarity IndexΒ Measure, which quantifies the visual - similarity between the simulation and experiment patch (enforces visual consistency in + - where d_SSIM(πœƒ) is the Structural Similarity IndexΒ Measure, which quantifies the visual + similarity between the simulation and experiment patch (enforces visual consistency in the content within the patches) - 𝛼, 𝛽 are weights to chosen to balance importance of the criteria """ - def __init__(self, sim_obj="", exp_img="", iter=0): + def __init__(self, sim_obj="", exp_img="", iter=0): """ Initialize a CongruityBuilder object with a ingrained.structure and experimental image. - + Args: sim_obj: (ingrained.structure) one of Bicrystal or PartialCharge stuctures exp_img: (np.array) the experimental image (consider preprocessing to enhance optimization solution) - """ + """ self.sim_obj = sim_obj self.exp_img = exp_img self.iter = iter - + def fit(self, sim_params = [], display=False, bias_x=0.0, bias_y=0.0): """ Find optimal correspondence between simulation and experiment for the specified parameters. Args: sim_params: (list) parameters required for the 'simulate_image' method in the sim_obj - display: (string) plt.show() for intermediate/final fit results + display: (string) plt.show() for intermediate/final fit results bias_x: (float) reduce search area in the x-direction by a fraction bias_y: (float) reduce search area in the y-direction by a fraction Returns: - Both 'fit' experimental and simulated image with the translation_distance between them and + Both 'fit' experimental and simulated image with the translation_distance between them and coodinates of the center pixel in the experiment - """ + """ # Simulate an image using specified parameters sim_img, sim_struct = self.sim_obj.simulate_image(sim_params=sim_params) - + # Display simulated image # self.sim_obj.display() - + # Read experimental image exp_img = self.exp_img.copy() - + # Downsample (with quantization) both simulated and experimental images # ds_by_factor = 2 ds_by_factor = 4 ds_sim_img = iop.apply_quantize_downsample(sim_img, factor=ds_by_factor) ds_exp_img = iop.apply_quantize_downsample(exp_img, factor=ds_by_factor) - + # Compute the Chi-squared (correlation) map for moving the downsampled simulation across downsampled experiment similarity_map = self.windowed_histogram_similarity(fixed=ds_exp_img, moving=ds_sim_img, bias_x=bias_x, bias_y=bias_y) # plt.imshow(similarity_map,cmap='hot'); plt.show() # Find pixel of fixed image where moving image can be positioned and not need further translation to improve similarity score ds_exp_patch, ds_stable_idxs, ds_shift_score = self.stabilize_inset_map(fixed=ds_exp_img, moving=ds_sim_img, similarity_map=similarity_map) - + # Quantizate both original simulated and experimental images (without downsampling, i.e. factor = 1) qt_sim_img = iop.apply_quantize_downsample(sim_img, factor=1) qt_exp_img = iop.apply_quantize_downsample(exp_img, factor=1) - - # Get upsampled idxs around 'ds_stable_idxs' - us_coords = iop.pixels_within_rectangle(np.array(ds_stable_idxs)[0]*ds_by_factor-ds_by_factor, - np.array(ds_stable_idxs)[1]*ds_by_factor-ds_by_factor, + + # Get upsampled idxs around 'ds_stable_idxs' + us_coords = iop.pixels_within_rectangle(np.array(ds_stable_idxs)[0]*ds_by_factor-ds_by_factor, + np.array(ds_stable_idxs)[1]*ds_by_factor-ds_by_factor, (2*ds_by_factor)+1, (2*ds_by_factor)+1) - + # Check if moving image is locked in place over upsampled (original spatial resolution) images us_exp_patch, us_stable_idxs, us_shift_score = self.stabilize_inset_coords(fixed=qt_exp_img, moving=qt_sim_img, critical_idxs=us_coords) - + # Get the shape of the moving image (size used as reference to extract exp_patch) nrows, ncols = np.shape(sim_img) - + # Get best fit experiment patch after search (us_exp_patch is quantized, we want patch from original image w/o quantization!) exp_patch = exp_img[int(us_stable_idxs[0]-(nrows-1)/2):int(us_stable_idxs[0]+(nrows-1)/2)+1,\ int(us_stable_idxs[1]-(ncols-1)/2):int(us_stable_idxs[1]+(ncols-1)/2)+1] - + # Convert to grayscale for viewing gs_exp_img = iop.scale_pixels(exp_patch,mode='grayscale') gs_sim_img = iop.scale_pixels(sim_img,mode='grayscale') - + if display: # Display matching pair after downsample stabilized plt.imshow(np.hstack([ds_exp_patch,15*np.ones((np.shape(ds_exp_patch)[0],5)),ds_sim_img]),cmap='hot'); plt.axis('off'); @@ -113,7 +113,7 @@ def fit(self, sim_params = [], display=False, bias_x=0.0, bias_y=0.0): plt.imshow(np.hstack([gs_exp_img,255*np.ones((np.shape(gs_exp_img)[0],5)),gs_sim_img]),cmap='hot'); plt.axis('off'); plt.show() return sim_img, sim_struct, exp_patch, us_shift_score, us_stable_idxs - + def fit_gb(self, sim_params = [], display=False, bias_x=0.0, bias_y=0.15): """ Find optimal correspondence between simulation and experiment for the specified parameters. @@ -125,40 +125,40 @@ def fit_gb(self, sim_params = [], display=False, bias_x=0.0, bias_y=0.15): bias_y: (float) reduce search area in the y-direction by a fraction Returns: - Both 'fit' experimental and simulated image with the translation_distance between them and + Both 'fit' experimental and simulated image with the translation_distance between them and coodinates of the center pixel in the experiment - """ + """ return self.fit(sim_params=sim_params, display=False, bias_x=bias_x, bias_y=bias_y) def windowed_histogram_similarity(self, fixed="", moving="", bias_x=0.0, bias_y=0.0): """ Compute the Chi squared distance metric between the moving image across all pixels in the fixed image, in an attempt to locate the moving image (or highly similar regions) within the original fixed image. - + Args: - fixed: (np.array) experimental image + fixed: (np.array) experimental image moving: (np.array) a simulated image patch (<= size of fixed image) bias_x: (float) reduce search area in the x-direction by a fraction bias_y: (float) reduce search area in the y-direction by a fraction - + Returns: A similarity map indicating regions of high correlation (reciprocal of Chi-squared similarity) between moving and fixed. """ # Compute histogram for simulated (moving) image, and normalize moving_hist, _ = np.histogram(moving.flatten(), bins=16, range=(0, 16)) moving_hist = moving_hist.astype(float) / np.sum(moving_hist) - + # Find appropriate size for a disk shaped mask that will define the shape of the sliding window radius = int((2/3) * np.max(np.shape(moving))) - selem = disk(radius) # A disk is (2 * radius + 1) wide - + selem = disk(radius) # A disk is (2 * radius + 1) wide + # Compute normalized windowed histogram feature vector for each pixel in the fixed image px_histograms = rank.windowed_histogram(fixed, selem=selem, n_bins=moving_hist.shape[0]) # Reshape moving histogram to (1,1,N) for broadcast when we want to use it in # arithmetic operations with the windowed histograms from the image moving_hist = moving_hist.reshape((1, 1) + moving_hist.shape) - + # Compute Chi-squared distance metric: sum((X-Y)^2 / (X+Y)); denom = px_histograms + moving_hist denom[denom == 0] = np.infty @@ -168,38 +168,38 @@ def windowed_histogram_similarity(self, fixed="", moving="", bias_x=0.0, bias_y= # Use reciprocal of Chi-squared similarity measure to create a full similarity map full_map = 1 / (chi_sqr + 1.0e-4) similarity_map = full_map.copy() - + # Get half length/width + 1 of moving image and use to define border size on fixed (where no sliding window comparisons permitted) pix_row = int(np.ceil(np.shape(moving)[0]/2) + 1) pix_col = int(np.ceil(np.shape(moving)[1]/2) + 1) - + # Bias values can further restrict area where comparisons are made pix_row += int(((np.shape(similarity_map)[0] - 2*pix_row) * bias_y)/2) pix_col += int(((np.shape(similarity_map)[1] - 2*pix_col) * bias_x)/2) - - # Construct the final correlation map + + # Construct the final correlation map brd_col = np.zeros((pix_row,np.shape(similarity_map)[1])) brd_row = np.zeros((np.shape(similarity_map)[0],pix_col)) - + similarity_map = np.vstack([brd_col,similarity_map[pix_row:-pix_row,::],brd_col]) similarity_map = np.hstack([brd_row,similarity_map[::,pix_col:-pix_col],brd_row]) return similarity_map - + def stabilize_inset_map(self, fixed="", moving="", similarity_map=""): """ - Find position (center pixel) on the fixed image, where the moving image can be placed without - further translation (or with minimal translation) to improve registration. The similarity map + Find position (center pixel) on the fixed image, where the moving image can be placed without + further translation (or with minimal translation) to improve registration. The similarity map defines the order that the center pixels are tested. - + Args: - fixed: (np.array) experimental image + fixed: (np.array) experimental image moving: (np.array) a simulated image patch (<= size of fixed image) similarity_map: (np.array) map of pixels measuring similarity between moving and fixed at a each center pixel - + Returns: - The stable center pixel and the correponding experimental patch (same size as moving image), - as well as the absolute value of the shift vector (in pixels) required to register the images. - A stable return is currently defined as having a shift error of zero. + The stable center pixel and the correponding experimental patch (same size as moving image), + as well as the absolute value of the shift vector (in pixels) required to register the images. + A stable return is currently defined as having a shift error of zero. """ # Get the shape of the moving image nrows, ncols = np.shape(moving) @@ -207,7 +207,7 @@ def stabilize_inset_map(self, fixed="", moving="", similarity_map=""): # Find the index (row, col) with the current highest similarity current_idx = np.unravel_index(similarity_map.argmax(), similarity_map.shape) - # Keep track of the current best index and shift score + # Keep track of the current best index and shift score current_best = None for i in range(np.sum(similarity_map > 0)): @@ -217,9 +217,9 @@ def stabilize_inset_map(self, fixed="", moving="", similarity_map=""): int(current_idx[1]-(ncols-1)/2):int(current_idx[1]+(ncols-1)/2)+1] # Find shift required for maximum cross correlation - shift, __, __ = phase_cross_correlation(moving,fixed_patch) + shift, __, __ = phase_cross_correlation(moving,fixed_patch) - # Calculate shift score + # Calculate shift score shift_score = np.sum(np.abs(shift)) # Check to see if shift score improved, if so, record the index! @@ -246,24 +246,24 @@ def stabilize_inset_map(self, fixed="", moving="", similarity_map=""): def stabilize_inset_coords(self, fixed="", moving="", critical_idxs=""): """ - Find position (center pixel) on the fixed image, where the moving image can be placed without - further translation (or with minimal translation) to improve registration. The critical indices + Find position (center pixel) on the fixed image, where the moving image can be placed without + further translation (or with minimal translation) to improve registration. The critical indices provided define the order that the center pixels are tested. - + Args: - fixed: (np.array) experimental image + fixed: (np.array) experimental image moving: (np.array) a simulated image patch (<= size of fixed image) - critical_idxs: (list) critical indices to test for stability - + critical_idxs: (list) critical indices to test for stability + Returns: - The stable center pixel and the correponding experimental patch (same size as moving image), - as well as the absolute value of the shift vector (in pixels) required to register the images. - A stable return is currently defined as having a shift error of zero. + The stable center pixel and the correponding experimental patch (same size as moving image), + as well as the absolute value of the shift vector (in pixels) required to register the images. + A stable return is currently defined as having a shift error of zero. """ # Get the shape of the moving image nrows, ncols = np.shape(moving) - # Keep track of the current best index and shift score + # Keep track of the current best index and shift score current_best = None for current_idx in critical_idxs: @@ -273,9 +273,9 @@ def stabilize_inset_coords(self, fixed="", moving="", critical_idxs=""): int(current_idx[1]-(ncols-1)/2):int(current_idx[1]+(ncols-1)/2)+1] # Find shift required for maximum cross correlation - shift, __, __ = phase_cross_correlation(moving,fixed_patch) + shift, __, __ = phase_cross_correlation(moving,fixed_patch) - # Calculate shift score + # Calculate shift score shift_score = np.sum(np.abs(shift)) # Check to see if shift score improved, if so, record the index! @@ -313,7 +313,7 @@ def display_panel(self, moving="", critical_idx="", iternum="", score="", cmap=' fixed = self.exp_img.copy() fixed = iop.scale_pixels(fixed, mode='grayscale') moving = iop.scale_pixels(moving, mode='grayscale') - + nrows, ncols = np.shape(moving) rcrds = iop.pixels_within_rectangle(int(critical_idx[1]-(ncols-1)/2), int(critical_idx[0]-(nrows-1)/2), ncols, nrows) @@ -323,7 +323,7 @@ def display_panel(self, moving="", critical_idx="", iternum="", score="", cmap=' for i in range(len(rcrds)): entry = rcrds[i] base_img[entry[1],entry[0]] = filler[i] - + # Likely to make simulation square (not always good). Can delete this! # moving = iop.apply_resize(moving, np.shape(fixed)) @@ -339,14 +339,14 @@ def display_panel(self, moving="", critical_idx="", iternum="", score="", cmap=' moving = np.hstack([brdx,moving,brdx]) except: pass - + # Ensure that all images are the same size (same as experiment) moving = iop.apply_resize(moving, np.shape(fixed)) base_img = iop.apply_resize(base_img, np.shape(fixed)) fig, axes = plt.subplots(nrows=1, ncols=3) fig.add_gridspec(nrows=1, ncols=3).update(wspace=0.0, hspace=0.0) - + axes[0].imshow(fixed, interpolation='quadric', cmap=cmap) axes[0].set_title(title_list[0],fontsize=8.5) axes[0].axis('off') @@ -358,7 +358,7 @@ def display_panel(self, moving="", critical_idx="", iternum="", score="", cmap=' axes[2].imshow(base_img, interpolation='quadric', cmap=cmap) axes[2].set_title(title_list[2],fontsize=9) axes[2].axis('off') - + pnt1 = np.min(rcrds,axis=0) pnt2 = np.max(rcrds,axis=0) axes[2].plot([pnt1[0],pnt1[0]],[pnt1[1],pnt2[1]],lw=0.5,color="#8AFF30") @@ -368,7 +368,7 @@ def display_panel(self, moving="", critical_idx="", iternum="", score="", cmap=' if score != "": axes[2].text(0.53*np.shape(base_img)[1],1.07*np.shape(base_img)[0],"FOM: "+"{:8.5f}".format(float(score)),va='center',ha='left',fontsize=9) - + if iternum != "": axes[2].text(0,1.07*np.shape(base_img)[0],"iteration: "+str(iternum),va='center',ha='left',fontsize=9) @@ -380,19 +380,19 @@ def display_panel(self, moving="", critical_idx="", iternum="", score="", cmap=' def taxicab_ssim_objective(self,x): """ - Objective function used to quantify how well a given set of input imaging paramters, x, - produce an image that can be arranged in such a way inside the experimental image that minimizes - the custom taxicab_ssim objective function: - - 𝐽(πœƒ)= 𝛼*d_TC(πœƒ) + 𝛽*d_𝑆𝑆𝐼𝑀(πœƒ) - - - where d_TC(πœƒ) is the taxicab distance required for optimal cross-correlation-based + Objective function used to quantify how well a given set of input imaging paramters, x, + produce an image that can be arranged in such a way inside the experimental image that minimizes + the custom taxicab_ssim objective function: + + 𝐽(πœƒ)= 𝛼*d_TC(πœƒ) + 𝛽*d_𝑆𝑆𝐼𝑀(πœƒ) + + - where d_TC(πœƒ) is the taxicab distance required for optimal cross-correlation-based registration after upsampling (enforces consistency in the pattern across boundaries) - - where d_SSIM(πœƒ) is the Structural Similarity IndexΒ Measure, which quantifies the visual - similarity between the simulation and experiment patch (enforces visual consistency in + - where d_SSIM(πœƒ) is the Structural Similarity IndexΒ Measure, which quantifies the visual + similarity between the simulation and experiment patch (enforces visual consistency in the content within the patches) - + Args: x: (np.array) set of parameters to test @@ -402,12 +402,12 @@ def taxicab_ssim_objective(self,x): xfit = [a for a in x[:-2]] + [int(a) for a in x[-2::]] try: sim_img, sim_struct, exp_patch, shift_score, __ = self.fit(xfit, display=False); - except Exception as e: + except Exception as e: print(e) sim_img = None - + if sim_img is not None: - match_ssim = iop.score_ssim(sim_img, exp_patch, win_size=35) + match_ssim = iop.score_ssim(sim_img, exp_patch) fom = 0.1*(shift_score)+match_ssim else: fom = 9999 @@ -422,19 +422,19 @@ def taxicab_ssim_objective(self,x): def taxicab_ssim_objective_gb(self,x): """ - Objective function used to quantify how well a given set of input imaging paramters, x, - produce an image that can be arranged in such a way inside the experimental image that minimizes - the custom taxicab_ssim objective function: - - 𝐽(πœƒ)= 𝛼*d_TC(πœƒ) + 𝛽*d_𝑆𝑆𝐼𝑀(πœƒ) - - - where d_TC(πœƒ) is the taxicab distance required for optimal cross-correlation-based + Objective function used to quantify how well a given set of input imaging paramters, x, + produce an image that can be arranged in such a way inside the experimental image that minimizes + the custom taxicab_ssim objective function: + + 𝐽(πœƒ)= 𝛼*d_TC(πœƒ) + 𝛽*d_𝑆𝑆𝐼𝑀(πœƒ) + + - where d_TC(πœƒ) is the taxicab distance required for optimal cross-correlation-based registration after upsampling (enforces consistency in the pattern across boundaries) - - where d_SSIM(πœƒ) is the Structural Similarity IndexΒ Measure, which quantifies the visual - similarity between the simulation and experiment patch (enforces visual consistency in + - where d_SSIM(πœƒ) is the Structural Similarity IndexΒ Measure, which quantifies the visual + similarity between the simulation and experiment patch (enforces visual consistency in the content within the patches) - + Uses special "fit_gb" which prioritizes seach closer to the center (horizontal) of the image. Args: @@ -446,12 +446,12 @@ def taxicab_ssim_objective_gb(self,x): xfit = [a for a in x[:-2]] + [int(a) for a in x[-2::]] try: sim_img, sim_struct, exp_patch, shift_score, __ = self.fit_gb(xfit, display=False); - except Exception as e: + except Exception as e: print(e) sim_img = None if sim_img is not None: - match_ssim = iop.score_ssim(sim_img, exp_patch, win_size=35) + match_ssim = iop.score_ssim(sim_img, exp_patch) fom = 0.1*(shift_score)+match_ssim else: fom = 9999 @@ -466,16 +466,16 @@ def taxicab_ssim_objective_gb(self,x): def find_correspondence(self,objective='taxicab_ssim', optimizer='Powell', initial_solution="", search_mode="stm"): """ - Wrapper around scipy.optimize.minimize solvers, to find optimal correspondence between simulation and experiment + Wrapper around scipy.optimize.minimize solvers, to find optimal correspondence between simulation and experiment by minimizing taxicab_ssim_objective. - + Args: - objective: (string) - optimizer: (string) + objective: (string) + optimizer: (string) initial_solution: (np.array) constraint_list: (list) search_mode: (string) Specify search mode ('gb' for STEM grain boundaries, and 'stm' for STM images) - + Returns: """ if os.path.isfile(os.getcwd()+"/progress.txt"): @@ -486,10 +486,10 @@ def find_correspondence(self,objective='taxicab_ssim', optimizer='Powell', initi os.remove(os.getcwd()+"/progress.txt") else: pass - + self.iter = 1 # Keep this in because its useful for constraints that need to be enforced by the optimizer. - # Currently, constraints are enforced by a combination of clamping (in sim_params), or throwing exceptions in the + # Currently, constraints are enforced by a combination of clamping (in sim_params), or throwing exceptions in the # image simulation steps constraints = [] # def make_constraint_list(c): @@ -504,7 +504,7 @@ def find_correspondence(self,objective='taxicab_ssim', optimizer='Powell', initi # constraints = make_constraint_list(constraint_list) if objective == "taxicab_ssim": return minimize(self.taxicab_ssim_objective, initial_solution, method='COBYLA',tol=1E-6,options={'disp': True, 'rhobeg': 0.25, 'catol': 0.01}, constraints=constraints) - + if optimizer == 'Powell': if objective == "taxicab_ssim": return minimize(self.taxicab_ssim_objective, initial_solution, method='Powell',tol=1E-6,options={'disp': True}) @@ -514,10 +514,10 @@ def find_correspondence(self,objective='taxicab_ssim', optimizer='Powell', initi # constraints = make_constraint_list(constraint_list) if objective == "taxicab_ssim": return minimize(self.taxicab_ssim_objective_gb, initial_solution, method='COBYLA',tol=1E-6,options={'disp': True, 'rhobeg': 0.25, 'catol': 0.01}, constraints=constraints) - + if optimizer == 'Powell': if objective == "taxicab_ssim": return minimize(self.taxicab_ssim_objective_gb, initial_solution, method='Powell',tol=1E-6,options={'disp': True}) else: - print("Search mode {} not understood! Select either 'gb' for STEM grain boundaries, or 'stm' for STM images") \ No newline at end of file + print("Search mode {} not understood! Select either 'gb' for STEM grain boundaries, or 'stm' for STM images") diff --git a/ingrained/structure.py b/ingrained/structure.py index 33e85d5..52f6f82 100644 --- a/ingrained/structure.py +++ b/ingrained/structure.py @@ -35,36 +35,36 @@ class PartialCharge(object): def __init__(self, config_file=""): """ Initialize a PartialCharge object with the path to a PARCHG file. - + Args: config_file: (string) path to a PARCHG file - """ + """ self.config_file = config_file - self.parchg = Chgcar.from_file(config_file) - - # Actual volumetric charge data from PARCHG - self.chgdata = self.parchg.data['total'] - + self.parchg = Chgcar.from_file(config_file) + + # Actual volumetric charge data from PARCHG + self.chgdata = self.parchg.data['total'] + # Structure associated w/ volumetric data - self.structure = self.parchg.structure - + self.structure = self.parchg.structure + # Fraction of cell occupied by structure - self.zmax = self.structure.cart_coords[:, 2].max()/self.structure.lattice.c - + self.zmax = self.structure.cart_coords[:, 2].max()/self.structure.lattice.c + # Get the real space sizes assocated with each voxel in volumetric data self.lwh = ( self.structure.lattice.a/self.parchg.dim[0] , \ self.structure.lattice.b/self.parchg.dim[1] , \ self.structure.lattice.c/self.parchg.dim[2] ) - + # Image fields self.cell = None self.image = None self.sim_params = None - + def _get_stm_vol(self, z_below="", z_above=""): - """ + """ Get the volumetric slab associated with the input depths. - + Args: z_below: (float) thickness or depth from top in Angstroms z_above: (float) distance above the surface to consider @@ -74,7 +74,7 @@ def _get_stm_vol(self, z_below="", z_above=""): Indice of the top slice (in reference to full slab) """ if self._parameter_check_z(z_below=z_below, z_above=z_above): - + # Get fractional coordinates for z_below and z_above z_below /= self.structure.lattice.c z_above /= self.structure.lattice.c @@ -86,18 +86,18 @@ def _get_stm_vol(self, z_below="", z_above=""): # Get volumetric slice (nzmin:nzmax,nrows,ncols) from charge data rhos = np.swapaxes(np.swapaxes(self.chgdata.copy(),0,2),1,2)[nzmin:nzmax,:,:]/self.structure.volume return rhos, nzmax - + def _shift_sites(self): """ - Translate all coordinates such that the maximum of the + Translate all coordinates such that the maximum of the maximum of the charge density is around the center of the cell. """ - + # Translate atomic positions so zmax is in the center of the cell - + self.structure.translate_sites(indices=range(self.structure.num_sites), \ vector=[0,0,self.structure.lattice.c*(.5-self.zmax)], \ frac_coords=False) @@ -118,22 +118,22 @@ def _shift_sites(self): new_dat[i].append(np.concatenate((dat_slice[zmax_index+int(ngzf/2):], dat_slice[:zmax_index], dat_slice[zmax_index:zmax_index+int(ngzf/2)]))) - - + + self.chgdata = new_dat self.zmax = self.structure.cart_coords[:, 2].max()/ \ - self.structure.lattice.c + self.structure.lattice.c self.structure.to('POSCAR','POSCAR_centered') - + def _get_stm_cell(self, z_below="", z_above="", r_val="", r_tol=""): """ - Calculate a single STM image cell from DFT-simulation. - - NOTE: Skew in the volumetric grid due to non-orthogonal lattice vectors + Calculate a single STM image cell from DFT-simulation. + + NOTE: Skew in the volumetric grid due to non-orthogonal lattice vectors in the imaging plane is NOT accounted for in the creation of the image_cell. Use get_stm_image (wrapper) to return an image that contains the proper skew! - + Args: z_below: (float) thickness or depth from top in Angstroms z_above: (float) distance above the surface to consider @@ -142,14 +142,14 @@ def _get_stm_cell(self, z_below="", z_above="", r_val="", r_tol=""): Returns: A numpy array (np.float64) of the simulated image tile. - + Credit: Chaitanya Kolluru """ - # Get volumetric slab corresponding to z_below and z_above parameters + # Get volumetric slab corresponding to z_below and z_above parameters rhos, nzmax = self._get_stm_vol(z_below=z_below, z_above=z_above) - + if self._parameter_check_r(r_val=r_val, r_tol=r_tol, rhos=rhos): - + # If absdiff between charge density value and r_val not is not within r_tol, switch pixels to off (i.e. -9E9) rhos[abs(rhos - r_val) >= r_tol] = -9E9 @@ -170,11 +170,11 @@ def _get_stm_cell(self, z_below="", z_above="", r_val="", r_tol=""): self.cell = image return image - + def simulate_image(self, sim_params=[]): """ Calculate a full displayable STM image from DFT-simulation. - + Args: sim_params: (list) specifying the following simulation parameters - z_below: (float) thickness or depth from top in (Γ…) @@ -190,16 +190,16 @@ def simulate_image(self, sim_params=[]): - sigma: (float) standard deviation for gaussian kernel used in postprocessing - crop_height: (int) final (cropped) image height in pixels - crop_width: (int) final (cropped) image width in pixels - + Returns: - A numpy array (np.float64) of the full simulated image. + A numpy array (np.float64) of the full simulated image. """ - # Enforce max stretch/squeeze and max/min shear value (both directions 0.30) - sim_params[4] = sorted((-0.50, sim_params[4], 0.50))[1] - sim_params[5] = sorted((-0.50, sim_params[5], 0.50))[1] - sim_params[6] = sorted((-0.50, sim_params[6], 0.50))[1] - sim_params[7] = sorted((-0.50, sim_params[7], 0.50))[1] - + # Enforce max stretch/squeeze and max/min shear value (both directions 0.30) + sim_params[4] = 0 #sorted((-0.50, sim_params[4], 0.50))[1] + sim_params[5] = 0 #sorted((-0.50, sim_params[5], 0.50))[1] + sim_params[6] = 0 #sorted((-0.50, sim_params[6], 0.50))[1] + sim_params[7] = 0 #sorted((-0.50, sim_params[7], 0.50))[1] + # Clamps the given pixel size to a range between 0.05 and 0.40 (Γ…) sim_params[8] = sorted((-360, sim_params[8], 360))[1] @@ -207,7 +207,7 @@ def simulate_image(self, sim_params=[]): sim_params[9] = sorted((0.05, sim_params[9], 0.40))[1] # Clamps the given sigma value to a range between 0 and 10 - sim_params[10] = sorted((0, sim_params[10], 10))[1] + sim_params[10] = 0 #sorted((0, sim_params[10], 10))[1] # sim_params[10] = 0 # In case we want to get rid of blur completely # Simulate the image cell @@ -215,14 +215,14 @@ def simulate_image(self, sim_params=[]): # Construct a larger cell by repeating image_cell in both directions, NREPS = 8 img_tiled = np.tile(self.cell,(16,16)) - + # Add shear to account for a non-orthogonal unit cell (if necessary) - img_tiled = self._apply_lattice_shear(img_tiled) + img_tiled = self._apply_lattice_shear(img_tiled) # Resize image based on pixel_size img_tiled = iop.apply_resize(img_tiled, np.array([int(a) for a in np.shape(img_tiled) * np.array(self.lwh[::-1][1::]) * (1./sim_params[9])])) - # Rotate tiled image according to rotation value + # Rotate tiled image according to rotation value img_tiled = iop.apply_rotation(img_tiled, sim_params[8]) # Apply any specified postprocessing shear @@ -236,79 +236,79 @@ def simulate_image(self, sim_params=[]): # Enforce odd sizes on width and length of cropped image (min_length = 25, max_length = current length) sim_params[-1] = sorted((37,int(2 * np.floor(sim_params[-1]/2) + 1),int(2 * np.floor((np.shape(img_tiled)[1]-2)/2) + 1)))[1] sim_params[-2] = sorted((37,int(2 * np.floor(sim_params[-2]/2) + 1),int(2 * np.floor((np.shape(img_tiled)[0]-2)/2) + 1)))[1] - + # Apply crop to image img_tiled = iop.apply_crop(img_tiled,sim_params[-1],sim_params[-2]) # Apply gaussian blur - image = iop.apply_blur(img_tiled, sigma=sim_params[10]) + image = iop.apply_blur(img_tiled, sigma=sim_params[10]) # If sucessful, record parameters! self.sim_params = sim_params self.image = image return image, None - - def simulate_image_random(self, rotation="", pix_size="", crop_height="", crop_width=""): + + def simulate_image_random(self, rotation="", pix_size="", crop_height="", crop_width=""): """ Randomly simulate a valid image. - + Args: rotation: (float) image rotation angle (in degrees CCW) pix_size: (float) real-space edge length of a single pixel (in Γ…) crop_height: (int) final (cropped) image height in pixels crop_width: (int) final (cropped) image width in pixels - + Returns: A numpy array (np.float64) of the full simulated image. """ sim_params = self._random_initialize(rotation=rotation, pix_size=pix_size, crop_height=crop_height, crop_width=crop_width) img_rand, __ = self.simulate_image(sim_params) return img_rand - + def _random_initialize(self, rotation="", pix_size="", crop_height="", crop_width=""): """ Randomly select a set of parameters to simulate a valid image - + Args: rotation: (float) image rotation angle (in degrees CCW) pix_size: (float) real-space edge length of a single pixel (in Γ…) crop_height: (int) final (cropped) image height in pixels crop_width: (int) final (cropped) image width in pixels - + Returns: A numpy array of sim_params (tuple) specifying the necessary simulation parameters for 'simulate_image' """ slab_thickness = self.structure.cart_coords[:, 2].max() - self.structure.cart_coords[:, 2].min() - + z_below = random.uniform(-1+1E-10,(0.5*slab_thickness)-1E-10) z_above = random.uniform(0,(1 - self.zmax) * self.structure.lattice.c) - + rhos, __ = self._get_stm_vol(z_below=z_below, z_above=z_above) - + r_val = random.uniform((rhos.max()/3)+1E-10, (rhos.max())-1E-10) r_tol = random.uniform(0.0001+1E-10,(0.999*r_val)-1E-10) - + x_shear = random.uniform(-0.3,0.3) - y_shear = random.uniform(-0.3,0.3) + y_shear = random.uniform(-0.3,0.3) x_stretch = random.uniform(-0.3,0.3) y_stretch = random.uniform(-0.3,0.3) - + sigma = random.uniform(0,4) if rotation != "": rotation = (rotation + random.choice(list(range(-360,360,90))) + random.uniform(-2,2))%360 else: rotation = random.choice(list(range(-360,360,1))) - + if pix_size != "": pix_size = pix_size + random.uniform(-0.015,0.015) else: pix_size = random.uniform(0.05,0.40) - + if crop_width == "": - crop_height = crop_width = int(2 * np.floor(random.uniform(33,121)/2) + 1) + crop_height = crop_width = int(2 * np.floor(random.uniform(33,121)/2) + 1) return [z_below, z_above, r_val, r_tol, x_shear, y_shear, x_stretch, y_stretch, rotation, pix_size, sigma, crop_height, crop_width] - + def _image_size_check(self, img): """ """ @@ -322,40 +322,40 @@ def _image_size_check(self, img): def _parameter_check_z(self, z_below="", z_above=""): """ Validate choice of depth parameters to ensure feasibility of STM image calculation - + Args: z_below: (float) thickness or depth from top in Angstroms z_above: (float) distance above the surface to consider - + Returns: A boolean (pass or fail) """ slab_thickness = self.structure.cart_coords[:, 2].max() - self.structure.cart_coords[:, 2].min() # Check to make sure zbelow is within valid depth range (in Angstroms) - if -1 <= z_below <= min([0.5*slab_thickness,4]): + if -1 <= z_below <= min([0.5*slab_thickness,2]): pass else: - raise Exception("ParameterError: z_below = {0} outside valid range [{1} to {2:.5}]".format(z_below,-1,min([0.5*slab_thickness,4]))) - + raise Exception("ParameterError: z_below = {0} outside valid range [{1} to {2:.5}]".format(z_below,-1,min([0.5*slab_thickness,2]))) + # Check to make sure z_above is within bounds of structure (within close proximity to surface, specifically) - if 0 < z_above <= ((1 - self.zmax) * 0.5 * self.structure.lattice.c): + if 1 < z_above <= ((1 - self.zmax) * self.structure.lattice.c) - 2: pass else: - raise Exception("ParameterError: z_above = {0} outside structure (must be <= {1:.5})".format(z_above,((1 - self.zmax) * 0.5 * self.structure.lattice.c))) + raise Exception("ParameterError: z_above = {0} outside structure (must be <= {1:.5})".format(z_above,((1 - self.zmax) * self.structure.lattice.c)) - 2) return True - + def _parameter_check_r(self, r_val="", r_tol="", rhos=""): """ Validate choice of charge density parameters to ensure feasibility of STM image calculation - + Args: r_val: (float) isosurface charge density plane r_tol: (float) tolerance to consider while determining isosurface rhos: (numpy array) the volumetric charge density slab - + Returns: - A boolean (pass or fail) + A boolean (pass or fail) """ # Check to make sure r_val is within valid isosurface (rmax/3 to rmax) if rhos.max()/3 < r_val < rhos.max(): @@ -369,23 +369,23 @@ def _parameter_check_r(self, r_val="", r_tol="", rhos=""): else: raise Exception("ParameterError: r_tol = {0} outside valid range ({1:.5} to {2:.5})".format(r_tol,0.0001,0.999*r_val)) return True - + def _apply_lattice_shear(self, img): """ - Apply centered shear of volumetric grid along lattice vectors. + Apply centered shear of volumetric grid along lattice vectors. Args: - img: (numpy array) - + img: (numpy array) + Returns: A square numpy array (np.float64) of the sheared image. - - NOTE: This function has limited use cases. Add more for robust handling! Currently assumes "a" - lattice vector in x-direction and "b" lattice vector resolves into x and y components. + + NOTE: This function has limited use cases. Add more for robust handling! Currently assumes "a" + lattice vector in x-direction and "b" lattice vector resolves into x and y components. """ # Get lattice vector a and b lva, lvb = self.parchg.structure.lattice.matrix[0], self.parchg.structure.lattice.matrix[1] - + # Use count of nonzero entries as method to detect vectors resolved along multiple directions if np.count_nonzero(lvb) == 1 and np.count_nonzero(lva) == 1: x_sh = 0 @@ -395,15 +395,15 @@ def _apply_lattice_shear(self, img): y_sh = 0 else: raise Exception('Shear detected, but not accounted for in image formed from volumetric grid!') - return iop.apply_shear(img, x_sh, y_sh) - + return iop.apply_shear(img, x_sh, y_sh) + def simulation_summary(self, iter="", verbose=False): """ Write the parameter summary to screen (or config_file). - + Args: config_file: (string) name of the write file - + Returns: None """ @@ -414,7 +414,7 @@ def simulation_summary(self, iter="", verbose=False): β€’ (r_val, r_tol) : {2} β€’ (x, y) shear (frac) : {3} β€’ (x, y) stretch (frac) : {4} - β€’ sigma (Gaussian blur) : {5} + β€’ sigma (Gaussian blur) : {5} β€’ rotation (deg CCW) : {6} β€’ pix_size (Γ…) : {7} β€’ img_size (pixels) : {8}""".format(iter,(z_below, z_above),(r_val, r_tol),(x_shear, y_shear),(x_stretch, y_stretch),sigma,rotation,pix_size,(crop_height,crop_width)) @@ -422,7 +422,7 @@ def simulation_summary(self, iter="", verbose=False): print(summary) else: return summary - + def display(self): """ Display the simulated STM image. @@ -433,17 +433,17 @@ def display(self): class Bicrystal(object): """ - This class contains two oriented supercell slabs prepared from structure queries to the Materials Project, - as well as a calculator to fuse the slabs together into a bicrystal, simulate a convolution HAADF STEM image, + This class contains two oriented supercell slabs prepared from structure queries to the Materials Project, + as well as a calculator to fuse the slabs together into a bicrystal, simulate a convolution HAADF STEM image, and perform postprocessing imaging manipulations. """ def __init__(self, config_file="", write_poscar=False, poscar_file=""): """ Initialize a Bicrystal object with the path to a bicrystal slabs (json) file. - + Args: config_file: (string) path to a bicrystal slabs file - """ + """ if poscar_file != "": poscar = Poscar.from_file(poscar_file) @@ -456,7 +456,7 @@ def __init__(self, config_file="", write_poscar=False, poscar_file=""): self.slab_1 = slab["slab_1"] self.slab_2 = slab["slab_2"] - + # Copy simulation folder from install directory to cwd try: shutil.copytree(os.path.dirname(__file__)+'/simulation', os.getcwd()+'/simulation') @@ -472,15 +472,15 @@ def __init__(self, config_file="", write_poscar=False, poscar_file=""): self.slab_1['uvw_project'], self.slab_1['uvw_upward'],\ self.slab_1['tilt_angle'], self.slab_1['max_dimension'],\ self.slab_1['flip_species']) - + # Actual structure for the bottom grain (slab positioned into bottom of bicrystal cell) bot_grain = BottomGrain(self.slab_2["chemical_formula"], self.slab_2["space_group"],\ self.slab_2['uvw_project'], self.slab_2['uvw_upward'],\ self.slab_2['tilt_angle'], self.slab_2['max_dimension'],\ - self.slab_2['flip_species']) - + self.slab_2['flip_species']) + structure, top_grain_fit, bot_grain_fit, strain_info = self._get_base_structure(top_grain, bot_grain, constraints) - + self.top_grain = top_grain_fit self.bot_grain = bot_grain_fit self.structure = structure @@ -503,7 +503,7 @@ def __init__(self, config_file="", write_poscar=False, poscar_file=""): # self.cell, __ = self._get_image_cell() # self.image = None # self.sim_params = None - + # self.lw = (self.structure.lattice.b/np.shape(self.cell)[1], # self.structure.lattice.c/np.shape(self.cell)[0]) @@ -518,33 +518,33 @@ def __init__(self, config_file="", write_poscar=False, poscar_file=""): def _get_base_structure(self, top_grain_structure, bot_grain_structure, constraints): """ Take grains and fuse them into a simple bicrystal, subject to input constraints - + Args: top_grain_structure (): an ingrained.contruct.TopGrain() object bot_grain_structure (): an ingrained.contruct.BottomGrain() object constraints (dict): dictionary of geometric constraints for bicrystal construction - + Returns: - The basic fused bicrystal (pymatgen structure) w/ corresponding ingrained.contruct.TopGrain() + The basic fused bicrystal (pymatgen structure) w/ corresponding ingrained.contruct.TopGrain() and ingrained.contruct.BottomGrain() objects after average expand/contract values applied """ structure, top_grain, bot_grain, strain_info= self._fuse_grains(top_grain_structure, bot_grain_structure, constraints) structure = self._adjust_interface_width(structure=structure, interface_width_1=constraints['interface_1_width'], interface_width_2=constraints['interface_2_width']) structure = self._remove_interface_collisions(structure=structure, collision_removal_1=constraints['collision_removal'][0], collision_removal_2=constraints['collision_removal'][1]) - + return structure, top_grain, bot_grain, strain_info - + def _fuse_grains(self, top_grain_structure, bot_grain_structure, constraints): """ - Take current grain structures and expand/contract in width/depth to minimize strain - between grains and assign same width/depth computed as the average of the ideal - expand/contract values for each grain. + Take current grain structures and expand/contract in width/depth to minimize strain + between grains and assign same width/depth computed as the average of the ideal + expand/contract values for each grain. Args: top_grain_structure (): an ingrained.contruct.TopGrain() object bot_grain_structure (): an ingrained.contruct.BottomGrain() object constraints (dict): dictionary of geometric constraints for bicrystal construction - + Return: A bicrystal (pymatgen structure) from fusion of top and bottom grains """ @@ -557,8 +557,8 @@ def _fuse_grains(self, top_grain_structure, bot_grain_structure, constraints): # Apply appropriate integer expansions to grain structures top_grain_structure.structure.make_supercell([top_n_width,1,top_n_depth]) bot_grain_structure.structure.make_supercell([bot_n_width,1,bot_n_depth]) - - # Get expanded width/depth values of the grains + + # Get expanded width/depth values of the grains widths = [top_grain_structure.structure.lattice.a,bot_grain_structure.structure.lattice.a] depths = [top_grain_structure.structure.lattice.c,bot_grain_structure.structure.lattice.c] @@ -578,32 +578,32 @@ def _fuse_grains(self, top_grain_structure, bot_grain_structure, constraints): bot_grain = AseAtomsAdaptor.get_atoms(bot_grain_structure.structure) # Strain to coincidence - top_grain.set_cell([width_bc,height_bc,depth_bc],scale_atoms=True) + top_grain.set_cell([width_bc,height_bc,depth_bc],scale_atoms=True) bot_grain.set_cell([width_bc,height_bc,depth_bc],scale_atoms=True) - # Set bottom grain as the gb structure + # Set bottom grain as the gb structure gb_full = bot_grain.copy() - + # Get chemical symbols and positions for top grain top_grain_sym = top_grain.get_chemical_symbols() top_grain_pos = top_grain.get_positions() - - # Insert top_grain info into gb structure (on top of the exisiting 'bottom' structure) + + # Insert top_grain info into gb structure (on top of the exisiting 'bottom' structure) for i in range(len(top_grain_sym)): gb_full.append(Atom(top_grain_sym[i],top_grain_pos[i])) # Get positions, cell information and shuffle axes gb_pos = gb_full.get_positions() gb_abc = gb_full.get_cell_lengths_and_angles()[0:3] - gb_full.set_positions(np.vstack([np.vstack([gb_pos[:,2],gb_pos[:,0]]),gb_pos[:,1]]).T) + gb_full.set_positions(np.vstack([np.vstack([gb_pos[:,2],gb_pos[:,0]]),gb_pos[:,1]]).T) gb_full.set_cell([gb_abc[2],gb_abc[0],gb_abc[1]]) - + # Convert back to pymatgen structure and resolve any boundary conflicts gb_full_structure = AseAtomsAdaptor.get_structure(gb_full) gb_full_structure = self._enforce_pb(gb_full_structure) gb_full_structure = self._resolve_pb_conflicts(gb_full_structure) return gb_full_structure, top_grain_structure, bot_grain_structure, [strain_top_width,strain_top_depth,strain_bot_width,strain_bot_depth] - + def _adjust_interface_width(self, structure="", interface_width_1=0, interface_width_2=0): """ Adjust spacing between both interfaces between grains, with option to remove collision conflicts @@ -615,9 +615,9 @@ def _adjust_interface_width(self, structure="", interface_width_1=0, interface_w Return: A bicrystal (pymatgen structure) that reflects specified interface widths - """ - # Get copy of current structure - bc = structure.copy() + """ + # Get copy of current structure + bc = structure.copy() if interface_width_1 != 0: # Update atom coordinate positions based on modified length in the c-direction from interface_width_1 new_crds, new_spec = [], [] @@ -628,112 +628,112 @@ def _adjust_interface_width(self, structure="", interface_width_1=0, interface_w else: new_crds.append([cx/bc.lattice.a,cy/bc.lattice.b,cz/(bc.lattice.c+interface_width_1)]) new_spec.append(str(bc.species[idx]).split('Element')[0]) - + # Create new 'Lattice' from parameters and a new pymatgen 'Structure' with updated positions lattice = Lattice.from_parameters(a=bc.lattice.a, b=bc.lattice.b, c=(bc.lattice.c+interface_width_1), alpha=90, beta=90, gamma=90) new_struct_1 = Structure(lattice, new_spec, new_crds) # Get copy of current structure with atom positions bc = new_struct_1.copy() - + if interface_width_2 != 0: crds = bc.cart_coords # Shift all coords up half distance of interface_width_2 crds[:,2] = (crds[:,2]-np.min(crds[:,2])) + interface_width_2/2 - + # Update atom coordinate positions based on modified length in the c-direction from interface_width_2 new_crds, new_spec = [], [] for idx in range(len(bc)): cx,cy,cz = crds[idx] new_crds.append([cx/bc.lattice.a,cy/bc.lattice.b,cz/(np.max(crds[:,2]) + interface_width_2/2)]) new_spec.append(str(bc.species[idx]).split('Element')[0]) - + # Create new 'Lattice' from parameters and a new pymatgen 'Structure' with updated positions lattice = Lattice.from_parameters(a=bc.lattice.a, b=bc.lattice.b, c=(np.max(crds[:,2]) + interface_width_2/2), alpha=90, beta=90, gamma=90) new_struct_2 = Structure(lattice, new_spec, new_crds) - + # Get copy of current structure with atom positions bc = new_struct_2.copy() - return bc - + return bc + def _remove_interface_collisions(self, structure="", collision_removal_1=False, collision_removal_2=False, collision_distance=1): """ Remove atoms at the interface that are within 'collision_distance' of another atom, starting removal with atoms at the bottom of the interface region - + Args: structure (pymatgen structure): bicrystal structure collision_removal_1 (boolean): remove atoms within 'collision_distance' in volume around interface_1 collision_removal_2 (boolean): remove atoms within 'collision_distance' in volume around interface_2 collision_distance (float): atoms less than or equal to this distance are considered collisions - + Return: A bicrystal (pymatgen structure) with collisions removed within the interface volumes """ # Get copy of current structure bc = structure.copy() - + # Interface 1 collision removal if collision_removal_1: - + # Get all atomic coordinates pos = bc.cart_coords - + # Define region surrounding interface where we check for conflict (3Γ… region) iw_zone = ((bc.lattice.c/2)-1.5,(bc.lattice.c/2)+1.5) - + # Retrieve all coordinates belonging to the iw_zone interface_coords = pos[np.where((pos[:,2]>=iw_zone[0]) & (pos[:,2]<=iw_zone[1]))] - + # Ensure they are sorted such that lower coordinates are tested first interface_coords = interface_coords[interface_coords[:,2].argsort()] - + for coord in interface_coords: # Get array of distances between coord and coords of current structure dist_check = distance_matrix([coord], bc.cart_coords)[0] - + # If an atom is within collision_distance, remove! if np.sum(dist_check <= collision_distance) > 1: indx = list(dist_check).index(0) - bc.remove_sites([indx]) - + bc.remove_sites([indx]) + # Interface 2 collision removal if collision_removal_2: - + # Get all atomic coordinates pos = bc.cart_coords - + # Define region surrounding interface where we check for conflict (3Γ… region) iw_zone = (0,1.5) - + # Retrieve all coordinates belonging to the iw_zone interface_coords = pos[np.where((pos[:,2]>=iw_zone[0]) & (pos[:,2]<=iw_zone[1]))] - + # Ensure they are sorted such that lower coordinates are tested first interface_coords = interface_coords[interface_coords[:,2].argsort()] - + for coord in interface_coords: # Move coordinate across PBC to check for conflict coord[2] = coord[2]+bc.lattice.c - + # Get array of distances between coord and coords of current structure dist_check = distance_matrix([coord], bc.cart_coords)[0] - + # If an atom is within collision_distance, remove! if np.sum(dist_check <= collision_distance) > 1: indx = list(dist_check).index(0) - bc.remove_sites([indx]) + bc.remove_sites([indx]) return bc - - def _get_image_cell(self, defocus=1, interface_width=0, pix_size=0.15, view=False): + + def _get_image_cell(self, defocus=1, interface_width=0, pix_size=0.15, view=False): """ - Calculate a single HAADF STEM image cell from bicrystal structure. - + Calculate a single HAADF STEM image cell from bicrystal structure. + Args: defocus: (float) controls degree to which edges blur in microscopy image (Γ…) interface_width: (float) spacing b/w max pos of bottom_grain and min pos of top grain (Γ…) pix_size: (float) real-space pixel size (Γ…) - view: (boolean) option to display cell after simulation + view: (boolean) option to display cell after simulation Returns: A numpy array (np.float64) of the simulated image tile. @@ -741,15 +741,15 @@ def _get_image_cell(self, defocus=1, interface_width=0, pix_size=0.15, view=Fals """ # Get copy of current structure - bc =self.structure.copy() - + bc =self.structure.copy() + # Apply additional interface width adjustments (on top of interface_1_width which is set during initialization) bc = self._adjust_interface_width(structure=bc, interface_width_1=interface_width) bc = self._remove_interface_collisions(structure=bc, collision_removal_1=True) - # Use pixel_size to define shape of the output image + # Use pixel_size to define shape of the output image pixx, pixy = np.round(np.array(bc.lattice.abc)/pix_size)[1::].astype(np.int) - + fmt = '% 4d', '% 8.4f', '% 9.4f', '% 9.4f', '% 4.2f', '% 4.3f' # Write input structure file required for Kirkland STEM simulation with open(os.getcwd()+'/simulation/SAMPLE.XYZ', "w") as sf: @@ -771,9 +771,9 @@ def _get_image_cell(self, defocus=1, interface_width=0, pix_size=0.15, view=Fals # Write parameter file required for Kirkland STEM simulation with open(os.getcwd()+'/simulation/params.txt', "w") as pf: - pf.write('SAMPLE.XYZ\n1 1 1\nSAMPLE.TIF\n'+str(pixx)+" "+str(pixy)+"\n") - pf.write("200 0 0 0 30\n100 150\nEND\n"+str(defocus)+"\n0") - + pf.write('SAMPLE.XYZ\n1 1 1\nSAMPLE.TIF\n'+str(pixx)+" "+str(pixy)+"\n") + pf.write("200 0 0 0 30\n100 150\nEND\n"+str(defocus)+"\n0") + # Simulate image with Kirkland incostem with cd(os.getcwd()+'/simulation'): subprocess.call("./incostem-osx", stdout=subprocess.PIPE) @@ -792,11 +792,11 @@ def _get_image_cell(self, defocus=1, interface_width=0, pix_size=0.15, view=Fals self.lw = (bc.lattice.b/np.shape(image)[1],bc.lattice.c/np.shape(image)[0]) return image, bc - + def simulate_image(self, sim_params=[]): """ Calculate a full displayable STM image from DFT-simulation. - + Args: sim_params: (list) specifying the following simulation parameters - pix_size: (float) real-space pixel size (Γ…) @@ -808,9 +808,9 @@ def simulate_image(self, sim_params=[]): - y_stretch: (float) fractional amt stretch (+) or compression (-) in y - crop_height: (int) final (cropped) image height in pixels - crop_width: (int) final (cropped) image width in pixels - + Returns: - A numpy array (np.float64) of the full simulated image. + A numpy array (np.float64) of the full simulated image. """ # Clamps the given pixel size to a range if self.pix_size != None: @@ -819,13 +819,13 @@ def simulate_image(self, sim_params=[]): # sim_params[0] = sorted((0.1, sim_params[0], 0.4))[1] sim_params[0] = sorted((0.05, sim_params[0], 0.4))[1] - # Clamp the interface_width to a range between -2 and 2 (Γ…) + # Clamp the interface_width to a range between -2 and 2 (Γ…) sim_params[1] = sorted((-2, sim_params[1], 2))[1] - - # Clamp the defocus to a range between 0.5 and 1.75 (Γ…) + + # Clamp the defocus to a range between 0.5 and 1.75 (Γ…) sim_params[2] = sorted((0.5, sim_params[2], 1.75))[1] - # Enforce max stretch/squeeze and max/min shear value (both directions 0.20) + # Enforce max stretch/squeeze and max/min shear value (both directions 0.20) sim_params[3] = sorted((-0.20, sim_params[3], 0.20))[1] sim_params[4] = sorted((-0.20, sim_params[4], 0.20))[1] sim_params[5] = sorted((-0.20, sim_params[5], 0.20))[1] @@ -834,9 +834,9 @@ def simulate_image(self, sim_params=[]): # Simulate the image cell old_struct = self.structure.copy() # print(self.lw) - + img_cell, new_struct = self._get_image_cell(defocus=sim_params[2], interface_width=sim_params[1], pix_size=sim_params[0], view=False) - + # if new_struct == old_struct: # print("No modifications were made to the structure!") # else: @@ -846,7 +846,7 @@ def simulate_image(self, sim_params=[]): # Construct a larger cell by repeating image_cell in both directions img_tiled = np.tile(img_cell,(1,4)) - + # Resize image based on pixel_size img_tiled = iop.apply_resize(img_tiled, np.array([int(a) for a in np.shape(img_tiled) * np.array(self.lw[::-1]) * (1./sim_params[0])])) @@ -859,40 +859,40 @@ def simulate_image(self, sim_params=[]): # Enforce odd sizes on width and length of cropped image (min_length = 25, max_length = current length) sim_params[-1] = sorted((35,int(2 * np.floor(sim_params[-1]/2) + 1),int(2 * np.floor((np.shape(img_tiled)[1]-2)/2) + 1)))[1] sim_params[-2] = sorted((35,int(2 * np.floor(sim_params[-2]/2) + 1),int(2 * np.floor((np.shape(img_tiled)[0]-2)/2) + 1)))[1] - + # Apply crop to image - image = iop.apply_crop(img_tiled,sim_params[-1],sim_params[-2]) + image = iop.apply_crop(img_tiled,sim_params[-1],sim_params[-2]) # If sucessful, record parameters! self.sim_params = sim_params self.image = image - return image, new_struct - - def simulate_image_random(self, pix_size="", crop_height="", crop_width=""): + return image, new_struct + + def simulate_image_random(self, pix_size="", crop_height="", crop_width=""): """ Randomly simulate a valid image. - + Args: pix_size: (float) real-space edge length of a single pixel (in Γ…) crop_height: (int) final (cropped) image height in pixels crop_width: (int) final (cropped) image width in pixels - + Returns: A numpy array (np.float64) of the full simulated image. """ sim_params = self._random_initialize(pix_size=pix_size, crop_height=crop_height, crop_width=crop_width) img_rand, __ = self.simulate_image(sim_params) return img_rand - + def _random_initialize(self, pix_size="", crop_height="", crop_width=""): """ Randomly select a set of parameters to simulate a valid image - + Args: pix_size: (float) real-space edge length of a single pixel (in Γ…) crop_height: (int) final (cropped) image height in pixels crop_width: (int) final (cropped) image width in pixels - + Returns: A numpy array of sim_params (tuple) specifying the necessary simulation parameters for 'simulate_image' """ @@ -901,15 +901,15 @@ def _random_initialize(self, pix_size="", crop_height="", crop_width=""): defocus = random.uniform(0.5+(1E-10),1.75-(1E-10)) x_shear = random.uniform(-0.1,0.1) - y_shear = random.uniform(-0.1,0.1) + y_shear = random.uniform(-0.1,0.1) x_stretch = random.uniform(-0.1,0.1) y_stretch = random.uniform(-0.1,0.1) - + if pix_size != "": pix_size = pix_size + random.uniform(-0.015,0.015) else: pix_size = random.uniform(0.05,0.40) - + if crop_width == "": if crop_height != "": crop_width = int(random.uniform(0.80,1.2)*crop_height) @@ -923,11 +923,11 @@ def _random_initialize(self, pix_size="", crop_height="", crop_width=""): crop_height = int(random.uniform(75,201)) return [pix_size, interface_width, defocus, x_shear, y_shear, x_stretch, y_stretch, crop_height, crop_width] - + @staticmethod def _find_optimal_expansion(x, y, min_len=5, max_len=100): """ - Given two floats, x and y, return integers a and b such that + Given two floats, x and y, return integers a and b such that |ax - by| β‰ˆ 0 subject to min/max constraints Args: @@ -937,23 +937,23 @@ def _find_optimal_expansion(x, y, min_len=5, max_len=100): max_length: max value of an approximate multiple Returns: - Smallest integers that best approximate the objective + Smallest integers that best approximate the objective """ alist = [x*i for i in range(1,100) if x*i < max_len and x*i > min_len] blist = [y*i for i in range(1,100) if y*i < max_len and y*i > min_len] dmat = distance.cdist(np.array(alist).reshape(-1,1),np.array(blist).reshape(-1,1),'euclidean') result = min((min((v, c) for c, v in enumerate(row)), r) for r, row in enumerate(dmat)) return int(alist[result[1]]/x), int(blist[result[0][1]]/y), np.min(dmat) - + @staticmethod def _enforce_pb(pymatgen_structure): """ - Wrap atoms that have fixed coords outside cell back to their coord in cell. + Wrap atoms that have fixed coords outside cell back to their coord in cell. Notice this sometimes happens when using AseAtomsAdaptor.get_structure() to convert ASE to pymatgen - + Args: pymatgen_structure - + Returns: pymatgen_structure """ @@ -961,15 +961,15 @@ def _enforce_pb(pymatgen_structure): for i in range(len(s.frac_coords)): s[i] = s[i].specie, [s[i].frac_coords[0]%1, s[i].frac_coords[1]%1, s[i].frac_coords[2]%1] return s - + @staticmethod def _resolve_pb_conflicts(pymatgen_structure): """ If atoms near the 0 bound (width or depth only) interfere w/ atoms @ max bounds w/ PBC's, delete them! - + Args: pymatgen_structure - + Returns: pymatgen_structure """ @@ -990,14 +990,14 @@ def _resolve_pb_conflicts(pymatgen_structure): i+=1 s.remove_sites(indel) return s - + def simulation_summary(self, iter="", verbose=False): """ Write the parameter summary to screen (or config_file). - + Args: config_file: (string) name of the write file - + Returns: None """ diff --git a/ingrained/utilities.py b/ingrained/utilities.py index fdd054c..8115bf4 100644 --- a/ingrained/utilities.py +++ b/ingrained/utilities.py @@ -1,6 +1,6 @@ import time import sys, os -sys.path.append('../../../') +#sys.path.append('../../../') import numpy as np from skimage.io import imsave import matplotlib.pyplot as plt @@ -39,9 +39,9 @@ def print_frames(config_file="", poscar_file="", exp_img="", exp_title="", progr 3,8,26 : multiple individual frames (#3, #8, #26) 4-74 : frames within a range (#4 through and including #74) 1:2:30 : start, step, stop sequence of frames - + >> Selection: """) - + else: decision = frame_selection @@ -74,8 +74,8 @@ def print_frames(config_file="", poscar_file="", exp_img="", exp_title="", progr progress = np.genfromtxt(progress_file, delimiter=',') save_path = os.path.dirname(os.path.realpath(progress_file)) - os.makedirs(save_path+"/frames/", exist_ok=True) - + os.makedirs(save_path+"/frames/", exist_ok=True) + if decision.lower() == 'all': for i in range(np.shape(progress)[0]): x = progress[i] @@ -185,7 +185,7 @@ def print_frames(config_file="", poscar_file="", exp_img="", exp_title="", progr congruity.display_panel(moving=sim_img, critical_idx=stable_idxs, score=x[-1], iternum=int(x[0]), cmap=cmap, title_list=describe_frames, savename=save_path+'/frames/ingrained'+'-'+str(idx).zfill(5)+'.png') enablePrint() Printer("Writing frame {} to file".format(idx)) - + else: Printer("Selection \"{}\" not understood!".format(decision)) print() @@ -195,14 +195,14 @@ def print_frames(config_file="", poscar_file="", exp_img="", exp_title="", progr def prepare_fantastx_input(config_file="",poscar_file="", exp_img="", progress_file=""): - + if poscar_file != "": # Initialize a Bicrystal object with the path to the initial poscar sim_obj = Bicrystal(poscar_file=poscar_file); else: # Initialize a Bicrystal object with the path to the slab json file sim_obj = Bicrystal(config_file=config_file); - + bias_y=0.15 # Initialize a ConguityBuilder with PARCHG and experimental image @@ -212,8 +212,8 @@ def prepare_fantastx_input(config_file="",poscar_file="", exp_img="", progress_f progress = np.genfromtxt(progress_file, delimiter=',') save_path = os.path.dirname(os.path.realpath(progress_file)) - os.makedirs(save_path+"/fantastx_start/", exist_ok=True) - + os.makedirs(save_path+"/fantastx_start/", exist_ok=True) + best_idx = int(np.argmin(progress[:,-1])) x = progress[best_idx] xfit = x[1:-1] @@ -234,4 +234,3 @@ def prepare_fantastx_input(config_file="",poscar_file="", exp_img="", progress_f np.savetxt(save_path+"/fantastx_start/"+'sim_params', sim_params) np.save(save_path+"/fantastx_start/"+"experiment.npy",exp_patch) - \ No newline at end of file From 459dcd41e29603715504f31bc1c88671dea5a3bd Mon Sep 17 00:00:00 2001 From: Chaitanya Kolluru Date: Wed, 21 Jul 2021 16:03:12 -0500 Subject: [PATCH 2/3] added setup.py and updated readme in dev_ch --- README.md | 85 +++++++++++++++++++++++++------- setup.py | 143 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 210 insertions(+), 18 deletions(-) create mode 100644 setup.py diff --git a/README.md b/README.md index 9c9dae5..dd4086f 100644 --- a/README.md +++ b/README.md @@ -7,41 +7,90 @@ Grain boundary structure initialization from atomic-resolution HAADF STEM imaging ([CdTe [110]-[110] tilt grain boundary](https://aip.scitation.org/doi/10.1063/1.5123169)) alt text - + STM structure validation from simulated images of partial charge densities volumes from DFT-simulation ([Cu2O(111)](https://pubs.rsc.org/en/content/articlelanding/2018/cp/c8cp06023a#!divAbstract)) alt text ## Getting started - + ### Requirements To install *ingrained*, you need: - * python 3.x + * python 3.x * a conda package manager to configure the basic installation environment * [Materials Project API key](https://materialsproject.org/open) ### Installation -The following will install a local copy of *ingrained* and create a custom [conda](https://docs.conda.io/en/latest/) development environment for both STEM and STM functionalities. +It is recommended to start installation on a new conda environment using Anaconda. + +Load Anaconda if available as a library +```sh +module load conda +``` +Or in some systems, Anaconda comes with python module. So +```sh +module load python +``` + +If above methods fail, download Anaconda for the system [here](https://docs.conda.io/en/latest/miniconda.html). Follow default instructions and install Anaconda. Restart the terminal so that the installation takes effect. + +Update conda and add conda-forge to the Anaconda channels & set higher priority. Skip this step if conda-forge is already added before. Install everything from conda-forge to be consistent & reduce compatibility issues across different pacakges. +```sh +conda update -n base -c defaults conda +conda config --add channels conda-forge +``` + +Create a new conda environment with path. Provide *path* to the new conda environment. Add the new environment (*fantastx*) to the system path. Install python 3. + +```sh +conda create -yp ~/miniconda3/envs/ingrained +conda activate ~/miniconda3/envs/ingrained +export PATH=~/miniconda3/envs/ingrained/bin:$PATH +conda install python=3 +``` + +Install Pymatgen (Installs Numpy, Scipy, Matplotlib) with pip instead of conda. + +```sh +pip install pymatgen +``` -Clone this repo and navigate to the root directory: +Get latest ingrained from github ```sh -git clone https://github.com/MaterialEyes/ingrained -cd ingrained +git clone https://github.com/MaterialEyes/ingrained.git ``` -Create the environment from the environment.yml file and activate it: + +Change branch to dev_ch (for the following installation) ```sh -conda env create -f utils/environment.yml -conda activate ingrained +git checkout dev_ch ``` + Set your [Materials Project API key](https://materialsproject.org/open) "MAPI_KEY" environment variable: -```python +```sh ./utils/set_key.sh mYaPiKeY ``` -Install additional (external) packages to parse standard experimental files (i.e. dm3, sxm, etc...) + +Install *pydm3* and *pySPM* readers (for .dm3 and .sxm files) with the following script ```sh ./utils/install_parsers.sh ``` + +Install ase +```sh +pip install --upgrade --user ase +``` + +Add the local bin (path) to sys +```sh +export PATH=~/.local/bin:$PATH +``` + +Install using python setup.py file with *develop* option. This approach allows installation of Ingrained in already existing conda environments as well. (For example, when using with Fantastx) +```sh +python setup.py develop +``` + ## *ingrained* workflow 1. Isolate a clean region of an experimental image, and preprocess if necessary (*i.e. Wiener filter, contrast enhancement, etc.*) 2. Initialize a *structure* object from an appropriate classes in [structure.py](https://github.com/MaterialEyes/ingrained/blob/master/ingrained/structure.py). For grain boundaries, a Bicrystal() object requires a path to a configuration file, config.json. For STM structure validation, a PartialCharge() object requires a path to a PARCHG file. @@ -54,7 +103,7 @@ The [clean_templates](https://github.com/MaterialEyes/ingrained/tree/master/exam ### Ex. 1: Grain boundary structure initialization from HAADF STEM -Here, we outline the contents of the [example](https://github.com/MaterialEyes/ingrained/tree/master/examples/completed_runs/gb/CdTe_110_CdTe_110_tilt) given for a [CdTe [110]-[110] tilt grain boundary](https://aip.scitation.org/doi/10.1063/1.5123169) structure +Here, we outline the contents of the [example](https://github.com/MaterialEyes/ingrained/tree/master/examples/completed_runs/gb/CdTe_110_CdTe_110_tilt) given for a [CdTe [110]-[110] tilt grain boundary](https://aip.scitation.org/doi/10.1063/1.5123169) structure (1) Experimental target image ([HAADF149.dm3](https://github.com/MaterialEyes/ingrained/blob/master/examples/completed_runs/gb/CdTe_110_CdTe_110_tilt/HAADF149.dm3)). @@ -179,7 +228,7 @@ Here, we outline the contents of the [example](https://github.com/MaterialEyes/i This file will need to be downloaded, as it is not stored in the repo. ```sh -bash parchg_download.sh +bash parchg_download.sh ``` (3) Python script to execute steps of the *ingrained* workflow ([run.py]()). @@ -245,9 +294,9 @@ Summary of optimization parameters: (4) Python script for final relaxation of z_above parameter ([relax_z_above.py](https://github.com/MaterialEyes/ingrained/tree/master/examples/completed_runs/stm/Cu2O_111/relax_z_above.py)) -Run only after find_correspondence has completed. Currently, z_above acts as a free variable -once it exceeds the top of the charge surface, and optimization may make this value arbitarily high when the best solution -is at the top of the surface. This function ensures that z_above represents a physically meaningful distance +Run only after find_correspondence has completed. Currently, z_above acts as a free variable +once it exceeds the top of the charge surface, and optimization may make this value arbitarily high when the best solution +is at the top of the surface. This function ensures that z_above represents a physically meaningful distance (Γ…) above the surface. ### Output @@ -257,7 +306,7 @@ Execution of the sequence outlined in [run.py]() will produce: * [strain_info.txt](https://github.com/MaterialEyes/ingrained/blob/master/examples/completed_runs/gb/CdTe_110_CdTe_110_tilt/strain_info) - (Ex. 1 only) record of the amount of strain in each bicrystal grain (given as % along width and depth) * [progress.txt](https://github.com/MaterialEyes/ingrained/blob/master/examples/completed_runs/gb/CdTe_110_CdTe_110_tilt/progress.txt) - record of the optimization solution and the respective figure-of-merit (FOM) at each optimization step. -Additional tools are included to view and write optimization progress to a movie +Additional tools are included to view and write optimization progress to a movie * [print_frames.py](https://github.com/MaterialEyes/ingrained/blob/master/examples/completed_runs/gb/CdTe_110_CdTe_110_tilt/print_frames.py) - writes specified optimization steps (frames) to custom image panels (.png) * [make_movie.sh](https://github.com/MaterialEyes/ingrained/blob/master/examples/completed_runs/gb/CdTe_110_CdTe_110_tilt/make_movie.sh) - wrapper around [FFmpeg](https://ffmpeg.org/) to create a movie from the sequence of panels created in [print_frames.py](https://github.com/MaterialEyes/ingrained/blob/master/examples/completed_runs/gb/CdTe_110_CdTe_110_tilt/print_frames.py) diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..d22563b --- /dev/null +++ b/setup.py @@ -0,0 +1,143 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +# Note: To use the 'upload' functionality of this file, you must: +# $ pipenv install twine --dev + +import io +import os +import sys +from shutil import rmtree + +from setuptools import find_packages, setup, Command +from setuptools.command.install import install + +# Package meta-data. +NAME = 'ingrained' +DESCRIPTION = 'Ingrained is a python toolkit for embedding STEM simulations of grain boundaries structures into experimental atomic resolution HAADF STEM images' +URL = 'https://github.com/MaterialEyes/ingrained' +EMAIL = 'eschwenk@u.northwestern.edu' +AUTHOR = 'Eric Schwenker' +REQUIRES_PYTHON = '>=3.6.0' +VERSION = '0.1.0' + +# What packages are required for this module to be executed? +REQUIRED = [ + 'opencv-python', + 'pillow>=2.3.1', + 'scipy', + 'numpy', + 'pymatgen', + 'pySPM', + 'sklearn', + 'scikit-image', + 'imutils', + 'jupyter' + ], + +# What packages are optional? +EXTRAS = { + # 'fancy feature': ['django'], +} + +# The rest you shouldn't have to touch too much :) +# ------------------------------------------------ +# Except, perhaps the License and Trove Classifiers! +# If you do change the License, remember to change the Trove Classifier for that! + +here = os.path.abspath(os.path.dirname(__file__)) + +# Import the README and use it as the long-description. +# Note: this will only work if 'README.md' is present in your MANIFEST.in file! +try: + with io.open(os.path.join(here, 'README.md'), encoding='utf-8') as f: + long_description = '\n' + f.read() +except FileNotFoundError: + long_description = DESCRIPTION + +# Load the package's __version__.py module as a dictionary. +about = {} +if not VERSION: + project_slug = NAME.lower().replace("-", "_").replace(" ", "_") + with open(os.path.join(here, project_slug, '__version__.py')) as f: + exec(f.read(), about) +else: + about['__version__'] = VERSION + +class UploadCommand(Command): + """Support setup.py upload.""" + + description = 'Build and publish the package.' + user_options = [] + + @staticmethod + def status(s): + """Prints things in bold.""" + print('\033[1m{0}\033[0m'.format(s)) + + def initialize_options(self): + pass + + def finalize_options(self): + pass + + def run(self): + try: + self.status('Removing previous builds...') + rmtree(os.path.join(here, 'dist')) + except OSError: + pass + + self.status('Building Source and Wheel (universal) distribution..') + os.system('{0} setup.py sdist bdist_wheel --universal'.format(sys.executable)) + + self.status('Uploading the package to PyPI via Twine..') + os.system('twine upload dist/*') + + self.status('Pushing git tags..') + os.system('git tag v{0}'.format(about['__version__'])) + os.system('git push --tags') + + sys.exit() + +# Where the magic happens: +setup( + name=NAME, + version=about['__version__'], + description=DESCRIPTION, + long_description=long_description, + long_description_content_type='text/markdown', + author=AUTHOR, + author_email=EMAIL, + python_requires=REQUIRES_PYTHON, + url=URL, + packages=find_packages(include=['ingrained', 'ingrained.*']), + # If your package is a single module, use this instead of 'packages': + # py_modules=['ingrained'], + + # entry_points={ + # 'console_scripts': ['mycli=mymodule:cli'], + # }, + install_requires=REQUIRED, + extras_require=EXTRAS, + include_package_data=True, + license='MIT', + classifiers=[ + # Trove classifiers + # Full list: https://pypi.python.org/pypi?%3Aaction=list_classifiers + 'License :: OSI Approved :: MIT License', + 'Programming Language :: Python', + 'Programming Language :: Python :: 3', + 'Programming Language :: Python :: 3.6', + 'Programming Language :: Python :: Implementation :: CPython', + 'Programming Language :: Python :: Implementation :: PyPy' + ], + #data_files=[(sys.prefix+'/ingrained/simulation/', ['ingrained/simulation/incostem-mac']), + # (sys.prefix+'/ingrained/simulation/', ['ingrained/simulation/incostem-linux']), + # ], + dependency_links=[], + # $ setup.py publish support. + cmdclass={ + 'upload': UploadCommand, + }, +) From 6fd6810b8acd9caadc17c4ab94eefc6d2fe85b4c Mon Sep 17 00:00:00 2001 From: Chaitanya Kolluru Date: Wed, 21 Jul 2021 16:11:59 -0500 Subject: [PATCH 3/3] updated instructions --- README.md | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index dd4086f..755d36a 100644 --- a/README.md +++ b/README.md @@ -22,7 +22,7 @@ To install *ingrained*, you need: * [Materials Project API key](https://materialsproject.org/open) ### Installation -It is recommended to start installation on a new conda environment using Anaconda. +It is recommended to start installation on a new conda environment using Anaconda. However, if needed to install Ingrained in an existing conda environment, skip this step and proceed to next step. Load Anaconda if available as a library ```sh @@ -50,6 +50,11 @@ export PATH=~/miniconda3/envs/ingrained/bin:$PATH conda install python=3 ``` +For installation in an existing conda environment, skip the above steps and start from here by activating the existing conda environment. +```sh +conda activate +``` + Install Pymatgen (Installs Numpy, Scipy, Matplotlib) with pip instead of conda. ```sh