Fix EPAM Indigo parsing so that the S groups needs to saved before collection block

core/indigo-core/molecule/src/molfile_saver.cpp

@@ -787,94 +787,6 @@ void MolfileSaver::_writeCtab(Output& output, BaseMolecule& mol, bool query)
 
     output.writeStringCR("M  V30 END BOND");
 
-    MoleculeStereocenters& stereocenters = mol.stereocenters;
-
-    if (stereocenters.begin() != stereocenters.end() || mol.hasHighlighting())
-    {
-        output.writeStringCR("M  V30 BEGIN COLLECTION");
-
-        QS_DEF(Array<int>, processed);
-
-        processed.clear_resize(mol.vertexEnd());
-        processed.zerofill();
-
-        for (i = mol.vertexBegin(); i != mol.vertexEnd(); i = mol.vertexNext(i))
-        {
-            if (processed[i])
-                continue;
-
-            ArrayOutput out(buf);
-            int j, type = stereocenters.getType(i);
-
-            if (type == MoleculeStereocenters::ATOM_ABS)
-                out.writeString("MDLV30/STEABS ATOMS=(");
-            else if (type == MoleculeStereocenters::ATOM_OR)
-                out.printf("MDLV30/STEREL%d ATOMS=(", stereocenters.getGroup(i));
-            else if (type == MoleculeStereocenters::ATOM_AND)
-                out.printf("MDLV30/STERAC%d ATOMS=(", stereocenters.getGroup(i));
-            else
-                continue;
-
-            QS_DEF(Array<int>, list);
-
-            list.clear();
-            list.push(i);
-
-            for (j = mol.vertexNext(i); j < mol.vertexEnd(); j = mol.vertexNext(j))
-                if (stereocenters.sameGroup(i, j))
-                {
-                    list.push(j);
-                    processed[j] = 1;
-                }
-
-            out.printf("%d", list.size());
-            for (j = 0; j < list.size(); j++)
-                out.printf(" %d", _atom_mapping[list[j]]);
-            out.writeChar(')');
-
-            _writeMultiString(output, buf.ptr(), buf.size());
-        }
-
-        if (mol.hasHighlighting())
-        {
-            if (mol.countHighlightedBonds() > 0)
-            {
-                ArrayOutput out(buf);
-
-                out.printf("MDLV30/HILITE BONDS=(%d", mol.countHighlightedBonds());
-
-                for (i = mol.edgeBegin(); i != mol.edgeEnd(); i = mol.edgeNext(i))
-                    if (mol.isBondHighlighted(i))
-                        out.printf(" %d", _bond_mapping[i]);
-                out.writeChar(')');
-
-                _writeMultiString(output, buf.ptr(), buf.size());
-            }
-            if (mol.countHighlightedAtoms() > 0)
-            {
-                ArrayOutput out(buf);
-                out.printf("MDLV30/HILITE ATOMS=(%d", mol.countHighlightedAtoms());
-                for (i = mol.vertexBegin(); i != mol.vertexEnd(); i = mol.vertexNext(i))
-                    if (mol.isAtomHighlighted(i))
-                        out.printf(" %d", _atom_mapping[i]);
-                out.writeChar(')');
-
-                _writeMultiString(output, buf.ptr(), buf.size());
-            }
-        }
-        if (mol.custom_collections.size() > 0)
-        {
-            for (i = mol.custom_collections.begin(); i != mol.custom_collections.end(); i = mol.custom_collections.next(i))
-            {
-                ArrayOutput out(buf);
-                out.printf("%s", mol.custom_collections.at(i));
-                _writeMultiString(output, buf.ptr(), buf.size());
-            }
-        }
-
-        output.writeStringCR("M  V30 END COLLECTION");
-    }
-
     QS_DEF(Array<int>, sgs_sorted);
     _checkSGroupIndices(mol, sgs_sorted);
 
@@ -1074,6 +986,94 @@ void MolfileSaver::_writeCtab(Output& output, BaseMolecule& mol, bool query)
         _removeImplicitSGroups(mol, implicit_sgroups_indexes);
     }
 
+    MoleculeStereocenters& stereocenters = mol.stereocenters;
+
+    if (stereocenters.begin() != stereocenters.end() || mol.hasHighlighting())
+    {
+        output.writeStringCR("M  V30 BEGIN COLLECTION");
+
+        QS_DEF(Array<int>, processed);
+
+        processed.clear_resize(mol.vertexEnd());
+        processed.zerofill();
+
+        for (i = mol.vertexBegin(); i != mol.vertexEnd(); i = mol.vertexNext(i))
+        {
+            if (processed[i])
+                continue;
+
+            ArrayOutput out(buf);
+            int j, type = stereocenters.getType(i);
+
+            if (type == MoleculeStereocenters::ATOM_ABS)
+                out.writeString("MDLV30/STEABS ATOMS=(");
+            else if (type == MoleculeStereocenters::ATOM_OR)
+                out.printf("MDLV30/STEREL%d ATOMS=(", stereocenters.getGroup(i));
+            else if (type == MoleculeStereocenters::ATOM_AND)
+                out.printf("MDLV30/STERAC%d ATOMS=(", stereocenters.getGroup(i));
+            else
+                continue;
+
+            QS_DEF(Array<int>, list);
+
+            list.clear();
+            list.push(i);
+
+            for (j = mol.vertexNext(i); j < mol.vertexEnd(); j = mol.vertexNext(j))
+                if (stereocenters.sameGroup(i, j))
+                {
+                    list.push(j);
+                    processed[j] = 1;
+                }
+
+            out.printf("%d", list.size());
+            for (j = 0; j < list.size(); j++)
+                out.printf(" %d", _atom_mapping[list[j]]);
+            out.writeChar(')');
+
+            _writeMultiString(output, buf.ptr(), buf.size());
+        }
+
+        if (mol.hasHighlighting())
+        {
+            if (mol.countHighlightedBonds() > 0)
+            {
+                ArrayOutput out(buf);
+
+                out.printf("MDLV30/HILITE BONDS=(%d", mol.countHighlightedBonds());
+
+                for (i = mol.edgeBegin(); i != mol.edgeEnd(); i = mol.edgeNext(i))
+                    if (mol.isBondHighlighted(i))
+                        out.printf(" %d", _bond_mapping[i]);
+                out.writeChar(')');
+
+                _writeMultiString(output, buf.ptr(), buf.size());
+            }
+            if (mol.countHighlightedAtoms() > 0)
+            {
+                ArrayOutput out(buf);
+                out.printf("MDLV30/HILITE ATOMS=(%d", mol.countHighlightedAtoms());
+                for (i = mol.vertexBegin(); i != mol.vertexEnd(); i = mol.vertexNext(i))
+                    if (mol.isAtomHighlighted(i))
+                        out.printf(" %d", _atom_mapping[i]);
+                out.writeChar(')');
+
+                _writeMultiString(output, buf.ptr(), buf.size());
+            }
+        }
+        if (mol.custom_collections.size() > 0)
+        {
+            for (i = mol.custom_collections.begin(); i != mol.custom_collections.end(); i = mol.custom_collections.next(i))
+            {
+                ArrayOutput out(buf);
+                out.printf("%s", mol.custom_collections.at(i));
+                _writeMultiString(output, buf.ptr(), buf.size());
+            }
+        }
+
+        output.writeStringCR("M  V30 END COLLECTION");
+    }
+
     output.writeStringCR("M  V30 END CTAB");
 
     int n_rgroups = mol.rgroups.getRGroupCount();
@@ -1212,6 +1212,7 @@ void MolfileSaver::_writeTGroup(Output& output, BaseMolecule& mol, int tg_idx)
 
 void MolfileSaver::_writeCtab2000(Output& output, BaseMolecule& mol, bool query)
 {
+    add_implicit_h = false;
     _handleCIP(mol);
     QueryMolecule* qmol = nullptr;
 

core/indigo-core/tests/common.cpp

@@ -1,4 +1,6 @@
 #include "common.h"
+#include "reaction/reaction.h"
+#include "reaction/reaction_auto_loader.h"
 #include "reaction/reaction_cml_saver.h"
 #include "reaction/reaction_json_saver.h"
 #include "reaction/rxnfile_saver.h"
@@ -25,6 +27,13 @@ void IndigoCoreTest::loadMolecule(const char* buf, Molecule& molecule)
     loader.loadMolecule(molecule);
 }
 
+void IndigoCoreTest::loadReaction(const char* buf, Reaction& reaction)
+{
+    BufferScanner scanner(buf);
+    ReactionAutoLoader loader(scanner);
+    loader.loadReaction(reaction);
+}
+
 void IndigoCoreTest::loadQueryMolecule(const char* buf, QueryMolecule& queryMolecule)
 {
     BufferScanner scanner(buf);

core/indigo-core/tests/tests/molecule.cpp

@@ -16,10 +16,17 @@
  * limitations under the License.
  ***************************************************************************/
 
+#ifdef _WIN32
+#include <direct.h>
+#define getcwd _getcwd
+#else
+#include <unistd.h>
+#endif
 #include <gtest/gtest.h>
 
 #include <base_cpp/output.h>
 #include <base_cpp/scanner.h>
+#include <fstream>
 #include <molecule/crippen.h>
 #include <molecule/hybridization.h>
 #include <molecule/lipinski.h>
@@ -29,6 +36,12 @@
 #include <molecule/tpsa.h>
 
 #include "common.h"
+#include "molecule/elements.h"
+#include "molecule/molfile_loader.h"
+#include "molecule/molfile_saver.h"
+#include "reaction/reaction.h"
+#include "reaction/reaction_auto_loader.h"
+#include "reaction/reaction_automapper.h"
 
 using namespace std;
 using namespace indigo;
@@ -352,3 +365,106 @@ TEST_F(IndigoCoreMoleculeTest, dearomatize_smarts)
     // printf("%s", sm.c_str());
     EXPECT_STREQ("c1-c=c-c=c-c=1", sm.c_str());
 }
+
+/*
+ * V3000 molfile roundtrip with add_implicit_h = true.
+ * Ensures SGROUP block is written before COLLECTION (CT file spec order) and
+ * that implicit H (e.g. MRV_IMPLICIT_H data sgroups) roundtrip correctly.
+ */
+TEST_F(IndigoCoreMoleculeTest, molfileV3000RoundtripImplicitH)
+{
+    Molecule mol;
+    loadMolecule("c1ccccc1O", mol); // phenol: aromatic carbons have implicit H
+
+    string molfileStr;
+    StringOutput out(molfileStr);
+    MolfileSaver saver(out);
+    saver.mode = MolfileSaver::MODE_3000;
+    saver.add_implicit_h = true;
+    saver.saveMolecule(mol);
+
+    // CT file order: SGROUP block must appear before COLLECTION block
+    size_t posSgroup = molfileStr.find("M  V30 BEGIN SGROUP");
+    size_t posCollection = molfileStr.find("M  V30 BEGIN COLLECTION");
+    if (posSgroup != string::npos && posCollection != string::npos)
+        EXPECT_LT(posSgroup, posCollection) << "SGROUP must be written before COLLECTION (CT file spec)";
+
+    BufferScanner scanner(molfileStr.c_str());
+    MolfileLoader loader(scanner);
+    Molecule mol2;
+    loader.loadMolecule(mol2);
+
+    EXPECT_EQ(mol.vertexCount(), mol2.vertexCount()) << "Roundtrip: atom count must match";
+    EXPECT_EQ(mol.edgeCount(), mol2.edgeCount()) << "Roundtrip: bond count must match";
+    // Same structure: roundtrip mol2 should match original as substructure and vice versa
+    EXPECT_TRUE(substructureMatch(smiles(mol2).c_str(), smiles(mol).c_str()));
+    EXPECT_TRUE(substructureMatch(smiles(mol).c_str(), smiles(mol2).c_str()));
+}
+
+/*
+ * V3000 molfile roundtrip with add_implicit_h = false (explicit H only).
+ * Ensures molecules without MRV_IMPLICIT_H sgroups roundtrip correctly.
+ */
+TEST_F(IndigoCoreMoleculeTest, molfileV3000RoundtripExplicitH)
+{
+    Molecule mol;
+    loadMolecule("CCO", mol);
+
+    string molfileStr;
+    StringOutput out(molfileStr);
+    MolfileSaver saver(out);
+    saver.mode = MolfileSaver::MODE_3000;
+    saver.add_implicit_h = false;
+    saver.saveMolecule(mol);
+
+    size_t posSgroup = molfileStr.find("M  V30 BEGIN SGROUP");
+    size_t posCollection = molfileStr.find("M  V30 BEGIN COLLECTION");
+    if (posSgroup != string::npos && posCollection != string::npos)
+        EXPECT_LT(posSgroup, posCollection) << "SGROUP must be written before COLLECTION (CT file spec)";
+
+    BufferScanner scanner(molfileStr.c_str());
+    MolfileLoader loader(scanner);
+    Molecule mol2;
+    loader.loadMolecule(mol2);
+
+    EXPECT_EQ(mol.vertexCount(), mol2.vertexCount()) << "Roundtrip: atom count must match";
+    EXPECT_EQ(mol.edgeCount(), mol2.edgeCount()) << "Roundtrip: bond count must match";
+    EXPECT_STREQ(smiles(mol).c_str(), smiles(mol2).c_str()) << "Roundtrip: SMILES must match";
+}
+
+TEST_F(IndigoCoreMoleculeTest, Reaction)
+{
+    Reaction reaction;
+    {
+        char dir[250];
+        getcwd(dir, 250);
+        cout << "Current directory: " << dir << "\n";
+        ifstream reactionFile("../../data/reactions/other/stereo_reaction.rxn");
+        string content;
+        string line;
+        while (getline(reactionFile, line))
+        {
+            content += line;
+            content.push_back('\n');
+        }
+        reactionFile.close();
+        //        cout << "The content is: \n" << content << "\n";
+        cout << "Loading reaction...";
+        loadReaction(content.c_str(), reaction);
+
+        reaction.aromatize(AromaticityOptions(AromaticityOptions::GENERIC));
+        ReactionAutomapper ram(reaction);
+        ram.automap(0);
+
+        for (int i = reaction.productBegin(); i < reaction.productEnd(); i = reaction.productNext(i))
+        {
+            cout << " *** Current Product Index is: " << i << " *** \n";
+            Molecule& mol = reaction.getMolecule(i);
+            string molOutStr;
+            StringOutput molOut(molOutStr);
+            MolfileSaver molSaver(molOut);
+            molSaver.saveMolecule(mol);
+            cout << molOutStr;
+        }
+    }
+}