def opt_with_frequency_flattener(cls,
                                     qchem_command,
                                     multimode="openmp",
                                     input_file="mol.qin",
                                     output_file="mol.qout",
                                     qclog_file="mol.qclog",
                                     max_iterations=10,
                                     max_molecule_perturb_scale=0.3,
                                     check_connectivity=True,
                                     **QCJob_kwargs):
        """
        Optimize a structure and calculate vibrational frequencies to check if the
        structure is in a true minima. If a frequency is negative, iteratively
        perturbe the geometry, optimize, and recalculate frequencies until all are
        positive, aka a true minima has been found.

        Args:
            qchem_command (str): Command to run QChem.
            multimode (str): Parallelization scheme, either openmp or mpi.
            input_file (str): Name of the QChem input file.
            output_file (str): Name of the QChem output file.
            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.
            check_connectivity (bool): Whether to check differences in connectivity
                introduced by structural perturbation. Defaults to True.
            **QCJob_kwargs: Passthrough kwargs to QCJob. See
                :class:`custodian.qchem.jobs.QCJob`.
        """

        min_molecule_perturb_scale = 0.1
        scale_grid = 10
        perturb_scale_grid = (
            max_molecule_perturb_scale - min_molecule_perturb_scale
        ) / scale_grid

        if not os.path.exists(input_file):
            raise AssertionError('Input file must be present!')
        orig_opt_input = QCInput.from_file(input_file)
        orig_opt_rem = copy.deepcopy(orig_opt_input.rem)
        orig_freq_rem = copy.deepcopy(orig_opt_input.rem)
        orig_freq_rem["job_type"] = "freq"
        first = True
        reversed_direction = False
        num_neg_freqs = []

        for ii in range(max_iterations):
            yield (QCJob(
                qchem_command=qchem_command,
                multimode=multimode,
                input_file=input_file,
                output_file=output_file,
                qclog_file=qclog_file,
                suffix=".opt_" + str(ii),
                backup=first,
                **QCJob_kwargs))
            first = False
            opt_outdata = QCOutput(output_file + ".opt_" + str(ii)).data
            if opt_outdata["structure_change"] == "unconnected_fragments" and not opt_outdata["completion"]:
                print("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)
                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!')
                if outdata.get('frequencies')[0] > 0.0:
                    print("All frequencies positive!")
                    break
                else:
                    num_neg_freqs += [sum(1 for freq in outdata.get('frequencies') if freq < 0)]
                    if len(num_neg_freqs) > 1:
                        if num_neg_freqs[-1] == num_neg_freqs[-2] and not reversed_direction:
                            reversed_direction = True
                        elif num_neg_freqs[-1] == num_neg_freqs[-2] and reversed_direction:
                            if len(num_neg_freqs) < 3:
                                raise AssertionError("ERROR: This should only be possible after at least three frequency flattening iterations! Exiting...")
                            else:
                                raise Exception("ERROR: Reversing the perturbation direction still could not flatten any frequencies. Exiting...")
                        elif num_neg_freqs[-1] != num_neg_freqs[-2] and reversed_direction:
                            reversed_direction = False

                    negative_freq_vecs = outdata.get("frequency_mode_vectors")[0]
                    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=outdata.get("initial_geometry"),
                            negative_freq_vecs=negative_freq_vecs,
                            molecule_perturb_scale=molecule_perturb_scale,
                            reversed_direction=reversed_direction)
                        new_molecule = Molecule(
                            species=outdata.get('species'),
                            coords=new_coords,
                            charge=outdata.get('charge'),
                            spin_multiplicity=outdata.get('multiplicity'))
                        if check_connectivity:
                            old_molgraph = MoleculeGraph.with_local_env_strategy(outdata.get("initial_molecule"),
                                                               OpenBabelNN(),
                                                               reorder=False,
                                                               extend_structure=False)
                            new_molgraph = MoleculeGraph.with_local_env_strategy(new_molecule,
                                                               OpenBabelNN(),
                                                               reorder=False,
                                                               extend_structure=False)
                            if old_molgraph.isomorphic_to(new_molgraph):
                                structure_successfully_perturbed = True
                                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)
                    new_opt_QCInput.write_file(input_file)