From 413a176ac09de055f24b117ba94d914aa976877c Mon Sep 17 00:00:00 2001 From: Mitchell Evan Lahm Date: Thu, 25 Jun 2026 20:22:48 -0400 Subject: [PATCH 1/2] First draft refactor of topological code --- concordantmodes/zmat.py | 481 +++++++++--------- .../molpro/CMA0A/No_Spec_Deloc/main.py | 4 +- .../molpro/CMA0A/nospec_deloc/main.py | 29 ++ .../molpro/CMA0A/nospec_deloc/templateA.dat | 12 + .../molpro/CMA0A/nospec_deloc/templateB.dat | 11 + .../5_benzene/molpro/CMA0A/nospec_deloc/zmat | 59 +++ 6 files changed, 365 insertions(+), 231 deletions(-) create mode 100644 examples/5_benzene/molpro/CMA0A/nospec_deloc/main.py create mode 100644 examples/5_benzene/molpro/CMA0A/nospec_deloc/templateA.dat create mode 100644 examples/5_benzene/molpro/CMA0A/nospec_deloc/templateB.dat create mode 100644 examples/5_benzene/molpro/CMA0A/nospec_deloc/zmat diff --git a/concordantmodes/zmat.py b/concordantmodes/zmat.py index 58aa0e98..ac5d6bd6 100644 --- a/concordantmodes/zmat.py +++ b/concordantmodes/zmat.py @@ -386,241 +386,20 @@ def zmat_process(self, zmat_output): ) self.torsion_variables.append("D" + str(i - first_index)) elif self.options.coords.upper() == "DELOCALIZED": - count = 0 - if self.options.covalent_radii: - c_r = CovalentRadii() - indices = [] - transdisp_inter = TransfDisp( - None, - self, - 1, - False, - np.array([]), - self.options, - indices, - ) - inter_atomic_len = np.zeros( - (len(self.cartesians_b), len(self.cartesians_b)) - ) - N = len(self.cartesians_b) - adj_mat = np.zeros((N, N)) - for i in range(len(self.cartesians_b)): - for j in range(i): - inter_atomic_len[j, i] = transdisp_inter.calc_bond( - self.cartesians_b[i], self.cartesians_b[j] - ) - if inter_atomic_len[j, i] < self.options.bond_threshold * ( - c_r.get(self.atom_list[i]) + c_r.get(self.atom_list[j]) - ): - count += 1 - adj_mat[i, j] = 1 - adj_mat[j, i] = 1 - self.bond_indices = np.append(self.bond_indices, 0) - self.bond_indices[-1] = np.array( - [str(j + 1), str(i + 1)], dtype=object - ) - self.bond_variables.append("R" + str(count)) - print("Interatomic Distance Matrix:") - print(inter_atomic_len) - print("Adjacency Matrix:") - print(adj_mat) - for i in range(len(adj_mat)): - print("Degree of vertex " + str(i)) - print(np.sum(adj_mat[i])) - adj_mat2 = np.dot(adj_mat, adj_mat) - for i in range(len(adj_mat2)): - adj_mat2[i, i] = 0 - print("Covalent Radius Scale Factor:") - print(self.options.bond_threshold) - print("Resulting bond indices:") - print(self.bond_indices) - else: - for i in range(len(zmat_output)): - if re.search(self.bond_regex, zmat_output[i]): - count += 1 - List = re.findall(self.bond_regex, zmat_output[i])[0] - self.bond_indices = np.append(self.bond_indices, 0) - self.bond_indices[-1] = np.array(List, dtype=object) - self.bond_variables.append("R" + str(count)) + + # Form all possible bonds, whether from qcelemental or manually input + self._build_bonds(zmat_output) # Form all possible angles from bonds - count = 0 - for i in range(len(self.bond_indices)): - for j in range(len(self.bond_indices) - i - 1): - a = np.setdiff1d(self.bond_indices[i], self.bond_indices[i + j + 1]) - b = np.intersect1d( - self.bond_indices[i], self.bond_indices[i + j + 1] - ) - c = np.setdiff1d(self.bond_indices[i + j + 1], self.bond_indices[i]) - if len(a) and len(b) and len(c): - d = np.array([a[0], b[0], c[0]]) - self.angle_indices = np.append(self.angle_indices, 0) - self.angle_indices[-1] = np.array(d, dtype=object) - count += 1 - self.angle_variables.append("A" + str(count)) + self._build_angles() # Form all possible torsions from angles - tor_count = 0 - oop_count = 0 - for i in range(len(self.angle_indices)): - for j in range(len(self.bond_indices)): - a = np.setdiff1d(self.angle_indices[i], self.bond_indices[j]) - b = np.intersect1d(self.angle_indices[i], self.bond_indices[j]) - c = np.setdiff1d(self.bond_indices[j], self.angle_indices[i]) - if len(c) == 1: - d = np.where(self.bond_indices[j] == c)[0][0] - f = 1 - d - g = self.bond_indices[j][f] - h = np.where(self.angle_indices[i] == g)[0][0] - if h == 1: - # This is an out of plane bend - oop = np.array( - [ - self.bond_indices[j][d], - self.bond_indices[j][f], - a[0], - a[1], - ], - dtype=object, - ) - - cont_bool = self.np_contains(self.oop_indices, oop) - if not cont_bool: - self.oop_indices = np.append(self.oop_indices, 0) - self.oop_indices[-1] = np.array(oop, dtype=object) - else: - # This is a torsion - if h: - tor = np.append( - self.angle_indices[i].copy(), - self.bond_indices[j][d], - ) - else: - tor = np.append( - self.bond_indices[j][d], - self.angle_indices[i].copy(), - ) - cont_bool = self.np_contains( - self.torsion_indices, tor, tor=True - ) - if not cont_bool: - self.torsion_indices = np.append( - self.torsion_indices, 0 - ) - self.torsion_indices[-1] = np.array(tor, dtype=object) - - for i in range(len(self.torsion_indices)): - self.torsion_variables.append("D" + str(i + 1)) - for i in range(len(self.oop_indices)): - self.oop_variables.append("O" + str(i + 1)) + self._build_torsions() + # Perform a topological analysis, this is the first step to the + # automatic generation of Natural Internal Coordinates if self.options.topo_analysis: - X_len_walks_dict = {} - X_len_walks_dict["2_length_walks"] = self.bond_indices.copy() - X_len_walks_dict["3_length_walks"] = self.angle_indices.copy() - X_len_walks_dict["4_length_walks"] = self.torsion_indices.copy() - - count = 4 - - prev_walks = self.torsion_indices.copy() - - while True: - new_walks = np.array([]) - for i in range(len(prev_walks)): - for j in range(len(self.bond_indices)): - a = np.where(prev_walks[i] == self.bond_indices[j][0])[0] - b = np.where(prev_walks[i] == self.bond_indices[j][1])[0] - if len(a) and not len(b): - if not a[0]: - new_walk = np.append( - self.bond_indices[j][1], prev_walks[i].copy() - ) - new_walks = np.append(new_walks, new_walk) - elif a[0] == count - 1: - new_walk = np.append( - prev_walks[i].copy(), self.bond_indices[j][1] - ) - new_walks = np.append(new_walks, new_walk) - if len(b) and not len(a): - if not b[0]: - new_walk = np.append( - self.bond_indices[j][0], prev_walks[i].copy() - ) - new_walks = np.append(new_walks, new_walk) - elif b[0] == count - 1: - new_walk = np.append( - prev_walks[i].copy(), self.bond_indices[j][0] - ) - new_walks = np.append(new_walks, new_walk) - - count += 1 - new_walks = new_walks.reshape((-1, count)) - print(new_walks) - new_walks = np.unique(new_walks, axis=0) - - del_list = np.array([]) - for i in range(len(new_walks)): - for j in range(len(new_walks) - i - 1): - a = np.array([new_walks[i], np.flip(new_walks[i + j + 1])]) - a = np.unique(a, axis=0) - if len(a) == 1: - del_list = np.append(del_list, [i + j + 1]) - - del_list = del_list.astype(int) - new_walks = np.delete(new_walks, del_list, axis=0) - prev_walks = new_walks - if count > self.options.topo_max_it or not len(new_walks): - print( - "Walk generator has terminated at walk lengths of " - + str(count) - ) - break - X_len_walks_dict[str(count) + "_length_walks"] = new_walks - print(str(count) + "_length_walks") - print(len(new_walks)) - print(new_walks) - - dict_len = len(X_len_walks_dict) - 1 - - cycles_dict = {} - for i in range(dict_len): - print(str(i + 3)) - cycles_dict[str(i + 3)] = np.array([]) - for j in range(len(X_len_walks_dict[str(i + 3) + "_length_walks"])): - a = np.array( - [ - X_len_walks_dict[str(i + 3) + "_length_walks"][j][0], - X_len_walks_dict[str(i + 3) + "_length_walks"][j][-1], - ] - ) - for k in range(len(self.bond_indices)): - b = np.intersect1d(a, self.bond_indices[k]) - if len(b) == 2: - cycle = X_len_walks_dict[str(i + 3) + "_length_walks"][ - j - ].copy() - cycles_dict[str(i + 3)] = np.append( - cycles_dict[str(i + 3)], cycle - ) - cycles_dict[str(i + 3)] = cycles_dict[str(i + 3)].reshape( - (-1, i + 3) - ) - del_array = np.array([]) - if len(cycles_dict[str(i + 3)]): - for j in range(len(cycles_dict[str(i + 3)])): - for k in range(len(cycles_dict[str(i + 3)]) - j - 1): - a = np.intersect1d( - cycles_dict[str(i + 3)][j], - cycles_dict[str(i + 3)][k + j + 1], - ) - if len(a) == len(cycles_dict[str(i + 3)][0]): - del_array = np.append(del_array, [k + j + 1]) - del_array = del_array.astype(int) - del_array = np.unique(del_array) - cycles_dict[str(i + 3)] = np.delete( - cycles_dict[str(i + 3)], del_array, axis=0 - ) - print(cycles_dict[str(i + 3)]) + self._generate_topology() elif self.options.coords.upper() == "CUSTOM": # This option will allow the user to specify a custom array of @@ -1025,3 +804,247 @@ def np_contains(self, array1, array2, tor=False): cont_bool = True return cont_bool + + def _build_bonds(self, zmat_output): + bond_count = 0 + if self.options.covalent_radii: + c_r = CovalentRadii() + indices = [] + transdisp_inter = TransfDisp( + None, + self, + 1, + False, + np.array([]), + self.options, + indices, + ) + inter_atomic_len = np.zeros( + (len(self.cartesians_b), len(self.cartesians_b)) + ) + N = len(self.cartesians_b) + adj_mat = np.zeros((N, N)) + for i in range(len(self.cartesians_b)): + for j in range(i): + inter_atomic_len[j, i] = transdisp_inter.calc_bond( + self.cartesians_b[i], self.cartesians_b[j] + ) + if inter_atomic_len[j, i] < self.options.bond_threshold * ( + c_r.get(self.atom_list[i]) + c_r.get(self.atom_list[j]) + ): + bond_count += 1 + adj_mat[i, j] = 1 + adj_mat[j, i] = 1 + self.bond_indices = np.append(self.bond_indices, 0) + self.bond_indices[-1] = np.array( + [str(j + 1), str(i + 1)], dtype=object + ) + self.bond_variables.append("R" + str(bond_count)) + print("Interatomic Distance Matrix:") + print(inter_atomic_len) + print("Adjacency Matrix:") + print(adj_mat) + for i in range(len(adj_mat)): + print("Degree of vertex " + str(i)) + print(np.sum(adj_mat[i])) + adj_mat2 = np.dot(adj_mat, adj_mat) + for i in range(len(adj_mat2)): + adj_mat2[i, i] = 0 + print("Covalent Radius Scale Factor:") + print(self.options.bond_threshold) + # print("Resulting bond indices:") + # print(self.bond_indices) + else: + for i in range(len(zmat_output)): + if re.search(self.bond_regex, zmat_output[i]): + bond_count += 1 + List = re.findall(self.bond_regex, zmat_output[i])[0] + self.bond_indices = np.append(self.bond_indices, 0) + self.bond_indices[-1] = np.array(List, dtype=object) + self.bond_variables.append("R" + str(bond_count)) + + def _build_angles(self): + ang_count = 0 + for i in range(len(self.bond_indices)): + for j in range(len(self.bond_indices) - i - 1): + a = np.setdiff1d(self.bond_indices[i], self.bond_indices[i + j + 1]) + b = np.intersect1d( + self.bond_indices[i], self.bond_indices[i + j + 1] + ) + c = np.setdiff1d(self.bond_indices[i + j + 1], self.bond_indices[i]) + if len(a) and len(b) and len(c): + d = np.array([a[0], b[0], c[0]]) + self.angle_indices = np.append(self.angle_indices, 0) + self.angle_indices[-1] = np.array(d, dtype=object) + ang_count += 1 + self.angle_variables.append("A" + str(ang_count)) + + def _build_torsions(self): + for i in range(len(self.angle_indices)): + for j in range(len(self.bond_indices)): + a = np.setdiff1d(self.angle_indices[i], self.bond_indices[j]) + b = np.intersect1d(self.angle_indices[i], self.bond_indices[j]) + c = np.setdiff1d(self.bond_indices[j], self.angle_indices[i]) + if len(c) == 1: + d = np.where(self.bond_indices[j] == c)[0][0] + f = 1 - d + g = self.bond_indices[j][f] + h = np.where(self.angle_indices[i] == g)[0][0] + if h == 1: + # This is an out of plane bend + oop = np.array( + [ + self.bond_indices[j][d], + self.bond_indices[j][f], + a[0], + a[1], + ], + dtype=object, + ) + + cont_bool = self.np_contains(self.oop_indices, oop) + if not cont_bool: + self.oop_indices = np.append(self.oop_indices, 0) + self.oop_indices[-1] = np.array(oop, dtype=object) + else: + # This is a torsion + if h: + tor = np.append( + self.angle_indices[i].copy(), + self.bond_indices[j][d], + ) + else: + tor = np.append( + self.bond_indices[j][d], + self.angle_indices[i].copy(), + ) + cont_bool = self.np_contains( + self.torsion_indices, tor, tor=True + ) + if not cont_bool: + self.torsion_indices = np.append( + self.torsion_indices, 0 + ) + self.torsion_indices[-1] = np.array(tor, dtype=object) + + for i in range(len(self.torsion_indices)): + self.torsion_variables.append("D" + str(i + 1)) + for i in range(len(self.oop_indices)): + self.oop_variables.append("O" + str(i + 1)) + + def _generate_topology(self): + X_len_walks_dict = {} + X_len_walks_dict["2_length_walks"] = self.bond_indices.copy() + X_len_walks_dict["3_length_walks"] = self.angle_indices.copy() + X_len_walks_dict["4_length_walks"] = self.torsion_indices.copy() + + count = 4 + + prev_walks = self.torsion_indices.copy() + + while True: + new_walks = np.array([]) + for i in range(len(prev_walks)): + for j in range(len(self.bond_indices)): + a = np.where(prev_walks[i] == self.bond_indices[j][0])[0] + b = np.where(prev_walks[i] == self.bond_indices[j][1])[0] + if len(a) and not len(b): + if not a[0]: + new_walk = np.append( + self.bond_indices[j][1], prev_walks[i].copy() + ) + new_walks = np.append(new_walks, new_walk) + elif a[0] == count - 1: + new_walk = np.append( + prev_walks[i].copy(), self.bond_indices[j][1] + ) + new_walks = np.append(new_walks, new_walk) + if len(b) and not len(a): + if not b[0]: + new_walk = np.append( + self.bond_indices[j][0], prev_walks[i].copy() + ) + new_walks = np.append(new_walks, new_walk) + elif b[0] == count - 1: + new_walk = np.append( + prev_walks[i].copy(), self.bond_indices[j][0] + ) + new_walks = np.append(new_walks, new_walk) + + count += 1 + # new_walks = new_walks.reshape((-1, count)) + new_walks = new_walks.reshape((-1, count)).tolist() + print(new_walks) + new_walks = np.unique(new_walks, axis=0) + + del_list = np.array([]) + for i in range(len(new_walks)): + for j in range(len(new_walks) - i - 1): + a = np.array([new_walks[i], np.flip(new_walks[i + j + 1])]) + a = np.unique(a, axis=0) + if len(a) == 1: + del_list = np.append(del_list, [i + j + 1]) + + del_list = del_list.astype(int) + new_walks = np.delete(new_walks, del_list, axis=0) + prev_walks = new_walks + if count > self.options.topo_max_it or not len(new_walks): + print( + "Walk generator has terminated at walk lengths of " + + str(count) + ) + break + X_len_walks_dict[str(count) + "_length_walks"] = new_walks + print(str(count) + "_length_walks") + print(len(new_walks)) + print(new_walks) + + dict_len = len(X_len_walks_dict) - 1 + + print('Looking for rings') + cycles_dict = {} + for i in range(dict_len): + # print(str(i + 3)) + cycles_dict[str(i + 3)] = np.array([]) + for j in range(len(X_len_walks_dict[str(i + 3) + "_length_walks"])): + a = np.array( + [ + X_len_walks_dict[str(i + 3) + "_length_walks"][j][0], + X_len_walks_dict[str(i + 3) + "_length_walks"][j][-1], + ] + ) + for k in range(len(self.bond_indices)): + b = np.intersect1d(a, self.bond_indices[k]) + if len(b) == 2: + cycle = X_len_walks_dict[str(i + 3) + "_length_walks"][ + j + ].copy() + cycles_dict[str(i + 3)] = np.append( + cycles_dict[str(i + 3)], cycle + ) + cycles_dict[str(i + 3)] = cycles_dict[str(i + 3)].reshape( + (-1, i + 3) + ) + del_array = np.array([]) + if len(cycles_dict[str(i + 3)]): + for j in range(len(cycles_dict[str(i + 3)])): + for k in range(len(cycles_dict[str(i + 3)]) - j - 1): + a = np.intersect1d( + cycles_dict[str(i + 3)][j], + cycles_dict[str(i + 3)][k + j + 1], + ) + if len(a) == len(cycles_dict[str(i + 3)][0]): + del_array = np.append(del_array, [k + j + 1]) + del_array = del_array.astype(int) + del_array = np.unique(del_array) + cycles_dict[str(i + 3)] = np.delete( + cycles_dict[str(i + 3)], del_array, axis=0 + ) + print(f"{i+3}-membered ring detected!" ) + print(cycles_dict[str(i + 3)]) + + def _generate_walks(self): + pass + + def _find_cycles(self): + pass diff --git a/examples/2_Methanol/molpro/CMA0A/No_Spec_Deloc/main.py b/examples/2_Methanol/molpro/CMA0A/No_Spec_Deloc/main.py index ac338843..690c5fab 100644 --- a/examples/2_Methanol/molpro/CMA0A/No_Spec_Deloc/main.py +++ b/examples/2_Methanol/molpro/CMA0A/No_Spec_Deloc/main.py @@ -11,9 +11,9 @@ "covalent_radii": True, "coords": "Delocalized", "topo_analysis": True, - # "calc_b" : False, + "calc_b" : False, # "calc_a" : False, - # "gen_disps_b" : False, + "gen_disps_b" : False, # "gen_disps_a" : False, "success_regex_b": r"Molpro calculation terminated", "success_regex_a": r"Molpro calculation terminated", diff --git a/examples/5_benzene/molpro/CMA0A/nospec_deloc/main.py b/examples/5_benzene/molpro/CMA0A/nospec_deloc/main.py new file mode 100644 index 00000000..c0da392f --- /dev/null +++ b/examples/5_benzene/molpro/CMA0A/nospec_deloc/main.py @@ -0,0 +1,29 @@ +from concordantmodes.options import Options + +options_kwargs = { + "cluster": "slurm", + "program_b": "molpro", + "program_a": "molpro", + "energy_regex_a": r"\(T\) total energy\s+(\-\d+\.\d+)", + "energy_regex_b": r"!RHF STATE 1.1 Energy\s*\s+(\-\d+\.\d+)", + "cart_insert_b": 7, + "cart_insert_a": 7, + "covalent_radii": True, + "coords": "Delocalized", + "topo_analysis": True, + "calc_b" : False, + # "calc_a" : False, + "gen_disps_b" : False, + # "gen_disps_a" : False, + "success_regex_a": r"Molpro calculation terminated", + "success_regex_b": r"Molpro calculation terminated", + # "symmetry": True, + # "autosalcs": True, +} +options_obj = Options(**options_kwargs) + +# 3. call Concordant Modes Program +from concordantmodes.cma import ConcordantModes + +CMA_obj = ConcordantModes(options_obj) +CMA_obj.run() diff --git a/examples/5_benzene/molpro/CMA0A/nospec_deloc/templateA.dat b/examples/5_benzene/molpro/CMA0A/nospec_deloc/templateA.dat new file mode 100644 index 00000000..84f082f2 --- /dev/null +++ b/examples/5_benzene/molpro/CMA0A/nospec_deloc/templateA.dat @@ -0,0 +1,12 @@ +*** Optimization of min3 +memory,1000,m +gthresh,energy=1.0d-12,orbital=1.0d-10,oneint=1.0d-16,twoint=1.0d-16,optgrad=1.0d-6,compress=1.0d-13 +geomtyp=xyz + + bohr + geometry = { +} + +basis=cc-pvtz +{rhf;maxit,200} +{ccsd(t);maxit,200} diff --git a/examples/5_benzene/molpro/CMA0A/nospec_deloc/templateB.dat b/examples/5_benzene/molpro/CMA0A/nospec_deloc/templateB.dat new file mode 100644 index 00000000..60d0b1c8 --- /dev/null +++ b/examples/5_benzene/molpro/CMA0A/nospec_deloc/templateB.dat @@ -0,0 +1,11 @@ +*** Optimization of min3 +memory,1000,m +gthresh,energy=1.0d-12,orbital=1.0d-10,oneint=1.0d-16,twoint=1.0d-16,optgrad=1.0d-6,compress=1.0d-13 +geomtyp=xyz + + bohr + geometry = { +} + +basis=cc-pvdz +{rhf;maxit,200} diff --git a/examples/5_benzene/molpro/CMA0A/nospec_deloc/zmat b/examples/5_benzene/molpro/CMA0A/nospec_deloc/zmat new file mode 100644 index 00000000..00ed2f0b --- /dev/null +++ b/examples/5_benzene/molpro/CMA0A/nospec_deloc/zmat @@ -0,0 +1,59 @@ +ZMAT begin +1 2 +2 3 +3 4 +4 5 +5 6 +6 1 + +1 7 +2 8 +3 9 +4 10 +5 11 +6 12 + +6 1 2 +1 2 3 +2 3 4 +3 4 5 +4 5 6 +5 6 1 + +7 1 2 +8 2 3 +9 3 4 + +10 4 5 +11 5 6 +12 6 1 + +1 2 3 4 T +2 3 4 5 T +3 4 5 6 T +4 5 6 1 T +5 6 1 2 T +6 1 2 3 T + +7 1 2 6 O +8 2 3 1 O +9 3 4 2 O +10 4 5 3 O +11 5 6 4 O +12 6 1 5 O +ZMAT end + +cart begin + C 0.00000000 2.64096559 -0.00000000 + C -2.28714329 1.32048279 -0.00000000 + C -2.28714329 -1.32048279 0.00000000 + C -0.00000000 -2.64096559 0.00000000 + C 2.28714329 -1.32048279 0.00000000 + C 2.28714329 1.32048279 -0.00000000 + H 0.00000000 4.68777554 -0.00000000 + H -4.05973271 2.34388777 -0.00000000 + H -4.05973271 -2.34388777 0.00000000 + H -0.00000000 -4.68777554 0.00000000 + H 4.05973271 -2.34388777 0.00000000 + H 4.05973271 2.34388777 -0.00000000 +cart end From a6757d84badee78d70ec17071629fc2c24568bd3 Mon Sep 17 00:00:00 2001 From: Mitchell Evan Lahm Date: Thu, 25 Jun 2026 20:24:58 -0400 Subject: [PATCH 2/2] N/A --- concordantmodes/zmat.py | 31 ++++++++++--------------------- 1 file changed, 10 insertions(+), 21 deletions(-) diff --git a/concordantmodes/zmat.py b/concordantmodes/zmat.py index ac5d6bd6..bc79c414 100644 --- a/concordantmodes/zmat.py +++ b/concordantmodes/zmat.py @@ -386,7 +386,7 @@ def zmat_process(self, zmat_output): ) self.torsion_variables.append("D" + str(i - first_index)) elif self.options.coords.upper() == "DELOCALIZED": - + # Form all possible bonds, whether from qcelemental or manually input self._build_bonds(zmat_output) @@ -868,9 +868,7 @@ def _build_angles(self): for i in range(len(self.bond_indices)): for j in range(len(self.bond_indices) - i - 1): a = np.setdiff1d(self.bond_indices[i], self.bond_indices[i + j + 1]) - b = np.intersect1d( - self.bond_indices[i], self.bond_indices[i + j + 1] - ) + b = np.intersect1d(self.bond_indices[i], self.bond_indices[i + j + 1]) c = np.setdiff1d(self.bond_indices[i + j + 1], self.bond_indices[i]) if len(a) and len(b) and len(c): d = np.array([a[0], b[0], c[0]]) @@ -922,9 +920,7 @@ def _build_torsions(self): self.torsion_indices, tor, tor=True ) if not cont_bool: - self.torsion_indices = np.append( - self.torsion_indices, 0 - ) + self.torsion_indices = np.append(self.torsion_indices, 0) self.torsion_indices[-1] = np.array(tor, dtype=object) for i in range(len(self.torsion_indices)): @@ -989,10 +985,7 @@ def _generate_topology(self): new_walks = np.delete(new_walks, del_list, axis=0) prev_walks = new_walks if count > self.options.topo_max_it or not len(new_walks): - print( - "Walk generator has terminated at walk lengths of " - + str(count) - ) + print("Walk generator has terminated at walk lengths of " + str(count)) break X_len_walks_dict[str(count) + "_length_walks"] = new_walks print(str(count) + "_length_walks") @@ -1000,8 +993,8 @@ def _generate_topology(self): print(new_walks) dict_len = len(X_len_walks_dict) - 1 - - print('Looking for rings') + + print("Looking for rings") cycles_dict = {} for i in range(dict_len): # print(str(i + 3)) @@ -1016,15 +1009,11 @@ def _generate_topology(self): for k in range(len(self.bond_indices)): b = np.intersect1d(a, self.bond_indices[k]) if len(b) == 2: - cycle = X_len_walks_dict[str(i + 3) + "_length_walks"][ - j - ].copy() + cycle = X_len_walks_dict[str(i + 3) + "_length_walks"][j].copy() cycles_dict[str(i + 3)] = np.append( cycles_dict[str(i + 3)], cycle ) - cycles_dict[str(i + 3)] = cycles_dict[str(i + 3)].reshape( - (-1, i + 3) - ) + cycles_dict[str(i + 3)] = cycles_dict[str(i + 3)].reshape((-1, i + 3)) del_array = np.array([]) if len(cycles_dict[str(i + 3)]): for j in range(len(cycles_dict[str(i + 3)])): @@ -1040,9 +1029,9 @@ def _generate_topology(self): cycles_dict[str(i + 3)] = np.delete( cycles_dict[str(i + 3)], del_array, axis=0 ) - print(f"{i+3}-membered ring detected!" ) + print(f"{i+3}-membered ring detected!") print(cycles_dict[str(i + 3)]) - + def _generate_walks(self): pass