diff --git a/arc/job/inputs.py b/arc/job/inputs.py index a529ecdc4d..94b5172c75 100644 --- a/arc/job/inputs.py +++ b/arc/job/inputs.py @@ -116,15 +116,13 @@ 'arkane_input_species': """#!/usr/bin/env python # -*- coding: utf-8 -*- -linear = {linear}{bonds} - -externalSymmetry = {symmetry} +{bonds}externalSymmetry = {symmetry} spinMultiplicity = {multiplicity} opticalIsomers = {optical} -energy = {{'{model_chemistry}': Log('{sp_path}')}} +energy = {{'{sp_level}': Log('{sp_path}')}} geometry = Log('{opt_path}') diff --git a/arc/main.py b/arc/main.py index a65199c2d1..498ecef516 100644 --- a/arc/main.py +++ b/arc/main.py @@ -67,8 +67,8 @@ class ARC(object): `output` ``dict`` Output dictionary with status and final QM file paths for all species Only used for restarting, the actual object used is in the Scheduler class `use_bac` ``bool`` Whether or not to use bond additivity corrections for thermo calculations - `model_chemistry` ``list`` The model chemistry in Arkane for energy corrections (AE, BAC). - This can be usually determined automatically. + `model_chemistry` ``str`` The model chemistry in Arkane for energy corrections (AE, BAC) + and frequencies/ZPE scaling factor. Can usually be determined automatically. `ess_settings` ``dict`` A dictionary of available ESS (keys) and a corresponding server list (values) `initial_trsh` ``dict`` Troubleshooting methods to try by default. Keys are ESS software, values are trshs 't0' ``float`` Initial time when the project was spawned @@ -772,66 +772,61 @@ def log_footer(self, level=logging.INFO): logging.log(level, 'ARC execution terminated on {0}'.format(time.asctime())) def determine_model_chemistry(self): - """Determine the model_chemistry used in Arkane""" + """Determine the model_chemistry to be used in Arkane + + Todo: + * Determine whether the model chemistry exists in Arkane automaticaly instead of hard coding + """ if self.model_chemistry: self.model_chemistry = self.model_chemistry.lower() - if self.model_chemistry not in ['cbs-qb3', 'cbs-qb3-paraskevas', 'ccsd(t)-f12/cc-pvdz-f12', - 'ccsd(t)-f12/cc-pvtz-f12', 'ccsd(t)-f12/cc-pvqz-f12', - 'b3lyp/cbsb7', 'b3lyp/6-311g(2d,d,p)', 'b3lyp/6-311+g(3df,2p)', - 'b3lyp/6-31g**']: + if self.model_chemistry.split('//')[0] not in [ + 'cbs-qb3', 'cbs-qb3-paraskevas', 'ccsd(t)-f12/cc-pvdz-f12', 'ccsd(t)-f12/cc-pvtz-f12', + 'ccsd(t)-f12/cc-pvqz-f12', 'b3lyp/cbsb7', 'b3lyp/6-311g(2d,d,p)', 'b3lyp/6-311+g(3df,2p)', + 'b3lyp/6-31g**']: logging.warn('No bond additivity corrections (BAC) are available in Arkane for "model chemistry"' ' {0}. As a result, thermodynamic parameters are expected to be inaccurate. Make sure that' ' atom energy corrections (AEC) were supplied or are available in Arkane to avoid' ' error.'.format(self.model_chemistry)) else: - # model chemistry was not given, try to determine it from the sp_level + # model chemistry was not given, try to determine it from the sp_level and freq_level model_chemistry = '' - if not self.composite_method: - sp_level = self.sp_level.lower() - else: - sp_level = self.composite_method - sp_level = sp_level.replace('f12a', 'f12').replace('f12b', 'f12') - if sp_level in ['ccsd(t)-f12/cc-pvdz', 'ccsd(t)-f12/cc-pvtz', 'ccsd(t)-f12/cc-pvqz']: - logging.warning('Using model chemistry {0} based on sp level {1}.'.format( - sp_level + '-f12', sp_level)) - model_chemistry = sp_level + '-f12' - elif not model_chemistry and sp_level in ['cbs-qb3', 'cbs-qb3-paraskevas', 'ccsd(t)-f12/cc-pvdz-f12', - 'ccsd(t)-f12/cc-pvtz-f12', 'ccsd(t)-f12/cc-pvqz-f12', - 'b3lyp/cbsb7', 'b3lyp/6-311g(2d,d,p)', 'b3lyp/6-311+g(3df,2p)', - 'b3lyp/6-31g**']: - model_chemistry = sp_level - elif self.use_bac: - logging.info('\n\n') - logging.warning('Could not determine appropriate Model Chemistry to be used in Arkane for' - ' thermochemical parameter calculations. Not using atom energy corrections and ' - 'bond additivity corrections!\n\n') - self.use_bac = False + if self.composite_method: + self.model_chemistry = self.composite_method.lower() else: - # use_bac is False, and no model chemistry was specified - if sp_level in ['m06-2x/cc-pvtz', 'g3', 'm08so/mg3s*', 'klip_1', 'klip_2', 'klip_3', 'klip_2_cc', - 'ccsd(t)-f12/cc-pvdz-f12_h-tz', 'ccsd(t)-f12/cc-pvdz-f12_h-qz', - 'ccsd(t)-f12/cc-pvdz-f12', 'ccsd(t)-f12/cc-pvtz-f12', 'ccsd(t)-f12/cc-pvqz-f12', - 'ccsd(t)-f12/cc-pcvdz-f12', 'ccsd(t)-f12/cc-pcvtz-f12', 'ccsd(t)-f12/cc-pcvqz-f12', - 'ccsd(t)-f12/cc-pvtz-f12(-pp)', 'ccsd(t)/aug-cc-pvtz(-pp)', 'ccsd(t)-f12/aug-cc-pvdz', - 'ccsd(t)-f12/aug-cc-pvtz', 'ccsd(t)-f12/aug-cc-pvqz', 'b-ccsd(t)-f12/cc-pvdz-f12', - 'b-ccsd(t)-f12/cc-pvtz-f12', 'b-ccsd(t)-f12/cc-pvqz-f12', 'b-ccsd(t)-f12/cc-pcvdz-f12', - 'b-ccsd(t)-f12/cc-pcvtz-f12', 'b-ccsd(t)-f12/cc-pcvqz-f12', 'b-ccsd(t)-f12/aug-cc-pvdz', - 'b-ccsd(t)-f12/aug-cc-pvtz', 'b-ccsd(t)-f12/aug-cc-pvqz', 'mp2_rmp2_pvdz', - 'mp2_rmp2_pvtz', 'mp2_rmp2_pvqz', 'ccsd-f12/cc-pvdz-f12', - 'ccsd(t)-f12/cc-pvdz-f12_noscale', 'g03_pbepbe_6-311++g_d_p', 'fci/cc-pvdz', - 'fci/cc-pvtz', 'fci/cc-pvqz', 'bmk/cbsb7', 'bmk/6-311g(2d,d,p)', 'b3lyp/6-31g**', - 'b3lyp/6-311+g(3df,2p)', 'MRCI+Davidson/aug-cc-pV(T+d)Z']: - model_chemistry = sp_level - self.model_chemistry = model_chemistry - logging.debug('Using {0} as model chemistry for energy corrections in Arkane'.format( - self.model_chemistry)) - if not self.model_chemistry: - logging.warn('Could not determine a Model Chemistry to be used in Arkane, NOT calculating thermodata') - for spc in self.arc_species_list: - spc.generate_thermo = False + sp_level = self.sp_level.replace('f12a', 'f12').replace('f12b', 'f12').lower() + freq_level = self.freq_level.replace('f12a', 'f12').replace('f12b', 'f12').lower() + if sp_level in ['ccsd(t)-f12/cc-pvdz', 'ccsd(t)-f12/cc-pvtz', 'ccsd(t)-f12/cc-pvqz']: + logging.warning('Using model chemistry {0} based on sp level {1}.'.format( + sp_level + '-f12', sp_level)) + sp_level += '-f12' + if sp_level not in ['ccsd(t)-f12/cc-pvdz-f12', 'ccsd(t)-f12/cc-pvtz-f12', 'ccsd(t)-f12/cc-pvqz-f12', + 'b3lyp/cbsb7', 'b3lyp/6-311g(2d,d,p)', 'b3lyp/6-311+g(3df,2p)', 'b3lyp/6-31g**']\ + and self.use_bac: + logging.info('\n\n') + logging.warning('Could not determine appropriate Model Chemistry to be used in Arkane for ' + 'thermochemical parameter calculations. Not using atom energy corrections and ' + 'bond additivity corrections!\n\n') + self.use_bac = False + elif sp_level not in ['m06-2x/cc-pvtz', 'g3', 'm08so/mg3s*', 'klip_1', 'klip_2', 'klip_3', 'klip_2_cc', + 'ccsd(t)-f12/cc-pvdz-f12_h-tz', 'ccsd(t)-f12/cc-pvdz-f12_h-qz', + 'ccsd(t)-f12/cc-pvdz-f12', 'ccsd(t)-f12/cc-pvtz-f12', 'ccsd(t)-f12/cc-pvqz-f12', + 'ccsd(t)-f12/cc-pcvdz-f12', 'ccsd(t)-f12/cc-pcvtz-f12', 'ccsd(t)-f12/cc-pcvqz-f12', + 'ccsd(t)-f12/cc-pvtz-f12(-pp)', 'ccsd(t)/aug-cc-pvtz(-pp)', 'ccsd(t)-f12/aug-cc-pvdz', + 'ccsd(t)-f12/aug-cc-pvtz', 'ccsd(t)-f12/aug-cc-pvqz', 'b-ccsd(t)-f12/cc-pvdz-f12', + 'b-ccsd(t)-f12/cc-pvtz-f12', 'b-ccsd(t)-f12/cc-pvqz-f12', 'b-ccsd(t)-f12/cc-pcvdz-f12', + 'b-ccsd(t)-f12/cc-pcvtz-f12', 'b-ccsd(t)-f12/cc-pcvqz-f12', 'b-ccsd(t)-f12/aug-cc-pvdz', + 'b-ccsd(t)-f12/aug-cc-pvtz', 'b-ccsd(t)-f12/aug-cc-pvqz', 'mp2_rmp2_pvdz', + 'mp2_rmp2_pvtz', 'mp2_rmp2_pvqz', 'ccsd-f12/cc-pvdz-f12', + 'ccsd(t)-f12/cc-pvdz-f12_noscale', 'g03_pbepbe_6-311++g_d_p', 'fci/cc-pvdz', + 'fci/cc-pvtz', 'fci/cc-pvqz', 'bmk/cbsb7', 'bmk/6-311g(2d,d,p)', 'b3lyp/6-31g**', + 'b3lyp/6-311+g(3df,2p)', 'MRCI+Davidson/aug-cc-pV(T+d)Z']: + logging.warn('Could not determine a Model Chemistry to be used in Arkane, ' + 'NOT calculating thermodata') + for spc in self.arc_species_list: + spc.generate_thermo = False + self.model_chemistry = sp_level + '//' + freq_level if self.model_chemistry: - logging.info('Using {0} as model chemistry for energy corrections in Arkane'.format( - self.model_chemistry)) + logging.info('Using {0} as a model chemistry in Arkane'.format(self.model_chemistry)) def determine_ess_settings(self, diagnostics=False): """ diff --git a/arc/mainTest.py b/arc/mainTest.py index 1b3ae12c59..b905f0fae2 100644 --- a/arc/mainTest.py +++ b/arc/mainTest.py @@ -53,13 +53,13 @@ def test_as_dict(self): arc_species_list=[spc1]) restart_dict = arc0.as_dict() expected_dict = {'composite_method': '', - 'conformer_level': 'b3lyp/6-31+g(d,p)', - 'ts_guess_level': 'b3lyp/6-31+g(d,p)', + 'conformer_level': 'b97d3/6-31+g(d,p)', + 'ts_guess_level': 'b97d3/6-31+g(d,p)', 'opt_level': 'wb97xd/6-311++g(d,p)', 'freq_level': 'wb97xd/6-311++g(d,p)', 'initial_trsh': 'scf=(NDump=30)', 'max_job_time': 120, - 'model_chemistry': 'ccsd(t)-f12/cc-pvtz-f12', + 'model_chemistry': 'ccsd(t)-f12/cc-pvtz-f12//wb97xd/6-311++g(d,p)', 'output': {}, 'project': 'arc_test', 'running_jobs': {}, @@ -266,9 +266,9 @@ def test_restart(self): spc2 = Species().fromSMILES(str('CC([O])=O')) spc2.generate_resonance_structures() spc2.thermo = db.thermo.getThermoData(spc2) - self.assertAlmostEqual(spc2.getEnthalpy(298), -178003.44650359568, 1) - self.assertAlmostEqual(spc2.getEntropy(298), 283.5983103176096, 1) - self.assertAlmostEqual(spc2.getHeatCapacity(1000), 118.99753808225603, 1) + self.assertAlmostEqual(spc2.getEnthalpy(298), -176074.01886272896, 1) + self.assertAlmostEqual(spc2.getEntropy(298), 283.2225158405262, 1) + self.assertAlmostEqual(spc2.getHeatCapacity(1000), 118.28356605714401, 1) self.assertTrue('arc_project_for_testing_delete_after_usage2' in spc2.thermo.comment) # delete the generated library from RMG-database @@ -332,6 +332,27 @@ def test_read_file(self): with self.assertRaises(InputError): read_file('nopath') + def test_determine_model_chemistry(self): + """Test determining the model chemistry""" + arc0 = ARC(project='arc_model_chemistry_test_0', level_of_theory='CBS-QB3') + # arc0.determine_model_chemistry() + self.assertEqual(arc0.model_chemistry, 'cbs-qb3') + + arc1 = ARC(project='arc_model_chemistry_test_1', model_chemistry='CBS-QB3-Paraskevas') + # arc1.determine_model_chemistry() + self.assertEqual(arc1.model_chemistry, 'cbs-qb3-paraskevas') + + arc2 = ARC(project='arc_model_chemistry_test_2', + level_of_theory='ccsd(t)-f12/cc-pvtz-f12//wb97xd/6-311++g(d,p)') + # arc2.determine_model_chemistry() + self.assertEqual(arc2.model_chemistry, 'ccsd(t)-f12/cc-pvtz-f12//wb97xd/6-311++g(d,p)') + + arc3 = ARC(project='arc_model_chemistry_test_2', + sp_level='ccsd(t)-f12/cc-pvtz-f12', opt_level='wb97xd/6-311++g(d,p)') + # arc3.determine_model_chemistry() + self.assertEqual(arc3.model_chemistry, 'ccsd(t)-f12/cc-pvtz-f12//wb97xd/6-311++g(d,p)') + + @classmethod def tearDownClass(cls): """ diff --git a/arc/parser.py b/arc/parser.py index 151bd6b899..c0b17934a6 100644 --- a/arc/parser.py +++ b/arc/parser.py @@ -77,7 +77,7 @@ def parse_t1(path): return t1 -def parse_e0(path): +def parse_e_elect(path, zpe_scale_factor=1.): """ Parse the zero K energy, E0, from an sp job output file """ @@ -85,11 +85,11 @@ def parse_e0(path): raise InputError('Could not find file {0}'.format(path)) log = determine_qm_software(fullpath=path) try: - e0 = log.loadEnergy(frequencyScaleFactor=1.) * 0.001 # convert to kJ/mol + e_elect = log.loadEnergy(zpe_scale_factor) * 0.001 # convert to kJ/mol except Exception: - logging.warning('Could not read E0 from {0}'.format(path)) - e0 = None - return e0 + logging.warning('Could not read e_elect from {0}'.format(path)) + e_elect = None + return e_elect def parse_xyz_from_file(path): diff --git a/arc/parserTest.py b/arc/parserTest.py index 5a85bef4d0..8efed8e13a 100644 --- a/arc/parserTest.py +++ b/arc/parserTest.py @@ -75,15 +75,15 @@ def test_parse_t1(self): t1 = parser.parse_t1(path) self.assertEqual(t1, 0.0086766) - def test_parse_e0(self): + def test_parse_e_elect(self): """Test parsing E0 from an sp job output file""" path1 = os.path.join(arc_path, 'arc', 'testing', 'mehylamine_CCSD(T).out') - e0 = parser.parse_e0(path1) - self.assertEqual(e0, -251377.49160993524) + e_elect = parser.parse_e_elect(path1) + self.assertEqual(e_elect, -251377.49160993524) path2 = os.path.join(arc_path, 'arc', 'testing', 'SO2OO_CBS-QB3.log') - e0 = parser.parse_e0(path2) - self.assertEqual(e0, -1833127.0939478774) + e_elect = parser.parse_e_elect(path2, zpe_scale_factor=0.99) + self.assertEqual(e_elect, -1833127.0939478774) def test_parse_dipole_moment(self): """Test parsing the dipole moment from an opt job output file""" diff --git a/arc/plotter.py b/arc/plotter.py index 4c97bda1e8..8ee23fedc8 100644 --- a/arc/plotter.py +++ b/arc/plotter.py @@ -40,7 +40,7 @@ def draw_3d(xyz=None, species=None, project_directory=None, save_only=False): """ Draws the molecule in a "3D-balls" style - If xyz ig given, it will be used, otherwise the function looks for species.final_xyz + If xyz is given, it will be used, otherwise the function looks for species.final_xyz Input coordinates are in string format Saves an image if a species and `project_directory` are provided If `save_only` is ``True``, then don't plot, only save the image diff --git a/arc/processor.py b/arc/processor.py index 626b0462c8..6f6857426c 100644 --- a/arc/processor.py +++ b/arc/processor.py @@ -12,7 +12,7 @@ from random import randint from arkane.input import species as arkane_input_species, transitionState as arkane_transition_state,\ - reaction as arkane_reaction + reaction as arkane_reaction, process_model_chemistry from arkane.statmech import StatMechJob, assign_frequency_scale_factor from arkane.thermo import ThermoJob from arkane.kinetics import KineticsJob @@ -40,8 +40,8 @@ class Processor(object): `rxn_list` ``list`` List of ARCReaction objects `output` ``dict`` Keys are labels, values are output file paths `use_bac` ``bool`` Whether or not to use bond additivity corrections for thermo calculations - `model_chemistry` ``list`` The model chemistry in Arkane for energy corrections (AE, BAC). - This can be usually determined automatically. + `sp_level` ``str`` The single point level of theory, used for atom and bond corrections in Arkane + `freq_level` ``str`` The frequency level of theory, used for the frequency scaling factor in Arkane `lib_long_desc` ``str`` A multiline description of levels of theory for the outputted RMG libraries `project_directory` ``str`` The path of the ARC project directory `t_min` ``tuple`` The minimum temperature for kinetics computations, e.g., (500, str('K')) @@ -58,7 +58,7 @@ def __init__(self, project, project_directory, species_dict, rxn_list, output, u self.rxn_list = rxn_list self.output = output self.use_bac = use_bac - self.model_chemistry = model_chemistry + self.sp_level, self.freq_level = process_model_chemistry(model_chemistry) self.lib_long_desc = lib_long_desc load_thermo_libs, load_kinetic_libs = False, False if any([species.is_ts and species.final_xyz for species in self.species_dict.values()])\ @@ -97,13 +97,6 @@ def _generate_arkane_species_file(self, species): species.arkane_file = species.yml_path return output_file_path - if not species.is_ts: - linear = species.is_linear() - if linear: - logging.info('Determined {0} to be a linear molecule'.format(species.label)) - species.long_thermo_description += 'Treated as a linear species\n' - else: - linear = False species.determine_symmetry() try: sp_path = self.output[species.label]['composite'] @@ -163,13 +156,13 @@ def _generate_arkane_species_file(self, species): input_file = input_files['arkane_input_species'] if self.use_bac and not species.is_ts: logging.info('Using the following BAC for {0}: {1}'.format(species.label, species.bond_corrections)) - bonds = '\n\nbonds = {0}'.format(species.bond_corrections) + bonds = 'bonds = {0}\n\n'.format(species.bond_corrections) else: logging.debug('NOT using BAC for {0}'.format(species.label)) bonds = '' - input_file = input_file.format(linear=linear, bonds=bonds, symmetry=species.external_symmetry, + input_file = input_file.format(bonds=bonds, symmetry=species.external_symmetry, multiplicity=species.multiplicity, optical=species.optical_isomers, - model_chemistry=self.model_chemistry, sp_path=sp_path, opt_path=opt_path, + sp_level=self.sp_level, sp_path=sp_path, opt_path=opt_path, freq_path=freq_path, rotors=rotors) with open(input_file_path, 'wb') as f: f.write(input_file) @@ -210,7 +203,7 @@ def process(self): species.thermo = arkane_spc.getThermoData() plotter.log_thermo(species.label, path=output_file_path[0]) species_for_thermo_lib.append(species) - if self.use_bac and self.model_chemistry: + if self.use_bac and self.sp_level: # If BAC was used, save another Arkane YAML file of this species with no BAC, so it can be used # for further rate calculations if needed (where the conformer.E0 has no BAC) statmech_success = self._run_statmech(arkane_spc, species.arkane_file, output_file_path[1], @@ -317,25 +310,29 @@ def process(self): f.write(str('\n')) def _run_statmech(self, arkane_spc, arkane_file, output_file_path=None, use_bac=False, kinetics=False, plot=False): - """ - A helper function for running an Arkane statmech job - `arkane_spc` is the species() function from Arkane's input.py - `arkane_file` is the Arkane species file (either .py or YAML form) - `output_file_path` is a path to the Arkane output.py file - `use_bac` is a bool flag indicating whether or not to use bond additivity corrections - `kinetics` is a bool flag indicating whether this specie sis part of a kinetics job, in which case..?? - `plot` is a bool flag indicating whether or not to plot a PDF of the calculated thermo properties + """A helper function for running an Arkane statmech job + + Args: + arkane_spc (str, unicode): An Arkane species() function representor. + arkane_file (str, unicode): The path to the Arkane species file (either in .py or YAML form). + output_file_path (str, unicode): The path to the Arkane output.py file. + use_bac (bool): A flag indicating whether or not to use bond additivity corrections (True to use). + kinetics (bool) A flag indicating whether this specie is part of a kinetics job. + plot (bool): A flag indicating whether to plot a PDF of the calculated thermo properties (True to plot) + + Returns: + bool: Whether the job was successful (True for successful). """ success = True stat_mech_job = StatMechJob(arkane_spc, arkane_file) - stat_mech_job.applyBondEnergyCorrections = use_bac and not kinetics and self.model_chemistry - if not kinetics or kinetics and self.model_chemistry: + stat_mech_job.applyBondEnergyCorrections = use_bac and not kinetics and self.sp_level + if not kinetics or (kinetics and self.sp_level): # currently we have to use a model chemistry for thermo - stat_mech_job.modelChemistry = self.model_chemistry + stat_mech_job.modelChemistry = self.sp_level else: - # if this is a klinetics computation and we don't have a valid model chemistry, don't bother about it + # if this is a kinetics computation and we don't have a valid model chemistry, don't bother about it stat_mech_job.applyAtomEnergyCorrections = False - stat_mech_job.frequencyScaleFactor = assign_frequency_scale_factor(self.model_chemistry) + stat_mech_job.frequencyScaleFactor = assign_frequency_scale_factor(self.freq_level) try: stat_mech_job.execute(outputFile=output_file_path, plot=plot) except Exception: diff --git a/arc/reaction.py b/arc/reaction.py index d9b404ca75..8ee4cbd2f3 100644 --- a/arc/reaction.py +++ b/arc/reaction.py @@ -310,24 +310,28 @@ def check_ts(self, log=True): Check that the TS E0 is above both reactants and products wells Return ``False`` if this test fails, else ``True`` """ - if any([spc.e0 is None for spc in self.r_species + self.p_species + [self.ts_species]]): + if any([spc.e_elect is None for spc in self.r_species + self.p_species + [self.ts_species]]): logging.error("Could not get E0's of all species participating in reaction {0}. Cannot check TS E0.".format( self.label)) - r_e0 = None if any([spc.e0 is None for spc in self.r_species]) else sum(spc.e0 for spc in self.r_species) - p_e0 = None if any([spc.e0 is None for spc in self.p_species]) else sum(spc.e0 for spc in self.p_species) - ts_e0 = self.ts_species.e0 - logging.error('Reactants E0: {0}\nProducts E0: {1}\nTS E0: {2}'.format(r_e0, p_e0, ts_e0)) + r_e_elect = None if any([spc.e_elect is None for spc in self.r_species])\ + else sum(spc.e_elect for spc in self.r_species) + p_e_elect = None if any([spc.e_elect is None for spc in self.p_species])\ + else sum(spc.e_elect for spc in self.p_species) + ts_e_elect = self.ts_species.e_elect + logging.error('Reactants E0: {0}\nProducts E0: {1}\nTS E0: {2}'.format(r_e_elect, p_e_elect, ts_e_elect)) return True - r_e0 = sum([spc.e0 for spc in self.r_species]) - p_e0 = sum([spc.e0 for spc in self.p_species]) - if self.ts_species.e0 < r_e0 or self.ts_species.e0 < p_e0: + r_e_elect = sum([spc.e_elect for spc in self.r_species]) + p_e_elect = sum([spc.e_elect for spc in self.p_species]) + if self.ts_species.e_elect < r_e_elect or self.ts_species.e_elect < p_e_elect: if log: logging.error('TS of reaction {0} has a lower E0 value than expected:\nReactants: {1} kJ/mol\nTS:' - ' {2} kJ/mol\nProducts: {3} kJ/mol'.format(self.label, r_e0, self.ts_species.e0, p_e0)) + ' {2} kJ/mol\nProducts: {3} kJ/mol'.format(self.label, r_e_elect, + self.ts_species.e_elect, p_e_elect)) return False if log: logging.info('Reaction {0} has the following path energies:\nReactants: {1} kJ/mol' - '\nTS: {2} kJ/mol\nProducts: {3} kJ/mol'.format(self.label, r_e0, self.ts_species.e0, p_e0)) + '\nTS: {2} kJ/mol\nProducts: {3} kJ/mol'.format(self.label, r_e_elect, + self.ts_species.e_elect, p_e_elect)) return True def check_attributes(self): diff --git a/arc/scheduler.py b/arc/scheduler.py index b5f58cac46..1ae777b040 100644 --- a/arc/scheduler.py +++ b/arc/scheduler.py @@ -1121,7 +1121,9 @@ def check_sp_job(self, label, job): self.output[label]['status'] += 'sp converged; ' self.output[label]['sp'] = os.path.join(job.local_path, 'output.out') self.species_dict[label].t1 = parser.parse_t1(self.output[label]['sp']) - self.species_dict[label].e0 = parser.parse_e0(self.output[label]['sp']) + zpe_scale_factor = 0.99 if self.composite_method.lower() == 'cbs-qb3' else 1.0 + self.species_dict[label].e_elect = parser.parse_e_elect(self.output[label]['sp'], + zpe_scale_factor=zpe_scale_factor) if self.species_dict[label].t1 is not None: txt = '' if self.species_dict[label].t1 > 0.02: diff --git a/arc/settings.py b/arc/settings.py index 4c1accb925..566a6f50f0 100644 --- a/arc/settings.py +++ b/arc/settings.py @@ -113,14 +113,14 @@ 'onedmin': 'output.out', } -default_levels_of_theory = {'conformer': 'b3lyp/6-31+g(d,p)', - 'ts_guesses': 'b3lyp/6-31+g(d,p)', # used for IRC as well +default_levels_of_theory = {'conformer': 'b97d3/6-31+g(d,p)', + 'ts_guesses': 'b97d3/6-31+g(d,p)', # used for IRC as well 'opt': 'wb97xd/6-311++g(d,p)', 'freq': 'wb97xd/6-311++g(d,p)', # should be the same level as opt (to calc freq at min E) - 'scan': 'wb97xd/6-311++g(d,p)', # should be the same level as freq (to project out rotors) + 'scan': 'b97d3/6-31+g(d,p)', # should be the same level as freq (to project out rotors) 'sp': 'ccsd(t)-f12/cc-pvtz-f12', # This should be a level for which BAC is available # 'sp': 'b3lyp/6-311+g(3df,2p)', - 'orbitals': 'b3lyp/6-311+g(d,p)', # save orbitals for visualization + 'orbitals': 'wb97x-d3/6-311++g(d,p)', # save orbitals for visualization 'scan_for_composite': 'B3LYP/CBSB7', # This is the frequency level of the CBS-QB3 method 'freq_for_composite': 'B3LYP/CBSB7', # This is the frequency level of the CBS-QB3 method } diff --git a/arc/species/converter.py b/arc/species/converter.py index 0a160651d4..a84e058886 100644 --- a/arc/species/converter.py +++ b/arc/species/converter.py @@ -392,12 +392,27 @@ def s_bonds_mol_from_xyz(xyz): def rdkit_conf_from_mol(mol, coordinates): - """A helper function generating an RDKit:Conformer object from an RMG:Molecule object""" - rd_mol, rd_inds = mol.toRDKitMol(removeHs=False, returnMapping=True) + """A helper function generating an RDKit Conformer object from an RMG Molecule object + + Args: + mol (Molecule): The RMG Molecule object + coordinates (list, str, unicode): Array of xyz coordinates for the conformer + + Returns: + Conformer: An RDKit Conformer object + RDMol: An RDKit Molecule object + dict: Atom index map. Keys are atom indices in the RMG Molecule, values are atom indices in the RDKit Molecule + """ + if isinstance(coordinates[0], (str, unicode)) and isinstance(coordinates, list): + raise InputError('The coordinates argument seem to be of wrong type. Got a list of strings:\n{0}'.format( + coordinates)) + if isinstance(coordinates, (str, unicode)): + coordinates = get_xyz_matrix(xyz=coordinates)[0] + rd_mol, rd_indices = mol.toRDKitMol(removeHs=False, returnMapping=True) Chem.AllChem.EmbedMolecule(rd_mol) # unfortunately, this mandatory embedding changes the coordinates indx_map = dict() for xyz_index, atom in enumerate(mol.atoms): # generate an atom index mapping dictionary - rd_index = rd_inds[atom] + rd_index = rd_indices[atom] indx_map[xyz_index] = rd_index conf = rd_mol.GetConformer(id=0) for i in range(rd_mol.GetNumAtoms()): # reset atom coordinates diff --git a/arc/species/species.py b/arc/species/species.py index 5dc0a8b6dd..1f75624076 100644 --- a/arc/species/species.py +++ b/arc/species/species.py @@ -18,7 +18,7 @@ import pybel as pyb from arkane.common import ArkaneSpecies, symbol_by_number -from arkane.statmech import determine_qm_software +from arkane.statmech import determine_qm_software, is_linear from rmgpy.molecule.converter import toOBMol from rmgpy.molecule.element import getElement from rmgpy.molecule.molecule import Atom, Molecule @@ -56,7 +56,8 @@ class ARCSpecies(object): it). Defaults to None. Important, e.g., if a Species is a bi-rad singlet, in which case the job should be unrestricted, but the multiplicity does not have the required information to make that decision (r vs. u) - `e0` ``float`` The total electronic energy E0 of the species at the chosen sp level (kJ/mol) + `e_elect` ``float`` The total electronic energy (without ZPE) of the species + at the chosen sp level, in kJ/mol `is_ts` ``bool`` Whether or not the species represents a transition state `number_of_rotors` ``int`` The number of potential rotors to scan `rotors_dict` ``dict`` A dictionary of rotors. structure given below. @@ -154,7 +155,7 @@ def __init__(self, is_ts=False, rmg_species=None, mol=None, label=None, xyz=None # Not reading from a dictionary self.is_ts = is_ts self.ts_conf_spawned = False - self.e0 = None + self.e_elect = None self.arkane_file = None if self.is_ts: if ts_methods is None: @@ -330,7 +331,7 @@ def as_dict(self): """A helper function for dumping this object as a dictionary in a YAML file for restarting ARC""" species_dict = dict() species_dict['is_ts'] = self.is_ts - species_dict['E0'] = self.e0 + species_dict['E0'] = self.e_elect species_dict['arkane_file'] = self.arkane_file if self.yml_path is not None: species_dict['yml_path'] = self.yml_path @@ -388,7 +389,7 @@ def from_dict(self, species_dict): raise InputError('All species must have a label') self.run_time = datetime.timedelta(seconds=species_dict['run_time']) if 'run_time' in species_dict else None self.t1 = species_dict['t1'] if 't1' in species_dict else None - self.e0 = species_dict['E0'] if 'E0' in species_dict else None + self.e_elect = species_dict['E electronic'] if 'E electronic' in species_dict else None self.arkane_file = species_dict['arkane_file'] if 'arkane_file' in species_dict else None self.yml_path = species_dict['yml_path'] if 'yml_path' in species_dict else None self.rxn_label = species_dict['rxn_label'] if 'rxn_label' in species_dict else None @@ -531,8 +532,8 @@ def from_yml_file(self, label=None): self.external_symmetry = external_symmetry_mode.symmetry if self.initial_xyz is not None: self.mol_from_xyz() - if self.e0 is None: - self.e0 = arkane_spc.conformer.E0.value_si * 0.001 # convert to kJ/mol + if self.e_elect is None: # TODO: this is actually the E0, not e_elect! be consistent! + self.e_elect = arkane_spc.conformer.E0.value_si * 0.001 # convert to kJ/mol def generate_conformers(self): """ @@ -767,114 +768,6 @@ def determine_multiplicity(self, smiles, adjlist, mol): if self.multiplicity is None: raise SpeciesError('Could not determine multiplicity for species {0}'.format(self.label)) - def is_linear(self): - """ - A helper function for determination of the species linearity from the .final_xyz attribute - """ - epsilon = 0.001 - if self.number_of_atoms == 1: - return False - if self.number_of_atoms == 2: - return True - xyz = self.final_xyz or self.initial_xyz or self.conformers[0] - if not xyz: - raise SpeciesError('Cannot determine linearity for {0} without the initial/final xyz coordinates'.format( - self.label)) - _, _, x, y, z = get_xyz_matrix(xyz) - # first check whether the problem can be reduced into two (or one) dimensions - reduced_coord = '' - if all([nearly_equal(element, x[0]) for element in x]): - reduced_coord += 'x' - if all([nearly_equal(element, y[0]) for element in y]): - reduced_coord += 'y' - if all([nearly_equal(element, z[0]) for element in z]): - reduced_coord += 'z' - if reduced_coord: - # The problem can (and should) be reduced to two dimensions. - # (might actually be buggy to deal with it in 3D due to divisions by zero) - # Use linear interpolation instead of a 3D line - # This is often the case for standard orientation of linear molecules - if len(reduced_coord) == 2: - # This is the trivial case where all atoms are on the same axis by definition - return True - if reduced_coord == 'x': - # don't use x - x = y - y = z - elif reduced_coord == 'y': - # keep x as it is, don't use y - y = z - elif reduced_coord == 'z': - # keep x, y as they are, don't use z - pass - else: - raise SpeciesError('Could not determine whether species {0} is linear, something must be wrong with its' - ' coordinates:\n{1}'.format(self.label, xyz)) - for i, _ in enumerate(x): - # Use a linear interpolation/extrapolation method, and compare extrapolated value and actual value - if i > 1: # only start the algorithm at the third atom (index 2) - if x[0] != x[1]: - y0 = y[0] + (x[i] - x[0]) * (y[1] - y[0]) / (x[1] - x[0]) - if y[i] != 0: - if nearly_equal(a=y0, b=y[i]): - return True - else: - return False - else: - # avoid dividing by 0, but check we do extrapolate to 0 within the tolerance - if abs(y0) < epsilon: - return True - else: - return False - elif y[0] != y[1]: - x0 = x[0] + (y[i] - y[0]) * (x[1] - x[0]) / (y[1] - y[0]) - if x[i] != 0: - if nearly_equal(a=x0, b=x[i]): - return True - else: - return False - else: - if abs(x0) < epsilon: - return True - else: - return False - else: - raise SpeciesError( - 'Could not determine whether species {0} is linear, something must be wrong with its' - ' coordinates:\n{1}'.format(self.label, xyz)) - else: - # this problem cannot be reduced, and must have a 3D treatment - vector_list = [] - for i, _ in enumerate(x): - # generate a vector list where each element in the list contains [v1, v2] - # which are vectors (subtracting two points in space) between atoms A-B and B-C (or i, i-1 and i-1, i-2) - if i > 1: # so only start this at the third atom (index 2) - v1 = [x[i] - x[i - 1], y[i] - y[i - 1], z[i] - z[i - 1]] # between points A, B - v2 = [x[i - 1] - x[i - 2], y[i - 1] - y[i - 2], z[i - 1] - z[i - 2]] # between points B, C - vector_list.append([v1, v2]) - for v1, v2 in vector_list: - if v1[0] and v2[0]: - ax = v1[0] / v2[0] # This is how much the x coord of point C is stretched relative to the AB vector - else: - ax = 0 - if v1[1] and v2[1]: - ay = v1[1] / v2[1] - else: - ay = 0 - if v1[2] and v2[2]: - az = v1[2] / v2[2] - else: - az = 0 - if ax and ay and not nearly_equal(ax, ay): - break - if ax and az and not nearly_equal(ax, az): - break - if ay and az and not nearly_equal(ay, az): - break - else: - return True - return False - def make_ts_report(self): """A helper function to write content into the .ts_report attribute""" self.ts_report = '' @@ -1013,11 +906,14 @@ def set_transport_data(self, lj_path, opt_path, bath_gas, opt_level, freq_path=' if self.number_of_atoms == 1: shape_index = 0 comment += '; The molecule is monoatomic' - elif self.is_linear(): - shape_index = 1 - comment += '; The molecule is linear' else: - shape_index = 2 + coordinates = get_xyz_matrix(self.final_xyz or self.initial_xyz or self.conformers[0])[0] + coordinates = np.array(coordinates) + if is_linear(coordinates): + shape_index = 1 + comment += '; The molecule is linear' + else: + shape_index = 2 if self.number_of_atoms > 1: dipole_moment = parse_dipole_moment(opt_path) or 0 if dipole_moment: diff --git a/arc/species/speciesTest.py b/arc/species/speciesTest.py index b53d739038..3f46ea3517 100644 --- a/arc/species/speciesTest.py +++ b/arc/species/speciesTest.py @@ -235,64 +235,6 @@ def test_xyz_format_conversion(self): self.assertEqual(y[1], 1.29155) self.assertEqual(z[-1], -0.932864) - def test_is_linear(self): - """Test determination of molecule linearity by xyz""" - xyz1 = """C 0.000000 0.000000 0.000000 - O 0.000000 0.000000 1.159076 - O 0.000000 0.000000 -1.159076""" # a trivial case - xyz2 = """S -0.06618943 -0.12360663 -0.07631983 - O -0.79539707 0.86755487 1.02675668 - O -0.68919931 0.25421823 -1.34830853 - N 0.01546439 -1.54297548 0.44580391 - C 1.59721519 0.47861334 0.00711000 - H 1.94428095 0.40772394 1.03719428 - H 2.20318015 -0.14715186 -0.64755729 - H 1.59252246 1.51178950 -0.33908352 - H -0.87856890 -2.02453514 0.38494433 - H -1.34135876 1.49608206 0.53295071""" # a non linear molecule - xyz3 = """N 0.0000000000 0.0000000000 0.3146069129 - O -1.0906813653 0.0000000000 -0.1376405244 - O 1.0906813653 0.0000000000 -0.1376405244""" # a non linear 3-atom molecule - xyz4 = """N 0.0000000000 0.0000000000 0.1413439534 - H -0.8031792912 0.0000000000 -0.4947038368 - H 0.8031792912 0.0000000000 -0.4947038368""" # a non linear 3-atom molecule - xyz5 = """S -0.5417345330 0.8208150346 0.0000000000 - O 0.9206183692 1.6432038228 0.0000000000 - H -1.2739176462 1.9692549926 0.0000000000""" # a non linear 3-atom molecule - xyz6 = """N 1.18784533 0.98526702 0.00000000 - C 0.04124533 0.98526702 0.00000000 - H -1.02875467 0.98526702 0.00000000""" # linear - xyz7 = """C -4.02394116 0.56169428 0.00000000 - H -5.09394116 0.56169428 0.00000000 - C -2.82274116 0.56169428 0.00000000 - H -1.75274116 0.56169428 0.00000000""" # linear - xyz8 = """C -1.02600933 2.12845307 0.00000000 - C -0.77966935 0.95278385 0.00000000 - H -1.23666197 3.17751246 0.00000000 - H -0.56023545 -0.09447399 0.00000000""" # just 0.5 degree off from linearity, so NOT linear... - xyz9 = """O -1.1998 0.1610 0.0275 - O -1.4021 0.6223 -0.8489 - O -1.48302 0.80682 -1.19946""" # just 3 points in space on a straight line (not a physical molecule) - spc1 = ARCSpecies(label=str('test_spc'), xyz=xyz1, multiplicity=1, charge=0, smiles=str('O=C=O')) - spc2 = ARCSpecies(label=str('test_spc'), xyz=xyz2, multiplicity=1, charge=0, smiles=str('[NH-][S+](=O)(O)C')) - spc3 = ARCSpecies(label=str('test_spc'), xyz=xyz3, multiplicity=2, charge=0, smiles=str('[O]N=O')) - spc4 = ARCSpecies(label=str('test_spc'), xyz=xyz4, multiplicity=2, charge=0, smiles=str('[NH2]')) - spc5 = ARCSpecies(label=str('test_spc'), xyz=xyz5, multiplicity=2, charge=0, smiles=str('[O]S')) - spc6 = ARCSpecies(label=str('test_spc'), xyz=xyz6, multiplicity=1, charge=0, smiles=str('C#N')) - spc7 = ARCSpecies(label=str('test_spc'), xyz=xyz7, multiplicity=1, charge=0, smiles=str('C#C')) - spc8 = ARCSpecies(label=str('test_spc'), xyz=xyz8, multiplicity=1, charge=0, smiles=str('C#C')) - spc9 = ARCSpecies(label=str('test_spc'), xyz=xyz9, multiplicity=1, charge=0, smiles=str('[O-][O+]=O')) - - self.assertTrue(spc1.is_linear()) - self.assertTrue(spc6.is_linear()) - self.assertTrue(spc7.is_linear()) - self.assertTrue(spc9.is_linear()) - self.assertFalse(spc2.is_linear()) - self.assertFalse(spc3.is_linear()) - self.assertFalse(spc4.is_linear()) - self.assertFalse(spc5.is_linear()) - self.assertFalse(spc8.is_linear()) - def test_charge_and_multiplicity(self): """Test determination of molecule charge and multiplicity""" spc1 = ARCSpecies(label='spc1', mol=Molecule(SMILES=str('C[CH]C')), generate_thermo=False)