Merge pull request #293 from OpenBioSim/fix_292

Fix issue #292

Author: lohedges

doc/source/changelog.rst

@@ -19,6 +19,8 @@ organisation on `GitHub <https://github.com/openbiosim/sire>`__.
 * Please add an item to this CHANGELOG for any new features or bug fixes when creating a PR.
 * Allow user to force fresh inference of stereochemistry when converting to RDKit format.
 * Fix setting of positive formal charge when reading SDF files.
+* Only use ``atom->setNoImplicit(True)`` inside custom RDKit sterochemistry inference function.
+* Fix redistribution of excess QM charge and make it the default behaviour.
 
 `2024.4.0 <https://github.com/openbiosim/sire/compare/2024.3.1...2024.4.0>`__ - Feb 2025
 ----------------------------------------------------------------------------------------

src/sire/qm/__init__.py

@@ -18,7 +18,7 @@ def create_engine(
     cutoff="7.5A",
     neighbour_list_frequency=0,
     mechanical_embedding=False,
-    redistribute_charge=False,
+    redistribute_charge=True,
     map=None,
 ):
     """
@@ -68,8 +68,11 @@ def create_engine(
 
     redistribute_charge : bool
         Whether to redistribute charge of the QM atoms to ensure that the total
-        charge of the QM region is an integer. Excess charge is redistributed
-        over the non QM atoms within the residues involved in the QM region.
+        charge of the QM region is an integer. If the QM region is a an entire
+        molecule, then the charge on each atom is shifted so that the total
+        charge is an integer. Alternatively, when the QM region is part of a
+        molecule, the excess charge is redistributed over the MM atoms within
+        the residues of the QM region.
 
     Returns
     -------

src/sire/qm/_emle.py

@@ -107,7 +107,7 @@ def emle(
     calculator,
     cutoff="7.5A",
     neighbour_list_frequency=0,
-    redistribute_charge=False,
+    redistribute_charge=True,
     map=None,
 ):
     """
@@ -137,8 +137,11 @@ def emle(
 
     redistribute_charge : bool
         Whether to redistribute charge of the QM atoms to ensure that the total
-        charge of the QM region is an integer. Excess charge is redistributed
-        over the non QM atoms within the residues involved in the QM region.
+        charge of the QM region is an integer. If the QM region is a an entire
+        molecule, then the charge on each atom is shifted so that the total
+        charge is an integer. Alternatively, when the QM region is part of a
+        molecule, the excess charge is redistributed over the MM atoms within
+        the residues of the QM region.
 
     Returns
     -------

src/sire/qm/_utils.py

@@ -102,6 +102,7 @@ def _check_charge(mols, qm_atoms, map, redistribute_charge=False, tol=1e-6):
     """
 
     import math as _math
+    import sys as _sys
 
     from sire.units import e_charge as _e_charge
 
@@ -127,8 +128,16 @@ def _check_charge(mols, qm_atoms, map, redistribute_charge=False, tol=1e-6):
             # Redistribute the charge over the QM atoms.
             qm_frac = excess_charge / len(qm_atoms)
 
-            # Redistribute the charge over the non QM atoms.
-            rem_frac = excess_charge / (residues.num_atoms() - len(qm_atoms))
+            # Work out the number of non-QM atoms.
+            num_mm = residues.num_atoms() - len(qm_atoms)
+
+            # Redistribute the charge over the MM atoms.
+            if num_mm > 0:
+                rem_frac = excess_charge / num_mm
+
+            print(
+                f"Redistributing excess QM charge of {excess_charge}", file=_sys.stderr
+            )
 
             # Loop over the residues.
             for res in residues:
@@ -140,10 +149,11 @@ def _check_charge(mols, qm_atoms, map, redistribute_charge=False, tol=1e-6):
 
                 # Loop over the atoms in the residue.
                 for atom in res:
-                    # Shift the charge.
+                    # Shift the charge of the QM atoms.
                     if atom in qm_atoms:
                         cursor[atom][charge_prop] -= qm_frac
-                    else:
+                    # Shift the charge of the MM atoms.
+                    elif num_mm > 0:
                         cursor[atom][charge_prop] += rem_frac
 
                 # Commit the changes.

tests/convert/test_rdkit.py

@@ -122,16 +122,11 @@ def test_rdkit_infer_bonds(ejm55_sdf, ejm55_gro):
     sdf = ejm55_sdf[0].molecule()
     gro = ejm55_gro["not (protein or water)"].molecule()
 
-    from rdkit import Chem
-
     assert sdf.smiles() == gro.smiles()
 
     match_sdf = sdf["smarts [NX3][CX3](=[OX1])[#6]"]
     match_gro = gro["smarts [NX3][CX3](=[OX1])[#6]"]
 
-    print(match_sdf)
-    print(match_gro)
-
     assert len(match_sdf) == 1
     assert len(match_gro) == 1
 
@@ -167,3 +162,20 @@ def test_rdkit_preserve_info(ala_mols, ejm55_gro):
 
             assert res0.name() == res1.name()
             assert res0.number() == res1.number()
+
+
+@pytest.mark.skipif(
+    "rdkit" not in sr.convert.supported_formats(),
+    reason="rdkit support is not available",
+)
+def test_rdkit_force_infer():
+    mol = sr.load_test_files("missing_cyanide_bond.sdf")[0]
+
+    rdmol = sr.convert.to(mol, "rdkit", map={"force_stereo_inference": False})
+    rdmol_infer = sr.convert.to(mol, "rdkit", map={"force_stereo_inference": True})
+
+    bond = rdmol.GetBonds()[0].GetBondType().name
+    bond_infer = rdmol_infer.GetBonds()[0].GetBondType().name
+
+    assert bond == "SINGLE"
+    assert bond_infer == "TRIPLE"

wrapper/Convert/SireRDKit/sire_rdkit.cpp

@@ -411,6 +411,7 @@ namespace SireRDKit
         {
             if (atom->getAtomicNum() > 1)
             {
+                atom->setNoImplicit(true);
                 atoms.append(std::make_pair(get_nb_unpaired_electrons(*atom),
                                             atom));
             }
@@ -646,9 +647,6 @@ namespace SireRDKit
 
             a->setAtomicNum(element.nProtons());
 
-            // don't automatically add hydrogens
-            a->setNoImplicit(true);
-
             elements.append(element);
 
             try
@@ -791,74 +789,11 @@ namespace SireRDKit
             infer_bond_info(molecule);
         }
 
-        // try each sanitisation step in turn, skipping failed
-        try
-        {
-            RDKit::MolOps::cleanUp(molecule);
-        }
-        catch (...)
-        {
-        }
-
-        try
-        {
-            molecule.updatePropertyCache();
-        }
-        catch (...)
-        {
-        }
-
-        try
-        {
-            RDKit::MolOps::symmetrizeSSSR(molecule);
-        }
-        catch (...)
-        {
-        }
-
-        try
-        {
-            RDKit::MolOps::Kekulize(molecule);
-        }
-        catch (...)
-        {
-        }
-
-        try
-        {
-            RDKit::MolOps::assignRadicals(molecule);
-        }
-        catch (...)
-        {
-        }
-
-        try
-        {
-            RDKit::MolOps::setAromaticity(molecule);
-        }
-        catch (...)
-        {
-        }
-
-        try
-        {
-            RDKit::MolOps::setConjugation(molecule);
-        }
-        catch (...)
-        {
-        }
-
+        // sanitze the molecule.
         try
         {
-            RDKit::MolOps::setHybridization(molecule);
-        }
-        catch (...)
-        {
-        }
+            RDKit::MolOps::sanitizeMol(molecule);
 
-        try
-        {
-            RDKit::MolOps::cleanupChirality(molecule);
         }
         catch (...)
         {