1 /*
2 * sketcherMinimizerAtom.h
3 *
4 * Created by Nicola Zonta on 19/10/2010.
5 * Copyright Schrodinger, LLC. All rights reserved.
6 *
7 */
8
9 #include "sketcherMinimizerAtom.h"
10 #include "sketcherMinimizer.h"
11 #include "sketcherMinimizerBond.h"
12 #include "sketcherMinimizerMaths.h"
13
14 #include <algorithm>
15 #include <numeric>
16 #include <queue>
17
18 using namespace std;
19
operator <(const CIPAtom & rhs) const20 bool CIPAtom::operator<(const CIPAtom& rhs) const
21 {
22 /* check if this has priority over rhs. An atom is better than another if it
23 has a parent in the chain that gets priority. Parents are evaluated
24 starting from the furthest to the closest.
25 Priority is assigned to the atom with highest atomic number, or has been
26 found to have priority due to its children in previous iterations of the
27 algorithm (scores) or in the present iteration (medals)
28 */
29 assert(allParents.size() == rhs.allParents.size());
30 for (size_t i = 0; i < allParents.size(); i++) {
31
32 if (allParents[i]->atomicNumber > rhs.allParents[i]->atomicNumber) {
33 return true;
34 }
35 if (allParents[i]->atomicNumber < rhs.allParents[i]->atomicNumber) {
36 return false;
37 }
38
39 if ((*scores)[allParents[i]] < (*rhs.scores)[rhs.allParents[i]]) {
40 return true;
41 }
42 if ((*scores)[allParents[i]] > (*rhs.scores)[rhs.allParents[i]]) {
43 return false;
44 }
45
46 vector<int> meds = (*medals)[allParents[i]];
47 vector<int> meds2 = (*rhs.medals)[rhs.allParents[i]];
48 size_t s = (meds.size() < meds2.size()) ? meds.size() : meds2.size();
49
50 for (size_t mm = 0; mm < s; mm++) {
51 if (meds[mm] > meds2[mm]) {
52 return true;
53 }
54 if (meds[mm] < meds2[mm]) {
55 return false;
56 }
57 }
58 if (meds.size() > meds2.size()) {
59 return true;
60 }
61 if (meds2.size() > meds.size()) {
62 return false;
63 }
64 }
65 size_t siz = theseAtoms.size();
66 if (rhs.theseAtoms.size() < siz) {
67 siz = rhs.theseAtoms.size();
68 }
69 for (size_t i = 0; i < siz; i++) {
70 if (theseAtoms[i].first > rhs.theseAtoms[i].first) {
71 return true;
72 }
73 if (theseAtoms[i].first < rhs.theseAtoms[i].first) {
74 return false;
75 }
76 }
77 if (theseAtoms.size() > rhs.theseAtoms.size()) {
78 return true;
79 }
80 if (theseAtoms.size() < rhs.theseAtoms.size()) {
81 return false;
82 }
83
84 return false;
85 }
86
operator ==(const CIPAtom & rhs) const87 bool CIPAtom::operator==(const CIPAtom& rhs) const
88 {
89 assert(allParents.size() == rhs.allParents.size());
90 for (size_t i = 0; i < allParents.size(); i++) {
91 if (allParents[i]->atomicNumber != rhs.allParents[i]->atomicNumber) {
92 return false;
93 }
94 if ((*scores)[allParents[i]] != (*rhs.scores)[rhs.allParents[i]]) {
95 return false;
96 }
97 }
98 if (theseAtoms.size() != rhs.theseAtoms.size()) {
99 return false;
100 }
101 for (size_t i = 0; i < theseAtoms.size(); i++) {
102 if (theseAtoms[i].first != rhs.theseAtoms[i].first) {
103 return false;
104 }
105 }
106 return true;
107 }
108
operator <<(std::ostream & os,const CIPAtom & a)109 std::ostream& operator<<(std::ostream& os, const CIPAtom& a)
110 {
111
112 for (size_t i = 0; i < a.allParents.size(); i++) {
113 os << a.allParents[i]->atomicNumber << "("
114 << (*a.scores)[a.allParents[i]] << ")";
115 if ((*a.medals)[a.allParents[i]].size()) {
116 cerr << "<";
117 for (int ii : (*a.medals)[a.allParents[i]]) {
118 cerr << ii << "|";
119 }
120 cerr << ">";
121 }
122 cerr << " ";
123 }
124 os << "-";
125 for (const auto& theseAtom : a.theseAtoms) {
126 os << " " << theseAtom.first;
127 }
128 return os;
129 }
130
isBetter(CIPAtom & rhs,map<sketcherMinimizerAtom *,unsigned int> * m) const131 bool CIPAtom::isBetter(CIPAtom& rhs,
132 map<sketcherMinimizerAtom*, unsigned int>* m) const
133 {
134 /*
135 Similar to the < operator, but in estabilishing priorities also considers
136 scores stored in m
137 */
138 assert(allParents.size() == rhs.allParents.size());
139 for (size_t i = 0; i < allParents.size(); i++) {
140
141 if ((*m)[allParents[i]] > (*m)[rhs.allParents[i]]) {
142 return true;
143 }
144 if ((*m)[allParents[i]] < (*m)[rhs.allParents[i]]) {
145 return false;
146 }
147
148 if (allParents[i]->atomicNumber > rhs.allParents[i]->atomicNumber) {
149 return true;
150 }
151 if (allParents[i]->atomicNumber < rhs.allParents[i]->atomicNumber) {
152 return false;
153 }
154
155 if ((*scores)[allParents[i]] < (*rhs.scores)[rhs.allParents[i]]) {
156 return true;
157 }
158 if ((*scores)[allParents[i]] > (*rhs.scores)[rhs.allParents[i]]) {
159 return false;
160 }
161
162 vector<int> meds = (*medals)[allParents[i]];
163 vector<int> meds2 = (*rhs.medals)[rhs.allParents[i]];
164 size_t s = (meds.size() < meds2.size()) ? meds.size() : meds2.size();
165
166 for (size_t mm = 0; mm < s; mm++) {
167 if (meds[mm] > meds2[mm]) {
168 return true;
169 }
170 if (meds[mm] < meds2[mm]) {
171 return false;
172 }
173 }
174 if (meds.size() > meds2.size()) {
175 return true;
176 }
177 if (meds2.size() > meds.size()) {
178 return false;
179 }
180 }
181 size_t siz = theseAtoms.size();
182 if (rhs.theseAtoms.size() < siz) {
183 siz = rhs.theseAtoms.size();
184 }
185 for (size_t i = 0; i < siz; i++) {
186 if (theseAtoms[i].first > rhs.theseAtoms[i].first) {
187 return true;
188 }
189 if (theseAtoms[i].first < rhs.theseAtoms[i].first) {
190 return false;
191 }
192 }
193 if (theseAtoms.size() > rhs.theseAtoms.size()) {
194 return true;
195 }
196 if (theseAtoms.size() < rhs.theseAtoms.size()) {
197 return false;
198 }
199
200 return false;
201 }
202
203 sketcherMinimizerAtom::~sketcherMinimizerAtom() = default;
204
sketcherMinimizerAtom()205 sketcherMinimizerAtom::sketcherMinimizerAtom()
206 : crossLayout(false), fixed(false), constrained(false), rigid(false),
207 isSharedAndInner(false), atomicNumber(6), charge(0), _valence(-10),
208 _generalUseN(-1), _generalUseN2(-1), m_chmN(-1),
209 _generalUseVisited(false), _generalUseVisited2(false), fragment(nullptr),
210 needsCheckForClashes(false), visited(false), coordinatesSet(false),
211 isR(true), hasStereochemistrySet(false), _hasRingChirality(false)
212 {
213 hidden = false;
214 m_pseudoZ = 0.f;
215 m_pocketDistance = 0.f;
216 m_x3D = INVALID_COORDINATES;
217 m_y3D = INVALID_COORDINATES;
218 m_z3D = INVALID_COORDINATES;
219 m_isClashing = false;
220 m_isLigand = false;
221 m_isWaterMap = false;
222 m_clockwiseInvert = false;
223 m_isStereogenic = false;
224 m_ignoreRingChirality = false;
225 }
226
227 sketcherMinimizerRing*
shareARing(const sketcherMinimizerAtom * atom1,const sketcherMinimizerAtom * atom2)228 sketcherMinimizerAtom::shareARing(const sketcherMinimizerAtom* atom1,
229 const sketcherMinimizerAtom* atom2)
230 {
231 /* return a ring shared by the two atoms. return a non-macrocycle if
232 * possible */
233 if (!atom1->rings.size()) {
234 return nullptr;
235 }
236 if (!atom2->rings.size()) {
237 return nullptr;
238 }
239
240 foreach (sketcherMinimizerRing* ring, atom1->rings) {
241 if (ring->isMacrocycle()) {
242 continue;
243 }
244 foreach (sketcherMinimizerRing* ring2, atom2->rings) {
245 if (ring == ring2) {
246 return ring;
247 }
248 }
249 }
250 foreach (sketcherMinimizerRing* ring, atom1->rings) {
251 foreach (sketcherMinimizerRing* ring2, atom2->rings) {
252 if (ring == ring2) {
253 return ring;
254 }
255 }
256 }
257 return nullptr;
258 }
259
findHsNumber() const260 int sketcherMinimizerAtom::findHsNumber() const
261 {
262 int valence = _valence;
263 if (valence == -10) {
264 valence = expectedValence(atomicNumber); // valence is not yet set
265 }
266 int nBondOrders = 0;
267 for (auto bond : bonds) {
268 nBondOrders += bond->bondOrder;
269 }
270 if (atomicNumber == 16) { // sulfite & sulfate
271 int nOs = 0;
272 for (size_t i = 0; i < neighbors.size(); ++i) {
273 if (neighbors[i]->atomicNumber == 8 && bonds[i]->bondOrder == 2) {
274 ++nOs;
275 }
276 }
277 if ((nOs) < 3) {
278 valence += nOs * 2;
279 }
280 }
281 if (atomicNumber == 15) { // P
282 int nOs = 0;
283 for (size_t i = 0; i < neighbors.size(); ++i) {
284 if (neighbors[i]->atomicNumber == 8 && bonds[i]->bondOrder == 2) {
285 ++nOs;
286 }
287 }
288 if (nOs < 2) {
289 valence += nOs * 2;
290 }
291 }
292 int out = valence - nBondOrders + charge;
293 if (out < 0) {
294 out = 0;
295 } else if (out > 4) {
296 out = 4;
297 }
298 return out;
299 }
300
301 unsigned int
expectedValence(unsigned int atomicNumber) const302 sketcherMinimizerAtom::expectedValence(unsigned int atomicNumber) const
303 {
304 switch (atomicNumber) {
305 case 1:
306 return 1;
307
308 case 5:
309 return 3;
310
311 case 6:
312 return 4;
313
314 case 7:
315 return 3;
316
317 case 8:
318 return 2;
319
320 case 9:
321 return 1;
322
323 case 14:
324 return 4;
325
326 case 15:
327 return 3;
328
329 case 16:
330 return 2;
331
332 case 17:
333 return 1;
334
335 case 34:
336 return 2;
337 case 35:
338 return 1;
339
340 case 53:
341 return 1;
342
343 default:
344 return 4;
345 }
346 return 4;
347 }
348
349 vector<sketcherMinimizerAtom*>
clockwiseOrderedNeighbors() const350 sketcherMinimizerAtom::clockwiseOrderedNeighbors() const
351 {
352 vector<sketcherMinimizerAtom*> orderedNeighs;
353 vector<sketcherMinimizerAtom*> neighs = neighbors;
354 int lastPoppedIndex = 0;
355 sketcherMinimizerAtom* lastPoppedAtom = neighs[lastPoppedIndex];
356 orderedNeighs.push_back(lastPoppedAtom);
357 neighs.erase(neighs.begin() + lastPoppedIndex);
358
359 while (neighs.size()) {
360 float smallestAngle = 361;
361 for (unsigned int i = 0; i < neighs.size(); i++) {
362 float newAngle = sketcherMinimizerMaths::signedAngle(
363 lastPoppedAtom->coordinates, coordinates,
364 neighs[i]->coordinates);
365 if (newAngle < 0) {
366 newAngle += 360;
367 }
368 if (newAngle < smallestAngle) {
369 smallestAngle = newAngle;
370 lastPoppedIndex = i;
371 }
372 }
373 lastPoppedAtom = neighs[lastPoppedIndex];
374 orderedNeighs.push_back(lastPoppedAtom);
375 neighs.erase(neighs.begin() + lastPoppedIndex);
376 }
377 return orderedNeighs;
378 }
379
writeStereoChemistry()380 void sketcherMinimizerAtom::writeStereoChemistry() // sets stereochemistry for
381 // this atom and from
382 // hasStereochemistrySet and
383 // isR
384 {
385
386 assert(neighbors.size() == bonds.size());
387
388 if (!hasStereochemistrySet) {
389 return;
390 }
391 if (!canBeChiral()) {
392 hasStereochemistrySet = false;
393 return;
394 } else {
395
396 sketcherMinimizerAtom dummyH;
397 dummyH.atomicNumber = 1;
398 dummyH.molecule = molecule;
399
400 sketcherMinimizerAtom dummyLP;
401 dummyLP.atomicNumber = 0;
402 dummyLP.molecule = molecule;
403
404 sketcherMinimizerAtom* dummy = &dummyH;
405
406 vector<sketcherMinimizerAtom*> neighs = neighbors;
407 vector<sketcherMinimizerAtom*> orderedNeighs;
408 vector<sketcherMinimizerBond*> bbonds = bonds;
409 vector<sketcherMinimizerBond*> orderedBonds;
410 vector<float> angles;
411
412 int lastPoppedIndex = 0;
413 sketcherMinimizerAtom* lastPoppedAtom = neighs[lastPoppedIndex];
414 orderedNeighs.push_back(lastPoppedAtom);
415 neighs.erase(neighs.begin() + lastPoppedIndex);
416 orderedBonds.push_back(bbonds[lastPoppedIndex]);
417 bbonds.erase(bbonds.begin() + lastPoppedIndex);
418
419 // TODO: consider using sketcherMinimizerAtom::clockwiseOrderedNeighbors
420 while (neighs.size()) { // order atoms
421 float smallestAngle = 361;
422 for (unsigned int i = 0; i < neighs.size(); i++) {
423 float newAngle = sketcherMinimizerMaths::signedAngle(
424 sketcherMinimizerPointF(lastPoppedAtom->coordinates.x(),
425 lastPoppedAtom->coordinates.y()),
426 sketcherMinimizerPointF(coordinates.x(), coordinates.y()),
427 sketcherMinimizerPointF(neighs[i]->coordinates.x(),
428 neighs[i]->coordinates.y()));
429 if (newAngle < 0) {
430 newAngle += 360;
431 }
432 if (newAngle < smallestAngle) {
433 smallestAngle = newAngle;
434 lastPoppedIndex = i;
435 }
436 }
437 angles.push_back(smallestAngle);
438 lastPoppedAtom = neighs[lastPoppedIndex];
439 orderedNeighs.push_back(lastPoppedAtom);
440 neighs.erase(neighs.begin() + lastPoppedIndex);
441 orderedBonds.push_back(bbonds[lastPoppedIndex]);
442 bbonds.erase(bbonds.begin() + lastPoppedIndex);
443 }
444 if ((atomicNumber == 7 || atomicNumber == 16) && _implicitHs == 0 &&
445 orderedBonds.size() == 3) {
446 dummy = &dummyLP;
447 }
448
449 bool four = true;
450 float totalAngle = std::accumulate(angles.begin(), angles.end(), 0.f);
451 angles.push_back(360.f - totalAngle);
452
453 vector<sketcherMinimizerAtomPriority> atomPriorities,
454 orderedAtomPriorities;
455 for (auto& orderedNeigh : orderedNeighs) {
456 sketcherMinimizerAtomPriority p;
457 p.a = orderedNeigh;
458 atomPriorities.push_back(p);
459 }
460 if (atomPriorities.size() == 3) {
461 four = false;
462 sketcherMinimizerAtomPriority p;
463 p.a = dummy;
464 atomPriorities.push_back(p);
465 }
466
467 bool isStereocenter = setCIPPriorities(atomPriorities, this);
468
469 if (!isStereocenter) {
470 if (!m_ignoreRingChirality) {
471 _hasRingChirality = true;
472 isStereocenter = setCIPPriorities(atomPriorities, this);
473 }
474 }
475 if (!isStereocenter) {
476 _hasRingChirality = false;
477 }
478
479 orderedAtomPriorities = atomPriorities;
480 orderAtomPriorities(orderedAtomPriorities, this);
481
482 if (isStereocenter) {
483 if (!four) {
484 if (orderedAtomPriorities[0].a == dummy) {
485 orderedAtomPriorities.push_back(orderedAtomPriorities[0]);
486 orderedAtomPriorities.insert(orderedAtomPriorities.begin(),
487 orderedAtomPriorities[3]);
488 orderedAtomPriorities.erase(orderedAtomPriorities.begin() +
489 1);
490 orderedAtomPriorities.erase(orderedAtomPriorities.begin() +
491 3);
492 }
493 }
494
495 unsigned int startIndex = 0;
496 for (unsigned int i = 0; i < atomPriorities.size(); i++) {
497 if (atomPriorities[i].a == orderedAtomPriorities[0].a) {
498 startIndex = i;
499 break;
500 }
501 }
502 for (unsigned int i = 0; i < startIndex; i++) {
503 atomPriorities.push_back(atomPriorities[0]);
504 atomPriorities.erase(atomPriorities.begin());
505 }
506
507 if (four) {
508 if (atomPriorities[1].a == orderedAtomPriorities[3].a) {
509 atomPriorities.push_back(atomPriorities[0]);
510 atomPriorities.erase(atomPriorities.begin());
511 }
512 }
513 sketcherMinimizerAtom* mainAtom = nullptr;
514 if (four) {
515 if (atomPriorities[0].a == orderedAtomPriorities[0].a) {
516 mainAtom = atomPriorities[0].a;
517 }
518 if (atomPriorities[3].a == orderedAtomPriorities[0].a) {
519 mainAtom = atomPriorities[3].a;
520 }
521 }
522 bool invert = false;
523 sketcherMinimizerBond* b1 = bondTo(atomPriorities[0].a);
524 if (b1) {
525 if (!four || mainAtom == atomPriorities[0].a ||
526 (!sketcherMinimizer::sameRing(this, atomPriorities[0].a) &&
527 !atomPriorities[0].a->hasStereochemistrySet)) {
528 bool reverse = false;
529 if (b1->startAtom != this) {
530 reverse = true;
531 }
532 b1->isWedge = !invert;
533 b1->hasStereochemistryDisplay = true;
534 b1->isReversed = reverse;
535 } else {
536 b1->hasStereochemistryDisplay = false;
537 }
538 }
539 if (four) {
540 sketcherMinimizerBond* b2 = bondTo(atomPriorities[3].a);
541 if (b2) {
542 if (mainAtom == atomPriorities[3].a ||
543 (!sketcherMinimizer::sameRing(this,
544 atomPriorities[3].a) &&
545 !atomPriorities[3].a->hasStereochemistrySet)) {
546 bool reverse = false;
547 if (b2->startAtom != this) {
548 reverse = true;
549 }
550 b2->isWedge = invert;
551 b2->hasStereochemistryDisplay = true;
552 b2->isReversed = reverse;
553 } else {
554 b2->hasStereochemistryDisplay = false;
555 }
556 }
557 }
558
559 sketcherMinimizerBond* b3 = bondTo(atomPriorities[1].a);
560 if (b3) {
561
562 b3->hasStereochemistryDisplay = false;
563 }
564 sketcherMinimizerBond* b4 = bondTo(atomPriorities[2].a);
565 if (b4) {
566
567 b4->hasStereochemistryDisplay = false;
568 }
569
570 int readS = readStereochemistry(true);
571 readS = -readS; // inverting stereochemistry cause isR is in
572 // sketcher coords and readStereo in coordgen coords
573 // (flipped y)
574
575 if ((readS == -1 && isR) || (readS == 1 && !isR)) { // inverting
576 invert = true;
577
578 if (b1) {
579 if (!four || mainAtom == atomPriorities[0].a ||
580 (!sketcherMinimizer::sameRing(this,
581 atomPriorities[0].a) &&
582 !atomPriorities[0].a->hasStereochemistrySet)) {
583 // cerr << "and setting it "<<endl;
584 bool reverse = false;
585 if (b1->startAtom != this) {
586 reverse = true;
587 }
588 b1->isWedge = !invert;
589 b1->hasStereochemistryDisplay = true;
590 b1->isReversed = reverse;
591 } else {
592 b1->hasStereochemistryDisplay = false;
593 }
594 }
595
596 if (four) {
597 sketcherMinimizerBond* b2 = bondTo(atomPriorities[3].a);
598
599 if (b2) {
600 if (mainAtom == atomPriorities[3].a ||
601 (!sketcherMinimizer::sameRing(
602 this, atomPriorities[3].a) &&
603 !atomPriorities[3].a->hasStereochemistrySet)) {
604 bool reverse = false;
605 if (b2->startAtom != this) {
606 reverse = true;
607 }
608 b2->isWedge = invert;
609 b2->hasStereochemistryDisplay = true;
610 b2->isReversed = reverse;
611 } else {
612 b2->hasStereochemistryDisplay = false;
613 }
614 }
615 }
616 }
617
618 }
619
620 else {
621 for (auto b : bonds) {
622 b->hasStereochemistryDisplay = false;
623 }
624 }
625 }
626 }
627
628 sketcherMinimizerAtomChiralityInfo::sketcherMinimizerChirality
getRelativeStereo(sketcherMinimizerAtom * lookingFrom,sketcherMinimizerAtom * atom1,sketcherMinimizerAtom * atom2)629 sketcherMinimizerAtom::getRelativeStereo(sketcherMinimizerAtom* lookingFrom,
630 sketcherMinimizerAtom* atom1,
631 sketcherMinimizerAtom* atom2)
632 {
633 readStereochemistry(); // to set m_RSPriorities
634 auto RSpriorities = m_RSPriorities;
635 if (RSpriorities.size() < 3) {
636 return sketcherMinimizerAtomChiralityInfo::unspecified;
637 }
638 vector<int> priorities(4, 3);
639
640 /* order the CIP priority of the atoms in the following order
641 atom1 - atom2 - atom3 - atomLookingFrom
642 */
643 for (unsigned int nn = 0; nn < neighbors.size(); nn++) {
644 sketcherMinimizerAtom* n = neighbors[nn];
645 if (n == atom1) {
646 priorities[0] = RSpriorities[nn];
647
648 } else if (n == atom2) {
649 priorities[1] = RSpriorities[nn];
650 } else if (n == lookingFrom) {
651 priorities[3] = RSpriorities[nn];
652 } else {
653 priorities[2] = RSpriorities[nn];
654 }
655 }
656 vector<int> can(4);
657 for (unsigned int i = 0; i < 4; i++) {
658 can[i] = i;
659 }
660 /*
661 this represents a molecule with
662 atom1 (priority 0 - highest)
663 atom2 (priority 1)
664 atom3 (priority 2)
665 atomLookingFrom (priority 3 -lowest)
666 which is the opposite of the CIP rules, where the the lowest priority atom
667 is AWAY from the observer. This is the reason why we return CCW for R and
668 CW for S.
669 */
670 bool match = sketcherMinimizerAtom::matchCIPSequence(priorities, can);
671 bool isClockWise = (match ? !isR : isR);
672 if (isClockWise) {
673 return sketcherMinimizerAtomChiralityInfo::clockwise;
674 }
675 return sketcherMinimizerAtomChiralityInfo::counterClockwise;
676 }
677
setAbsoluteStereoFromChiralityInfo()678 bool sketcherMinimizerAtom::setAbsoluteStereoFromChiralityInfo()
679 {
680 auto info = m_chiralityInfo;
681 if (info.direction == sketcherMinimizerAtomChiralityInfo::unspecified) {
682 return true;
683 }
684 readStereochemistry(); // to set m_RSPriorities
685 auto RSpriorities = m_RSPriorities;
686 if (RSpriorities.size() < 3) {
687 cerr << "CHMMol-> sketcher stereo error: wrong number for RSpriorities"
688 << endl;
689 return false;
690 }
691
692 vector<int> priorities(4, 5);
693
694 bool at3 = false;
695 for (unsigned int nn = 0; nn < neighbors.size(); nn++) {
696 sketcherMinimizerAtom* n = neighbors[nn];
697 if (n == info.atom1) {
698 priorities[0] = RSpriorities[nn];
699
700 } else if (n == info.atom2) {
701
702 priorities[1] = RSpriorities[nn];
703 } else if (n == info.lookingFrom) {
704 priorities[3] = RSpriorities[nn];
705 } else {
706 if (at3) {
707 cerr << "CHMMol-> sketcher stereo error: more than 1 atom not "
708 "matching"
709 << endl;
710 return false;
711
712 } else {
713 at3 = true;
714 priorities[2] = RSpriorities[nn];
715 }
716 }
717 }
718 int addingHN = 0;
719 if (priorities[0] == 5) {
720 priorities[0] = 3;
721 addingHN++;
722 }
723 if (priorities[1] == 5) {
724 priorities[1] = 3;
725 addingHN++;
726 }
727 if (priorities[2] == 5) {
728 priorities[2] = 3;
729 addingHN++;
730 }
731 if (priorities[3] == 5) {
732 priorities[3] = 3;
733 addingHN++;
734 }
735 if (addingHN > 1) {
736 cerr << "CHMMol-> sketcher stereo error: more than 1 H on chiral center"
737 << endl;
738 return false;
739 }
740
741 bool invert = false;
742
743 vector<int> can(4);
744 for (unsigned int i = 0; i < 4; i++) {
745 can[i] = i;
746 }
747 if (!sketcherMinimizerAtom::matchCIPSequence(priorities, can)) {
748 invert = !invert;
749 }
750 bool isRBool = true;
751 if (info.direction == sketcherMinimizerAtomChiralityInfo::clockwise) {
752 isRBool = false;
753 }
754 if (invert) {
755 isRBool = !isRBool;
756 }
757 isR = isRBool;
758 hasStereochemistrySet = true;
759 return true;
760 }
761
matchCIPSequence(vector<int> & v1,vector<int> & v2)762 bool sketcherMinimizerAtom::matchCIPSequence(vector<int>& v1, vector<int>& v2)
763
764 {
765 if (v1.size() < v2.size()) {
766 v1.push_back(3);
767 } else if (v2.size() < v1.size()) {
768 v2.push_back(3);
769 }
770
771 int outofPlaceNs = 0;
772 for (unsigned int i = 0; i < v1.size(); i++) {
773 if (v1[i] != v2[i]) {
774 outofPlaceNs++;
775 }
776 }
777 if (outofPlaceNs == 2) {
778 return false;
779 } else if (outofPlaceNs == 4) {
780 int n1 = v1[0];
781 int index2 = 0;
782 for (unsigned int j = 0; j < v2.size(); j++) {
783 if (v2[j] == n1) {
784 index2 = j;
785 break;
786 }
787 }
788 if (v1[index2] != v2[0]) {
789 return false;
790 }
791 }
792 return true;
793 }
794
setCoordinates(sketcherMinimizerPointF coords)795 void sketcherMinimizerAtom::setCoordinates(sketcherMinimizerPointF coords)
796 {
797 coordinates = std::move(coords);
798 coordinates.round();
799 coordinatesSet = true;
800 }
801
hasNoStereoActiveBonds() const802 bool sketcherMinimizerAtom::hasNoStereoActiveBonds() const
803 {
804 for (auto bond : bonds) {
805 if (bond->isStereo()) {
806 return false;
807 }
808 }
809 return true;
810 }
811
orderAtomPriorities(vector<sketcherMinimizerAtomPriority> & atomPriorities,sketcherMinimizerAtom * center)812 void sketcherMinimizerAtom::orderAtomPriorities(
813 vector<sketcherMinimizerAtomPriority>& atomPriorities,
814 sketcherMinimizerAtom* center) // orders trying to keep long chains in
815 // position 2 and 3 and side substituents in
816 // 1 and 4
817 {
818 assert(atomPriorities.size() == 4);
819 vector<float> weights(4);
820 for (unsigned int i = 0; i < 4; i++) {
821 queue<sketcherMinimizerAtom*> q;
822 for (auto& _atom : center->molecule->_atoms) {
823 _atom->_generalUseVisited = false;
824 }
825
826 q.push(atomPriorities[i].a);
827
828 center->_generalUseVisited = true;
829 atomPriorities[i].a->_generalUseVisited = true;
830 int counter = 0;
831 while (q.size()) {
832 counter++;
833 sketcherMinimizerAtom* at = q.front();
834 q.pop();
835 for (auto n : at->neighbors) {
836 if (!n->_generalUseVisited) {
837 q.push(n);
838 n->_generalUseVisited = true;
839 }
840 }
841 }
842 weights[i] = static_cast<float>(counter);
843 sketcherMinimizerBond* b = center->bondTo(atomPriorities[i].a);
844 if (b) {
845 if (b->bondOrder == 2) {
846 weights[i] -=
847 0.25; // so that =O get lower priority than -OH in phosphate
848 }
849 if (center->atomicNumber == 16 && b->bondOrder == 2) {
850 weights[i] += 2000; // forcing the wedge away from double bond
851 }
852 // in sulphoxide
853
854 if (sketcherMinimizer::sameRing(b->startAtom, b->endAtom)) {
855 weights[i] +=
856 500; // force same ring atoms to be in position 3 and 4
857 }
858 }
859 if (atomPriorities[i].a->atomicNumber == 6) {
860 weights[i] += 0.5; // carbons get priority over other heavy atoms
861 }
862 if (atomPriorities[i].a->atomicNumber == 1) {
863 weights[i] -= 0.5;
864 }
865 if (atomPriorities[i].a->isSharedAndInner &&
866 !center->isSharedAndInner) {
867 weights[i] -= 2000; // forced bond to shared and inner
868 }
869
870 if (center->crossLayout) {
871 if (atomPriorities[i].a->neighbors.size() > 1) {
872 weights[i] += 200;
873 }
874 }
875 if (/* atomPriorities[i].a->isStereogenic && */ atomPriorities[i]
876 .a->hasStereochemistrySet) {
877 weights[i] += 10000; // to avoid problems with wedges when 2
878 }
879 // stereocenters are near
880 for (unsigned int j = 0; j < atomPriorities[i].a->bonds.size(); j++) {
881 if (atomPriorities[i].a->bonds[j]->bondOrder == 2) {
882 weights[i] += 100;
883 break;
884 }
885 }
886 }
887
888 float lowestWeight = weights[0];
889 int index = 0;
890 sketcherMinimizerAtomPriority firstAtom = atomPriorities[0];
891 for (unsigned int i = 1; i < 4; i++) {
892 if (weights[i] < lowestWeight) {
893 lowestWeight = weights[i];
894 index = i;
895 }
896 }
897 firstAtom = atomPriorities[index];
898 atomPriorities.erase(atomPriorities.begin() + index);
899 weights.erase(weights.begin() + index);
900
901 index = 0;
902 lowestWeight = weights[0];
903 sketcherMinimizerAtomPriority secondAtom = atomPriorities[0];
904 for (unsigned int i = 1; i < 3; i++) {
905 if (weights[i] < lowestWeight) {
906 lowestWeight = weights[i];
907 index = i;
908 }
909 }
910
911 secondAtom = atomPriorities[index];
912 atomPriorities.erase(atomPriorities.begin() + index);
913
914 if ((center->atomicNumber != 16 && center->atomicNumber != 15) ||
915 center->neighbors.size() != 4) {
916 atomPriorities.push_back(secondAtom);
917 atomPriorities.insert(atomPriorities.begin(), firstAtom);
918 } else {
919 atomPriorities.insert(atomPriorities.begin() + 1, secondAtom);
920 atomPriorities.insert(atomPriorities.begin(), firstAtom);
921 }
922 }
923
setCIPPriorities(vector<sketcherMinimizerAtomPriority> & atomPriorities,sketcherMinimizerAtom * center)924 bool sketcherMinimizerAtom::setCIPPriorities(
925 vector<sketcherMinimizerAtomPriority>& atomPriorities,
926 sketcherMinimizerAtom* center)
927 {
928 for (auto& atomPrioritie : atomPriorities) {
929 atomPrioritie.priority = 3;
930 }
931 if (atomPriorities.size() != 4) {
932 // cerr << "coordgen: stereo error. (wrong number of atom priorities:
933 // "<< atomPriorities.size () << ")"<<endl; // commented for
934 // Ev:134037
935 return false;
936 }
937 for (unsigned int i = 0; i < atomPriorities.size() - 1; i++) {
938 for (unsigned int j = i + 1; j < atomPriorities.size(); j++) {
939
940 sketcherMinimizerAtom* res =
941 CIPPriority(atomPriorities[i].a, atomPriorities[j].a, center);
942
943 if (res == atomPriorities[i].a) {
944 atomPriorities[i].priority--;
945 } else if (res == atomPriorities[j].a) {
946 atomPriorities[j].priority--;
947 }
948 }
949 }
950 vector<bool> found(4, false);
951
952 for (auto& atomPrioritie : atomPriorities) {
953 if (found[atomPrioritie.priority]) {
954 return false; // two atoms have the same priority
955 }
956 found[atomPrioritie.priority] = true;
957 }
958 return true;
959 }
960
961 sketcherMinimizerAtom*
CIPPriority(sketcherMinimizerAtom * at1,sketcherMinimizerAtom * at2,sketcherMinimizerAtom * center)962 sketcherMinimizerAtom::CIPPriority(sketcherMinimizerAtom* at1,
963 sketcherMinimizerAtom* at2,
964 sketcherMinimizerAtom* center)
965 {
966 assert(center->molecule);
967 assert(at1->molecule == center->molecule);
968 assert(at2->molecule == center->molecule);
969 assert(center->molecule->_atoms.size());
970 assert(at1);
971 assert(at2);
972
973 if (at1->atomicNumber > at2->atomicNumber) {
974 return at1;
975 } else if (at2->atomicNumber > at1->atomicNumber) {
976 return at2;
977 }
978 if (center->_hasRingChirality && !center->m_ignoreRingChirality) {
979 if (sketcherMinimizer::sameRing(center, at1, at2)) {
980 if (at1->_generalUseN > at2->_generalUseN) {
981 return at1;
982 } else {
983 return at2;
984 }
985 }
986 }
987
988 vector<CIPAtom> AN1, AN2;
989
990 map<sketcherMinimizerAtom *, int> score1,
991 score2; // used to keep track if a parent atom has been found to have
992 // priority over another
993 map<sketcherMinimizerAtom *, vector<int>> medals1,
994 medals2; // marks if an atom is a parent of the atoms being evaluated in
995 // the current iteration
996 map<sketcherMinimizerAtom *, int> visited1,
997 visited2; // marks at which iteration this atom was evaluated
998
999 visited1[center] = 1;
1000 visited2[center] = 1;
1001 visited1[at1] = 2;
1002 visited2[at2] = 2;
1003
1004 vector<pair<int, sketcherMinimizerAtom *>> v1, v2;
1005 v1.emplace_back(at1->atomicNumber, at1);
1006 v2.emplace_back(at2->atomicNumber, at2);
1007
1008 vector<sketcherMinimizerAtom *> parents1, parents2;
1009 parents1.push_back(center);
1010 parents1.push_back(at1);
1011 parents2.push_back(center);
1012 parents2.push_back(at2);
1013
1014 AN1.emplace_back(v1, center, parents1, &score1, &medals1, &visited1);
1015 AN2.emplace_back(v2, center, parents2, &score2, &medals2, &visited2);
1016
1017 int level = 1;
1018
1019 while (AN1.size() || AN2.size()) {
1020 level++;
1021
1022 stable_sort(AN1.begin(), AN1.end());
1023 stable_sort(AN2.begin(), AN2.end());
1024
1025 sketcherMinimizerAtom::assignMedals(AN1);
1026 sketcherMinimizerAtom::assignMedals(AN2);
1027
1028 sketcherMinimizerAtom::chooseFirstAndSortAccordingly(AN1);
1029 sketcherMinimizerAtom::chooseFirstAndSortAccordingly(AN2);
1030
1031 sketcherMinimizerAtom::finalizeScores(AN1);
1032 sketcherMinimizerAtom::finalizeScores(AN2);
1033
1034 size_t nn = AN1.size();
1035 if (AN2.size() < nn) {
1036 nn = AN2.size();
1037 }
1038
1039 for (size_t i = 0; i < nn; ++i) {
1040 if (AN1[i] < AN2[i]) {
1041 return at1;
1042 }
1043 if (AN2[i] < AN1[i]) {
1044 return at2;
1045 }
1046 }
1047
1048 if (AN1.size() > AN2.size()) {
1049 return at1;
1050 }
1051 if (AN2.size() > AN1.size()) {
1052 return at2;
1053 }
1054
1055 AN1 = sketcherMinimizerAtom::expandOneLevel(AN1);
1056 AN2 = sketcherMinimizerAtom::expandOneLevel(AN2);
1057 }
1058 return nullptr;
1059 }
1060
chooseFirstAndSortAccordingly(vector<CIPAtom> & V)1061 void sketcherMinimizerAtom::chooseFirstAndSortAccordingly(vector<CIPAtom>& V)
1062 {
1063 if (V.size() < 2) {
1064 return;
1065 }
1066 vector<CIPAtom> copyV = V;
1067 V.clear();
1068 map<sketcherMinimizerAtom*, unsigned int> friendsMask;
1069 while (copyV.size()) {
1070 int bestI = 0;
1071 for (unsigned int i = 1; i < copyV.size(); i++) {
1072 if (copyV[i].isBetter(copyV[bestI], &friendsMask)) {
1073 bestI = i;
1074 }
1075 }
1076 CIPAtom newBest = copyV[bestI];
1077 copyV.erase(copyV.begin() + bestI);
1078 V.push_back(newBest);
1079
1080 for (auto allParent : newBest.allParents) {
1081 friendsMask[allParent] |= (1 << copyV.size());
1082 }
1083 }
1084 }
1085
expandOneLevel(vector<CIPAtom> & oldV)1086 vector<CIPAtom> sketcherMinimizerAtom::expandOneLevel(vector<CIPAtom>& oldV)
1087 {
1088 // we need to keep the bound atoms together, one by one will not work.
1089
1090 map<sketcherMinimizerAtom*, bool> visitedThisRound;
1091 vector<CIPAtom> newV;
1092 for (auto& an : oldV) {
1093 for (unsigned int aa = 0; aa < an.theseAtoms.size(); aa++) {
1094 sketcherMinimizerAtom* a = an.theseAtoms[aa].second;
1095 if (a == nullptr) {
1096 continue; // dummy atom
1097 }
1098 // if (visitedThisRound[a]) continue; // a is present twice
1099 // because closing a ring and has already been dealt with
1100 visitedThisRound[a] = true;
1101 map<sketcherMinimizerAtom*, int>* visited = an.visited;
1102 map<sketcherMinimizerAtom*, vector<int>>* medals = an.medals;
1103 map<sketcherMinimizerAtom*, int>* scores = an.scores;
1104
1105 vector<sketcherMinimizerAtom*> allParents = an.allParents;
1106 allParents.push_back(a);
1107 vector<pair<int, sketcherMinimizerAtom*>> theseAts;
1108
1109 for (unsigned int n = 0; n < a->bonds.size(); n++) {
1110 sketcherMinimizerAtom* neigh = a->neighbors[n];
1111 if (neigh != an.parent) { // if n is not the direct parent of a
1112 bool ghost =
1113 (neigh->atomicNumber == 1 ||
1114 ((*visited)[neigh] &&
1115 (*visited)[neigh] !=
1116 (*visited)[a] + 1) // closing a ring to an atom
1117 // already visited in a
1118 // previous cycle
1119 );
1120 theseAts.emplace_back(
1121 neigh->atomicNumber,
1122 ghost ? ((sketcherMinimizerAtom*) nullptr)
1123 : neigh); // put a ghost for hydrogens and atoms
1124 // closing a ring
1125 if (!ghost) {
1126 (*visited)[neigh] = (*visited)[a] + 1;
1127 }
1128 }
1129 if (a->bonds[n]->bondOrder == 2) { // put ghosts for multiple
1130 // order bonds, even to the
1131 // parent
1132 theseAts.emplace_back(neigh->atomicNumber,
1133 (sketcherMinimizerAtom*) nullptr);
1134 }
1135
1136 else if (a->bonds[n]->bondOrder == 3) {
1137 theseAts.emplace_back(neigh->atomicNumber,
1138 (sketcherMinimizerAtom*) nullptr);
1139 theseAts.emplace_back(neigh->atomicNumber,
1140 (sketcherMinimizerAtom*) nullptr);
1141 }
1142 }
1143
1144 for (int counter = 0; counter < a->_implicitHs;
1145 counter++) { // put ghosts for implicit Hs
1146 theseAts.emplace_back(1, (sketcherMinimizerAtom*) nullptr);
1147 }
1148 stable_sort(theseAts.begin(), theseAts.end());
1149 reverse(theseAts.begin(), theseAts.end());
1150 newV.emplace_back(theseAts, a, allParents, scores, medals, visited);
1151 }
1152 }
1153 return newV;
1154 }
1155
assignMedals(vector<CIPAtom> & v)1156 void sketcherMinimizerAtom::assignMedals(vector<CIPAtom>& v)
1157 {
1158
1159 if (v.size() < 1) {
1160 return;
1161 }
1162 map<sketcherMinimizerAtom*, vector<int>>* medals = v[0].medals;
1163
1164 vector<bool> isEqualToPrevious(v.size());
1165 for (unsigned int i = 1; i < v.size();
1166 i++) { // need to be done before assigning the medals because they are
1167 // considered in the == operator
1168
1169 isEqualToPrevious[i] = (v[i] == v[i - 1]);
1170 }
1171 unsigned int medalLvl = 0;
1172 for (unsigned int i = 0; i < v.size(); i++) {
1173 if (i > 0) {
1174 if (isEqualToPrevious[i]) {
1175 assert(medalLvl > 0);
1176 medalLvl--;
1177 }
1178 }
1179 for (auto allParent : v[i].allParents) {
1180 vector<int> medalsV = (*medals)[allParent];
1181 while (medalsV.size() < medalLvl) {
1182 medalsV.push_back(0);
1183 }
1184 if (medalsV.size() > medalLvl) {
1185 medalsV[medalLvl] = medalsV[medalLvl] + 1;
1186 } else {
1187 medalsV.push_back(1);
1188 }
1189 (*medals)[allParent] = medalsV;
1190 }
1191 medalLvl++;
1192 }
1193 }
1194
finalizeScores(vector<CIPAtom> & v)1195 void sketcherMinimizerAtom::finalizeScores(vector<CIPAtom>& v)
1196 {
1197
1198 if (v.size() < 1) {
1199 return;
1200 }
1201 vector<bool> isEqualToPrevious(v.size());
1202 for (unsigned int i = 1; i < v.size();
1203 i++) { // need to be done before assigning the scores because they are
1204 // considered in the == operator
1205
1206 isEqualToPrevious[i] = (v[i] == v[i - 1]);
1207 }
1208 map<sketcherMinimizerAtom*, vector<int>>* medals = v[0].medals;
1209 map<sketcherMinimizerAtom*, int>* scores = v[0].scores;
1210 scores->clear();
1211 int score = 1;
1212 for (unsigned int i = 0; i < v.size(); i++) {
1213 if (i > 0) {
1214 if (isEqualToPrevious[i]) {
1215 score--;
1216 }
1217 }
1218 /* write the score */
1219 for (auto allParent : v[i].allParents) {
1220 if ((*scores)[allParent] == 0) {
1221 (*scores)[allParent] = score;
1222 }
1223 }
1224 score++;
1225 }
1226 medals->clear();
1227 }
1228
1229 sketcherMinimizerBond*
bondTo(sketcherMinimizerAtom * at) const1230 sketcherMinimizerAtom::bondTo(sketcherMinimizerAtom* at) const
1231 {
1232 for (unsigned int i = 0; i < neighbors.size(); i++) {
1233 if (neighbors[i] == at) {
1234 return bonds[i];
1235 }
1236 }
1237 return nullptr;
1238 }
1239
isResidue() const1240 bool sketcherMinimizerAtom::isResidue() const
1241 {
1242 return false;
1243 }
1244
readStereochemistry(bool readOnly)1245 int sketcherMinimizerAtom::readStereochemistry(
1246 bool readOnly) // 0 if not assigned, 1 if R, -1 if S
1247 {
1248 if (!readOnly) {
1249 _hasRingChirality = false;
1250 m_isStereogenic = false;
1251 }
1252 if (!canBeChiral()) {
1253 return 0;
1254 }
1255 // if (neighbors.size () != 4 && neighbors.size () != 3) return 0;
1256
1257 sketcherMinimizerAtom dummyH;
1258 dummyH.atomicNumber = 1;
1259 dummyH.molecule = molecule;
1260 sketcherMinimizerAtom dummyLP;
1261 dummyLP.atomicNumber = 0;
1262 dummyLP.molecule = molecule;
1263
1264 vector<sketcherMinimizerAtom*> neighs = neighbors;
1265 vector<sketcherMinimizerAtom*> orderedNeighs;
1266 vector<sketcherMinimizerBond*> bnds = bonds;
1267 vector<sketcherMinimizerBond*> orderedBonds;
1268 vector<float> angles;
1269 int lastPoppedIndex = 0;
1270
1271 sketcherMinimizerAtom* lastPoppedAtom = neighs[lastPoppedIndex];
1272 orderedNeighs.push_back(lastPoppedAtom);
1273 neighs.erase(neighs.begin() + lastPoppedIndex);
1274 orderedBonds.push_back(bnds[lastPoppedIndex]);
1275 bnds.erase(bnds.begin() + lastPoppedIndex);
1276
1277 // TODO: consider using sketcherMinimizerAtom::clockwiseOrderedNeighbors
1278 while (neighs.size()) { // order atoms
1279 float smallestAngle = 361;
1280 for (unsigned int i = 0; i < neighs.size(); i++) {
1281 float newAngle = sketcherMinimizerMaths::signedAngle(
1282 lastPoppedAtom->coordinates, coordinates,
1283 neighs[i]->coordinates);
1284 if (newAngle < 0) {
1285 newAngle += 360;
1286 }
1287 if (newAngle < smallestAngle) {
1288 smallestAngle = newAngle;
1289 lastPoppedIndex = i;
1290 }
1291 }
1292 angles.push_back(smallestAngle);
1293 lastPoppedAtom = neighs[lastPoppedIndex];
1294 orderedNeighs.push_back(lastPoppedAtom);
1295 neighs.erase(neighs.begin() + lastPoppedIndex);
1296 orderedBonds.push_back(bnds[lastPoppedIndex]);
1297 bnds.erase(bnds.begin() + lastPoppedIndex);
1298 }
1299
1300 float totalAngle = 0;
1301 for (float angle : angles) {
1302 totalAngle += angle;
1303 }
1304 angles.push_back(360.f - totalAngle);
1305
1306 bool semiplane = false;
1307 sketcherMinimizerAtom* centralAtom = nullptr;
1308 if (angles.size() == 3) {
1309 for (unsigned int i = 0; i < angles.size(); i++) {
1310 if (angles[i] > 180.f) {
1311 semiplane = true;
1312 size_t precI;
1313 if (i == 0) {
1314 precI = angles.size() - 1;
1315 } else {
1316 precI = i - 1;
1317 }
1318 centralAtom = orderedNeighs[precI];
1319 }
1320 }
1321 }
1322
1323 int wedgeN = -1;
1324 int dashedN = -1;
1325 bool giveUp = false;
1326 assert(orderedBonds.size() == 4 || orderedBonds.size() == 3);
1327 for (unsigned int i = 0; i < orderedBonds.size(); i++) {
1328 if (orderedBonds[i]->hasStereochemistryDisplay) {
1329 if (!orderedBonds[i]->isReversed &&
1330 orderedBonds[i]->startAtom != this) {
1331 continue;
1332 }
1333 if (orderedBonds[i]->isReversed &&
1334 orderedBonds[i]->endAtom != this) {
1335 continue;
1336 }
1337
1338 bool wedge = orderedBonds[i]->isWedge;
1339
1340 if (orderedBonds[i]->isReversed) {
1341 wedge = !wedge;
1342 }
1343
1344 if (orderedBonds[i]->startAtom != this) {
1345 wedge = !wedge;
1346 }
1347
1348 if (wedge) {
1349 if (wedgeN == -1) {
1350 wedgeN = i;
1351 } else {
1352 giveUp = true;
1353 }
1354 } else {
1355 if (dashedN == -1) {
1356 dashedN = i;
1357 } else {
1358 giveUp = true;
1359 }
1360 }
1361 }
1362 }
1363
1364 int startIndex = 0;
1365 bool invert = false;
1366 if (dashedN == -1 && wedgeN == -1) {
1367 giveUp = true;
1368 } else if (dashedN == -1) {
1369 startIndex = wedgeN;
1370
1371 } else if (wedgeN == -1) {
1372 startIndex = dashedN;
1373 invert = true;
1374
1375 } else {
1376 if (orderedBonds.size() == 3) {
1377 return 0;
1378 }
1379 if (wedgeN - dashedN == 1 || wedgeN - dashedN == -3) {
1380 startIndex = wedgeN;
1381
1382 } else if (wedgeN - dashedN == -1 || wedgeN - dashedN == 3) {
1383 startIndex = dashedN;
1384 invert = true;
1385 }
1386
1387 else {
1388 giveUp = true;
1389 }
1390 }
1391
1392 int nOfHs = _implicitHs;
1393 size_t totalSubstituentsN = neighbors.size() + nOfHs;
1394
1395 if (orderedNeighs.size() == 3 && nOfHs == 1) {
1396 if (semiplane) {
1397 if (centralAtom == orderedNeighs[startIndex]) {
1398 invert = !invert;
1399 }
1400 }
1401
1402 orderedNeighs.insert(orderedNeighs.begin() + startIndex, &dummyH);
1403 invert = !invert;
1404 }
1405 if (orderedNeighs.size() == 3 &&
1406 (atomicNumber == 7 ||
1407 atomicNumber == 16)) { // nitrogen or sulphoxide // should we check for
1408 // the =O on the sulphoxide?
1409
1410 orderedNeighs.insert(orderedNeighs.begin() + startIndex, &dummyLP);
1411 invert = !invert;
1412 }
1413
1414 vector<sketcherMinimizerAtomPriority> atomPriorities;
1415 for (auto& orderedNeigh : orderedNeighs) {
1416 sketcherMinimizerAtomPriority p;
1417 p.a = orderedNeigh;
1418 atomPriorities.push_back(p);
1419 }
1420 if (atomPriorities.size() != 4) {
1421 return 0;
1422 }
1423
1424 vector<int> canonical;
1425 canonical.push_back(0);
1426 canonical.push_back(1);
1427 canonical.push_back(2);
1428 canonical.push_back(3);
1429
1430 m_RSPriorities.clear();
1431 bool isStereocenter = setCIPPriorities(atomPriorities, this);
1432
1433 if (!isStereocenter) {
1434 if (!readOnly) {
1435 if (!m_ignoreRingChirality) {
1436 _hasRingChirality = true;
1437 isStereocenter = setCIPPriorities(atomPriorities, this);
1438 }
1439 }
1440 }
1441
1442 if (!isStereocenter) {
1443 if (!readOnly) {
1444 _hasRingChirality = false;
1445 }
1446 }
1447
1448 if (!isStereocenter) {
1449 giveUp = true;
1450 } else {
1451 if (!readOnly) {
1452 m_isStereogenic = true;
1453 }
1454 }
1455 if (totalSubstituentsN < 4 && (atomicNumber != 7 && atomicNumber != 16)) {
1456 if (!readOnly) {
1457 m_isStereogenic = false;
1458 }
1459 }
1460
1461 if (!m_isStereogenic) {
1462 if (!readOnly) {
1463 _hasRingChirality = false;
1464 }
1465 }
1466 for (auto n : neighbors) {
1467 for (auto& atomPrioritie : atomPriorities) {
1468 if (atomPrioritie.a == n) {
1469 m_RSPriorities.push_back(atomPrioritie.priority);
1470 break;
1471 }
1472 }
1473 }
1474
1475 if (!matchCIPSequence(canonical, m_RSPriorities)) {
1476 m_clockwiseInvert = true;
1477 } else {
1478 m_clockwiseInvert = false;
1479 }
1480
1481 if (!giveUp) {
1482 int outofPlaceAtoms = 0;
1483 for (unsigned int i = 0; i < atomPriorities.size(); i++) {
1484 // cerr <<i<<" "<<atomPriorities[i].a->atomicNumber<<endl;
1485 int n = startIndex + i;
1486 if (n > 3) {
1487 n -= 4;
1488 }
1489 if (atomPriorities[n].priority != i) {
1490 outofPlaceAtoms++;
1491 }
1492 }
1493 if (outofPlaceAtoms == 2) {
1494 invert = !invert;
1495 } else if (outofPlaceAtoms == 4) {
1496 int n2 = atomPriorities[startIndex].priority + startIndex;
1497 if (n2 > 3) {
1498 n2 -= 4;
1499 }
1500 if (atomPriorities[n2].priority != 0) {
1501 invert = !invert; // if I swapped position 0 with another will
1502 }
1503 // that settle both atoms?
1504 }
1505
1506 if (invert) {
1507 return -1;
1508 } else {
1509 return 1;
1510 }
1511 }
1512 return 0;
1513 }
1514
canBeChiral() const1515 bool sketcherMinimizerAtom::canBeChiral() const
1516 {
1517 if (atomicNumber == 16) {
1518 if (neighbors.size() == 3) {
1519 return true;
1520 }
1521 }
1522 if (atomicNumber == 7) {
1523 if (neighbors.size() == 3 || neighbors.size() == 4) {
1524 return true;
1525 }
1526 }
1527 if (neighbors.size() != 3 && neighbors.size() != 4) {
1528 return false;
1529 }
1530 if ((neighbors.size() + _implicitHs != 4)) {
1531 return false;
1532 }
1533
1534 return true;
1535 }
1536
1537 /* get a vector that points out towards the most free region of space around the
1538 * atom. Used to determined where an arrow to the atom will point from */
getSingleAdditionVector(const vector<sketcherMinimizerAtom * > & ats)1539 sketcherMinimizerPointF sketcherMinimizerAtom::getSingleAdditionVector(
1540 const vector<sketcherMinimizerAtom*>& ats)
1541 {
1542 sketcherMinimizerPointF out(0.f, 0.f);
1543 int n = 0;
1544 for (unsigned int i = 0; i < ats.size(); i++) {
1545 sketcherMinimizerAtom* at = ats[i];
1546 for (unsigned int j = 0; j < at->neighbors.size(); j++) {
1547 sketcherMinimizerAtom* ne = at->neighbors[j];
1548 if (ne->neighbors.size() > 1 &&
1549 find(ats.begin(), ats.end(), ne) == ats.end()) {
1550 out += ne->coordinates - at->coordinates;
1551 n++;
1552 }
1553 }
1554 }
1555 if (n > 0) {
1556 out /= n;
1557 } else {
1558 return sketcherMinimizerPointF(50, 0);
1559 }
1560 out *= -1;
1561 return out;
1562 }
1563
getSingleAdditionVector() const1564 sketcherMinimizerPointF sketcherMinimizerAtom::getSingleAdditionVector() const
1565 {
1566 sketcherMinimizerPointF out(0.f, 0.f);
1567 float totalf = 0.f;
1568 if (neighbors.size()) {
1569 for (auto neighbor : neighbors) {
1570 float f = 1.f;
1571 if (sketcherMinimizer::sameRing(this, neighbor)) {
1572 f = 4.f;
1573 }
1574 out += f * (neighbor->coordinates - coordinates);
1575 totalf += f;
1576 }
1577 out /= totalf;
1578 }
1579 out *= -1;
1580 return out;
1581 }
1582
isMetal(const unsigned int atomicNumber)1583 bool sketcherMinimizerAtom::isMetal(const unsigned int atomicNumber)
1584 {
1585 if (atomicNumber >= 3 && atomicNumber <= 4) {
1586 return true;
1587 }
1588 if (atomicNumber >= 11 && atomicNumber <= 12) {
1589 return true;
1590 }
1591 if (atomicNumber >= 19 && atomicNumber <= 20) {
1592 return true;
1593 }
1594 if (atomicNumber >= 37 && atomicNumber <= 38) {
1595 return true;
1596 }
1597 if (atomicNumber >= 55 && atomicNumber <= 56) {
1598 return true;
1599 }
1600 if (atomicNumber >= 87 && atomicNumber <= 88) {
1601 return true;
1602 }
1603
1604 if (atomicNumber >= 21 && atomicNumber <= 30) {
1605 return true;
1606 }
1607 if (atomicNumber >= 39 && atomicNumber <= 48) {
1608 return true;
1609 }
1610 if (atomicNumber >= 72 && atomicNumber <= 80) {
1611 return true;
1612 }
1613 if (atomicNumber >= 104 && atomicNumber <= 112) {
1614 return true;
1615 }
1616 if (atomicNumber >= 57 && atomicNumber <= 71) {
1617 return true; // lanthanoids
1618 }
1619 if (atomicNumber >= 89 && atomicNumber <= 103) {
1620 return true; // actinoids
1621 }
1622 if (atomicNumber == 13) {
1623 return true;
1624 }
1625 if (atomicNumber == 31) {
1626 return true;
1627 }
1628 if (atomicNumber == 49) {
1629 return true;
1630 }
1631 if (atomicNumber == 81) {
1632 return true;
1633 }
1634 if (atomicNumber == 32) {
1635 return true;
1636 }
1637 if (atomicNumber == 50) {
1638 return true;
1639 }
1640 if (atomicNumber == 82) {
1641 return true;
1642 }
1643 if (atomicNumber == 51) {
1644 return true;
1645 }
1646 if (atomicNumber == 83) {
1647 return true;
1648 }
1649 if (atomicNumber == 84) {
1650 return true;
1651 }
1652 return false;
1653 }
1654
mirrorCoordinates(sketcherMinimizerAtom * at,const sketcherMinimizerBond * bond)1655 void sketcherMinimizerAtom::mirrorCoordinates(sketcherMinimizerAtom* at,
1656 const sketcherMinimizerBond* bond)
1657 {
1658 at->setCoordinates(sketcherMinimizerMaths::mirrorPoint(
1659 at->getCoordinates(), bond->getStartAtom()->getCoordinates(),
1660 bond->getEndAtom()->getCoordinates()));
1661 }
1662
hasValid3DCoordinates() const1663 bool sketcherMinimizerAtom::hasValid3DCoordinates() const
1664 {
1665 return (m_x3D < INVALID_COORDINATES && m_y3D < INVALID_COORDINATES &&
1666 m_z3D < INVALID_COORDINATES);
1667 }
1668
1669 vector<sketcherMinimizerAtom*>
getSubmolecule(sketcherMinimizerAtom * excludedAtom)1670 sketcherMinimizerAtom::getSubmolecule(sketcherMinimizerAtom* excludedAtom)
1671 {
1672 vector<sketcherMinimizerAtom*> subMolecule;
1673 queue<sketcherMinimizerAtom*> q;
1674 map<sketcherMinimizerAtom*, bool> isVisited;
1675 isVisited[excludedAtom] = true;
1676 q.push(this);
1677 isVisited[this] = true;
1678 while (q.size()) {
1679 sketcherMinimizerAtom* atom = q.front();
1680 subMolecule.push_back(atom);
1681 q.pop();
1682 foreach (sketcherMinimizerAtom* neighbor, atom->neighbors) {
1683 if (!isVisited[neighbor]) {
1684 q.push(neighbor);
1685 isVisited[neighbor] = true;
1686 }
1687 }
1688 }
1689 return subMolecule;
1690 }
1691