Add multiwfn QTAIM parsing capabilities (#3926)

* Initial commit; thank you, Santiago! * Refactoring to clean up parsing; add rings and cages * More progress * Almost done with first draft; just need to put all of the pieces together * Some cleaning up; some work still to be done here * First draft, done * First steps towards tests * Progress on tests; some small bugfixes * More tests * Finished! * Remove unnecessary files * Trying to fight linters and type checkers and all that nonsense * pre-commit auto-fixes * Please, ruff, be kind * Please no more linting, please no more linting * pre-commit auto-fixes * (hopefully) just one more mypy issue * Whoops * Opinionated linters hate Union and Optional * pre-commit auto-fixes * Someday, I will meet the creators of mypy, and they will feel my rage * Small type change suggested by @DanielYang59 * pre-commit auto-fixes * Added wavefunction output as a feature of Q-Chem sets * Add back in qtaim-based bond definition (rather than distance-based definition); seems to fail for metals, but might be better elsewhere * Somehow, the linter was mad at me for things that I never touched. Injustice! * Added a third way to define bonds * mypy, of course --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
materialsproject · Jul 30, 2024 · 720fdf8 · 720fdf8
1 parent 5256fce
commit 720fdf8
Show file tree

Hide file tree

Showing 15 changed files with 21,014 additions and 3 deletions.
diff --git a/src/pymatgen/core/structure.py b/src/pymatgen/core/structure.py
@@ -1874,7 +1874,7 @@ def get_symmetric_neighbor_list(
                         redundant.append(jdx)
 
             # Delete the redundant neighbors
-            m = ~np.in1d(np.arange(len(bonds[0])), redundant)
+            m = ~np.isin(np.arange(len(bonds[0])), redundant)
             idcs_dist = np.argsort(bonds[3][m])
             bonds = (bonds[0][m][idcs_dist], bonds[1][m][idcs_dist], bonds[2][m][idcs_dist], bonds[3][m][idcs_dist])
 

diff --git a/src/pymatgen/io/multiwfn.py b/src/pymatgen/io/multiwfn.py
diff --git a/src/pymatgen/io/qchem/sets.py b/src/pymatgen/io/qchem/sets.py
@@ -150,6 +150,7 @@ def __init__(
         max_scf_cycles: int = 100,
         geom_opt_max_cycles: int = 200,
         plot_cubes: bool = False,
+        output_wavefunction: bool = False,
         nbo_params: dict | None = None,
         geom_opt: dict | None = None,
         cdft_constraints: list[list[dict]] | None = None,
@@ -229,6 +230,8 @@ def __init__(
             max_scf_cycles (int): Maximum number of SCF iterations. (Default: 100)
             geom_opt_max_cycles (int): Maximum number of geometry optimization iterations. (Default: 200)
             plot_cubes (bool): Whether to write CUBE files of the electron density. (Default: False)
+            output_wavefunction (bool): Whether to write a wavefunction (*.wfn) file of the electron density
+                (Default: False)
             nbo_params (dict): A dict containing the desired NBO params. Note that a key:value pair of
                 "version":7 will trigger NBO7 analysis. Otherwise, NBO5 analysis will be performed,
                 including if an empty dict is passed. Besides a key of "version", all other key:value
@@ -368,6 +371,7 @@ def __init__(
         self.max_scf_cycles = max_scf_cycles
         self.geom_opt_max_cycles = geom_opt_max_cycles
         self.plot_cubes = plot_cubes
+        self.output_wavefunction = output_wavefunction
         self.nbo_params = nbo_params
         self.geom_opt = geom_opt
         self.cdft_constraints = cdft_constraints
@@ -430,6 +434,10 @@ def __init__(
         if self.job_type.lower() in ["opt", "ts", "pes_scan"]:
             rem["geom_opt_max_cycles"] = str(self.geom_opt_max_cycles)
 
+        # To keep things simpler on the analysis side, don't give user option to change *.wfn file name
+        if self.output_wavefunction:
+            rem["write_wfn"] = "wavefunction"
+
         solvent_def = 0
         for a in [self.pcm_dielectric, self.isosvp_dielectric, self.smd_solvent, self.cmirs_solvent]:
             if a is not None:
@@ -644,6 +652,7 @@ def __init__(
         custom_smd: str | None = None,
         max_scf_cycles: int = 100,
         plot_cubes: bool = False,
+        output_wavefunction: bool = False,
         nbo_params: dict | None = None,
         vdw_mode: Literal["atomic", "sequential"] = "atomic",
         cdft_constraints: list[list[dict]] | None = None,
@@ -713,6 +722,8 @@ def __init__(
                 Refer to the QChem manual for further details.
             max_scf_cycles (int): Maximum number of SCF iterations. (Default: 100)
             plot_cubes (bool): Whether to write CUBE files of the electron density. (Default: False)
+            output_wavefunction (bool): Whether to write a wavefunction (*.wfn) file of the electron density
+                (Default: False)
             cdft_constraints (list of lists of dicts):
                 A list of lists of dictionaries, where each dictionary represents a charge
                 constraint in the cdft section of the QChem input file.
@@ -842,6 +853,7 @@ def __init__(
             qchem_version=qchem_version,
             max_scf_cycles=self.max_scf_cycles,
             plot_cubes=plot_cubes,
+            output_wavefunction=output_wavefunction,
             nbo_params=nbo_params,
             vdw_mode=vdw_mode,
             cdft_constraints=cdft_constraints,
@@ -868,6 +880,7 @@ def __init__(
         custom_smd: str | None = None,
         max_scf_cycles: int = 100,
         plot_cubes: bool = False,
+        output_wavefunction: bool = False,
         nbo_params: dict | None = None,
         opt_variables: dict[str, list] | None = None,
         geom_opt_max_cycles: int = 200,
@@ -940,6 +953,8 @@ def __init__(
                 explicitly requested by passing in a dictionary (empty or otherwise) for this input parameter.
                 (Default: False)
             plot_cubes (bool): Whether to write CUBE files of the electron density. (Default: False)
+            output_wavefunction (bool): Whether to write a wavefunction (*.wfn) file of the electron density
+                (Default: False)
             vdw_mode ('atomic' | 'sequential'): Method of specifying custom van der Waals radii. Applies
                 only if you are using overwrite_inputs to add a $van_der_waals section to the input.
                 In 'atomic' mode (default), dict keys represent the atomic number associated with each
@@ -1053,6 +1068,7 @@ def __init__(
             max_scf_cycles=self.max_scf_cycles,
             geom_opt_max_cycles=self.geom_opt_max_cycles,
             plot_cubes=plot_cubes,
+            output_wavefunction=output_wavefunction,
             nbo_params=nbo_params,
             geom_opt=geom_opt,
             cdft_constraints=cdft_constraints,
@@ -1077,6 +1093,7 @@ def __init__(
         custom_smd: str | None = None,
         max_scf_cycles: int = 100,
         plot_cubes: bool = False,
+        output_wavefunction: bool = False,
         nbo_params: dict | None = None,
         opt_variables: dict[str, list] | None = None,
         geom_opt_max_cycles: int = 200,
@@ -1146,6 +1163,8 @@ def __init__(
                 explicitly requested by passing in a dictionary (empty or otherwise) for this input parameter.
                 (Default: False)
             plot_cubes (bool): Whether to write CUBE files of the electron density. (Default: False)
+            output_wavefunction (bool): Whether to write a wavefunction (*.wfn) file of the electron density
+                (Default: False)
             overwrite_inputs (dict): Dictionary of QChem input sections to add or overwrite variables.
                 The currently available sections (keys) are rem, pcm,
                 solvent, smx, opt, scan, van_der_waals, and plots. The value of each key is a
@@ -1187,6 +1206,7 @@ def __init__(
             max_scf_cycles=self.max_scf_cycles,
             geom_opt_max_cycles=self.geom_opt_max_cycles,
             plot_cubes=plot_cubes,
+            output_wavefunction=output_wavefunction,
             nbo_params=nbo_params,
             geom_opt=geom_opt,
             overwrite_inputs=overwrite_inputs,
@@ -1211,6 +1231,7 @@ def __init__(
         custom_smd: str | None = None,
         max_scf_cycles: int = 100,
         plot_cubes: bool = False,
+        output_wavefunction: bool = False,
         nbo_params: dict | None = None,
         vdw_mode: Literal["atomic", "sequential"] = "atomic",
         cdft_constraints: list[list[dict]] | None = None,
@@ -1271,6 +1292,8 @@ def __init__(
                 Refer to the QChem manual for further details.
             max_scf_cycles (int): Maximum number of SCF iterations. (Default: 100)
             plot_cubes (bool): Whether to write CUBE files of the electron density. (Default: False)
+            output_wavefunction (bool): Whether to write a wavefunction (*.wfn) file of the electron density
+                (Default: False)
             vdw_mode ('atomic' | 'sequential'): Method of specifying custom van der Waals radii. Applies
                 only if you are using overwrite_inputs to add a $van_der_waals section to the input.
                 In 'atomic' mode (default), dict keys represent the atomic number associated with each
@@ -1376,6 +1399,7 @@ def __init__(
             qchem_version=qchem_version,
             max_scf_cycles=self.max_scf_cycles,
             plot_cubes=plot_cubes,
+            output_wavefunction=output_wavefunction,
             nbo_params=nbo_params,
             vdw_mode=vdw_mode,
             cdft_constraints=cdft_constraints,
@@ -1400,6 +1424,7 @@ def __init__(
         custom_smd: str | None = None,
         max_scf_cycles: int = 100,
         plot_cubes: bool = False,
+        output_wavefunction: bool = False,
         nbo_params: dict | None = None,
         vdw_mode: Literal["atomic", "sequential"] = "atomic",
         cdft_constraints: list[list[dict]] | None = None,
@@ -1460,6 +1485,8 @@ def __init__(
                 Refer to the QChem manual for further details.
             max_scf_cycles (int): Maximum number of SCF iterations. (Default: 100)
             plot_cubes (bool): Whether to write CUBE files of the electron density. (Default: False)
+            output_wavefunction (bool): Whether to write a wavefunction (*.wfn) file of the electron density
+                (Default: False)
             vdw_mode ('atomic' | 'sequential'): Method of specifying custom van der Waals radii. Applies
                 only if you are using overwrite_inputs to add a $van_der_waals section to the input.
                 In 'atomic' mode (default), dict keys represent the atomic number associated with each
@@ -1565,6 +1592,7 @@ def __init__(
             qchem_version=qchem_version,
             max_scf_cycles=self.max_scf_cycles,
             plot_cubes=plot_cubes,
+            output_wavefunction=output_wavefunction,
             nbo_params=nbo_params,
             vdw_mode=vdw_mode,
             cdft_constraints=cdft_constraints,
@@ -1597,6 +1625,7 @@ def __init__(
         custom_smd: str | None = None,
         max_scf_cycles: int = 100,
         plot_cubes: bool = False,
+        output_wavefunction: bool = False,
         nbo_params: dict | None = None,
         opt_variables: dict[str, list] | None = None,
         scan_variables: dict[str, list] | None = None,
@@ -1670,6 +1699,8 @@ def __init__(
                 Refer to the QChem manual for further details.
             max_scf_cycles (int): Maximum number of SCF iterations. (Default: 100)
             plot_cubes (bool): Whether to write CUBE files of the electron density. (Default: False)
+            output_wavefunction (bool): Whether to write a wavefunction (*.wfn) file of the electron density
+                (Default: False)
             overwrite_inputs (dict): Dictionary of QChem input sections to add or overwrite variables.
                 The currently available sections (keys) are rem, pcm,
                 solvent, smx, opt, scan, van_der_waals, and plots. The value of each key is a
@@ -1714,6 +1745,7 @@ def __init__(
             qchem_version=qchem_version,
             max_scf_cycles=self.max_scf_cycles,
             plot_cubes=plot_cubes,
+            output_wavefunction=output_wavefunction,
             nbo_params=nbo_params,
             overwrite_inputs=overwrite_inputs,
             vdw_mode=vdw_mode,

diff --git a/src/pymatgen/io/vasp/inputs.py b/src/pymatgen/io/vasp/inputs.py
@@ -2867,9 +2867,9 @@ def from_directory(
                 full_zpath = zpath(os.path.join(input_dir, fname))
                 sub_dct[fname.lower()] = ftype.from_file(full_zpath)  # type: ignore[attr-defined]
             except FileNotFoundError:  # handle the case where there is no KPOINTS file
-                sub_dct[fname.lower()] = None
+                sub_dct[fname.lower()] = None  # type: ignore[assignment]
 
-        sub_dct["optional_files"] = {
+        sub_dct["optional_files"] = {  # type: ignore[assignment]
             fname: ftype.from_file(os.path.join(input_dir, fname)) for fname, ftype in (optional_files or {}).items()
         }