rdkit/Code/GraphMol/ROMol.cpp

//
//  Copyright (C) 2003-2024 Greg Landrum and other RDKit contributors
//
//   @@ All Rights Reserved @@
//  This file is part of the RDKit.
//  The contents are covered by the terms of the BSD license
//  which is included in the file license.txt, found at the root
//  of the RDKit source tree.
//

// our stuff
#include <RDGeneral/Invariant.h>
#include <RDGeneral/RDLog.h>
#include "ROMol.h"
#include "Atom.h"
#include "QueryAtom.h"
#include "Bond.h"
#include "QueryBond.h"
#include "MolPickler.h"
#include "Conformer.h"
#include "SubstanceGroup.h"

#ifdef RDK_USE_BOOST_SERIALIZATION
#include <RDGeneral/BoostStartInclude.h>
#include <boost/archive/text_oarchive.hpp>
#include <boost/archive/text_iarchive.hpp>
#include <RDGeneral/BoostEndInclude.h>
#endif

namespace RDKit {
class QueryAtom;
class QueryBond;

const int ci_RIGHTMOST_ATOM = -0xBADBEEF;
const int ci_LEADING_BOND = -0xBADBEEF + 1;
const int ci_ATOM_HOLDER = -0xDEADD06;

void ROMol::destroy() {
  d_atomBookmarks.clear();
  d_bondBookmarks.clear();

  auto atItP = boost::vertices(d_graph);
  while (atItP.first != atItP.second) {
    delete (d_graph)[*(atItP.first++)];
  }

  auto bondItP = boost::edges(d_graph);
  while (bondItP.first != bondItP.second) {
    delete (d_graph)[*(bondItP.first++)];
  }

  d_graph.clear();

  delete dp_ringInfo;

  d_sgroups.clear();
  d_stereo_groups.clear();
}

ROMol::ROMol(const std::string &pickle) : RDProps() {
  initMol();
  numBonds = 0;
  MolPickler::molFromPickle(pickle, *this);
  numBonds = rdcast<unsigned int>(boost::num_edges(d_graph));
}

ROMol::ROMol(const std::string &pickle, unsigned int propertyFlags)
    : RDProps() {
  initMol();
  numBonds = 0;
  MolPickler::molFromPickle(pickle, *this, propertyFlags);
  numBonds = rdcast<unsigned int>(boost::num_edges(d_graph));
}

void ROMol::initFromOther(const ROMol &other, bool quickCopy, int confId) {
  if (this == &other) {
    return;
  }
  numBonds = 0;
  // std::cerr<<"    init from other: "<<this<<" "<<&other<<std::endl;
  // copy over the atoms
  // Avoid repeated reallocations when copying: for MolGraph's vecS vertex
  // container, reserving upfront can reduce allocation churn.
  d_graph.m_vertices.reserve(other.getNumAtoms());
  for (const auto oatom : other.atoms()) {
    constexpr bool updateLabel = false;
    constexpr bool takeOwnership = true;
    addAtom(oatom->copy(), updateLabel, takeOwnership);
  }

  // and the bonds:
  for (const auto obond : other.bonds()) {
    addBond(obond->copy(), true);
  }

  // ring information
  delete dp_ringInfo;
  if (other.dp_ringInfo) {
    dp_ringInfo = new RingInfo(*(other.dp_ringInfo));
  } else {
    dp_ringInfo = new RingInfo();
  }

  // enhanced stereochemical information
  d_stereo_groups.clear();
  d_stereo_groups.reserve(other.d_stereo_groups.size());
  for (auto &otherGroup : other.d_stereo_groups) {
    std::vector<Atom *> atoms;
    for (auto &otherAtom : otherGroup.getAtoms()) {
      atoms.push_back(getAtomWithIdx(otherAtom->getIdx()));
    }
    std::vector<Bond *> bonds;
    for (auto &otherBond : otherGroup.getBonds()) {
      bonds.push_back(getBondWithIdx(otherBond->getIdx()));
    }
    d_stereo_groups.emplace_back(otherGroup.getGroupType(), std::move(atoms),
                                 std::move(bonds), otherGroup.getReadId());
    d_stereo_groups.back().setWriteId(otherGroup.getWriteId());
  }

  if (other.dp_delAtoms) {
    dp_delAtoms.reset(new boost::dynamic_bitset<>(*other.dp_delAtoms));
  } else {
    dp_delAtoms.reset(nullptr);
  }
  if (other.dp_delBonds) {
    dp_delBonds.reset(new boost::dynamic_bitset<>(*other.dp_delBonds));
  } else {
    dp_delBonds.reset(nullptr);
  }

  if (!quickCopy) {
    // copy conformations
    for (const auto &conf : other.d_confs) {
      if (confId < 0 || rdcast<int>(conf->getId()) == confId) {
        this->addConformer(new Conformer(*conf));
      }
    }

    // Copy sgroups
    for (const auto &sg : getSubstanceGroups(other)) {
      addSubstanceGroup(*this, sg);
    }

    d_props = other.d_props;

    // Bookmarks should be copied as well:
    for (auto abmI : other.d_atomBookmarks) {
      for (const auto *aptr : abmI.second) {
        setAtomBookmark(getAtomWithIdx(aptr->getIdx()), abmI.first);
      }
    }
    for (auto bbmI : other.d_bondBookmarks) {
      for (const auto *bptr : bbmI.second) {
        setBondBookmark(getBondWithIdx(bptr->getIdx()), bbmI.first);
      }
    }
  } else {
    d_props.reset();
    STR_VECT computed;
    d_props.setVal(RDKit::detail::computedPropName, computed);
  }

  // std::cerr<<"---------    done init from other: "<<this<<"
  // "<<&other<<std::endl;
}

void ROMol::initMol() {
  d_props.reset();
  dp_ringInfo = new RingInfo();
}

unsigned int ROMol::getAtomDegree(const Atom *at) const {
  PRECONDITION(at, "no atom");
  PRECONDITION(&at->getOwningMol() == this,
               "atom not associated with this molecule");
  return rdcast<unsigned int>(boost::out_degree(at->getIdx(), d_graph));
};

unsigned int ROMol::getNumAtoms(bool onlyExplicit) const {
  int res = rdcast<int>(boost::num_vertices(d_graph));
  if (!onlyExplicit) {
    // if we are interested in hydrogens as well add them up from
    // each
    for (const auto atom : atoms()) {
      res += atom->getTotalNumHs();
    }
  }
  return res;
};
unsigned int ROMol::getNumHeavyAtoms() const {
  unsigned int res = 0;
  for (const auto atom : atoms()) {
    if (atom->getAtomicNum() > 1) {
      ++res;
    }
  }
  return res;
};

Atom *ROMol::getAtomWithIdx(unsigned int idx) {
  URANGE_CHECK(idx, getNumAtoms());

  auto vd = boost::vertex(idx, d_graph);
  auto res = d_graph[vd];
  POSTCONDITION(res, "");
  return res;
}

const Atom *ROMol::getAtomWithIdx(unsigned int idx) const {
  URANGE_CHECK(idx, getNumAtoms());

  auto vd = boost::vertex(idx, d_graph);
  const auto res = d_graph[vd];

  POSTCONDITION(res, "");
  return res;
}

// returns the first inserted atom with the given bookmark
Atom *ROMol::getAtomWithBookmark(int mark) {
  auto lu = d_atomBookmarks.find(mark);
  PRECONDITION((lu != d_atomBookmarks.end() && !lu->second.empty()),
               "atom bookmark not found");
  return lu->second.front();
};

// returns all atoms with the given bookmark
ROMol::ATOM_PTR_LIST &ROMol::getAllAtomsWithBookmark(int mark) {
  auto lu = d_atomBookmarks.find(mark);
  PRECONDITION(lu != d_atomBookmarks.end(), "atom bookmark not found");
  return lu->second;
};

// returns the unique atom with the given bookmark
Atom *ROMol::getUniqueAtomWithBookmark(int mark) {
  auto lu = d_atomBookmarks.find(mark);
  PRECONDITION((lu != d_atomBookmarks.end()), "bookmark not found");
  return lu->second.front();
}

// returns the first inserted bond with the given bookmark
Bond *ROMol::getBondWithBookmark(int mark) {
  auto lu = d_bondBookmarks.find(mark);
  PRECONDITION((lu != d_bondBookmarks.end() && !lu->second.empty()),
               "bond bookmark not found");
  return lu->second.front();
};

// returns all bonds with the given bookmark
ROMol::BOND_PTR_LIST &ROMol::getAllBondsWithBookmark(int mark) {
  auto lu = d_bondBookmarks.find(mark);
  PRECONDITION(lu != d_bondBookmarks.end(), "bond bookmark not found");
  return lu->second;
};

// returns the unique bond with the given bookmark
Bond *ROMol::getUniqueBondWithBookmark(int mark) {
  auto lu = d_bondBookmarks.find(mark);
  PRECONDITION((lu != d_bondBookmarks.end()), "bookmark not found");
  return lu->second.front();
}

void ROMol::clearAtomBookmark(const int mark) { d_atomBookmarks.erase(mark); }

void ROMol::clearAtomBookmark(int mark, const Atom *atom) {
  PRECONDITION(atom, "no atom");
  auto lu = d_atomBookmarks.find(mark);
  if (lu != d_atomBookmarks.end()) {
    auto &marks = lu->second;
    unsigned int tgtIdx = atom->getIdx();
    auto entry = std::find_if(marks.begin(), marks.end(), [&tgtIdx](auto ptr) {
      return ptr->getIdx() == tgtIdx;
    });
    if (entry != marks.end()) {
      marks.erase(entry);
    }
    if (marks.empty()) {
      d_atomBookmarks.erase(mark);
    }
  }
}

void ROMol::clearBondBookmark(int mark) { d_bondBookmarks.erase(mark); }
void ROMol::clearBondBookmark(int mark, const Bond *bond) {
  PRECONDITION(bond, "no bond");
  auto lu = d_bondBookmarks.find(mark);
  if (lu != d_bondBookmarks.end()) {
    auto &marks = lu->second;
    unsigned int tgtIdx = bond->getIdx();
    auto entry = std::find_if(marks.begin(), marks.end(), [&tgtIdx](auto ptr) {
      return ptr->getIdx() == tgtIdx;
    });
    if (entry != marks.end()) {
      marks.erase(entry);
    }
    if (marks.empty()) {
      d_bondBookmarks.erase(mark);
    }
  }
}

unsigned int ROMol::getNumBonds(bool onlyHeavy) const {
  // By default return the bonds that connect only the heavy atoms
  // hydrogen connecting bonds are ignores
  auto res = numBonds;
  if (!onlyHeavy) {
    // If we need hydrogen connecting bonds add them up
    for (const auto atom : atoms()) {
      res += atom->getTotalNumHs();
    }
  }
  return res;
}

Bond *ROMol::getBondWithIdx(unsigned int idx) {
  return const_cast<Bond *>(static_cast<const ROMol *>(this)->getBondWithIdx(
      idx));  // avoid code duplication
}

const Bond *ROMol::getBondWithIdx(unsigned int idx) const {
  URANGE_CHECK(idx, getNumBonds());

  // boost::graph doesn't give us random-access to edges,
  // so we have to iterate to it
  auto [iter, end] = getEdges();
  for (unsigned int i = 0; i < idx; i++) {
    ++iter;
  }
  const Bond *res = d_graph[*iter];

  POSTCONDITION(res != nullptr, "Invalid bond requested");
  return res;
}

Bond *ROMol::getBondBetweenAtoms(unsigned int idx1, unsigned int idx2) {
  return const_cast<Bond *>(
      static_cast<const ROMol *>(this)->getBondBetweenAtoms(
          idx1, idx2));  // avoid code duplication
}

const Bond *ROMol::getBondBetweenAtoms(unsigned int idx1,
                                       unsigned int idx2) const {
  URANGE_CHECK(idx1, getNumAtoms());
  URANGE_CHECK(idx2, getNumAtoms());
  const Bond *res = nullptr;

  auto [edge, found] = boost::edge(boost::vertex(idx1, d_graph),
                                   boost::vertex(idx2, d_graph), d_graph);
  if (found) {
    res = d_graph[edge];
  }
  return res;
}

ROMol::ADJ_ITER_PAIR ROMol::getAtomNeighbors(Atom const *at) const {
  PRECONDITION(at, "no atom");
  PRECONDITION(&at->getOwningMol() == this,
               "atom not associated with this molecule");
  return boost::adjacent_vertices(at->getIdx(), d_graph);
};

ROMol::OBOND_ITER_PAIR ROMol::getAtomBonds(Atom const *at) const {
  PRECONDITION(at, "no atom");
  PRECONDITION(&at->getOwningMol() == this,
               "atom not associated with this molecule");
  return boost::out_edges(at->getIdx(), d_graph);
}

ROMol::ATOM_ITER_PAIR ROMol::getVertices() { return boost::vertices(d_graph); }
ROMol::BOND_ITER_PAIR ROMol::getEdges() { return boost::edges(d_graph); }
ROMol::ATOM_ITER_PAIR ROMol::getVertices() const {
  return boost::vertices(d_graph);
}
ROMol::BOND_ITER_PAIR ROMol::getEdges() const { return boost::edges(d_graph); }

unsigned int ROMol::addAtom(Atom *atom_pin, bool updateLabel,
                            bool takeOwnership) {
  PRECONDITION(atom_pin, "null atom passed in");
  PRECONDITION(!takeOwnership || !atom_pin->hasOwningMol() ||
                   &atom_pin->getOwningMol() == this,
               "cannot take ownership of an atom which already has an owner");
  Atom *atom_p;
  if (!takeOwnership) {
    atom_p = atom_pin->copy();
  } else {
    atom_p = atom_pin;
  }

  atom_p->setOwningMol(this);
  auto which = boost::add_vertex(d_graph);
  d_graph[which] = atom_p;
  atom_p->setIdx(which);
  if (updateLabel) {
    replaceAtomBookmark(atom_p, ci_RIGHTMOST_ATOM);
  }
  for (auto &conf : d_confs) {
    conf->setAtomPos(which, RDGeom::Point3D(0.0, 0.0, 0.0));
  }
  return rdcast<unsigned int>(which);
};

unsigned int ROMol::addBond(Bond *bond_pin, bool takeOwnership) {
  PRECONDITION(bond_pin, "null bond passed in");
  PRECONDITION(!takeOwnership || !bond_pin->hasOwningMol() ||
                   &bond_pin->getOwningMol() == this,
               "cannot take ownership of an bond which already has an owner");
  URANGE_CHECK(bond_pin->getBeginAtomIdx(), getNumAtoms());
  URANGE_CHECK(bond_pin->getEndAtomIdx(), getNumAtoms());
  PRECONDITION(bond_pin->getBeginAtomIdx() != bond_pin->getEndAtomIdx(),
               "attempt to add self-bond");
  PRECONDITION(!(boost::edge(bond_pin->getBeginAtomIdx(),
                             bond_pin->getEndAtomIdx(), d_graph)
                     .second),
               "bond already exists");

  Bond *bond_p;
  if (!takeOwnership) {
    bond_p = bond_pin->copy();
  } else {
    bond_p = bond_pin;
  }

  bond_p->setOwningMol(this);
  auto [which, ok] = boost::add_edge(bond_p->getBeginAtomIdx(),
                                     bond_p->getEndAtomIdx(), d_graph);
  CHECK_INVARIANT(ok, "bond could not be added");
  d_graph[which] = bond_p;
  bond_p->setIdx(numBonds);
  numBonds++;
  return numBonds;
}

void ROMol::setStereoGroups(std::vector<StereoGroup> stereo_groups) {
  auto is_abs = [](const auto &sg) {
    return sg.getGroupType() == StereoGroupType::STEREO_ABSOLUTE;
  };

  // if there's more than one ABS group, merge them
  if (auto num_abs = std::ranges::count_if(stereo_groups, is_abs);
      num_abs <= 1) {
    d_stereo_groups = std::move(stereo_groups);
  } else {
    std::vector<Atom *> abs_atoms;
    std::vector<Bond *> abs_bonds;
    std::vector<StereoGroup> new_stereo_groups;
    new_stereo_groups.reserve(stereo_groups.size() - num_abs + 1);
    for (auto &&sg : stereo_groups) {
      if (is_abs(sg)) {
        auto &other_atoms = sg.getAtoms();
        auto &other_bonds = sg.getBonds();
        abs_atoms.insert(abs_atoms.begin(), other_atoms.begin(),
                         other_atoms.end());
        abs_bonds.insert(abs_bonds.begin(), other_bonds.begin(),
                         other_bonds.end());
      } else {
        new_stereo_groups.push_back(std::move(sg));
      }
    }
    new_stereo_groups.emplace_back(StereoGroupType::STEREO_ABSOLUTE,
                                   std::move(abs_atoms), std::move(abs_bonds));
    d_stereo_groups = std::move(new_stereo_groups);
  }
}

void ROMol::debugMol(std::ostream &str) const {
  str << "Atoms:" << std::endl;
  for (const auto atom : atoms()) {
    str << "\t" << *atom << std::endl;
  }

  str << "Bonds:" << std::endl;
  for (const auto bond : bonds()) {
    str << "\t" << *bond << std::endl;
  }

  const auto &sgs = getSubstanceGroups(*this);
  if (!sgs.empty()) {
    str << "Substance Groups:" << std::endl;
    for (const auto &sg : sgs) {
      str << "\t" << sg << std::endl;
    }
  }

  const auto &stgs = getStereoGroups();
  if (!stgs.empty()) {
    unsigned idx = 0;
    str << "Stereo Groups:" << std::endl;
    for (const auto &stg : stgs) {
      str << "\t" << idx << ' ' << stg << std::endl;
      ++idx;
    }
  }
}

// --------------------------------------------
//
//  Iterators
//
// --------------------------------------------
ROMol::AtomIterator ROMol::beginAtoms() { return AtomIterator(this); }
ROMol::ConstAtomIterator ROMol::beginAtoms() const {
  return ConstAtomIterator(this);
}
ROMol::AtomIterator ROMol::endAtoms() {
  return AtomIterator(this, getNumAtoms());
}
ROMol::ConstAtomIterator ROMol::endAtoms() const {
  return ConstAtomIterator(this, getNumAtoms());
}

ROMol::AromaticAtomIterator ROMol::beginAromaticAtoms() {
  return AromaticAtomIterator(this);
}
ROMol::ConstAromaticAtomIterator ROMol::beginAromaticAtoms() const {
  return ConstAromaticAtomIterator(this);
}
ROMol::AromaticAtomIterator ROMol::endAromaticAtoms() {
  return AromaticAtomIterator(this, getNumAtoms());
}
ROMol::ConstAromaticAtomIterator ROMol::endAromaticAtoms() const {
  return ConstAromaticAtomIterator(this, getNumAtoms());
}

ROMol::HeteroatomIterator ROMol::beginHeteros() {
  return HeteroatomIterator(this);
}
ROMol::ConstHeteroatomIterator ROMol::beginHeteros() const {
  return ConstHeteroatomIterator(this);
}
ROMol::HeteroatomIterator ROMol::endHeteros() {
  return HeteroatomIterator(this, getNumAtoms());
}
ROMol::ConstHeteroatomIterator ROMol::endHeteros() const {
  return ConstHeteroatomIterator(this, getNumAtoms());
}

bool ROMol::hasQuery() const {
  for (auto atom : atoms()) {
    if (atom->hasQuery()) {
      return true;
    }
  }
  for (auto bond : bonds()) {
    if (bond->hasQuery()) {
      return true;
    }
  }
  return false;
}

ROMol::QueryAtomIterator ROMol::beginQueryAtoms(QueryAtom const *what) {
  return QueryAtomIterator(this, what);
}
ROMol::ConstQueryAtomIterator ROMol::beginQueryAtoms(
    QueryAtom const *what) const {
  return ConstQueryAtomIterator(this, what);
}
ROMol::QueryAtomIterator ROMol::endQueryAtoms() {
  return QueryAtomIterator(this, getNumAtoms());
}
ROMol::ConstQueryAtomIterator ROMol::endQueryAtoms() const {
  return ConstQueryAtomIterator(this, getNumAtoms());
}
ROMol::MatchingAtomIterator ROMol::beginMatchingAtoms(bool (*what)(Atom *)) {
  return MatchingAtomIterator(this, what);
}
ROMol::ConstMatchingAtomIterator ROMol::beginMatchingAtoms(
    bool (*what)(const Atom *)) const {
  return ConstMatchingAtomIterator(this, what);
}
ROMol::MatchingAtomIterator ROMol::endMatchingAtoms() {
  return MatchingAtomIterator(this, getNumAtoms());
}
ROMol::ConstMatchingAtomIterator ROMol::endMatchingAtoms() const {
  return ConstMatchingAtomIterator(this, getNumAtoms());
}

ROMol::BondIterator ROMol::beginBonds() { return BondIterator(this); }
ROMol::ConstBondIterator ROMol::beginBonds() const {
  return ConstBondIterator(this);
}
ROMol::BondIterator ROMol::endBonds() {
  auto [beg, end] = getEdges();
  return BondIterator(this, end);
}
ROMol::ConstBondIterator ROMol::endBonds() const {
  auto [beg, end] = getEdges();
  return ConstBondIterator(this, end);
}

void ROMol::clearComputedProps(bool includeRings) const {
  // the SSSR information:
  if (includeRings) {
    this->dp_ringInfo->reset();
  }

  RDProps::clearComputedProps();

  for (auto atom : atoms()) {
    atom->clearComputedProps();
  }

  for (auto bond : bonds()) {
    bond->clearComputedProps();
  }
}

void ROMol::updatePropertyCache(bool strict) {
  for (auto atom : atoms()) {
    atom->updatePropertyCache(strict);
  }
  for (auto bond : bonds()) {
    bond->updatePropertyCache(strict);
  }
}

bool ROMol::needsUpdatePropertyCache() const {
  for (const auto atom : atoms()) {
    if (atom->needsUpdatePropertyCache()) {
      return true;
    }
  }
  // there is no test for bonds yet since they do not obtain a valence property
  return false;
}

void ROMol::clearPropertyCache() {
  for (auto atom : atoms()) {
    atom->clearPropertyCache();
  }
}

const Conformer &ROMol::getConformer(int id) const {
  // make sure we have more than one conformation
  if (d_confs.size() == 0) {
    throw ConformerException("No conformations available on the molecule");
  }

  if (id < 0) {
    return *(d_confs.front());
  }
  auto cid = (unsigned int)id;
  for (auto conf : d_confs) {
    if (conf->getId() == cid) {
      return *conf;
    }
  }
  // we did not find a conformation with the specified ID
  std::string mesg = "Can't find conformation with ID: ";
  mesg += id;
  throw ConformerException(mesg);
}

Conformer &ROMol::getConformer(int id) {
  return const_cast<Conformer &>(
      static_cast<const ROMol *>(this)->getConformer(id));
}

void ROMol::removeConformer(unsigned int id) {
  for (auto ci = d_confs.begin(); ci != d_confs.end(); ++ci) {
    if ((*ci)->getId() == id) {
      d_confs.erase(ci);
      return;
    }
  }
}

unsigned int ROMol::addConformer(Conformer *conf, bool assignId) {
  PRECONDITION(conf, "bad conformer");
  PRECONDITION(conf->getNumAtoms() == this->getNumAtoms(),
               "Number of atom mismatch");
  if (assignId) {
    int maxId = -1;
    for (auto cptr : d_confs) {
      maxId = std::max((int)(cptr->getId()), maxId);
    }
    maxId++;
    conf->setId((unsigned int)maxId);
  }
  conf->setOwningMol(this);
  CONFORMER_SPTR nConf(conf);
  d_confs.push_back(nConf);
  return conf->getId();
}

#ifdef RDK_USE_BOOST_SERIALIZATION
template <class Archive>
void ROMol::save(Archive &ar, const unsigned int) const {
  std::string pkl;
  MolPickler::pickleMol(*this, pkl, PicklerOps::AllProps);
  ar << pkl;
}

template <class Archive>
void ROMol::load(Archive &ar, const unsigned int) {
  std::string pkl;
  ar >> pkl;

  delete dp_ringInfo;
  initMol();

  numBonds = 0;
  MolPickler::molFromPickle(pkl, *this, PicklerOps::AllProps);
  numBonds = rdcast<unsigned int>(boost::num_edges(d_graph));
}

template RDKIT_GRAPHMOL_EXPORT void ROMol::save<boost::archive::text_oarchive>(
    boost::archive::text_oarchive &, const unsigned int) const;
template RDKIT_GRAPHMOL_EXPORT void ROMol::load<boost::archive::text_iarchive>(
    boost::archive::text_iarchive &, const unsigned int);
#endif

}  // namespace RDKit