source: src/molecule_graph.cpp@ eb33c4

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

Moved files bondgraph.?pp -> Graph/BondGraph.?pp.

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