Merge pull request #593 from espottesmith/qchem_update

Qchem update

Author: itsduowang

atomate/qchem/drones.py

@@ -215,10 +215,10 @@ def generate_doc(self, dir_name, qcinput_files, qcoutput_files, multirun):
             elif "ESP" in d_calc_final:
                 d["output"]["esp"] = d_calc_final["ESP"][-1]
 
-            if (
-                d["output"]["job_type"] == "opt"
-                or d["output"]["job_type"] == "optimization"
-            ):
+            if "nbo_data" in d_calc_final:
+                d["output"]["nbo"] = d_calc_final["nbo_data"]
+
+            if d["output"]["job_type"] in ["opt", "optimization", "ts"]:
                 if "molecule_from_optimized_geometry" in d_calc_final:
                     d["output"]["optimized_molecule"] = d_calc_final[
                         "molecule_from_optimized_geometry"
@@ -231,28 +231,44 @@ def generate_doc(self, dir_name, qcinput_files, qcoutput_files, multirun):
                         d_calc_final["opt_constraint"][0],
                         float(d_calc_final["opt_constraint"][6]),
                     ]
-            if (
-                d["output"]["job_type"] == "freq"
-                or d["output"]["job_type"] == "frequency"
-            ):
+            if d["output"]["job_type"] in ["freq", "frequency"]:
                 d["output"]["frequencies"] = d_calc_final["frequencies"]
+                d["output"]["frequency_modes"] = d_calc_final["frequency_mode_vectors"]
                 d["output"]["enthalpy"] = d_calc_final["total_enthalpy"]
                 d["output"]["entropy"] = d_calc_final["total_entropy"]
-                if (
-                    d["input"]["job_type"] == "opt"
-                    or d["input"]["job_type"] == "optimization"
-                ):
+                if d["input"]["job_type"] in ["opt", "optimization", "ts"]:
                     d["output"]["optimized_molecule"] = d_calc_final["initial_molecule"]
                     d["output"]["final_energy"] = d["calcs_reversed"][1]["final_energy"]
+                # For frequency-first calcs
+                else:
+                    if len(d["calcs_reversed"]) > 1:
+                        first_calc = d["calcs_reversed"][-1]["input"]["rem"]["job_type"]
+                        second_calc = d["calcs_reversed"][-2]["input"]["rem"]["job_type"]
+                        if first_calc in ["freq", "frequency"] and second_calc in ["opt", "optimization", "ts"]:
+                            d["output"]["optimized_molecule"] = d_calc_final["initial_molecule"]
+                            d["output"]["final_energy"] = d["calcs_reversed"][1]["final_energy"]
+
+            if d["output"]["job_type"] == "pes_scan":
+                d["output"]["scan_energies"] = d_calc_final.get("scan_energies")
+                d["input"]["scan_variables"] = d_calc_final.get("scan_variables")
+                d["output"]["scan_geometries"] = d_calc_final.get("optimized_geometries")
+                d["output"]["scan_molecules"] = d_calc_final.get("molecules_from_optimized_geometries")
+
+            if d["output"]["job_type"] == "force":
+                d["output"]["gradients"] = d_calc_final["gradients"][0]
+                if d_calc_final["pcm_gradients"] is not None:
+                    d["output"]["pcm_gradients"] = d_calc_final["pcm_gradients"][0]
+                if d_calc_final["CDS_gradients"] is not None:
+                    d["output"]["CDS_gradients"] = d_calc_final["CDS_gradients"][0]
 
             opt_trajectory = []
             calcs = copy.deepcopy(d["calcs_reversed"])
             calcs.reverse()
             for calc in calcs:
                 job_type = calc["input"]["rem"]["job_type"]
-                if job_type == "opt" or job_type == "optimization":
-                    for ii, geom in enumerate(calc["geometries"]):
-                        site_properties = {"Mulliken": calc["Mulliken"][ii]}
+                if job_type in ["opt", "optimization", "ts"]:
+                    for ii,geom in enumerate(calc["geometries"]):
+                        site_properties = {"Mulliken":calc["Mulliken"][ii]}
                         if "RESP" in calc:
                             site_properties["RESP"] = calc["RESP"][ii]
                         mol = Molecule(
@@ -310,30 +326,36 @@ def generate_doc(self, dir_name, qcinput_files, qcoutput_files, multirun):
 
             d["state"] = "successful" if d_calc_final["completion"] else "unsuccessful"
             if "special_run_type" in d:
-                if d["special_run_type"] == "frequency_flattener":
+                if d["special_run_type"] in ["frequency_flattener", "ts_frequency_flattener"]:
                     if d["state"] == "successful":
-                        orig_num_neg_freq = sum(
-                            1
-                            for freq in d["calcs_reversed"][-2]["frequencies"]
-                            if freq < 0
-                        )
-                        orig_energy = d_calc_init["final_energy"]
-                        final_num_neg_freq = sum(
-                            1 for freq in d_calc_final["frequencies"] if freq < 0
-                        )
+                        orig_num_neg_freq = None
+                        for calc in d["calcs_reversed"][::-1]:
+                            if "frequencies" in calc:
+                                orig_num_neg_freq = sum(1 for freq in calc["frequencies"] if freq < 0)
+                                break
+                        orig_energy = None
+                        for calc in d["calcs_reversed"][::-1]:
+                            if "final_energy" in calc:
+                                if calc["final_energy"] is not None:
+                                    orig_energy = calc["final_energy"]
+                                    break
+                        final_num_neg_freq = sum(1 for freq in d_calc_final["frequencies"] if freq < 0)
                         final_energy = d["calcs_reversed"][1]["final_energy"]
-                        d["num_frequencies_flattened"] = (
-                            orig_num_neg_freq - final_num_neg_freq
-                        )
-                        if final_num_neg_freq > 0:  # If a negative frequency remains,
-                            # and it's too large to ignore,
-                            if (
-                                final_num_neg_freq > 1
-                                or abs(d["output"]["frequencies"][0]) >= 15.0
-                            ):
-                                d[
-                                    "state"
-                                ] = "unsuccessful"  # then the flattening was unsuccessful
+                        if orig_num_neg_freq is None:
+                            d["num_frequencies_flattened"] = 0
+                        else:
+                            d["num_frequencies_flattened"] = orig_num_neg_freq - final_num_neg_freq
+
+                        if d["special_run_type"] == "frequency_flattener":
+                            if final_num_neg_freq > 0: # If a negative frequency remains,
+                                # and it's too large to ignore,
+                                if final_num_neg_freq > 1 or abs(d["output"]["frequencies"][0]) >= 15.0:
+                                    d["state"] = "unsuccessful" # then the flattening was unsuccessful
+                        else:
+                            if final_num_neg_freq > 1: # If a negative frequency remains,
+                                # and it's too large to ignore,
+                                if final_num_neg_freq > 2 or abs(d["output"]["frequencies"][1]) >= 15.0:
+                                    d["state"] = "unsuccessful" # then the flattening was unsuccessful
                         if final_energy > orig_energy:
                             d["warnings"]["energy_increased"] = True
 

atomate/qchem/firetasks/geo_transformations.py

@@ -1,4 +1,9 @@
-# This module defines firetasks for writing QChem input files
+# coding: utf-8
+
+
+# This module defines firetasks for modifying molecular geometries
+
+import copy
 
 import numpy as np
 from fireworks import FiretaskBase, FWAction, explicit_serialize
@@ -58,3 +63,36 @@ def run_task(self, fw_spec):
         update_spec = {"prev_calc_molecule": rotated_mol}
 
         return FWAction(update_spec=update_spec)
+
+
+@explicit_serialize
+class PerturbGeometry(FiretaskBase):
+    """
+    Modify a Molecule geometry.
+    """
+
+    required_params = list()
+    optional_params = ["molecule", "mode", "scale"]
+
+    def run_task(self, fw_spec):
+        if self.get("molecule"):
+            mol = self.get("molecule")
+        elif fw_spec.get("prev_calc_molecule"):
+            mol = fw_spec.get("prev_calc_molecule")
+        else:
+            raise KeyError("No molecule present; add as an optional param or check fw_spec")
+
+        if self.get("mode") is not None:
+            mode = self.get("mode")
+        else:
+            raise KeyError("No mode present; add as an optional param or check fw_spec")
+
+        mol_copy = copy.deepcopy(mol)
+
+        for ii in range(len(mol)):
+            vec = np.array(mode[ii])
+            mol_copy.translate_sites(indices=[ii], vector=vec * self.get("scale", 1.0))
+
+        update_spec = {"prev_calc_molecule": mol_copy}
+
+        return FWAction(update_spec=update_spec)

atomate/qchem/firetasks/parse_outputs.py

@@ -41,17 +41,14 @@ class QChemToDb(FiretaskBase):
             of this key in the fw_spec.
         multirun (bool): Whether the job to parse includes multiple
             calculations in one input / output pair.
+        runs (list): Series of file suffixes that the Drone should look for
+            when parsing output.
     """
 
     optional_params = [
-        "calc_dir",
-        "calc_loc",
-        "input_file",
-        "output_file",
-        "additional_fields",
-        "db_file",
-        "fw_spec_field",
-        "multirun",
+        "calc_dir", "calc_loc", "input_file", "output_file",
+        "additional_fields", "db_file", "fw_spec_field", "multirun",
+        "runs"
     ]
 
     def run_task(self, fw_spec):
@@ -64,13 +61,14 @@ def run_task(self, fw_spec):
         input_file = self.get("input_file", "mol.qin")
         output_file = self.get("output_file", "mol.qout")
         multirun = self.get("multirun", False)
+        runs = self.get("runs", None)
 
         # parse the QChem directory
         logger.info(f"PARSING DIRECTORY: {calc_dir}")
 
         additional_fields = self.get("additional_fields", {})
 
-        drone = QChemDrone(additional_fields=additional_fields)
+        drone = QChemDrone(runs=runs, additional_fields=additional_fields)
 
         # assimilate (i.e., parse)
         task_doc = drone.assimilate(

atomate/qchem/firetasks/run_calc.py

@@ -14,7 +14,7 @@
 
 from atomate.utils.utils import env_chk, get_logger
 
-__author__ = "Samuel Blau"
+__author__ = "Samuel Blau, Evan Spotte-Smith"
 __copyright__ = "Copyright 2018, The Materials Project"
 __version__ = "0.1"
 __maintainer__ = "Samuel Blau"
@@ -91,22 +91,11 @@ class RunQChemCustodian(FiretaskBase):
 
     required_params = ["qchem_cmd"]
     optional_params = [
-        "multimode",
-        "input_file",
-        "output_file",
-        "max_cores",
-        "qclog_file",
-        "suffix",
-        "calc_loc",
-        "save_scratch",
-        "max_errors",
-        "job_type",
-        "handler_group",
-        "gzipped_output",
-        "backup",
-        "linked",
-        "max_iterations",
-        "max_molecule_perturb_scale",
+        "multimode", "input_file", "output_file", "max_cores", "qclog_file",
+        "suffix", "calc_loc", "save_scratch", "max_errors", "job_type",
+        "handler_group", "gzipped_output", "backup", "linked",
+        "max_iterations", "max_molecule_perturb_scale", "freq_before_opt",
+        "transition_state"
     ]
 
     def run_task(self, fw_spec):
@@ -135,6 +124,8 @@ def run_task(self, fw_spec):
         max_molecule_perturb_scale = self.get("max_molecule_perturb_scale", 0.3)
         job_type = self.get("job_type", "normal")
         gzipped_output = self.get("gzipped_output", True)
+        transition_state = self.get("transition_state", False)
+        freq_before_opt = self.get("freq_before_opt", False)
 
         handler_groups = {
             "default": [
@@ -169,6 +160,8 @@ def run_task(self, fw_spec):
                     qclog_file=qclog_file,
                     max_iterations=max_iterations,
                     linked=linked,
+                    freq_before_opt=freq_before_opt,
+                    transition_state=transition_state,
                     save_final_scratch=save_scratch,
                     max_cores=max_cores,
                     calc_loc=calc_loc,
@@ -183,6 +176,8 @@ def run_task(self, fw_spec):
                     max_iterations=max_iterations,
                     max_molecule_perturb_scale=max_molecule_perturb_scale,
                     linked=linked,
+                    freq_before_opt=freq_before_opt,
+                    transition_state=transition_state,
                     save_final_scratch=save_scratch,
                     max_cores=max_cores,
                     calc_loc=calc_loc,

atomate/qchem/firetasks/tests/test_geo_transformations.py

@@ -1,13 +1,14 @@
 import os
 import unittest
 
-import numpy as np
-from pymatgen.core import Molecule
-
-from atomate.qchem.firetasks.geo_transformations import RotateTorsion
+from atomate.qchem.firetasks.geo_transformations import RotateTorsion, PerturbGeometry
 from atomate.utils.testing import AtomateTest
+from pymatgen.core.structure import Molecule
+import numpy as np
+from pymatgen.io.qchem.outputs import QCOutput
+import numpy as np
 
-__author__ = "Brandon Wood"
+__author__ = "Brandon Wood, Evan Spotte-Smith"
 __email__ = "[email protected]"
 
 module_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)))
@@ -48,5 +49,39 @@ def test_rotate_torsion(self):
         )
 
 
+class TestPerturbGeometry(AtomateTest):
+    @classmethod
+    def setUpClass(cls):
+
+        cls.ts_init = Molecule.from_file(
+            os.path.join(module_dir, "..", "..", "test_files",
+                         "ts_init.xyz"))
+        cls.ts_perturbed = Molecule.from_file(
+            os.path.join(module_dir, "..", "..", "test_files",
+                         "ts_perturbed.xyz"))
+        cls.mode = QCOutput(os.path.join(
+            module_dir, "..", "..", "test_files", "ts.out"
+        )).data["frequency_mode_vectors"][0]
+
+    def setUp(self, lpad=False):
+        super(TestPerturbGeometry, self).setUp(lpad=False)
+
+    def tearDown(self):
+        pass
+
+    def test_perturb(self):
+        ft = PerturbGeometry({
+            "molecule": self.ts_init,
+            "mode": self.mode,
+            "scale": 1.0
+        })
+        pert_mol = ft.run_task({})
+        test_mol = Molecule.from_dict(
+            pert_mol.as_dict()["update_spec"]["prev_calc_molecule"])
+        np.testing.assert_equal(self.ts_perturbed.species, test_mol.species)
+        np.testing.assert_allclose(
+            self.ts_perturbed.cart_coords, test_mol.cart_coords, atol=0.0001)
+
+
 if __name__ == "__main__":
     unittest.main()