source: src/molecule_fragmentation.cpp@ ba9f5b

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

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

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

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

  • Property mode set to 100644
File size: 79.2 KB
Line 
1/*
2 * molecule_fragmentation.cpp
3 *
4 * Created on: Oct 5, 2009
5 * Author: heber
6 */
7
8#include "Helpers/MemDebug.hpp"
9
10#include <cstring>
11
12#include "World.hpp"
13#include "atom.hpp"
14#include "bond.hpp"
15#include "config.hpp"
16#include "element.hpp"
17#include "helpers.hpp"
18#include "lists.hpp"
19#include "log.hpp"
20#include "memoryallocator.hpp"
21#include "molecule.hpp"
22#include "periodentafel.hpp"
23#include "World.hpp"
24
25/************************************* Functions for class molecule *********************************/
26
27
28/** Estimates by educated guessing (using upper limit) the expected number of fragments.
29 * The upper limit is
30 * \f[
31 * n = N \cdot C^k
32 * \f]
33 * where \f$C=2^c\f$ and c is the maximum bond degree over N number of atoms.
34 * \param *out output stream for debugging
35 * \param order bond order k
36 * \return number n of fragments
37 */
38int molecule::GuesstimateFragmentCount(int order)
39{
40 size_t c = 0;
41 int FragmentCount;
42 // get maximum bond degree
43 for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
44 c = ((*iter)->ListOfBonds.size() > c) ? (*iter)->ListOfBonds.size() : c;
45 }
46 FragmentCount = NoNonHydrogen*(1 << (c*order));
47 DoLog(1) && (Log() << Verbose(1) << "Upper limit for this subgraph is " << FragmentCount << " for " << NoNonHydrogen << " non-H atoms with maximum bond degree of " << c << "." << endl);
48 return FragmentCount;
49};
50
51/** Scans a single line for number and puts them into \a KeySet.
52 * \param *out output stream for debugging
53 * \param *buffer buffer to scan
54 * \param &CurrentSet filled KeySet on return
55 * \return true - at least one valid atom id parsed, false - CurrentSet is empty
56 */
57bool ScanBufferIntoKeySet(char *buffer, KeySet &CurrentSet)
58{
59 stringstream line;
60 int AtomNr;
61 int status = 0;
62
63 line.str(buffer);
64 while (!line.eof()) {
65 line >> AtomNr;
66 if (AtomNr >= 0) {
67 CurrentSet.insert(AtomNr); // insert at end, hence in same order as in file!
68 status++;
69 } // else it's "-1" or else and thus must not be added
70 }
71 DoLog(1) && (Log() << Verbose(1) << "The scanned KeySet is ");
72 for(KeySet::iterator runner = CurrentSet.begin(); runner != CurrentSet.end(); runner++) {
73 DoLog(0) && (Log() << Verbose(0) << (*runner) << "\t");
74 }
75 DoLog(0) && (Log() << Verbose(0) << endl);
76 return (status != 0);
77};
78
79/** Parses the KeySet file and fills \a *FragmentList from the known molecule structure.
80 * Does two-pass scanning:
81 * -# Scans the keyset file and initialises a temporary graph
82 * -# Scans TEFactors file and sets the TEFactor of each key set in the temporary graph accordingly
83 * Finally, the temporary graph is inserted into the given \a FragmentList for return.
84 * \param &path path to file
85 * \param *FragmentList empty, filled on return
86 * \return true - parsing successfully, false - failure on parsing (FragmentList will be NULL)
87 */
88bool ParseKeySetFile(std::string &path, Graph *&FragmentList)
89{
90 bool status = true;
91 ifstream InputFile;
92 stringstream line;
93 GraphTestPair testGraphInsert;
94 int NumberOfFragments = 0;
95 string filename;
96
97 if (FragmentList == NULL) { // check list pointer
98 FragmentList = new Graph;
99 }
100
101 // 1st pass: open file and read
102 DoLog(1) && (Log() << Verbose(1) << "Parsing the KeySet file ... " << endl);
103 filename = path + KEYSETFILE;
104 InputFile.open(filename.c_str());
105 if (InputFile.good()) {
106 // each line represents a new fragment
107 char buffer[MAXSTRINGSIZE];
108 // 1. parse keysets and insert into temp. graph
109 while (!InputFile.eof()) {
110 InputFile.getline(buffer, MAXSTRINGSIZE);
111 KeySet CurrentSet;
112 if ((strlen(buffer) > 0) && (ScanBufferIntoKeySet(buffer, CurrentSet))) { // if at least one valid atom was added, write config
113 testGraphInsert = FragmentList->insert(GraphPair (CurrentSet,pair<int,double>(NumberOfFragments++,1))); // store fragment number and current factor
114 if (!testGraphInsert.second) {
115 DoeLog(0) && (eLog()<< Verbose(0) << "KeySet file must be corrupt as there are two equal key sets therein!" << endl);
116 performCriticalExit();
117 }
118 }
119 }
120 // 2. Free and done
121 InputFile.close();
122 InputFile.clear();
123 DoLog(1) && (Log() << Verbose(1) << "\t ... done." << endl);
124 } else {
125 DoLog(1) && (Log() << Verbose(1) << "\t ... File " << filename << " not found." << endl);
126 status = false;
127 }
128
129 return status;
130};
131
132/** Parses the TE factors file and fills \a *FragmentList from the known molecule structure.
133 * -# Scans TEFactors file and sets the TEFactor of each key set in the temporary graph accordingly
134 * \param *out output stream for debugging
135 * \param *path path to file
136 * \param *FragmentList graph whose nodes's TE factors are set on return
137 * \return true - parsing successfully, false - failure on parsing
138 */
139bool ParseTEFactorsFile(char *path, Graph *FragmentList)
140{
141 bool status = true;
142 ifstream InputFile;
143 stringstream line;
144 GraphTestPair testGraphInsert;
145 int NumberOfFragments = 0;
146 double TEFactor;
147 char filename[MAXSTRINGSIZE];
148
149 if (FragmentList == NULL) { // check list pointer
150 FragmentList = new Graph;
151 }
152
153 // 2nd pass: open TEFactors file and read
154 DoLog(1) && (Log() << Verbose(1) << "Parsing the TEFactors file ... " << endl);
155 sprintf(filename, "%s/%s%s", path, FRAGMENTPREFIX, TEFACTORSFILE);
156 InputFile.open(filename);
157 if (InputFile != NULL) {
158 // 3. add found TEFactors to each keyset
159 NumberOfFragments = 0;
160 for(Graph::iterator runner = FragmentList->begin();runner != FragmentList->end(); runner++) {
161 if (!InputFile.eof()) {
162 InputFile >> TEFactor;
163 (*runner).second.second = TEFactor;
164 DoLog(2) && (Log() << Verbose(2) << "Setting " << ++NumberOfFragments << " fragment's TEFactor to " << (*runner).second.second << "." << endl);
165 } else {
166 status = false;
167 break;
168 }
169 }
170 // 4. Free and done
171 InputFile.close();
172 DoLog(1) && (Log() << Verbose(1) << "done." << endl);
173 } else {
174 DoLog(1) && (Log() << Verbose(1) << "File " << filename << " not found." << endl);
175 status = false;
176 }
177
178 return status;
179};
180
181/** Stores key sets to file.
182 * \param KeySetList Graph with Keysets
183 * \param &path path to file
184 * \return true - file written successfully, false - writing failed
185 */
186bool StoreKeySetFile(Graph &KeySetList, std::string &path)
187{
188 bool status = true;
189 string line = path + KEYSETFILE;
190 ofstream output(line.c_str());
191
192 // open KeySet file
193 DoLog(1) && (Log() << Verbose(1) << "Saving key sets of the total graph ... ");
194 if(output.good()) {
195 for(Graph::iterator runner = KeySetList.begin(); runner != KeySetList.end(); runner++) {
196 for (KeySet::iterator sprinter = (*runner).first.begin();sprinter != (*runner).first.end(); sprinter++) {
197 if (sprinter != (*runner).first.begin())
198 output << "\t";
199 output << *sprinter;
200 }
201 output << endl;
202 }
203 DoLog(0) && (Log() << Verbose(0) << "done." << endl);
204 } else {
205 DoeLog(0) && (eLog()<< Verbose(0) << "Unable to open " << line << " for writing keysets!" << endl);
206 performCriticalExit();
207 status = false;
208 }
209 output.close();
210 output.clear();
211
212 return status;
213};
214
215
216/** Stores TEFactors to file.
217 * \param *out output stream for debugging
218 * \param KeySetList Graph with factors
219 * \param *path path to file
220 * \return true - file written successfully, false - writing failed
221 */
222bool StoreTEFactorsFile(Graph &KeySetList, char *path)
223{
224 ofstream output;
225 bool status = true;
226 string line;
227
228 // open TEFactors file
229 line = path;
230 line.append("/");
231 line += FRAGMENTPREFIX;
232 line += TEFACTORSFILE;
233 output.open(line.c_str(), ios::out);
234 DoLog(1) && (Log() << Verbose(1) << "Saving TEFactors of the total graph ... ");
235 if(output != NULL) {
236 for(Graph::iterator runner = KeySetList.begin(); runner != KeySetList.end(); runner++)
237 output << (*runner).second.second << endl;
238 DoLog(1) && (Log() << Verbose(1) << "done." << endl);
239 } else {
240 DoLog(1) && (Log() << Verbose(1) << "failed to open " << line << "." << endl);
241 status = false;
242 }
243 output.close();
244
245 return status;
246};
247
248/** For a given graph, sorts KeySets into a (index, keyset) map.
249 * \param *GlobalKeySetList list of keysets with global ids (valid in "this" molecule) needed for adaptive increase
250 * \return map from index to keyset
251 */
252map<int,KeySet> * GraphToIndexedKeySet(Graph *GlobalKeySetList)
253{
254 map<int,KeySet> *IndexKeySetList = new map<int,KeySet>;
255 for(Graph::iterator runner = GlobalKeySetList->begin(); runner != GlobalKeySetList->end(); runner++) {
256 IndexKeySetList->insert( pair<int,KeySet>(runner->second.first,runner->first) );
257 }
258 return IndexKeySetList;
259};
260
261/** Inserts a (\a No, \a value) pair into the list, overwriting present one.
262 * Note if values are equal, No will decided on which is first
263 * \param *out output stream for debugging
264 * \param &AdaptiveCriteriaList list to insert into
265 * \param &IndexedKeySetList list to find key set for a given index \a No
266 * \param FragOrder current bond order of fragment
267 * \param No index of keyset
268 * \param value energy value
269 */
270void InsertIntoAdaptiveCriteriaList(map<int, pair<double,int> > *AdaptiveCriteriaList, map<int,KeySet> &IndexKeySetList, int FragOrder, int No, double Value)
271{
272 map<int,KeySet>::iterator marker = IndexKeySetList.find(No); // find keyset to Frag No.
273 if (marker != IndexKeySetList.end()) { // if found
274 Value *= 1 + MYEPSILON*(*((*marker).second.begin())); // in case of equal energies this makes them not equal without changing anything actually
275 // as the smallest number in each set has always been the root (we use global id to keep the doubles away), seek smallest and insert into AtomMask
276 pair <map<int, pair<double,int> >::iterator, bool> InsertedElement = AdaptiveCriteriaList->insert( make_pair(*((*marker).second.begin()), pair<double,int>( fabs(Value), FragOrder) ));
277 map<int, pair<double,int> >::iterator PresentItem = InsertedElement.first;
278 if (!InsertedElement.second) { // this root is already present
279 if ((*PresentItem).second.second < FragOrder) // if order there is lower, update entry with higher-order term
280 //if ((*PresentItem).second.first < (*runner).first) // as higher-order terms are not always better, we skip this part (which would always include this site into adaptive increase)
281 { // if value is smaller, update value and order
282 (*PresentItem).second.first = fabs(Value);
283 (*PresentItem).second.second = FragOrder;
284 DoLog(2) && (Log() << Verbose(2) << "Updated element (" << (*PresentItem).first << ",[" << (*PresentItem).second.first << "," << (*PresentItem).second.second << "])." << endl);
285 } else {
286 DoLog(2) && (Log() << Verbose(2) << "Did not update element " << (*PresentItem).first << " as " << FragOrder << " is less than or equal to " << (*PresentItem).second.second << "." << endl);
287 }
288 } else {
289 DoLog(2) && (Log() << Verbose(2) << "Inserted element (" << (*PresentItem).first << ",[" << (*PresentItem).second.first << "," << (*PresentItem).second.second << "])." << endl);
290 }
291 } else {
292 DoLog(1) && (Log() << Verbose(1) << "No Fragment under No. " << No << "found." << endl);
293 }
294};
295
296/** Scans the adaptive order file and insert (index, value) into map.
297 * \param &path path to ENERGYPERFRAGMENT file (may be NULL if Order is non-negative)
298 * \param &IndexedKeySetList list to find key set for a given index \a No
299 * \return adaptive criteria list from file
300 */
301map<int, pair<double,int> > * ScanAdaptiveFileIntoMap(std::string &path, map<int,KeySet> &IndexKeySetList)
302{
303 map<int, pair<double,int> > *AdaptiveCriteriaList = new map<int, pair<double,int> >;
304 int No = 0, FragOrder = 0;
305 double Value = 0.;
306 char buffer[MAXSTRINGSIZE];
307 string filename = path + ENERGYPERFRAGMENT;
308 ifstream InputFile(filename.c_str());
309
310 if (InputFile.fail()) {
311 DoeLog(1) && (eLog() << Verbose(1) << "Cannot find file " << filename << "." << endl);
312 return AdaptiveCriteriaList;
313 }
314
315 if (CountLinesinFile(InputFile) > 0) {
316 // each line represents a fragment root (Atom::nr) id and its energy contribution
317 InputFile.getline(buffer, MAXSTRINGSIZE); // skip comment lines
318 InputFile.getline(buffer, MAXSTRINGSIZE);
319 while(!InputFile.eof()) {
320 InputFile.getline(buffer, MAXSTRINGSIZE);
321 if (strlen(buffer) > 2) {
322 //Log() << Verbose(2) << "Scanning: " << buffer << endl;
323 stringstream line(buffer);
324 line >> FragOrder;
325 line >> ws >> No;
326 line >> ws >> Value; // skip time entry
327 line >> ws >> Value;
328 No -= 1; // indices start at 1 in file, not 0
329 //Log() << Verbose(2) << " - yields (" << No << "," << Value << ", " << FragOrder << ")" << endl;
330
331 // clean the list of those entries that have been superceded by higher order terms already
332 InsertIntoAdaptiveCriteriaList(AdaptiveCriteriaList, IndexKeySetList, FragOrder, No, Value);
333 }
334 }
335 // close and done
336 InputFile.close();
337 InputFile.clear();
338 }
339
340 return AdaptiveCriteriaList;
341};
342
343/** Maps adaptive criteria list back onto (Value, (Root Nr., Order))
344 * (i.e. sorted by value to pick the highest ones)
345 * \param *out output stream for debugging
346 * \param &AdaptiveCriteriaList list to insert into
347 * \param *mol molecule with atoms
348 * \return remapped list
349 */
350map<double, pair<int,int> > * ReMapAdaptiveCriteriaListToValue(map<int, pair<double,int> > *AdaptiveCriteriaList, molecule *mol)
351{
352 atom *Walker = NULL;
353 map<double, pair<int,int> > *FinalRootCandidates = new map<double, pair<int,int> > ;
354 DoLog(1) && (Log() << Verbose(1) << "Root candidate list is: " << endl);
355 for(map<int, pair<double,int> >::iterator runner = AdaptiveCriteriaList->begin(); runner != AdaptiveCriteriaList->end(); runner++) {
356 Walker = mol->FindAtom((*runner).first);
357 if (Walker != NULL) {
358 //if ((*runner).second.second >= Walker->AdaptiveOrder) { // only insert if this is an "active" root site for the current order
359 if (!Walker->MaxOrder) {
360 DoLog(2) && (Log() << Verbose(2) << "(" << (*runner).first << ",[" << (*runner).second.first << "," << (*runner).second.second << "])" << endl);
361 FinalRootCandidates->insert( make_pair( (*runner).second.first, pair<int,int>((*runner).first, (*runner).second.second) ) );
362 } else {
363 DoLog(2) && (Log() << Verbose(2) << "Excluding (" << *Walker << ", " << (*runner).first << ",[" << (*runner).second.first << "," << (*runner).second.second << "]), as it has reached its maximum order." << endl);
364 }
365 } else {
366 DoeLog(0) && (eLog()<< Verbose(0) << "Atom No. " << (*runner).second.first << " was not found in this molecule." << endl);
367 performCriticalExit();
368 }
369 }
370 return FinalRootCandidates;
371};
372
373/** Marks all candidate sites for update if below adaptive threshold.
374 * Picks a given number of highest values and set *AtomMask to true.
375 * \param *out output stream for debugging
376 * \param *AtomMask defines true/false per global Atom::nr to mask in/out each nuclear site, used to activate given number of site to increment order adaptively
377 * \param FinalRootCandidates list candidates to check
378 * \param Order desired order
379 * \param *mol molecule with atoms
380 * \return true - if update is necessary, false - not
381 */
382bool MarkUpdateCandidates(bool *AtomMask, map<double, pair<int,int> > &FinalRootCandidates, int Order, molecule *mol)
383{
384 atom *Walker = NULL;
385 int No = -1;
386 bool status = false;
387 for(map<double, pair<int,int> >::iterator runner = FinalRootCandidates.upper_bound(pow(10.,Order)); runner != FinalRootCandidates.end(); runner++) {
388 No = (*runner).second.first;
389 Walker = mol->FindAtom(No);
390 //if (Walker->AdaptiveOrder < MinimumRingSize[Walker->nr]) {
391 DoLog(2) && (Log() << Verbose(2) << "Root " << No << " is still above threshold (10^{" << Order <<"}: " << runner->first << ", setting entry " << No << " of Atom mask to true." << endl);
392 AtomMask[No] = true;
393 status = true;
394 //} else
395 //Log() << Verbose(2) << "Root " << No << " is still above threshold (10^{" << Order <<"}: " << runner->first << ", however MinimumRingSize of " << MinimumRingSize[Walker->nr] << " does not allow further adaptive increase." << endl;
396 }
397 return status;
398};
399
400/** print atom mask for debugging.
401 * \param *out output stream for debugging
402 * \param *AtomMask defines true/false per global Atom::nr to mask in/out each nuclear site, used to activate given number of site to increment order adaptively
403 * \param AtomCount number of entries in \a *AtomMask
404 */
405void PrintAtomMask(bool *AtomMask, int AtomCount)
406{
407 DoLog(2) && (Log() << Verbose(2) << " ");
408 for(int i=0;i<AtomCount;i++)
409 DoLog(0) && (Log() << Verbose(0) << (i % 10));
410 DoLog(0) && (Log() << Verbose(0) << endl);
411 DoLog(2) && (Log() << Verbose(2) << "Atom mask is: ");
412 for(int i=0;i<AtomCount;i++)
413 DoLog(0) && (Log() << Verbose(0) << (AtomMask[i] ? "t" : "f"));
414 DoLog(0) && (Log() << Verbose(0) << endl);
415};
416
417/** Checks whether the OrderAtSite is still below \a Order at some site.
418 * \param *AtomMask defines true/false per global Atom::nr to mask in/out each nuclear site, used to activate given number of site to increment order adaptively
419 * \param *GlobalKeySetList list of keysets with global ids (valid in "this" molecule) needed for adaptive increase
420 * \param Order desired Order if positive, desired exponent in threshold criteria if negative (0 is single-step)
421 * \param *MinimumRingSize array of max. possible order to avoid loops
422 * \param path path to ENERGYPERFRAGMENT file (may be NULL if Order is non-negative)
423 * \return true - needs further fragmentation, false - does not need fragmentation
424 */
425bool molecule::CheckOrderAtSite(bool *AtomMask, Graph *GlobalKeySetList, int Order, int *MinimumRingSize, std::string path)
426{
427 bool status = false;
428
429 // initialize mask list
430 for(int i=getAtomCount();i--;)
431 AtomMask[i] = false;
432
433 if (Order < 0) { // adaptive increase of BondOrder per site
434 if (AtomMask[getAtomCount()] == true) // break after one step
435 return false;
436
437 // transmorph graph keyset list into indexed KeySetList
438 if (GlobalKeySetList == NULL) {
439 DoeLog(1) && (eLog()<< Verbose(1) << "Given global key set list (graph) is NULL!" << endl);
440 return false;
441 }
442 map<int,KeySet> *IndexKeySetList = GraphToIndexedKeySet(GlobalKeySetList);
443
444 // parse the EnergyPerFragment file
445 map<int, pair<double,int> > *AdaptiveCriteriaList = ScanAdaptiveFileIntoMap(path, *IndexKeySetList); // (Root No., (Value, Order)) !
446 if (AdaptiveCriteriaList->empty()) {
447 DoeLog(2) && (eLog()<< Verbose(2) << "Unable to parse file, incrementing all." << endl);
448 for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
449 #ifdef ADDHYDROGEN
450 if ((*iter)->type->Z != 1) // skip hydrogen
451 #endif
452 {
453 AtomMask[(*iter)->nr] = true; // include all (non-hydrogen) atoms
454 status = true;
455 }
456 }
457 }
458 // then map back onto (Value, (Root Nr., Order)) (i.e. sorted by value to pick the highest ones)
459 map<double, pair<int,int> > *FinalRootCandidates = ReMapAdaptiveCriteriaListToValue(AdaptiveCriteriaList, this);
460
461 // pick the ones still below threshold and mark as to be adaptively updated
462 MarkUpdateCandidates(AtomMask, *FinalRootCandidates, Order, this);
463
464 delete[](IndexKeySetList);
465 delete[](AdaptiveCriteriaList);
466 delete[](FinalRootCandidates);
467 } else { // global increase of Bond Order
468 for(molecule::const_iterator iter = begin(); iter != end(); ++iter) {
469 #ifdef ADDHYDROGEN
470 if ((*iter)->type->Z != 1) // skip hydrogen
471 #endif
472 {
473 AtomMask[(*iter)->nr] = true; // include all (non-hydrogen) atoms
474 if ((Order != 0) && ((*iter)->AdaptiveOrder < Order)) // && ((*iter)->AdaptiveOrder < MinimumRingSize[(*iter)->nr]))
475 status = true;
476 }
477 }
478 if ((!Order) && (!AtomMask[getAtomCount()])) // single stepping, just check
479 status = true;
480
481 if (!status) {
482 if (Order == 0)
483 DoLog(1) && (Log() << Verbose(1) << "Single stepping done." << endl);
484 else
485 DoLog(1) && (Log() << Verbose(1) << "Order at every site is already equal or above desired order " << Order << "." << endl);
486 }
487 }
488
489 PrintAtomMask(AtomMask, getAtomCount()); // for debugging
490
491 return status;
492};
493
494/** Create a SortIndex to map from atomic labels to the sequence in which the atoms are given in the config file.
495 * \param *out output stream for debugging
496 * \param *&SortIndex Mapping array of size molecule::AtomCount
497 * \return true - success, false - failure of SortIndex alloc
498 */
499bool molecule::CreateMappingLabelsToConfigSequence(int *&SortIndex)
500{
501 if (SortIndex != NULL) {
502 DoLog(1) && (Log() << Verbose(1) << "SortIndex is " << SortIndex << " and not NULL as expected." << endl);
503 return false;
504 }
505 SortIndex = new int[getAtomCount()];
506 for(int i=getAtomCount();i--;)
507 SortIndex[i] = -1;
508
509 int AtomNo = 0;
510 SetIndexedArrayForEachAtomTo( SortIndex, &atom::nr, &IncrementalAbsoluteValue, AtomNo );
511
512 return true;
513};
514
515
516
517/** Creates a lookup table for true father's Atom::Nr -> atom ptr.
518 * \param *start begin of list (STL iterator, i.e. first item)
519 * \paran *end end of list (STL iterator, i.e. one past last item)
520 * \param **Lookuptable pointer to return allocated lookup table (should be NULL on start)
521 * \param count optional predetermined size for table (otherwise we set the count to highest true father id)
522 * \return true - success, false - failure
523 */
524bool molecule::CreateFatherLookupTable(atom **&LookupTable, int count)
525{
526 bool status = true;
527 int AtomNo;
528
529 if (LookupTable != NULL) {
530 Log() << Verbose(0) << "Pointer for Lookup table is not NULL! Aborting ..." <<endl;
531 return false;
532 }
533
534 // count them
535 if (count == 0) {
536 for (molecule::iterator iter = begin(); iter != end(); ++iter) { // create a lookup table (Atom::nr -> atom) used as a marker table lateron
537 count = (count < (*iter)->GetTrueFather()->nr) ? (*iter)->GetTrueFather()->nr : count;
538 }
539 }
540 if (count <= 0) {
541 Log() << Verbose(0) << "Count of lookup list is 0 or less." << endl;
542 return false;
543 }
544
545 // allocate and fill
546 LookupTable = new atom *[count];
547 if (LookupTable == NULL) {
548 eLog() << Verbose(0) << "LookupTable memory allocation failed!" << endl;
549 performCriticalExit();
550 status = false;
551 } else {
552 for (int i=0;i<count;i++)
553 LookupTable[i] = NULL;
554 for (molecule::iterator iter = begin(); iter != end(); ++iter) {
555 AtomNo = (*iter)->GetTrueFather()->nr;
556 if ((AtomNo >= 0) && (AtomNo < count)) {
557 //*out << "Setting LookupTable[" << AtomNo << "] to " << *(*iter) << endl;
558 LookupTable[AtomNo] = (*iter);
559 } else {
560 Log() << Verbose(0) << "Walker " << *(*iter) << " exceeded range of nuclear ids [0, " << count << ")." << endl;
561 status = false;
562 break;
563 }
564 }
565 }
566
567 return status;
568};
569
570/** Performs a many-body bond order analysis for a given bond order.
571 * -# parses adjacency, keysets and orderatsite files
572 * -# performs DFS to find connected subgraphs (to leave this in was a design decision: might be useful later)
573 * -# RootStack is created for every subgraph (here, later we implement the "update 10 sites with highest energ
574y contribution", and that's why this consciously not done in the following loop)
575 * -# in a loop over all subgraphs
576 * -# calls FragmentBOSSANOVA with this RootStack and within the subgraph molecule structure
577 * -# creates molecule (fragment)s from the returned keysets (StoreFragmentFromKeySet)
578 * -# combines the generated molecule lists from all subgraphs
579 * -# saves to disk: fragment configs, adjacency, orderatsite, keyset files
580 * Note that as we split "this" molecule up into a list of subgraphs, i.e. a MoleculeListClass, we have two sets
581 * of vertex indices: Global always means the index in "this" molecule, whereas local refers to the molecule or
582 * subgraph in the MoleculeListClass.
583 * \param Order up to how many neighbouring bonds a fragment contains in BondOrderScheme::BottumUp scheme
584 * \param &prefix path and prefix of the bond order configs to be written
585 * \return 1 - continue, 2 - stop (no fragmentation occured)
586 */
587int molecule::FragmentMolecule(int Order, std::string &prefix)
588{
589 MoleculeListClass *BondFragments = NULL;
590 int *MinimumRingSize = new int[getAtomCount()];
591 int FragmentCounter;
592 MoleculeLeafClass *MolecularWalker = NULL;
593 MoleculeLeafClass *Subgraphs = NULL; // list of subgraphs from DFS analysis
594 fstream File;
595 bool FragmentationToDo = true;
596 class StackClass<bond *> *BackEdgeStack = NULL, *LocalBackEdgeStack = NULL;
597 bool CheckOrder = false;
598 Graph **FragmentList = NULL;
599 Graph *ParsedFragmentList = NULL;
600 Graph TotalGraph; // graph with all keysets however local numbers
601 int TotalNumberOfKeySets = 0;
602 atom **ListOfAtoms = NULL;
603 atom ***ListOfLocalAtoms = NULL;
604 bool *AtomMask = NULL;
605
606 DoLog(0) && (Log() << Verbose(0) << endl);
607#ifdef ADDHYDROGEN
608 DoLog(0) && (Log() << Verbose(0) << "I will treat hydrogen special and saturate dangling bonds with it." << endl);
609#else
610 DoLog(0) && (Log() << Verbose(0) << "Hydrogen is treated just like the rest of the lot." << endl);
611#endif
612
613 // ++++++++++++++++++++++++++++ INITIAL STUFF: Bond structure analysis, file parsing, ... ++++++++++++++++++++++++++++++++++++++++++
614
615 // ===== 1. Check whether bond structure is same as stored in files ====
616
617 // create lookup table for Atom::nr
618 FragmentationToDo = FragmentationToDo && CreateFatherLookupTable(ListOfAtoms, getAtomCount());
619
620 // === compare it with adjacency file ===
621 FragmentationToDo = FragmentationToDo && CheckAdjacencyFileAgainstMolecule(prefix, ListOfAtoms);
622 delete[](ListOfAtoms);
623
624 // ===== 2. perform a DFS analysis to gather info on cyclic structure and a list of disconnected subgraphs =====
625 Subgraphs = DepthFirstSearchAnalysis(BackEdgeStack);
626
627 // analysis of the cycles (print rings, get minimum cycle length) for each subgraph
628 for(int i=getAtomCount();i--;)
629 MinimumRingSize[i] = getAtomCount();
630 MolecularWalker = Subgraphs;
631 const int LeafCount = Subgraphs->next->Count();
632 FragmentCounter = 0;
633 while (MolecularWalker->next != NULL) {
634 MolecularWalker = MolecularWalker->next;
635 // fill the bond structure of the individually stored subgraphs
636 ListOfAtoms = NULL;
637 MolecularWalker->FillBondStructureFromReference(this, ListOfAtoms, false); // we want to keep the created ListOfLocalAtoms
638 DoLog(0) && (Log() << Verbose(0) << "Analysing the cycles of subgraph " << MolecularWalker->Leaf << " with nr. " << FragmentCounter << "." << endl);
639 LocalBackEdgeStack = new StackClass<bond *> (MolecularWalker->Leaf->BondCount);
640// // check the list of local atoms for debugging
641// Log() << Verbose(0) << "ListOfLocalAtoms for this subgraph is:" << endl;
642// for (int i=0;i<getAtomCount();i++)
643// if (ListOfLocalAtoms[FragmentCounter][i] == NULL)
644// Log() << Verbose(0) << "\tNULL";
645// else
646// Log() << Verbose(0) << "\t" << ListOfLocalAtoms[FragmentCounter][i]->Name;
647 DoLog(0) && (Log() << Verbose(0) << "Gathering local back edges for subgraph " << MolecularWalker->Leaf << " with nr. " << FragmentCounter << "." << endl);
648 MolecularWalker->Leaf->PickLocalBackEdges(ListOfAtoms, BackEdgeStack, LocalBackEdgeStack);
649 DoLog(0) && (Log() << Verbose(0) << "Analysing the cycles of subgraph " << MolecularWalker->Leaf << " with nr. " << FragmentCounter << "." << endl);
650 MolecularWalker->Leaf->CyclicStructureAnalysis(LocalBackEdgeStack, MinimumRingSize);
651 DoLog(0) && (Log() << Verbose(0) << "Done with Analysing the cycles of subgraph " << MolecularWalker->Leaf << " with nr. " << FragmentCounter << "." << endl);
652 delete(LocalBackEdgeStack);
653 delete(ListOfAtoms);
654 FragmentCounter++;
655 }
656 delete(BackEdgeStack);
657
658 // ===== 3. if structure still valid, parse key set file and others =====
659 FragmentationToDo = FragmentationToDo && ParseKeySetFile(prefix, ParsedFragmentList);
660
661 // ===== 4. check globally whether there's something to do actually (first adaptivity check)
662 FragmentationToDo = FragmentationToDo && ParseOrderAtSiteFromFile(prefix);
663
664 // =================================== Begin of FRAGMENTATION ===============================
665 // ===== 6a. assign each keyset to its respective subgraph =====
666 ListOfLocalAtoms = new atom **[LeafCount];
667 for (int i=0;i<LeafCount;i++)
668 ListOfLocalAtoms[i] = NULL;
669 FragmentCounter = 0;
670 Subgraphs->next->AssignKeySetsToFragment(this, ParsedFragmentList, ListOfLocalAtoms, FragmentList, FragmentCounter, true);
671 delete[](ListOfLocalAtoms);
672
673 // ===== 6b. prepare and go into the adaptive (Order<0), single-step (Order==0) or incremental (Order>0) cycle
674 KeyStack *RootStack = new KeyStack[Subgraphs->next->Count()];
675 AtomMask = new bool[getAtomCount()+1];
676 AtomMask[getAtomCount()] = false;
677 FragmentationToDo = false; // if CheckOrderAtSite just ones recommends fragmentation, we will save fragments afterwards
678 while ((CheckOrder = CheckOrderAtSite(AtomMask, ParsedFragmentList, Order, MinimumRingSize, prefix))) {
679 FragmentationToDo = FragmentationToDo || CheckOrder;
680 AtomMask[getAtomCount()] = true; // last plus one entry is used as marker that we have been through this loop once already in CheckOrderAtSite()
681 // ===== 6b. fill RootStack for each subgraph (second adaptivity check) =====
682 Subgraphs->next->FillRootStackForSubgraphs(RootStack, AtomMask, (FragmentCounter = 0));
683
684 // ===== 7. fill the bond fragment list =====
685 FragmentCounter = 0;
686 MolecularWalker = Subgraphs;
687 while (MolecularWalker->next != NULL) {
688 MolecularWalker = MolecularWalker->next;
689 DoLog(1) && (Log() << Verbose(1) << "Fragmenting subgraph " << MolecularWalker << "." << endl);
690 //MolecularWalker->Leaf->OutputListOfBonds(out); // output atom::ListOfBonds for debugging
691 if (MolecularWalker->Leaf->hasBondStructure()) {
692 // call BOSSANOVA method
693 DoLog(0) && (Log() << Verbose(0) << endl << " ========== BOND ENERGY of subgraph " << FragmentCounter << " ========================= " << endl);
694 MolecularWalker->Leaf->FragmentBOSSANOVA(FragmentList[FragmentCounter], RootStack[FragmentCounter], MinimumRingSize);
695 } else {
696 DoeLog(1) && (eLog()<< Verbose(1) << "Subgraph " << MolecularWalker << " has no atoms!" << endl);
697 }
698 FragmentCounter++; // next fragment list
699 }
700 }
701 DoLog(2) && (Log() << Verbose(2) << "CheckOrder is " << CheckOrder << "." << endl);
702 delete[](RootStack);
703 delete[](AtomMask);
704 delete(ParsedFragmentList);
705 delete[](MinimumRingSize);
706
707 // ==================================== End of FRAGMENTATION ============================================
708
709 // ===== 8a. translate list into global numbers (i.e. ones that are valid in "this" molecule, not in MolecularWalker->Leaf)
710 Subgraphs->next->TranslateIndicesToGlobalIDs(FragmentList, (FragmentCounter = 0), TotalNumberOfKeySets, TotalGraph);
711
712 // free subgraph memory again
713 FragmentCounter = 0;
714 if (Subgraphs != NULL) {
715 while (Subgraphs->next != NULL) {
716 Subgraphs = Subgraphs->next;
717 delete(FragmentList[FragmentCounter++]);
718 delete(Subgraphs->previous);
719 }
720 delete(Subgraphs);
721 }
722 delete[](FragmentList);
723
724 // ===== 8b. gather keyset lists (graphs) from all subgraphs and transform into MoleculeListClass =====
725 //if (FragmentationToDo) { // we should always store the fragments again as coordination might have changed slightly without changing bond structure
726 // allocate memory for the pointer array and transmorph graphs into full molecular fragments
727 BondFragments = new MoleculeListClass(World::getPointer());
728 int k=0;
729 for(Graph::iterator runner = TotalGraph.begin(); runner != TotalGraph.end(); runner++) {
730 KeySet test = (*runner).first;
731 DoLog(0) && (Log() << Verbose(0) << "Fragment No." << (*runner).second.first << " with TEFactor " << (*runner).second.second << "." << endl);
732 BondFragments->insert(StoreFragmentFromKeySet(test, World::getInstance().getConfig()));
733 k++;
734 }
735 DoLog(0) && (Log() << Verbose(0) << k << "/" << BondFragments->ListOfMolecules.size() << " fragments generated from the keysets." << endl);
736
737 // ===== 9. Save fragments' configuration and keyset files et al to disk ===
738 if (BondFragments->ListOfMolecules.size() != 0) {
739 // create the SortIndex from BFS labels to order in the config file
740 int *SortIndex = NULL;
741 CreateMappingLabelsToConfigSequence(SortIndex);
742
743 DoLog(1) && (Log() << Verbose(1) << "Writing " << BondFragments->ListOfMolecules.size() << " possible bond fragmentation configs" << endl);
744 if (BondFragments->OutputConfigForListOfFragments(prefix, SortIndex))
745 DoLog(1) && (Log() << Verbose(1) << "All configs written." << endl);
746 else
747 DoLog(1) && (Log() << Verbose(1) << "Some config writing failed." << endl);
748
749 // store force index reference file
750 BondFragments->StoreForcesFile(prefix, SortIndex);
751
752 // store keysets file
753 StoreKeySetFile(TotalGraph, prefix);
754
755 {
756 // store Adjacency file
757 std::string filename = prefix + ADJACENCYFILE;
758 StoreAdjacencyToFile(filename);
759 }
760
761 // store Hydrogen saturation correction file
762 BondFragments->AddHydrogenCorrection(prefix);
763
764 // store adaptive orders into file
765 StoreOrderAtSiteFile(prefix);
766
767 // restore orbital and Stop values
768 //CalculateOrbitals(*configuration);
769
770 // free memory for bond part
771 DoLog(1) && (Log() << Verbose(1) << "Freeing bond memory" << endl);
772 delete[](SortIndex);
773 } else {
774 DoLog(1) && (Log() << Verbose(1) << "FragmentList is zero on return, splitting failed." << endl);
775 }
776 delete(BondFragments);
777 DoLog(0) && (Log() << Verbose(0) << "End of bond fragmentation." << endl);
778
779 return ((int)(!FragmentationToDo)+1); // 1 - continue, 2 - stop (no fragmentation occured)
780};
781
782
783/** Stores pairs (Atom::nr, Atom::AdaptiveOrder) into file.
784 * Atoms not present in the file get "-1".
785 * \param &path path to file ORDERATSITEFILE
786 * \return true - file writable, false - not writable
787 */
788bool molecule::StoreOrderAtSiteFile(std::string &path)
789{
790 string line;
791 ofstream file;
792
793 line = path + ORDERATSITEFILE;
794 file.open(line.c_str());
795 DoLog(1) && (Log() << Verbose(1) << "Writing OrderAtSite " << ORDERATSITEFILE << " ... " << endl);
796 if (file.good()) {
797 ActOnAllAtoms( &atom::OutputOrder, &file );
798 file.close();
799 DoLog(1) && (Log() << Verbose(1) << "done." << endl);
800 return true;
801 } else {
802 DoLog(1) && (Log() << Verbose(1) << "failed to open file " << line << "." << endl);
803 return false;
804 }
805};
806
807/** Parses pairs(Atom::nr, Atom::AdaptiveOrder) from file and stores in molecule's Atom's.
808 * Atoms not present in the file get "0".
809 * \param &path path to file ORDERATSITEFILEe
810 * \return true - file found and scanned, false - file not found
811 * \sa ParseKeySetFile() and CheckAdjacencyFileAgainstMolecule() as this is meant to be used in conjunction with the two
812 */
813bool molecule::ParseOrderAtSiteFromFile(std::string &path)
814{
815 unsigned char *OrderArray = new unsigned char[getAtomCount()];
816 bool *MaxArray = new bool[getAtomCount()];
817 bool status;
818 int AtomNr, value;
819 string line;
820 ifstream file;
821
822 for(int i=0;i<getAtomCount();i++) {
823 OrderArray[i] = 0;
824 MaxArray[i] = false;
825 }
826
827 DoLog(1) && (Log() << Verbose(1) << "Begin of ParseOrderAtSiteFromFile" << endl);
828 line = path + ORDERATSITEFILE;
829 file.open(line.c_str());
830 if (file.good()) {
831 while (!file.eof()) { // parse from file
832 AtomNr = -1;
833 file >> AtomNr;
834 if (AtomNr != -1) { // test whether we really parsed something (this is necessary, otherwise last atom is set twice and to 0 on second time)
835 file >> value;
836 OrderArray[AtomNr] = value;
837 file >> value;
838 MaxArray[AtomNr] = value;
839 //Log() << Verbose(2) << "AtomNr " << AtomNr << " with order " << (int)OrderArray[AtomNr] << " and max order set to " << (int)MaxArray[AtomNr] << "." << endl;
840 }
841 }
842 file.close();
843
844 // set atom values
845 SetAtomValueToIndexedArray( OrderArray, &atom::nr, &atom::AdaptiveOrder );
846 SetAtomValueToIndexedArray( MaxArray, &atom::nr, &atom::MaxOrder );
847
848 DoLog(1) && (Log() << Verbose(1) << "\t ... done." << endl);
849 status = true;
850 } else {
851 DoLog(1) && (Log() << Verbose(1) << "\t ... failed to open file " << line << "." << endl);
852 status = false;
853 }
854 delete[](OrderArray);
855 delete[](MaxArray);
856
857 DoLog(1) && (Log() << Verbose(1) << "End of ParseOrderAtSiteFromFile" << endl);
858 return status;
859};
860
861
862
863/** Looks through a StackClass<atom *> and returns the likeliest removal candiate.
864 * \param *out output stream for debugging messages
865 * \param *&Leaf KeySet to look through
866 * \param *&ShortestPathList list of the shortest path to decide which atom to suggest as removal candidate in the end
867 * \param index of the atom suggested for removal
868 */
869int molecule::LookForRemovalCandidate(KeySet *&Leaf, int *&ShortestPathList)
870{
871 atom *Runner = NULL;
872 int SP, Removal;
873
874 DoLog(2) && (Log() << Verbose(2) << "Looking for removal candidate." << endl);
875 SP = -1; //0; // not -1, so that Root is never removed
876 Removal = -1;
877 for (KeySet::iterator runner = Leaf->begin(); runner != Leaf->end(); runner++) {
878 Runner = FindAtom((*runner));
879 if (Runner->type->Z != 1) { // skip all those added hydrogens when re-filling snake stack
880 if (ShortestPathList[(*runner)] > SP) { // remove the oldest one with longest shortest path
881 SP = ShortestPathList[(*runner)];
882 Removal = (*runner);
883 }
884 }
885 }
886 return Removal;
887};
888
889/** Initializes some value for putting fragment of \a *mol into \a *Leaf.
890 * \param *mol total molecule
891 * \param *Leaf fragment molecule
892 * \param &Leaflet pointer to KeySet structure
893 * \param **SonList calloc'd list which atom of \a *Leaf is a son of which atom in \a *mol
894 * \return number of atoms in fragment
895 */
896int StoreFragmentFromKeySet_Init(molecule *mol, molecule *Leaf, KeySet &Leaflet, atom **SonList)
897{
898 atom *FatherOfRunner = NULL;
899
900 Leaf->BondDistance = mol->BondDistance;
901
902 // first create the minimal set of atoms from the KeySet
903 int size = 0;
904 for(KeySet::iterator runner = Leaflet.begin(); runner != Leaflet.end(); runner++) {
905 FatherOfRunner = mol->FindAtom((*runner)); // find the id
906 SonList[FatherOfRunner->nr] = Leaf->AddCopyAtom(FatherOfRunner);
907 size++;
908 }
909 return size;
910};
911
912/** Creates an induced subgraph out of a fragmental key set, adding bonds and hydrogens (if treated specially).
913 * \param *out output stream for debugging messages
914 * \param *mol total molecule
915 * \param *Leaf fragment molecule
916 * \param IsAngstroem whether we have Ansgtroem or bohrradius
917 * \param **SonList list which atom of \a *Leaf is a son of which atom in \a *mol
918 */
919void CreateInducedSubgraphOfFragment(molecule *mol, molecule *Leaf, atom **SonList, bool IsAngstroem)
920{
921 bool LonelyFlag = false;
922 atom *OtherFather = NULL;
923 atom *FatherOfRunner = NULL;
924
925#ifdef ADDHYDROGEN
926 molecule::const_iterator runner;
927#endif
928 // we increment the iter just before skipping the hydrogen
929 for (molecule::const_iterator iter = Leaf->begin(); iter != Leaf->end();) {
930 LonelyFlag = true;
931 FatherOfRunner = (*iter)->father;
932 ASSERT(FatherOfRunner,"Atom without father found");
933 if (SonList[FatherOfRunner->nr] != NULL) { // check if this, our father, is present in list
934 // create all bonds
935 for (BondList::const_iterator BondRunner = FatherOfRunner->ListOfBonds.begin(); BondRunner != FatherOfRunner->ListOfBonds.end(); (++BondRunner)) {
936 OtherFather = (*BondRunner)->GetOtherAtom(FatherOfRunner);
937// Log() << Verbose(2) << "Father " << *FatherOfRunner << " of son " << *SonList[FatherOfRunner->nr] << " is bound to " << *OtherFather;
938 if (SonList[OtherFather->nr] != NULL) {
939// Log() << Verbose(0) << ", whose son is " << *SonList[OtherFather->nr] << "." << endl;
940 if (OtherFather->nr > FatherOfRunner->nr) { // add bond (nr check is for adding only one of both variants: ab, ba)
941// Log() << Verbose(3) << "Adding Bond: ";
942// Log() << Verbose(0) <<
943 Leaf->AddBond((*iter), SonList[OtherFather->nr], (*BondRunner)->BondDegree);
944// Log() << Verbose(0) << "." << endl;
945 //NumBonds[(*iter)->nr]++;
946 } else {
947// Log() << Verbose(3) << "Not adding bond, labels in wrong order." << endl;
948 }
949 LonelyFlag = false;
950 } else {
951// Log() << Verbose(0) << ", who has no son in this fragment molecule." << endl;
952#ifdef ADDHYDROGEN
953 //Log() << Verbose(3) << "Adding Hydrogen to " << (*iter)->Name << " and a bond in between." << endl;
954 if(!Leaf->AddHydrogenReplacementAtom((*BondRunner), (*iter), FatherOfRunner, OtherFather, IsAngstroem))
955 exit(1);
956#endif
957 //NumBonds[(*iter)->nr] += Binder->BondDegree;
958 }
959 }
960 } else {
961 DoeLog(1) && (eLog()<< Verbose(1) << "Son " << (*iter)->getName() << " has father " << FatherOfRunner->getName() << " but its entry in SonList is " << SonList[FatherOfRunner->nr] << "!" << endl);
962 }
963 if ((LonelyFlag) && (Leaf->getAtomCount() > 1)) {
964 DoLog(0) && (Log() << Verbose(0) << **iter << "has got bonds only to hydrogens!" << endl);
965 }
966 ++iter;
967#ifdef ADDHYDROGEN
968 while ((iter != Leaf->end()) && ((*iter)->type->Z == 1)){ // skip added hydrogen
969 iter++;
970 }
971#endif
972 }
973};
974
975/** Stores a fragment from \a KeySet into \a molecule.
976 * First creates the minimal set of atoms from the KeySet, then creates the bond structure from the complete
977 * molecule and adds missing hydrogen where bonds were cut.
978 * \param *out output stream for debugging messages
979 * \param &Leaflet pointer to KeySet structure
980 * \param IsAngstroem whether we have Ansgtroem or bohrradius
981 * \return pointer to constructed molecule
982 */
983molecule * molecule::StoreFragmentFromKeySet(KeySet &Leaflet, bool IsAngstroem)
984{
985 atom **SonList = new atom*[getAtomCount()];
986 molecule *Leaf = World::getInstance().createMolecule();
987
988 for(int i=0;i<getAtomCount();i++)
989 SonList[i] = NULL;
990
991// Log() << Verbose(1) << "Begin of StoreFragmentFromKeyset." << endl;
992 StoreFragmentFromKeySet_Init(this, Leaf, Leaflet, SonList);
993 // create the bonds between all: Make it an induced subgraph and add hydrogen
994// Log() << Verbose(2) << "Creating bonds from father graph (i.e. induced subgraph creation)." << endl;
995 CreateInducedSubgraphOfFragment(this, Leaf, SonList, IsAngstroem);
996
997 //Leaflet->Leaf->ScanForPeriodicCorrection(out);
998 delete[](SonList);
999// Log() << Verbose(1) << "End of StoreFragmentFromKeyset." << endl;
1000 return Leaf;
1001};
1002
1003
1004/** Clears the touched list
1005 * \param *out output stream for debugging
1006 * \param verbosity verbosity level
1007 * \param *&TouchedList touched list
1008 * \param SubOrder current suborder
1009 * \param TouchedIndex currently touched
1010 */
1011void SPFragmentGenerator_ClearingTouched(int verbosity, int *&TouchedList, int SubOrder, int &TouchedIndex)
1012{
1013 Log() << Verbose(1+verbosity) << "Clearing touched list." << endl;
1014 for (TouchedIndex=SubOrder+1;TouchedIndex--;) // empty touched list
1015 TouchedList[TouchedIndex] = -1;
1016 TouchedIndex = 0;
1017
1018}
1019
1020/** Adds the current combination of the power set to the snake stack.
1021 * \param *out output stream for debugging
1022 * \param verbosity verbosity level
1023 * \param CurrentCombination
1024 * \param SetDimension maximum number of bits in power set
1025 * \param *FragmentSet snake stack to remove from
1026 * \param *&TouchedList touched list
1027 * \param TouchedIndex currently touched
1028 * \return number of set bits
1029 */
1030int AddPowersetToSnakeStack(int verbosity, int CurrentCombination, int SetDimension, KeySet *FragmentSet, bond **BondsSet, int *&TouchedList, int &TouchedIndex)
1031{
1032 atom *OtherWalker = NULL;
1033 bool bit = false;
1034 KeySetTestPair TestKeySetInsert;
1035
1036 int Added = 0;
1037 for (int j=0;j<SetDimension;j++) { // pull out every bit by shifting
1038 bit = ((CurrentCombination & (1 << j)) != 0); // mask the bit for the j-th bond
1039 if (bit) { // if bit is set, we add this bond partner
1040 OtherWalker = BondsSet[j]->rightatom; // rightatom is always the one more distant, i.e. the one to add
1041 //Log() << Verbose(1+verbosity) << "Current Bond is " << BondsSet[j] << ", checking on " << *OtherWalker << "." << endl;
1042 Log() << Verbose(2+verbosity) << "Adding " << *OtherWalker << " with nr " << OtherWalker->nr << "." << endl;
1043 TestKeySetInsert = FragmentSet->insert(OtherWalker->nr);
1044 if (TestKeySetInsert.second) {
1045 TouchedList[TouchedIndex++] = OtherWalker->nr; // note as added
1046 Added++;
1047 } else {
1048 Log() << Verbose(2+verbosity) << "This was item was already present in the keyset." << endl;
1049 }
1050 } else {
1051 Log() << Verbose(2+verbosity) << "Not adding." << endl;
1052 }
1053 }
1054 return Added;
1055};
1056
1057/** Counts the number of elements in a power set.
1058 * \param *SetFirst
1059 * \param *SetLast
1060 * \param *&TouchedList touched list
1061 * \param TouchedIndex currently touched
1062 * \return number of elements
1063 */
1064int CountSetMembers(bond *SetFirst, bond *SetLast, int *&TouchedList, int TouchedIndex)
1065{
1066 int SetDimension = 0;
1067 bond *Binder = SetFirst; // start node for this level
1068 while (Binder->next != SetLast) { // compare to end node of this level
1069 Binder = Binder->next;
1070 for (int k=TouchedIndex;k--;) {
1071 if (Binder->Contains(TouchedList[k])) // if we added this very endpiece
1072 SetDimension++;
1073 }
1074 }
1075 return SetDimension;
1076};
1077
1078/** Counts the number of elements in a power set.
1079 * \param *BondsList bonds list to fill
1080 * \param *SetFirst
1081 * \param *SetLast
1082 * \param *&TouchedList touched list
1083 * \param TouchedIndex currently touched
1084 * \return number of elements
1085 */
1086int FillBondsList(bond **BondsList, bond *SetFirst, bond *SetLast, int *&TouchedList, int TouchedIndex)
1087{
1088 int SetDimension = 0;
1089 bond *Binder = SetFirst; // start node for this level
1090 while (Binder->next != SetLast) { // compare to end node of this level
1091 Binder = Binder->next;
1092 for (int k=0;k<TouchedIndex;k++) {
1093 if (Binder->leftatom->nr == TouchedList[k]) // leftatom is always the close one
1094 BondsList[SetDimension++] = Binder;
1095 }
1096 }
1097 return SetDimension;
1098};
1099
1100/** Remove all items that were added on this SP level.
1101 * \param *out output stream for debugging
1102 * \param verbosity verbosity level
1103 * \param *FragmentSet snake stack to remove from
1104 * \param *&TouchedList touched list
1105 * \param TouchedIndex currently touched
1106 */
1107void RemoveAllTouchedFromSnakeStack(int verbosity, KeySet *FragmentSet, int *&TouchedList, int &TouchedIndex)
1108{
1109 int Removal = 0;
1110 for(int j=0;j<TouchedIndex;j++) {
1111 Removal = TouchedList[j];
1112 Log() << Verbose(2+verbosity) << "Removing item nr. " << Removal << " from snake stack." << endl;
1113 FragmentSet->erase(Removal);
1114 TouchedList[j] = -1;
1115 }
1116 DoLog(2) && (Log() << Verbose(2) << "Remaining local nr.s on snake stack are: ");
1117 for(KeySet::iterator runner = FragmentSet->begin(); runner != FragmentSet->end(); runner++)
1118 DoLog(0) && (Log() << Verbose(0) << (*runner) << " ");
1119 DoLog(0) && (Log() << Verbose(0) << endl);
1120 TouchedIndex = 0; // set Index to 0 for list of atoms added on this level
1121};
1122
1123/** From a given set of Bond sorted by Shortest Path distance, create all possible fragments of size \a SetDimension.
1124 * -# loops over every possible combination (2^dimension of edge set)
1125 * -# inserts current set, if there's still space left
1126 * -# yes: calls SPFragmentGenerator with structure, created new edge list and size respective to root dist
1127ance+1
1128 * -# no: stores fragment into keyset list by calling InsertFragmentIntoGraph
1129 * -# removes all items added into the snake stack (in UniqueFragments structure) added during level (root
1130distance) and current set
1131 * \param *out output stream for debugging
1132 * \param FragmentSearch UniqueFragments structure with all values needed
1133 * \param RootDistance current shortest path level, whose set of edges is represented by **BondsSet
1134 * \param SetDimension Number of possible bonds on this level (i.e. size of the array BondsSet[])
1135 * \param SubOrder remaining number of allowed vertices to add
1136 */
1137void molecule::SPFragmentGenerator(struct UniqueFragments *FragmentSearch, int RootDistance, bond **BondsSet, int SetDimension, int SubOrder)
1138{
1139 int verbosity = 0; //FragmentSearch->ANOVAOrder-SubOrder;
1140 int NumCombinations;
1141 int bits, TouchedIndex, SubSetDimension, SP, Added;
1142 int SpaceLeft;
1143 int *TouchedList = new int[SubOrder + 1];
1144 KeySetTestPair TestKeySetInsert;
1145
1146 NumCombinations = 1 << SetDimension;
1147
1148 // here for all bonds of Walker all combinations of end pieces (from the bonds)
1149 // have to be added and for the remaining ANOVA order GraphCrawler be called
1150 // recursively for the next level
1151
1152 Log() << Verbose(1+verbosity) << "Begin of SPFragmentGenerator." << endl;
1153 Log() << Verbose(1+verbosity) << "We are " << RootDistance << " away from Root, which is " << *FragmentSearch->Root << ", SubOrder is " << SubOrder << ", SetDimension is " << SetDimension << " and this means " << NumCombinations-1 << " combination(s)." << endl;
1154
1155 // initialised touched list (stores added atoms on this level)
1156 SPFragmentGenerator_ClearingTouched(verbosity, TouchedList, SubOrder, TouchedIndex);
1157
1158 // create every possible combination of the endpieces
1159 Log() << Verbose(1+verbosity) << "Going through all combinations of the power set." << endl;
1160 for (int i=1;i<NumCombinations;i++) { // sweep through all power set combinations (skip empty set!)
1161 // count the set bit of i
1162 bits = 0;
1163 for (int j=SetDimension;j--;)
1164 bits += (i & (1 << j)) >> j;
1165
1166 Log() << Verbose(1+verbosity) << "Current set is " << Binary(i | (1 << SetDimension)) << ", number of bits is " << bits << "." << endl;
1167 if (bits <= SubOrder) { // if not greater than additional atoms allowed on stack, continue
1168 // --1-- add this set of the power set of bond partners to the snake stack
1169 Added = AddPowersetToSnakeStack(verbosity, i, SetDimension, FragmentSearch->FragmentSet, BondsSet, TouchedList, TouchedIndex);
1170
1171 SpaceLeft = SubOrder - Added ;// SubOrder - bits; // due to item's maybe being already present, this does not work anymore
1172 if (SpaceLeft > 0) {
1173 Log() << Verbose(1+verbosity) << "There's still some space left on stack: " << SpaceLeft << "." << endl;
1174 if (SubOrder > 1) { // Due to Added above we have to check extra whether we're not already reaching beyond the desired Order
1175 // --2-- look at all added end pieces of this combination, construct bond subsets and sweep through a power set of these by recursion
1176 SP = RootDistance+1; // this is the next level
1177
1178 // first count the members in the subset
1179 SubSetDimension = CountSetMembers(FragmentSearch->BondsPerSPList[2*SP], FragmentSearch->BondsPerSPList[2*SP+1], TouchedList, TouchedIndex);
1180
1181 // then allocate and fill the list
1182 bond *BondsList[SubSetDimension];
1183 SubSetDimension = FillBondsList(BondsList, FragmentSearch->BondsPerSPList[2*SP], FragmentSearch->BondsPerSPList[2*SP+1], TouchedList, TouchedIndex);
1184
1185 // then iterate
1186 Log() << Verbose(2+verbosity) << "Calling subset generator " << SP << " away from root " << *FragmentSearch->Root << " with sub set dimension " << SubSetDimension << "." << endl;
1187 SPFragmentGenerator(FragmentSearch, SP, BondsList, SubSetDimension, SubOrder-bits);
1188 }
1189 } else {
1190 // --2-- otherwise store the complete fragment
1191 Log() << Verbose(1+verbosity) << "Enough items on stack for a fragment!" << endl;
1192 // store fragment as a KeySet
1193 DoLog(2) && (Log() << Verbose(2) << "Found a new fragment[" << FragmentSearch->FragmentCounter << "], local nr.s are: ");
1194 for(KeySet::iterator runner = FragmentSearch->FragmentSet->begin(); runner != FragmentSearch->FragmentSet->end(); runner++)
1195 DoLog(0) && (Log() << Verbose(0) << (*runner) << " ");
1196 DoLog(0) && (Log() << Verbose(0) << endl);
1197 //if (!CheckForConnectedSubgraph(FragmentSearch->FragmentSet))
1198 //DoeLog(1) && (eLog()<< Verbose(1) << "The found fragment is not a connected subgraph!" << endl);
1199 InsertFragmentIntoGraph(FragmentSearch);
1200 }
1201
1202 // --3-- remove all added items in this level from snake stack
1203 Log() << Verbose(1+verbosity) << "Removing all items that were added on this SP level " << RootDistance << "." << endl;
1204 RemoveAllTouchedFromSnakeStack(verbosity, FragmentSearch->FragmentSet, TouchedList, TouchedIndex);
1205 } else {
1206 Log() << Verbose(2+verbosity) << "More atoms to add for this set (" << bits << ") than space left on stack " << SubOrder << ", skipping this set." << endl;
1207 }
1208 }
1209 delete[](TouchedList);
1210 Log() << Verbose(1+verbosity) << "End of SPFragmentGenerator, " << RootDistance << " away from Root " << *FragmentSearch->Root << " and SubOrder is " << SubOrder << "." << endl;
1211};
1212
1213/** Allocates memory for UniqueFragments::BondsPerSPList.
1214 * \param *out output stream
1215 * \param Order bond order (limits BFS exploration and "number of digits" in power set generation
1216 * \param FragmentSearch UniqueFragments
1217 * \sa FreeSPList()
1218 */
1219void InitialiseSPList(int Order, struct UniqueFragments &FragmentSearch)
1220{
1221 FragmentSearch.BondsPerSPList = new bond* [Order * 2];
1222 FragmentSearch.BondsPerSPCount = new int[Order];
1223 for (int i=Order;i--;) {
1224 FragmentSearch.BondsPerSPList[2*i] = new bond(); // start node
1225 FragmentSearch.BondsPerSPList[2*i+1] = new bond(); // end node
1226 FragmentSearch.BondsPerSPList[2*i]->next = FragmentSearch.BondsPerSPList[2*i+1]; // intertwine these two
1227 FragmentSearch.BondsPerSPList[2*i+1]->previous = FragmentSearch.BondsPerSPList[2*i];
1228 FragmentSearch.BondsPerSPCount[i] = 0;
1229 }
1230};
1231
1232/** Free's memory for for UniqueFragments::BondsPerSPList.
1233 * \param *out output stream
1234 * \param Order bond order (limits BFS exploration and "number of digits" in power set generation
1235 * \param FragmentSearch UniqueFragments\
1236 * \sa InitialiseSPList()
1237 */
1238void FreeSPList(int Order, struct UniqueFragments &FragmentSearch)
1239{
1240 delete[](FragmentSearch.BondsPerSPCount);
1241 for (int i=Order;i--;) {
1242 delete(FragmentSearch.BondsPerSPList[2*i]);
1243 delete(FragmentSearch.BondsPerSPList[2*i+1]);
1244 }
1245 delete[](FragmentSearch.BondsPerSPList);
1246};
1247
1248/** Sets FragmenSearch to initial value.
1249 * Sets UniqueFragments::ShortestPathList entries to zero, UniqueFragments::BondsPerSPCount to zero (except zero level to 1) and
1250 * adds initial bond UniqueFragments::Root to UniqueFragments::Root to UniqueFragments::BondsPerSPList
1251 * \param *out output stream
1252 * \param Order bond order (limits BFS exploration and "number of digits" in power set generation
1253 * \param FragmentSearch UniqueFragments
1254 * \sa FreeSPList()
1255 */
1256void SetSPList(int Order, struct UniqueFragments &FragmentSearch)
1257{
1258 // prepare Label and SP arrays of the BFS search
1259 FragmentSearch.ShortestPathList[FragmentSearch.Root->nr] = 0;
1260
1261 // prepare root level (SP = 0) and a loop bond denoting Root
1262 for (int i=Order;i--;)
1263 FragmentSearch.BondsPerSPCount[i] = 0;
1264 FragmentSearch.BondsPerSPCount[0] = 1;
1265 bond *Binder = new bond(FragmentSearch.Root, FragmentSearch.Root);
1266 add(Binder, FragmentSearch.BondsPerSPList[1]);
1267};
1268
1269/** Resets UniqueFragments::ShortestPathList and cleans bonds from UniqueFragments::BondsPerSPList.
1270 * \param *out output stream
1271 * \param Order bond order (limits BFS exploration and "number of digits" in power set generation
1272 * \param FragmentSearch UniqueFragments
1273 * \sa InitialiseSPList()
1274 */
1275void ResetSPList(int Order, struct UniqueFragments &FragmentSearch)
1276{
1277 bond *Binder = NULL;
1278 DoLog(0) && (Log() << Verbose(0) << "Free'ing all found lists. and resetting index lists" << endl);
1279 for(int i=Order;i--;) {
1280 DoLog(1) && (Log() << Verbose(1) << "Current SP level is " << i << ": ");
1281 Binder = FragmentSearch.BondsPerSPList[2*i];
1282 while (Binder->next != FragmentSearch.BondsPerSPList[2*i+1]) {
1283 Binder = Binder->next;
1284 // Log() << Verbose(0) << "Removing atom " << Binder->leftatom->nr << " and " << Binder->rightatom->nr << "." << endl; // make sure numbers are local
1285 FragmentSearch.ShortestPathList[Binder->leftatom->nr] = -1;
1286 FragmentSearch.ShortestPathList[Binder->rightatom->nr] = -1;
1287 }
1288 // delete added bonds
1289 cleanup(FragmentSearch.BondsPerSPList[2*i], FragmentSearch.BondsPerSPList[2*i+1]);
1290 // also start and end node
1291 DoLog(0) && (Log() << Verbose(0) << "cleaned." << endl);
1292 }
1293};
1294
1295
1296/** Fills the Bonds per Shortest Path List and set the vertex labels.
1297 * \param *out output stream
1298 * \param Order bond order (limits BFS exploration and "number of digits" in power set generation
1299 * \param FragmentSearch UniqueFragments
1300 * \param *mol molecule with atoms and bonds
1301 * \param RestrictedKeySet Restricted vertex set to use in context of molecule
1302 */
1303void FillSPListandLabelVertices(int Order, struct UniqueFragments &FragmentSearch, molecule *mol, KeySet RestrictedKeySet)
1304{
1305 // Actually, we should construct a spanning tree vom the root atom and select all edges therefrom and put them into
1306 // according shortest path lists. However, we don't. Rather we fill these lists right away, as they do form a spanning
1307 // tree already sorted into various SP levels. That's why we just do loops over the depth (CurrentSP) and breadth
1308 // (EdgeinSPLevel) of this tree ...
1309 // In another picture, the bonds always contain a direction by rightatom being the one more distant from root and hence
1310 // naturally leftatom forming its predecessor, preventing the BFS"seeker" from continuing in the wrong direction.
1311 int AtomKeyNr = -1;
1312 atom *Walker = NULL;
1313 atom *OtherWalker = NULL;
1314 atom *Predecessor = NULL;
1315 bond *CurrentEdge = NULL;
1316 bond *Binder = NULL;
1317 int RootKeyNr = FragmentSearch.Root->GetTrueFather()->nr;
1318 int RemainingWalkers = -1;
1319 int SP = -1;
1320
1321 DoLog(0) && (Log() << Verbose(0) << "Starting BFS analysis ..." << endl);
1322 for (SP = 0; SP < (Order-1); SP++) {
1323 DoLog(1) && (Log() << Verbose(1) << "New SP level reached: " << SP << ", creating new SP list with " << FragmentSearch.BondsPerSPCount[SP] << " item(s)");
1324 if (SP > 0) {
1325 DoLog(0) && (Log() << Verbose(0) << ", old level closed with " << FragmentSearch.BondsPerSPCount[SP-1] << " item(s)." << endl);
1326 FragmentSearch.BondsPerSPCount[SP] = 0;
1327 } else
1328 DoLog(0) && (Log() << Verbose(0) << "." << endl);
1329
1330 RemainingWalkers = FragmentSearch.BondsPerSPCount[SP];
1331 CurrentEdge = FragmentSearch.BondsPerSPList[2*SP]; /// start of this SP level's list
1332 while (CurrentEdge->next != FragmentSearch.BondsPerSPList[2*SP+1]) { /// end of this SP level's list
1333 CurrentEdge = CurrentEdge->next;
1334 RemainingWalkers--;
1335 Walker = CurrentEdge->rightatom; // rightatom is always the one more distant
1336 Predecessor = CurrentEdge->leftatom; // ... and leftatom is predecessor
1337 AtomKeyNr = Walker->nr;
1338 DoLog(0) && (Log() << Verbose(0) << "Current Walker is: " << *Walker << " with nr " << Walker->nr << " and SP of " << SP << ", with " << RemainingWalkers << " remaining walkers on this level." << endl);
1339 // check for new sp level
1340 // go through all its bonds
1341 DoLog(1) && (Log() << Verbose(1) << "Going through all bonds of Walker." << endl);
1342 for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
1343 OtherWalker = (*Runner)->GetOtherAtom(Walker);
1344 if ((RestrictedKeySet.find(OtherWalker->nr) != RestrictedKeySet.end())
1345 #ifdef ADDHYDROGEN
1346 && (OtherWalker->type->Z != 1)
1347 #endif
1348 ) { // skip hydrogens and restrict to fragment
1349 DoLog(2) && (Log() << Verbose(2) << "Current partner is " << *OtherWalker << " with nr " << OtherWalker->nr << " in bond " << *(*Runner) << "." << endl);
1350 // set the label if not set (and push on root stack as well)
1351 if ((OtherWalker != Predecessor) && (OtherWalker->GetTrueFather()->nr > RootKeyNr)) { // only pass through those with label bigger than Root's
1352 FragmentSearch.ShortestPathList[OtherWalker->nr] = SP+1;
1353 DoLog(3) && (Log() << Verbose(3) << "Set Shortest Path to " << FragmentSearch.ShortestPathList[OtherWalker->nr] << "." << endl);
1354 // add the bond in between to the SP list
1355 Binder = new bond(Walker, OtherWalker); // create a new bond in such a manner, that bond::rightatom is always the one more distant
1356 add(Binder, FragmentSearch.BondsPerSPList[2*(SP+1)+1]);
1357 FragmentSearch.BondsPerSPCount[SP+1]++;
1358 DoLog(3) && (Log() << Verbose(3) << "Added its bond to SP list, having now " << FragmentSearch.BondsPerSPCount[SP+1] << " item(s)." << endl);
1359 } else {
1360 if (OtherWalker != Predecessor)
1361 DoLog(3) && (Log() << Verbose(3) << "Not passing on, as index of " << *OtherWalker << " " << OtherWalker->GetTrueFather()->nr << " is smaller than that of Root " << RootKeyNr << "." << endl);
1362 else
1363 DoLog(3) && (Log() << Verbose(3) << "This is my predecessor " << *Predecessor << "." << endl);
1364 }
1365 } else Log() << Verbose(2) << "Is not in the restricted keyset or skipping hydrogen " << *OtherWalker << "." << endl;
1366 }
1367 }
1368 }
1369};
1370
1371/** prints the Bonds per Shortest Path list in UniqueFragments.
1372 * \param *out output stream
1373 * \param Order bond order (limits BFS exploration and "number of digits" in power set generation
1374 * \param FragmentSearch UniqueFragments
1375 */
1376void OutputSPList(int Order, struct UniqueFragments &FragmentSearch)
1377{
1378 bond *Binder = NULL;
1379 DoLog(0) && (Log() << Verbose(0) << "Printing all found lists." << endl);
1380 for(int i=1;i<Order;i++) { // skip the root edge in the printing
1381 Binder = FragmentSearch.BondsPerSPList[2*i];
1382 DoLog(1) && (Log() << Verbose(1) << "Current SP level is " << i << "." << endl);
1383 while (Binder->next != FragmentSearch.BondsPerSPList[2*i+1]) {
1384 Binder = Binder->next;
1385 DoLog(2) && (Log() << Verbose(2) << *Binder << endl);
1386 }
1387 }
1388};
1389
1390/** Simply counts all bonds in all UniqueFragments::BondsPerSPList lists.
1391 * \param *out output stream
1392 * \param Order bond order (limits BFS exploration and "number of digits" in power set generation
1393 * \param FragmentSearch UniqueFragments
1394 */
1395int CountNumbersInBondsList(int Order, struct UniqueFragments &FragmentSearch)
1396{
1397 bond *Binder = NULL;
1398 int SP = -1; // the Root <-> Root edge must be subtracted!
1399 for(int i=Order;i--;) { // sum up all found edges
1400 Binder = FragmentSearch.BondsPerSPList[2*i];
1401 while (Binder->next != FragmentSearch.BondsPerSPList[2*i+1]) {
1402 Binder = Binder->next;
1403 SP++;
1404 }
1405 }
1406 return SP;
1407};
1408
1409/** Creates a list of all unique fragments of certain vertex size from a given graph \a Fragment for a given root vertex in the context of \a this molecule.
1410 * -# initialises UniqueFragments structure
1411 * -# fills edge list via BFS
1412 * -# creates the fragment by calling recursive function SPFragmentGenerator with UniqueFragments structure, 0 as
1413 root distance, the edge set, its dimension and the current suborder
1414 * -# Free'ing structure
1415 * Note that we may use the fact that the atoms are SP-ordered on the atomstack. I.e. when popping always the last, we first get all
1416 * with SP of 2, then those with SP of 3, then those with SP of 4 and so on.
1417 * \param *out output stream for debugging
1418 * \param Order bond order (limits BFS exploration and "number of digits" in power set generation
1419 * \param FragmentSearch UniqueFragments structure containing TEFactor, root atom and so on
1420 * \param RestrictedKeySet Restricted vertex set to use in context of molecule
1421 * \return number of inserted fragments
1422 * \note ShortestPathList in FragmentSearch structure is probably due to NumberOfAtomsSPLevel and SP not needed anymore
1423 */
1424int molecule::PowerSetGenerator(int Order, struct UniqueFragments &FragmentSearch, KeySet RestrictedKeySet)
1425{
1426 int Counter = FragmentSearch.FragmentCounter; // mark current value of counter
1427
1428 DoLog(0) && (Log() << Verbose(0) << endl);
1429 DoLog(0) && (Log() << Verbose(0) << "Begin of PowerSetGenerator with order " << Order << " at Root " << *FragmentSearch.Root << "." << endl);
1430
1431 SetSPList(Order, FragmentSearch);
1432
1433 // do a BFS search to fill the SP lists and label the found vertices
1434 FillSPListandLabelVertices(Order, FragmentSearch, this, RestrictedKeySet);
1435
1436 // outputting all list for debugging
1437 OutputSPList(Order, FragmentSearch);
1438
1439 // creating fragments with the found edge sets (may be done in reverse order, faster)
1440 int SP = CountNumbersInBondsList(Order, FragmentSearch);
1441 DoLog(0) && (Log() << Verbose(0) << "Total number of edges is " << SP << "." << endl);
1442 if (SP >= (Order-1)) {
1443 // start with root (push on fragment stack)
1444 DoLog(0) && (Log() << Verbose(0) << "Starting fragment generation with " << *FragmentSearch.Root << ", local nr is " << FragmentSearch.Root->nr << "." << endl);
1445 FragmentSearch.FragmentSet->clear();
1446 DoLog(0) && (Log() << Verbose(0) << "Preparing subset for this root and calling generator." << endl);
1447
1448 // prepare the subset and call the generator
1449 bond* BondsList[FragmentSearch.BondsPerSPCount[0]];
1450 for(int i=0;i<FragmentSearch.BondsPerSPCount[0];i++)
1451 BondsList[i] = NULL;
1452 BondsList[0] = FragmentSearch.BondsPerSPList[0]->next; // on SP level 0 there's only the root bond
1453
1454 SPFragmentGenerator(&FragmentSearch, 0, BondsList, FragmentSearch.BondsPerSPCount[0], Order);
1455 } else {
1456 DoLog(0) && (Log() << Verbose(0) << "Not enough total number of edges to build " << Order << "-body fragments." << endl);
1457 }
1458
1459 // as FragmentSearch structure is used only once, we don't have to clean it anymore
1460 // remove root from stack
1461 DoLog(0) && (Log() << Verbose(0) << "Removing root again from stack." << endl);
1462 FragmentSearch.FragmentSet->erase(FragmentSearch.Root->nr);
1463
1464 // free'ing the bonds lists
1465 ResetSPList(Order, FragmentSearch);
1466
1467 // return list
1468 DoLog(0) && (Log() << Verbose(0) << "End of PowerSetGenerator." << endl);
1469 return (FragmentSearch.FragmentCounter - Counter);
1470};
1471
1472bool KeyCompare::operator() (const KeySet SubgraphA, const KeySet SubgraphB) const
1473{
1474 //Log() << Verbose(0) << "my check is used." << endl;
1475 if (SubgraphA.size() < SubgraphB.size()) {
1476 return true;
1477 } else {
1478 if (SubgraphA.size() > SubgraphB.size()) {
1479 return false;
1480 } else {
1481 KeySet::iterator IteratorA = SubgraphA.begin();
1482 KeySet::iterator IteratorB = SubgraphB.begin();
1483 while ((IteratorA != SubgraphA.end()) && (IteratorB != SubgraphB.end())) {
1484 if ((*IteratorA) < (*IteratorB))
1485 return true;
1486 else if ((*IteratorA) > (*IteratorB)) {
1487 return false;
1488 } // else, go on to next index
1489 IteratorA++;
1490 IteratorB++;
1491 } // end of while loop
1492 }// end of check in case of equal sizes
1493 }
1494 return false; // if we reach this point, they are equal
1495};
1496
1497
1498/** Combines all KeySets from all orders into single ones (with just unique entries).
1499 * \param *out output stream for debugging
1500 * \param *&FragmentList list to fill
1501 * \param ***FragmentLowerOrdersList
1502 * \param &RootStack stack with all root candidates (unequal to each atom in complete molecule if adaptive scheme is applied)
1503 * \param *mol molecule with atoms and bonds
1504 */
1505int CombineAllOrderListIntoOne(Graph *&FragmentList, Graph ***FragmentLowerOrdersList, KeyStack &RootStack, molecule *mol)
1506{
1507 int RootNr = 0;
1508 int RootKeyNr = 0;
1509 int StartNr = 0;
1510 int counter = 0;
1511 int NumLevels = 0;
1512 atom *Walker = NULL;
1513
1514 DoLog(0) && (Log() << Verbose(0) << "Combining the lists of all orders per order and finally into a single one." << endl);
1515 if (FragmentList == NULL) {
1516 FragmentList = new Graph;
1517 counter = 0;
1518 } else {
1519 counter = FragmentList->size();
1520 }
1521
1522 StartNr = RootStack.back();
1523 do {
1524 RootKeyNr = RootStack.front();
1525 RootStack.pop_front();
1526 Walker = mol->FindAtom(RootKeyNr);
1527 NumLevels = 1 << (Walker->AdaptiveOrder - 1);
1528 for(int i=0;i<NumLevels;i++) {
1529 if (FragmentLowerOrdersList[RootNr][i] != NULL) {
1530 InsertGraphIntoGraph(*FragmentList, (*FragmentLowerOrdersList[RootNr][i]), &counter);
1531 }
1532 }
1533 RootStack.push_back(Walker->nr);
1534 RootNr++;
1535 } while (RootKeyNr != StartNr);
1536 return counter;
1537};
1538
1539/** Free's memory allocated for all KeySets from all orders.
1540 * \param *out output stream for debugging
1541 * \param ***FragmentLowerOrdersList
1542 * \param &RootStack stack with all root candidates (unequal to each atom in complete molecule if adaptive scheme is applied)
1543 * \param *mol molecule with atoms and bonds
1544 */
1545void FreeAllOrdersList(Graph ***FragmentLowerOrdersList, KeyStack &RootStack, molecule *mol)
1546{
1547 DoLog(1) && (Log() << Verbose(1) << "Free'ing the lists of all orders per order." << endl);
1548 int RootNr = 0;
1549 int RootKeyNr = 0;
1550 int NumLevels = 0;
1551 atom *Walker = NULL;
1552 while (!RootStack.empty()) {
1553 RootKeyNr = RootStack.front();
1554 RootStack.pop_front();
1555 Walker = mol->FindAtom(RootKeyNr);
1556 NumLevels = 1 << (Walker->AdaptiveOrder - 1);
1557 for(int i=0;i<NumLevels;i++) {
1558 if (FragmentLowerOrdersList[RootNr][i] != NULL) {
1559 delete(FragmentLowerOrdersList[RootNr][i]);
1560 }
1561 }
1562 delete[](FragmentLowerOrdersList[RootNr]);
1563 RootNr++;
1564 }
1565 delete[](FragmentLowerOrdersList);
1566};
1567
1568
1569/** Performs BOSSANOVA decomposition at selected sites, increasing the cutoff by one at these sites.
1570 * -# constructs a complete keyset of the molecule
1571 * -# In a loop over all possible roots from the given rootstack
1572 * -# increases order of root site
1573 * -# calls PowerSetGenerator with this order, the complete keyset and the rootkeynr
1574 * -# for all consecutive lower levels PowerSetGenerator is called with the suborder, the higher order keyset
1575as the restricted one and each site in the set as the root)
1576 * -# these are merged into a fragment list of keysets
1577 * -# All fragment lists (for all orders, i.e. from all destination fields) are merged into one list for return
1578 * Important only is that we create all fragments, it is not important if we create them more than once
1579 * as these copies are filtered out via use of the hash table (KeySet).
1580 * \param *out output stream for debugging
1581 * \param Fragment&*List list of already present keystacks (adaptive scheme) or empty list
1582 * \param &RootStack stack with all root candidates (unequal to each atom in complete molecule if adaptive scheme is applied)
1583 * \param *MinimumRingSize minimum ring size for each atom (molecule::Atomcount)
1584 * \return pointer to Graph list
1585 */
1586void molecule::FragmentBOSSANOVA(Graph *&FragmentList, KeyStack &RootStack, int *MinimumRingSize)
1587{
1588 Graph ***FragmentLowerOrdersList = NULL;
1589 int NumLevels = 0;
1590 int NumMolecules = 0;
1591 int TotalNumMolecules = 0;
1592 int *NumMoleculesOfOrder = NULL;
1593 int Order = 0;
1594 int UpgradeCount = RootStack.size();
1595 KeyStack FragmentRootStack;
1596 int RootKeyNr = 0;
1597 int RootNr = 0;
1598 struct UniqueFragments FragmentSearch;
1599
1600 DoLog(0) && (Log() << Verbose(0) << "Begin of FragmentBOSSANOVA." << endl);
1601
1602 // FragmentLowerOrdersList is a 2D-array of pointer to MoleculeListClass objects, one dimension represents the ANOVA expansion of a single order (i.e. 5)
1603 // with all needed lower orders that are subtracted, the other dimension is the BondOrder (i.e. from 1 to 5)
1604 NumMoleculesOfOrder = new int[UpgradeCount];
1605 FragmentLowerOrdersList = new Graph**[UpgradeCount];
1606
1607 for(int i=0;i<UpgradeCount;i++) {
1608 NumMoleculesOfOrder[i] = 0;
1609 FragmentLowerOrdersList[i] = NULL;
1610 }
1611
1612 // initialise the fragments structure
1613 FragmentSearch.FragmentCounter = 0;
1614 FragmentSearch.FragmentSet = new KeySet;
1615 FragmentSearch.Root = FindAtom(RootKeyNr);
1616 FragmentSearch.ShortestPathList = new int[getAtomCount()];
1617 for (int i=getAtomCount();i--;) {
1618 FragmentSearch.ShortestPathList[i] = -1;
1619 }
1620
1621 // Construct the complete KeySet which we need for topmost level only (but for all Roots)
1622 KeySet CompleteMolecule;
1623 for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
1624 CompleteMolecule.insert((*iter)->GetTrueFather()->nr);
1625 }
1626
1627 // this can easily be seen: if Order is 5, then the number of levels for each lower order is the total sum of the number of levels above, as
1628 // each has to be split up. E.g. for the second level we have one from 5th, one from 4th, two from 3th (which in turn is one from 5th, one from 4th),
1629 // hence we have overall four 2th order levels for splitting. This also allows for putting all into a single array (FragmentLowerOrdersList[])
1630 // with the order along the cells as this: 5433222211111111 for BondOrder 5 needing 16=pow(2,5-1) cells (only we use bit-shifting which is faster)
1631 RootNr = 0; // counts through the roots in RootStack
1632 while ((RootNr < UpgradeCount) && (!RootStack.empty())) {
1633 RootKeyNr = RootStack.front();
1634 RootStack.pop_front();
1635 atom *Walker = FindAtom(RootKeyNr);
1636 // check cyclic lengths
1637 //if ((MinimumRingSize[Walker->GetTrueFather()->nr] != -1) && (Walker->GetTrueFather()->AdaptiveOrder+1 > MinimumRingSize[Walker->GetTrueFather()->nr])) {
1638 // Log() << Verbose(0) << "Bond order " << Walker->GetTrueFather()->AdaptiveOrder << " of Root " << *Walker << " greater than or equal to Minimum Ring size of " << MinimumRingSize << " found is not allowed." << endl;
1639 //} else
1640 {
1641 // increase adaptive order by one
1642 Walker->GetTrueFather()->AdaptiveOrder++;
1643 Order = Walker->AdaptiveOrder = Walker->GetTrueFather()->AdaptiveOrder;
1644
1645 // initialise Order-dependent entries of UniqueFragments structure
1646 InitialiseSPList(Order, FragmentSearch);
1647
1648 // allocate memory for all lower level orders in this 1D-array of ptrs
1649 NumLevels = 1 << (Order-1); // (int)pow(2,Order);
1650 FragmentLowerOrdersList[RootNr] = new Graph*[NumLevels];
1651 for (int i=0;i<NumLevels;i++)
1652 FragmentLowerOrdersList[RootNr][i] = NULL;
1653
1654 // create top order where nothing is reduced
1655 DoLog(0) && (Log() << Verbose(0) << "==============================================================================================================" << endl);
1656 DoLog(0) && (Log() << Verbose(0) << "Creating KeySets of Bond Order " << Order << " for " << *Walker << ", " << (RootStack.size()-RootNr) << " Roots remaining." << endl); // , NumLevels is " << NumLevels << "
1657
1658 // Create list of Graphs of current Bond Order (i.e. F_{ij})
1659 FragmentLowerOrdersList[RootNr][0] = new Graph;
1660 FragmentSearch.TEFactor = 1.;
1661 FragmentSearch.Leaflet = FragmentLowerOrdersList[RootNr][0]; // set to insertion graph
1662 FragmentSearch.Root = Walker;
1663 NumMoleculesOfOrder[RootNr] = PowerSetGenerator(Walker->AdaptiveOrder, FragmentSearch, CompleteMolecule);
1664
1665 // output resulting number
1666 DoLog(1) && (Log() << Verbose(1) << "Number of resulting KeySets is: " << NumMoleculesOfOrder[RootNr] << "." << endl);
1667 if (NumMoleculesOfOrder[RootNr] != 0) {
1668 NumMolecules = 0;
1669 } else {
1670 Walker->GetTrueFather()->MaxOrder = true;
1671 }
1672 // now, we have completely filled each cell of FragmentLowerOrdersList[] for the current Walker->AdaptiveOrder
1673 //NumMoleculesOfOrder[Walker->AdaptiveOrder-1] = NumMolecules;
1674 TotalNumMolecules += NumMoleculesOfOrder[RootNr];
1675// Log() << Verbose(1) << "Number of resulting molecules for Order " << (int)Walker->GetTrueFather()->AdaptiveOrder << " is: " << NumMoleculesOfOrder[RootNr] << "." << endl;
1676 RootStack.push_back(RootKeyNr); // put back on stack
1677 RootNr++;
1678
1679 // free Order-dependent entries of UniqueFragments structure for next loop cycle
1680 FreeSPList(Order, FragmentSearch);
1681 }
1682 }
1683 DoLog(0) && (Log() << Verbose(0) << "==============================================================================================================" << endl);
1684 DoLog(1) && (Log() << Verbose(1) << "Total number of resulting molecules is: " << TotalNumMolecules << "." << endl);
1685 DoLog(0) && (Log() << Verbose(0) << "==============================================================================================================" << endl);
1686
1687 // cleanup FragmentSearch structure
1688 delete[](FragmentSearch.ShortestPathList);
1689 delete(FragmentSearch.FragmentSet);
1690
1691 // now, FragmentLowerOrdersList is complete, it looks - for BondOrder 5 - as this (number is the ANOVA Order of the terms therein)
1692 // 5433222211111111
1693 // 43221111
1694 // 3211
1695 // 21
1696 // 1
1697
1698 // Subsequently, we combine all into a single list (FragmentList)
1699 CombineAllOrderListIntoOne(FragmentList, FragmentLowerOrdersList, RootStack, this);
1700 FreeAllOrdersList(FragmentLowerOrdersList, RootStack, this);
1701 delete[](NumMoleculesOfOrder);
1702
1703 DoLog(0) && (Log() << Verbose(0) << "End of FragmentBOSSANOVA." << endl);
1704};
1705
1706/** Corrects the nuclei position if the fragment was created over the cell borders.
1707 * Scans all bonds, checks the distance, if greater than typical, we have a candidate for the correction.
1708 * We remove the bond whereafter the graph probably separates. Then, we translate the one component periodically
1709 * and re-add the bond. Looping on the distance check.
1710 * \param *out ofstream for debugging messages
1711 */
1712void molecule::ScanForPeriodicCorrection()
1713{
1714 bond *Binder = NULL;
1715 bond *OtherBinder = NULL;
1716 atom *Walker = NULL;
1717 atom *OtherWalker = NULL;
1718 double * const cell_size = World::getInstance().getDomain();
1719 double *matrix = ReturnFullMatrixforSymmetric(cell_size);
1720 enum Shading *ColorList = NULL;
1721 double tmp;
1722 Vector Translationvector;
1723 //class StackClass<atom *> *CompStack = NULL;
1724 class StackClass<atom *> *AtomStack = new StackClass<atom *>(getAtomCount());
1725 bool flag = true;
1726
1727 DoLog(2) && (Log() << Verbose(2) << "Begin of ScanForPeriodicCorrection." << endl);
1728
1729 ColorList = new enum Shading[getAtomCount()];
1730 for (int i=0;i<getAtomCount();i++)
1731 ColorList[i] = (enum Shading)0;
1732 while (flag) {
1733 // remove bonds that are beyond bonddistance
1734 Translationvector.Zero();
1735 // scan all bonds
1736 flag = false;
1737 for(molecule::iterator AtomRunner = begin(); (!flag) && (AtomRunner != end()); ++AtomRunner)
1738 for(BondList::iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); (!flag) && (BondRunner != (*AtomRunner)->ListOfBonds.end()); ++BondRunner) {
1739 Binder = (*BondRunner);
1740 for (int i=NDIM;i--;) {
1741 tmp = fabs(Binder->leftatom->x[i] - Binder->rightatom->x[i]);
1742 //Log() << Verbose(3) << "Checking " << i << "th distance of " << *Binder->leftatom << " to " << *Binder->rightatom << ": " << tmp << "." << endl;
1743 if (tmp > BondDistance) {
1744 OtherBinder = Binder->next; // note down binding partner for later re-insertion
1745 unlink(Binder); // unlink bond
1746 DoLog(2) && (Log() << Verbose(2) << "Correcting at bond " << *Binder << "." << endl);
1747 flag = true;
1748 break;
1749 }
1750 }
1751 }
1752 if (flag) {
1753 // create translation vector from their periodically modified distance
1754 for (int i=NDIM;i--;) {
1755 tmp = Binder->leftatom->x[i] - Binder->rightatom->x[i];
1756 if (fabs(tmp) > BondDistance)
1757 Translationvector[i] = (tmp < 0) ? +1. : -1.;
1758 }
1759 Translationvector.MatrixMultiplication(matrix);
1760 //Log() << Verbose(3) << "Translation vector is ";
1761 Log() << Verbose(0) << Translationvector << endl;
1762 // apply to all atoms of first component via BFS
1763 for (int i=getAtomCount();i--;)
1764 ColorList[i] = white;
1765 AtomStack->Push(Binder->leftatom);
1766 while (!AtomStack->IsEmpty()) {
1767 Walker = AtomStack->PopFirst();
1768 //Log() << Verbose (3) << "Current Walker is: " << *Walker << "." << endl;
1769 ColorList[Walker->nr] = black; // mark as explored
1770 Walker->x += Translationvector; // translate
1771 for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
1772 if ((*Runner) != Binder) {
1773 OtherWalker = (*Runner)->GetOtherAtom(Walker);
1774 if (ColorList[OtherWalker->nr] == white) {
1775 AtomStack->Push(OtherWalker); // push if yet unexplored
1776 }
1777 }
1778 }
1779 }
1780 // re-add bond
1781 link(Binder, OtherBinder);
1782 } else {
1783 DoLog(3) && (Log() << Verbose(3) << "No corrections for this fragment." << endl);
1784 }
1785 //delete(CompStack);
1786 }
1787 // free allocated space from ReturnFullMatrixforSymmetric()
1788 delete(AtomStack);
1789 delete[](ColorList);
1790 delete[](matrix);
1791 DoLog(2) && (Log() << Verbose(2) << "End of ScanForPeriodicCorrection." << endl);
1792};
Note: See TracBrowser for help on using the repository browser.