1. Fix commands with list

aqme/csearch/base.py

@@ -753,29 +753,38 @@ def summ_search(
                     )
 
         if self.args.program.lower() in ['crest']:
+            stop_xtb_opt = False
             if not complex_ts:
                 # mol_crest is the RDKit-optimized mol object
-                rdmolfiles.MolToXYZFile(mol_crest, name + "_crest.xyz")
+                if mol_crest is not None:
+                    rdmolfiles.MolToXYZFile(mol_crest, name + "_crest.xyz")
+                else:
+                    stop_xtb_opt = True
+                    status = -1
             else:
                 # mol is the raw mol object (no optimization with RDKit to avoid problems when using
                 # noncovalent complexes and TSs)
-                rdmolfiles.MolToXYZFile(mol, name + "_crest.xyz")
-
-            status = xtb_opt_main(
-                    f"{name}_{self.args.program.lower()}",
-                    dup_data,
-                    dup_data_idx,
-                    self,
-                    charge,
-                    mult,
-                    constraints_atoms,
-                    constraints_dist,
-                    constraints_angle,
-                    constraints_dihedral,
-                    'crest',
-                    complex_ts=complex_ts,
-                    mol=mol # this is necessary for CREST calculations with constraints
-                )
+                if mol is not None:
+                    rdmolfiles.MolToXYZFile(mol, name + "_crest.xyz")
+                else:
+                    stop_xtb_opt = True
+                    status = -1
+            if not stop_xtb_opt:
+                status = xtb_opt_main(
+                        f"{name}_{self.args.program.lower()}",
+                        dup_data,
+                        dup_data_idx,
+                        self,
+                        charge,
+                        mult,
+                        constraints_atoms,
+                        constraints_dist,
+                        constraints_angle,
+                        constraints_dihedral,
+                        'crest',
+                        complex_ts=complex_ts,
+                        mol=mol # this is necessary for CREST calculations with constraints
+                    )
 
         return status
 
@@ -1189,8 +1198,13 @@ def rdkit_to_sdf(
         Conversion from RDKit to SDF
         """
 
-        Chem.SanitizeMol(mol)
-        mol = Chem.AddHs(mol)
+        try:
+            Chem.SanitizeMol(mol)
+            mol = Chem.AddHs(mol)
+        except Chem.AtomValenceException: # this happens sometimes with complex metals when substituting the metal with an I atom
+            self.args.log.write(f'\nx  The species provided could not be converted into a mol object wth RDKit. It normally happens with tricky metal complexes and might be fixed with different structures (i.e., changing a single bond + positive charge with a double bond).')
+            return -1, None, None, None
+
         mol.SetProp("_Name", name)
 
         # detects and applies auto-detection of initial number of conformers

aqme/csearch/utils.py

@@ -501,8 +501,8 @@ def smi_to_mol(
         try:
             mol = Chem.MolFromSmiles(smi[0], params)
         except Chem.AtomValenceException:
-            log.write(f"\nx  The SMILES string provided ( {smi[0]} ) contains errors. For example, N atoms from ligands of metal complexes should be N+ since they're drawn with four bonds in ChemDraw, same for O atoms in carbonyl ligands, etc.\n")
-            sys.exit()
+            log.write(f"\nx  The SMILES string provided ( {smi[0]} ) contains errors or the molecule needs to be drawn in a different way. For example, N atoms from ligands of metal complexes should be N+ since they're drawn with four bonds in ChemDraw, same for O atoms in carbonyl ligands, etc.\n")
+            mol = None
 
     return (
         mol,

aqme/utils.py

@@ -240,14 +240,15 @@ def rules_get_charge(mol, args):
 
     M_ligands, N_carbenes, bridge_atoms, C_accounted, neighbours = [], [], [], [], []
     charge_rules = np.zeros(len(mol.GetAtoms()), dtype=int)
-    neighbours, metal_found, sanit_step = [], False, True
+    neighbours, metal_found = [], False
+    try:
+        Chem.SanitizeMol(mol)
+    except Chem.AtomValenceException: # this happens sometimes with complex metals when substituting the metal with an I atom
+        args.log.write(f'\nx  The charge can not be safely calculated for the system provided. If the charge is not right, you can assign it manually with charge=INT.')
     for i, atom in enumerate(mol.GetAtoms()):
         # get the neighbours of metal atom and calculate the charge of metal center + ligands
         if atom.GetIdx() in args.metal_idx:
             # a sanitation step is needed to ensure that metals and ligands show correct valences
-            if sanit_step:
-                Chem.SanitizeMol(mol)
-                sanit_step = False
             metal_found = True
             charge_idx = args.metal_idx.index(atom.GetIdx())
             neighbours = atom.GetNeighbors()
@@ -473,7 +474,11 @@ def command_line_args():
                         try:
                             value = ast.literal_eval(value)
                         except (SyntaxError, ValueError):
-                            pass
+                            # this line fixes issues when using "[X]" or ["X"] instead of "['X']" when using lists
+                            if arg_name.lower() in ["files", "metal_oxi", "metal_atoms", "gen_atoms"]:
+                                value = value.replace('[',']').replace(',',']').split(']')
+                                while('' in value):
+                                    value.remove('')
                 kwargs[arg_name] = value
 
     # Second, load all the default variables as an "add_option" object

tests/test_csearch.py

@@ -439,6 +439,15 @@ def test_double_bond_chrg(
     file = str("CSEARCH/" + name + "_" + program + ".sdf")
     mols = rdkit.Chem.SDMolSupplier(file, removeHs=False)
     assert charge == int(mols[0].GetProp("Real charge"))
+    # check that H atoms are included
+    outfile = open(file,"r")
+    outlines = outfile.readlines()
+    outfile.close()
+    Hatoms_found = False
+    for line in outlines:
+        if "H   0" in line:
+            Hatoms_found = True
+    assert Hatoms_found
     os.chdir(w_dir_main)
 
 
@@ -579,6 +588,15 @@ def test_csearch_rdkit_summ_parameters(
     assert len(mols) == output_nummols
     assert charge == int(mols[0].GetProp("Real charge"))
     assert mult == int(mols[0].GetProp("Mult"))
+    # check that H atoms are included
+    outfile = open(file,"r")
+    outlines = outfile.readlines()
+    outfile.close()
+    Hatoms_found = False
+    for line in outlines:
+        if "H   0" in line:
+            Hatoms_found = True
+    assert Hatoms_found    
     os.chdir(w_dir_main)
 
 

tests/test_qcorr.py

@@ -865,6 +865,7 @@ def test_QCORR_analysis(init_folder, file, command_line, target_folder, restore_
             ]
             subprocess.run(cmd_aqme)
 
+        # the energy of the json file was changed to avoid the duplicate filter
         if file.split(".")[0] == "CH4":
             assert path.exists(f"{w_dir_QCORR_5b}/{target_folder}/{file}")
             assert path.exists(