#include "utils.h" #include #include #include #include #include namespace RDKit { namespace ChemDraw { std::string NodeType(CDXNodeType nodetype) { switch (nodetype) { case kCDXNodeType_Unspecified: return "Unspecified"; case kCDXNodeType_Element: return "Element"; case kCDXNodeType_ElementList: return "ElementList"; case kCDXNodeType_ElementListNickname: return "ElementListNickname"; case kCDXNodeType_Nickname: return "Nickname"; case kCDXNodeType_Fragment: return "Fragment"; case kCDXNodeType_Formula: return "Formula"; case kCDXNodeType_GenericNickname: return "GenericNickname"; case kCDXNodeType_AnonymousAlternativeGroup: return "Anonymous Alternative Group"; case kCDXNodeType_NamedAlternativeGroup: return "Named Alternative Group"; case kCDXNodeType_MultiAttachment: return "MultiAttachment"; case kCDXNodeType_VariableAttachment: return "Variable Attachment"; case kCDXNodeType_ExternalConnectionPoint: return "ExternalConnectionPoint"; case kCDXNodeType_LinkNode: return "LinkNode"; case kCDXNodeType_Monomer: return "Monomer"; default: return "?"; } } void scaleBonds(const ROMol &mol, Conformer &conf, double targetBondLength, double bondLength) { double avg_bond_length = 0.0; if (bondLength < 0) { // If we don't have a bond length for any reason, just scale the average // bond length for (auto &bond : mol.bonds()) { avg_bond_length += (conf.getAtomPos(bond->getBeginAtomIdx()) - conf.getAtomPos(bond->getEndAtomIdx())) .length(); } avg_bond_length /= mol.getNumBonds(); } else { avg_bond_length = bondLength; } if (avg_bond_length > 0) { double scale = targetBondLength / avg_bond_length; for (auto &pos : conf.getPositions()) { pos *= scale; } } } unsigned int get_fuse_label(Atom *atm) { // return atm->getAtomMapNum(); easier debugging unsigned int label = 0; // default is no label atm->getPropIfPresent(FUSE_LABEL, label); return label; } void set_fuse_label(Atom *atm, unsigned int idx) { // atm->setAtomMapNum(idx); //for debugging if (idx) { atm->setProp(FUSE_LABEL, idx); } else { atm->clearProp(FUSE_LABEL); } } struct FragmentReplacement { // R = Replacement // F = Fragment // C = Connection // C R C F F // N=*=C.*=CCC=* // label 1 1 1 // has bond ordering // // goal replace the atom R with the connections unsigned int label = 0; Atom *replacement_atom = nullptr; std::vector replacement_connection_atoms; std::vector fragment_atoms; bool replace(RWMol &mol) { if (!replacement_atom) { return true; } auto bond_ordering = replacement_atom->getProp>(CDX_BOND_ORDERING); // The "addBond" lower in the loop potentially modifies the atomBonds // iterator. To ensure safety, we copy the bonds first. std::vector replacement_bonds( mol.atomBonds(replacement_atom).begin(), mol.atomBonds(replacement_atom).end()); std::vector xbonds; // Reuse this vector to reduce repeated allocations // Find the connecting atoms and and do the replacement for (auto bond : replacement_bonds) { // find the position of the attachment bonds in the bond ordering unsigned bond_id = 0; if (!bond->getPropIfPresent(CDX_BOND_ID, bond_id)) { BOOST_LOG(rdWarningLog) << "bond missing internal CDX BOND id, can't attach fragment at bond:" << std::endl; return false; } auto it = std::find(bond_ordering.begin(), bond_ordering.end(), bond_id); if (it == bond_ordering.end()) { return false; } auto pos = std::distance(bond_ordering.begin(), it); if (pos < 0 || (size_t)pos >= fragment_atoms.size()) { BOOST_LOG(rdWarningLog) << "bond ordering and number of atoms in fragment mismatch, can't attach fragment at bond:" << bond_id << std::endl; return false; } auto &xatom = fragment_atoms[pos]; // The "addBond" lower in the loop potentially modifies the atomBonds // iterator. To ensure safety, we copy the bonds first. xbonds.assign(mol.atomBonds(xatom).begin(), mol.atomBonds(xatom).end()); for (auto &xbond : xbonds) { // xatom is the fragment dummy atom // xbond is the fragment bond if (bond->getBeginAtom() == replacement_atom) { mol.addBond(xbond->getOtherAtom(xatom), bond->getEndAtom(), bond->getBondType()); } else { mol.addBond(bond->getBeginAtom(), xbond->getOtherAtom(xatom), bond->getBondType()); } } } mol.removeAtom(replacement_atom); for (auto &atom : fragment_atoms) { mol.removeAtom(atom); } return true; } }; // Replace fragments that are not possible with molzip bool replaceFragments(RWMol &mol) { // Anything with a single atom that is supposed to be replaced via a fragment // is here std::map replacements; for (auto &atom : mol.atoms()) { auto label = get_fuse_label(atom); if (label) { if (atom->hasProp(CDX_BOND_ORDERING)) { auto &frag = replacements[label]; frag.label = label; frag.replacement_atom = atom; } else { // The is the fragment attachment atoms that need to // be attached to the ones connected to the atom being replaced auto &frag = replacements[label]; frag.fragment_atoms.push_back(atom); } } } mol.beginBatchEdit(); for (auto &replacement : replacements) { replacement.second.replace(mol); } mol.commitBatchEdit(); return true; } namespace { Atom::ChiralType getChirality(ROMol &mol, Atom *center_atom, Conformer &conf) { if (center_atom->hasProp(CDX_BOND_ORDERING)) { auto bond_ordering = center_atom->getProp>(CDX_BOND_ORDERING); if (bond_ordering.size() < 3) { return Atom::ChiralType::CHI_UNSPECIFIED; } std::vector atoms; std::vector> angles; auto center = conf.getAtomPos(center_atom->getIdx()); for (auto cdx_id : bond_ordering) { if (cdx_id == 0) { continue; } for (auto bond : mol.atomBonds(center_atom)) { int bond_id; if (bond->getPropIfPresent(CDX_BOND_ID, bond_id)) { } else { return Atom::ChiralType::CHI_UNSPECIFIED; } if (bond_id == cdx_id) { auto atom = bond->getOtherAtom(center_atom); if (!atom) { // something went really wrong return Atom::ChiralType::CHI_UNSPECIFIED; } auto pos = conf.getAtomPos(atom->getIdx()) - center; double angle = atan2(pos.x, pos.y); angles.push_back(std::make_pair(angle, bond->getIdx())); } } } std::sort(angles.begin(), angles.end()); // angles are now sorted in a clockwise rotation INT_LIST bonds; for (auto &angle : angles) { bonds.push_back(angle.second); } if (bonds.size() < 3) { return Atom::ChiralType::CHI_UNSPECIFIED; } auto nswaps = center_atom->getPerturbationOrder(bonds); if (bonds.size() == 3 && center_atom->getTotalNumHs() == 1) { ++nswaps; } // This is supports the HDot and HDash available in chemdraw // one is an implicit wedged hydrogen and one is a dashed hydrogen if (center_atom->hasProp(CDX_IMPLICIT_HYDROGEN_STEREO) && center_atom->getProp(CDX_IMPLICIT_HYDROGEN_STEREO) == 'w') { nswaps++; } if (nswaps % 2) { return Atom::ChiralType::CHI_TETRAHEDRAL_CCW; } return Atom::ChiralType::CHI_TETRAHEDRAL_CW; } return Atom::ChiralType::CHI_UNSPECIFIED; } } // namespace void checkChemDrawTetrahedralGeometries(RWMol &mol) { std::vector> unsetTetrahedralAtoms; Conformer *conf = nullptr; if (mol.getNumConformers()) { conf = &mol.getConformer(); } bool chiralityChanged = false; for (auto atom : mol.atoms()) { // only deal with unspecified chiralities if (atom->getChiralTag() != Atom::ChiralType::CHI_UNSPECIFIED) { atom->clearProp(CDX_CIP); continue; } if (conf && !conf->is3D()) { atom->setChiralTag(getChirality(mol, atom, *conf)); if (atom->getChiralTag() != Atom::ChiralType::CHI_UNSPECIFIED) { chiralityChanged = true; } } // If we have a cip code, might as well check it too CDXAtomCIPType cip; if (atom->getPropIfPresent(CDX_CIP, cip)) { // assign, possibly wrong, initial stereo. // note: we can probably deduce this through CDX_BOND_ORDERING, but // I currently don't understand that well enough. switch (cip) { case kCDXCIPAtom_R: if (!chiralityChanged) { atom->setChiralTag(Atom::ChiralType::CHI_TETRAHEDRAL_CW); } unsetTetrahedralAtoms.push_back(std::make_pair('R', atom)); break; case kCDXCIPAtom_r: if (!chiralityChanged) { atom->setChiralTag(Atom::ChiralType::CHI_TETRAHEDRAL_CW); } unsetTetrahedralAtoms.push_back(std::make_pair('r', atom)); break; case kCDXCIPAtom_S: if (!chiralityChanged) { atom->setChiralTag(Atom::ChiralType::CHI_TETRAHEDRAL_CW); } unsetTetrahedralAtoms.push_back(std::make_pair('S', atom)); break; case kCDXCIPAtom_s: if (!chiralityChanged) { atom->setChiralTag(Atom::ChiralType::CHI_TETRAHEDRAL_CCW); } unsetTetrahedralAtoms.push_back(std::make_pair('s', atom)); break; default: break; } } } // Now that we have missing chiralities, let's check the CIP codes and reset // if necessary. // This is an expensive way of doing this, but we only have stereo->cip not // cip->stereo implemented currently for (auto cipatom : unsetTetrahedralAtoms) { try { CIPLabeler::assignCIPLabels(mol); } catch (...) { // can throw std::runtime error? break; } std::string cipcode; if (cipatom.second->getPropIfPresent( common_properties::_CIPCode, cipcode)) { if (cipcode.size() && cipcode[0] != cipatom.first) { // need to swap if (cipatom.second->getChiralTag() == Atom::ChiralType::CHI_TETRAHEDRAL_CW) { cipatom.second->setChiralTag(Atom::ChiralType::CHI_TETRAHEDRAL_CCW); cipatom.second->updatePropertyCache(); chiralityChanged = true; } else if (cipatom.second->getChiralTag() == Atom::ChiralType::CHI_TETRAHEDRAL_CCW) { cipatom.second->setChiralTag(Atom::ChiralType::CHI_TETRAHEDRAL_CW); cipatom.second->updatePropertyCache(); chiralityChanged = true; } } } } if (chiralityChanged) { const bool cleanIt = true; const bool force = true; MolOps::assignStereochemistry(mol, cleanIt, force); } } } // namespace ChemDraw } // namespace RDKit