From 0d7e6cd3b7800664e850938a7ebfe1240a3bd319 Mon Sep 17 00:00:00 2001 From: LiuHanyu <41718895+lhycms@users.noreply.github.com> Date: Thu, 18 Jan 2024 19:21:41 +0800 Subject: [PATCH] Add `pymatgen.io.pwmat` module (#3512) * Add pymatgen.io.pwmat module * rm unitest inside pymatgen.io.pwmat dir * remove try,except expression: remove format expression: remove format expression: codespell * Add pwmat fmt within Structure.from_file() and Structure.from_str() * Remove abtractclass LineLocaterBase * Add output file of PWmat: Report, OUT.FERMI * Add output file of PWmat: DOS.totalspin, DOS.spinup, DOS.spindown * Rename part to part_upper to impove code quality to pass automatically test * add OutFermi, Report, Dosspin object to __init__.py in pwmat * Add input file named gen.kpt * Add high_symmetry_point * Add test for all new features * snake case variables, fix doc strings and return types * Revise comments and compress the test file. to pass test * Revise annotation again * Further modify the annotations * Further modify the annotations * Fix mispelling according to codespell * Modify annotation according to mypy * Modify annotation of filename: str -> PathLike * Add pwmat in FileFormats * Change the type annotation of filename to PathLike * Change type annotation: np.array -> np.ndarray * Utilize mypy check locally to pass 21 tests for pymatgen.io.pwmat * fix typos * del tests/io/pwmat/__init__.py * more informative unrecognized file extension error in Structure.from_file() * fnmatch compressed *.config* or *.pwmat* files * compress test files * snake_case dict keys and method names * don't write tmp test files to git repo * re.escape err msg * drop filename from expected err msg * Add test for Structure pwmat IO format * Add test for Structure pwmat IO format * Try to pass test in win * Decompress the test file and reduce the number of lines. * Add test for Structure.from_file() and Structure.to() in pwmat format * improve test_inputs.py assertions --------- Signed-off-by: LiuHanyu <41718895+lhycms@users.noreply.github.com> Co-authored-by: Janosh Riebesell --- pymatgen/core/structure.py | 18 +- pymatgen/io/pwmat/__init__.py | 5 + pymatgen/io/pwmat/inputs.py | 679 +++++++++++++++++++++++++ pymatgen/io/pwmat/outputs.py | 377 ++++++++++++++ tests/core/test_structure.py | 13 +- tests/files/pwmat/DOS.spinup_projected | 21 + tests/files/pwmat/MOVEMENT.lzma | Bin 0 -> 6367 bytes tests/files/pwmat/OUT.FERMI.lzma | Bin 0 -> 63 bytes tests/files/pwmat/REPORT | 179 +++++++ tests/files/pwmat/atom.config | 14 + tests/io/pwmat/test_inputs.py | 95 ++++ tests/io/pwmat/test_outputs.py | 74 +++ 12 files changed, 1470 insertions(+), 5 deletions(-) create mode 100644 pymatgen/io/pwmat/__init__.py create mode 100644 pymatgen/io/pwmat/inputs.py create mode 100644 pymatgen/io/pwmat/outputs.py create mode 100644 tests/files/pwmat/DOS.spinup_projected create mode 100644 tests/files/pwmat/MOVEMENT.lzma create mode 100644 tests/files/pwmat/OUT.FERMI.lzma create mode 100644 tests/files/pwmat/REPORT create mode 100644 tests/files/pwmat/atom.config create mode 100644 tests/io/pwmat/test_inputs.py create mode 100644 tests/io/pwmat/test_outputs.py diff --git a/pymatgen/core/structure.py b/pymatgen/core/structure.py index 6a772913201..1368f60f888 100644 --- a/pymatgen/core/structure.py +++ b/pymatgen/core/structure.py @@ -60,7 +60,7 @@ from pymatgen.util.typing import CompositionLike, SpeciesLike -FileFormats = Literal["cif", "poscar", "cssr", "json", "yaml", "yml", "xsf", "mcsqs", "res", ""] +FileFormats = Literal["cif", "poscar", "cssr", "json", "yaml", "yml", "xsf", "mcsqs", "res", "pwmat", ""] class Neighbor(Site): @@ -2671,7 +2671,7 @@ def to(self, filename: str | Path = "", fmt: FileFormats = "", **kwargs) -> str: fmt (str): Format to output to. Defaults to JSON unless filename is provided. If fmt is specifies, it overrides whatever the filename is. Options include "cif", "poscar", "cssr", "json", - "xsf", "mcsqs", "prismatic", "yaml", "yml", "fleur-inpgen". + "xsf", "mcsqs", "prismatic", "yaml", "yml", "fleur-inpgen", "pwmat". Non-case sensitive. **kwargs: Kwargs passthru to relevant methods. E.g., This allows the passing of parameters like symprec to the @@ -2752,6 +2752,10 @@ def to(self, filename: str | Path = "", fmt: FileFormats = "", **kwargs) -> str: with zopen(filename, mode="wt", encoding="utf8") as file: file.write(res_str) return res_str + elif fmt == "pwmat" or fnmatch(filename.lower(), "*.pwmat") or fnmatch(filename.lower(), "*.config"): + from pymatgen.io.pwmat import AtomConfig + + writer = AtomConfig(self, **kwargs) else: if fmt == "": raise ValueError(f"Format not specified and could not infer from {filename=}") @@ -2832,6 +2836,10 @@ def from_str( # type: ignore[override] from pymatgen.io.res import ResIO struct = ResIO.structure_from_str(input_string, **kwargs) + elif fmt == "pwmat": + from pymatgen.io.pwmat import AtomConfig + + struct = AtomConfig.from_str(input_string, **kwargs).structure else: raise ValueError(f"Invalid {fmt=}, valid options are {get_args(FileFormats)}") @@ -2910,8 +2918,12 @@ def from_file( # type: ignore[override] from pymatgen.io.res import ResIO struct = ResIO.structure_from_file(filename, **kwargs) + elif fnmatch(fname.lower(), "*.config*") or fnmatch(fname.lower(), "*.pwmat*"): + from pymatgen.io.pwmat import AtomConfig + + struct = AtomConfig.from_file(filename, **kwargs).structure else: - raise ValueError("Unrecognized file extension!") + raise ValueError(f"Unrecognized extension in {filename=}") if sort: struct = struct.get_sorted_structure() if merge_tol: diff --git a/pymatgen/io/pwmat/__init__.py b/pymatgen/io/pwmat/__init__.py new file mode 100644 index 00000000000..344ed9f9875 --- /dev/null +++ b/pymatgen/io/pwmat/__init__.py @@ -0,0 +1,5 @@ +"""This package implements modules for input and output to and from PWmat.""" +from __future__ import annotations + +from .inputs import AtomConfig +from .outputs import DosSpin, Movement, OutFermi, Report diff --git a/pymatgen/io/pwmat/inputs.py b/pymatgen/io/pwmat/inputs.py new file mode 100644 index 00000000000..3d48add1bbe --- /dev/null +++ b/pymatgen/io/pwmat/inputs.py @@ -0,0 +1,679 @@ +from __future__ import annotations + +import linecache +from abc import ABC, abstractmethod +from collections import Counter +from typing import TYPE_CHECKING + +import numpy as np +from monty.io import zopen +from monty.json import MSONable + +from pymatgen.core import Lattice, Structure +from pymatgen.symmetry.kpath import KPathSeek + +if TYPE_CHECKING: + from pymatgen.util.typing import PathLike + +__author__ = "Hanyu Liu" +__email__ = "domainofbuaa@gmail.com" +__date__ = "2024-1-16" + + +class LineLocator(MSONable): + """Find the line indices (starts from 1) of a certain paragraph of text from the file.""" + + @staticmethod + def locate_all_lines(file_path: PathLike, content: str) -> list[int]: + """Locate the line in file where a certain paragraph of text is located (return all indices) + + Args: + file_path (PathLike): Absolute path to file. + content (str): Certain paragraph of text that needs to be located. + """ + row_idxs: list[int] = [] # starts from 1 to be compatible with linecache package + row_no: int = 0 + with zopen(file_path, mode="rt") as f: + for row_content in f: + row_no += 1 + if content.upper() in row_content.upper(): + row_idxs.append(row_no) + return row_idxs + + +class ListLocator(MSONable): + """Find the element indices (starts from 0) of a certain paragraph of text from the list.""" + + @staticmethod + def locate_all_lines(strs_lst: list[str], content: str) -> list[int]: + """Locate the elements in list where a certain paragraph of text is located (return all indices) + + Args: + strs_lst (list[str]): List of strings. + content (str): Certain paragraph of text that needs to be located. + """ + str_idxs: list[int] = [] # starts from 0 to be compatible with list + str_no: int = -1 + for tmp_str in strs_lst: + str_no += 1 + if content.upper() in tmp_str.upper(): + str_idxs.append(str_no) + return str_idxs + + +class ACExtractorBase(ABC): + """A parent class of ACExtractor and ACstrExtractor, ensuring that they are as consistent as possible""" + + @abstractmethod + def get_n_atoms(self) -> int: + """Get the number of atoms in structure defined by atom.config file.""" + + @abstractmethod + def get_lattice(self) -> np.ndarray: + """Get the lattice of structure defined by atom.config file""" + + @abstractmethod + def get_types(self) -> np.ndarray: + """Get atomic number of atoms in structure defined by atom.config file""" + + @abstractmethod + def get_coords(self) -> np.ndarray: + """Get fractional coordinates of atoms in structure defined by atom.config file""" + + @abstractmethod + def get_magmoms(self) -> np.ndarray: + """Get atomic magmoms of atoms in structure defined by atom.config file""" + + +class ACExtractor(ACExtractorBase): + """Extract information contained in atom.config : number of atoms, lattice, types, frac_coords, magmoms""" + + def __init__(self, file_path: PathLike) -> None: + """Initialization function + + Args + file_path (str): The absolute path of atom.config file. + """ + self.atom_config_path = file_path + self.n_atoms = self.get_n_atoms() + self.lattice = self.get_lattice() + self.types = self.get_types() + self.coords = self.get_coords() + self.magmoms = self.get_magmoms() + + def get_n_atoms(self) -> int: + """Return the number of atoms in the structure.""" + first_row = linecache.getline(str(self.atom_config_path), 1) + return int(first_row.split()[0]) + + def get_lattice(self) -> np.ndarray: + """Return the lattice of structure. + + Returns: + lattice: np.ndarray, shape = (9,) + """ + basis_vectors: list[float] = [] + content: str = "LATTICE" + idx_row: int = LineLocator.locate_all_lines(file_path=self.atom_config_path, content=content)[0] + for row_idx in [idx_row + 1, idx_row + 2, idx_row + 3]: + row_content: list[str] = linecache.getline(str(self.atom_config_path), row_idx).split()[:3] + for value in row_content: + basis_vectors.append(float(value)) + + return np.array(basis_vectors) + + def get_types(self) -> np.ndarray: + """Return the atomic number of atoms in structure. + + Returns: + np.ndarray: Atomic numbers in order corresponding to sites + """ + content = "POSITION" + idx_row = LineLocator.locate_all_lines(file_path=self.atom_config_path, content=content)[0] + with open(self.atom_config_path) as file: + atom_config_content = file.readlines() + atomic_numbers_content = atom_config_content[idx_row : idx_row + self.n_atoms] + atomic_numbers_lst = [int(row.split()[0]) for row in atomic_numbers_content] # convert str to int + return np.array(atomic_numbers_lst) + + def get_coords(self) -> np.ndarray: + """Return the fractional coordinates in structure. + + Returns: + np.ndarray: Fractional coordinates. + """ + coords_lst: list[np.ndarray] = [] + content: str = "POSITION" + idx_row: int = LineLocator.locate_all_lines(file_path=self.atom_config_path, content=content)[0] + with open(self.atom_config_path) as f: + atom_config_content = f.readlines() + """ + row_content: + '29 0.377262291145329 0.128590184800933 0.257759805813488 1 1 1' + """ + for row_content in atom_config_content[idx_row : idx_row + self.n_atoms]: + row_content_lst = row_content.split() + coord_tmp = [float(value) for value in row_content_lst[1:4]] # convert str to float. + coords_lst.append(np.array(coord_tmp)) + return np.array(coords_lst).reshape(-1) + + def get_magmoms(self) -> np.ndarray: + """Return the magenetic moments of atoms in structure. + + Returns: + np.ndarray: The magnetic moments of individual atoms. + """ + content: str = "MAGNETIC" + magnetic_moments_lst: list[float] = [] + try: # Error: not containing magmoms info. + idx_row = LineLocator.locate_all_lines(file_path=self.atom_config_path, content=content)[-1] + + with open(self.atom_config_path) as f: + atom_config_content = f.readlines() + + magnetic_moments_content = atom_config_content[idx_row : idx_row + self.n_atoms] + # MAGNETIC + # 3 0.0 # atomic_number magmom + # ... + magnetic_moments_lst = [ + float(tmp_magnetic_moment.split()[-1]) for tmp_magnetic_moment in magnetic_moments_content + ] + except Exception: + magnetic_moments_lst = [0 for _ in range(self.n_atoms)] + return np.array(magnetic_moments_lst) + + +class ACstrExtractor(ACExtractorBase): + """Extract information from atom.config file. You can get str by slicing the MOVEMENT.""" + + def __init__(self, atom_config_str: str): + """Initialization function. + + Args: + atom_config_str (str): A string describing the structure in atom.config file. + """ + self.atom_config_str = atom_config_str + self.strs_lst = self.atom_config_str.split("\n") + self.num_atoms = self.get_n_atoms() + + def get_n_atoms(self) -> int: + """Return the number of atoms in structure. + + Returns: + int: The number of atoms + """ + return int(self.strs_lst[0].split()[0].strip()) + + def get_lattice(self) -> np.ndarray: + """Return the lattice of structure. + + Returns: + np.ndarray: Lattice basis vectors of shape=(9,) + """ + basis_vectors_lst = [] + aim_content = "LATTICE" + aim_idx = ListLocator.locate_all_lines(strs_lst=self.strs_lst, content=aim_content)[0] + for idx_str in [aim_idx + 1, aim_idx + 2, aim_idx + 3]: + # ['0.8759519000E+01', '0.0000000000E+00', '0.0000000000E+00'] + str_lst = self.strs_lst[idx_str].split()[:3] + for tmp_str in str_lst: + basis_vectors_lst.append(float(tmp_str)) # convert str to float + return np.array(basis_vectors_lst) + + def get_types(self) -> np.ndarray: + """Return the atomic number of atoms in structure. + + Returns: + np.ndarray: Types of elements. + """ + aim_content = "POSITION" + aim_idx = ListLocator.locate_all_lines(strs_lst=self.strs_lst, content=aim_content)[0] + strs_lst = self.strs_lst[aim_idx + 1 : aim_idx + self.num_atoms + 1] + atomic_numbers_lst = [int(entry.split()[0]) for entry in strs_lst] + return np.array(atomic_numbers_lst) + + def get_coords(self) -> np.ndarray: + """Return the fractional coordinate of atoms in structure. + + Returns: + np.ndarray: Fractional coordinates of atoms of shape=(num_atoms*3,) + """ + coords_lst = [] + aim_content = "POSITION" + aim_idx = ListLocator.locate_all_lines(strs_lst=self.strs_lst, content=aim_content)[0] + for tmp_str in self.strs_lst[aim_idx + 1 : aim_idx + self.num_atoms + 1]: + # ['14', '0.751401861790384', '0.501653718883189', '0.938307102003243', '1', '1', '1'] + tmp_strs_lst = tmp_str.split() + tmp_coord = [float(value) for value in tmp_strs_lst[1:4]] # convert str to float + coords_lst.append(tmp_coord) + return np.array(coords_lst).reshape(-1) + + def get_magmoms(self) -> np.ndarray: + """Return the magnetic moments of atoms in structure. + + Returns: + np.ndarray: Atomic magnetic moments. + """ + magnetic_moments_lst: list[float] = [] + aim_content: str = "MAGNETIC" + aim_idxs: list[int] = ListLocator.locate_all_lines(strs_lst=self.strs_lst, content=aim_content) + if len(aim_idxs) == 0: + magnetic_moments_lst = [0.0 for _ in range(self.num_atoms)] + else: + aim_idx = aim_idxs[0] + magnetic_moments_content = self.strs_lst[aim_idx + 1 : aim_idx + self.num_atoms + 1] + magnetic_moments_lst = [ + float(tmp_magnetic_moment.split()[-1]) for tmp_magnetic_moment in magnetic_moments_content + ] + return np.array(magnetic_moments_lst) + + def get_e_tot(self) -> np.ndarray: + """Return the total energy of structure. + + Returns: + np.ndarray: The total energy of the material system. + """ + # strs_lst: + # [' 216 atoms', 'Iteration (fs) = 0.3000000000E+01', + # ' Etot', 'Ep', 'Ek (eV) = -0.2831881714E+05 -0.2836665392E+05 0.4783678177E+02', + # ' SCF = 7'] + strs_lst = self.strs_lst[0].split(",") + aim_index = ListLocator.locate_all_lines(strs_lst=strs_lst, content="EK (EV) =")[0] + # strs_lst[aim_index].split() : + # ['Ek', '(eV)', '=', '-0.2831881714E+05', '-0.2836665392E+05', '0.4783678177E+02'] + return np.array([float(strs_lst[aim_index].split()[3].strip())]) + + def get_atom_energies(self) -> np.ndarray | None: + """Return the energies of individual atoms in material system. + + Returns: + np.ndarray | None : The energies of individual atoms within the material system. + + Description: + When turn on `ENERGY DEPOSITION`, PWmat will output energy per atom. + """ + energies = [] + aim_content = "Atomic-Energy, ".upper() + aim_idxs = ListLocator.locate_all_lines(strs_lst=self.strs_lst, content=aim_content) + if len(aim_idxs) == 0: + return None + aim_idx = aim_idxs[0] + for tmp_str in self.strs_lst[aim_idx + 1 : aim_idx + self.num_atoms + 1]: + """ + Atomic-Energy, Etot(eV),E_nonloc(eV),Q_atom:dE(eV)= -0.1281163115E+06 + 14 0.6022241483E+03 0.2413350871E+02 0.3710442365E+01 + """ + energies.append(float(tmp_str.split()[1])) + return np.array(energies) + + def get_atom_forces(self) -> np.ndarray: + """Return the force on atoms in material system. + + Returns: + np.ndarray: Forces acting on individual atoms of shape=(num_atoms*3,) + """ + forces = [] + aim_content = "Force".upper() + aim_idx = ListLocator.locate_all_lines(strs_lst=self.strs_lst, content=aim_content)[0] + + for line in self.strs_lst[aim_idx + 1 : aim_idx + self.num_atoms + 1]: + # ['14', '0.089910342901203', '0.077164252174742', '0.254144099204679'] + forces.append([float(val) for val in line.split()[1:4]]) + return -np.array(forces).reshape(-1) + + def get_virial(self) -> np.ndarray | None: + """Return the virial tensor of material system. + + Returns: + np.ndarray | None: Virial tensor of shape=(9,) + """ + virial_tensor = [] + aim_content = "LATTICE" + aim_idx = ListLocator.locate_all_lines(strs_lst=self.strs_lst, content=aim_content)[0] + + for tmp_idx in [aim_idx + 1, aim_idx + 2, aim_idx + 3]: + # tmp_strs_lst = + # ['0.8759519000E+01', '0.0000000000E+00', '0.0000000000E+00', + # 'stress', '(eV):', '0.115558E+02', '0.488108E+01', '0.238778E+01'] + tmp_strs_lst = self.strs_lst[tmp_idx].split() + tmp_aim_row_lst = ListLocator.locate_all_lines(strs_lst=tmp_strs_lst, content="STRESS") + if len(tmp_aim_row_lst) == 0: + return None + + for tmp_idx in [aim_idx + 1, aim_idx + 2, aim_idx + 3]: + # tmp_str_lst = ['0.120972E+02', '0.483925E+01', '0.242063E+01'] + tmp_str_lst = self.strs_lst[tmp_idx].split()[-3:] + virial_tensor.append(float(tmp_str_lst[0])) + virial_tensor.append(float(tmp_str_lst[1])) + virial_tensor.append(float(tmp_str_lst[2])) + + return np.array(virial_tensor) + + +class AtomConfig(MSONable): + """Object for representing the data in a atom.config or final.config file.""" + + def __init__(self, structure: Structure, sort_structure: bool = False): + """Initialization function. + + Args: + structure (Structure): Structure object + sort_structure (bool, optional): Whether to sort the structure. Useful if species + are not grouped properly together. Defaults to False. + """ + self.structure: Structure = structure + if sort_structure: + self.structure = self.structure.get_sorted_structure() + elements_counter = dict(sorted(Counter(self.structure.species).items())) + true_names = [f"{tmp_key}{tmp_value}" for (tmp_key, tmp_value) in elements_counter.items()] + self.true_names = "".join(true_names) + + def __repr__(self): + return self.get_str() + + def __str__(self): + return self.get_str() + + @classmethod + def from_str(cls, data: str, mag: bool = False) -> AtomConfig: + """Reads a atom.config from a string. + + Args: + data (str): string containing atom.config data + mag (bool, optional): Whether to read magnetic moment information. + + Returns: + AtomConfig object + """ + ac_extractor = ACstrExtractor(atom_config_str=data) + properties: dict[str, float] = {} + structure = Structure( + lattice=ac_extractor.get_lattice(), + species=ac_extractor.get_types(), + coords=ac_extractor.get_coords().reshape(-1, 3), + coords_are_cartesian=False, + properties=properties, + ) + if mag: + magmoms = ac_extractor.get_magmoms() + for idx, tmp_site in enumerate(structure): + tmp_site.properties.update({"magmom": magmoms[idx]}) + + return cls(structure) + + @classmethod + def from_file(cls, filename: PathLike, mag: bool = False) -> AtomConfig: + """Returns a AtomConfig from a file + + Args: + filename (PathLike): File name containing AtomConfig data + mag (bool, optional): Whether to read magnetic moments. Defaults to True. + + Returns: + AtomConfig object. + """ + with zopen(filename, "rt") as file: + return cls.from_str(data=file.read(), mag=mag) + + @classmethod + def from_dict(cls, dct: dict) -> AtomConfig: + """Returns a AtomConfig object from a dictionary. + + Args: + dct: dict containing atom.config data + + Returns: + AtomConfig object. + """ + return cls(Structure.from_dict(dct["structure"])) + + def get_str(self) -> str: + """Return a string describing the structure in atom.config format. + + Returns: + str: String representation of atom.config + """ + # This corrects for VASP really annoying bug of crashing on lattices + # which have triple product < 0. We will just invert the lattice + # vectors. + latt = self.structure.lattice + if np.linalg.det(latt.matrix) < 0: + latt = Lattice(-latt.matrix) + + lines: list[str] = [] + lines.append(f"\t{self.structure.num_sites} atoms\n") + lines.append("Lattice vector\n") + for ii in range(3): + lines.append(f"{latt.matrix[ii][0]:>15f}{latt.matrix[ii][1]:>15f}{latt.matrix[ii][2]:>15f}\n") + lines.append("Position, move_x, move_y, move_z\n") + for ii in range(self.structure.num_sites): + lines.append(f"{int(self.structure.species[ii].Z):>4d}") + lines.append(f"{self.structure.frac_coords[ii][0]:>15f}") + lines.append(f"{self.structure.frac_coords[ii][1]:>15f}") + lines.append(f"{self.structure.frac_coords[ii][2]:>15f}") + lines.append(" 1 1 1\n") + if "magmom" in self.structure.sites[0].properties: + lines.append("MAGNETIC\n") + for _, tmp_site in enumerate(self.structure.sites): + lines.append(f"{int(tmp_site.specie.Z):>4d}{tmp_site.properties['magmom']:>15f}\n") + return "".join(lines) + + def write_file(self, filename: PathLike, **kwargs): + """Writes AtomConfig to a file.""" + with zopen(filename, "wt") as f: + f.write(self.get_str(**kwargs)) + + def as_dict(self): + """ + Returns: + dict + """ + return { + "@module": type(self).__module__, + "@class": type(self).__name__, + "structure": self.structure.as_dict(), + "true_names": self.true_names, + } + + +class GenKpt(MSONable): + """GenKpt object for reading and writing gen.kpt. This file just generate line-mode kpoints.""" + + def __init__( + self, + reciprocal_lattice: np.ndarray, + kpoints: dict[str, np.ndarray], + path: list[list[str]], + density: float = 0.01, + ): + """Initialization function. + + Args: + reciprocal_lattice (np.array): Reciprocal lattice with factor of 2*pi. + kpoints (dict[str, np.array]): Kpoints and their corresponding fractional coordinates. + kpath (list[list[str]]): All kpaths, with each list representing one kpath. + density (float): The density of kpoints mesh with factor of 2*pi. + """ + self.reciprocal_lattice: np.ndarray = reciprocal_lattice + self.kpath: dict = {} + self.kpath.update({"kpoints": kpoints}) + self.kpath.update({"path": path}) + self.density = density + + @staticmethod + def from_structure(structure: Structure, dim: int, density: float = 0.01) -> GenKpt: + """Obtain a AtomConfig object from Structure object. + + Args: + strutcure (Structure): A structure object. + dim (int): The dimension of the material system (2 or 3). + density (float): Kpoints mesh without factor with 2*pi. Program will + automatically convert it with 2*pi. + """ + kpath_set = KPathSeek(structure) + if dim == 2: + kpts_2d: dict[str, np.ndarray] = {} + for tmp_name, tmp_kpt in kpath_set.kpath["kpoints"].items(): + if (tmp_kpt[2]) == 0: + kpts_2d.update({tmp_name: tmp_kpt}) + + path_2d: list[list[str]] = [] + for tmp_path in kpath_set.kpath["path"]: + tmp_path_2d: list[str] = [] + for tmp_hsp in tmp_path: + if tmp_hsp in kpts_2d: + tmp_path_2d.append(tmp_hsp) + if len(tmp_path_2d) > 1: + path_2d.append(tmp_path_2d) + kpts: dict[str, np.ndarray] = kpts_2d + path: list[list[str]] = path_2d + else: + kpts = kpath_set.kpath["kpoints"] + path = kpath_set.kpath["path"] + rec_lattice: np.ndarray = structure.lattice.reciprocal_lattice.matrix # with 2*pi + return GenKpt(rec_lattice, kpts, path, density * 2 * np.pi) + + def get_str(self): + """Returns a string to be written as a gen.kpt file.""" + + def calc_distance(hsp1: str, hsp2: str) -> float: + """Calculate the distance between two high symmetry points. + + Args: + hsp1 (str): The name of the first high symmetry point. + hsp2 (str): The name of the second high symmetry point. + + Returns: + distance (float): Distance between two high symmetry points.With factor of 2*pi + """ + hsp1_coord: np.ndarray = np.dot( + np.array(self.kpath["kpoints"][hsp1]).reshape(1, 3), self.reciprocal_lattice + ) + hsp2_coord: np.ndarray = np.dot( + np.array(self.kpath["kpoints"][hsp2]).reshape(1, 3), self.reciprocal_lattice + ) + return float(np.linalg.norm(hsp2_coord - hsp1_coord)) + + discontinue_pairs: list[list[int]] = [] + for ii in range(len(self.kpath["path"]) - 1): + discontinue_pairs.append([self.kpath["path"][ii][-1], self.kpath["path"][ii + 1][0]]) + + flatten_paths: list[str] = [tmp_hsp for tmp_path in self.kpath["path"] for tmp_hsp in tmp_path] + genkpt_str: str = f"Generated by pymatgen. density={self.density/(2*np.pi)}, " + genkpt_str += "fractional coordinates in reciprocal lattice.\n" + for ii in range(len(flatten_paths) - 1): + if [flatten_paths[ii], flatten_paths[ii + 1]] not in discontinue_pairs: + genkpt_str += f"{np.ceil(calc_distance(flatten_paths[ii], flatten_paths[ii+1])/self.density)}\n" + genkpt_str += f" {self.kpath['kpoints'][flatten_paths[ii]][0]:>12.6f}\t" + genkpt_str += f"{self.kpath['kpoints'][flatten_paths[ii]][1]:>12.6f}\t" + genkpt_str += f"{self.kpath['kpoints'][flatten_paths[ii]][2]:>12.6f}\t" + genkpt_str += f"{flatten_paths[ii]}\n" + genkpt_str += f" {self.kpath['kpoints'][flatten_paths[ii+1]][0]:>12.6f}\t" + genkpt_str += f"{self.kpath['kpoints'][flatten_paths[ii+1]][1]:>12.6f}\t" + genkpt_str += f"{self.kpath['kpoints'][flatten_paths[ii+1]][2]:>12.6f}\t" + genkpt_str += f"{flatten_paths[ii+1]}\n" + return genkpt_str + + def write_file(self, filename: PathLike): + """Writes gen.kpt to a file. + + Args: + filename (PathLike): The absolute path of file to be written. + """ + with zopen(filename, "wt") as file: + file.write(self.get_str()) + + +class HighSymmetryPoint(MSONable): + """HighSymmetryPoint object for reading and writing HIGH_SYMMETRY_POINTS file which generate line-mode kpoints.""" + + def __init__(self, reciprocal_lattice: np.ndarray, kpts: dict[str, list], path: list[list[str]], density: float): + """Initialization function. + + Args: + reciprocal_lattice (np.array): Reciprocal lattice. + kpts (dict[str, list[float]]): Kpoints and their corresponding fractional coordinates. + path (list[list[str]]): All k-paths, with each list representing one k-path. + density (float): Density of kpoints mesh with factor of 2*pi. + """ + self.reciprocal_lattice: np.ndarray = reciprocal_lattice + self.kpath: dict = {} + self.kpath.update({"kpoints": kpts}) + self.kpath.update({"path": path}) + self.density = density + + @staticmethod + def from_structure(structure: Structure, dim: int, density: float = 0.01) -> HighSymmetryPoint: + """Obtain HighSymmetry object from Structure object. + + Args: + structure (Structure): A structure object. + dim (int): Dimension of the material system (2 or 3). + density (float, optional): Density of kpoints mesh without factor of 2*pi. Defaults to 0.01. + The program will automatically convert it to with factor of 2*pi. + """ + reciprocal_lattice: np.ndarray = structure.lattice.reciprocal_lattice.matrix + gen_kpt = GenKpt.from_structure(structure=structure, dim=dim, density=density) + return HighSymmetryPoint( + reciprocal_lattice, gen_kpt.kpath["kpoints"], gen_kpt.kpath["path"], density * 2 * np.pi + ) + + def get_str(self) -> str: + """Returns a string describing high symmetry points in HIGH_SYMMETRY_POINTS format.""" + + def calc_distance(hsp1: str, hsp2: str) -> float: + """Calculate the distance of two high symmetry points. + + Returns: + distance (float): Calculate the distance of two high symmetry points. With factor of 2*pi. + """ + hsp1_coord: np.ndarray = np.dot( + np.array(self.kpath["kpoints"][hsp1]).reshape(1, 3), self.reciprocal_lattice + ) + hsp2_coord: np.ndarray = np.dot( + np.array(self.kpath["kpoints"][hsp2]).reshape(1, 3), self.reciprocal_lattice + ) + return float(np.linalg.norm(hsp2_coord - hsp1_coord)) + + def get_hsp_row_str(label: str, index: int, coordinate: float) -> str: + """ + Return string containing name, index, coordinate of the certain high symmetry point + in HIGH_SYMMETRY_POINTS format. + + Args: + label (str): Name of the high symmetry point. + index (int): Index of the high symmetry point. + coordinate (float): Coordinate in bandstructure of the high symmetry point. + + Returns: + str: String containing name, index, coordinate of the certain high symmetry point + in HIGH_SYMMETRY_POINTS format. + """ + if label == "GAMMA": + return f"G {index:>4d} {coordinate:>.6f}\n" + return f"{label} {index:>4d} {coordinate:>.6f}\n" + + discontinue_pairs: list[list[str]] = [] + for ii in range(len(self.kpath["path"]) - 1): + discontinue_pairs.append([self.kpath["path"][ii][-1], self.kpath["path"][ii + 1][0]]) + flatten_paths: list[str] = [tmp_hsp for tmp_path in self.kpath["path"] for tmp_hsp in tmp_path] + # flatten_paths = [hsp.replace("GAMMA", "G") for hsp in flatten_paths] + + index: int = 1 + coordinate: float = 0.0 + hsp_str: str = "Label Index Coordinate\n" + hsp_str += get_hsp_row_str(flatten_paths[0], index, coordinate) + for ii in range(1, len(flatten_paths)): + if [flatten_paths[ii - 1], flatten_paths[ii]] not in discontinue_pairs: + coordinate += calc_distance(flatten_paths[ii - 1], flatten_paths[ii]) / (2 * np.pi) + index += int(np.ceil(calc_distance(flatten_paths[ii - 1], flatten_paths[ii]) / self.density + 1)) + else: + coordinate += 0.0 + index += 1 + hsp_str += get_hsp_row_str(flatten_paths[ii], index, coordinate) + return hsp_str + + def write_file(self, filename: PathLike): + """Write HighSymmetryPoint to a file.""" + with zopen(filename, "wt") as file: + file.write(self.get_str()) diff --git a/pymatgen/io/pwmat/outputs.py b/pymatgen/io/pwmat/outputs.py new file mode 100644 index 00000000000..59897878dca --- /dev/null +++ b/pymatgen/io/pwmat/outputs.py @@ -0,0 +1,377 @@ +from __future__ import annotations + +import copy +import linecache +from io import StringIO +from typing import TYPE_CHECKING + +import numpy as np +from monty.io import zopen +from monty.json import MSONable + +from pymatgen.io.pwmat.inputs import ACstrExtractor, AtomConfig, LineLocator + +if TYPE_CHECKING: + from pymatgen.core import Structure + from pymatgen.util.typing import PathLike + + +__author__ = "Hanyu Liu" +__email__ = "domainofbuaa@gmail.com" +__date__ = "2024-1-16" + + +class Movement(MSONable): + """Parser for data in MOVEMENT which records trajectory during MD.""" + + def __init__(self, filename: PathLike, ionic_step_skip: int | None = None, ionic_step_offset: int | None = None): + """Initialization function. + + Args: + filename (PathLike): The path of MOVEMENT + ionic_step_skip (int | None, optional): If ionic_step_skip is a number > 1, + only every ionic_step_skip ionic steps will be read for + structure and energies. This is very useful if you are parsing + very large MOVEMENT files. Defaults to None. + ionic_step_offset (int | None, optional): Used together with ionic_step_skip. + If set, the first ionic step read will be offset by the amount of + ionic_step_offset. Defaults to None. + """ + self.filename: PathLike = filename + self.ionic_step_skip: int | None = ionic_step_skip + self.ionic_step_offset: int | None = ionic_step_offset + self.split_mark: str = "--------------------------------------" + + self.chunk_sizes, self.chunk_starts = self._get_chunk_info() + self.n_ionic_steps: int = len(self.chunk_sizes) + + self.ionic_steps: list[dict] = self._parse_sefv() + if self.ionic_step_offset and self.ionic_step_skip: + self.ionic_steps = self.ionic_steps[self.ionic_step_offset :: self.ionic_step_skip] + + def _get_chunk_info(self) -> tuple[list[int], list[int]]: + """Split MOVEMENT into many chunks, so that program process it chunk by chunk. + + Returns: + tuple[list[int], list[int]]: + chunk_sizes (list[int]): The number of lines occupied by structural + information in each step. + chunk_starts (list[int]): The starting line number for structural + information in each step. + """ + chunk_sizes: list[int] = [] + row_idxs: list[int] = LineLocator.locate_all_lines(self.filename, self.split_mark) + chunk_sizes.append(row_idxs[0]) + for ii in range(1, len(row_idxs)): + chunk_sizes.append(row_idxs[ii] - row_idxs[ii - 1]) + chunk_sizes_bak: list[int] = copy.deepcopy(chunk_sizes) + chunk_sizes_bak.insert(0, 0) + chunk_starts: list[int] = np.cumsum(chunk_sizes_bak).tolist() + chunk_starts.pop(-1) + return chunk_sizes, chunk_starts + + @property + def atom_configs(self) -> list[Structure]: + """Returns AtomConfig object for structures contained in MOVEMENT. + + Returns: + list[Structure]: List of Structure objects for the structure at each ionic step. + """ + return [step["atom_config"] for _, step in enumerate(self.ionic_steps)] + + @property + def e_tots(self) -> np.ndarray: + """Returns total energies of each ionic step structures contained in MOVEMENT. + + Returns: + np.ndarray: Total energy of of each ionic step structure, + with shape of (n_ionic_steps,). + """ + return np.array([step["e_tot"] for _, step in enumerate(self.ionic_steps)]) + + @property + def atom_forces(self) -> np.ndarray: + """Returns forces on atoms in each structures contained in MOVEMENT. + + Returns: + np.ndarray: The forces on atoms of each ionic step structure, + with shape of (n_ionic_steps, n_atoms, 3). + """ + return np.array([step["atom_forces"] for _, step in enumerate(self.ionic_steps)]) + + @property + def e_atoms(self) -> np.ndarray: + """ + Returns individual energies of atoms in each ionic step structures + contained in MOVEMENT. + + Returns: + np.ndarray: The individual energy of atoms in each ionic step structure, + with shape of (n_ionic_steps, n_atoms). + """ + return np.array([step["eatoms"] for _, step in enumerate(self.ionic_steps) if ("eatoms" in step)]) + + @property + def virials(self) -> np.ndarray: + """Returns virial tensor of each ionic step structure contained in MOVEMENT. + + Returns: + np.ndarray: The virial tensor of each ionic step structure, + with shape of (n_ionic_steps, 3, 3) + """ + return np.array([step["virial"] for _, step in enumerate(self.ionic_steps) if ("virial" in step)]) + + def _parse_sefv(self) -> list[dict]: + """ + Parse the MOVEMENT file, return information ionic step structure containing + structures, energies, forces on atoms and virial tensor. + + Returns: + list[dict]: Structure containing structures, energies, forces on atoms + and virial tensor. The corresponding keys are 'atom_config', 'e_tot', + 'atom_forces' and 'virial'. + """ + ionic_steps: list[dict] = [] + with zopen(self.filename, "rt") as mvt: + tmp_step: dict = {} + for ii in range(self.n_ionic_steps): + tmp_chunk: str = "" + for _ in range(self.chunk_sizes[ii]): + tmp_chunk += mvt.readline() + tmp_step.update({"atom_config": AtomConfig.from_str(tmp_chunk)}) + tmp_step.update({"e_tot": ACstrExtractor(tmp_chunk).get_e_tot()[0]}) + tmp_step.update({"atom_forces": ACstrExtractor(tmp_chunk).get_atom_forces().reshape(-1, 3)}) + e_atoms: np.ndarray | None = ACstrExtractor(tmp_chunk).get_atom_forces() + if e_atoms is not None: + tmp_step.update({"atom_energies": ACstrExtractor(tmp_chunk).get_atom_energies()}) + else: + print(f"Ionic step #{ii} : Energy deposition is turn down.") + virial: np.ndarray | None = ACstrExtractor(tmp_chunk).get_virial() + if virial is not None: + tmp_step.update({"virial": virial.reshape(3, 3)}) + else: + print(f"Ionic step #{ii} : No virial information.") + ionic_steps.append(tmp_step) + return ionic_steps + + +class OutFermi(MSONable): + """Extract fermi energy (eV) from OUT.FERMI""" + + def __init__(self, filename: PathLike): + """Initialization function + + Args: + filename (PathLike): The absolute path of OUT.FERMI file. + """ + self.filename: PathLike = filename + with zopen(self.filename, "rt") as f: + self._e_fermi: float = np.round(float(f.readline().split()[-2].strip()), 3) + + @property + def e_fermi(self) -> float: + """Returns the fermi energy level. + + Returns: + float: Fermi energy level. + """ + return self._e_fermi + + +class Report(MSONable): + """Extract information of spin, kpoints, bands, eigenvalues from REPORT file.""" + + def __init__(self, filename: PathLike): + """Initialization function. + + Args: + filename (PathLike): The absolute path of REPORT file. + """ + self.filename: PathLike = filename + self._spin, self._num_kpts, self._num_bands = self._parse_band() + self._eigenvalues = self._parse_eigen() + self._kpts, self._kpts_weight, self._hsps = self._parse_kpt() + + def _parse_band(self) -> tuple[int, int, int]: + """ + Parse REPORT file to obtain spin switches, the number of kpoints + and the number of bands. + + Returns: + spin (int): Whether turn on spin or not + 1: turn down the spin + 2: turn on the spin. + num_kpts (int): The number of kpoints. + num_bands (int): The number of bands. + """ + content: str = "SPIN" + row_idx: int = LineLocator.locate_all_lines(file_path=self.filename, content=content)[0] + spin = int(linecache.getline(str(self.filename), row_idx).split()[-1].strip()) + + content = "NUM_KPT" + row_idx = LineLocator.locate_all_lines(file_path=self.filename, content=content)[0] + num_kpts = int(linecache.getline(str(self.filename), row_idx).split()[-1].strip()) + + content = "NUM_BAND" + row_idx = LineLocator.locate_all_lines(file_path=self.filename, content=content)[0] + num_bands = int(linecache.getline(str(self.filename), row_idx).split()[-1].strip()) + return spin, num_kpts, num_bands + + def _parse_eigen(self) -> np.ndarray: + """Parse REPORT file to obtain information about eigenvalues. + + Returns: + np.ndarray: Eigenvalues with shape of (1 or 2, n_kpoints, n_bands). + The first index represents spin, the second index represents + kpoints, the third index represents band. + """ + num_rows: int = int(np.ceil(self._num_bands / 5)) + content: str = "eigen energies, in eV" + rows_lst: list[int] = LineLocator.locate_all_lines(file_path=self.filename, content=content) + rows_array: np.ndarray = np.array(rows_lst).reshape(self._spin, -1) + eigenvalues: np.ndarray = np.zeros((self._spin, self._num_kpts, self._num_bands)) + + for ii in range(self._spin): + for jj in range(self._num_kpts): + tmp_eigenvalues_str = "" + for kk in range(num_rows): + tmp_eigenvalues_str += linecache.getline(str(self.filename), rows_array[ii][jj] + kk + 1) + tmp_eigenvalues_array = np.array([float(eigen_value) for eigen_value in tmp_eigenvalues_str.split()]) + for kk in range(self._num_bands): + eigenvalues[ii][jj][kk] = tmp_eigenvalues_array[kk] + return eigenvalues + + def _parse_kpt(self) -> tuple[np.ndarray, np.ndarray, dict[str, np.ndarray]]: + """Parse REPORT file to obtain information about kpoints. + + Returns: + 3-tuple containing: + kpts (np.ndarray): The fractional coordinates of kpoints. + kpts_weight (np.ndarray): The weight of kpoints. + hsps (dict[str, np.ndarray]): The name and coordinates of high symmetric points. + """ + num_rows: int = int(self._num_kpts) + content: str = "total number of K-point:" + row_idx: int = LineLocator.locate_all_lines(self.filename, content)[0] + kpts: np.array = np.zeros((self._num_kpts, 3)) + kpts_weight: np.array = np.zeros(self._num_kpts) + hsps: dict[str, np.array] = {} + for ii in range(num_rows): + # 0.00000 0.00000 0.00000 0.03704 G + tmp_row_lst: list[str] = linecache.getline(str(self.filename), row_idx + ii + 1).split() + for jj in range(3): + kpts[ii][jj] = float(tmp_row_lst[jj].strip()) + kpts_weight[ii] = float(tmp_row_lst[3].strip()) + + if len(tmp_row_lst) == 5: + hsps.update( + {tmp_row_lst[4]: np.array([float(tmp_row_lst[0]), float(tmp_row_lst[1]), float(tmp_row_lst[2])])} + ) + return kpts, kpts_weight, hsps + + @property + def spin(self) -> int: + """Return the spin switches. + + Returns: + int: Spin switches. 1 represents turn on spin, 2 represents turn down spin. + """ + return self._spin + + @property + def n_kpoints(self) -> int: + """Returns the number of k-points.""" + return self._num_kpts + + @property + def n_bands(self) -> int: + """Returns the number of bands.""" + return self._num_bands + + @property + def eigenvalues(self) -> np.ndarray: + """Returns the eigenvalues. + + Returns: + np.ndarray: The first index represents spin, the second index + represents kpoint, the third index represents band. + """ + return self._eigenvalues + + @property + def kpoints(self) -> np.ndarray: + """Returns the fractional coordinates of kpoints.""" + return self._kpts + + @property + def kpoints_weight(self) -> np.ndarray: + """Returns the weight of kpoints.""" + return self._kpts_weight + + @property + def hsps(self) -> dict[str, np.ndarray]: + """Return the high symmetry points. + + Returns: + dict[str, np.ndarray]: The label and fractional coordinate of + high symmetry points. Return empty dict when task is not + line-mode kpath. + """ + return self._hsps + + +class DosSpin(MSONable): + """Extract information of DOS from DOS_SPIN file: + - DOS.totalspin, DOS.totalspin_projected + - DOS.spinup, DOS.spinup_projected + - DOS.spindown, DOS.spindown_projected + """ + + def __init__(self, filename: PathLike): + self.filename: PathLike = filename + self._labels, self._dos = self._parse() + + def _parse(self): + """Parse the DOS_SPIN file to get name and values of partial DOS. + + Returns: + labels (list[str]): The label of DOS, e.g. Total, Cr-3S, ... + dos (np.array): Value of density of state. + """ + labels: list[str] = [] + labels = linecache.getline(str(self.filename), 1).split()[1:] + dos_str: str = "" + with zopen(self.filename, mode="rt") as file: + file.readline() + dos_str = file.read() + dos: np.array = np.loadtxt(StringIO(dos_str)) + return labels, dos + + @property + def labels(self) -> list[str]: + """Returns the name of the partial density of states""" + return self._labels + + @property + def dos(self) -> np.ndarray: + """Returns value of density of state.""" + return self._dos + + def get_partial_dos(self, part: str) -> np.ndarray: + """Get partial dos for give element or orbital. + + Args: + part (str): The name of partial dos. + e.g. 'Energy', 'Total', 'Cr-3S', 'Cr-3P', + 'Cr-4S', 'Cr-3D', 'I-4D', 'I-5S', 'I-5P', 'Cr-3S', 'Cr-3Pz', + 'Cr-3Px', 'Cr-3Py', 'Cr-4S', 'Cr-3Dz2','Cr-3Dxz', 'Cr-3Dyz', + 'Cr-3D(x^2-y^2)', 'Cr-3Dxy', 'I-4Dz2', 'I-4Dxz', 'I-4Dyz', + 'I-4D(x^2-y^2)', 'I-4Dxy', 'I-5S', 'I-5Pz', 'I-5Px', 'I-5Py' + + Returns: + partial_dos: np.array + """ + part_upper: str = part.upper() + labels_upper: list[str] = [tmp_label.upper() for tmp_label in self._labels] + idx_dos = labels_upper.index(part_upper) + return self._dos[:, idx_dos] diff --git a/tests/core/test_structure.py b/tests/core/test_structure.py index 778a380b35f..18e279ee995 100644 --- a/tests/core/test_structure.py +++ b/tests/core/test_structure.py @@ -818,7 +818,7 @@ def test_get_dist_matrix(self): assert_allclose(self.struct.distance_matrix, ans) def test_to_from_file_and_string(self): - for fmt in ("cif", "json", "poscar", "cssr"): + for fmt in ("cif", "json", "poscar", "cssr", "pwmat"): struct = self.struct.to(fmt=fmt) assert struct is not None ss = IStructure.from_str(struct, fmt=fmt) @@ -849,6 +849,10 @@ def test_to_from_file_and_string(self): os.replace(yaml_path, yml_path) assert Structure.from_file(yml_path) == self.struct + atom_config_path = f"{self.tmp_path}/atom-test.config" + self.struct.to(filename=atom_config_path) + assert Structure.from_file(atom_config_path) == self.struct + with pytest.raises(ValueError, match="Format not specified and could not infer from filename='whatever'"): self.struct.to(filename="whatever") with pytest.raises(ValueError, match="Invalid fmt='badformat'"): @@ -1284,7 +1288,7 @@ def test_to_from_abivars(self): def test_to_from_file_str(self): # to/from string - for fmt in ("cif", "json", "poscar", "cssr", "yaml", "yml", "xsf", "res"): + for fmt in ("cif", "json", "poscar", "cssr", "yaml", "yml", "xsf", "res", "pwmat"): struct = self.struct.to(fmt=fmt) assert struct is not None ss = Structure.from_str(struct, fmt=fmt) @@ -1301,6 +1305,11 @@ def test_to_from_file_str(self): assert os.path.isfile(f"json-struct{ext}") assert Structure.from_file(f"json-struct{ext}") == self.struct + # test Structure.from_file with unsupported file extension (using tmp JSON file with wrong ext) + Path(filename := f"{self.tmp_path}/bad.extension").write_text(self.struct.to(fmt="json")) + with pytest.raises(ValueError, match="Unrecognized extension in filename="): + self.struct.from_file(filename=filename) + def test_from_spacegroup(self): s1 = Structure.from_spacegroup("Fm-3m", Lattice.cubic(3), ["Li", "O"], [[0.25, 0.25, 0.25], [0, 0, 0]]) assert s1.formula == "Li8 O4" diff --git a/tests/files/pwmat/DOS.spinup_projected b/tests/files/pwmat/DOS.spinup_projected new file mode 100644 index 00000000000..1483168ad7b --- /dev/null +++ b/tests/files/pwmat/DOS.spinup_projected @@ -0,0 +1,21 @@ + # Energy Total Cr-3S Cr-3Pz Cr-3Px Cr-3Py Cr-4S Cr-3Dz2 Cr-3Dxz Cr-3Dyz Cr-3D(x^2-y^2) Cr-3Dxy I-4Dz2 I-4Dxz I-4Dyz I-4D(x^2-y^2) I-4Dxy I-5S I-5Pz I-5Px I-5Py + -0.13461E+01 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 + -0.13246E+01 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 + -0.13031E+01 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 + -0.12816E+01 0.34071E+00 0.35900E-09 0.10133E-07 0.25147E-07 0.23790E-07 0.12139E-03 0.20441E-02 0.84560E-01 0.91495E-01 0.37174E-01 0.44490E-01 0.14927E-06 0.10245E-06 0.12008E-06 0.18766E-06 0.18648E-06 0.20205E-02 0.45595E-01 0.18111E-01 0.15096E-01 + -0.12601E+01 0.31440E+01 0.33128E-08 0.93505E-07 0.23205E-06 0.21953E-06 0.11202E-02 0.18863E-01 0.78031E+00 0.84430E+00 0.34304E+00 0.41054E+00 0.13775E-05 0.94543E-06 0.11081E-05 0.17317E-05 0.17208E-05 0.18645E-01 0.42075E+00 0.16712E+00 0.13931E+00 + -0.12386E+01 0.12168E+02 0.12652E-07 0.36904E-06 0.88627E-06 0.83881E-06 0.42823E-02 0.74463E-01 0.30229E+01 0.32783E+01 0.13227E+01 0.15761E+01 0.53079E-05 0.36590E-05 0.43012E-05 0.67135E-05 0.66721E-05 0.71811E-01 0.16209E+01 0.64645E+00 0.54987E+00 + -0.12171E+01 0.20942E+02 0.20787E-07 0.68535E-06 0.14385E-05 0.13863E-05 0.71044E-02 0.14140E+00 0.52372E+01 0.57049E+01 0.22484E+01 0.26086E+01 0.89334E-05 0.62860E-05 0.75092E-05 0.11640E-04 0.11609E-04 0.12075E+00 0.27294E+01 0.11156E+01 0.10282E+01 + -0.11956E+01 0.21697E+02 0.19394E-07 0.91468E-06 0.12197E-05 0.13576E-05 0.70079E-02 0.20340E+00 0.55807E+01 0.60513E+01 0.22707E+01 0.23191E+01 0.84012E-05 0.64825E-05 0.81936E-05 0.12344E-04 0.12560E-04 0.11428E+00 0.25881E+01 0.12005E+01 0.13616E+01 + -0.11741E+01 0.22704E+02 0.23226E-07 0.14172E-05 0.12329E-05 0.18420E-05 0.90313E-02 0.30569E+00 0.60199E+01 0.61867E+01 0.24599E+01 0.20584E+01 0.73862E-05 0.70120E-05 0.90919E-05 0.13119E-04 0.13807E-04 0.10694E+00 0.23466E+01 0.14553E+01 0.17553E+01 + -0.11526E+01 0.23745E+02 0.29542E-07 0.25670E-05 0.19648E-05 0.26289E-05 0.10844E-01 0.35531E+00 0.58583E+01 0.62687E+01 0.26746E+01 0.25069E+01 0.69076E-05 0.87350E-05 0.94940E-05 0.13565E-04 0.13608E-04 0.11375E+00 0.21809E+01 0.18277E+01 0.19477E+01 + -0.11311E+01 0.19620E+02 0.21256E-07 0.33466E-05 0.22220E-05 0.21729E-05 0.64687E-02 0.26815E+00 0.40402E+01 0.55248E+01 0.20412E+01 0.26279E+01 0.56193E-05 0.89586E-05 0.75819E-05 0.11121E-04 0.95211E-05 0.10137E+00 0.16247E+01 0.16614E+01 0.17237E+01 + -0.11096E+01 0.10324E+02 0.79401E-08 0.21544E-05 0.12752E-05 0.94299E-06 0.16730E-02 0.11774E+00 0.18007E+01 0.32235E+01 0.92294E+00 0.15281E+01 0.29676E-05 0.51970E-05 0.39460E-05 0.56898E-05 0.42689E-05 0.54951E-01 0.77733E+00 0.88606E+00 0.10112E+01 + -0.10881E+01 0.44444E+01 0.19455E-08 0.86787E-06 0.40237E-06 0.22697E-06 0.18111E-03 0.25819E-01 0.70127E+00 0.17218E+01 0.24756E+00 0.51230E+00 0.10482E-05 0.18348E-05 0.19250E-05 0.18140E-05 0.15609E-05 0.19726E-01 0.22449E+00 0.36723E+00 0.62392E+00 + -0.10666E+01 0.36435E+01 0.12693E-08 0.56964E-06 0.18634E-06 0.10056E-06 0.14838E-03 0.33384E-02 0.57736E+00 0.16163E+01 0.11903E+00 0.26118E+00 0.62096E-06 0.10494E-05 0.17908E-05 0.93458E-06 0.11694E-05 0.12073E-01 0.95928E-01 0.30205E+00 0.65601E+00 + -0.10451E+01 0.33347E+01 0.43016E-08 0.33979E-06 0.21867E-06 0.38270E-06 0.14544E-02 0.82830E-02 0.65831E+00 0.11043E+01 0.32135E+00 0.22903E+00 0.40064E-06 0.93506E-06 0.15780E-05 0.92715E-06 0.11646E-05 0.90440E-02 0.94345E-01 0.40197E+00 0.50658E+00 + -0.10236E+01 0.58741E+01 0.16215E-07 0.19114E-06 0.58394E-06 0.14115E-05 0.60329E-02 0.33000E-01 0.15510E+01 0.86911E+00 0.11230E+01 0.41715E+00 0.30213E-06 0.17047E-05 0.24358E-05 0.19353E-05 0.21688E-05 0.10903E-01 0.19493E+00 0.10710E+01 0.59803E+00 + -0.10021E+01 0.10835E+02 0.31232E-07 0.40408E-06 0.11559E-05 0.24150E-05 0.11305E-01 0.57235E-01 0.30726E+01 0.15089E+01 0.19684E+01 0.70274E+00 0.38109E-06 0.29833E-05 0.40935E-05 0.32807E-05 0.34284E-05 0.15621E-01 0.29270E+00 0.20991E+01 0.11066E+01 + -0.98065E+00 0.10815E+02 0.29570E-07 0.61451E-06 0.11825E-05 0.19044E-05 0.10200E-01 0.45669E-01 0.31833E+01 0.17801E+01 0.16335E+01 0.61788E+00 0.30659E-06 0.27595E-05 0.37117E-05 0.27968E-05 0.26573E-05 0.11797E-01 0.19789E+00 0.21163E+01 0.12185E+01 + -0.95915E+00 0.56892E+01 0.13959E-07 0.39339E-06 0.65418E-06 0.78874E-06 0.45785E-02 0.17782E-01 0.16670E+01 0.10485E+01 0.73084E+00 0.33344E+00 0.13235E-06 0.14037E-05 0.18132E-05 0.12331E-05 0.10482E-05 0.46764E-02 0.66384E-01 0.11086E+01 0.70733E+00 + -0.93765E+00 0.26331E+01 0.33024E-08 0.10313E-06 0.38272E-06 0.38138E-06 0.99164E-03 0.34258E-02 0.57600E+00 0.44439E+00 0.39186E+00 0.31456E+00 0.56428E-07 0.81113E-06 0.89247E-06 0.43818E-06 0.37850E-06 0.22301E-02 0.37830E-01 0.47422E+00 0.38761E+00 \ No newline at end of file diff --git a/tests/files/pwmat/MOVEMENT.lzma b/tests/files/pwmat/MOVEMENT.lzma new file mode 100644 index 0000000000000000000000000000000000000000..c4151eb0ce0b8bb3ccf732ef311ec1a5eaa2fae2 GIT binary patch literal 6367 zcmV<57$E0e004jh|NsC0|NsC001#l^#g`+d~G<-$_mM zE%v6NH*{w@yPZ7nF?w5%jTRXU?geMDRJ{ARfM7>9wVho_wBV#B5z-Q_GV!TbE5D~q zv>WHMW5NHLwLn4N{>?IXVF=v*PZ=3}@S;5!LS8uF^7Nf|k+A;zhi3GYae3EN_hRhy zDkIlGIoFKeUS??{wr{LM^77B-_?8bVp#Dfn*#fz7Q{}nK1g3TP;g>)H9O}w8%JmF1 zUT;_gcbll>P$w^L#`&rc`4q|1H}@DZ^q~YIFUk+pg~?V;Vw7-h5NmLGHI^h~f1lJG zDL3F!Sm%1U)x~uKM}W6UyiaU%<;KtsdUQuCzf5ec`AB2l`btcf#$tKu$c@QI>e=XW4ZP{!(c_5D;2j=Hb<%)T9+Z zajGGNnLrCvWZi_?+%O=pv}!v^G1GB)4p%ISvj$Xp$*q(AB7Jsnr&#U=P9~W(pCp*z zZ6I>-hK*t8Pz%3S-pDN;bxt1`dAY`JI5MJn*yH_u4LfIy)hZ6@r|%Gp08nJTma#cw zJCe=q?&yFnKyeL}%ef6VZaNg=^pEX1Lc^9Vi6r@u37mNk6l*Ns``8#I zr-K6*Y`S7LIY@XBj-}f9z9p)X1NPi*nwlCQ?Xt+w=dlpQzeevj9~>TJYoir&j5-@K z4OBbqudQOZ(K0o~3*FP`^GYw{yo3K5TbD~rs-@Mbw&$H~_Z*~SGmqGynxzr;~bh-u{#|`{J->dl5XI zPpCNS=L82ODk4m@%BjNa8i?M<=gGh$lna7n)j|9Y)d1QAsYiIUKs|`qh-5T-#MneE zV|y6^rZh>rdBDr1v{EojI9~Pwwc5Z5J&Cgcrb;MD^|(UA7TxXahEhECHoXGW7@$Qh zyUKPv3_MUDq_&!@s5ddrbi=ZWu50Xw z9Kig9KDuIc0lBXd??xMwo~3hsNWataV-n&*{wFR>p9i6eaDUE096ML}*LjAOppso1 z47Y>frI9pU%ruV|2H0LRWqZfcZsQu@$rYz#iU;%tb*DR7kSxH1cFK+JC7b)*oSPP; z>&F1DgsskO2hz=E<0f<+@px50CGMc@R;S1yAFgT@CrB1t6O>=u zz^B$cLyZaK!Q+XO02G<(pJJ;Nm4l-YH8-K$l6a~<2W#;48SY$WIAUuKTs3lr@jPaJ ztKRT?rQ)#PPlO9uTk@N(#nO;E(Hmjn zQ#p&X?FLtQXIcA=7HuyFb4ZL6K+(CYN{xrJL=ZwZBrf-TZtgPQWMGEq!jRaJvBsTM z&28@hLC=mI4SKb`E#mPo`x^`$x3=FB6TQ6VO^*HH0=$vad(HUfG!X1QHdl19@s9so zfvv=Xk*%(;SZoZmj_^RuvxFrmq{^2orxl*=cJ7DSq3g98HEVGCzUCMG>nRW730Qk- zuGfn$1b;Alx+c>JN&q#SP=0D42UaH24>+ovoXt3iU^N9iS}@ps&V#B%1q9^Z^vGQi zWch+Cl5Vune+~R2t9Nn|y4?wwr9M`I^1t@lVu72{qA4U8?owuvMpi)kc>x)*yY>OUZwhh9-Tq1Fhy(#iVAxxyEHZ0BT` zd`Nh-!IV>V^QOz8{0A&7r;Q9Xm)n34GbX&%qzF-ZiupB$Nz)6%XN{ac3AjLjxdrkM zTYGLOC4b(e4 zFI907^qPJf>o+E~ySb1{8gI-RPsUi7jX9$F5ypx6bnq)5MirrCy>aMA*rq=nY`XXt z%ZAiJCyHS)iaUNskhwSnR4c~ki~#?mTf+|UkHN?j?M0)mI>*Md8@bu?A~0TPP*&L` zL3yUDdHXq~Z9Oah@IP=*h2$ZdJSwO@$W{gHH$)=P_)w>F(5PI2K;mz$=-Oa< zkZxWSZlZZ!Biqf{LoPCxCt2R>qmNwYhMShBBytDJ8Q>J1qh_89ZjS5itFERCs}uY? zm+$j93BsznsQ$72cNMW)$2gXLbO}YZk4axTx=~1sw9*evL+afHajKT*3|Gz8sFCw) zxM1n4a~S5@n!Z|i=S|`r1JEk>1Cm`ZhD4!^e8(Fp!FuyZfrge) zQeXxQ5{e;xH7kE+5KIz()WQOaC&VSnUI_lGzsD6g;o10FFlxbEhamY4J_q-^Xb?QT zw$~{|X;r;P>b|wqg({+Jw=+3RAWZ`L97#dLDPYh_CJ`;kUhq^eALuI3KDXB&NPtZA zWT!+@NZZ1tI=R-JE6xRGrdT|8gYA?kTXqWY*^7w?eN0y6)K&aE=%xi8`$5seSQMZ}q2)|iWwI)qw>TpS62FD4A=!}n#Z)>cg5b7tf2^BG{|Y=!msrG@g?Qq@1uA-;DGHEOxwQXu zA`#bm>>p55O$#V5gkQ>P!X8^&GFXP{bDOP3R)KBYSE8wS5;ktG~RxcnmXr~o=A2JkA+aAm~RpwmA*P7zJRs1Ocm;Up85C&}QS=;hy zf#%6as9PlpDZ=dKbuT{MsypG2QnLa zgnu))3M9Ru7?Eizb=y43%+E0sAix*VMZZ!&2hg_F(kwj=blb-zCM+o>HTv%+MF}CE z^tBWLfgt$n<~6cb&>GX%)FB=m7~B3EUE5-<)cVOk&HAIz>Ybffx^~U7bvP?ZkXq|O zsB00Wv?3~*A)ozhdS*exzM z+cQ~zZYKyTE!f*4P+n9bI1xWH9h_$-66*OWqmNxl{y?(>)3g`Z0}T zAtMwqq5N%Zu2 z(J2FgC4wQdVzg9~$oZYP^FB7)C~0MDg0s}M+o$k0A9Nrluow2u?c;MIT+i?&(@|df z4EIf+@i>P^qb`*v1%5P0d}5ZHoXg2Pkd<&zN*gRxZfY4=sGSMftM#I_^vv9IB^h%K z*Q8QL#Vb7)TQW?ehS&Jff6nlZa?v%=3(V*i#F<(6A}lH&?H!bbJ>+zc14IADj1ymW zOy0BpOGB8dhj?Q(DyDx8xh1XwV;ciUU5pLy9}VA_pdnh+PR-b`zyeqyLG(8MmeQSm z95^^_IjZe?=zaVI@x()oBjICLw>kJ93X|qKx2&w}ZgZpT5*8s0&u$Q1?pkIqXL&02BS5o=P8)`-<6#7;yMz||-YapVY~o4S|r;u|w! zkdfXJTIsAWbpQ-=)~VenS7H=V+waH6ma7a@=!=>2e5)L;%~Dq=TJ!WsE4;><<`Ix@ zFd&ID-@;u761h*;AM;n4+j)V<^FIs6+B@LEK_o1E3FERVOd<{7YblAwbB*Yj@s?%h zFlh3=p5!4N?h%}7m~Ps2vfz22XZlybUggXbT!q0IUzeyatC_9oEK%XaTjP9Jg*nrY zQygDvZD9AYV8bUY6gk`QW)>}mr_^F;9h+syvy(KL1#Co2`sSlk&3fZ8e`i&`?QOV1 zjWw??8c~|PbX-vUI9~H?$BRZFXpb4i7)>!zl5J94DA6#@fH%YY-WWLMF^dCPSuwUt zBqF;ZDEa;EnWH7x9;@=8# z{uM7>KF9LCfl^6(O5duVx<96j&WC0Qt4Ly+CLAB4$Cg}x4uKE9X?MCL7$ zZzNU_=`(9|ca9$gbYz7=sHB*nAbai{mZ7}AfC#%o^P_ocN9lv;(v z8Lt5bhOtjfQREDo^=biZztMUpR?l6~Ez@Z4OL6VW+!f`;Z4*1`zrH{%jH?-H>j?K- zRurFS-cM{ntg#56i=-L0c26+#(}By?)t8ya_s~+YaI3r^8N+LSNSMYUCt@4Ekc_wG zvj}if*N?~@W?As?`IT_3%_jkz=VjD-^hkx%@Eg9R@MJoR@Dts4gQ<19O zX*0f!9z5V3fT`5OXGOO+bG1L*2|LF{X^`+rEx#7Wx)RR#0$M}U>%IEo6&+01#Uj`N#%F!Uy{B?M8?pF?w930g>Wf^^=_<7TPituyoL{Idv;`O=0n zcbf<=ZtwW8pZ%E?*2g%&n&>ok!`w#kJ&fMdcdT8SmRj9l`s9#my{9bX_(|qs76A1k z)}t%M@N_yeFd+OwGf?2W3FFbh3kBR&QEo!SRb(+P{^w;`-)zb^EJrvZk6l+2+dd#d zS0$?>Ev*8&4>9mZXlN`N8_4)9RanbSb@&4T!6QzShy{XXS6t3DWy5s;u}|)#4tSIp zQER+v=r1PoyFy`3@dsStOQ?1N%9p3>5EIlqeHx`k_1pqia(<_pTSZ|?a^8Pz)je1W5))&u7K@q-`$_k>-!6pHojIM#&5x7E6wPVjX9E~s zZk+P>k5_|7Gqk$bNam$08NF|84UdgYG}vDKaG0nDkFJ14M}q+j^oE-&d}{lC2^QkNKFFtKF1sZ z33d~2A&YBbJKR&6H%9~