Extract unique puckering mode and make plots

Author: JintaoWu98

arc/plotter.py

@@ -1014,12 +1014,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:
@@ -1048,6 +1051,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:
@@ -1121,6 +1126,76 @@ def plot_torsion_angles(torsion_angles,
     return num_comb
 
 
+def plot_ring_torsion_angles(conformers, plot_path=None, tolerance=15):
+    """
+    Plot the torsion angles of the generated conformers for each ring,
+    considering torsion similarity within a given tolerance.
+
+    Args:
+        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.
+    """
+    if 'puckering' not in conformers[0]:
+        return
+
+    # Dictionary to store unique angle sets for each ring
+    ring_angle_data = {}
+    
+    # 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))
+
+    
+    # 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)
+    
+    return ring_angle_data
+
+
 def plot_1d_rotor_scan(angles: Optional[Union[list, tuple, np.array]] = None,
                        energies: Optional[Union[list, tuple, np.array]] = None,
                        results: Optional[dict] = None,

arc/species/conformers.py

@@ -348,18 +348,25 @@ 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)
+        angles_sorted = arc.plotter.plot_ring_torsion_angles(conformers=conformers, plot_path=plot_path)
+    # Process the ring angle data
+    if angles_sorted is not None:
+        for ring, angle_counts in angles_sorted.items():
+            angles = [list(angle) for angle, _ in angle_counts]
+            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)
@@ -370,10 +377,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]))
@@ -536,7 +543,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)
@@ -552,7 +560,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)