source: src/molecule_graph.cpp@ 76c0d6

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 76c0d6 was 0de7e8, checked in by Frederik Heber <heber@…>, 15 years ago

FIX: Test Fragmentation is at MaxOrder was broken.

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