source: src/molecule_graph.cpp@ ba9f5b

Action_Thermostats Add_AtomRandomPerturbation Add_FitFragmentPartialChargesAction Add_RotateAroundBondAction Add_SelectAtomByNameAction Added_ParseSaveFragmentResults AddingActions_SaveParseParticleParameters Adding_Graph_to_ChangeBondActions Adding_MD_integration_tests Adding_ParticleName_to_Atom Adding_StructOpt_integration_tests AtomFragments Automaking_mpqc_open AutomationFragmentation_failures Candidate_v1.5.4 Candidate_v1.6.0 Candidate_v1.6.1 ChangeBugEmailaddress ChangingTestPorts ChemicalSpaceEvaluator CombiningParticlePotentialParsing Combining_Subpackages Debian_Package_split Debian_package_split_molecuildergui_only Disabling_MemDebug Docu_Python_wait EmpiricalPotential_contain_HomologyGraph EmpiricalPotential_contain_HomologyGraph_documentation Enable_parallel_make_install Enhance_userguide Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Example_ManyWaysToTranslateAtom Exclude_Hydrogens_annealWithBondGraph FitPartialCharges_GlobalError Fix_BoundInBox_CenterInBox_MoleculeActions Fix_ChargeSampling_PBC Fix_ChronosMutex Fix_FitPartialCharges Fix_FitPotential_needs_atomicnumbers Fix_ForceAnnealing Fix_IndependentFragmentGrids Fix_ParseParticles Fix_ParseParticles_split_forward_backward_Actions Fix_PopActions Fix_QtFragmentList_sorted_selection Fix_Restrictedkeyset_FragmentMolecule Fix_StatusMsg Fix_StepWorldTime_single_argument Fix_Verbose_Codepatterns Fix_fitting_potentials Fixes ForceAnnealing_goodresults ForceAnnealing_oldresults ForceAnnealing_tocheck ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_continued ForceAnnealing_with_BondGraph_continued_betteresults ForceAnnealing_with_BondGraph_contraction-expansion FragmentAction_writes_AtomFragments FragmentMolecule_checks_bonddegrees GeometryObjects Gui_Fixes Gui_displays_atomic_force_velocity ImplicitCharges IndependentFragmentGrids IndependentFragmentGrids_IndividualZeroInstances IndependentFragmentGrids_IntegrationTest IndependentFragmentGrids_Sole_NN_Calculation JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool JobMarket_unresolvable_hostname_fix MoreRobust_FragmentAutomation ODR_violation_mpqc_open PartialCharges_OrthogonalSummation PdbParser_setsAtomName PythonUI_with_named_parameters QtGui_reactivate_TimeChanged_changes Recreated_GuiChecks Rewrite_FitPartialCharges RotateToPrincipalAxisSystem_UndoRedo SaturateAtoms_findBestMatching SaturateAtoms_singleDegree StoppableMakroAction Subpackage_CodePatterns Subpackage_JobMarket Subpackage_LinearAlgebra Subpackage_levmar Subpackage_mpqc_open Subpackage_vmg Switchable_LogView ThirdParty_MPQC_rebuilt_buildsystem TrajectoryDependenant_MaxOrder TremoloParser_IncreasedPrecision TremoloParser_MultipleTimesteps TremoloParser_setsAtomName Ubuntu_1604_changes stable
Last change on this file since ba9f5b was c27778, checked in by Frederik Heber <heber@…>, 14 years ago

MEMFIXES: ListOfLocalAtoms in molecule::FragmentMolecule() was not free'd correctly.

NOTE: All of these lists and maps are hard to understand and make the code very confusing. It's really high time for refactoring.

Signed-off-by: Frederik Heber <heber@…>

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