Fix d2k function (#3932)

* Fix d2k function

Need to use the transpose of the inverse matrix for recipmatrix

* breaking: snake_case method args kpt_density and recip_cell (better to do early while user base small)

* d2k_recip_cell fix kpt_density type anno and doc str

* Add test for aims density conversions

* Fix mypy errors

* Add print line to see what the actual errors look like

All tests pass locally

* for supercell round coords

should fix error in tests

* remove extra print line

* refactor test_static_si_no_kgrid

---------

Co-authored-by: Janosh Riebesell <[email protected]>

Author: tpurcell90

src/pymatgen/io/aims/sets/base.py

@@ -5,7 +5,6 @@
 import copy
 import json
 import logging
-from collections.abc import Iterable
 from dataclasses import dataclass, field
 from typing import TYPE_CHECKING, Any
 from warnings import warn
@@ -19,7 +18,7 @@
 from pymatgen.io.core import InputFile, InputGenerator, InputSet
 
 if TYPE_CHECKING:
-    from collections.abc import Sequence
+    from collections.abc import Iterable, Sequence
 
     from pymatgen.util.typing import PathLike
 
@@ -234,7 +233,7 @@ def _read_previous(
             prev_dir (str or Path): The previous directory for the calculation
         """
         prev_structure: Structure | Molecule | None = None
-        prev_parameters = {}
+        prev_params = {}
         prev_results: dict[str, Any] = {}
 
         if prev_dir:
@@ -243,7 +242,7 @@ def _read_previous(
             # jobflow_remote)
             split_prev_dir = str(prev_dir).split(":")[-1]
             with open(f"{split_prev_dir}/parameters.json") as param_file:
-                prev_parameters = json.load(param_file, cls=MontyDecoder)
+                prev_params = json.load(param_file, cls=MontyDecoder)
 
             try:
                 aims_output: Sequence[Structure | Molecule] = read_aims_output(
@@ -256,7 +255,7 @@ def _read_previous(
             except (IndexError, AimsParseError):
                 pass
 
-        return prev_structure, prev_parameters, prev_results
+        return prev_structure, prev_params, prev_results
 
     @staticmethod
     def _get_properties(
@@ -308,12 +307,9 @@ def _get_input_parameters(
         Returns:
             dict: The input object
         """
-        # Get the default configuration
-        # FHI-aims recommends using their defaults so bare-bones default parameters
-        parameters: dict[str, Any] = {
-            "xc": "pbe",
-            "relativistic": "atomic_zora scalar",
-        }
+        # Get the default config
+        # FHI-aims recommends using their defaults so bare-bones default params
+        params: dict[str, Any] = {"xc": "pbe", "relativistic": "atomic_zora scalar"}
 
         # Override default parameters with previous parameters
         prev_parameters = {} if prev_parameters is None else copy.deepcopy(prev_parameters)
@@ -327,25 +323,25 @@ def _get_input_parameters(
                 kpt_settings["density"] = density
 
         parameter_updates = self.get_parameter_updates(structure, prev_parameters)
-        parameters = recursive_update(parameters, parameter_updates)
+        params = recursive_update(params, parameter_updates)
 
         # Override default parameters with user_params
-        parameters = recursive_update(parameters, self.user_params)
-        if ("k_grid" in parameters) and ("density" in kpt_settings):
+        params = recursive_update(params, self.user_params)
+        if ("k_grid" in params) and ("density" in kpt_settings):
             warn(
                 "WARNING: the k_grid is set in user_params and in the kpt_settings,"
                 " using the one passed in user_params.",
                 stacklevel=1,
             )
-        elif isinstance(structure, Structure) and ("k_grid" not in parameters):
+        elif isinstance(structure, Structure) and ("k_grid" not in params):
             density = kpt_settings.get("density", 5.0)
             even = kpt_settings.get("even", True)
-            parameters["k_grid"] = self.d2k(structure, density, even)
-        elif isinstance(structure, Molecule) and "k_grid" in parameters:
+            params["k_grid"] = self.d2k(structure, density, even)
+        elif isinstance(structure, Molecule) and "k_grid" in params:
             warn("WARNING: removing unnecessary k_grid information", stacklevel=1)
-            del parameters["k_grid"]
+            del params["k_grid"]
 
-        return parameters
+        return params
 
     def get_parameter_updates(
         self,
@@ -366,7 +362,7 @@ def get_parameter_updates(
     def d2k(
         self,
         structure: Structure,
-        kptdensity: float | list[float] = 5.0,
+        kpt_density: float | tuple[float, float, float] = 5.0,
         even: bool = True,
     ) -> Iterable[float]:
         """Convert k-point density to Monkhorst-Pack grid size.
@@ -376,15 +372,15 @@ def d2k(
         Args:
             structure (Structure): Contains unit cell and
                 information about boundary conditions.
-            kptdensity (float | list[float]): Required k-point
+            kpt_density (float | list[float]): Required k-point
                 density.  Default value is 5.0 point per Ang^-1.
             even (bool): Round up to even numbers.
 
         Returns:
             dict: Monkhorst-Pack grid size in all directions
         """
-        recipcell = structure.lattice.inv_matrix
-        return self.d2k_recipcell(recipcell, structure.lattice.pbc, kptdensity, even)
+        recip_cell = structure.lattice.inv_matrix.transpose()
+        return self.d2k_recip_cell(recip_cell, structure.lattice.pbc, kpt_density, even)
 
     def k2d(self, structure: Structure, k_grid: np.ndarray[int]):
         """Generate the kpoint density in each direction from given k_grid.
@@ -398,36 +394,36 @@ def k2d(self, structure: Structure, k_grid: np.ndarray[int]):
         Returns:
             dict: Density of kpoints in each direction. result.mean() computes average density
         """
-        recipcell = structure.lattice.inv_matrix
-        densities = k_grid / (2 * np.pi * np.sqrt((recipcell**2).sum(axis=1)))
+        recip_cell = structure.lattice.inv_matrix.transpose()
+        densities = k_grid / (2 * np.pi * np.sqrt((recip_cell**2).sum(axis=1)))
         return np.array(densities)
 
     @staticmethod
-    def d2k_recipcell(
-        recipcell: np.ndarray,
+    def d2k_recip_cell(
+        recip_cell: np.ndarray,
         pbc: Sequence[bool],
-        kptdensity: float | Sequence[float] = 5.0,
+        kpt_density: float | tuple[float, float, float] = 5.0,
         even: bool = True,
     ) -> Sequence[int]:
         """Convert k-point density to Monkhorst-Pack grid size.
 
         Args:
-            recipcell (Cell): The reciprocal cell
+            recip_cell (Cell): The reciprocal cell
             pbc (Sequence[bool]): If element of pbc is True
                 then system is periodic in that direction
-            kptdensity (float or list[floats]): Required k-point
-                density.  Default value is 3.5 point per Ang^-1.
+            kpt_density (float or list[floats]): Required k-point
+                density. Default value is 5 points per Ang^-1.
         even(bool): Round up to even numbers.
 
         Returns:
             dict: Monkhorst-Pack grid size in all directions
         """
-        if not isinstance(kptdensity, Iterable):
-            kptdensity = 3 * [float(kptdensity)]
+        if isinstance(kpt_density, float):
+            kpt_density = (kpt_density, kpt_density, kpt_density)
         kpts: list[int] = []
         for i in range(3):
             if pbc[i]:
-                k = 2 * np.pi * np.sqrt((recipcell[i] ** 2).sum()) * float(kptdensity[i])
+                k = 2 * np.pi * np.sqrt((recip_cell[i] ** 2).sum()) * float(kpt_density[i])
                 if even:
                     kpts.append(2 * int(np.ceil(k / 2)))
                 else:

src/pymatgen/util/testing/aims.py

@@ -84,6 +84,7 @@ def comp_system(
     generator_cls: type,
     properties: list[str] | None = None,
     prev_dir: str | None | Path = None,
+    user_kpt_settings: dict[str, Any] | None = None,
 ) -> None:
     """Compare files generated by tests with ones in reference directories.
 
@@ -96,16 +97,24 @@ def comp_system(
         generator_cls (type): The class of the generator
         properties (list[str] | None): The list of properties to calculate
         prev_dir (str | Path | None): The previous directory to pull outputs from
+        user_kpt_settings (dict[str, Any] | None): settings for k-point density in FHI-aims
 
     Raises:
         AssertionError: If the input files are not the same
     """
+    if user_kpt_settings is None:
+        user_kpt_settings = {}
+
     k_point_density = user_params.pop("k_point_density", 20)
 
     try:
-        generator = generator_cls(user_params=user_params, k_point_density=k_point_density)
+        generator = generator_cls(
+            user_params=user_params,
+            k_point_density=k_point_density,
+            user_kpoints_settings=user_kpt_settings,
+        )
     except TypeError:
-        generator = generator_cls(user_params=user_params)
+        generator = generator_cls(user_params=user_params, user_kpoints_settings=user_kpt_settings)
 
     input_set = generator.get_input_set(structure, prev_dir, properties)
     input_set.write_input(work_dir / test_name)
@@ -117,7 +126,7 @@ def compare_single_files(ref_file: str | Path, test_file: str | Path) -> None:
     """Compare single files generated by tests with ones in reference directories.
 
     Args:
-        ref_file (str | Path): The reference file to cmpare against
+        ref_file (str | Path): The reference file to compare against
         test_file (str | Path): The file to compare against the reference
 
     Raises:

tests/io/aims/aims_input_generator_ref/static-no-kgrid-si/control.in.gz

@@ -4,7 +4,7 @@
   "species_dir": "/home/tpurcell/git/atomate2/tests/aims/species_dir/light",
   "k_grid": [
     12,
-    12,
-    12
+    6,
+    4
   ]
 }

tests/io/aims/aims_input_generator_ref/static-no-kgrid-si/geometry.in.gz

@@ -19,7 +19,11 @@ def test_static_si(tmp_path):
 
 def test_static_si_no_kgrid(tmp_path):
     parameters = {"species_dir": "light"}
-    comp_system(Si, parameters, "static-no-kgrid-si", tmp_path, ref_path, StaticSetGenerator)
+    Si_supercell = Si.make_supercell([1, 2, 3], in_place=False)
+    for site in Si_supercell:
+        # round site.coords to ignore floating point errors
+        site.coords = [round(x, 15) for x in site.coords]
+    comp_system(Si_supercell, parameters, "static-no-kgrid-si", tmp_path, ref_path, StaticSetGenerator)
 
 
 def test_static_o2(tmp_path):