Merge pull request #406 from ReactionMechanismGenerator/solvation_sp_2

Follow the solvation scheme if one was requested

arc/common.py

@@ -368,7 +368,7 @@ def save_yaml_file(path: str,
     yaml.add_representer(str, string_representer)
     logger.debug('Creating a restart file...')
     content = yaml.dump(data=content)
-    if os.path.dirname(path) and not os.path.exists(os.path.dirname(path)):
+    if '/' in path and os.path.dirname(path) and not os.path.exists(os.path.dirname(path)):
         os.makedirs(os.path.dirname(path))
     with open(path, 'w') as f:
         f.write(content)

arc/job/inputs.py

@@ -181,6 +181,26 @@
 {rotors}
 
 """,
+
+    'arkane_input_species_explicit_e': """#!/usr/bin/env python3
+# encoding: utf-8
+
+{bonds}externalSymmetry = {symmetry}
+
+spinMultiplicity = {multiplicity}
+
+opticalIsomers = {optical}
+
+energy = {e_elect}
+
+geometry = Log('{opt_path}')
+
+frequencies = Log('{freq_path}')
+
+{rotors}
+
+""",
+
     'arkane_hindered_rotor':
         """HinderedRotor(scanLog=Log('{rotor_path}'), pivots={pivots}, top={top}, symmetry={symmetry}, fit='fourier')""",
 

arc/level.py

@@ -144,6 +144,15 @@ def __str__(self) -> str:
             str_ += f' ({self.method_type})'
         return str_
 
+    def copy(self):
+        """
+        A method to create a copy of the object.
+
+        Returns:
+            Level: A copy of the object.
+        """
+        return Level(repr=self.as_dict())
+
     def simple(self) -> str:
         """
         Return a simple humane-readable string representation of the object.

arc/levelTest.py

@@ -151,6 +151,13 @@ def test_ess_methods_yml(self):
                 self.assertIsInstance(method, str)
         self.assertIn('gaussian', list(ess_methods.keys()))
 
+    def test_copy(self):
+        """Test copying the object"""
+        level_1 = Level(repr='wB97xd/def2-tzvp')
+        level_2 = level_1.copy()
+        self.assertIsNot(level_1, level_2)
+        self.assertEqual(level_1.as_dict(), level_2.as_dict())
+
 
 if __name__ == '__main__':
     unittest.main(testRunner=unittest.TextTestRunner(verbosity=2))

arc/scheduler.py

@@ -1035,30 +1035,42 @@ def run_freq_job(self, label):
             self.run_job(label=label, xyz=self.species_dict[label].get_xyz(generate=False),
                          level_of_theory=self.freq_level, job_type='freq')
 
-    def run_sp_job(self, label):
+    def run_sp_job(self,
+                   label: str,
+                   level: Optional[Level] = None,
+                   ):
         """
         Spawn a single point job using 'final_xyz' for species ot TS 'label'.
         If the method is MRCI, first spawn a simple CCSD job, and use orbital determination to run the MRCI job.
 
         Args:
             label (str): The species label.
+            level (Level): An alternative level of theory to run at. If ``None``, self.sp_level will be used.
         """
+        level = level or self.sp_level
+
         # determine_occ(xyz=self.xyz, charge=self.charge)
-        if self.sp_level == self.opt_level and not self.composite_method \
+        if level == self.opt_level and not self.composite_method \
                 and 'paths' in self.output[label] and 'geo' in self.output[label]['paths'] \
                 and self.output[label]['paths']['geo']:
-            logger.info(f'Not running an sp job for {label} at {self.sp_level} since the optimization was done at the '
+            logger.info(f'Not running an sp job for {label} at {level} since the optimization was done at the '
                         f'same level of theory. Using the optimization output to parse the sp energy.')
             recent_opt_job_name, recent_opt_job = 'opt_a0', None
             if 'opt' in self.job_dict[label]:
                 for opt_job_name, opt_job in self.job_dict[label]['opt'].items():
                     if int(opt_job_name.split('_a')[-1]) > int(recent_opt_job_name.split('_a')[-1]):
                         recent_opt_job_name, recent_opt_job = opt_job_name, opt_job
-                self.post_sp_actions(label=label, sp_path=os.path.join(recent_opt_job.local_path, 'output.out'))
+                self.post_sp_actions(label=label,
+                                     sp_path=os.path.join(recent_opt_job.local_path, 'output.out'),
+                                     level=level,
+                                     )
 
             # If opt is not in the job dictionary, the likely explanation is this job has been restarted
             elif 'geo' in self.output[label]['paths']:  # Then just use this path directly
-                self.post_sp_actions(label=label, sp_path=self.output[label]['paths']['geo'])
+                self.post_sp_actions(label=label,
+                                     sp_path=self.output[label]['paths']['geo'],
+                                     level=level,
+                                     )
 
             else:
                 raise RuntimeError(f'Unable to set the path for the sp job for species {label}')
@@ -1069,7 +1081,7 @@ def run_sp_job(self, label):
             self.job_dict[label]['sp'] = dict()
         if self.composite_method:
             raise SchedulerError(f'run_sp_job() was called for {label} which has a composite method level of theory')
-        if 'mrci' in self.sp_level.method:
+        if 'mrci' in level.method:
             if self.job_dict[label]['sp']:
                 # Parse orbital information from the CCSD job, then run MRCI
                 job0 = None
@@ -1107,8 +1119,9 @@ def run_sp_job(self, label):
         if self.job_types['sp']:
             self.run_job(label=label,
                          xyz=self.species_dict[label].get_xyz(generate=False),
-                         level_of_theory=self.sp_level,
-                         job_type='sp')
+                         level_of_theory=level,
+                         job_type='sp',
+                         )
 
     def run_scan_jobs(self, label):
         """
@@ -1315,6 +1328,7 @@ def spawn_post_opt_jobs(self,
                     # this is an 'optfreq' job type, don't run freq
                     self.check_freq_job(label=label, job=self.job_dict[label]['optfreq'][job_name])
             if not composite:
+                # don't use a solvation correction for this sp job if a solvation_scheme_level was specified
                 self.run_sp_job(label)
             if self.species_dict[label].mol is None:
                 # useful for TS species where xyz might not be given to perceive a .mol attribute,
@@ -2233,7 +2247,10 @@ def check_sp_job(self, label, job):
             # This is a CCSD job ran before MRCI. Spawn MRCI
             self.run_sp_job(label)
         elif job.job_status[1]['status'] == 'done':
-            self.post_sp_actions(label, sp_path=os.path.join(job.local_path, 'output.out'))
+            self.post_sp_actions(label,
+                                 sp_path=os.path.join(job.local_path, 'output.out'),
+                                 level=job.level,
+                                 )
             # Update restart dictionary and save the yaml restart file:
             self.save_restart_dict()
             if self.species_dict[label].number_of_atoms == 1:
@@ -2244,15 +2261,20 @@ def check_sp_job(self, label, job):
                                   job=job,
                                   level_of_theory=job.level)
 
-    def post_sp_actions(self, label, sp_path):
+    def post_sp_actions(self,
+                        label: str,
+                        sp_path: str,
+                        level: Optional[Level] = None,
+                        ):
         """
         Perform post-sp actions.
 
         Args:
             label (str): The species label.
-            sp_path (str): The path to 'output.out' for the single point job
+            sp_path (str): The path to 'output.out' for the single point job.
+            level (Level, optional): The level of theory used for the sp job.
         """
-        self.output[label]['job_types']['sp'] = True
+        original_sp_path = self.output[label]['paths']['sp'] if 'sp' in self.output[label]['paths'] else None
         self.output[label]['paths']['sp'] = sp_path
         if self.sp_level is not None and 'ccsd' in self.sp_level.method:
             self.species_dict[label].t1 = parser.parse_t1(self.output[label]['paths']['sp'])
@@ -2269,6 +2291,27 @@ def post_sp_actions(self, label, sp_path):
             logger.info(f'Species {label} has a T1 diagnostic parameter of {self.species_dict[label].t1}{txt}')
             self.output[label]['info'] += f'T1 = {self.species_dict[label].t1}; '
 
+        if self.sp_level.solvation_scheme_level is not None:
+            # a complex solvation correction behavior was requested for the single-point energy value
+            if not self.output[label]['job_types']['sp']:
+                # this is the first "original" sp job, spawn two more at the sp_level.solvation_scheme_level level,
+                # with and without solvation corrections
+                solvation_sp_level = self.sp_level.solvation_scheme_level.copy()
+                solvation_sp_level.solvation_method = self.sp_level.solvation_method
+                solvation_sp_level.solvent = self.sp_level.solvent
+                self.run_sp_job(label=label, level=solvation_sp_level)
+                self.run_sp_job(label=label, level=self.sp_level.solvation_scheme_level)
+            else:
+                # this is one of the additional sp jobs spawned by the above previously
+                if level is not None and level.solvation_method is not None:
+                    self.output[label]['paths']['sp_sol'] = sp_path
+                else:
+                    self.output[label]['paths']['sp_no_sol'] = sp_path
+                self.output[label]['paths']['sp'] = original_sp_path  # restore the original path
+
+        # set *at the end* to differentiate between sp jobs when using complex solvation corrections
+        self.output[label]['job_types']['sp'] = True
+
     def check_irc_job(self, label, job):
         """
         Check an IRC job and perform post-job tasks.

arc/species/converter.py

@@ -1282,8 +1282,9 @@ def order_atoms_in_mol_list(ref_mol, mol_list):
             try:  # TODO: flag as unordered (or solve)
                 order_atoms(ref_mol, mol)
             except SanitizationError as e:
-                logger.warning(f'Could not order atoms in\n{mol.copy(deep=True).to_adjacency_list}'
-                               f'\nGot the following error:\n{e}')
+                logger.warning(f'Could not order atoms in\n'
+                               f'{mol.copy(deep=True).to_adjacency_list()}\n'
+                               f'Got the following error:\n{e}')
                 return False
     else:
         logger.warning('Could not order atoms')

arc/species/species.py

@@ -1629,6 +1629,8 @@ def _scissors(self, indices: tuple) -> list:
         if not mol_copy.has_bond(atom1, atom2):
             raise SpeciesError('Attempted to remove a nonexistent bond.')
         bond = mol_copy.get_bond(atom1, atom2)
+        if not bond.is_single():
+            logger.warning(f'Scissors were requested to remove a non-single bond in {self.label}.')
         mol_copy.remove_bond(bond)
         mol_splits = mol_copy.split()
         if len(mol_splits) == 2:

arc/statmech/arkane.py

@@ -16,6 +16,9 @@
 from arkane.kinetics import KineticsJob
 from arkane.statmech import StatMechJob
 from arkane.thermo import ThermoJob
+from arkane.ess import ess_factory
+
+import rmgpy.constants as constants
 
 import arc.plotter as plotter
 from arc.common import get_logger
@@ -166,7 +169,7 @@ def compute_high_p_rate_coefficient(self) -> None:
         Populates the reaction.kinetics attribute.
         """
         ts_species = self.species_dict[self.reaction.ts_label]
-        if self.output_dict[ts_species.label]['convergence'] and self.reaction.check_ts():
+        if self.output_dict[ts_species.label]['convergence']:
             success = True
             arkane_output_path = self.generate_arkane_species_file(species=ts_species,
                                                                    bac_type=None)
@@ -182,6 +185,7 @@ def compute_high_p_rate_coefficient(self) -> None:
                              f'of reaction {self.reaction.label}')
             else:
                 ts_species.e0 = arkane_ts_species.conformer.E0.value_si * 0.001  # convert to kJ/mol
+                self.reaction.check_ts()
                 self.reaction.dh_rxn298 = \
                     sum([product.thermo.get_enthalpy(298) * self.reaction.get_species_count(product, well=1)
                          for product in self.reaction.p_species]) \
@@ -386,22 +390,48 @@ def generate_arkane_species_file(self,
         # write the Arkane species input file
         bac_txt = '' if bac_type is not None else '_no_BAC'
         input_file_path = os.path.join(species_folder_path, f'{species.label}_arkane_input{bac_txt}.py')
-        input_file = input_files['arkane_input_species']
+        input_file = input_files['arkane_input_species'] if 'sp_sol' not in self.output_dict[species.label]['paths'] \
+            else input_files['arkane_input_species_explicit_e']
         if bac_type is not None and not species.is_ts:
             logger.info(f'Using the following BAC (type {bac_type}) for {species.label}: {species.bond_corrections}')
             bonds = f'bonds = {species.bond_corrections}\n\n'
         else:
             logger.debug(f'NOT using BAC for {species.label}')
             bonds = ''
 
-        input_file = input_file.format(bonds=bonds,
-                                       symmetry=species.external_symmetry,
-                                       multiplicity=species.multiplicity,
-                                       optical=species.optical_isomers,
-                                       sp_path=sp_path,
-                                       opt_path=opt_path,
-                                       freq_path=freq_path,
-                                       rotors=rotors)
+        if 'sp_sol' not in self.output_dict[species.label]['paths']:
+            input_file = input_file.format(bonds=bonds,
+                                           symmetry=species.external_symmetry,
+                                           multiplicity=species.multiplicity,
+                                           optical=species.optical_isomers,
+                                           sp_path=sp_path,
+                                           opt_path=opt_path,
+                                           freq_path=freq_path,
+                                           rotors=rotors)
+        else:
+            # e_elect = e_original + sp_e_sol_corrected - sp_e_uncorrected
+            original_log = ess_factory(self.output_dict[species.label]['paths']['sp'])
+            e_original = original_log.load_energy()
+            e_sol_log = ess_factory(self.output_dict[species.label]['paths']['sp_sol'])
+            e_sol = e_sol_log.load_energy()
+            e_no_sol_log = ess_factory(self.output_dict[species.label]['paths']['sp_no_sol'])
+            e_no_sol = e_no_sol_log.load_energy()
+            e_elect = (e_original + e_sol - e_no_sol) / (constants.E_h * constants.Na)  # convert J/mol to Hartree
+            logger.info(f'\nSolvation correction scheme for {species.label}:\n'
+                        f'Original electronic energy: {e_original * 0.001} kJ/mol\n'
+                        f'Solvation correction: {(e_sol - e_no_sol) * 0.001} kJ/mol\n'
+                        f'New electronic energy: {(e_original + e_sol - e_no_sol) * 0.001} kJ/mol\n\n')
+            print(f'e_elect final: {(e_original + e_sol - e_no_sol) * 0.001} kJ/mol\n\n')
+            input_file = input_files['arkane_input_species_explicit_e']
+            input_file = input_file.format(bonds=bonds,
+                                           symmetry=species.external_symmetry,
+                                           multiplicity=species.multiplicity,
+                                           optical=species.optical_isomers,
+                                           sp_level=self.sp_level,
+                                           e_elect=e_elect,
+                                           opt_path=opt_path,
+                                           freq_path=freq_path,
+                                           rotors=rotors)
 
         if freq_path:
             with open(input_file_path, 'w') as f: