Merge pull request #17 from CCQC/MitchSandbox

Mitch sandbox
CCQC · Mar 20, 2023 · 4e8e77d · 4e8e77d
2 parents f5b51d6 + 9db1a94
commit 4e8e77d
Show file tree

Hide file tree

Showing 135 changed files with 7,569 additions and 207 deletions.
diff --git a/Example/2_Methanol/molpro/CMA_0B/.template.dat.swp b/Example/2_Methanol/molpro/CMA_0B/.template.dat.swp
diff --git a/Example/2_Methanol/molpro/CMA_0B/main.py b/Example/2_Methanol/molpro/CMA_0B/main.py
@@ -1,7 +1,7 @@
 from concordantmodes.options import Options
 
 options_kwargs = {
-    "queue": "gen4.q,gen6.q,debug.q",
+    "queue": "gen4.q,gen6.q",
     "program": "[email protected]+mpi",
     "energy_regex": r"\(T\) total energy\s+(\-\d+\.\d+)",
     "cart_insert": 9,

diff --git a/README.md b/README.md
@@ -37,6 +37,10 @@ User must have at least:
 Python version 3.5.4
 Numpy version 1.13.1
 
+Developers must have:
+Pytest 7.2.0 (-c condaforge)
+Pytest-xdist 3.0.2 (-c condaforge)
+
 A simple way to ensure this program may be run is by creating a conda environment with the following command:
 
 conda create --name CMA python=3.9 sympy numpy scipy

diff --git a/__pycache__/__init__.cpython-35.pyc b/__pycache__/__init__.cpython-35.pyc
diff --git a/__pycache__/__init__.cpython-38.pyc b/__pycache__/__init__.cpython-38.pyc
diff --git a/__pycache__/__init__.cpython-39.pyc b/__pycache__/__init__.cpython-39.pyc
diff --git a/__pycache__/algorithm.cpython-35.pyc b/__pycache__/algorithm.cpython-35.pyc
diff --git a/__pycache__/algorithm.cpython-38.pyc b/__pycache__/algorithm.cpython-38.pyc
diff --git a/__pycache__/algorithm.cpython-39.pyc b/__pycache__/algorithm.cpython-39.pyc
diff --git a/__pycache__/cma.cpython-35.pyc b/__pycache__/cma.cpython-35.pyc
diff --git a/__pycache__/cma.cpython-38.pyc b/__pycache__/cma.cpython-38.pyc
diff --git a/__pycache__/cma.cpython-39.pyc b/__pycache__/cma.cpython-39.pyc
diff --git a/__pycache__/directory_tree.cpython-35.pyc b/__pycache__/directory_tree.cpython-35.pyc
diff --git a/__pycache__/directory_tree.cpython-38.pyc b/__pycache__/directory_tree.cpython-38.pyc
diff --git a/__pycache__/directory_tree.cpython-39.pyc b/__pycache__/directory_tree.cpython-39.pyc
diff --git a/__pycache__/f_convert.cpython-35.pyc b/__pycache__/f_convert.cpython-35.pyc
diff --git a/__pycache__/f_convert.cpython-38.pyc b/__pycache__/f_convert.cpython-38.pyc
diff --git a/__pycache__/f_convert.cpython-39.pyc b/__pycache__/f_convert.cpython-39.pyc
diff --git a/__pycache__/f_read.cpython-35.pyc b/__pycache__/f_read.cpython-35.pyc
diff --git a/__pycache__/f_read.cpython-38.pyc b/__pycache__/f_read.cpython-38.pyc
diff --git a/__pycache__/f_read.cpython-39.pyc b/__pycache__/f_read.cpython-39.pyc
diff --git a/__pycache__/force_constant.cpython-35.pyc b/__pycache__/force_constant.cpython-35.pyc
diff --git a/__pycache__/force_constant.cpython-38.pyc b/__pycache__/force_constant.cpython-38.pyc
diff --git a/__pycache__/force_constant.cpython-39.pyc b/__pycache__/force_constant.cpython-39.pyc
diff --git a/__pycache__/g_matrix.cpython-35.pyc b/__pycache__/g_matrix.cpython-35.pyc
diff --git a/__pycache__/g_matrix.cpython-38.pyc b/__pycache__/g_matrix.cpython-38.pyc
diff --git a/__pycache__/g_matrix.cpython-39.pyc b/__pycache__/g_matrix.cpython-39.pyc
diff --git a/__pycache__/gf_method.cpython-35.pyc b/__pycache__/gf_method.cpython-35.pyc
diff --git a/__pycache__/gf_method.cpython-38.pyc b/__pycache__/gf_method.cpython-38.pyc
diff --git a/__pycache__/gf_method.cpython-39.pyc b/__pycache__/gf_method.cpython-39.pyc
diff --git a/__pycache__/int2cart.cpython-35.pyc b/__pycache__/int2cart.cpython-35.pyc
diff --git a/__pycache__/int2cart.cpython-38.pyc b/__pycache__/int2cart.cpython-38.pyc
diff --git a/__pycache__/int2cart.cpython-39.pyc b/__pycache__/int2cart.cpython-39.pyc
diff --git a/__pycache__/masses.cpython-35.pyc b/__pycache__/masses.cpython-35.pyc
diff --git a/__pycache__/masses.cpython-38.pyc b/__pycache__/masses.cpython-38.pyc
diff --git a/__pycache__/masses.cpython-39.pyc b/__pycache__/masses.cpython-39.pyc
diff --git a/__pycache__/options.cpython-35.pyc b/__pycache__/options.cpython-35.pyc
diff --git a/__pycache__/options.cpython-38.pyc b/__pycache__/options.cpython-38.pyc
diff --git a/__pycache__/options.cpython-39.pyc b/__pycache__/options.cpython-39.pyc
diff --git a/__pycache__/reap.cpython-35.pyc b/__pycache__/reap.cpython-35.pyc
diff --git a/__pycache__/reap.cpython-38.pyc b/__pycache__/reap.cpython-38.pyc
diff --git a/__pycache__/reap.cpython-39.pyc b/__pycache__/reap.cpython-39.pyc
diff --git a/__pycache__/rmsd.cpython-38.pyc b/__pycache__/rmsd.cpython-38.pyc
diff --git a/__pycache__/rmsd.cpython-39.pyc b/__pycache__/rmsd.cpython-39.pyc
diff --git a/__pycache__/s_vectors.cpython-38.pyc b/__pycache__/s_vectors.cpython-38.pyc
diff --git a/__pycache__/s_vectors.cpython-39.pyc b/__pycache__/s_vectors.cpython-39.pyc
diff --git a/__pycache__/sapelo_template.cpython-38.pyc b/__pycache__/sapelo_template.cpython-38.pyc
diff --git a/__pycache__/sapelo_template.cpython-39.pyc b/__pycache__/sapelo_template.cpython-39.pyc
diff --git a/__pycache__/submit.cpython-38.pyc b/__pycache__/submit.cpython-38.pyc
diff --git a/__pycache__/submit.cpython-39.pyc b/__pycache__/submit.cpython-39.pyc
diff --git a/__pycache__/ted.cpython-38.pyc b/__pycache__/ted.cpython-38.pyc
diff --git a/__pycache__/ted.cpython-39.pyc b/__pycache__/ted.cpython-39.pyc
diff --git a/__pycache__/transf_disp.cpython-38.pyc b/__pycache__/transf_disp.cpython-38.pyc
diff --git a/__pycache__/transf_disp.cpython-39.pyc b/__pycache__/transf_disp.cpython-39.pyc
diff --git a/__pycache__/vulcan_template.cpython-38.pyc b/__pycache__/vulcan_template.cpython-38.pyc
diff --git a/__pycache__/vulcan_template.cpython-39.pyc b/__pycache__/vulcan_template.cpython-39.pyc
diff --git a/__pycache__/zmat.cpython-38.pyc b/__pycache__/zmat.cpython-38.pyc
diff --git a/__pycache__/zmat.cpython-39.pyc b/__pycache__/zmat.cpython-39.pyc
diff --git a/algorithm.py b/algorithm.py
@@ -61,8 +61,8 @@ def loop(self, a, off_diag, lim):
     def coupling_diagnostic(self, a, initial_fc, tolerance):
         indices = []
         diag = np.zeros((a, a))
-        with open("S_p.npy", "rb") as x:
-            S = np.load(x)
+        # with open("S_p.npy", "rb") as x:
+            # S = np.load(x)
         for x in range(a):
             for y in range(a):
                 if x == y:
@@ -72,8 +72,8 @@ def coupling_diagnostic(self, a, initial_fc, tolerance):
         # print(diag, "printing diag")
         diag = np.absolute(diag)
         print(diag)
-        with open("D.npy", "wb") as q:
-            np.save(q, diag)
+        # with open("D.npy", "wb") as q:
+            # np.save(q, diag)
         # diag[np.abs(diag) < 1e-31] = 1e-30
         data = np.abs(diag)
         hist, bin_edges = np.histogram(data, bins=a)
@@ -85,7 +85,7 @@ def coupling_diagnostic(self, a, initial_fc, tolerance):
             if index[1] > index[0]:
                 indices_new.append(index)
         indices = indices_new
-        if self.options.clean_house:
-            os.system("rm D.npy  S_p.npy")
+        # if self.options.clean_house:
+            # os.system("rm D.npy  S_p.npy")
 
         return indices
diff --git a/cma.py b/cma.py
@@ -21,7 +21,7 @@
 from concordantmodes.s_vectors import SVectors
 from concordantmodes.submit import Submit
 from concordantmodes.ted import TED
-from concordantmodes.trans_disp import TransDisp
+from concordantmodes.transf_disp import TransfDisp
 from concordantmodes.vulcan_template import VulcanTemplate
 from concordantmodes.sapelo_template import SapeloTemplate
 from concordantmodes.zmat import Zmat
@@ -62,6 +62,8 @@ def run(self):
         )
         if self.options.man_proj:
             proj = self.proj
+            print(proj.shape)
+            print(proj)
             s_vec.run(
                 self.zmat_obj.cartesians_init,
                 True,
@@ -107,7 +109,7 @@ def run(self):
             else:
                 indices = np.arange(len(eigs_init))
 
-            init_disp = TransDisp(
+            init_disp = TransfDisp(
                 s_vec,
                 self.zmat_obj,
                 self.options.disp,
@@ -281,7 +283,7 @@ def run(self):
             self.options.proj_tol,
             self.zmat_obj,
             self.TED_obj,
-            "init",
+            cma="init",
         )
         init_GF.run()
 
@@ -299,7 +301,7 @@ def run(self):
             self.options.proj_tol,
             self.zmat_obj,
             self.TED_obj,
-            False,
+            cma=False,
         )
         TED_GF.run()
 
@@ -320,7 +322,7 @@ def run(self):
         )
         s_vec.run(self.zmat_obj.cartesians_final, False, proj=self.TED_obj.proj)
 
-        transdisp = TransDisp(
+        transf_disp = TransfDisp(
             s_vec,
             self.zmat_obj,
             self.options.disp,
@@ -337,9 +339,9 @@ def run(self):
         #    self.zmat_obj, self.zmat_obj.cartesians_final, B_tensor, self.options
         # )
         # second_B.run()
-        transdisp.run()
-        p_disp = transdisp.p_disp
-        m_disp = transdisp.m_disp
+        transf_disp.run()
+        p_disp = transf_disp.p_disp
+        m_disp = transf_disp.m_disp
         if self.options.disp_check:
             raise RuntimeError
 
@@ -356,7 +358,7 @@ def run(self):
             dir_obj = DirectoryTree(
                 progname,
                 self.zmat_obj,
-                transdisp,
+                transf_disp,
                 self.options.cart_insert,
                 p_disp,
                 m_disp,
@@ -420,16 +422,16 @@ def run(self):
         reap_obj = Reap(
             progname,
             self.zmat_obj,
-            transdisp.disp_cart,
+            transf_disp.disp_cart,
             self.options,
-            transdisp.n_coord,
+            transf_disp.n_coord,
             eigs,
             algo.indices,
             self.options.energy_regex,
             self.options.gradient_regex,
             self.options.molly_regex_init,
             self.options.success_regex,
-            #disp_sym = transdisp.disp_sym
+            #disp_sym = transf_disp.disp_sym
         )
         reap_obj.run()
 
@@ -438,7 +440,7 @@ def run(self):
         ref_en = reap_obj.ref_en
 
         fc = ForceConstant(
-            transdisp, p_en_array, m_en_array, ref_en, self.options, algo.indices
+            transf_disp, p_en_array, m_en_array, ref_en, self.options, algo.indices
         )
         fc.run()
         print("Computed Force Constants:")
@@ -460,7 +462,7 @@ def run(self):
         if self.options.coords != "ZMAT":
             g_mat.G = np.dot(self.TED_obj.proj.T, np.dot(g_mat.G, self.TED_obj.proj))
 
-        self.G = np.dot(np.dot(transdisp.eig_inv, g_mat.G), transdisp.eig_inv.T)
+        self.G = np.dot(np.dot(transf_disp.eig_inv, g_mat.G), transf_disp.eig_inv.T)
         self.G[np.abs(self.G) < self.options.tol] = 0
 
         if self.options.benchmark_full:
@@ -477,7 +479,7 @@ def run(self):
             self.options.proj_tol,
             self.zmat_obj,
             self.TED_obj,
-            cma,
+            cma=cma,
         )
         final_GF.run()
 
@@ -504,7 +506,7 @@ def run(self):
         # This code converts the force constants back into cartesian
         # coordinates and writes out an "output.default.hess" file, which
         # is of the same format as FCMFINAL of CFOUR.
-        self.F = np.dot(np.dot(transdisp.eig_inv.T, self.F), transdisp.eig_inv)
+        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))
         cart_conv = FcConv(
@@ -563,7 +565,7 @@ def run(self):
                     self.options.proj_tol,
                     self.zmat_obj,
                     self.TED_obj,
-                    cma,
+                    cma=cma,
                 )
                 final_GF.run()
 
@@ -590,7 +592,7 @@ def run(self):
                 # This code converts the force constants back into cartesian
                 # coordinates and writes out an "output.default.hess" file, which
                 # is of the same format as FCMFINAL of CFOUR.
-                # self.F = np.dot(np.dot(transdisp.eig_inv.T, self.F), transdisp.eig_inv)
+                # 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))
                 # cart_conv = FcConv(
@@ -623,15 +625,15 @@ def run(self):
             algo.run()
             print("printing algo indices", algo.indices)
             diagnostic_indices = algo.indices
-            with open("L_full.npy", "rb") as w:
-                L_full = np.load(w)
-            with open("L_0.npy", "rb") as x:
-                L_0 = np.load(x)
-
-            L_inv = LA.inv(L_full)
-            S = np.dot(L_inv, L_0)
-            print("printing overlap")
-            print(S)
+            # with open("L_full.npy", "rb") as w:
+                # L_full = np.load(w)
+            # with open("L_0.npy", "rb") as x:
+                # L_0 = np.load(x)
+
+            # L_inv = LA.inv(L_full)
+            # S = np.dot(L_inv, L_0)
+            # print("printing overlap")
+            # print(S)
             # with open("S.npy", "wb") as y:
             #    np.save(y, S)
             self.options.mode_coupling_check = False
@@ -656,7 +658,7 @@ def run(self):
                 self.options.proj_tol,
                 self.zmat_obj,
                 self.TED_obj,
-                cma,
+                cma=cma,
             )
             final_GF.run()
 
@@ -667,8 +669,8 @@ def run(self):
             print("//{:^40s}//".format(" Final TED"))
             print("////////////////////////////////////////////")
             self.TED_obj.run(np.dot(init_GF.L, final_GF.L), final_GF.freq)
-            if self.options.clean_house:
-                os.system("rm L_full.npy")
+            # if self.options.clean_house:
+                # os.system("rm L_full.npy")
 
             # This code prints out the frequencies in order of energy as well
             # as the ZPVE in several different units.
@@ -681,10 +683,10 @@ def run(self):
                 + "{:6.2f}".format(0.5 * np.sum(final_GF.freq) / 219474.6313708)
                 + " (hartrees) "
             )
-            if self.options.clean_house:
-                os.system("rm L_0.npy L_full.npy S.npy Full_fc_levelB.npy")
-        if self.options.clean_house:
-            os.system("rm S_p.npy")
-            path = os.getcwd() + "/L_full.npy"
-            if os.path.isfile(path):
-                os.system("rm L_full.npy")
+            # if self.options.clean_house:
+                # os.system("rm L_0.npy L_full.npy S.npy Full_fc_levelB.npy")
+        # if self.options.clean_house:
+            # os.system("rm S_p.npy")
+            # path = os.getcwd() + "/L_full.npy"
+            # if os.path.isfile(path):
+                # os.system("rm L_full.npy")
diff --git a/directory_tree.py b/directory_tree.py
@@ -25,7 +25,7 @@ def __init__(
         self.zmat = zmat
         self.insertion_index = insertion_index
         self.dir_name = dir_name
-        self.disps = disps  # This should be the 'TransDisp' object
+        self.disps = disps  # This should be the 'TransfDisp' object
         # nate
         self.p_disp = p_disp
         self.m_disp = m_disp
@@ -41,7 +41,7 @@ def make_input(self, data, dispp, n_at, at, index):
         if index == -1:
             print(
                 "The user needs to specify a different value for the \
-                   cartInsert keyword."
+                   cart_insert keyword."
             )
             raise RuntimeError
         else:
@@ -72,6 +72,7 @@ def run(self):
         if prog_name == "molpro" or prog_name == "psi4" or prog_name == "cfour":
             with open(self.template, "r") as file:
                 data = file.readlines()
+                # print(data)
         else:
             print("Specified program not supported: " + prog_name)
             raise RuntimeError

diff --git a/gf_method.py b/gf_method.py
@@ -11,7 +11,7 @@ class GFMethod(object):
     TODO: Insert standard uncertainties of amu_elmass and HARTREE_WAVENUM
     """
 
-    def __init__(self, G, F, tol, proj_tol, zmat, ted, cma):
+    def __init__(self, G, F, tol, proj_tol, zmat, ted, cma=None):
         self.G = G
         self.F = F
         self.tol = tol
@@ -34,21 +34,21 @@ def run(self):
         self.L = np.real(self.L)
         L = np.absolute(self.L)
         L_p = np.real(self.L_p)
-        S_p = np.dot(np.absolute(LA.inv(L_p)), np.absolute(L_p))
+        # S_p = np.dot(np.absolute(LA.inv(L_p)), np.absolute(L_p))
         self.L_p = L_p
         # make sure the correct "primed" or unprimed version is being saved in each instance
-        if self.cma:
-            with open("L_full.npy", "wb") as z:
-                np.save(z, L)
-        # if self.cma == False:
-        #    with open("L_intermediate.npy", "wb") as z:
-        #        np.save(z, L)
-        if self.cma is None:
-            with open("L_0.npy", "wb") as z:
-                np.save(z, L)
-        if self.cma == "init":
-            with open("S_p.npy", "wb") as z:
-                np.save(z, S_p)
+        # if self.cma:
+            # with open("L_full.npy", "wb") as z:
+                # np.save(z, L)
+        # # if self.cma == False:
+        # #    with open("L_intermediate.npy", "wb") as z:
+        # #        np.save(z, L)
+        # if self.cma is None:
+            # with open("L_0.npy", "wb") as z:
+                # np.save(z, L)
+        # if self.cma == "init":
+            # with open("S_p.npy", "wb") as z:
+                # np.save(z, S_p)
         # Construct the normal mode overlap matrix. Will be useful for off-diagonal diagnostics.
         L = np.absolute(np.real(self.L))
         L_inv = LA.inv(L)
@@ -57,13 +57,18 @@ def run(self):
 
         self.L[np.abs(self.L) < self.tol] = 0
         # Compute the frequencies by the square root of the eigenvalues.
-        self.freq = np.sqrt(self.eig_v, dtype=np.complex)
+        self.freq = np.sqrt(self.eig_v, dtype=complex)
         # Filter for imaginary modes.
         for i in range(len(self.freq)):
             if np.real(self.freq[i]) > 0.0:
                 self.freq[i] = np.real(self.freq[i])
             else:
                 self.freq[i] = -np.imag(self.freq[i])
+
+        # # This seems to fix the constant warning I get about casting out imaginary values from the complex numbers, 
+        # # however I will need to test it on transition states to see if it captures the imaginary frequencies.
+        # self.freq = np.real(self.freq)
+        # print(self.freq)
         self.freq = self.freq.astype(float)
         # Convert from Hartrees to wavenumbers.
         self.freq *= self.HARTREE_WAVENUM

diff --git a/int2cart.py b/int2cart.py
@@ -2,7 +2,6 @@
 from numpy import linalg as LA
 from . import masses
 
-
 class Int2Cart(object):
     def __init__(self, zmat, var_dict):
         self.zmat = zmat

diff --git a/options.py b/options.py
@@ -1,7 +1,7 @@
 class Options(object):
     def __init__(self, **kwargs):
         self.benchmark_full = kwargs.pop("benchmark_full", False)
-        self.bond_tol = kwargs.pop("inter_tol", 3.778)
+        self.bond_threshold = kwargs.pop("bond_threshold", 1.2)
         self.calc = kwargs.pop("calc", True)
         self.cart_coords = kwargs.pop("cart_coords", "Bohr")
         self.calc_init = kwargs.pop("calc_init", True)
@@ -10,6 +10,7 @@ def __init__(self, **kwargs):
         self.clean_house = kwargs.pop("clean_house", True)
         self.cluster = kwargs.pop("cluster", "vulcan")
         self.coords = kwargs.pop("coords", "ZMAT")
+        self.covalent_radii = kwargs.pop("covalent_radii", False)
         self.deriv_level = kwargs.pop("deriv_level", 0)
         self.deriv_level_init = kwargs.pop("deriv_level_init", 0)
         self.dir_reap = kwargs.pop("dir_reap", True)
@@ -22,7 +23,6 @@ def __init__(self, **kwargs):
         self.gen_disps_init = kwargs.pop("gen_disps_init",True)
         self.geom_check = kwargs.pop("geom_check", False)
         self.gradient_regex = kwargs.pop("gradient_regex", "")
-        self.interatomic_distance = kwargs.pop("interatomic_distance", False)
         self.man_proj = kwargs.pop("man_proj", False)
         self.mode_coupling_check = kwargs.pop("mode_coupling_check", False)
         self.molly_regex_init = kwargs.pop("molly_regex_init", "")
@@ -43,4 +43,6 @@ def __init__(self, **kwargs):
         self.success_regex_init = kwargs.pop("success_regex_init", "")
         self.tight_disp = kwargs.pop("tight_disp", False)
         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)
         self.units = kwargs.pop("units", "HartreeBohr")