Fix redistribution of excess QM charge. [closes #294]

Author: lohedges

doc/source/changelog.rst

@@ -20,6 +20,7 @@ organisation on `GitHub <https://github.com/openbiosim/sire>`__.
 * 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.