diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 90b2bbfd3a6..762fc45e001 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -56,3 +56,9 @@ repos: # MD041: first line in a file should be a top-level heading # MD025: single title args: [--disable, MD013, MD024, MD025, MD033, MD041, "--"] + + - repo: https://github.com/kynan/nbstripout + rev: 0.6.1 + hooks: + - id: nbstripout + args: [--drop-empty-cells, --keep-output] diff --git a/examples/aims_io/FHI-aims-example.ipynb b/examples/aims_io/FHI-aims-example.ipynb new file mode 100644 index 00000000000..af4d207dc30 --- /dev/null +++ b/examples/aims_io/FHI-aims-example.ipynb @@ -0,0 +1,161 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "8e09cccf-4335-4c59-bfa7-d3e3caeef404", + "metadata": {}, + "outputs": [], + "source": [ + "from pymatgen.io.aims.inputs import AimsGeometryIn, AimsCube, AimsControlIn\n", + "from pymatgen.io.aims.outputs import AimsOutput\n", + "\n", + "from pymatgen.core import Structure, Lattice\n", + "\n", + "import numpy as np\n", + "\n", + "from pathlib import Path\n", + "from subprocess import check_call\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "536af143-f892-4549-86d7-ff426ba265ed", + "metadata": {}, + "outputs": [], + "source": [ + "# AIMS_CMD should be modified to match your system\n", + "AIMS_CMD = \"aims.x\"\n", + "AIMS_OUTPUT = \"aims.out\"\n", + "AIMS_SD = \"species_dir\"\n", + "AIMS_TEST_DIR = \"../../tests/io/aims/species_directory/light/\"\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4d8d9fda-af37-45eb-971c-56fed59f3a27", + "metadata": {}, + "outputs": [], + "source": [ + "# Create test structure\n", + "structure = Structure(\n", + " lattice=Lattice(\n", + " np.array([[0, 2.715, 2.715],[2.715, 0, 2.715], [2.715, 2.715, 0]])\n", + " ),\n", + " species=[\"Si\", \"Si\"],\n", + " coords=np.array([np.zeros(3), np.ones(3) * 0.25])\n", + ")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8e67e134-84f8-4c35-afe4-87cd66e2e781", + "metadata": {}, + "outputs": [], + "source": [ + "# Create the geometry file from the structure\n", + "geo_in = AimsGeometryIn.from_structure(structure)\n", + "\n", + "# Create the control.in file\n", + "cont_in = AimsControlIn(\n", + " {\n", + " \"xc\": \"pw-lda\",\n", + " \"relax_geometry\": \"trm 0.01\",\n", + " \"relax_unit_cell\": \"full\",\n", + " \"species_dir\": AIMS_SD,\n", + " }\n", + ")\n", + "\n", + "# Add new parameters as if AimsControl\n", + "cont_in[\"k_grid\"] = [1, 1, 1]\n", + "\n", + "# Output options to control in automatically append the list\n", + "cont_in[\"output\"] = \"hirshfeld\"\n", + "cont_in[\"output\"] = [\"eigenvectors\"]\n", + "\n", + "# Cube file output controlled by the AimsCube class\n", + "cont_in[\"cubes\"] = [\n", + " AimsCube(\"total_density\", origin=[0,0,0], points=[11, 11, 11]),\n", + " AimsCube(\"eigenstate_density 1\", origin=[0, 0, 0], points=[11, 11, 11]),\n", + "]\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "096337d6-871a-48dc-b4b3-a3c7c6fd812e", + "metadata": {}, + "outputs": [], + "source": [ + "# Write the input files\n", + "workdir = Path.cwd() / \"workdir/\"\n", + "workdir.mkdir(exist_ok=True)\n", + "\n", + "geo_in.write_file(workdir, overwrite=True)\n", + "cont_in.write_file(structure, workdir, overwrite=True)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c9994cd2-5e45-4071-ab87-b6b3e3af1174", + "metadata": {}, + "outputs": [], + "source": [ + "# Run the calculation\n", + "with open(f\"{workdir}/{AIMS_OUTPUT}\", \"w\") as outfile:\n", + " aims_run = check_call([AIMS_CMD], cwd=workdir, stdout=outfile)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0c42a6c6-bb45-472f-a2cb-190cdd922047", + "metadata": {}, + "outputs": [], + "source": [ + "# Read the aims output file and the final relaxed geometry\n", + "outputs = AimsOutput.from_outfile(f\"{workdir}/{AIMS_OUTPUT}\")\n", + "relaxed_structure = AimsGeometryIn.from_file(f\"{workdir}/geometry.in.next_step\")\n", + "\n", + "# Check the results\n", + "assert outputs.get_results_for_image(-1).lattice == relaxed_structure.structure.lattice\n", + "assert np.all(outputs.get_results_for_image(-1).frac_coords == relaxed_structure.structure.frac_coords)\n", + "\n", + "assert np.allclose(\n", + " outputs.get_results_for_image(-1).properties[\"stress\"],\n", + " outputs.stress\n", + ")\n", + "\n", + "assert np.allclose(\n", + " outputs.get_results_for_image(-1).site_properties[\"force\"],\n", + " outputs.forces\n", + ")\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.16" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/pymatgen/io/aims/__init__.py b/pymatgen/io/aims/__init__.py new file mode 100644 index 00000000000..d0d04312bee --- /dev/null +++ b/pymatgen/io/aims/__init__.py @@ -0,0 +1 @@ +"""IO interface for FHI-aims.""" diff --git a/pymatgen/io/aims/inputs.py b/pymatgen/io/aims/inputs.py new file mode 100644 index 00000000000..209a3a7498a --- /dev/null +++ b/pymatgen/io/aims/inputs.py @@ -0,0 +1,619 @@ +"""Classes for reading/manipulating/writing FHI-aims input files. + +Works for aims cube objects, geometry.in and control.in +""" + +from __future__ import annotations + +import gzip +import time +from copy import deepcopy +from dataclasses import dataclass, field +from pathlib import Path +from typing import TYPE_CHECKING, Any + +import numpy as np +from monty.json import MontyDecoder, MSONable + +from pymatgen.core import Lattice, Molecule, Structure + +if TYPE_CHECKING: + from collections.abc import Sequence + + +@dataclass +class AimsGeometryIn(MSONable): + """Representation of an aims geometry.in file + + Attributes: + _content (str): The content of the input file + _structure (Structure or Molecule): The structure or molecule + representation of the file + """ + + _content: str + _structure: Structure | Molecule + + @classmethod + def from_str(cls, contents: str) -> AimsGeometryIn: + """Create an input from the content of an input file + + Args: + contents (str): The content of the file + + Returns: + The AimsGeometryIn file for the string contents + """ + content_lines = [ + line.strip() for line in contents.split("\n") if len(line.strip()) > 0 and line.strip()[0] != "#" + ] + + species = [] + coords = [] + is_frac = [] + lattice_vectors = [] + charges_dct = {} + moments_dct = {} + + for line in content_lines: + inp = line.split() + if (inp[0] == "atom") or (inp[0] == "atom_frac"): + coords.append([float(ii) for ii in line.split()[1:4]]) + species.append(inp[4]) + is_frac.append(inp[0] == "atom_frac") + if inp[0] == "lattice_vector": + lattice_vectors.append([float(ii) for ii in line.split()[1:4]]) + if inp[0] == "initial_moment": + moments_dct[len(coords) - 1] = float(inp[1]) + if inp[0] == "initial_charge": + charges_dct[len(coords) - 1] = float(inp[1]) + + charge = np.zeros(len(coords)) + for key, val in charges_dct.items(): + charge[key] = val + + magmom = np.zeros(len(coords)) + for key, val in moments_dct.items(): + magmom[key] = val + + if len(lattice_vectors) == 3: + lattice = Lattice(lattice_vectors) + for cc in range(len(coords)): + if is_frac[cc]: + coords[cc] = lattice.get_cartesian_coords(np.array(coords[cc]).reshape(1, 3)).flatten() + elif len(lattice_vectors) == 0: + lattice = None + if any(is_frac): + raise ValueError("Fractional coordinates given in file with no lattice vectors.") + else: + raise ValueError("Incorrect number of lattice vectors passed.") + + site_props = {"magmom": magmom, "charge": charge} + if lattice is None: + structure = Molecule(species, coords, np.sum(charge), site_properties=site_props) + else: + structure = Structure( + lattice, species, coords, np.sum(charge), coords_are_cartesian=True, site_properties=site_props + ) + + return cls(_content="\n".join(content_lines), _structure=structure) + + @classmethod + def from_file(cls, filepath: str | Path) -> AimsGeometryIn: + """Create an AimsGeometryIn from an input file. + + Args: + filepath (str | Path): The path to the input file (either plain text of gzipped) + + Returns: + AimsGeometryIn: The input object represented in the file + """ + if str(filepath).endswith(".gz"): + with gzip.open(filepath, "rt") as infile: + content = infile.read() + else: + with open(filepath) as infile: + content = infile.read() + return cls.from_str(content) + + @classmethod + def from_structure(cls, structure: Structure | Molecule) -> AimsGeometryIn: + """Construct an input file from an input structure. + + Args: + structure (Structure or Molecule): The structure for the file + + Returns: + AimsGeometryIn: The input object for the structure + """ + content_lines = [] + + if isinstance(structure, Structure): + for lv in structure.lattice.matrix: + content_lines.append(f"lattice_vector {lv[0]: .12e} {lv[1]: .12e} {lv[2]: .12e}") + + charges = structure.site_properties.get("charge", np.zeros(len(structure.species))) + magmoms = structure.site_properties.get("magmom", np.zeros(len(structure.species))) + for species, coord, charge, magmom in zip(structure.species, structure.cart_coords, charges, magmoms): + content_lines.append(f"atom {coord[0]: .12e} {coord[1]: .12e} {coord[2]: .12e} {species}") + if charge != 0: + content_lines.append(f" initial_charge {charge:.12e}") + if magmom != 0: + content_lines.append(f" initial_moment {magmom:.12e}") + + return cls(_content="\n".join(content_lines), _structure=structure) + + @property + def structure(self) -> Structure | Molecule: + """Access structure for the file""" + return self._structure + + @property + def content(self) -> str: + """Access the contents of the file""" + return self._content + + def write_file(self, directory: str | Path | None = None, overwrite: bool = False) -> None: + """Write the geometry.in file + + Args: + directory (str | Path | None): The directory to write the geometry.in file + overwrite (bool): If True allow to overwrite existing files + """ + directory = directory or Path.cwd() + + if not overwrite and (Path(directory) / "geometry.in").exists(): + raise ValueError(f"geometry.in file exists in {directory}") + + with open(f"{directory}/geometry.in", "w") as fd: + fd.write("#" + "=" * 72 + "\n") + fd.write(f"# FHI-aims geometry file: {directory}/geometry.in\n") + fd.write("# File generated from pymatgen\n") + fd.write(f"# {time.asctime()}\n") + fd.write("#" + "=" * 72 + "\n") + fd.write(self.content) + fd.write("\n") + + def as_dict(self) -> dict[str, Any]: + """Get a dictionary representation of the geometry.in file.""" + dct = {} + dct["@module"] = type(self).__module__ + dct["@class"] = type(self).__name__ + dct["content"] = self.content + dct["structure"] = self.structure + return dct + + @classmethod + def from_dict(cls, dct: dict[str, Any]) -> AimsGeometryIn: + """Initialize from dictionary. + + Args: + dct (dict[str, Any]): The MontyEncoded dictionary of the AimsGeometryIn object + + Returns: + The input object represented by the dict + """ + decoded = {key: MontyDecoder().process_decoded(val) for key, val in dct.items() if not key.startswith("@")} + + return cls(_content=decoded["content"], _structure=decoded["structure"]) + + +ALLOWED_AIMS_CUBE_TYPES = ( + "delta_density", + "spin_density", + "stm", + "total_density", + "total_density_integrable", + "long_range_potential", + "hartree_potential", + "xc_potential", + "delta_v", + "ion_dens", + "dielec_func", + "elf", +) + +ALLOWED_AIMS_CUBE_TYPES_STATE = ( + "first_order_density", + "eigenstate", + "eigenstate_imag", + "eigenstate_density", +) + +ALLOWED_AIMS_CUBE_FORMATS = ( + "cube", + "gOpenMol", + "xsf", +) + + +@dataclass +class AimsCube(MSONable): + """Class representing the FHI-aims cubes + + Attributes: + type (str): The value to be outputted as a cube file + origin (Sequence[float] or tuple[float, float, float]): The origin of the cube + edges (Sequence[Sequence[float]]): Specifies the edges of a cube: dx, dy, dz + dx (float): The length of the step in the x direction + dy (float): The length of the step in the y direction + dx (float): The length of the step in the x direction + points (Sequence[int] or tuple[int, int, int]): The number of points + along each edge + spin_state (int): The spin-channel to use either 1 or 2 + kpoint (int): The k-point to use (the index of the list printed from + `output k_point_list`) + filename (str): The filename to use + format (str): The format to output the cube file in: cube, gOpenMol, or xsf + elf_type (int): The type of electron localization function to use ( + see FHI-aims manual) + """ + + type: str = field(default_factory=str) + origin: Sequence[float] | tuple[float, float, float] = field(default_factory=lambda: [0.0, 0.0, 0.0]) + edges: Sequence[Sequence[float]] = field(default_factory=lambda: 0.1 * np.eye(3)) + points: Sequence[int] | tuple[int, int, int] = field(default_factory=lambda: [0, 0, 0]) + format: str = "cube" + spin_state: int | None = None + kpoint: int | None = None + filename: str | None = None + elf_type: int | None = None + + def __eq__(self, other: object) -> bool: + """Check if two cubes are equal to each other""" + if not isinstance(other, AimsCube): + return NotImplemented + + if self.type != other.type: + return False + + if not np.allclose(self.origin, other.origin): + return False + + if not np.allclose(self.edges, other.edges): + return False + + if not np.allclose(self.points, other.points): + return False + + if self.format != other.format: + return False + + if self.spin_state != other.spin_state: + return False + + if self.kpoint != other.kpoint: + return False + + if self.filename != other.filename: + return False + + if self.elf_type != other.elf_type: + return False + + return True + + def __post_init__(self) -> None: + """Check the inputted variables to make sure they are correct + + Raises: + ValueError: If any of the inputs is invalid + """ + split_type = self.type.split() + if split_type[0] in ALLOWED_AIMS_CUBE_TYPES: + if len(split_type) > 1: + msg = f"Cube of type {split_type[0]} can not have a state associated with it" + raise ValueError(msg) + elif split_type[0] in ALLOWED_AIMS_CUBE_TYPES_STATE: + if len(split_type) != 2: + msg = f"Cube of type {split_type[0]} must have a state associated with it" + raise ValueError(msg) + else: + raise ValueError("Cube type undefined") + + if self.format not in ALLOWED_AIMS_CUBE_FORMATS: + raise ValueError(f"{self.format} is invalid. Cube files must have a format of {ALLOWED_AIMS_CUBE_FORMATS}") + + valid_spins = (1, 2, None) + if self.spin_state not in valid_spins: + raise ValueError(f"Spin state must be one of {valid_spins}") + + if len(self.origin) != 3: + raise ValueError("The cube origin must have 3 components") + + if len(self.points) != 3: + raise ValueError("The number of points per edge must have 3 components") + + if len(self.edges) != 3: + raise ValueError("Only three cube edges can be passed") + + for edge in self.edges: + if len(edge) != 3: + raise ValueError("Each cube edge must have 3 components") + + if self.elf_type is not None and self.type != "elf": + raise ValueError("elf_type is only used when the cube type is elf. Otherwise it must be None") + + @property + def control_block(self) -> str: + """Get the block of text for the control.in file of the Cube""" + cb = f"output cube {self.type}\n" + cb += f" cube origin {self.origin[0]: .12e} {self.origin[1]: .12e} {self.origin[2]: .12e}\n" + for idx in range(3): + cb += f" cube edge {self.points[idx]} " + cb += f"{self.edges[idx][0]: .12e} " + cb += f"{self.edges[idx][1]: .12e} " + cb += f"{self.edges[idx][2]: .12e}\n" + cb += f" cube format {self.format}\n" + if self.spin_state is not None: + cb += f" cube spinstate {self.spin_state}\n" + if self.kpoint is not None: + cb += f" cube kpoint {self.kpoint}\n" + if self.filename is not None: + cb += f" cube filename {self.filename}\n" + if self.elf_type is not None: + cb += f" cube elf_type {self.elf_type}\n" + + return cb + + def as_dict(self) -> dict[str, Any]: + """Get a dictionary representation of the geometry.in file.""" + dct: dict[str, Any] = {} + dct["@module"] = type(self).__module__ + dct["@class"] = type(self).__name__ + dct["type"] = self.type + dct["origin"] = self.origin + dct["edges"] = self.edges + dct["points"] = self.points + dct["format"] = self.format + dct["spin_state"] = self.spin_state + dct["kpoint"] = self.kpoint + dct["filename"] = self.filename + dct["elf_type"] = self.elf_type + return dct + + @classmethod + def from_dict(cls, dct: dict[str, Any]) -> AimsCube: + """Initialize from dictionary. + + Args: + dct (dict[str, Any]): The MontyEncoded dictionary + + Returns: + AimsCube + """ + attrs = ("type", "origin", "edges", "points", "format", "spin_state", "kpoint", "filename", "elf_type") + decoded = {key: MontyDecoder().process_decoded(dct[key]) for key in attrs} + + return cls(**decoded) + + +@dataclass +class AimsControlIn(MSONable): + """Class representing and FHI-aims control.in file + + Attributes: + _parameters (dict[str, Any]): The parameters dictionary containing all input + flags (key) and values for the control.in file + """ + + _parameters: dict[str, Any] = field(default_factory=dict) + + def __post_init__(self) -> None: + """Initialize the output list of _parameters""" + if "output" not in self._parameters: + self._parameters["output"] = [] + + def __getitem__(self, key: str) -> Any: + """Get an input parameter + + Args: + key (str): The parameter to get + + Returns: + The setting for that parameter + + Raises: + KeyError: If the key is not in self._parameters + """ + if key not in self._parameters: + raise KeyError(f"{key} not set in AimsControlIn") + return self._parameters[key] + + def __setitem__(self, key: str, value: Any) -> None: + """Set an attribute of the class + + Args: + key (str): The parameter to get + value (Any): The value for that parameter + """ + if key == "output": + if isinstance(value, str): + value = [value] + self._parameters[key] += value + else: + self._parameters[key] = value + + def __delitem__(self, key: str) -> Any: + """Delete a parameter from the input object + + Args: + key (str): The key in the parameter to remove + + Returns: + Either the value of the deleted parameter or None if key is + not in self._parameters + """ + return self._parameters.pop(key, None) + + @property + def parameters(self) -> dict[str, Any]: + """The dictionary of input parameters for control.in""" + return self._parameters + + @parameters.setter + def parameters(self, parameters: dict[str, Any]) -> None: + """Reset a control.in inputs from a parameters dictionary + + Args: + parameters (dict[str, Any]): The new set of parameters to use + + """ + self._parameters = parameters + if "output" not in self._parameters: + self._parameters["output"] = [] + + def get_aims_control_parameter_str(self, key: str, value: Any, format: str) -> str: + """Get the string needed to add a parameter to the control.in file + + Args: + key(str): The name of the input flag + value(Any): The value to be set for the flag + format(str): The format string to apply to the value + + Returns: + The line to add to the control.in file + """ + return f"{key :35s}" + (format % value) + "\n" + + def write_file( + self, + structure: Structure | Molecule, + directory: str | Path | None = None, + verbose_header: bool = False, + overwrite: bool = False, + ) -> None: + """Writes the control.in file + + Args: + structure (Structure or Molecule): The structure to write the input + file for + directory (str or Path): The directory to write the control.in file. + If None use cwd + verbose_header (bool): If True print the input option dictionary + overwrite (bool): If True allow to overwrite existing files + + Raises: + ValueError: If a file must be overwritten and overwrite is False + ValueError: If k-grid is not provided for the periodic structures + """ + directory = directory or Path.cwd() + + if (Path(directory) / "control.in").exists() and not overwrite: + raise ValueError(f"control.in file already in {directory}") + + lim = "#" + "=" * 79 + + if isinstance(structure, Structure) and ( + "k_grid" not in self._parameters and "k_grid_density" not in self._parameters + ): + raise ValueError("k-grid must be defined for periodic systems") + + parameters = deepcopy(self._parameters) + + with open(f"{directory}/control.in", "w") as fd: + fd.write("#" + "=" * 72 + "\n") + fd.write(f"# FHI-aims geometry file: {directory}/geometry.in\n") + fd.write("# File generated from pymatgen\n") + fd.write(f"# {time.asctime()}\n") + fd.write("#" + "=" * 72 + "\n") + + if parameters["xc"] == "LDA": + parameters["xc"] = "pw-lda" + + cubes = parameters.pop("cubes", None) + + if verbose_header: + fd.write("# \n# List of parameters used to initialize the calculator:") + for p, v in parameters.items(): + s = f"# {p}:{v}\n" + fd.write(s) + fd.write(lim + "\n") + + assert not ("smearing" in parameters and "occupation_type" in parameters) + + for key, value in parameters.items(): + if key in ["species_dir", "plus_u"]: + continue + if key == "smearing": + name = parameters["smearing"][0].lower() + if name == "fermi-dirac": + name = "fermi" + width = parameters["smearing"][1] + if name == "methfessel-paxton": + order = parameters["smearing"][2] + order = " %d" % order + else: + order = "" + + fd.write(self.get_aims_control_parameter_str("occupation_type", (name, width, order), "%s %f%s")) + elif key == "output": + for output_type in value: + fd.write(self.get_aims_control_parameter_str(key, output_type, "%s")) + elif key == "vdw_correction_hirshfeld" and value: + fd.write(self.get_aims_control_parameter_str(key, "", "%s")) + elif isinstance(value, bool): + fd.write(self.get_aims_control_parameter_str(key, str(value).lower(), ".%s.")) + elif isinstance(value, (tuple, list)): + fd.write(self.get_aims_control_parameter_str(key, " ".join([str(x) for x in value]), "%s")) + elif isinstance(value, str): + fd.write(self.get_aims_control_parameter_str(key, value, "%s")) + else: + fd.write(self.get_aims_control_parameter_str(key, value, "%r")) + + if cubes: + for cube in cubes: + fd.write(cube.control_block) + + fd.write(lim + "\n\n") + fd.write(self.get_species_block(structure, self._parameters["species_dir"])) + + def get_species_block(self, structure: Structure | Molecule, species_dir: str | Path) -> str: + """Get the basis set information for a structure + + Args: + structure (Molecule or Structure): The structure to get the basis set information for + species_dir (str or Pat:): The directory to find the species files in + + Returns: + The block to add to the control.in file for the species + + Raises: + ValueError: If a file for the species is not found + """ + sb = "" + species = np.unique(structure.species) + for sp in species: + filename = f"{species_dir}/{sp.Z:02d}_{sp.symbol}_default" + if Path(filename).exists(): + with open(filename) as sf: + sb += "".join(sf.readlines()) + elif Path(f"{filename}.gz").exists(): + with gzip.open(f"{filename}.gz", "rt") as sf: + sb += "".join(sf.readlines()) + else: + raise ValueError("Species file for {sp.symbol} not found.") + + return sb + + def as_dict(self) -> dict[str, Any]: + """Get a dictionary representation of the geometry.in file.""" + dct: dict[str, Any] = {} + dct["@module"] = type(self).__module__ + dct["@class"] = type(self).__name__ + dct["parameters"] = self.parameters + return dct + + @classmethod + def from_dict(cls, dct: dict[str, Any]) -> AimsControlIn: + """Initialize from dictionary. + + Args: + dct (dict[str, Any]): The MontyEncoded dictionary + + Returns: + The AimsControlIn for dct + """ + decoded = {key: MontyDecoder().process_decoded(val) for key, val in dct.items() if not key.startswith("@")} + + return cls(_parameters=decoded["parameters"]) diff --git a/pymatgen/io/aims/outputs.py b/pymatgen/io/aims/outputs.py new file mode 100644 index 00000000000..b07cea8a33c --- /dev/null +++ b/pymatgen/io/aims/outputs.py @@ -0,0 +1,217 @@ +"""A representation of FHI-aims output (based on ASE output parser).""" +from __future__ import annotations + +from typing import TYPE_CHECKING, Any + +import numpy as np +from monty.json import MontyDecoder, MSONable + +from pymatgen.io.aims.parsers import ( + read_aims_header_info, + read_aims_header_info_from_content, + read_aims_output, + read_aims_output_from_content, +) + +if TYPE_CHECKING: + from collections.abc import Sequence + from pathlib import Path + + from emmet.core.math import Matrix3D, Vector3D + + from pymatgen.core import Molecule, Structure + + +class AimsOutput(MSONable): + """The main output file for FHI-aims.""" + + def __init__( + self, + results: Molecule | Structure | Sequence[Molecule | Structure], + metadata: dict[str, Any], + structure_summary: dict[str, Any], + ) -> None: + """AimsOutput object constructor. + + Args: + results (Molecule or Structure or Sequence[Molecule or Structure]): A list + of all images in an output file + metadata (Dict[str, Any]): The metadata of the executable used to perform + the calculation + structure_summary (Dict[str, Any]): The summary of the starting + atomic structure + """ + self._results = results + self._metadata = metadata + self._structure_summary = structure_summary + + def as_dict(self) -> dict[str, Any]: + """Create a dict representation of the outputs for MSONable.""" + d: dict[str, Any] = { + "@module": self.__class__.__module__, + "@class": self.__class__.__name__, + } + + d["results"] = self._results + d["metadata"] = self._metadata + d["structure_summary"] = self._structure_summary + return d + + @classmethod + def from_outfile(cls, outfile: str | Path) -> AimsOutput: + """Construct an AimsOutput from an output file. + + Args: + outfile: str | Path: The aims.out file to parse + + Returns: + The AimsOutput object for the output file + """ + metadata, structure_summary = read_aims_header_info(outfile) + results = read_aims_output(outfile, index=slice(0, None)) + + return cls(results, metadata, structure_summary) + + @classmethod + def from_str(cls, content: str) -> AimsOutput: + """Construct an AimsOutput from an output file. + + Args: + content (str): The content of the aims.out file + + Returns: + The AimsOutput for the output file content + """ + metadata, structure_summary = read_aims_header_info_from_content(content) + results = read_aims_output_from_content(content, index=slice(0, None)) + + return cls(results, metadata, structure_summary) + + @classmethod + def from_dict(cls, d: dict[str, Any]) -> AimsOutput: + """Construct an AimsOutput from a dictionary. + + Args: + d (dict[str, Any]): The dictionary used to create AimsOutput + + Returns: + The AimsOutput for d + """ + decoded = {k: MontyDecoder().process_decoded(v) for k, v in d.items() if not k.startswith("@")} + for struct in decoded["results"]: + struct.properties = {k: MontyDecoder().process_decoded(v) for k, v in struct.properties.items()} + + return cls( + decoded["results"], + decoded["metadata"], + decoded["structure_summary"], + ) + + def get_results_for_image(self, image_ind: int) -> Structure | Molecule: + """Get the results dictionary for a particular image or slice of images. + + Args: + image_ind (int): The index of the image to get the results for + + Returns: + The results of the image with index images_ind + """ + return self._results[image_ind] + + @property + def structure_summary(self) -> dict[str, Any]: + """The summary of the material/molecule that the calculations represent.""" + return self._structure_summary + + @property + def metadata(self) -> dict[str, Any]: + """The system metadata.""" + return self._metadata + + @property + def n_images(self) -> int: + """The number of images in results.""" + return len(self._results) + + @property + def initial_structure(self) -> Structure | Molecule: + """The initial structure for the calculations.""" + return self._structure_summary["initial_structure"] + + @property + def final_structure(self) -> Structure | Molecule: + """The final structure for the calculation.""" + return self._results[-1] + + @property + def structures(self) -> Sequence[Structure | Molecule]: + """All images in the output file.""" + return self._results + + @property + def fermi_energy(self) -> float: + """The Fermi energy for the final structure in the calculation.""" + return self.get_results_for_image(-1).properties["fermi_energy"] + + @property + def vbm(self) -> float: + """The HOMO level for the final structure in the calculation.""" + return self.get_results_for_image(-1).properties["vbm"] + + @property + def cbm(self) -> float: + """The LUMO level for the final structure in the calculation.""" + return self.get_results_for_image(-1).properties["cbm"] + + @property + def band_gap(self) -> float: + """The band gap for the final structure in the calculation.""" + return self.get_results_for_image(-1).properties["gap"] + + @property + def direct_band_gap(self) -> float: + """The direct band gap for the final structure in the calculation.""" + return self.get_results_for_image(-1).properties["direct_gap"] + + @property + def final_energy(self) -> float: + """The total energy for the final structure in the calculation.""" + return self.get_results_for_image(-1).properties["energy"] + + @property + def completed(self) -> bool: + """Did the calculation complete.""" + return len(self._results) > 0 + + @property + def aims_version(self) -> str: + """The version of FHI-aims used for the calculation.""" + return self._metadata["version_number"] + + @property + def forces(self) -> Sequence[Vector3D] | None: + """The forces for the final image of the calculation.""" + force_array = self.get_results_for_image(-1).site_properties.get("force", None) + if isinstance(force_array, np.ndarray): + return force_array.tolist() + + return force_array + + @property + def stress(self) -> Matrix3D: + """The stress for the final image of the calculation.""" + return self.get_results_for_image(-1).properties.get("stress", None) + + @property + def stresses(self) -> Sequence[Matrix3D] | None: + """The atomic virial stresses for the final image of the calculation.""" + stresses_array = self.get_results_for_image(-1).site_properties.get("atomic_virial_stress", None) + if isinstance(stresses_array, np.ndarray): + return stresses_array.tolist() + return stresses_array + + @property + def all_forces(self) -> list[list[Vector3D]]: + """The forces for all images in the calculation.""" + all_forces_array = [res.site_properties.get("force", None) for res in self._results] + return [af.tolist() if isinstance(af, np.ndarray) else af for af in all_forces_array] diff --git a/pymatgen/io/aims/parsers.py b/pymatgen/io/aims/parsers.py new file mode 100644 index 00000000000..c35d4874cba --- /dev/null +++ b/pymatgen/io/aims/parsers.py @@ -0,0 +1,1161 @@ +"""AIMS output parser, taken from ASE with modifications.""" +from __future__ import annotations + +import gzip +from dataclasses import dataclass, field +from pathlib import Path +from typing import TYPE_CHECKING, Any + +import numpy as np + +from pymatgen.core import Lattice, Molecule, Structure + +if TYPE_CHECKING: + from collections.abc import Generator, Sequence + from io import TextIOWrapper + + from emmet.core.math import Matrix3D, Vector3D + +# TARP: Originally an object, but type hinting needs this to be an int +LINE_NOT_FOUND = -1000 +EV_PER_A3_TO_KBAR = 1.60217653e-19 * 1e22 + + +class ParseError(Exception): + """Parse error during reading of a file""" + + +class AimsParseError(Exception): + """Exception raised if an error occurs when parsing an Aims output file.""" + + def __init__(self, message: str) -> None: + """Initialize the error with the message, message""" + self.message = message + super().__init__(self.message) + + +# Read aims.out files +scalar_property_to_line_key = { + "free_energy": ["| Electronic free energy"], + "number_of_iterations": ["| Number of self-consistency cycles"], + "magnetic_moment": ["N_up - N_down"], + "n_atoms": ["| Number of atoms"], + "n_bands": [ + "Number of Kohn-Sham states", + "Reducing total number of Kohn-Sham states", + "Reducing total number of Kohn-Sham states", + ], + "n_electrons": ["The structure contains"], + "n_kpts": ["| Number of k-points"], + "n_spins": ["| Number of spin channels"], + "electronic_temp": ["Occupation type:"], + "fermi_energy": ["| Chemical potential (Fermi level)"], +} + + +@dataclass +class AimsOutChunk: + """Base class for AimsOutChunks. + + Attributes: + lines (list[str]): The list of all lines in the chunk + """ + + lines: list[str] = field(default_factory=list) + + def reverse_search_for(self, keys: list[str], line_start: int = 0) -> int: + """Find the last time one of the keys appears in self.lines. + + Args: + keys (list[str]): The key strings to search for in self.lines + line_start (int): The lowest index to search for in self.lines + + Returns: + The last time one of the keys appears in self.lines + """ + for ll, line in enumerate(self.lines[line_start:][::-1]): + if any(key in line for key in keys): + return len(self.lines) - ll - 1 + + return LINE_NOT_FOUND + + def search_for_all(self, key: str, line_start: int = 0, line_end: int = -1) -> list[int]: + """Find the all times the key appears in self.lines. + + Args: + key (str): The key string to search for in self.lines + line_start (int): The first line to start the search from + line_end (int): The last line to end the search at + + Returns: + All times the key appears in the lines + """ + line_index = [] + for ll, line in enumerate(self.lines[line_start:line_end]): + if key in line: + line_index.append(ll + line_start) + return line_index + + def parse_scalar(self, property: str) -> float | None: + """Parse a scalar property from the chunk. + + Args: + property (str): The property key to parse + + Returns: + The scalar value of the property or None if not found + """ + line_start = self.reverse_search_for(scalar_property_to_line_key[property]) + + if line_start == LINE_NOT_FOUND: + return None + + line = self.lines[line_start] + return float(line.split(":")[-1].strip().split()[0]) + + +@dataclass +class AimsOutHeaderChunk(AimsOutChunk): + """The header of the aims.out file containing general information.""" + + lines: list[str] = field(default_factory=list) + _cache: dict[str, Any] = field(default_factory=dict) + + @property + def commit_hash(self) -> str: + """The commit hash for the FHI-aims version.""" + line_start = self.reverse_search_for(["Commit number"]) + if line_start == LINE_NOT_FOUND: + raise AimsParseError("This file does not appear to be an aims-output file") + + return self.lines[line_start].split(":")[1].strip() + + @property + def aims_uuid(self) -> str: + """The aims-uuid for the calculation.""" + line_start = self.reverse_search_for(["aims_uuid"]) + if line_start == LINE_NOT_FOUND: + raise AimsParseError("This file does not appear to be an aims-output file") + + return self.lines[line_start].split(":")[1].strip() + + @property + def version_number(self) -> str: + """The commit hash for the FHI-aims version.""" + line_start = self.reverse_search_for(["FHI-aims version"]) + if line_start == LINE_NOT_FOUND: + raise AimsParseError("This file does not appear to be an aims-output file") + + return self.lines[line_start].split(":")[1].strip() + + @property + def fortran_compiler(self) -> str | None: + """The fortran compiler used to make FHI-aims.""" + line_start = self.reverse_search_for(["Fortran compiler :"]) + if line_start == LINE_NOT_FOUND: + raise AimsParseError("This file does not appear to be an aims-output file") + + return self.lines[line_start].split(":")[1].split("/")[-1].strip() + + @property + def c_compiler(self) -> str | None: + """The C compiler used to make FHI-aims.""" + line_start = self.reverse_search_for(["C compiler :"]) + if line_start == LINE_NOT_FOUND: + return None + + return self.lines[line_start].split(":")[1].split("/")[-1].strip() + + @property + def fortran_compiler_flags(self) -> str | None: + """The fortran compiler flags used to make FHI-aims.""" + line_start = self.reverse_search_for(["Fortran compiler flags"]) + if line_start == LINE_NOT_FOUND: + raise AimsParseError("This file does not appear to be an aims-output file") + + return self.lines[line_start].split(":")[1].strip() + + @property + def c_compiler_flags(self) -> str | None: + """The C compiler flags used to make FHI-aims.""" + line_start = self.reverse_search_for(["C compiler flags"]) + if line_start == LINE_NOT_FOUND: + return None + + return self.lines[line_start].split(":")[1].strip() + + @property + def build_type(self) -> list[str]: + """The optional build flags passed to cmake.""" + line_end = self.reverse_search_for(["Linking against:"]) + line_inds = self.search_for_all("Using", line_end=line_end) + + return [" ".join(self.lines[ind].split()[1:]).strip() for ind in line_inds] + + @property + def linked_against(self) -> list[str]: + """Get all libraries used to link the FHI-aims executable.""" + line_start = self.reverse_search_for(["Linking against:"]) + if line_start == LINE_NOT_FOUND: + return [] + + linked_libs = [self.lines[line_start].split(":")[1].strip()] + line_start += 1 + while "lib" in self.lines[line_start]: + linked_libs.append(self.lines[line_start].strip()) + line_start += 1 + + return linked_libs + + @property + def initial_lattice(self) -> Lattice | None: + """The initial lattice vectors from the aims.out file.""" + line_start = self.reverse_search_for(["| Unit cell:"]) + if line_start == LINE_NOT_FOUND: + return None + + return Lattice( + np.array( + [[float(inp) for inp in line.split()[-3:]] for line in self.lines[line_start + 1 : line_start + 4]] + ) + ) + + @property + def initial_structure(self) -> Structure | Molecule: + """The initial structure + + Using the FHI-aims output file recreate the initial structure for + the calculation. + """ + lattice = self.initial_lattice + + line_start = self.reverse_search_for(["Atomic structure:"]) + if line_start == LINE_NOT_FOUND: + raise AimsParseError("No information about the structure in the chunk.") + + line_start += 2 + + coords = np.zeros((self.n_atoms, 3)) + species = [""] * self.n_atoms + for ll, line in enumerate(self.lines[line_start : line_start + self.n_atoms]): + inp = line.split() + coords[ll, :] = [float(pos) for pos in inp[4:7]] + species[ll] = inp[3] + + site_properties = {"charge": self.initial_charges} + if self.initial_magnetic_moments is not None: + site_properties["magmoms"] = self.initial_magnetic_moments + + if lattice: + return Structure( + lattice, + species, + coords, + np.sum(self.initial_charges), + coords_are_cartesian=True, + site_properties=site_properties, + ) + + return Molecule( + species, + coords, + np.sum(self.initial_charges), + site_properties=site_properties, + ) + + @property + def initial_charges(self) -> Sequence[float]: + """The initial charges for the structure""" + if "initial_charges" not in self._cache: + self._parse_initial_charges_and_moments() + return self._cache["initial_charges"] + + @property + def initial_magnetic_moments(self) -> Sequence[float]: + """The initial magnetic Moments""" + if "initial_magnetic_moments" not in self._cache: + self._parse_initial_charges_and_moments() + return self._cache["initial_magnetic_moments"] + + def _parse_initial_charges_and_moments(self) -> None: + """Parse the initial charges and magnetic moments from a file""" + charges = np.zeros(self.n_atoms) + magmoms = None + line_start = self.reverse_search_for(["Initial charges", "Initial moments and charges"]) + if line_start != LINE_NOT_FOUND: + line_start += 2 + magmoms = np.zeros(self.n_atoms) + for ll, line in enumerate(self.lines[line_start : line_start + self.n_atoms]): + inp = line.split() + if len(inp) == 4: + charges[ll] = float(inp[2]) + magmoms = None + else: + charges[ll] = float(inp[3]) + magmoms[ll] = float(inp[2]) + + self._cache["initial_charges"] = charges + self._cache["initial_magnetic_moments"] = magmoms + + @property + def is_md(self) -> bool: + """Is the output for a molecular dynamics calculation?""" + return self.reverse_search_for(["Complete information for previous time-step:"]) != LINE_NOT_FOUND + + @property + def is_relaxation(self) -> bool: + """Is the output for a relaxation?""" + return self.reverse_search_for(["Geometry relaxation:"]) != LINE_NOT_FOUND + + def _parse_k_points(self) -> None: + """Parse the list of k-points used in the calculation.""" + n_kpts = self.parse_scalar("n_kpts") + if n_kpts is None: + self._cache.update( + { + "k_points": None, + "k_point_weights": None, + } + ) + return + n_kpts = int(n_kpts) + + line_start = self.reverse_search_for(["| K-points in task"]) + line_end = self.reverse_search_for(["| k-point:"]) + if (line_start == LINE_NOT_FOUND) or (line_end == LINE_NOT_FOUND) or (line_end - line_start != n_kpts): + self._cache.update( + { + "k_points": None, + "k_point_weights": None, + } + ) + return + + k_points = np.zeros((n_kpts, 3)) + k_point_weights = np.zeros(n_kpts) + for kk, line in enumerate(self.lines[line_start + 1 : line_end + 1]): + k_points[kk] = [float(inp) for inp in line.split()[4:7]] + k_point_weights[kk] = float(line.split()[-1]) + + self._cache.update( + { + "k_points": k_points, + "k_point_weights": k_point_weights, + } + ) + + @property + def n_atoms(self) -> int: + """The number of atoms for the material.""" + n_atoms = self.parse_scalar("n_atoms") + if n_atoms is None: + raise AimsParseError("No information about the number of atoms in the header.") + return int(n_atoms) + + @property + def n_bands(self) -> int | None: + """The number of Kohn-Sham states for the chunk.""" + line_start = self.reverse_search_for(scalar_property_to_line_key["n_bands"]) + + if line_start == LINE_NOT_FOUND: + raise AimsParseError("No information about the number of Kohn-Sham states in the header.") + + line = self.lines[line_start] + if "| Number of Kohn-Sham states" in line: + return int(line.split(":")[-1].strip().split()[0]) + + return int(line.split()[-1].strip()[:-1]) + + @property + def n_electrons(self) -> int | None: + """The number of electrons for the chunk.""" + line_start = self.reverse_search_for(scalar_property_to_line_key["n_electrons"]) + + if line_start == LINE_NOT_FOUND: + raise AimsParseError("No information about the number of electrons in the header.") + + line = self.lines[line_start] + return int(float(line.split()[-2])) + + @property + def n_k_points(self) -> int | None: + """The number of k_ppoints for the calculation.""" + n_kpts = self.parse_scalar("n_kpts") + if n_kpts is None: + return None + + return int(n_kpts) + + @property + def n_spins(self) -> int | None: + """The number of spin channels for the chunk.""" + n_spins = self.parse_scalar("n_spins") + if n_spins is None: + raise AimsParseError("No information about the number of spin channels in the header.") + return int(n_spins) + + @property + def electronic_temperature(self) -> float: + """The electronic temperature for the chunk.""" + line_start = self.reverse_search_for(scalar_property_to_line_key["electronic_temp"]) + # TARP: Default FHI-aims value + if line_start == LINE_NOT_FOUND: + return 0.00 + + line = self.lines[line_start] + return float(line.split("=")[-1].strip().split()[0]) + + @property + def k_points(self) -> Sequence[Vector3D]: + """All k-points listed in the calculation.""" + if "k_points" not in self._cache: + self._parse_k_points() + + return self._cache["k_points"] + + @property + def k_point_weights(self) -> Sequence[float]: + """The k-point weights for the calculation.""" + if "k_point_weights" not in self._cache: + self._parse_k_points() + + return self._cache["k_point_weights"] + + @property + def header_summary(self) -> dict[str, Any]: + """Dictionary summarizing the information inside the header.""" + return { + "initial_structure": self.initial_structure, + "initial_lattice": self.initial_lattice, + "is_relaxation": self.is_relaxation, + "is_md": self.is_md, + "n_atoms": self.n_atoms, + "n_bands": self.n_bands, + "n_electrons": self.n_electrons, + "n_spins": self.n_spins, + "electronic_temperature": self.electronic_temperature, + "n_k_points": self.n_k_points, + "k_points": self.k_points, + "k_point_weights": self.k_point_weights, + } + + @property + def metadata_summary(self) -> dict[str, list[str] | str | None]: + """Dictionary containing all metadata for FHI-aims build.""" + return { + "commit_hash": self.commit_hash, + "aims_uuid": self.aims_uuid, + "version_number": self.version_number, + "fortran_compiler": self.fortran_compiler, + "c_compiler": self.c_compiler, + "fortran_compiler_flags": self.fortran_compiler_flags, + "c_compiler_flags": self.c_compiler_flags, + "build_type": self.build_type, + "linked_against": self.linked_against, + } + + +class AimsOutCalcChunk(AimsOutChunk): + """A part of the aims.out file corresponding to a single structure.""" + + def __init__(self, lines: list[str], header: AimsOutHeaderChunk) -> None: + """Construct the AimsOutCalcChunk. + + Args: + lines (list[str]): The lines used for the structure + header (.AimsOutHeaderChunk): A summary of the relevant information from + the aims.out header + """ + super().__init__(lines) + self._header = header.header_summary + self._cache: dict[str, Any] = {} + + def _parse_structure(self) -> Structure | Molecule: + """Parse a structure object from the file. + + For the given section of the aims output file generate the + calculated structure. + + Returns: + The structure or molecule for the calculation + """ + species, coords, velocities, lattice = self._parse_lattice_atom_pos() + + site_properties: dict[str, Sequence[Any]] = dict() + if len(velocities) > 0: + site_properties["velocity"] = np.array(velocities) + + results = self.results + site_prop_keys = { + "forces": "force", + "stresses": "atomic_virial_stress", + "hirshfeld_charges": "hirshfeld_charge", + "hirshfeld_volumes": "hirshfeld_volume", + "hirshfeld_atomic_dipoles": "hirshfeld_atomic_dipole", + } + properties = {prop: results[prop] for prop in results if prop not in site_prop_keys} + for prop, site_key in site_prop_keys.items(): + if prop in results: + site_properties[site_key] = results[prop] + + if lattice is not None: + return Structure( + lattice, + species, + coords, + site_properties=site_properties, + properties=properties, + coords_are_cartesian=True, + ) + return Molecule( + species, + coords, + site_properties=site_properties, + properties=properties, + ) + + def _parse_lattice_atom_pos( + self, + ) -> tuple[list[str], list[Vector3D], list[Vector3D], Lattice | None]: + """Parse the lattice and atomic positions of the structure + + Returns: + list[str]: The species symbols for the atoms in the structure + list[Vector3D]: The Cartesian coordinates of the atoms + list[Vector3D]: The velocities of the atoms + Lattice or None: The Lattice for the structure + """ + lattice_vectors = [] + velocities: list[Vector3D] = [] + species: list[str] = [] + coords: list[Vector3D] = [] + + start_keys = [ + "Atomic structure (and velocities) as used in the preceding time step", + "Updated atomic structure", + "Atomic structure that was used in the preceding time step of the wrapper", + ] + line_start = self.reverse_search_for(start_keys) + if line_start == LINE_NOT_FOUND: + species = [sp.symbol for sp in self.initial_structure.species] + coords = self.initial_structure.cart_coords.tolist() + velocities = list(self.initial_structure.site_properties.get("velocity", [])) + lattice = self.initial_lattice + + return (species, coords, velocities, lattice) + + line_start += 1 + + line_end = self.reverse_search_for( + ['Writing the current geometry to file "geometry.in.next_step"'], + line_start, + ) + if line_end == LINE_NOT_FOUND: + line_end = len(self.lines) + + for line in self.lines[line_start:line_end]: + if "lattice_vector " in line: + lattice_vectors.append([float(inp) for inp in line.split()[1:]]) + elif "atom " in line: + line_split = line.split() + species.append(line_split[4]) + coords.append([float(inp) for inp in line_split[1:4]]) + elif "velocity " in line: + velocities.append([float(inp) for inp in line.split()[1:]]) + + lattice = Lattice(lattice_vectors) if len(lattice_vectors) == 3 else None + return species, coords, velocities, lattice + + @property + def species(self) -> list[str]: + """The list of atomic symbols for all atoms in the structure""" + if "species" not in self._cache: + ( + self._cache["species"], + self._cache["coords"], + self._cache["velocities"], + self._cache["lattice"], + ) = self._parse_lattice_atom_pos() + return self._cache["species"] + + @property + def coords(self) -> list[Vector3D]: + """The cartesian coordinates of the atoms""" + if "coords" not in self._cache: + ( + self._cache["species"], + self._cache["coords"], + self._cache["velocities"], + self._cache["lattice"], + ) = self._parse_lattice_atom_pos() + return self._cache["coords"] + + @property + def velocities(self) -> list[Vector3D]: + """The velocities of the atoms""" + if "velocities" not in self._cache: + ( + self._cache["species"], + self._cache["coords"], + self._cache["velocities"], + self._cache["lattice"], + ) = self._parse_lattice_atom_pos() + return self._cache["velocities"] + + @property + def lattice(self) -> Lattice: + """The Lattice object for the structure""" + if "lattice" not in self._cache: + ( + self._cache["species"], + self._cache["coords"], + self._cache["velocities"], + self._cache["lattice"], + ) = self._parse_lattice_atom_pos() + return self._cache["lattice"] + + @property + def forces(self) -> list[Vector3D] | None: + """The forces from the aims.out file.""" + line_start = self.reverse_search_for(["Total atomic forces"]) + if line_start == LINE_NOT_FOUND: + return None + + line_start += 1 + + return np.array( + [[float(inp) for inp in line.split()[-3:]] for line in self.lines[line_start : line_start + self.n_atoms]] + ) + + @property + def stresses(self) -> list[Matrix3D] | None: + """The stresses from the aims.out file and convert to kbar.""" + line_start = self.reverse_search_for(["Per atom stress (eV) used for heat flux calculation"]) + if line_start == LINE_NOT_FOUND: + return None + line_start += 3 + stresses = [] + for line in self.lines[line_start : line_start + self.n_atoms]: + xx, yy, zz, xy, xz, yz = (float(d) for d in line.split()[2:8]) + stresses.append([xx, yy, zz, yz, xz, xy]) + + return np.array(stresses) * EV_PER_A3_TO_KBAR + + @property + def stress(self) -> Matrix3D | None: + """The stress from the aims.out file and convert to kbar.""" + from ase.stress import full_3x3_to_voigt_6_stress + + line_start = self.reverse_search_for( + [ + "Analytical stress tensor - Symmetrized", + "Numerical stress tensor", + ] + ) # Offset to relevant lines + if line_start == LINE_NOT_FOUND: + return None + + stress = [[float(inp) for inp in line.split()[2:5]] for line in self.lines[line_start + 5 : line_start + 8]] + return full_3x3_to_voigt_6_stress(stress) * EV_PER_A3_TO_KBAR + + @property + def is_metallic(self) -> bool: + """Is the system is metallic.""" + line_start = self.reverse_search_for( + ["material is metallic within the approximate finite broadening function (occupation_type)"] + ) + return line_start != LINE_NOT_FOUND + + @property + def energy(self) -> float: + """The energy from the aims.out file.""" + if self.initial_lattice is not None and self.is_metallic: + line_ind = self.reverse_search_for(["Total energy corrected"]) + else: + line_ind = self.reverse_search_for(["Total energy uncorrected"]) + if line_ind == LINE_NOT_FOUND: + raise AimsParseError("No energy is associated with the structure.") + + return float(self.lines[line_ind].split()[5]) + + @property + def dipole(self) -> Vector3D | None: + """The electric dipole moment from the aims.out file.""" + line_start = self.reverse_search_for(["Total dipole moment [eAng]"]) + if line_start == LINE_NOT_FOUND: + return None + + line = self.lines[line_start] + return np.array([float(inp) for inp in line.split()[6:9]]) + + @property + def dielectric_tensor(self) -> Matrix3D | None: + """The dielectric tensor from the aims.out file.""" + line_start = self.reverse_search_for(["PARSE DFPT_dielectric_tensor"]) + if line_start == LINE_NOT_FOUND: + return None + + # we should find the tensor in the next three lines: + lines = self.lines[line_start + 1 : line_start + 4] + + # make ndarray and return + return np.array([np.fromstring(line, sep=" ") for line in lines]) + + @property + def polarization(self) -> Vector3D | None: + """The polarization vector from the aims.out file.""" + line_start = self.reverse_search_for(["| Cartesian Polarization"]) + if line_start == LINE_NOT_FOUND: + return None + line = self.lines[line_start] + return np.array([float(s) for s in line.split()[-3:]]) + + def _parse_homo_lumo(self) -> dict[str, float]: + """Parse the HOMO/LUMO values and get band gap if periodic.""" + line_start = self.reverse_search_for(["Highest occupied state (VBM)"]) + homo = float(self.lines[line_start].split(" at ")[1].split("eV")[0].strip()) + + line_start = self.reverse_search_for(["Lowest unoccupied state (CBM)"]) + lumo = float(self.lines[line_start].split(" at ")[1].split("eV")[0].strip()) + + line_start = self.reverse_search_for(["verall HOMO-LUMO gap"]) + homo_lumo_gap = float(self.lines[line_start].split(":")[1].split("eV")[0].strip()) + + line_start = self.reverse_search_for(["Smallest direct gap"]) + if line_start == LINE_NOT_FOUND: + return { + "vbm": homo, + "cbm": lumo, + "gap": homo_lumo_gap, + "direct_gap": homo_lumo_gap, + } + + direct_gap = float(self.lines[line_start].split(":")[1].split("eV")[0].strip()) + return { + "vbm": homo, + "cbm": lumo, + "gap": homo_lumo_gap, + "direct_gap": direct_gap, + } + + def _parse_hirshfeld( + self, + ) -> None: + """Parse the Hirshfled charges volumes, and dipole moments.""" + line_start = self.reverse_search_for(["Performing Hirshfeld analysis of fragment charges and moments."]) + if line_start == LINE_NOT_FOUND: + self._cache.update( + { + "hirshfeld_charges": None, + "hirshfeld_volumes": None, + "hirshfeld_atomic_dipoles": None, + "hirshfeld_dipole": None, + } + ) + return + + line_inds = self.search_for_all("Hirshfeld charge", line_start, -1) + hirshfeld_charges = np.array([float(self.lines[ind].split(":")[1]) for ind in line_inds]) + + line_inds = self.search_for_all("Hirshfeld volume", line_start, -1) + hirshfeld_volumes = np.array([float(self.lines[ind].split(":")[1]) for ind in line_inds]) + + line_inds = self.search_for_all("Hirshfeld dipole vector", line_start, -1) + hirshfeld_atomic_dipoles = np.array( + [[float(inp) for inp in self.lines[ind].split(":")[1].split()] for ind in line_inds] + ) + + if self.lattice is None: + hirshfeld_dipole = np.sum( + hirshfeld_charges.reshape((-1, 1)) * self.coords, + axis=1, + ) + else: + hirshfeld_dipole = None + + self._cache.update( + { + "hirshfeld_charges": hirshfeld_charges, + "hirshfeld_volumes": hirshfeld_volumes, + "hirshfeld_atomic_dipoles": hirshfeld_atomic_dipoles, + "hirshfeld_dipole": hirshfeld_dipole, + } + ) + + @property + def structure(self) -> Structure | Molecule: + """The pytmagen SiteCollection of the chunk.""" + if "structure" not in self._cache: + self._cache["structure"] = self._parse_structure() + return self._cache["structure"] + + @property + def results(self) -> dict[str, Any]: + """Convert an AimsOutChunk to a Results Dictionary.""" + results = { + "energy": self.energy, + "free_energy": self.free_energy, + "forces": self.forces, + "stress": self.stress, + "stresses": self.stresses, + "magmom": self.magmom, + "dipole": self.dipole, + "fermi_energy": self.E_f, + "n_iter": self.n_iter, + "hirshfeld_charges": self.hirshfeld_charges, + "hirshfeld_dipole": self.hirshfeld_dipole, + "hirshfeld_volumes": self.hirshfeld_volumes, + "hirshfeld_atomic_dipoles": self.hirshfeld_atomic_dipoles, + "dielectric_tensor": self.dielectric_tensor, + "polarization": self.polarization, + "vbm": self.vbm, + "cbm": self.cbm, + "gap": self.gap, + "direct_gap": self.direct_gap, + } + + return {key: value for key, value in results.items() if value is not None} + + # Properties from the aims.out header + @property + def initial_structure(self) -> Structure | Molecule: + """The initial structure for the calculation""" + return self._header["initial_structure"] + + @property + def initial_lattice(self) -> Lattice | None: + """The initial Lattice of the structure""" + return self._header["initial_lattice"] + + @property + def n_atoms(self) -> int: + """The number of atoms in the structure""" + return self._header["n_atoms"] + + @property + def n_bands(self) -> int: + """The number of Kohn-Sham states for the chunk.""" + return self._header["n_bands"] + + @property + def n_electrons(self) -> int: + """The number of electrons for the chunk.""" + return self._header["n_electrons"] + + @property + def n_spins(self) -> int: + """The number of spin channels for the chunk.""" + return self._header["n_spins"] + + @property + def electronic_temperature(self) -> float: + """The electronic temperature for the chunk.""" + return self._header["electronic_temperature"] + + @property + def n_k_points(self) -> int: + """The number of k_ppoints for the calculation.""" + return self._header["n_k_points"] + + @property + def k_points(self) -> Sequence[Vector3D]: + """All k-points listed in the calculation.""" + return self._header["k_points"] + + @property + def k_point_weights(self) -> Sequence[float]: + """The k-point weights for the calculation.""" + return self._header["k_point_weights"] + + @property + def free_energy(self) -> float | None: + """The free energy of the calculation""" + return self.parse_scalar("free_energy") + + @property + def n_iter(self) -> int | None: + """The number of steps needed to converge the SCF cycle for the chunk.""" + val = self.parse_scalar("number_of_iterations") + if val is not None: + return int(val) + return None + + @property + def magmom(self) -> float | None: + """The magnetic moment of the structure""" + return self.parse_scalar("magnetic_moment") + + @property + def E_f(self) -> float | None: + """The Fermi energy""" + return self.parse_scalar("fermi_energy") + + @property + def converged(self) -> bool: + """True if the calculation is converged""" + return (len(self.lines) > 0) and ("Have a nice day." in self.lines[-5:]) + + @property + def hirshfeld_charges(self) -> Sequence[float] | None: + """The Hirshfeld charges of the system""" + if "hirshfeld_charges" not in self._cache: + self._parse_hirshfeld() + return self._cache["hirshfeld_charges"] + + @property + def hirshfeld_atomic_dipoles(self) -> Sequence[Vector3D] | None: + """The Hirshfeld atomic dipoles of the system""" + if "hirshfeld_atomic_dipoles" not in self._cache: + self._parse_hirshfeld() + return self._cache["hirshfeld_atomic_dipoles"] + + @property + def hirshfeld_volumes(self) -> Sequence[float] | None: + """The Hirshfeld atomic dipoles of the system""" + if "hirshfeld_volumes" not in self._cache: + self._parse_hirshfeld() + return self._cache["hirshfeld_volumes"] + + @property + def hirshfeld_dipole(self) -> None | Vector3D: + """The Hirshfeld dipole of the system""" + if "hirshfeld_dipole" not in self._cache: + self._parse_hirshfeld() + + return self._cache["hirshfeld_dipole"] + + @property + def vbm(self) -> float: + """The valance band maximum""" + return self._parse_homo_lumo()["vbm"] + + @property + def cbm(self) -> float: + """The conduction band minimnum""" + return self._parse_homo_lumo()["cbm"] + + @property + def gap(self) -> float: + """The band gap""" + return self._parse_homo_lumo()["gap"] + + @property + def direct_gap(self) -> float: + """The direct bandgap""" + return self._parse_homo_lumo()["direct_gap"] + + +def get_lines(content: str | TextIOWrapper) -> list[str]: + """Get a list of lines from a str or file of content + + Args: + content: the content of the file to parse + + Returns: + The list of lines + """ + if isinstance(content, str): + return [line.strip() for line in content.split("\n")] + return [line.strip() for line in content.readlines()] + + +def get_header_chunk(content: str | TextIOWrapper) -> AimsOutHeaderChunk: + """Get the header chunk for an output + + Args: + content (str or TextIOWrapper): the content to parse + + Returns: + The AimsHeaderChunk of the file + """ + lines = get_lines(content) + header = [] + stopped = False + # Stop the header once the first SCF cycle begins + for line in lines: + header.append(line) + if ( + "Convergence: q app. | density | eigen (eV) | Etot (eV)" in line + or "Begin self-consistency iteration #" in line + ): + stopped = True + break + + if not stopped: + raise ParseError("No SCF steps present, calculation failed at setup.") + + return AimsOutHeaderChunk(header) + + +def get_aims_out_chunks(content: str | TextIOWrapper, header_chunk: AimsOutHeaderChunk) -> Generator: + """Yield unprocessed chunks (header, lines) for each AimsOutChunk image. + + Args: + content (str or TextIOWrapper): the content to parse + header_chunk (AimsOutHeaderChunk): The AimsOutHeader for the calculation + + Yields: + The next AimsOutChunk + """ + lines = list(filter(lambda x: x not in header_chunk.lines, get_lines(content))) + if len(lines) == 0: + return + + # If the calculation is relaxation the updated structural information + # occurs before the re-initialization + if header_chunk.is_relaxation: + chunk_end_line = "Geometry optimization: Attempting to predict improved coordinates." + else: + chunk_end_line = "Begin self-consistency loop: Re-initialization" + + # If SCF is not converged then do not treat the next chunk_end_line as a + # new chunk until after the SCF is re-initialized + ignore_chunk_end_line = False + line_iter = lines.__iter__() + while True: + try: + line = next(line_iter).strip() # Raises StopIteration on empty file + except StopIteration: + break + + chunk_lines = [] + while chunk_end_line not in line or ignore_chunk_end_line: + chunk_lines.append(line) + # If SCF cycle not converged or numerical stresses are requested, + # don't end chunk on next Re-initialization + patterns = [ + ("Self-consistency cycle not yet converged - restarting mixer to attempt better convergence."), + ( + "Components of the stress tensor (for mathematical " + "background see comments in numerical_stress.f90)." + ), + "Calculation of numerical stress completed", + ] + if any(pattern in line for pattern in patterns): + ignore_chunk_end_line = True + elif "Begin self-consistency loop: Re-initialization" in line: + ignore_chunk_end_line = False + + try: + line = next(line_iter).strip() + except StopIteration: + break + yield AimsOutCalcChunk(chunk_lines, header_chunk) + + +def check_convergence(chunks: list[AimsOutCalcChunk], non_convergence_ok: bool = False) -> bool: + """Check if the aims output file is for a converged calculation. + + Args: + chunks(list[.AimsOutCalcChunk]): The list of chunks for the aims calculations + non_convergence_ok(bool): True if it is okay for the calculation to not be converged + chunks: list[AimsOutCalcChunk]: + non_convergence_ok: bool: (Default value = False) + + Returns: + True if the calculation is converged + """ + if not non_convergence_ok and not chunks[-1].converged: + raise ParseError("The calculation did not complete successfully") + return True + + +def read_aims_header_info_from_content( + content: str, +) -> tuple[dict[str, list[str] | None | str], dict[str, Any]]: + """Read the FHI-aims header information. + + Args: + content (str): The content of the output file to check + + Returns: + The metadata for the header of the aims calculation + """ + header_chunk = get_header_chunk(content) + return header_chunk.metadata_summary, header_chunk.header_summary + + +def read_aims_header_info( + filename: str | Path, +) -> tuple[dict[str, None | list[str] | str], dict[str, Any]]: + """Read the FHI-aims header information. + + Args: + filename(str or Path): The file to read + + Returns: + The metadata for the header of the aims calculation + """ + content = None + for path in [Path(filename), Path(f"{filename}.gz")]: + if not path.exists(): + continue + if path.suffix == ".gz": + with gzip.open(filename, "rt") as fd: + content = fd.read() + else: + with open(filename) as fd: + content = fd.read() + + if content is None: + raise FileNotFoundError(f"The requested output file {filename} does not exist.") + + return read_aims_header_info_from_content(content) + + +def read_aims_output_from_content( + content: str, index: int | slice = -1, non_convergence_ok: bool = False +) -> Structure | Molecule | Sequence[Structure | Molecule]: + """Read and aims output file from the content of a file + + Args: + content (str): The content of the file to read + index: int | slice: (Default value = -1) + non_convergence_ok: bool: (Default value = False) + + Returns: + The list of images to get + """ + header_chunk = get_header_chunk(content) + chunks = list(get_aims_out_chunks(content, header_chunk)) + + check_convergence(chunks, non_convergence_ok) + # Relaxations have an additional footer chunk due to how it is split + images = ( + [chunk.structure for chunk in chunks[:-1]] if header_chunk.is_relaxation else [chunk.atoms for chunk in chunks] + ) + return images[index] + + +def read_aims_output( + filename: str | Path, + index: int | slice = -1, + non_convergence_ok: bool = False, +) -> Structure | Molecule | Sequence[Structure | Molecule]: + """Import FHI-aims output files with all data available. + + Includes all structures for relaxations and MD runs with FHI-aims + + Args: + filename(str or Path): The file to read + index(int or slice): The index of the images to read + non_convergence_ok(bool): True if the calculations do not have to be converged + + Returns: + The list of images to get + """ + content = None + for path in [Path(filename), Path(f"{filename}.gz")]: + if not path.exists(): + continue + if path.suffix == ".gz": + with gzip.open(path, "rt") as fd: + content = fd.read() + else: + with open(path) as fd: + content = fd.read() + + if content is None: + raise FileNotFoundError(f"The requested output file {filename} does not exist.") + + return read_aims_output_from_content(content, index, non_convergence_ok) diff --git a/tests/io/aims/input_files/control.in.h2o.gz b/tests/io/aims/input_files/control.in.h2o.gz new file mode 100644 index 00000000000..208aa4ce02d Binary files /dev/null and b/tests/io/aims/input_files/control.in.h2o.gz differ diff --git a/tests/io/aims/input_files/control.in.si.gz b/tests/io/aims/input_files/control.in.si.gz new file mode 100644 index 00000000000..828dc548be4 Binary files /dev/null and b/tests/io/aims/input_files/control.in.si.gz differ diff --git a/tests/io/aims/input_files/geometry.in.h2o.gz b/tests/io/aims/input_files/geometry.in.h2o.gz new file mode 100644 index 00000000000..ab64943e765 Binary files /dev/null and b/tests/io/aims/input_files/geometry.in.h2o.gz differ diff --git a/tests/io/aims/input_files/geometry.in.h2o.ref.gz b/tests/io/aims/input_files/geometry.in.h2o.ref.gz new file mode 100644 index 00000000000..9bb21423adb Binary files /dev/null and b/tests/io/aims/input_files/geometry.in.h2o.ref.gz differ diff --git a/tests/io/aims/input_files/geometry.in.si.gz b/tests/io/aims/input_files/geometry.in.si.gz new file mode 100644 index 00000000000..48b159d8215 Binary files /dev/null and b/tests/io/aims/input_files/geometry.in.si.gz differ diff --git a/tests/io/aims/input_files/geometry.in.si.ref.gz b/tests/io/aims/input_files/geometry.in.si.ref.gz new file mode 100644 index 00000000000..c5b46dfae44 Binary files /dev/null and b/tests/io/aims/input_files/geometry.in.si.ref.gz differ diff --git a/tests/io/aims/input_files/h2o_ref.json.gz b/tests/io/aims/input_files/h2o_ref.json.gz new file mode 100644 index 00000000000..4030391d9ff Binary files /dev/null and b/tests/io/aims/input_files/h2o_ref.json.gz differ diff --git a/tests/io/aims/input_files/si_ref.json.gz b/tests/io/aims/input_files/si_ref.json.gz new file mode 100644 index 00000000000..29148f3f399 Binary files /dev/null and b/tests/io/aims/input_files/si_ref.json.gz differ diff --git a/tests/io/aims/output_files/h2o.out.gz b/tests/io/aims/output_files/h2o.out.gz new file mode 100644 index 00000000000..f4e8ba9ca0f Binary files /dev/null and b/tests/io/aims/output_files/h2o.out.gz differ diff --git a/tests/io/aims/output_files/h2o_ref.json.gz b/tests/io/aims/output_files/h2o_ref.json.gz new file mode 100644 index 00000000000..36e451bab83 Binary files /dev/null and b/tests/io/aims/output_files/h2o_ref.json.gz differ diff --git a/tests/io/aims/output_files/si.out.gz b/tests/io/aims/output_files/si.out.gz new file mode 100644 index 00000000000..20bfa80681b Binary files /dev/null and b/tests/io/aims/output_files/si.out.gz differ diff --git a/tests/io/aims/output_files/si_ref.json.gz b/tests/io/aims/output_files/si_ref.json.gz new file mode 100644 index 00000000000..edc0f16180a Binary files /dev/null and b/tests/io/aims/output_files/si_ref.json.gz differ diff --git a/tests/io/aims/parser_checks/calc_chunk.out.gz b/tests/io/aims/parser_checks/calc_chunk.out.gz new file mode 100644 index 00000000000..4205879470f Binary files /dev/null and b/tests/io/aims/parser_checks/calc_chunk.out.gz differ diff --git a/tests/io/aims/parser_checks/header_chunk.out.gz b/tests/io/aims/parser_checks/header_chunk.out.gz new file mode 100644 index 00000000000..4cb108c84e4 Binary files /dev/null and b/tests/io/aims/parser_checks/header_chunk.out.gz differ diff --git a/tests/io/aims/parser_checks/molecular_calc_chunk.out.gz b/tests/io/aims/parser_checks/molecular_calc_chunk.out.gz new file mode 100644 index 00000000000..2fcb4ba0409 Binary files /dev/null and b/tests/io/aims/parser_checks/molecular_calc_chunk.out.gz differ diff --git a/tests/io/aims/parser_checks/molecular_header_chunk.out.gz b/tests/io/aims/parser_checks/molecular_header_chunk.out.gz new file mode 100644 index 00000000000..c8d980c27be Binary files /dev/null and b/tests/io/aims/parser_checks/molecular_header_chunk.out.gz differ diff --git a/tests/io/aims/parser_checks/numerical_stress.out.gz b/tests/io/aims/parser_checks/numerical_stress.out.gz new file mode 100644 index 00000000000..b68472dd2cd Binary files /dev/null and b/tests/io/aims/parser_checks/numerical_stress.out.gz differ diff --git a/tests/io/aims/species_directory/light/01_H_default.gz b/tests/io/aims/species_directory/light/01_H_default.gz new file mode 100644 index 00000000000..986db78288f Binary files /dev/null and b/tests/io/aims/species_directory/light/01_H_default.gz differ diff --git a/tests/io/aims/species_directory/light/08_O_default.gz b/tests/io/aims/species_directory/light/08_O_default.gz new file mode 100644 index 00000000000..d9ca5514c8e Binary files /dev/null and b/tests/io/aims/species_directory/light/08_O_default.gz differ diff --git a/tests/io/aims/species_directory/light/14_Si_default.gz b/tests/io/aims/species_directory/light/14_Si_default.gz new file mode 100644 index 00000000000..c60e230eeab Binary files /dev/null and b/tests/io/aims/species_directory/light/14_Si_default.gz differ diff --git a/tests/io/aims/test_aims_inputs.py b/tests/io/aims/test_aims_inputs.py new file mode 100644 index 00000000000..9725d870d37 --- /dev/null +++ b/tests/io/aims/test_aims_inputs.py @@ -0,0 +1,230 @@ +from __future__ import annotations + +import gzip +import json +from pathlib import Path + +import numpy as np +import pytest +from monty.json import MontyDecoder, MontyEncoder + +from pymatgen.io.aims.inputs import ( + ALLOWED_AIMS_CUBE_TYPES, + ALLOWED_AIMS_CUBE_TYPES_STATE, + AimsControlIn, + AimsCube, + AimsGeometryIn, +) + +infile_dir = Path(__file__).parent / "input_files" + + +def compare_files(ref_file, test_file): + with open(test_file) as tf: + test_lines = tf.readlines()[5:] + + with gzip.open(f"{ref_file}.gz", "rt") as rf: + ref_lines = rf.readlines()[5:] + + for test_line, ref_line in zip(test_lines, ref_lines): + if "species_dir" in ref_line: + continue + + assert test_line.strip() == ref_line.strip() + + +def test_read_write_si_in(tmp_path): + si = AimsGeometryIn.from_file(infile_dir / "geometry.in.si.gz") + + in_lattice = np.array([[0.0, 2.715, 2.716], [2.717, 0.0, 2.718], [2.719, 2.720, 0.0]]) + in_coords = np.array([[0.0, 0.0, 0.0], [0.25, 0.24, 0.26]]) + + assert all(sp.symbol == "Si" for sp in si.structure.species) + assert np.allclose(si.structure.lattice.matrix, in_lattice) + assert np.allclose(si.structure.frac_coords.flatten(), in_coords.flatten()) + + si_test_from_struct = AimsGeometryIn.from_structure(si.structure) + assert si.structure == si_test_from_struct.structure + + si_test_from_struct.write_file(directory=tmp_path, overwrite=True) + with pytest.raises( + ValueError, + match="geometry.in file exists in ", + ): + si_test_from_struct.write_file(directory=tmp_path, overwrite=False) + + compare_files(infile_dir / "geometry.in.si.ref", f"{tmp_path}/geometry.in") + + with gzip.open(f"{infile_dir}/si_ref.json.gz", "rt") as si_ref_json: + si_from_dct = json.load(si_ref_json, cls=MontyDecoder) + + assert si.structure == si_from_dct.structure + + +def test_read_h2o_in(tmp_path): + h2o = AimsGeometryIn.from_file(infile_dir / "geometry.in.h2o.gz") + + in_coords = np.array( + [ + [0.0, 0.0, 0.119262], + [0.0, 0.763239, -0.477047], + [0.0, -0.763239, -0.477047], + ] + ) + + assert all(sp.symbol == symb for sp, symb in zip(h2o.structure.species, ["O", "H", "H"])) + assert np.allclose(h2o.structure.cart_coords.flatten(), in_coords.flatten()) + + h2o_test_from_struct = AimsGeometryIn.from_structure(h2o.structure) + assert h2o.structure == h2o_test_from_struct.structure + + h2o_test_from_struct.write_file(directory=tmp_path, overwrite=True) + + with pytest.raises( + ValueError, + match="geometry.in file exists in ", + ): + h2o_test_from_struct.write_file(directory=tmp_path, overwrite=False) + + compare_files(infile_dir / "geometry.in.h2o.ref", f"{tmp_path}/geometry.in") + + with gzip.open(f"{infile_dir}/h2o_ref.json.gz", "rt") as h2o_ref_json: + h2o_from_dct = json.load(h2o_ref_json, cls=MontyDecoder) + + assert h2o.structure == h2o_from_dct.structure + + +def check_wrong_type_aims_cube(type, exp_err): + with pytest.raises(ValueError, match=exp_err): + AimsCube(type=type) + + +def test_aims_cube(): + check_wrong_type_aims_cube(type="INCORRECT_TYPE", exp_err="Cube type undefined") + + for type in ALLOWED_AIMS_CUBE_TYPES_STATE: + check_wrong_type_aims_cube( + type=type, + exp_err=f"Cube of type {type} must have a state associated with it", + ) + + for type in ALLOWED_AIMS_CUBE_TYPES: + check_wrong_type_aims_cube( + type=f"{type} 1", + exp_err=f"Cube of type {type} can not have a state associated with it", + ) + + with pytest.raises( + ValueError, + match="TEST_ERR is invalid. Cube files must have a format of", + ): + AimsCube(type=ALLOWED_AIMS_CUBE_TYPES[0], format="TEST_ERR") + + with pytest.raises(ValueError, match=r"Spin state must be one of \(1, 2, None\)"): + AimsCube(type=ALLOWED_AIMS_CUBE_TYPES[0], spin_state=3) + + with pytest.raises(ValueError, match="The cube origin must have 3 components"): + AimsCube(type=ALLOWED_AIMS_CUBE_TYPES[0], origin=[0]) + + with pytest.raises(ValueError, match="Only three cube edges can be passed"): + AimsCube(type=ALLOWED_AIMS_CUBE_TYPES[0], edges=[[0.0, 0.0, 0.1]]) + + with pytest.raises(ValueError, match="Each cube edge must have 3 components"): + AimsCube( + type=ALLOWED_AIMS_CUBE_TYPES[0], + edges=[[0.0, 0.0, 0.1], [0.1, 0.0, 0.0], [0.1, 0.0]], + ) + + with pytest.raises(ValueError, match="elf_type is only used when the cube type is elf. Otherwise it must be None"): + AimsCube(type=ALLOWED_AIMS_CUBE_TYPES[0], elf_type=1) + + with pytest.raises(ValueError, match="The number of points per edge must have 3 components"): + AimsCube(type=ALLOWED_AIMS_CUBE_TYPES[0], points=[100, 100, 100, 100]) + + test_cube = AimsCube( + type="elf", + origin=[0, 0, 0], + edges=[[0.01, 0, 0], [0.0, 0.01, 0], [0.0, 0, 0.01]], + points=[100, 100, 100], + spin_state=1, + kpoint=1, + filename="test.cube", + format="cube", + elf_type=1, + ) + + test_cube_block = [ + "output cube elf", + " cube origin 0.000000000000e+00 0.000000000000e+00 0.000000000000e+00", + " cube edge 100 1.000000000000e-02 0.000000000000e+00 0.000000000000e+00", + " cube edge 100 0.000000000000e+00 1.000000000000e-02 0.000000000000e+00", + " cube edge 100 0.000000000000e+00 0.000000000000e+00 1.000000000000e-02", + " cube format cube", + " cube spinstate 1", + " cube kpoint 1", + " cube filename test.cube", + " cube elf_type 1", + "", + ] + assert test_cube.control_block == "\n".join(test_cube_block) + + test_cube_from_dict = json.loads(json.dumps(test_cube.as_dict(), cls=MontyEncoder), cls=MontyDecoder) + assert test_cube_from_dict.control_block == test_cube.control_block + + +def test_aims_control_in(tmp_path): + parameters = { + "cubes": [ + AimsCube(type="eigenstate 1", points=[10, 10, 10]), + AimsCube(type="total_density", points=[10, 10, 10]), + ], + "xc": "LDA", + "smearing": ["fermi-dirac", 0.01], + "vdw_correction_hirshfeld": True, + "compute_forces": True, + "relax_geometry": ["trm", "1e-3"], + "batch_size_limit": 200, + "species_dir": str(infile_dir.parent / "species_directory/light"), + } + + aims_control = AimsControlIn(parameters.copy()) + + for key, val in parameters.items(): + assert aims_control[key] == val + + del aims_control["xc"] + assert "xc" not in aims_control.parameters + aims_control.parameters = parameters + + h2o = AimsGeometryIn.from_file(infile_dir / "geometry.in.h2o.gz").structure + + si = AimsGeometryIn.from_file(infile_dir / "geometry.in.si.gz").structure + aims_control.write_file(h2o, directory=tmp_path, overwrite=True) + + compare_files(infile_dir / "control.in.h2o", f"{tmp_path}/control.in") + + with pytest.raises( + ValueError, + match="k-grid must be defined for periodic systems", + ): + aims_control.write_file(si, directory=tmp_path, overwrite=True) + aims_control["k_grid"] = [1, 1, 1] + + with pytest.raises( + ValueError, + match="control.in file already in ", + ): + aims_control.write_file(si, directory=tmp_path, overwrite=False) + + aims_control["output"] = "band 0 0 0 0.5 0 0.5 10 G X" + aims_control["output"] = "band 0 0 0 0.5 0.5 0.5 10 G L" + + aims_control_from_dict = json.loads(json.dumps(aims_control.as_dict(), cls=MontyEncoder), cls=MontyDecoder) + for key, val in aims_control.parameters.items(): + print("\n\n", key, "\n", val, "\n", aims_control_from_dict[key], "\n\n") + if key in ["output", "cubes"]: + np.all(aims_control_from_dict[key] == val) + assert aims_control_from_dict[key] == val + + aims_control_from_dict.write_file(si, directory=tmp_path, verbose_header=True, overwrite=True) + compare_files(infile_dir / "control.in.si", f"{tmp_path}/control.in") diff --git a/tests/io/aims/test_aims_outputs.py b/tests/io/aims/test_aims_outputs.py new file mode 100644 index 00000000000..7c174f09e1f --- /dev/null +++ b/tests/io/aims/test_aims_outputs.py @@ -0,0 +1,112 @@ +from __future__ import annotations + +import gzip +import json +from pathlib import Path + +import numpy as np +from monty.json import MontyDecoder, MontyEncoder + +from pymatgen.core import Structure +from pymatgen.io.aims.outputs import AimsOutput + +outfile_dir = Path(__file__).parent / "output_files" + + +def comp_images(test, ref): + assert test.species == ref.species + if isinstance(test, Structure): + assert np.allclose(test.lattice.matrix, ref.lattice.matrix) + assert np.allclose(test.cart_coords, ref.cart_coords) + + for key, val in test.site_properties.items(): + assert np.allclose(val, ref.site_properties[key]) + + for key, val in test.properties.items(): + assert np.allclose(val, ref.properties[key]) + + +def test_aims_output_si(): + si = AimsOutput.from_outfile(f"{outfile_dir}/si.out.gz") + with gzip.open(f"{outfile_dir}/si_ref.json.gz") as ref_file: + si_ref = json.load(ref_file, cls=MontyDecoder) + + assert si_ref.metadata == si.metadata + assert si_ref.structure_summary == si.structure_summary + + assert si_ref.n_images == si.n_images + for ii in range(si.n_images): + comp_images(si.get_results_for_image(ii), si_ref.get_results_for_image(ii)) + + +def test_aims_output_h2o(): + h2o = AimsOutput.from_outfile(f"{outfile_dir}/h2o.out.gz") + with gzip.open(f"{outfile_dir}/h2o_ref.json.gz", "rt") as ref_file: + h2o_ref = json.load(ref_file, cls=MontyDecoder) + + assert h2o_ref.metadata == h2o.metadata + assert h2o_ref.structure_summary == h2o.structure_summary + + assert h2o_ref.n_images == h2o.n_images + for ii in range(h2o.n_images): + comp_images(h2o.get_results_for_image(ii), h2o_ref.get_results_for_image(ii)) + + +def test_aims_output_si_str(): + with gzip.open(f"{outfile_dir}/si.out.gz", "rt") as si_out: + si = AimsOutput.from_str(si_out.read()) + + with gzip.open(f"{outfile_dir}/si_ref.json.gz", "rt") as ref_file: + si_ref = json.load(ref_file, cls=MontyDecoder) + + assert si_ref.metadata == si.metadata + assert si_ref.structure_summary == si.structure_summary + + assert si_ref.n_images == si.n_images + for ii in range(si.n_images): + comp_images(si.get_results_for_image(ii), si_ref.get_results_for_image(ii)) + + +def test_aims_output_h2o_str(): + with gzip.open(f"{outfile_dir}/h2o.out.gz", "rt") as h2o_out: + h2o = AimsOutput.from_str(h2o_out.read()) + + with gzip.open(f"{outfile_dir}/h2o_ref.json.gz", "rt") as ref_file: + h2o_ref = json.load(ref_file, cls=MontyDecoder) + + assert h2o_ref.metadata == h2o.metadata + assert h2o_ref.structure_summary == h2o.structure_summary + + assert h2o_ref.n_images == h2o.n_images + for ii in range(h2o.n_images): + comp_images(h2o.get_results_for_image(ii), h2o_ref.get_results_for_image(ii)) + + +def test_aims_output_si_dict(): + si = AimsOutput.from_outfile(f"{outfile_dir}/si.out.gz") + si = json.loads(json.dumps(si.as_dict(), cls=MontyEncoder), cls=MontyDecoder) + + with gzip.open(f"{outfile_dir}/si_ref.json.gz") as ref_file: + si_ref = json.load(ref_file, cls=MontyDecoder) + + assert si_ref.metadata == si.metadata + assert si_ref.structure_summary == si.structure_summary + + assert si_ref.n_images == si.n_images + for ii in range(si.n_images): + comp_images(si.get_results_for_image(ii), si_ref.get_results_for_image(ii)) + + +def test_aims_output_h2o_dict(): + h2o = AimsOutput.from_outfile(f"{outfile_dir}/h2o.out.gz") + h2o = json.loads(json.dumps(h2o.as_dict(), cls=MontyEncoder), cls=MontyDecoder) + + with gzip.open(f"{outfile_dir}/h2o_ref.json.gz", "rt") as ref_file: + h2o_ref = json.load(ref_file, cls=MontyDecoder) + + assert h2o_ref.metadata == h2o.metadata + assert h2o_ref.structure_summary == h2o.structure_summary + + assert h2o_ref.n_images == h2o.n_images + for ii in range(h2o.n_images): + comp_images(h2o.get_results_for_image(ii), h2o_ref.get_results_for_image(ii)) diff --git a/tests/io/aims/test_aims_parsers.py b/tests/io/aims/test_aims_parsers.py new file mode 100644 index 00000000000..9916eca1feb --- /dev/null +++ b/tests/io/aims/test_aims_parsers.py @@ -0,0 +1,593 @@ +# flake8: noqa +import numpy as np +from pymatgen.io.aims.parsers import ( + AimsParseError, + AimsOutChunk, + AimsOutHeaderChunk, + AimsOutCalcChunk, + EV_PER_A3_TO_KBAR, + LINE_NOT_FOUND, +) +import gzip +from pathlib import Path +from ase.stress import full_3x3_to_voigt_6_stress + +import pytest + +eps_hp = 1e-15 # The espsilon value used to compare numbers that are high-precision +eps_lp = 1e-7 # The espsilon value used to compare numbers that are low-precision + +parser_file_dir = Path(__file__).parent / "parser_checks" + + +@pytest.fixture +def default_chunk(): + lines = ["TEST", "A", "TEST", "| Number of atoms: 200 atoms"] + return AimsOutChunk(lines) + + +def test_reverse_search_for(default_chunk): + assert default_chunk.reverse_search_for(["TEST"]) == 2 + assert default_chunk.reverse_search_for(["TEST"], 1) == 2 + + assert default_chunk.reverse_search_for(["TEST A"]) == LINE_NOT_FOUND + + assert default_chunk.reverse_search_for(["A"]) == 1 + assert default_chunk.reverse_search_for(["A"], 2) == LINE_NOT_FOUND + + +def test_search_for_all(default_chunk): + assert default_chunk.search_for_all("TEST") == [0, 2] + assert default_chunk.search_for_all("TEST", 0, 1) == [0] + assert default_chunk.search_for_all("TEST", 1, -1) == [2] + + +def test_search_parse_scalar(default_chunk): + assert default_chunk.parse_scalar("n_atoms") == 200 + assert default_chunk.parse_scalar("n_electrons") is None + + +@pytest.fixture +def empty_header_chunk(): + return AimsOutHeaderChunk([]) + + +@pytest.mark.parametrize("attrname", ["n_atoms", "n_bands", "n_electrons", "n_spins", "initial_structure"]) +def test_missing_parameter(attrname, empty_header_chunk): + with pytest.raises(AimsParseError, match="No information about"): + getattr(empty_header_chunk, attrname) + + +def test_default_header_electronic_temperature(empty_header_chunk): + assert empty_header_chunk.electronic_temperature == 0.0 + + +def test_default_header_initial_lattice(empty_header_chunk): + assert empty_header_chunk.initial_lattice is None + + +def test_default_header_is_md(empty_header_chunk): + assert not empty_header_chunk.is_md + + +def test_default_header_is_relaxation(empty_header_chunk): + assert not empty_header_chunk.is_relaxation + + +def test_default_header_n_k_points(empty_header_chunk): + assert empty_header_chunk.n_k_points is None + + +def test_default_header_k_points(empty_header_chunk): + assert empty_header_chunk.k_points is None + + +def test_default_header_k_point_weights(empty_header_chunk): + assert empty_header_chunk.k_point_weights is None + + +@pytest.fixture +def initial_lattice(): + return np.array( + [ + [1.00000000, 2.70300000, 3.70300000], + [4.70300000, 2.00000000, 6.70300000], + [8.70300000, 7.70300000, 3.00000000], + ] + ) + + +@pytest.fixture +def header_chunk(): + with gzip.open(f"{parser_file_dir}/header_chunk.out.gz", "rt") as hc_file: + lines = hc_file.readlines() + + for ll, line in enumerate(lines): + lines[ll] = line.strip() + + return AimsOutHeaderChunk(lines) + + +def test_header_n_atoms(header_chunk): + assert header_chunk.n_atoms == 2 + + +def test_header_n_bands(header_chunk): + assert header_chunk.n_bands == 3 + + +def test_header_n_electrons(header_chunk): + assert header_chunk.n_electrons == 28 + + +def test_header_n_spins(header_chunk): + assert header_chunk.n_spins == 2 + + +def test_header_initial_structure(header_chunk, initial_lattice): + initial_positions = np.array([[0.000, 1.000, 2.000], [2.703, 3.703, 4.703]]) + assert len(header_chunk.initial_structure) == 2 + assert np.allclose( + header_chunk.initial_structure.lattice.matrix, + initial_lattice, + ) + assert np.allclose(header_chunk.initial_structure.cart_coords, initial_positions) + assert np.all(["Na", "Cl"] == [sp.symbol for sp in header_chunk.initial_structure.species]) + + +def test_header_initial_lattice(header_chunk, initial_lattice): + assert np.allclose(header_chunk.initial_lattice.matrix, initial_lattice) + + +def test_header_electronic_temperature(header_chunk): + assert header_chunk.electronic_temperature == 0.05 + + +def test_header_is_md(header_chunk): + assert header_chunk.is_md + + +def test_header_is_relaxation(header_chunk): + assert header_chunk.is_relaxation + + +def test_header_n_k_points(header_chunk): + assert header_chunk.n_k_points == 8 + + +@pytest.fixture +def k_points(): + return np.array( + [ + [0.000, 0.000, 0.000], + [0.000, 0.000, 0.500], + [0.000, 0.500, 0.000], + [0.000, 0.500, 0.500], + [0.500, 0.000, 0.000], + [0.500, 0.000, 0.500], + [0.500, 0.500, 0.000], + [0.500, 0.500, 0.500], + ] + ) + + +def test_header_k_point_weights( + header_chunk, +): + assert np.allclose(header_chunk.k_point_weights, np.full((8), 0.125)) + + +def test_header_k_points(header_chunk, k_points): + assert np.allclose(header_chunk.k_points, k_points) + + +def test_header_header_summary(header_chunk, k_points): + header_summary = { + "initial_structure": header_chunk.initial_structure, + "initial_lattice": header_chunk.initial_lattice, + "is_relaxation": True, + "is_md": True, + "n_atoms": 2, + "n_bands": 3, + "n_electrons": 28, + "n_spins": 2, + "electronic_temperature": 0.05, + "n_k_points": 8, + "k_points": k_points, + "k_point_weights": np.full((8), 0.125), + } + for key, val in header_chunk.header_summary.items(): + if isinstance(val, np.ndarray): + assert np.allclose(val, header_summary[key]) + else: + assert val == header_summary[key] + + +@pytest.fixture +def empty_calc_chunk(header_chunk): + return AimsOutCalcChunk([], header_chunk) + + +def test_header_transfer_n_atoms(empty_calc_chunk): + assert empty_calc_chunk.n_atoms == 2 + + +def test_header_transfer_n_bands(empty_calc_chunk): + assert empty_calc_chunk.n_bands == 3 + + +def test_header_transfer_n_electrons(empty_calc_chunk): + assert empty_calc_chunk.n_electrons == 28 + + +def test_header_transfer_n_spins(empty_calc_chunk): + assert empty_calc_chunk.n_spins == 2 + + +def test_header_transfer_initial_lattice(empty_calc_chunk, initial_lattice): + assert np.allclose(empty_calc_chunk.initial_lattice.matrix, initial_lattice) + + +def test_header_transfer_initial_structure(empty_calc_chunk, initial_lattice): + initial_positions = np.array([[0.000, 1.000, 2.000], [2.703, 3.703, 4.703]]) + + assert np.allclose(empty_calc_chunk.initial_structure.lattice.matrix, empty_calc_chunk.initial_lattice.matrix) + assert len(empty_calc_chunk.initial_structure) == 2 + assert np.allclose(empty_calc_chunk.initial_structure.lattice.matrix, initial_lattice) + assert np.allclose(empty_calc_chunk.initial_structure.cart_coords, initial_positions) + assert np.all(["Na", "Cl"] == [sp.symbol for sp in empty_calc_chunk.initial_structure.species]) + + +def test_header_transfer_electronic_temperature(empty_calc_chunk): + assert empty_calc_chunk.electronic_temperature == 0.05 + + +def test_header_transfer_n_k_points(empty_calc_chunk): + assert empty_calc_chunk.n_k_points == 8 + + +def test_header_transfer_k_point_weights(empty_calc_chunk): + assert np.allclose(empty_calc_chunk.k_point_weights, np.full((8), 0.125)) + + +def test_header_transfer_k_points(empty_calc_chunk, k_points): + assert np.allclose(empty_calc_chunk.k_points, k_points) + + +def test_default_calc_energy_raises_error(empty_calc_chunk): + with pytest.raises(AimsParseError, match="No energy is associated with the structure."): + getattr(empty_calc_chunk, "energy") + + +@pytest.mark.parametrize( + "attrname", + [ + "forces", + "stresses", + "stress", + "free_energy", + "n_iter", + "magmom", + "E_f", + "dipole", + "hirshfeld_charges", + "hirshfeld_volumes", + "hirshfeld_atomic_dipoles", + "hirshfeld_dipole", + ], +) +def test_chunk_defaults_none(attrname, empty_calc_chunk): + assert getattr(empty_calc_chunk, attrname) is None + + +def test_default_calc_is_metallic(empty_calc_chunk): + assert not empty_calc_chunk.is_metallic + + +def test_default_calc_converged(empty_calc_chunk): + assert not empty_calc_chunk.converged + + +@pytest.fixture +def calc_chunk(header_chunk): + with gzip.open(f"{parser_file_dir}/calc_chunk.out.gz", "rt") as fd: + lines = fd.readlines() + + for ll, line in enumerate(lines): + lines[ll] = line.strip() + return AimsOutCalcChunk(lines, header_chunk) + + +@pytest.fixture +def numerical_stress_chunk(header_chunk): + with gzip.open(f"{parser_file_dir}/numerical_stress.out.gz", "rt") as fd: + lines = fd.readlines() + + for ll, line in enumerate(lines): + lines[ll] = line.strip() + return AimsOutCalcChunk(lines, header_chunk) + + +def test_calc_structure(calc_chunk, initial_lattice): + initial_positions = np.array([[0.000, 1.000, 2.000], [2.703, 3.703, 4.703]]) + + assert len(calc_chunk.structure.species) == 2 + assert np.allclose(calc_chunk.structure.lattice.matrix, initial_lattice) + assert np.allclose(calc_chunk.structure.cart_coords, initial_positions) + assert np.all(["Na", "Cl"] == [sp.symbol for sp in calc_chunk.structure.species]) + + +def test_calc_forces(calc_chunk): + forces = np.array([[1.0, 2.0, 3.0], [6.0, 5.0, 4.0]]) + assert np.allclose(calc_chunk.forces, forces) + + # Different because of the constraints + assert np.allclose(calc_chunk.structure.site_properties["force"], forces) + assert np.allclose(calc_chunk.results["forces"], forces) + + +def test_calc_stresses(calc_chunk): + stresses = EV_PER_A3_TO_KBAR * np.array( + [ + [-10.0, -20.0, -30.0, -60.0, -50.0, -40.0], + [10.0, 20.0, 30.0, 60.0, 50.0, 40.0], + ] + ) + assert np.allclose(calc_chunk.stresses, stresses) + assert np.allclose(calc_chunk.structure.site_properties["atomic_virial_stress"], stresses) + assert np.allclose(calc_chunk.results["stresses"], stresses) + + +def test_calc_stress(calc_chunk): + stress = EV_PER_A3_TO_KBAR * full_3x3_to_voigt_6_stress( + np.array( + [ + [1.00000000, 2.00000000, 3.00000000], + [2.00000000, 5.00000000, 6.00000000], + [3.00000000, 6.00000000, 7.00000000], + ] + ) + ) + assert np.allclose(calc_chunk.stress, stress) + assert np.allclose(calc_chunk.structure.properties["stress"], stress) + assert np.allclose(calc_chunk.results["stress"], stress) + + +def test_calc_num_stress(numerical_stress_chunk): + stress = EV_PER_A3_TO_KBAR * full_3x3_to_voigt_6_stress( + np.array( + [ + [1.00000000, 2.00000000, 3.00000000], + [2.00000000, 5.00000000, 6.00000000], + [3.00000000, 6.00000000, 7.00000000], + ] + ) + ) + assert np.allclose(numerical_stress_chunk.stress, stress) + assert np.allclose(numerical_stress_chunk.structure.properties["stress"], stress) + assert np.allclose(numerical_stress_chunk.results["stress"], stress) + + +def test_calc_free_energy(calc_chunk): + free_energy = -3.169503986610555e05 + assert np.abs(calc_chunk.free_energy - free_energy) < eps_hp + assert np.abs(calc_chunk.structure.properties["free_energy"] - free_energy) < eps_hp + assert np.abs(calc_chunk.results["free_energy"] - free_energy) < eps_hp + + +def test_calc_energy(calc_chunk): + energy = -2.169503986610555e05 + assert np.abs(calc_chunk.energy - energy) < eps_hp + assert np.abs(calc_chunk.structure.properties["energy"] - energy) < eps_hp + assert np.abs(calc_chunk.results["energy"] - energy) < eps_hp + + +def test_calc_magnetic_moment(calc_chunk): + magmom = 0 + assert calc_chunk.magmom == magmom + assert calc_chunk.structure.properties["magmom"] == magmom + assert calc_chunk.results["magmom"] == magmom + + +def test_calc_n_iter(calc_chunk): + n_iter = 58 + assert calc_chunk.n_iter == n_iter + assert calc_chunk.results["n_iter"] == n_iter + + +def test_calc_fermi_energy(calc_chunk): + Ef = -8.24271207 + assert np.abs(calc_chunk.E_f - Ef) < eps_lp + assert np.abs(calc_chunk.results["fermi_energy"] - Ef) < eps_lp + + +def test_calc_dipole(calc_chunk): + assert calc_chunk.dipole is None + + +def test_calc_is_metallic(calc_chunk): + assert calc_chunk.is_metallic + + +def test_calc_converged(calc_chunk): + assert calc_chunk.converged + + +def test_calc_hirshfeld_charges(calc_chunk): + hirshfeld_charges = [0.20898543, -0.20840994] + assert np.allclose(calc_chunk.hirshfeld_charges, hirshfeld_charges) + assert np.allclose(calc_chunk.results["hirshfeld_charges"], hirshfeld_charges) + + +def test_calc_hirshfeld_volumes(calc_chunk): + hirshfeld_volumes = [73.39467444, 62.86011074] + assert np.allclose(calc_chunk.hirshfeld_volumes, hirshfeld_volumes) + assert np.allclose(calc_chunk.results["hirshfeld_volumes"], hirshfeld_volumes) + + +def test_calc_hirshfeld_atomic_dipoles(calc_chunk): + hirshfeld_atomic_dipoles = np.zeros((2, 3)) + assert np.allclose(calc_chunk.hirshfeld_atomic_dipoles, hirshfeld_atomic_dipoles) + assert np.allclose(calc_chunk.results["hirshfeld_atomic_dipoles"], hirshfeld_atomic_dipoles) + + +def test_calc_hirshfeld_dipole(calc_chunk): + assert calc_chunk.hirshfeld_dipole is None + + +@pytest.fixture +def molecular_header_chunk(): + with gzip.open(f"{parser_file_dir}/molecular_header_chunk.out.gz", "rt") as fd: + lines = fd.readlines() + + for ll, line in enumerate(lines): + lines[ll] = line.strip() + + return AimsOutHeaderChunk(lines) + + +@pytest.mark.parametrize( + "attrname", + [ + "k_points", + "k_point_weights", + "initial_lattice", + "n_k_points", + ], +) +def test_chunk_molecular_header_defaults_none(attrname, molecular_header_chunk): + assert getattr(molecular_header_chunk, attrname) is None + + +def test_molecular_header_n_bands(molecular_header_chunk): + assert molecular_header_chunk.n_bands == 7 + + +def test_molecular_header_initial_structure(molecular_header_chunk, molecular_positions): + assert len(molecular_header_chunk.initial_structure) == 3 + assert np.all(["O", "H", "H"] == [sp.symbol for sp in molecular_header_chunk.initial_structure.species]) + assert np.allclose( + molecular_header_chunk.initial_structure.cart_coords, + np.array( + [ + [0.00000000, 0.00000000, 0.00000000], + [0.95840000, 0.00000000, 0.00000000], + [-0.24000000, 0.92790000, 0.00000000], + ] + ), + ) + + +@pytest.fixture +def molecular_calc_chunk(molecular_header_chunk): + with gzip.open(f"{parser_file_dir}/molecular_calc_chunk.out.gz", "rt") as fd: + lines = fd.readlines() + + for ll, line in enumerate(lines): + lines[ll] = line.strip() + return AimsOutCalcChunk(lines, molecular_header_chunk) + + +@pytest.fixture +def molecular_positions(): + return np.array( + [ + [-0.00191785, -0.00243279, 0.00000000], + [0.97071531, -0.00756333, 0.00000000], + [-0.25039746, 0.93789612, 0.00000000], + ] + ) + + +def test_molecular_calc_atoms(molecular_calc_chunk, molecular_positions): + assert len(molecular_calc_chunk.structure.species) == 3 + assert np.allclose(molecular_calc_chunk.structure.cart_coords, molecular_positions) + assert np.all(["O", "H", "H"] == [sp.symbol for sp in molecular_calc_chunk.structure.species]) + + +def test_molecular_calc_forces(molecular_calc_chunk): + forces = np.array( + [ + [0.502371357164392e-03, 0.518627676606471e-03, 0.000000000000000e00], + [-0.108826758257187e-03, -0.408128912334209e-03, -0.649037698626122e-27], + [-0.393544598907207e-03, -0.110498764272267e-03, -0.973556547939183e-27], + ] + ) + assert np.allclose(molecular_calc_chunk.forces, forces) + assert np.allclose(molecular_calc_chunk.structure.site_properties["force"], forces) + assert np.allclose(molecular_calc_chunk.results["forces"], forces) + + +@pytest.mark.parametrize("attrname", ["stresses", "stress", "magmom", "E_f"]) +def test_chunk_molecular_defaults_none(attrname, molecular_calc_chunk): + assert getattr(molecular_calc_chunk, attrname) is None + + +def test_molecular_calc_free_energy(molecular_calc_chunk): + free_energy = -2.206778551123339e04 + assert np.abs(molecular_calc_chunk.free_energy - free_energy) < eps_hp + assert np.abs(molecular_calc_chunk.results["free_energy"] - free_energy) < eps_hp + assert np.abs(molecular_calc_chunk.structure.properties["free_energy"] - free_energy) < eps_hp + + +def test_molecular_calc_energy(molecular_calc_chunk): + energy = -0.206778551123339e04 + assert np.abs(molecular_calc_chunk.energy - energy) < eps_hp + assert np.abs(molecular_calc_chunk.structure.properties["energy"] - energy) < eps_hp + assert np.abs(molecular_calc_chunk.results["energy"] - energy) < eps_hp + + +def test_molecular_calc_n_iter(molecular_calc_chunk): + n_iter = 7 + assert molecular_calc_chunk.n_iter == n_iter + assert molecular_calc_chunk.results["n_iter"] == n_iter + + +def test_molecular_calc_dipole(molecular_calc_chunk): + dipole = [0.260286493869765, 0.336152447755231, 0.470003778119121e-15] + assert np.allclose(molecular_calc_chunk.dipole, dipole) + assert np.allclose(molecular_calc_chunk.structure.properties["dipole"], dipole) + assert np.allclose(molecular_calc_chunk.results["dipole"], dipole) + + +def test_molecular_calc_is_metallic(molecular_calc_chunk): + assert not molecular_calc_chunk.is_metallic + + +def test_molecular_calc_converged(molecular_calc_chunk): + assert molecular_calc_chunk.converged + + +def test_molecular_calc_hirshfeld_charges(molecular_calc_chunk): + molecular_hirshfeld_charges = np.array([-0.32053200, 0.16022630, 0.16020375]) + assert np.allclose(molecular_calc_chunk.hirshfeld_charges, molecular_hirshfeld_charges) + assert np.allclose(molecular_calc_chunk.results["hirshfeld_charges"], molecular_hirshfeld_charges) + + +def test_molecular_calc_hirshfeld_volumes(molecular_calc_chunk): + hirshfeld_volumes = np.array([21.83060659, 6.07674041, 6.07684447]) + assert np.allclose(molecular_calc_chunk.hirshfeld_volumes, hirshfeld_volumes) + assert np.allclose(molecular_calc_chunk.results["hirshfeld_volumes"], hirshfeld_volumes) + + +def test_molecular_calc_hirshfeld_atomic_dipoles(molecular_calc_chunk): + hirshfeld_atomic_dipoles = np.array( + [ + [0.04249319, 0.05486053, 0.00000000], + [0.13710134, -0.00105126, 0.00000000], + [-0.03534982, 0.13248706, 0.00000000], + ] + ) + assert np.allclose(molecular_calc_chunk.hirshfeld_atomic_dipoles, hirshfeld_atomic_dipoles) + assert np.allclose( + molecular_calc_chunk.results["hirshfeld_atomic_dipoles"], + hirshfeld_atomic_dipoles, + ) + + +def test_molecular_calc_hirshfeld_dipole(molecular_calc_chunk, molecular_positions): + molecular_hirshfeld_charges = np.array([-0.32053200, 0.16022630, 0.16020375]) + hirshfeld_dipole = np.sum(molecular_hirshfeld_charges.reshape((-1, 1)) * molecular_positions, axis=1) + + assert np.allclose(molecular_calc_chunk.hirshfeld_dipole, hirshfeld_dipole) + assert np.allclose(molecular_calc_chunk.results["hirshfeld_dipole"], hirshfeld_dipole)