Repository: materialsproject/custodian
Pull request #163 (2020-12-23 20:20:19)

Qchem updates: improved handling of scratch files and scratch directories
## Summary

Resolves previous problems with scratch directories and saving scratch files. Now explicitly sets both QCSCRATCH and QCLOCALSCR correctly and employs a better procedure for saving scratch files of interest. Also standardizes scratch directory naming and placement to avoid unnecessary confusion. 


Patch:
diff --git a/custodian/qchem/__init__.py b/custodian/qchem/__init__.py
index 148cd9d2..69c63b65 100644
--- a/custodian/qchem/__init__.py
+++ b/custodian/qchem/__init__.py
@@ -1,3 +1,9 @@
+"""
+This package implements various QChem Jobs and Error Handlers.
+"""
+
+
 from __future__ import unicode_literals
 
-__author__ = "xiaohuiqu"
+__author__ = "Samuel Blau, Brandon Wood, Shyam Dwaraknath"
+__credits__ = "Xiaohui Qu"
diff --git a/custodian/qchem/handlers.py b/custodian/qchem/handlers.py
index e583d508..81c848ea 100644
--- a/custodian/qchem/handlers.py
+++ b/custodian/qchem/handlers.py
@@ -54,7 +54,9 @@ class QChemErrorHandler(ErrorHandler):
         self.opt_error_history = []
 
     def check(self):
-        # Checks output file for errors.
+        """
+        Checks output file for errors
+        """
         self.outdata = QCOutput(self.output_file).data
         self.errors = self.outdata.get("errors")
         self.warnings = self.outdata.get("warnings")
@@ -70,6 +72,9 @@ class QChemErrorHandler(ErrorHandler):
         return len(self.errors) > 0
 
     def correct(self):
+        """
+        Perform corrections
+        """
         backup({self.input_file, self.output_file})
         actions = []
         self.qcinp = QCInput.from_file(self.input_file)
@@ -139,7 +144,7 @@ class QChemErrorHandler(ErrorHandler):
 
         elif "unable_to_determine_lamda" in self.errors:
             # Set last geom as new starting geom and rerun. If no opt cycles,
-            # use diff SCF strat? Diff initial guess? Change basis?
+            # use diff SCF strat? Diff initial guess? Change basis? Unclear.
             if len(self.outdata.get("energy_trajectory")) > 1:
                 self.qcinp.molecule = self.outdata.get("molecule_from_last_geometry")
                 actions.append({"molecule": "molecule_from_last_geometry"})
@@ -266,11 +271,16 @@ class QChemSCFErrorHandler(ErrorHandler):
         self.errors = None
 
     def check(self):
-        # Checks output file for errors.
+        """
+        Checks output file for errors
+        """
         self.outdata = QCOutput(self.output_file).data
         self.errors = self.outdata.get("errors")
         return len(self.errors) > 0
 
     def correct(self):
+        """
+        Corrects errors, but it hasn't been implemented yet
+        """
         print("This hasn't been implemented yet!")
         return {"errors": self.errors, "actions": None}
diff --git a/custodian/qchem/jobs.py b/custodian/qchem/jobs.py
index 3481aada..121af765 100644
--- a/custodian/qchem/jobs.py
+++ b/custodian/qchem/jobs.py
@@ -1,11 +1,14 @@
 # coding: utf-8
 
+"""
+This module implements basic kinds of jobs for QChem runs.
+"""
+
 import math
 import os
 import shutil
 import copy
 import subprocess
-import random
 import numpy as np
 from pymatgen.core import Molecule
 from pymatgen.io.qchem.inputs import QCInput
@@ -36,10 +39,9 @@ class QCJob(Job):
         output_file="mol.qout",
         qclog_file="mol.qclog",
         suffix="",
-        scratch_dir=os.getcwd(),
+        calc_loc=None,
         save_scratch=False,
-        save_name="saved_scratch",
-        backup=True,
+        backup=True
     ):
         """
         Args:
@@ -51,11 +53,10 @@ class QCJob(Job):
             qclog_file (str): Name of the file to redirect the standard output
                 to. None means not to record the standard output.
             suffix (str): String to append to the file in postprocess.
-            scratch_dir (str): QCSCRATCH directory. Defaults to current directory.
-            save_scratch (bool): Whether to save basic scratch directory contents.
+            calc_loc (str): Path where Q-Chem should run. Defaults to None, in
+                which case Q-Chem will run in the system-defined QCLOCALSCR.
+            save_scratch (bool): Whether to save full scratch directory contents.
                 Defaults to False.
-            save_name (str): Name of the saved scratch directory. Defaults to
-                to "saved_scratch".
             backup (bool): Whether to backup the initial input file. If True, the
                 input will be copied with a ".orig" appended. Defaults to True.
         """
@@ -66,13 +67,15 @@ class QCJob(Job):
         self.max_cores = max_cores
         self.qclog_file = qclog_file
         self.suffix = suffix
-        self.scratch_dir = scratch_dir
+        self.calc_loc = calc_loc
         self.save_scratch = save_scratch
-        self.save_name = save_name
         self.backup = backup
 
     @property
     def current_command(self):
+        """
+        The command to run QChem
+        """
         multi = {"openmp": "-nt", "mpi": "-np"}
         if self.multimode not in multi:
             raise RuntimeError("ERROR: Multimode should only be set to openmp or mpi")
@@ -81,9 +84,8 @@ class QCJob(Job):
             str(self.max_cores),
             self.input_file,
             self.output_file,
+            "scratch"
         ]
-        if self.save_scratch:
-            command.append(self.save_name)
         command = self.qchem_command + command
         com_str = ""
         for part in command:
@@ -91,18 +93,36 @@ class QCJob(Job):
         return com_str
 
     def setup(self):
+        """
+        Sets up environment variables necessary to efficiently run QChem
+        """
         if self.backup:
             shutil.copy(self.input_file, "{}.orig".format(self.input_file))
-        os.environ["QCSCRATCH"] = self.scratch_dir
-        if self.multimode == "openmp":
-            os.environ["QCTHREADS"] = str(self.max_cores)
-            os.environ["OMP_NUM_THREADS"] = str(self.max_cores)
+        if self.multimode == 'openmp':
+            os.environ['QCTHREADS'] = str(self.max_cores)
+            os.environ['OMP_NUM_THREADS'] = str(self.max_cores)
+        os.environ["QCSCRATCH"] = os.getcwd()
+        if self.calc_loc is not None:
+            os.environ["QCLOCALSCR"] = self.calc_loc
 
     def postprocess(self):
+        """
+        Renames and removes scratch files after running QChem
+        """
+        scratch_dir = os.path.join(os.environ["QCSCRATCH"], "scratch")
+        for file in ["HESS", "GRAD", "plots/dens.0.cube"]:
+            file_path = os.path.join(scratch_dir, file)
+            if os.path.exists(file_path):
+                shutil.copy(file_path, os.getcwd())
         if self.suffix != "":
             shutil.move(self.input_file, self.input_file + self.suffix)
             shutil.move(self.output_file, self.output_file + self.suffix)
             shutil.move(self.qclog_file, self.qclog_file + self.suffix)
+            for file in ["HESS", "GRAD", "dens.0.cube"]:
+                if os.path.exists(file):
+                    shutil.move(file, file + self.suffix)
+        if not self.save_scratch:
+            shutil.rmtree(scratch_dir)
 
     def run(self):
         """
@@ -111,11 +131,10 @@ class QCJob(Job):
         Returns:
             (subprocess.Popen) Used for monitoring.
         """
-        myrand = str(random.randint(1, 1000000000))
-        mydir = os.path.join("/tmp", "qchem" + myrand)
-        os.mkdir(mydir)
-        os.environ["QCLOCALSCR"] = mydir
-        qclog = open(self.qclog_file, "w")
+        local_scratch = os.path.join(os.environ["QCLOCALSCR"], "scratch")
+        if os.path.exists(local_scratch):
+            shutil.rmtree(local_scratch)
+        qclog = open(self.qclog_file, 'w')
         p = subprocess.Popen(self.current_command, stdout=qclog, shell=True)
         return p
 
@@ -131,6 +150,7 @@ class QCJob(Job):
         max_molecule_perturb_scale=0.3,
         check_connectivity=True,
         linked=True,
+        save_final_scratch=False,
         **QCJob_kwargs
     ):
         """
@@ -150,6 +170,9 @@ class QCJob(Job):
                 can be applied to the molecule. Defaults to 0.3.
             check_connectivity (bool): Whether to check differences in connectivity
                 introduced by structural perturbation. Defaults to True.
+            linked (bool): Whether or not to use the linked flattener. Defaults to True.
+            save_final_scratch (bool): Whether to save full scratch directory contents
+                at the end of the flattening. Defaults to False.
             **QCJob_kwargs: Passthrough kwargs to QCJob. See
                 :class:`custodian.qchem.jobs.QCJob`.
         """
@@ -178,14 +201,16 @@ class QCJob(Job):
                         output_file=output_file,
                         qclog_file=qclog_file,
                         suffix=".opt_" + str(ii),
-                        scratch_dir=os.getcwd(),
                         save_scratch=True,
-                        save_name="chain_scratch",
                         backup=first,
                         **QCJob_kwargs
                     )
                 )
                 opt_outdata = QCOutput(output_file + ".opt_" + str(ii)).data
+                opt_indata = QCInput.from_file(input_file + ".opt_" + str(ii))
+                if opt_indata.rem["scf_algorithm"] != freq_rem["scf_algorithm"]:
+                    freq_rem["scf_algorithm"] = opt_indata.rem["scf_algorithm"]
+                    opt_rem["scf_algorithm"] = opt_indata.rem["scf_algorithm"]
                 first = False
                 if (
                     opt_outdata["structure_change"] == "unconnected_fragments"
@@ -195,70 +220,70 @@ class QCJob(Job):
                         "Unstable molecule broke into unconnected fragments which failed to optimize! Exiting..."
                     )
                     break
-                else:
-                    energy_history.append(opt_outdata.get("final_energy"))
-                    freq_QCInput = QCInput(
-                        molecule=opt_outdata.get("molecule_from_optimized_geometry"),
-                        rem=freq_rem,
-                        opt=orig_input.opt,
-                        pcm=orig_input.pcm,
-                        solvent=orig_input.solvent,
-                        smx=orig_input.smx,
+                energy_history.append(opt_outdata.get("final_energy"))
+                freq_QCInput = QCInput(
+                    molecule=opt_outdata.get("molecule_from_optimized_geometry"),
+                    rem=freq_rem,
+                    opt=orig_input.opt,
+                    pcm=orig_input.pcm,
+                    solvent=orig_input.solvent,
+                    smx=orig_input.smx,
+                )
+                freq_QCInput.write_file(input_file)
+                yield (
+                    QCJob(
+                        qchem_command=qchem_command,
+                        multimode=multimode,
+                        input_file=input_file,
+                        output_file=output_file,
+                        qclog_file=qclog_file,
+                        suffix=".freq_" + str(ii),
+                        save_scratch=True,
+                        backup=first,
+                        **QCJob_kwargs
                     )
-                    freq_QCInput.write_file(input_file)
-                    yield (
-                        QCJob(
-                            qchem_command=qchem_command,
-                            multimode=multimode,
-                            input_file=input_file,
-                            output_file=output_file,
-                            qclog_file=qclog_file,
-                            suffix=".freq_" + str(ii),
-                            scratch_dir=os.getcwd(),
-                            save_scratch=True,
-                            save_name="chain_scratch",
-                            backup=first,
-                            **QCJob_kwargs
-                        )
+                )
+                outdata = QCOutput(output_file + ".freq_" + str(ii)).data
+                indata = QCInput.from_file(input_file + ".freq_" + str(ii))
+                if indata.rem["scf_algorithm"] != freq_rem["scf_algorithm"]:
+                    freq_rem["scf_algorithm"] = indata.rem["scf_algorithm"]
+                    opt_rem["scf_algorithm"] = indata.rem["scf_algorithm"]
+                errors = outdata.get("errors")
+                if len(errors) != 0:
+                    raise AssertionError(
+                        "No errors should be encountered while flattening frequencies!"
                     )
-                    outdata = QCOutput(output_file + ".freq_" + str(ii)).data
-                    errors = outdata.get("errors")
-                    if len(errors) != 0:
-                        raise AssertionError(
-                            "No errors should be encountered while flattening frequencies!"
-                        )
-                    if outdata.get("frequencies")[0] > 0.0:
-                        print("All frequencies positive!")
-                        break
-                    elif (
-                        abs(outdata.get("frequencies")[0]) < 15.0
-                        and outdata.get("frequencies")[1] > 0.0
+                if outdata.get("frequencies")[0] > 0.0:
+                    print("All frequencies positive!")
+                    break
+                if (
+                    abs(outdata.get("frequencies")[0]) < 15.0
+                    and outdata.get("frequencies")[1] > 0.0
+                ):
+                    print(
+                        "One negative frequency smaller than 15.0 - not worth further flattening!"
+                    )
+                    break
+                if len(energy_history) > 1:
+                    if (
+                        abs(energy_history[-1] - energy_history[-2])
+                        < energy_diff_cutoff
                     ):
-                        print(
-                            "One negative frequency smaller than 15.0 - not worth further flattening!"
-                        )
+                        print("Energy change below cutoff!")
                         break
-                    else:
-                        if len(energy_history) > 1:
-                            if (
-                                abs(energy_history[-1] - energy_history[-2])
-                                < energy_diff_cutoff
-                            ):
-                                print("Energy change below cutoff!")
-                                break
-                        opt_QCInput = QCInput(
-                            molecule=opt_outdata.get(
-                                "molecule_from_optimized_geometry"
-                            ),
-                            rem=opt_rem,
-                            opt=orig_input.opt,
-                            pcm=orig_input.pcm,
-                            solvent=orig_input.solvent,
-                            smx=orig_input.smx,
-                        )
-                        opt_QCInput.write_file(input_file)
-            if os.path.exists(os.path.join(os.getcwd(), "chain_scratch")):
-                shutil.rmtree(os.path.join(os.getcwd(), "chain_scratch"))
+                opt_QCInput = QCInput(
+                    molecule=opt_outdata.get(
+                        "molecule_from_optimized_geometry"
+                    ),
+                    rem=opt_rem,
+                    opt=orig_input.opt,
+                    pcm=orig_input.pcm,
+                    solvent=orig_input.solvent,
+                    smx=orig_input.smx,
+                )
+                opt_QCInput.write_file(input_file)
+            if not save_final_scratch:
+                shutil.rmtree(os.path.join(os.getcwd(), "scratch"))
 
         else:
             if not os.path.exists(input_file):
@@ -298,212 +323,213 @@ class QCJob(Job):
                         "Unstable molecule broke into unconnected fragments which failed to optimize! Exiting..."
                     )
                     break
-                else:
-                    freq_QCInput = QCInput(
-                        molecule=opt_outdata.get("molecule_from_optimized_geometry"),
-                        rem=orig_freq_rem,
-                        opt=orig_opt_input.opt,
-                        pcm=orig_opt_input.pcm,
-                        solvent=orig_opt_input.solvent,
-                        smx=orig_opt_input.smx,
+                freq_QCInput = QCInput(
+                    molecule=opt_outdata.get("molecule_from_optimized_geometry"),
+                    rem=orig_freq_rem,
+                    opt=orig_opt_input.opt,
+                    pcm=orig_opt_input.pcm,
+                    solvent=orig_opt_input.solvent,
+                    smx=orig_opt_input.smx,
+                )
+                freq_QCInput.write_file(input_file)
+                yield (
+                    QCJob(
+                        qchem_command=qchem_command,
+                        multimode=multimode,
+                        input_file=input_file,
+                        output_file=output_file,
+                        qclog_file=qclog_file,
+                        suffix=".freq_" + str(ii),
+                        backup=first,
+                        **QCJob_kwargs
                     )
-                    freq_QCInput.write_file(input_file)
-                    yield (
-                        QCJob(
-                            qchem_command=qchem_command,
-                            multimode=multimode,
-                            input_file=input_file,
-                            output_file=output_file,
-                            qclog_file=qclog_file,
-                            suffix=".freq_" + str(ii),
-                            backup=first,
-                            **QCJob_kwargs
-                        )
+                )
+                outdata = QCOutput(output_file + ".freq_" + str(ii)).data
+                errors = outdata.get("errors")
+                if len(errors) != 0:
+                    raise AssertionError(
+                        "No errors should be encountered while flattening frequencies!"
                     )
-                    outdata = QCOutput(output_file + ".freq_" + str(ii)).data
-                    errors = outdata.get("errors")
-                    if len(errors) != 0:
-                        raise AssertionError(
-                            "No errors should be encountered while flattening frequencies!"
-                        )
-                    if outdata.get("frequencies")[0] > 0.0:
-                        print("All frequencies positive!")
-                        if opt_outdata.get("final_energy") > orig_energy:
-                            print(
-                                "WARNING: Energy increased during frequency flattening!"
-                            )
-                        break
-                    else:
-                        hist = {}
-                        hist["molecule"] = copy.deepcopy(
-                            outdata.get("initial_molecule")
-                        )
-                        hist["geometry"] = copy.deepcopy(
-                            outdata.get("initial_geometry")
-                        )
-                        hist["frequencies"] = copy.deepcopy(outdata.get("frequencies"))
-                        hist["frequency_mode_vectors"] = copy.deepcopy(
-                            outdata.get("frequency_mode_vectors")
-                        )
-                        hist["num_neg_freqs"] = sum(
-                            1 for freq in outdata.get("frequencies") if freq < 0
+                if outdata.get("frequencies")[0] > 0.0:
+                    print("All frequencies positive!")
+                    if opt_outdata.get("final_energy") > orig_energy:
+                        print(
+                            "WARNING: Energy increased during frequency flattening!"
                         )
-                        hist["energy"] = copy.deepcopy(opt_outdata.get("final_energy"))
-                        hist["index"] = len(history)
-                        hist["children"] = []
-                        history.append(hist)
+                    break
+                hist = {}
+                hist["molecule"] = copy.deepcopy(
+                    outdata.get("initial_molecule")
+                )
+                hist["geometry"] = copy.deepcopy(
+                    outdata.get("initial_geometry")
+                )
+                hist["frequencies"] = copy.deepcopy(outdata.get("frequencies"))
+                hist["frequency_mode_vectors"] = copy.deepcopy(
+                    outdata.get("frequency_mode_vectors")
+                )
+                hist["num_neg_freqs"] = sum(
+                    1 for freq in outdata.get("frequencies") if freq < 0
+                )
+                hist["energy"] = copy.deepcopy(opt_outdata.get("final_energy"))
+                hist["index"] = len(history)
+                hist["children"] = []
+                history.append(hist)
 
-                        ref_mol = history[-1]["molecule"]
-                        geom_to_perturb = history[-1]["geometry"]
-                        negative_freq_vecs = history[-1]["frequency_mode_vectors"][0]
-                        reversed_direction = False
-                        standard = True
+                ref_mol = history[-1]["molecule"]
+                geom_to_perturb = history[-1]["geometry"]
+                negative_freq_vecs = history[-1]["frequency_mode_vectors"][0]
+                reversed_direction = False
+                standard = True
 
-                        # If we've found one or more negative frequencies in two consecutive iterations, let's dig in
-                        # deeper:
-                        if len(history) > 1:
-                            # Start by finding the latest iteration's parent:
-                            if history[-1]["index"] in history[-2]["children"]:
-                                parent_hist = history[-2]
-                                history[-1]["parent"] = parent_hist["index"]
-                            elif history[-1]["index"] in history[-3]["children"]:
-                                parent_hist = history[-3]
-                                history[-1]["parent"] = parent_hist["index"]
-                            else:
-                                raise AssertionError(
-                                    "ERROR: your parent should always be one or two iterations behind you! Exiting..."
-                                )
+                # If we've found one or more negative frequencies in two consecutive iterations, let's dig in
+                # deeper:
+                if len(history) > 1:
+                    # Start by finding the latest iteration's parent:
+                    if history[-1]["index"] in history[-2]["children"]:
+                        parent_hist = history[-2]
+                        history[-1]["parent"] = parent_hist["index"]
+                    elif history[-1]["index"] in history[-3]["children"]:
+                        parent_hist = history[-3]
+                        history[-1]["parent"] = parent_hist["index"]
+                    else:
+                        raise AssertionError(
+                            "ERROR: your parent should always be one or two iterations behind you! Exiting..."
+                        )
 
-                            # if the number of negative frequencies has remained constant or increased from parent to
-                            # child,
-                            if (
-                                history[-1]["num_neg_freqs"]
-                                >= parent_hist["num_neg_freqs"]
-                            ):
-                                # check to see if the parent only has one child, aka only the positive perturbation has
-                                # been tried,
-                                # in which case just try the negative perturbation from the same parent
-                                if len(parent_hist["children"]) == 1:
-                                    ref_mol = parent_hist["molecule"]
-                                    geom_to_perturb = parent_hist["geometry"]
-                                    negative_freq_vecs = parent_hist[
-                                        "frequency_mode_vectors"
-                                    ][0]
-                                    reversed_direction = True
-                                    standard = False
-                                    parent_hist["children"].append(len(history))
-                                # If the parent has two children, aka both directions have been tried, then we have to
-                                # get creative:
-                                elif len(parent_hist["children"]) == 2:
-                                    # If we're dealing with just one negative frequency,
-                                    if parent_hist["num_neg_freqs"] == 1:
-                                        make_good_child_next_parent = False
-                                        if (
-                                            history[parent_hist["children"][0]][
-                                                "energy"
-                                            ]
-                                            < history[-1]["energy"]
-                                        ):
-                                            good_child = copy.deepcopy(
-                                                history[parent_hist["children"][0]]
-                                            )
-                                        else:
-                                            good_child = copy.deepcopy(history[-1])
-                                        if good_child["num_neg_freqs"] > 1:
-                                            raise Exception(
-                                                "ERROR: Child with lower energy has more negative frequencies! "
-                                                "Exiting..."
-                                            )
-                                        elif (
-                                            good_child["energy"] < parent_hist["energy"]
-                                        ):
-                                            make_good_child_next_parent = True
-                                        elif (
-                                            vector_list_diff(
-                                                good_child["frequency_mode_vectors"][0],
-                                                parent_hist["frequency_mode_vectors"][
-                                                    0
-                                                ],
-                                            )
-                                            > 0.2
-                                        ):
-                                            make_good_child_next_parent = True
-                                        else:
-                                            raise Exception(
-                                                "ERROR: Good child not good enough! Exiting..."
-                                            )
-                                        if make_good_child_next_parent:
-                                            good_child["index"] = len(history)
-                                            history.append(good_child)
-                                            ref_mol = history[-1]["molecule"]
-                                            geom_to_perturb = history[-1]["geometry"]
-                                            negative_freq_vecs = history[-1][
-                                                "frequency_mode_vectors"
-                                            ][0]
-                                    else:
-                                        raise Exception(
-                                            "ERROR: Can't deal with multiple neg frequencies yet! Exiting..."
-                                        )
+                    # if the number of negative frequencies has remained constant or increased from parent to
+                    # child,
+                    if (
+                        history[-1]["num_neg_freqs"]
+                        >= parent_hist["num_neg_freqs"]
+                    ):
+                        # check to see if the parent only has one child, aka only the positive perturbation has
+                        # been tried,
+                        # in which case just try the negative perturbation from the same parent
+                        if len(parent_hist["children"]) == 1:
+                            ref_mol = parent_hist["molecule"]
+                            geom_to_perturb = parent_hist["geometry"]
+                            negative_freq_vecs = parent_hist[
+                                "frequency_mode_vectors"
+                            ][0]
+                            reversed_direction = True
+                            standard = False
+                            parent_hist["children"].append(len(history))
+                        # If the parent has two children, aka both directions have been tried, then we have to
+                        # get creative:
+                        elif len(parent_hist["children"]) == 2:
+                            # If we're dealing with just one negative frequency,
+                            if parent_hist["num_neg_freqs"] == 1:
+                                make_good_child_next_parent = False
+                                if (
+                                    history[parent_hist["children"][0]][
+                                        "energy"
+                                    ]
+                                    < history[-1]["energy"]
+                                ):
+                                    good_child = copy.deepcopy(
+                                        history[parent_hist["children"][0]]
+                                    )
                                 else:
-                                    raise AssertionError(
-                                        "ERROR: Parent cannot have more than two childen! Exiting..."
+                                    good_child = copy.deepcopy(history[-1])
+                                if good_child["num_neg_freqs"] > 1:
+                                    raise Exception(
+                                        "ERROR: Child with lower energy has more negative frequencies! "
+                                        "Exiting..."
                                     )
-                            # Implicitly, if the number of negative frequencies decreased from parent to child,
-                            # continue normally.
-                        if standard:
-                            history[-1]["children"].append(len(history))
-
-                        min_molecule_perturb_scale = 0.1
-                        scale_grid = 10
-                        perturb_scale_grid = (
-                            max_molecule_perturb_scale - min_molecule_perturb_scale
-                        ) / scale_grid
-
-                        structure_successfully_perturbed = False
-                        for molecule_perturb_scale in np.arange(
-                            max_molecule_perturb_scale,
-                            min_molecule_perturb_scale,
-                            -perturb_scale_grid,
-                        ):
-                            new_coords = perturb_coordinates(
-                                old_coords=geom_to_perturb,
-                                negative_freq_vecs=negative_freq_vecs,
-                                molecule_perturb_scale=molecule_perturb_scale,
-                                reversed_direction=reversed_direction,
-                            )
-                            new_molecule = Molecule(
-                                species=orig_species,
-                                coords=new_coords,
-                                charge=orig_charge,
-                                spin_multiplicity=orig_multiplicity,
-                            )
-                            if check_connectivity:
-                                structure_successfully_perturbed = (
-                                    check_for_structure_changes(ref_mol, new_molecule)
-                                    == "no_change"
+                                if (
+                                    good_child["energy"] < parent_hist["energy"]
+                                ):
+                                    make_good_child_next_parent = True
+                                elif (
+                                    vector_list_diff(
+                                        good_child["frequency_mode_vectors"][0],
+                                        parent_hist["frequency_mode_vectors"][
+                                            0
+                                        ],
+                                    )
+                                    > 0.2
+                                ):
+                                    make_good_child_next_parent = True
+                                else:
+                                    raise Exception(
+                                        "ERROR: Good child not good enough! Exiting..."
+                                    )
+                                if make_good_child_next_parent:
+                                    good_child["index"] = len(history)
+                                    history.append(good_child)
+                                    ref_mol = history[-1]["molecule"]
+                                    geom_to_perturb = history[-1]["geometry"]
+                                    negative_freq_vecs = history[-1][
+                                        "frequency_mode_vectors"
+                                    ][0]
+                            else:
+                                raise Exception(
+                                    "ERROR: Can't deal with multiple neg frequencies yet! Exiting..."
                                 )
-                                if structure_successfully_perturbed:
-                                    break
-                        if not structure_successfully_perturbed:
-                            raise Exception(
-                                "ERROR: Unable to perturb coordinates to remove negative frequency without changing "
-                                "the connectivity! Exiting..."
+                        else:
+                            raise AssertionError(
+                                "ERROR: Parent cannot have more than two childen! Exiting..."
                             )
+                    # Implicitly, if the number of negative frequencies decreased from parent to child,
+                    # continue normally.
+                if standard:
+                    history[-1]["children"].append(len(history))
 
-                        new_opt_QCInput = QCInput(
-                            molecule=new_molecule,
-                            rem=orig_opt_rem,
-                            opt=orig_opt_input.opt,
-                            pcm=orig_opt_input.pcm,
-                            solvent=orig_opt_input.solvent,
-                            smx=orig_opt_input.smx,
+                min_molecule_perturb_scale = 0.1
+                scale_grid = 10
+                perturb_scale_grid = (
+                    max_molecule_perturb_scale - min_molecule_perturb_scale
+                ) / scale_grid
+
+                structure_successfully_perturbed = False
+                for molecule_perturb_scale in np.arange(
+                    max_molecule_perturb_scale,
+                    min_molecule_perturb_scale,
+                    -perturb_scale_grid,
+                ):
+                    new_coords = perturb_coordinates(
+                        old_coords=geom_to_perturb,
+                        negative_freq_vecs=negative_freq_vecs,
+                        molecule_perturb_scale=molecule_perturb_scale,
+                        reversed_direction=reversed_direction,
+                    )
+                    new_molecule = Molecule(
+                        species=orig_species,
+                        coords=new_coords,
+                        charge=orig_charge,
+                        spin_multiplicity=orig_multiplicity,
+                    )
+                    if check_connectivity:
+                        structure_successfully_perturbed = (
+                            check_for_structure_changes(ref_mol, new_molecule)
+                            == "no_change"
                         )
-                        new_opt_QCInput.write_file(input_file)
+                        if structure_successfully_perturbed:
+                            break
+                if not structure_successfully_perturbed:
+                    raise Exception(
+                        "ERROR: Unable to perturb coordinates to remove negative frequency without changing "
+                        "the connectivity! Exiting..."
+                    )
+
+                new_opt_QCInput = QCInput(
+                    molecule=new_molecule,
+                    rem=orig_opt_rem,
+                    opt=orig_opt_input.opt,
+                    pcm=orig_opt_input.pcm,
+                    solvent=orig_opt_input.solvent,
+                    smx=orig_opt_input.smx,
+                )
+                new_opt_QCInput.write_file(input_file)
 
 
 def perturb_coordinates(
     old_coords, negative_freq_vecs, molecule_perturb_scale, reversed_direction
 ):
+    """
+    Perturbs a structure along the imaginary mode vibrational frequency vectors
+    """
     max_dis = max([math.sqrt(sum([x ** 2 for x in vec])) for vec in negative_freq_vecs])
     scale = molecule_perturb_scale / max_dis
     normalized_vecs = [[x * scale for x in vec] for vec in negative_freq_vecs]
@@ -517,6 +543,9 @@ def perturb_coordinates(
 
 
 def vector_list_diff(vecs1, vecs2):
+    """
+    Calculates the summed difference of two vectors
+    """
     diff = 0.0
     if len(vecs1) != len(vecs2):
         raise AssertionError("ERROR: Vectors must be of equal length! Exiting...")


Test patch:
diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml
index e1fd05f2..4428a41b 100644
--- a/.github/workflows/pytest.yml
+++ b/.github/workflows/pytest.yml
@@ -14,18 +14,24 @@ jobs:
 
     steps:
     - uses: actions/checkout@v2
-    - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v2
+
+    - uses: conda-incubator/setup-miniconda@v2
       with:
         python-version: ${{ matrix.python-version }}
+
+    - name: Install OpenBabel
+      shell: bash -l {0}
+      run: conda install -c conda-forge openbabel
+
     - name: Install dependencies
+      shell: bash -l {0}
       run: |
-        sudo apt-get install python-openbabel
         python -m pip install --upgrade pip
         pip install -r requirements.txt
         pip install -r requirements-ci.txt
         pip install -e .
     - name: pytest
+      shell: bash -l {0}
       env:
         PMG_MAPI_KEY: ${{ secrets.PMG_MAPI_KEY }}
         MPLBACKEND: "Agg"
diff --git a/custodian/qchem/tests/__init__.py b/custodian/qchem/tests/__init__.py
index 148cd9d2..075a2a56 100644
--- a/custodian/qchem/tests/__init__.py
+++ b/custodian/qchem/tests/__init__.py
@@ -1,3 +1,8 @@
+"""
+Tests for QChem Custodian
+"""
+
 from __future__ import unicode_literals
 
-__author__ = "xiaohuiqu"
+__author__ = "Samuel Blau, Brandon Wood, Shyam Dwaraknath"
+__credits__ = "Xiaohui Qu"
diff --git a/custodian/qchem/tests/test_jobs.py b/custodian/qchem/tests/test_jobs.py
index 82d89c60..7afca8dd 100644
--- a/custodian/qchem/tests/test_jobs.py
+++ b/custodian/qchem/tests/test_jobs.py
@@ -36,7 +36,7 @@ class QCJobTest(TestCase):
     def test_defaults(self):
         with patch("custodian.qchem.jobs.shutil.copy") as copy_patch:
             myjob = QCJob(qchem_command="qchem", max_cores=32)
-            self.assertEqual(myjob.current_command, " qchem -nt 32 mol.qin mol.qout")
+            self.assertEqual(myjob.current_command, " qchem -nt 32 mol.qin mol.qout scratch")
             myjob.setup()
             self.assertEqual(copy_patch.call_args_list[0][0][0], "mol.qin")
             self.assertEqual(copy_patch.call_args_list[0][0][1], "mol.qin.orig")
@@ -51,27 +51,27 @@ class QCJobTest(TestCase):
             input_file="different.qin",
             output_file="not_default.qout",
             max_cores=12,
-            scratch_dir="/not/default/scratch/",
+            calc_loc="/not/default/",
             backup=False,
         )
         self.assertEqual(
-            myjob.current_command, " qchem -slurm -np 12 different.qin not_default.qout"
+            myjob.current_command, " qchem -slurm -np 12 different.qin not_default.qout scratch"
         )
         myjob.setup()
-        self.assertEqual(os.environ["QCSCRATCH"], "/not/default/scratch/")
+        self.assertEqual(os.environ["QCSCRATCH"], os.getcwd())
+        self.assertEqual(os.environ["QCLOCALSCR"], "/not/default/")
 
     def test_save_scratch(self):
         with patch("custodian.qchem.jobs.shutil.copy") as copy_patch:
             myjob = QCJob(
                 qchem_command="qchem -slurm",
                 max_cores=32,
-                scratch_dir=os.getcwd(),
+                calc_loc="/tmp/scratch",
                 save_scratch=True,
-                save_name="freq_scratch",
             )
             self.assertEqual(
                 myjob.current_command,
-                " qchem -slurm -nt 32 mol.qin mol.qout freq_scratch",
+                " qchem -slurm -nt 32 mol.qin mol.qout scratch",
             )
             myjob.setup()
             self.assertEqual(copy_patch.call_args_list[0][0][0], "mol.qin")
@@ -79,6 +79,7 @@ class QCJobTest(TestCase):
             self.assertEqual(os.environ["QCSCRATCH"], os.getcwd())
             self.assertEqual(os.environ["QCTHREADS"], "32")
             self.assertEqual(os.environ["OMP_NUM_THREADS"], "32")
+            self.assertEqual(os.environ["QCLOCALSCR"], "/tmp/scratch")
 
 
 class OptFFTest(TestCase):
@@ -829,10 +830,15 @@ class OptFFTest5690(TestCase):
 class OptFF_small_neg_freq(TestCase):
     def setUp(self):
         os.makedirs(scr_dir)
+        os.makedirs(os.path.join(scr_dir,"scratch"))
         shutil.copyfile(
             os.path.join(test_dir, "small_neg_freq/mol.qin.orig"),
             os.path.join(scr_dir, "mol.qin"),
         )
+        shutil.copyfile(
+            os.path.join(test_dir, "small_neg_freq/mol.qin.opt_0"),
+            os.path.join(scr_dir, "mol.qin.opt_0"),
+        )
         shutil.copyfile(
             os.path.join(test_dir, "small_neg_freq/mol.qout.opt_0"),
             os.path.join(scr_dir, "mol.qout.opt_0"),
@@ -878,9 +884,7 @@ class OptFF_small_neg_freq(TestCase):
             input_file="mol.qin",
             output_file="mol.qout",
             suffix=".opt_0",
-            scratch_dir=os.getcwd(),
             save_scratch=True,
-            save_name="chain_scratch",
             backup=True,
         ).as_dict()
         self.assertEqual(next(myjob).as_dict(), expected_next)
@@ -891,9 +895,7 @@ class OptFF_small_neg_freq(TestCase):
             input_file="mol.qin",
             output_file="mol.qout",
             suffix=".freq_0",
-            scratch_dir=os.getcwd(),
             save_scratch=True,
-            save_name="chain_scratch",
             backup=False,
         ).as_dict()
         self.assertEqual(next(myjob).as_dict(), expected_next)
@@ -903,6 +905,10 @@ class OptFF_small_neg_freq(TestCase):
             ).as_dict(),
             QCInput.from_file(os.path.join(scr_dir, "mol.qin")).as_dict(),
         )
+        shutil.copyfile(
+            os.path.join(scr_dir, "mol.qin"),
+            os.path.join(scr_dir, "mol.qin.freq_0"),
+        )
         expected_next = QCJob(
             qchem_command="qchem -slurm",
             max_cores=32,
@@ -910,9 +916,7 @@ class OptFF_small_neg_freq(TestCase):
             input_file="mol.qin",
             output_file="mol.qout",
             suffix=".opt_1",
-            scratch_dir=os.getcwd(),
             save_scratch=True,
-            save_name="chain_scratch",
             backup=False,
         ).as_dict()
         self.assertEqual(next(myjob).as_dict(), expected_next)
@@ -922,6 +926,10 @@ class OptFF_small_neg_freq(TestCase):
             ).as_dict(),
             QCInput.from_file(os.path.join(scr_dir, "mol.qin")).as_dict(),
         )
+        shutil.copyfile(
+            os.path.join(scr_dir, "mol.qin"),
+            os.path.join(scr_dir, "mol.qin.opt_1"),
+        )
         expected_next = QCJob(
             qchem_command="qchem -slurm",
             max_cores=32,
@@ -929,9 +937,7 @@ class OptFF_small_neg_freq(TestCase):
             input_file="mol.qin",
             output_file="mol.qout",
             suffix=".freq_1",
-            scratch_dir=os.getcwd(),
             save_scratch=True,
-            save_name="chain_scratch",
             backup=False,
         ).as_dict()
         self.assertEqual(next(myjob).as_dict(), expected_next)
@@ -941,6 +947,10 @@ class OptFF_small_neg_freq(TestCase):
             ).as_dict(),
             QCInput.from_file(os.path.join(scr_dir, "mol.qin")).as_dict(),
         )
+        shutil.copyfile(
+            os.path.join(scr_dir, "mol.qin"),
+            os.path.join(scr_dir, "mol.qin.freq_1"),
+        )
         expected_next = QCJob(
             qchem_command="qchem -slurm",
             max_cores=32,
@@ -948,9 +958,7 @@ class OptFF_small_neg_freq(TestCase):
             input_file="mol.qin",
             output_file="mol.qout",
             suffix=".opt_2",
-            scratch_dir=os.getcwd(),
             save_scratch=True,
-            save_name="chain_scratch",
             backup=False,
         ).as_dict()
         self.assertEqual(next(myjob).as_dict(), expected_next)
@@ -960,6 +968,10 @@ class OptFF_small_neg_freq(TestCase):
             ).as_dict(),
             QCInput.from_file(os.path.join(scr_dir, "mol.qin")).as_dict(),
         )
+        shutil.copyfile(
+            os.path.join(scr_dir, "mol.qin"),
+            os.path.join(scr_dir, "mol.qin.opt_2"),
+        )
         expected_next = QCJob(
             qchem_command="qchem -slurm",
             max_cores=32,
@@ -967,9 +979,7 @@ class OptFF_small_neg_freq(TestCase):
             input_file="mol.qin",
             output_file="mol.qout",
             suffix=".freq_2",
-            scratch_dir=os.getcwd(),
             save_scratch=True,
-            save_name="chain_scratch",
             backup=False,
         ).as_dict()
         self.assertEqual(next(myjob).as_dict(), expected_next)
@@ -979,6 +989,10 @@ class OptFF_small_neg_freq(TestCase):
             ).as_dict(),
             QCInput.from_file(os.path.join(scr_dir, "mol.qin")).as_dict(),
         )
+        shutil.copyfile(
+            os.path.join(scr_dir, "mol.qin"),
+            os.path.join(scr_dir, "mol.qin.freq_2"),
+        )
         self.assertRaises(StopIteration, myjob.__next__)