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>

Author: espottesmith

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])
 

src/pymatgen/io/multiwfn.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,

src/pymatgen/io/qchem/sets.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()
         }