diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 0e03a85..559ff36 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -6,16 +6,13 @@ jobs:
     strategy:
       matrix:
         cfg:
-          - conda-env: base
-            python-version: 3.9
-          - conda-env: base
-            python-version: "3.10"
           - conda-env: base
             python-version: "3.11"
           - conda-env: base
             python-version: "3.12"
           - conda-env: base
             python-version: "3.13"
+            pytest-version: "9.0.1"
     env:
       PYVER: ${{ matrix.cfg.python-version }}
       CONDA_ENV: ${{ matrix.cfg.conda-env }}
diff --git a/concordantmodes/cma.py b/concordantmodes/cma.py
index 0b7f00f..555f65e 100644
--- a/concordantmodes/cma.py
+++ b/concordantmodes/cma.py
@@ -26,8 +26,9 @@
 from concordantmodes.symmetry import Symmetry
 from concordantmodes.ted import TED
 from concordantmodes.transf_disp import TransfDisp
-from concordantmodes.vulcan_template import VulcanTemplate
-from concordantmodes.sapelo_template import SapeloTemplate
+
+# from concordantmodes.vulcan_template import VulcanTemplate
+# from concordantmodes.sapelo_template import SapeloTemplate
 from concordantmodes.zmat import Zmat
 
 
@@ -81,22 +82,27 @@ def run(self, sym_sort=[]):
         s_vec = SVectors(
             self.zmat_obj, self.options, self.zmat_obj.variable_dictionary_b
         )
-        if self.options.man_proj:
-            proj = self.proj
-            np.set_printoptions(precision=4, linewidth=240)
-            s_vec.run(
-                self.zmat_obj.cartesians_b,
-                True,
-                proj=proj,
-                second_order=self.options.second_order,
-            )
-            s_vec.proj = proj
-        else:
-            s_vec.run(
-                self.zmat_obj.cartesians_b,
-                True,
-                second_order=self.options.second_order,
-            )
+        s_vec.run(
+            self.zmat_obj.cartesians_b,
+            True,
+            proj=self.proj,
+            second_order=self.options.second_order,
+        )
+        # if self.options.man_proj:
+        # proj = self.proj
+        # s_vec.run(
+        # self.zmat_obj.cartesians_b,
+        # True,
+        # proj=proj,
+        # second_order=self.options.second_order,
+        # )
+        # # s_vec.proj = proj
+        # else:
+        # s_vec.run(
+        # self.zmat_obj.cartesians_b,
+        # True,
+        # second_order=self.options.second_order,
+        # )
 
         if self.options.molsym_symmetry:
             self.symm_obj.make_proj(s_vec)
@@ -105,10 +111,10 @@ def run(self, sym_sort=[]):
         self.TED_obj = TED(s_vec.proj, self.zmat_obj, self.options)
 
         # Print out the percentage composition of the projected coordinates
-        if self.options.coords != "ZMAT":
-            self.TED_obj.run(
-                np.eye(len(self.TED_obj.proj.T)), np.zeros(len(self.TED_obj.proj.T))
-            )
+        # if self.options.coords != "ZMAT":
+        # self.TED_obj.run(
+        # np.eye(len(self.TED_obj.proj.T)), np.zeros(len(self.TED_obj.proj.T))
+        # )
 
         # Compute G-Matrix
         g_mat = GMatrix(self.zmat_obj, s_vec, self.options)
@@ -119,7 +125,7 @@ def run(self, sym_sort=[]):
         if os.path.exists(rootdir + "/fc_b.grad"):
             g_read_obj = GrRead("fc_b.grad")
             g_read_obj.run(self.zmat_obj.cartesians_b)
-            print(self.zmat_obj.cartesians_b)
+            # print(self.zmat_obj.cartesians_b)
 
         num_deg_free = s_vec.proj.shape[1]
         # Read in FC matrix in cartesians, then convert to internals.
@@ -127,6 +133,7 @@ def run(self, sym_sort=[]):
         self.options.init_bool = False
         if os.path.exists(rootdir + "/fc_b.dat"):
             f_read_obj = FcRead("fc_b.dat")
+            self.options.cart_fc_b = True
         elif os.path.exists(rootdir + "/FCMFINAL"):
             f_read_obj = FcRead("FCMFINAL")
         else:
@@ -138,10 +145,18 @@ def run(self, sym_sort=[]):
             if self.options.cart_fc_b:
                 eigs_b = np.eye(len(self.zmat_obj.cartesians_b.flatten()))
                 coord_type = "cartesian"
+
+            algo = Algorithm(
+                num_deg_free,
+                cma_level,
+                self.options,
+                proj_irreps=self.symm_obj.proj_irreps,
+            )
+            algo.run()
             if not self.options.deriv_level_b:
                 if not self.options.molsym_symmetry:
-                    algo = Algorithm(num_deg_free, cma_level, self.options)
-                    algo.run()
+                    # algo = Algorithm(num_deg_free, cma_level, self.options)
+                    # algo.run()
                     # Sym_sort doesn't seem to be working
                     print("symmetric displacements:")
                     if len(sym_sort) > 1:
@@ -150,23 +165,23 @@ def run(self, sym_sort=[]):
                         )
                         # indices = self.symm_obj.create_sym_sort_disps(sym_sort, algo.indices)
                 else:
-                    algo = Algorithm(
-                        num_deg_free,
-                        cma_level,
-                        self.options,
-                        proj_irreps=self.symm_obj.proj_irreps,
-                    )
-                    algo.run()
+                    # algo = Algorithm(
+                    # num_deg_free,
+                    # cma_level,
+                    # self.options,
+                    # proj_irreps=self.symm_obj.proj_irreps,
+                    # )
+                    # algo.run()
                     self.symm_obj.indices_by_irrep = algo.indices_by_irrep
             else:
-                self.symm_obj.proj_irreps = None
-                algo = Algorithm(
-                    num_deg_free,
-                    cma_level,
-                    self.options,
-                    proj_irreps=self.symm_obj.proj_irreps,
-                )
-                algo.run()
+                # self.symm_obj.proj_irreps = None
+                # algo = Algorithm(
+                # num_deg_free,
+                # cma_level,
+                # self.options,
+                # proj_irreps=self.symm_obj.proj_irreps,
+                # )
+                # algo.run()
                 self.symm_obj.indices_by_irrep = algo.indices_by_irrep
 
             b_disp = TransfDisp(
@@ -188,11 +203,11 @@ def run(self, sym_sort=[]):
             prog_name_b = prog_b.split("@")[0]
 
             if self.options.gen_disps_b:
-                if os.path.exists(rootdir + "/DispsB"):
-                    if os.path.exists(rootdir + "/oldDispsB"):
-                        shutil.rmtree(rootdir + "/oldDispsB")
-                    shutil.copytree(rootdir + "/DispsB", rootdir + "/oldDispsB")
-                    shutil.rmtree(rootdir + "/DispsB")
+                # if os.path.exists(rootdir + "/DispsB"):
+                # if os.path.exists(rootdir + "/oldDispsB"):
+                # shutil.rmtree(rootdir + "/oldDispsB")
+                # shutil.copytree(rootdir + "/DispsB", rootdir + "/oldDispsB")
+                # shutil.rmtree(rootdir + "/DispsB")
 
                 ref_geom_b = b_disp.disp_cart["ref"]
 
@@ -227,37 +242,8 @@ def run(self, sym_sort=[]):
             os.chdir(rootdir + "/DispsB")
 
             if self.options.calc_b:
-                disp_list = []
-                for i in os.listdir(rootdir + "/DispsB"):
-                    disp_list.append(i)
-
-                if self.options.cluster != "sapelo":
-                    v_template = VulcanTemplate(
-                        self.options, len(disp_list), prog_name_b, prog_b
-                    )
-                    out = v_template.run()
-                    with open("displacements.sh", "w") as file:
-                        file.write(out)
-
-                    # Submits an array, then checks if all jobs have finished every
-                    # 30 seconds.
-                    sub = Submit(disp_list, self.options)
-                    sub.run()
-                else:
-                    s_template = SapeloTemplate(
-                        self.options, len(disp_list), prog_name_b, prog_b
-                    )
-                    out = s_template.run()
-                    with open("optstep.sh", "w") as file:
-                        file.write(out)
-                    for z in range(0, len(disp_list)):
-                        source = os.getcwd() + "/optstep.sh"
-                        os.chdir("./" + str(z + 1))
-                        destination = os.getcwd()
-                        shutil.copy2(source, destination)
-                        os.chdir("../")
-                    sub = Submit(disp_list, self.options)
-                    sub.run()
+                sub = Submit(self.options, cma_level, rootdir, prog_name_b, prog_b)
+                sub.run()
 
             reap_obj_b = Reap(
                 self.options,
@@ -269,6 +255,7 @@ def run(self, sym_sort=[]):
             )
             reap_obj_b.run()
 
+            # can this be folded into reap obj?
             if not self.options.deriv_level_b:
                 p_array_b = reap_obj_b.p_en_array
                 m_array_b = reap_obj_b.m_en_array
@@ -282,6 +269,9 @@ def run(self, sym_sort=[]):
                 ref_en_b = None
 
                 # Need to convert this array here from cartesians to internals using projected A-tensor
+                # Surely this can be folded into reap_obj. Then we simply define p_array and m_array so
+                # so force_constant can read them in and determine what do do with them based off of
+                # deriv_level.
                 for i in range(len(algo.indices)):
                     grad_s_vec = SVectors(
                         self.zmat_obj,
@@ -311,9 +301,35 @@ def run(self, sym_sort=[]):
                 cma_level=cma_level,
             )
             fc_b.run()
-            print("Computed Force Constants:")
-            print(fc_b.FC)
-            if self.options.second_order:
+            # print("Computed Force Constants:")
+            # print(fc_b.FC.shape)
+            # print(fc_b.FC)
+            # np.set_printoptions(precision=6, linewidth=240)
+            # print("Computed Gradient:")
+            # print(fc_b.gradient.shape)
+            # print(fc_b.gradient)
+
+            f_conv_test = FcConv(
+                fc_b.FC,
+                s_vec,
+                self.zmat_obj,
+                "cartesian",
+                False,
+                self.TED_obj,
+                self.options,
+            )
+            f_conv_test.run(grad=fc_b.gradient)
+            # self.F_cart_b = f_conv_test.F
+            # self.grad_cart_b = f_conv_test.grad
+            # print("Cartesian Force Constants:")
+            # print(f_conv_test.F.shape)
+            # print(f_conv_test.F)
+            # print("Computed Gradient:")
+            # print(f_conv_test.grad.shape)
+            # print(f_conv_test.grad.reshape(-1,3))
+            # raise RuntimeError
+
+            if self.options.second_order and self.options.cart_fc_b:
                 f_conv_obj = FcConv(
                     fc_b.FC,
                     s_vec,
@@ -326,10 +342,6 @@ def run(self, sym_sort=[]):
                 f_conv_obj.run(grad=fc_b.gradient)
                 fc_b.FC = f_conv_obj.F
 
-        # Temporary code to ensure nothing breaks
-        self.options.deriv_level_b = 0
-        self.options.deriv_level_a = 0
-
         if not self.options.init_bool:
             f_read_obj.run()
             f_conv_obj = FcConv(
@@ -349,8 +361,10 @@ def run(self, sym_sort=[]):
         else:
             F = fc_b.FC
 
-        if self.options.coords != "ZMAT" and not self.options.init_bool:
-            F = np.dot(self.TED_obj.proj.T, np.dot(F, self.TED_obj.proj))
+        # Folded into F_conv
+        # if self.options.coords != "ZMAT" and not self.options.init_bool:
+        # F = np.dot(self.TED_obj.proj.T, np.dot(F, self.TED_obj.proj))
+        # Fold into g_mat
         if self.options.coords != "ZMAT":
             g_mat.G = np.dot(self.TED_obj.proj.T, np.dot(g_mat.G, self.TED_obj.proj))
 
@@ -375,6 +389,7 @@ def run(self, sym_sort=[]):
             abs_row = np.abs(row)
             row[abs_row < np.max(abs_row) * del_tol] = 0
         F = F.T
+        # raise RuntimeError
 
         if len(sym_sort) > 1:
             F, g_mat.G = self.symm_obj.GF_sym_sort(F, g_mat, sym_sort)
@@ -391,10 +406,6 @@ def run(self, sym_sort=[]):
         )
         b_GF.run()
 
-        # For a quick test, to be deleted later
-        g_b = g_mat.G.copy()
-        f_b = F.copy()
-
         ted_b = b_GF.ted.TED
         # more sym_sort stuff
         if len(sym_sort):
@@ -406,9 +417,6 @@ def run(self, sym_sort=[]):
             flat_sym_modes_b = [x for xs in self.irreps_b for x in xs]
             print(flat_sym_modes_b)
             del_list = []
-            # for i in range(len(self.irreps_b)):
-            # if len(self.irreps_b[i]) == 1:
-            # del_list.append(self.irreps_b[i][0])
             flat_sym_modes_b = np.delete(np.array(flat_sym_modes_b), del_list)
             ted_b = ted_b.T
             ted_b = ted_b[flat_sym_modes_b]
@@ -443,7 +451,6 @@ def run(self, sym_sort=[]):
                     "No extra computations, running back through level B single points"
                 )
                 os.chdir("DispsB")
-                print(os.getcwd())
                 reap_obj_c = Reap(
                     self.options,
                     len(eigs_b),
@@ -506,7 +513,12 @@ def run(self, sym_sort=[]):
         s_vec = SVectors(
             self.zmat_obj, self.options, self.zmat_obj.variable_dictionary_a
         )
-        s_vec.run(self.zmat_obj.cartesians_a, False, proj=self.TED_obj.proj)
+        s_vec.run(
+            self.zmat_obj.cartesians_a,
+            False,
+            proj=self.TED_obj.proj,
+            second_order=self.options.second_order,
+        )
 
         transf_disp = TransfDisp(
             s_vec,
@@ -526,18 +538,13 @@ def run(self, sym_sort=[]):
             raise RuntimeError
 
         # The displacements have been generated, now we have to run them!
-        prog = self.options.program_a
-        progname = prog.split("@")[0]
+        prog_a = self.options.program_a
+        prog_name_a = prog_a.split("@")[0]
         if self.options.gen_disps_a:
-            if os.path.exists(rootdir + "DispsA"):
-                if os.path.exists(rootdir + "/oldDispsA"):
-                    shutil.rmtree(rootdir + "/oldDispsA")
-                shutil.copytree(rootdir + "/DispsA", rootdir + "/oldDispsA")
-                shutil.rmtree(rootdir + "/DispsA")
 
             ref_geom_a = transf_disp.disp_cart["ref"]
             dir_obj_a = DirectoryTree(
-                progname,
+                prog_name_a,
                 self.zmat_obj,
                 ref_geom_a,
                 cma_level,
@@ -566,38 +573,8 @@ def run(self, sym_sort=[]):
         os.chdir(rootdir + "/DispsA")
 
         if self.options.calc_a:
-            disp_list = []
-            for i in os.listdir(rootdir + "/DispsA"):
-                disp_list.append(i)
-
-            # Generates the submit script for the displacements.
-            if self.options.cluster != "sapelo":
-                v_template = VulcanTemplate(
-                    self.options, len(disp_list), progname, prog
-                )
-                out = v_template.run()
-                with open("displacements.sh", "w") as file:
-                    file.write(out)
-
-                # Submits an array, then checks if all jobs have finished every
-                # 10 seconds.
-                sub = Submit(disp_list, self.options)
-                sub.run()
-            else:
-                s_template = SapeloTemplate(
-                    self.options, len(disp_list), progname, prog
-                )
-                out = s_template.run()
-                with open("optstep.sh", "w") as file:
-                    file.write(out)
-                for z in range(0, len(disp_list)):
-                    source = os.getcwd() + "/optstep.sh"
-                    os.chdir("./" + str(z + 1))
-                    destination = os.getcwd()
-                    shutil.copy2(source, destination)
-                    os.chdir("../")
-                sub = Submit(disp_list, self.options)
-                sub.run()
+            sub = Submit(self.options, cma_level, rootdir, prog_name_a, prog_a)
+            sub.run()
 
         # After this point, all of the jobs have finished, and it's time
         # to reap the energies as well as checking for sucesses
@@ -626,15 +603,13 @@ def run(self, sym_sort=[]):
         fc_a.run()
 
         np.set_printoptions(precision=4, linewidth=240)
-        print("Computed Force Constants:")
-        print(fc_a.FC)
+        # print("Computed Force Constants:")
+        # print(fc_a.FC)
+
+        # print("Computed Normal Mode Gradient:")
+        # print(fc_a.gradient)
 
         self.F = fc_a.FC
-        # f_b = np.dot(
-        # np.dot(LA.inv(transf_disp.eig_inv).T, f_b), LA.inv(transf_disp.eig_inv)
-        # )
-        # print("Level B FC in same basis:")
-        # print(f_b)
 
         # Recompute the G-matrix with the new geometry, and then transform
         # the G-matrix using the lower level of theory eigenvalue matrix.
@@ -693,12 +668,9 @@ def run(self, sym_sort=[]):
         # coordinates and writes out an "output.default.hess" file, which
         # is of the same format as FCMFINAL of CFOUR.
 
-        # To be implemented: logic for when the final force constants are computed
-        # atop a non-stationary point.
-
         self.F = np.dot(np.dot(transf_disp.eig_inv.T, self.F), transf_disp.eig_inv)
-        if self.options.coords != "ZMAT":
-            self.F = np.dot(self.TED_obj.proj, np.dot(self.F, self.TED_obj.proj.T))
+        self.gradient = np.dot(fc_a.gradient, transf_disp.eig_inv)
+
         cart_conv = FcConv(
             self.F,
             s_vec,
@@ -708,7 +680,13 @@ def run(self, sym_sort=[]):
             self.TED_obj,
             self.options,
         )
-        cart_conv.run()
+        if self.options.second_order:
+            cart_conv.run(grad=self.gradient)
+            self.grad_cart = cart_conv.grad
+        else:
+            cart_conv.run()
+
+        self.F_cart = cart_conv.F
 
         # if mol2.size != 0:
         #    if self.options.rmsd:
diff --git a/concordantmodes/directory_tree.py b/concordantmodes/directory_tree.py
index b931390..5cf3212 100644
--- a/concordantmodes/directory_tree.py
+++ b/concordantmodes/directory_tree.py
@@ -74,6 +74,11 @@ def run(self):
         data_buff = data.copy()
         if os.path.exists(os.getcwd() + "/old" + self.dir_name):
             shutil.rmtree("old" + self.dir_name, ignore_errors=True)
+        if os.path.exists(os.getcwd() + "/" + self.dir_name):
+            shutil.copytree(
+                os.getcwd() + "/" + self.dir_name, os.getcwd() + "/old" + self.dir_name
+            )
+            shutil.rmtree(os.getcwd() + "/" + self.dir_name)
         inp = ""
         if self.prog_name == "cfour":
             inp = "ZMAT"
diff --git a/concordantmodes/f_convert.py b/concordantmodes/f_convert.py
index e3683a5..4cc2a58 100644
--- a/concordantmodes/f_convert.py
+++ b/concordantmodes/f_convert.py
@@ -36,36 +36,67 @@ def run(self, grad=np.array([])):
                 B = self.s_vec.B
             else:
                 B = np.dot(self.ted.proj.T, self.s_vec.B)
-            BT = B.T
-            G = np.dot(B, BT)
+            G = np.dot(B, B.T)
             self.A_T = np.dot(LA.inv(G), B)
-            if self.ted.proj is not None:
-                self.A_T = np.dot(self.ted.proj, self.A_T)
             if self.options.units == "MdyneAng":
                 self.F /= self.BOHR_ANG
                 self.F *= self.MDYNE_HART
             self.F = np.einsum("pi,rj,ij->pr", self.A_T, self.A_T, self.F)
             # Non-stationary, gradient correction to internal coordinate force constants
             V2 = self.F.copy() * 0
-            if len(grad) and self.options.second_order:
+            if len(grad) and self.options.second_order and self.options.cart_fc_b:
                 # Note: I may want to use projected A_T to reduce the size of the
                 # internal coordinate basis, thus speeding up the computation.
                 # The basis will eventually need to be projected anyways.
                 np.set_printoptions(precision=6, linewidth=240)
                 self.v_q = np.dot(self.A_T, grad)
-                C2 = np.einsum("rij,pi,qj->rpq", self.s_vec.B2, self.A_T, self.A_T)
+                if self.ted.proj is None:
+                    B2 = self.s_vec.B2
+                else:
+                    B2 = np.einsum("rp,pij->rij", self.ted.proj.T, self.s_vec.B2)
+                C2 = np.einsum("rij,pi,qj->rpq", B2, self.A_T, self.A_T)
                 V2 = np.einsum("q,qpr->pr", self.v_q, C2)
 
-            self.F = self.F - V2
+                grad = np.dot(grad, self.A_T.T)
+
+            self.F -= V2
 
             if self.print_f:
-                self.print_const(fc_name="fc_int.dat")
+                self.print_const(fc_name="fc_int.dat", grad=grad)
         elif self.coord.lower() == "cartesian":
-            self.F = np.einsum("pi,rj,pr->ij", self.s_vec.B, self.s_vec.B, self.F)
+
+            if self.ted.proj is None:
+                B = self.s_vec.B
+            else:
+                B = np.dot(self.ted.proj.T, self.s_vec.B)
+
+            self.F = np.einsum("pi,rj,pr->ij", B, B, self.F)
+
+            V2 = self.F.copy() * 0
+
+            if len(grad) and self.options.second_order:
+                if self.ted.proj is None:
+                    B2 = self.s_vec.B2
+                else:
+                    B2 = np.einsum("rp,pij->rij", self.ted.proj.T, self.s_vec.B2)
+                V2 = np.einsum("rij,r->ij", B2, grad)
+
+                grad = np.dot(grad, B)
+                self.grad = grad
+
+            self.F += V2
+
+            # print("The Bs:")
+            # print(self.s_vec.B.shape)
+            # print(self.s_vec.B)
+            # print(B.shape)
+            # print(B)
+            # raise RuntimeError
+
             if self.print_f:
-                self.print_const()
+                self.print_const(grad=grad)
 
-    def print_const(self, fc_name="fc_a.dat"):
+    def print_const(self, fc_name="fc_a.dat", grad=np.array([])):
         self.N = len(self.F)
         fc_output = ""
         fc_output += "{:5d}{:5d}\n".format(len(self.zmat.atom_list), self.N)
@@ -85,3 +116,22 @@ def print_const(self, fc_name="fc_a.dat"):
             fc_output += "\n"
         with open(fc_name, "w+") as file:
             file.write(fc_output)
+
+        if len(grad):
+            g_print = grad
+            # g_print = g_print.flatten()
+            gr_output = ""
+            for i in range(len(g_print) // 3):
+                gr_output += "{:20.10f}".format(g_print[3 * i])
+                gr_output += "{:20.10f}".format(g_print[3 * i + 1])
+                gr_output += "{:20.10f}".format(g_print[3 * i + 2])
+                gr_output += "\n"
+            if len(g_print) % 3:
+                for i in range(len(g_print) % 3):
+                    gr_output += "{:20.10f}".format(
+                        g_print[3 * (len(g_print) // 3) + i]
+                    )
+                fc_output += "\n"
+            # Introduce grad_name
+            with open("fc_a.grad", "w+") as file:
+                file.write(gr_output)
diff --git a/concordantmodes/options.py b/concordantmodes/options.py
index b19184e..1cfce62 100644
--- a/concordantmodes/options.py
+++ b/concordantmodes/options.py
@@ -63,7 +63,7 @@ def __init__(self, **kwargs):
         self.success_regex_c = kwargs.pop("success_regex_c", "")
         self.molsym_symmetry = kwargs.pop("molsym_symmetry", False)
         self.tight_disp = kwargs.pop("tight_disp", False)
-        self.time_limit = kwargs.pop("time_limit", "4:00:00")
+        self.time_limit = kwargs.pop("time_limit", "10:00:00")
         self.tol = kwargs.pop("tol", 1.0e-14)
         self.topo_analysis = kwargs.pop("topo_analysis", False)
         self.topo_max_it = kwargs.pop("topo_max_it", 20)
diff --git a/concordantmodes/reap.py b/concordantmodes/reap.py
index c57fd9b..47b78f6 100644
--- a/concordantmodes/reap.py
+++ b/concordantmodes/reap.py
@@ -56,8 +56,11 @@ def run(self):
             grad_regex2 = re.compile(self.gradient_regex[1])
         size = self.num_deg_free
 
-        if self.options.second_order:
-            size = self.indices[-1][0] + 1
+        # if self.options.second_order:
+        # print(self.indices)
+        # print(self.indices[-1])
+        # raise RuntimeError
+        # size = self.indices[-1][0] + 1
 
         self.fail_list = []
         if not self.deriv_level:
@@ -98,6 +101,9 @@ def run(self):
 
             Sum = 1
             h = 0
+            # print(p_en_array.shape)
+            # print(p_en_array)
+            # print(indices)
             for index in indices:
                 i, j = index[0], index[1]
                 p_en_array[i, j] = energy = self.reap_energies(
diff --git a/concordantmodes/s_vectors.py b/concordantmodes/s_vectors.py
index bc0c940..9d58a12 100644
--- a/concordantmodes/s_vectors.py
+++ b/concordantmodes/s_vectors.py
@@ -563,6 +563,7 @@ def run(self, carts, B_proj, proj=None, second_order=False):
 
         # Run numerical second order B-tensor here.
         if second_order:
+            B = self.B.copy()
             Proj = self.proj.copy()
             self.B2 = self.second_order_B()
             self.B2[np.abs(self.B2) < 1e-14] = 0
@@ -578,6 +579,11 @@ def run(self, carts, B_proj, proj=None, second_order=False):
                         self.B2[i, k, j + 1] = self.B2[i, j + 1, k]
             self.B2 = self.B2.astype(float)
             self.proj = Proj
+            # np.set_printoptions(precision=6, linewidth=240)
+            # print(self.B)
+            self.B = B
+            # print(self.B)
+            # raise RuntimeError
 
     def compute_STRE(self, x1, x2):
         s = (x1 - x2) / self.compute_r(x1, x2)
@@ -778,15 +784,15 @@ def second_order_B(self):
 
         B_disp.run()
         self.run(B_disp.p_disp[0], False, proj=self.proj, second_order=False)
-        B_list_p = np.array([self.B], dtype=object)
+        B_list_p = np.array([self.B.copy()], dtype=object)
         self.run(B_disp.m_disp[0], False, proj=self.proj, second_order=False)
-        B_list_m = np.array([self.B], dtype=object)
+        B_list_m = np.array([self.B.copy()], dtype=object)
 
         for i in range(len(B_disp.p_disp) - 1):
             self.run(B_disp.p_disp[i + 1], False, second_order=False)
-            B_list_p = np.append(B_list_p, [self.B], axis=0)
+            B_list_p = np.append(B_list_p, [self.B.copy()], axis=0)
             self.run(B_disp.m_disp[i + 1], False, second_order=False)
-            B_list_m = np.append(B_list_m, [self.B], axis=0)
+            B_list_m = np.append(B_list_m, [self.B.copy()], axis=0)
 
         # And differentiate
         B2 = self.num_differentiate(B_list_p, B_list_m)
diff --git a/concordantmodes/submit.py b/concordantmodes/submit.py
index 317722e..698d4b4 100644
--- a/concordantmodes/submit.py
+++ b/concordantmodes/submit.py
@@ -1,17 +1,35 @@
 import re
+import shutil
 import subprocess
 import time
 import os
 from subprocess import Popen
+from concordantmodes.vulcan_template import VulcanTemplate
+from concordantmodes.sapelo_template import SapeloTemplate
 
 
 class Submit(object):
-    def __init__(self, disp_list, options):
-        self.disp_list = disp_list
+    def __init__(self, options, cma_level, rootdir, prog_name, prog):
+        self.cma_level = cma_level
         self.options = options
+        self.prog = prog
+        self.prog_name = prog_name
+        self.rootdir = rootdir
 
     def run(self):
+        disp_list = []
+        for i in os.listdir(self.rootdir + "/Disps" + self.cma_level):
+            disp_list.append(i)
+
         if self.options.cluster != "sapelo":
+            v_template = VulcanTemplate(
+                self.options, len(disp_list), self.prog_name, self.prog
+            )
+            out = v_template.run()
+
+            with open("displacements.sh", "w") as file:
+                file.write(out)
+
             pipe = subprocess.PIPE
 
             process = subprocess.run(
@@ -27,7 +45,7 @@ def run(self):
                 )
                 qacct_string = str(qacct_proc.stdout)
                 job_match = re.findall(self.job_fin_regex, qacct_string)
-                if len(job_match) == len(self.disp_list):
+                if len(job_match) == len(disp_list):
                     break
                 time.sleep(30)
 
@@ -36,11 +54,23 @@ def run(self):
             pass
 
         else:
-            from subprocess import Popen
+            s_template = SapeloTemplate(
+                self.options, len(disp_list), self.prog_name, self.prog
+            )
+            out = s_template.run()
+
+            with open("optstep.sh", "w") as file:
+                file.write(out)
+            for z in range(0, len(disp_list)):
+                source = os.getcwd() + "/optstep.sh"
+                os.chdir("./" + str(z + 1))
+                destination = os.getcwd()
+                shutil.copy2(source, destination)
+                os.chdir("../")
 
             processes = []
             print(os.getcwd())
-            for z in range(len(self.disp_list)):
+            for z in range(len(disp_list)):
                 path = str(z + 1) + "/"
                 pipe = subprocess.PIPE
                 job = subprocess.run(
@@ -49,8 +79,9 @@ def run(self):
                 processes.append(job)
                 time.sleep(3)
 
-            print("First Nap")
-            time.sleep(10)
+            # print("First Nap")
+            # print("Napping")
+            # time.sleep(10)
             for q in range(len(processes)):
                 while True:
                     job = processes[q]
@@ -70,5 +101,6 @@ def run(self):
                             + str(q)
                         )
                         break
-            print("Second Nap")
-            time.sleep(10)
+            # print("Second Nap")
+            print("Napping")
+            time.sleep(15)
diff --git a/concordantmodes/symmetry.py b/concordantmodes/symmetry.py
index 37a4d8a..e4450f6 100644
--- a/concordantmodes/symmetry.py
+++ b/concordantmodes/symmetry.py
@@ -22,6 +22,7 @@ def __init__(self, zmat, options, proj):
         self.zmat = zmat
         self.options = options
         self.proj = proj
+        self.proj_irreps = None
 
     def dummy_obj(self):
         pass
@@ -315,21 +316,25 @@ def mode_symmetry_sort(self, TED, sym_sort, freqs):
         sym_freqs = copy.deepcopy(sym_modes)
         del_list = []
         for i in range(len(sym_modes)):
-            if len(sym_modes[i]) > 1:
-                for j in range(len(sym_modes[i])):
-                    index = sym_modes[i][j]
-                    sym_freqs[i][j] = freqs[index].copy()
-                sym_freqs[i].reverse()
-            elif len(sym_modes[i]) == 1:
-                del_list.append(i)
-            else:
-                pass
+            # if len(sym_modes[i]) > 1:
+            # for j in range(len(sym_modes[i])):
+            # index = sym_modes[i][j]
+            # sym_freqs[i][j] = freqs[index].copy()
+            # sym_freqs[i].reverse()
+            # elif len(sym_modes[i]) == 1:
+            # del_list.append(i)
+            # else:
+            # pass
+            for j in range(len(sym_modes[i])):
+                index = sym_modes[i][j]
+                sym_freqs[i][j] = freqs[index].copy()
+            sym_freqs[i].reverse()
         del_list.reverse()
         if len(del_list):
             for i in del_list:
                 print(freqs[sym_modes[i][0]])
-        for i in del_list:
-            del sym_freqs[i]
+        # for i in del_list:
+        # del sym_freqs[i]
         flat_sym_freqs = [x for xs in sym_freqs for x in xs]
         flat_sym_freqs = np.array(flat_sym_freqs)
 
diff --git a/concordantmodes/tests/suite_execute.py b/concordantmodes/tests/suite_execute.py
index 494a5cc..4844ba9 100644
--- a/concordantmodes/tests/suite_execute.py
+++ b/concordantmodes/tests/suite_execute.py
@@ -61,7 +61,8 @@ def run(self):
         )
         self.f_conv.run()
 
-        self.F = np.dot(self.TED_obj.proj.T, np.dot(self.f_conv.F, self.TED_obj.proj))
+        self.F = self.f_conv.F
+        # self.F = np.dot(self.TED_obj.proj.T, np.dot(self.f_conv.F, self.TED_obj.proj))
         self.G = np.dot(self.TED_obj.proj.T, np.dot(self.g_mat.G, self.TED_obj.proj))
         self.GF = GFMethod(
             self.G, self.F, self.ZMAT, self.TED_obj, self.options, self.symm_obj.symtext
diff --git a/concordantmodes/tests/test_f_convert.py b/concordantmodes/tests/test_f_convert.py
index 42d3a7c..0fa40e9 100644
--- a/concordantmodes/tests/test_f_convert.py
+++ b/concordantmodes/tests/test_f_convert.py
@@ -32,6 +32,7 @@ def test_f_convert2int():
         suite.options,
     )
     FCint.run()
+    FCint.F = np.dot(suite.TED_obj.proj, np.dot(FCint.F, suite.TED_obj.proj.T))
 
     FCintR = FcRead(suite.path + "/fc_int.dat")
     FCintR.run()
diff --git a/concordantmodes/tests/_test_transf_disp.py b/concordantmodes/tests/test_transf_disp.py
similarity index 89%
rename from concordantmodes/tests/_test_transf_disp.py
rename to concordantmodes/tests/test_transf_disp.py
index 0a7bd07..20fc270 100644
--- a/concordantmodes/tests/_test_transf_disp.py
+++ b/concordantmodes/tests/test_transf_disp.py
@@ -49,7 +49,8 @@ def test_transf_disp():
 
     os.chdir(suite.root)
 
-    assert np.allclose(np.asarray(disp_ref), disps.p_disp[3][3], rtol=0.0, atol=1e-5)
+    assert np.allclose(np.asarray(disp_ref), disps.m_disp[3][3], rtol=0.0, atol=1e-5)
+    # assert np.allclose(np.asarray(disp_ref), disps.p_disp[3][3], rtol=0.0, atol=1e-5)
 
 
 test_transf_disp()
diff --git a/concordantmodes/transf_disp.py b/concordantmodes/transf_disp.py
index 9f35b27..6eb7dad 100644
--- a/concordantmodes/transf_disp.py
+++ b/concordantmodes/transf_disp.py
@@ -89,16 +89,22 @@ def run(self, fc=np.array([])):
                 u[i] = 1.0 / u[i]
         u = np.repeat(u, 3)
         u = np.diag(u)
-        self.A = self.compute_A(
-            self.B.copy(), self.ted.proj, self.eig_inv, u, cma_level=self.cma_level
-        )
+        # print(self.coord_type)
+        # raise RuntimeError
+        if self.coord_type != "cartesian":
+            # print(self.coord_type)
+            self.A = self.compute_A(
+                self.B.copy(), self.ted.proj, self.eig_inv, u, cma_level=self.cma_level
+            )
 
-        # Generate the normal coordinate values for the reference structure
-        self.n_coord = self.int_c(self.ref_carts, self.eig_inv, self.ted.proj)
+            # Generate the normal coordinate values for the reference structure
+            self.n_coord = self.int_c(self.ref_carts, self.eig_inv, self.ted.proj)
 
-        print("Normal Coordinate Values:")
-        for i in range(len(self.n_coord)):
-            print("Normal Coordinate #{:<4n}: {: 3.5f}".format(i + 1, self.n_coord[i]))
+            print("Normal Coordinate Values:")
+            for i in range(len(self.n_coord)):
+                print(
+                    "Normal Coordinate #{:<4n}: {: 3.5f}".format(i + 1, self.n_coord[i])
+                )
 
         # Next, we will have to specify our Normal mode internal coordinate
         # displacement sizes
@@ -311,6 +317,8 @@ def run(self, fc=np.array([])):
                         (len(self.ref_carts.flatten()), len(self.ref_carts.flatten())),
                         dtype=object,
                     )
+                    # print(self.indices)
+                    # raise RuntimeError
                     for index in self.indices:
                         i, j = index[0], index[1]
                         if i == j:
diff --git a/examples/1_H2O/molpro/CMA0A/No_Spec_Deloc/zmat b/examples/1_H2O/molpro/CMA0A/No_Spec_Deloc/zmat
index b3d3a44..62f62a7 100644
--- a/examples/1_H2O/molpro/CMA0A/No_Spec_Deloc/zmat
+++ b/examples/1_H2O/molpro/CMA0A/No_Spec_Deloc/zmat
@@ -1,8 +1,3 @@
-ZMAT begin
-1 2
-2 3
-ZMAT end
-
 cart begin
      H        -0.00000000    -1.42462540     0.99592408
      O         0.00000000     0.00000000    -0.12550454
diff --git a/examples/6_Psi4_script/Psi4_CMA.py b/examples/6_Psi4_script/Psi4_CMA.py
new file mode 100644
index 0000000..d982f0c
--- /dev/null
+++ b/examples/6_Psi4_script/Psi4_CMA.py
@@ -0,0 +1,392 @@
+import psi4
+import concordantmodes
+import time
+import numpy as np
+from numpy import linalg as LA
+import os
+import re
+import sys
+import shutil
+import subprocess
+from concordantmodes.algorithm import Algorithm
+from concordantmodes.f_convert import FcConv
+from concordantmodes.f_read import FcRead
+from concordantmodes.force_constant import ForceConstant
+from concordantmodes.gf_method import GFMethod
+from concordantmodes.g_matrix import GMatrix
+from concordantmodes.g_read import GrRead
+from concordantmodes import masses
+from concordantmodes.molden_writer import MoldenWriter
+from concordantmodes.options import Options
+from concordantmodes.s_vectors import SVectors
+from concordantmodes.symmetry import Symmetry
+from concordantmodes.ted import TED
+from concordantmodes.transf_disp import TransfDisp
+from concordantmodes.zmat import Zmat
+
+t1 = time.time()
+
+# First things first. Let's get a Psi4 hessian running.
+print(psi4.__version__)
+
+psi4.set_output_file("output.dat", False)
+# psi4.core.be_quiet()
+
+# Set psi4 options here
+psi4.set_memory("10 GB")
+psi4.set_options({"freeze_core": True})
+psi4.set_options({"normal_modes_write": True})
+psi4.set_options({"e_convergence": 10})
+psi4.set_options({"d_convergence": 10})
+psi4.set_options({"r_convergence": 8})
+psi4.set_options({"cc_type": "conv"})
+# psi4.set_options({'scf_type': 'pk'})
+
+mol = psi4.geometry(
+    """
+     C         0.00000000     0.00000000     1.14423616
+     H         0.00000000    -1.76987484     2.24620597
+     H         0.00000000     1.76987484     2.24620597
+     O        -0.00000000     0.00000000    -1.14151276
+        """
+)
+
+# psi4 units set as bohr, if not then CMA units will need to be set as angstrom
+mol.set_units(psi4.core.GeometryUnits.Bohr)
+mol.fix_orientation(True)
+mol.fix_com(True)
+
+# Set the CMA variables here:
+level_a = "ccsd(t)/cc-pvtz"
+level_b = "mp2/cc-pvtz"
+
+options_kwargs = {
+    "man_proj": False,
+    "coords": "Delocalized",
+    "reduced_disp": True,
+    "cart_fc_b": True,
+    "covalent_radii": True,
+    "second_order": True,
+    # "reduced_disp"           : 0.002,
+    # "molsym_symmetry" : True,
+    # "autosalcs" : True,
+}
+
+options = Options(**options_kwargs)
+print("Disp:")
+print(options.disp)
+
+# Change this if man_proj=True
+proj = None
+# modify this for custom symmetry
+sym_sort = []
+
+
+# string value to access "state" of CMA program. Are we in level A, B, C, etc...
+cma_level = "B"
+
+# Read in Psi4 Molecule information into our zmat object.
+zmat_obj = Zmat(options)
+zmat_obj.divider = False
+
+zmat_obj.cartesians_b = mol.geometry().to_array()
+zmat_obj.cartesians_a = mol.geometry().to_array()
+
+atom_list = []
+for i in range(len(zmat_obj.cartesians_a)):
+    atom_list.append(mol.fsymbol(i))
+zmat_obj.atom_list = atom_list
+
+zmat_obj.masses = [masses.get_mass(label) for label in zmat_obj.atom_list]
+for i in range(len(zmat_obj.masses)):
+    zmat_obj.masses[i] = zmat_obj.masses[i] / zmat_obj.amu_elMass
+zmat_obj.mass_weight = np.diag(np.array(zmat_obj.masses).repeat(3))
+
+zmat_obj.zmat_process("")
+zmat_obj.zmat_calc()
+zmat_obj.zmat_compile()
+zmat_obj.zmat_print()
+
+
+# Do we want to use molsym_symmetry or "manual" symmetry via sym_sort?
+symm_obj = Symmetry(zmat_obj, options, proj)
+if options.molsym_symmetry:
+    symm_obj.run()
+else:
+    """
+    We won't run the symmetry code, but we'll create a dummy object to be passed as an argument.
+    #TODO: This is a hacky way to do this, but it's a quick fix for now. Maybe reincorporate symmetry as a s_vector obj?
+    """
+    symm_obj.dummy_obj()
+    symm_obj.symtext = None
+    # check if sym_sort object was passed in. If so, intialize sym_sort objects
+    if len(sym_sort) > 1:
+        symm_obj.create_flat_sym_sort(sym_sort)
+
+
+# Compute the initial s-vectors
+s_vec = SVectors(zmat_obj, options, zmat_obj.variable_dictionary_b)
+s_vec.run(
+    zmat_obj.cartesians_b,
+    True,
+    second_order=options.second_order,
+)
+
+if options.molsym_symmetry:
+    symm_obj.make_proj(s_vec)
+    s_vec.proj = copy.deepcopy(symm_obj.salc_proj)
+
+# Define the Total Energy Distribution object
+TED_obj = TED(s_vec.proj, zmat_obj, options)
+
+# Print out the percentage composition of the projected coordinates
+if options.coords != "ZMAT":
+    TED_obj.run(np.eye(len(TED_obj.proj.T)), np.zeros(len(TED_obj.proj.T)))
+
+# Compute G-Matrix
+g_mat = GMatrix(zmat_obj, s_vec, options)
+g_mat.run()
+
+psi4.energy(level_b)
+
+np.set_printoptions(precision=4, linewidth=240)
+
+grad_e, grad_wfn = psi4.gradient(level_b, molecule=mol, return_wfn=True)
+# print(grad_wfn.gradient().to_array().flatten())
+cart_grad = grad_wfn.gradient().to_array().flatten()
+
+if options.second_order:
+    g_read_obj = GrRead("")
+    g_read_obj.cart_grad = cart_grad
+
+freq_e, freq_wfn = psi4.frequency(level_b, molecule=mol, return_wfn=True)
+print(freq_wfn.hessian().to_array())
+cart_hess = freq_wfn.hessian().to_array()
+
+# raise RuntimeError
+
+# This might not be necessary, just holding onto it for now
+options.deriv_level_b = 0
+options.deriv_level_a = 0
+
+f_conv_obj = FcConv(
+    cart_hess,
+    # f_read_obj.fc_mat,
+    s_vec,
+    zmat_obj,
+    "internal",
+    False,
+    TED_obj,
+    options,
+)
+if options.second_order:
+    f_conv_obj.run(grad=g_read_obj.cart_grad)
+else:
+    f_conv_obj.run()
+
+F = f_conv_obj.F
+
+g_mat.G = np.dot(TED_obj.proj.T, np.dot(g_mat.G, TED_obj.proj))
+
+options.second_order = False
+# Numerical trimming
+F[np.abs(F) < options.tol] = 0
+g_mat.G[np.abs(g_mat.G) < options.tol] = 0
+
+# Run the GF matrix method with the internal F-Matrix and computed G-Matrix!
+print("Initial Frequencies:")
+b_GF = GFMethod(
+    g_mat.G.copy(),
+    F.copy(),
+    zmat_obj,
+    TED_obj,
+    options,
+    symm_obj.symtext,
+    cma="init",
+)
+b_GF.run()
+
+# raise RuntimeError
+
+ted_b = b_GF.ted.TED
+if len(sym_sort):
+    irreps_b, flat_sym_freqs = symm_obj.mode_symmetry_sort(
+        b_GF.ted.TED, sym_sort, b_GF.freq
+    )
+    ref_b = np.array(flat_sym_freqs)
+    #### this block could probably be moved inside the symmetry.py module?
+    flat_sym_modes_b = [x for xs in irreps_b for x in xs]
+    print(flat_sym_modes_b)
+    del_list = []
+    for i in range(len(irreps_b)):
+        if len(irreps_b[i]) == 1:
+            del_list.append(irreps_b[i][0])
+    flat_sym_modes_b = np.delete(np.array(flat_sym_modes_b), del_list2)
+    ted_b = ted_b.T
+    ted_b = ted_b[flat_sym_modes_b]
+    ted_b = ted_b.T
+    #### end of block that could probably be moved inside the symmetry.py module?
+
+
+# Now for the TED check.
+G = np.dot(np.dot(LA.inv(b_GF.L), g_mat.G), LA.inv(b_GF.L).T)
+G[np.abs(G) < options.tol] = 0
+F = np.dot(np.dot(b_GF.L.T, F), b_GF.L)
+F[np.abs(F) < options.tol] = 0
+
+print("TED Frequencies:")
+TED_GF = GFMethod(
+    G,
+    F,
+    zmat_obj,
+    TED_obj,
+    options,
+    symm_obj.symtext,
+    cma=False,
+)
+TED_GF.run()
+
+initial_fc = TED_GF.eig_v
+eigs = len(TED_GF.L)
+
+
+cma_level = "A"
+num_deg_free = s_vec.proj.shape[1]
+if options.molsym_symmetry:
+    algo = Algorithm(num_deg_free, cma_level, options, symm_obj.proj_irreps)
+else:
+    algo = Algorithm(num_deg_free, cma_level, options, None)
+# Only generate CMA-0 indices for now
+algo.run()
+if options.molsym_symmetry:
+    symm_obj.indices_by_irrep = algo.indices_by_irrep
+
+transf_disp = TransfDisp(
+    s_vec, zmat_obj, b_GF.L, True, TED_obj, options, algo.indices, symm_obj=symm_obj
+)
+transf_disp.run(fc=F)
+p_disp = transf_disp.p_disp
+m_disp = transf_disp.m_disp
+
+print("Computing Ref Energy:")
+sub_t1 = time.time()
+ref_en, first_wfn = psi4.energy(level_a, molecule=mol, return_wfn=True)
+print(ref_en)
+psi4.core.clean()
+sub_t2 = time.time()
+print("This energy took " + str(sub_t2 - sub_t1) + " seconds to run.")
+
+mol_disp = mol.clone()
+
+p_en_array = []
+print("Positive disps ({length:d} total):".format(length=len(p_disp)))
+for i in range(len(p_disp)):
+    sub_t1 = time.time()
+    disp = p_disp[i, i]
+    disp_shape = disp.shape
+    a, b = disp_shape[0], disp_shape[1]
+    psi_disp = psi4.core.Matrix(a, b).from_array(disp)  # initialize psi4.Matrix
+    mol_disp.set_geometry(psi_disp)
+    e = psi4.energy(level_a, molecule=mol_disp)
+    p_en_array.append(e)
+    psi4.core.clean()
+    print(i)
+    print(e)
+    sub_t2 = time.time()
+    print("This energy took " + str(sub_t2 - sub_t1) + " seconds to run.")
+
+p_en_array = np.diag(p_en_array)
+print(p_en_array)
+
+m_en_array = []
+print("Negative disps ({length:d} total):".format(length=len(p_disp)))
+for i in range(len(m_disp)):
+    sub_t1 = time.time()
+    disp = m_disp[i, i]
+    disp_shape = disp.shape
+    a, b = disp_shape[0], disp_shape[1]
+    psi_disp = psi4.core.Matrix(a, b).from_array(disp)  # initialize psi4.Matrix
+    mol_disp.set_geometry(psi_disp)
+    e = psi4.energy(level_a, molecule=mol_disp)
+    m_en_array.append(e)
+    psi4.core.clean()
+    print(i)
+    print(e)
+    sub_t2 = time.time()
+    print("This energy took " + str(sub_t2 - sub_t1) + " seconds to run.")
+
+m_en_array = np.diag(m_en_array)
+print(m_en_array)
+
+f_c = ForceConstant(transf_disp, p_en_array, m_en_array, ref_en, options, algo.indices)
+f_c.run()
+print("Computed Force Constants:")
+print(f_c.FC)
+F = f_c.FC
+
+g_mat = GMatrix(zmat_obj, s_vec, options)
+g_mat.run()
+g_mat.G = np.dot(TED_obj.proj.T, np.dot(g_mat.G, TED_obj.proj))
+
+g_mat.G = np.dot(np.dot(transf_disp.eig_inv, g_mat.G), transf_disp.eig_inv.T)
+g_mat.G[np.abs(g_mat.G) < options.tol] = 0
+G = g_mat.G
+
+# Final GF Matrix run
+print("Final Harmonic Frequencies:")
+a_GF = GFMethod(
+    G,
+    F,
+    zmat_obj,
+    TED_obj,
+    options,
+    symm_obj.symtext,
+)
+a_GF.run()
+
+
+print("////////////////////////////////////////////")
+print("//{:^40s}//".format(" Final TED"))
+print("////////////////////////////////////////////")
+TED_obj.run(np.dot(b_GF.L, a_GF.L), a_GF.freq)
+
+# This code prints out the frequencies in order of energy as well
+# as the ZPVE in several different units.
+print(
+    "Final Harmonic ZPVE in: "
+    + "{:6.2f}".format(np.sum(a_GF.freq) / 2)
+    + " (cm^-1) "
+    + "{:6.2f}".format(0.5 * np.sum(a_GF.freq) / 349.7550881133)
+    + " (kcal mol^-1) "
+    + "{:6.2f}".format(0.5 * np.sum(a_GF.freq) / 219474.6313708)
+    + " (hartrees) "
+)
+
+F = np.dot(np.dot(transf_disp.eig_inv.T, F), transf_disp.eig_inv)
+if options.coords != "ZMAT":
+    F = np.dot(TED_obj.proj, np.dot(F, TED_obj.proj.T))
+options.second_order = False
+cart_conv = FcConv(
+    F,
+    s_vec,
+    zmat_obj,
+    "cartesian",
+    True,
+    TED_obj,
+    options,
+)
+cart_conv.run()
+
+molden = MoldenWriter(zmat_obj, transf_disp, a_GF.freq)
+molden.run()
+
+t2 = time.time()
+print("This program took " + str(t2 - t1) + " seconds to run.")
+
+# Bonus, run the psi4 TZ frequencies to see how well we did
+t1 = time.time()
+psi4_e, psi4_wfn = psi4.frequencies("ccsd(t)/cc-pvtz", molecule=mol, return_wfn=True)
+t2 = time.time()
+print("Psi4 freq - CMA freq:")
+print(psi4_wfn.frequencies().to_array() - a_GF.freq)
+print("Psi4 took " + str(t2 - t1) + " seconds to run.")