diff --git a/aqme/csearch/base.py b/aqme/csearch/base.py index 28b68fc9..536d2c2d 100644 --- a/aqme/csearch/base.py +++ b/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 diff --git a/aqme/csearch/utils.py b/aqme/csearch/utils.py index 317dbba4..134b041d 100644 --- a/aqme/csearch/utils.py +++ b/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, diff --git a/aqme/utils.py b/aqme/utils.py index d41ad544..6e04b1cb 100644 --- a/aqme/utils.py +++ b/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 diff --git a/tests/test_csearch.py b/tests/test_csearch.py index f5980843..a6569df9 100644 --- a/tests/test_csearch.py +++ b/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) diff --git a/tests/test_qcorr.py b/tests/test_qcorr.py index ba991916..db810aa9 100644 --- a/tests/test_qcorr.py +++ b/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(