/*
Contributors: Nicola Zonta
Copyright Schrodinger, LLC. All rights reserved
*/
#include "CoordgenFragmentBuilder.h"
#include "CoordgenMinimizer.h"
#include "sketcherMinimizerFragment.h"
#include "sketcherMinimizer.h"
#include <numeric>
#include <algorithm>
using namespace std;
const int bondLength = BONDLENGTH;
const int PERFECTLY_PLANAR_SYSTEM_SCORE = 50;
const int NON_PLANAR_SYSTEM_SCORE = 1000;
const int UNTREATABLE_SYSTEM_PLANARITY_SCORE = 200000;
void CoordgenFragmentBuilder::initializeCoordinates(
sketcherMinimizerFragment* fragment) const
{
assert(fragment->getAtoms().size());
buildFragment(fragment);
fragment->storeCoordinateInformation();
}
void CoordgenFragmentBuilder::rotateMainFragment(
sketcherMinimizerFragment* f) const
{
if (f->fixed)
return;
if (f->isTemplated)
return;
if (!f->constrained)
return;
sketcherMinimizerPointF constrainOldCenter(0.f, 0.f);
sketcherMinimizerPointF constrainNewCenter(0.f, 0.f);
vector<sketcherMinimizerAtom*> constrainedAtoms;
for (const auto& atom : f->getAtoms()) {
if (atom->constrained) {
constrainedAtoms.push_back(atom);
}
}
for (const auto& child : f->_children) {
sketcherMinimizerAtom* atom = child->_bondToParent->endAtom;
if (atom->constrained) {
constrainedAtoms.push_back(atom);
}
}
for (const auto& a : constrainedAtoms) {
constrainOldCenter += a->templateCoordinates;
constrainNewCenter += a->coordinates;
}
if (constrainedAtoms.size() > 0) {
constrainOldCenter /= constrainedAtoms.size();
constrainNewCenter /= constrainedAtoms.size();
}
vector<sketcherMinimizerPointF> v1, v2;
for (const auto& a : constrainedAtoms) {
v2.push_back(a->coordinates - constrainNewCenter);
v1.push_back(a->templateCoordinates - constrainOldCenter);
}
float rotMat[4];
sketcherMinimizer::alignmentMatrix(v1, v2, rotMat);
vector<sketcherMinimizerPointF> rotatedV2;
for (auto p : v2) {
auto rotatedPoint =
sketcherMinimizerPointF(p.x() * rotMat[0] + p.y() * rotMat[1],
p.x() * rotMat[2] + p.y() * rotMat[3]);
rotatedV2.push_back(rotatedPoint);
}
for (sketcherMinimizerAtom* a : f->getAtoms()) {
sketcherMinimizerPointF v = a->getCoordinates() - constrainNewCenter;
v = sketcherMinimizerPointF(v.x() * rotMat[0] + v.y() * rotMat[1],
v.x() * rotMat[2] + v.y() * rotMat[3]);
a->setCoordinates(v + constrainOldCenter);
}
for (sketcherMinimizerFragment* child : f->_children) {
sketcherMinimizerAtom* at = child->_bondToParent->endAtom;
sketcherMinimizerPointF v = at->getCoordinates() - constrainNewCenter;
v = sketcherMinimizerPointF(v.x() * rotMat[0] + v.y() * rotMat[1],
v.x() * rotMat[2] + v.y() * rotMat[3]);
at->setCoordinates(v + constrainOldCenter);
at->coordinatesSet = false;
}
}
bool CoordgenFragmentBuilder::findTemplate(
const vector<sketcherMinimizerRing*> rings) const
{
vector<sketcherMinimizerAtom*> allAtoms =
rings[0]->_atoms[0]->fragment->getAtoms();
bool foundTemplate = false;
sketcherMinimizer::loadTemplates();
vector<sketcherMinimizerAtom*> atoms;
vector<sketcherMinimizerBond*> bonds;
map<sketcherMinimizerAtom*, bool> isVisited;
vector<int> oldIndices;
for (sketcherMinimizerRing* r : rings) {
for (sketcherMinimizerAtom* a : r->_atoms) {
if (!isVisited[a]) {
isVisited[a] = true;
oldIndices.push_back(a->_generalUseN);
a->_generalUseN = atoms.size();
atoms.push_back(a);
}
}
}
for (sketcherMinimizerRing* r : rings) {
for (sketcherMinimizerBond* b : r->_bonds) {
bool found = false;
for (sketcherMinimizerBond* b2 : bonds) {
if (b2 == b) {
found = true;
break;
}
}
if (!found) {
bonds.push_back(b);
}
}
}
vector<unsigned int> mapping;
for (auto& temp : sketcherMinimizer::m_templates.getTemplates()) {
foundTemplate = sketcherMinimizer::compare(atoms, bonds, temp, mapping);
if (foundTemplate) {
if (atoms.size() > 0)
atoms[0]->fragment->isTemplated = true;
for (unsigned int i = 0; i < atoms.size(); i++) {
atoms[i]->setCoordinates(temp->_atoms[mapping[i]]->coordinates *
bondLength);
atoms[i]->rigid = true;
}
for (sketcherMinimizerRing* r : rings) {
r->coordinatesGenerated = true;
}
break;
}
}
for (unsigned int i = 0; i < oldIndices.size(); i++) {
atoms[i]->_generalUseN = oldIndices[i];
}
return foundTemplate;
}
sketcherMinimizerRing* CoordgenFragmentBuilder::findCentralRingOfSystem(
const vector<sketcherMinimizerRing*> rings) const
{
sketcherMinimizerRing* highest = rings.at(0);
int high_score = -1;
for (sketcherMinimizerRing* r : rings) {
int priority = 0;
if (r->isMacrocycle())
priority += 1000;
if (r->_atoms.size() == 6)
priority += 100;
priority += r->_atoms.size();
priority += 10 * (r->fusedWith.size());
if (priority > high_score) {
highest = r;
high_score = priority;
}
}
return highest;
}
void CoordgenFragmentBuilder::generateCoordinatesCentralRings(
vector<sketcherMinimizerRing*> rings) const
{
if (rings.size() == 1) { // linearfusedRings
buildRing(rings[0]);
} else {
bool foundTemplate = findTemplate(rings);
if (!foundTemplate) {
float planarityScore = newScorePlanarity(rings);
if (planarityScore < NON_PLANAR_SYSTEM_SCORE) {
bool needsTemplate =
planarityScore > PERFECTLY_PLANAR_SYSTEM_SCORE;
if (needsTemplate) {
findTemplate(rings);
}
while (rings.size() != 0) {
auto centralRing = findCentralRingOfSystem(rings);
buildRing(centralRing);
rings.erase(
std::remove(rings.begin(), rings.end(), centralRing),
rings.end());
}
CoordgenMinimizer::maybeMinimizeRings(rings);
} else if (planarityScore > UNTREATABLE_SYSTEM_PLANARITY_SCORE) {
return;
} else {
sketcherMinimizerRing* firstRing =
findCentralRingOfSystem(rings);
m_macrocycleBuilder.openCycleAndGenerateCoords(firstRing);
firstRing->getAtoms()
.at(0)
->getMolecule()
->requireMinimization();
}
}
}
}
float CoordgenFragmentBuilder::newScorePlanarity(
const vector<sketcherMinimizerRing*> rings)
const // if score > 1000 then it is not planar
{
float score = 0.f;
for (const auto& ring : rings) {
if (ring->isMacrocycle() &&
m_macrocycleBuilder.findBondToOpen(ring) == NULL) {
continue;
}
if (ring->isMacrocycle()) {
for (const auto& otherRing : ring->fusedWith) {
if (otherRing->isMacrocycle()) {
score += NON_PLANAR_SYSTEM_SCORE;
}
}
}
for (const auto& bond : ring->_bonds) {
if (bond->rings.size() > 2) {
score += NON_PLANAR_SYSTEM_SCORE * (bond->rings.size() - 2);
}
}
for (const auto& atom : ring->getAtoms()) {
if (atom->neighbors.size() > 3) {
float angle = 0;
for (sketcherMinimizerRing* r : atom->rings) {
angle += M_PI - (2 * M_PI / r->_atoms.size());
}
if (angle >= 1.99 * M_PI) {
score += NON_PLANAR_SYSTEM_SCORE;
}
}
}
}
return score;
}
void CoordgenFragmentBuilder::generateCoordinatesSideRings(
stack<sketcherMinimizerRing*> sideRings) const
{
while (sideRings.size()) {
sketcherMinimizerRing* ring = sideRings.top();
sideRings.pop();
buildRing(ring);
}
}
sketcherMinimizerRing*
CoordgenFragmentBuilder::getSharedAtomsWithAlreadyDrawnRing(
const sketcherMinimizerRing* ring,
vector<sketcherMinimizerAtom*>& fusionAtoms,
sketcherMinimizerBond*& fusionBond) const
{
sketcherMinimizerRing* parent = NULL;
for (unsigned int i = 0; i < ring->fusedWith.size(); i++) {
if (ring->fusedWith[i]->coordinatesGenerated) {
if (!parent)
parent = ring->fusedWith[i];
else {
if (ring->fusedWith[i]->_atoms.size() > parent->_atoms.size())
parent = ring->fusedWith[i];
}
}
}
if (parent) {
for (unsigned int i = 0; i < parent->fusedWith.size(); i++) {
if (parent->fusedWith[i] == ring) {
fusionAtoms = parent->fusionAtoms[i];
}
}
for (sketcherMinimizerBond* b : parent->fusionBonds) {
if (ring->containsAtom(b->startAtom) ||
ring->containsAtom(b->endAtom)) {
fusionBond = b;
break;
}
}
}
return parent;
}
vector<sketcherMinimizerAtom*> CoordgenFragmentBuilder::orderChainOfAtoms(
const vector<sketcherMinimizerAtom*> atoms,
sketcherMinimizerAtom* startAtom)
{
vector<sketcherMinimizerAtom*> orderedAtoms;
map<sketcherMinimizerAtom*, bool> stillToAdd;
for (sketcherMinimizerAtom* a : atoms) {
stillToAdd[a] = true;
}
sketcherMinimizerAtom* atomToAdd = startAtom;
while (orderedAtoms.size() < atoms.size()) {
orderedAtoms.push_back(atomToAdd);
stillToAdd[atomToAdd] = false;
if (orderedAtoms.size() >= atoms.size())
break;
for (sketcherMinimizerAtom* neighbor : atomToAdd->neighbors) {
if (stillToAdd[neighbor]) {
atomToAdd = neighbor;
break;
}
}
}
return orderedAtoms;
}
vector<sketcherMinimizerAtom*>
CoordgenFragmentBuilder::orderRingAtoms(const sketcherMinimizerRing* ring)
{
vector<sketcherMinimizerAtom*> ringAtoms = ring->getAtoms();
assert(!ringAtoms.empty());
return orderChainOfAtoms(ringAtoms, ringAtoms.at(0));
}
void CoordgenFragmentBuilder::buildRing(sketcherMinimizerRing* ring) const
{
if (ring->coordinatesGenerated)
return;
vector<sketcherMinimizerAtom*> atoms;
sketcherMinimizerRing* parent = NULL;
sketcherMinimizerBond* fusionBond = NULL;
vector<sketcherMinimizerAtom*> fusionAtoms;
parent = getSharedAtomsWithAlreadyDrawnRing(ring, fusionAtoms, fusionBond);
atoms = orderRingAtoms(ring);
vector<sketcherMinimizerPointF> coords;
if (atoms.size() >= MACROCYCLE) {
coords = m_macrocycleBuilder.newMacrocycle(ring, atoms);
} else {
coords = listOfCoordinatesFromListofRingAtoms(atoms);
}
if (parent) {
if (fusionAtoms.size() < 2) {
sketcherMinimizerPointF coordinatesOfPivotAtom;
sketcherMinimizerAtom* pivotAtom = NULL;
sketcherMinimizerAtom* pivotAtomOnParent = NULL;
if (!fusionAtoms.empty()) {
pivotAtom = fusionAtoms.at(0);
pivotAtomOnParent = fusionAtoms.at(0);
} else {
assert(fusionBond != NULL);
if (find(atoms.begin(), atoms.end(),
fusionBond->getStartAtom()) != atoms.end()) {
pivotAtom = fusionBond->getStartAtom();
pivotAtomOnParent = fusionBond->getEndAtom();
} else {
pivotAtom = fusionBond->getEndAtom();
pivotAtomOnParent = fusionBond->getStartAtom();
}
}
assert(find(atoms.begin(), atoms.end(), pivotAtom) != atoms.end());
assert(parent->containsAtom(pivotAtomOnParent));
for (unsigned int i = 0; i < atoms.size(); i++) {
if (atoms[i] == pivotAtom) {
coordinatesOfPivotAtom = coords[i];
break;
}
}
for (unsigned int i = 0; i < coords.size(); i++) {
coords[i] -= coordinatesOfPivotAtom;
}
sketcherMinimizerPointF center =
accumulate(coords.begin(), coords.end(),
sketcherMinimizerPointF(0, 0)) /
coords.size();
center.normalize();
int neighborsN = 0;
sketcherMinimizerPointF neighborsMean(0.f, 0.f);
for (sketcherMinimizerAtom* neighbor :
pivotAtomOnParent->neighbors) {
if (parent->containsAtom(neighbor)) {
neighborsMean += neighbor->getCoordinates();
++neighborsN;
}
}
assert(neighborsN == 2);
neighborsMean *= 0.5;
sketcherMinimizerPointF pivotAtomCoordinates =
pivotAtomOnParent->getCoordinates();
sketcherMinimizerPointF direction =
pivotAtomCoordinates - neighborsMean;
direction.normalize();
sketcherMinimizerPointF bondVector(0.f, 0.f);
if (fusionBond != NULL) {
bondVector = direction * bondLength;
}
for (unsigned int i = 0; i < coords.size(); i++) {
sketcherMinimizerPointF coordinates = coords[i];
coordinates.rotate(center.y(), center.x());
coordinates.rotate(-direction.y(), direction.x());
coords[i] = coordinates + pivotAtomCoordinates + bondVector;
}
for (unsigned int i = 0; i < coords.size(); i++) {
atoms[i]->setCoordinates(coords[i]);
}
if (fusionBond) {
// check ZE inversions
if (!fusionBond->checkStereoChemistry()) {
for (sketcherMinimizerAtom* atom : atoms) {
sketcherMinimizerAtom::mirrorCoordinates(atom,
fusionBond);
}
}
}
}
else {
// build 2 lists of coordinates, mirror of each other.
vector<sketcherMinimizerPointF> coords2 = coords;
for (unsigned int i = 0; i < coords2.size(); i++) {
coords2[i] =
sketcherMinimizerPointF(coords2[i].x(), -coords2[i].y());
}
map<sketcherMinimizerAtom *, sketcherMinimizerPointF> map1, map2;
for (unsigned int i = 0; i < atoms.size(); i++) {
map1[atoms[i]] = coords[i];
map2[atoms[i]] = coords2[i];
}
int lastI = fusionAtoms.size() - 1;
sketcherMinimizerPointF mean = (fusionAtoms[0]->coordinates +
fusionAtoms[lastI]->coordinates) *
0.5;
sketcherMinimizerPointF newMean =
(map1[fusionAtoms[0]] + map1[fusionAtoms[lastI]]) * 0.5;
sketcherMinimizerPointF v =
fusionAtoms[lastI]->coordinates - fusionAtoms[0]->coordinates;
float distance = v.length();
if (distance < SKETCHER_EPSILON) {
distance = SKETCHER_EPSILON;
}
float sine = v.y() / distance;
float cosine = v.x() / distance;
sketcherMinimizerPointF newV =
(map1[fusionAtoms[lastI]] - map1[fusionAtoms[0]]);
float newDistance = newV.length();
float newSine = newV.y() / newDistance;
float newCosine = newV.x() / newDistance;
for (unsigned int i = 0; i < atoms.size();
i++) { // rotate and translate coords
sketcherMinimizerPointF d = coords[i] - newMean;
d.rotate(newSine, newCosine);
d.rotate(-sine, cosine);
coords[i] = mean + d;
}
newV = (map2[fusionAtoms[lastI]] - map2[fusionAtoms[0]]);
newMean = (map2[fusionAtoms[0]] + map2[fusionAtoms[lastI]]) * 0.5;
newDistance = newV.length();
newSine = newV.y() / newDistance;
newCosine = newV.x() / newDistance;
for (unsigned int i = 0; i < atoms.size();
i++) { // rotate and translate coords2
sketcherMinimizerPointF d = coords2[i] - newMean;
d.rotate(newSine, newCosine);
d.rotate(-sine, cosine);
coords2[i] = mean + d;
}
// pick the list of coordinates that will place the ring outside the
// parent i.e further away from its center.
sketcherMinimizerPointF center = parent->findCenter();
map<sketcherMinimizerAtom*, bool> fusionMap;
for (unsigned int i = 0; i < fusionAtoms.size(); i++)
fusionMap[fusionAtoms[i]] = true;
float distance1 = 0.f, distance2 = 0.f;
for (unsigned int i = 0; i < atoms.size(); i++) {
if (fusionMap[atoms[i]])
continue;
sketcherMinimizerPointF v1 = coords[i] - center;
sketcherMinimizerPointF v2 = coords2[i] - center;
distance1 += v1.length();
distance2 += v2.length();
}
vector<sketcherMinimizerPointF>* targetCoords = &coords;
if (distance1 < distance2)
targetCoords = &coords2;
for (unsigned int i = 0; i < atoms.size(); i++) {
atoms[i]->setCoordinates((*targetCoords)[i]);
}
// TODO forbid if double bond to ring
if (ring->fusedWith.size() == 1 &&
(*fusionAtoms.begin())->hasNoStereoActiveBonds() &&
(*fusionAtoms.rbegin())->hasNoStereoActiveBonds()) {
CoordgenFlipRingDOF* dof =
new CoordgenFlipRingDOF(ring, fusionAtoms);
ring->getAtoms().at(0)->fragment->addDof(dof);
}
}
} else {
for (unsigned int i = 0; i < coords.size(); i++) {
atoms[i]->setCoordinates(coords[i]);
}
}
ring->coordinatesGenerated = true;
}
vector<sketcherMinimizerPointF>
CoordgenFragmentBuilder::listOfCoordinatesFromListofRingAtoms(
const vector<sketcherMinimizerAtom*> atoms)
{
vector<sketcherMinimizerPointF> out;
assert(atoms.size());
float a = 2 * M_PI / atoms.size();
sketcherMinimizerPointF coords(0.f, 0.f);
float angle = 0;
for (unsigned int n = 0; n < atoms.size(); n++) {
out.push_back(coords);
angle = a * n;
coords += sketcherMinimizerPointF(cos(angle) * bondLength,
-sin(angle) * bondLength);
}
return out;
}
void CoordgenFragmentBuilder::buildNonRingAtoms(
sketcherMinimizerFragment* fragment) const
{
set<sketcherMinimizerAtom*> isAtomVisited;
queue<sketcherMinimizerAtom*> q;
for (sketcherMinimizerAtom* atom : fragment->getAtoms()) {
if (!atom->getRings().empty()) {
q.push(atom);
isAtomVisited.insert(atom);
}
}
if (q.empty()) {
if (fragment->getParent()) {
fragment->_bondToParent->getStartAtom()->setCoordinates(
sketcherMinimizerPointF(-bondLength, 0));
}
sketcherMinimizerAtom* atomToStartFrom =
(fragment->getParent() ? fragment->_bondToParent->getEndAtom()
: fragment->getAtoms()[0]);
atomToStartFrom->setCoordinates(sketcherMinimizerPointF(0, 0));
q.push(atomToStartFrom);
isAtomVisited.insert(atomToStartFrom);
}
while (q.size()) {
generateCoordinatesNeighborsOfFirstAtomInQueue(q, isAtomVisited,
fragment);
}
}
vector<float> CoordgenFragmentBuilder::neighborsAnglesAtCenter(
const sketcherMinimizerAtom* atom) const
{
int angleDivision = atom->neighbors.size();
vector<float> angles;
if (!m_evenAngles) {
if (atom->neighbors.size() == 2) {
if (atom->getAtomicNumber() == 6 ||
!atom->neighbors.at(0)->crossLayout ||
!atom->neighbors.at(1)->crossLayout) {
angleDivision = 3;
}
if (atom->bonds[0]->getBondOrder() +
atom->bonds[1]->getBondOrder() >=
4)
angleDivision = 2;
} else if (atom->neighbors.size() == 4 && !atom->crossLayout) {
angles.push_back(M_PI / 3);
angles.push_back(M_PI * 0.5);
angles.push_back(2 * M_PI / 3);
angles.push_back(M_PI * 0.5);
}
}
if (!angles.size()) {
for (unsigned int i = 0; i < atom->neighbors.size(); i++) {
angles.push_back(2 * M_PI / angleDivision);
}
}
return angles;
}
void CoordgenFragmentBuilder::initializeVariablesForNeighboursCoordinates(
sketcherMinimizerAtom* atom, set<sketcherMinimizerAtom*>& isAtomVisited,
sketcherMinimizerPointF& startCoordinates,
vector<sketcherMinimizerAtom*>& orderedNeighbors,
vector<float>& angles) const
{
if (atom->rings.size()) {
return initializeVariablesForNeighboursCoordinatesRingAtom(
atom, isAtomVisited, startCoordinates, orderedNeighbors, angles);
}
orderedNeighbors.clear();
if (atom->neighbors.size() != 4) {
orderedNeighbors = atom->neighbors;
} else {
vector<sketcherMinimizerAtomPriority> atomPriorities;
for (unsigned int i = 0; i < atom->neighbors.size(); i++) {
sketcherMinimizerAtomPriority p;
p.a = atom->neighbors[i];
atomPriorities.push_back(p);
}
sketcherMinimizerAtom::orderAtomPriorities(atomPriorities, atom);
for (sketcherMinimizerAtomPriority ap : atomPriorities) {
orderedNeighbors.push_back(ap.a);
}
}
unsigned int startN = 0;
for (unsigned int i = 0; i < orderedNeighbors.size(); i++) {
if (isAtomVisited.find(orderedNeighbors[i]) != isAtomVisited.end()) {
startCoordinates =
orderedNeighbors[i]->getCoordinates() - atom->getCoordinates();
startN = i;
break;
}
}
for (unsigned int i = 0; i < startN; i++) {
orderedNeighbors.push_back(orderedNeighbors.at(0));
orderedNeighbors.erase(orderedNeighbors.begin());
}
angles = neighborsAnglesAtCenter(atom);
}
void CoordgenFragmentBuilder::
initializeVariablesForNeighboursCoordinatesRingAtom(
const sketcherMinimizerAtom* atom, set<sketcherMinimizerAtom*>&,
sketcherMinimizerPointF& startCoordinates,
vector<sketcherMinimizerAtom*>& orderedNeighbors,
vector<float>& angles) const
{
vector<pair<float, sketcherMinimizerAtom*>> ringNeighboursAndAngles;
orderedNeighbors.clear();
for (sketcherMinimizerAtom* neigh : atom->neighbors) {
if (sketcherMinimizer::sameRing(neigh, atom)) {
float ang = atan2(neigh->coordinates.y() - atom->coordinates.y(),
neigh->coordinates.x() - atom->coordinates.x());
if (ang < 0)
ang += 2 * M_PI;
pair<float, sketcherMinimizerAtom*> pairToAdd(ang, neigh);
ringNeighboursAndAngles.push_back(pairToAdd);
} else {
orderedNeighbors.push_back(neigh);
}
}
stable_sort(ringNeighboursAndAngles.begin(), ringNeighboursAndAngles.end());
vector<float> gaps;
vector<float> scaledGaps;
for (unsigned int i = 0; i < ringNeighboursAndAngles.size(); i++) {
int next = (i + 1) % ringNeighboursAndAngles.size();
float gap = ringNeighboursAndAngles[next].first -
ringNeighboursAndAngles[i].first;
if (gap < 0)
gap += 2 * M_PI;
gaps.push_back(gap);
bool rin = false;
float mid = ringNeighboursAndAngles[i].first + gap * 0.5;
float sine = sin(-mid);
float cosine = cos(-mid);
sketcherMinimizerPointF p(bondLength * 0.1, 0);
p.rotate(sine, cosine);
sketcherMinimizerPointF point = atom->coordinates + p;
vector<sketcherMinimizerRing*> rings = atom->getFragment()->getRings();
for (sketcherMinimizerRing* r : rings) {
if (r->isMacrocycle())
continue;
if (r->contains(point)) {
rin = true;
break;
}
}
float scaledGap = gap;
if (gap > M_PI)
scaledGap *= 10;
else if (rin)
scaledGap *= 0.2;
scaledGaps.push_back(scaledGap);
}
int bestI = 0;
for (unsigned int i = 0; i < gaps.size(); i++) {
if (scaledGaps[i] > scaledGaps[bestI]) {
bestI = i;
}
}
angles.clear();
float biggestGap = gaps[bestI];
startCoordinates = ringNeighboursAndAngles[bestI].second->getCoordinates() -
atom->getCoordinates();
for (unsigned int i = 0; i < orderedNeighbors.size(); i++) {
angles.push_back(-biggestGap / (orderedNeighbors.size() + 1));
}
}
void CoordgenFragmentBuilder::generateCoordinatesNeighborsOfFirstAtomInQueue(
queue<sketcherMinimizerAtom*>& atomQueue,
set<sketcherMinimizerAtom*>& isAtomVisited,
const sketcherMinimizerFragment* fragment) const
{
sketcherMinimizerAtom* at = atomQueue.front();
atomQueue.pop();
sketcherMinimizerPointF initialCoords(-bondLength, 0.f);
vector<sketcherMinimizerAtom*> orderedNeighbors;
vector<float> angles;
initializeVariablesForNeighboursCoordinates(
at, isAtomVisited, initialCoords, orderedNeighbors, angles);
for (unsigned int i = 0; i < orderedNeighbors.size(); i++) {
int n = i;
sketcherMinimizerAtom* neigh = orderedNeighbors[n];
if (isAtomVisited.find(neigh) != isAtomVisited.end())
continue;
float s = sin(angles[n]);
float c = cos(angles[n]);
initialCoords.rotate(s, c);
neigh->setCoordinates(at->coordinates + initialCoords);
if (neigh->fragment == fragment) {
isAtomVisited.insert(neigh);
atomQueue.push(neigh);
} else {
neigh->coordinatesSet = false;
}
if (at->needsCheckForClashes) {
neigh->needsCheckForClashes = true;
CoordgenMinimizer::checkForClashes(neigh);
}
for (auto& dof : at->fragment->getDofsOfAtom(at)) {
if (dof->getFragment() == fragment) {
dof->addAtom(neigh);
}
}
}
avoidZEInversions(at, isAtomVisited);
maybeAddMacrocycleDOF(at);
for (auto& neighbor : at->neighbors) {
if (!sketcherMinimizerAtom::shareARing(at, neighbor) &&
at->getFragment() == neighbor->getFragment()) {
CoordgenScaleAtomsDOF* dof = new CoordgenScaleAtomsDOF(at);
dof->addAtom(neighbor);
at->fragment->addDof(dof);
}
}
}
void CoordgenFragmentBuilder::maybeAddMacrocycleDOF(
sketcherMinimizerAtom* atom) const
{
if (atom->getRings().size() == 1 &&
atom->getRings().at(0)->isMacrocycle() && atom->neighbors.size() == 3) {
for (auto& bond : atom->getBonds()) {
if (bond->isStereo() && !bond->isTerminal())
return;
}
for (auto& neighbor : atom->neighbors) {
if (!sketcherMinimizerAtom::shareARing(atom, neighbor)) {
CoordgenInvertBondDOF* dof =
new CoordgenInvertBondDOF(atom, neighbor);
atom->fragment->addDof(dof);
}
}
}
}
void CoordgenFragmentBuilder::avoidZEInversions(
const sketcherMinimizerAtom* at,
set<sketcherMinimizerAtom*>& isAtomVisited) const
{
if (!at->getRings().empty())
return;
sketcherMinimizerBond* doubleBond = NULL;
vector<sketcherMinimizerAtom*> atomsToMirror;
for (unsigned int i = 0; i < at->bonds.size(); i++) {
if (at->bonds[i]->isStereo() &&
isAtomVisited.find(at->neighbors[i]) != isAtomVisited.end()) {
doubleBond = at->bonds[i];
} else {
atomsToMirror.push_back(at->neighbors[i]);
}
}
if (doubleBond == NULL)
return;
if (atomsToMirror.size() && doubleBond) {
sketcherMinimizerAtom* firstCIPNeighborStart =
doubleBond->startAtomCIPFirstNeighbor();
if (firstCIPNeighborStart == NULL)
return;
sketcherMinimizerAtom* firstCIPNeighborEnd =
doubleBond->endAtomCIPFirstNeighbor();
if (firstCIPNeighborEnd == NULL)
return;
if (!doubleBond->checkStereoChemistry()) {
for (sketcherMinimizerAtom* a : atomsToMirror) {
sketcherMinimizerAtom::mirrorCoordinates(a, doubleBond);
}
}
}
}
void CoordgenFragmentBuilder::buildFragment(
sketcherMinimizerFragment* fragment) const
{
buildRings(fragment);
buildNonRingAtoms(fragment);
CoordgenMinimizer::avoidInternalClashes(fragment);
fallbackIfNanCoordinates(fragment);
if (!fragment->getParent() && fragment->constrained)
rotateMainFragment(fragment);
// reset coordinates for fixed fragments:
if (fragment->fixed) {
fragment->setAllCoordinatesToTemplate();
}
}
void CoordgenFragmentBuilder::fallbackIfNanCoordinates(
sketcherMinimizerFragment* fragment) const
{
vector<sketcherMinimizerAtom*> fragmentAtoms = fragment->getAtoms();
if (CoordgenMinimizer::hasNaNCoordinates(fragmentAtoms) &&
CoordgenMinimizer::hasValid3DCoordinates(fragmentAtoms)) {
CoordgenMinimizer::fallbackOn3DCoordinates(fragmentAtoms);
}
}
void CoordgenFragmentBuilder::initializeFusedRingInformation(
sketcherMinimizerFragment* fragment) const
{
/*
for each ring in the fragment find out which other rings it shares atoms
with,
or if they are connected in other ways (i.e. via a double bond).
*/
if (fragment->getRings().size() < 2)
return;
vector<sketcherMinimizerAtom*> fragmentAtoms = fragment->getAtoms();
for (sketcherMinimizerAtom* atom : fragmentAtoms) {
if (atom->rings.size() > 1) {
for (unsigned int i = 0; i < atom->rings.size(); i++) {
for (unsigned int j = i + 1; j < atom->rings.size(); j++) {
bool found = false;
for (unsigned int r = 0;
r < atom->rings[i]->fusedWith.size(); r++) {
if (atom->rings[i]->fusedWith[r] == atom->rings[j]) {
atom->rings[i]->fusionAtoms[r].push_back(atom);
for (unsigned int rr = 0;
rr < atom->rings[j]->fusedWith.size(); rr++) {
if (atom->rings[j]->fusedWith[rr] ==
atom->rings[i]) {
atom->rings[j]->fusionAtoms[rr].push_back(
atom);
break;
}
}
found = true;
break;
}
}
if (!found) {
atom->rings[i]->fusedWith.push_back(atom->rings[j]);
atom->rings[j]->fusedWith.push_back(atom->rings[i]);
vector<sketcherMinimizerAtom*> ats;
ats.push_back(atom);
atom->rings[i]->fusionAtoms.push_back(ats);
atom->rings[j]->fusionAtoms.push_back(ats);
}
}
}
}
}
vector<sketcherMinimizerBond*> bonds = fragment->getBonds();
for (sketcherMinimizerBond* b : bonds) {
if ((b->bondOrder != 1) && (b->startAtom->rings.size()) &&
(b->endAtom->rings.size()) &&
(!sketcherMinimizer::sameRing(b->startAtom, b->endAtom))) {
for (sketcherMinimizerRing* r : b->startAtom->rings) {
for (sketcherMinimizerRing* r2 : b->endAtom->rings) {
r->fusedWith.push_back(r2);
r2->fusedWith.push_back(r);
vector<sketcherMinimizerAtom*> ats;
r->fusionAtoms.push_back(ats);
r2->fusionAtoms.push_back(ats);
r->fusionBonds.push_back(b);
r2->fusionBonds.push_back(b);
}
}
}
}
vector<sketcherMinimizerRing*> fragmentRings = fragment->getRings();
for (sketcherMinimizerRing* r : fragmentRings) {
for (unsigned int i = 0; i < r->fusedWith.size(); i++) {
vector<sketcherMinimizerAtom*> fusionAtoms = r->fusionAtoms[i];
if (fusionAtoms.size() <= 2)
continue;
// find an atom with only one neighbor in the vector
sketcherMinimizerAtom* startAtom = NULL;
for (sketcherMinimizerAtom* a : fusionAtoms) {
int counter = 0;
for (sketcherMinimizerAtom* n : a->neighbors) {
if (find(fusionAtoms.begin(), fusionAtoms.end(), n) !=
fusionAtoms.end()) {
counter++;
}
}
if (counter == 1) {
startAtom = a;
break;
}
}
assert(startAtom != NULL);
r->fusionAtoms[i] = orderChainOfAtoms(fusionAtoms, startAtom);
}
}
for (sketcherMinimizerRing* r : fragmentRings) {
for (unsigned int ii = 0; ii < r->fusedWith.size(); ii++) {
if (r->fusionAtoms[ii].size() > 2) {
for (sketcherMinimizerAtom* a :
r->fusionAtoms[ii]) { // looking for atoms with no
// neighbors in r but not in
// r->fusedWith[ii]
bool found = false;
for (sketcherMinimizerAtom* n : a->neighbors) {
int nn = 0;
for (sketcherMinimizerRing* nr : n->rings) {
if (nr == r)
nn++;
else if (nr == r->fusedWith[ii])
nn++;
}
if (nn == 1) {
found = true;
break;
}
}
if (!found) {
a->isSharedAndInner = true;
}
}
}
}
}
}
void CoordgenFragmentBuilder::simplifyRingSystem(
const vector<sketcherMinimizerRing*> allRings,
stack<sketcherMinimizerRing*>& sideRings,
vector<sketcherMinimizerRing*>& centralRings) const
{
map<sketcherMinimizerRing*, bool> alreadyChosenAsSide;
bool found = true;
/* simplify the structure by iteratively removing rings that have only one
fusion partner. These side rings will have coordinates generated
separately,
one at a time, starting from the closest to the core.
If we think of rings as atoms and fusion between two rings as bonds, this
is the same as stripping away non-ring atoms from a fragment to generate
coordinate
for the ring system first
*/
while (found) {
found = false;
for (sketcherMinimizerRing* r : allRings) {
if (alreadyChosenAsSide[r])
continue;
if (r->isMacrocycle())
continue;
int n = 0;
for (auto ringCounter = 0u; ringCounter < r->fusedWith.size();
++ringCounter) {
auto fusedRing = r->fusedWith.at(ringCounter);
if (!alreadyChosenAsSide[fusedRing]) {
n++;
if (fusedRing->isMacrocycle()) {
/* disable macrocycles */
n++;
}
if (r->fusionAtoms.at(ringCounter).size() > 3) {
/* disable rings that share too many atoms */
n++;
}
}
}
if (n == 1) /* ring is fused with only one ring,
it will be removed from the core of the system
and treated as a side ring
*/
{
sideRings.push(r);
found = true;
alreadyChosenAsSide[r] = true;
r->side = true;
}
}
}
/*
everything that is left after removing the side rings is the core of the
ring system
*/
for (sketcherMinimizerRing* ring : allRings) {
if (!alreadyChosenAsSide[ring])
centralRings.push_back(ring);
}
}
void CoordgenFragmentBuilder::buildRings(
sketcherMinimizerFragment* fragment) const
{
if (fragment->getRings().empty())
return;
initializeFusedRingInformation(fragment);
vector<sketcherMinimizerRing*> fragmentRings = fragment->getRings();
if (fragmentRings.size() == 1) {
buildRing(fragmentRings[0]);
} else { // more than 1 ring
stack<sketcherMinimizerRing*> sideRings;
vector<sketcherMinimizerRing*> centralRings;
simplifyRingSystem(fragmentRings, sideRings, centralRings);
generateCoordinatesCentralRings(centralRings);
generateCoordinatesSideRings(sideRings);
}
}