diff --git a/pymatgen/command_line/chargemol_caller.py b/pymatgen/command_line/chargemol_caller.py index 82bd19c25ca..4f9a2cf694d 100644 --- a/pymatgen/command_line/chargemol_caller.py +++ b/pymatgen/command_line/chargemol_caller.py @@ -42,12 +42,15 @@ from __future__ import annotations +import hashlib import os import subprocess +import urllib.request import warnings from glob import glob from shutil import which from typing import TYPE_CHECKING +from zipfile import ZipFile import numpy as np from monty.tempfile import ScratchDir @@ -101,9 +104,18 @@ def __init__( " Please download the library at https://sourceforge.net/projects/ddec/files" "and follow the instructions." ) - if atomic_densities_path == "": - atomic_densities_path = os.getcwd() - self._atomic_densities_path = atomic_densities_path + self._atomic_densities_path = atomic_densities_path or os.getenv("DDEC6_ATOMIC_DENSITIES_DIR") or os.getcwd() + + if not os.path.isdir(self._atomic_densities_path): + try: + self._download_and_unzip_atomic_densities() + except Exception as exc: + raise FileNotFoundError( + "Atomic densities not found and auto-download failed. " + "Please set the DDEC6_ATOMIC_DENSITIES_DIR environment variable to the " + "directory containing the atomic densities from Chargemol at " + "https://sourceforge.net/projects/ddec/files" + ) from exc self._chgcar_path = self._get_filepath(path, "CHGCAR") self._potcar_path = self._get_filepath(path, "POTCAR") @@ -132,6 +144,30 @@ def __init__( else: self._from_data_dir(chargemol_output_path=path) + def _download_and_unzip_atomic_densities(self, version: str = "latest", verbose: bool = True) -> None: + download_url = f"https://sourceforge.net/projects/ddec/files/{version}/download" + if verbose: + print(f"Downloading atomic densities from {download_url}...") + download_path = f"/tmp/chargemol_{version}.zip" + extraction_path = "~/.cache/pymatgen/ddec" + + urllib.request.urlretrieve(download_url, download_path) + + sha1 = "3246a2f3ff7c6f88a8ef0796bde9f73b1da78e32" + # verify sha1 + with open(download_path, "rb") as f: + calculated_sha1 = hashlib.sha1(f.read()).hexdigest() + print(f"Calculated SHA-1: {calculated_sha1}") + assert calculated_sha1 == sha1 + + with ZipFile(download_path, "r") as zip_ref: + # extract to only the atomic_densities directory + zip_ref.extractall(extraction_path, members=["atomic_densities/*"]) + + self._atomic_densities_path = os.path.expanduser(extraction_path) + + os.remove(download_path) # clean up downloaded file + @staticmethod def _get_filepath(path, filename, suffix=""): """Returns the full path to the filename in the path. Works even if the file has @@ -163,10 +199,9 @@ def _execute_chargemol(self, **job_control_kwargs): """Internal function to run Chargemol. Args: - atomic_densities_path (str): Path to the atomic densities directory - required by Chargemol. If None, Pymatgen assumes that this is - defined in a "DDEC6_ATOMIC_DENSITIES_DIR" environment variable. - Default: None. + atomic_densities_path (str): Path to the atomic densities directory required by Chargemol. + If None, Pymatgen assumes that this is defined in a "DDEC6_ATOMIC_DENSITIES_DIR" + environment variable. Default: None. job_control_kwargs: Keyword arguments for _write_jobscript_for_chargemol. """ with ScratchDir("."): @@ -531,16 +566,15 @@ def _get_data_from_xyz(xyz_path): list[float]: site-specific properties """ props = [] - if os.path.exists(xyz_path): - with open(xyz_path) as r: - for i, line in enumerate(r): - if i <= 1: - continue - if line.strip() == "": - break - props.append(float(line.split()[-1])) - else: + if not os.path.exists(xyz_path): raise FileNotFoundError(f"{xyz_path} not found") + with open(xyz_path) as file: + for idx, line in enumerate(file): + if idx <= 1: + continue + if line.strip() == "": + break + props.append(float(line.split()[-1])) return props diff --git a/pymatgen/core/structure.py b/pymatgen/core/structure.py index 3efa3902025..1fb0779fc35 100644 --- a/pymatgen/core/structure.py +++ b/pymatgen/core/structure.py @@ -493,7 +493,9 @@ def remove_site_property(self, property_name: str): del site.properties[property_name] def replace_species( - self, species_mapping: dict[SpeciesLike, SpeciesLike | dict[SpeciesLike, float]], in_place: bool = True + self, + species_mapping: dict[SpeciesLike, SpeciesLike | dict[SpeciesLike, float]], + in_place: bool = True, ) -> SiteCollection: """Swap species. @@ -599,7 +601,11 @@ def add_spin_by_element(self, spins: dict[str, float]) -> None: for sp, occu in site.species.items(): sym = sp.symbol oxi_state = getattr(sp, "oxi_state", None) - species = Species(sym, oxidation_state=oxi_state, spin=spins.get(str(sp), spins.get(sym))) + species = Species( + sym, + oxidation_state=oxi_state, + spin=spins.get(str(sp), spins.get(sym)), + ) new_species[species] = occu site.species = Composition(new_species) @@ -730,7 +736,10 @@ def _relax( opt_kwargs = opt_kwargs or {} is_molecule = isinstance(self, Molecule) # UIP=universal interatomic potential - run_uip = isinstance(calculator, str) and calculator.lower() in ("m3gnet", "chgnet") + run_uip = isinstance(calculator, str) and calculator.lower() in ( + "m3gnet", + "chgnet", + ) calc_params = dict(stress_weight=stress_weight) if not is_molecule else {} calculator = self._prep_calculator(calculator, **calc_params) @@ -1198,7 +1207,13 @@ def from_magnetic_spacegroup( all_site_properties["magmom"] = all_magmoms - return cls(latt, all_sp, all_coords, site_properties=all_site_properties, labels=all_labels) + return cls( + latt, + all_sp, + all_coords, + site_properties=all_site_properties, + labels=all_labels, + ) def unset_charge(self): """Reset the charge to None. E.g. to compute it dynamically based on oxidation states.""" @@ -1281,7 +1296,10 @@ def get_space_group_info(self, symprec: float = 1e-2, angle_tolerance: float = 5 from pymatgen.symmetry.analyzer import SpacegroupAnalyzer spg_analyzer = SpacegroupAnalyzer(self, symprec=symprec, angle_tolerance=angle_tolerance) - return spg_analyzer.get_space_group_symbol(), spg_analyzer.get_space_group_number() + return ( + spg_analyzer.get_space_group_symbol(), + spg_analyzer.get_space_group_number(), + ) def matches(self, other: IStructure | Structure, anonymous: bool = False, **kwargs) -> bool: """Check whether this structure is similar to another structure. @@ -1609,7 +1627,12 @@ def get_neighbor_list( if exclude_self: self_pair = (center_indices == points_indices) & (distances <= numerical_tol) cond = ~self_pair - return (center_indices[cond], points_indices[cond], images[cond], distances[cond]) + return ( + center_indices[cond], + points_indices[cond], + images[cond], + distances[cond], + ) def get_symmetric_neighbor_list( self, @@ -1695,7 +1718,12 @@ def get_symmetric_neighbor_list( # delete the redundant neighbors m = ~np.in1d(np.arange(len(bonds[0])), redundant) idcs_dist = np.argsort(bonds[3][m]) - bonds = (bonds[0][m][idcs_dist], bonds[1][m][idcs_dist], bonds[2][m][idcs_dist], bonds[3][m][idcs_dist]) + bonds = ( + bonds[0][m][idcs_dist], + bonds[1][m][idcs_dist], + bonds[2][m][idcs_dist], + bonds[3][m][idcs_dist], + ) # expand the output tuple by symmetry_indices and symmetry_ops. n_bonds = len(bonds[0]) @@ -1725,23 +1753,32 @@ def get_symmetric_neighbor_list( to_b = self[bonds[1][jdx]].frac_coords r_b = bonds[2][jdx] for op in ops: - are_related, is_reversed = op.are_symmetrically_related_vectors( - from_a, to_a, r_a, from_b, to_b, r_b - ) + ( + are_related, + is_reversed, + ) = op.are_symmetrically_related_vectors(from_a, to_a, r_a, from_b, to_b, r_b) if are_related and not is_reversed: symmetry_indices[jdx] = symmetry_index symmetry_ops[jdx] = op elif are_related and is_reversed: symmetry_indices[jdx] = symmetry_index symmetry_ops[jdx] = op - bonds[0][jdx], bonds[1][jdx] = bonds[1][jdx], bonds[0][jdx] + bonds[0][jdx], bonds[1][jdx] = ( + bonds[1][jdx], + bonds[0][jdx], + ) bonds[2][jdx] = -bonds[2][jdx] symmetry_index += 1 # the bonds are ordered by their symmetry index idcs_symid = np.argsort(symmetry_indices) - bonds = (bonds[0][idcs_symid], bonds[1][idcs_symid], bonds[2][idcs_symid], bonds[3][idcs_symid]) + bonds = ( + bonds[0][idcs_symid], + bonds[1][idcs_symid], + bonds[2][idcs_symid], + bonds[3][idcs_symid], + ) symmetry_indices = symmetry_indices[idcs_symid] symmetry_ops = symmetry_ops[idcs_symid] @@ -1754,7 +1791,10 @@ def get_symmetric_neighbor_list( first_idx = np.argmax(symmetry_indices == symmetry_idx) for second_idx in identity_idcs: if symmetry_indices[second_idx] == symmetry_idx: - idcs_symop[first_idx], idcs_symop[second_idx] = idcs_symop[second_idx], idcs_symop[first_idx] + idcs_symop[first_idx], idcs_symop[second_idx] = ( + idcs_symop[second_idx], + idcs_symop[first_idx], + ) return ( bonds[0][idcs_symop], @@ -2008,7 +2048,12 @@ def get_all_neighbors_old(self, r, include_index=False, include_image=False, inc return neighbors def get_neighbors_in_shell( - self, origin: ArrayLike, r: float, dr: float, include_index: bool = False, include_image: bool = False + self, + origin: ArrayLike, + r: float, + dr: float, + include_index: bool = False, + include_image: bool = False, ) -> list[PeriodicNeighbor]: """Returns all sites in a shell centered on origin (coords) between radii r-dr and r+dr. @@ -2233,7 +2278,13 @@ def interpolate( lattice = self.lattice fcoords = start_coords + x * vec structs.append( - self.__class__(lattice, sp, fcoords, site_properties=self.site_properties, labels=self.labels) + self.__class__( + lattice, + sp, + fcoords, + site_properties=self.site_properties, + labels=self.labels, + ) ) return structs @@ -2264,7 +2315,10 @@ def get_miller_index_from_site_indexes(self, site_ids, round_dp=4, verbose=True) ) def get_primitive_structure( - self, tolerance: float = 0.25, use_site_props: bool = False, constrain_latt: list | dict | None = None + self, + tolerance: float = 0.25, + use_site_props: bool = False, + constrain_latt: list | dict | None = None, ): """This finds a smaller unit cell than the input. Sometimes it doesn"t find the smallest possible one, so this method is recursively called @@ -2449,7 +2503,9 @@ def factors(n: int): # Default behavior p = struct.get_primitive_structure( - tolerance=tolerance, use_site_props=use_site_props, constrain_latt=constrain_latt + tolerance=tolerance, + use_site_props=use_site_props, + constrain_latt=constrain_latt, ).get_reduced_structure() if not constrain_latt: return p @@ -2662,13 +2718,11 @@ def to(self, filename: str = "", fmt: str = "", **kwargs) -> str: """ fmt = fmt.lower() - if fmt == "cif" or fnmatch(filename.lower(), "*.cif*"): - from pymatgen.io.cif import CifWriter + from pymatgen.io.cif import CifWriter + if fmt == "cif" or fnmatch(filename.lower(), "*.cif*"): writer = CifWriter(self, **kwargs) elif fmt == "mcif" or fnmatch(filename.lower(), "*.mcif*"): - from pymatgen.io.cif import CifWriter - writer = CifWriter(self, write_magmoms=True, **kwargs) elif fmt == "poscar" or fnmatch(filename, "*POSCAR*"): from pymatgen.io.vasp import Poscar @@ -2822,7 +2876,12 @@ def from_str( @classmethod def from_file( - cls, filename: str | Path, primitive: bool = False, sort: bool = False, merge_tol: float = 0.0, **kwargs + cls, + filename: str | Path, + primitive: bool = False, + sort: bool = False, + merge_tol: float = 0.0, + **kwargs, ) -> Structure | IStructure: """Reads a structure from a file. For example, anything ending in a "cif" is assumed to be a Crystallographic Information Format file. @@ -2859,26 +2918,75 @@ def from_file( with zopen(filename, "rt", errors="replace") as f: contents = f.read() if fnmatch(fname.lower(), "*.cif*") or fnmatch(fname.lower(), "*.mcif*"): - return cls.from_str(contents, fmt="cif", primitive=primitive, sort=sort, merge_tol=merge_tol, **kwargs) + return cls.from_str( + contents, + fmt="cif", + primitive=primitive, + sort=sort, + merge_tol=merge_tol, + **kwargs, + ) if fnmatch(fname, "*POSCAR*") or fnmatch(fname, "*CONTCAR*") or fnmatch(fname, "*.vasp"): - struct = cls.from_str(contents, fmt="poscar", primitive=primitive, sort=sort, merge_tol=merge_tol, **kwargs) + struct = cls.from_str( + contents, + fmt="poscar", + primitive=primitive, + sort=sort, + merge_tol=merge_tol, + **kwargs, + ) elif fnmatch(fname, "CHGCAR*") or fnmatch(fname, "LOCPOT*"): struct = Chgcar.from_file(filename, **kwargs).structure elif fnmatch(fname, "vasprun*.xml*"): struct = Vasprun(filename, **kwargs).final_structure elif fnmatch(fname.lower(), "*.cssr*"): - return cls.from_str(contents, fmt="cssr", primitive=primitive, sort=sort, merge_tol=merge_tol, **kwargs) + return cls.from_str( + contents, + fmt="cssr", + primitive=primitive, + sort=sort, + merge_tol=merge_tol, + **kwargs, + ) elif fnmatch(fname, "*.json*") or fnmatch(fname, "*.mson*"): - return cls.from_str(contents, fmt="json", primitive=primitive, sort=sort, merge_tol=merge_tol, **kwargs) + return cls.from_str( + contents, + fmt="json", + primitive=primitive, + sort=sort, + merge_tol=merge_tol, + **kwargs, + ) elif fnmatch(fname, "*.yaml*") or fnmatch(fname, "*.yml*"): - return cls.from_str(contents, fmt="yaml", primitive=primitive, sort=sort, merge_tol=merge_tol, **kwargs) + return cls.from_str( + contents, + fmt="yaml", + primitive=primitive, + sort=sort, + merge_tol=merge_tol, + **kwargs, + ) elif fnmatch(fname, "*.xsf"): - return cls.from_str(contents, fmt="xsf", primitive=primitive, sort=sort, merge_tol=merge_tol, **kwargs) + return cls.from_str( + contents, + fmt="xsf", + primitive=primitive, + sort=sort, + merge_tol=merge_tol, + **kwargs, + ) elif fnmatch(fname, "input*.xml"): return ExcitingInput.from_file(fname, **kwargs).structure elif fnmatch(fname, "*rndstr.in*") or fnmatch(fname, "*lat.in*") or fnmatch(fname, "*bestsqs*"): - return cls.from_str(contents, fmt="mcsqs", primitive=primitive, sort=sort, merge_tol=merge_tol, **kwargs) + return cls.from_str( + contents, + fmt="mcsqs", + primitive=primitive, + sort=sort, + merge_tol=merge_tol, + **kwargs, + ) elif fnmatch(fname, "CTRL*"): return LMTOCtrl.from_file(filename=filename, **kwargs).structure elif fnmatch(fname, "inp*.xml") or fnmatch(fname, "*.in*") or fnmatch(fname, "inp_*"): @@ -3222,7 +3330,13 @@ def get_zmatrix(self): angle = self.get_angle(idx, nn[0], nn[1]) dih = self.get_dihedral(idx, nn[0], nn[1], nn[2]) output.append(f"{self[idx].specie} {nn[0] + 1} B{idx} {nn[1] + 1} A{idx} {nn[2] + 1} D{idx}") - output_var.extend((f"B{idx}={bond_length:.6f}", f"A{idx}={angle:.6f}", f"D{idx}={dih:.6f}")) + output_var.extend( + ( + f"B{idx}={bond_length:.6f}", + f"A{idx}={angle:.6f}", + f"D{idx}={dih:.6f}", + ) + ) return "\n".join(output) + "\n\n" + "\n".join(output_var) def _find_nn_pos_before_site(self, site_idx): @@ -3280,7 +3394,12 @@ def from_dict(cls, dct) -> IMolecule | Molecule: charge = dct.get("charge", 0) spin_multiplicity = dct.get("spin_multiplicity") properties = dct.get("properties") - return cls.from_sites(sites, charge=charge, spin_multiplicity=spin_multiplicity, properties=properties) + return cls.from_sites( + sites, + charge=charge, + spin_multiplicity=spin_multiplicity, + properties=properties, + ) def get_distance(self, i: int, j: int) -> float: """Get distance between site i and j. @@ -3308,7 +3427,16 @@ def get_sites_in_sphere(self, pt: ArrayLike, r: float) -> list[Neighbor]: for idx, site in enumerate(self._sites): dist = site.distance_from_point(pt) if dist <= r: - neighbors.append(Neighbor(site.species, site.coords, site.properties, dist, idx, label=site.label)) + neighbors.append( + Neighbor( + site.species, + site.coords, + site.properties, + dist, + idx, + label=site.label, + ) + ) return neighbors def get_neighbors(self, site: Site, r: float) -> list[Neighbor]: @@ -3535,7 +3663,9 @@ def to(self, filename: str = "", fmt: str = "") -> str | None: @classmethod def from_str( - cls, input_string: str, fmt: Literal["xyz", "gjf", "g03", "g09", "com", "inp", "json", "yaml"] + cls, + input_string: str, + fmt: Literal["xyz", "gjf", "g03", "g09", "com", "inp", "json", "yaml"], ) -> IMolecule | Molecule: """Reads the molecule from a string. @@ -3960,7 +4090,13 @@ def substitute(self, index: int, func_group: IMolecule | Molecule | str, bond_or # group. del self[index] for site in fgroup[1:]: - s_new = PeriodicSite(site.species, site.coords, self.lattice, coords_are_cartesian=True, label=site.label) + s_new = PeriodicSite( + site.species, + site.coords, + self.lattice, + coords_are_cartesian=True, + label=site.label, + ) self._sites.append(s_new) def remove_species(self, species: Sequence[SpeciesLike]) -> None: @@ -4080,7 +4216,11 @@ def sort(self, key: Callable | None = None, reverse: bool = False) -> None: self._sites.sort(key=key, reverse=reverse) def translate_sites( - self, indices: int | Sequence[int], vector: ArrayLike, frac_coords: bool = True, to_unit_cell: bool = True + self, + indices: int | Sequence[int], + vector: ArrayLike, + frac_coords: bool = True, + to_unit_cell: bool = True, ) -> None: """Translate specific sites by some vector, keeping the sites within the unit cell. @@ -4181,7 +4321,12 @@ def get_rand_vec(): for idx in range(len(self._sites)): self.translate_sites([idx], get_rand_vec(), frac_coords=False) - def make_supercell(self, scaling_matrix: ArrayLike, to_unit_cell: bool = True, in_place: bool = True) -> Structure: + def make_supercell( + self, + scaling_matrix: ArrayLike, + to_unit_cell: bool = True, + in_place: bool = True, + ) -> Structure: """Create a supercell. Args: @@ -4353,33 +4498,54 @@ def from_prototype(cls, prototype: str, species: Sequence, **kwargs) -> Structur return Structure.from_spacegroup("Im-3m", Lattice.cubic(kwargs["a"]), species, [[0, 0, 0]]) if prototype == "hcp": return Structure.from_spacegroup( - "P6_3/mmc", Lattice.hexagonal(kwargs["a"], kwargs["c"]), species, [[1 / 3, 2 / 3, 1 / 4]] + "P6_3/mmc", + Lattice.hexagonal(kwargs["a"], kwargs["c"]), + species, + [[1 / 3, 2 / 3, 1 / 4]], ) if prototype == "diamond": return Structure.from_spacegroup("Fd-3m", Lattice.cubic(kwargs["a"]), species, [[0, 0, 0]]) if prototype == "rocksalt": return Structure.from_spacegroup( - "Fm-3m", Lattice.cubic(kwargs["a"]), species, [[0, 0, 0], [0.5, 0.5, 0]] + "Fm-3m", + Lattice.cubic(kwargs["a"]), + species, + [[0, 0, 0], [0.5, 0.5, 0]], ) if prototype == "perovskite": return Structure.from_spacegroup( - "Pm-3m", Lattice.cubic(kwargs["a"]), species, [[0, 0, 0], [0.5, 0.5, 0.5], [0.5, 0.5, 0]] + "Pm-3m", + Lattice.cubic(kwargs["a"]), + species, + [[0, 0, 0], [0.5, 0.5, 0.5], [0.5, 0.5, 0]], ) if prototype in ("cscl"): return Structure.from_spacegroup( - "Pm-3m", Lattice.cubic(kwargs["a"]), species, [[0, 0, 0], [0.5, 0.5, 0.5]] + "Pm-3m", + Lattice.cubic(kwargs["a"]), + species, + [[0, 0, 0], [0.5, 0.5, 0.5]], ) if prototype in ("fluorite", "caf2"): return Structure.from_spacegroup( - "Fm-3m", Lattice.cubic(kwargs["a"]), species, [[0, 0, 0], [1 / 4, 1 / 4, 1 / 4]] + "Fm-3m", + Lattice.cubic(kwargs["a"]), + species, + [[0, 0, 0], [1 / 4, 1 / 4, 1 / 4]], ) if prototype in ("antifluorite"): return Structure.from_spacegroup( - "Fm-3m", Lattice.cubic(kwargs["a"]), species, [[1 / 4, 1 / 4, 1 / 4], [0, 0, 0]] + "Fm-3m", + Lattice.cubic(kwargs["a"]), + species, + [[1 / 4, 1 / 4, 1 / 4], [0, 0, 0]], ) if prototype in ("zincblende"): return Structure.from_spacegroup( - "F-43m", Lattice.cubic(kwargs["a"]), species, [[0, 0, 0], [1 / 4, 1 / 4, 3 / 4]] + "F-43m", + Lattice.cubic(kwargs["a"]), + species, + [[0, 0, 0], [1 / 4, 1 / 4, 3 / 4]], ) except KeyError as exc: raise ValueError(f"Required parameter {exc} not specified as a kwargs!") from exc @@ -4448,7 +4614,9 @@ def __init__( self._sites: list[Site] = list(self._sites) # type: ignore def __setitem__( # type: ignore - self, idx: int | slice | Sequence[int] | SpeciesLike, site: SpeciesLike | Site | Sequence + self, + idx: int | slice | Sequence[int] | SpeciesLike, + site: SpeciesLike | Site | Sequence, ) -> None: """Modify a site in the molecule. @@ -4586,7 +4754,14 @@ def remove_species(self, species: Sequence[SpeciesLike]): for site in self: new_sp_occu = {sp: amt for sp, amt in site.species.items() if sp not in species} if len(new_sp_occu) > 0: - new_sites.append(Site(new_sp_occu, site.coords, properties=site.properties, label=site.label)) + new_sites.append( + Site( + new_sp_occu, + site.coords, + properties=site.properties, + label=site.label, + ) + ) self.sites = new_sites def remove_sites(self, indices: Sequence[int]): @@ -4612,7 +4787,12 @@ def translate_sites(self, indices: Sequence[int] | None = None, vector: ArrayLik vector = [0, 0, 0] for idx in indices: site = self[idx] - new_site = Site(site.species, site.coords + vector, properties=site.properties, label=site.label) + new_site = Site( + site.species, + site.coords + vector, + properties=site.properties, + label=site.label, + ) self[idx] = new_site def rotate_sites( @@ -4742,7 +4922,10 @@ def substitute(self, index: int, func_group: IMolecule | Molecule | str, bond_or functional_group = functional_group.copy() vec = functional_group[0].coords - functional_group[1].coords vec /= np.linalg.norm(vec) - functional_group[0] = "X", functional_group[1].coords + float(bond_len) * vec + functional_group[0] = ( + "X", + functional_group[1].coords + float(bond_len) * vec, + ) # Align X to the origin. x = functional_group[0] diff --git a/pymatgen/transformations/site_transformations.py b/pymatgen/transformations/site_transformations.py index 71d294a4d7f..9f08a274905 100644 --- a/pymatgen/transformations/site_transformations.py +++ b/pymatgen/transformations/site_transformations.py @@ -18,6 +18,7 @@ from pymatgen.analysis.ewald import EwaldMinimizer, EwaldSummation from pymatgen.analysis.local_env import MinimumDistanceNN from pymatgen.symmetry.analyzer import SpacegroupAnalyzer +from pymatgen.transformations.advanced_transformations import EnumerateStructureTransformation from pymatgen.transformations.transformation_abc import AbstractTransformation if TYPE_CHECKING: @@ -428,9 +429,8 @@ def _enumerate_ordering(self, structure: Structure): for ind in indices: new_sp = {sp: occu * fraction for sp, occu in structure[ind].species.items()} struct[ind] = new_sp - # Perform enumeration - from pymatgen.transformations.advanced_transformations import EnumerateStructureTransformation + # Perform enumeration trans = EnumerateStructureTransformation() return trans.apply_transformation(struct, 10000) diff --git a/tests/command_line/test_chargemol_caller.py b/tests/command_line/test_chargemol_caller.py index 5f7fb8e5557..3c59f0ccdfa 100644 --- a/tests/command_line/test_chargemol_caller.py +++ b/tests/command_line/test_chargemol_caller.py @@ -1,11 +1,19 @@ from __future__ import annotations +import os import unittest +import zipfile +from typing import TYPE_CHECKING + +import pytest from pymatgen.command_line.chargemol_caller import ChargemolAnalysis from pymatgen.core import Element from pymatgen.util.testing import TEST_FILES_DIR +if TYPE_CHECKING: + from pathlib import Path + class TestChargemolAnalysis(unittest.TestCase): def test_parse_chargemol(self): @@ -25,15 +33,16 @@ def test_parse_chargemol(self): assert ca.ddec_rcubed_moments == [14.496002, 88.169236] assert ca.ddec_rfourth_moments == [37.648248, 277.371929] assert ca.cm5_charges == [0.420172, -0.420172] - assert ca.summary["ddec"]["partial_charges"] == ca.ddec_charges - assert ca.summary["ddec"]["dipoles"] == ca.dipoles - assert ca.summary["ddec"]["bond_order_sums"] == ca.bond_order_sums - assert ca.summary["ddec"]["rsquared_moments"] == ca.ddec_rsquared_moments - assert ca.summary["ddec"]["rcubed_moments"] == ca.ddec_rcubed_moments - assert ca.summary["ddec"]["rfourth_moments"] == ca.ddec_rfourth_moments + summary = ca.summary["ddec"] + assert summary["partial_charges"] == ca.ddec_charges + assert summary["dipoles"] == ca.dipoles + assert summary["bond_order_sums"] == ca.bond_order_sums + assert summary["rsquared_moments"] == ca.ddec_rsquared_moments + assert summary["rcubed_moments"] == ca.ddec_rcubed_moments + assert summary["rfourth_moments"] == ca.ddec_rfourth_moments assert ca.summary["cm5"]["partial_charges"] == ca.cm5_charges - assert ca.summary["ddec"]["bond_order_dict"] == ca.bond_order_dict - assert ca.summary["ddec"].get("spin_moments") is None + assert summary["bond_order_dict"] == ca.bond_order_dict + assert summary.get("spin_moments") is None assert ca.natoms == [1, 1] assert ca.structure is not None assert len(ca.bond_order_dict) == 2 @@ -55,3 +64,82 @@ def test_parse_chargemol2(self): assert ca.summary["ddec"]["spin_moments"] == ca.ddec_spin_moments assert ca.natoms is None assert ca.structure is None + + +@pytest.fixture() +def mock_xyz(tmp_path: Path): + test_dir = tmp_path / "chargemol" / "spin_unpolarized" + test_dir.mkdir(parents=True) + xyz_file_content = """2 + +C 0.00000 0.00000 0.00000 0.8432 +O 1.20000 0.00000 0.00000 -0.8432""" + + (test_dir / "DDEC6_even_tempered_net_atomic_charges.xyz").write_text(xyz_file_content) + return test_dir + + +def test_parse_chargemol(monkeypatch, mock_xyz): + monkeypatch.setattr( + ChargemolAnalysis, + "_download_and_unzip_atomic_densities", + lambda self, version, verbose: fake_download(self), + ) + + ca = ChargemolAnalysis(path=str(mock_xyz), run_chargemol=False) + ca._atomic_densities_path = os.path.expanduser("~/.cache/pymatgen/ddec/fake_file") + + assert ca.ddec_charges == [0.8432, -0.8432] + + +def fake_download(self, version: str = "latest", verbose: bool = True) -> None: + extraction_path = "~/.cache/pymatgen/ddec" + os.makedirs(os.path.expanduser(extraction_path), exist_ok=True) + + # Create a fake zipfile with a folder named "atomic_densities" + fake_zip_path = os.path.expanduser(f"{extraction_path}/fake_file.zip") + fake_folder_path = os.path.expanduser(f"{extraction_path}/atomic_densities") + + # Create a fake file inside the "atomic_densities" folder + fake_file_path = os.path.join(fake_folder_path, "fake_file.txt") + os.makedirs(fake_folder_path, exist_ok=True) + with open(fake_file_path, "w") as f: + f.write("fake content") + + # Create a zipfile containing the "atomic_densities" folder + with zipfile.ZipFile(fake_zip_path, "w") as fake_zip: + fake_zip.write(fake_file_path, arcname="atomic_densities/fake_file.txt") + + self._atomic_densities_path = os.path.expanduser(extraction_path) + + +def test_fake_download_and_modify_path(monkeypatch): + # mock os.path.exists to simulate atomic densities are not found + # monkeypatch.setattr(os.path, "exists", lambda *args, **kwargs: False) + + # mock subprocess.Popen to simulate a successful Chargemol execution + monkeypatch.setattr("subprocess.Popen", lambda *args, **kwargs: type("", (), {"returncode": 0})()) + + # monkeypatch urllib.request.urlretrieve to fake download + import urllib.request + + monkeypatch.setattr(urllib.request, "urlretrieve", lambda url, path: None) + + # monkeypatch the _download_and_unzip_atomic_densities method + monkeypatch.setattr( + ChargemolAnalysis, + "_download_and_unzip_atomic_densities", + lambda self: fake_download(self), + ) + + # Create an instance of ChargemolAnalysis + test_dir = f"{TEST_FILES_DIR}/chargemol/spin_unpolarized" + + ca = ChargemolAnalysis( + path=test_dir, + atomic_densities_path="fake_path", # force _download_and_unzip_atomic_densities to be called + run_chargemol=False, + ) + + # Your assertions + assert ca._atomic_densities_path == os.path.expanduser("~/.cache/pymatgen/ddec")