Add Firework for Multiwfn/QTAIM analysis of wavefunctions from Q-Chem

atomate/qchem/drones.py

@@ -9,11 +9,14 @@
 
 from monty.io import zopen
 from monty.json import jsanitize
+from monty.shutil import compress_file, decompress_file
+
 from pymatgen.apps.borg.hive import AbstractDrone
 from pymatgen.core import Molecule
 from pymatgen.io.babel import BabelMolAdaptor
 from pymatgen.io.qchem.inputs import QCInput
 from pymatgen.io.qchem.outputs import QCOutput
+from pymatgen.io.multiwfn import process_multiwfn_qtaim
 from pymatgen.symmetry.analyzer import PointGroupAnalyzer
 
 from atomate import __version__ as atomate_version
@@ -502,6 +505,28 @@ def post_process(dir_name, d):
             if len(filenames) >= 1:
                 with zopen(filenames[0], "rt") as f:
                     d["critic2"]["bonding"] = json.load(f)
+        filenames = glob.glob(os.path.join(fullpath, "CPprop.txt*"))
+        if len(filenames) >= 1:
+            filename = filenames[0]
+            recompress = False
+            if filename[-3:] == ".gz":
+                recompress = True
+                decompress_file(os.path.join(fullpath, filename))
+                filename = filename[:-3]
+
+            if "optimized_molecule" in d["output"]:
+                mol = d["output"]["optimized_molecule"]
+            else:
+                mol = d["output"]["initial_molecule"]
+
+            if not isinstance(mol, Molecule):
+                mol = Molecule.from_dict(mol)
+
+            d["output"]["qtaim"] = process_multiwfn_qtaim(mol, os.path.join(fullpath, filename))
+
+            if recompress:
+                compress_file(os.path.join(fullpath, filename))
+
 
     def validate_doc(self, d):
         """

atomate/qchem/firetasks/multiwfn.py

@@ -0,0 +1,93 @@
+# This module defines Firetask that run Multiwfn to analyze a wavefunction (*.wfn) file produced by e.g. Q-Chem.
+
+import os
+from pathlib import Path
+import subprocess
+
+from fireworks import FiretaskBase, explicit_serialize
+from monty.serialization import dumpfn, loadfn
+from monty.shutil import compress_file, decompress_file
+
+from atomate.utils.utils import env_chk, get_logger
+
+__author__ = "Evan Spotte-Smith"
+__copyright__ = "Copyright 2024, The Materials Project"
+__version__ = "0.1"
+__maintainer__ = "Evan Spotte-Smith"
+__email__ = "espottesmith@gmail.com"
+__status__ = "Alpha"
+__date__ = "07/17/2024"
+
+
+logger = get_logger(__name__)
+
+
+@explicit_serialize
+class RunMultiwfn_QTAIM(FiretaskBase):
+    """
+    Run the Multiwfn package on an electron density wavefunction (*.wfn) file produced by a Q-Chem calculation
+    to generate a CPprop file for quantum theory of atoms in molecules (QTAIM) analysis.
+
+    Required params:
+        molecule (Molecule): Molecule object of the molecule whose electron density is being analyzed
+                             Note that if prev_calc_molecule is set in the firework spec it will override
+                             the molecule required param.
+        multiwfn_command (str): Shell command to run Multiwfn
+        wfn_file (str): Name of the wavefunction file being analyzed
+    """
+
+    required_params = ["molecule", "multiwfn_command", "wfn_file"]
+
+    def run_task(self, fw_spec):
+        if fw_spec.get("prev_calc_molecule"):
+            molecule = fw_spec.get("prev_calc_molecule")
+        else:
+            molecule = self.get("molecule")
+        if molecule is None:
+            raise ValueError(
+                "No molecule passed and no prev_calc_molecule found in spec! Exiting..."
+            )
+
+        compress_at_end = False
+
+        wfn = self.get("wfn_file")
+
+        # File might be compressed
+        if not os.path.exists(wfn) and not wfn.endswith(".gz"):
+            wfn += ".gz"
+
+        if wfn[-3:] == ".gz":
+            compress_at_end = True
+            decompress_file(wfn)
+            wfn = wfn[:-3]
+
+        # This will run through an interactive Multiwfn dialogue and select the necessary options for QTAIM
+        input_script = """
+2
+2
+3
+4
+5
+6
+-1
+-9
+8
+7
+0
+-10
+q
+        """
+
+        with open("multiwfn_options.txt", "w") as file:
+            file.write(input_script)
+
+        base_command = env_chk(self.get("multiwfn_command"), fw_spec)
+
+        cmd = f"{base_command} {wfn} < multiwfn_options.txt"
+
+        logger.info(f"Running command: {cmd}")
+        return_code = subprocess.call(cmd, shell=True)
+        logger.info(f"Command {cmd} finished running with return code: {return_code}")
+
+        if compress_at_end:
+            compress_file(wfn)

atomate/qchem/firetasks/tests/test_multiwfn.py

@@ -0,0 +1,83 @@
+import json
+import os
+import shutil
+import unittest
+from shutil import which
+
+from pymatgen.io.qchem.outputs import QCOutput
+from pymatgen.io.multiwfn import process_multiwfn_qtaim
+
+from atomate.qchem.firetasks.multiwfn import RunMultiwfn_QTAIM
+from atomate.utils.testing import AtomateTest
+
+__author__ = "Evan Spotte-Smith"
+__email__ = "espottesmith@gmail.com"
+
+module_dir = os.path.dirname(os.path.abspath(__file__))
+
+
+@unittest.skipIf(not which("Multiwfn_noGUI"), "Multiwfn executable not present")
+class TestRunMultiwfn_QTAIM(AtomateTest):
+
+    def setUp(self, lpad=False):
+        os.chdir(
+            os.path.join(
+                module_dir,
+                "..",
+                "..",
+                "test_files",
+                "multiwfn_example",
+            )
+        )
+        out_file = "mol.qout.gz"
+        qc_out = QCOutput(filename=out_file)
+        self.mol = qc_out.data["initial_molecule"]
+        self.wavefunction = "WAVEFUNCTION.wfn.gz"
+        super().setUp(lpad=False)
+
+    def tearDown(self):
+        os.remove("multiwfn_options.txt")
+        os.remove("CPprop.txt")
+
+    def test_run(self):
+        os.chdir(
+            os.path.join(
+                module_dir,
+                "..",
+                "..",
+                "test_files",
+                "multiwfn_example"
+            )
+        )
+        firetask = RunMultiwfn_QTAIM(
+            molecule=self.mol,
+            multiwfn_command="Multiwfn_noGUI",
+            wfn_file="WAVEFUNCTION.wfn.gz"
+        )
+        firetask.run_task(fw_spec={})
+
+        reference = process_multiwfn_qtaim(
+            self.mol,
+            "CPprop_correct.txt"
+        )
+
+        this_output = process_multiwfn_qtaim(
+            self.mol,
+            "CPprop.txt"
+        )
+
+        for root in ["atom", "bond", "ring", "cage"]:
+            assert len(reference[root]) == len(this_output[root])
+
+            for k, v in reference[root].items():
+                assert k in this_output[root]
+
+                for kk, vv in reference[root][k].items():
+                    output_val = this_output[root][k].get(kk)
+                    if isinstance(vv, list):
+                        assert isinstance(output_val, list)
+                        assert len(vv) == len(output_val)
+                        for index, vvelem in enumerate(vv):
+                            self.assertAlmostEqual(vvelem, output_val[index])
+
+                    self.assertAlmostEqual(vv, output_val)

atomate/qchem/fireworks/core.py

@@ -11,6 +11,7 @@
 from atomate.qchem.firetasks.critic2 import ProcessCritic2, RunCritic2
 from atomate.qchem.firetasks.fragmenter import FragmentMolecule
 from atomate.qchem.firetasks.geo_transformations import PerturbGeometry
+from atomate.qchem.firetasks.multiwfn import RunMultiwfn_QTAIM
 from atomate.qchem.firetasks.parse_outputs import ProtCalcToDb, QChemToDb
 from atomate.qchem.firetasks.run_calc import RunQChemCustodian
 from atomate.qchem.firetasks.write_inputs import WriteInputFromIOSet
@@ -1036,3 +1037,92 @@ def __init__(
             )
         )
         super().__init__(t, parents=parents, name=name, **kwargs)
+
+
+class WfnAndQTAIMFW(Firework):
+    def __init__(
+        self,
+        molecule=None,
+        name="Wavefunction and QTAIM",
+        qchem_cmd=">>qchem_cmd<<",
+        multimode=">>multimode<<",
+        max_cores=">>max_cores<<",
+        multiwfn_command=">>multiwfn_command<<",
+        qchem_input_params=None,
+        db_file=None,
+        parents=None,
+        max_errors=5,
+        **kwargs
+    ):
+        """
+        Perform a Q-Chem single point calculation in order to generate a wavefunction file of the electron density
+        and then analyze the electron density critical points with the Multiwfn package.
+
+        Args:
+            molecule (Molecule): Input molecule.
+            name (str): Name for the Firework.
+            qchem_cmd (str): Command to run QChem. Supports env_chk.
+            multimode (str): Parallelization scheme, either openmp or mpi. Supports env_chk.
+            max_cores (int): Maximum number of cores to parallelize over. Supports env_chk.
+            multiwfn_command (str): Terminal command for Multiwfn. Supports env_chk.
+            qchem_input_params (dict): Specify kwargs for instantiating the input set parameters.
+                Basic uses would be to modify the default inputs of the set, such as dft_rung,
+                basis_set, pcm_dielectric, scf_algorithm, or max_scf_cycles. See
+                pymatgen/io/qchem/sets.py for default values of all input parameters. For
+                instance, if a user wanted to use a more advanced DFT functional, include a pcm
+                with a dielectric of 30, and use a larger basis, the user would set
+                qchem_input_params = {"dft_rung": 5, "pcm_dielectric": 30, "basis_set":
+                "6-311++g**"}. However, more advanced customization of the input is also
+                possible through the overwrite_inputs key which allows the user to directly
+                modify the rem, pcm, smd, and solvent dictionaries that QChemDictSet passes to
+                inputs.py to print an actual input file. For instance, if a user wanted to set
+                the sym_ignore flag in the rem section of the input file to true, then they
+                would set qchem_input_params = {"overwrite_inputs": "rem": {"sym_ignore":
+                "true"}}. Of course, overwrite_inputs could be used in conjunction with more
+                typical modifications, as seen in the test_double_FF_opt workflow test.
+            db_file (str): Path to file specifying db credentials to place output parsing.
+            parents ([Firework]): Parents of this particular Firework.
+            **kwargs: Other kwargs that are passed to Firework.__init__.
+        """
+
+        qchem_input_params = copy.deepcopy(qchem_input_params) or {}
+        qchem_input_params["output_wavefunction"] = True
+
+        input_file = "mol.qin"
+        output_file = "mol.qout"
+        t = []
+        t.append(
+            WriteInputFromIOSet(
+                molecule=molecule,
+                qchem_input_set="SinglePointSet",
+                input_file=input_file,
+                qchem_input_params=qchem_input_params,
+            )
+        )
+        t.append(
+            RunQChemCustodian(
+                qchem_cmd=qchem_cmd,
+                multimode=multimode,
+                input_file=input_file,
+                output_file=output_file,
+                max_cores=max_cores,
+                max_errors=max_errors,
+                job_type="normal",
+            )
+        )
+        t.append(
+            RunMultiwfn_QTAIM(
+                molecule=molecule,
+                multiwfn_command=multiwfn_command,
+                wfn_file="WAVEFUNCTION.wfn"
+            )
+        )
+        t.append(
+            QChemToDb(
+                db_file=db_file,
+                input_file=input_file,
+                output_file=output_file,
+                additional_fields={"task_label": name},
+            )
+        )
+        super().__init__(t, parents=parents, name=name, **kwargs)