source: src/Fragmentation/Fragmentation.cpp@ 9eb71b3

Action_Thermostats Add_AtomRandomPerturbation Add_RotateAroundBondAction Add_SelectAtomByNameAction Adding_Graph_to_ChangeBondActions Adding_MD_integration_tests Adding_StructOpt_integration_tests AutomationFragmentation_failures Candidate_v1.6.1 ChangeBugEmailaddress ChangingTestPorts ChemicalSpaceEvaluator Disabling_MemDebug Docu_Python_wait EmpiricalPotential_contain_HomologyGraph_documentation Enhance_userguide Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Example_ManyWaysToTranslateAtom Exclude_Hydrogens_annealWithBondGraph Fix_ChronosMutex Fix_StatusMsg Fix_StepWorldTime_single_argument Fix_Verbose_Codepatterns ForceAnnealing_goodresults ForceAnnealing_oldresults ForceAnnealing_tocheck ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_continued ForceAnnealing_with_BondGraph_continued_betteresults ForceAnnealing_with_BondGraph_contraction-expansion GeometryObjects Gui_displays_atomic_force_velocity IndependentFragmentGrids_IntegrationTest JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool PythonUI_with_named_parameters QtGui_reactivate_TimeChanged_changes Recreated_GuiChecks RotateToPrincipalAxisSystem_UndoRedo StoppableMakroAction TremoloParser_IncreasedPrecision TremoloParser_MultipleTimesteps Ubuntu_1604_changes
Last change on this file since 9eb71b3 was 9eb71b3, checked in by Frederik Heber <frederik.heber@…>, 8 years ago

Commented out MemDebug include and Memory::ignore.

  • MemDebug clashes with various allocation operators that use a specific placement in memory. It is so far not possible to wrap new/delete fully. Hence, we stop this effort which so far has forced us to put ever more includes (with clashes) into MemDebug and thereby bloat compilation time.
  • MemDebug does not add that much usefulness which is not also provided by valgrind.
  • Property mode set to 100644
File size: 26.0 KB
Line 
1/*
2 * Project: MoleCuilder
3 * Description: creates and alters molecular systems
4 * Copyright (C) 2010-2012 University of Bonn. All rights reserved.
5 * Copyright (C) 2013 Frederik Heber. All rights reserved.
6 *
7 *
8 * This file is part of MoleCuilder.
9 *
10 * MoleCuilder is free software: you can redistribute it and/or modify
11 * it under the terms of the GNU General Public License as published by
12 * the Free Software Foundation, either version 2 of the License, or
13 * (at your option) any later version.
14 *
15 * MoleCuilder is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU General Public License for more details.
19 *
20 * You should have received a copy of the GNU General Public License
21 * along with MoleCuilder. If not, see <http://www.gnu.org/licenses/>.
22 */
23
24/*
25 * Fragmentation.cpp
26 *
27 * Created on: Oct 18, 2011
28 * Author: heber
29 */
30
31#ifdef HAVE_CONFIG_H
32#include <config.h>
33#endif
34
35#include <boost/bimap.hpp>
36
37//#include "CodePatterns/MemDebug.hpp"
38
39#include "Fragmentation.hpp"
40
41#include "CodePatterns/Assert.hpp"
42#include "CodePatterns/Info.hpp"
43#include "CodePatterns/IteratorAdaptors.hpp"
44#include "CodePatterns/Log.hpp"
45
46#include "Atom/atom.hpp"
47#include "Bond/bond.hpp"
48#include "Descriptors/MoleculeDescriptor.hpp"
49#include "Element/element.hpp"
50#include "Element/periodentafel.hpp"
51#include "Fragmentation/AdaptivityMap.hpp"
52#include "Fragmentation/AtomMask.hpp"
53#include "Fragmentation/fragmentation_helpers.hpp"
54#include "Fragmentation/Graph.hpp"
55#include "Fragmentation/helpers.hpp"
56#include "Fragmentation/KeySet.hpp"
57#include "Fragmentation/PowerSetGenerator.hpp"
58#include "Fragmentation/UniqueFragments.hpp"
59#include "Graph/BondGraph.hpp"
60#include "Graph/AdjacencyList.hpp"
61#include "Graph/ListOfLocalAtoms.hpp"
62#include "molecule.hpp"
63#include "World.hpp"
64
65
66/** Constructor of class Fragmentation.
67 *
68 * \param _mol molecule which we currently fragment
69 * \param _FileChecker instance contains adjacency parsed from elsewhere
70 * \param _treatment whether to treat hydrogen special and saturate dangling bonds or not
71 */
72Fragmentation::Fragmentation(molecule *_mol, AdjacencyList &_FileChecker, const enum HydrogenTreatment _treatment) :
73 mol(_mol),
74 treatment(_treatment),
75 FileChecker(_FileChecker)
76{}
77
78/** Destructor of class Fragmentation.
79 *
80 */
81Fragmentation::~Fragmentation()
82{}
83
84
85/** Performs a many-body bond order analysis for a given bond order.
86 *
87 * \note during fragmentation we switch to so-called local ids, atomic ids
88 * that are valid only for the specific molecule (representing a connected
89 * subgraph of the molecular system).
90 *
91 * -# create the local id to global id mapping
92 * -# parse the adjacency file and require the above mapping for translation
93 * -# initialize a mask for the molecule's atoms, telling which atoms are
94 * treated and which atoms neglected during fragmentation
95 * -# parse and orderatsite file and check whether there's something to do. This
96 allows for iterative calls to fragmentation
97 * -# fragments from the last fragmentation stored to file are converted into
98 * keysets (sets of atomic indices that describe one fragment)
99 * -# in a loop as long as order at site is not correct yet
100 * -# prepare a stack containing the initial ids where the fragmentation or
101 * rather the graph algorithms start from
102 * -# call fragmentBOSSANOVA()
103 * -# afterwards in case we don't saturate we remove single-atom fragments
104 * -# translate the local ids of the keysets into global ids.
105 * -# updated order at site file is written
106 *
107 * note that all created fragments or rather their describing key sets are
108 * contained in the Graph Fragmentation::FragmentList.
109 *
110 * \param atomids atomic ids (local to Fragmentation::mol) to fragment, used in AtomMask
111 * \param Order up to how many neighbouring bonds a fragment contains in BondOrderScheme::BottumUp scheme
112 * \param prefix prefix string for every fragment file name (may include path)
113 * \param ParsedFragmentList all already created key sets
114 * \return 1 - continue, 2 - stop (no fragmentation occured)
115 */
116int Fragmentation::FragmentMolecule(
117 const std::vector<atomId_t> &atomids,
118 int Order,
119 const std::string &prefix,
120 const Graph &ParsedFragmentList,
121 const bool _ParseStateFile)
122{
123 std::fstream File;
124 bool CheckOrder = false;
125 int TotalNumberOfKeySets = 0;
126
127 LOG(0, std::endl);
128 switch (treatment) {
129 case ExcludeHydrogen:
130 LOG(1, "INFO: I will treat hydrogen special.");
131 break;
132 case IncludeHydrogen:
133 LOG(1, "INFO: Hydrogen is treated just like the rest of the lot.");
134 break;
135 default:
136 ASSERT(0, "Fragmentation::FragmentMolecule() - there is a HydrogenTreatment setting which I have no idea about.");
137 break;
138 }
139
140 // ++++++++++++++++++++++++++++ INITIAL STUFF: Bond structure analysis, file parsing, ... ++++++++++++++++++++++++++++++++++++++++++
141 bool FragmentationToDo = true;
142
143 // ===== 1. Check whether bond structure is same as stored in files ====
144
145 // create a lookup global <-> local id for atomic ids valid in World and in molecule
146 Global_local_bimap_t Global_local_bimap;
147 for (std::vector<local_t>::const_iterator iter = atomids.begin();
148 iter != atomids.end();
149 ++iter) {
150 const atom * _atom = mol->FindAtom(*iter);
151 ASSERT( _atom != NULL,
152 "Fragmentation::FragmentMolecule() - could not find atom "+toString(*iter)+".");
153 Global_local_bimap.insert(
154 idpair_t(
155 (global_t)(_atom->getId()), (local_t)(*iter)
156 )
157 );
158 }
159
160 // === compare it with adjacency file ===
161 {
162 const std::vector<atomId_t> globalids(
163 MapKeyConstIterator<Global_local_bimap_t::left_const_iterator>(Global_local_bimap.left.begin()),
164 MapKeyConstIterator<Global_local_bimap_t::left_const_iterator>(Global_local_bimap.left.end())
165 );
166 AdjacencyList WorldAdjacency(globalids);
167 FragmentationToDo = FragmentationToDo && (FileChecker > WorldAdjacency);
168 }
169
170 // ===== 2. create AtomMask that takes Saturation condition into account
171 AtomMask_t AtomMask(atomids);
172 for (molecule::const_iterator iter = const_cast<const molecule *>(mol)->begin();
173 iter != const_cast<const molecule *>(mol)->end();
174 ++iter) {
175 // remove in hydrogen and we do saturate
176 if ((treatment == ExcludeHydrogen) && ((*iter)->getType()->getAtomicNumber() == 1)) // skip hydrogen
177 AtomMask.setFalse((*iter)->getNr());
178 }
179
180 // ===== 4. check globally whether there's something to do actually (first adaptivity check)
181 if (_ParseStateFile)
182 FragmentationToDo = FragmentationToDo && ParseOrderAtSiteFromFile(atomids, prefix, Global_local_bimap);
183
184 // =================================== Begin of FRAGMENTATION ===============================
185 // ===== 6a. assign each keyset to its respective subgraph =====
186 ListOfLocalAtoms_t ListOfLocalAtoms;
187 Graph FragmentList;
188 AssignKeySetsToFragment(ParsedFragmentList, ListOfLocalAtoms, FragmentList, true);
189
190 // ===== 6b. prepare and go into the adaptive (Order<0), single-step (Order==0) or incremental (Order>0) cycle
191 KeyStack RootStack;
192 FragmentationToDo = false; // if CheckOrderAtSite just ones recommends fragmentation, we will save fragments afterwards
193 bool LoopDoneAlready = false;
194 while ((CheckOrder = CheckOrderAtSite(AtomMask, ParsedFragmentList, Order, prefix, LoopDoneAlready))) {
195 FragmentationToDo = FragmentationToDo || CheckOrder;
196 LoopDoneAlready = true; // last plus one entry is used as marker that we have been through this loop once already in CheckOrderAtSite()
197 // ===== 6b. fill RootStack for each subgraph (second adaptivity check) =====
198 FillRootStackForSubgraphs(RootStack, AtomMask);
199
200 // call BOSSANOVA method
201 FragmentBOSSANOVA(mol, FragmentList, RootStack);
202 }
203 LOG(3, "DEBUG: CheckOrder is " << CheckOrder << ".");
204
205 // ==================================== End of FRAGMENTATION ============================================
206
207 // if hydrogen is not treated special, we may have single hydrogens and other
208 // fragments which are note single-determinant. These need to be removed
209 if (treatment == IncludeHydrogen) {
210 // remove all single atom fragments from FragmentList
211 Graph::iterator iter = FragmentList.begin();
212 while ( iter != FragmentList.end()) {
213 if ((*iter).first.size() == 1) {
214 LOG(1, "INFO: Removing KeySet " << (*iter).first << " as is not single-determinant.");
215 Graph::iterator eraseiter = iter++;
216 FragmentList.erase(eraseiter);
217 } else
218 ++iter;
219 }
220 }
221
222 // ===== 8a. translate list into global numbers (i.e. ones that are valid in "this" molecule, not in MolecularWalker->Leaf)
223 TranslateIndicesToGlobalIDs(FragmentList, TotalNumberOfKeySets, TotalGraph);
224
225 LOG(1, "STATUS: We have created " << TotalGraph.size() << " fragments.");
226
227
228 // store adaptive orders into file
229 StoreOrderAtSiteFile(prefix);
230
231 return ((int)(!FragmentationToDo)+1); // 1 - continue, 2 - stop (no fragmentation occured)
232};
233
234
235/** Performs BOSSANOVA decomposition at selected sites, increasing the cutoff by one at these sites.
236 * -# constructs a complete keyset of the molecule
237 * -# In a loop over all possible roots from the given rootstack
238 * -# increases order of root site
239 * -# calls PowerSetGenerator with this order, the complete keyset and the rootkeynr
240 * -# for all consecutive lower levels PowerSetGenerator is called with the suborder, the higher order keyset
241as the restricted one and each site in the set as the root)
242 * -# these are merged into a fragment list of keysets
243 * -# All fragment lists (for all orders, i.e. from all destination fields) are merged into one list for return
244 * Important only is that we create all fragments, it is not important if we create them more than once
245 * as these copies are filtered out via use of the hash table (KeySet).
246 * \param *out output stream for debugging
247 * \param Fragment&*List list of already present keystacks (adaptive scheme) or empty list
248 * \param &RootStack stack with all root candidates (unequal to each atom in complete molecule if adaptive scheme is applied)
249 * \return pointer to Graph list
250 */
251void Fragmentation::FragmentBOSSANOVA(molecule *mol, Graph &FragmentList, KeyStack &RootStack)
252{
253 Info FunctionInfo(__func__);
254 std::vector<Graph*> *FragmentLowerOrdersList = NULL;
255 size_t NumLevels = 0;
256// size_t NumMolecules = 0;
257 size_t TotalNumMolecules = 0;
258 int *NumMoleculesOfOrder = NULL;
259 int Order = 0;
260 int UpgradeCount = RootStack.size();
261 KeyStack FragmentRootStack;
262 int RootKeyNr = 0;
263 int RootNr = 0;
264
265 // FragmentLowerOrdersList is a 2D-array of pointer to vector of molecule objects, one dimension represents the ANOVA expansion of a single order (i.e. 5)
266 // with all needed lower orders that are subtracted, the other dimension is the BondOrder (i.e. from 1 to 5)
267 NumMoleculesOfOrder = new int[UpgradeCount];
268 FragmentLowerOrdersList = new std::vector<Graph*>[UpgradeCount];
269
270 for(int i=0;i<UpgradeCount;i++)
271 NumMoleculesOfOrder[i] = 0;
272
273 // Construct the restricted key set (contained in the root stack), i.e. which atoms take part
274 KeySet RestrictedKeyset;
275 for (KeyStack::const_iterator iter = RootStack.begin(); iter != RootStack.end(); ++iter)
276 RestrictedKeyset.insert(*iter);
277
278 // 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
279 // 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),
280 // hence we have overall four 2th order levels for splitting. This also allows for putting all into a single array (FragmentLowerOrdersList[])
281 // 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)
282 RootNr = 0; // counts through the roots in RootStack
283 while ((RootNr < UpgradeCount) && (!RootStack.empty())) {
284 RootKeyNr = RootStack.front();
285 RootStack.pop_front();
286 atom *Walker = mol->FindAtom(RootKeyNr);
287 // check cyclic lengths
288 //if ((MinimumRingSize[Walker->GetTrueFather()->getNr()] != -1) && (Walker->GetTrueFather()->getAdaptiveOrder()+1 > MinimumRingSize[Walker->GetTrueFather()->getNr()])) {
289 // LOG(0, "Bond order " << Walker->GetTrueFather()->getAdaptiveOrder() << " of Root " << *Walker << " greater than or equal to Minimum Ring size of " << MinimumRingSize << " found is not allowed.");
290 //} else
291 {
292 // set adaptive order to desired max order
293 Walker->GetTrueFather()->setAdaptiveOrder(Walker->GetTrueFather()->getMaxOrder());
294 Order = Walker->GetTrueFather()->getAdaptiveOrder();
295 Walker->setAdaptiveOrder(Order);
296
297 // allocate memory for all lower level orders
298 NumLevels = Order;
299 FragmentLowerOrdersList[RootNr].resize(NumLevels, NULL);
300 for( size_t i=0;i<NumLevels;++i)
301 FragmentLowerOrdersList[RootNr][i] = new Graph;
302
303 // initialise Order-dependent entries of UniqueFragments structure
304 UniqueFragments FragmentSearch(1., FragmentLowerOrdersList[RootNr], Walker);
305 PowerSetGenerator PSG(&FragmentSearch, Walker->getAdaptiveOrder());
306
307 // create top order where nothing is reduced
308 LOG(0, "==============================================================================================================");
309 LOG(0, "Creating KeySets up till Bond Order " << Order << " for " << *Walker << ", " << (RootStack.size()-RootNr) << " Roots remaining."); // , NumLevels is " << NumLevels << "
310
311 // Create list of Graphs of current Bond Order (i.e. F_{ij})
312 NumMoleculesOfOrder[RootNr] = PSG(RestrictedKeyset, treatment);
313
314 // output resulting number
315 LOG(1, "INFO: Number of resulting KeySets is: " << NumMoleculesOfOrder[RootNr] << ".");
316// if (NumMoleculesOfOrder[RootNr] != 0) {
317// NumMolecules = 0;
318// }
319 // now, we have completely filled each cell of FragmentLowerOrdersList[] for the current Walker->getAdaptiveOrder()
320 //NumMoleculesOfOrder[Walker->getAdaptiveOrder()-1] = NumMolecules;
321 TotalNumMolecules += NumMoleculesOfOrder[RootNr];
322// LOG(1, "Number of resulting molecules for Order " << (int)Walker->GetTrueFather()->getAdaptiveOrder() << " is: " << NumMoleculesOfOrder[RootNr] << ".");
323 RootStack.push_back(RootKeyNr); // put back on stack
324 RootNr++;
325 }
326 }
327 LOG(0, "==============================================================================================================");
328 LOG(0, "\tTotal number of resulting fragments is: " << TotalNumMolecules << ".");
329 LOG(0, "==============================================================================================================");
330
331 // now, FragmentLowerOrdersList is complete, it looks - for BondOrder 5 - as this (number is the ANOVA Order of the terms therein)
332 // 5433222211111111
333 // 43221111
334 // 3211
335 // 21
336 // 1
337
338 // Subsequently, we combine all into a single list (FragmentList)
339 CombineAllOrderListIntoOne(FragmentList, FragmentLowerOrdersList, RootStack, mol);
340 FreeAllOrdersList(FragmentLowerOrdersList, RootStack, mol);
341 delete[](NumMoleculesOfOrder);
342};
343
344/** Estimates by educated guessing (using upper limit) the expected number of fragments.
345 * The upper limit is
346 * \f[
347 * n = N \cdot C^k
348 * \f]
349 * where \f$C=2^c\f$ and c is the maximum bond degree over N number of atoms.
350 * \param *out output stream for debugging
351 * \param order bond order k
352 * \return number n of fragments
353 */
354int Fragmentation::GuesstimateFragmentCount(int order)
355{
356 size_t c = 0;
357 int FragmentCount;
358 // get maximum bond degree
359 for (molecule::const_iterator iter = const_cast<const molecule *>(mol)->begin();
360 iter != const_cast<const molecule *>(mol)->end();
361 ++iter) {
362 const BondList& ListOfBonds = (*iter)->getListOfBonds();
363 c = (ListOfBonds.size() > c) ? ListOfBonds.size() : c;
364 }
365 FragmentCount = (treatment == ExcludeHydrogen ? mol->getNoNonHydrogen() : mol->getAtomCount()) *(1 << (c*order));
366 LOG(1, "INFO: Upper limit for this subgraph is " << FragmentCount << " for "
367 << mol->getNoNonHydrogen() << " non-H atoms with maximum bond degree of " << c << ".");
368 return FragmentCount;
369};
370
371
372/** Checks whether the OrderAtSite is still below \a Order at some site.
373 * \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
374 * \param *GlobalKeySetList list of keysets with global ids (valid in "this" molecule) needed for adaptive increase
375 * \param Order desired Order if positive, desired exponent in threshold criteria if negative (0 is single-step)
376 * \param path path to ENERGYPERFRAGMENT file (may be NULL if Order is non-negative)
377 * \param LoopDoneAlready indicate whether we have done a fragmentation loop already
378 * \return true - needs further fragmentation, false - does not need fragmentation
379 */
380bool Fragmentation::CheckOrderAtSite(AtomMask_t &AtomMask, const Graph &GlobalKeySetList, int Order, const std::string &path, bool LoopDoneAlready)
381{
382 bool status = false;
383
384 if (Order < 0) { // adaptive increase of BondOrder per site
385 if (LoopDoneAlready) // break after one step
386 return false;
387
388 // transmorph graph keyset list into indexed KeySetList
389 AdaptivityMap * IndexKeySetList = GlobalKeySetList.GraphToAdaptivityMap();
390
391 // parse the EnergyPerFragment file
392 IndexKeySetList->ScanAdaptiveFileIntoMap(path); // (Root No., (Value, Order)) !
393 // then map back onto (Value, (Root Nr., Order)) (i.e. sorted by value to pick the highest ones)
394 IndexKeySetList->ReMapAdaptiveCriteriaListToValue(mol);
395
396 // pick the ones still below threshold and mark as to be adaptively updated
397 if (IndexKeySetList->IsAdaptiveCriteriaListEmpty()) {
398 ELOG(2, "Unable to parse file, incrementing all.");
399 status = true;
400 } else {
401 // mark as false all sites that are below threshold already
402 status = IndexKeySetList->MarkUpdateCandidates(AtomMask, Order, mol);
403 }
404
405 delete[](IndexKeySetList);
406 } else { // global increase of Bond Order
407 for(molecule::iterator iter = mol->begin(); iter != mol->end(); ++iter) {
408 atom * const Walker = *iter;
409 if (AtomMask.isTrue(Walker->getNr())) { // skip masked out
410 Walker->setMaxOrder((Order != 0 ? Order : Walker->getMaxOrder()+1));
411 // remove all that have reached desired order
412 if (Walker->getAdaptiveOrder() >= Walker->getMaxOrder()) // && (Walker->getAdaptiveOrder() < MinimumRingSize[Walker->getNr()]))
413 AtomMask.setFalse(Walker->getNr());
414 else
415 status = true;
416 }
417 }
418 if ((!Order) && (!LoopDoneAlready)) // single stepping, just check
419 status = true;
420
421 if (!status) {
422 if (Order == 0)
423 LOG(1, "INFO: Single stepping done.");
424 else
425 LOG(1, "INFO: Order at every site is already equal or above desired order " << Order << ".");
426 }
427 }
428
429 PrintAtomMask(AtomMask, AtomMask.size()); // for debugging
430
431 return status;
432};
433
434/** Stores pairs (Atom::Nr, Atom::AdaptiveOrder) into file.
435 * Atoms not present in the file get "-1".
436 * \param &path path to file ORDERATSITEFILE
437 * \return true - file writable, false - not writable
438 */
439bool Fragmentation::StoreOrderAtSiteFile(
440 const std::string &path)
441{
442 string line;
443 ofstream file;
444
445 line = path + ORDERATSITEFILE;
446 file.open(line.c_str(), std::ofstream::out | std::ofstream::app);
447 std::stringstream output;
448 output << "INFO: Writing OrderAtSite " << ORDERATSITEFILE << " ... ";
449 if (file.good()) {
450 for (molecule::const_iterator iter = const_cast<const molecule *>(mol)->begin();
451 iter != const_cast<const molecule *>(mol)->end();
452 ++iter) {
453 file << (*iter)->getId()
454 << "\t" << (int)(*iter)->getAdaptiveOrder()
455 << "\t" << (int)(*iter)->getMaxOrder() << std::endl;
456 }
457 file.close();
458 output << "done.";
459 return true;
460 } else {
461 output << "failed to open file " << line << ".";
462 return false;
463 }
464 LOG(1, output.str());
465 return true;
466};
467
468
469/** Parses pairs(Atom::Nr, Atom::AdaptiveOrder) from file and stores in molecule's Atom's.
470 * Atoms not present in the file get "0".
471 * \param atomids atoms to fragment, used in AtomMask
472 * \param &path path to file ORDERATSITEFILE
473 * \param global_local_bimap translate global to local id
474 * \return true - file found and scanned, false - file not found
475 * \sa ParseKeySetFile() and CheckAdjacencyFileAgainstMolecule() as this is meant to be used in conjunction with the two
476 */
477bool Fragmentation::ParseOrderAtSiteFromFile(
478 const std::vector<atomId_t> &atomids,
479 const std::string &path,
480 const Global_local_bimap_t &global_local_bimap)
481{
482// Info FunctionInfo(__func__);
483 typedef unsigned char order_t;
484 typedef std::map<atomId_t, order_t> OrderArray_t;
485 OrderArray_t OrderArray;
486 AtomMask_t MaxArray(atomids);
487 bool status;
488 int AtomNr, ordervalue, maxvalue;
489 string line;
490 ifstream file;
491
492 line = path + ORDERATSITEFILE;
493 file.open(line.c_str());
494 if (file.good()) {
495 while (!file.eof()) { // parse from file
496 AtomNr = -1;
497 file >> AtomNr;
498 file >> ordervalue;
499 file >> maxvalue;
500 if (AtomNr != -1) { // test whether we really parsed something (this is necessary, otherwise last atom is set twice and to 0 on second time)
501 // parsed id is global, must be translated to local id
502 Global_local_bimap_t::left_const_iterator iter = global_local_bimap.left.find(AtomNr);
503 // skip global ids we don't know about, must be in other molecule
504 if (iter != global_local_bimap.left.end()) {
505 const int LocalId = iter->second;
506 OrderArray[LocalId] = ordervalue;
507 MaxArray.setValue(LocalId, (bool)maxvalue);
508 //LOG(2, "AtomNr " << LocalId << " with order " << (int)OrderArray[LocalId] << " and max order set to " << (int)MaxArray[LocalId] << ".");
509 }
510 }
511 }
512 file.close();
513
514 // set atom values
515 for(molecule::iterator iter=mol->begin();iter!=mol->end();++iter){
516 (*iter)->setAdaptiveOrder(OrderArray[(*iter)->getNr()]);
517 (*iter)->setMaxOrder(OrderArray[(*iter)->getNr()]); //MaxArray.isTrue((*iter)->getNr());
518 }
519 //SetAtomValueToIndexedArray( OrderArray, &atom::getNr(), &atom::AdaptiveOrder );
520 //SetAtomValueToIndexedArray( MaxArray, &atom::getNr(), &atom::MaxOrder );
521
522 status = true;
523 } else {
524 ELOG(1, "Failed to open OrdersAtSite file " << line << ".");
525 status = false;
526 }
527
528 return status;
529};
530
531/** Fills the root stack for sites to be used as root in fragmentation depending on order or adaptivity criteria
532 * Again, as in \sa FillBondStructureFromReference steps recursively through each Leaf in this chain list of molecule's.
533 * \param &RootStack stack to be filled
534 * \param AtomMask defines true/false per global Atom::Nr to mask in/out each nuclear site
535 * \return true - stack is non-empty, fragmentation necessary, false - stack is empty, no more sites to update
536 */
537void Fragmentation::FillRootStackForSubgraphs(KeyStack &RootStack, const AtomMask_t &AtomMask)
538{
539 for(molecule::const_iterator iter = const_cast<const molecule *>(mol)->begin();
540 iter != const_cast<const molecule *>(mol)->end();
541 ++iter) {
542 const atom * const Father = (*iter)->GetTrueFather();
543 if (AtomMask.isTrue(Father->getNr())) // apply mask
544 if ((treatment == IncludeHydrogen) || ((*iter)->getType()->getAtomicNumber() != 1)) // skip hydrogen
545 RootStack.push_front((*iter)->getNr());
546 }
547}
548
549/** The indices per keyset are compared to the respective father's Atom::Nr in each subgraph and thus put into \a **&FragmentList.
550 * \param *KeySetList list with all keysets
551 * \param ListOfLocalAtoms Lookup table for each subgraph and index of each atom in global molecule, may be NULL on start, then it is filled
552 * \param **&FragmentList list to be allocated and returned
553 * \param FreeList true - ***ListOfLocalAtoms is free'd before return, false - it is not
554 * \retuen true - success, false - failure
555 */
556bool Fragmentation::AssignKeySetsToFragment(const Graph &KeySetList, ListOfLocalAtoms_t &ListOfLocalAtoms, Graph &FragmentList, bool FreeList)
557{
558// Info FunctionInfo(__func__);
559 bool status = true;
560 size_t KeySetCounter = 0;
561
562 // fill ListOfLocalAtoms if NULL was given
563 if (!mol->FillListOfLocalAtoms(ListOfLocalAtoms, mol->getAtomCount())) {
564 ELOG(1, "Filling of ListOfLocalAtoms failed.");
565 return false;
566 }
567
568 if (KeySetList.size() != 0) { // if there are some scanned keysets at all
569 // assign scanned keysets
570 KeySet TempSet;
571 for (Graph::const_iterator runner = KeySetList.begin(); runner != KeySetList.end(); runner++) { // key sets contain global numbers!
572 if (ListOfLocalAtoms[mol->FindAtom(*((*runner).first.begin()))->getNr()] != NULL) {// as we may assume that that bond structure is unchanged, we only test the first key in each set
573 // translate keyset to local numbers
574 for (KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++)
575 TempSet.insert(ListOfLocalAtoms[mol->FindAtom(*sprinter)->getNr()]->getNr());
576 // insert into FragmentList
577 FragmentList.insert(GraphPair(TempSet, pair<int, double> (KeySetCounter++, (*runner).second.second)));
578 }
579 TempSet.clear();
580 }
581 } else
582 ELOG(2, "KeySetList is NULL or empty.");
583
584 if (FreeList) {
585 // free the index lookup list
586 ListOfLocalAtoms.clear();
587 }
588 return status;
589}
590
591/** Translate list into global numbers (i.e. ones that are valid in "this" molecule, not in MolecularWalker->Leaf)
592 * \param &FragmentList Graph with local numbers per fragment
593 * \param &TotalNumberOfKeySets global key set counter
594 * \param &TotalGraph Graph to be filled with global numbers
595 */
596void Fragmentation::TranslateIndicesToGlobalIDs(Graph &FragmentList, int &TotalNumberOfKeySets, Graph &TotalGraph)
597{
598// Info FunctionInfo(__func__);
599 for (Graph::iterator runner = FragmentList.begin(); runner != FragmentList.end(); runner++) {
600 KeySet TempSet;
601 for (KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++)
602 TempSet.insert((mol->FindAtom(*sprinter))->GetTrueFather()->getId());
603 TotalGraph.insert(GraphPair(TempSet, pair<int, double> (TotalNumberOfKeySets++, (*runner).second.second)));
604 }
605}
606
Note: See TracBrowser for help on using the repository browser.