source: src/molecule_graph.cpp@ 5f8660a

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 5f8660a was bcf653, checked in by Frederik Heber <heber@…>, 14 years ago

Added copyright note to each .cpp file and an extensive one to builder.cpp.

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