Fix distance geometry stereo issues (#2158)

* MDLFormat: Set tetrahedral stereo to unspecified for atom stereo parity
3. This was already the default behaviour for cis/trans stereochemistry.

* OBStereoFacade: Add functions to get all stereo objects of a certain type.

* Gen3DStereoHelper: Helper class for 3D coord generation

* OBDistanceGeometry:

* Fix handling of implicit stereo references
* Fix Calculate13Angle to work for triple bonds
* Use new Gen3DStereoHelper class

* OpGen3D: Fix stereochemistry issue (check after energy minimization)

* CisTransFrom3D:

* Fix issue with implicit stereo refs (L2236+L2260)
* Use torsion angles (fixes issue with non-planar geometry)

* Add and enable more distance geometry roundtrip test cases.
(Some disabled to ensure the test runs fast enough)

* CisTransFrom3D: Change log messages from warning to info

* OBDistanceGeometry/OpGen3D:

- Increase number of iterations
- Improve error handling

include/openbabel/distgeom.h

@@ -84,6 +84,13 @@ namespace OpenBabel {
     void AddConformer();
     void GetConformers(OBMol &mol);
 
+    /**
+     * Check if last call to AddConformer was successful.
+     *
+     * \return Success or failure.
+     */
+    bool WasSuccessful() const;
+
     /**
      * Convenience method to set up this molecule, generate a geometry and return it
      *

include/openbabel/stereo/stereo.h

@@ -402,6 +402,10 @@ namespace OpenBabel {
        * Get the number of tetrahedral stereocenters.
        */
       unsigned int NumTetrahedralStereo();
+      /**
+       * Get all the OBTetrahedralStereo objects.
+       */
+      std::vector<OBTetrahedralStereo*> GetAllTetrahedralStereo();
       /**
        * Check if atom with @p id is a tetrahedral center.
        * @return True if the atom with @p id has tetrahedral stereochemistry.
@@ -421,6 +425,10 @@ namespace OpenBabel {
        * Get the number of cis/trans stereocenters.
        */
       unsigned int NumCisTransStereo();
+      /**
+       * Get all the OBCisTransStereo objects.
+       */
+      std::vector<OBCisTransStereo*> GetAllCisTransStereo();
       /**
        * Check if bond with @p id is a stereogenic cis/trans double bond.
        * @return True if the bond with @p id has cis/trans stereochemistry.
@@ -440,6 +448,10 @@ namespace OpenBabel {
        * Get the number of square-planar stereocenters.
        */
       unsigned int NumSquarePlanarStereo();
+      /**
+       * Get all the OBSquarePlanarStereo objects.
+       */
+      std::vector<OBSquarePlanarStereo*> GetAllSquarePlanarStereo();
       /**
        * Check if atom with @p id is a stereogenic square-planar atom.
        * @return True if the atom with @p id has square-planar stereochemistry.

src/CMakeLists.txt

@@ -106,6 +106,7 @@ set(stereo_srcs
   stereo/tetrahedral.cpp
   stereo/perception.cpp
   stereo/facade.cpp
+  stereo/gen3dstereohelper.cpp
 )
 
 set(openbabel_srcs

src/distgeom.cpp

@@ -30,6 +30,7 @@ GNU General Public License for more details.
 #include <openbabel/generic.h>
 #include "rand.h"
 #include <LBFGS.h>
+#include "stereo/gen3dstereohelper.h"
 
 #include <openbabel/stereo/stereo.h>
 #include <openbabel/stereo/cistrans.h>
@@ -62,6 +63,7 @@ namespace OpenBabel {
       bounds = Eigen::MatrixXf(static_cast<int>(N), static_cast<int>(N));
       preMet = Eigen::MatrixXf(bounds);
       debug = false;
+      success = false;
     }
     ~DistanceGeometryPrivate()
     { }
@@ -102,7 +104,10 @@ namespace OpenBabel {
     }
 
     Eigen::MatrixXf bounds, preMet;
-    bool debug; double maxBoxSize; };
+    bool debug; double maxBoxSize;
+    OBGen3DStereoHelper stereoHelper;
+    bool success;
+  };
 
 
   OBDistanceGeometry::OBDistanceGeometry(): _d(NULL) {}
@@ -133,13 +138,12 @@ namespace OpenBabel {
     dim = 4;
     _mol = mol;
 
-    OBConversion conv;
-    conv.SetOutFormat("can");
-    input_smiles = conv.WriteString(&_mol, true);
+    _d = new DistanceGeometryPrivate(mol.NumAtoms());
+
+    _d->stereoHelper.Setup(&_mol);
 
     _mol.SetDimension(3);
     _vdata = _mol.GetAllData(OBGenericDataType::StereoData);
-    _d = new DistanceGeometryPrivate(mol.NumAtoms());
 
     SetUpperBounds();
     // Do we use the current geometry for default 1-2 and 1-3 bounds?
@@ -180,15 +184,18 @@ namespace OpenBabel {
         OBTetrahedralStereo::Config config = ts->GetConfig();
         vector<unsigned long> nbrs;
 
-        nbrs.push_back(_mol.GetAtomById(config.from)->GetIdx()-1);
-        for(size_t i=0; i<config.refs.size(); i++) {
-          nbrs.push_back(_mol.GetAtomById(config.refs[i])->GetIdx()-1);
-        }
-
-        // This is required to avoid segfault (why?)
         unsigned long centerIdx = _mol.GetAtomById(config.center)->GetIdx()-1;
-        for(size_t i=0; i<nbrs.size(); i++) {
-          if (nbrs[i] > _mol.NumAtoms()) nbrs[i] = centerIdx;
+
+        if (config.from == OBStereo::ImplicitRef)
+          nbrs.push_back(centerIdx);
+        else
+          nbrs.push_back(_mol.GetAtomById(config.from)->GetIdx()-1);
+
+        for(size_t i=0; i<config.refs.size(); i++) {
+          if (config.refs[i] == OBStereo::ImplicitRef)
+            nbrs.push_back(centerIdx);
+          else
+            nbrs.push_back(_mol.GetAtomById(config.refs[i])->GetIdx()-1);
         }
 
         if(config.winding == OBStereo::Clockwise) {
@@ -287,7 +294,15 @@ namespace OpenBabel {
 
   inline double Calculate13Angle(double a, double b, double c)
   {
-    return acos((SQUARE(a) + SQUARE(b) - SQUARE(c)) / (2.0*a*b));
+    double cosine = (SQUARE(a) + SQUARE(b) - SQUARE(c)) / (2.0*a*b);
+
+    // Handle cases where cosine is outside [-1, 1] range.
+    if (cosine > 1.0)
+      return 0.0;
+    if (cosine < -1.0)
+      return M_PI;
+
+    return acos(cosine);
   }
 
   // When atoms i and j are in a 1-3 relationship, the distance
@@ -911,12 +926,16 @@ namespace OpenBabel {
 
   bool OBDistanceGeometry::CheckStereoConstraints()
   {
+    return _d->stereoHelper.Check(&_mol);
+
+    /*
     // Check stereo by canonical SMILES
     StereoFrom3D(&_mol, true);
     OBConversion conv;
     conv.SetOutFormat("can");
     std::string predicted_smiles = conv.WriteString(&_mol, true);
     return input_smiles == predicted_smiles;
+    */
 
     // Check all stereo constraints
     // First, gather the known, specified stereochemistry
@@ -1164,22 +1183,24 @@ namespace OpenBabel {
       cerr << " max box size: " << _d->maxBoxSize << endl;
     }
 
-    bool success = false;
-    unsigned int maxIter = 1 * _mol.NumAtoms();
+    _d->success = false;
+    unsigned int maxIter = 10 * _mol.NumAtoms();
     for (unsigned int trial = 0; trial < maxIter; trial++) {
       generateInitialCoords();
       firstMinimization();
       if (dim == 4) minimizeFourthDimension();
       if (CheckStereoConstraints() && CheckBounds()) {
-        success = true;
+        _d->success = true;
         break;
       }
-      if (_d->debug && !success)
+      if (_d->debug && !_d->success)
         cerr << "Stereo unsatisfied, trying again" << endl;
     }
-    if(!success) {
-      obErrorLog.ThrowError(__FUNCTION__, "Distance Geometry failed.", obWarning);
-    }
+  }
+
+  bool OBDistanceGeometry::WasSuccessful() const
+  {
+    return _d->success;
   }
 
   bool OBDistanceGeometry::CheckBounds()
@@ -1262,7 +1283,7 @@ namespace OpenBabel {
     AddConformer();
     GetConformers(mol);
 
-    return true;
+    return _d->success;
   }
 
   double DistGeomFunc::operator() (const Eigen::VectorXd& x, Eigen::VectorXd& grad) {

src/formats/mdlformat.cpp

@@ -857,14 +857,28 @@ namespace OpenBabel
         mol.SetDimension(3);
       // use 3D coordinates to determine stereochemistry
       StereoFrom3D(&mol);
+      OpenBabel::OBStereoFacade facade(&mol);
+
       if (pConv->IsOption("s", OBConversion::INOPTIONS)) { // Use the parities for tet stereo instead
         TetStereoFromParity(mol, parities, true); // True means "delete existing TetStereo first"
+      } else {
+        // Set stereo to unspecified for atom stereo parity 3 (1 & 2 determined from 3D coords)
+        for (std::size_t i = 0; i < parities.size(); ++i) {
+          if (parities[i] != Unknown)
+            continue;
+          unsigned long atomId = mol.GetAtom(i+1)->GetId();
+          OBTetrahedralStereo *ts = facade.GetTetrahedralStereo(atomId);
+          if (!ts)
+            continue;
+          OBTetrahedralStereo::Config config = ts->GetConfig();
+          config.specified = false;
+          ts->SetConfig(config);
+        }
       }
 
       // For unspecified cis/trans stereos, set their Configs to unspecified
       // This should really be done in CisTransFrom3D like in CisTransFrom2D but can't change the API now :-/
       map<OBBond*, OBStereo::BondDirection>::const_iterator bd_it;
-      OpenBabel::OBStereoFacade facade(&mol);
       for(bd_it=updown.begin(); bd_it!=updown.end(); ++bd_it) {
         OBBond* bond = bd_it->first;
         if (bond->GetBondOrder()!=2 || bd_it->second != OBStereo::UnknownDir)