1 | /*
|
---|
2 | * Project: MoleCuilder
|
---|
3 | * Description: creates and alters molecular systems
|
---|
4 | * Copyright (C) 2010 University of Bonn. All rights reserved.
|
---|
5 | * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
|
---|
6 | */
|
---|
7 |
|
---|
8 | /*
|
---|
9 | * molecule_graph.cpp
|
---|
10 | *
|
---|
11 | * Created on: Oct 5, 2009
|
---|
12 | * Author: heber
|
---|
13 | */
|
---|
14 |
|
---|
15 | // include config.h
|
---|
16 | #ifdef HAVE_CONFIG_H
|
---|
17 | #include <config.h>
|
---|
18 | #endif
|
---|
19 |
|
---|
20 | #include "Helpers/MemDebug.hpp"
|
---|
21 |
|
---|
22 | #include "atom.hpp"
|
---|
23 | #include "bond.hpp"
|
---|
24 | #include "bondgraph.hpp"
|
---|
25 | #include "config.hpp"
|
---|
26 | #include "defs.hpp"
|
---|
27 | #include "element.hpp"
|
---|
28 | #include "Helpers/helpers.hpp"
|
---|
29 | #include "Helpers/Info.hpp"
|
---|
30 | #include "linkedcell.hpp"
|
---|
31 | #include "lists.hpp"
|
---|
32 | #include "Helpers/Verbose.hpp"
|
---|
33 | #include "Helpers/Log.hpp"
|
---|
34 | #include "molecule.hpp"
|
---|
35 | #include "World.hpp"
|
---|
36 | #include "Helpers/fast_functions.hpp"
|
---|
37 | #include "Helpers/Assert.hpp"
|
---|
38 | #include "LinearAlgebra/Matrix.hpp"
|
---|
39 | #include "Box.hpp"
|
---|
40 | #include "stackclass.hpp"
|
---|
41 |
|
---|
42 | struct BFSAccounting
|
---|
43 | {
|
---|
44 | atom **PredecessorList;
|
---|
45 | int *ShortestPathList;
|
---|
46 | enum Shading *ColorList;
|
---|
47 | class StackClass<atom *> *BFSStack;
|
---|
48 | class StackClass<atom *> *TouchedStack;
|
---|
49 | int AtomCount;
|
---|
50 | int BondOrder;
|
---|
51 | atom *Root;
|
---|
52 | bool BackStepping;
|
---|
53 | int CurrentGraphNr;
|
---|
54 | int ComponentNr;
|
---|
55 | };
|
---|
56 |
|
---|
57 | /** Accounting data for Depth First Search.
|
---|
58 | */
|
---|
59 | struct DFSAccounting
|
---|
60 | {
|
---|
61 | class StackClass<atom *> *AtomStack;
|
---|
62 | class StackClass<bond *> *BackEdgeStack;
|
---|
63 | int CurrentGraphNr;
|
---|
64 | int ComponentNumber;
|
---|
65 | atom *Root;
|
---|
66 | bool BackStepping;
|
---|
67 | };
|
---|
68 |
|
---|
69 | /************************************* Functions for class molecule *********************************/
|
---|
70 |
|
---|
71 | /** Creates an adjacency list of the molecule.
|
---|
72 | * We obtain an outside file with the indices of atoms which are bondmembers.
|
---|
73 | */
|
---|
74 | void molecule::CreateAdjacencyListFromDbondFile(ifstream *input)
|
---|
75 | {
|
---|
76 | Info FunctionInfo(__func__);
|
---|
77 | // 1 We will parse bonds out of the dbond file created by tremolo.
|
---|
78 | int atom1, atom2;
|
---|
79 | atom *Walker, *OtherWalker;
|
---|
80 | char line[MAXSTRINGSIZE];
|
---|
81 |
|
---|
82 | if (input->fail()) {
|
---|
83 | DoeLog(0) && (eLog() << Verbose(0) << "Opening of bond file failed \n");
|
---|
84 | performCriticalExit();
|
---|
85 | };
|
---|
86 | doCountAtoms();
|
---|
87 |
|
---|
88 | // skip header
|
---|
89 | input->getline(line,MAXSTRINGSIZE);
|
---|
90 | DoLog(1) && (Log() << Verbose(1) << "Scanning file ... \n");
|
---|
91 | while (!input->eof()) // Check whether we read everything already
|
---|
92 | {
|
---|
93 | input->getline(line,MAXSTRINGSIZE);
|
---|
94 | stringstream zeile(line);
|
---|
95 | zeile >> atom1;
|
---|
96 | zeile >> atom2;
|
---|
97 |
|
---|
98 | DoLog(2) && (Log() << Verbose(2) << "Looking for atoms " << atom1 << " and " << atom2 << "." << endl);
|
---|
99 | if (atom2 < atom1) //Sort indices of atoms in order
|
---|
100 | flip(atom1, atom2);
|
---|
101 | Walker = FindAtom(atom1);
|
---|
102 | ASSERT(Walker,"Could not find an atom with the ID given in dbond file");
|
---|
103 | OtherWalker = FindAtom(atom2);
|
---|
104 | ASSERT(OtherWalker,"Could not find an atom with the ID given in dbond file");
|
---|
105 | AddBond(Walker, OtherWalker); //Add the bond between the two atoms with respective indices.
|
---|
106 | }
|
---|
107 | }
|
---|
108 | ;
|
---|
109 |
|
---|
110 | /** Creates an adjacency list of the molecule.
|
---|
111 | * Generally, we use the CSD approach to bond recognition, that is the the distance
|
---|
112 | * between two atoms A and B must be within [Rcov(A)+Rcov(B)-t,Rcov(A)+Rcov(B)+t] with
|
---|
113 | * a threshold t = 0.4 Angstroem.
|
---|
114 | * To make it O(N log N) the function uses the linked-cell technique as follows:
|
---|
115 | * The procedure is step-wise:
|
---|
116 | * -# Remove every bond in list
|
---|
117 | * -# Count the atoms in the molecule with CountAtoms()
|
---|
118 | * -# partition cell into smaller linked cells of size \a bonddistance
|
---|
119 | * -# put each atom into its corresponding cell
|
---|
120 | * -# go through every cell, check the atoms therein against all possible bond partners in the 27 adjacent cells, add bond if true
|
---|
121 | * -# correct the bond degree iteratively (single->double->triple bond)
|
---|
122 | * -# finally print the bond list to \a *out if desired
|
---|
123 | * \param *out out stream for printing the matrix, NULL if no output
|
---|
124 | * \param bonddistance length of linked cells (i.e. maximum minimal length checked)
|
---|
125 | * \param IsAngstroem whether coordinate system is gauged to Angstroem or Bohr radii
|
---|
126 | * \param *minmaxdistance function to give upper and lower bound on whether particle is bonded to some other
|
---|
127 | * \param *BG BondGraph with the member function above or NULL, if just standard covalent should be used.
|
---|
128 | */
|
---|
129 | void molecule::CreateAdjacencyList(double bonddistance, bool IsAngstroem, void (BondGraph::*minmaxdistance)(BondedParticle * const , BondedParticle * const , double &, double &, bool), BondGraph *BG)
|
---|
130 | {
|
---|
131 | atom *Walker = NULL;
|
---|
132 | atom *OtherWalker = NULL;
|
---|
133 | int n[NDIM];
|
---|
134 | double MinDistance, MaxDistance;
|
---|
135 | LinkedCell *LC = NULL;
|
---|
136 | bool free_BG = false;
|
---|
137 | Box &domain = World::getInstance().getDomain();
|
---|
138 |
|
---|
139 | if (BG == NULL) {
|
---|
140 | BG = new BondGraph(IsAngstroem);
|
---|
141 | free_BG = true;
|
---|
142 | }
|
---|
143 |
|
---|
144 | BondDistance = bonddistance; // * ((IsAngstroem) ? 1. : 1./AtomicLengthToAngstroem);
|
---|
145 | DoLog(0) && (Log() << Verbose(0) << "Begin of CreateAdjacencyList." << endl);
|
---|
146 | // remove every bond from the list
|
---|
147 | for(molecule::iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner)
|
---|
148 | for(BondList::iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); !(*AtomRunner)->ListOfBonds.empty(); BondRunner = (*AtomRunner)->ListOfBonds.begin())
|
---|
149 | if ((*BondRunner)->leftatom == *AtomRunner)
|
---|
150 | delete((*BondRunner));
|
---|
151 | BondCount = 0;
|
---|
152 |
|
---|
153 | // count atoms in molecule = dimension of matrix (also give each unique name and continuous numbering)
|
---|
154 | DoLog(1) && (Log() << Verbose(1) << "AtomCount " << getAtomCount() << " and bonddistance is " << bonddistance << "." << endl);
|
---|
155 |
|
---|
156 | if ((getAtomCount() > 1) && (bonddistance > 0.1)) {
|
---|
157 | DoLog(2) && (Log() << Verbose(2) << "Creating Linked Cell structure ... " << endl);
|
---|
158 | LC = new LinkedCell(*this, bonddistance);
|
---|
159 |
|
---|
160 | // create a list to map Tesselpoint::nr to atom *
|
---|
161 | DoLog(2) && (Log() << Verbose(2) << "Creating TesselPoint to atom map ... " << endl);
|
---|
162 |
|
---|
163 | // set numbers for atoms that can later be used
|
---|
164 | int i=0;
|
---|
165 | for(internal_iterator iter = atoms.begin();iter!= atoms.end(); ++iter){
|
---|
166 | (*iter)->nr = i++;
|
---|
167 | }
|
---|
168 |
|
---|
169 | // 3a. go through every cell
|
---|
170 | DoLog(2) && (Log() << Verbose(2) << "Celling ... " << endl);
|
---|
171 | for (LC->n[0] = 0; LC->n[0] < LC->N[0]; LC->n[0]++)
|
---|
172 | for (LC->n[1] = 0; LC->n[1] < LC->N[1]; LC->n[1]++)
|
---|
173 | for (LC->n[2] = 0; LC->n[2] < LC->N[2]; LC->n[2]++) {
|
---|
174 | const LinkedCell::LinkedNodes *List = LC->GetCurrentCell();
|
---|
175 | // Log() << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << " containing " << List->size() << " points." << endl;
|
---|
176 | if (List != NULL) {
|
---|
177 | for (LinkedCell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
|
---|
178 | Walker = dynamic_cast<atom*>(*Runner);
|
---|
179 | ASSERT(Walker,"Tesselpoint that was not an atom retrieved from LinkedNode");
|
---|
180 | //Log() << Verbose(0) << "Current Atom is " << *Walker << "." << endl;
|
---|
181 | // 3c. check for possible bond between each atom in this and every one in the 27 cells
|
---|
182 | for (n[0] = -1; n[0] <= 1; n[0]++)
|
---|
183 | for (n[1] = -1; n[1] <= 1; n[1]++)
|
---|
184 | for (n[2] = -1; n[2] <= 1; n[2]++) {
|
---|
185 | const LinkedCell::LinkedNodes *OtherList = LC->GetRelativeToCurrentCell(n);
|
---|
186 | // Log() << Verbose(2) << "Current relative cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << " containing " << List->size() << " points." << endl;
|
---|
187 | if (OtherList != NULL) {
|
---|
188 | for (LinkedCell::LinkedNodes::const_iterator OtherRunner = OtherList->begin(); OtherRunner != OtherList->end(); OtherRunner++) {
|
---|
189 | if ((*OtherRunner)->nr > Walker->nr) {
|
---|
190 | OtherWalker = dynamic_cast<atom*>(*OtherRunner);
|
---|
191 | ASSERT(OtherWalker,"TesselPoint that was not an atom retrieved from LinkedNode");
|
---|
192 | //Log() << Verbose(1) << "Checking distance " << OtherWalker->x.PeriodicDistanceSquared(&(Walker->x), cell_size) << " against typical bond length of " << bonddistance*bonddistance << "." << endl;
|
---|
193 | (BG->*minmaxdistance)(Walker, OtherWalker, MinDistance, MaxDistance, IsAngstroem);
|
---|
194 | const double distance = domain.periodicDistanceSquared(OtherWalker->getPosition(),Walker->getPosition());
|
---|
195 | const bool status = (distance <= MaxDistance * MaxDistance) && (distance >= MinDistance * MinDistance);
|
---|
196 | // Log() << Verbose(1) << "MinDistance is " << MinDistance << " and MaxDistance is " << MaxDistance << "." << endl;
|
---|
197 | if (OtherWalker->father->nr > Walker->father->nr) {
|
---|
198 | if (status) { // create bond if distance is smaller
|
---|
199 | // Log() << Verbose(1) << "Adding Bond between " << *Walker << " and " << *OtherWalker << " in distance " << sqrt(distance) << "." << endl;
|
---|
200 | AddBond(Walker->father, OtherWalker->father, 1); // also increases molecule::BondCount
|
---|
201 | } else {
|
---|
202 | // Log() << Verbose(1) << "Not Adding: distance too great." << endl;
|
---|
203 | }
|
---|
204 | } else {
|
---|
205 | // Log() << Verbose(1) << "Not Adding: Wrong order of labels." << endl;
|
---|
206 | }
|
---|
207 | }
|
---|
208 | }
|
---|
209 | }
|
---|
210 | }
|
---|
211 | }
|
---|
212 | }
|
---|
213 | }
|
---|
214 | delete (LC);
|
---|
215 | DoLog(1) && (Log() << Verbose(1) << "I detected " << BondCount << " bonds in the molecule with distance " << BondDistance << "." << endl);
|
---|
216 |
|
---|
217 | // correct bond degree by comparing valence and bond degree
|
---|
218 | DoLog(2) && (Log() << Verbose(2) << "Correcting bond degree ... " << endl);
|
---|
219 | CorrectBondDegree();
|
---|
220 |
|
---|
221 | // output bonds for debugging (if bond chain list was correctly installed)
|
---|
222 | for_each(atoms.begin(),atoms.end(),mem_fun(&atom::OutputBondOfAtom));
|
---|
223 | } else
|
---|
224 | DoLog(1) && (Log() << Verbose(1) << "AtomCount is " << getAtomCount() << ", thus no bonds, no connections!." << endl);
|
---|
225 | DoLog(0) && (Log() << Verbose(0) << "End of CreateAdjacencyList." << endl);
|
---|
226 | if (free_BG)
|
---|
227 | delete(BG);
|
---|
228 | }
|
---|
229 | ;
|
---|
230 |
|
---|
231 | /** Checks for presence of bonds within atom list.
|
---|
232 | * TODO: more sophisticated check for bond structure (e.g. connected subgraph, ...)
|
---|
233 | * \return true - bonds present, false - no bonds
|
---|
234 | */
|
---|
235 | bool molecule::hasBondStructure() const
|
---|
236 | {
|
---|
237 | for(molecule::const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner)
|
---|
238 | if (!(*AtomRunner)->ListOfBonds.empty())
|
---|
239 | return true;
|
---|
240 | return false;
|
---|
241 | }
|
---|
242 |
|
---|
243 | /** Counts the number of present bonds.
|
---|
244 | * \return number of bonds
|
---|
245 | */
|
---|
246 | unsigned int molecule::CountBonds() const
|
---|
247 | {
|
---|
248 | unsigned int counter = 0;
|
---|
249 | for(molecule::const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner)
|
---|
250 | for(BondList::const_iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); BondRunner != (*AtomRunner)->ListOfBonds.end(); ++BondRunner)
|
---|
251 | if ((*BondRunner)->leftatom == *AtomRunner)
|
---|
252 | counter++;
|
---|
253 | return counter;
|
---|
254 | }
|
---|
255 |
|
---|
256 | /** Prints a list of all bonds to \a *out.
|
---|
257 | * \param output stream
|
---|
258 | */
|
---|
259 | void molecule::OutputBondsList() const
|
---|
260 | {
|
---|
261 | DoLog(1) && (Log() << Verbose(1) << endl << "From contents of bond chain list:");
|
---|
262 | for(molecule::const_iterator AtomRunner = molecule::begin(); AtomRunner != molecule::end(); ++AtomRunner)
|
---|
263 | for(BondList::const_iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); BondRunner != (*AtomRunner)->ListOfBonds.end(); ++BondRunner)
|
---|
264 | if ((*BondRunner)->leftatom == *AtomRunner) {
|
---|
265 | DoLog(0) && (Log() << Verbose(0) << *(*BondRunner) << "\t" << endl);
|
---|
266 | }
|
---|
267 | DoLog(0) && (Log() << Verbose(0) << endl);
|
---|
268 | }
|
---|
269 | ;
|
---|
270 |
|
---|
271 | /** correct bond degree by comparing valence and bond degree.
|
---|
272 | * correct Bond degree of each bond by checking both bond partners for a mismatch between valence and current sum of bond degrees,
|
---|
273 | * iteratively increase the one first where the other bond partner has the fewest number of bonds (i.e. in general bonds oxygene
|
---|
274 | * preferred over carbon bonds). Beforehand, we had picked the first mismatching partner, which lead to oxygenes with single instead of
|
---|
275 | * double bonds as was expected.
|
---|
276 | * \param *out output stream for debugging
|
---|
277 | * \return number of bonds that could not be corrected
|
---|
278 | */
|
---|
279 | int molecule::CorrectBondDegree() const
|
---|
280 | {
|
---|
281 | int No = 0, OldNo = -1;
|
---|
282 |
|
---|
283 | if (BondCount != 0) {
|
---|
284 | DoLog(1) && (Log() << Verbose(1) << "Correcting Bond degree of each bond ... " << endl);
|
---|
285 | do {
|
---|
286 | OldNo = No;
|
---|
287 | No=0;
|
---|
288 | BOOST_FOREACH(atom *atom,atoms){
|
---|
289 | No+=atom->CorrectBondDegree();
|
---|
290 | }
|
---|
291 | } while (OldNo != No);
|
---|
292 | DoLog(0) && (Log() << Verbose(0) << " done." << endl);
|
---|
293 | } else {
|
---|
294 | DoLog(1) && (Log() << Verbose(1) << "BondCount is " << BondCount << ", no bonds between any of the " << getAtomCount() << " atoms." << endl);
|
---|
295 | }
|
---|
296 | DoLog(0) && (Log() << Verbose(0) << No << " bonds could not be corrected." << endl);
|
---|
297 |
|
---|
298 | return (No);
|
---|
299 | }
|
---|
300 | ;
|
---|
301 |
|
---|
302 | /** Counts all cyclic bonds and returns their number.
|
---|
303 | * \note Hydrogen bonds can never by cyclic, thus no check for that
|
---|
304 | * \param *out output stream for debugging
|
---|
305 | * \return number opf cyclic bonds
|
---|
306 | */
|
---|
307 | int molecule::CountCyclicBonds()
|
---|
308 | {
|
---|
309 | NoCyclicBonds = 0;
|
---|
310 | int *MinimumRingSize = NULL;
|
---|
311 | MoleculeLeafClass *Subgraphs = NULL;
|
---|
312 | class StackClass<bond *> *BackEdgeStack = NULL;
|
---|
313 | for(molecule::iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner)
|
---|
314 | if ((!(*AtomRunner)->ListOfBonds.empty()) && ((*(*AtomRunner)->ListOfBonds.begin())->Type == Undetermined)) {
|
---|
315 | DoLog(0) && (Log() << Verbose(0) << "No Depth-First-Search analysis performed so far, calling ..." << endl);
|
---|
316 | Subgraphs = DepthFirstSearchAnalysis(BackEdgeStack);
|
---|
317 | while (Subgraphs->next != NULL) {
|
---|
318 | Subgraphs = Subgraphs->next;
|
---|
319 | delete (Subgraphs->previous);
|
---|
320 | }
|
---|
321 | delete (Subgraphs);
|
---|
322 | delete[] (MinimumRingSize);
|
---|
323 | break;
|
---|
324 | }
|
---|
325 | for(molecule::iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner)
|
---|
326 | for(BondList::iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); BondRunner != (*AtomRunner)->ListOfBonds.end(); ++BondRunner)
|
---|
327 | if ((*BondRunner)->leftatom == *AtomRunner)
|
---|
328 | if ((*BondRunner)->Cyclic)
|
---|
329 | NoCyclicBonds++;
|
---|
330 | delete (BackEdgeStack);
|
---|
331 | return NoCyclicBonds;
|
---|
332 | }
|
---|
333 | ;
|
---|
334 |
|
---|
335 | /** Returns Shading as a char string.
|
---|
336 | * \param color the Shading
|
---|
337 | * \return string of the flag
|
---|
338 | */
|
---|
339 | string molecule::GetColor(enum Shading color) const
|
---|
340 | {
|
---|
341 | switch (color) {
|
---|
342 | case white:
|
---|
343 | return "white";
|
---|
344 | break;
|
---|
345 | case lightgray:
|
---|
346 | return "lightgray";
|
---|
347 | break;
|
---|
348 | case darkgray:
|
---|
349 | return "darkgray";
|
---|
350 | break;
|
---|
351 | case black:
|
---|
352 | return "black";
|
---|
353 | break;
|
---|
354 | default:
|
---|
355 | return "uncolored";
|
---|
356 | break;
|
---|
357 | };
|
---|
358 | }
|
---|
359 | ;
|
---|
360 |
|
---|
361 | /** Sets atom::GraphNr and atom::LowpointNr to BFSAccounting::CurrentGraphNr.
|
---|
362 | * \param *out output stream for debugging
|
---|
363 | * \param *Walker current node
|
---|
364 | * \param &BFS structure with accounting data for BFS
|
---|
365 | */
|
---|
366 | void DepthFirstSearchAnalysis_SetWalkersGraphNr(atom *&Walker, struct DFSAccounting &DFS)
|
---|
367 | {
|
---|
368 | if (!DFS.BackStepping) { // if we don't just return from (8)
|
---|
369 | Walker->GraphNr = DFS.CurrentGraphNr;
|
---|
370 | Walker->LowpointNr = DFS.CurrentGraphNr;
|
---|
371 | DoLog(1) && (Log() << Verbose(1) << "Setting Walker[" << Walker->getName() << "]'s number to " << Walker->GraphNr << " with Lowpoint " << Walker->LowpointNr << "." << endl);
|
---|
372 | DFS.AtomStack->Push(Walker);
|
---|
373 | DFS.CurrentGraphNr++;
|
---|
374 | }
|
---|
375 | }
|
---|
376 | ;
|
---|
377 |
|
---|
378 | /** During DFS goes along unvisited bond and touches other atom.
|
---|
379 | * Sets bond::type, if
|
---|
380 | * -# BackEdge: set atom::LowpointNr and push on \a BackEdgeStack
|
---|
381 | * -# TreeEgde: set atom::Ancestor and continue with Walker along this edge
|
---|
382 | * Continue until molecule::FindNextUnused() finds no more unused bonds.
|
---|
383 | * \param *out output stream for debugging
|
---|
384 | * \param *mol molecule with atoms and finding unused bonds
|
---|
385 | * \param *&Binder current edge
|
---|
386 | * \param &DFS DFS accounting data
|
---|
387 | */
|
---|
388 | void DepthFirstSearchAnalysis_ProbeAlongUnusedBond(const molecule * const mol, atom *&Walker, bond *&Binder, struct DFSAccounting &DFS)
|
---|
389 | {
|
---|
390 | atom *OtherAtom = NULL;
|
---|
391 |
|
---|
392 | do { // (3) if Walker has no unused egdes, go to (5)
|
---|
393 | DFS.BackStepping = false; // reset backstepping flag for (8)
|
---|
394 | if (Binder == NULL) // if we don't just return from (11), Binder is already set to next unused
|
---|
395 | Binder = mol->FindNextUnused(Walker);
|
---|
396 | if (Binder == NULL)
|
---|
397 | break;
|
---|
398 | DoLog(2) && (Log() << Verbose(2) << "Current Unused Bond is " << *Binder << "." << endl);
|
---|
399 | // (4) Mark Binder used, ...
|
---|
400 | Binder->MarkUsed(black);
|
---|
401 | OtherAtom = Binder->GetOtherAtom(Walker);
|
---|
402 | DoLog(2) && (Log() << Verbose(2) << "(4) OtherAtom is " << OtherAtom->getName() << "." << endl);
|
---|
403 | if (OtherAtom->GraphNr != -1) {
|
---|
404 | // (4a) ... if "other" atom has been visited (GraphNr != 0), set lowpoint to minimum of both, go to (3)
|
---|
405 | Binder->Type = BackEdge;
|
---|
406 | DFS.BackEdgeStack->Push(Binder);
|
---|
407 | Walker->LowpointNr = (Walker->LowpointNr < OtherAtom->GraphNr) ? Walker->LowpointNr : OtherAtom->GraphNr;
|
---|
408 | DoLog(3) && (Log() << Verbose(3) << "(4a) Visited: Setting Lowpoint of Walker[" << Walker->getName() << "] to " << Walker->LowpointNr << "." << endl);
|
---|
409 | } else {
|
---|
410 | // (4b) ... otherwise set OtherAtom as Ancestor of Walker and Walker as OtherAtom, go to (2)
|
---|
411 | Binder->Type = TreeEdge;
|
---|
412 | OtherAtom->Ancestor = Walker;
|
---|
413 | Walker = OtherAtom;
|
---|
414 | DoLog(3) && (Log() << Verbose(3) << "(4b) Not Visited: OtherAtom[" << OtherAtom->getName() << "]'s Ancestor is now " << OtherAtom->Ancestor->getName() << ", Walker is OtherAtom " << OtherAtom->getName() << "." << endl);
|
---|
415 | break;
|
---|
416 | }
|
---|
417 | Binder = NULL;
|
---|
418 | } while (1); // (3)
|
---|
419 | }
|
---|
420 | ;
|
---|
421 |
|
---|
422 | /** Checks whether we have a new component.
|
---|
423 | * if atom::LowpointNr of \a *&Walker is greater than atom::GraphNr of its atom::Ancestor, we have a new component.
|
---|
424 | * Meaning that if we touch upon a node who suddenly has a smaller atom::LowpointNr than its ancestor, then we
|
---|
425 | * have a found a new branch in the graph tree.
|
---|
426 | * \param *out output stream for debugging
|
---|
427 | * \param *mol molecule with atoms and finding unused bonds
|
---|
428 | * \param *&Walker current node
|
---|
429 | * \param &DFS DFS accounting data
|
---|
430 | */
|
---|
431 | void DepthFirstSearchAnalysis_CheckForaNewComponent(const molecule * const mol, atom *&Walker, struct DFSAccounting &DFS, MoleculeLeafClass *&LeafWalker)
|
---|
432 | {
|
---|
433 | atom *OtherAtom = NULL;
|
---|
434 |
|
---|
435 | // (5) if Ancestor of Walker is ...
|
---|
436 | DoLog(1) && (Log() << Verbose(1) << "(5) Number of Walker[" << Walker->getName() << "]'s Ancestor[" << Walker->Ancestor->getName() << "] is " << Walker->Ancestor->GraphNr << "." << endl);
|
---|
437 |
|
---|
438 | if (Walker->Ancestor->GraphNr != DFS.Root->GraphNr) {
|
---|
439 | // (6) (Ancestor of Walker is not Root)
|
---|
440 | if (Walker->LowpointNr < Walker->Ancestor->GraphNr) {
|
---|
441 | // (6a) set Ancestor's Lowpoint number to minimum of of its Ancestor and itself, go to Step(8)
|
---|
442 | Walker->Ancestor->LowpointNr = (Walker->Ancestor->LowpointNr < Walker->LowpointNr) ? Walker->Ancestor->LowpointNr : Walker->LowpointNr;
|
---|
443 | DoLog(2) && (Log() << Verbose(2) << "(6) Setting Walker[" << Walker->getName() << "]'s Ancestor[" << Walker->Ancestor->getName() << "]'s Lowpoint to " << Walker->Ancestor->LowpointNr << "." << endl);
|
---|
444 | } else {
|
---|
445 | // (7) (Ancestor of Walker is a separating vertex, remove all from stack till Walker (including), these and Ancestor form a component
|
---|
446 | Walker->Ancestor->SeparationVertex = true;
|
---|
447 | DoLog(2) && (Log() << Verbose(2) << "(7) Walker[" << Walker->getName() << "]'s Ancestor[" << Walker->Ancestor->getName() << "]'s is a separating vertex, creating component." << endl);
|
---|
448 | mol->SetNextComponentNumber(Walker->Ancestor, DFS.ComponentNumber);
|
---|
449 | DoLog(3) && (Log() << Verbose(3) << "(7) Walker[" << Walker->getName() << "]'s Ancestor's Compont is " << DFS.ComponentNumber << "." << endl);
|
---|
450 | mol->SetNextComponentNumber(Walker, DFS.ComponentNumber);
|
---|
451 | DoLog(3) && (Log() << Verbose(3) << "(7) Walker[" << Walker->getName() << "]'s Compont is " << DFS.ComponentNumber << "." << endl);
|
---|
452 | do {
|
---|
453 | OtherAtom = DFS.AtomStack->PopLast();
|
---|
454 | LeafWalker->Leaf->AddCopyAtom(OtherAtom);
|
---|
455 | mol->SetNextComponentNumber(OtherAtom, DFS.ComponentNumber);
|
---|
456 | DoLog(3) && (Log() << Verbose(3) << "(7) Other[" << OtherAtom->getName() << "]'s Compont is " << DFS.ComponentNumber << "." << endl);
|
---|
457 | } while (OtherAtom != Walker);
|
---|
458 | DFS.ComponentNumber++;
|
---|
459 | }
|
---|
460 | // (8) Walker becomes its Ancestor, go to (3)
|
---|
461 | DoLog(2) && (Log() << Verbose(2) << "(8) Walker[" << Walker->getName() << "] is now its Ancestor " << Walker->Ancestor->getName() << ", backstepping. " << endl);
|
---|
462 | Walker = Walker->Ancestor;
|
---|
463 | DFS.BackStepping = true;
|
---|
464 | }
|
---|
465 | }
|
---|
466 | ;
|
---|
467 |
|
---|
468 | /** Cleans the root stack when we have found a component.
|
---|
469 | * If we are not DFSAccounting::BackStepping, then we clear the root stack by putting everything into a
|
---|
470 | * component down till we meet DFSAccounting::Root.
|
---|
471 | * \param *out output stream for debugging
|
---|
472 | * \param *mol molecule with atoms and finding unused bonds
|
---|
473 | * \param *&Walker current node
|
---|
474 | * \param *&Binder current edge
|
---|
475 | * \param &DFS DFS accounting data
|
---|
476 | */
|
---|
477 | void DepthFirstSearchAnalysis_CleanRootStackDownTillWalker(const molecule * const mol, atom *&Walker, bond *&Binder, struct DFSAccounting &DFS, MoleculeLeafClass *&LeafWalker)
|
---|
478 | {
|
---|
479 | atom *OtherAtom = NULL;
|
---|
480 |
|
---|
481 | if (!DFS.BackStepping) { // coming from (8) want to go to (3)
|
---|
482 | // (9) remove all from stack till Walker (including), these and Root form a component
|
---|
483 | //DFS.AtomStack->Output(out);
|
---|
484 | mol->SetNextComponentNumber(DFS.Root, DFS.ComponentNumber);
|
---|
485 | DoLog(3) && (Log() << Verbose(3) << "(9) Root[" << DFS.Root->getName() << "]'s Component is " << DFS.ComponentNumber << "." << endl);
|
---|
486 | mol->SetNextComponentNumber(Walker, DFS.ComponentNumber);
|
---|
487 | DoLog(3) && (Log() << Verbose(3) << "(9) Walker[" << Walker->getName() << "]'s Component is " << DFS.ComponentNumber << "." << endl);
|
---|
488 | do {
|
---|
489 | OtherAtom = DFS.AtomStack->PopLast();
|
---|
490 | LeafWalker->Leaf->AddCopyAtom(OtherAtom);
|
---|
491 | mol->SetNextComponentNumber(OtherAtom, DFS.ComponentNumber);
|
---|
492 | DoLog(3) && (Log() << Verbose(3) << "(7) Other[" << OtherAtom->getName() << "]'s Compont is " << DFS.ComponentNumber << "." << endl);
|
---|
493 | } while (OtherAtom != Walker);
|
---|
494 | DFS.ComponentNumber++;
|
---|
495 |
|
---|
496 | // (11) Root is separation vertex, set Walker to Root and go to (4)
|
---|
497 | Walker = DFS.Root;
|
---|
498 | Binder = mol->FindNextUnused(Walker);
|
---|
499 | DoLog(1) && (Log() << Verbose(1) << "(10) Walker is Root[" << DFS.Root->getName() << "], next Unused Bond is " << Binder << "." << endl);
|
---|
500 | if (Binder != NULL) { // Root is separation vertex
|
---|
501 | DoLog(1) && (Log() << Verbose(1) << "(11) Root is a separation vertex." << endl);
|
---|
502 | Walker->SeparationVertex = true;
|
---|
503 | }
|
---|
504 | }
|
---|
505 | }
|
---|
506 | ;
|
---|
507 |
|
---|
508 | /** Initializes DFSAccounting structure.
|
---|
509 | * \param *out output stream for debugging
|
---|
510 | * \param &DFS accounting structure to allocate
|
---|
511 | * \param *mol molecule with AtomCount, BondCount and all atoms
|
---|
512 | */
|
---|
513 | void DepthFirstSearchAnalysis_Init(struct DFSAccounting &DFS, const molecule * const mol)
|
---|
514 | {
|
---|
515 | DFS.AtomStack = new StackClass<atom *> (mol->getAtomCount());
|
---|
516 | DFS.CurrentGraphNr = 0;
|
---|
517 | DFS.ComponentNumber = 0;
|
---|
518 | DFS.BackStepping = false;
|
---|
519 | mol->ResetAllBondsToUnused();
|
---|
520 | DFS.BackEdgeStack->ClearStack();
|
---|
521 | }
|
---|
522 | ;
|
---|
523 |
|
---|
524 | /** Free's DFSAccounting structure.
|
---|
525 | * \param *out output stream for debugging
|
---|
526 | * \param &DFS accounting structure to free
|
---|
527 | */
|
---|
528 | void DepthFirstSearchAnalysis_Finalize(struct DFSAccounting &DFS)
|
---|
529 | {
|
---|
530 | delete (DFS.AtomStack);
|
---|
531 | // delete (DFS.BackEdgeStack); // DON'T free, see DepthFirstSearchAnalysis(), is returned as allocated
|
---|
532 | }
|
---|
533 | ;
|
---|
534 |
|
---|
535 | void molecule::init_DFS(struct DFSAccounting &DFS) const{
|
---|
536 | DepthFirstSearchAnalysis_Init(DFS, this);
|
---|
537 | for_each(atoms.begin(),atoms.end(),mem_fun(&atom::resetGraphNr));
|
---|
538 | for_each(atoms.begin(),atoms.end(),mem_fun(&atom::InitComponentNr));
|
---|
539 | }
|
---|
540 |
|
---|
541 | /** Performs a Depth-First search on this molecule.
|
---|
542 | * Marks bonds in molecule as cyclic, bridge, ... and atoms as
|
---|
543 | * articulations points, ...
|
---|
544 | * We use the algorithm from [Even, Graph Algorithms, p.62].
|
---|
545 | * \param *out output stream for debugging
|
---|
546 | * \param *&BackEdgeStack NULL pointer to StackClass with all the found back edges, allocated and filled on return
|
---|
547 | * \return list of each disconnected subgraph as an individual molecule class structure
|
---|
548 | */
|
---|
549 | MoleculeLeafClass * molecule::DepthFirstSearchAnalysis(class StackClass<bond *> *&BackEdgeStack) const
|
---|
550 | {
|
---|
551 | struct DFSAccounting DFS;
|
---|
552 | BackEdgeStack = new StackClass<bond *> (BondCount);
|
---|
553 | DFS.BackEdgeStack = BackEdgeStack;
|
---|
554 | MoleculeLeafClass *SubGraphs = new MoleculeLeafClass(NULL);
|
---|
555 | MoleculeLeafClass *LeafWalker = SubGraphs;
|
---|
556 | int OldGraphNr = 0;
|
---|
557 | atom *Walker = NULL;
|
---|
558 | bond *Binder = NULL;
|
---|
559 |
|
---|
560 | if (getAtomCount() == 0)
|
---|
561 | return SubGraphs;
|
---|
562 | DoLog(0) && (Log() << Verbose(0) << "Begin of DepthFirstSearchAnalysis" << endl);
|
---|
563 | init_DFS(DFS);
|
---|
564 |
|
---|
565 | for (molecule::const_iterator iter = begin(); iter != end();) {
|
---|
566 | DFS.Root = *iter;
|
---|
567 | // (1) mark all edges unused, empty stack, set atom->GraphNr = -1 for all
|
---|
568 | DFS.AtomStack->ClearStack();
|
---|
569 |
|
---|
570 | // put into new subgraph molecule and add this to list of subgraphs
|
---|
571 | LeafWalker = new MoleculeLeafClass(LeafWalker);
|
---|
572 | LeafWalker->Leaf = World::getInstance().createMolecule();
|
---|
573 | LeafWalker->Leaf->AddCopyAtom(DFS.Root);
|
---|
574 |
|
---|
575 | OldGraphNr = DFS.CurrentGraphNr;
|
---|
576 | Walker = DFS.Root;
|
---|
577 | do { // (10)
|
---|
578 | do { // (2) set number and Lowpoint of Atom to i, increase i, push current atom
|
---|
579 | DepthFirstSearchAnalysis_SetWalkersGraphNr(Walker, DFS);
|
---|
580 |
|
---|
581 | DepthFirstSearchAnalysis_ProbeAlongUnusedBond(this, Walker, Binder, DFS);
|
---|
582 |
|
---|
583 | if (Binder == NULL) {
|
---|
584 | DoLog(2) && (Log() << Verbose(2) << "No more Unused Bonds." << endl);
|
---|
585 | break;
|
---|
586 | } else
|
---|
587 | Binder = NULL;
|
---|
588 | } while (1); // (2)
|
---|
589 |
|
---|
590 | // if we came from backstepping, yet there were no more unused bonds, we end up here with no Ancestor, because Walker is Root! Then we are finished!
|
---|
591 | if ((Walker == DFS.Root) && (Binder == NULL))
|
---|
592 | break;
|
---|
593 |
|
---|
594 | DepthFirstSearchAnalysis_CheckForaNewComponent(this, Walker, DFS, LeafWalker);
|
---|
595 |
|
---|
596 | DepthFirstSearchAnalysis_CleanRootStackDownTillWalker(this, Walker, Binder, DFS, LeafWalker);
|
---|
597 |
|
---|
598 | } while ((DFS.BackStepping) || (Binder != NULL)); // (10) halt only if Root has no unused edges
|
---|
599 |
|
---|
600 | // From OldGraphNr to CurrentGraphNr ranges an disconnected subgraph
|
---|
601 | DoLog(0) && (Log() << Verbose(0) << "Disconnected subgraph ranges from " << OldGraphNr << " to " << DFS.CurrentGraphNr << "." << endl);
|
---|
602 | LeafWalker->Leaf->Output((ofstream *)&(Log() << Verbose(0)));
|
---|
603 | DoLog(0) && (Log() << Verbose(0) << endl);
|
---|
604 |
|
---|
605 | // step on to next root
|
---|
606 | while ((iter != end()) && ((*iter)->GraphNr != -1)) {
|
---|
607 | //Log() << Verbose(1) << "Current next subgraph root candidate is " << (*iter)->Name << "." << endl;
|
---|
608 | if ((*iter)->GraphNr != -1) // if already discovered, step on
|
---|
609 | iter++;
|
---|
610 | }
|
---|
611 | }
|
---|
612 | // set cyclic bond criterium on "same LP" basis
|
---|
613 | CyclicBondAnalysis();
|
---|
614 |
|
---|
615 | OutputGraphInfoPerAtom();
|
---|
616 |
|
---|
617 | OutputGraphInfoPerBond();
|
---|
618 |
|
---|
619 | // free all and exit
|
---|
620 | DepthFirstSearchAnalysis_Finalize(DFS);
|
---|
621 | DoLog(0) && (Log() << Verbose(0) << "End of DepthFirstSearchAnalysis" << endl);
|
---|
622 | return SubGraphs;
|
---|
623 | }
|
---|
624 | ;
|
---|
625 |
|
---|
626 | /** Scans through all bonds and set bond::Cyclic to true where atom::LowpointNr of both ends is equal: LP criterion.
|
---|
627 | */
|
---|
628 | void molecule::CyclicBondAnalysis() const
|
---|
629 | {
|
---|
630 | NoCyclicBonds = 0;
|
---|
631 | for(molecule::const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner)
|
---|
632 | for(BondList::const_iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); BondRunner != (*AtomRunner)->ListOfBonds.end(); ++BondRunner)
|
---|
633 | if ((*BondRunner)->leftatom == *AtomRunner)
|
---|
634 | if ((*BondRunner)->rightatom->LowpointNr == (*BondRunner)->leftatom->LowpointNr) { // cyclic ??
|
---|
635 | (*BondRunner)->Cyclic = true;
|
---|
636 | NoCyclicBonds++;
|
---|
637 | }
|
---|
638 | }
|
---|
639 | ;
|
---|
640 |
|
---|
641 | /** Output graph information per atom.
|
---|
642 | * \param *out output stream
|
---|
643 | */
|
---|
644 | void molecule::OutputGraphInfoPerAtom() const
|
---|
645 | {
|
---|
646 | DoLog(1) && (Log() << Verbose(1) << "Final graph info for each atom is:" << endl);
|
---|
647 | for_each(atoms.begin(),atoms.end(),mem_fun(&atom::OutputGraphInfo));
|
---|
648 | }
|
---|
649 | ;
|
---|
650 |
|
---|
651 | /** Output graph information per bond.
|
---|
652 | * \param *out output stream
|
---|
653 | */
|
---|
654 | void molecule::OutputGraphInfoPerBond() const
|
---|
655 | {
|
---|
656 | bond *Binder = NULL;
|
---|
657 | DoLog(1) && (Log() << Verbose(1) << "Final graph info for each bond is:" << endl);
|
---|
658 | for(molecule::const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner)
|
---|
659 | for(BondList::const_iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); BondRunner != (*AtomRunner)->ListOfBonds.end(); ++BondRunner)
|
---|
660 | if ((*BondRunner)->leftatom == *AtomRunner) {
|
---|
661 | Binder = *BondRunner;
|
---|
662 | if (DoLog(2)) {
|
---|
663 | ostream &out = (Log() << Verbose(2));
|
---|
664 | out << ((Binder->Type == TreeEdge) ? "TreeEdge " : "BackEdge ") << *Binder << ": <";
|
---|
665 | out << ((Binder->leftatom->SeparationVertex) ? "SP," : "") << "L" << Binder->leftatom->LowpointNr << " G" << Binder->leftatom->GraphNr << " Comp.";
|
---|
666 | Binder->leftatom->OutputComponentNumber(&out);
|
---|
667 | out << " === ";
|
---|
668 | out << ((Binder->rightatom->SeparationVertex) ? "SP," : "") << "L" << Binder->rightatom->LowpointNr << " G" << Binder->rightatom->GraphNr << " Comp.";
|
---|
669 | Binder->rightatom->OutputComponentNumber(&out);
|
---|
670 | out << ">." << endl;
|
---|
671 | }
|
---|
672 | if (Binder->Cyclic) // cyclic ??
|
---|
673 | DoLog(3) && (Log() << Verbose(3) << "Lowpoint at each side are equal: CYCLIC!" << endl);
|
---|
674 | }
|
---|
675 | }
|
---|
676 | ;
|
---|
677 |
|
---|
678 | /** Initialise each vertex as white with no predecessor, empty queue, color Root lightgray.
|
---|
679 | * \param *out output stream for debugging
|
---|
680 | * \param &BFS accounting structure
|
---|
681 | * \param AtomCount number of entries in the array to allocate
|
---|
682 | */
|
---|
683 | void InitializeBFSAccounting(struct BFSAccounting &BFS, int AtomCount)
|
---|
684 | {
|
---|
685 | BFS.AtomCount = AtomCount;
|
---|
686 | BFS.PredecessorList = new atom*[AtomCount];
|
---|
687 | BFS.ShortestPathList = new int[AtomCount];
|
---|
688 | BFS.ColorList = new enum Shading[AtomCount];
|
---|
689 | BFS.BFSStack = new StackClass<atom *> (AtomCount);
|
---|
690 | BFS.TouchedStack = new StackClass<atom *> (AtomCount);
|
---|
691 |
|
---|
692 | for (int i = AtomCount; i--;) {
|
---|
693 | BFS.ShortestPathList[i] = -1;
|
---|
694 | BFS.PredecessorList[i] = 0;
|
---|
695 | BFS.ColorList[i] = white;
|
---|
696 | }
|
---|
697 | };
|
---|
698 |
|
---|
699 | /** Free's accounting structure.
|
---|
700 | * \param *out output stream for debugging
|
---|
701 | * \param &BFS accounting structure
|
---|
702 | */
|
---|
703 | void FinalizeBFSAccounting(struct BFSAccounting &BFS)
|
---|
704 | {
|
---|
705 | delete[](BFS.PredecessorList);
|
---|
706 | delete[](BFS.ShortestPathList);
|
---|
707 | delete[](BFS.ColorList);
|
---|
708 | delete (BFS.BFSStack);
|
---|
709 | delete (BFS.TouchedStack);
|
---|
710 | BFS.AtomCount = 0;
|
---|
711 | };
|
---|
712 |
|
---|
713 | /** Clean the accounting structure.
|
---|
714 | * \param *out output stream for debugging
|
---|
715 | * \param &BFS accounting structure
|
---|
716 | */
|
---|
717 | void CleanBFSAccounting(struct BFSAccounting &BFS)
|
---|
718 | {
|
---|
719 | atom *Walker = NULL;
|
---|
720 | while (!BFS.TouchedStack->IsEmpty()) {
|
---|
721 | Walker = BFS.TouchedStack->PopFirst();
|
---|
722 | BFS.PredecessorList[Walker->nr] = NULL;
|
---|
723 | BFS.ShortestPathList[Walker->nr] = -1;
|
---|
724 | BFS.ColorList[Walker->nr] = white;
|
---|
725 | }
|
---|
726 | };
|
---|
727 |
|
---|
728 | /** Resets shortest path list and BFSStack.
|
---|
729 | * \param *out output stream for debugging
|
---|
730 | * \param *&Walker current node, pushed onto BFSAccounting::BFSStack and BFSAccounting::TouchedStack
|
---|
731 | * \param &BFS accounting structure
|
---|
732 | */
|
---|
733 | void ResetBFSAccounting(atom *&Walker, struct BFSAccounting &BFS)
|
---|
734 | {
|
---|
735 | BFS.ShortestPathList[Walker->nr] = 0;
|
---|
736 | BFS.BFSStack->ClearStack(); // start with empty BFS stack
|
---|
737 | BFS.BFSStack->Push(Walker);
|
---|
738 | BFS.TouchedStack->Push(Walker);
|
---|
739 | };
|
---|
740 |
|
---|
741 | /** Performs a BFS from \a *Root, trying to find the same node and hence a cycle.
|
---|
742 | * \param *out output stream for debugging
|
---|
743 | * \param *&BackEdge the edge from root that we don't want to move along
|
---|
744 | * \param &BFS accounting structure
|
---|
745 | */
|
---|
746 | void CyclicStructureAnalysis_CyclicBFSFromRootToRoot(bond *&BackEdge, struct BFSAccounting &BFS)
|
---|
747 | {
|
---|
748 | atom *Walker = NULL;
|
---|
749 | atom *OtherAtom = NULL;
|
---|
750 | do { // look for Root
|
---|
751 | Walker = BFS.BFSStack->PopFirst();
|
---|
752 | DoLog(2) && (Log() << Verbose(2) << "Current Walker is " << *Walker << ", we look for SP to Root " << *BFS.Root << "." << endl);
|
---|
753 | for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
|
---|
754 | if ((*Runner) != BackEdge) { // only walk along DFS spanning tree (otherwise we always find SP of one being backedge Binder)
|
---|
755 | OtherAtom = (*Runner)->GetOtherAtom(Walker);
|
---|
756 | #ifdef ADDHYDROGEN
|
---|
757 | if (OtherAtom->getType()->getAtomicNumber() != 1) {
|
---|
758 | #endif
|
---|
759 | DoLog(2) && (Log() << Verbose(2) << "Current OtherAtom is: " << OtherAtom->getName() << " for bond " << *(*Runner) << "." << endl);
|
---|
760 | if (BFS.ColorList[OtherAtom->nr] == white) {
|
---|
761 | BFS.TouchedStack->Push(OtherAtom);
|
---|
762 | BFS.ColorList[OtherAtom->nr] = lightgray;
|
---|
763 | BFS.PredecessorList[OtherAtom->nr] = Walker; // Walker is the predecessor
|
---|
764 | BFS.ShortestPathList[OtherAtom->nr] = BFS.ShortestPathList[Walker->nr] + 1;
|
---|
765 | DoLog(2) && (Log() << Verbose(2) << "Coloring OtherAtom " << OtherAtom->getName() << " lightgray, its predecessor is " << Walker->getName() << " and its Shortest Path is " << BFS.ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl);
|
---|
766 | //if (BFS.ShortestPathList[OtherAtom->nr] < MinimumRingSize[Walker->GetTrueFather()->nr]) { // Check for maximum distance
|
---|
767 | DoLog(3) && (Log() << Verbose(3) << "Putting OtherAtom into queue." << endl);
|
---|
768 | BFS.BFSStack->Push(OtherAtom);
|
---|
769 | //}
|
---|
770 | } else {
|
---|
771 | DoLog(3) && (Log() << Verbose(3) << "Not Adding, has already been visited." << endl);
|
---|
772 | }
|
---|
773 | if (OtherAtom == BFS.Root)
|
---|
774 | break;
|
---|
775 | #ifdef ADDHYDROGEN
|
---|
776 | } else {
|
---|
777 | DoLog(2) && (Log() << Verbose(2) << "Skipping hydrogen atom " << *OtherAtom << "." << endl);
|
---|
778 | BFS.ColorList[OtherAtom->nr] = black;
|
---|
779 | }
|
---|
780 | #endif
|
---|
781 | } else {
|
---|
782 | DoLog(2) && (Log() << Verbose(2) << "Bond " << *(*Runner) << " not Visiting, is the back edge." << endl);
|
---|
783 | }
|
---|
784 | }
|
---|
785 | BFS.ColorList[Walker->nr] = black;
|
---|
786 | DoLog(1) && (Log() << Verbose(1) << "Coloring Walker " << Walker->getName() << " black." << endl);
|
---|
787 | if (OtherAtom == BFS.Root) { // if we have found the root, check whether this cycle wasn't already found beforehand
|
---|
788 | // step through predecessor list
|
---|
789 | while (OtherAtom != BackEdge->rightatom) {
|
---|
790 | if (!OtherAtom->GetTrueFather()->IsCyclic) // if one bond in the loop is not marked as cyclic, we haven't found this cycle yet
|
---|
791 | break;
|
---|
792 | else
|
---|
793 | OtherAtom = BFS.PredecessorList[OtherAtom->nr];
|
---|
794 | }
|
---|
795 | if (OtherAtom == BackEdge->rightatom) { // if each atom in found cycle is cyclic, loop's been found before already
|
---|
796 | DoLog(3) && (Log() << Verbose(3) << "This cycle was already found before, skipping and removing seeker from search." << endl);
|
---|
797 | do {
|
---|
798 | OtherAtom = BFS.TouchedStack->PopLast();
|
---|
799 | if (BFS.PredecessorList[OtherAtom->nr] == Walker) {
|
---|
800 | DoLog(4) && (Log() << Verbose(4) << "Removing " << *OtherAtom << " from lists and stacks." << endl);
|
---|
801 | BFS.PredecessorList[OtherAtom->nr] = NULL;
|
---|
802 | BFS.ShortestPathList[OtherAtom->nr] = -1;
|
---|
803 | BFS.ColorList[OtherAtom->nr] = white;
|
---|
804 | BFS.BFSStack->RemoveItem(OtherAtom);
|
---|
805 | }
|
---|
806 | } while ((!BFS.TouchedStack->IsEmpty()) && (BFS.PredecessorList[OtherAtom->nr] == NULL));
|
---|
807 | BFS.TouchedStack->Push(OtherAtom); // last was wrongly popped
|
---|
808 | OtherAtom = BackEdge->rightatom; // set to not Root
|
---|
809 | } else
|
---|
810 | OtherAtom = BFS.Root;
|
---|
811 | }
|
---|
812 | } while ((!BFS.BFSStack->IsEmpty()) && (OtherAtom != BFS.Root) && (OtherAtom != NULL)); // || (ShortestPathList[OtherAtom->nr] < MinimumRingSize[Walker->GetTrueFather()->nr])));
|
---|
813 | };
|
---|
814 |
|
---|
815 | /** Climb back the BFSAccounting::PredecessorList and find cycle members.
|
---|
816 | * \param *out output stream for debugging
|
---|
817 | * \param *&OtherAtom
|
---|
818 | * \param *&BackEdge denotes the edge we did not want to travel along when doing CyclicBFSFromRootToRoot()
|
---|
819 | * \param &BFS accounting structure
|
---|
820 | * \param *&MinimumRingSize minimum distance from this node possible without encountering oneself, set on return for each atom
|
---|
821 | * \param &MinRingSize global minimum distance from one node without encountering oneself, set on return
|
---|
822 | */
|
---|
823 | void CyclicStructureAnalysis_RetrieveCycleMembers(atom *&OtherAtom, bond *&BackEdge, struct BFSAccounting &BFS, int *&MinimumRingSize, int &MinRingSize)
|
---|
824 | {
|
---|
825 | atom *Walker = NULL;
|
---|
826 | int NumCycles = 0;
|
---|
827 | int RingSize = -1;
|
---|
828 |
|
---|
829 | if (OtherAtom == BFS.Root) {
|
---|
830 | // now climb back the predecessor list and thus find the cycle members
|
---|
831 | NumCycles++;
|
---|
832 | RingSize = 1;
|
---|
833 | BFS.Root->GetTrueFather()->IsCyclic = true;
|
---|
834 | DoLog(1) && (Log() << Verbose(1) << "Found ring contains: ");
|
---|
835 | Walker = BFS.Root;
|
---|
836 | while (Walker != BackEdge->rightatom) {
|
---|
837 | DoLog(0) && (Log() << Verbose(0) << Walker->getName() << " <-> ");
|
---|
838 | Walker = BFS.PredecessorList[Walker->nr];
|
---|
839 | Walker->GetTrueFather()->IsCyclic = true;
|
---|
840 | RingSize++;
|
---|
841 | }
|
---|
842 | DoLog(0) && (Log() << Verbose(0) << Walker->getName() << " with a length of " << RingSize << "." << endl << endl);
|
---|
843 | // walk through all and set MinimumRingSize
|
---|
844 | Walker = BFS.Root;
|
---|
845 | MinimumRingSize[Walker->GetTrueFather()->nr] = RingSize;
|
---|
846 | while (Walker != BackEdge->rightatom) {
|
---|
847 | Walker = BFS.PredecessorList[Walker->nr];
|
---|
848 | if (RingSize < MinimumRingSize[Walker->GetTrueFather()->nr])
|
---|
849 | MinimumRingSize[Walker->GetTrueFather()->nr] = RingSize;
|
---|
850 | }
|
---|
851 | if ((RingSize < MinRingSize) || (MinRingSize == -1))
|
---|
852 | MinRingSize = RingSize;
|
---|
853 | } else {
|
---|
854 | DoLog(1) && (Log() << Verbose(1) << "No ring containing " << *BFS.Root << " with length equal to or smaller than " << MinimumRingSize[BFS.Root->GetTrueFather()->nr] << " found." << endl);
|
---|
855 | }
|
---|
856 | };
|
---|
857 |
|
---|
858 | /** From a given node performs a BFS to touch the next cycle, for whose nodes \a *&MinimumRingSize is set and set it accordingly.
|
---|
859 | * \param *out output stream for debugging
|
---|
860 | * \param *&Root node to look for closest cycle from, i.e. \a *&MinimumRingSize is set for this node
|
---|
861 | * \param *&MinimumRingSize minimum distance from this node possible without encountering oneself, set on return for each atom
|
---|
862 | * \param AtomCount number of nodes in graph
|
---|
863 | */
|
---|
864 | void CyclicStructureAnalysis_BFSToNextCycle(atom *&Root, atom *&Walker, int *&MinimumRingSize, int AtomCount)
|
---|
865 | {
|
---|
866 | struct BFSAccounting BFS;
|
---|
867 | atom *OtherAtom = Walker;
|
---|
868 |
|
---|
869 | InitializeBFSAccounting(BFS, AtomCount);
|
---|
870 |
|
---|
871 | ResetBFSAccounting(Walker, BFS);
|
---|
872 | while (OtherAtom != NULL) { // look for Root
|
---|
873 | Walker = BFS.BFSStack->PopFirst();
|
---|
874 | //Log() << Verbose(2) << "Current Walker is " << *Walker << ", we look for SP to Root " << *Root << "." << endl;
|
---|
875 | for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
|
---|
876 | // "removed (*Runner) != BackEdge) || " from next if, is u
|
---|
877 | if ((Walker->ListOfBonds.size() == 1)) { // only walk along DFS spanning tree (otherwise we always find SP of 1 being backedge Binder), but terminal hydrogens may be connected via backedge, hence extra check
|
---|
878 | OtherAtom = (*Runner)->GetOtherAtom(Walker);
|
---|
879 | //Log() << Verbose(2) << "Current OtherAtom is: " << OtherAtom->Name << " for bond " << *Binder << "." << endl;
|
---|
880 | if (BFS.ColorList[OtherAtom->nr] == white) {
|
---|
881 | BFS.TouchedStack->Push(OtherAtom);
|
---|
882 | BFS.ColorList[OtherAtom->nr] = lightgray;
|
---|
883 | BFS.PredecessorList[OtherAtom->nr] = Walker; // Walker is the predecessor
|
---|
884 | BFS.ShortestPathList[OtherAtom->nr] = BFS.ShortestPathList[Walker->nr] + 1;
|
---|
885 | //Log() << Verbose(2) << "Coloring OtherAtom " << OtherAtom->Name << " lightgray, its predecessor is " << Walker->Name << " and its Shortest Path is " << ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl;
|
---|
886 | if (OtherAtom->GetTrueFather()->IsCyclic) { // if the other atom is connected to a ring
|
---|
887 | MinimumRingSize[Root->GetTrueFather()->nr] = BFS.ShortestPathList[OtherAtom->nr] + MinimumRingSize[OtherAtom->GetTrueFather()->nr];
|
---|
888 | OtherAtom = NULL; //break;
|
---|
889 | break;
|
---|
890 | } else
|
---|
891 | BFS.BFSStack->Push(OtherAtom);
|
---|
892 | } else {
|
---|
893 | //Log() << Verbose(3) << "Not Adding, has already been visited." << endl;
|
---|
894 | }
|
---|
895 | } else {
|
---|
896 | //Log() << Verbose(3) << "Not Visiting, is a back edge." << endl;
|
---|
897 | }
|
---|
898 | }
|
---|
899 | BFS.ColorList[Walker->nr] = black;
|
---|
900 | //Log() << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl;
|
---|
901 | }
|
---|
902 | //CleanAccountingLists(TouchedStack, PredecessorList, ShortestPathList, ColorList);
|
---|
903 |
|
---|
904 | FinalizeBFSAccounting(BFS);
|
---|
905 | }
|
---|
906 | ;
|
---|
907 |
|
---|
908 | /** All nodes that are not in cycles get assigned a \a *&MinimumRingSizeby BFS to next cycle.
|
---|
909 | * \param *out output stream for debugging
|
---|
910 | * \param *&MinimumRingSize array with minimum distance without encountering onself for each atom
|
---|
911 | * \param &MinRingSize global minium distance
|
---|
912 | * \param &NumCyles number of cycles in graph
|
---|
913 | * \param *mol molecule with atoms
|
---|
914 | */
|
---|
915 | void CyclicStructureAnalysis_AssignRingSizetoNonCycleMembers(int *&MinimumRingSize, int &MinRingSize, int &NumCycles, const molecule * const mol)
|
---|
916 | {
|
---|
917 | atom *Root = NULL;
|
---|
918 | atom *Walker = NULL;
|
---|
919 | if (MinRingSize != -1) { // if rings are present
|
---|
920 | // go over all atoms
|
---|
921 | for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
|
---|
922 | Root = *iter;
|
---|
923 |
|
---|
924 | if (MinimumRingSize[Root->GetTrueFather()->nr] == mol->getAtomCount()) { // check whether MinimumRingSize is set, if not BFS to next where it is
|
---|
925 | Walker = Root;
|
---|
926 |
|
---|
927 | //Log() << Verbose(1) << "---------------------------------------------------------------------------------------------------------" << endl;
|
---|
928 | CyclicStructureAnalysis_BFSToNextCycle(Root, Walker, MinimumRingSize, mol->getAtomCount());
|
---|
929 |
|
---|
930 | }
|
---|
931 | DoLog(1) && (Log() << Verbose(1) << "Minimum ring size of " << *Root << " is " << MinimumRingSize[Root->GetTrueFather()->nr] << "." << endl);
|
---|
932 | }
|
---|
933 | DoLog(1) && (Log() << Verbose(1) << "Minimum ring size is " << MinRingSize << ", over " << NumCycles << " cycles total." << endl);
|
---|
934 | } else
|
---|
935 | DoLog(1) && (Log() << Verbose(1) << "No rings were detected in the molecular structure." << endl);
|
---|
936 | }
|
---|
937 | ;
|
---|
938 |
|
---|
939 | /** Analyses the cycles found and returns minimum of all cycle lengths.
|
---|
940 | * We begin with a list of Back edges found during DepthFirstSearchAnalysis(). We go through this list - one end is the Root,
|
---|
941 | * the other our initial Walker - and do a Breadth First Search for the Root. We mark down each Predecessor and as soon as
|
---|
942 | * we have found the Root via BFS, we may climb back the closed cycle via the Predecessors. Thereby we mark atoms and bonds
|
---|
943 | * as cyclic and print out the cycles.
|
---|
944 | * \param *out output stream for debugging
|
---|
945 | * \param *BackEdgeStack stack with all back edges found during DFS scan. Beware: This stack contains the bonds from the total molecule, not from the subgraph!
|
---|
946 | * \param *&MinimumRingSize contains smallest ring size in molecular structure on return or -1 if no rings were found, if set is maximum search distance
|
---|
947 | * \todo BFS from the not-same-LP to find back to starting point of tributary cycle over more than one bond
|
---|
948 | */
|
---|
949 | void molecule::CyclicStructureAnalysis(class StackClass<bond *> * BackEdgeStack, int *&MinimumRingSize) const
|
---|
950 | {
|
---|
951 | struct BFSAccounting BFS;
|
---|
952 | atom *Walker = NULL;
|
---|
953 | atom *OtherAtom = NULL;
|
---|
954 | bond *BackEdge = NULL;
|
---|
955 | int NumCycles = 0;
|
---|
956 | int MinRingSize = -1;
|
---|
957 |
|
---|
958 | InitializeBFSAccounting(BFS, getAtomCount());
|
---|
959 |
|
---|
960 | //Log() << Verbose(1) << "Back edge list - ";
|
---|
961 | //BackEdgeStack->Output(out);
|
---|
962 |
|
---|
963 | DoLog(1) && (Log() << Verbose(1) << "Analysing cycles ... " << endl);
|
---|
964 | NumCycles = 0;
|
---|
965 | while (!BackEdgeStack->IsEmpty()) {
|
---|
966 | BackEdge = BackEdgeStack->PopFirst();
|
---|
967 | // this is the target
|
---|
968 | BFS.Root = BackEdge->leftatom;
|
---|
969 | // this is the source point
|
---|
970 | Walker = BackEdge->rightatom;
|
---|
971 |
|
---|
972 | ResetBFSAccounting(Walker, BFS);
|
---|
973 |
|
---|
974 | DoLog(1) && (Log() << Verbose(1) << "---------------------------------------------------------------------------------------------------------" << endl);
|
---|
975 | OtherAtom = NULL;
|
---|
976 | CyclicStructureAnalysis_CyclicBFSFromRootToRoot(BackEdge, BFS);
|
---|
977 |
|
---|
978 | CyclicStructureAnalysis_RetrieveCycleMembers(OtherAtom, BackEdge, BFS, MinimumRingSize, MinRingSize);
|
---|
979 |
|
---|
980 | CleanBFSAccounting(BFS);
|
---|
981 | }
|
---|
982 | FinalizeBFSAccounting(BFS);
|
---|
983 |
|
---|
984 | CyclicStructureAnalysis_AssignRingSizetoNonCycleMembers(MinimumRingSize, MinRingSize, NumCycles, this);
|
---|
985 | };
|
---|
986 |
|
---|
987 | /** Sets the next component number.
|
---|
988 | * This is O(N) as the number of bonds per atom is bound.
|
---|
989 | * \param *vertex atom whose next atom::*ComponentNr is to be set
|
---|
990 | * \param nr number to use
|
---|
991 | */
|
---|
992 | void molecule::SetNextComponentNumber(atom *vertex, int nr) const
|
---|
993 | {
|
---|
994 | size_t i = 0;
|
---|
995 | if (vertex != NULL) {
|
---|
996 | for (; i < vertex->ListOfBonds.size(); i++) {
|
---|
997 | if (vertex->ComponentNr[i] == -1) { // check if not yet used
|
---|
998 | vertex->ComponentNr[i] = nr;
|
---|
999 | break;
|
---|
1000 | } else if (vertex->ComponentNr[i] == nr) // if number is already present, don't add another time
|
---|
1001 | break; // breaking here will not cause error!
|
---|
1002 | }
|
---|
1003 | if (i == vertex->ListOfBonds.size()) {
|
---|
1004 | DoeLog(0) && (eLog()<< Verbose(0) << "Error: All Component entries are already occupied!" << endl);
|
---|
1005 | performCriticalExit();
|
---|
1006 | }
|
---|
1007 | } else {
|
---|
1008 | DoeLog(0) && (eLog()<< Verbose(0) << "Error: Given vertex is NULL!" << endl);
|
---|
1009 | performCriticalExit();
|
---|
1010 | }
|
---|
1011 | }
|
---|
1012 | ;
|
---|
1013 |
|
---|
1014 | /** Returns next unused bond for this atom \a *vertex or NULL of none exists.
|
---|
1015 | * \param *vertex atom to regard
|
---|
1016 | * \return bond class or NULL
|
---|
1017 | */
|
---|
1018 | bond * molecule::FindNextUnused(atom *vertex) const
|
---|
1019 | {
|
---|
1020 | for (BondList::const_iterator Runner = vertex->ListOfBonds.begin(); Runner != vertex->ListOfBonds.end(); (++Runner))
|
---|
1021 | if ((*Runner)->IsUsed() == white)
|
---|
1022 | return ((*Runner));
|
---|
1023 | return NULL;
|
---|
1024 | }
|
---|
1025 | ;
|
---|
1026 |
|
---|
1027 | /** Resets bond::Used flag of all bonds in this molecule.
|
---|
1028 | * \return true - success, false - -failure
|
---|
1029 | */
|
---|
1030 | void molecule::ResetAllBondsToUnused() const
|
---|
1031 | {
|
---|
1032 | for(molecule::const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner)
|
---|
1033 | for(BondList::const_iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); BondRunner != (*AtomRunner)->ListOfBonds.end(); ++BondRunner)
|
---|
1034 | if ((*BondRunner)->leftatom == *AtomRunner)
|
---|
1035 | (*BondRunner)->ResetUsed();
|
---|
1036 | }
|
---|
1037 | ;
|
---|
1038 |
|
---|
1039 | /** Output a list of flags, stating whether the bond was visited or not.
|
---|
1040 | * \param *out output stream for debugging
|
---|
1041 | * \param *list
|
---|
1042 | */
|
---|
1043 | void OutputAlreadyVisited(int *list)
|
---|
1044 | {
|
---|
1045 | DoLog(4) && (Log() << Verbose(4) << "Already Visited Bonds:\t");
|
---|
1046 | for (int i = 1; i <= list[0]; i++)
|
---|
1047 | DoLog(0) && (Log() << Verbose(0) << list[i] << " ");
|
---|
1048 | DoLog(0) && (Log() << Verbose(0) << endl);
|
---|
1049 | }
|
---|
1050 | ;
|
---|
1051 |
|
---|
1052 | /** Storing the bond structure of a molecule to file.
|
---|
1053 | * Simply stores Atom::nr and then the Atom::nr of all bond partners per line.
|
---|
1054 | * \param &filename name of file
|
---|
1055 | * \param path path to file, defaults to empty
|
---|
1056 | * \return true - file written successfully, false - writing failed
|
---|
1057 | */
|
---|
1058 | bool molecule::StoreAdjacencyToFile(std::string filename, std::string path)
|
---|
1059 | {
|
---|
1060 | ofstream AdjacencyFile;
|
---|
1061 | string line;
|
---|
1062 | bool status = true;
|
---|
1063 |
|
---|
1064 | if (path != "")
|
---|
1065 | line = path + "/" + filename;
|
---|
1066 | else
|
---|
1067 | line = filename;
|
---|
1068 | AdjacencyFile.open(line.c_str(), ios::out);
|
---|
1069 | DoLog(1) && (Log() << Verbose(1) << "Saving adjacency list ... " << endl);
|
---|
1070 | if (AdjacencyFile.good()) {
|
---|
1071 | AdjacencyFile << "m\tn" << endl;
|
---|
1072 | for_each(atoms.begin(),atoms.end(),bind2nd(mem_fun(&atom::OutputAdjacency),&AdjacencyFile));
|
---|
1073 | AdjacencyFile.close();
|
---|
1074 | DoLog(1) && (Log() << Verbose(1) << "\t... done." << endl);
|
---|
1075 | } else {
|
---|
1076 | DoLog(1) && (Log() << Verbose(1) << "\t... failed to open file " << line << "." << endl);
|
---|
1077 | status = false;
|
---|
1078 | }
|
---|
1079 |
|
---|
1080 | return status;
|
---|
1081 | }
|
---|
1082 | ;
|
---|
1083 |
|
---|
1084 | /** Storing the bond structure of a molecule to file.
|
---|
1085 | * Simply stores Atom::nr and then the Atom::nr of all bond partners, one per line.
|
---|
1086 | * \param &filename name of file
|
---|
1087 | * \param path path to file, defaults to empty
|
---|
1088 | * \return true - file written successfully, false - writing failed
|
---|
1089 | */
|
---|
1090 | bool molecule::StoreBondsToFile(std::string filename, std::string path)
|
---|
1091 | {
|
---|
1092 | ofstream BondFile;
|
---|
1093 | string line;
|
---|
1094 | bool status = true;
|
---|
1095 |
|
---|
1096 | if (path != "")
|
---|
1097 | line = path + "/" + filename;
|
---|
1098 | else
|
---|
1099 | line = filename;
|
---|
1100 | BondFile.open(line.c_str(), ios::out);
|
---|
1101 | DoLog(1) && (Log() << Verbose(1) << "Saving adjacency list ... " << endl);
|
---|
1102 | if (BondFile.good()) {
|
---|
1103 | BondFile << "m\tn" << endl;
|
---|
1104 | for_each(atoms.begin(),atoms.end(),bind2nd(mem_fun(&atom::OutputBonds),&BondFile));
|
---|
1105 | BondFile.close();
|
---|
1106 | DoLog(1) && (Log() << Verbose(1) << "\t... done." << endl);
|
---|
1107 | } else {
|
---|
1108 | DoLog(1) && (Log() << Verbose(1) << "\t... failed to open file " << line << "." << endl);
|
---|
1109 | status = false;
|
---|
1110 | }
|
---|
1111 |
|
---|
1112 | return status;
|
---|
1113 | }
|
---|
1114 | ;
|
---|
1115 |
|
---|
1116 | bool CheckAdjacencyFileAgainstMolecule_Init(std::string &path, ifstream &File, int *&CurrentBonds)
|
---|
1117 | {
|
---|
1118 | string filename;
|
---|
1119 | filename = path + ADJACENCYFILE;
|
---|
1120 | File.open(filename.c_str(), ios::out);
|
---|
1121 | DoLog(1) && (Log() << Verbose(1) << "Looking at bond structure stored in adjacency file and comparing to present one ... " << endl);
|
---|
1122 | if (File.fail())
|
---|
1123 | return false;
|
---|
1124 |
|
---|
1125 | // allocate storage structure
|
---|
1126 | CurrentBonds = new int[8]; // contains parsed bonds of current atom
|
---|
1127 | for(int i=0;i<8;i++)
|
---|
1128 | CurrentBonds[i] = 0;
|
---|
1129 | return true;
|
---|
1130 | }
|
---|
1131 | ;
|
---|
1132 |
|
---|
1133 | void CheckAdjacencyFileAgainstMolecule_Finalize(ifstream &File, int *&CurrentBonds)
|
---|
1134 | {
|
---|
1135 | File.close();
|
---|
1136 | File.clear();
|
---|
1137 | delete[](CurrentBonds);
|
---|
1138 | }
|
---|
1139 | ;
|
---|
1140 |
|
---|
1141 | void CheckAdjacencyFileAgainstMolecule_CompareBonds(bool &status, int &NonMatchNumber, atom *&Walker, size_t &CurrentBondsOfAtom, int AtomNr, int *&CurrentBonds, atom **ListOfAtoms)
|
---|
1142 | {
|
---|
1143 | size_t j = 0;
|
---|
1144 | int id = -1;
|
---|
1145 |
|
---|
1146 | //Log() << Verbose(2) << "Walker is " << *Walker << ", bond partners: ";
|
---|
1147 | if (CurrentBondsOfAtom == Walker->ListOfBonds.size()) {
|
---|
1148 | for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
|
---|
1149 | id = (*Runner)->GetOtherAtom(Walker)->nr;
|
---|
1150 | j = 0;
|
---|
1151 | for (; (j < CurrentBondsOfAtom) && (CurrentBonds[j++] != id);)
|
---|
1152 | ; // check against all parsed bonds
|
---|
1153 | if (CurrentBonds[j - 1] != id) { // no match ? Then mark in ListOfAtoms
|
---|
1154 | ListOfAtoms[AtomNr] = NULL;
|
---|
1155 | NonMatchNumber++;
|
---|
1156 | status = false;
|
---|
1157 | DoeLog(2) && (eLog() << Verbose(2) << id << " can not be found in list." << endl);
|
---|
1158 | } else {
|
---|
1159 | //Log() << Verbose(0) << "[" << id << "]\t";
|
---|
1160 | }
|
---|
1161 | }
|
---|
1162 | //Log() << Verbose(0) << endl;
|
---|
1163 | } else {
|
---|
1164 | DoLog(0) && (Log() << Verbose(0) << "Number of bonds for Atom " << *Walker << " does not match, parsed " << CurrentBondsOfAtom << " against " << Walker->ListOfBonds.size() << "." << endl);
|
---|
1165 | status = false;
|
---|
1166 | }
|
---|
1167 | }
|
---|
1168 | ;
|
---|
1169 |
|
---|
1170 | /** Checks contents of adjacency file against bond structure in structure molecule.
|
---|
1171 | * \param *out output stream for debugging
|
---|
1172 | * \param *path path to file
|
---|
1173 | * \param **ListOfAtoms allocated (molecule::AtomCount) and filled lookup table for ids (Atom::nr) to *Atom
|
---|
1174 | * \return true - structure is equal, false - not equivalence
|
---|
1175 | */
|
---|
1176 | bool molecule::CheckAdjacencyFileAgainstMolecule(std::string &path, atom **ListOfAtoms)
|
---|
1177 | {
|
---|
1178 | ifstream File;
|
---|
1179 | bool status = true;
|
---|
1180 | atom *Walker = NULL;
|
---|
1181 | int *CurrentBonds = NULL;
|
---|
1182 | int NonMatchNumber = 0; // will number of atoms with differing bond structure
|
---|
1183 | size_t CurrentBondsOfAtom = -1;
|
---|
1184 | const int AtomCount = getAtomCount();
|
---|
1185 |
|
---|
1186 | if (!CheckAdjacencyFileAgainstMolecule_Init(path, File, CurrentBonds)) {
|
---|
1187 | DoLog(1) && (Log() << Verbose(1) << "Adjacency file not found." << endl);
|
---|
1188 | return true;
|
---|
1189 | }
|
---|
1190 |
|
---|
1191 | char buffer[MAXSTRINGSIZE];
|
---|
1192 | // Parse the file line by line and count the bonds
|
---|
1193 | while (!File.eof()) {
|
---|
1194 | File.getline(buffer, MAXSTRINGSIZE);
|
---|
1195 | stringstream line;
|
---|
1196 | line.str(buffer);
|
---|
1197 | int AtomNr = -1;
|
---|
1198 | line >> AtomNr;
|
---|
1199 | CurrentBondsOfAtom = -1; // we count one too far due to line end
|
---|
1200 | // parse into structure
|
---|
1201 | if ((AtomNr >= 0) && (AtomNr < AtomCount)) {
|
---|
1202 | Walker = ListOfAtoms[AtomNr];
|
---|
1203 | while (!line.eof())
|
---|
1204 | line >> CurrentBonds[++CurrentBondsOfAtom];
|
---|
1205 | // compare against present bonds
|
---|
1206 | CheckAdjacencyFileAgainstMolecule_CompareBonds(status, NonMatchNumber, Walker, CurrentBondsOfAtom, AtomNr, CurrentBonds, ListOfAtoms);
|
---|
1207 | } else {
|
---|
1208 | if (AtomNr != -1)
|
---|
1209 | DoeLog(2) && (eLog() << Verbose(2) << AtomNr << " is not valid in the range of ids [" << 0 << "," << AtomCount << ")." << endl);
|
---|
1210 | }
|
---|
1211 | }
|
---|
1212 | CheckAdjacencyFileAgainstMolecule_Finalize(File, CurrentBonds);
|
---|
1213 |
|
---|
1214 | if (status) { // if equal we parse the KeySetFile
|
---|
1215 | DoLog(1) && (Log() << Verbose(1) << "done: Equal." << endl);
|
---|
1216 | } else
|
---|
1217 | DoLog(1) && (Log() << Verbose(1) << "done: Not equal by " << NonMatchNumber << " atoms." << endl);
|
---|
1218 | return status;
|
---|
1219 | }
|
---|
1220 | ;
|
---|
1221 |
|
---|
1222 | /** Picks from a global stack with all back edges the ones in the fragment.
|
---|
1223 | * \param *out output stream for debugging
|
---|
1224 | * \param **ListOfLocalAtoms array of father atom::nr to local atom::nr (reverse of atom::father)
|
---|
1225 | * \param *ReferenceStack stack with all the back egdes
|
---|
1226 | * \param *LocalStack stack to be filled
|
---|
1227 | * \return true - everything ok, false - ReferenceStack was empty
|
---|
1228 | */
|
---|
1229 | bool molecule::PickLocalBackEdges(atom **ListOfLocalAtoms, class StackClass<bond *> *&ReferenceStack, class StackClass<bond *> *&LocalStack) const
|
---|
1230 | {
|
---|
1231 | bool status = true;
|
---|
1232 | if (ReferenceStack->IsEmpty()) {
|
---|
1233 | DoLog(1) && (Log() << Verbose(1) << "ReferenceStack is empty!" << endl);
|
---|
1234 | return false;
|
---|
1235 | }
|
---|
1236 | bond *Binder = ReferenceStack->PopFirst();
|
---|
1237 | bond *FirstBond = Binder; // mark the first bond, so that we don't loop through the stack indefinitely
|
---|
1238 | atom *Walker = NULL, *OtherAtom = NULL;
|
---|
1239 | ReferenceStack->Push(Binder);
|
---|
1240 |
|
---|
1241 | do { // go through all bonds and push local ones
|
---|
1242 | Walker = ListOfLocalAtoms[Binder->leftatom->nr]; // get one atom in the reference molecule
|
---|
1243 | if (Walker != NULL) // if this Walker exists in the subgraph ...
|
---|
1244 | for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
|
---|
1245 | OtherAtom = (*Runner)->GetOtherAtom(Walker);
|
---|
1246 | if (OtherAtom == ListOfLocalAtoms[(*Runner)->rightatom->nr]) { // found the bond
|
---|
1247 | LocalStack->Push((*Runner));
|
---|
1248 | DoLog(3) && (Log() << Verbose(3) << "Found local edge " << *(*Runner) << "." << endl);
|
---|
1249 | break;
|
---|
1250 | }
|
---|
1251 | }
|
---|
1252 | Binder = ReferenceStack->PopFirst(); // loop the stack for next item
|
---|
1253 | DoLog(3) && (Log() << Verbose(3) << "Current candidate edge " << Binder << "." << endl);
|
---|
1254 | ReferenceStack->Push(Binder);
|
---|
1255 | } while (FirstBond != Binder);
|
---|
1256 |
|
---|
1257 | return status;
|
---|
1258 | }
|
---|
1259 | ;
|
---|
1260 |
|
---|
1261 | void BreadthFirstSearchAdd_Init(struct BFSAccounting &BFS, atom *&Root, int AtomCount, int BondOrder, atom **AddedAtomList = NULL)
|
---|
1262 | {
|
---|
1263 | BFS.AtomCount = AtomCount;
|
---|
1264 | BFS.BondOrder = BondOrder;
|
---|
1265 | BFS.PredecessorList = new atom*[AtomCount];
|
---|
1266 | BFS.ShortestPathList = new int[AtomCount];
|
---|
1267 | BFS.ColorList = new enum Shading[AtomCount];
|
---|
1268 | BFS.BFSStack = new StackClass<atom *> (AtomCount);
|
---|
1269 |
|
---|
1270 | BFS.Root = Root;
|
---|
1271 | BFS.BFSStack->ClearStack();
|
---|
1272 | BFS.BFSStack->Push(Root);
|
---|
1273 |
|
---|
1274 | // initialise each vertex as white with no predecessor, empty queue, color Root lightgray
|
---|
1275 | for (int i = AtomCount; i--;) {
|
---|
1276 | BFS.PredecessorList[i] = NULL;
|
---|
1277 | BFS.ShortestPathList[i] = -1;
|
---|
1278 | if ((AddedAtomList != NULL) && (AddedAtomList[i] != NULL)) // mark already present atoms (i.e. Root and maybe others) as visited
|
---|
1279 | BFS.ColorList[i] = lightgray;
|
---|
1280 | else
|
---|
1281 | BFS.ColorList[i] = white;
|
---|
1282 | }
|
---|
1283 | //BFS.ShortestPathList[Root->nr] = 0; // done by Calloc
|
---|
1284 | }
|
---|
1285 | ;
|
---|
1286 |
|
---|
1287 | void BreadthFirstSearchAdd_Free(struct BFSAccounting &BFS)
|
---|
1288 | {
|
---|
1289 | delete[](BFS.PredecessorList);
|
---|
1290 | delete[](BFS.ShortestPathList);
|
---|
1291 | delete[](BFS.ColorList);
|
---|
1292 | delete (BFS.BFSStack);
|
---|
1293 | BFS.AtomCount = 0;
|
---|
1294 | }
|
---|
1295 | ;
|
---|
1296 |
|
---|
1297 | void BreadthFirstSearchAdd_UnvisitedNode(molecule *Mol, struct BFSAccounting &BFS, atom *&Walker, atom *&OtherAtom, bond *&Binder, bond *&Bond, atom **&AddedAtomList, bond **&AddedBondList, bool IsAngstroem)
|
---|
1298 | {
|
---|
1299 | if (Binder != Bond) // let other atom white if it's via Root bond. In case it's cyclic it has to be reached again (yet Root is from OtherAtom already black, thus no problem)
|
---|
1300 | BFS.ColorList[OtherAtom->nr] = lightgray;
|
---|
1301 | BFS.PredecessorList[OtherAtom->nr] = Walker; // Walker is the predecessor
|
---|
1302 | BFS.ShortestPathList[OtherAtom->nr] = BFS.ShortestPathList[Walker->nr] + 1;
|
---|
1303 | DoLog(2) && (Log() << Verbose(2) << "Coloring OtherAtom " << OtherAtom->getName() << " " << ((BFS.ColorList[OtherAtom->nr] == white) ? "white" : "lightgray") << ", its predecessor is " << Walker->getName() << " and its Shortest Path is " << BFS.ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl);
|
---|
1304 | if ((((BFS.ShortestPathList[OtherAtom->nr] < BFS.BondOrder) && (Binder != Bond)))) { // Check for maximum distance
|
---|
1305 | DoLog(3) && (Log() << Verbose(3));
|
---|
1306 | if (AddedAtomList[OtherAtom->nr] == NULL) { // add if it's not been so far
|
---|
1307 | AddedAtomList[OtherAtom->nr] = Mol->AddCopyAtom(OtherAtom);
|
---|
1308 | DoLog(0) && (Log() << Verbose(0) << "Added OtherAtom " << OtherAtom->getName());
|
---|
1309 | AddedBondList[Binder->nr] = Mol->CopyBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder);
|
---|
1310 | DoLog(0) && (Log() << Verbose(0) << " and bond " << *(AddedBondList[Binder->nr]) << ", ");
|
---|
1311 | } else { // this code should actually never come into play (all white atoms are not yet present in BondMolecule, that's why they are white in the first place)
|
---|
1312 | DoLog(0) && (Log() << Verbose(0) << "Not adding OtherAtom " << OtherAtom->getName());
|
---|
1313 | if (AddedBondList[Binder->nr] == NULL) {
|
---|
1314 | AddedBondList[Binder->nr] = Mol->CopyBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder);
|
---|
1315 | DoLog(0) && (Log() << Verbose(0) << ", added Bond " << *(AddedBondList[Binder->nr]));
|
---|
1316 | } else
|
---|
1317 | DoLog(0) && (Log() << Verbose(0) << ", not added Bond ");
|
---|
1318 | }
|
---|
1319 | DoLog(0) && (Log() << Verbose(0) << ", putting OtherAtom into queue." << endl);
|
---|
1320 | BFS.BFSStack->Push(OtherAtom);
|
---|
1321 | } else { // out of bond order, then replace
|
---|
1322 | if ((AddedAtomList[OtherAtom->nr] == NULL) && (Binder->Cyclic))
|
---|
1323 | BFS.ColorList[OtherAtom->nr] = white; // unmark if it has not been queued/added, to make it available via its other bonds (cyclic)
|
---|
1324 | if (Binder == Bond)
|
---|
1325 | DoLog(3) && (Log() << Verbose(3) << "Not Queueing, is the Root bond");
|
---|
1326 | else if (BFS.ShortestPathList[OtherAtom->nr] >= BFS.BondOrder)
|
---|
1327 | DoLog(3) && (Log() << Verbose(3) << "Not Queueing, is out of Bond Count of " << BFS.BondOrder);
|
---|
1328 | if (!Binder->Cyclic)
|
---|
1329 | DoLog(0) && (Log() << Verbose(0) << ", is not part of a cyclic bond, saturating bond with Hydrogen." << endl);
|
---|
1330 | if (AddedBondList[Binder->nr] == NULL) {
|
---|
1331 | if ((AddedAtomList[OtherAtom->nr] != NULL)) { // .. whether we add or saturate
|
---|
1332 | AddedBondList[Binder->nr] = Mol->CopyBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder);
|
---|
1333 | } else {
|
---|
1334 | #ifdef ADDHYDROGEN
|
---|
1335 | if (!Mol->AddHydrogenReplacementAtom(Binder, AddedAtomList[Walker->nr], Walker, OtherAtom, IsAngstroem))
|
---|
1336 | exit(1);
|
---|
1337 | #endif
|
---|
1338 | }
|
---|
1339 | }
|
---|
1340 | }
|
---|
1341 | }
|
---|
1342 | ;
|
---|
1343 |
|
---|
1344 | void BreadthFirstSearchAdd_VisitedNode(molecule *Mol, struct BFSAccounting &BFS, atom *&Walker, atom *&OtherAtom, bond *&Binder, bond *&Bond, atom **&AddedAtomList, bond **&AddedBondList, bool IsAngstroem)
|
---|
1345 | {
|
---|
1346 | DoLog(3) && (Log() << Verbose(3) << "Not Adding, has already been visited." << endl);
|
---|
1347 | // This has to be a cyclic bond, check whether it's present ...
|
---|
1348 | if (AddedBondList[Binder->nr] == NULL) {
|
---|
1349 | if ((Binder != Bond) && (Binder->Cyclic) && (((BFS.ShortestPathList[Walker->nr] + 1) < BFS.BondOrder))) {
|
---|
1350 | AddedBondList[Binder->nr] = Mol->CopyBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder);
|
---|
1351 | } else { // if it's root bond it has to broken (otherwise we would not create the fragments)
|
---|
1352 | #ifdef ADDHYDROGEN
|
---|
1353 | if(!Mol->AddHydrogenReplacementAtom(Binder, AddedAtomList[Walker->nr], Walker, OtherAtom, IsAngstroem))
|
---|
1354 | exit(1);
|
---|
1355 | #endif
|
---|
1356 | }
|
---|
1357 | }
|
---|
1358 | }
|
---|
1359 | ;
|
---|
1360 |
|
---|
1361 | /** Adds atoms up to \a BondCount distance from \a *Root and notes them down in \a **AddedAtomList.
|
---|
1362 | * Gray vertices are always enqueued in an StackClass<atom *> FIFO queue, the rest is usual BFS with adding vertices found was
|
---|
1363 | * white and putting into queue.
|
---|
1364 | * \param *out output stream for debugging
|
---|
1365 | * \param *Mol Molecule class to add atoms to
|
---|
1366 | * \param **AddedAtomList list with added atom pointers, index is atom father's number
|
---|
1367 | * \param **AddedBondList list with added bond pointers, index is bond father's number
|
---|
1368 | * \param *Root root vertex for BFS
|
---|
1369 | * \param *Bond bond not to look beyond
|
---|
1370 | * \param BondOrder maximum distance for vertices to add
|
---|
1371 | * \param IsAngstroem lengths are in angstroem or bohrradii
|
---|
1372 | */
|
---|
1373 | void molecule::BreadthFirstSearchAdd(molecule *Mol, atom **&AddedAtomList, bond **&AddedBondList, atom *Root, bond *Bond, int BondOrder, bool IsAngstroem)
|
---|
1374 | {
|
---|
1375 | struct BFSAccounting BFS;
|
---|
1376 | atom *Walker = NULL, *OtherAtom = NULL;
|
---|
1377 | bond *Binder = NULL;
|
---|
1378 |
|
---|
1379 | // add Root if not done yet
|
---|
1380 | if (AddedAtomList[Root->nr] == NULL) // add Root if not yet present
|
---|
1381 | AddedAtomList[Root->nr] = Mol->AddCopyAtom(Root);
|
---|
1382 |
|
---|
1383 | BreadthFirstSearchAdd_Init(BFS, Root, BondOrder, getAtomCount(), AddedAtomList);
|
---|
1384 |
|
---|
1385 | // and go on ... Queue always contains all lightgray vertices
|
---|
1386 | while (!BFS.BFSStack->IsEmpty()) {
|
---|
1387 | // we have to pop the oldest atom from stack. This keeps the atoms on the stack always of the same ShortestPath distance.
|
---|
1388 | // e.g. if current atom is 2, push to end of stack are of length 3, but first all of length 2 would be popped. They again
|
---|
1389 | // append length of 3 (their neighbours). Thus on stack we have always atoms of a certain length n at bottom of stack and
|
---|
1390 | // followed by n+1 till top of stack.
|
---|
1391 | Walker = BFS.BFSStack->PopFirst(); // pop oldest added
|
---|
1392 | DoLog(1) && (Log() << Verbose(1) << "Current Walker is: " << Walker->getName() << ", and has " << Walker->ListOfBonds.size() << " bonds." << endl);
|
---|
1393 | for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
|
---|
1394 | if ((*Runner) != NULL) { // don't look at bond equal NULL
|
---|
1395 | Binder = (*Runner);
|
---|
1396 | OtherAtom = (*Runner)->GetOtherAtom(Walker);
|
---|
1397 | DoLog(2) && (Log() << Verbose(2) << "Current OtherAtom is: " << OtherAtom->getName() << " for bond " << *(*Runner) << "." << endl);
|
---|
1398 | if (BFS.ColorList[OtherAtom->nr] == white) {
|
---|
1399 | BreadthFirstSearchAdd_UnvisitedNode(Mol, BFS, Walker, OtherAtom, Binder, Bond, AddedAtomList, AddedBondList, IsAngstroem);
|
---|
1400 | } else {
|
---|
1401 | BreadthFirstSearchAdd_VisitedNode(Mol, BFS, Walker, OtherAtom, Binder, Bond, AddedAtomList, AddedBondList, IsAngstroem);
|
---|
1402 | }
|
---|
1403 | }
|
---|
1404 | }
|
---|
1405 | BFS.ColorList[Walker->nr] = black;
|
---|
1406 | DoLog(1) && (Log() << Verbose(1) << "Coloring Walker " << Walker->getName() << " black." << endl);
|
---|
1407 | }
|
---|
1408 | BreadthFirstSearchAdd_Free(BFS);
|
---|
1409 | }
|
---|
1410 | ;
|
---|
1411 |
|
---|
1412 | /** Adds a bond as a copy to a given one
|
---|
1413 | * \param *left leftatom of new bond
|
---|
1414 | * \param *right rightatom of new bond
|
---|
1415 | * \param *CopyBond rest of fields in bond are copied from this
|
---|
1416 | * \return pointer to new bond
|
---|
1417 | */
|
---|
1418 | bond * molecule::CopyBond(atom *left, atom *right, bond *CopyBond)
|
---|
1419 | {
|
---|
1420 | bond *Binder = AddBond(left, right, CopyBond->BondDegree);
|
---|
1421 | Binder->Cyclic = CopyBond->Cyclic;
|
---|
1422 | Binder->Type = CopyBond->Type;
|
---|
1423 | return Binder;
|
---|
1424 | }
|
---|
1425 | ;
|
---|
1426 |
|
---|
1427 | void BuildInducedSubgraph_Init(atom **&ParentList, int AtomCount)
|
---|
1428 | {
|
---|
1429 | // reset parent list
|
---|
1430 | ParentList = new atom*[AtomCount];
|
---|
1431 | for (int i=0;i<AtomCount;i++)
|
---|
1432 | ParentList[i] = NULL;
|
---|
1433 | DoLog(3) && (Log() << Verbose(3) << "Resetting ParentList." << endl);
|
---|
1434 | }
|
---|
1435 | ;
|
---|
1436 |
|
---|
1437 | void BuildInducedSubgraph_FillParentList(const molecule *mol, const molecule *Father, atom **&ParentList)
|
---|
1438 | {
|
---|
1439 | // fill parent list with sons
|
---|
1440 | DoLog(3) && (Log() << Verbose(3) << "Filling Parent List." << endl);
|
---|
1441 | for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
|
---|
1442 | ParentList[(*iter)->father->nr] = (*iter);
|
---|
1443 | // Outputting List for debugging
|
---|
1444 | DoLog(4) && (Log() << Verbose(4) << "Son[" << (*iter)->father->nr << "] of " << (*iter)->father << " is " << ParentList[(*iter)->father->nr] << "." << endl);
|
---|
1445 | }
|
---|
1446 | };
|
---|
1447 |
|
---|
1448 | void BuildInducedSubgraph_Finalize(atom **&ParentList)
|
---|
1449 | {
|
---|
1450 | delete[](ParentList);
|
---|
1451 | }
|
---|
1452 | ;
|
---|
1453 |
|
---|
1454 | bool BuildInducedSubgraph_CreateBondsFromParent(molecule *mol, const molecule *Father, atom **&ParentList)
|
---|
1455 | {
|
---|
1456 | bool status = true;
|
---|
1457 | atom *OtherAtom = NULL;
|
---|
1458 | // check each entry of parent list and if ok (one-to-and-onto matching) create bonds
|
---|
1459 | DoLog(3) && (Log() << Verbose(3) << "Creating bonds." << endl);
|
---|
1460 | for (molecule::const_iterator iter = Father->begin(); iter != Father->end(); ++iter) {
|
---|
1461 | if (ParentList[(*iter)->nr] != NULL) {
|
---|
1462 | if (ParentList[(*iter)->nr]->father != (*iter)) {
|
---|
1463 | status = false;
|
---|
1464 | } else {
|
---|
1465 | for (BondList::const_iterator Runner = (*iter)->ListOfBonds.begin(); Runner != (*iter)->ListOfBonds.end(); (++Runner)) {
|
---|
1466 | OtherAtom = (*Runner)->GetOtherAtom((*iter));
|
---|
1467 | if (ParentList[OtherAtom->nr] != NULL) { // if otheratom is also a father of an atom on this molecule, create the bond
|
---|
1468 | DoLog(4) && (Log() << Verbose(4) << "Endpoints of Bond " << (*Runner) << " are both present: " << ParentList[(*iter)->nr]->getName() << " and " << ParentList[OtherAtom->nr]->getName() << "." << endl);
|
---|
1469 | mol->AddBond(ParentList[(*iter)->nr], ParentList[OtherAtom->nr], (*Runner)->BondDegree);
|
---|
1470 | }
|
---|
1471 | }
|
---|
1472 | }
|
---|
1473 | }
|
---|
1474 | }
|
---|
1475 | return status;
|
---|
1476 | }
|
---|
1477 | ;
|
---|
1478 |
|
---|
1479 | /** Adds bond structure to this molecule from \a Father molecule.
|
---|
1480 | * This basically causes this molecule to become an induced subgraph of the \a Father, i.e. for every bond in Father
|
---|
1481 | * with end points present in this molecule, bond is created in this molecule.
|
---|
1482 | * Special care was taken to ensure that this is of complexity O(N), where N is the \a Father's molecule::AtomCount.
|
---|
1483 | * \param *out output stream for debugging
|
---|
1484 | * \param *Father father molecule
|
---|
1485 | * \return true - is induced subgraph, false - there are atoms with fathers not in \a Father
|
---|
1486 | * \todo not checked, not fully working probably
|
---|
1487 | */
|
---|
1488 | bool molecule::BuildInducedSubgraph(const molecule *Father)
|
---|
1489 | {
|
---|
1490 | bool status = true;
|
---|
1491 | atom **ParentList = NULL;
|
---|
1492 | DoLog(2) && (Log() << Verbose(2) << "Begin of BuildInducedSubgraph." << endl);
|
---|
1493 | BuildInducedSubgraph_Init(ParentList, Father->getAtomCount());
|
---|
1494 | BuildInducedSubgraph_FillParentList(this, Father, ParentList);
|
---|
1495 | status = BuildInducedSubgraph_CreateBondsFromParent(this, Father, ParentList);
|
---|
1496 | BuildInducedSubgraph_Finalize(ParentList);
|
---|
1497 | DoLog(2) && (Log() << Verbose(2) << "End of BuildInducedSubgraph." << endl);
|
---|
1498 | return status;
|
---|
1499 | }
|
---|
1500 | ;
|
---|
1501 |
|
---|
1502 | /** For a given keyset \a *Fragment, checks whether it is connected in the current molecule.
|
---|
1503 | * \param *out output stream for debugging
|
---|
1504 | * \param *Fragment Keyset of fragment's vertices
|
---|
1505 | * \return true - connected, false - disconnected
|
---|
1506 | * \note this is O(n^2) for it's just a bug checker not meant for permanent use!
|
---|
1507 | */
|
---|
1508 | bool molecule::CheckForConnectedSubgraph(KeySet *Fragment)
|
---|
1509 | {
|
---|
1510 | atom *Walker = NULL, *Walker2 = NULL;
|
---|
1511 | bool BondStatus = false;
|
---|
1512 | int size;
|
---|
1513 |
|
---|
1514 | DoLog(1) && (Log() << Verbose(1) << "Begin of CheckForConnectedSubgraph" << endl);
|
---|
1515 | DoLog(2) && (Log() << Verbose(2) << "Disconnected atom: ");
|
---|
1516 |
|
---|
1517 | // count number of atoms in graph
|
---|
1518 | size = 0;
|
---|
1519 | for (KeySet::iterator runner = Fragment->begin(); runner != Fragment->end(); runner++)
|
---|
1520 | size++;
|
---|
1521 | if (size > 1)
|
---|
1522 | for (KeySet::iterator runner = Fragment->begin(); runner != Fragment->end(); runner++) {
|
---|
1523 | Walker = FindAtom(*runner);
|
---|
1524 | BondStatus = false;
|
---|
1525 | for (KeySet::iterator runners = Fragment->begin(); runners != Fragment->end(); runners++) {
|
---|
1526 | Walker2 = FindAtom(*runners);
|
---|
1527 | for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
|
---|
1528 | if ((*Runner)->GetOtherAtom(Walker) == Walker2) {
|
---|
1529 | BondStatus = true;
|
---|
1530 | break;
|
---|
1531 | }
|
---|
1532 | if (BondStatus)
|
---|
1533 | break;
|
---|
1534 | }
|
---|
1535 | }
|
---|
1536 | if (!BondStatus) {
|
---|
1537 | DoLog(0) && (Log() << Verbose(0) << (*Walker) << endl);
|
---|
1538 | return false;
|
---|
1539 | }
|
---|
1540 | }
|
---|
1541 | else {
|
---|
1542 | DoLog(0) && (Log() << Verbose(0) << "none." << endl);
|
---|
1543 | return true;
|
---|
1544 | }
|
---|
1545 | DoLog(0) && (Log() << Verbose(0) << "none." << endl);
|
---|
1546 |
|
---|
1547 | DoLog(1) && (Log() << Verbose(1) << "End of CheckForConnectedSubgraph" << endl);
|
---|
1548 |
|
---|
1549 | return true;
|
---|
1550 | }
|
---|