Merge pull request #190 from pierrelb/checkQMGeom

Check QMTP optimized molecular geometries are the desired species.

It is possible that the geometry optimization "reacted" the species to a different isomer.
This checks that the final geometry is isomorphic with the desired molecule.
Bonds are guessed based on distances between atoms. Bond orders (single, double, etc) are not checked,
only connectivity.

Author: rwest

rmgpy/molecule/element.pxd

@@ -30,5 +30,6 @@ cdef class Element:
     cdef public str name
     cdef public str symbol
     cdef public float mass
+    cdef public float covRadius
 
 cpdef Element getElement(value)

rmgpy/molecule/element.py

@@ -41,6 +41,7 @@
 """
 
 import cython
+from rdkit.Chem import GetPeriodicTable
 
 ################################################################################
 
@@ -53,7 +54,7 @@ class ElementError(Exception):
     pass
 
 ################################################################################
-
+_rdkit_periodic_table = GetPeriodicTable()
 class Element:
     """
     A chemical element. The attributes are:
@@ -65,6 +66,7 @@ class Element:
     `symbol`    ``str``         The symbol used for the element
     `name`      ``str``         The IUPAC name of the element
     `mass`      ``float``       The mass of the element in kg/mol
+    `covRadius` ``float``       Covalent bond radius in Angstrom
     =========== =============== ================================================
     
     This class is specifically for properties that all atoms of the same element
@@ -76,6 +78,7 @@ def __init__(self, number, symbol, name, mass):
         self.symbol = intern(symbol)
         self.name = name
         self.mass = mass
+        self.covRadius = _rdkit_periodic_table.GetRcovalent(symbol)
 
     def __str__(self):
         """

rmgpy/molecule/molecule.pxd

@@ -29,6 +29,7 @@ from .atomtype cimport AtomType
 from .group cimport GroupAtom, GroupBond, Group
 from .element cimport Element
 cimport rmgpy.constants as constants
+cimport numpy
 
 ################################################################################
 cdef public dict _known_smiles
@@ -41,7 +42,7 @@ cdef class Atom(Vertex):
     cdef public short charge
     cdef public str label
     cdef public AtomType atomType
-    cdef public list coords
+    cdef public numpy.ndarray coords
     cdef public short lonePairs
 
     cpdef bint equivalent(self, Vertex other) except -2
@@ -170,6 +171,8 @@ cdef class Molecule(Graph):
 
     cpdef fromAdjacencyList(self, str adjlist, bint saturateH=?)
 
+    cpdef fromXYZ(self, numpy.ndarray atomicNums, numpy.ndarray coordinates)
+    
     cpdef str toInChI(self)
 
     cpdef str toAugmentedInChI(self)

rmgpy/molecule/molecule.py

@@ -98,7 +98,7 @@ class Atom(Vertex):
     `spinMultiplicity`  ``short``           The spin multiplicity of the atom
     `charge`            ``short``           The formal charge of the atom
     `label`             ``str``             A string label that can be used to tag individual atoms
-    `coords`
+    `coords`            ``numpy array``     The (x,y,z) coordinates in Angstrom
     `lonePairs`         ``short``           The number of lone electron pairs
     =================== =================== ====================================
 
@@ -107,7 +107,7 @@ class Atom(Vertex):
     e.g. ``atom.symbol`` instead of ``atom.element.symbol``.
     """
 
-    def __init__(self, element=None, radicalElectrons=0, spinMultiplicity=1, charge=0, label='', lonePairs=0):
+    def __init__(self, element=None, radicalElectrons=0, spinMultiplicity=1, charge=0, label='', lonePairs=0, coords=None):
         Vertex.__init__(self)
         if isinstance(element, str):
             self.element = elements.__dict__[element]
@@ -119,7 +119,7 @@ def __init__(self, element=None, radicalElectrons=0, spinMultiplicity=1, charge=
         self.label = label
         self.atomType = None
         self.lonePairs = lonePairs
-        self.coords = list()
+        self.coords = coords
 
     def __str__(self):
         """
@@ -829,6 +829,51 @@ def deleteHydrogens(self):
         for atom in hydrogens:
             self.removeAtom(atom)
 
+    def connectTheDots(self):
+        """
+        Delete all bonds, and set them again based on the Atoms' coords.
+        Does not detect bond type.
+        """
+        cython.declare(criticalDistance=float, i=int, atom1=Atom, atom2=Atom,
+                       bond=Bond, atoms=list, zBoundary=float)
+                       # groupBond=GroupBond, 
+        self._fingerprint = None
+        
+        atoms = self.vertices
+        
+        # Ensure there are coordinates to work with
+        for atom in atoms:
+            assert atom.coords != None
+        
+        # If there are any bonds, remove them
+        for atom1 in atoms:
+            for bond in self.getBonds(atom1):
+                self.removeEdge(bond)
+        
+        # Sort atoms by distance on the z-axis
+        sortedAtoms = sorted(atoms, key=lambda x: x.coords[2])
+        
+        for i, atom1 in enumerate(sortedAtoms):
+            for atom2 in sortedAtoms[i+1:]:
+                # Set upper limit for bond distance
+                criticalDistance = (atom1.element.covRadius + atom2.element.covRadius + 0.45)**2
+                
+                # First atom that is more than 4.0 Anstroms away in the z-axis, break the loop
+                # Atoms are sorted along the z-axis, so all following atoms should be even further
+                zBoundary = (atom1.coords[2] - atom2.coords[2])**2
+                if zBoundary > 16.0:
+                    break
+                
+                distanceSquared = sum((atom1.coords - atom2.coords)**2)
+                
+                if distanceSquared > criticalDistance or distanceSquared < 0.40:
+                    continue
+                else:
+                    # groupBond = GroupBond(atom1, atom2, ['S','D','T','B'])
+                    bond = Bond(atom1, atom2, 'S')
+                    self.addBond(bond)
+        self.updateAtomTypes()
+        
     def updateAtomTypes(self):
         """
         Iterate through the atoms in the structure, checking their atom types
@@ -1153,6 +1198,44 @@ def fromAdjacencyList(self, adjlist, saturateH=False):
         self.updateAtomTypes()
 
         return self
+        
+    def fromXYZ(self, atomicNums, coordinates):
+        """
+        Create an RMG molecule from a list of coordinates and a corresponding
+        list of atomic numbers. These are typically received from CCLib and the
+        molecule is sent to `ConnectTheDots` so will only contain single bonds.
+        """
+        
+        _rdkit_periodic_table = elements.GetPeriodicTable()
+        
+        for i, atNum in enumerate(atomicNums):
+            atom = Atom(_rdkit_periodic_table.GetElementSymbol(int(atNum)))
+            atom.coords = coordinates[i]
+            self.addAtom(atom)
+        return self.connectTheDots()
+    
+    def toSingleBonds(self):
+        """
+        Returns a copy of the current molecule, consisting of only single bonds.
+        
+        This is useful for isomorphism comparison against something that was made
+        via fromXYZ, which does not attempt to perceive bond orders
+        """
+        cython.declare(atom1=Atom, atom2=Atom, bond=Bond, newMol=Molecule, atoms=list, mapping=dict)
+    
+        newMol = Molecule()
+        atoms = self.atoms
+        mapping = {}
+        for atom1 in atoms:
+            atom2 = newMol.addAtom(Atom(atom1.element))
+            mapping[atom1] = atom2
+    
+        for atom1 in atoms:
+            for atom2 in atom1.bonds:
+                bond = Bond(mapping[atom1], mapping[atom2], 'S')
+                newMol.addBond(bond)
+        newMol.updateAtomTypes()
+        return newMol
 
     def toInChI(self):
         """

rmgpy/qm/gaussian.py

@@ -5,6 +5,7 @@
 from subprocess import Popen, PIPE
 import re
 
+from rmgpy.molecule import Molecule
 from qmdata import CCLibData
 from molecule import QMMolecule
 
@@ -116,9 +117,19 @@ def verifyOutputFile(self):
             return False
 
         if InChIMatch:
-            logging.info("Successful Gaussian quantum result found in {0}".format(self.outputFilePath))
-            # " + self.molfile.name + " ("+self.molfile.InChIAug+") has been found. This log file will be used.")
-            return True
+            # Compare the optimized geometry to the original molecule
+            qmData = self.parse()
+            
+            cclibMol = Molecule()
+            cclibMol.fromXYZ(qmData.atomicNumbers, qmData.atomCoords.value)
+            testMol = self.molecule.toSingleBonds()
+            
+            if cclibMol.isIsomorphic(testMol):
+                logging.info("Successful MOPAC quantum result found in {0}".format(self.outputFilePath))
+                return True
+            else:
+                logging.info("Incorrect connectivity for optimized geometry in file {0}".format(self.outputFilePath))
+                return False
         else:
             return False # until the next line works
 

rmgpy/qm/molecule.py

@@ -68,7 +68,7 @@ def generateRDKitGeometries(self, boundsMatrix=None):
             distGeomAttempts = 5*(atoms-3) #number of conformer attempts is just a linear scaling with molecule size, due to time considerations in practice, it is probably more like 3^(n-3) or something like that
 
         rdmol, minEid = self.rd_embed(rdmol, distGeomAttempts, boundsMatrix)
-        self.save_coordinates(rdmol, minEid, rdAtIdx)
+        self.saveCoordinatesFromRDMol(rdmol, minEid, rdAtIdx)
 
     def rd_build(self):
         """
@@ -104,11 +104,17 @@ def rd_embed(self, rdmol, numConfAttempts, boundsMatrix=None):
 
         return rdmol, minEid
 
-    def save_coordinates(self, rdmol, minEid, rdAtIdx):
+    def saveCoordinatesFromRDMol(self, rdmol, minEid, rdAtIdx):
         # Save xyz coordinates on each atom in molecule ****
         for atom in self.molecule.atoms:
             point = rdmol.GetConformer(minEid).GetAtomPosition(atom.sortingLabel)
-            atom.coords = [point.x, point.y, point.z]
+            atom.coords = numpy.array([point.x, point.y, point.z])
+    
+    def saveCoordinatesFromQMData(self, qmdata):
+        """
+        Save geometry info from QMData (eg CCLibData)
+        """
+        raise NotImplementedError
 
 def loadThermoDataFile(filePath):
     """

rmgpy/qm/mopac.py

@@ -4,6 +4,7 @@
 import logging
 from subprocess import Popen, PIPE
 
+from rmgpy.molecule import Molecule
 from qmdata import CCLibData
 from molecule import QMMolecule
 
@@ -127,9 +128,19 @@ def verifyOutputFile(self):
             return False
 
         if InChIMatch:
-            logging.info("Successful MOPAC quantum result found in {0}".format(self.outputFilePath))
-            # " + self.molfile.name + " ("+self.molfile.InChIAug+") has been found. This log file will be used.")
-            return True
+            # Compare the optimized geometry to the original molecule
+            qmData = self.parse()
+            
+            cclibMol = Molecule()
+            cclibMol.fromXYZ(qmData.atomicNumbers, qmData.atomCoords.value)
+            testMol = self.molecule.toSingleBonds()
+            
+            if cclibMol.isIsomorphic(testMol):
+                logging.info("Successful MOPAC quantum result found in {0}".format(self.outputFilePath))
+                return True
+            else:
+                logging.info("Incorrect connectivity for optimized geometry in file {0}".format(self.outputFilePath))
+                return False
 
         #InChIs do not match (most likely due to limited name length mirrored in log file (240 characters), but possibly due to a collision)
         return self.checkForInChiKeyCollision(logFileInChI) # Not yet implemented!