1 /*
2  Contributors: Nicola Zonta
3  Copyright Schrodinger, LLC. All rights reserved
4  */
5 
6 #include "CoordgenMinimizer.h"
7 #include "sketcherMinimizer.h" // should be removed at the end of refactoring
8 #include "sketcherMinimizerAtom.h"
9 #include "sketcherMinimizerBendInteraction.h"
10 #include "sketcherMinimizerBond.h"
11 #include "sketcherMinimizerClashInteraction.h"
12 #include "sketcherMinimizerConstraintInteraction.h"
13 #include "sketcherMinimizerEZConstrainInteraction.h"
14 #include "sketcherMinimizerFragment.h"
15 #include "sketcherMinimizerMaths.h"
16 #include "sketcherMinimizerResidue.h"
17 #include "sketcherMinimizerResidueInteraction.h"
18 #include "sketcherMinimizerRing.h"
19 #include "sketcherMinimizerStretchInteraction.h"
20 #include <algorithm>
21 #include <limits>
22 #include <queue>
23 
24 using namespace std;
25 static const float bondLength = BONDLENGTH;
26 static const float clashEnergyThreshold = 10;
27 
28 #define SAME_SIDE_DPR_PENALTY 100
29 #define SAME_SIDE_DPR_PENALTY_2 50
30 static const float FORCE_MULTIPLIER = 0.3f;
31 
32 static const float STANDARD_CROSSING_BOND_PENALTY = 2500.f;
33 static const float TERMINAL_BOND_CROSSING_MULTIPLIER = 0.5f;
34 static const float MACROCYCLE_BOND_CROSSING_MULTIPLIER = 8.f;
35 static const float RING_BOND_CROSSING_MULTIPLIER = 2.f;
36 
37 static const unsigned int MAXIMUM_NUMBER_OF_SCORED_SOLUTIONS = 10000;
38 static const float REJECTED_SOLUTION_SCORE = 99999999.f;
39 
40 static const unsigned int ITERATION_HISTORY_SIZE = 100;
41 static const float MAX_NET_ENERGY_CHANGE = 20.f;
42 
CoordgenMinimizer()43 CoordgenMinimizer::CoordgenMinimizer()
44 {
45     m_maxIterations = 1000;
46     skipMinimization = false;
47     skipFlipFragments = false;
48     skipAvoidClashes = false;
49     m_scoreResidueInteractions = true;
50     m_precision = 1.f;
51     energy_list = {};
52     all_coordinates = {};
53 }
54 
~CoordgenMinimizer()55 CoordgenMinimizer::~CoordgenMinimizer()
56 {
57     clearInteractions();
58 }
59 
clearInteractions()60 void CoordgenMinimizer::clearInteractions()
61 {
62     for (auto& _interaction : _interactions) {
63         delete _interaction;
64     }
65     _interactions.clear();
66     _intramolecularClashInteractions.clear();
67     _extraInteractions.clear();
68     _stretchInteractions.clear();
69     _bendInteractions.clear();
70 }
71 
run()72 void CoordgenMinimizer::run()
73 {
74     if (skipMinimization) {
75         return;
76     }
77     if (!_interactions.size()) {
78         setupInteractions();
79     }
80 
81 #ifdef DEBUG_MINIMIZATION_COORDINATES
82     // to seperate energy and DOF minimization
83     energy_list.push_back(-1.f);
84     all_coordinates.push_back({sketcherMinimizerPointF(0,0)});
85 #endif
86 
87     std::vector<float> local_energy_list(m_maxIterations);
88     std::vector<sketcherMinimizerPointF> lowest_energy_coords(m_atoms.size());
89     float min_energy = std::numeric_limits<float>::max();
90     for (unsigned int iterations = 0; iterations < m_maxIterations; ++iterations) {
91         local_energy_list[iterations] = scoreInteractions();
92         // track coordinates with lowest energy
93         if (local_energy_list[iterations] < min_energy) {
94             for (size_t i = 0; i < m_atoms.size(); ++i) {
95                 lowest_energy_coords[i] = m_atoms[i]->coordinates;
96             }
97         }
98 #ifdef DEBUG_MINIMIZATION_COORDINATES
99         // store data from this minimization step to be written to a file later
100         energy_list.push_back(local_energy_list[iterations]);
101         std::vector<sketcherMinimizerPointF> these_coordinates;
102         for (auto atom : m_atoms) {
103             these_coordinates.push_back(atom->coordinates);
104         }
105         all_coordinates.push_back(these_coordinates);
106 #endif
107         if (!applyForces(0.1f)) {
108             break;
109         }
110         if (iterations < 2 * ITERATION_HISTORY_SIZE) {
111             continue;
112         }
113         if (local_energy_list[iterations - ITERATION_HISTORY_SIZE] - local_energy_list[iterations] < MAX_NET_ENERGY_CHANGE) {
114             break;
115         }
116     }
117     // set coordinates back to lowest energy state
118     if (min_energy < std::numeric_limits<float>::max()) {
119         for (size_t i = 0; i < m_atoms.size(); ++i) {
120             m_atoms[i]->coordinates = lowest_energy_coords[i];
121         }
122     }
123 }
124 
applyForces(float maxd)125 bool CoordgenMinimizer::applyForces(float maxd)
126 {
127     float delta = 0.001f; // minimum squared displacement
128     float distance = 0.f;
129     for (auto atom : m_atoms) {
130         if (atom->fixed) {
131             continue;
132         }
133         sketcherMinimizerPointF displacement = atom->force * FORCE_MULTIPLIER;
134         if (displacement.x() != displacement.x() ||
135             displacement.y() != displacement.y()) {
136             displacement = sketcherMinimizerPointF(0.f, 0.f);
137         }
138         float dsquare = displacement.x() * displacement.x() +
139                         displacement.y() * displacement.y();
140         if (dsquare < SKETCHER_EPSILON) {
141             dsquare = SKETCHER_EPSILON;
142         }
143         if (dsquare > maxd * maxd) {
144             displacement *= maxd / sqrt(dsquare);
145         }
146         atom->coordinates += displacement;
147         distance += displacement.squareLength();
148         atom->force = sketcherMinimizerPointF(0, 0);
149     }
150     return distance >= delta;
151 }
152 
153 /* store extra interaction to be used when minimizing molecule.
154  cis amides constraints are an example as they need 3d coordinates
155  to be detected */
addExtraInteraction(sketcherMinimizerMolecule * molecule,sketcherMinimizerInteraction * interaction)156 void CoordgenMinimizer::addExtraInteraction(
157     sketcherMinimizerMolecule* molecule,
158     sketcherMinimizerInteraction* interaction)
159 {
160     _extraInteractionsOfMolecule[molecule].push_back(interaction);
161 }
162 
addClashInteractionsOfMolecule(sketcherMinimizerMolecule * molecule,bool intrafragmentClashes)163 void CoordgenMinimizer::addClashInteractionsOfMolecule(
164     sketcherMinimizerMolecule* molecule, bool intrafragmentClashes)
165 {
166     vector<sketcherMinimizerAtom*> atoms = molecule->getAtoms();
167     vector<sketcherMinimizerBond*> bonds = molecule->getBonds();
168 
169     if (atoms.size() > 1) {
170         for (sketcherMinimizerAtom* atom : atoms) {
171             if (atom->isResidue()) {
172                 continue;
173             }
174 
175             for (sketcherMinimizerBond* bond : bonds) {
176                 if (bond->isResidueInteraction()) {
177                     continue;
178                 }
179                 sketcherMinimizerAtom* at2 = atom;
180                 sketcherMinimizerAtom* at1 = bond->startAtom;
181                 sketcherMinimizerAtom* at3 = bond->endAtom;
182                 if (at1 == at2 || at1 == at3 || at2 == at3) {
183                     continue;
184                 }
185                 if (at1->fragment->getDofsOfAtom(at1).empty() &&
186                     at2->fragment->getDofsOfAtom(at2).empty() &&
187                     at3->fragment->getDofsOfAtom(at3).empty() &&
188                     !intrafragmentClashes) {
189                     if (at1->fragment == at2->fragment) {
190                         continue;
191                     }
192                     if (at3->fragment == at2->fragment) {
193                         continue;
194                     }
195                 }
196                 if (at2->fixed && at1->fixed && at3->fixed) {
197                     continue;
198                 }
199 
200                 if (at1->isNeighborOf(at2)) {
201                     continue;
202                 }
203                 for (sketcherMinimizerAtom* n : at1->neighbors) {
204                     if (n->isNeighborOf(at2)) {
205                         continue;
206                     }
207                 }
208                 if (at3->isNeighborOf(at2)) {
209                     continue;
210                 }
211                 for (sketcherMinimizerAtom* n : at3->neighbors) {
212                     if (n->isNeighborOf(at2)) {
213                         continue;
214                     }
215                 }
216 
217                 if (!(at1->rigid && at2->rigid && at3->rigid)) {
218 
219                     //                }
220                     auto* interaction =
221                         new sketcherMinimizerClashInteraction(at1, at2, at3);
222                     float restVK = 0.8f;
223                     if (at2->atomicNumber == 6 && at2->charge == 0) {
224                         restVK -= 0.1f;
225                     }
226                     if (at1->atomicNumber == 6 && at1->charge == 0 &&
227                         at3->atomicNumber == 6 && at3->charge == 0) {
228                         restVK -= 0.1f;
229                     }
230 
231                     interaction->restV =
232                         (bondLength * restVK) * (bondLength * restVK);
233                     _intramolecularClashInteractions.push_back(interaction);
234                     _interactions.push_back(interaction);
235                 }
236             }
237         }
238     }
239 }
240 
addStretchInteractionsOfMolecule(sketcherMinimizerMolecule * molecule)241 void CoordgenMinimizer::addStretchInteractionsOfMolecule(
242                                 sketcherMinimizerMolecule* molecule)
243 {
244     vector<sketcherMinimizerBond*> bonds = molecule->getBonds();
245     for (sketcherMinimizerBond* bo : bonds) {
246         if (bo->isResidueInteraction()) {
247             continue;
248         }
249         sketcherMinimizerAtom* at1 = bo->startAtom;
250         sketcherMinimizerAtom* at2 = bo->endAtom;
251         auto* interaction = new sketcherMinimizerStretchInteraction(at1, at2);
252         interaction->k *= 0.1f;
253         interaction->restV = bondLength;
254         if (at1->rigid && at2->rigid) {
255             sketcherMinimizerPointF v = at2->coordinates - at1->coordinates;
256             interaction->restV = v.length();
257         }
258         auto sharedRing = sketcherMinimizer::sameRing(at1, at2);
259         if (sharedRing && !sharedRing->isMacrocycle()) {
260             interaction->k *= 50;
261         }
262         _interactions.push_back(interaction);
263         _stretchInteractions.push_back(interaction);
264     }
265 }
266 
267 /* return a set of all carbons part of a carbonyl group, i.e. doubly bonded to
268  * an oxygen. */
getChetoCs(const std::vector<sketcherMinimizerAtom * > & allAtoms) const269 std::set<sketcherMinimizerAtom*> CoordgenMinimizer::getChetoCs(
270     const std::vector<sketcherMinimizerAtom*>& allAtoms) const
271 {
272     std::set<sketcherMinimizerAtom*> chetoCs;
273     for (auto atom : allAtoms) {
274         if (atom->atomicNumber != 6) {
275             continue;
276         }
277         for (auto bondedAtom : atom->neighbors) {
278             if (bondedAtom->atomicNumber == 8) {
279                 auto bond = sketcherMinimizer::getBond(atom, bondedAtom);
280                 if (bond && bond->bondOrder == 2) {
281                     chetoCs.insert(atom);
282                     continue;
283                 }
284             }
285         }
286     }
287     return chetoCs;
288 }
289 
290 /* return a set of all amino nitrogens. Not chemically accurate, doesn't filter
291  * out nitro Ns for instance. */
getAminoNs(const std::vector<sketcherMinimizerAtom * > & allAtoms) const292 std::set<sketcherMinimizerAtom*> CoordgenMinimizer::getAminoNs(
293     const std::vector<sketcherMinimizerAtom*>& allAtoms) const
294 {
295     std::set<sketcherMinimizerAtom*> aminoNs;
296     for (auto atom : allAtoms) {
297         if (atom->atomicNumber == 7) {
298             aminoNs.insert(atom);
299         }
300     }
301     return aminoNs;
302 }
303 
304 /* return a set of all aminoacid alpha carbon, i.e. a carbon that is bound to a
305  nitrogen
306  and a cheto carbon. */
getAlphaCs(const std::vector<sketcherMinimizerAtom * > & allAtoms,const std::set<sketcherMinimizerAtom * > & chetoCs,const std::set<sketcherMinimizerAtom * > & aminoNs) const307 std::set<sketcherMinimizerAtom*> CoordgenMinimizer::getAlphaCs(
308     const std::vector<sketcherMinimizerAtom*>& allAtoms,
309     const std::set<sketcherMinimizerAtom*>& chetoCs,
310     const std::set<sketcherMinimizerAtom*>& aminoNs) const
311 {
312     std::set<sketcherMinimizerAtom*> alphaCs;
313     for (auto atom : allAtoms) {
314         bool bondedToCheto = false;
315         bool bondedToAminoN = false;
316         if (atom->atomicNumber != 6) {
317             continue;
318         }
319         if (chetoCs.find(atom) != chetoCs.end()) {
320             continue;
321         }
322         for (auto bondedAtom : atom->neighbors) {
323             if (chetoCs.find(bondedAtom) != chetoCs.end()) {
324                 bondedToCheto = true;
325             }
326             if (aminoNs.find(bondedAtom) != aminoNs.end()) {
327                 bondedToAminoN = true;
328             }
329         }
330         if (bondedToCheto && bondedToAminoN) {
331             alphaCs.insert(atom);
332         }
333     }
334     return alphaCs;
335 }
336 
337 /* add constraints to keep all backbone atoms of a peptide in a straight trans
338  * line. */
addPeptideBondInversionConstraintsOfMolecule(sketcherMinimizerMolecule * molecule)339 void CoordgenMinimizer::addPeptideBondInversionConstraintsOfMolecule(
340     sketcherMinimizerMolecule* molecule)
341 {
342     auto atoms = molecule->getAtoms();
343     auto chetoCs = getChetoCs(atoms);
344     if (chetoCs.size() < 2) {
345         return;
346     }
347     auto aminoNs = getAminoNs(atoms);
348     if (aminoNs.size() < 2) {
349         return;
350     }
351     auto alphaCs = getAlphaCs(atoms, chetoCs, aminoNs);
352     if (alphaCs.size() < 2) {
353         return;
354     }
355     std::vector<std::vector<sketcherMinimizerAtom*>> consecutiveAtomsGroups;
356     getFourConsecutiveAtomsThatMatchSequence(consecutiveAtomsGroups, chetoCs,
357                                              aminoNs, alphaCs, chetoCs);
358     getFourConsecutiveAtomsThatMatchSequence(consecutiveAtomsGroups, aminoNs,
359                                              alphaCs, chetoCs, aminoNs);
360     getFourConsecutiveAtomsThatMatchSequence(consecutiveAtomsGroups, alphaCs,
361                                              chetoCs, aminoNs, alphaCs);
362     for (auto torsionAtoms : consecutiveAtomsGroups) {
363         bool cis = false;
364         auto interaction = new sketcherMinimizerEZConstrainInteraction(
365             torsionAtoms[0], torsionAtoms[1], torsionAtoms[2], torsionAtoms[3],
366             cis);
367         _extraInteractions.push_back(interaction);
368         _interactions.push_back(interaction);
369     }
370 }
371 
372 /* find chains of four bound atoms that are part of the four provided sets.
373  Useful to detect portions of a peptide backbone for instance. */
getFourConsecutiveAtomsThatMatchSequence(std::vector<std::vector<sketcherMinimizerAtom * >> & consecutiveAtomsGroups,const std::set<sketcherMinimizerAtom * > & firstSet,const std::set<sketcherMinimizerAtom * > & secondSet,const std::set<sketcherMinimizerAtom * > & thirdSet,const std::set<sketcherMinimizerAtom * > & fourthSet) const374 void CoordgenMinimizer::getFourConsecutiveAtomsThatMatchSequence(
375     std::vector<std::vector<sketcherMinimizerAtom*>>& consecutiveAtomsGroups,
376     const std::set<sketcherMinimizerAtom*>& firstSet,
377     const std::set<sketcherMinimizerAtom*>& secondSet,
378     const std::set<sketcherMinimizerAtom*>& thirdSet,
379     const std::set<sketcherMinimizerAtom*>& fourthSet) const
380 {
381     for (auto firstAtom : firstSet) {
382         for (auto secondAtom : firstAtom->neighbors) {
383             if (secondSet.find(secondAtom) == secondSet.end()) {
384                 continue;
385             }
386             for (auto thirdAtom : secondAtom->neighbors) {
387                 if (thirdSet.find(thirdAtom) == thirdSet.end()) {
388                     continue;
389                 }
390                 for (auto fourthAtom : thirdAtom->neighbors) {
391                     if (fourthSet.find(fourthAtom) == fourthSet.end()) {
392                         continue;
393                     }
394                     std::vector<sketcherMinimizerAtom*> fourMatchingAtoms(4);
395                     fourMatchingAtoms.at(0) = firstAtom;
396                     fourMatchingAtoms.at(1) = secondAtom;
397                     fourMatchingAtoms.at(2) = thirdAtom;
398                     fourMatchingAtoms.at(3) = fourthAtom;
399                     consecutiveAtomsGroups.push_back(fourMatchingAtoms);
400                 }
401             }
402         }
403     }
404 }
405 
addConstrainedInteractionsOfMolecule(sketcherMinimizerMolecule * molecule)406 void CoordgenMinimizer::addConstrainedInteractionsOfMolecule(
407     sketcherMinimizerMolecule* molecule)
408 {
409     for (auto atom : molecule->getAtoms()) {
410         if (atom->constrained) {
411             auto interaction = new sketcherMinimizerConstraintInteraction(
412                 atom, atom->templateCoordinates);
413             _intramolecularClashInteractions.push_back(interaction);
414             _interactions.push_back(interaction);
415         }
416     }
417 }
418 
addChiralInversionConstraintsOfMolecule(sketcherMinimizerMolecule * molecule)419 void CoordgenMinimizer::addChiralInversionConstraintsOfMolecule(
420     sketcherMinimizerMolecule* molecule)
421 {
422     for (auto ring : molecule->getRings()) {
423         if (ring->isMacrocycle()) {
424             vector<sketcherMinimizerAtom*> atoms =
425                 CoordgenFragmentBuilder::orderRingAtoms(ring);
426             for (unsigned int i = 0; i < atoms.size(); i++) {
427                 int size = static_cast<int>(atoms.size());
428                 int a1 = (i - 1 + size) % size;
429                 int a11 = (i - 2 + size) % size;
430                 int a2 = (i + 1) % size;
431 
432                 sketcherMinimizerBond* bond =
433                     sketcherMinimizer::getBond(atoms[a1], atoms[i]);
434                 if (bond->isStereo()) {
435                     bool cis = bond->markedAsCis(atoms[a11], atoms[a2]);
436                     auto* ezint = new sketcherMinimizerEZConstrainInteraction(
437                         atoms[a11], atoms[a1], atoms[i], atoms[a2], cis);
438                     _interactions.push_back(ezint);
439                 }
440             }
441         }
442     }
443 }
444 
addBendInteractionsOfMolecule(sketcherMinimizerMolecule * molecule)445 void CoordgenMinimizer::addBendInteractionsOfMolecule(
446     sketcherMinimizerMolecule* molecule)
447 {
448     vector<sketcherMinimizerAtom*> atoms = molecule->getAtoms();
449     vector<sketcherMinimizerBond*> bonds = molecule->getBonds();
450     for (sketcherMinimizerAtom* at : atoms) {
451         vector<sketcherMinimizerBendInteraction*> interactions;
452         vector<sketcherMinimizerBendInteraction*> ringInteractions;
453         vector<sketcherMinimizerBendInteraction*> nonRingInteractions;
454         int nbonds = static_cast<int>(at->neighbors.size());
455         bool invertedMacrocycleBond = false;
456         if (nbonds > 1) {
457             // order bonds so that they appear in clockwise order.
458             vector<sketcherMinimizerAtom*> orderedNeighs =
459                 at->clockwiseOrderedNeighbors();
460             float angle = 120.f;
461             if (nbonds == 2) {
462                 if (at->bonds[0]->bondOrder + at->bonds[1]->bondOrder > 3) {
463                     angle = 180.f;
464                 }
465             }
466             if (nbonds > 2) {
467                 angle = 360.f / nbonds;
468             }
469             for (int i = 0; i < nbonds; i++) {
470                 int j = (i - 1 + nbonds) % nbonds;
471                 if (nbonds == 2 && i == 1) {
472                     continue; // first and last interaction are the same if
473                 }
474                 // there is just one interaction
475                 sketcherMinimizerAtom* at1 = orderedNeighs[i];
476                 sketcherMinimizerAtom* at2 = at;
477                 sketcherMinimizerAtom* at3 = orderedNeighs[j];
478                 auto* interaction =
479                     new sketcherMinimizerBendInteraction(at1, at2, at3);
480                 interactions.push_back(interaction);
481                 interaction->restV = angle;
482                 sketcherMinimizerRing* r = sketcherMinimizer::sameRing(
483                     at, orderedNeighs[i], orderedNeighs[j]);
484                 if (r) {
485                     if (!r->isMacrocycle()) {
486                         int extraAtoms = 0;
487                         /* if the rings are to be drawn as fused, they will
488                          * result in a bigger ring */
489                         for (unsigned int i = 0; i < r->fusedWith.size(); i++) {
490                             if (r->fusedWith[i]->isMacrocycle()) {
491                                 continue;
492                             }
493                             if (r->fusionAtoms[i].size() > 2) {
494                                 extraAtoms += static_cast<int>(
495                                     r->fusedWith[i]->_atoms.size() -
496                                     r->fusionAtoms[i].size());
497                             }
498                         }
499                         interaction->isRing = true;
500                         interaction->k *= 100;
501                         interaction->restV = static_cast<float>(
502                             180. - (360. / (r->size() + extraAtoms)));
503                         ringInteractions.push_back(interaction);
504                     } else {
505                         if (nbonds == 3) {
506                             sketcherMinimizerAtom* otherAtom = nullptr;
507                             for (auto atom : orderedNeighs) {
508                                 if (atom != at1 && atom != at3) {
509                                     otherAtom = atom;
510                                     break;
511                                 }
512                             }
513                             if (otherAtom) {
514                                 if (sketcherMinimizerMaths::sameSide(
515                                         at3->coordinates,
516                                         otherAtom->coordinates,
517                                         at1->coordinates, at2->coordinates)) {
518                                     invertedMacrocycleBond = true;
519                                 }
520                             }
521                         }
522                         bool fusedToRing = false;
523                         if (orderedNeighs.size() > 2) {
524                             fusedToRing = true;
525                         }
526                         for (auto atom : orderedNeighs) {
527                             if (!sketcherMinimizer::sameRing(at, atom)) {
528                                 fusedToRing = false;
529                                 break;
530                             }
531                         }
532                         if (fusedToRing || invertedMacrocycleBond) {
533                             ringInteractions.push_back(interaction);
534                         } else {
535                             /* macrocycles are treated as non rings */
536                             nonRingInteractions.push_back(interaction);
537                         }
538                     }
539                 } else {
540                     nonRingInteractions.push_back(interaction);
541                 }
542                 if (interaction->atom1->rigid && interaction->atom2->rigid &&
543                     interaction->atom3->rigid) {
544                     interaction->restV = sketcherMinimizerMaths::unsignedAngle(
545                         interaction->atom1->coordinates,
546                         interaction->atom2->coordinates,
547                         interaction->atom3->coordinates);
548                 }
549             }
550         }
551         if (ringInteractions.size() != 1 || nonRingInteractions.size() != 2) {
552             invertedMacrocycleBond = false;
553         }
554         if (ringInteractions.size()) { // subtract all the ring angles from 360
555                                        // and divide the remaining equally
556                                        // between the other interactions
557             float totalAngleInRings = 0;
558             for (sketcherMinimizerBendInteraction* i : ringInteractions) {
559                 totalAngleInRings += i->restV;
560             }
561             if (invertedMacrocycleBond) {
562                 totalAngleInRings = 360 - totalAngleInRings;
563             }
564             for (sketcherMinimizerBendInteraction* i : nonRingInteractions) {
565                 i->restV =
566                     (360 - totalAngleInRings) / nonRingInteractions.size();
567             }
568         } else { // do nothing if 1 or 3 interactions (defaults to 120 degrees)
569                  // if 4 or more set angles accordingly
570 
571             if (nonRingInteractions.size() == 4) {
572                 if (at->crossLayout || m_evenAngles) {
573                     nonRingInteractions[0]->restV = 90;
574                     nonRingInteractions[1]->restV = 90;
575                     nonRingInteractions[2]->restV = 90;
576                     nonRingInteractions[3]->restV = 90;
577                 } else {
578                     int indexOfBiggestAngle = 0;
579                     float biggestAngle = 0;
580                     int counter = 0;
581                     for (auto interaction : nonRingInteractions) {
582                         float angle = sketcherMinimizerMaths::unsignedAngle(
583                             interaction->atom1->coordinates,
584                             interaction->atom2->coordinates,
585                             interaction->atom3->coordinates);
586                         if (angle > biggestAngle) {
587                             biggestAngle = angle;
588                             indexOfBiggestAngle = counter;
589                         }
590                         counter++;
591                     }
592                     nonRingInteractions[indexOfBiggestAngle]->restV = 120;
593                     nonRingInteractions[(indexOfBiggestAngle + 1) % 4]->restV =
594                         90;
595                     nonRingInteractions[(indexOfBiggestAngle + 2) % 4]->restV =
596                         60;
597                     nonRingInteractions[(indexOfBiggestAngle + 3) % 4]->restV =
598                         90;
599                 }
600             } else if (nonRingInteractions.size() > 4) {
601                 for (sketcherMinimizerBendInteraction* i :
602                      nonRingInteractions) {
603                     i->restV =
604                         static_cast<float>(360. / nonRingInteractions.size());
605                 }
606             }
607         }
608         for (auto interaction : interactions) {
609             if (!(interaction->atom1->fixed && interaction->atom2->fixed &&
610                   interaction->atom3->fixed)) {
611                 _interactions.push_back(interaction);
612                 _bendInteractions.push_back(interaction);
613             } else {
614                 delete interaction;
615             }
616         }
617     }
618 }
619 
minimizeMolecule(sketcherMinimizerMolecule * molecule)620 void CoordgenMinimizer::minimizeMolecule(sketcherMinimizerMolecule* molecule)
621 {
622     std::map<sketcherMinimizerAtom*, sketcherMinimizerPointF>
623         previousCoordinates;
624     for (auto atom : molecule->getAtoms()) {
625         previousCoordinates[atom] = atom->getCoordinates();
626     }
627     clearInteractions();
628     addInteractionsOfMolecule(molecule, true);
629     run();
630     for (auto bond : molecule->getBonds()) {
631         if (!bond->checkStereoChemistry()) {
632             for (auto atom : molecule->getAtoms()) {
633                 atom->setCoordinates(previousCoordinates[atom]);
634             }
635             break;
636         }
637     }
638 }
639 
minimizeResidues()640 void CoordgenMinimizer::minimizeResidues()
641 {
642     setupInteractionsOnlyResidues();
643     run();
644 }
645 
minimizeProteinOnlyLID(const std::map<std::string,std::vector<sketcherMinimizerResidue * >> & chains)646 void CoordgenMinimizer::minimizeProteinOnlyLID(
647     const std::map<std::string, std::vector<sketcherMinimizerResidue*>>& chains)
648 {
649     setupInteractionsProteinOnly(chains);
650     run();
651 }
652 
minimizeAll()653 void CoordgenMinimizer::minimizeAll()
654 {
655     setupInteractions(true);
656     run();
657 }
658 
addInteractionsOfMolecule(sketcherMinimizerMolecule * molecule,bool intrafragmentClashes)659 void CoordgenMinimizer::addInteractionsOfMolecule(
660                                                   sketcherMinimizerMolecule* molecule, bool intrafragmentClashes)
661 {
662     addClashInteractionsOfMolecule(molecule, intrafragmentClashes);
663     addStretchInteractionsOfMolecule(molecule);
664     addBendInteractionsOfMolecule(molecule);
665     addChiralInversionConstraintsOfMolecule(molecule);
666 }
667 
setupInteractionsOnlyResidues()668 void CoordgenMinimizer::setupInteractionsOnlyResidues()
669 {
670     const float CLASH_DISTANCE = bondLength * 1.5f;
671     for (auto res : m_residues) {
672         for (auto res2 : m_residues) {
673             if (res2 >= res) {
674                 continue;
675             }
676             auto* minimizerInteraction =
677                 new sketcherMinimizerClashInteraction(res, res2, res);
678             minimizerInteraction->restV = CLASH_DISTANCE * CLASH_DISTANCE;
679             _interactions.push_back(minimizerInteraction);
680         }
681     }
682 }
683 
setupInteractionsProteinOnly(const std::map<std::string,std::vector<sketcherMinimizerResidue * >> & chains)684 void CoordgenMinimizer::setupInteractionsProteinOnly(
685     const std::map<std::string, std::vector<sketcherMinimizerResidue*>>& chains)
686 {
687     clearInteractions();
688     std::set<sketcherMinimizerBond*> interactions;
689     std::set<sketcherMinimizerResidue*> residues;
690     for (const auto& chain : chains) {
691         for (auto res : chain.second) {
692             residues.insert(res);
693             for (auto interaction : res->residueInteractions) {
694                 interactions.insert(interaction);
695             }
696         }
697     }
698     for (auto res : residues) {
699         for (auto interaction : interactions) {
700             if (res == interaction->startAtom || res == interaction->endAtom) {
701                 continue;
702             }
703             auto* minimizerInteraction = new sketcherMinimizerClashInteraction(
704                 interaction->startAtom, res, interaction->endAtom);
705             minimizerInteraction->restV = bondLength * bondLength;
706             _interactions.push_back(minimizerInteraction);
707         }
708     }
709 }
710 
setupInteractions(bool intrafragmentClashes)711 void CoordgenMinimizer::setupInteractions(bool intrafragmentClashes)
712 {
713     clearInteractions();
714     for (sketcherMinimizerMolecule* molecule : m_molecules) {
715         addInteractionsOfMolecule(molecule, intrafragmentClashes);
716     }
717 }
718 
scoreInteractions()719 float CoordgenMinimizer::scoreInteractions()
720 {
721     float totalEnergy = 0.f;
722     for (auto interaction : _interactions) {
723         interaction->score(totalEnergy);
724     }
725     return totalEnergy;
726 }
727 
728 // returns true if the two molecules have a atom-atoms or atom-bond pair with
729 // distance < threshold or crossing bonds
findIntermolecularClashes(sketcherMinimizerMolecule * mol1,sketcherMinimizerMolecule * mol2,float threshold)730 bool CoordgenMinimizer::findIntermolecularClashes(
731     sketcherMinimizerMolecule* mol1, sketcherMinimizerMolecule* mol2,
732     float threshold)
733 {
734     // could be made faster for instance checking the molecules bounding boxes
735     // first
736     if (mol1 == mol2) {
737         return false;
738     }
739     float threshold2 = threshold * threshold;
740     for (sketcherMinimizerAtom* a : mol1->_atoms) {
741         for (sketcherMinimizerAtom* a2 : mol2->_atoms) {
742             if (sketcherMinimizerMaths::squaredDistance(
743                     a->coordinates, a2->coordinates) < threshold2) {
744                 return true;
745             }
746         }
747     }
748 
749     for (sketcherMinimizerAtom* a : mol1->_atoms) {
750         for (sketcherMinimizerBond* b : mol2->_bonds) {
751             if (sketcherMinimizerMaths::squaredDistancePointSegment(
752                     a->coordinates, b->startAtom->coordinates,
753                     b->endAtom->coordinates) < threshold2) {
754                 return true;
755             }
756         }
757     }
758     for (sketcherMinimizerAtom* a : mol2->_atoms) {
759         for (sketcherMinimizerBond* b : mol1->_bonds) {
760             if (sketcherMinimizerMaths::squaredDistancePointSegment(
761                     a->coordinates, b->startAtom->coordinates,
762                     b->endAtom->coordinates) < threshold2) {
763                 return true;
764             }
765         }
766     }
767     for (sketcherMinimizerBond* b : mol1->_bonds) {
768         for (sketcherMinimizerBond* b2 : mol2->_bonds) {
769             if (sketcherMinimizerMaths::intersectionOfSegments(
770                     b->startAtom->coordinates, b->endAtom->coordinates,
771                     b2->startAtom->coordinates, b2->endAtom->coordinates)) {
772                 return true;
773             }
774         }
775     }
776 
777     return false;
778 }
779 
findIntermolecularClashes(const vector<sketcherMinimizerMolecule * > & mols,float threshold)780 bool CoordgenMinimizer::findIntermolecularClashes(
781     const vector<sketcherMinimizerMolecule*>& mols, float threshold)
782 {
783     for (unsigned int i = 0; i < mols.size(); i++) {
784         for (unsigned int j = i + 1; j < mols.size(); j++) {
785             if (findIntermolecularClashes(mols[i], mols[j], threshold)) {
786                 return true;
787             }
788         }
789     }
790     return false;
791 }
792 
scoreClashes(sketcherMinimizerMolecule * molecule,bool residueInteractions,bool scoreProximityRelationsOnOppositeSid) const793 float CoordgenMinimizer::scoreClashes(
794     sketcherMinimizerMolecule* molecule, bool residueInteractions,
795     bool scoreProximityRelationsOnOppositeSid) const
796 {
797     float E = 0.f;
798     for (auto i : _intramolecularClashInteractions) {
799         i->score(E, true);
800     }
801     for (sketcherMinimizerInteraction* i : _extraInteractions) {
802         i->score(E, true);
803     }
804     E += scoreDofs(molecule);
805     E += scoreCrossBonds(molecule, residueInteractions);
806     E += scoreAtomsInsideRings();
807     if (scoreProximityRelationsOnOppositeSid) {
808         E += scoreProximityRelationsOnOppositeSides();
809     }
810     return E;
811 }
812 
scoreDofs(sketcherMinimizerMolecule * molecule) const813 float CoordgenMinimizer::scoreDofs(sketcherMinimizerMolecule* molecule) const
814 {
815     float E = 0.f;
816     for (const auto& fragment : molecule->getFragments()) {
817         for (const auto& dof : fragment->getDofs()) {
818             E += dof->getCurrentPenalty();
819         }
820     }
821     return E;
822 }
823 
scoreCrossBonds(sketcherMinimizerMolecule * molecule,bool residueInteractions) const824 float CoordgenMinimizer::scoreCrossBonds(sketcherMinimizerMolecule* molecule,
825                                          bool residueInteractions) const
826 {
827     if (!m_scoreResidueInteractions) {
828         residueInteractions = false;
829     }
830 
831     float out = 0.f;
832     vector<sketcherMinimizerBond*> bonds = molecule->getBonds();
833     if (molecule->getBonds().size() > 2) {
834 
835         for (unsigned int b = 0; b < bonds.size() - 1; b++) {
836             sketcherMinimizerBond* b1 = bonds[b];
837             if (b1->isResidueInteraction()) {
838                 continue;
839             }
840             for (unsigned int bb = b + 1; bb < bonds.size(); bb++) {
841                 sketcherMinimizerBond* b2 = bonds[bb];
842                 if (b2->isResidueInteraction()) {
843                     continue;
844                 }
845                 if (b2->startAtom->molecule != b1->startAtom->molecule) {
846                     continue;
847                 }
848                 if (bondsClash(b1, b2)) {
849                     float penalty = STANDARD_CROSSING_BOND_PENALTY *
850                         b1->crossingBondPenaltyMultiplier *
851                         b2->crossingBondPenaltyMultiplier;
852                     if (b1->isTerminal() || b2->isTerminal()) {
853                         penalty *= TERMINAL_BOND_CROSSING_MULTIPLIER;
854                     }
855                     if (b1->isInMacrocycle() || b2->isInMacrocycle()) {
856                         penalty *= MACROCYCLE_BOND_CROSSING_MULTIPLIER;
857                     }
858                     if (b1->isInSmallRing() || b2->isInSmallRing()) {
859                         penalty *= RING_BOND_CROSSING_MULTIPLIER;
860                     }
861                     out += penalty;
862                 }
863             }
864         }
865     }
866     if (m_residueInteractions.size() && residueInteractions) {
867         for (auto r : m_residues) {
868             if (r->residueInteractions.size() > 1) {
869                 for (unsigned int ri1 = 0;
870                      ri1 < r->residueInteractions.size() - 1; ri1++) {
871                     for (unsigned int ri2 = 1;
872                          ri2 < r->residueInteractions.size(); ri2++) {
873 
874                         sketcherMinimizerAtom* a1 =
875                             r->residueInteractions[ri1]->endAtom;
876                         sketcherMinimizerAtom* a2 =
877                             r->residueInteractions[ri2]->endAtom;
878                         if (sketcherMinimizerMaths::intersectionOfSegments(
879                                 a1->coordinates +
880                                     a1->getSingleAdditionVector() * 0.2f,
881                                 a2->coordinates +
882                                     a2->getSingleAdditionVector() * 0.2f,
883                                 a1->coordinates, a2->coordinates)) {
884                             out += 15.f;
885                         }
886 
887                         for (auto b2 : m_bonds) {
888                             if (b2->startAtom ==
889                                 r->residueInteractions[ri1]->endAtom) {
890                                 continue;
891                             }
892                             if (b2->endAtom ==
893                                 r->residueInteractions[ri1]->endAtom) {
894                                 continue;
895                             }
896                             if (b2->startAtom ==
897                                 r->residueInteractions[ri2]->endAtom) {
898                                 continue;
899                             }
900                             if (b2->endAtom ==
901                                 r->residueInteractions[ri2]->endAtom) {
902                                 continue;
903                             }
904 
905                             if (sketcherMinimizerMaths::intersectionOfSegments(
906                                     a1->coordinates, a2->coordinates,
907                                     b2->startAtom->coordinates,
908                                     b2->endAtom->coordinates)) {
909                                 out += 10.f;
910                             }
911                         }
912                     }
913                 }
914             }
915         }
916     }
917     return out;
918 }
919 
scoreAtomsInsideRings() const920 float CoordgenMinimizer::scoreAtomsInsideRings() const
921 {
922     float out = 0.f;
923     float cutOff = bondLength;
924     for (sketcherMinimizerMolecule* m : m_molecules) {
925         for (sketcherMinimizerRing* r : m->_rings) {
926             if (r->_atoms.size() > MACROCYCLE) {
927                 continue;
928             }
929             if (r->_atoms.size() < 3) {
930                 continue;
931             }
932             sketcherMinimizerPointF c = r->findCenter();
933             for (sketcherMinimizerAtom* a : m->_atoms) {
934                 if (a->fragment == r->_atoms[0]->fragment) {
935                     continue;
936                 }
937                 sketcherMinimizerPointF d = c - a->coordinates;
938                 if (d.x() > cutOff) {
939                     continue;
940                 }
941                 if (d.y() > cutOff) {
942                     continue;
943                 }
944                 if (d.x() < -cutOff) {
945                     continue;
946                 }
947                 if (d.y() < -cutOff) {
948                     continue;
949                 }
950                 float sq = d.squareLength();
951                 if (sq > cutOff * cutOff) {
952                     continue;
953                 }
954                 float dist = d.length();
955                 if (dist < cutOff) {
956                     out += 50 + 100 * (1 - (dist / cutOff));
957                 }
958             }
959         }
960     }
961     return out;
962 }
963 
scoreProximityRelationsOnOppositeSides() const964 float CoordgenMinimizer::scoreProximityRelationsOnOppositeSides() const
965 {
966     float out = 0.f;
967     for (sketcherMinimizerMolecule* m : m_molecules) {
968         if (m->_atoms.size() < 2) {
969             continue;
970         }
971         for (unsigned int i = 0; i < m->m_proximityRelations.size(); i++) {
972             sketcherMinimizerPointF v1, v2;
973             sketcherMinimizerMolecule* otherMol1 = nullptr;
974             sketcherMinimizerBond* pr1 = m->m_proximityRelations[i];
975             sketcherMinimizerFragment* f1 = nullptr;
976             if (pr1->startAtom->molecule == m) {
977                 f1 = pr1->startAtom->fragment;
978                 v1 = pr1->startAtom->getSingleAdditionVector();
979                 otherMol1 = pr1->endAtom->molecule;
980             } else {
981                 f1 = pr1->endAtom->fragment;
982                 v1 = pr1->endAtom->getSingleAdditionVector();
983                 otherMol1 = pr1->startAtom->molecule;
984             }
985             if (otherMol1 == m) {
986                 continue;
987             }
988 
989             for (unsigned int j = i + 1; j < m->m_proximityRelations.size();
990                  j++) {
991                 sketcherMinimizerMolecule* otherMol2 = nullptr;
992 
993                 sketcherMinimizerBond* pr2 = m->m_proximityRelations[j];
994                 if (pr2->startAtom->molecule == m) {
995                     if (pr2->startAtom->fragment == f1) {
996                         continue;
997                     }
998                     v2 = pr2->startAtom->getSingleAdditionVector();
999                     otherMol2 = pr2->endAtom->molecule;
1000 
1001                 } else {
1002                     if (pr2->endAtom->fragment == f1) {
1003                         continue;
1004                     }
1005                     v2 = pr2->endAtom->getSingleAdditionVector();
1006                     otherMol2 = pr2->startAtom->molecule;
1007                 }
1008                 if (otherMol2 == m) {
1009                     continue;
1010                 }
1011                 if (otherMol1 != otherMol2) {
1012                     continue;
1013                 }
1014                 float angle = sketcherMinimizerMaths::unsignedAngle(
1015                     v1, sketcherMinimizerPointF(0.f, 0.f), v2);
1016                 if (angle > 90) {
1017                     out += SAME_SIDE_DPR_PENALTY +
1018                            SAME_SIDE_DPR_PENALTY_2 * (angle - 90);
1019                 }
1020             }
1021         }
1022     }
1023     return out;
1024 }
1025 
runExhaustiveSearch(sketcherMinimizerMolecule * molecule,vector<CoordgenFragmentDOF * > dofs,float & clashE,CoordgenDOFSolutions & solutions)1026 bool CoordgenMinimizer::runExhaustiveSearch(sketcherMinimizerMolecule* molecule,
1027                                             vector<CoordgenFragmentDOF*> dofs,
1028                                             float& clashE,
1029                                             CoordgenDOFSolutions& solutions)
1030 {
1031     float bestResult = clashE;
1032     bool abort = false;
1033     runExhaustiveSearchLevel(molecule, dofs.begin(), dofs, bestResult, abort,
1034                              solutions);
1035     for (auto dof : dofs) {
1036         dof->setToOptimalValue();
1037     }
1038     clashE = bestResult;
1039     return (bestResult < clashEnergyThreshold);
1040 }
1041 
runExhaustiveSearchLevel(sketcherMinimizerMolecule * molecule,const vector<CoordgenFragmentDOF * >::iterator & iterator,vector<CoordgenFragmentDOF * > & dofs,float & bestResult,bool & abort,CoordgenDOFSolutions & solutions)1042 void CoordgenMinimizer::runExhaustiveSearchLevel(
1043     sketcherMinimizerMolecule* molecule,
1044     const vector<CoordgenFragmentDOF*>::iterator& iterator,
1045     vector<CoordgenFragmentDOF*>& dofs, float& bestResult, bool& abort,
1046     CoordgenDOFSolutions& solutions)
1047 {
1048     if (abort) {
1049         return;
1050     }
1051     if (iterator == dofs.end()) {
1052         float result = solutions.scoreCurrentSolution();
1053         if (result < clashEnergyThreshold) {
1054             for (auto dof : dofs) {
1055                 dof->storeCurrentValueAsOptimal();
1056             }
1057             abort = true;
1058         } else if (result < bestResult - SKETCHER_EPSILON) {
1059             bestResult = result;
1060             for (auto dof : dofs) {
1061                 dof->storeCurrentValueAsOptimal();
1062             }
1063         }
1064     } else {
1065         vector<CoordgenFragmentDOF*>::iterator nextIter = iterator;
1066         ++nextIter;
1067         for (int i = 0; i < (*iterator)->numberOfStates(); ++i) {
1068             runExhaustiveSearchLevel(molecule, nextIter, dofs, bestResult,
1069                                      abort, solutions);
1070             (*iterator)->changeState();
1071         }
1072     }
1073 }
1074 
1075 std::vector<std::vector<CoordgenFragmentDOF*>>
buildTuplesOfDofs(const vector<CoordgenFragmentDOF * > & dofs,unsigned int order) const1076 CoordgenMinimizer::buildTuplesOfDofs(const vector<CoordgenFragmentDOF*>& dofs,
1077                                      unsigned int order) const
1078 {
1079     std::vector<std::vector<CoordgenFragmentDOF *>> growingVector,
1080         lastOrderVector;
1081     for (auto dof : dofs) {
1082         std::vector<CoordgenFragmentDOF*> tuple;
1083         tuple.push_back(dof);
1084         growingVector.push_back(tuple);
1085     }
1086     for (unsigned int i = 1; i < order; ++i) {
1087         lastOrderVector = growingVector;
1088         growingVector.clear();
1089         for (auto lastOrderTuple : lastOrderVector) {
1090             bool copy = false;
1091             for (auto dof : dofs) {
1092                 if (copy) {
1093                     auto newTuple = lastOrderTuple;
1094                     newTuple.push_back(dof);
1095                     growingVector.push_back(newTuple);
1096                 } else if (dof == *(lastOrderTuple.rbegin())) {
1097                     copy = true;
1098                 }
1099             }
1100         }
1101     }
1102     return growingVector;
1103 }
1104 
growSolutions(std::set<std::vector<short unsigned int>> & allScoredSolutions,int & currentTier,std::map<std::vector<short unsigned int>,float> & growingSolutions,CoordgenDOFSolutions & solutions,float & bestScore)1105 bool CoordgenMinimizer::growSolutions(
1106     std::set<std::vector<short unsigned int>>& allScoredSolutions,
1107     int& currentTier,
1108     std::map<std::vector<short unsigned int>, float>& growingSolutions,
1109     CoordgenDOFSolutions& solutions, float& bestScore)
1110 {
1111     std::map<std::vector<short unsigned int>, float> oldGrowingSolutions =
1112         growingSolutions;
1113     float bestScoreForRun = bestScore;
1114     std::vector<std::pair<float, std::vector<short unsigned int>>>
1115         bestSolutions;
1116     bestSolutions.reserve(growingSolutions.size());
1117     for (const auto& solution : growingSolutions) {
1118         bestSolutions.emplace_back(solution.second, solution.first);
1119     }
1120     sort(bestSolutions.begin(), bestSolutions.end());
1121     growingSolutions.clear();
1122     int maxN = static_cast<int>(6 * getPrecision());
1123     if (maxN < 1) {
1124         maxN = 1;
1125     }
1126     int n = 0;
1127 
1128     for (const auto& solution : bestSolutions) {
1129         if (n > maxN) {
1130             break;
1131         }
1132         for (auto dof : solutions.getAllDofs()) {
1133             if (dof->tier() > currentTier) {
1134                 continue;
1135             }
1136             solutions.loadSolution(solution.second);
1137             for (int i = 1; i < dof->numberOfStates(); ++i) {
1138                 dof->changeState();
1139 
1140                 auto newSolution = solutions.getCurrentSolution();
1141                 if (allScoredSolutions.find(newSolution) ==
1142                     allScoredSolutions.end()) {
1143                     float score = solutions.scoreCurrentSolution();
1144                     if (score == REJECTED_SOLUTION_SCORE) {
1145                         return false;
1146                     }
1147                     allScoredSolutions.insert(newSolution);
1148                     if (score < bestScore) {
1149                         bestScore = score;
1150                     }
1151                     if (score < bestScoreForRun &&
1152                         score < REJECTED_SOLUTION_SCORE) {
1153                         growingSolutions[newSolution] = score;
1154                     }
1155                 }
1156             }
1157         }
1158         n++;
1159     }
1160     if (growingSolutions.empty() && currentTier < 5) {
1161         currentTier += 3;
1162         growingSolutions = oldGrowingSolutions;
1163     }
1164     return true;
1165 }
1166 
runSearch(int tier,CoordgenDOFSolutions & solutions)1167 bool CoordgenMinimizer::runSearch(int tier, CoordgenDOFSolutions& solutions)
1168 {
1169     std::map<std::vector<short unsigned int>, float> growingSolutions;
1170     std::set<std::vector<short unsigned int>> allScoredSolutions;
1171     float bestScore = solutions.scoreCurrentSolution();
1172     growingSolutions[solutions.getCurrentSolution()] = bestScore;
1173     int i = 0;
1174     bool hasValidSolution = true;
1175     do {
1176 #ifdef DEBUG_MINIMIZATION_COORDINATES
1177         // store data from this minimization step to be written to a file later
1178         energy_list.push_back(solutions.scoreCurrentSolution());
1179         std::vector<sketcherMinimizerPointF> these_coordinates;
1180         for (auto atom : _atoms) {
1181             these_coordinates.push_back(atom->coordinates);
1182         }
1183         all_coordinates.push_back(these_coordinates);
1184 #endif
1185         ++i;
1186         hasValidSolution = growSolutions(
1187             allScoredSolutions, tier, growingSolutions, solutions, bestScore);
1188     } while ((hasValidSolution && !growingSolutions.empty()) && i < 100);
1189     std::pair<std::vector<short unsigned int>, float> bestSolution =
1190         solutions.findBestSolution();
1191     solutions.loadSolution(bestSolution.first);
1192     return bestSolution.second < clashEnergyThreshold;
1193 }
1194 
runLocalSearch(sketcherMinimizerMolecule * molecule,const vector<CoordgenFragmentDOF * > & dofs,int levels,float & clashE,CoordgenDOFSolutions & solutions)1195 bool CoordgenMinimizer::runLocalSearch(sketcherMinimizerMolecule* molecule,
1196                                        const vector<CoordgenFragmentDOF*>& dofs,
1197                                        int levels, float& clashE,
1198                                        CoordgenDOFSolutions& solutions)
1199 {
1200     bool downhill = false;
1201     auto combinationsOfDofs = buildTuplesOfDofs(dofs, levels);
1202     do {
1203         downhill = false;
1204         for (const auto& combinationOfDofs : combinationsOfDofs) {
1205             float lastResult = clashE;
1206             bool foundOptimalPosition = runExhaustiveSearch(
1207                 molecule, combinationOfDofs, clashE, solutions);
1208             if (foundOptimalPosition) {
1209                 return true;
1210             } else if (clashE < lastResult - SKETCHER_EPSILON) {
1211                 downhill = true;
1212             }
1213         }
1214     } while (downhill);
1215     return false;
1216 }
1217 
flipFragments(sketcherMinimizerMolecule * molecule,float & clashE)1218 bool CoordgenMinimizer::flipFragments(sketcherMinimizerMolecule* molecule,
1219                                       float& clashE)
1220 {
1221     float bestResult = clashE;
1222     if (skipFlipFragments) {
1223         return true;
1224     }
1225     if (bestResult < clashEnergyThreshold) {
1226         return true;
1227     }
1228     vector<CoordgenFragmentDOF *> dofs, onlyFlipDofs;
1229     vector<sketcherMinimizerFragment*> fragments = molecule->getFragments();
1230     reverse(fragments.begin(), fragments.end());
1231     for (auto fragment : fragments) {
1232         if (!fragment->fixed) {
1233             for (auto dof : fragment->getDofs()) {
1234                 if (dof->numberOfStates() > 1) {
1235                     dofs.push_back(dof);
1236                     if (dof == *(fragment->getDofs().begin())) {
1237                         onlyFlipDofs.push_back(dof);
1238                     }
1239                 }
1240             }
1241         }
1242     }
1243     CoordgenDOFSolutions solutions(this, molecule, dofs);
1244     bool cleanPose = runSearch(0, solutions);
1245     //  if (!cleanPose) cleanPose = runSearch(6, solutions);
1246     buildMoleculeFromFragments(molecule, false);
1247     return cleanPose;
1248 }
1249 
avoidClashesOfMolecule(sketcherMinimizerMolecule * molecule,const std::vector<sketcherMinimizerInteraction * > & extraInteractions)1250 bool CoordgenMinimizer::avoidClashesOfMolecule(
1251     sketcherMinimizerMolecule* molecule,
1252     const std::vector<sketcherMinimizerInteraction*>& extraInteractions)
1253 {
1254     clearInteractions();
1255     addClashInteractionsOfMolecule(molecule, false);
1256     addPeptideBondInversionConstraintsOfMolecule(molecule);
1257     for (sketcherMinimizerInteraction* interaction : extraInteractions) {
1258         _interactions.push_back(interaction);
1259         _extraInteractions.push_back(interaction);
1260     }
1261     for (auto interaction : _extraInteractionsOfMolecule[molecule]) {
1262         _extraInteractions.push_back(interaction);
1263         _interactions.push_back(interaction);
1264     }
1265     bool scoreResidueInteractions = true;
1266     bool doNotComputeForces = true;
1267     float clashE =
1268         scoreClashes(molecule, scoreResidueInteractions, doNotComputeForces);
1269     bool cleanPose = flipFragments(molecule, clashE);
1270     if (!cleanPose) {
1271         avoidTerminalClashes(molecule, clashE);
1272         molecule->requireMinimization();
1273     }
1274     if (molecule->minimizationIsRequired()) {
1275         minimizeMolecule(molecule);
1276     }
1277     return cleanPose;
1278 }
1279 
avoidClashes()1280 bool CoordgenMinimizer::avoidClashes()
1281 {
1282     bool allCleanPoses = true;
1283     if (skipAvoidClashes) {
1284         return true;
1285     }
1286     for (sketcherMinimizerMolecule* molecule : m_molecules) {
1287         auto cleanPose = avoidClashesOfMolecule(molecule);
1288         allCleanPoses = allCleanPoses && cleanPose;
1289     }
1290     return allCleanPoses;
1291 }
1292 
avoidInternalClashes(sketcherMinimizerFragment * fragment)1293 void CoordgenMinimizer::avoidInternalClashes(
1294     sketcherMinimizerFragment* fragment)
1295 {
1296     // avoid intraFragmentClashes
1297     vector<sketcherMinimizerAtom*> fragmentAtoms = fragment->getAtoms();
1298     for (sketcherMinimizerAtom* a : fragmentAtoms) {
1299         if (a->neighbors.size() != 1) {
1300             continue;
1301         }
1302         if (a->needsCheckForClashes) {
1303             continue;
1304         }
1305         if (a->fixed) {
1306             continue;
1307         }
1308         if (!fragment->getDofsOfAtom(a).empty()) {
1309             continue;
1310         }
1311         for (sketcherMinimizerAtom* a2 : fragmentAtoms) {
1312             if (a == a2) {
1313                 continue;
1314             }
1315             if (!fragment->getDofsOfAtom(a2).empty()) {
1316                 continue;
1317             }
1318             if (sketcherMinimizer::getBond(a, a2)) {
1319                 continue;
1320             }
1321             float dx = a2->coordinates.x() - a->coordinates.x();
1322             if (dx > bondLength * 0.5f) {
1323                 continue;
1324             }
1325             if (dx < -bondLength * 0.5f) {
1326                 continue;
1327             }
1328             float dy = a2->coordinates.y() - a->coordinates.y();
1329             if (dy > bondLength * 0.5f) {
1330                 continue;
1331             }
1332             if (dy < -bondLength * 0.5f) {
1333                 continue;
1334             }
1335             float squareD = dx * dx + dy * dy;
1336             if (squareD > bondLength * 0.5f * bondLength * 0.5f) {
1337                 continue;
1338             }
1339 
1340             sketcherMinimizerPointF vec =
1341                 a->coordinates - a->neighbors[0]->coordinates;
1342             vec *= 0.3f;
1343             a->coordinates -= vec;
1344             if (a2->neighbors.size() == 1) {
1345                 a2->coordinates += vec;
1346                 a2->coordinates.round();
1347             }
1348         }
1349     }
1350 }
1351 
bondsClash(sketcherMinimizerBond * bond,sketcherMinimizerBond * bond2) const1352 bool CoordgenMinimizer::bondsClash(sketcherMinimizerBond* bond,
1353                                    sketcherMinimizerBond* bond2) const
1354 {
1355     if (bond == bond2) {
1356         return false;
1357     }
1358     if (bond->getStartAtom() == bond2->getStartAtom() ||
1359         bond->getStartAtom() == bond2->getEndAtom() ||
1360         bond->getEndAtom() == bond2->getStartAtom() ||
1361         bond->getEndAtom() == bond2->getEndAtom()) {
1362         return false;
1363     }
1364     auto& start1 = bond->getStartAtom()->coordinates;
1365     auto& start2 = bond2->getStartAtom()->coordinates;
1366     auto& end1 = bond->getEndAtom()->coordinates;
1367     auto& end2 = bond2->getEndAtom()->coordinates;
1368     // coincidence and intersection calculations are expensive. Often bonds
1369     // are nowhere near each other, so skip the remaining work if a bond is
1370     // strictly to the left or right of another bond.
1371     if (max(start1.x(), end1.x()) < min(start2.x(), end2.x()) ||
1372         max(start1.y(), end1.y()) < min(start2.y(), end2.y()) ||
1373         min(start1.x(), end1.x()) > max(start2.x(), end2.x()) ||
1374         min(start1.y(), end1.y()) > max(start2.y(), end2.y())) {
1375         return false;
1376     }
1377 
1378     if (sketcherMinimizerMaths::pointsCoincide(
1379             bond->getStartAtom()->coordinates,
1380             bond2->getStartAtom()->coordinates) ||
1381         sketcherMinimizerMaths::pointsCoincide(
1382             bond->getStartAtom()->coordinates,
1383             bond2->getEndAtom()->coordinates) ||
1384         sketcherMinimizerMaths::pointsCoincide(
1385             bond->getEndAtom()->coordinates,
1386             bond2->getStartAtom()->coordinates) ||
1387         sketcherMinimizerMaths::pointsCoincide(
1388             bond->getEndAtom()->coordinates,
1389             bond2->getEndAtom()->coordinates)) {
1390         return true;
1391     }
1392     return (sketcherMinimizerMaths::intersectionOfSegments(
1393         bond->startAtom->coordinates, bond->endAtom->coordinates,
1394         bond2->startAtom->coordinates, bond2->endAtom->coordinates));
1395 }
1396 
avoidTerminalClashes(sketcherMinimizerMolecule * molecule,float & clashE)1397 void CoordgenMinimizer::avoidTerminalClashes(
1398     sketcherMinimizerMolecule* molecule, float& clashE)
1399 {
1400     if (clashE < 0.1) {
1401         return;
1402     }
1403     for (auto bond : molecule->getBonds()) {
1404         if (bond->isResidueInteraction()) {
1405             continue;
1406         }
1407         if (!bond->isTerminal()) {
1408             continue;
1409         }
1410         sketcherMinimizerAtom* terminalAtom = bond->getEndAtom();
1411         sketcherMinimizerAtom* rootAtom = bond->getStartAtom();
1412         if (terminalAtom->getBonds().size() != 1) {
1413             terminalAtom = bond->getStartAtom();
1414             rootAtom = bond->getEndAtom();
1415         }
1416         if (terminalAtom->fixed) {
1417             continue;
1418         }
1419         for (auto bond2 : molecule->getBonds()) {
1420             if (bond2->isResidueInteraction()) {
1421                 continue;
1422             }
1423             if (bondsClash(bond, bond2)) {
1424                 terminalAtom->setCoordinates(rootAtom->getCoordinates() +
1425                                              (terminalAtom->getCoordinates() -
1426                                               rootAtom->getCoordinates()) *
1427                                                  0.1);
1428             }
1429         }
1430     }
1431     clashE = scoreClashes(molecule);
1432 }
1433 
maybeMinimizeRings(const vector<sketcherMinimizerRing * > & rings)1434 void CoordgenMinimizer::maybeMinimizeRings(
1435     const vector<sketcherMinimizerRing*>& rings)
1436 {
1437     bool found = false;
1438     for (auto r : rings) {
1439         if (r->_atoms.size() == 5) {
1440             for (auto& _atom : r->_atoms) {
1441                 if (_atom->rings.size() > 2) {
1442                     found = true;
1443                 }
1444             }
1445         }
1446         if (r->isMacrocycle() && r->_atoms.size() % 2 != 0) {
1447             for (auto& _atom : r->_atoms) {
1448                 if (_atom->rings.size() > 2) {
1449                     found = true;
1450                 }
1451             }
1452         }
1453     }
1454 
1455     if (!found) {
1456         return;
1457     }
1458     rings.at(0)->getAtoms().at(0)->molecule->requireMinimization();
1459 }
1460 
buildMoleculeFromFragments(sketcherMinimizerMolecule * molecule,bool firstTime) const1461 void CoordgenMinimizer::buildMoleculeFromFragments(
1462     sketcherMinimizerMolecule* molecule, bool firstTime) const
1463 {
1464     for (auto fragment : molecule->getFragments()) {
1465         float angle = 0;
1466         sketcherMinimizerPointF position(0.f, 0.f);
1467         if (fragment->getParent()) {
1468             sketcherMinimizerPointF p1 =
1469                 fragment->_bondToParent->startAtom->coordinates;
1470             sketcherMinimizerPointF p2 =
1471                 fragment->_bondToParent->endAtom->coordinates;
1472             sketcherMinimizerPointF p = p2 - p1;
1473             angle = atan2(-p.y(), p.x());
1474             position = fragment->_bondToParent->endAtom->coordinates;
1475             if (firstTime) {
1476                 sketcherMinimizer::alignWithParentDirection(fragment, position,
1477                                                             angle);
1478             }
1479         }
1480         fragment->setCoordinates(position, angle);
1481     }
1482 }
1483 
buildFromFragments(bool firstTime) const1484 void CoordgenMinimizer::buildFromFragments(bool firstTime) const
1485 {
1486     for (sketcherMinimizerMolecule* molecule : m_molecules) {
1487         buildMoleculeFromFragments(molecule, firstTime);
1488     }
1489 }
1490 
hasValid3DCoordinates(const vector<sketcherMinimizerAtom * > & atoms)1491 bool CoordgenMinimizer::hasValid3DCoordinates(
1492     const vector<sketcherMinimizerAtom*>& atoms)
1493 {
1494     for (sketcherMinimizerAtom* atom : atoms) {
1495         if (!atom->hasValid3DCoordinates()) {
1496             return false;
1497         }
1498     }
1499     return true;
1500 }
1501 
fallbackOn3DCoordinates(const vector<sketcherMinimizerAtom * > & atoms)1502 void CoordgenMinimizer::fallbackOn3DCoordinates(
1503     const vector<sketcherMinimizerAtom*>& atoms)
1504 {
1505     float scale = 35.f; // ratio between average bond length and 2d bond length
1506     /* TODO find best projection */
1507     for (sketcherMinimizerAtom* atom : atoms) {
1508         atom->setCoordinates(
1509             sketcherMinimizerPointF(atom->m_x3D * scale, -atom->m_y3D * scale));
1510     }
1511 }
1512 
hasNaNCoordinates(const std::vector<sketcherMinimizerAtom * > & atoms)1513 bool CoordgenMinimizer::hasNaNCoordinates(
1514     const std::vector<sketcherMinimizerAtom*>& atoms)
1515 {
1516     for (sketcherMinimizerAtom* a : atoms)
1517         if (a->coordinates.x() != a->coordinates.x() ||
1518             a->coordinates.y() != a->coordinates.y()) {
1519             return true;
1520         }
1521     return false;
1522 }
1523 
hasNaNCoordinates()1524 bool CoordgenMinimizer::hasNaNCoordinates()
1525 {
1526     return hasNaNCoordinates(m_atoms);
1527 }
1528 
checkForClashes(sketcherMinimizerAtom * a)1529 void CoordgenMinimizer::checkForClashes(sketcherMinimizerAtom* a)
1530 {
1531 
1532     if (a->fixed) {
1533         return;
1534     }
1535 
1536     sketcherMinimizerPointF oldCoordinates = a->coordinates;
1537     vector<sketcherMinimizerPointF> coordsVect;
1538     coordsVect.push_back(oldCoordinates);
1539     coordsVect.push_back(oldCoordinates +
1540                          sketcherMinimizerPointF(bondLength * 0.25f, 0.f));
1541     coordsVect.push_back(oldCoordinates +
1542                          sketcherMinimizerPointF(-bondLength * 0.25f, 0.f));
1543     coordsVect.push_back(oldCoordinates +
1544                          sketcherMinimizerPointF(0.f, bondLength * 0.25f));
1545     coordsVect.push_back(oldCoordinates +
1546                          sketcherMinimizerPointF(0.f, -bondLength * 0.25f));
1547     coordsVect.push_back(oldCoordinates + sketcherMinimizerPointF(
1548                                               bondLength * 0.25f * 0.7071f,
1549                                               -bondLength * 0.25f * 0.7071f));
1550     coordsVect.push_back(oldCoordinates + sketcherMinimizerPointF(
1551                                               -bondLength * 0.25f * 0.7071f,
1552                                               -bondLength * 0.25f * 0.7071f));
1553     coordsVect.push_back(oldCoordinates +
1554                          sketcherMinimizerPointF(-bondLength * 0.25f * 0.7071f,
1555                                                  bondLength * 0.25f * 0.7071f));
1556     coordsVect.push_back(oldCoordinates +
1557                          sketcherMinimizerPointF(bondLength * 0.25f * 0.7071f,
1558                                                  bondLength * 0.25f * 0.7071f));
1559     float bestE = 999999.f;
1560     int bestI = 0;
1561     for (unsigned int i = 0; i < coordsVect.size(); i++) {
1562         a->coordinates = coordsVect[i];
1563 
1564         // solves intrafragment clashes by shifting the atomic coordinates up
1565         // down left right or diagonally
1566         sketcherMinimizerClashInteraction clashI(a, a, a);
1567         clashI.restV = 300;
1568         float clashE = 0;
1569         vector<sketcherMinimizerBond*> bonds = a->getFragment()->getBonds();
1570         for (sketcherMinimizerBond* b : bonds) {
1571             if (!b->startAtom->coordinatesSet) {
1572                 continue;
1573             }
1574             if (!b->endAtom->coordinatesSet) {
1575                 continue;
1576             }
1577             if (b->startAtom == a) {
1578                 continue;
1579             }
1580             if (b->endAtom == a) {
1581                 continue;
1582             }
1583             clashI.atom1 = b->startAtom;
1584             clashI.atom2 = a;
1585             clashI.atom3 = b->endAtom;
1586             clashI.energy(clashE);
1587         }
1588 
1589         for (sketcherMinimizerBond* b : a->bonds) {
1590             vector<sketcherMinimizerAtom*> atoms = a->getFragment()->getAtoms();
1591             for (sketcherMinimizerAtom* atom : atoms) {
1592                 if (atom == a) {
1593                     continue;
1594                 }
1595                 if (!b->startAtom->coordinatesSet) {
1596                     continue;
1597                 }
1598                 if (!b->endAtom->coordinatesSet) {
1599                     continue;
1600                 }
1601                 if (b->startAtom == atom) {
1602                     continue;
1603                 }
1604                 if (b->endAtom == atom) {
1605                     continue;
1606                 }
1607                 clashI.atom1 = b->startAtom;
1608                 clashI.atom2 = atom;
1609                 clashI.atom3 = b->endAtom;
1610                 clashI.energy(clashE);
1611             }
1612         }
1613 
1614         if (clashE < SKETCHER_EPSILON) {
1615             return;
1616         }
1617         if (i == 0) {
1618             bestE = clashE;
1619         }
1620         if (clashE < bestE) {
1621             bestE = clashE;
1622             bestI = i;
1623         }
1624     }
1625     a->setCoordinates(coordsVect[bestI]);
1626 }
1627 
getPrecision() const1628 float CoordgenMinimizer::getPrecision() const
1629 {
1630     return m_precision;
1631 }
1632 
setPrecision(float f)1633 void CoordgenMinimizer::setPrecision(float f)
1634 {
1635     m_precision = f;
1636 }
1637 
1638 std::pair<std::vector<short unsigned int>, float>
findBestSolution() const1639 CoordgenDOFSolutions::findBestSolution() const
1640 {
1641     std::pair<std::vector<short unsigned int>, float> bestSolution =
1642         *m_solutions.begin();
1643     for (auto solution : m_solutions) {
1644         if (solution.second < bestSolution.second) {
1645             bestSolution = solution;
1646         }
1647     }
1648     return bestSolution;
1649 }
1650 
getCurrentSolution()1651 std::vector<short unsigned int> CoordgenDOFSolutions::getCurrentSolution()
1652 {
1653     std::vector<short unsigned int> solution;
1654     for (auto dof : m_allDofs) {
1655         solution.push_back(dof->getCurrentState());
1656     }
1657     return solution;
1658 }
1659 
loadSolution(const std::vector<short unsigned int> & solution)1660 void CoordgenDOFSolutions::loadSolution(
1661     const std::vector<short unsigned int>& solution)
1662 {
1663     assert(solution.size() == m_allDofs.size());
1664     for (unsigned int i = 0; i < solution.size(); ++i) {
1665         m_allDofs.at(i)->setState(solution.at(i));
1666     }
1667 }
1668 
hasSolution(const std::vector<short unsigned int> & solution)1669 bool CoordgenDOFSolutions::hasSolution(
1670     const std::vector<short unsigned int>& solution)
1671 {
1672     return m_solutions.find(solution) != m_solutions.end();
1673 }
1674 
scoreCurrentSolution()1675 float CoordgenDOFSolutions::scoreCurrentSolution()
1676 {
1677     std::vector<short unsigned int> solution;
1678     for (auto dof : m_allDofs) {
1679         solution.push_back(dof->getCurrentState());
1680     }
1681     //   for (auto dof : solution) cerr <<dof;
1682     //   cerr << endl;
1683     auto position = m_solutions.find(solution);
1684     if (position != m_solutions.end()) {
1685         return position->second;
1686     } else {
1687         if (m_solutions.size() >
1688             MAXIMUM_NUMBER_OF_SCORED_SOLUTIONS * m_minimizer->getPrecision()) {
1689             return REJECTED_SOLUTION_SCORE;
1690         }
1691         m_minimizer->buildMoleculeFromFragments(m_molecule, false);
1692         float result = m_minimizer->scoreClashes(m_molecule, true);
1693         m_solutions[solution] = result;
1694         return result;
1695     }
1696 }
1697