temp2

arc/plotter.py

@@ -1015,12 +1015,15 @@ def plot_torsion_angles(torsion_angles,
     """
     num_comb = None
     torsions = list(torsion_angles.keys()) if torsions_sampling_points is None \
-        else list(torsions_sampling_points.keys())
+        else list(torsions_sampling_points.keys())    
+    torsions = [torsion for torsion in torsions if not (isinstance(torsion, tuple) and all(isinstance(sub_t, tuple) for sub_t in torsion))]
+
     ticks = [0, 60, 120, 180, 240, 300, 360]
     sampling_points = dict()
     if torsions_sampling_points is not None:
         for tor, points in torsions_sampling_points.items():
-            sampling_points[tor] = [point if point <= 360 else point - 360 for point in points]
+            if not isinstance(points[0],list):
+                sampling_points[tor] = [point if point <= 360 else point - 360 for point in points]
     if not torsions:
         return
     if len(torsions) == 1:
@@ -1049,6 +1052,8 @@ def plot_torsion_angles(torsion_angles,
         fig.dpi = 120
         num_comb = 1
         for i, torsion in enumerate(torsions):
+            if tuple(torsion) not in torsion_angles:
+                continue
             axs[i].plot(np.array(torsion_angles[tuple(torsion)]),
                         np.zeros_like(np.arange(len(torsion_angles[tuple(torsion)]))), 'g.')
             if wells_dict is not None:
@@ -1122,88 +1127,75 @@ def plot_torsion_angles(torsion_angles,
     return num_comb
 
 
-def plot_ring_torsion_angles(conformers, plot_path=None):
+def plot_ring_torsion_angles(conformers, plot_path=None, tolerance=15):
     """
-    Plot the torsion angles of the generated conformers.
+    Plot the torsion angles of the generated conformers for each ring,
+    considering torsion similarity within a given tolerance.
 
     Args:
-        torsion_angles (dict): Keys are torsions, values are lists of corresponding angles.
-        torsions_sampling_points (dict, optional): Keys are torsions, values are sampling points.
-        wells_dict (dict, optional): Keys are torsions, values are lists of wells.
-                                     Each entry in such a list is a well dictionary with the following keys:
-                                     ``start_idx``, ``end_idx``, ``start_angle``, ``end_angle``, and ``angles``.
-        e_conformers (list, optional): Entries are conformers corresponding to the sampling points with FF energies.
-        de_threshold (float, optional): Energy threshold, plotted as a dashed horizontal line.
+        conformers (list): A list of conformers, each containing a 'puckering' key with ring torsion angles.
         plot_path (str, optional): The path for saving the plot.
+        tolerance (float, optional): The angular tolerance to consider two torsion angle sets as similar.
     """
-    num_comb = None
     if 'puckering' not in conformers[0]:
         return
-    ring_num = len(conformers[0]['puckering'])
-    # fig, axs = plt.subplots(nrows=ring_num, ncols=1, sharex=True, sharey=True, gridspec_kw={'hspace': 0})
-    # fig.dpi = 120
-    num_comb = 1
-    for i in range(ring_num):
-        unique_angle_sets = {}
-        for conformer in conformers:
-            angles = conformer['puckering']
-            for ring, angle_set in angles.items():
-                # Round angles to the nearest integer
-                rounded_angle_set = tuple(round(angle) for angle in angle_set)
-                found = False
-                for existing_angle_set in unique_angle_sets.keys():
-                    if almost_equal_lists(existing_angle_set, rounded_angle_set):
-                        unique_angle_sets[existing_angle_set] += 1
-                        found = True
-                        break
-                if not found:
-                    unique_angle_sets[rounded_angle_set] = 1
-    angles = [list(angles) for angles in unique_angle_sets.keys()]
-    counts = list(unique_angle_sets.values())
-
-    # Combine angles and counts into a single list of tuples for sorting
-    angles_counts = list(zip(angles, counts))
-
-    # Sort by counts in descending order
-    angles_counts_sorted = sorted(angles_counts, key=lambda x: x[1], reverse=True)
-
-    # Unzip the sorted tuples back into separate lists
-    angles_sorted, counts_sorted = zip(*angles_counts_sorted)
+
+    # Dictionary to store unique angle sets for each ring
+    ring_angle_data = {}
 
-    fig, ax = plt.subplots(figsize=(10, 5))
-    x = np.arange(len(angles_sorted))  # the label locations
+    # Process each conformer
+    for conformer in conformers:
+        rings = conformer['puckering']  # Retrieve the puckering angles for rings
+        for torsions, angle_set in rings.items():
+            rounded_angle_set = tuple(round(angle) for angle in angle_set)  # Round angles
+            if torsions not in ring_angle_data:
+                ring_angle_data[torsions] = []
+
+            # Check for similarity within the current ring
+            is_similar = False
+            for i, (existing_set, count) in enumerate(ring_angle_data[torsions]):
+                if all(abs(a1 - a2) <= tolerance for a1, a2 in zip(rounded_angle_set, existing_set)):
+                    # If similar, increment count
+                    ring_angle_data[torsions][i] = (existing_set, count + 1)
+                    is_similar = True
+                    break
+            if not is_similar:
+                # Add unique angle set with a count
+                ring_angle_data[torsions].append((rounded_angle_set, 1))
+
 
-    ax.bar(x, counts_sorted, color='blue')
-    ax.set_xlabel('Base Angles Set (Rounded to Nearest Integer)')
-    ax.set_ylabel('Frequency')
-    ax.set_title('Frequency of Different Angle Sets Across Conformers')
-    ax.set_xticks(x)
-    ax.set_xticklabels([f'{angle}' for angle in angles_sorted], rotation=45, ha='right')
+    # Plot data for each ring
+    for ring, angle_counts in ring_angle_data.items():
+        # Extract and sort data
+        angles, counts = zip(*angle_counts)
+        angles_counts_sorted = sorted(zip(angles, counts), key=lambda x: x[1], reverse=True)
+        angles_sorted, counts_sorted = zip(*angles_counts_sorted)
+
+        # Create bar plot for this ring
+        fig, ax = plt.subplots(figsize=(10, 5))
+        x = np.arange(len(angles_sorted))  # Label positions
+        ax.bar(x, counts_sorted, color='blue')
+        ax.set_xlabel('Rounded Angle Sets (Nearest Integer)')
+        ax.set_ylabel('Frequency')
+        ax.set_title(f'Frequency of Different Angle Sets for Ring {ring}')
+        ax.set_xticks(x)
+        ax.set_xticklabels([f'{angle}' for angle in angles_sorted], rotation=45, ha='right')
+
+        # Save or display the plot
+        if plot_path is not None:
+            ring_plot_path = os.path.join(plot_path, f'conformer_ring_torsions_{ring}.png')
+            if not os.path.isdir(plot_path):
+                os.makedirs(plot_path)
+            try:
+                plt.savefig(ring_plot_path, bbox_inches='tight')
+            except FileNotFoundError:
+                pass
+        if is_notebook():
+            plt.show()
+        plt.close(fig)
 
-    # plt.tight_layout()
-    # plt.show()
-    # plt.setp(axs, xticks=ticks)  # set the x ticks of all subplots
-    # fig.set_size_inches(8, len(torsions) * 1.5)
-    if plot_path is not None:
-        if not os.path.isdir(plot_path):
-            os.makedirs(plot_path)
-        file_names = list()
-        for (_, _, files) in os.walk(plot_path):
-            file_names.extend(files)
-            break  # don't continue to explore subdirectories
-        i = 0
-        for file_ in file_names:
-            if 'conformer torsions' in file_:
-                i += 1
-        image_path = os.path.join(plot_path, f'conformer ring torsions {i}.png')
-        try:
-            plt.savefig(image_path, bbox_inches='tight')
-        except FileNotFoundError:
-            pass
-    if is_notebook():
-        plt.show()
-    plt.close(fig)
-    return num_comb
+    # Return the collected data
+    return ring_angle_data
 
 
 def plot_1d_rotor_scan(angles: Optional[Union[list, tuple, np.array]] = None,

arc/species/conformers.py

@@ -352,19 +352,27 @@ def deduce_new_conformers(label, conformers, torsions, tops, mol_list, smeared_s
         symmetries[tuple(torsion)] = symmetry
     logger.debug(f'Identified {len([s for s in symmetries.values() if s > 1])} symmetric wells for {label}')
 
-    torsions_sampling_points, wells_dict = dict(), dict()
+    torsions_sampling_points, wells_dict, ring_sampling_points = dict(), dict(), dict()
     for tor, tor_angles in torsion_angles.items():
         torsions_sampling_points[tor], wells_dict[tor] = \
             determine_torsion_sampling_points(label, tor_angles, smeared_scan_res=smeared_scan_res,
                                               symmetry=symmetries[tor])
-
     if plot_path is not None:
         arc.plotter.plot_torsion_angles(torsion_angles, torsions_sampling_points, wells_dict=wells_dict,
                                         plot_path=plot_path)
-        arc.plotter.plot_ring_torsion_angles(conformers=conformers, plot_path=plot_path)
+        angles_sorted = arc.plotter.plot_ring_torsion_angles(conformers=conformers, plot_path=plot_path)
+    # Process the ring angle data
+    for ring, angle_counts in angles_sorted.items():
+        # Extract unique angles only
+        angles = [list(angle) for angle, _ in angle_counts]  # Ignore counts, only keep the angles
+
+        # Flatten the structure and add to torsions_sampling_points
+        ring_sampling_points[tuple(ring)] = angles
+
+    combined_sampling_points = {**torsions_sampling_points, **ring_sampling_points}
 
     hypothetical_num_comb = 1
-    for points in torsions_sampling_points.values():
+    for points in combined_sampling_points.values():
         hypothetical_num_comb *= len(points)
     number_of_chiral_centers = get_number_of_chiral_centers(label, mol, conformer=conformers[0],
                                                             just_get_the_number=True)
@@ -375,10 +383,10 @@ def deduce_new_conformers(label, conformers, torsions, tops, mol_list, smeared_s
         hypothetical_num_comb_str = str(hypothetical_num_comb)
     logger.info(f'\nHypothetical number of conformer combinations for {label}: {hypothetical_num_comb_str}')
 
-    # split torsions_sampling_points into two lists, use combinations only for those with multiple sampling points
+    # split combined_sampling_points into two lists, use combinations only for those with multiple sampling points
     single_tors, multiple_tors, single_sampling_point, multiple_sampling_points = list(), list(), list(), list()
     multiple_sampling_points_dict = dict()  # used for plotting an energy "scan"
-    for tor, points in torsions_sampling_points.items():
+    for tor, points in combined_sampling_points.items():
         if len(points) == 1:
             single_tors.append(tor)
             single_sampling_point.append((points[0]))
@@ -541,7 +549,8 @@ def conformers_combinations_by_lowest_conformer(label, mol, base_xyz, multiple_t
                                  'FF energy': round(energy, 3),
                                  'source': f'Changing dihedrals on most stable conformer, iteration {i}',
                                  'torsion': tor,
-                                 'dihedral': round(dihedral, 2),
+                                 'dihedral': round(dihedral, 2) if isinstance(dihedral, float)
+                                             else [round(angle, 2) for angle in dihedral],
                                  'dmat': dmat1,
                                  'fl_distance': fl_distance1}
                     newest_conformers_dict[tor].append(conformer)
@@ -557,7 +566,7 @@ def conformers_combinations_by_lowest_conformer(label, mol, base_xyz, multiple_t
                          'FF energy': None,
                          'source': f'Changing dihedrals on most stable conformer, iteration {i}, but FF energy is None',
                          'torsion': tor,
-                         'dihedral': round(dihedral, 2),
+                         'dihedral':  round(dihedral, 2) if isinstance(dihedral, float) else [round(angle, 2) for angle in dihedral],
                          'dmat': dmat1,
                          'fl_distance': fl_distance1})
         new_conformers.extend(newest_conformer_list)
@@ -611,7 +620,6 @@ def generate_all_combinations(label, mol, base_xyz, multiple_tors, multiple_samp
     # Generate sampling points combinations.
     product_combinations = list(product(*multiple_sampling_points))
     new_conformers = list()
-
     if multiple_tors:
         xyzs, energies = change_dihedrals_and_force_field_it(label, mol, xyz=base_xyz, torsions=multiple_tors,
                                                              new_dihedrals=product_combinations, optimize=True,
@@ -756,12 +764,38 @@ def change_dihedrals_and_force_field_it(label, mol, xyz, torsions, new_dihedrals
         new_dihedrals = [new_dihedrals]
 
     for dihedrals in new_dihedrals:
+        skip_to_next_dihedral = False  # Initialize a flag
         xyz_dihedrals = xyz
         for torsion, dihedral in zip(torsions, dihedrals):
             conf, rd_mol = converter.rdkit_conf_from_mol(mol, xyz_dihedrals)
             if conf is not None:
-                torsion_0_indexed = [tor - 1 for tor in torsion]
-                xyz_dihedrals = converter.set_rdkit_dihedrals(conf, rd_mol, torsion_0_indexed, deg_abs=dihedral)
+                # if isinstance(torsion[0], int):         # Assume torsion is a flat tuple
+                #     torsion_0_indexed = [tor - 1 for tor in torsion]
+                # elif isinstance(torsion[0], tuple):     # Check if torsion is a tuple of tuples
+                #     torsion_0_indexed = [[tor - 0 for tor in sub_torsion] for sub_torsion in torsion]
+                if not isinstance(dihedral, list):
+                    torsion_0_indexed = [tor - 1 for tor in torsion]
+                    xyz_dihedrals = converter.set_rdkit_dihedrals(conf, rd_mol, torsion_0_indexed, deg_abs=dihedral)
+                elif isinstance(dihedral, list):
+                    try:
+                        torsion_0_indexed = [[tor - 0 for tor in sub_torsion] for sub_torsion in torsion]
+                        # Calculate the head and tail
+                        ring_length = len(torsion_0_indexed)
+                        head = torsion_0_indexed[0][0]
+                        for torsion in torsion_0_indexed:
+                            if head == torsion[-1]:
+                                tail = torsion[-2]
+                                break
+                        xyz_dihedrals = converter.set_rdkit_ring_dihedrals(
+                            conf, rd_mol, head, tail, torsion_0_indexed[0:ring_length - 3], dihedral[0:ring_length - 3]
+                        )
+                    except Exception as e:  # Catch exceptions specifically from set_rdkit_ring_dihedrals
+                        print(f"Skipping ring dihedral due to an error: {e}")
+                        skip_to_next_dihedral = True  # Set the flag to skip this dihedral set
+                        break  # Exit the inner loop for this dihedral set
+        if skip_to_next_dihedral:
+            continue  # Skip the rest of this `dihedrals` iteration
+
         xyz_, energy = get_force_field_energies(label, mol=mol, xyz=xyz_dihedrals, optimize=True,
                                                 force_field=force_field, suppress_warning=True)
         if energy and xyz_:
@@ -926,17 +960,14 @@ def determine_puckering(conformers, rings_indexes):
         list: Entries are conformer dictionaries.
     """
     for conformer in conformers:
-        # print(f'conformer:{conformer}')
         if isinstance(conformer['xyz'], str):
             xyz = converter.str_to_xyz(conformer['xyz'])
         else:
             xyz = conformer['xyz']
         if 'puckering' not in conformer or not conformer['puckering']:
             conformer['puckering'] = dict()
             for i, ring in enumerate(rings_indexes):
-                # print(f'torsion:{ring}')
                 theta = vectors.calculate_ring_dihedral_angles(coords=xyz['coords'], ring=ring, index=0)
-                print('line932, ring:',ring)
                 conformer['puckering'][tuple((ring[i%len(ring)], ring[(i + 1)%len(ring)], ring[(i + 2)%len(ring)], ring[(i + 3)%len(ring)]) for i in range(len(ring)))] = theta
                 # z = vectors.calculate_puckering_coord(coords=xyz['coords'], ring=ring, index=0)
                 # conformer['puckering'][tuple(ring)] = z

arc/species/converter.py

@@ -1913,8 +1913,8 @@ def set_rdkit_ring_dihedrals(conf_original, rd_mol, ring_head, ring_tail, torsio
         conf_mod.SetAtomPosition(i, pos)
 
     # Visualize before removing hydrogens
-    print("Before removing hydrogens:")
-    show_3d_molecule(rd_mol_mod, title="Before removing hydrogens")
+    # print("Before removing hydrogens:")
+    # show_3d_molecule(rd_mol_mod, title="Before removing hydrogens")
 
     # Remove hydrogens
     rd_mol_noH = Chem.RemoveHs(rd_mol_mod)
@@ -1934,17 +1934,17 @@ def set_rdkit_ring_dihedrals(conf_original, rd_mol, ring_head, ring_tail, torsio
     rd_mol_noH.AddConformer(conf_noH)
 
     # Visualize after removing hydrogens
-    print("After removing hydrogens:")
-    show_3d_molecule(rd_mol_noH, title="After removing hydrogens")
+    # print("After removing hydrogens:")
+    # show_3d_molecule(rd_mol_noH, title="After removing hydrogens")
 
     # Remove the bond to open the ring
     rd_mol_noH = Chem.RWMol(rd_mol_noH)
     rd_mol_noH.RemoveBond(atom_map[ring_head], atom_map[ring_tail])
     Chem.SanitizeMol(rd_mol_noH)
 
     # Visualize after opening the ring
-    print("After opening the ring:")
-    show_3d_molecule(rd_mol_noH, title="After opening the ring")
+    # print("After opening the ring:")
+    # show_3d_molecule(rd_mol_noH, title="After opening the ring")
 
     # Set the specified dihedral angles
     conf_noH = rd_mol_noH.GetConformer()
@@ -1953,8 +1953,8 @@ def set_rdkit_ring_dihedrals(conf_original, rd_mol, ring_head, ring_tail, torsio
         rdMolTransforms.SetDihedralDeg(conf_noH, *torsion_noH, dihedral)
 
     # Visualize after setting dihedral angles
-    print("After setting dihedral angles:")
-    show_3d_molecule(rd_mol_noH, title="After setting dihedral angles")
+    # print("After setting dihedral angles:")
+    # show_3d_molecule(rd_mol_noH, title="After setting dihedral angles")
 
     # Re-add the bond to close the ring
     rd_mol_noH.AddBond(
@@ -1963,8 +1963,8 @@ def set_rdkit_ring_dihedrals(conf_original, rd_mol, ring_head, ring_tail, torsio
     Chem.SanitizeMol(rd_mol_noH)
 
     # Visualize after closing the ring
-    print("After closing the ring:")
-    show_3d_molecule(rd_mol_noH, title="After closing the ring")
+    # print("After closing the ring:")
+    # show_3d_molecule(rd_mol_noH, title="After closing the ring")
 
     # Optimize the molecule
     uff_ff = AllChem.UFFGetMoleculeForceField(rd_mol_noH)
@@ -1987,8 +1987,8 @@ def set_rdkit_ring_dihedrals(conf_original, rd_mol, ring_head, ring_tail, torsio
     conf_noH = rd_mol_noH.GetConformer()
 
     # Visualize after optimization
-    print("After optimization:")
-    show_3d_molecule(rd_mol_noH, title="After optimization")
+    # print("After optimization:")
+    # show_3d_molecule(rd_mol_noH, title="After optimization")
 
     # Add hydrogens back to the optimized molecule
     rd_mol_opt_H = Chem.AddHs(rd_mol_noH)
@@ -1997,8 +1997,8 @@ def set_rdkit_ring_dihedrals(conf_original, rd_mol, ring_head, ring_tail, torsio
     AllChem.UFFOptimizeMolecule(rd_mol_opt_H)
 
     # Visualize after adding hydrogens back
-    print("After adding hydrogens back:")
-    show_3d_molecule(rd_mol_opt_H, title="After adding hydrogens back")
+    # print("After adding hydrogens back:")
+    # show_3d_molecule(rd_mol_opt_H, title="After adding hydrogens back")
 
     # Extract updated coordinates
     conf_opt_H = rd_mol_opt_H.GetConformer()