File: Renumber.cpp

package info (click to toggle)
rdkit 202503.1-5
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid
  • size: 220,160 kB
  • sloc: cpp: 399,240; python: 77,453; ansic: 25,517; java: 8,173; javascript: 4,005; sql: 2,389; yacc: 1,565; lex: 1,263; cs: 1,081; makefile: 580; xml: 229; fortran: 183; sh: 105
file content (134 lines) | stat: -rw-r--r-- 4,324 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
//
//  Copyright (C) 2013 Greg Landrum
//
//   @@ 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.
//
#include <GraphMol/GraphMol.h>
#include <GraphMol/MolOps.h>
#include <RDGeneral/Exceptions.h>
#include <GraphMol/AtomIterators.h>
#include <GraphMol/BondIterators.h>

namespace RDKit {
namespace MolOps {

ROMol *renumberAtoms(const ROMol &mol,
                     const std::vector<unsigned int> &newOrder) {
  unsigned int nAts = mol.getNumAtoms();
  PRECONDITION(newOrder.size() == nAts, "bad newOrder size");

  std::vector<unsigned int> revOrder(nAts);
  for (unsigned int nIdx = 0; nIdx < nAts; ++nIdx) {
    unsigned int oIdx = newOrder[nIdx];
    if (oIdx > nAts) {
      throw ValueErrorException("idx value exceeds numAtoms");
    }
    revOrder[oIdx] = nIdx;
  }

  // ------
  // newOrder[i] : which atom should be in position i of the new mol
  // revOrder[i] : where atom i of the original mol landed in the new mol
  auto *res = new RWMol();

  // copy over the atoms:
  for (unsigned int nIdx = 0; nIdx < nAts; ++nIdx) {
    unsigned int oIdx = newOrder[nIdx];
    const Atom *oAtom = mol.getAtomWithIdx(oIdx);
    Atom *nAtom = oAtom->copy();
    res->addAtom(nAtom, false, true);

    // take care of atom-numbering-dependent properties:
    INT_VECT nAtoms;
    if (nAtom->getPropIfPresent(common_properties::_ringStereoAtoms, nAtoms)) {
      // FIX: ought to be able to avoid this copy.
      for (auto &val : nAtoms) {
        if (val < 0) {
          val = -1 * (revOrder[(-val - 1)] + 1);
        } else {
          val = revOrder[val - 1] + 1;
        }
      }
      nAtom->setProp(common_properties::_ringStereoAtoms, nAtoms, true);
    }

    unsigned int otherAtom;
    if (nAtom->getPropIfPresent(common_properties::_ringStereoOtherAtom,
                                otherAtom)) {
      otherAtom = revOrder[otherAtom];
      nAtom->setProp(common_properties::_ringStereoOtherAtom, otherAtom, true);
    }
  }

  // now the bonds:
  for (ROMol::ConstBondIterator bi = mol.beginBonds(); bi != mol.endBonds();
       ++bi) {
    const Bond *oBond = (*bi);
    Bond *nBond = oBond->copy();
    nBond->setBeginAtomIdx(revOrder[oBond->getBeginAtomIdx()]);
    nBond->setEndAtomIdx(revOrder[oBond->getEndAtomIdx()]);
    res->addBond(nBond, true);
    // take care of atom-numbering-dependent properties:
    for (auto &idx : nBond->getStereoAtoms()) {
      idx = revOrder[idx];
    }
  }

  // Conformers:
  for (auto oConf = mol.beginConformers(); oConf != mol.endConformers();
       ++oConf) {
    auto *nConf = new Conformer(nAts);
    for (unsigned int i = 0; i < nAts; ++i) {
      nConf->setAtomPos(i, (*oConf)->getAtomPos(newOrder[i]));
    }
    nConf->setId((*oConf)->getId());
    nConf->set3D((*oConf)->is3D());
    res->addConformer(nConf);
  }

  // update the ring info:
  const RingInfo *oRings = mol.getRingInfo();
  if (oRings && oRings->isInitialized()) {
    RingInfo *nRings = res->getRingInfo();
    nRings->reset();
    nRings->initialize();
    for (unsigned int i = 0; i < oRings->numRings(); ++i) {
      const INT_VECT &oRing = oRings->atomRings()[i];
      INT_VECT nRing(oRing.size());
      for (unsigned int j = 0; j < oRing.size(); ++j) {
        nRing[j] = revOrder[oRing[j]];
      }
      nRings->addRing(nRing, oRings->bondRings()[i]);
    }
  }

  if (mol.getStereoGroups().size()) {
    std::vector<StereoGroup> nsgs;
    nsgs.reserve(mol.getStereoGroups().size());
    for (const auto &osg : mol.getStereoGroups()) {
      std::vector<Atom *> ats;
      std::vector<Bond *> bds;
      ats.reserve(osg.getAtoms().size());
      bds.reserve(osg.getBonds().size());
      for (const auto aptr : osg.getAtoms()) {
        ats.push_back(res->getAtomWithIdx(revOrder[aptr->getIdx()]));
      }
      for (const auto bptr : osg.getBonds()) {
        bds.push_back(
            res->getBondWithIdx(bptr->getIdx()));  // bonds do not change order
      }

      nsgs.emplace_back(osg.getGroupType(), ats, bds, osg.getReadId());
    }
    res->setStereoGroups(std::move(nsgs));
  }

  return dynamic_cast<ROMol *>(res);
}  // namespace

};  // end of namespace MolOps
};  // end of namespace RDKit