Allow ARCSpecies to be initialized with a conformer list file

arc/job/job.py

@@ -234,7 +234,10 @@ def __init__(self, project, settings, species_name, xyz, job_type, level_of_theo
         self.submit = ''
         self.input = ''
         self.server_nodes = list()
-        self._write_initiated_job_to_csv_file()
+        if job_num is None:
+            # this checks jon_num and not self.job_num on purpose
+            # if job_num was given, then don't save as initiated jobs, this is a restarted job
+            self._write_initiated_job_to_csv_file()
 
     def as_dict(self):
         """A helper function for dumping this object as a dictionary in a YAML file for restarting ARC"""

arc/job/submit.py

@@ -31,12 +31,12 @@
 SubmitDir=`pwd`
 
 GAUSS_SCRDIR=/scratch/users/{un}/g09/$SLURM_JOB_NAME-$SLURM_JOB_ID
-export  GAUSS_SCRDIR
+export GAUSS_SCRDIR
 
 mkdir -p $GAUSS_SCRDIR
 mkdir -p $WorkDir
 
-cd  $WorkDir
+cd $WorkDir
 . $g09root/g09/bsd/g09.profile
 
 cp $SubmitDir/input.gjf .
@@ -79,7 +79,7 @@
 SubmitDir=`pwd`
 
 mkdir -p $WorkDir
-cd  $WorkDir
+cd $WorkDir
 
 cp $SubmitDir/input.inp .
 

arc/main.py

@@ -3,6 +3,7 @@
 
 from __future__ import (absolute_import, division, print_function, unicode_literals)
 import logging
+import warnings
 import sys
 import os
 import time
@@ -776,6 +777,8 @@ def log_header(self, level=logging.INFO):
             logging.log(level, 'The current git HEAD for ARC is:')
             logging.log(level, '    {0}\n    {1}\n'.format(head, date))
         logging.info('Starting project {0}'.format(self.project))
+        # ignore Paramiko warnings:
+        warnings.filterwarnings(action='ignore', module='.*paramiko.*')
 
     def log_footer(self, level=logging.INFO):
         """

arc/mainTest.py

@@ -75,7 +75,6 @@ def test_as_dict(self):
                                       'multiplicity': 1,
                                       'neg_freqs_trshed': [],
                                       'number_of_rotors': 0,
-                                      'opt_level': '',
                                       'rotors_dict': {},
                                       't1': None}],
                          }

arc/scheduler.py

@@ -308,7 +308,7 @@ def __init__(self, project, settings, species_list, composite_method, conformer_
                             # restart-related check are performed in run_scan_jobs()
                             self.run_scan_jobs(species.label)
                 elif not self.species_dict[species.label].is_ts and self.generate_conformers\
-                        and 'geo' not in self.output[species.label]:
+                        and 'geo' not in self.output[species.label] and not species.conformers:
                     self.species_dict[species.label].generate_conformers()
             else:
                 # Species is loaded from a YAML file
@@ -325,6 +325,13 @@ def schedule_jobs(self):
         """
         The main job scheduling block
         """
+        for species in self.species_dict.values():
+            if not species.initial_xyz and not species.final_xyz and species.conformers and species.conformer_energies:
+                self.determine_most_stable_conformer(species.label)
+                if self.composite_method:
+                    self.run_composite_job(species.label)
+                else:
+                    self.run_opt_job(species.label)
         if self.generate_conformers:
             self.run_conformer_jobs()
         while self.running_jobs != {}:  # loop while jobs are still running
@@ -442,7 +449,8 @@ def schedule_jobs(self):
                             folder_name = 'rxns' if self.species_dict[label].is_ts else 'Species'
                             orbitals_path = os.path.join(self.project_directory, 'output', folder_name, label, 'geometry',
                                                          'orbitals.fchk')
-                            shutil.copyfile(job.local_path_to_orbitals_file, orbitals_path)
+                            if os.path.isfile(job.local_path_to_orbitals_file):
+                                shutil.copyfile(job.local_path_to_orbitals_file, orbitals_path)
                         self.timer = False
                         break
 
@@ -549,15 +557,8 @@ def run_conformer_jobs(self):
         """
         for label in self.unique_species_labels:
             if not self.species_dict[label].is_ts and 'opt converged' not in self.output[label]['status']\
-                        and 'opt' not in self.job_dict[label]:
-                geo_dir = os.path.join(self.project_directory, 'output', 'Species', label, 'geometry')
-                if not os.path.exists(geo_dir):
-                    os.makedirs(geo_dir)
-                conf_path = os.path.join(geo_dir, 'conformers_before_optimization.txt')
-                with open(conf_path, 'w') as f:
-                    for conf in self.species_dict[label].conformers:
-                        f.write(conf)
-                        f.write('\n\n')
+                        and 'opt' not in self.job_dict[label] and not self.species_dict[label].conformer_energies:
+                self.save_conformers_file(label)
                 if not self.testing:
                     if len(self.species_dict[label].conformers) > 1:
                         self.job_dict[label]['conformers'] = dict()
@@ -714,13 +715,13 @@ def run_orbitals_job(self, label):
 
     def parse_conformer_energy(self, job, label, i):
         """
-        Parse E0 (Hartree) from the conformer opt output file, and save it in the 'conformer_energies' attribute.
+        Parse E0 (J/mol) from the conformer opt output file, and save it in the 'conformer_energies' attribute.
         """
         if job.job_status[1] == 'done':
             log = Log(path='')
             log.determine_qm_software(fullpath=job.local_path_to_output_file)
-            e0 = log.software_log.loadEnergy()
-            self.species_dict[label].conformer_energies[i] = e0  # in J/mol
+            self.species_dict[label].conformer_energies[i] = log.software_log.loadEnergy()  # in J/mol
+            logging.debug('energy for {0} is {1}'.format(i, self.species_dict[label].conformer_energies[i]))
         else:
             logging.warn('Conformer {i} for {label} did not converge!'.format(i=i, label=label))
 
@@ -733,7 +734,7 @@ def determine_most_stable_conformer(self, label):
         if self.species_dict[label].is_ts:
             raise SchedulerError('The determine_most_stable_conformer() method does not deal with transition'
                                  ' state guesses.')
-        if all(e == 0.0 for e in self.species_dict[label].conformer_energies):
+        if 'conformers' in self.job_dict[label] and all(e is None for e in self.species_dict[label].conformer_energies):
             logging.error('No conformer converged for species {0}! Trying to troubleshoot conformer jobs...'.format(
                 label))
             for i, job in self.job_dict[label]['conformers'].items():
@@ -743,34 +744,19 @@ def determine_most_stable_conformer(self, label):
             conformer_xyz = None
             xyzs = list()
             log = Log(path='')
-            for job in self.job_dict[label]['conformers'].values():
-                log.determine_qm_software(fullpath=job.local_path_to_output_file)
-                try:
-                    coord, number, _ = log.software_log.loadGeometry()
-                except RMGInputError:
-                    xyzs.append(None)
-                else:
-                    xyzs.append(get_xyz_string(xyz=coord, number=number))
-            energies, xyzs = (list(t) for t in zip(*sorted(zip(self.species_dict[label].conformer_energies, xyzs))))
-            smiles_list = list()
-            for xyz in xyzs:
-                b_mol = molecules_from_xyz(xyz, multiplicity=self.species_dict[label].multiplicity,
-                                           charge=self.species_dict[label].charge)[1]
-                smiles = b_mol.toSMILES() if b_mol is not None else 'no 2D structure'
-                smiles_list.append(smiles)
-            geo_dir = os.path.join(self.project_directory, 'output', 'Species', label, 'geometry')
-            if not os.path.exists(geo_dir):
-                os.makedirs(geo_dir)
-            conf_path = os.path.join(geo_dir, 'conformers_after_optimization.txt')
-            with open(conf_path, 'w') as f:
-                for i, xyz in enumerate(xyzs):
-                    f.write('conformer {0}:\n'.format(i))
-                    if xyz is not None:
-                        f.write(xyz + '\n')
-                        f.write('SMILES: ' + smiles_list[i] + '\n')
-                        f.write('Relative Energy: {0} kJ/mol\n\n\n'.format((energies[i] - min(energies)) * 0.001))
+            if self.species_dict[label].conformer_energies:
+                xyzs = self.species_dict[label].conformers
+            else:
+                for job in self.job_dict[label]['conformers'].values():
+                    log.determine_qm_software(fullpath=job.local_path_to_output_file)
+                    try:
+                        coord, number, _ = log.software_log.loadGeometry()
+                    except RMGInputError:
+                        xyzs.append(None)
                     else:
-                        f.write('Failed to converge')
+                        xyzs.append(get_xyz_string(xyz=coord, number=number))
+            energies, xyzs = (list(t) for t in zip(*sorted(zip(self.species_dict[label].conformer_energies, xyzs))))
+            self.save_conformers_file(label, xyzs=xyzs, energies=energies)
             # Run isomorphism checks if a 2D representation is available
             if self.species_dict[label].mol is not None:
                 for i, xyz in enumerate(xyzs):
@@ -781,11 +767,14 @@ def determine_most_stable_conformer(self, label):
                             is_isomorphic = check_isomorphism(self.species_dict[label].mol, b_mol)
                         except ValueError as e:
                             if self.species_dict[label].charge:
-                                logging.error('Could not determine isomorphism for charged species. Got the '
-                                              'following error:\n{0}'.format(e.message))
+                                logging.error('Could not determine isomorphism for charged species {0}. '
+                                              'Optimizing the most stable conformer anyway. Got the '
+                                              'following error:\n{1}'.format(label, e))
                             else:
-                                logging.error('Could not determine isomorphism for (non-charged) species. Got the '
-                                              'following error:\n{0}'.format(e.message))
+                                logging.error('Could not determine isomorphism for (non-charged) species {0}. '
+                                              'Optimizing the most stable conformer anyway. Got the '
+                                              'following error:\n{1}'.format(label, e))
+                            conformer_xyz = xyzs[0]
                             break
                         if is_isomorphic:
                             if i == 0:
@@ -818,7 +807,7 @@ def determine_most_stable_conformer(self, label):
                         smiles_list.append(molecules_from_xyz(xyz, multiplicity=self.species_dict[label].multiplicity,
                                                               charge=self.species_dict[label].charge)[1])
                     if self.allow_nonisomorphic_2d or self.species_dict[label].charge:
-                        # we'll optimize the most stable conformer even if it not isomorphic to the 2D graph
+                        # we'll optimize the most stable conformer even if it is not isomorphic to the 2D graph
                         logging.error('No conformer for {0} was found to be isomorphic with the 2D graph representation'
                                       ' {1} (got: {2}). Optimizing the most stable conformer anyway.'.format(
                                        label, self.species_dict[label].mol.toSMILES(), smiles_list))
@@ -846,7 +835,7 @@ def determine_most_likely_ts_conformer(self, label):
         if not self.species_dict[label].is_ts:
             raise SchedulerError('The determine_most_likely_ts_conformer() method only deals with transition'
                                  ' state guesses.')
-        if all(e == 0.0 for e in self.species_dict[label].conformer_energies):
+        if all(e is None for e in self.species_dict[label].conformer_energies):
             logging.error('No guess converged for TS {0}!')
             # for i, job in self.job_dict[label]['conformers'].items():
             #     self.troubleshoot_ess(label, job, level_of_theory=job.level_of_theory, job_type='conformer',
@@ -893,6 +882,7 @@ def parse_composite_geo(self, label, job):
             self.species_dict[label].final_xyz = get_xyz_string(xyz=coord, number=number)
             self.output[label]['status'] += 'composite converged; '
             self.output[label]['composite'] = os.path.join(job.local_path, 'output.out')
+            self.species_dict[label].opt_level = self.composite_method
             rxn_str = ''
             if self.species_dict[label].is_ts:
                 rxn_str = ' of reaction {0}'.format(self.species_dict[label].rxn_label)
@@ -914,7 +904,6 @@ def parse_composite_geo(self, label, job):
                 return True  # run freq / scan jobs on this optimized geometry
             elif not self.species_dict[label].is_ts:
                 self.troubleshoot_negative_freq(label=label, job=job)
-            self.species_dict[label].opt_level = self.composite_method
         if job.job_status[1] != 'done' or not freq_ok:
             self.troubleshoot_ess(label=label, job=job, level_of_theory=job.level_of_theory, job_type='composite')
         return False  # return ``False``, so no freq / scan jobs are initiated for this unoptimized geometry
@@ -1314,7 +1303,7 @@ def troubleshoot_negative_freq(self, label, job):
             xyz2 = atomcoords - factor * displacement
             self.species_dict[label].conformers.append(get_xyz_string(xyz=xyz1, number=atomnos))
             self.species_dict[label].conformers.append(get_xyz_string(xyz=xyz2, number=atomnos))
-            self.species_dict[label].conformer_energies.extend([0.0, 0.0])  # a placeholder (lists are synced)
+            self.species_dict[label].conformer_energies.extend([None, None])  # a placeholder (lists are synced)
         self.job_dict[label]['conformers'] = dict()  # initialize the conformer job dictionary
         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)
@@ -1736,16 +1725,53 @@ def make_reaction_labels_info_file(self):
             f.write('Reaction labels and respective TS labels:\n\n')
         return rxn_info_path
 
+    def save_conformers_file(self, label, xyzs=None, energies=None):
+        """
+        A helper function for saving the conformers for species `label` before and after optimization
+        If `xyzs` is given, then it is used as the conformers xyz list, otherwise it is taken from the species.conformer
+        attribute. If `energies`, a list of respective conformer energies in J/mol, is given, it will also be reported.
+        """
+        geo_dir = os.path.join(self.project_directory, 'output', 'Species', label, 'geometry')
+        if not os.path.exists(geo_dir):
+            os.makedirs(geo_dir)
+        smiles_list = list()
+        xyzs = xyzs or self.species_dict[label].conformers
+        for xyz in xyzs:
+            b_mol = molecules_from_xyz(xyz, multiplicity=self.species_dict[label].multiplicity,
+                                       charge=self.species_dict[label].charge)[1]
+            smiles = b_mol.toSMILES() if b_mol is not None else 'no 2D structure'
+            smiles_list.append(smiles)
+        if energies is not None:
+            conf_path = os.path.join(geo_dir, 'conformers_after_optimization.txt')
+        else:
+            conf_path = os.path.join(geo_dir, 'conformers_before_optimization.txt')
+        with open(conf_path, 'w') as f:
+            if energies is not None:
+                f.write('conformers optimized at {0}\n\n'.format(self.conformer_level))
+            for i, conf in enumerate(xyzs):
+                f.write('conformer {0}:\n'.format(i))
+                if conf is not None:
+                    f.write(conf + '\n')
+                    f.write('SMILES: ' + smiles_list[i] + '\n')
+                    if energies is not None:
+                        if energies[i] == min(energies):
+                            f.write('Relative Energy: 0 kJ/mol (lowest)')
+                        else:
+                            f.write('Relative Energy: {0:.3f} kJ/mol'.format((energies[i] - min(energies)) * 0.001))
+                else:
+                    f.write('Failed to converge')
+                f.write('\n\n\n')
+
 
-# Add a custom string representer to use block literals for multiline strings
 def string_representer(dumper, data):
+    """Add a custom string representer to use block literals for multiline strings"""
     if len(data.splitlines()) > 1:
         return dumper.represent_scalar(tag='tag:yaml.org,2002:str', value=data, style='|')
     return dumper.represent_scalar(tag='tag:yaml.org,2002:str', value=data)
 
 
-# Add a custom unicode representer to use block literals for multiline strings
 def unicode_representer(dumper, data):
+    """Add a custom unicode representer to use block literals for multiline strings"""
     if len(data.splitlines()) > 1:
         return yaml.ScalarNode(tag='tag:yaml.org,2002:str', value=data, style='|')
     return yaml.ScalarNode(tag='tag:yaml.org,2002:str', value=data)

arc/schedulerTest.py

@@ -63,10 +63,14 @@ def test_conformers(self):
         self.sched1.job_dict[label]['conformers'] = dict()
         self.sched1.job_dict[label]['conformers'][0] = self.job1
         self.sched1.job_dict[label]['conformers'][1] = self.job2
+        self.sched1.species_dict[label].conformer_energies = [None, None]
+        self.sched1.species_dict[label].conformers = [None, None]
         self.sched1.parse_conformer_energy(job=self.job1, label=label, i=0)
         self.sched1.parse_conformer_energy(job=self.job2, label=label, i=1)
         expecting = [-251596443.5088726, -254221943.3698632]
         self.assertEqual(self.sched1.species_dict[label].conformer_energies, expecting)
+        self.sched1.species_dict[label].conformers[0] = parser.parse_xyz_from_file(self.job1.local_path_to_output_file)
+        self.sched1.species_dict[label].conformers[1] = parser.parse_xyz_from_file(self.job2.local_path_to_output_file)
 
         self.sched1.determine_most_stable_conformer(label=label)
         expecting = """N      -0.75555952   -0.12937106    0.00000000
@@ -83,20 +87,21 @@ def test_conformers(self):
         self.assertTrue(os.path.isfile(methylamine_conf_path))
         with open(methylamine_conf_path, 'r') as f:
             lines = f.readlines()
-        self.assertEqual(lines[0], 'conformer 0:\n')
-        self.assertEqual(lines[9], 'SMILES: CN\n')
-        self.assertTrue('Relative Energy:' in lines[10])
-        self.assertEqual(lines[14][0], 'N')
+        self.assertTrue('conformers optimized at' in lines[0])
+        self.assertEqual(lines[11], 'SMILES: CN\n')
+        self.assertTrue('Relative Energy:' in lines[12])
+        self.assertEqual(lines[16][0], 'N')
 
         self.sched1.run_conformer_jobs()
+        self.sched1.save_conformers_file(label='C2H6')
         c2h6_conf_path = os.path.join(self.sched1.project_directory, 'output', 'Species', 'C2H6', 'geometry',
                                       'conformers_before_optimization.txt')
         self.assertTrue(os.path.isfile(c2h6_conf_path))
         with open(c2h6_conf_path, 'r') as f:
             lines = f.readlines()
-        self.assertEqual(lines[0][0], 'C')
+        self.assertEqual(lines[1][0], 'C')
         self.assertEqual(lines[9], '\n')
-        self.assertEqual(lines[12][0], 'H')
+        self.assertEqual(lines[17][0], 'H')
 
     def test_check_negative_freq(self):
         """Test the check_negative_freq() method"""