diff --git a/pymatgen/analysis/chemenv/coordination_environments/coordination_geometry_finder.py b/pymatgen/analysis/chemenv/coordination_environments/coordination_geometry_finder.py index b8398ac5e4b..165978b3f97 100644 --- a/pymatgen/analysis/chemenv/coordination_environments/coordination_geometry_finder.py +++ b/pymatgen/analysis/chemenv/coordination_environments/coordination_geometry_finder.py @@ -120,9 +120,9 @@ def __init__( if central_site is None: raise ValueError("The centroid includes the central site but no central site is given") total += self.bare_centre - self.centre = total / (np.float_(len(bare_coords)) + 1.0) + self.centre = total / (np.float64(len(bare_coords)) + 1.0) else: - self.centre = total / np.float_(len(bare_coords)) + self.centre = total / np.float64(len(bare_coords)) elif centering_type == "central_site": if include_central_site_in_centroid: raise ValueError( @@ -138,9 +138,9 @@ def __init__( if central_site is None: raise ValueError("The centroid includes the central site but no central site is given") total += self.bare_centre - self.centre = total / (np.float_(len(bare_coords)) + 1.0) + self.centre = total / (np.float64(len(bare_coords)) + 1.0) else: - self.centre = total / np.float_(len(bare_coords)) + self.centre = total / np.float64(len(bare_coords)) self._bare_coords = self.bare_points_without_centre self._coords = self._bare_coords - self.centre self.central_site = self.bare_central_site - self.centre diff --git a/pymatgen/analysis/chemenv/utils/coordination_geometry_utils.py b/pymatgen/analysis/chemenv/utils/coordination_geometry_utils.py index 37f0784da32..e47f3bd6636 100644 --- a/pymatgen/analysis/chemenv/utils/coordination_geometry_utils.py +++ b/pymatgen/analysis/chemenv/utils/coordination_geometry_utils.py @@ -631,12 +631,12 @@ def __init__(self, coefficients, p1=None, p2=None, p3=None): raise ValueError("Normal vector is equal to 0.0") if self.normal_vector[non_zeros[0]] < 0.0: self.normal_vector = -self.normal_vector - dd = -np.float_(coefficients[3]) / norm_v + dd = -np.float64(coefficients[3]) / norm_v else: - dd = np.float_(coefficients[3]) / norm_v + dd = np.float64(coefficients[3]) / norm_v self._coefficients = np.array( [self.normal_vector[0], self.normal_vector[1], self.normal_vector[2], dd], - np.float_, + np.float64, ) self._crosses_origin = np.isclose(dd, 0.0, atol=1e-7, rtol=0.0) self.p1 = p1 diff --git a/pymatgen/analysis/elasticity/elastic.py b/pymatgen/analysis/elasticity/elastic.py index 0df68712e08..c512e6521cc 100644 --- a/pymatgen/analysis/elasticity/elastic.py +++ b/pymatgen/analysis/elasticity/elastic.py @@ -1080,7 +1080,7 @@ def get_diff_coeff(hvec, n=1): hvec (1D array-like): sampling stencil n (int): degree of derivative to find """ - hvec = np.array(hvec, dtype=np.float_) + hvec = np.array(hvec, dtype=np.float64) acc = len(hvec) exp = np.column_stack([np.arange(acc)] * acc) a = np.vstack([hvec] * acc) ** exp diff --git a/pymatgen/analysis/ewald.py b/pymatgen/analysis/ewald.py index ed7712176dd..e2e37127dba 100644 --- a/pymatgen/analysis/ewald.py +++ b/pymatgen/analysis/ewald.py @@ -317,8 +317,8 @@ def _calc_recip(self): """ n_sites = len(self._struct) prefactor = 2 * pi / self._vol - e_recip = np.zeros((n_sites, n_sites), dtype=np.float_) - forces = np.zeros((n_sites, 3), dtype=np.float_) + e_recip = np.zeros((n_sites, n_sites), dtype=np.float64) + forces = np.zeros((n_sites, 3), dtype=np.float64) coords = self._coords rcp_latt = self._struct.lattice.reciprocal_lattice recip_nn = rcp_latt.get_points_in_sphere([[0, 0, 0]], [0, 0, 0], self._gmax) @@ -363,9 +363,9 @@ def _calc_real_and_point(self): force_pf = 2 * self._sqrt_eta / sqrt(pi) coords = self._coords n_sites = len(self._struct) - e_real = np.empty((n_sites, n_sites), dtype=np.float_) + e_real = np.empty((n_sites, n_sites), dtype=np.float64) - forces = np.zeros((n_sites, 3), dtype=np.float_) + forces = np.zeros((n_sites, 3), dtype=np.float64) qs = np.array(self._oxi_states) diff --git a/pymatgen/analysis/graphs.py b/pymatgen/analysis/graphs.py index 00630572924..fcbf0bcbdf8 100644 --- a/pymatgen/analysis/graphs.py +++ b/pymatgen/analysis/graphs.py @@ -654,8 +654,8 @@ def map_indices(grp): atoms = len(grp) - 1 offset = len(self.structure) - atoms - for i in range(atoms): - grp_map[i] = i + offset + for idx in range(atoms): + grp_map[idx] = idx + offset return grp_map @@ -1011,7 +1011,7 @@ def types_of_coordination_environments(self, anonymous=False): if anonymous: mapping = {centre_sp: "A"} - available_letters = [chr(66 + i) for i in range(25)] + available_letters = [chr(66 + idx) for idx in range(25)] for label in labels: sp = label[1] if sp not in mapping: @@ -2164,8 +2164,8 @@ def map_indices(grp): atoms = len(grp) - 1 offset = len(self.molecule) - atoms - for i in range(atoms): - grp_map[i] = i + offset + for idx in range(atoms): + grp_map[idx] = idx + offset return grp_map @@ -2298,9 +2298,9 @@ def replace_group( sizes[neighbor[2]] = len(nx.descendants(disconnected, neighbor[2])) keep = max(sizes, key=lambda x: sizes[x]) - for i in sizes: - if i != keep: - to_remove.add(i) + for idx in sizes: + if idx != keep: + to_remove.add(idx) self.remove_nodes(list(to_remove)) self.substitute_group( diff --git a/pymatgen/analysis/interface_reactions.py b/pymatgen/analysis/interface_reactions.py index cbd2c980d70..c8a19222d5c 100644 --- a/pymatgen/analysis/interface_reactions.py +++ b/pymatgen/analysis/interface_reactions.py @@ -155,30 +155,30 @@ def get_kinks(self) -> list[tuple[int, float, float, Reaction, float]]: react_kink = [self._get_reaction(x) for x in x_kink] num_atoms = [(x * self.comp1.num_atoms + (1 - x) * self.comp2.num_atoms) for x in x_kink] energy_per_rxt_formula = [ - energy_kink[i] - * self._get_elem_amt_in_rxn(react_kink[i]) - / num_atoms[i] + energy_kink[idx] + * self._get_elem_amt_in_rxn(react_kink[idx]) + / num_atoms[idx] * InterfacialReactivity.EV_TO_KJ_PER_MOL - for i in range(2) + for idx in range(2) ] else: - for i in reversed(critical_comp): + for idx in reversed(critical_comp): # Gets mixing ratio x at kinks. - c = self.pd.pd_coords(i) - x = float(np.linalg.norm(c - c2_coord) / np.linalg.norm(c1_coord - c2_coord)) + coords = self.pd.pd_coords(idx) + mixing_ratio = float(np.linalg.norm(coords - c2_coord) / np.linalg.norm(c1_coord - c2_coord)) # Modifies mixing ratio in case compositions self.comp1 and # self.comp2 are not normalized. - x = x * n2 / (n1 + x * (n2 - n1)) - n_atoms = x * self.comp1.num_atoms + (1 - x) * self.comp2.num_atoms + mixing_ratio = mixing_ratio * n2 / (n1 + mixing_ratio * (n2 - n1)) + n_atoms = mixing_ratio * self.comp1.num_atoms + (1 - mixing_ratio) * self.comp2.num_atoms # Converts mixing ratio in comp1 - comp2 tie line to that in # c1 - c2 tie line. - x_converted = self._convert(x, self.factor1, self.factor2) + x_converted = self._convert(mixing_ratio, self.factor1, self.factor2) x_kink.append(x_converted) # Gets reaction energy at kinks - normalized_energy = self._get_energy(x) + normalized_energy = self._get_energy(mixing_ratio) energy_kink.append(normalized_energy) # Gets balanced reaction at kinks - rxt = self._get_reaction(x) + rxt = self._get_reaction(mixing_ratio) react_kink.append(rxt) rxt_energy = normalized_energy * self._get_elem_amt_in_rxn(rxt) / n_atoms energy_per_rxt_formula.append(rxt_energy * self.EV_TO_KJ_PER_MOL) diff --git a/pymatgen/analysis/local_env.py b/pymatgen/analysis/local_env.py index a07884b17db..677f46fdde4 100644 --- a/pymatgen/analysis/local_env.py +++ b/pymatgen/analysis/local_env.py @@ -751,8 +751,6 @@ def get_voronoi_polyhedra(self, structure: Structure, n: int): targets = structure.elements if self.targets is None else self.targets center = structure[n] - cutoff = self.cutoff - # max cutoff is the longest diagonal of the cell + room for noise corners = [[1, 1, 1], [-1, 1, 1], [1, -1, 1], [1, 1, -1]] d_corners = [np.linalg.norm(structure.lattice.get_cartesian_coords(c)) for c in corners] @@ -760,11 +758,11 @@ def get_voronoi_polyhedra(self, structure: Structure, n: int): while True: try: - neighbors = structure.get_sites_in_sphere(center.coords, cutoff) - neighbors = [i[0] for i in sorted(neighbors, key=lambda s: s[1])] + neighbors = structure.get_sites_in_sphere(center.coords, self.cutoff) + neighbors = [ngbr[0] for ngbr in sorted(neighbors, key=lambda s: s[1])] # Run the Voronoi tessellation - qvoronoi_input = [s.coords for s in neighbors] + qvoronoi_input = [site.coords for site in neighbors] voro = Voronoi(qvoronoi_input) # can give seg fault if cutoff is too small @@ -773,12 +771,12 @@ def get_voronoi_polyhedra(self, structure: Structure, n: int): break except RuntimeError as exc: - if cutoff >= max_cutoff: + if self.cutoff >= max_cutoff: if exc.args and "vertex" in exc.args[0]: # pass through the error raised by _extract_cell_info raise exc raise RuntimeError("Error in Voronoi neighbor finding; max cutoff exceeded") - cutoff = min(cutoff * 2, max_cutoff + 0.001) + self.cutoff = min(self.cutoff * 2, max_cutoff + 0.001) return cell_info def get_all_voronoi_polyhedra(self, structure: Structure): @@ -826,7 +824,7 @@ def get_all_voronoi_polyhedra(self, structure: Structure): # Get the non-duplicates (using the site indices for numerical stability) indices = np.array(indices, dtype=int) # type: ignore indices, uniq_inds = np.unique(indices, return_index=True, axis=0) # type: ignore[assignment] - sites = [sites[i] for i in uniq_inds] + sites = [sites[idx] for idx in uniq_inds] # Sort array such that atoms in the root image are first # Exploit the fact that the array is sorted by the unique operation such that @@ -841,7 +839,8 @@ def get_all_voronoi_polyhedra(self, structure: Structure): # Get the information for each neighbor return [ - self._extract_cell_info(i, sites, targets, voro, self.compute_adj_neighbors) for i in root_images.tolist() + self._extract_cell_info(idx, sites, targets, voro, self.compute_adj_neighbors) + for idx in root_images.tolist() ] def _extract_cell_info(self, site_idx, sites, targets, voro, compute_adj_neighbors=False): @@ -889,7 +888,7 @@ def _extract_cell_info(self, site_idx, sites, targets, voro, compute_adj_neighbo raise RuntimeError("This structure is pathological, infinite vertex in the Voronoi construction") # Get the solid angle of the face - facets = [all_vertices[i] for i in vind] + facets = [all_vertices[idx] for idx in vind] angle = solid_angle(center_coords, facets) # Compute the volume of associated with this face @@ -950,7 +949,7 @@ def _extract_cell_info(self, site_idx, sites, targets, voro, compute_adj_neighbo # If desired, determine which neighbors are adjacent if compute_adj_neighbors: # Initialize storage for the adjacent neighbors - adj_neighbors = {i: [] for i in result_weighted} + adj_neighbors = {idx: [] for idx in result_weighted} # Find the neighbors that are adjacent by finding those # that contain exactly two vertices @@ -1926,22 +1925,22 @@ def solid_angle(center, coords): The solid angle. """ # Compute the displacement from the center - r = [np.subtract(c, center) for c in coords] + disp = [np.subtract(c, center) for c in coords] # Compute the magnitude of each vector - r_norm = [np.linalg.norm(i) for i in r] + r_norm = [np.linalg.norm(idx) for idx in disp] # Compute the solid angle for each tetrahedron that makes up the facet # Following: https://wikipedia.org/wiki/Solid_angle#Tetrahedron angle = 0 - for i in range(1, len(r) - 1): - j = i + 1 - tp = np.abs(np.dot(r[0], np.cross(r[i], r[j]))) + for ii in range(1, len(disp) - 1): + jj = ii + 1 + tp = np.abs(np.dot(disp[0], np.cross(disp[ii], disp[jj]))) de = ( - r_norm[0] * r_norm[i] * r_norm[j] - + r_norm[j] * np.dot(r[0], r[i]) - + r_norm[i] * np.dot(r[0], r[j]) - + r_norm[0] * np.dot(r[i], r[j]) + r_norm[0] * r_norm[ii] * r_norm[jj] + + r_norm[jj] * np.dot(disp[0], disp[ii]) + + r_norm[ii] * np.dot(disp[0], disp[jj]) + + r_norm[0] * np.dot(disp[ii], disp[jj]) ) _angle = (0.5 * pi if tp > 0 else -0.5 * pi) if de == 0 else np.arctan(tp / de) angle += (_angle if _angle > 0 else _angle + np.pi) * 2 @@ -2877,7 +2876,7 @@ def get_order_parameters( neighsites = vnn.get_nn(structure, n) else: # Structure.get_sites_in_sphere --> also other periodic images - neighsitestmp = [i[0] for i in structure.get_sites_in_sphere(centsite.coords, self._cutoff)] + neighsitestmp = [idx[0] for idx in structure.get_sites_in_sphere(centsite.coords, self._cutoff)] neighsites = [] if centsite not in neighsitestmp: raise ValueError("Could not find center site!") @@ -2923,15 +2922,15 @@ def get_order_parameters( # norms = [[[] for j in range(nneigh)] for t in self._types] # First, coordination number and distance-based OPs. - for i, t in enumerate(self._types): - if t == "cn": - ops[i] = nneigh / self._params[i]["norm"] - elif t == "sgl_bd": + for idx, typ in enumerate(self._types): + if typ == "cn": + ops[idx] = nneigh / self._params[idx]["norm"] + elif typ == "sgl_bd": dist_sorted = sorted(dist) if len(dist_sorted) == 1: - ops[i] = 1 + ops[idx] = 1 elif len(dist_sorted) > 1: - ops[i] = 1 - dist_sorted[0] / dist_sorted[1] + ops[idx] = 1 - dist_sorted[0] / dist_sorted[1] # Then, bond orientational OPs based on spherical harmonics # according to Steinhardt et al., Phys. Rev. B, 28, 784-805, 1983. @@ -2961,13 +2960,13 @@ def get_order_parameters( # Note that None flags that we have too few neighbors # for calculating BOOPS. - for i, t in enumerate(self._types): - if t == "q2": - ops[i] = self.get_q2(thetas, phis) if len(thetas) > 0 else None - elif t == "q4": - ops[i] = self.get_q4(thetas, phis) if len(thetas) > 0 else None - elif t == "q6": - ops[i] = self.get_q6(thetas, phis) if len(thetas) > 0 else None + for idx, typ in enumerate(self._types): + if typ == "q2": + ops[idx] = self.get_q2(thetas, phis) if len(thetas) > 0 else None + elif typ == "q4": + ops[idx] = self.get_q4(thetas, phis) if len(thetas) > 0 else None + elif typ == "q6": + ops[idx] = self.get_q6(thetas, phis) if len(thetas) > 0 else None # Then, deal with the Peters-style OPs that are tailor-made # to recognize common structural motifs @@ -2986,9 +2985,9 @@ def get_order_parameters( kc = 0 for k in range(nneigh): # From neighbor k, we construct if j != k: # the prime meridian. - for i in range(len(self._types)): - qsp_theta[i][j].append(0.0) - norms[i][j].append(0) + for idx in range(len(self._types)): + qsp_theta[idx][j].append(0.0) + norms[idx][j].append(0) tmp = max(-1.0, min(np.inner(zaxis, rij_norm[k]), 1.0)) thetak = acos(tmp) xaxis = gramschmidt(rij_norm[k], zaxis) @@ -3006,33 +3005,33 @@ def get_order_parameters( # Contributions of j-i-k angles, where i represents the # central atom and j and k two of the neighbors. - for i, t in enumerate(self._types): - if t in ["bent", "sq_pyr_legacy"]: - tmp = self._params[i]["IGW_TA"] * (thetak * ipi - self._params[i]["TA"]) - qsp_theta[i][j][kc] += exp(-0.5 * tmp * tmp) - norms[i][j][kc] += 1 - elif t in ["tri_plan", "tri_plan_max", "tet", "tet_max"]: - tmp = self._params[i]["IGW_TA"] * (thetak * ipi - self._params[i]["TA"]) - gaussthetak[i] = exp(-0.5 * tmp * tmp) - if t in ["tri_plan_max", "tet_max"]: - qsp_theta[i][j][kc] += gaussthetak[i] - norms[i][j][kc] += 1 - elif t in ["T", "tri_pyr", "sq_pyr", "pent_pyr", "hex_pyr"]: - tmp = self._params[i]["IGW_EP"] * (thetak * ipi - 0.5) - qsp_theta[i][j][kc] += exp(-0.5 * tmp * tmp) - norms[i][j][kc] += 1 - elif t in [ + for idx, typ in enumerate(self._types): + if typ in ["bent", "sq_pyr_legacy"]: + tmp = self._params[idx]["IGW_TA"] * (thetak * ipi - self._params[idx]["TA"]) + qsp_theta[idx][j][kc] += exp(-0.5 * tmp * tmp) + norms[idx][j][kc] += 1 + elif typ in ["tri_plan", "tri_plan_max", "tet", "tet_max"]: + tmp = self._params[idx]["IGW_TA"] * (thetak * ipi - self._params[idx]["TA"]) + gaussthetak[idx] = exp(-0.5 * tmp * tmp) + if typ in ["tri_plan_max", "tet_max"]: + qsp_theta[idx][j][kc] += gaussthetak[idx] + norms[idx][j][kc] += 1 + elif typ in ["T", "tri_pyr", "sq_pyr", "pent_pyr", "hex_pyr"]: + tmp = self._params[idx]["IGW_EP"] * (thetak * ipi - 0.5) + qsp_theta[idx][j][kc] += exp(-0.5 * tmp * tmp) + norms[idx][j][kc] += 1 + elif typ in [ "sq_plan", "oct", "oct_legacy", "cuboct", "cuboct_max", ]: - if thetak >= self._params[i]["min_SPP"]: - tmp = self._params[i]["IGW_SPP"] * (thetak * ipi - 1.0) - qsp_theta[i][j][kc] += self._params[i]["w_SPP"] * exp(-0.5 * tmp * tmp) - norms[i][j][kc] += self._params[i]["w_SPP"] - elif t in [ + if thetak >= self._params[idx]["min_SPP"]: + tmp = self._params[idx]["IGW_SPP"] * (thetak * ipi - 1.0) + qsp_theta[idx][j][kc] += self._params[idx]["w_SPP"] * exp(-0.5 * tmp * tmp) + norms[idx][j][kc] += self._params[idx]["w_SPP"] + elif typ in [ "see_saw_rect", "tri_bipyr", "sq_bipyr", @@ -3042,31 +3041,31 @@ def get_order_parameters( "sq_plan_max", "hex_plan_max", ]: - if thetak < self._params[i]["min_SPP"]: + if thetak < self._params[idx]["min_SPP"]: tmp = ( - self._params[i]["IGW_EP"] * (thetak * ipi - 0.5) - if t != "hex_plan_max" - else self._params[i]["IGW_TA"] - * (fabs(thetak * ipi - 0.5) - self._params[i]["TA"]) + self._params[idx]["IGW_EP"] * (thetak * ipi - 0.5) + if typ != "hex_plan_max" + else self._params[idx]["IGW_TA"] + * (fabs(thetak * ipi - 0.5) - self._params[idx]["TA"]) ) - qsp_theta[i][j][kc] += exp(-0.5 * tmp * tmp) - norms[i][j][kc] += 1 - elif t in ["pent_plan", "pent_plan_max"]: - tmp = 0.4 if thetak <= self._params[i]["TA"] * pi else 0.8 - tmp2 = self._params[i]["IGW_TA"] * (thetak * ipi - tmp) - gaussthetak[i] = exp(-0.5 * tmp2 * tmp2) - if t == "pent_plan_max": - qsp_theta[i][j][kc] += gaussthetak[i] - norms[i][j][kc] += 1 - elif t == "bcc" and j < k: - if thetak >= self._params[i]["min_SPP"]: - tmp = self._params[i]["IGW_SPP"] * (thetak * ipi - 1.0) - qsp_theta[i][j][kc] += self._params[i]["w_SPP"] * exp(-0.5 * tmp * tmp) - norms[i][j][kc] += self._params[i]["w_SPP"] - elif t == "sq_face_cap_trig_pris" and thetak < self._params[i]["TA3"]: - tmp = self._params[i]["IGW_TA1"] * (thetak * ipi - self._params[i]["TA1"]) - qsp_theta[i][j][kc] += exp(-0.5 * tmp * tmp) - norms[i][j][kc] += 1 + qsp_theta[idx][j][kc] += exp(-0.5 * tmp * tmp) + norms[idx][j][kc] += 1 + elif typ in ["pent_plan", "pent_plan_max"]: + tmp = 0.4 if thetak <= self._params[idx]["TA"] * pi else 0.8 + tmp2 = self._params[idx]["IGW_TA"] * (thetak * ipi - tmp) + gaussthetak[idx] = exp(-0.5 * tmp2 * tmp2) + if typ == "pent_plan_max": + qsp_theta[idx][j][kc] += gaussthetak[idx] + norms[idx][j][kc] += 1 + elif typ == "bcc" and j < k: + if thetak >= self._params[idx]["min_SPP"]: + tmp = self._params[idx]["IGW_SPP"] * (thetak * ipi - 1.0) + qsp_theta[idx][j][kc] += self._params[idx]["w_SPP"] * exp(-0.5 * tmp * tmp) + norms[idx][j][kc] += self._params[idx]["w_SPP"] + elif typ == "sq_face_cap_trig_pris" and thetak < self._params[idx]["TA3"]: + tmp = self._params[idx]["IGW_TA1"] * (thetak * ipi - self._params[idx]["TA1"]) + qsp_theta[idx][j][kc] += exp(-0.5 * tmp * tmp) + norms[idx][j][kc] += 1 for m in range(nneigh): if (m != j) and (m != k) and (not flag_xaxis): @@ -3087,7 +3086,7 @@ def get_order_parameters( ) # South pole contributions of m. if ( - t + typ in [ "tri_bipyr", "sq_bipyr", @@ -3098,60 +3097,68 @@ def get_order_parameters( "hex_plan_max", "see_saw_rect", ] - and thetam >= self._params[i]["min_SPP"] + and thetam >= self._params[idx]["min_SPP"] ): - tmp = self._params[i]["IGW_SPP"] * (thetam * ipi - 1.0) - qsp_theta[i][j][kc] += exp(-0.5 * tmp * tmp) - norms[i][j][kc] += 1 + tmp = self._params[idx]["IGW_SPP"] * (thetam * ipi - 1.0) + qsp_theta[idx][j][kc] += exp(-0.5 * tmp * tmp) + norms[idx][j][kc] += 1 # Contributions of j-i-m angle and # angles between plane j-i-k and i-m vector. if not flag_xaxis and not flag_xtwoaxis: - for i, t in enumerate(self._types): - if t in [ + for idx, typ in enumerate(self._types): + if typ in [ "tri_plan", "tri_plan_max", "tet", "tet_max", ]: - tmp = self._params[i]["IGW_TA"] * (thetam * ipi - self._params[i]["TA"]) - tmp2 = cos(self._params[i]["fac_AA"] * phi) ** self._params[i]["exp_cos_AA"] - tmp3 = 1 if t in ["tri_plan_max", "tet_max"] else gaussthetak[i] - qsp_theta[i][j][kc] += tmp3 * exp(-0.5 * tmp * tmp) * tmp2 - norms[i][j][kc] += 1 - elif t in ["pent_plan", "pent_plan_max"]: - tmp = 0.4 if thetam <= self._params[i]["TA"] * pi else 0.8 - tmp2 = self._params[i]["IGW_TA"] * (thetam * ipi - tmp) + tmp = self._params[idx]["IGW_TA"] * (thetam * ipi - self._params[idx]["TA"]) + tmp2 = ( + cos(self._params[idx]["fac_AA"] * phi) + ** self._params[idx]["exp_cos_AA"] + ) + tmp3 = 1 if typ in ["tri_plan_max", "tet_max"] else gaussthetak[idx] + qsp_theta[idx][j][kc] += tmp3 * exp(-0.5 * tmp * tmp) * tmp2 + norms[idx][j][kc] += 1 + elif typ in ["pent_plan", "pent_plan_max"]: + tmp = 0.4 if thetam <= self._params[idx]["TA"] * pi else 0.8 + tmp2 = self._params[idx]["IGW_TA"] * (thetam * ipi - tmp) tmp3 = cos(phi) - tmp4 = 1 if t == "pent_plan_max" else gaussthetak[i] - qsp_theta[i][j][kc] += tmp4 * exp(-0.5 * tmp2 * tmp2) * tmp3 * tmp3 - norms[i][j][kc] += 1 - elif t in [ + tmp4 = 1 if typ == "pent_plan_max" else gaussthetak[idx] + qsp_theta[idx][j][kc] += tmp4 * exp(-0.5 * tmp2 * tmp2) * tmp3 * tmp3 + norms[idx][j][kc] += 1 + elif typ in [ "T", "tri_pyr", "sq_pyr", "pent_pyr", "hex_pyr", ]: - tmp = cos(self._params[i]["fac_AA"] * phi) ** self._params[i]["exp_cos_AA"] - tmp3 = self._params[i]["IGW_EP"] * (thetam * ipi - 0.5) - qsp_theta[i][j][kc] += tmp * exp(-0.5 * tmp3 * tmp3) - norms[i][j][kc] += 1 - elif t in ["sq_plan", "oct", "oct_legacy"]: + tmp = ( + cos(self._params[idx]["fac_AA"] * phi) + ** self._params[idx]["exp_cos_AA"] + ) + tmp3 = self._params[idx]["IGW_EP"] * (thetam * ipi - 0.5) + qsp_theta[idx][j][kc] += tmp * exp(-0.5 * tmp3 * tmp3) + norms[idx][j][kc] += 1 + elif typ in ["sq_plan", "oct", "oct_legacy"]: if ( - thetak < self._params[i]["min_SPP"] - and thetam < self._params[i]["min_SPP"] + thetak < self._params[idx]["min_SPP"] + and thetam < self._params[idx]["min_SPP"] ): tmp = ( - cos(self._params[i]["fac_AA"] * phi) - ** self._params[i]["exp_cos_AA"] + cos(self._params[idx]["fac_AA"] * phi) + ** self._params[idx]["exp_cos_AA"] ) - tmp2 = self._params[i]["IGW_EP"] * (thetam * ipi - 0.5) - qsp_theta[i][j][kc] += tmp * exp(-0.5 * tmp2 * tmp2) - if t == "oct_legacy": - qsp_theta[i][j][kc] -= tmp * self._params[i][6] * self._params[i][7] - norms[i][j][kc] += 1 - elif t in [ + tmp2 = self._params[idx]["IGW_EP"] * (thetam * ipi - 0.5) + qsp_theta[idx][j][kc] += tmp * exp(-0.5 * tmp2 * tmp2) + if typ == "oct_legacy": + qsp_theta[idx][j][kc] -= ( + tmp * self._params[idx][6] * self._params[idx][7] + ) + norms[idx][j][kc] += 1 + elif typ in [ "tri_bipyr", "sq_bipyr", "pent_bipyr", @@ -3161,99 +3168,99 @@ def get_order_parameters( "hex_plan_max", ]: if ( - thetam < self._params[i]["min_SPP"] - and thetak < self._params[i]["min_SPP"] + thetam < self._params[idx]["min_SPP"] + and thetak < self._params[idx]["min_SPP"] ): tmp = ( - cos(self._params[i]["fac_AA"] * phi) - ** self._params[i]["exp_cos_AA"] + cos(self._params[idx]["fac_AA"] * phi) + ** self._params[idx]["exp_cos_AA"] ) tmp2 = ( - self._params[i]["IGW_EP"] * (thetam * ipi - 0.5) - if t != "hex_plan_max" - else self._params[i]["IGW_TA"] - * (fabs(thetam * ipi - 0.5) - self._params[i]["TA"]) + self._params[idx]["IGW_EP"] * (thetam * ipi - 0.5) + if typ != "hex_plan_max" + else self._params[idx]["IGW_TA"] + * (fabs(thetam * ipi - 0.5) - self._params[idx]["TA"]) ) - qsp_theta[i][j][kc] += tmp * exp(-0.5 * tmp2 * tmp2) - norms[i][j][kc] += 1 - elif t == "bcc" and j < k: - if thetak < self._params[i]["min_SPP"]: + qsp_theta[idx][j][kc] += tmp * exp(-0.5 * tmp2 * tmp2) + norms[idx][j][kc] += 1 + elif typ == "bcc" and j < k: + if thetak < self._params[idx]["min_SPP"]: fac = 1 if thetak > piover2 else -1 tmp = (thetam - piover2) / asin(1 / 3) - qsp_theta[i][j][kc] += ( + qsp_theta[idx][j][kc] += ( fac * cos(3 * phi) * fac_bcc * tmp * exp(-0.5 * tmp * tmp) ) - norms[i][j][kc] += 1 - elif t == "see_saw_rect": + norms[idx][j][kc] += 1 + elif typ == "see_saw_rect": if ( - thetam < self._params[i]["min_SPP"] - and thetak < self._params[i]["min_SPP"] + thetam < self._params[idx]["min_SPP"] + and thetak < self._params[idx]["min_SPP"] and phi < 0.75 * pi ): tmp = ( - cos(self._params[i]["fac_AA"] * phi) - ** self._params[i]["exp_cos_AA"] + cos(self._params[idx]["fac_AA"] * phi) + ** self._params[idx]["exp_cos_AA"] ) - tmp2 = self._params[i]["IGW_EP"] * (thetam * ipi - 0.5) - qsp_theta[i][j][kc] += tmp * exp(-0.5 * tmp2 * tmp2) - norms[i][j][kc] += 1.0 - elif t in ["cuboct", "cuboct_max"]: + tmp2 = self._params[idx]["IGW_EP"] * (thetam * ipi - 0.5) + qsp_theta[idx][j][kc] += tmp * exp(-0.5 * tmp2 * tmp2) + norms[idx][j][kc] += 1.0 + elif typ in ["cuboct", "cuboct_max"]: if ( - thetam < self._params[i]["min_SPP"] - and self._params[i][4] < thetak < self._params[i][2] + thetam < self._params[idx]["min_SPP"] + and self._params[idx][4] < thetak < self._params[idx][2] ): - if self._params[i][4] < thetam < self._params[i][2]: + if self._params[idx][4] < thetam < self._params[idx][2]: tmp = cos(phi) - tmp2 = self._params[i][5] * (thetam * ipi - 0.5) - qsp_theta[i][j][kc] += tmp * tmp * exp(-0.5 * tmp2 * tmp2) - norms[i][j][kc] += 1.0 - elif thetam < self._params[i][4]: + tmp2 = self._params[idx][5] * (thetam * ipi - 0.5) + qsp_theta[idx][j][kc] += tmp * tmp * exp(-0.5 * tmp2 * tmp2) + norms[idx][j][kc] += 1.0 + elif thetam < self._params[idx][4]: tmp = 0.0556 * (cos(phi - 0.5 * pi) - 0.81649658) - tmp2 = self._params[i][6] * (thetam * ipi - onethird) - qsp_theta[i][j][kc] += exp(-0.5 * tmp * tmp) * exp( + tmp2 = self._params[idx][6] * (thetam * ipi - onethird) + qsp_theta[idx][j][kc] += exp(-0.5 * tmp * tmp) * exp( -0.5 * tmp2 * tmp2 ) - norms[i][j][kc] += 1.0 - elif thetam > self._params[i][2]: + norms[idx][j][kc] += 1.0 + elif thetam > self._params[idx][2]: tmp = 0.0556 * (cos(phi - 0.5 * pi) - 0.81649658) - tmp2 = self._params[i][6] * (thetam * ipi - twothird) - qsp_theta[i][j][kc] += exp(-0.5 * tmp * tmp) * exp( + tmp2 = self._params[idx][6] * (thetam * ipi - twothird) + qsp_theta[idx][j][kc] += exp(-0.5 * tmp * tmp) * exp( -0.5 * tmp2 * tmp2 ) - norms[i][j][kc] += 1.0 + norms[idx][j][kc] += 1.0 elif ( - t == "sq_face_cap_trig_pris" + typ == "sq_face_cap_trig_pris" and not flag_yaxis - and thetak < self._params[i]["TA3"] + and thetak < self._params[idx]["TA3"] ): - if thetam < self._params[i]["TA3"]: + if thetam < self._params[idx]["TA3"]: tmp = ( - cos(self._params[i]["fac_AA1"] * phi2) - ** self._params[i]["exp_cos_AA1"] + cos(self._params[idx]["fac_AA1"] * phi2) + ** self._params[idx]["exp_cos_AA1"] ) - tmp2 = self._params[i]["IGW_TA1"] * ( - thetam * ipi - self._params[i]["TA1"] + tmp2 = self._params[idx]["IGW_TA1"] * ( + thetam * ipi - self._params[idx]["TA1"] ) else: tmp = ( cos( - self._params[i]["fac_AA2"] - * (phi2 + self._params[i]["shift_AA2"]) + self._params[idx]["fac_AA2"] + * (phi2 + self._params[idx]["shift_AA2"]) ) - ** self._params[i]["exp_cos_AA2"] + ** self._params[idx]["exp_cos_AA2"] ) - tmp2 = self._params[i]["IGW_TA2"] * ( - thetam * ipi - self._params[i]["TA2"] + tmp2 = self._params[idx]["IGW_TA2"] * ( + thetam * ipi - self._params[idx]["TA2"] ) - qsp_theta[i][j][kc] += tmp * exp(-0.5 * tmp2 * tmp2) - norms[i][j][kc] += 1 + qsp_theta[idx][j][kc] += tmp * exp(-0.5 * tmp2 * tmp2) + norms[idx][j][kc] += 1 kc += 1 # Normalize Peters-style OPs. - for i, t in enumerate(self._types): - if t in [ + for idx, typ in enumerate(self._types): + if typ in [ "tri_plan", "tet", "bent", @@ -3263,12 +3270,12 @@ def get_order_parameters( "cuboct", "pent_plan", ]: - ops[i] = tmp_norm = 0.0 + ops[idx] = tmp_norm = 0.0 for j in range(nneigh): - ops[i] += sum(qsp_theta[i][j]) - tmp_norm += float(sum(norms[i][j])) - ops[i] = ops[i] / tmp_norm if tmp_norm > 1.0e-12 else None # type: ignore - elif t in [ + ops[idx] += sum(qsp_theta[idx][j]) + tmp_norm += float(sum(norms[idx][j])) + ops[idx] = ops[idx] / tmp_norm if tmp_norm > 1.0e-12 else None # type: ignore + elif typ in [ "T", "tri_pyr", "see_saw_rect", @@ -3288,35 +3295,35 @@ def get_order_parameters( "hex_plan_max", "sq_face_cap_trig_pris", ]: - ops[i] = None # type: ignore + ops[idx] = None # type: ignore if nneigh > 1: for j in range(nneigh): - for k in range(len(qsp_theta[i][j])): - qsp_theta[i][j][k] = ( - qsp_theta[i][j][k] / norms[i][j][k] if norms[i][j][k] > 1.0e-12 else 0.0 + for k in range(len(qsp_theta[idx][j])): + qsp_theta[idx][j][k] = ( + qsp_theta[idx][j][k] / norms[idx][j][k] if norms[idx][j][k] > 1.0e-12 else 0.0 ) - ops[i] = max(qsp_theta[i][j]) if j == 0 else max(ops[i], *qsp_theta[i][j]) - elif t == "bcc": - ops[i] = 0.0 + ops[idx] = max(qsp_theta[idx][j]) if j == 0 else max(ops[idx], *qsp_theta[idx][j]) + elif typ == "bcc": + ops[idx] = 0.0 for j in range(nneigh): - ops[i] += sum(qsp_theta[i][j]) + ops[idx] += sum(qsp_theta[idx][j]) if nneigh > 3: - ops[i] = ops[i] / float(0.5 * float(nneigh * (6 + (nneigh - 2) * (nneigh - 3)))) + ops[idx] = ops[idx] / float(0.5 * float(nneigh * (6 + (nneigh - 2) * (nneigh - 3)))) else: - ops[i] = None # type: ignore - elif t == "sq_pyr_legacy": + ops[idx] = None # type: ignore + elif typ == "sq_pyr_legacy": if nneigh > 1: dmean = np.mean(dist) acc = 0.0 for d in dist: - tmp = self._params[i][2] * (d - dmean) + tmp = self._params[idx][2] * (d - dmean) acc = acc + exp(-0.5 * tmp * tmp) for j in range(nneigh): - ops[i] = max(qsp_theta[i][j]) if j == 0 else max(ops[i], *qsp_theta[i][j]) - ops[i] = acc * ops[i] / float(nneigh) + ops[idx] = max(qsp_theta[idx][j]) if j == 0 else max(ops[idx], *qsp_theta[idx][j]) + ops[idx] = acc * ops[idx] / float(nneigh) # nneigh * (nneigh - 1)) else: - ops[i] = None # type: ignore + ops[idx] = None # type: ignore # Then, deal with the new-style OPs that require vectors between # neighbors. @@ -3338,20 +3345,20 @@ def get_order_parameters( b = min(distjk_unique) if len(distjk_unique) > 0 else 0 dhalf = max(distjk_unique) / 2 if len(distjk_unique) > 0 else 0 - for i, t in enumerate(self._types): - if t in ("reg_tri", "sq"): + for idx, typ in enumerate(self._types): + if typ in ("reg_tri", "sq"): if nneigh < 3: - ops[i] = None # type: ignore + ops[idx] = None # type: ignore else: - ops[i] = 1.0 - if t == "reg_tri": + ops[idx] = 1.0 + if typ == "reg_tri": a = 2 * asin(b / (2 * sqrt(h * h + (b / (2 * cos(3 * pi / 18))) ** 2))) # type: ignore nmax = 3 - elif t == "sq": + elif typ == "sq": a = 2 * asin(b / (2 * sqrt(h * h + dhalf * dhalf))) # type: ignore nmax = 4 for j in range(min([nneigh, nmax])): - ops[i] = ops[i] * exp(-0.5 * ((aijs[j] - a) * self._params[i][0]) ** 2) + ops[idx] = ops[idx] * exp(-0.5 * ((aijs[j] - a) * self._params[idx][0]) ** 2) return ops @@ -3407,9 +3414,9 @@ def get_nn_info(self, structure: Structure, n: int): """ site = structure[n] neighs_dists = structure.get_neighbors(site, self.cutoff) - ds = sorted(i.nn_distance for i in neighs_dists) + ds = sorted(idx.nn_distance for idx in neighs_dists) - ns = [1 / ds[i] - 1 / ds[i + 1] for i in range(len(ds) - 1)] + ns = [1 / ds[idx] - 1 / ds[idx + 1] for idx in range(len(ds) - 1)] d_max = ds[ns.index(max(ns))] siw = [] @@ -3479,9 +3486,9 @@ def get_nn_info(self, structure: Structure, n: int): """ site = structure[n] neighs_dists = structure.get_neighbors(site, self.cutoff) - ds = sorted(i.nn_distance for i in neighs_dists) + ds = sorted(idx.nn_distance for idx in neighs_dists) - ns = [ds[i + 1] / ds[i] for i in range(len(ds) - 1)] + ns = [ds[idx + 1] / ds[idx] for idx in range(len(ds) - 1)] d_max = ds[ns.index(max(ns))] siw = [] @@ -3551,9 +3558,9 @@ def get_nn_info(self, structure: Structure, n: int): """ site = structure[n] neighs_dists = structure.get_neighbors(site, self.cutoff) - ds = sorted(i.nn_distance for i in neighs_dists) + ds = sorted(idx.nn_distance for idx in neighs_dists) - ns = [ds[i + 1] - ds[i] for i in range(len(ds) - 1)] + ns = [ds[idx + 1] - ds[idx] for idx in range(len(ds) - 1)] d_max = ds[ns.index(max(ns))] siw = [] diff --git a/pymatgen/analysis/molecule_matcher.py b/pymatgen/analysis/molecule_matcher.py index dbd2fb02fac..64e8dc922f3 100644 --- a/pymatgen/analysis/molecule_matcher.py +++ b/pymatgen/analysis/molecule_matcher.py @@ -154,8 +154,8 @@ def uniform_labels(self, mol1, mol2): if elements1 != elements2: continue vmol2 = openbabel.OBMol() - for i in label2: - vmol2.AddAtom(ob_mol2.GetAtom(i)) + for idx in label2: + vmol2.AddAtom(ob_mol2.GetAtom(idx)) aligner.SetTargetMol(vmol2) aligner.Align() rmsd = aligner.GetRMSD() @@ -255,11 +255,11 @@ def _inchi_labels(mol): inchi = match.group("inchi") label_text = match.group("labels") eq_atom_text = match.group("eq_atoms") - heavy_atom_labels = tuple(int(i) for i in label_text.replace(";", ",").split(",")) + heavy_atom_labels = tuple(int(idx) for idx in label_text.replace(";", ",").split(",")) eq_atoms = [] if eq_atom_text is not None: eq_tokens = re.findall(r"\(((?:[0-9]+,)+[0-9]+)\)", eq_atom_text.replace(";", ",")) - eq_atoms = tuple(tuple(int(i) for i in t.split(",")) for t in eq_tokens) + eq_atoms = tuple(tuple(int(idx) for idx in t.split(",")) for t in eq_tokens) return heavy_atom_labels, eq_atoms, inchi @staticmethod @@ -275,8 +275,8 @@ def _group_centroid(mol, ilabels, group_atoms): Centroid. Tuple (x, y, z) """ c1x, c1y, c1z = 0.0, 0.0, 0.0 - for i in group_atoms: - orig_idx = ilabels[i - 1] + for idx in group_atoms: + orig_idx = ilabels[idx - 1] oa1 = mol.GetAtom(orig_idx) c1x += float(oa1.x()) c1y += float(oa1.y()) @@ -309,8 +309,8 @@ def _virtual_molecule(self, mol, ilabels, eq_atoms): unique_atom_labels = sorted(all_atoms - non_unique_atoms) # try to align molecules using unique atoms - for i in unique_atom_labels: - orig_idx = ilabels[i - 1] + for idx in unique_atom_labels: + orig_idx = ilabels[idx - 1] oa1 = mol.GetAtom(orig_idx) a1 = vmol.NewAtom() a1.SetAtomicNum(oa1.GetAtomicNum()) @@ -321,8 +321,8 @@ def _virtual_molecule(self, mol, ilabels, eq_atoms): for symm in eq_atoms: c1x, c1y, c1z = self._group_centroid(mol, ilabels, symm) min_distance = float("inf") - for i in range(1, vmol.NumAtoms() + 1): - va = vmol.GetAtom(i) + for idx in range(1, vmol.NumAtoms() + 1): + va = vmol.GetAtom(idx) distance = math.sqrt((c1x - va.x()) ** 2 + (c1y - va.y()) ** 2 + (c1z - va.z()) ** 2) if distance < min_distance: min_distance = distance @@ -434,14 +434,14 @@ def _align_hydrogen_atoms(mol1, mol2, heavy_indices1, heavy_indices2): label2 = heavy_indices2 + tuple(hydrogen_atoms2) cmol1 = openbabel.OBMol() - for i in label1: - oa1 = mol1.GetAtom(i) + for idx in label1: + oa1 = mol1.GetAtom(idx) a1 = cmol1.NewAtom() a1.SetAtomicNum(oa1.GetAtomicNum()) a1.SetVector(oa1.GetVector()) cmol2 = openbabel.OBMol() - for i in label2: - oa2 = mol2.GetAtom(i) + for idx in label2: + oa2 = mol2.GetAtom(idx) a2 = cmol2.NewAtom() a2.SetAtomicNum(oa2.GetAtomicNum()) a2.SetVector(oa2.GetVector()) @@ -489,7 +489,7 @@ def _get_elements(mol, label): Returns: Elements. List of integers. """ - return [int(mol.GetAtom(i).GetAtomicNum()) for i in label] + return [int(mol.GetAtom(idx).GetAtomicNum()) for idx in label] def _is_molecule_linear(self, mol): """ @@ -505,8 +505,8 @@ def _is_molecule_linear(self, mol): return True a1 = mol.GetAtom(1) a2 = mol.GetAtom(2) - for i in range(3, mol.NumAtoms() + 1): - angle = float(mol.GetAtom(i).GetAngle(a2, a1)) + for idx in range(3, mol.NumAtoms() + 1): + angle = float(mol.GetAtom(idx).GetAngle(a2, a1)) if angle < 0.0: angle = -angle if angle > 90.0: @@ -633,14 +633,14 @@ def _calc_rms(mol1, mol2, clabel1, clabel2): ob_mol2 = BabelMolAdaptor(mol2).openbabel_mol cmol1 = openbabel.OBMol() - for i in clabel1: - oa1 = ob_mol1.GetAtom(i) + for idx in clabel1: + oa1 = ob_mol1.GetAtom(idx) a1 = cmol1.NewAtom() a1.SetAtomicNum(oa1.GetAtomicNum()) a1.SetVector(oa1.GetVector()) cmol2 = openbabel.OBMol() - for i in clabel2: - oa2 = ob_mol2.GetAtom(i) + for idx in clabel2: + oa2 = ob_mol2.GetAtom(idx) a2 = cmol2.NewAtom() a2.SetAtomicNum(oa2.GetAtomicNum()) a2.SetVector(oa2.GetVector()) @@ -691,7 +691,7 @@ def group_molecules(self, mol_list): group_indices.append(sorted(current_group)) group_indices.sort(key=lambda x: (len(x), -x[0]), reverse=True) - return [[mol_list[i] for i in g] for g in group_indices] + return [[mol_list[idx] for idx in g] for g in group_indices] def as_dict(self): """ @@ -925,7 +925,7 @@ def fit(self, p: Molecule, ignore_warning=False): """ inds, U, V, rmsd = self.match(p, ignore_warning=ignore_warning) - p_prime = Molecule.from_sites([p[i] for i in inds]) + p_prime = Molecule.from_sites([p[idx] for idx in inds]) for site in p_prime: site.coords = np.dot(site.coords, U) + V @@ -1016,7 +1016,7 @@ def fit(self, p: Molecule): inds, U, V, rmsd = self.match(p) # Translate and rotate `mol1` unto `mol2` using Kabsch algorithm. - p_prime = Molecule.from_sites([p[i] for i in inds]) + p_prime = Molecule.from_sites([p[idx] for idx in inds]) for site in p_prime: site.coords = np.dot(site.coords, U) + V @@ -1025,7 +1025,7 @@ def fit(self, p: Molecule): @staticmethod def permutations(p_atoms, p_centroid, p_weights, q_atoms, q_centroid, q_weights): """Generates two possible permutations of atom order. This method uses the principle component - of the inertia tensor to prealign the molecules and hungarian method to determine the order. + of the inertia tensor to pre-align the molecules and hungarian method to determine the order. There are always two possible permutation depending on the way to pre-aligning the molecules. Args: @@ -1204,7 +1204,7 @@ def match(self, p: Molecule): out = [] for inds in self.permutations(p): p_prime = p.copy() - p_prime._sites = [p_prime[i] for i in inds] + p_prime._sites = [p_prime[idx] for idx in inds] U, V, rmsd = super().match(p_prime) @@ -1220,14 +1220,14 @@ def fit(self, p: Molecule): p: a `Molecule` object what will be matched with the target one. Returns: - Array of the possible matches where the elements are: + list[tuple[Molecule, float]]: possible matches where the elements are: p_prime: Rotated and translated of the `p` `Molecule` object rmsd: Root-mean-square-deviation between `p_prime` and the `target` """ - out = [] + out: list[tuple[Molecule, float]] = [] for inds in self.permutations(p): p_prime = p.copy() - p_prime._sites = [p_prime[i] for i in inds] + p_prime._sites = [p_prime[idx] for idx in inds] U, V, rmsd = super().match(p_prime) @@ -1236,7 +1236,7 @@ def fit(self, p: Molecule): for site in p_prime: site.coords = np.dot(site.coords, U) + V - out.append((p_prime, rmsd)) + out += [(p_prime, rmsd)] return out diff --git a/pymatgen/analysis/molecule_structure_comparator.py b/pymatgen/analysis/molecule_structure_comparator.py index 1a74e992c88..541cfdc1fc2 100644 --- a/pymatgen/analysis/molecule_structure_comparator.py +++ b/pymatgen/analysis/molecule_structure_comparator.py @@ -216,7 +216,7 @@ def _get_bonds(self, mol): num_atoms = len(mol) # index starting from 0 if self.ignore_ionic_bond: - covalent_atoms = [i for i in range(num_atoms) if mol.species[i].symbol not in self.ionic_element_list] + covalent_atoms = [idx for idx in range(num_atoms) if mol.species[idx].symbol not in self.ionic_element_list] else: covalent_atoms = list(range(num_atoms)) all_pairs = list(itertools.combinations(covalent_atoms, 2)) diff --git a/pymatgen/analysis/phase_diagram.py b/pymatgen/analysis/phase_diagram.py index 8fab50fb146..81345ff178d 100644 --- a/pymatgen/analysis/phase_diagram.py +++ b/pymatgen/analysis/phase_diagram.py @@ -391,7 +391,7 @@ def __init__( self.el_refs = dict(computed_data["el_refs"]) self.qhull_entries = tuple(computed_data["qhull_entries"]) self._qhull_spaces = tuple(frozenset(e.elements) for e in self.qhull_entries) - self._stable_entries = tuple({self.qhull_entries[i] for i in set(itertools.chain(*self.facets))}) + self._stable_entries = tuple({self.qhull_entries[idx] for idx in set(itertools.chain(*self.facets))}) self._stable_spaces = tuple(frozenset(e.elements) for e in self._stable_entries) def as_dict(self): @@ -458,7 +458,7 @@ def _compute(self) -> dict[str, Any]: # Add the elemental references idx.extend([min_entries.index(el) for el in el_refs.values()]) - qhull_entries = [min_entries[i] for i in idx] + qhull_entries = [min_entries[idx] for idx in idx] qhull_data = data[idx][:, 1:] # Add an extra point to enforce full dimensionality. @@ -652,8 +652,8 @@ def _get_facet_chempots(self, facet): Returns: {element: chempot} for all elements in the phase diagram. """ - comp_list = [self.qhull_entries[i].composition for i in facet] - energy_list = [self.qhull_entries[i].energy_per_atom for i in facet] + comp_list = [self.qhull_entries[idx].composition for idx in facet] + energy_list = [self.qhull_entries[idx].energy_per_atom for idx in facet] atom_frac_mat = [[c.get_atomic_fraction(e) for e in self.elements] for c in comp_list] chempots = np.linalg.solve(atom_frac_mat, energy_list) @@ -1189,8 +1189,10 @@ def get_chempot_range_map( data2 = self.facets[combi[1]] common_ent_ind = set(data1).intersection(set(data2)) if len(common_ent_ind) == len(elements): - common_entries = [self.qhull_entries[i] for i in common_ent_ind] - data = np.array([[all_chempots[i][j] - el_energies[self.elements[j]] for j in inds] for i in combi]) + common_entries = [self.qhull_entries[idx] for idx in common_ent_ind] + data = np.array( + [[all_chempots[ii][jj] - el_energies[self.elements[jj]] for jj in inds] for ii in combi] + ) sim = Simplex(data) for entry in common_entries: chempot_ranges[entry].append(sim) @@ -1647,7 +1649,7 @@ def __init__( # Add the elemental references inds.extend([min_entries.index(el) for el in el_refs.values()]) - self.qhull_entries = tuple(min_entries[i] for i in inds) + self.qhull_entries = tuple(min_entries[idx] for idx in inds) # make qhull spaces frozensets since they become keys to self.pds dict and frozensets are hashable # prevent repeating elements in chemical space and avoid the ordering problem (i.e. Fe-O == O-Fe automatically) self._qhull_spaces = tuple(frozenset(e.elements) for e in self.qhull_entries) @@ -1661,14 +1663,14 @@ def __init__( systems = set() # NOTE reduce the number of comparisons by only comparing to larger sets - for i in range(2, max_size + 1): - test = (s for s in spaces if len(s) == i) - refer = (s for s in spaces if len(s) > i) + for idx in range(2, max_size + 1): + test = (s for s in spaces if len(s) == idx) + refer = (s for s in spaces if len(s) > idx) systems |= {t for t in test if not any(t.issubset(r) for r in refer)} spaces = systems - # TODO comprhys: refactor to have self._compute method to allow serialisation + # TODO comprhys: refactor to have self._compute method to allow serialization self.spaces = sorted(spaces, key=len, reverse=False) # Calculate pds for smaller dimension spaces first self.pds = dict(self._get_pd_patch_for_space(s) for s in tqdm(self.spaces, disable=not verbose)) self.all_entries = all_entries @@ -1930,7 +1932,7 @@ def fmt(fl): for facet in pd.facets: for face in itertools.combinations(facet, len(facet) - 1): - face_entries = [pd.qhull_entries[i] for i in face] + face_entries = [pd.qhull_entries[idx] for idx in face] if any(e.reduced_formula in terminal_formulas for e in face_entries): continue @@ -2455,19 +2457,19 @@ class (pymatgen.analysis.chempot_diagram). if not (is_x and is_y): if is_x: coords = sorted(coords, key=lambda c: c[1]) - for i in [0, -1]: - x = [min(xlim), coords[i][0]] - y = [coords[i][1], coords[i][1]] + for idx in [0, -1]: + x = [min(xlim), coords[idx][0]] + y = [coords[idx][1], coords[idx][1]] plt.plot(x, y, "k") center_x += min(xlim) - center_y += coords[i][1] + center_y += coords[idx][1] elif is_y: coords = sorted(coords, key=lambda c: c[0]) - for i in [0, -1]: - x = [coords[i][0], coords[i][0]] - y = [coords[i][1], min(ylim)] + for idx in [0, -1]: + x = [coords[idx][0], coords[idx][0]] + y = [coords[idx][1], min(ylim)] plt.plot(x, y, "k") - center_x += coords[i][0] + center_x += coords[idx][0] center_y += min(ylim) xy = (center_x / (n + 2), center_y / (n + 2)) else: diff --git a/pymatgen/command_line/bader_caller.py b/pymatgen/command_line/bader_caller.py index d737b4daa28..3fd69cf2235 100644 --- a/pymatgen/command_line/bader_caller.py +++ b/pymatgen/command_line/bader_caller.py @@ -203,7 +203,11 @@ def __init__( ) def _parse_atomic_densities(self): # Deprecation tracker - if datetime(2025, 2, 26) < datetime.now() and "CI" in os.environ: + if ( + datetime(2025, 2, 26) < datetime.now() + and os.getenv("CI") + and os.getenv("GITHUB_REPOSITORY") == "materialsproject/pymatgen" + ): raise RuntimeError("This method should have been removed, see #3656.") def slice_from_center(data, x_width, y_width, z_width): diff --git a/pymatgen/core/structure.py b/pymatgen/core/structure.py index d23deb80ea5..3dcca940673 100644 --- a/pymatgen/core/structure.py +++ b/pymatgen/core/structure.py @@ -1131,7 +1131,7 @@ def from_spacegroup( raise ValueError(f"Supplied species and coords lengths ({len(species)} vs {len(coords)}) are different!") frac_coords = ( - np.array(coords, dtype=np.float_) if not coords_are_cartesian else latt.get_fractional_coords(coords) + np.array(coords, dtype=np.float64) if not coords_are_cartesian else latt.get_fractional_coords(coords) ) props = {} if site_properties is None else site_properties @@ -1754,7 +1754,7 @@ def get_symmetric_neighbor_list( # expand the output tuple by symmetry_indices and symmetry_ops. n_bonds = len(bonds[0]) symmetry_indices = np.empty(n_bonds) - symmetry_indices[:] = np.NaN + symmetry_indices[:] = np.nan symmetry_ops = np.empty(len(symmetry_indices), dtype=object) symmetry_identity = SymmOp.from_rotation_and_translation(np.eye(3), np.zeros(3)) symmetry_index = 0 diff --git a/pymatgen/entries/computed_entries.py b/pymatgen/entries/computed_entries.py index d74d256b34a..263bbac1aa8 100644 --- a/pymatgen/entries/computed_entries.py +++ b/pymatgen/entries/computed_entries.py @@ -416,7 +416,7 @@ def normalize(self, mode: Literal["formula_unit", "atom"] = "formula_unit") -> C def __repr__(self) -> str: n_atoms = self.composition.num_atoms output = [ - f"{self.entry_id} {type(self).__name__:<10} " f"- {self.formula:<12} ({self.reduced_formula})", + f"{self.entry_id} {type(self).__name__:<10} - {self.formula:<12} ({self.reduced_formula})", f"{'Energy (Uncorrected)':<24} = {self._energy:<9.4f} eV ({self._energy / n_atoms:<8.4f} eV/atom)", f"{'Correction':<24} = {self.correction:<9.4f} eV ({self.correction / n_atoms:<8.4f} eV/atom)", f"{'Energy (Final)':<24} = {self.energy:<9.4f} eV ({self.energy_per_atom:<8.4f} eV/atom)", diff --git a/pymatgen/io/cif.py b/pymatgen/io/cif.py index 0ff60e099d2..182888ada6d 100644 --- a/pymatgen/io/cif.py +++ b/pymatgen/io/cif.py @@ -17,6 +17,7 @@ from typing import TYPE_CHECKING, Any, Literal import numpy as np +from monty.dev import deprecated from monty.io import zopen from monty.serialization import loadfn @@ -34,11 +35,11 @@ __author__ = "Shyue Ping Ong, Will Richards, Matthew Horton" -sub_spgrp = partial(re.sub, r"[\s_]", "") +sub_space_group = partial(re.sub, r"[\s_]", "") -space_groups = {sub_spgrp(key): key for key in SYMM_DATA["space_group_encoding"]} # type: ignore +space_groups = {sub_space_group(key): key for key in SYMM_DATA["space_group_encoding"]} # type: ignore -space_groups.update({sub_spgrp(key): key for key in SYMM_DATA["space_group_encoding"]}) # type: ignore +space_groups.update({sub_space_group(key): key for key in SYMM_DATA["space_group_encoding"]}) # type: ignore class CifBlock: @@ -711,7 +712,7 @@ def get_symops(self, data): msg_template = "No _symmetry_equiv_pos_as_xyz type key found. Spacegroup from {} used." if sg: - sg = sub_spgrp(sg) + sg = sub_space_group(sg) try: spg = space_groups.get(sg) if spg: @@ -1136,7 +1137,7 @@ def get_matching_coord(coord): return struct return None - @np.deprecate( + @deprecated( message="get_structures is deprecated and will be removed in 2024. Use parse_structures instead." "The only difference is that primitive defaults to False in the new parse_structures method." "So parse_structures(primitive=True) is equivalent to the old behavior of get_structures().", @@ -1184,8 +1185,12 @@ def parse_structures( Returns: list[Structure]: All structures in CIF file. """ - if os.getenv("CI") and datetime.now() > datetime(2024, 3, 1): # March 2024 seems long enough # pragma: no cover - raise RuntimeError("remove the change of default primitive=True to False made on 2023-10-24") + if ( + os.getenv("CI") + and os.getenv("GITHUB_REPOSITORY") == "materialsproject/pymatgen" + and datetime.now() > datetime(2024, 10, 1) + ): # pragma: no cover + raise RuntimeError("remove the warning about changing default primitive=True to False on 2023-10-24") if primitive is None: primitive = False warnings.warn( diff --git a/pymatgen/io/cp2k/utils.py b/pymatgen/io/cp2k/utils.py index 8e49ab9e609..430946e5e41 100644 --- a/pymatgen/io/cp2k/utils.py +++ b/pymatgen/io/cp2k/utils.py @@ -46,7 +46,7 @@ def postprocessor(data: str) -> str | float | bool | None: except ValueError as exc: raise ValueError(f"Error parsing {data!r} as float in CP2K file.") from exc if re.match(r"\*+", data): - return np.NaN + return np.nan return data diff --git a/pymatgen/io/lobster/outputs.py b/pymatgen/io/lobster/outputs.py index a670a993cd1..ad284bf0a34 100644 --- a/pymatgen/io/lobster/outputs.py +++ b/pymatgen/io/lobster/outputs.py @@ -1003,9 +1003,9 @@ def _get_spillings(data, number_of_spins): splitrow = row.split() if len(splitrow) > 2 and splitrow[2] == "spilling:": if splitrow[1] == "charge": - charge_spilling.append(np.float_(splitrow[3].replace("%", "")) / 100.0) + charge_spilling.append(np.float64(splitrow[3].replace("%", "")) / 100.0) if splitrow[1] == "total": - total_spilling.append(np.float_(splitrow[3].replace("%", "")) / 100.0) + total_spilling.append(np.float64(splitrow[3].replace("%", "")) / 100.0) if len(charge_spilling) == number_of_spins and len(total_spilling) == number_of_spins: break diff --git a/pymatgen/io/vasp/optics.py b/pymatgen/io/vasp/optics.py index dc7f71f0f91..efc17ad6ada 100644 --- a/pymatgen/io/vasp/optics.py +++ b/pymatgen/io/vasp/optics.py @@ -388,7 +388,7 @@ def epsilon_imag( min_band1, max_band1 = np.min(np.where(cderm)[1]), np.max(np.where(cderm)[1]) except ValueError as exc: if "zero-size array" in str(exc): - return egrid, np.zeros_like(egrid, dtype=np.complex_) + return egrid, np.zeros_like(egrid, dtype=np.complex128) raise exc _, _, nk, nspin = cderm.shape[:4] iter_idx = [ diff --git a/pymatgen/io/vasp/outputs.py b/pymatgen/io/vasp/outputs.py index 23fceb0cbd4..2c7489986f9 100644 --- a/pymatgen/io/vasp/outputs.py +++ b/pymatgen/io/vasp/outputs.py @@ -3792,7 +3792,7 @@ def __init__(self, filename): phase_factors = defaultdict( lambda: np.full( (n_kpoints, n_bands, n_ions, len(headers)), - np.NaN, + np.nan, dtype=np.complex128, ) ) @@ -4731,7 +4731,7 @@ def fft_mesh(self, kpoint: int, band: int, spin: int = 0, spinor: int = 0, shift else: tcoeffs = self.coeffs[kpoint][band] - mesh = np.zeros(tuple(self.ng), dtype=np.complex_) + mesh = np.zeros(tuple(self.ng), dtype=np.complex128) for gp, coeff in zip(self.Gpoints[kpoint], tcoeffs): t = tuple(gp.astype(int) + (self.ng / 2).astype(int)) mesh[t] = coeff @@ -5081,8 +5081,8 @@ def read_data(dtype): return np.frombuffer(data, dtype=dtype) nbands, nelect, nk, ispin = read_data(np.int32) - _ = read_data(np.float_) # nodes_in_dielectric_function - _ = read_data(np.float_) # wplasmon + _ = read_data(np.float64) # nodes_in_dielectric_function + _ = read_data(np.float64) # wplasmon me_datatype = np.dtype(data_type) cder = read_data(me_datatype) diff --git a/pymatgen/io/vasp/sets.py b/pymatgen/io/vasp/sets.py index ba4713acf17..165fb2f6b12 100644 --- a/pymatgen/io/vasp/sets.py +++ b/pymatgen/io/vasp/sets.py @@ -1814,9 +1814,8 @@ def __post_init__(self): and not isinstance(self.structure[0].magmom, list) ): raise ValueError( - "The structure must have the 'magmom' site " - "property and each magnetic moment value must have 3 " - "components. eg:- magmom = [0,0,2]" + "The structure must have the 'magmom' site property and each magnetic " + "moment value must have 3 components. e.g. magmom = [0,0,2]" ) @property @@ -1866,6 +1865,7 @@ def structure(self, structure: Structure | None) -> None: # magmom has to be 3D for SOC calculation. if hasattr(structure[0], "magmom"): if not isinstance(structure[0].magmom, list): + # project magmom to z-axis structure = structure.copy(site_properties={"magmom": [[0, 0, site.magmom] for site in structure]}) else: raise ValueError("Neither the previous structure has magmom property nor magmom provided") @@ -2755,7 +2755,7 @@ def get_vasprun_outcar(path: str | Path, parse_dos: bool = True, parse_eigen: bo ) -def get_structure_from_prev_run(vasprun, outcar=None): +def get_structure_from_prev_run(vasprun, outcar=None) -> Structure: """ Process structure from previous run. @@ -2764,8 +2764,8 @@ def get_structure_from_prev_run(vasprun, outcar=None): outcar (Outcar): Outcar that contains the magnetization info from previous run. Returns: - The magmom-decorated structure that can be passed to get VASP input files, e.g. - get_kpoints. + Structure: The magmom-decorated structure that can be passed to get VASP input files, e.g. + get_kpoints(). """ structure = vasprun.final_structure @@ -2776,7 +2776,7 @@ def get_structure_from_prev_run(vasprun, outcar=None): site_properties["magmom"] = [i["tot"] for i in outcar.magnetization] else: site_properties["magmom"] = vasprun.parameters["MAGMOM"] - # ldau + # LDAU if vasprun.parameters.get("LDAU", False): for k in ("LDAUU", "LDAUJ", "LDAUL"): vals = vasprun.incar[k] diff --git a/pymatgen/symmetry/analyzer.py b/pymatgen/symmetry/analyzer.py index f6420a1f71e..d692f226c2d 100644 --- a/pymatgen/symmetry/analyzer.py +++ b/pymatgen/symmetry/analyzer.py @@ -459,17 +459,17 @@ def get_conventional_to_primitive_transformation_matrix(self, international_mono if abs(lengths[0] - lengths[2]) < 0.0001: transf = np.eye else: - transf = np.array([[-1, 1, 1], [2, 1, 1], [-1, -2, 1]], dtype=np.float_) / 3 + transf = np.array([[-1, 1, 1], [2, 1, 1], [-1, -2, 1]], dtype=np.float64) / 3 elif "I" in self.get_space_group_symbol(): - transf = np.array([[-1, 1, 1], [1, -1, 1], [1, 1, -1]], dtype=np.float_) / 2 + transf = np.array([[-1, 1, 1], [1, -1, 1], [1, 1, -1]], dtype=np.float64) / 2 elif "F" in self.get_space_group_symbol(): - transf = np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]], dtype=np.float_) / 2 + transf = np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]], dtype=np.float64) / 2 elif "C" in self.get_space_group_symbol() or "A" in self.get_space_group_symbol(): if self.get_crystal_system() == "monoclinic": - transf = np.array([[1, 1, 0], [-1, 1, 0], [0, 0, 2]], dtype=np.float_) / 2 + transf = np.array([[1, 1, 0], [-1, 1, 0], [0, 0, 2]], dtype=np.float64) / 2 else: - transf = np.array([[1, -1, 0], [1, 1, 0], [0, 0, 2]], dtype=np.float_) / 2 + transf = np.array([[1, -1, 0], [1, 1, 0], [0, 0, 2]], dtype=np.float64) / 2 else: transf = np.eye(3) diff --git a/pymatgen/transformations/advanced_transformations.py b/pymatgen/transformations/advanced_transformations.py index 180a0f2917e..610c4098bad 100644 --- a/pymatgen/transformations/advanced_transformations.py +++ b/pymatgen/transformations/advanced_transformations.py @@ -1518,8 +1518,8 @@ def __init__( please use max_atoms angle_tolerance: tolerance to determine the 90 degree angles. """ - self.min_atoms = min_atoms or -np.Inf - self.max_atoms = max_atoms or np.Inf + self.min_atoms = min_atoms or -np.inf + self.max_atoms = max_atoms or np.inf self.min_length = min_length self.force_diagonal = force_diagonal self.force_90_degrees = force_90_degrees diff --git a/pymatgen/vis/structure_vtk.py b/pymatgen/vis/structure_vtk.py index 30a6deedb19..33d6f28aad8 100644 --- a/pymatgen/vis/structure_vtk.py +++ b/pymatgen/vis/structure_vtk.py @@ -601,7 +601,7 @@ def add_faces(self, faces, color, opacity=0.35): center = np.zeros(3, float) for site in face: center += site - center /= np.float_(len(face)) + center /= np.float64(len(face)) for ii, f in enumerate(face): points = vtk.vtkPoints() triangle = vtk.vtkTriangle() diff --git a/tests/analysis/test_interface_reactions.py b/tests/analysis/test_interface_reactions.py index 3679c7f9504..361a35bc214 100644 --- a/tests/analysis/test_interface_reactions.py +++ b/tests/analysis/test_interface_reactions.py @@ -7,7 +7,7 @@ from matplotlib.figure import Figure as mpl_figure from numpy.testing import assert_allclose from pandas import DataFrame -from plotly.graph_objects import Figure as plotly_figure +from plotly.graph_objects import Figure from scipy.spatial import ConvexHull from pymatgen.analysis.interface_reactions import GrandPotentialInterfacialReactivity, InterfacialReactivity @@ -391,7 +391,7 @@ def test_plot(self): assert fig, isinstance(fig, mpl_figure) fig = ir.plot(backend="plotly") - assert isinstance(fig, plotly_figure) + assert isinstance(fig, Figure) def test_get_dataframe(self): for ir in self.irs: @@ -415,14 +415,8 @@ def test_minimum(self): (0.3333333, -3.333333), (0.3333333, -4.0), ] - for i, j in zip(self.irs, answer): - ( - assert_allclose(i.minimum, j, atol=1e-7), - ( - f"minimum: the system with {i.c1_original.reduced_formula} and {i.c2_original.reduced_formula} " - f"gets error!{j} expected, but gets {i.minimum}" - ), - ) + for inter_react, expected in zip(self.irs, answer): + assert_allclose(inter_react.minimum, expected, atol=1e-7) def test_get_no_mixing_energy(self): answer = [ diff --git a/tests/command_line/test_bader_caller.py b/tests/command_line/test_bader_caller.py index 1abab791cd4..d8a298352eb 100644 --- a/tests/command_line/test_bader_caller.py +++ b/tests/command_line/test_bader_caller.py @@ -48,8 +48,8 @@ def test_init(self): 1.021523, 1.024357, ] - for i in range(14): - assert ans[i] == approx(analysis.get_charge_transfer(i), abs=1e-3) + for idx in range(14): + assert ans[idx] == approx(analysis.get_charge_transfer(idx), abs=1e-3) assert analysis.get_partial_charge(0) == -analysis.get_charge_transfer(0) struct = analysis.get_oxidation_state_decorated_structure() assert struct[0].specie.oxi_state == approx(1.3863218, abs=1e-3) diff --git a/tests/io/vasp/test_outputs.py b/tests/io/vasp/test_outputs.py index 19a05340ee4..019298fb32f 100644 --- a/tests/io/vasp/test_outputs.py +++ b/tests/io/vasp/test_outputs.py @@ -1767,13 +1767,13 @@ def test_n2_spin(self): assert len(w.coeffs[0][0]) == w.nb assert len(w.band_energy[0]) == w.nk - temp_ggp = Wavecar._generate_G_points + orig_gen_g_points = Wavecar._generate_G_points try: - Wavecar._generate_G_points = lambda x, y, gamma: [] + Wavecar._generate_G_points = lambda _x, _y, gamma: [] with pytest.raises(ValueError, match=r"not enough values to unpack \(expected 3, got 0\)"): Wavecar(f"{TEST_FILES_DIR}/WAVECAR.N2") finally: - Wavecar._generate_G_points = temp_ggp + Wavecar._generate_G_points = orig_gen_g_points def test_generate_nbmax(self): self.wavecar._generate_nbmax()