Fix incorrect rendering of radicals

Rdkit has a complicated system for including Hydrogen atoms. It
is discussed here:
https://sourceforge.net/p/rdkit/mailman/message/36699970/

In brief, there are three levels of hydrogen inclusion: Explicit
and in graph, explicit, and implicit. This is meant to model the
difference between hydrogens that are explicit in a SMILES string,
and those that are inferred by the toolkit, but it unfortunately
goes even further than that. For instance, if hydrogens are removed
from a molecule with the rdkit.Chem.rdmolops.RemoveHs() function,
which is the recommended way to make explicit hydrogens implicit
for visualization, RdKit can then imply the existence of hydrogens
using its valence model. This can, for example, convert a radical
to the stable species of equivalent charge, even if this means
adding hydrogens. In fact, before this commit, simply calling the
OFFTK Molecule.to_rdkit method could convert, for example, the
hydroxyl radical to water:

>>> from openff.toolkit.topology import Molecule
>>> from rdkit import Chem
>>>
>>> radical = Molecule.from_smiles("[O][H]")
>>> rd_radical = radical.to_rdkit()
>>> Chem.MolToSmiles(rd_radical)
'[H]O'

Visualizing rd_radical clearly demonstrates that this is water;
the oxygen atom has an implicit hydrogen. This can also be seen
by building the radical up from individual atoms and bonds.

OFFTK has a different philosophy: all the hydrogens are either
provided or inferred when the Molecule is created, and then they are
explicitly represented in the molecular graph forever. This commit
sets the `NoImplicit` property for all atoms in a rdkit molecule
created by to_rdkit to True. Since all the hydrogens in a Molecule
are already explicitly represented, this more faithfully represents
the desired molecule in the rdkit ecosystem. The molecular species
is even now stable over sanitization:

>>> radical = Molecule.from_smiles("[O][H]")
>>> rd_hydroxyl_rad = radical.to_rdkit()
>>> Chem.SanitizeMol(rd_hydroxyl_rad) # Modifies in place
>>> Chem.MolToSmiles(rd_hydroxyl_rad)
[H][O]

Visualization even includes the radical dot!

openff/toolkit/tests/test_molecule.py

@@ -3724,6 +3724,14 @@ def test_visualize_nglview(self):
         # Ensure an NGLView widget is returned
         assert isinstance(mol.visualize(backend="nglview"), nglview.NGLWidget)
 
+        # Providing other arguments is an error
+        with pytest.raises(ValueError):
+            mol.visualize(backend="nglview", width=100)
+        with pytest.raises(ValueError):
+            mol.visualize(backend="nglview", height=100)
+        with pytest.raises(ValueError):
+            mol.visualize(backend="nglview", show_all_hydrogens=True)
+
     @requires_pkg("ipython")
     @requires_openeye
     def test_visualize_openeye(self):

openff/toolkit/topology/molecule.py

@@ -6298,7 +6298,9 @@ def add_conformer(self, coordinates):
 
         return self._add_conformer(coordinates)
 
-    def visualize(self, backend="rdkit", width=500, height=300, show_hydrogens=True):
+    def visualize(
+        self, backend="rdkit", width=None, height=None, show_all_hydrogens=None
+    ):
         """
         Render a visualization of the molecule in Jupyter
 
@@ -6317,8 +6319,8 @@ def visualize(self, backend="rdkit", width=500, height=300, show_hydrogens=True)
         height : int, optional, default=300
             Width of the generated representation (only applicable to
             ``backend=openeye`` or ``backend=rdkit``)
-        show_hydrogens : bool, optional, default=True
-            Whether to explicitly depict hydrogen atoms (only applicable to
+        show_all_hydrogens : bool, optional, default=True
+            Whether to explicitly depict all hydrogen atoms. (only applicable to
             ``backend=openeye`` or ``backend=rdkit``)
 
         Returns
@@ -6340,6 +6342,16 @@ def visualize(self, backend="rdkit", width=500, height=300, show_hydrogens=True)
                 import nglview as nv
             except ImportError:
                 raise MissingDependencyError("nglview")
+
+            if (
+                width is not None
+                or height is not None
+                or show_all_hydrogens is not None
+            ):
+                raise ValueError(
+                    "The width, height, and show_hydrogens arguments do not apply to the nglview backend."
+                )
+
             if self.conformers:
                 from openff.toolkit.utils.viz import _OFFTrajectoryNGLView
 
@@ -6351,18 +6363,27 @@ def visualize(self, backend="rdkit", width=500, height=300, show_hydrogens=True)
                     "Visualizing with NGLview requires that the molecule has "
                     "conformers."
                 )
+
+        width = 500 if width is None else width
+        height = 300 if height is None else height
+        show_all_hydrogens = True if show_all_hydrogens is None else show_all_hydrogens
+
         if backend == "rdkit":
             if RDKIT_AVAILABLE:
                 from IPython.display import SVG
                 from rdkit.Chem.Draw import rdDepictor, rdMolDraw2D
                 from rdkit.Chem.rdmolops import RemoveHs
 
                 rdmol = self.to_rdkit()
-                if not show_hydrogens:
-                    rdmol = RemoveHs(rdmol)
+
+                if not show_all_hydrogens:
+                    # updateExplicitCount: Keep a record of the hydrogens we remove.
+                    # This is used in visualization to distinguish eg radicals from normal species
+                    rdmol = RemoveHs(rdmol, updateExplicitCount=True)
 
                 rdDepictor.SetPreferCoordGen(True)
                 rdDepictor.Compute2DCoords(rdmol)
+                rdmol = rdMolDraw2D.PrepareMolForDrawing(rdmol)
 
                 drawer = rdMolDraw2D.MolDraw2DSVG(width, height)
                 drawer.DrawMolecule(rdmol)
@@ -6387,7 +6408,7 @@ def visualize(self, backend="rdkit", width=500, height=300, show_hydrogens=True)
                     width, height, oedepict.OEScale_AutoScale
                 )
 
-                if show_hydrogens:
+                if show_all_hydrogens:
                     opts.SetHydrogenStyle(oedepict.OEHydrogenStyle_ImplicitAll)
 
                 oedepict.OEPrepareDepiction(oemol)

openff/toolkit/utils/rdkit_wrapper.py

@@ -1688,6 +1688,9 @@ def to_rdkit(cls, molecule, aromaticity_model=DEFAULT_AROMATICITY_MODEL):
             elif atom.stereochemistry == "R":
                 rdatom.SetChiralTag(Chem.CHI_TETRAHEDRAL_CCW)
 
+            # Stop rdkit from adding implicit hydrogens
+            rdatom.SetNoImplicit(True)
+
             rd_index = rdmol.AddAtom(rdatom)
 
             # Let's make sure al the atom indices in the two molecules