From 3b09d18112f089231b4d0fc3691555de41da3471 Mon Sep 17 00:00:00 2001 From: Evan Walter Clark Spotte-Smith Date: Wed, 4 Sep 2019 15:45:56 -0400 Subject: [PATCH 1/4] Added FrequencyFlatteningTransitionStatteFW --- atomate/qchem/fireworks/core.py | 120 ++++++++++++++++++++++++++++++-- 1 file changed, 115 insertions(+), 5 deletions(-) diff --git a/atomate/qchem/fireworks/core.py b/atomate/qchem/fireworks/core.py index c68c7b4f9..d3bf8d23d 100644 --- a/atomate/qchem/fireworks/core.py +++ b/atomate/qchem/fireworks/core.py @@ -226,9 +226,10 @@ def __init__(self, """ qchem_input_params = qchem_input_params or {} - input_file="mol.qin" - output_file="mol.qout" - t = [] + input_file = "mol.qin" + output_file = "mol.qout" + + t = list() t.append( WriteInputFromIOSet( molecule=molecule, @@ -262,6 +263,117 @@ def __init__(self, **kwargs) +#TODO: test (actually with Q-Chem and by writing tests) +class FrequencyFlatteningTransitionStateFW(Firework): + def __init__(self, + molecule=None, + name="frequency flattening transition state optimization", + qchem_cmd=">>qchem_cmd<<", + multimode=">>multimode<<", + max_cores=">>max_cores<<", + directory=None, + qchem_input_params=None, + max_iterations=10, + max_molecule_perturb_scale=0.3, + reversed_direction=False, + db_file=None, + parents=None, + **kwargs): + """ + First, perform a search over the potential energy surface between reactants and products in order to determine + a guess for the transition state. Then, iteratively optimize the transition state structure and flatten + imaginary frequencies to ensure that the resulting structure is a true transition state. + + Args: + molecule (dict of Molecules): Input molecules. The dict should have two entries, "reactants" and "products". + Note that the order of the molecules, and the atoms of the molecules, is important. The FSM and GSM + methods will only work properly if the atoms in the reactants and products correspond exactly. + 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. + directory (str): Location where calculation should take place. Default is current working directory. + 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 conjuction with more typical modifications, + as seen in the test_double_FF_opt workflow test. + max_iterations (int): Number of perturbation -> optimization -> frequency + iterations to perform. Defaults to 10. + max_molecule_perturb_scale (float): The maximum scaled perturbation that can be + applied to the molecule. Defaults to 0.3. + reversed_direction (bool): Whether to reverse the direction of the vibrational + frequency vectors. Defaults to False. + 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 = qchem_input_params or {} + input_file = "mol.qin" + output_file = "mol.qout" + + t = list() + t.append( + WriteInputFromIOSet( + molecule=molecule, + qchem_input_set="FreezingStringSet", + 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, + job_type="ts_with_frequency_flattener", + max_iterations=max_iterations, + max_molecule_perturb_scale=max_molecule_perturb_scale, + reversed_direction=reversed_direction, + gzipped_output=False)) + if directory is None: + t.append( + QChemToDb( + db_file=db_file, + input_file=input_file, + output_file=output_file, + additional_fields={ + "task_label": name, + "special_run_type": "ts_frequency_flattener" + })) + else: + t.append( + QChemToDb( + db_file=db_file, + calc_dir=directory, + input_file="mol.qin", + output_file="mol.qout", + additional_fields={ + "task_label": name, + "special_run_type": "ts_frequency_flattener" + })) + + super(FrequencyFlatteningTransitionStateFW, self).__init__( + t, + parents=parents, + name=name, + **kwargs) + + class FragmentFW(Firework): def __init__(self, molecule=None, @@ -270,8 +382,6 @@ def __init__(self, additional_charges=None, do_triplets=True, name="fragment and optimize", - qchem_cmd=">>qchem_cmd<<", - multimode=">>multimode<<", max_cores=">>max_cores<<", qchem_input_params=None, db_file=None, From 99d59e2449d9787817274e9f4ddfdcbce557e25a Mon Sep 17 00:00:00 2001 From: Evan Walter Clark Spotte-Smith Date: Thu, 30 Jan 2020 10:39:07 -0800 Subject: [PATCH 2/4] Get up-to-date with recent changes in pymatgen --- atomate/qchem/firetasks/fragmenter.py | 6 ++---- atomate/qchem/firetasks/write_inputs.py | 20 ++++++-------------- atomate/qchem/tests/test_drones.py | 16 ++++------------ 3 files changed, 12 insertions(+), 30 deletions(-) diff --git a/atomate/qchem/firetasks/fragmenter.py b/atomate/qchem/firetasks/fragmenter.py index 8edc8ae63..084cedfb0 100644 --- a/atomate/qchem/firetasks/fragmenter.py +++ b/atomate/qchem/firetasks/fragmenter.py @@ -208,13 +208,11 @@ def _in_database(self, molecule): # otherwise, look through the docs for an entry with an isomorphic molecule with # equivalent charge and multiplicity else: - new_mol_graph = MoleculeGraph.with_local_env_strategy(molecule, OpenBabelNN(), - reorder=False, extend_structure=False) + new_mol_graph = MoleculeGraph.with_local_env_strategy(molecule, OpenBabelNN()) for doc in self.all_relevant_docs: if molecule.composition.reduced_formula == doc["formula_pretty"]: old_mol = Molecule.from_dict(doc["input"]["initial_molecule"]) - old_mol_graph = MoleculeGraph.with_local_env_strategy(old_mol, OpenBabelNN(), - reorder=False, extend_structure=False) + old_mol_graph = MoleculeGraph.with_local_env_strategy(old_mol, OpenBabelNN()) # If such an equivalent molecule is found, return true if new_mol_graph.isomorphic_to(old_mol_graph) and molecule.charge == old_mol_graph.molecule.charge and molecule.spin_multiplicity == old_mol_graph.molecule.spin_multiplicity: return True diff --git a/atomate/qchem/firetasks/write_inputs.py b/atomate/qchem/firetasks/write_inputs.py index e196722f1..3cf598525 100644 --- a/atomate/qchem/firetasks/write_inputs.py +++ b/atomate/qchem/firetasks/write_inputs.py @@ -77,13 +77,9 @@ def run_task(self, fw_spec): mol = self.get("molecule") # check if mol and prev_calc_mol are isomorphic mol_graph = MoleculeGraph.with_local_env_strategy(mol, - OpenBabelNN(), - reorder=False, - extend_structure=False) - prev_mol_graph = MoleculeGraph.with_local_env_strategy(prev_calc_molecule, - OpenBabelNN(), - reorder=False, - extend_structure=False) + OpenBabelNN()) + prev_mol_graph = MoleculeGraph.with_local_env_strategy(prev_calc_mol, + OpenBabelNN()) # If they are isomorphic, aka a previous FW has not changed bonding, # then we will use prev_calc_mol. If bonding has changed, we will use mol. if mol_graph.isomorphic_to(prev_mol_graph): @@ -148,13 +144,9 @@ def run_task(self, fw_spec): mol = self.get("molecule") # check if mol and prev_calc_mol are isomorphic mol_graph = MoleculeGraph.with_local_env_strategy(mol, - OpenBabelNN(), - reorder=False, - extend_structure=False) - prev_mol_graph = MoleculeGraph.with_local_env_strategy(prev_calc_molecule, - OpenBabelNN(), - reorder=False, - extend_structure=False) + OpenBabelNN()) + prev_mol_graph = MoleculeGraph.with_local_env_strategy(prev_calc_mol, + OpenBabelNN()) if mol_graph.isomorphic_to(prev_mol_graph): mol = prev_calc_mol else: diff --git a/atomate/qchem/tests/test_drones.py b/atomate/qchem/tests/test_drones.py index 970f53b53..252f05b1f 100644 --- a/atomate/qchem/tests/test_drones.py +++ b/atomate/qchem/tests/test_drones.py @@ -208,13 +208,9 @@ def test_assimilate_unstable_opt(self): self.assertEqual(doc["orig"]["rem"], doc["calcs_reversed"][-1]["input"]["rem"]) self.assertEqual(doc["orig"]["molecule"], doc["calcs_reversed"][-1]["input"]["molecule"]) orig_molgraph = MoleculeGraph.with_local_env_strategy(Molecule.from_dict(doc["orig"]["molecule"]), - OpenBabelNN(), - reorder=False, - extend_structure=False) + OpenBabelNN()) initial_molgraph = MoleculeGraph.with_local_env_strategy(Molecule.from_dict(doc["input"]["initial_molecule"]), - OpenBabelNN(), - reorder=False, - extend_structure=False) + OpenBabelNN()) self.assertEqual(orig_molgraph.isomorphic_to(initial_molgraph), True) def test_assimilate_opt_with_hidden_changes_from_handler(self): @@ -238,13 +234,9 @@ def test_assimilate_opt_with_hidden_changes_from_handler(self): self.assertEqual(doc["pointgroup"], "C1") self.assertEqual(doc["orig"]["rem"], doc["calcs_reversed"][-1]["input"]["rem"]) orig_molgraph = MoleculeGraph.with_local_env_strategy(Molecule.from_dict(doc["orig"]["molecule"]), - OpenBabelNN(), - reorder=False, - extend_structure=False) + OpenBabelNN()) initial_molgraph = MoleculeGraph.with_local_env_strategy(Molecule.from_dict(doc["input"]["initial_molecule"]), - OpenBabelNN(), - reorder=False, - extend_structure=False) + OpenBabelNN()) self.assertEqual(orig_molgraph.isomorphic_to(initial_molgraph), False) def test_assimilate_disconnected_opt(self): From 68fd39fdd071d02bab72e115e48dc7729e46e52d Mon Sep 17 00:00:00 2001 From: Evan Walter Clark Spotte-Smith Date: Mon, 4 Dec 2023 13:58:48 -0800 Subject: [PATCH 3/4] Add hashes to drone --- atomate/qchem/drones.py | 49 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 49 insertions(+) diff --git a/atomate/qchem/drones.py b/atomate/qchem/drones.py index 08c416b4b..28e54d06b 100644 --- a/atomate/qchem/drones.py +++ b/atomate/qchem/drones.py @@ -11,10 +11,14 @@ from monty.json import jsanitize from pymatgen.apps.borg.hive import AbstractDrone from pymatgen.core import Molecule +from pymatgen.core.periodic_table import Element from pymatgen.io.babel import BabelMolAdaptor from pymatgen.io.qchem.inputs import QCInput from pymatgen.io.qchem.outputs import QCOutput +from pymatgen.analysis.graphs import MoleculeGraph +from pymatgen.analysis.local_env import OpenBabelNN, metal_edge_extender from pymatgen.symmetry.analyzer import PointGroupAnalyzer +from pymatgen.util.graph_hashing import weisfeiler_lehman_graph_hash from atomate import __version__ as atomate_version from atomate.utils.utils import get_logger @@ -31,6 +35,13 @@ logger = get_logger(__name__) +METALS = { + str(e) + for e in [Element.from_Z(i) for i in range(1, 119)] + if e.is_metal +} + + class QChemDrone(AbstractDrone): """ A QChem drone to parse QChem calculations and insert an organized, searchable entry into the database. @@ -53,6 +64,9 @@ class QChemDrone(AbstractDrone): "chemsys", "pointgroup", "formula_alphabetical", + "species_hash", + "coord_hash", + "species_hash_nometal" }, "input": {"initial_molecule", "job_type"}, "output": {"initial_molecule", "job_type", "final_energy"}, @@ -346,6 +360,41 @@ def generate_doc(self, dir_name, qcinput_files, qcoutput_files, multirun): smiles = pbmol.write("smi").split()[0] d["smiles"] = smiles + # Add graph hashes + # This is primarily for emmet builders + if "optimized_molecule" in d["output"]: + hash_mol = d["output"]["optimized_molecule"] + else: + hash_mol = d["output"]["initial_molecule"] + + hash_mg = MoleculeGraph.with_local_env_strategy(hash_mol, OpenBabelNN()) + hash_mg = metal_edge_extender(hash_mg) + undir_mg = hash_mg.graph.to_undirected() + + metal_inds = [i for i, e in enumerate(hash_mol.species) if str(e) in METALS] + + to_delete = list() + for bond in hash_mg.graph.edges(): + if bond[0] in metal_inds or bond[1] in metal_inds: + to_delete.append((bond[0], bond[1])) + + mg_nometal = copy.deepcopy(hash_mg) + for b in to_delete: + mg_nometal.break_edge(b[0], b[1], allow_reverse=True) + + d["coord_hash"] = weisfeiler_lehman_graph_hash( + undir_mg, + node_attr="coords" + ) + d["species_hash"] = weisfeiler_lehman_graph_hash( + undir_mg, + node_attr="specie" + ) + d["species_hash_nometal"] = weisfeiler_lehman_graph_hash( + mg_nometal.graph.to_undirected(), + node_attr="specie" + ) + d["state"] = "successful" if d_calc_final["completion"] else "unsuccessful" if "special_run_type" in d: if d["special_run_type"] in [ From adffd135ac3a2cd9b5fe235ee06c34a6aa4c1d15 Mon Sep 17 00:00:00 2001 From: Evan Walter Clark Spotte-Smith Date: Mon, 4 Dec 2023 14:20:46 -0800 Subject: [PATCH 4/4] Tests added; tests pass --- atomate/qchem/tests/test_drones.py | 33 ++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/atomate/qchem/tests/test_drones.py b/atomate/qchem/tests/test_drones.py index 64c8aff6d..464de18fe 100644 --- a/atomate/qchem/tests/test_drones.py +++ b/atomate/qchem/tests/test_drones.py @@ -41,6 +41,9 @@ def test_assimilate_opt(self): self.assertEqual(doc["smiles"], "O1[C](O[Li])OC=C1") self.assertEqual(doc["formula_pretty"], "LiH2(CO)3") self.assertEqual(doc["formula_anonymous"], "AB2C3D3") + self.assertEqual(doc["species_hash"], "44fb6a8c8e99aed63f23e573e720018d") + self.assertEqual(doc["coord_hash"], "acbb408648e992fea44acb87e912fd5f") + self.assertEqual(doc["species_hash_nometal"], "48ba8b7456a39b5ee1f8d76772d9f4c8") self.assertEqual(doc["chemsys"], "C-H-Li-O") self.assertEqual(doc["pointgroup"], "Cs") self.assertIn("custodian", doc) @@ -77,6 +80,9 @@ def test_assimilate_pes_scan(self): self.assertEqual(doc["smiles"], "[O]C(=O)[O]") self.assertEqual(doc["formula_pretty"], "CO3") self.assertEqual(doc["formula_anonymous"], "AB3") + self.assertEqual(doc["species_hash"], "75e7a4125709cb5a14d1ce2b84c3cdbd") + self.assertEqual(doc["coord_hash"], "65b93a5088773337b9372c4ce65aeb37") + self.assertEqual(doc["species_hash_nometal"], "75e7a4125709cb5a14d1ce2b84c3cdbd") self.assertEqual(doc["chemsys"], "C-O") self.assertEqual(doc["pointgroup"], "C2v") self.assertIn("custodian", doc) @@ -130,6 +136,9 @@ def test_assimilate_freq(self): self.assertEqual(doc["smiles"], "O1[C](O[Li])OC=C1") self.assertEqual(doc["formula_pretty"], "LiH2(CO)3") self.assertEqual(doc["formula_anonymous"], "AB2C3D3") + self.assertEqual(doc["species_hash"], "44fb6a8c8e99aed63f23e573e720018d") + self.assertEqual(doc["coord_hash"], "acbb408648e992fea44acb87e912fd5f") + self.assertEqual(doc["species_hash_nometal"], "48ba8b7456a39b5ee1f8d76772d9f4c8") self.assertEqual(doc["chemsys"], "C-H-Li-O") self.assertEqual(doc["pointgroup"], "Cs") self.assertIn("custodian", doc) @@ -203,6 +212,9 @@ def test_assimilate_FF(self): self.assertEqual(doc["smiles"], "O1[C](O[Li])OC=C1") self.assertEqual(doc["formula_pretty"], "LiH2(CO)3") self.assertEqual(doc["formula_anonymous"], "AB2C3D3") + self.assertEqual(doc["species_hash"], "44fb6a8c8e99aed63f23e573e720018d") + self.assertEqual(doc["coord_hash"], "acbb408648e992fea44acb87e912fd5f") + self.assertEqual(doc["species_hash_nometal"], "48ba8b7456a39b5ee1f8d76772d9f4c8") self.assertEqual(doc["chemsys"], "C-H-Li-O") self.assertEqual(doc["pointgroup"], "Cs") self.assertIn("custodian", doc) @@ -313,6 +325,9 @@ def test_assimilate_ffts(self): self.assertEqual(doc["smiles"], "O(C(=O)[O])[Li].[CH2]COC(=O)O[Li]") self.assertEqual(doc["formula_pretty"], "LiH2C2O3") self.assertEqual(doc["formula_anonymous"], "AB2C2D3") + self.assertEqual(doc["species_hash"], "b58892da682cac0193cf85f25fe8c25b") + self.assertEqual(doc["coord_hash"], "ba40774a9d7a39f8354d0ca7efaff6d0") + self.assertEqual(doc["species_hash_nometal"], "d7ab8a26c0d207bad6bd1316a368c8d5") self.assertEqual(doc["chemsys"], "C-H-Li-O") self.assertEqual(doc["pointgroup"], "C1") self.assertIn("custodian", doc) @@ -396,6 +411,9 @@ def test_multirun(self): self.assertEqual(doc["smiles"], "O1[C](O[Li])OC=C1") self.assertEqual(doc["formula_pretty"], "LiH2(CO)3") self.assertEqual(doc["formula_anonymous"], "AB2C3D3") + self.assertEqual(doc["species_hash"], "44fb6a8c8e99aed63f23e573e720018d") + self.assertEqual(doc["coord_hash"], "d11f6abeef573141250f38af1388ca0c") + self.assertEqual(doc["species_hash_nometal"], "48ba8b7456a39b5ee1f8d76772d9f4c8") self.assertEqual(doc["chemsys"], "C-H-Li-O") self.assertEqual(doc["pointgroup"], "C2") self.assertIn("calcs_reversed", doc) @@ -437,6 +455,9 @@ def test_assimilate_unstable_opt(self): self.assertEqual(doc["cputime"], None) self.assertEqual(doc["formula_pretty"], "CS2NO") self.assertEqual(doc["formula_anonymous"], "ABCD2") + self.assertEqual(doc["species_hash"], "1559ce7584cf8c27f1c6044a6af76dd1") + self.assertEqual(doc["coord_hash"], "8698b987cdb70eed57bd0a7e77b7e00c") + self.assertEqual(doc["species_hash_nometal"], "1559ce7584cf8c27f1c6044a6af76dd1") self.assertEqual(doc["chemsys"], "C-N-O-S") self.assertEqual(doc["pointgroup"], "C1") self.assertEqual(doc["orig"]["rem"], doc["calcs_reversed"][-1]["input"]["rem"]) @@ -471,6 +492,9 @@ def test_assimilate_opt_with_hidden_changes_from_handler(self): self.assertEqual(doc["cputime"], 7471.17) self.assertEqual(doc["formula_pretty"], "HC2O") self.assertEqual(doc["formula_anonymous"], "ABC2") + self.assertEqual(doc["species_hash"], "6dc4aca792bcd6bd45bc5176f42f6aee") + self.assertEqual(doc["coord_hash"], "7cd547f71ddf74efcfb00743161d07f2") + self.assertEqual(doc["species_hash_nometal"], "6dc4aca792bcd6bd45bc5176f42f6aee") self.assertEqual(doc["chemsys"], "C-H-O") self.assertEqual(doc["pointgroup"], "C1") self.assertEqual(doc["orig"]["rem"], doc["calcs_reversed"][-1]["input"]["rem"]) @@ -504,6 +528,9 @@ def test_assimilate_disconnected_opt(self): self.assertEqual(doc["cputime"], 8825.76) self.assertEqual(doc["formula_pretty"], "H2C2O3") self.assertEqual(doc["formula_anonymous"], "A2B2C3") + self.assertEqual(doc["species_hash"], "c87c6b5a4bb8632cdb934e400a0237fb") + self.assertEqual(doc["coord_hash"], "6adeeb1d55585a35a6bcf9e1513218f2") + self.assertEqual(doc["species_hash_nometal"], "c87c6b5a4bb8632cdb934e400a0237fb") self.assertEqual(doc["chemsys"], "C-H-O") self.assertEqual(doc["pointgroup"], "C1") self.assertEqual(doc["orig"]["rem"], doc["calcs_reversed"][-1]["input"]["rem"]) @@ -527,6 +554,9 @@ def test_assimilate_sp(self): self.assertEqual(doc["smiles"], "[O]") self.assertEqual(doc["formula_pretty"], "O2") self.assertEqual(doc["formula_anonymous"], "A") + self.assertEqual(doc["species_hash"], "d6bfd1bfa289860fea9c71dd64fab914") + self.assertEqual(doc["coord_hash"], "d1e604b6971a7e2d889e3172456da7db") + self.assertEqual(doc["species_hash_nometal"], "d6bfd1bfa289860fea9c71dd64fab914") self.assertEqual(doc["chemsys"], "O") self.assertEqual(doc["pointgroup"], "Kh") self.assertIn("custodian", doc) @@ -572,6 +602,9 @@ def test_sp_with_orig(self): self.assertEqual(doc["smiles"], "[O]") self.assertEqual(doc["formula_pretty"], "O2") self.assertEqual(doc["formula_anonymous"], "A") + self.assertEqual(doc["species_hash"], "d6bfd1bfa289860fea9c71dd64fab914") + self.assertEqual(doc["coord_hash"], "4b44e5b5c47ec269779254ae49ca0b51") + self.assertEqual(doc["species_hash_nometal"], "d6bfd1bfa289860fea9c71dd64fab914") self.assertEqual(doc["chemsys"], "O") self.assertEqual(doc["pointgroup"], "Kh") self.assertIn("custodian", doc)