source: src/Fragmentation/Fragmentation.cpp@ d760bb

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

REFACTOR: PowerSetGenerator now creates all orders in one go, not limited to one by one increments.

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