From 5f8bca9922e2c46e9d39f6637cf52384b333eea7 Mon Sep 17 00:00:00 2001 From: Jintao Date: Sun, 8 Dec 2024 11:19:59 +0800 Subject: [PATCH] Extract unique puckering mode and make plots --- arc/plotter.py | 79 ++++++++++++++++++++++++++++++++++++++- arc/species/conformers.py | 21 +++++++---- 2 files changed, 91 insertions(+), 9 deletions(-) diff --git a/arc/plotter.py b/arc/plotter.py index 316f4d3058..ce75e6bb02 100644 --- a/arc/plotter.py +++ b/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, diff --git a/arc/species/conformers.py b/arc/species/conformers.py index 8279e682a2..dbcf630c0d 100644 --- a/arc/species/conformers.py +++ b/arc/species/conformers.py @@ -348,18 +348,24 @@ 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 + 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 +376,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 +542,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 +559,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)