Solve problems in arc's directed scans. (#602)

A few functions returned a `np.float64`, which caused errors when dumped
into yaml files, and several other small corrections.

arc/job/adapters/scripts/ob_script.py

@@ -1,3 +1,6 @@
+#!/usr/bin/env python3
+# encoding: utf-8
+
 from openbabel import openbabel
 import argparse
 import os

arc/job/adapters/scripts/tani_script.py

@@ -329,9 +329,9 @@ def main():
     device = input_dict["device"]
     model = None
     if input_dict["model"].lower() == 'ani1ccx':
-        model = torchani.models.ani1ccx(periodic_table_index=True).to(device)
+        model = torchani.models.ANI1ccx(periodic_table_index=True).to(device)
     elif input_dict["model"].lower() =='ani1x':
-        model = torchani.models.ani1x(periodic_table_index=True).to(device)
+        model = torchani.models.ANI1x(periodic_table_index=True).to(device)
     elif input_dict["model"].lower() == 'ani2x' or model is None:
         model = torchani.models.ANI2x(periodic_table_index=True).to(device)
 
@@ -345,7 +345,7 @@ def main():
         forces = run_force(xyz=xyz, device=device, model=model)
         save_output_file(path = str(args.yml_path), key="force", val=forces)
 
-    elif job_type in ['opt', 'conformers', 'optfreq']:
+    elif job_type in ['opt', 'conformers', 'directed_scan', 'optfreq']:
         constraints = input_dict["constraints"] if "constraints" in input_dict.keys() else None
         opt_xyz = run_opt(xyz=xyz, constraints=constraints, fmax=input_dict["fmax"], model=model,
                           steps=input_dict["steps"] if "steps" in input_dict.keys() else None, engine=input_dict["engine"])

arc/job/adapters/torch_ani.py

@@ -237,8 +237,12 @@ def execute_incore(self):
         self._log_job_execution()
 
         if self.level is not None and self.level.args is not None and any(key in self.level.args for key in ["model", "device", "fmax", "engine", "steps"]):
-                    for key in self.level.args:    
+                    for key in self.level.args:  
                         tani_default_options_dict.update({key : self.level.args[key]})
+        if self.args is not None and self.args["keyword"] is not None \
+            and any(key in self.args["keyword"] for key in ["model", "device", "fmax", "engine", "steps"]):
+                for key in self.args["keyword"]:
+                    tani_default_options_dict.update({key : self.args["keyword"][key]})
 
         if not self.check_settings(tani_default_options_dict):
             logger.error(f'given wrong settings in input: {self.level.args} or '

arc/main.py

@@ -996,10 +996,11 @@ def set_levels_of_theory(self):
                     default_flag = ' (default)'
                 else:
                     default_flag = ''
-                self.opt_level = Level(repr=self.opt_level)
                 logger.info(f'Geometry optimization:{default_flag} {self.opt_level}')
             else:
                 logger.warning("Not performing geometry optimization, since it was not requested by the user.")
+            if self.opt_level:
+                self.opt_level = Level(repr=self.opt_level)
 
             if self.job_types['freq']:
                 if not self.freq_level:
@@ -1011,21 +1012,23 @@ def set_levels_of_theory(self):
                         info = ' (default)'
                 else:
                     info = ''
-                self.freq_level = Level(repr=self.freq_level)
                 logger.info(f'Frequencies:{info} {self.freq_level}')
             else:
                 logger.warning("Not performing frequency calculation, since it was not requested by the user.")
+            if self.freq_level:
+                self.freq_level = Level(repr=self.freq_level)
 
             if self.job_types['sp']:
                 if not self.sp_level:
                     self.sp_level = default_levels_of_theory['sp']
                     default_flag = ' (default)'
                 else:
                     default_flag = ''
-                self.sp_level = Level(repr=self.sp_level)
                 logger.info(f'Energy:{default_flag} {self.sp_level}')
             else:
                 logger.warning("Not performing single point energy calculation, since it was not requested by the user.")
+            if self.sp_level:
+                self.sp_level = Level(repr=self.sp_level)
 
             if self.job_types['rotors']:
                 if not self.scan_level:
@@ -1037,35 +1040,38 @@ def set_levels_of_theory(self):
                         info = ' (default)'
                 else:
                     info = ''
-                self.scan_level = Level(repr=self.scan_level)
                 logger.info(f'Rotor scans:{info} {self.scan_level}')
             else:
                 logger.warning("Not performing rotor scans, since it was not requested by the user. This might "
                                "compromise finding the best conformer, as dihedral angles won't be corrected. "
                                "Also, the calculated thermodynamic properties and rate coefficients "
                                "will be less accurate.")
+            if self.scan_level:
+                self.scan_level = Level(repr=self.scan_level)
 
             if self.job_types['irc']:
                 if not self.irc_level:
                     self.irc_level = default_levels_of_theory['irc']
                     default_flag = ' (default)'
                 else:
                     default_flag = ''
-                self.irc_level = Level(repr=self.irc_level)
                 logger.info(f'IRC:{default_flag} {self.irc_level}')
             else:
                 logger.warning("Not running IRC computations, since it was not requested by the user.")
+            if self.irc_level:
+                self.irc_level = Level(repr=self.irc_level)
 
             if self.job_types['orbitals']:
                 if not self.orbitals_level:
                     self.orbitals_level = default_levels_of_theory['orbitals']
                     default_flag = ' (default)'
                 else:
                     default_flag = ''
-                self.orbitals_level = Level(repr=self.orbitals_level)
                 logger.info(f'Orbitals:{default_flag} {self.orbitals_level}')
             else:
                 logger.debug("Not running molecular orbitals visualization, since it was not requested by the user.")
+            if self.orbitals_level:
+                self.orbitals_level = Level(repr=self.orbitals_level)
 
         if not self.job_types['fine'] and (self.opt_level is not None and self.opt_level.method_type == 'dft'
                                            or self.composite_method is not None):

arc/main_test.py

@@ -12,6 +12,7 @@
 from arc.common import ARC_PATH
 from arc.exceptions import InputError
 from arc.imports import settings
+from arc.level import Level
 from arc.main import ARC, StatmechEnum, process_adaptive_levels
 from arc.species.species import ARCSpecies
 
@@ -440,6 +441,40 @@ def test_process_adaptive_levels(self):
                                      (15, 'inf'): {('opt', 'freq'): 'b3lyp/6-31g(d,p)',
                                                    'sp': 'b3lyp/6-311+g(d,p)'}})
 
+    def test_process_level_of_theory(self):
+        """
+        Tests the process_level_of_theory function.
+        """
+        arc0 = ARC(project='test_0', level_of_theory='ccsd(t)-f12/cc-pvdz-f12//b3lyp/6-311+g(3df,2p)', bac_type=False, freq_scale_factor=1,)
+        arc1 = ARC(project='test_1', level_of_theory='wb97xd/6-311+g(2d,2p)',
+                                    arkane_level_of_theory="b3lyp/6-311+g(3df,2p)",
+                                    bac_type=False,
+                                    freq_scale_factor=1,
+                                    job_types= {"freq": True,
+                                                "sp"  : True,
+                                                "opt" : False,})
+        arc2 = ARC(project='test_2', sp_level='wb97xd/6-311+g(2d,2p)',
+                                     opt_level='wb97xd/6-311+g(2d,2p)',
+                                     arkane_level_of_theory="b3lyp/6-311+g(3df,2p)",
+                                     bac_type=False,
+                                     freq_scale_factor=1,
+                                     job_types= {"freq": True,
+                                                 "sp"  : False,
+                                                 "opt" : False,})
+        arc3 = ARC(project='test_3', sp_level='wb97xd/6-311+g(2d,2p)',
+                                     opt_level='wb97xd/6-311+g(2d,2p)',
+                                     arkane_level_of_theory="b3lyp/6-311+g(3df,2p)",
+                                     bac_type=False,
+                                     freq_scale_factor=1,
+                                     job_types= {"opt" : False,})
+
+        arc0.process_level_of_theory(), arc1.process_level_of_theory(), arc2.process_level_of_theory(), arc3.process_level_of_theory()
+        for arc in [arc0, arc1, arc2, arc3]:
+            print("arc")
+            self.assertIsInstance(arc.sp_level, Level)
+            self.assertIsInstance(arc.opt_level, Level)
+            self.assertIsInstance(arc.freq_level, Level)
+
     @classmethod
     def tearDownClass(cls):
         """

arc/parser.py

@@ -233,6 +233,8 @@ def parse_geometry(path: str) -> Optional[Dict[str, tuple]]:
         content = read_yaml_file(path)
         if isinstance(content, dict) and 'xyz' in content.keys():
             return content['xyz']
+        if isinstance(content, dict) and 'opt_xyz' in content.keys():
+            return content['opt_xyz']
     software = identify_ess(path)
     xyz_str = ''
     if software == 'xtb':
@@ -737,6 +739,8 @@ def parse_xyz_from_file(path: str) -> Optional[Dict[str, tuple]]:
             xyz = content["xyz"]
         elif "opt_xyz" in content.keys():
             xyz = content["opt_xyz"]
+        if isinstance(xyz, str):
+            xyz = str_to_xyz(xyz)
     else:
         record = False
         for line in lines:

arc/plotter.py

@@ -32,8 +32,9 @@
                         get_close_tuple,
                         get_logger,
                         is_notebook,
+                        is_str_float,
                         save_yaml_file,
-                        sort_two_lists_by_the_first)
+                        sort_two_lists_by_the_first,)
 from arc.exceptions import InputError, SanitizationError
 from arc.level import Level
 from arc.parser import parse_trajectory
@@ -1145,7 +1146,13 @@ def plot_1d_rotor_scan(angles: Optional[Union[list, tuple, np.array]] = None,
     plt.xticks(np.arange(min_angle, min_angle + 361, step=60))
     plt.ylabel('Electronic energy (kJ/mol)')
     if label is not None:
-        original_dihedral_str = f' from {original_dihedral:.1f}$^o$' if original_dihedral is not None else ''
+        if isinstance(original_dihedral, str):
+            if is_str_float(original_dihedral):
+                original_dihedral = float(original_dihedral)
+        if isinstance(original_dihedral, float):
+            original_dihedral_str = f' from {original_dihedral:.1f}$^o$' if original_dihedral is not None else ''
+        else:
+            original_dihedral_str = ""
         plt.title(f'{label} 1D scan of {scan}{original_dihedral_str}')
     plt.tight_layout()
 

arc/scheduler.py

@@ -424,7 +424,10 @@ def __init__(self,
                     if self.job_types['sp']:
                         self.run_sp_job(species.label)
                     if self.job_types['rotors']:
-                        self.run_sp_job(species.label)
+                        self.run_sp_job(species.label, level = self.scan_level)
+                        if not self.job_types['opt']: # The user provided an optimized coordinets
+                            self.run_scan_jobs(species.label)
+
                 elif ((species.initial_xyz is not None or species.final_xyz is not None)
                         or species.is_ts and species.rxn_label is None) and not self.testing:
                     # For restarting purposes: check before running jobs whether they were already terminated
@@ -764,6 +767,8 @@ def run_job(self,
             args['shift'] = shift
         if scan_trsh:
             args['keyword']['scan_trsh'] = scan_trsh
+        if level_of_theory is not None and level_of_theory.args is not None:
+            args.update(level_of_theory.args)
 
         job = job_factory(job_adapter=job_adapter,
                           project=self.project,

arc/species/converter.py

@@ -1901,7 +1901,7 @@ def get_center_of_mass(xyz):
     cm_x /= sum(masses)
     cm_y /= sum(masses)
     cm_z /= sum(masses)
-    return cm_x, cm_y, cm_z
+    return float(cm_x), float(cm_y), float(cm_z)
 
 
 def translate_to_center_of_mass(xyz):

arc/species/converter_test.py

@@ -4009,6 +4009,10 @@ def test_get_center_of_mass(self):
         self.assertAlmostEqual(cm_y, -0.4880, 3)
         self.assertAlmostEqual(cm_z, -0.1603, 3)
 
+        self.assertIsInstance(cm_x, float)
+        self.assertIsInstance(cm_y, float)
+        self.assertIsInstance(cm_z, float)
+
         xyz = """C	1.1714680	-0.4048940	0.0000000
 C	0.0000000	0.5602500	0.0000000
 O	-1.1945070	-0.2236470	0.0000000
@@ -4673,6 +4677,27 @@ def test_modify_coords(self):
                                           fragments=fragments,
                                           )
         self.assertAlmostEqual(calculate_dihedral_angle(coords=new_xyz, torsion=indices, index=1), 200, places=3)
+
+        coords={'coords': ((-0.7862825353221515, -0.28824023055636216, 0.4782944637692894),
+                           (0.21968869054702736, 0.40094256193652866, -0.2919820499085219),
+                           (-0.07796443595084417, 0.5692847962524797, -1.6621913220858304),
+                           (-1.102200211589376, -1.1132157833188596, -0.01879031191901484),
+                           (-1.5973749070505925, 0.29546848172306867, 0.6474145668621136),
+                           (0.4237940503863438, 1.3660724867336205, 0.19101403432872205),
+                           (1.1352054736534014, -0.1980893380251006, -0.2652264470061931),
+                           (-0.7497944593402266, 1.258221857416732, -1.7507029654486272)),
+                'isotopes': (14, 12, 16, 1, 1, 1, 1, 1),
+                'symbols': ('N', 'C', 'O', 'H', 'H', 'H', 'H', 'H')}
+        indices=[3, 0, 1, 2]
+        new_value=53.76
+        modification_type="groups"
+        mol=Molecule(smiles="NCO")
+        new_xyz = converter.modify_coords(coords=coords,
+                                          indices=indices,
+                                          new_value=new_value,
+                                          modification_type=modification_type,
+                                          mol=mol)
+        self.assertTrue(type(new_xyz["coords"][0][0] is float))
 
     def test_relocate_zmat_dummy_atoms_to_the_end(self):
         """Test the relocating dummy atoms in a zmat to the end of the cartesian coordinates atom list"""

arc/species/zmat.py

@@ -1113,7 +1113,7 @@ def _add_nth_atom_to_coords(zmat: dict,
                       cd_length * math.sin(bcd_angle) * math.sin(abcd_dihedral)])
         d = m.dot(d)  # Rotate the coordinate system into the reference frame of orientation defined by A, B, C.
         # Add the coordinates of atom C to the resulting atom D:
-        coords.append((d[0] + coords[c_index][0], d[1] + coords[c_index][1], d[2] + coords[c_index][2]))
+        coords.append((float(d[0] + coords[c_index][0]), float(d[1] + coords[c_index][1]), float(d[2] + coords[c_index][2])))
     return coords