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