Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
45 changes: 45 additions & 0 deletions src/schrodinger/rdkit_extensions/stereochemistry.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -72,9 +72,43 @@ void wedgeMolBonds(RDKit::ROMol& mol, const RDKit::Conformer* conf)
}
}

// Save wiggly bonds (unknown stereochemistry) before
// ClearSingleBondDirFlags clears them. We check both the properties (set by
// MDL input) and BondDir (set when CXSMILES wiggly bonds are copied in
// mol_model.cpp).
std::vector<bool> had_wiggly_bond;
had_wiggly_bond.reserve(mol.getNumBonds());
for (auto bond : mol.bonds()) {
int wiggly_bond_v2000{0};
int wiggly_bond_v3000{0};
bond->getPropIfPresent(RDKit::common_properties::_MolFileBondStereo,
wiggly_bond_v2000);
bond->getPropIfPresent(RDKit::common_properties::_MolFileBondCfg,
wiggly_bond_v3000);
// MDL V2000 uses _MolFileBondStereo=4 for wiggly bonds
// MDL V3000 uses _MolFileBondCfg=2 for wiggly bonds
// (These values are from the MDL molfile specification, not defined as
// constants in RDKit)
bool is_wiggly = (wiggly_bond_v2000 == 4 || wiggly_bond_v3000 == 2 ||
bond->getBondDir() == RDKit::Bond::BondDir::UNKNOWN);
had_wiggly_bond.push_back(is_wiggly);
}

RDKit::ClearSingleBondDirFlags(mol);
RDKit::Chirality::clearMolBlockWedgingInfo(mol);

// Temporarily change wiggly bonds to OTHER type so WedgeMolBonds skips
// them. WedgeMolBonds only considers Bond::SINGLE bonds, so this prevents
// it from stealing a wiggly bond to express the stereo of an adjacent
// chiral atom — which would be overwritten when we restore UNKNOWN,
// losing that atom's stereo representation. This mirrors the same technique
// used above for attachment dummy bonds.
for (size_t bond_idx = 0; bond_idx < mol.getNumBonds(); ++bond_idx) {
if (had_wiggly_bond[bond_idx]) {
mol.getBondWithIdx(bond_idx)->setBondType(RDKit::Bond::OTHER);
}
}

try {
// Temporarily silence RDKit's loggers
RDLog::LogStateSetter silence_rdkit_logging;
Expand All @@ -88,6 +122,17 @@ void wedgeMolBonds(RDKit::ROMol& mol, const RDKit::Conformer* conf)
for (auto bond : attachment_dummy_bonds) {
bond->setBondType(RDKit::Bond::SINGLE);
}

// Restore wiggly bonds: return to SINGLE type and set UNKNOWN direction.
// WedgeMolBonds skipped them (OTHER type), so neighboring chiral atoms had
// their stereo assigned to other bonds.
for (size_t bond_idx = 0; bond_idx < mol.getNumBonds(); ++bond_idx) {
if (had_wiggly_bond[bond_idx]) {
auto bond = mol.getBondWithIdx(bond_idx);
bond->setBondType(RDKit::Bond::SINGLE);
bond->setBondDir(RDKit::Bond::BondDir::UNKNOWN);
}
}
}

} // namespace rdkit_extensions
Expand Down
15 changes: 15 additions & 0 deletions src/schrodinger/sketcher/model/mol_model.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -2671,6 +2671,21 @@ void MolModel::addMolCommandFunc(RDKit::ROMol mol)
++bond_index) {
RDKit::Bond* bond = m_mol.getBondWithIdx(bond_index);
setTagForBond(bond, m_next_bond_tag++);

// Preserve wiggly bonds (BondDir::UNKNOWN) from the source molecule.
// RDKit's insertMol preserves BondDir via bond copy, but the explicit
// call ensures correctness across all code paths. Additionally, for
// CXSMILES input RDKit sets BondDir::UNKNOWN but not _MolFileBondCfg;
// set _MolFileBondCfg=2 so MolToCXSmiles can detect the wiggly bond
// after clearMolBlockWedgingInfo clears the CXSMILES wiggly annotation.
// (2 = wiggly bond in MDL V3000 format)
auto* src_bond = mol.getBondWithIdx(bond_index - old_num_bonds);
bond->setBondDir(src_bond->getBondDir());
if (src_bond->getBondDir() == RDKit::Bond::BondDir::UNKNOWN &&
!src_bond->hasProp(RDKit::common_properties::_MolFileBondCfg)) {
bond->setProp(RDKit::common_properties::_MolFileBondCfg, 2);
}

if (contains_two_monomer_linkages(bond)) {
setSecondaryConnectionTagForBond(bond, m_next_bond_tag++);
}
Expand Down
12 changes: 12 additions & 0 deletions test/schrodinger/rdkit_extensions/test_convert.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -1257,3 +1257,15 @@ M END
BOOST_TEST(cxsmiles == "CC[C@H](C)[C@H](C)C[C@H](C)N |a:4,7,o1:2|");
}
}

BOOST_AUTO_TEST_CASE(test_wiggly_bond_cxsmiles_roundtrip)
{
// Test that wiggly bonds are preserved when round-tripping through
// to_rdkit and to_string with EXTENDED_SMILES format
const std::string input_smiles = "CCC(C)N |w:2.3|";
auto mol = to_rdkit(input_smiles, Format::EXTENDED_SMILES);
BOOST_REQUIRE(mol != nullptr);

std::string output_smiles = to_string(*mol, Format::EXTENDED_SMILES);
BOOST_TEST(output_smiles == input_smiles);
}
Original file line number Diff line number Diff line change
@@ -0,0 +1,69 @@
/* -------------------------------------------------------------------------
* Tests demonstrating RDKit behaviors with wiggly bonds (unknown
* stereochemistry). These tests document the specific RDKit behaviors that
* necessitate wiggly bond preservation code in the sketcher.
*
* Copyright Schrodinger LLC, All Rights Reserved.
--------------------------------------------------------------------------- */

#define BOOST_TEST_MODULE rdkit_extensions_wiggly_bond_behaviors

#include <boost/test/unit_test.hpp>

#include <rdkit/GraphMol/RWMol.h>
#include <rdkit/GraphMol/SmilesParse/SmilesParse.h>
#include <rdkit/GraphMol/SmilesParse/SmilesWrite.h>

/**
* Demonstrate that CXSMILES wiggly bonds lack _MolFileBondCfg property
* (which prepare_mol checks to identify wiggly bonds for preservation).
* This test also shows that insertMol preserves BondDir but not the property.
*/
BOOST_AUTO_TEST_CASE(test_cxsmiles_wiggly_bonds_lack_molfile_property)
{
// Create a molecule with a wiggly bond using CXSMILES
const std::string cxsmiles = "CCC(C)N |w:2.3|";
RDKit::SmilesParserParams params;
params.sanitize = true;
params.removeHs = false;
std::unique_ptr<RDKit::RWMol> mol(RDKit::SmilesToMol(cxsmiles, params));
BOOST_REQUIRE(mol != nullptr);

// Verify CXSMILES roundtrip preserves wiggly bond
std::string cxsmiles_output = RDKit::MolToCXSmiles(*mol);
BOOST_TEST(cxsmiles_output == cxsmiles);

// Get the wiggly bond (bond index 2: between atom 2 and atom 3)
auto* bond = mol->getBondWithIdx(2);
BOOST_REQUIRE(bond != nullptr);

// The key issue: CXSMILES input does NOT set _MolFileBondCfg property
// (unlike MDL input which sets both BondDir AND the property)
int cfg{-1};
bond->getPropIfPresent(RDKit::common_properties::_MolFileBondCfg, cfg);
BOOST_TEST(cfg == -1); // property not set for CXSMILES input

// RDKit's CXSMILES parser also doesn't set BondDir to UNKNOWN initially.
// The wiggly bond information is only preserved in the CXSMILES extension.
// When coordinates are generated and the molecule is processed, BondDir
// gets set.

// However, we can manually set BondDir to UNKNOWN (like mol_model.cpp does)
// to simulate what happens when the mol is copied
bond->setBondDir(RDKit::Bond::BondDir::UNKNOWN);
BOOST_TEST(bond->getBondDir() == RDKit::Bond::BondDir::UNKNOWN);

// Test insertMol preserves BondDir
RDKit::RWMol target_mol;
target_mol.insertMol(*mol);

auto* dest_bond = target_mol.getBondWithIdx(2);
BOOST_REQUIRE(dest_bond != nullptr);
BOOST_TEST(dest_bond->getBondDir() == RDKit::Bond::BondDir::UNKNOWN);

// But the property still isn't there - this is why we need to set it in
// mol_model.cpp
cfg = -1;
dest_bond->getPropIfPresent(RDKit::common_properties::_MolFileBondCfg, cfg);
BOOST_TEST(cfg == -1);
}
54 changes: 54 additions & 0 deletions test/schrodinger/sketcher/model/test_mol_model.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -4289,5 +4289,59 @@ BOOST_AUTO_TEST_CASE(test_stereo_labels_update_on_atom_movement)
<< initial_label);
}

/**
* Test that wiggly bonds (unknown stereochemistry) are preserved when
* round-tripping through add_text_to_mol_model and EXTENDED_SMILES export.
*/
BOOST_AUTO_TEST_CASE(test_wiggly_bond_preservation)
{
// Test MDL V3000 format roundtrip
QUndoStack undo_stack_mdl;
TestMolModel model_mdl(&undo_stack_mdl);

const std::string mdl_wiggly_bond = R"(
RDKit 2D

0 0 0 0 0 0 0 0 0 0999 V3000
M V30 BEGIN CTAB
M V30 COUNTS 5 4 0 0 0
M V30 BEGIN ATOM
M V30 1 C -2.078434 0.000000 0.000000 0
M V30 2 C -0.779413 -0.750012 0.000000 0
M V30 3 C 0.519609 0.000000 0.000000 0
M V30 4 C 1.818630 -0.750012 0.000000 0
M V30 5 N 0.519609 1.500025 0.000000 0
M V30 END ATOM
M V30 BEGIN BOND
M V30 1 1 1 2
M V30 2 1 2 3
M V30 3 1 3 4
M V30 4 1 3 5 CFG=2
M V30 END BOND
M V30 END CTAB
M END
$$$$
)";

add_text_to_mol_model(model_mdl, mdl_wiggly_bond);

// Verify EXTENDED_SMILES export shows wiggly bond from MDL input
std::string exported_smiles_from_mdl =
get_mol_text(&model_mdl, Format::EXTENDED_SMILES);
BOOST_TEST(exported_smiles_from_mdl == "CCC(C)N |w:2.3|");

// Test EXTENDED_SMILES format roundtrip
QUndoStack undo_stack_smiles;
TestMolModel model_smiles(&undo_stack_smiles);

// Add a molecule with a wiggly bond using add_text_to_mol_model
add_text_to_mol_model(model_smiles, "CCC(C)N |w:2.3|");

// Extract it as EXTENDED_SMILES and verify the wiggly bond is preserved
std::string exported_smiles =
get_mol_text(&model_smiles, Format::EXTENDED_SMILES);
BOOST_TEST(exported_smiles == "CCC(C)N |w:2.3|");
}

} // namespace sketcher
} // namespace schrodinger
83 changes: 83 additions & 0 deletions test/schrodinger/sketcher/rdkit/test_stereochemistry.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -101,5 +101,88 @@ BOOST_AUTO_TEST_CASE(test_wedgeMolBonds)
false);
}

/**
* Verify that wedgeMolBonds preserves wiggly bonds (unknown stereochemistry).
* Without this preservation, wiggly bonds would be lost when recalculating
* wedging, since ClearSingleBondDirFlags clears BondDir::UNKNOWN.
*/
BOOST_AUTO_TEST_CASE(test_wedgeMolBonds_preserves_wiggly_bonds)
{
// Create a molecule with a wiggly bond using CXSMILES
const std::string cxsmiles = "CCC(C)N |w:2.3|";
auto mol = rdkit_extensions::to_rdkit(
cxsmiles, rdkit_extensions::Format::EXTENDED_SMILES);
BOOST_REQUIRE(mol != nullptr);

// Generate 2D coordinates
rdkit_extensions::compute2DCoords(*mol);

// Set BondDir to UNKNOWN to simulate a wiggly bond
auto* bond = mol->getBondWithIdx(2);
BOOST_REQUIRE(bond != nullptr);
bond->setBondDir(RDKit::Bond::BondDir::UNKNOWN);
BOOST_TEST(bond->getBondDir() == RDKit::Bond::BondDir::UNKNOWN);

// Call wedgeMolBonds - it should preserve the wiggly bond
rdkit_extensions::wedgeMolBonds(*mol, &mol->getConformer());

// Verify that BondDir::UNKNOWN was preserved
BOOST_TEST(bond->getBondDir() == RDKit::Bond::BondDir::UNKNOWN);
}

/**
* Verify that wedgeMolBonds does not steal a wiggly bond to express the stereo
* of an adjacent chiral atom.
*
* WedgeMolBonds prefers terminal-neighbor bonds for wedging. If the wiggly bond
* connects to a terminal atom (N in this case), WedgeMolBonds would naively
* steal it to express the adjacent chiral atom's stereo. When we then restore
* it to UNKNOWN, the chiral atom loses its wedge — the stereo is no longer
* represented.
*
* The fix: before calling WedgeMolBonds, temporarily change wiggly bonds to
* Bond::OTHER type so WedgeMolBonds skips them entirely (it only considers
* Bond::SINGLE bonds), then restore them to SINGLE + UNKNOWN afterward.
*/
BOOST_AUTO_TEST_CASE(
test_wedgeMolBonds_does_not_steal_wiggly_bond_from_adjacent_chiral_atom)
{
// CC[C@@H](CC)N: atom 2 (C@@H) has defined stereo;
// the bond to N (atom 5, degree-1 terminal) is the wiggly bond.
// WedgeMolBonds prefers terminal-neighbor bonds, so without the fix it
// picks this bond for atom 2's stereo — which we then overwrite with
// UNKNOWN.
auto mol = rdkit_extensions::to_rdkit("CC[C@@H](CC)N",
rdkit_extensions::Format::SMILES);
BOOST_REQUIRE(mol != nullptr);
rdkit_extensions::compute2DCoords(*mol);

// Simulate wiggly bond state (as mol_model.cpp does when adding a molecule)
auto* wiggly_bond = mol->getBondBetweenAtoms(2, 5);
BOOST_REQUIRE(wiggly_bond != nullptr);
wiggly_bond->setBondDir(RDKit::Bond::BondDir::UNKNOWN);

rdkit_extensions::wedgeMolBonds(*mol, &mol->getConformer());

// The wiggly bond must remain UNKNOWN
BOOST_TEST(wiggly_bond->getBondDir() == RDKit::Bond::BondDir::UNKNOWN);

// The chiral atom (C@@H, atom 2) must still have its stereo represented
// by at least one other bond with a wedge/dash direction.
auto* stereo_atom = mol->getAtomWithIdx(2);
bool has_stereo_bond = false;
for (auto* b : mol->atomBonds(stereo_atom)) {
if (b != wiggly_bond &&
(b->getBondDir() == RDKit::Bond::BondDir::BEGINWEDGE ||
b->getBondDir() == RDKit::Bond::BondDir::BEGINDASH)) {
has_stereo_bond = true;
break;
}
}
BOOST_TEST(has_stereo_bond,
"The chiral atom adjacent to the wiggly bond should still have "
"its stereo expressed by a wedge or dash on another bond");
}

} // namespace sketcher
} // namespace schrodinger