Skip to content

Commit

Permalink
Merge pull request #425 from ReactionMechanismGenerator/scan_fixes
Browse files Browse the repository at this point in the history 
Scan fixes
  • Loading branch information
@alongd
alongd committed Sep 27, 2020
2 parents 648d74f + cf2e940 commit fba1e60
Show file tree
Hide file tree
Showing 4 changed files with 79 additions and 64 deletions.
23 changes: 23 additions & 0 deletions arc/common.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1054,3 +1054,26 @@ def get_ordered_intersection_of_two_lists(l1: list,
l3 = [v for v in l3 if v not in lookup and lookup.add(v) is None]

return l3


def get_angle_in_180_range(angle: float,
round_to: Optional[int] = 2,
) -> float:
"""
Get the corresponding angle in the -180 to +180 degree range.
Args:
angle (float): An angle in degrees.
round_to (int, optional): The number of decimal figures to round the result to.
``None`` to not round. Default: 2.
Returns:
float: The corresponding angle in the -180 to +180 degree range.
"""
angle = float(angle)
while not (-180 <= angle <= 180):
factor = 360 if angle < -180 else -360
angle += factor
if round_to is not None:
return round(angle, round_to)
return angle
19 changes: 18 additions & 1 deletion arc/commonTest.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -592,7 +592,7 @@ def test_is_same_sequence_sublist(self):
self.assertFalse(common.is_same_sequence_sublist([4, 3, 2], [1, 2, 3, 4]))

def test_get_ordered_intersection_of_two_lists(self):
"""Test get ordered intersection of two lists."""
"""Test get ordered intersection of two lists"""
l1 = [1, 2, 3, 3, 5, 6]
l2 = [6, 3, 5, 5, 1]

Expand All @@ -616,6 +616,23 @@ def test_get_ordered_intersection_of_two_lists(self):
l2 = [7]
self.assertEqual(common.get_ordered_intersection_of_two_lists(l1, l2), list())

def test_get_angle_in_180_range(self):
"""Test the getting a corresponding angle in the -180 to +180 range"""
self.assertEqual(common.get_angle_in_180_range(0), 0)
self.assertEqual(common.get_angle_in_180_range(10), 10)
self.assertEqual(common.get_angle_in_180_range(-5.364589, round_to=None), -5.364589)
self.assertEqual(common.get_angle_in_180_range(-5.364589), -5.36)
self.assertEqual(common.get_angle_in_180_range(-120), -120)
self.assertEqual(common.get_angle_in_180_range(179.999), 180)
self.assertEqual(common.get_angle_in_180_range(180), 180)
self.assertEqual(common.get_angle_in_180_range(-180), -180)
self.assertEqual(common.get_angle_in_180_range(181), -179)
self.assertEqual(common.get_angle_in_180_range(360), 0)
self.assertEqual(common.get_angle_in_180_range(362), 2)
self.assertEqual(common.get_angle_in_180_range(1000), -80)
self.assertEqual(common.get_angle_in_180_range(-1000), 80)
self.assertEqual(common.get_angle_in_180_range(247.62), -112.38)

@classmethod
def tearDownClass(cls):
"""
Expand Down
60 changes: 17 additions & 43 deletions arc/plotter.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@
from rmgpy.species import Species

from arc.common import (extermum_list,
get_angle_in_180_range,
get_logger,
is_notebook,
save_yaml_file,
Expand Down Expand Up @@ -190,7 +191,7 @@ def plot_3d_mol_as_scatter(xyz, path=None, plot_h=True, show_plot=True, name='',
ax.scatter(xs=x, ys=y, zs=z, s=sizes, c=colors, depthshade=True)
for i, symbol in enumerate(symbols):
text = symbol if index is None else symbol + ' ' + str(i + index)
ax.text(x[i]+0.01, y[i]+0.01, z[i]+0.01, text, size=10)
ax.text(x[i] + 0.01, y[i] + 0.01, z[i] + 0.01, text, size=10)
plt.axis('off')
if show_plot:
plt.show()
Expand Down Expand Up @@ -363,7 +364,7 @@ def draw_thermo_parity_plots(species_list: list,
thermo_sources = '\nSources of thermoproperties determined by RMG for the parity plots:\n'
max_label_len = max([len(label) for label in labels])
for i, label in enumerate(labels):
thermo_sources += ' {0}: {1}{2}\n'.format(label, ' '*(max_label_len - len(label)), comments[i])
thermo_sources += ' {0}: {1}{2}\n'.format(label, ' ' * (max_label_len - len(label)), comments[i])
logger.info(thermo_sources)
if path is not None:
with open(os.path.join(path, 'thermo.info'), 'w') as f:
Expand Down Expand Up @@ -556,8 +557,8 @@ def text_plotter(x_data, y_data, labels, text_positions, axis, txt_width, txt_he
for x, y, lab, t in zip(x_data, y_data, labels, text_positions):
axis.text(x - .03, 1.02 * t, f'{lab}', rotation=0, color='black', fontsize=10)
if y != t:
axis.arrow(x, t + 20, 0, y-t, color='blue', alpha=0.2, width=txt_width*0.0,
head_width=.02, head_length=txt_height*0.5,
axis.arrow(x, t + 20, 0, y - t, color='blue', alpha=0.2, width=txt_width * 0.0,
head_width=.02, head_length=txt_height * 0.5,
zorder=0, length_includes_head=True)


Expand Down Expand Up @@ -903,10 +904,10 @@ def plot_torsion_angles(torsion_angles, torsions_sampling_points=None, wells_dic
axs.frameon = False
axs.set_ylabel(str(torsion), labelpad=10)
axs.set_yticklabels(['' for _ in range(len(torsions))])
axs.tick_params(axis='y', # changes apply to the x-axis
which='both', # both major and minor ticks are affected
left=False, # ticks along the bottom edge are off
right=False, # ticks along the top edge are off
axs.tick_params(axis='y', # changes apply to the x-axis
which='both', # both major and minor ticks are affected
left=False, # ticks along the bottom edge are off
right=False, # ticks along the top edge are off
labelleft=False) # labels along the bottom edge are off
axs.set_title('Dihedral angle (degrees)')
axs.axes.xaxis.set_ticks(ticks=ticks)
Expand Down Expand Up @@ -955,10 +956,10 @@ def plot_torsion_angles(torsion_angles, torsions_sampling_points=None, wells_dic
# axs[i].yaxis.label.set_rotation(0)
if e_conformers is None:
axs[i].set_yticklabels(['' for _ in range(len(torsions))])
axs[i].tick_params(axis='y', # changes apply to the x-axis
which='both', # both major and minor ticks are affected
left=False, # ticks along the bottom edge are off
right=False, # ticks along the top edge are off
axs[i].tick_params(axis='y', # changes apply to the x-axis
which='both', # both major and minor ticks are affected
left=False, # ticks along the bottom edge are off
right=False, # ticks along the top edge are off
labelleft=False) # labels along the bottom edge are off
axs[0].set_title('Dihedral angle (degrees)')
# Hide x labels and tick labels for all but bottom plot.
Expand Down Expand Up @@ -1139,38 +1140,13 @@ def plot_2d_rotor_scan(results, path=None, label='', cmap='Blues', resolution=90
phis1 = np.append(phis1, phis1[-1] + res1)
zero_phi0, zero_phi1 = list(), list()
energies = np.zeros(shape=(phis0.size, phis1.size), dtype=np.float64)
keys_list = list(results['directed_scan'].keys())
e_min = None
for i, phi0 in enumerate(phis0):
for j, phi1 in enumerate(phis1):
key = tuple(f'{dihedral:.2f}' for dihedral in [phi0, phi1])
if key in keys_list:
energies[i, j] = results['directed_scan'][key]['energy']
if e_min is None or energies[i, j] < e_min:
e_min = energies[i, j]
else:
keys = list()
keys.append(tuple(f'{dihedral:.2f}' for dihedral in [360.0 - phi0, phi1]))
keys.append(tuple(f'{dihedral:.2f}' for dihedral in [phi0, 360.0 - phi1]))
keys.append(tuple(f'{dihedral:.2f}' for dihedral in [360.0 + phi0, phi1]))
keys.append(tuple(f'{dihedral:.2f}' for dihedral in [phi0, 360.0 + phi1]))
keys.append(tuple(f'{dihedral:.2f}' for dihedral in [360.0 + phi0, 360.0 + phi1]))
keys.append(tuple(f'{dihedral:.2f}' for dihedral in [360.0 - phi0, 360.0 - phi1]))
keys.append(tuple(f'{dihedral:.2f}' for dihedral in [-phi0, phi1]))
keys.append(tuple(f'{dihedral:.2f}' for dihedral in [phi0, -phi1]))
keys.append(tuple(f'{dihedral:.2f}' for dihedral in [-phi0, -phi1]))
keys.append(tuple(f'{dihedral:.2f}' for dihedral in [phi0 - 360.0, phi1]))
keys.append(tuple(f'{dihedral:.2f}' for dihedral in [phi0, phi1 - 360.0]))
keys.append(tuple(f'{dihedral:.2f}' for dihedral in [phi0 - 360.0, phi1 - 360.0]))
keys.append(tuple(f'{dihedral:.2f}' for dihedral in [phi0 - 360.0, phi1 + 360.0]))
keys.append(tuple(f'{dihedral:.2f}' for dihedral in [phi0 + 360.0, phi1 - 360.0]))
for key_ in keys:
if key_ in keys_list:
energies[i, j] = results['directed_scan'][key_]['energy']
break
else:
logger.warning(f'Could not replace {key} when plotting 2D scan for {label} ({path}).\n'
f'Tried: {keys}.')
key = tuple(f'{get_angle_in_180_range(dihedral):.2f}' for dihedral in [phi0, phi1])
energies[i, j] = results['directed_scan'][key]['energy']
if e_min is None or energies[i, j] < e_min:
e_min = energies[i, j]
if mark_lowest_conformations and energies[i, j] == 0:
zero_phi0.append(phi0)
zero_phi1.append(phi1)
Expand All @@ -1190,11 +1166,9 @@ def plot_2d_rotor_scan(results, path=None, label='', cmap='Blues', resolution=90
plt.ylabel(f'Dihedral 2 for {results["scans"][1]} (degrees)')
label = ' for ' + label if label else ''
plt.title(f'2D scan energies (kJ/mol){label}')
# min_x = int(np.ceil(np.min(phis0) / 10.0)) * 10
min_x = min_y = -180
plt.xlim = (min_x, min_x + 360)
plt.xticks(np.arange(min_x, min_x + 361, step=60))
# min_y = int(np.ceil(np.min(phis1) / 10.0)) * 10
plt.ylim = (min_y, min_y + 360)
plt.yticks(np.arange(min_y, min_y + 361, step=60))

Expand Down
41 changes: 21 additions & 20 deletions arc/scheduler.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@

from arc import parser, plotter
from arc.common import (extermum_list,
get_angle_in_180_range,
get_logger,
get_ordinal_indicator,
read_yaml_file,
Expand Down Expand Up @@ -219,7 +220,7 @@ def __init__(self,
restart_dict: Optional[dict] = None,
max_job_time: Optional[float] = None,
allow_nonisomorphic_2d: Optional[bool] = False,
memory: Optional[float] = None,
memory: Optional[float] = None,
testing: Optional[bool] = False,
dont_gen_confs: Optional[list] = None,
n_confs: Optional[int] = 10,
Expand Down Expand Up @@ -419,7 +420,7 @@ def __init__(self,
# composite is done; do other jobs
if not self.output[species.label]['job_types']['freq'] \
and 'freq' not in list(self.job_dict[species.label].keys()) \
and (self.species_dict[species.label].is_ts
and (self.species_dict[species.label].is_ts
or self.species_dict[species.label].number_of_atoms > 1):
self.run_freq_job(species.label)
if self.job_types['rotors']:
Expand Down Expand Up @@ -1088,7 +1089,7 @@ def run_sp_job(self,

else:
raise RuntimeError(f'Unable to set the path for the sp job for species {label}')

return
if 'sp' not in self.job_dict[label]: # Check whether or not single point jobs have been spawned yet
# we're spawning the first sp job for this species
Expand Down Expand Up @@ -1197,7 +1198,7 @@ def run_scan_jobs(self, label):
else:
self.spawn_directed_scan_jobs(label, rotor_index=i)
else:
self.spawn_directed_scan_jobs(label, rotor_index=i)
self.spawn_directed_scan_jobs(label, rotor_index=i)
else:
# this is a "normal" scan (not directed)
# check this job isn't already running on the server(from a restarted project)
Expand Down Expand Up @@ -1424,12 +1425,13 @@ def spawn_directed_scan_jobs(self, label, rotor_index, xyz=None):
elif 'brute' in directed_scan_type:
# spawn jobs all at once
dihedrals = dict()

for scan in scans:
original_dihedral = calculate_dihedral_angle(coords=xyz['coords'], torsion=scan)
dihedrals[tuple(scan)] = [round(original_dihedral + i * increment
if original_dihedral + i * increment <= 180.0
else original_dihedral + i * increment - 360.0, 2)
for i in range(int(360 / increment) + 1)]
original_dihedral = get_angle_in_180_range(calculate_dihedral_angle(coords=xyz['coords'],
torsion=scan,
index=1))
dihedrals[tuple(scan)] = [get_angle_in_180_range(original_dihedral + i * increment) for i in
range(int(360 / increment) + 1)]
modified_xyz = xyz
if 'diagonal' not in directed_scan_type:
# increment dihedrals one by one (resulting in an ND scan)
Expand Down Expand Up @@ -1484,8 +1486,7 @@ def spawn_directed_scan_jobs(self, label, rotor_index, xyz=None):
max_num = 360 / increment + 1 # dihedral angles per scan
original_dihedrals = list()
for dihedral in rotor_dict['original_dihedrals']:
f_dihedral = float(dihedral)
original_dihedrals.append(f_dihedral if f_dihedral < 180.0 else f_dihedral - 360.0)
original_dihedrals.append(get_angle_in_180_range(dihedral))
if not any(self.species_dict[label].rotors_dict[rotor_index]['cont_indices']):
# this is the first call for this cont_opt directed rotor, spawn the first job w/o changing dihedrals
self.run_job(label=label,
Expand Down Expand Up @@ -1516,10 +1517,10 @@ def spawn_directed_scan_jobs(self, label, rotor_index, xyz=None):
dihedrals = list()
for index, (original_dihedral, scan_) in enumerate(zip(original_dihedrals, scans)):
dihedral = original_dihedral + \
self.species_dict[label].rotors_dict[rotor_index]['cont_indices'][index] * increment
self.species_dict[label].rotors_dict[rotor_index]['cont_indices'][index] * increment
# change the original dihedral so we won't end up with two calcs for 180.0, but none for -180.0
# (it only matters for plotting, the geometry is of course the same)
dihedral = dihedral if dihedral <= 180.0 else dihedral - 360.0
dihedral = get_angle_in_180_range(dihedral)
dihedrals.append(dihedral)
# Only change the dihedrals in the xyz if this torsion corresponds to the current index,
# or if this is a diagonal scan.
Expand Down Expand Up @@ -1549,7 +1550,7 @@ def spawn_directed_scan_jobs(self, label, rotor_index, xyz=None):
self.species_dict[label].rotors_dict[rotor_index]['cont_indices'][index] += 1
break
elif (self.species_dict[label].rotors_dict[rotor_index]['cont_indices'][index] == max_num - 1
and index < len(scans) - 1):
and index < len(scans) - 1):
self.species_dict[label].rotors_dict[rotor_index]['cont_indices'][index] = 0

def spawn_md_jobs(self, label, prev_conf_list=None, num_confs=None):
Expand Down Expand Up @@ -2417,7 +2418,7 @@ def check_scan_job(self, label: str, job: Job) -> None:
pivots=job.pivots,
energies=energies,
scan_res=job.scan_res,
used_methods=self.species_dict[label].rotors_dict[i]['trsh_methods'],
used_methods=self.species_dict[label].rotors_dict[i]['trsh_methods'],
log_file=job.local_path_to_output_file,
species=self.species_dict[label],
preserve_params=self.species_dict[label].preserve_param_in_scan,
Expand Down Expand Up @@ -2660,7 +2661,7 @@ def check_md_job(self, label, job, max_iterations=10):
# spawn a new MD simulation
ordinal = get_ordinal_indicator(self.species_dict[label].recent_md_conformer[2] + 1)
logger.info(f'{self.species_dict[label].recent_md_conformer[2] + 1}{ordinal} conformer for '
f'{label}:\n{ lowest_conf[0]}')
f'{label}:\n{lowest_conf[0]}')
plotter.draw_structure(xyz=lowest_conf[0], species=self.species_dict[label])
ordinal = get_ordinal_indicator(self.species_dict[label].recent_md_conformer[2] + 2)
logger.info(f'Spawning the {self.species_dict[label].recent_md_conformer[2] + 2}{ordinal} round of MD '
Expand Down Expand Up @@ -2704,8 +2705,8 @@ def check_all_done(self, label):
if any([job_type not in ['conformers', 'opt', 'composite']
for job_type in self.job_dict[label].keys()]) else zero_delta
self.species_dict[label].run_time = self.species_dict[label].run_time \
or (conf_time or zero_delta) + (opt_time or zero_delta) \
+ (comp_time or zero_delta) + (other_time or zero_delta)
or (conf_time or zero_delta) + (opt_time or zero_delta) \
+ (comp_time or zero_delta) + (other_time or zero_delta)
logger.info(f'\nAll jobs for species {label} successfully converged. '
f'Run time: {self.species_dict[label].run_time}')
else:
Expand Down Expand Up @@ -2968,7 +2969,7 @@ def troubleshoot_ess(self,
self.species_dict[label].checkfile = job.checkfile
# determine if the species is a hydrogen (or its isotope) atom
is_h = self.species_dict[label].number_of_atoms == 1 and \
self.species_dict[label].mol.atoms[0].element.symbol in ['H', 'D', 'T']
self.species_dict[label].mol.atoms[0].element.symbol in ['H', 'D', 'T']
output_errors, ess_trsh_methods, remove_checkfile, level_of_theory, \
software, job_type, fine, trsh_keyword, memory, shift, cpu_cores, dont_rerun = \
trsh_ess_job(label=label,
Expand Down Expand Up @@ -3009,7 +3010,7 @@ def troubleshoot_ess(self,
directed_dihedrals=job.directed_dihedrals,
directed_scans=job.directed_scans,
directed_scan_type=job.directed_scan_type,
rotor_index = job.rotor_index,
rotor_index=job.rotor_index,
cpu_cores=cpu_cores,
)
self.save_restart_dict()
Expand Down

0 comments on commit fba1e60

Please sign in to comment.