diff --git a/arc/job/job.py b/arc/job/job.py index f6ba97057f..1bb7458f6c 100644 --- a/arc/job/job.py +++ b/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""" diff --git a/arc/job/submit.py b/arc/job/submit.py index a05fc8cb9b..431c2080cf 100644 --- a/arc/job/submit.py +++ b/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 . diff --git a/arc/main.py b/arc/main.py index 85e831e9af..89873b2632 100644 --- a/arc/main.py +++ b/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): """ diff --git a/arc/mainTest.py b/arc/mainTest.py index fb7cbb53c7..ee39dd3f0f 100644 --- a/arc/mainTest.py +++ b/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}], } diff --git a/arc/scheduler.py b/arc/scheduler.py index 71b2606a5c..6930c5c296 100644 --- a/arc/scheduler.py +++ b/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) diff --git a/arc/schedulerTest.py b/arc/schedulerTest.py index 0807e80f95..008b6f1b41 100644 --- a/arc/schedulerTest.py +++ b/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""" diff --git a/arc/species/converter.py b/arc/species/converter.py index 022f5d7199..d25627508d 100644 --- a/arc/species/converter.py +++ b/arc/species/converter.py @@ -303,7 +303,7 @@ def order_atoms_in_mol_list(ref_mol, mol_list): def order_atoms(ref_mol, mol): """Order the atoms in `mol` by the atom order in ref_mol""" - if mol is not None: + if ref_mol is not None and mol is not None: ref_mol_is_iso_copy = ref_mol.copy(deep=True) mol_is_iso_copy = mol.copy(deep=True) ref_mol_find_iso_copy = ref_mol.copy(deep=True) diff --git a/arc/species/species.py b/arc/species/species.py index 144209ec8e..c9d909891e 100644 --- a/arc/species/species.py +++ b/arc/species/species.py @@ -50,7 +50,7 @@ class ARCSpecies(object): `number_of_rotors` ``int`` The number of potential rotors to scan `rotors_dict` ``dict`` A dictionary of rotors. structure given below. `conformers` ``list`` A list of selected conformers XYZs - `conformer_energies` ``list`` A list of conformers E0 (Hartree) + `conformer_energies` ``list`` A list of conformers E0 (in J/mol) `initial_xyz` ``string`` The initial geometry guess `final_xyz` ``string`` The optimized species geometry `opt_level` ``string`` Level of theory for geometry optimization. Saved for archiving. @@ -106,11 +106,13 @@ class ARCSpecies(object): def __init__(self, is_ts=False, rmg_species=None, mol=None, label=None, xyz=None, multiplicity=None, charge=None, smiles='', adjlist='', inchi='', bond_corrections=None, generate_thermo=True, species_dict=None, yml_path=None, ts_methods=None, ts_number=None, rxn_label=None, external_symmetry=None, - optical_isomers=None, run_time=None): + optical_isomers=None, run_time=None, conformers_path=None): self.t1 = None self.ts_number = ts_number self.conformers = list() self.conformer_energies = list() + if conformers_path is not None: + self.append_conformers(conformers_path) self.xyzs = list() # used for conformer search self.thermo = None self.rmg_thermo = None @@ -152,7 +154,7 @@ def __init__(self, is_ts=False, rmg_species=None, mol=None, label=None, xyz=None self.chosen_ts = None self.generate_thermo = generate_thermo if not self.is_ts else False self.long_thermo_description = '' - self.opt_level = '' + self.opt_level = None self.ts_report = '' self.yml_path = yml_path self.final_xyz = '' @@ -240,7 +242,7 @@ def __init__(self, is_ts=False, rmg_species=None, mol=None, label=None, xyz=None # consider the initial guess as one of the conformers if generating others. # otherwise, just consider it as the conformer. self.conformers.append(self.initial_xyz) - self.conformer_energies.append(0.0) # dummy + self.conformer_energies.append(None) # dummy self.neg_freqs_trshed = list() @@ -259,7 +261,7 @@ def __init__(self, is_ts=False, rmg_species=None, mol=None, label=None, xyz=None if not isinstance(self.charge, int): raise SpeciesError('Charge for species {0} is not an integer (got {1}, a {2})'.format( self.label, self.charge, type(self.charge))) - if not self.is_ts and self.initial_xyz is None and self.mol is None: + if not self.is_ts and self.initial_xyz is None and self.mol is None and not self.conformers: raise SpeciesError('No structure (xyz, SMILES, adjList, RMG:Species, or RMG:Molecule) was given for' ' species {0}'.format(self.label)) if self.label is None: @@ -296,7 +298,8 @@ def as_dict(self): species_dict['multiplicity'] = self.multiplicity species_dict['charge'] = self.charge species_dict['generate_thermo'] = self.generate_thermo - species_dict['opt_level'] = self.opt_level + if self.opt_level is not None: + species_dict['opt_level'] = self.opt_level species_dict['final_xyz'] = self.final_xyz species_dict['number_of_rotors'] = self.number_of_rotors species_dict['rotors_dict'] = self.rotors_dict @@ -371,7 +374,7 @@ def from_dict(self, species_dict): self.generate_thermo = False else: self.generate_thermo = species_dict['generate_thermo'] if 'generate_thermo' in species_dict else True - self.opt_level = species_dict['opt_level'] if 'opt_level' in species_dict else '' + self.opt_level = species_dict['opt_level'] if 'opt_level' in species_dict else None self.number_of_rotors = species_dict['number_of_rotors'] if 'number_of_rotors' in species_dict else 0 self.rotors_dict = species_dict['rotors_dict'] if 'rotors_dict' in species_dict else dict() self.external_symmetry = species_dict['external_symmetry'] if 'external_symmetry' in species_dict else None @@ -407,13 +410,15 @@ def from_dict(self, species_dict): if not self.charge: self.mol_list = self.mol.generate_resonance_structures(keep_isomorphic=False, filter_structures=True) - if self.mol is None and self.initial_xyz is None and not self.final_xyz: + if 'conformers_path' in species_dict: + self.append_conformers(species_dict['conformers_path']) + if self.mol is None and self.initial_xyz is None and not self.final_xyz and not self.conformers: raise SpeciesError('Must have either mol or xyz for species {0}'.format(self.label)) if self.initial_xyz is not None and not self.final_xyz: # consider the initial guess as one of the conformers if generating others. # otherwise, just consider it as the conformer. self.conformers.append(self.initial_xyz) - self.conformer_energies.append(0.0) # dummy + self.conformer_energies.append(None) # dummy def from_yml_file(self, label=None): """ @@ -470,7 +475,7 @@ def generate_conformers(self): self.find_conformers(self.mol) for xyz in self.xyzs: self.conformers.append(xyz) - self.conformer_energies.append(0.0) # a placeholder (lists are synced) + self.conformer_energies.append(None) # a placeholder (lists are synced) else: # generate "conformers" from the results of the different TS guess methods, if successful tsg_index = 0 @@ -479,7 +484,7 @@ def generate_conformers(self): for tsg in self.ts_guesses: if tsg.success: self.conformers.append(tsg.xyz) - self.conformer_energies.append(0.0) # a placeholder (lists are synced) + self.conformer_energies.append(None) # a placeholder (lists are synced) tsg.index = tsg_index tsg_index += 1 self.successful_methods.append(tsg.method) @@ -845,6 +850,30 @@ def mol_from_xyz(self, xyz=None): filter_structures=True) order_atoms_in_mol_list(ref_mol=self.mol, mol_list=self.mol_list) + def append_conformers(self, conformers_path): + """ + Populate the conformers and conformer energies lists with data from an ARC's conformers file. + The `conformers_path` should direct to either a "conformers_before_optimization" or + a "conformers_after_optimization" ARC file. + """ + if not os.path.isfile(conformers_path): + raise ValueError('Conformers file {0} could not be opened'.format(conformers_path)) + with open(conformers_path, 'r') as f: + lines = f.readlines() + conformer = '' + for line in lines: + if 'conformer' in line or 'SMILES' in line or 'Failed to converge' in line or line in ['\r', '\n', '\r\n']: + continue_reading_conformer = False + elif 'Relative Energy' in line: + self.conformer_energies.append(float(line.split()[2]) * 1000) # convert kJ/mol to J/mol + continue_reading_conformer = False + else: + conformer += line + continue_reading_conformer = True + if not continue_reading_conformer and conformer: + self.conformers.append(conformer) + conformer = '' + class TSGuess(object): """ diff --git a/arc/species/speciesTest.py b/arc/species/speciesTest.py index 4da7b8fced..49a6469604 100644 --- a/arc/species/speciesTest.py +++ b/arc/species/speciesTest.py @@ -8,6 +8,7 @@ from __future__ import (absolute_import, division, print_function, unicode_literals) import unittest import os +import shutil from rmgpy.molecule.molecule import Molecule from rmgpy.species import Species @@ -16,7 +17,9 @@ from arc.species.species import ARCSpecies, TSGuess, get_min_energy_conformer,\ determine_rotor_type, determine_rotor_symmetry, check_species_xyz from arc.species.converter import get_xyz_string, get_xyz_matrix, molecules_from_xyz -from arc.settings import arc_path +from arc.settings import arc_path, default_levels_of_theory +from arc.rmgdb import make_rmg_database_object +from arc.scheduler import Scheduler ################################################################################ @@ -364,7 +367,6 @@ def test_as_dict(self): 'charge': 0, 'is_ts': False, 'final_xyz': '', - 'opt_level': '', 't1': None, 'bond_corrections': {'C-H': 3, 'C-N': 1, 'H-N': 2}, 'rotors_dict': {}} @@ -512,6 +514,77 @@ def test_preserving_multiplicity(self): self.assertEqual(spc.mol.multiplicity, multiplicity_list[i]) self.assertTrue(all([structure.multiplicity == spc.multiplicity for structure in spc.mol_list])) + def test_append_conformers(self): + """Test that ARC correctly parses its own conformer files""" + settings = {'gaussian': 'server1', 'molpro': 'server2', 'qchem': 'server1', 'ssh': False} + project_directory = os.path.join(arc_path, 'Projects', 'arc_project_for_testing_delete_after_usage4') + spc1 = ARCSpecies(label=str('vinoxy'), smiles=str('C=C[O]')) + rmgdb = make_rmg_database_object() + sched1 = Scheduler(project='project_test', settings=settings, species_list=[spc1], + composite_method='', conformer_level=default_levels_of_theory['conformer'], + opt_level=default_levels_of_theory['opt'], freq_level=default_levels_of_theory['freq'], + sp_level=default_levels_of_theory['sp'], scan_level=default_levels_of_theory['scan'], + ts_guess_level=default_levels_of_theory['ts_guesses'], rmgdatabase=rmgdb, + project_directory=project_directory, generate_conformers=True, testing=True, + orbitals_level=default_levels_of_theory['orbitals']) + xyzs = ["""O 1.09068700 0.26516800 -0.16706300 +C 2.92204100 -1.18335700 -0.38884900 +C 2.27655500 -0.00373900 0.08543500 +H 2.36544800 -1.88781000 -0.99914600 +H 3.96112000 -1.38854500 -0.14958800 +H 2.87813500 0.68828400 0.70399400 +""", + """O 1.19396100 -0.06003700 0.03890100 +C 3.18797000 0.77061300 -0.87352700 +C 2.43591200 -0.04439300 0.02171600 +H 4.27370000 0.76090200 -0.86286100 +H 2.66641700 1.41155700 -1.57757300 +H 3.00398000 -0.68336800 0.72359800 +""", + """O 1.35241100 -1.02956000 -0.24056200 +C -0.72084300 0.01308200 0.09573000 +C 0.69217700 0.01185100 -0.09044300 +H -1.25803800 -0.93018100 0.10926800 +H -1.26861200 0.94177100 0.22420100 +H 1.20290400 0.99303700 -0.09819400 +""", + """O -1.40102900 -0.98575100 -0.11588500 +C 0.72457000 -0.01076700 0.06448800 +C -0.69494600 0.03450000 -0.06206300 +H 1.22539000 -0.97248000 0.11741200 +H 1.31277400 0.90087100 0.10878400 +H -1.16675800 1.03362600 -0.11273700"""] + energies = [0, 5, 5, 5] # J/mol + + sched1.save_conformers_file(label='vinoxy', xyzs=xyzs) + self.assertTrue(os.path.isfile(os.path.join(project_directory, 'output', 'Species', 'vinoxy', 'geometry', + 'conformers_before_optimization.txt'))) + + sched1.save_conformers_file(label='vinoxy', xyzs=xyzs, energies=energies) + self.assertTrue(os.path.isfile(os.path.join(project_directory, 'output', 'Species', 'vinoxy', 'geometry', + 'conformers_after_optimization.txt'))) + + spc2 = ARCSpecies(label=str('vinoxy'), smiles=str('C=C[O]'), conformers_path=os.path.join(project_directory, + 'output', 'Species', 'vinoxy', 'geometry', 'conformers_before_optimization.txt')) + + spc3 = ARCSpecies(label=str('vinoxy'), smiles=str('C=C[O]'), conformers_path=os.path.join(project_directory, + 'output', 'Species', 'vinoxy', 'geometry', 'conformers_after_optimization.txt')) + + self.assertEqual(spc2.conformers[2], xyzs[2]) + self.assertEqual(spc3.conformers[2], xyzs[2]) + self.assertEqual(spc3.conformer_energies[2], energies[2]) + + @classmethod + def tearDownClass(cls): + """ + A function that is run ONCE after all unit tests in this class. + Delete all project directories created during these unit tests + """ + projects = ['arc_project_for_testing_delete_after_usage4'] + for project in projects: + project_directory = os.path.join(arc_path, 'Projects', project) + shutil.rmtree(project_directory) + class TestTSGuess(unittest.TestCase): """