diff --git a/arc/plotter.py b/arc/plotter.py index bf08bb5fd7..9d2caa84c4 100644 --- a/arc/plotter.py +++ b/arc/plotter.py @@ -1016,11 +1016,14 @@ 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()) + print("torsions_sampling_points:", torsions_sampling_points) + print("torsion_angles:", torsion_angles) 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,11 @@ def plot_torsion_angles(torsion_angles, fig.dpi = 120 num_comb = 1 for i, torsion in enumerate(torsions): + torsion_key = tuple(torsion) + print(f"Checking torsion key: {torsion_key}") + if torsion_key not in torsion_angles: + print(f"Key {torsion_key} not found in torsion_angles. Skipping...") + 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 +1130,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, diff --git a/arc/species/conformers.py b/arc/species/conformers.py index d999592cc8..bfbabb1ecf 100644 --- a/arc/species/conformers.py +++ b/arc/species/conformers.py @@ -98,7 +98,7 @@ MAX_COMBINATION_ITERATIONS = 25 # A threshold below which all combinations will be generated. Above it just samples of the entire search space. -COMBINATION_THRESHOLD = 1000 +COMBINATION_THRESHOLD = 10 # Consolidation tolerances for Z matrices CONSOLIDATION_TOLS = {'R': 1e-2, 'A': 1e-2, 'D': 1e-2} @@ -352,19 +352,29 @@ 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]) - + print('wells_dict:', wells_dict) 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 + print('ring_sampling_points:', ring_sampling_points) + + 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 +385,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 +551,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 +568,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) @@ -572,6 +583,7 @@ def conformers_combinations_by_lowest_conformer(label, mol, base_xyz, multiple_t if plot_path is not None: logger.info(converter.xyz_to_str(lowest_conf_i['xyz'])) arc.plotter.draw_structure(xyz=lowest_conf_i['xyz']) + print('hhhh') num_comb = arc.plotter.plot_torsion_angles(torsion_angles, multiple_sampling_points_dict, wells_dict=wells_dict, e_conformers=newest_conformers_dict, de_threshold=de_threshold, plot_path=plot_path) @@ -760,8 +772,26 @@ def change_dihedrals_and_force_field_it(label, mol, xyz, torsions, new_dihedrals 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] + print("change_torsion", torsion) + print("change_dihedral", dihedral) + if not isinstance(dihedral, list): + xyz_dihedrals = converter.set_rdkit_dihedrals(conf, rd_mol, torsion_0_indexed, deg_abs=dihedral) + elif isinstance(dihedral, list): + # Calculate the head and tail + ring_length = len(torsion_0_indexed) # Length of the ring + head = torsion_0_indexed[0][0] # First element of the first torsion + # Find the tail + for torsion in torsion_0_indexed: + if head == torsion[-1]: # Check if the head appears in the torsion + print("torsion", torsion) + tail = torsion[-2] # The last element in this torsion becomes the tail candidate + break + print(f'we are here ring twist, head {head}, tail {tail}') + xyz_dihedrals = converter.set_rdkit_ring_dihedrals(conf, rd_mol, head, tail, torsion_0_indexed[0:ring_length-3], dihedral[0:ring_length-3]) 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,7 +956,6 @@ 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: diff --git a/arc/species/vectors.py b/arc/species/vectors.py index d675811a80..89f7eb56a4 100644 --- a/arc/species/vectors.py +++ b/arc/species/vectors.py @@ -252,106 +252,106 @@ def calculate_ring_dihedral_angles(coords: Union[list, tuple, dict], return angles -def calculate_plane_angles(coords: Union[list, tuple, dict], - ring: list, - index: int = 0 - ) -> list: - if isinstance(coords, dict) and 'coords' in coords: - coords = coords['coords'] - if not isinstance(coords, (list, tuple)): - raise TypeError(f'coords must be a list or a tuple, got {type(coords)}') - - coords = np.array(coords, dtype=np.float32) - ring = [atom - index for atom in ring] # Adjusting for zero-indexed - print('ring:', ring) - vectors = [coords[ring[i]] - coords[ring[i - 1]] for i in range(len(ring))] - normals = [np.cross(vectors[i], vectors[(i + 1) % len(vectors)]) for i in range(len(vectors))] - normals = [n / np.linalg.norm(n) if np.linalg.norm(n) != 0 else n for n in normals] - print('normals:', normals) - angles = [] - for i in range(len(normals)): - n1 = normals[i] - n2 = normals[(i + 1) % len(normals)] - dot_product = np.dot(n1, n2) - dot_product = max(min(dot_product, 1.0), -1.0) - print('dot pdt:', dot_product) - angle = np.arccos(dot_product) - angle_deg = np.degrees(angle) - angles.append(angle_deg) - - return angles - - -def calculate_puckering_coord(coords: Union[list, tuple, dict], - ring: list, - index: int = 0, - units: str = 'degs', - ) -> float: - """ - Calculate a dihedral angle. - - Args: - coords (list, tuple, dict): The array-format or tuple-format coordinates, or the xyz dict. - torsion (list): The 4 atoms defining the dihedral angle. - index (int, optional): Whether ``torsion`` is 0-indexed or 1-indexed (values are 0 or 1). - units (str, optional): The desired units, either 'rads' for radians, or 'degs' for degrees. - - Raises: - VectorsError: If ``index`` is out of range, or ``torsion`` is of wrong length or has repeating indices. - TypeError: If ``coords`` is of wrong type. - - Returns: float - The dihedral angle in a 0-360 degrees range. - """ - if isinstance(coords, dict) and 'coords' in coords: - coords = coords['coords'] - if not isinstance(coords, (list, tuple)): - raise TypeError(f'coords must be a list or a tuple, got\n{coords}\nwhich is a {type(coords)}') - if index not in [0, 1]: - raise VectorsError(f'index must be either 0 or 1, got {index}') - if ring is None: - raise VectorsError(f'torsion cannot be None') - # if len(ring) != 4: - # raise VectorsError(f'torsion atom list must be of length four, got {len(ring)}') - # if len(set(ring)) < 4: - # raise VectorsError(f'some indices in torsion are repetitive: {ring}') - new_ring = list() - for atoms in ring: - if isinstance(atoms, str) and 'X' in atom: - new_ring.append(int(atom[1:])) - else: - new_ring.append(atoms) - # if not all([isinstance(t, int) for t in new_ring]): - # raise VectorsError(f'all entries in torsion must be integers, got: {new_ring} ' - # f'({[type(t) for t in new_ring]})') - print('ring:',ring) - new_ring = [t - index for t in ring] # convert 1-index to 0-index if needed - print('new_ring:',new_ring) - coords = np.asarray(coords, dtype=np.float32) - center = 0 - for i, atom in enumerate(new_ring): - v1 = coords[new_ring[i]] - center += v1 +# def calculate_plane_angles(coords: Union[list, tuple, dict], +# ring: list, +# index: int = 0 +# ) -> list: +# if isinstance(coords, dict) and 'coords' in coords: +# coords = coords['coords'] +# if not isinstance(coords, (list, tuple)): +# raise TypeError(f'coords must be a list or a tuple, got {type(coords)}') + +# coords = np.array(coords, dtype=np.float32) +# ring = [atom - index for atom in ring] # Adjusting for zero-indexed +# print('ring:', ring) +# vectors = [coords[ring[i]] - coords[ring[i - 1]] for i in range(len(ring))] +# normals = [np.cross(vectors[i], vectors[(i + 1) % len(vectors)]) for i in range(len(vectors))] +# normals = [n / np.linalg.norm(n) if np.linalg.norm(n) != 0 else n for n in normals] +# print('normals:', normals) +# angles = [] +# for i in range(len(normals)): +# n1 = normals[i] +# n2 = normals[(i + 1) % len(normals)] +# dot_product = np.dot(n1, n2) +# dot_product = max(min(dot_product, 1.0), -1.0) +# print('dot pdt:', dot_product) +# angle = np.arccos(dot_product) +# angle_deg = np.degrees(angle) +# angles.append(angle_deg) + +# return angles + + +# def calculate_puckering_coord(coords: Union[list, tuple, dict], +# ring: list, +# index: int = 0, +# units: str = 'degs', +# ) -> float: +# """ +# Calculate a dihedral angle. + +# Args: +# coords (list, tuple, dict): The array-format or tuple-format coordinates, or the xyz dict. +# torsion (list): The 4 atoms defining the dihedral angle. +# index (int, optional): Whether ``torsion`` is 0-indexed or 1-indexed (values are 0 or 1). +# units (str, optional): The desired units, either 'rads' for radians, or 'degs' for degrees. + +# Raises: +# VectorsError: If ``index`` is out of range, or ``torsion`` is of wrong length or has repeating indices. +# TypeError: If ``coords`` is of wrong type. + +# Returns: float +# The dihedral angle in a 0-360 degrees range. +# """ +# if isinstance(coords, dict) and 'coords' in coords: +# coords = coords['coords'] +# if not isinstance(coords, (list, tuple)): +# raise TypeError(f'coords must be a list or a tuple, got\n{coords}\nwhich is a {type(coords)}') +# if index not in [0, 1]: +# raise VectorsError(f'index must be either 0 or 1, got {index}') +# if ring is None: +# raise VectorsError(f'torsion cannot be None') +# # if len(ring) != 4: +# # raise VectorsError(f'torsion atom list must be of length four, got {len(ring)}') +# # if len(set(ring)) < 4: +# # raise VectorsError(f'some indices in torsion are repetitive: {ring}') +# new_ring = list() +# for atoms in ring: +# if isinstance(atoms, str) and 'X' in atom: +# new_ring.append(int(atom[1:])) +# else: +# new_ring.append(atoms) +# # if not all([isinstance(t, int) for t in new_ring]): +# # raise VectorsError(f'all entries in torsion must be integers, got: {new_ring} ' +# # f'({[type(t) for t in new_ring]})') +# print('ring:',ring) +# new_ring = [t - index for t in ring] # convert 1-index to 0-index if needed +# print('new_ring:',new_ring) +# coords = np.asarray(coords, dtype=np.float32) +# center = 0 +# for i, atom in enumerate(new_ring): +# v1 = coords[new_ring[i]] +# center += v1 - center = center/len(new_ring) - N = len(new_ring) - R_list = coords - center - print('R_list:',R_list) - R_prime = np.zeros(3) - R_double_prime = np.zeros(3) - - for j, Rj in enumerate(new_ring, start=1): - angle = 2 * np.pi * (j - 1) / N - R_prime += np.sin(angle) * R_list[Rj] - R_double_prime += np.cos(angle) * R_list[Rj] - - # Calculate the unit normal vector to the mean plane - n = np.cross(R_prime, R_double_prime) - n /= np.linalg.norm(n) - - # Compute displacements from the mean plane - displacements = [np.dot(R_list[Rj], n) for Rj in new_ring] - return displacements +# center = center/len(new_ring) +# N = len(new_ring) +# R_list = coords - center +# print('R_list:',R_list) +# R_prime = np.zeros(3) +# R_double_prime = np.zeros(3) + +# for j, Rj in enumerate(new_ring, start=1): +# angle = 2 * np.pi * (j - 1) / N +# R_prime += np.sin(angle) * R_list[Rj] +# R_double_prime += np.cos(angle) * R_list[Rj] + +# # Calculate the unit normal vector to the mean plane +# n = np.cross(R_prime, R_double_prime) +# n /= np.linalg.norm(n) + +# # Compute displacements from the mean plane +# displacements = [np.dot(R_list[Rj], n) for Rj in new_ring] +# return displacements def calculate_param(coords: Union[list, tuple, dict],