diff --git a/arc/plotter.py b/arc/plotter.py index bf08bb5fd7..70665717ae 100644 --- a/arc/plotter.py +++ b/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, diff --git a/arc/species/conformers.py b/arc/species/conformers.py index d999592cc8..419ae03728 100644 --- a/arc/species/conformers.py +++ b/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,7 +960,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: @@ -934,9 +967,7 @@ def determine_puckering(conformers, rings_indexes): 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 diff --git a/arc/species/converter.py b/arc/species/converter.py index 4659c7ea65..a8dd861ed1 100644 --- a/arc/species/converter.py +++ b/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,8 +1934,8 @@ 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) @@ -1943,8 +1943,8 @@ def set_rdkit_ring_dihedrals(conf_original, rd_mol, ring_head, ring_tail, torsio 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() diff --git a/arc/species/vectors.py b/arc/species/vectors.py index d675811a80..eaa1ea628d 100644 --- a/arc/species/vectors.py +++ b/arc/species/vectors.py @@ -243,7 +243,6 @@ def calculate_ring_dihedral_angles(coords: Union[list, tuple, dict], coords = np.array(coords, dtype=np.float32) ring = [atom - index for atom in ring] # Adjusting for zero-indexed - print('ring:', ring) angles = [] for i in range(len(ring)): angle_deg = calculate_dihedral_angle(coords, [ring[i%len(ring)], ring[(i + 1)%len(ring)], ring[(i + 2)%len(ring)], ring[(i + 3)%len(ring)]]) @@ -252,106 +251,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],