Merge pull request #123 from jvalegre/jv_branch

v1.4.3

aqme/csearch/base.py

@@ -694,9 +694,6 @@ def conformer_generation(
             error_message = "\nx  ERROR: The structure is not valid or no conformers were obtained from this SMILES string"
             self.args.log.write(error_message)
 
-        # if self.args.program.lower() == "crest" and valid_structure:
-        #     shutil.rmtree(f"{self.csearch_folder}/../crest_xyz")
-
         n_seconds = round(time.time() - start_time, 2)
         dup_data.at[dup_data_idx, "CSEARCH time (seconds)"] = n_seconds
 
@@ -935,18 +932,21 @@ def genConformer_r(
 
                 # setting the metal back instead of I
                 set_metal_atomic_number(mol, self.args.metal_idx, self.args.metal_sym)
-            try:
-                sdwriter.write(mol, conf)
-            except (TypeError):
-                raise
-
+
             # if CREST is used, this RDKit preoptimzed mol object will be employed to initializethe the trajectories
             if self.args.program.lower() in ["crest"]:
                 return mol
             else:
+                try:
+                    sdwriter.write(mol, conf)
+                except (TypeError):
+                    raise
                 return 1
 
-        if self.args.program.lower() in ["crest"]:
+        elif self.args.program.lower() in ["crest"]:
+            # setting the metal back instead of I
+            set_metal_atomic_number(mol, self.args.metal_idx, self.args.metal_sym)
+
             return mol
 
         # when SUMM is selected, this cycle generates conformers based on rotation of dihedral angles

aqme/utils.py

@@ -24,7 +24,7 @@
 J_TO_AU = 4.184 * 627.509541 * 1000.0  # UNIT CONVERSION
 T = 298.15
 
-aqme_version = "1.4.2"
+aqme_version = "1.4.3"
 time_run = time.strftime("%Y/%m/%d %H:%M:%S", time.localtime())
 aqme_ref = f"AQME v {aqme_version}, Alegre-Requena, J. V.; Sowndarya, S.; Perez-Soto, R.; Alturaifi, T. M.; Paton, R. S., 2022. https://github.com/jvalegre/aqme"
 
@@ -232,13 +232,13 @@ def rules_get_charge(mol, args):
     Automatically sets the charge for metal complexes
     """
 
-    C_group = ["C", "Se", "Ge"]
-    N_group = ["N", "P", "As"]
-    O_group = ["O", "S", "Se"]
+    C_group = ["C", "Si", "Ge", "Sn"]
+    N_group = ["N", "P", "As", "Sb"]
+    O_group = ["O", "S", "Se", "Te"]
     F_group = ["F", "Cl", "Br", "I"]
 
 
-    M_ligands, N_carbenes, bridge_atoms, neighbours = [], [], [], []
+    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
     for i, atom in enumerate(mol.GetAtoms()):
@@ -253,12 +253,32 @@ def rules_get_charge(mol, args):
             neighbours = atom.GetNeighbors()
             charge_rules[i] = args.metal_oxi[charge_idx]
             for neighbour in neighbours:
+                double_bond = False
                 M_ligands.append(neighbour.GetIdx())
                 if neighbour.GetTotalValence() == 4:
                     if neighbour.GetSymbol() in C_group:
+                        # correct for C=C interacting with M with pi interactions
+                        # when drawing these rings in ChemDraw, all the aromatic C atoms will be linked
+                        # to the M atom as part of 3 member rings
+                        atom_rings = mol.GetRingInfo().AtomRings()
+                        for ring in atom_rings:
+                            if neighbour.GetIdx() in list(ring) and args.metal_idx[0] in list(ring):
+                                # for each double bond
+                                if len(ring) == 3:
+                                    C_count, C_accounted_indiv = 0,[]
+                                    for ring_member in ring:
+                                        if mol.GetAtoms()[ring_member].GetSymbol() == "C" and ring_member not in C_accounted:
+                                            C_accounted_indiv.append(ring_member)
+                                            C_count += 1
+                                    if C_count == 2:
+                                        for ele in C_accounted_indiv:
+                                            C_accounted.append(ele)
+                                        break
                         # first, detects carbenes to adjust charge
                         carbene_like = False
                         bridge_ligand = False
+                        if neighbour.GetIdx() in C_accounted:
+                            double_bond = True
                         for inside_neighbour in neighbour.GetNeighbors():
                             if inside_neighbour.GetSymbol() in N_group:
                                 if inside_neighbour.GetTotalValence() == 4:
@@ -271,7 +291,7 @@ def rules_get_charge(mol, args):
                                     if not bridge_ligand:
                                         carbene_like = True
                                         N_carbenes.append(inside_neighbour.GetIdx())
-                        if not carbene_like:
+                        if not carbene_like and not double_bond:
                             charge_rules[i] = charge_rules[i] - 1
                 elif neighbour.GetTotalValence() == 3:
                     if neighbour.GetSymbol() in N_group and neighbour.GetFormalCharge() == 0:
@@ -288,7 +308,7 @@ def rules_get_charge(mol, args):
                         charge_rules[i] = charge_rules[i] - 2
                     if neighbour.GetSymbol() in F_group:
                         charge_rules[i] = charge_rules[i] - 1
-
+                                            
     # for charges not in the metal, neighbours or exceptions (i.e., C=N+ from carbenes or CN from bridge atoms)
     invalid_charged_atoms = M_ligands + N_carbenes + bridge_atoms + args.metal_idx
     for i, atom in enumerate(mol.GetAtoms()):

docs/Misc/versions.rst

@@ -4,8 +4,15 @@
 Versions
 ========
 
+Version 1.4.3 [`url <https://github.com/jvalegre/aqme/releases/tag/1.4.2>`__]
+   -  Return metal into RDKit mol object when using the metal_atoms option with CSEARCH-CREST
+   -  Doubles bonds do not add extra charges in metal complexes when using the automated charge 
+      calculation from SMILES
+   -  Deprotonated SiR3 groups add -1 charge to metal complexes when using the automated charge 
+      calculation from SMILES
+
 Version 1.4.2 [`url <https://github.com/jvalegre/aqme/releases/tag/1.4.2>`__]
-   -  Fixed an error that raised when using CSEARCH-CREST with organic molecules.
+   -  Fixed an error that raised when using CSEARCH-CREST with organic molecules
    -  Adding more information printed when running CSEARCH
    -  Updated README with citations from external programs
    -  Fixed a bug during filtering of xTB conformers in CMIN (using kcal/mol instead of Hartree

meta.yaml

@@ -1,5 +1,5 @@
 {% set name = "aqme" %}
-{% set version = "1.4.2" %}
+{% set version = "1.4.3" %}
 
 package:
   name: {{ name|lower }}

setup.py

@@ -12,7 +12,7 @@
     name="aqme",
     packages=find_packages(exclude=["tests"]),
     package_data={"aqme": ["templates/*"]},
-    version="1.4.2",
+    version="1.4.3",
     license="MIT",
     description="Automated Quantum Mechanical Environments",
     long_description="Automated Quantum Mechanical Environments",
@@ -29,7 +29,7 @@
         "automated",
     ],
     url="https://github.com/jvalegre/aqme",
-    download_url="https://github.com/jvalegre/aqme/archive/refs/tags/1.4.2.tar.gz",
+    download_url="https://github.com/jvalegre/aqme/archive/refs/tags/1.4.3.tar.gz",
     classifiers=[
         "Development Status :: 5 - Production/Stable",  # Chose either "3 - Alpha", "4 - Beta" or "5 - Production/Stable" as the current state of your package
         "Intended Audience :: Developers",  # Define that your audience are developers

tests/test_csearch.py

@@ -409,6 +409,38 @@ def test_csearch_fullmonte_parameters(
     assert mult == int(mols[0].GetProp("Mult"))
     os.chdir(w_dir_main)
 
+# tests for parameters of metals with double bonds
+@pytest.mark.parametrize(
+    "program, smi, name, charge",
+    [
+        ("rdkit", "N[SiH](N)[Cu]1CC1", "Cu_ethene", 0),
+        ("rdkit", "N[SiH](N)[Cu]1234C5C1C2C3C54", "Cu_Cp", -1),
+        ("rdkit", "N[SiH](N)[Cu]12345C6C1C2C3C4C65", "Cu_Ph", 0),
+    ],
+)
+def test_double_bond_chrg(
+    program,
+    smi,
+    name,
+    charge,
+):
+    os.chdir(csearch_methods_dir)
+    # runs the program with the different tests
+    csearch(
+        program=program,
+        smi=smi,
+        name=name,
+        sample=10,
+        metal_atoms=['Cu'],
+        metal_oxi=[1]
+    )
+
+    # tests here
+    file = str("CSEARCH/" + name + "_" + program + ".sdf")
+    mols = rdkit.Chem.SDMolSupplier(file, removeHs=False)
+    assert charge == int(mols[0].GetProp("Real charge"))
+    os.chdir(w_dir_main)
+
 
 # tests for parameters of csearch rdkit
 @pytest.mark.parametrize(
@@ -609,6 +641,25 @@ def test_csearch_rdkit_summ_parameters(
             False,
             4,
         ),
+        # metal atoms
+        (
+            "rdkit",
+            "I[Pd]([PH3+])(F)Cl",
+            "Pd_metal_only",
+            False,
+            True,
+            ["Pd"],
+            [2],
+            None,
+            None,
+            None,
+            None,
+            -1,
+            1,
+            None,
+            False,
+            1
+        ),
         # multiple templates
         (
             "rdkit",
@@ -786,6 +837,7 @@ def test_csearch_methods(
             smi=smi,
             name=name
         )
+
     elif not complex and not metal_complex:
         csearch(
             w_dir_main=csearch_methods_dir,
@@ -795,17 +847,29 @@ def test_csearch_methods(
         )
 
     elif metal_complex is True:
-        csearch(
-            w_dir_main=csearch_methods_dir,
-            program=program,
-            smi=smi,
-            name=name,
-            metal_atoms=metal,
-            metal_oxi=metal_oxi,
-            complex_type=complex_type,
-            mult=mult,
-            sample=10
-        )
+        if name == 'Pd_metal_only':
+            csearch(
+                w_dir_main=csearch_methods_dir,
+                program=program,
+                smi=smi,
+                name=name,
+                metal_atoms=metal,
+                metal_oxi=metal_oxi,
+                mult=mult,
+                sample=10
+            )
+        else:
+            csearch(
+                w_dir_main=csearch_methods_dir,
+                program=program,
+                smi=smi,
+                name=name,
+                metal_atoms=metal,
+                metal_oxi=metal_oxi,
+                complex_type=complex_type,
+                mult=mult,
+                sample=10
+            )
 
     elif complex is True:
         csearch(
@@ -821,7 +885,7 @@ def test_csearch_methods(
 
     if destination:
         file = str(csearch_methods_dir+"/Et_sdf_files/" + name + "_" + program + ".sdf")
-    elif metal_complex is False or name in ['Ag_complex_crest','Cu_trigonal']:
+    elif metal_complex is False or name in ['Ag_complex_crest','Cu_trigonal','Pd_metal_only']:
         file = str(csearch_methods_dir+"/CSEARCH/" + name + "_" + program + ".sdf")
     else:
         file = str(
@@ -845,6 +909,25 @@ def test_csearch_methods(
     mols = rdkit.Chem.SDMolSupplier(file, removeHs=False, sanitize=False)
     assert charge == int(mols[-1].GetProp("Real charge"))
     assert mult == int(mols[-1].GetProp("Mult"))
+
+    # check that the metal is added back to the RDKit mol objects
+    metal_found = False
+    if name in ['Pd_complex','Pd_metal_only']:
+        outfile = open(file, "r")
+        outlines_sdf = outfile.readlines()
+        outfile.close()
+        for line in outlines_sdf:
+            if 'Pd  0' in line:
+                metal_found = True
+        assert metal_found
+    if name == 'Ag_complex_crest':
+        outfile = open(file, "r")
+        outlines_sdf = outfile.readlines()
+        outfile.close()
+        for line in outlines_sdf:
+            if 'Ag  0' in line:
+                metal_found = True
+        assert metal_found
     if name == 'nci':
         assert len(mols) > 350
     # the n of conformers decreases when --nci is used