Repository: openforcefield/openff-toolkit

Issue:
Molecule.from_rdkit() allows input of implicit-hydrogen molecules
**Describe the bug**
The `Molecule.from_rdkit(rdmol)` method will allow implicit-hydrogen molecules to be imported without converting them to explicit-hydrogen molecules.

**To Reproduce**
```python
from rdkit import Chem
molOrig = Chem.MolFromSmiles("CCCCCCCCCC")
ofmol = Molecule.from_rdkit(molOrig)
print(ofmol.atoms)
from rdkit import Chem
molOrig = Chem.MolFromSmiles("CCCCCCCCCC")
ofmol = Molecule.from_rdkit(molOrig)
print(ofmol.atoms)
```
returns a molecule without explicit protons that works in our pipeline just fine, except will always fail to produce anything useful:
**Output**
```
[Atom(name=, atomic number=6), Atom(name=, atomic number=6), Atom(name=, atomic number=6), Atom(name=, atomic number=6), Atom(name=, atomic number=6), Atom(name=, atomic number=6), Atom(name=, atomic number=6), Atom(name=, atomic number=6), Atom(name=, atomic number=6), Atom(name=, atomic number=6)]
```

**Computing environment (please complete the following information):**
 - `osx`
 - `openforcefield==0.6.0` and `rdkit==2019.09.03`

**Additional Context**
cc: https://github.com/openmm/openmmforcefields/issues/107

Remedies are either
* throw an exception on import
* add hydrogens explicitly
* don't worry about anything until `Molecule._validate()` is called just before we need a valid molecule (most permissive, hardest to debug issues)

Discussion:
mattwthompson: The last commit involved manually fixing several merge conflicts, since this has lingered for months. Passing locally and the diffs against `upstream/master` look the same before and after, so I think it went smoothly.

Do note that **some** tests are still going to fail while #889 is unresolved, but they are not necessarily related to this patch.
lilyminium: I ran into something similar when I was trying to add caps to residues. I'm not sure why this is default, but RDKit doesn't think you have implicit hydrogens unless you set atom.SetNoImplicit(False) (i.e., allow implicit hydrogens). 

When I add this to your `from_smiles()` [function](https://github.com/openforcefield/openff-toolkit/blob/ae2a812a4be4184d9b5db71e370341346a3c0343/openff/toolkit/utils/toolkits.py#L3355-L3360), @j-wags's tests give me the desired 8 atoms. 

```diff
        for atom in rdmol.GetAtoms():
+++            atom.SetNoImplicit(False)
            if atom.GetAtomMapNum() != 0:
                # set the map back to zero but hide the index in the atom prop data
                atom.SetProp("_map_idx", str(atom.GetAtomMapNum()))
                # set it back to zero
                atom.SetAtomMapNum(0)
+++        rdmol.UpdatePropertyCache(strict=False)     
```

> the RDKit backend silently places them at (0,0,0)

I don't know how to query positions from OFF Molecules, but does this fix that?

```diff
--- rdmol = Chem.AddHs(rdmol)
+++ rdmol = Chem.AddHs(rdmol, addCoords=True)
```

If you run into complaints about not calling calcImplicitValence, try `rdmol.UpdatePropertyCache(strict=False)`. Hydrogens are added in both `from_smiles` and `from_rdkit`.

Incidentally you might want to update [this line](https://github.com/openforcefield/openff-toolkit/blob/ae2a812a4be4184d9b5db71e370341346a3c0343/openff/toolkit/utils/toolkits.py#L3395-L3397):

```diff
		molecule = self.from_rdkit(
            rdmol, _cls=_cls, allow_undefined_stereo=allow_undefined_stereo,
+++			hydrogens_are_explicit=hydrogens_are_explicit,
        )
```
j-wags: @lilyminium Your changes above are great. I've put in the `addCoords` kwarg and tested it. The resulting behavior seems good. 

The other suggested changes gave the right behavior for `from_smiles` when given partial atom maps, but then I realized that OpenEyeToolkitWrapper has the same problem. I don't want to expand the scope for this PR to also modify `from_openeye`, so I'll open a new issue for those. 

Patch:
diff --git a/docs/environment.yml b/docs/environment.yml
index 3c2ef1d9..470cc184 100644
--- a/docs/environment.yml
+++ b/docs/environment.yml
@@ -9,6 +9,7 @@ dependencies:
     - sphinx_bootstrap_theme
     - myst-parser
     - testpath ==0.3.1
+    - docutils ==0.16 # Temporary https://github.com/executablebooks/MyST-Parser/issues/343
     # conda build dependencies
     - python
     - setuptools
diff --git a/docs/releasehistory.md b/docs/releasehistory.md
index ba0a2882..fac25475 100644
--- a/docs/releasehistory.md
+++ b/docs/releasehistory.md
@@ -8,8 +8,16 @@ Releases follow the `major.minor.micro` scheme recommended by [PEP440](https://w
 
 ## Changes since last release
 
-### Behaviors changed
-
+### New features and behaviors changed
+
+- [PR #762](https://github.com/openforcefield/openforcefield/pull/762): `Molecule.from_rdkit` now converts
+  implicit hydrogens into explicit hydrogens by default. This change may affect 
+  `RDKitToolkitWrapper/Molecule.from_smiles`, 
+  `from_mapped_smiles`, `from_file`, `from_file_obj`, `from_inchi`, and `from_qcschema`. 
+  This new behavior can be disabled using the
+  `hydrogens_are_explicit=True` keyword argument to `from_smiles`, or loading the molecule into
+  the desired protonation state in RDKit, and calling `from_rdkit` on the RDKit molecule with 
+  `hydrogens_are_explicit=True`.
 - [PR #894](https://github.com/openforcefield/openforcefield/pull/894): Calls to `Molecule.from_openeye`, 
   `Molecule.from_rdkit`, `Molecule.from_smiles`, `OpenEyeToolkitWrapper.from_smiles`, and 
   `RDKitToolkitWrapper.from_smiles` will now load atom maps into the the resulting 
diff --git a/openff/toolkit/topology/molecule.py b/openff/toolkit/topology/molecule.py
index 16a88068..62a57aa1 100644
--- a/openff/toolkit/topology/molecule.py
+++ b/openff/toolkit/topology/molecule.py
@@ -4495,7 +4495,9 @@ class FrozenMolecule(Serializable):
 
     @classmethod
     @RDKitToolkitWrapper.requires_toolkit()
-    def from_rdkit(cls, rdmol, allow_undefined_stereo=False):
+    def from_rdkit(
+        cls, rdmol, allow_undefined_stereo=False, hydrogens_are_explicit=False
+    ):
         """
         Create a Molecule from an RDKit molecule.
 
@@ -4507,6 +4509,8 @@ class FrozenMolecule(Serializable):
             An RDKit molecule
         allow_undefined_stereo : bool, default=False
             If false, raises an exception if oemol contains undefined stereochemistry.
+        hydrogens_are_explicit : bool, default=False
+            If False, RDKit will perform hydrogen addition using Chem.AddHs
 
         Returns
         -------
@@ -4528,6 +4532,7 @@ class FrozenMolecule(Serializable):
         molecule = toolkit.from_rdkit(
             rdmol,
             allow_undefined_stereo=allow_undefined_stereo,
+            hydrogens_are_explicit=hydrogens_are_explicit,
             _cls=cls,
         )
         return molecule
diff --git a/openff/toolkit/utils/toolkits.py b/openff/toolkit/utils/toolkits.py
index 2deb22bc..d2d8525e 100644
--- a/openff/toolkit/utils/toolkits.py
+++ b/openff/toolkit/utils/toolkits.py
@@ -2885,6 +2885,7 @@ class RDKitToolkitWrapper(ToolkitWrapper):
         pdbmol = self.from_rdkit(
             Chem.MolFromPDBFile(file_path, removeHs=False),
             allow_undefined_stereo=True,
+            hydrogens_are_explicit=True,
             _cls=_cls,
         )
 
@@ -3393,7 +3394,10 @@ class RDKitToolkitWrapper(ToolkitWrapper):
                     )
 
         molecule = self.from_rdkit(
-            rdmol, _cls=_cls, allow_undefined_stereo=allow_undefined_stereo
+            rdmol,
+            _cls=_cls,
+            allow_undefined_stereo=allow_undefined_stereo,
+            hydrogens_are_explicit=hydrogens_are_explicit,
         )
 
         return molecule
@@ -3949,7 +3953,13 @@ class RDKitToolkitWrapper(ToolkitWrapper):
 
         molecule._conformers = diverse_conformers
 
-    def from_rdkit(self, rdmol, allow_undefined_stereo=False, _cls=None):
+    def from_rdkit(
+        self,
+        rdmol,
+        allow_undefined_stereo=False,
+        hydrogens_are_explicit=False,
+        _cls=None,
+    ):
         """
         Create a Molecule from an RDKit molecule.
 
@@ -3963,6 +3973,8 @@ class RDKitToolkitWrapper(ToolkitWrapper):
             An RDKit molecule
         allow_undefined_stereo : bool, default=False
             If false, raises an exception if rdmol contains undefined stereochemistry.
+        hydrogens_are_explicit : bool, default=False
+            If False, RDKit will perform hydrogen addition using Chem.AddHs
         _cls : class
             Molecule constructor
 
@@ -3994,6 +4006,9 @@ class RDKitToolkitWrapper(ToolkitWrapper):
         # Make a copy of the RDKit Mol as we'll need to change it (e.g. assign stereo).
         rdmol = Chem.Mol(rdmol)
 
+        if not hydrogens_are_explicit:
+            rdmol = Chem.AddHs(rdmol, addCoords=True)
+
         # Sanitizing the molecule. We handle aromaticity and chirality manually.
         # This SanitizeMol(...) calls cleanUp, updatePropertyCache, symmetrizeSSSR,
         # assignRadicals, setConjugation, and setHybridization.


Test patch:
diff --git a/openff/toolkit/tests/test_toolkits.py b/openff/toolkit/tests/test_toolkits.py
index c78d8a9d..5ad0f166 100644
--- a/openff/toolkit/tests/test_toolkits.py
+++ b/openff/toolkit/tests/test_toolkits.py
@@ -198,20 +198,6 @@ rdkit_inchi_stereochemistry_lost = [
     "DrugBank_6865",
 ]
 
-rdkit_inchi_isomorphic_fails = [
-    "DrugBank_178",
-    "DrugBank_246",
-    "DrugBank_5847",
-    "DrugBank_700",
-    "DrugBank_1564",
-    "DrugBank_1700",
-    "DrugBank_4662",
-    "DrugBank_2052",
-    "DrugBank_2077",
-    "DrugBank_2082",
-    "DrugBank_2210",
-    "DrugBank_2642",
-]
 rdkit_inchi_roundtrip_mangled = ["DrugBank_2684"]
 
 openeye_iupac_bad_stereo = [
@@ -1955,17 +1941,9 @@ class TestRDKitToolkitWrapper:
 
             # compare the full molecule excluding the properties dictionary
             # turn of the bond order matching as this could move in the aromatic rings
-            if molecule.name in rdkit_inchi_isomorphic_fails:
-                # Some molecules graphs change during the round trip testing
-                # we test quite strict isomorphism here
-                with pytest.raises(AssertionError):
-                    assert molecule.is_isomorphic_with(
-                        mol2, bond_order_matching=False, toolkit_registry=toolkit
-                    )
-            else:
-                assert molecule.is_isomorphic_with(
-                    mol2, bond_order_matching=False, toolkit_registry=toolkit
-                )
+            assert molecule.is_isomorphic_with(
+                mol2, bond_order_matching=False, toolkit_registry=toolkit
+            )
 
     def test_smiles_charged(self):
         """Test RDKitWrapper functions for reading/writing charged SMILES"""
@@ -2135,6 +2113,21 @@ class TestRDKitToolkitWrapper:
             == expected_output_smiles
         )
 
+    def test_from_rdkit_implicit_hydrogens(self):
+        """
+        Test that hydrogens are inferred from hydrogen-less RDKit molecules,
+        unless the option is turned off.
+        """
+        from rdkit import Chem
+
+        rdmol = Chem.MolFromSmiles("CC")
+        offmol = Molecule.from_rdkit(rdmol)
+
+        assert any([a.atomic_number == 1 for a in offmol.atoms])
+
+        offmol_no_h = Molecule.from_rdkit(rdmol, hydrogens_are_explicit=True)
+        assert not any([a.atomic_number == 1 for a in offmol_no_h.atoms])
+
     @pytest.mark.parametrize(
         "smiles, expected_map", [("[Cl:1][Cl]", {0: 1}), ("[Cl:1][Cl:2]", {0: 1, 1: 2})]
     )