diff --git a/.github/workflows/linting.yml b/.github/workflows/linting.yml index a690cc8..49a23da 100644 --- a/.github/workflows/linting.yml +++ b/.github/workflows/linting.yml @@ -8,7 +8,7 @@ jobs: strategy: max-parallel: 4 matrix: - python-version: [3.8] + python-version: [3.11] steps: - uses: actions/checkout@v2 @@ -27,5 +27,5 @@ jobs: mypy --namespace-packages --explicit-package-bases pymatgen - name: ruff run: | - ruff pymatgen + ruff check pymatgen ruff format pymatgen diff --git a/pymatgen/analysis/diffusion/aimd/pathway.py b/pymatgen/analysis/diffusion/aimd/pathway.py index 950e145..9ddc21a 100644 --- a/pymatgen/analysis/diffusion/aimd/pathway.py +++ b/pymatgen/analysis/diffusion/aimd/pathway.py @@ -47,9 +47,7 @@ def __init__(self, structure, trajectories, interval=0.5, species=("Li", "Na")): assert np.all(trajectories >= 0) assert np.all(trajectories <= 1) - indices = [ - j for j, site in enumerate(structure) if site.specie.symbol in species - ] + indices = [j for j, site in enumerate(structure) if site.specie.symbol in species] lattice = structure.lattice frac_interval = [interval / length for length in lattice.abc] nsteps = len(trajectories) @@ -83,35 +81,20 @@ def __init__(self, structure, trajectories, interval=0.5, species=("Li", "Na")): for i in range(3): next_i[i] = corner_i[i] + 1 if corner_i[i] < lens[i] - 1 else 0 - agrid = ( - np.array([corner_i[0], next_i[0]])[:, None] - * np.array([1, 0, 0])[None, :] - ) - bgrid = ( - np.array([corner_i[1], next_i[1]])[:, None] - * np.array([0, 1, 0])[None, :] - ) - cgrid = ( - np.array([corner_i[2], next_i[2]])[:, None] - * np.array([0, 0, 1])[None, :] - ) + agrid = np.array([corner_i[0], next_i[0]])[:, None] * np.array([1, 0, 0])[None, :] + bgrid = np.array([corner_i[1], next_i[1]])[:, None] * np.array([0, 1, 0])[None, :] + cgrid = np.array([corner_i[2], next_i[2]])[:, None] * np.array([0, 0, 1])[None, :] - grid_indices = ( - agrid[:, None, None] + bgrid[None, :, None] + cgrid[None, None, :] - ) + grid_indices = agrid[:, None, None] + bgrid[None, :, None] + cgrid[None, None, :] grid_indices = grid_indices.reshape(8, 3) mini_grid = [grid[indx[0], indx[1], indx[2]] for indx in grid_indices] dist_matrix = lattice.get_all_distances(mini_grid, fcoord) - indx = np.where(dist_matrix == np.min(dist_matrix, axis=0)[None, :])[0][ - 0 - ] + indx = np.where(dist_matrix == np.min(dist_matrix, axis=0)[None, :])[0][0] # 3-index label mapping to single index min_indx = ( - grid_indices[indx][0] * len(rb) * len(rc) - + grid_indices[indx][1] * len(rc) - + grid_indices[indx][2] + grid_indices[indx][0] * len(rb) * len(rc) + grid_indices[indx][1] * len(rc) + grid_indices[indx][2] ) # make sure the index does not go out of bound. @@ -133,9 +116,7 @@ def __init__(self, structure, trajectories, interval=0.5, species=("Li", "Na")): self.stable_sites = None @classmethod - def from_diffusion_analyzer( - cls, diffusion_analyzer, interval=0.5, species=("Li", "Na") - ): + def from_diffusion_analyzer(cls, diffusion_analyzer, interval=0.5, species=("Li", "Na")): """ Create a ProbabilityDensityAnalysis from a diffusion_analyzer object. @@ -154,9 +135,7 @@ def from_diffusion_analyzer( trajectories = np.array(trajectories) - return ProbabilityDensityAnalysis( - structure, trajectories, interval=interval, species=species - ) + return ProbabilityDensityAnalysis(structure, trajectories, interval=interval, species=species) def generate_stable_sites(self, p_ratio=0.25, d_cutoff=1.0): """ @@ -319,14 +298,10 @@ def __init__(self, structure, coords_ref, trajectories, species=("Li", "Na")): trajectories = np.array(trajectories) count = Counter() - indices = [ - i for i, site in enumerate(structure) if site.specie.symbol in species - ] + indices = [i for i, site in enumerate(structure) if site.specie.symbol in species] for it in range(len(trajectories)): - dist_matrix = lattice.get_all_distances( - coords_ref, trajectories[it, indices, :] - ) + dist_matrix = lattice.get_all_distances(coords_ref, trajectories[it, indices, :]) labels = np.where(dist_matrix == np.min(dist_matrix, axis=0)[None, :])[0] count.update(labels) @@ -344,15 +319,11 @@ def __init__(self, structure, coords_ref, trajectories, species=("Li", "Na")): self.site_occ = site_occ def get_average_site_occupancy(self, indices): - """ - Get the average site occupancy over a subset of reference sites. - """ + """Get the average site occupancy over a subset of reference sites.""" return np.sum(self.site_occ[indices]) / len(indices) @classmethod - def from_diffusion_analyzer( - cls, coords_ref, diffusion_analyzer, species=("Li", "Na") - ): + def from_diffusion_analyzer(cls, coords_ref, diffusion_analyzer, species=("Li", "Na")): """ Create a SiteOccupancyAnalyzer object using a diffusion_analyzer object. diff --git a/pymatgen/analysis/diffusion/neb/full_path_mapper.py b/pymatgen/analysis/diffusion/neb/full_path_mapper.py index 22fe029..3c64559 100644 --- a/pymatgen/analysis/diffusion/neb/full_path_mapper.py +++ b/pymatgen/analysis/diffusion/neb/full_path_mapper.py @@ -1,8 +1,5 @@ -# Copyright (c) Materials Virtual Lab. -# Distributed under the terms of the BSD License. -""" -Migraiton Graph Analysis -""" +"""Migration Graph Analysis.""" + from __future__ import annotations __author__ = "Jimmy Shen" @@ -39,7 +36,8 @@ def generic_groupby(list_in: list, comp: Callable = operator.eq): """ - Group a list of unsortable objects + Group a list of unsortable objects. + Args: list_in: A list of generic objects comp: (Default value = operator.eq) The comparator @@ -112,22 +110,16 @@ def __init__( @property def only_sites(self) -> Structure: - """ - A structure that only contains the migrating species - """ + """A structure that only contains the migrating species.""" return self.m_graph.structure @property def host_structure(self) -> Structure: - """ - A structure that only contains the non-migrating species - """ + """A structure that only contains the non-migrating species.""" host_struct = self.structure.copy() rm_sites = set() for isite in self.only_sites: - neighbors_ = host_struct.get_neighbors_in_shell( - isite.coords, r=0.0, dr=0.05, include_index=True - ) + neighbors_ = host_struct.get_neighbors_in_shell(isite.coords, r=0.0, dr=0.05, include_index=True) if len(neighbors_) == 0: continue for n_ in neighbors_: @@ -137,9 +129,7 @@ def host_structure(self) -> Structure: @property def symm_structure(self) -> SymmetrizedStructure: - """ - The symmetrized structure with the present item's symprec value - """ + """The symmetrized structure with the present item's symprec value.""" a = SpacegroupAnalyzer(self.structure, symprec=self.symprec) sym_struct = a.get_symmetrized_structure() if not isinstance(sym_struct, SymmetrizedStructure): @@ -148,9 +138,7 @@ def symm_structure(self) -> SymmetrizedStructure: @property def unique_hops(self): - """ - The unique hops dictionary keyed by the hop label - """ + """The unique hops dictionary keyed by the hop label.""" # reversed so that the first instance represents the group of distinct hops ihop_data = list(reversed(list(self.m_graph.graph.edges(data=True)))) for u, v, d in ihop_data: @@ -160,13 +148,12 @@ def unique_hops(self): return {d["hop_label"]: d for u, v, d in ihop_data} @classmethod - def with_base_structure( - cls, base_structure: Structure, m_graph: StructureGraph, **kwargs - ) -> MigrationGraph: + def with_base_structure(cls, base_structure: Structure, m_graph: StructureGraph, **kwargs) -> MigrationGraph: """ Args: - base_structure: base framework structure that does not contain any - migrating sites. + base_structure: base framework structure that does not contain any migrating sites. + m_graph: The StructureGraph object that defines the migration network. + **kwargs: Passthrough for kwargs. Returns: A constructed MigrationGraph object @@ -186,6 +173,7 @@ def with_local_env_strategy( structure: Input structure that contains all sites. migrating_specie: The specie that migrates. E.g. "Li". nn: The specific local environment object used to connect the migrating ion sites. + **kwargs: Passthrough for kwargs. Returns: A constructed MigrationGraph object @@ -202,9 +190,12 @@ def with_distance( Using a specific nn strategy to get the connectivity graph between all the migrating ion sites. Args: + structure: Input structure that contains all sites. + migrating_specie: The specie that migrates. E.g. "Li". max_distance: Maximum length of NEB path in the unit of Angstrom. Defaults to None, which means you are setting the value to the min cutoff until finding 1D or >1D percolating paths. + **kwargs: Passthrough for kwargs. Returns: A constructed MigrationGraph object @@ -229,6 +220,7 @@ def get_structure_from_entries( Args: entries: list of entries, must contain a mixture of inserted and empty structures. migrating_ion_entry: The metallic phase of the working ion, used to calculate insertion energies. + **kwargs: Passthrough for kwargs. Additional Kwargs: symprec: symmetry parameter for SpacegroupAnalyzer @@ -296,7 +288,7 @@ def _get_pos_and_migration_hop(self, u, v, w): Args: u (int): index of initial node v (int): index of final node - w (int): index for multiple edges that share the same two nodes + w (int): index for multiple edges that share the same two nodes. """ edge = self.m_graph.graph[u][v][w] i_site = self.only_sites.sites[u] @@ -311,26 +303,17 @@ def _get_pos_and_migration_hop(self, u, v, w): edge["ipos_cart"] = np.dot(i_site.frac_coords, self.only_sites.lattice.matrix) edge["epos_cart"] = np.dot(e_site.frac_coords, self.only_sites.lattice.matrix) - edge["hop"] = MigrationHop( - i_site, e_site, self.symm_structure, symprec=self.symprec - ) + edge["hop"] = MigrationHop(i_site, e_site, self.symm_structure, symprec=self.symprec) def _populate_edges_with_migration_hops(self): - """ - Populate the edges with the data for the Migration Paths - """ + """Populate the edges with the data for the Migration Paths.""" list(starmap(self._get_pos_and_migration_hop, self.m_graph.graph.edges)) def _group_and_label_hops(self): - """ - Group the MigrationHop objects together and label all the symmetrically equlivaelnt hops with the same label - """ + """Group the MigrationHop objects together and label all the symmetrically equlivaelnt hops with the same label.""" hops = list(nx.get_edge_attributes(self.m_graph.graph, "hop").items()) labs = generic_groupby(hops, comp=lambda x, y: x[1] == y[1]) - new_attr = { - g_index: {"hop_label": labs[edge_index]} - for edge_index, (g_index, _) in enumerate(hops) - } + new_attr = {g_index: {"hop_label": labs[edge_index]} for edge_index, (g_index, _) in enumerate(hops)} nx.set_edge_attributes(self.m_graph.graph, new_attr) return new_attr @@ -346,24 +329,19 @@ def add_data_to_similar_edges( target_label: The edge uniqueness label are adding data data: The data to passed to the different edges m_hop: If the data is an array, and m_hop is set, it uses the reference migration path to - determine whether the data needs to be flipped so that 0-->1 is different from 1-->0 + determine whether the data needs to be flipped so that 0-->1 is different from 1-->0. """ for _u, _v, d in self.m_graph.graph.edges(data=True): if d["hop_label"] == target_label: d.update(data) # Try to override the data. - if ( - m_hop is not None - and not m_hop.symm_structure.spacegroup.are_symmetrically_equivalent( - [m_hop.isite], [d["hop"].isite] - ) + if m_hop is not None and not m_hop.symm_structure.spacegroup.are_symmetrically_equivalent( + [m_hop.isite], [d["hop"].isite] ): # "The data going to this edge needs to be flipped" for k in data: if isinstance(data[k], (np.ndarray, np.generic)): - raise Warning( - "The data provided will only be flipped if it a list" - ) + raise Warning("The data provided will only be flipped if it a list") if not isinstance(data[k], list): continue d[k] = d[k][::-1] # flip the data in the array @@ -373,7 +351,7 @@ def assign_cost_to_graph(self, cost_keys=None): Read the data dict on each add and populate a cost key Args: cost_keys: a list of keys for data on each edge. - The SC Graph is decorated with a "cost" key that is the product of the different keys here + The SC Graph is decorated with a "cost" key that is the product of the different keys here. """ if cost_keys is None: cost_keys = ["hop_distance"] @@ -401,9 +379,7 @@ def get_path(self, max_val=100000, flip_hops=True): Each dict contains the information of a hop """ if len(self.unique_hops) != len(self.unique_hops): - logger.error( - f"There are {len(self.unique_hops)} SC hops but {len(self.unique_hops)} UC hops in {self}" - ) + logger.error(f"There are {len(self.unique_hops)} SC hops but {len(self.unique_hops)} UC hops in {self}") # for u, v, k, d in self.m_graph.graph.edges(data=True, keys=True): for u in self.m_graph.graph.nodes(): @@ -417,9 +393,7 @@ def get_path(self, max_val=100000, flip_hops=True): for tmp_u, tmp_v, tmp_k in cut_edges: path_graph.remove_edge(tmp_u, tmp_v, key=tmp_k) # populate the entire graph with multiple images - best_ans, path_parent = periodic_dijkstra( - path_graph, sources={u}, weight="cost", max_image=2 - ) + best_ans, path_parent = periodic_dijkstra(path_graph, sources={u}, weight="cost", max_image=2) # find a way to a u site that is not in the (0,0,0) image all_paths = [] for idx, jimage in path_parent: @@ -442,9 +416,7 @@ def get_path(self, max_val=100000, flip_hops=True): # displacement +/- (jimage1 - jimage2) is present on of of the edges # Note: there should only ever be one valid to_jimage for a u->v pair i1_, i2_ = sorted((idx1, idx2)) - all_edge_data = [ - *path_graph.get_edge_data(i1_, i2_, default={}).items() - ] + all_edge_data = [*path_graph.get_edge_data(i1_, i2_, default={}).items()] image_diff = np.subtract(jimage2, jimage1) found_ = 0 for _k, tmp_d in all_edge_data: @@ -460,7 +432,13 @@ def get_path(self, max_val=100000, flip_hops=True): def get_summary_dict(self, added_keys: list[str] | None = None) -> dict: """ - Dictionary format, for saving to database + Dictionary format, for saving to database. + + Args: + added_keys: a list of keys for data on each edge. + + Returns: + Dict. """ hops = [] keys = ["hop_label", "to_jimage", "ipos", "epos", "ipos_cart", "epos_cart"] @@ -479,8 +457,8 @@ def get_keys(d): unique_hops = [] for d in self.unique_hops.values(): - new_hop["iindex"] = d["iindex"] - new_hop["eindex"] = d["eindex"] + new_hop["iindex"] = d["iindex"] # type: ignore + new_hop["eindex"] = d["eindex"] # type: ignore unique_hops.append(get_keys(d)) unique_hops = sorted(unique_hops, key=lambda x: x["hop_label"]) @@ -495,9 +473,7 @@ def get_keys(d): class ChargeBarrierGraph(MigrationGraph): - """ - A Migration graph with additional charge density analysis on the charge density of the host material - """ + """A Migration graph with additional charge density analysis on the charge density of the host material.""" def __init__( self, @@ -512,11 +488,12 @@ def __init__( The graph is constructed using the structure, and cost values are assigned based on charge density analysis. Args: + structure (Structure): Input structure. + m_graph (StructureGraph): Input structure graph. potential_field: Input VolumetricData object that describes the field does not have to contains all the metastable sites. - migrating_specie (Specie-like): The specie that migrates. E.g., - "Li". - symprec (float): Symmetry precision to determine equivalence. + potential_data_key (str): Key for potential data. + **kwargs: Passthru for kwargs. """ self.potential_field = potential_field self.potential_data_key = potential_data_key @@ -524,17 +501,11 @@ def __init__( self._setup_grids() def _setup_grids(self): - """Populate the internal variables used for defining the grid points in the charge density analysis""" + """Populate the internal variables used for defining the grid points in the charge density analysis.""" # set up the grid - aa = np.linspace( - 0, 1, len(self.potential_field.get_axis_grid(0)), endpoint=False - ) - bb = np.linspace( - 0, 1, len(self.potential_field.get_axis_grid(1)), endpoint=False - ) - cc = np.linspace( - 0, 1, len(self.potential_field.get_axis_grid(2)), endpoint=False - ) + aa = np.linspace(0, 1, len(self.potential_field.get_axis_grid(0)), endpoint=False) + bb = np.linspace(0, 1, len(self.potential_field.get_axis_grid(1)), endpoint=False) + cc = np.linspace(0, 1, len(self.potential_field.get_axis_grid(2)), endpoint=False) # move the grid points to the center aa, bb, dd = map(_shift_grid, [aa, bb, cc]) @@ -553,15 +524,9 @@ def _setup_grids(self): def _dist_mat(self, pos_frac): # return a matrix that contains the distances to pos_frac - aa = np.linspace( - 0, 1, len(self.potential_field.get_axis_grid(0)), endpoint=False - ) - bb = np.linspace( - 0, 1, len(self.potential_field.get_axis_grid(1)), endpoint=False - ) - cc = np.linspace( - 0, 1, len(self.potential_field.get_axis_grid(2)), endpoint=False - ) + aa = np.linspace(0, 1, len(self.potential_field.get_axis_grid(0)), endpoint=False) + bb = np.linspace(0, 1, len(self.potential_field.get_axis_grid(1)), endpoint=False) + cc = np.linspace(0, 1, len(self.potential_field.get_axis_grid(2)), endpoint=False) aa, bb, cc = map(_shift_grid, [aa, bb, cc]) AA, BB, CC = np.meshgrid(aa, bb, cc, indexing="ij") dist_from_pos = self.potential_field.structure.lattice.get_all_distances( @@ -581,15 +546,9 @@ def _get_pathfinder_from_hop(self, migration_hop: MigrationHop, n_images=20): mid_struct = self.potential_field.structure.copy() # the moving ion is always inserted on the zero index - start_struct.insert( - 0, migration_hop.isite.species_string, ipos, properties=dict(magmom=0) - ) - end_struct.insert( - 0, migration_hop.isite.species_string, epos, properties=dict(magmom=0) - ) - mid_struct.insert( - 0, migration_hop.isite.species_string, mpos, properties=dict(magmom=0) - ) + start_struct.insert(0, migration_hop.isite.species_string, ipos, properties=dict(magmom=0)) + end_struct.insert(0, migration_hop.isite.species_string, epos, properties=dict(magmom=0)) + mid_struct.insert(0, migration_hop.isite.species_string, mpos, properties=dict(magmom=0)) chgpot = ChgcarPotential(self.potential_field, normalize=False) return NEBPathfinder( @@ -601,9 +560,7 @@ def _get_pathfinder_from_hop(self, migration_hop: MigrationHop, n_images=20): mid_struct=mid_struct, ) - def _get_avg_chg_at_max( - self, migration_hop, radius=None, chg_along_path=False, output_positions=False - ): + def _get_avg_chg_at_max(self, migration_hop, radius=None, chg_along_path=False, output_positions=False): """Obtain the maximum average charge along the path Args: migration_hop (MigrationHop): MigrationPath object that represents a given hop @@ -619,9 +576,8 @@ def _get_avg_chg_at_max( Returns: [float]: maximum of the charge density, (optional: entire list of charge density) """ - if radius is None: - rr = self._tube_radius - if rr <= 0: + rr = radius or self._tube_radius + if rr <= 0: # type: ignore raise ValueError("The integration radius must be positive.") npf = self._get_pathfinder_from_hop(migration_hop) @@ -631,9 +587,7 @@ def _get_avg_chg_at_max( for ict in centers: dist_mat = self._dist_mat(ict) mask = dist_mat < rr - vol_sphere = self.potential_field.structure.volume * ( - mask.sum() / self.potential_field.ngridpts - ) + vol_sphere = self.potential_field.structure.volume * (mask.sum() / self.potential_field.ngridpts) avg_chg.append( np.sum(self.potential_field.data[self.potential_data_key] * mask) / self.potential_field.ngridpts @@ -651,7 +605,7 @@ def _get_chg_between_sites_tube(self, migration_hop, mask_file_seedname=None): Args: migration_hop: MigrationHop object that represents a given hop mask_file_seedname(string): seed name for output of the migration path masks (for debugging and - visualization) (Default value = None) + visualization) (Default value = None). Returns: float: The total charge density in a tube that connects two sites of a given edges of the graph @@ -659,9 +613,7 @@ def _get_chg_between_sites_tube(self, migration_hop, mask_file_seedname=None): try: _ = self._tube_radius except NameError: - logger.warning( - "The radius of the tubes for charge analysis need to be defined first." - ) + logger.warning("The radius of the tubes for charge analysis need to be defined first.") ipos = migration_hop.isite.frac_coords epos = migration_hop.esite.frac_coords @@ -669,15 +621,10 @@ def _get_chg_between_sites_tube(self, migration_hop, mask_file_seedname=None): cart_epos = np.dot(epos, self.potential_field.structure.lattice.matrix) pbc_mask = np.zeros(self._uc_grid_shape, dtype=bool).flatten() for img in self._images: - grid_pos = np.dot( - self._fcoords + img, self.potential_field.structure.lattice.matrix - ) - proj_on_line = np.dot(grid_pos - cart_ipos, cart_epos - cart_ipos) / ( - np.linalg.norm(cart_epos - cart_ipos) - ) + grid_pos = np.dot(self._fcoords + img, self.potential_field.structure.lattice.matrix) + proj_on_line = np.dot(grid_pos - cart_ipos, cart_epos - cart_ipos) / (np.linalg.norm(cart_epos - cart_ipos)) dist_to_line = np.linalg.norm( - np.cross(grid_pos - cart_ipos, cart_epos - cart_ipos) - / (np.linalg.norm(cart_epos - cart_ipos)), + np.cross(grid_pos - cart_ipos, cart_epos - cart_ipos) / (np.linalg.norm(cart_epos - cart_ipos)), axis=-1, ) @@ -725,8 +672,7 @@ def populate_edges_with_chg_density_info(self, tube_radius=1): v["hop"], chg_along_path=True, output_positions=True ) images = [ - {"position": ifrac, "average_charge": ichg} - for ifrac, ichg in zip(frac_coords_list, avg_chg_list) + {"position": ifrac, "average_charge": ichg} for ifrac, ichg in zip(frac_coords_list, avg_chg_list) ] v.update( dict( @@ -741,7 +687,7 @@ def get_least_chg_path(self): """ obtain an intercolating pathway through the material that has the least amount of charge Returns: - list of hops + list of hops. """ min_chg = 100000000 min_path = [] @@ -756,9 +702,7 @@ def get_least_chg_path(self): return min_path def get_summary_dict(self, add_keys: list[str] | None = None): - """ - Dictionary format, for saving to database - """ + """Dictionary format, for saving to database.""" a_keys = ["max_avg_chg", "chg_total"] if add_keys is not None: a_keys += add_keys @@ -766,9 +710,7 @@ def get_summary_dict(self, add_keys: list[str] | None = None): # Utility functions -def get_only_sites_from_structure( - structure: Structure, migrating_specie: str -) -> Structure: +def get_only_sites_from_structure(structure: Structure, migrating_specie: str) -> Structure: """ Get a copy of the structure with only the migrating sites. @@ -791,15 +733,13 @@ def _shift_grid(vv): """ Move the grid points by half a step so that they sit in the center Args: - vv: equally space grid points in 1-D + vv: equally space grid points in 1-D. """ step = vv[1] - vv[0] return vv + step / 2.0 -def get_hop_site_sequence( - hop_list: list[dict], start_u: int | str, key: str | None = None -) -> list: +def get_hop_site_sequence(hop_list: list[dict], start_u: int | str, key: str | None = None) -> list: """ Read in a list of hop dictionaries and print the sequence of sites (and relevant property values if any). @@ -813,11 +753,7 @@ def get_hop_site_sequence( """ hops = iter(hop_list) ihop = next(hops) - site_seq = ( - [ihop["eindex"], ihop["iindex"]] - if ihop["eindex"] == start_u - else [ihop["iindex"], ihop["eindex"]] - ) + site_seq = [ihop["eindex"], ihop["iindex"]] if ihop["eindex"] == start_u else [ihop["iindex"], ihop["eindex"]] for ihop in hops: if ihop["iindex"] == site_seq[-1]: @@ -909,9 +845,7 @@ def order_path(hop_list: list[dict], start_u: int | str) -> list[dict]: def almost(a, b): - """ - return true if the values are almost equal - """ + """Return true if the values are almost equal.""" SMALL_VAL = 1e-4 try: return all(almost(i, j) for i, j in zip(list(a), list(b))) @@ -924,7 +858,7 @@ def almost(a, b): def check_uc_hop(sc_hop, uc_hop): """ See if hop in the 2X2X2 supercell and a unit cell hop - are equivalent under lattice translation + are equivalent under lattice translation. Args: sc_hop: MigrationHop object form pymatgen-diffusion. diff --git a/pymatgen/analysis/diffusion/neb/io.py b/pymatgen/analysis/diffusion/neb/io.py index e686940..5550892 100644 --- a/pymatgen/analysis/diffusion/neb/io.py +++ b/pymatgen/analysis/diffusion/neb/io.py @@ -1,9 +1,8 @@ # Copyright (c) Materials Virtual Lab. # Distributed under the terms of the BSD License. -""" -Generate input fields for NEB calculations. -""" +"""Generate input fields for NEB calculations.""" + from __future__ import annotations import copy @@ -15,15 +14,13 @@ class MVLCINEBEndPointSet(MITRelaxSet): - """ - Class for writing NEB end points relaxation inputs. - """ + """Class for writing NEB end points relaxation inputs.""" def __init__(self, structure, **kwargs): r""" Args: structure: Structure - \*\*kwargs: Keyword args supported by VaspInputSets. + **kwargs: Keyword args supported by VaspInputSets. """ user_incar_settings = kwargs.get("user_incar_settings", {}) defaults = { @@ -48,14 +45,14 @@ class MVLCINEBSet(MITNEBSet): """ MAVRL-tested settings for CI-NEB calculations. Note that these parameters requires the VTST modification of VASP from the Henkelman group. See - http://theory.cm.utexas.edu/vtsttools/ + http://theory.cm.utexas.edu/vtsttools/. """ def __init__(self, structures, **kwargs): r""" Args: - structure: Structure - \*\*kwargs: Keyword args supported by VaspInputSets. + structures: Input structures. + **kwargs: Keyword args supported by VaspInputSets. """ user_incar_settings = kwargs.get("user_incar_settings", {}) diff --git a/pymatgen/analysis/diffusion/neb/pathfinder.py b/pymatgen/analysis/diffusion/neb/pathfinder.py index 1119651..c6a64ab 100644 --- a/pymatgen/analysis/diffusion/neb/pathfinder.py +++ b/pymatgen/analysis/diffusion/neb/pathfinder.py @@ -1,9 +1,8 @@ # Copyright (c) Materials Virtual Lab. # Distributed under the terms of the BSD License. -""" -Algorithms for NEB migration path analysis. -""" +"""Algorithms for NEB migration path analysis.""" + from __future__ import annotations import itertools @@ -86,9 +85,7 @@ def __init__(self, structures): if ni not in [0, nimages + 1]: for j in range(i + 1, natoms): - img = latt.get_distance_and_image( - frac_coords, structures[ni][j].frac_coords - )[1] + img = latt.get_distance_and_image(frac_coords, structures[ni][j].frac_coords)[1] translations[ni - 1, i, j] = latt.get_cartesian_coords(img) translations[ni - 1, j, i] = -latt.get_cartesian_coords(img) @@ -143,9 +140,7 @@ def run( indices = list(range(len(self.structures[0]))) else: species = [get_el_sp(sp) for sp in species] - indices = [ - i for i, site in enumerate(self.structures[0]) if site.specie in species - ] + indices = [i for i, site in enumerate(self.structures[0]) if site.specie in species] if len(indices) == 0: raise ValueError("The given species are not in the system!") @@ -155,15 +150,11 @@ def run( # Get the sets of objective functions, true and total force # matrices. funcs, true_forces = self._get_funcs_and_forces(coords) - tot_forces = self._get_total_forces( - coords, true_forces, spring_const=spring_const - ) + tot_forces = self._get_total_forces(coords, true_forces, spring_const=spring_const) # Each atom is allowed to move up to max_disp disp_mat = step_size * tot_forces[:, indices, :] - disp_mat = np.where( - np.abs(disp_mat) > max_disp, np.sign(disp_mat) * max_disp, disp_mat - ) + disp_mat = np.where(np.abs(disp_mat) > max_disp, np.sign(disp_mat) * max_disp, disp_mat) coords[1 : (self.nimages + 1), indices] += disp_mat max_force = np.abs(tot_forces[:, indices, :]).max() @@ -175,9 +166,7 @@ def run( old_funcs = funcs else: - warnings.warn( - "Maximum iteration number is reached without convergence!", UserWarning - ) + warnings.warn("Maximum iteration number is reached without convergence!", UserWarning) for ni in range(self.nimages): # generate the improved image structure @@ -202,7 +191,11 @@ def run( @classmethod def from_endpoints( - cls, endpoints, nimages=5, sort_tol=1.0, interpolate_lattices=False + cls, + endpoints, + nimages: int = 5, + sort_tol: float = 1.0, + interpolate_lattices: bool = False, ): """ A class method that starts with end-point structures instead. The @@ -215,6 +208,7 @@ def from_endpoints( sort_tol (float): Distance tolerance (in Angstrom) used to match the atomic indices between start and end structures. Need to increase the value in some cases. + interpolate_lattices (bool): Whether to interpolate lattices between the start and end structures. """ try: images = endpoints[0].interpolate( @@ -226,8 +220,7 @@ def from_endpoints( except Exception as e: if "Unable to reliably match structures " in str(e): warnings.warn( - "Auto sorting is turned off because it is unable" - " to match the end-point structures!", + "Auto sorting is turned off because it is unable" " to match the end-point structures!", UserWarning, ) images = endpoints[0].interpolate( @@ -244,7 +237,7 @@ def from_endpoints( def _get_funcs_and_forces(self, x): """ Calculate the set of objective functions as well as their gradients, - i.e. "effective true forces" + i.e. "effective true forces". """ funcs = [] funcs_prime = [] @@ -257,11 +250,7 @@ def _get_funcs_and_forces(self, x): vec = [x[ni + 1, i] - x[ni + 1] - trans[ni, i] for i in range(natoms)] trial_dist = np.linalg.norm(vec, axis=2) - aux = ( - (trial_dist - target_dists[ni]) - * weights[ni] - / (trial_dist + np.eye(natoms, dtype=np.float64)) - ) + aux = (trial_dist - target_dists[ni]) * weights[ni] / (trial_dist + np.eye(natoms, dtype=np.float64)) # Objective function func = np.sum((trial_dist - target_dists[ni]) ** 2 * weights[ni]) @@ -303,24 +292,18 @@ def _get_total_forces(self, x, true_forces, spring_const): tangent = self.get_unit_vector(tangent) # Spring force - spring_force = ( - spring_const * (np.linalg.norm(vec1) - np.linalg.norm(vec2)) * tangent - ) + spring_force = spring_const * (np.linalg.norm(vec1) - np.linalg.norm(vec2)) * tangent # Total force flat_ft = true_forces[ni - 1].copy().flatten() - total_force = true_forces[ni - 1] + ( - spring_force - np.dot(flat_ft, tangent) * tangent - ).reshape(natoms, 3) + total_force = true_forces[ni - 1] + (spring_force - np.dot(flat_ft, tangent) * tangent).reshape(natoms, 3) total_forces.append(total_force) return np.array(total_forces) class MigrationHop(MSONable): - """ - A convenience container representing a migration path. - """ + """A convenience container representing a migration path.""" def __init__( self, @@ -337,7 +320,7 @@ def __init__( symm_structure: SymmetrizedStructure host_symm_struct: SymmetrizedStructure of the host structure, used to for its spacegroup - symprec: used to determine equivalence + symprec: used to determine equivalence. """ self.isite = isite self.esite = esite @@ -346,9 +329,7 @@ def __init__( self.symm_structure = symm_structure self.host_symm_struct = host_symm_struct self.symprec = symprec - self.msite = PeriodicSite( - esite.specie, (isite.frac_coords + esite.frac_coords) / 2, esite.lattice - ) + self.msite = PeriodicSite(esite.specie, (isite.frac_coords + esite.frac_coords) / 2, esite.lattice) sg = ( self.host_symm_struct.spacegroup # type: ignore[union-attr] @@ -382,13 +363,9 @@ def __init__( break if self.iindex is None: - raise RuntimeError( - f"No symmetrically equivalent site was found for {isite}" - ) + raise RuntimeError(f"No symmetrically equivalent site was found for {isite}") if self.eindex is None: - raise RuntimeError( - f"No symmetrically equivalent site was found for {esite}" - ) + raise RuntimeError(f"No symmetrically equivalent site was found for {esite}") def __repr__(self): ifc = self.isite.frac_coords @@ -451,27 +428,19 @@ def get_structures(self, nimages=5, vac_mode=True, idpp=False, **idpp_kwargs): idpp (bool): Defaults to False. If True, the generated structures will be run through the IDPPSolver to generate a better guess for the minimum energy path. - \*\*idpp_kwargs: Passthrough kwargs for the IDPPSolver.run. + **idpp_kwargs: Passthrough kwargs for the IDPPSolver.run. Returns: [Structure] Note that the first site of each structure is always the migrating ion. This makes it easier to perform subsequent analysis. """ - migrating_specie_sites, other_sites = self._split_migrating_and_other_sites( - vac_mode - ) + migrating_specie_sites, other_sites = self._split_migrating_and_other_sites(vac_mode) - start_structure = Structure.from_sites( - [self.isite, *migrating_specie_sites, *other_sites] - ) - end_structure = Structure.from_sites( - [self.esite, *migrating_specie_sites, *other_sites] - ) + start_structure = Structure.from_sites([self.isite, *migrating_specie_sites, *other_sites]) + end_structure = Structure.from_sites([self.esite, *migrating_specie_sites, *other_sites]) - structures = start_structure.interpolate( - end_structure, nimages=nimages + 1, pbc=False - ) + structures = start_structure.interpolate(end_structure, nimages=nimages + 1, pbc=False) if idpp: solver = IDPPSolver(structures) @@ -486,10 +455,7 @@ def _split_migrating_and_other_sites(self, vac_mode): if site.specie != self.isite.specie: other_sites.append(site) else: - if ( - self.isite.distance(site) <= 1e-8 - or self.esite.distance(site) <= 1e-8 - ): + if self.isite.distance(site) <= 1e-8 or self.esite.distance(site) <= 1e-8: migrating_specie_sites.append(site) continue @@ -523,9 +489,7 @@ def get_sc_structures( If in vacancy mode, the base structure is the fully intercalated structure """ - migrating_specie_sites, other_sites = self._split_migrating_and_other_sites( - vac_mode - ) + migrating_specie_sites, other_sites = self._split_migrating_and_other_sites(vac_mode) if vac_mode: base_struct = Structure.from_sites(other_sites + migrating_specie_sites) else: @@ -552,7 +516,7 @@ def write_path(self, fname, **kwargs): Args: fname (str): File name. - \*\*kwargs: Kwargs supported by NEBPath.get_structures. + **kwargs: Kwargs supported by NEBPath.get_structures. """ sites = [] for st in self.get_structures(**kwargs): @@ -643,9 +607,7 @@ def get_paths(self): for sites in self.symm_structure.equivalent_sites: if sites[0].specie == self.migrating_specie: site0 = sites[0] - for nn in self.symm_structure.get_neighbors( - site0, r=round(self.max_path_length, 3) + 0.01 - ): + for nn in self.symm_structure.get_neighbors(site0, r=round(self.max_path_length, 3) + 0.01): if nn.specie == self.migrating_specie: path = MigrationHop(site0, nn, self.symm_structure) paths.add(path) @@ -661,18 +623,16 @@ def write_all_paths(self, fname, nimages=5, **kwargs): Args: fname (str): Filename nimages (int): Number of images per path. - \*\*kwargs: Passthrough kwargs to path.get_structures. + **kwargs: Passthrough kwargs to path.get_structures. """ sites = [] for p in self.get_paths(): - structures = p.get_structures( - nimages=nimages, species=[self.migrating_specie], **kwargs - ) + structures = p.get_structures(nimages=nimages, species=[self.migrating_specie], **kwargs) sites.append(structures[0][0]) sites.append(structures[-1][0]) for s in structures[1:-1]: sites.append(PeriodicSite("H", s[0].frac_coords, s.lattice)) - sites.extend(structures[0].sites[1:]) + sites.extend(structures[0].sites[1:]) # type: ignore Structure.from_sites(sites).to(filename=fname) @@ -690,12 +650,11 @@ class NEBPathfinder: Ceder, The Journal of Chemical Physics 145 (7), 074112 """ - def __init__( - self, start_struct, end_struct, relax_sites, v, n_images=20, mid_struct=None - ): + def __init__(self, start_struct, end_struct, relax_sites, v, n_images=20, mid_struct=None): """ Args: - start_struct, end_struct: Endpoint structures to interpolate + start_struct: Starting structure + end_struct: End structure to interpolate relax_sites: List of site indices whose interpolation paths should be relaxed v: Static potential field to use for the elastic band relaxation @@ -727,18 +686,12 @@ def interpolate(self): # to make arithmetic easier we will do the interpolation in two parts with n # images each then just take every other image at the end, this results in # exactly n images - images_0 = self.__s1.interpolate( - self.__mid, nimages=self.__n_images, interpolate_lattices=False - )[:-1] - images_1 = self.__mid.interpolate( - self.__s2, nimages=self.__n_images, interpolate_lattices=False - ) + images_0 = self.__s1.interpolate(self.__mid, nimages=self.__n_images, interpolate_lattices=False)[:-1] + images_1 = self.__mid.interpolate(self.__s2, nimages=self.__n_images, interpolate_lattices=False) images = images_0 + images_1 images = images[::2] else: - images = self.__s1.interpolate( - self.__s2, nimages=self.__n_images, interpolate_lattices=False - ) + images = self.__s1.interpolate(self.__s2, nimages=self.__n_images, interpolate_lattices=False) for site_i in self.__relax_sites: start_f = images[0].sites[site_i].frac_coords end_f = images[-1].sites[site_i].frac_coords @@ -757,8 +710,7 @@ def interpolate(self): for image_i, image in enumerate(images): image.translate_sites( site_i, - NEBPathfinder.__d2f(path[image_i], self.__v) - - image.sites[site_i].frac_coords, + NEBPathfinder.__d2f(path[image_i], self.__v) - image.sites[site_i].frac_coords, frac_coords=True, to_unit_cell=True, ) @@ -776,7 +728,9 @@ def plot_images(self, outfile): """ Generates a POSCAR with the calculated diffusion path with respect to the first endpoint. - :param outfile: Output file for the POSCAR. + + Args: + outfile: Output file for the POSCAR. """ sum_struct = self.__images[0].sites for image in self.__images: @@ -847,11 +801,7 @@ def string_relax( # logger.debug(f"Getting path from {start} to {end} (coords wrt V grid)") # Set parameters - dr = ( - np.array([1 / V.shape[0], 1 / V.shape[1], 1 / V.shape[2]]) - if not dr - else np.array(dr, dtype=float) - ) + dr = np.array([1 / V.shape[0], 1 / V.shape[1], 1 / V.shape[2]]) if not dr else np.array(dr, dtype=float) keff = k * dr * n_images h0 = h @@ -880,23 +830,16 @@ def string_relax( # Evolve string for step in range(max_iter): # Gradually decay step size to prevent oscillations - h = ( - h0 * np.exp(-2 * (step - min_iter) / max_iter) - if step > min_iter - else h0 - ) + h = h0 * np.exp(-2 * (step - min_iter) / max_iter) if step > min_iter else h0 # Calculate forces acting on string d = V.shape s0 = s.copy() # store copy for endpoint fixing below (fixes GH 2732) edV = np.array( [ [ - dV[0][int(pt[0]) % d[0]][int(pt[1]) % d[1]][int(pt[2]) % d[2]] - / dr[0], - dV[1][int(pt[0]) % d[0]][int(pt[1]) % d[1]][int(pt[2]) % d[2]] - / dr[0], - dV[2][int(pt[0]) % d[0]][int(pt[1]) % d[1]][int(pt[2]) % d[2]] - / dr[0], + dV[0][int(pt[0]) % d[0]][int(pt[1]) % d[1]][int(pt[2]) % d[2]] / dr[0], + dV[1][int(pt[0]) % d[0]][int(pt[1]) % d[1]][int(pt[2]) % d[2]] / dr[0], + dV[2][int(pt[0]) % d[0]][int(pt[1]) % d[1]][int(pt[2]) % d[2]] / dr[0], ] for pt in s ] @@ -910,16 +853,8 @@ def string_relax( ds_plus[0] = ds_plus[0] - ds_plus[0] ds_minus[-1] = ds_minus[-1] - ds_minus[-1] Fpot = edV - Fel = ( - keff - * (la.norm(ds_plus) - la.norm(ds0_plus)) - * (ds_plus / la.norm(ds_plus)) - ) - Fel += ( - keff - * (la.norm(ds_minus) - la.norm(ds0_minus)) - * (ds_minus / la.norm(ds_minus)) - ) + Fel = keff * (la.norm(ds_plus) - la.norm(ds0_plus)) * (ds_plus / la.norm(ds_plus)) + Fel += keff * (la.norm(ds_minus) - la.norm(ds0_minus)) * (ds_minus / la.norm(ds_minus)) s -= h * (Fpot + Fel) # Fix endpoints @@ -937,10 +872,7 @@ def string_relax( tol = la.norm((s - s0) * dr) / n_images / h if tol > 1e10: - raise ValueError( - "Pathfinding failed, path diverged! Consider reducing h to avoid " - "divergence." - ) + raise ValueError("Pathfinding failed, path diverged! Consider reducing h to avoid " "divergence.") if step > min_iter and tol < max_tol: logger.debug(f"Converged at {step=}") @@ -1011,24 +943,17 @@ def rescale_field(self, new_dim): """ v_dim = self.__v.shape padded_v = np.lib.pad(self.__v, ((0, 1), (0, 1), (0, 1)), mode="wrap") - ogrid_list = np.array( - [ - list(c) - for c in list(np.ndindex(v_dim[0] + 1, v_dim[1] + 1, v_dim[2] + 1)) - ] - ) - v_ogrid = padded_v.reshape( - ((v_dim[0] + 1) * (v_dim[1] + 1) * (v_dim[2] + 1), -1) - ) + ogrid_list = np.array([list(c) for c in list(np.ndindex(v_dim[0] + 1, v_dim[1] + 1, v_dim[2] + 1))]) + v_ogrid = padded_v.reshape(((v_dim[0] + 1) * (v_dim[1] + 1) * (v_dim[2] + 1), -1)) ngrid_a, ngrid_b, ngrid_c = np.mgrid[ 0 : v_dim[0] : v_dim[0] / new_dim[0], 0 : v_dim[1] : v_dim[1] / new_dim[1], 0 : v_dim[2] : v_dim[2] / new_dim[2], ] - v_ngrid = scipy.interpolate.griddata( - ogrid_list, v_ogrid, (ngrid_a, ngrid_b, ngrid_c), method="linear" - ).reshape((new_dim[0], new_dim[1], new_dim[2])) + v_ngrid = scipy.interpolate.griddata(ogrid_list, v_ogrid, (ngrid_a, ngrid_b, ngrid_c), method="linear").reshape( + (new_dim[0], new_dim[1], new_dim[2]) + ) self.__v = v_ngrid def gaussian_smear(self, r): @@ -1065,9 +990,9 @@ def gaussian_smear(self, r): for g_b in np.arange(-2 * r_disc[1], 2 * r_disc[1] + 1, 1): for g_c in np.arange(-2 * r_disc[2], 2 * r_disc[2] + 1, 1): g = np.array([g_a / v_dim[0], g_b / v_dim[1], g_c / v_dim[2]]).T - gauss_dist[int(g_a + r_disc[0])][int(g_b + r_disc[1])][ - int(g_c + r_disc[2]) - ] = la.norm(np.dot(self.__s.lattice.matrix, g)) / r + gauss_dist[int(g_a + r_disc[0])][int(g_b + r_disc[1])][int(g_c + r_disc[2])] = ( + la.norm(np.dot(self.__s.lattice.matrix, g)) / r + ) gauss = scipy.stats.norm.pdf(gauss_dist) gauss = gauss / np.sum(gauss, dtype=float) padded_v = np.pad( @@ -1131,9 +1056,7 @@ def __add_gaussians(s, dim, r=1.5): for b_d in np.arange(0, dim[1], 1): for c_d in np.arange(0, dim[2], 1): coords_f = np.array([a_d / dim[0], b_d / dim[1], c_d / dim[2]]) - d_f = sorted( - s.get_sites_in_sphere(coords_f, s.lattice.a), key=lambda x: x[1] - )[0][1] + d_f = sorted(s.get_sites_in_sphere(coords_f, s.lattice.a), key=lambda x: x[1])[0][1] # logger.debug(d_f) gauss_dist[int(a_d)][int(b_d)][int(c_d)] = d_f / r return scipy.stats.norm.pdf(gauss_dist) diff --git a/pymatgen/analysis/diffusion/neb/periodic_dijkstra.py b/pymatgen/analysis/diffusion/neb/periodic_dijkstra.py index d5646d9..4bafd9e 100644 --- a/pymatgen/analysis/diffusion/neb/periodic_dijkstra.py +++ b/pymatgen/analysis/diffusion/neb/periodic_dijkstra.py @@ -1,8 +1,7 @@ # Copyright (c) Materials Virtual Lab. # Distributed under the terms of the BSD License. -""" -Dijkstra's path search on a graph where the nodes are on a periodic graph -""" +"""Dijkstra's path search on a graph where the nodes are on a periodic graph.""" + from __future__ import annotations __author__ = "Jimmy Shen" @@ -30,6 +29,9 @@ def _get_adjacency_with_images(G: Graph) -> dict: Note: the current implementation assumes that the original "to_jimage" value always corresponds to a an edge u -> v where u <= v. + Args: + G (pymatgen.analysis.graphs.StructureGraph): Structure graph. + Returns: dict: Nested dictionary with [start][end][edge_key][data_field] """ @@ -75,7 +77,7 @@ def periodic_dijkstra( G (Graph): The graph object with additional "to_jimage" fields to indicate edges across periodic images. sources (set): the index of the source node - target (int, optional): The index of of target node, if None populate all nodes. Defaults to None. + weight: the weight of the edges. max_image (int, optional): Defaults to 3. target_reached (callable, optional): A function of (site_index, jimage) used to check for stop iteration. This function is always called on the top of heap so it might miss the optimal path but @@ -134,7 +136,7 @@ def periodic_dijkstra_on_sgraph( Args: sgraph (Graph): The StructureGraph object used for path searching sources (set): the index of the source node - target (int, optional): The index of of target node, if None populate all nodes. Defaults to None. + weight: the weight of the edges. max_image (int, optional): Defaults to 3. target_reached (callable, optional): A function of (site_index, jimage) used to check for stop iteration. This function is always called on the top of heap so it might miss the optimal path but @@ -156,9 +158,7 @@ def periodic_dijkstra_on_sgraph( def get_optimal_pathway_rev(path_parent: dict, leaf_node: tuple): - """ - follow a leaf node all the way up to source. - """ + """Follow a leaf node all the way up to source.""" cur = leaf_node while cur in path_parent: yield cur diff --git a/pymatgen/analysis/diffusion/utils/__init__.py b/pymatgen/analysis/diffusion/utils/__init__.py index 88e19df..72c34ee 100644 --- a/pymatgen/analysis/diffusion/utils/__init__.py +++ b/pymatgen/analysis/diffusion/utils/__init__.py @@ -1,3 +1 @@ -""" -Utils Modules -""" +"""Utils Modules.""" diff --git a/pymatgen/analysis/diffusion/utils/edge_data_from_sc.py b/pymatgen/analysis/diffusion/utils/edge_data_from_sc.py index acbe135..b7c84f7 100644 --- a/pymatgen/analysis/diffusion/utils/edge_data_from_sc.py +++ b/pymatgen/analysis/diffusion/utils/edge_data_from_sc.py @@ -1,9 +1,8 @@ # Copyright (c) Materials Virtual Lab. # Distributed under the terms of the BSD License. -""" -Function to add edge data to MigrationGraph through 2 SC structures -""" +"""Function to add edge data to MigrationGraph through 2 SC structures.""" + from __future__ import annotations __author__ = "Haoming Li" @@ -39,12 +38,12 @@ def add_edge_data_from_sc( supercell structures. Args: + mg: MigrationGraph object. i_sc: Supercell structure containing working ion at initial position e_sc: Supercell structure containing working ion at ending position data_array: The data to be added to the edges key: Key of the edge attribute to be added - use_host_sg: Flag t whether or not to use the host structure's spacegroup to - initiate MigrationHop + use_host_sg: Flag whether to use the host structure's spacegroup to initiate MigrationHop Returns: None @@ -53,15 +52,11 @@ def add_edge_data_from_sc( i_wi = [x for x in i_sc.sites if x.species_string == wi] e_wi = [x for x in e_sc.sites if x.species_string == wi] if len(i_wi) != 1 or len(e_wi) != 1: - raise ValueError( - "The number of working ions in each supercell structure should be one" - ) + raise ValueError("The number of working ions in each supercell structure should be one") isite, esite = i_wi[0], e_wi[0] uhop_index, mh_from_sc = get_unique_hop(mg, i_sc, isite, esite, use_host_sg) add_dict = {key: data_array} - mg.add_data_to_similar_edges( - target_label=uhop_index, data=add_dict, m_hop=mh_from_sc - ) + mg.add_data_to_similar_edges(target_label=uhop_index, data=add_dict, m_hop=mh_from_sc) def get_uc_pos( @@ -71,7 +66,7 @@ def get_uc_pos( sc: Structure, sm: StructureMatcher, ) -> tuple[PeriodicSite, PeriodicSite, PeriodicSite]: - """Take positions in the supercell and transform into the unit cell positions + """Take positions in the supercell and transform into the unit cell positions. Args: isite: initial site in the SC @@ -85,10 +80,7 @@ def get_uc_pos( """ mapping = get_matched_structure_mapping(base=uc, inserted=sc, sm=sm) if mapping is None: - raise ValueError( - "Cannot obtain inverse mapping, consider lowering tolerances " - "in StructureMatcher" - ) + raise ValueError("Cannot obtain inverse mapping, consider lowering tolerances " "in StructureMatcher") sc_m, total_t = mapping sc_ipos = isite.frac_coords sc_ipos_t = sc_ipos - total_t @@ -123,9 +115,7 @@ def get_uc_pos( def _get_first_close_site(frac_coord, structure, stol=0.1): for site in structure.sites: - dist, image = structure.lattice.get_distance_and_image( - frac_coord, site.frac_coords - ) + dist, image = structure.lattice.get_distance_and_image(frac_coord, site.frac_coords) if dist < stol: return np.add(site.frac_coords, image) return frac_coord @@ -133,7 +123,7 @@ def _get_first_close_site(frac_coord, structure, stol=0.1): def mh_eq(mh1, mh2): """ - Allow for symmetric matching of MigrationPath objects with variable precession + Allow for symmetric matching of MigrationPath objects with variable precession. Args: mh1: MigrationHop object @@ -153,7 +143,7 @@ def get_unique_hop( esite: PeriodicSite, use_host_sg: bool = True, ) -> tuple[int, MigrationHop]: - """Get the unique hop label that correspond to two end positions in the SC + """Get the unique hop label that correspond to two end positions in the SC. Args: mg: Object containing the migration analysis @@ -166,16 +156,10 @@ def get_unique_hop( Returns: The index of the unique hop, the MigrationHop object transformed from the SC """ - sm = StructureMatcher( - ignored_species=[ - next(iter(mg.m_graph.graph.edges(data=True)))[2]["hop"].isite.specie.name - ] - ) + sm = StructureMatcher(ignored_species=[next(iter(mg.m_graph.graph.edges(data=True)))[2]["hop"].isite.specie.name]) uc_isite, uc_msite, uc_esite = get_uc_pos(isite, esite, mg.symm_structure, sc, sm) if use_host_sg: - base_ss = SpacegroupAnalyzer( - mg.host_structure, symprec=mg.symprec - ).get_symmetrized_structure() + base_ss = SpacegroupAnalyzer(mg.host_structure, symprec=mg.symprec).get_symmetrized_structure() mh_from_sc = MigrationHop( uc_isite, uc_esite, @@ -184,9 +168,7 @@ def get_unique_hop( symprec=mg.symprec, ) else: - mh_from_sc = MigrationHop( - uc_isite, uc_esite, symm_structure=mg.symm_structure, symprec=mg.symprec - ) + mh_from_sc = MigrationHop(uc_isite, uc_esite, symm_structure=mg.symm_structure, symprec=mg.symprec) result = [] for k, v in mg.unique_hops.items(): # tolerance may be changed here @@ -198,7 +180,5 @@ def get_unique_hop( if len(result) == 0: raise ValueError("No matches between UC and SC") - assert mg.symm_structure.spacegroup.are_symmetrically_equivalent( - [uc_msite], [mh_from_sc.msite] - ) + assert mg.symm_structure.spacegroup.are_symmetrically_equivalent([uc_msite], [mh_from_sc.msite]) return result[0], mh_from_sc diff --git a/pymatgen/analysis/diffusion/utils/maggma.py b/pymatgen/analysis/diffusion/utils/maggma.py index 62a0035..a49b305 100644 --- a/pymatgen/analysis/diffusion/utils/maggma.py +++ b/pymatgen/analysis/diffusion/utils/maggma.py @@ -3,8 +3,9 @@ """ Functions for querying Materials Project style MongoStores that contains cathode materials The functions are isolated from the rest of the package so -that the rest of the package will not depend on Maggma +that the rest of the package will not depend on Maggma. """ + from __future__ import annotations __author__ = "Jimmy Shen" @@ -42,7 +43,7 @@ def get_entries_from_dbs( material_store: Material documents one per each similar structure ( multiple tasks) migrating_ion: The name of the migrating ion - material_ids list with topotactic structures + material_id: Material id """ with structure_group_store as store: sg_doc = store.query_one({structure_group_store.key: material_id}) diff --git a/pymatgen/analysis/diffusion/utils/parse_entries.py b/pymatgen/analysis/diffusion/utils/parse_entries.py index dcbd779..e289ec7 100644 --- a/pymatgen/analysis/diffusion/utils/parse_entries.py +++ b/pymatgen/analysis/diffusion/utils/parse_entries.py @@ -1,9 +1,8 @@ # Copyright (c) Materials Virtual Lab. # Distributed under the terms of the BSD License. -""" -Functions for combining many ComputedEntry objects into MigrationGraph objects. -""" +"""Functions for combining many ComputedEntry objects into MigrationGraph objects.""" + from __future__ import annotations import logging @@ -79,23 +78,16 @@ def process_entries( # grouping of inserted structures with base structures all_sga = [ - SpacegroupAnalyzer( - itr_base_ent.structure, symprec=symprec, angle_tolerance=angle_tol - ) + SpacegroupAnalyzer(itr_base_ent.structure, symprec=symprec, angle_tolerance=angle_tol) for itr_base_ent in base_entries ] entries_with_num_symmetry_ops = [ - (ient, len(all_sga[itr_ent].get_space_group_operations())) - for itr_ent, ient in enumerate(base_entries) + (ient, len(all_sga[itr_ent].get_space_group_operations())) for itr_ent, ient in enumerate(base_entries) ] - entries_with_num_symmetry_ops = sorted( - entries_with_num_symmetry_ops, key=lambda x: x[0].energy_per_atom - ) - entries_with_num_symmetry_ops = sorted( - entries_with_num_symmetry_ops, key=lambda x: x[1], reverse=True - ) + entries_with_num_symmetry_ops = sorted(entries_with_num_symmetry_ops, key=lambda x: x[0].energy_per_atom) + entries_with_num_symmetry_ops = sorted(entries_with_num_symmetry_ops, key=lambda x: x[1], reverse=True) entries_with_num_symmetry_ops = sorted( entries_with_num_symmetry_ops, key=lambda x: x[0].structure.num_sites, @@ -131,9 +123,7 @@ def _meta_stable_sites(base_ent, inserted_ent): ) continue - struct_sym = get_sym_migration_ion_sites( - base_ent.structure, struct_wo_sym_ops, working_ion - ) + struct_sym = get_sym_migration_ion_sites(base_ent.structure, struct_wo_sym_ops, working_ion) results.append( { "base": base_ent.structure, @@ -142,14 +132,10 @@ def _meta_stable_sites(base_ent, inserted_ent): ) results = filter(lambda x: len(x["inserted"]) != 0, results) # type: ignore - return sorted( - results, key=lambda x: x["inserted"].composition[working_ion], reverse=True - ) + return sorted(results, key=lambda x: x["inserted"].composition[working_ion], reverse=True) -def get_matched_structure_mapping( - base: Structure, inserted: Structure, sm: StructureMatcher -): +def get_matched_structure_mapping(base: Structure, inserted: Structure, sm: StructureMatcher): """ Get the mapping from the inserted structure onto the base structure, assuming that the inserted structure sans the working ion is some kind @@ -167,15 +153,11 @@ def get_matched_structure_mapping( s1, s2 = sm._process_species([base, inserted]) fu, _ = sm._get_supercell_size(s1, s2) try: - val, dist, sc_m, total_t, mapping = sm._strict_match( - s1, s2, fu=fu, s1_supercell=True - ) + val, dist, sc_m, total_t, mapping = sm._strict_match(s1, s2, fu=fu, s1_supercell=True) except TypeError: return None sc = s1 * sc_m - sc.lattice = Lattice.from_parameters( - *sc.lattice.abc, *sc.lattice.angles, vesta=True - ) # type: ignore + sc.lattice = Lattice.from_parameters(*sc.lattice.abc, *sc.lattice.angles, vesta=True) # type: ignore return sc_m, total_t @@ -201,9 +183,7 @@ def get_inserted_on_base( Returns: List of entries for each working ion in the list of """ - mapped_result = get_matched_structure_mapping( - base_ent.structure, inserted_ent.structure, sm - ) + mapped_result = get_matched_structure_mapping(base_ent.structure, inserted_ent.structure, sm) if mapped_result is None: return None @@ -239,9 +219,9 @@ def get_sym_migration_ion_sites( from the base and inserted entries. Args: - inserted_entry: entry that contains cation - base_struct_entry: the entry containing the base structure - migrating_ion_entry: the name of the migrating species + base_struct: the base structure. + inserted_struct: inserted structure. + migrating_ion: Ion that migrates. symprec: the symprec tolerance for the space group analysis angle_tol: the angle tolerance for the space group analysis @@ -274,9 +254,7 @@ def get_sym_migration_ion_sites( # must clean up as you go or the number of sites explodes if len(sym_migration_struct) > 1: - sym_migration_struct.merge_sites( - tol=SITE_MERGE_R, mode="average" - ) # keeps removing duplicates + sym_migration_struct.merge_sites(tol=SITE_MERGE_R, mode="average") # keeps removing duplicates return sym_migration_struct @@ -290,9 +268,7 @@ def _filter_and_merge(inserted_structure: Structure) -> Structure | None: migration_sites = [] base_sites = [] for i_site in inserted_structure: - if "insertion_energy" in i_site.properties and isinstance( - i_site.properties["insertion_energy"], float - ): + if "insertion_energy" in i_site.properties and isinstance(i_site.properties["insertion_energy"], float): migration_sites.append(i_site) else: base_sites.append(i_site) @@ -327,7 +303,7 @@ def get_insertion_energy( (E[inserted] - (E[Base] + n * E[working_ion])) / n Where n is the number of working ions and E[inserted]. Additionally, and E[base] and E[inserted] are for structures of the same size - (sans working ion) + (sans working ion). """ wi_ = str(migrating_ion_entry.composition.elements[0]) comp_inserted_no_wi = inserted_entry.composition.as_dict() diff --git a/pymatgen/analysis/diffusion/utils/supercells.py b/pymatgen/analysis/diffusion/utils/supercells.py index ec59185..8b49e79 100644 --- a/pymatgen/analysis/diffusion/utils/supercells.py +++ b/pymatgen/analysis/diffusion/utils/supercells.py @@ -1,8 +1,5 @@ -# Copyright (c) Materials Virtual Lab. -# Distributed under the terms of the BSD License. -""" -Functions for creating supercells for NEB calculations -""" +"""Functions for creating supercells for NEB calculations.""" + from __future__ import annotations import logging @@ -22,7 +19,7 @@ logger = logging.getLogger(__name__) -# Helper functions for MigraionHop.get_sc_struture +# Helper functions for MigrationHop.get_sc_struture def get_sc_fromstruct( @@ -62,7 +59,7 @@ def _get_sc_from_struct_pmg( min_length: float = 10.0, ) -> list[list[int]] | None: """ - Generate the best supercell from a unitcell using the pymatgen CubicSupercellTransformation + Generate the best supercell from a unitcell using the pymatgen CubicSupercellTransformation. Args: base_struct: structure of the unit cell @@ -74,9 +71,7 @@ def _get_sc_from_struct_pmg( 3x3 matrix: supercell matrix """ - cst = CubicSupercellTransformation( - min_atoms=min_atoms, max_atoms=max_atoms, min_length=min_length - ) + cst = CubicSupercellTransformation(min_atoms=min_atoms, max_atoms=max_atoms, min_length=min_length) try: cst.apply_transformation(base_struct) @@ -85,54 +80,6 @@ def _get_sc_from_struct_pmg( return cst.transformation_matrix -# Something is broken with how ASE generates supercells -# now many calls to `find_optimal_cell_shape` results in -# `Failed to find a transformation matrix` -# So remove this functionality now in favour of decreasing the max_length -# -# def _get_sc_from_struct_ase( -# base_struct: Structure, min_atoms: int=80, max_atoms: int=240 -# ) -> List[List[int]]: -# """generate the best supercell from a unit cell using ASE's method, is slower but -# more exhaustive -# -# Args: -# base_struct (pymatgen.Structure): unit cell -# min_atoms (int, optional): Minimum number of atoms in the desired supercell. Defaults to 80. -# -# Returns: -# 3x3 matrix: supercell matrix -# """ -# -# num_cells_min = int(np.ceil(min_atoms / base_struct.num_sites)) -# num_cells_max = int(np.ceil(max_atoms / base_struct.num_sites)) -# -# atoms = AseAtomsAdaptor().get_atoms(base_struct) -# res = [] -# for icell in range(num_cells_min, num_cells_max): -# if icell % 2 != 0 and icell % 3 != 0: -# continue # cells with many factors are more likely to be square -# logger.info(f"Getting cell shape {icell} x unit cell") -# sc_mat = find_optimal_cell_shape(atoms.cell, icell, "sc") -# if sc_mat is None: -# continue -# deviation = get_deviation_from_optimal_cell_shape(np.dot(sc_mat, atoms.cell)) -# dd = {"deviation": deviation, "P1": sc_mat} -# # base_struct_sc.to('poscar', f'endpoints_{base_struct_sc.num_sites}.vasp') -# res.append(dd) -# if deviation < 0.3: -# logger.debug("Found good cell") -# return sc_mat -# -# else: -# best_case = min(res, key=lambda x: x["deviation"]) -# logger.warning( -# f"Could not find case with deviation from cubic was less than 0.3 using the ASE cubic supercell finder \ -# \nThe best one had a deviation of {best_case['deviation']}" -# ) -# return best_case["P1"] - - def get_start_end_structures( isite: PeriodicSite, esite: PeriodicSite, @@ -146,10 +93,14 @@ def get_start_end_structures( Obtain the starting and terminating structures in a supercell for NEB calculations. Args: + isite: Initial site index. + esite: End site index. hop: object presenting the migration event base_struct: unit cell representation of the structure sc_mat: supercell transformation to create the simulation cell for the NEB calc - tol: toleranace for identifying isite/esite within base_struct + vac_mode: Vacuum mode. + debug: debug mode. + tol: toleranace for identifying isite/esite within base_struct. Returns: initial structure, final structure, empty structure all in the supercell @@ -191,19 +142,9 @@ def remove_site_at_pos(structure: Structure, site: PeriodicSite, tol: float): if debug: icart = base_sc.lattice.get_cartesian_coords(ipos_sc) ecart = base_sc.lattice.get_cartesian_coords(epos_sc) - assert ( - abs( - np.linalg.norm(icart - ecart) - - np.linalg.norm(isite.coords - esite.coords) - ) - < 1e-5 - ) - i_ref_ = PeriodicSite( - species=esite.species_string, coords=ipos_sc, lattice=base_sc.lattice - ) - e_ref_ = PeriodicSite( - species=esite.species_string, coords=epos_sc, lattice=base_sc.lattice - ) + assert abs(np.linalg.norm(icart - ecart) - np.linalg.norm(isite.coords - esite.coords)) < 1e-5 + i_ref_ = PeriodicSite(species=esite.species_string, coords=ipos_sc, lattice=base_sc.lattice) + e_ref_ = PeriodicSite(species=esite.species_string, coords=epos_sc, lattice=base_sc.lattice) start_struct = remove_site_at_pos(start_struct, e_ref_, tol) end_struct = remove_site_at_pos(end_struct, i_ref_, tol) return start_struct, end_struct, base_sc