Skip to content

Commit

Permalink
Extract unique puckering mode and make plots
Browse files Browse the repository at this point in the history
  • Loading branch information
JintaoWu98 committed Dec 8, 2024
1 parent b4b7b08 commit 5f8bca9
Show file tree
Hide file tree
Showing 2 changed files with 91 additions and 9 deletions.
79 changes: 77 additions & 2 deletions arc/plotter.py
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -1121,6 +1126,76 @@ def plot_torsion_angles(torsion_angles,
return num_comb


def plot_ring_torsion_angles(conformers, plot_path=None, tolerance=15):

Check notice

Code scanning / CodeQL

Explicit returns mixed with implicit (fall through) returns Note

Mixing implicit and explicit returns may indicate an error as implicit returns always return None.
"""
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:

Check notice

Code scanning / CodeQL

Empty except Note

'except' clause does nothing but pass and there is no explanatory comment.
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,
Expand Down
21 changes: 14 additions & 7 deletions arc/species/conformers.py
Original file line number Diff line number Diff line change
Expand Up @@ -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():

Check failure

Code scanning / CodeQL

Potentially uninitialized local variable Error

Local variable 'angles_sorted' may be used before it is initialized.
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)
Expand All @@ -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]))
Expand Down Expand Up @@ -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)
Expand All @@ -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)
Expand Down

0 comments on commit 5f8bca9

Please sign in to comment.