Troubleshoot rotor scans

arc/job/inputs.py

@@ -42,7 +42,7 @@
 
 {charge} {multiplicity}
 {xyz}
-{scan}
+{scan}(scan_trsh}
 
 
 """,

arc/job/job.py

@@ -41,6 +41,7 @@ class Job(object):
  (e.g., "2 1 3 5" as a string or [2, 1, 3, 5] as a list of integers)
  `pivots` ``list`` The rotor scan pivots, if the job type is scan. Not used directly in these
  methods, but used to identify the rotor.
+ `scan_res` ``int`` The rotor scan resolution in degrees
  `software` ``str`` The electronic structure software to be used
  `server_nodes` ``list`` A list of nodes this job was submitted to (for troubleshooting)
  `memory` ``int`` The allocated memory (1500 mb by default)
@@ -63,8 +64,9 @@ class Job(object):
  `server` ``str`` Server's name. Determined automatically
  'trsh' ''str'' A troubleshooting handle to be appended to input files
  'ess_trsh_methods' ``list`` A list of troubleshooting methods already tried out for ESS convergence
- `occ` ``int`` The number of occupied orbitals (core + val) from a molpro CCSD sp calc
  `initial_trsh` ``dict`` Troubleshooting methods to try by default. Keys are server names, values are trshs
+ `scan_trsh` ``str`` A troubleshooting method for rotor scans
+ `occ` ``int`` The number of occupied orbitals (core + val) from a molpro CCSD sp calc
  `project_directory` ``str`` The path to the project directory
  `max_job_time` ``int`` The maximal allowed job time on the server in hours
  ================ =================== ===============================================================================
@@ -76,9 +78,9 @@ class Job(object):
  """
  def __init__(self, project, settings, species_name, xyz, job_type, level_of_theory, multiplicity, project_directory,
  charge=0, conformer=-1, fine=False, shift='', software=None, is_ts=False, scan='', pivots=None,
- memory=1500, comments='', trsh='', ess_trsh_methods=None, occ=None, initial_trsh=None, job_num=None,
- job_server_name=None, job_name=None, job_id=None, server=None, date_time=None, run_time=None,
- max_job_time=5):
+ memory=1500, comments='', trsh='', scan_trsh='', ess_trsh_methods=None, initial_trsh=None, job_num=None,
+ job_server_name=None, job_name=None, job_id=None, server=None, date_time=None, run_time=None, occ=None,
+ max_job_time=120, scan_res=None):
  self.project = project
  self.settings=settings
  self.date_time = date_time if date_time is not None else datetime.datetime.now()
@@ -94,6 +96,8 @@ def __init__(self, project, settings, species_name, xyz, job_type, level_of_theo
  self.ess_trsh_methods = ess_trsh_methods if ess_trsh_methods is not None else list()
  self.trsh = trsh
  self.initial_trsh = initial_trsh if initial_trsh is not None else dict()
+ self.scan_trsh = scan_trsh
+ self.scan_res = scan_res if scan_res is not None else rotor_scan_resolution
  self.max_job_time = max_job_time
  job_types = ['conformer', 'opt', 'freq', 'optfreq', 'sp', 'composite', 'scan', 'gsm', 'irc', 'ts_guess']
  # the 'conformer' job type is identical to 'opt', but we differentiate them to be identifiable in Scheduler
@@ -319,7 +323,7 @@ def write_completed_job_to_csv_file(self):
  def write_submit_script(self):
  un = servers[self.server]['un'] # user name
  if self.max_job_time > 9999 or self.max_job_time == 0:
- self.max_job_time = 9999
+ self.max_job_time = 120
  if t_max_format[servers[self.server]['cluster_soft']] == 'days':
  # e.g., 5-0:00:00
  d, h = divmod(self.max_job_time, 24)
@@ -400,7 +404,8 @@ def write_input_file(self):
  else:
  job_type_1 = 'opt=(calcfc, noeigentest)'
  if self.fine:
- fine = 'scf=(tight, direct) int=(ultrafine, Acc2E=12)'
+ # Note that the Acc2E argument is not available in Gaussian03
+ fine = 'scf=(tight, direct) integral=(grid=ultrafine, Acc2E=12)'
  if self.is_ts:
  job_type_1 = 'opt=(ts, calcfc, noeigentest, tight)'
  else:
@@ -425,7 +430,7 @@ def write_input_file(self):
 
  elif self.job_type == 'freq':
  if self.software == 'gaussian':
- job_type_2 = 'freq iop(7/33=1) scf=(tight, direct) int=(ultrafine, Acc2E=12)'
+ job_type_2 = 'freq iop(7/33=1) scf=(tight, direct) integral=(grid=ultrafine, Acc2E=12)'
  elif self.software == 'qchem':
  job_type_1 = 'freq'
  elif self.software == 'molpro':
@@ -439,7 +444,7 @@ def write_input_file(self):
  job_type_1 = 'opt=(calcfc, noeigentest)'
  job_type_2 = 'freq iop(7/33=1)'
  if self.fine:
- fine = 'scf=(tight, direct) int=(ultrafine, Acc2E=12)'
+ fine = 'scf=(tight, direct) integral=(grid=ultrafine, Acc2E=12)'
  if self.is_ts:
  job_type_1 = 'opt=(ts, calcfc, noeigentest, tight)'
  else:
@@ -476,7 +481,7 @@ def write_input_file(self):
 
  if self.job_type == 'sp':
  if self.software == 'gaussian':
- job_type_1 = 'scf=(tight, direct) int=(ultrafine, Acc2E=12)'
+ job_type_1 = 'scf=(tight, direct) integral=(grid=ultrafine, Acc2E=12)'
  elif self.software == 'qchem':
  job_type_1 = 'sp'
  elif self.software == 'molpro':
@@ -485,7 +490,7 @@ def write_input_file(self):
  if self.job_type == 'composite':
  if self.software == 'gaussian':
  if self.fine:
- fine = 'scf=(tight, direct) int=(ultrafine, Acc2E=12)'
+ fine = 'scf=(tight, direct) integral=(grid=ultrafine, Acc2E=12)'
  if self.is_ts:
  job_type_1 = 'opt=(ts, calcfc, noeigentest, tight)'
  else:
@@ -496,14 +501,16 @@ def write_input_file(self):
  if self.job_type == 'scan':
  if self.software == 'gaussian':
  if self.is_ts:
- job_type_1 = 'opt=(ts, modredundant, calcfc, noeigentest)'
+ job_type_1 = 'opt=(ts, modredundant, calcfc, noeigentest, maxStep=5) scf=(tight, direct)' \
+ ' integral=(grid=ultrafine, Acc2E=12)'
  else:
- job_type_1 = 'opt=(modredundant, calcfc, noeigentest)'
+ job_type_1 = 'opt=(modredundant, calcfc, noeigentest, maxStep=5) scf=(tight, direct)' \
+ ' integral=(grid=ultrafine, Acc2E=12)'
  scan_string = ''.join([str(num) + ' ' for num in self.scan])
- if not divmod(360, rotor_scan_resolution):
- raise JobError('Scan job got an illegal rotor scan resolution of {0}'.format(rotor_scan_resolution))
- scan_string = 'D ' + scan_string + 'S ' + str(int(360 / rotor_scan_resolution)) + ' ' +\
- '{0:10}'.format(float(rotor_scan_resolution))
+ if not divmod(360, self.scan_res):
+ raise JobError('Scan job got an illegal rotor scan resolution of {0}'.format(self.scan_res))
+ scan_string = 'D ' + scan_string + 'S ' + str(int(360 / self.scan_res)) + ' ' +\
+ '{0:10}'.format(float(self.scan_res))
  else:
  raise ValueError('Currently rotor scan is only supported in gaussian. Got: {0} using the {1} level of'
  ' theory'.format(self.software, self.method + '/' + self.basis_set))
@@ -538,7 +545,7 @@ def write_input_file(self):
  basis=self.basis_set, charge=self.charge, multiplicity=self.multiplicity,
  spin=self.spin, xyz=self.xyz, job_type_1=job_type_1,
  job_type_2=job_type_2, scan=scan_string, restricted=restricted, fine=fine,
- shift=self.shift, trsh=self.trsh)
+ shift=self.shift, trsh=self.trsh, scan_trsh=self.scan_trsh)
  except KeyError as e:
  logging.error('Could not interpret all input file keys in\n{0}'.format(self.input))
  raise e

arc/scheduler.py

@@ -106,7 +106,7 @@ class Scheduler(object):
  """
  def __init__(self, project, settings, species_list, composite_method, conformer_level, opt_level, freq_level,
  sp_level, scan_level, ts_guess_level, project_directory, rmgdatabase, fine=False, scan_rotors=True,
- generate_conformers=True, initial_trsh=None, rxn_list=None, restart_dict=None, max_job_time=5,
+ generate_conformers=True, initial_trsh=None, rxn_list=None, restart_dict=None, max_job_time=120,
  testing=False):
  self.rmgdb = rmgdatabase
  self.restart_dict = restart_dict
@@ -458,7 +458,7 @@ def schedule_jobs(self):
  self.species_dict[spc.label] = spc
 
  def run_job(self, label, xyz, level_of_theory, job_type, fine=False, software=None, shift='', trsh='', memory=1500,
- conformer=-1, ess_trsh_methods=None, scan='', pivots=None, occ=None):
+ conformer=-1, ess_trsh_methods=None, scan='', pivots=None, occ=None, scan_trsh='', scan_res=None):
  """
  A helper function for running (all) jobs
  """
@@ -471,7 +471,8 @@ def run_job(self, label, xyz, level_of_theory, job_type, fine=False, software=No
  level_of_theory=level_of_theory, multiplicity=species.multiplicity, charge=species.charge, fine=fine,
  shift=shift, software=software, is_ts=species.is_ts, memory=memory, trsh=trsh, conformer=conformer,
  ess_trsh_methods=ess_trsh_methods, scan=scan, pivots=pivots, occ=occ, initial_trsh=self.initial_trsh,
- project_directory=self.project_directory, max_job_time=self.max_job_time)
+ project_directory=self.project_directory, max_job_time=self.max_job_time, scan_trsh=scan_trsh,
+ scan_res=scan_res)
  if conformer < 0:
  # this is NOT a conformer job
  self.running_jobs[label].append(job.job_name) # mark as a running job
@@ -1001,21 +1002,20 @@ def check_scan_job(self, label, job):
  for i in range(self.species_dict[label].number_of_rotors):
  message = ''
  invalidation_reason = ''
+ trsh = False
  if self.species_dict[label].rotors_dict[i]['pivots'] == job.pivots:
  invalidate = False
  if job.job_status[1] == 'done':
- # ESS converged. Get PES using Arkane:
+ # ESS converged. Get PES scan using Arkane:
  log = Log(path='')
  log.determine_qm_software(fullpath=job.local_path_to_output_file)
- plot_scan = True
  try:
  v_list, angle = log.software_log.loadScanEnergies()
  except ZeroDivisionError:
  logging.error('Energies from rotor scan of {label} between pivots {pivots} could not'
  'be read. Invalidating rotor.'.format(label=label, pivots=job.pivots))
  invalidate = True
  invalidation_reason = 'could not read energies'
- plot_scan = False
  else:
  v_list = np.array(v_list, np.float64)
  v_list = v_list * 0.001 # convert to kJ/mol
@@ -1033,6 +1033,10 @@ def check_scan_job(self, label, job):
  invalidate = True
  invalidation_reason = 'initial and final points are inconsistent by more than {0}' \
  ' kJ/mol'.format(inconsistency_az)
+ if not job.scan_trsh:
+ logging.info('Trying to troubleshoot rotor {0} of {1}...'.format(job.pivots, label))
+ trsh = True
+ self.troubleshoot_scan_job(job=job)
  if not invalidate:
  v_last = v_list[-1]
  for v in v_list:
@@ -1048,6 +1052,11 @@ def check_scan_job(self, label, job):
  invalidate = True
  invalidation_reason = 'two consecutive points are inconsistent by more than {0}' \
  ' kJ/mol'.format(inconsistency_ab)
+ if not job.scan_trsh:
+ logging.info('Trying to troubleshoot rotor {0} of {1}...'.format(
+ job.pivots, label))
+ trsh = True
+ self.troubleshoot_scan_job(job=job)
  break
  if abs(v - v_list[0]) > maximum_barrier:
  # The barrier for the hinderd rotor is higher than `maximum_barrier` kJ/mol.
@@ -1064,7 +1073,7 @@ def check_scan_job(self, label, job):
  v_last = v
  # 2. Check conformation:
  invalidated = ''
- if not invalidate:
+ if not invalidate and not trsh:
  v_diff = (v_list[0] - np.min(v_list))
  if v_diff >= 2 or v_diff > 0.5 * (max(v_list) - min(v_list)):
  self.species_dict[label].rotors_dict[i]['success'] = False
@@ -1085,15 +1094,14 @@ def check_scan_job(self, label, job):
  self.run_opt_job(label) # run opt on new initial_xyz with the desired dihedral
  else:
  self.species_dict[label].rotors_dict[i]['success'] = True
- else:
+ elif invalidate:
  invalidated = '*INVALIDATED* '
- symmetry = ''
  if self.species_dict[label].rotors_dict[i]['success']:
  self.species_dict[label].rotors_dict[i]['symmetry'], _ = determine_rotor_symmetry(
  rotor_path=job.local_path_to_output_file, label=label,
  pivots=self.species_dict[label].rotors_dict[i]['pivots'])
  symmetry = ' has symmetry {0}'.format(self.species_dict[label].rotors_dict[i]['symmetry'])
- if plot_scan:
+
  logging.info('{invalidated}Rotor scan {scan} between pivots {pivots}'
  ' for {label}{symmetry}'.format(invalidated=invalidated,
  scan=self.species_dict[label].rotors_dict[i]['scan'],
@@ -1107,13 +1115,15 @@ def check_scan_job(self, label, job):
  else:
  # scan job crashed
  invalidate = True
- if invalidate:
- self.species_dict[label].rotors_dict[i]['success'] = False
- else:
- self.species_dict[label].rotors_dict[i]['success'] = True
- self.species_dict[label].rotors_dict[i]['scan_path'] = job.local_path_to_output_file
- self.species_dict[label].rotors_dict[i]['invalidation_reason'] = invalidation_reason
- break # `job` has only one pivot. Break if found, otherwise raise an error.
+ invalidation_reason = 'scan job crashed'
+ if not trsh:
+ if invalidate:
+ self.species_dict[label].rotors_dict[i]['success'] = False
+ else:
+ self.species_dict[label].rotors_dict[i]['success'] = True
+ self.species_dict[label].rotors_dict[i]['scan_path'] = job.local_path_to_output_file
+ self.species_dict[label].rotors_dict[i]['invalidation_reason'] = invalidation_reason
+ break # A job object has only one pivot. Break if found, otherwise raise an error.
  else:
  raise SchedulerError('Could not match rotor with pivots {0} in species {1}'.format(job.pivots, label))
  self.save_restart_dict()
@@ -1201,9 +1211,30 @@ def troubleshoot_negative_freq(self, label, job):
  for i, xyz in enumerate(self.species_dict[label].conformers):
  self.run_job(label=label, xyz=xyz, level_of_theory=self.conformer_level, job_type='conformer', conformer=i)
 
+ def troubleshoot_scan_job(self, job):
+ """
+ Try freezing all dihedrals other than the scan's pivots for this job
+ """
+ label = job.species_name
+ species_scan_lists = [rotor_dict['scan'] for rotor_dict in self.species_dict[label].rotors_dict.values()]
+ if job.scan not in species_scan_lists:
+ raise SchedulerError('Could not find the dihedral to troubleshoot for in the dcan list of species'
+ ' {0}'.format(label))
+ species_scan_lists.pop(species_scan_lists.index(job.scan))
+ if len(species_scan_lists):
+ scan_trsh = '\n'
+ for scan in species_scan_lists:
+ scan_trsh += 'D ' + ''.join([str(num) + ' ' for num in scan]) + 'F\n'
+ scan_res = min(4, int(job.scan_res / 2))
+ # make sure mod(360, scan res) is 0:
+ if scan_res not in [4, 2, 1]:
+ scan_res = min([4, 2, 1], key=lambda x:abs(x - scan_res))
+ self.run_job(label=label, xyz=job.xyz, level_of_theory=job.level_of_theory, job_type='scan',
+ scan=job.scan, pivots=job.pivots, scan_trsh=scan_trsh, scan_res=4)
+
  def troubleshoot_opt_jobs(self, label):
  """
- we're troubleshooting for opt jobs.
+ We're troubleshooting for opt jobs.
  First check for server status and troubleshoot if needed. Then check for ESS status and troubleshoot
  if needed. Finally, check whether or not the last job had fine=True, add if it didn't run with fine.
  """