Merge pull request #765 from rdguha1995/main

Correcting the Proton Energy FW

Author: Zhuoying

atomate/qchem/firetasks/parse_outputs.py

@@ -126,13 +126,13 @@ class ProtCalcToDb(FiretaskBase):
 
     Optional params:
         calc_dir (str): path to dir (on current filesystem) that contains QChem
-            input and output files for both the H0 and H2+. Default: use current working directory.
+            input and output files for both the H2O and H3O+ calculations. Default: use current working directory.
         calc_loc (str OR bool): if True will set most recent calc_loc. If str
             search for the most recent calc_loc with the matching name
-        input_file_H0 (str): name of the QChem input file for H0 calculation
-        output_file_H0 (str): name of the QChem output file for H0 calculation
-        input_file_H2 (str): name of the QChem input file for H2 calculation
-        output_file_H2 (str): name of the QChem output file for H2 calculation
+        input_file_H2O (str): name of the QChem input file for H2O calculation
+        output_file_H2O (str): name of the QChem output file for H2O calculation
+        input_file_H3O (str): name of the QChem input file for H3O+ calculation
+        output_file_H3O (str): name of the QChem output file for H3O+ calculation
         additional_fields (dict): dict of additional fields to add
         db_file (str): path to file containing the database credentials.
             Supports env_chk. Default: write data to JSON file.
@@ -143,10 +143,10 @@ class ProtCalcToDb(FiretaskBase):
     optional_params = [
         "calc_dir",
         "calc_loc",
-        "input_file_H0",
-        "output_file_H0",
-        "input_file_H2",
-        "output_file_H2",
+        "input_file_H2O",
+        "output_file_H2O",
+        "input_file_H3O",
+        "output_file_H3O",
         "additional_fields",
         "db_file",
         "runs",
@@ -159,10 +159,10 @@ def run_task(self, fw_spec):
             calc_dir = self["calc_dir"]
         elif self.get("calc_loc"):
             calc_dir = get_calc_loc(self["calc_loc"], fw_spec["calc_locs"])["path"]
-        input_file_H0 = self.get("input_file_H0", "H0.qin")
-        output_file_H0 = self.get("output_file_H0", "H0.qout")
-        input_file_H2 = self.get("input_file_H2", "H2_plus.qin")
-        output_file_H2 = self.get("output_file_H2", "H2_plus.qout")
+        input_file_H2O = self.get("input_file_H2O", "water.qin")
+        output_file_H2O = self.get("output_file_H2O", "water.qout")
+        input_file_H3O = self.get("input_file_H3O", "hydronium.qin")
+        output_file_H3O = self.get("output_file_H3O", "hydronium.qout")
         runs = self.get("runs", None)
 
         # parse the QChem directory
@@ -175,15 +175,15 @@ def run_task(self, fw_spec):
         # assimilate (i.e., parse)
         task_doc_1 = drone.assimilate(
             path=calc_dir,
-            input_file=input_file_H0,
-            output_file=output_file_H0,
+            input_file=input_file_H2O,
+            output_file=output_file_H2O,
             multirun=False,
         )
 
         task_doc_2 = drone.assimilate(
             path=calc_dir,
-            input_file=input_file_H2,
-            output_file=output_file_H2,
+            input_file=input_file_H3O,
+            output_file=output_file_H3O,
             multirun=False,
         )
 
@@ -195,6 +195,12 @@ def run_task(self, fw_spec):
         task_doc_clean["orig"]["molecule"]["spin_multiplicity"] = 1
         task_doc_clean["output"]["initial_molecule"]["charge"] = 1
         task_doc_clean["output"]["initial_molecule"]["spin_multiplicity"] = 1
+        task_doc_clean["output"]["initial_molecule"]["sites"] = [{'name': 'H', 'species': [{'element': 'H', 'occu': 1}], 'xyz': [0.0, 0.0, 0.0], 'properties': {}}]
+        task_doc_clean["output"]["mulliken"] = [+1.000000]
+        task_doc_clean["output"]["resp"] = [+1.000000]
+        task_doc_clean["output"]["optimized_molecule"]["charge"] = 1
+        task_doc_clean["output"]["optimized_molecule"]["spin_multiplicity"] = 1
+        task_doc_clean["output"]["optimized_molecule"]["sites"] = [{'name': 'H', 'species': [{'element': 'H', 'occu': 1}], 'xyz': [0.0, 0.0, 0.0], 'properties': {}}]
         task_doc_clean["output"]["final_energy"] = (
             task_doc_2["output"]["final_energy"] - task_doc_1["output"]["final_energy"]
         )

atomate/qchem/fireworks/core.py

@@ -117,11 +117,11 @@ def __init__(
         """
         For this custom Firework the electronic energy of a proton in a specific solvent environment is approximated.
         Since a proton has 0 electrons,running a QChem job would yield an error. The energy can be approximated by
-        calculating the electronic energy of a proton and a hydrogen atom at infinite separation and then subtracting
-        the electronic energy of a hydrogen atom. This Firework combines these two calculations and adds a task doc to the
+        calculating the electronic energy of a hydronium ion and a water molecule and then subtracting
+        the respective electronic energies. This Firework combines these two calculations and adds a task doc to the
         DB with the separate calculation details and the effective energy after subtraction.
 
-        Args:
+        Arg
             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.
@@ -150,21 +150,31 @@ def __init__(
 
         qchem_input_params = qchem_input_params or {}
 
-        H_site = Site("H", [0.0, 0.0, 0.0])
-        H_site_inf = Site("H", [100000.0, 0.0, 0.0])
-        H0_atom = Molecule.from_sites([H_site])
-        H2_plus_mol = Molecule.from_sites([H_site, H_site_inf])
-        H0_atom.set_charge_and_spin(0, 2)
-        H2_plus_mol.set_charge_and_spin(1, 2)
-        input_file_1 = "H0.qin"
-        output_file_1 = "H0.qout"
-        input_file_2 = "H2_plus.qin"
-        output_file_2 = "H2_plus.qout"
+        H_site_1_H2O = Site("H", [0.18338, 2.20176, 0.01351])
+        H_site_2_H2O = Site("H", [-1.09531, 1.61602, 0.70231])
+        O_site_H2O = Site("O", [-0.80595, 2.22952, -0.01914])
+        H2O_molecule = Molecule.from_sites([H_site_1_H2O, H_site_2_H2O, O_site_H2O])
+
+        H_site_1_H3O = Site("H", [0.11550, 2.34733, 0.00157])
+        H_site_2_H3O = Site("H", [-1.17463, 1.77063, 0.67652])
+        H_site_3_H3O = Site("H", [-1.29839, 2.78012, -0.51436])
+        O_site_H3O = Site("O", [-0.78481, 1.99137, -0.20661])
+        H3O_ion = Molecule.from_sites(
+            [H_site_1_H3O, H_site_2_H3O, H_site_3_H3O, O_site_H3O]
+        )
+
+        H2O_molecule.set_charge_and_spin(0, 1)
+        H3O_ion.set_charge_and_spin(1, 1)
+
+        input_file_1 = "water.qin"
+        output_file_1 = "water.qout"
+        input_file_2 = "hydronium.qin"
+        output_file_2 = "hydronium.qout"
         t = []
         t.append(
             WriteInputFromIOSet(
-                molecule=H0_atom,
-                qchem_input_set="SinglePointSet",
+                molecule=H2O_molecule,
+                qchem_input_set="OptSet",
                 input_file=input_file_1,
                 qchem_input_params=qchem_input_params,
             )
@@ -184,8 +194,8 @@ def __init__(
 
         t.append(
             WriteInputFromIOSet(
-                molecule=H2_plus_mol,
-                qchem_input_set="SinglePointSet",
+                molecule=H3O_ion,
+                qchem_input_set="OptSet",
                 input_file=input_file_2,
                 qchem_input_params=qchem_input_params,
             )
@@ -205,10 +215,10 @@ def __init__(
         t.append(
             ProtCalcToDb(
                 db_file=db_file,
-                input_file_H0=input_file_1,
-                output_file_H0=output_file_1,
-                input_file_H2=input_file_2,
-                output_file_H2=output_file_2,
+                input_file_H2O=input_file_1,
+                output_file_H2O=output_file_1,
+                input_file_H3O=input_file_2,
+                output_file_H3O=output_file_2,
                 additional_fields={"task_label": name},
             )
         )

atomate/qchem/fireworks/tests/test_core.py

@@ -137,19 +137,30 @@ def test_SinglePointFW_not_defaults(self):
         self.assertEqual(firework.name, "special single point")
 
     def test_ProtonEnergyFW(self):
-        H_site = Site("H", [0.0, 0.0, 0.0])
-        H_site_inf = Site("H", [100000.0, 0.0, 0.0])
-        H0_atom = Molecule.from_sites([H_site])
-        H2_plus_mol = Molecule.from_sites([H_site, H_site_inf])
-        H0_atom.set_charge_and_spin(0, 2)
-        H2_plus_mol.set_charge_and_spin(1, 2)
+
+        H_site_1_H2O = Site("H", [0.18338, 2.20176, 0.01351])
+        H_site_2_H2O = Site("H", [-1.09531, 1.61602, 0.70231])
+        O_site_H2O = Site("O", [-0.80595, 2.22952, -0.01914])
+        H2O_molecule = Molecule.from_sites([H_site_1_H2O, H_site_2_H2O, O_site_H2O])
+
+        H_site_1_H3O = Site("H", [0.11550, 2.34733, 0.00157])
+        H_site_2_H3O = Site("H", [-1.17463, 1.77063, 0.67652])
+        H_site_3_H3O = Site("H", [-1.29839, 2.78012, -0.51436])
+        O_site_H3O = Site("O", [-0.78481, 1.99137, -0.20661])
+        H3O_ion = Molecule.from_sites(
+            [H_site_1_H3O, H_site_2_H3O, H_site_3_H3O, O_site_H3O]
+        )
+
+        H2O_molecule.set_charge_and_spin(0, 1)
+        H3O_ion.set_charge_and_spin(1, 1)
+
         firework = ProtonEnergyFW(qchem_input_params={"smd_solvent": "water"})
         self.assertEqual(
             firework.tasks[0].as_dict(),
             WriteInputFromIOSet(
-                molecule=H0_atom,
-                qchem_input_set="SinglePointSet",
-                input_file="H0.qin",
+                molecule=H2O_molecule,
+                qchem_input_set="OptSet",
+                input_file="water.qin",
                 qchem_input_params={"smd_solvent": "water"},
             ).as_dict(),
         )
@@ -158,8 +169,8 @@ def test_ProtonEnergyFW(self):
             RunQChemCustodian(
                 qchem_cmd=">>qchem_cmd<<",
                 multimode=">>multimode<<",
-                input_file="H0.qin",
-                output_file="H0.qout",
+                input_file="water.qin",
+                output_file="water.qout",
                 max_cores=">>max_cores<<",
                 max_errors=5,
                 job_type="normal",
@@ -169,9 +180,9 @@ def test_ProtonEnergyFW(self):
         self.assertEqual(
             firework.tasks[2].as_dict(),
             WriteInputFromIOSet(
-                molecule=H2_plus_mol,
-                qchem_input_set="SinglePointSet",
-                input_file="H2_plus.qin",
+                molecule=H3O_ion,
+                qchem_input_set="OptSet",
+                input_file="hydronium.qin",
                 qchem_input_params={"smd_solvent": "water"},
             ).as_dict(),
         )
@@ -180,8 +191,8 @@ def test_ProtonEnergyFW(self):
             RunQChemCustodian(
                 qchem_cmd=">>qchem_cmd<<",
                 multimode=">>multimode<<",
-                input_file="H2_plus.qin",
-                output_file="H2_plus.qout",
+                input_file="hydronium.qin",
+                output_file="hydronium.qout",
                 max_cores=">>max_cores<<",
                 max_errors=5,
                 job_type="normal",
@@ -192,10 +203,10 @@ def test_ProtonEnergyFW(self):
             firework.tasks[4].as_dict(),
             ProtCalcToDb(
                 db_file=None,
-                input_file_H0="H0.qin",
-                output_file_H0="H0.qout",
-                input_file_H2="H2_plus.qin",
-                output_file_H2="H2_plus.qout",
+                input_file_H2O="water.qin",
+                output_file_H2O="water.qout",
+                input_file_H3O="hydronium.qin",
+                output_file_H3O="hydronium.qout",
                 additional_fields={"task_label": "proton electronic energy"},
             ).as_dict(),
         )