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