File: utils.cpp

package info (click to toggle)
rdkit 202503.6-3
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid
  • size: 222,000 kB
  • sloc: cpp: 411,111; python: 78,482; ansic: 26,181; java: 8,285; javascript: 4,404; sql: 2,393; yacc: 1,626; lex: 1,267; cs: 1,090; makefile: 581; xml: 229; fortran: 183; sh: 121
file content (325 lines) | stat: -rw-r--r-- 10,562 bytes parent folder | download
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
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
#include "utils.h"
#include <GraphMol/MolOps.h>
#include <GraphMol/CIPLabeler/CIPLabeler.h>

#include <GraphMol/SmilesParse/SmilesParse.h>
#include <GraphMol/SmilesParse/SmilesWrite.h>
#include <GraphMol/FileParsers/FileWriters.h>
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 "Forumla";
    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 avgerage
    // 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<unsigned int>(FUSE_LABEL, label);
  return label;
}

void set_fuse_label(Atom *atm, unsigned int idx) {
  // atm->setAtomMapNum(idx); //for debugging
  if (idx) {
    atm->setProp<unsigned int>(FUSE_LABEL, idx);
  } else {
    atm->clearProp(FUSE_LABEL);
  }
}

struct FragmentReplacement {
  // R = Replacement
  // F = Fragment
  // C = Conneciton
  //                    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<Atom *> replacement_connection_atoms;
  std::vector<Atom *> fragment_atoms;

  bool replace(RWMol &mol) {
    if (!replacement_atom) return true;

    auto bond_ordering =
        replacement_atom->getProp<std::vector<int>>(CDX_BOND_ORDERING);

    // Find the connecting atoms and and do the replacement
    for (auto bond : mol.atomBonds(replacement_atom)) {
      // find the position of the attachement bonds in the bond ordering
      auto bond_id = bond->getProp<unsigned int>(CDX_BOND_ID);
      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);

      auto &xatom = fragment_atoms[pos];

      for (auto &xbond : mol.atomBonds(xatom)) {
        // 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<int, FragmentReplacement> 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)) {
    std::vector<int> bond_ordering =
        center_atom->getProp<std::vector<int>>(CDX_BOND_ORDERING);
    if (bond_ordering.size() < 3) {
      return Atom::ChiralType::CHI_UNSPECIFIED;
    }
    std::vector<Atom *> atoms;

    std::vector<std::pair<double, unsigned int>> 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<int>(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<char>(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<std::pair<char, Atom *>> 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<CDXAtomCIPType>(CDX_CIP, cip)) {
        // assign, possibly wrong, initial stereo.
        // note: we can probably deduce this through CDX_BOND_ORDERING, but
        //  I currenlty 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<std::string>(
            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 RDKit