Merge pull request #76 from molssi-seamm/dev

Rotated molecule from SMILES, InChI, or InChIKey to standard orientation

HISTORY.rst

@@ -1,6 +1,11 @@
 =======
 History
 =======
+2024.5.6 -- Rotated molecule from SMILES, InChI, or InChIKey to standard orientation
+    * Molecules created from line notation are created in an random orientation. This
+      enhancement rotates them to the standard orientation, which will look nice for
+      small, symmetric molecules.
+
 2024.5.5 -- Bugfix: bonds in RDKit
     * There was an indexing bug translating bonds back from RDKit to SEAMM. The famous
       0/1 problem!

molsystem/inchi.py

@@ -74,7 +74,7 @@ def to_inchi(self, key=False, openbabel=False):
 
         return inchi.strip()
 
-    def from_inchi(self, inchi, name=None, openbabel=True):
+    def from_inchi(self, inchi, name=None, reorient=True, openbabel=True):
         """Create the system from a InChI string.
 
         Parameters
@@ -83,6 +83,8 @@ def from_inchi(self, inchi, name=None, openbabel=True):
             The InChI string
         name : str = None
             The name of the molecule
+        reorient : bool = True
+            Whether to reorient to the standard orientation
         openbabel : bool = False
             Whether to use Openbabel rather than default of RDKit
 
@@ -116,15 +118,18 @@ def from_inchi(self, inchi, name=None, openbabel=True):
             AllChem.EmbedMolecule(mol)
             self.from_RDKMol(mol)
 
-            print(f"{self.n_atoms=}")
-            print(f"{self.atoms.symbols=}")
+        # Rotate to standard orientation
+        rdkMol = self.to_RDKMol()
+        rdkConf = rdkMol.GetConformers()[0]
+        Chem.rdMolTransforms.CanonicalizeConformer(rdkConf)
+        self.from_RDKMol(rdkMol)
 
         if name is not None:
             self.name = name
         else:
             self.name = save
 
-    def from_inchikey(self, inchikey, name=None):
+    def from_inchikey(self, inchikey, name=None, reorient=True):
         """Create the system from an InChIKey string.
 
         Parameters
@@ -133,33 +138,15 @@ def from_inchikey(self, inchikey, name=None):
             The InChIKey string
         name : str = None
             The name of the molecule
+        reorient : bool = True
+            Whether to reorient to the standard orientation
 
         Returns
         -------
         None
         """
         inchi = self._get_inchi(inchikey)
-
-        save = self.name
-
-        obConversion = OB.OBConversion()
-        obConversion.SetInAndOutFormats("inchi", "mdl")
-        mol = OB.OBMol()
-        obConversion.ReadString(mol, inchi)
-
-        # Add hydrogens
-        mol.AddHydrogens()
-
-        # Get coordinates for a 3-D structure
-        builder = OB.OBBuilder()
-        builder.Build(mol)
-
-        self.from_OBMol(mol)
-
-        if name is not None:
-            self.name = name
-        else:
-            self.name = save
+        self.from_inchi(inchi, reorient=reorient)
 
     def _get_inchi(self, inchikey):
         """Get the InChI from PubChem given the InChIKey."""

molsystem/rdkit_.py

@@ -56,6 +56,7 @@ def to_RDKMol(self):
             1: Chem.BondType.SINGLE,
             2: Chem.BondType.DOUBLE,
             3: Chem.BondType.TRIPLE,
+            5: Chem.BondType.AROMATIC,
         }
         for row in self.bonds.bonds():
             rdk_mol.AddBond(
@@ -108,6 +109,7 @@ def from_RDKMol(self, rdk_mol):
             Chem.BondType.SINGLE: 1,
             Chem.BondType.DOUBLE: 2,
             Chem.BondType.TRIPLE: 3,
+            Chem.BondType.AROMATIC: 5,
         }
         for rdk_bond in rdk_mol.GetBonds():
             i = rdk_bond.GetBeginAtom().GetIdx()

molsystem/smiles.py

@@ -94,7 +94,7 @@ def to_smiles(self, canonical=False, hydrogens=False, isomeric=True, rdkit=False
 
         return smiles.strip()
 
-    def from_smiles(self, smiles, name=None, rdkit=False):
+    def from_smiles(self, smiles, name=None, reorient=True, rdkit=False):
         """Create the system from a SMILES string.
 
         Parameters
@@ -103,6 +103,8 @@ def from_smiles(self, smiles, name=None, rdkit=False):
             The SMILES string
         name : str = None
             The name of the molecule
+        reorient : bool = True
+            Whether to reorient to the standard orientation
         rdkit : bool = False
             Whether to use RDKit rather than default of OpenBabel
 
@@ -136,6 +138,13 @@ def from_smiles(self, smiles, name=None, rdkit=False):
 
             self.from_OBMol(mol)
 
+        # Rotate to standard orientation
+        if reorient:
+            rdkMol = self.to_RDKMol()
+            rdkConf = rdkMol.GetConformers()[0]
+            Chem.rdMolTransforms.CanonicalizeConformer(rdkConf)
+            self.from_RDKMol(rdkMol)
+
         if name is not None:
             self.name = name
         else:

tests/test_angiotensin.py

@@ -118,7 +118,7 @@ def testdb():
     # Angiotensin II from SMILES
     system = db.create_system("angiotensin")
     configuration = system.create_configuration("SMILES")
-    configuration.from_smiles(SMILES)
+    configuration.from_smiles(SMILES, reorient=False)
 
     # And read in
     system.read_cif_file(data_path / "1n9v.cif")

tests/test_openbabel.py

@@ -458,7 +458,7 @@ def test_all_residue_search(configuration):
         "PHE_LL": [(64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75)],
     }
 
-    configuration.from_smiles(SMILES)
+    configuration.from_smiles(SMILES, reorient=False)
 
     for name, sc in sidechains.items():
         smarts = f"[NH1X3][C@@H]({sc})[CX3]=[OX1]"