Changeset 5034e1 for src/molecule_fragmentation.cpp

Timestamp:

Oct 12, 2009, 12:10:43 PM (16 years ago)

Author:

Frederik Heber <heber@…>

Branches:

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

Children:

d2943b

Parents:

4a7776a

git-author:

Frederik Heber <heber@…> (10/12/09 12:04:44)

git-committer:

Frederik Heber <heber@…> (10/12/09 12:10:43)

Message:

First half of molecule_fragmentation.cpp refactoring.

• new template helper function: IncrementalAbsoluteValue() which also increments itself (fixed by one)
• new template molecule function: SetAtomValueToIndexedArray() which sets an attribute of atom class according to some array
• Store/ParseKeySetFile() removed from class molecule, as there was no need
• molecule::StoreFragmentFromKeySet() split up into: StoreFragmentFromKeySet_Init() and CreateInducedSubgraphOfFragment()
• molecule::CreateMappingLabelsToConfigSequence() uses SetIndexedArrayForEachAtomTo() with new IncrementalAbsoluteValue()
• molecule::ParseOrderAtSiteFromFile() uses SetAtomValueToIndexedArray() on MaxOrder and Adaptive Order
• atom has new member function: OutputOrder()
• molecule::StoreOrderAtSiteFile() uses ActOnAllAtoms with above
• molecule::CheckOrderAtSite() rewritten, new helper functions are: GraphToIndexedKeySet(), ScanAdaptiveFileIntoMap() (therein InsertIntoAdaptiveCriteriaList()), ReMapAdaptiveCriteriaListToValue(), MarkUpdateCandidates(), PrintAtomMask()
• StoreKeySetFile split up, TEFactors are now stored in StoreTEFactorsFile, same goes with Parse...()

File:

• src/molecule_fragmentation.cpp (modified) (22 diffs)

src/molecule_fragmentation.cpp

@@ -50 +50 @@
  * \return true - at least one valid atom id parsed, false - CurrentSet is empty
  */
-bool molecule::ScanBufferIntoKeySet(ofstream *out, char *buffer, KeySet &CurrentSet)
+bool ScanBufferIntoKeySet(ofstream *out, char *buffer, KeySet &CurrentSet)
 {
   stringstream line;
@@ -59 +59 @@
   while (!line.eof()) {
     line >> AtomNr;
-    if ((AtomNr >= 0) && (AtomNr < AtomCount)) {
+    if (AtomNr >= 0) {
       CurrentSet.insert(AtomNr);  // insert at end, hence in same order as in file!
       status++;
@@ -82 +82 @@
  * \return true - parsing successfully, false - failure on parsing (FragmentList will be NULL)
  */
-bool molecule::ParseKeySetFile(ofstream *out, char *path, Graph *&FragmentList)
+bool ParseKeySetFile(ofstream *out, char *path, Graph *&FragmentList)
 {
   bool status = true;
@@ -89 +89 @@
   GraphTestPair testGraphInsert;
   int NumberOfFragments = 0;
-  double TEFactor;
   char *filename = Malloc<char>(MAXSTRINGSIZE, "molecule::ParseKeySetFile - filename");
 
@@ -124 +123 @@
   }
 
+  return status;
+};
+
+/** Parses the TE factors file and fills \a *FragmentList from the known molecule structure.
+ * -# Scans TEFactors file and sets the TEFactor of each key set in the temporary graph accordingly
+ * \param *out output stream for debugging
+ * \param *path path to file
+ * \param *FragmentList graph whose nodes's TE factors are set on return
+ * \return true - parsing successfully, false - failure on parsing
+ */
+bool ParseTEFactorsFile(ofstream *out, char *path, Graph *FragmentList)
+{
+  bool status = true;
+  ifstream InputFile;
+  stringstream line;
+  GraphTestPair testGraphInsert;
+  int NumberOfFragments = 0;
+  double TEFactor;
+  char *filename = Malloc<char>(MAXSTRINGSIZE, "molecule::ParseTEFactorsFile - filename");
+
+  if (FragmentList == NULL) { // check list pointer
+    FragmentList = new Graph;
+  }
+
   // 2nd pass: open TEFactors file and read
   *out << Verbose(1) << "Parsing the TEFactors file ... " << endl;
@@ -155 +178 @@
 };
 
-/** Stores keysets and TEFactors to file.
- * \param *out output stream for debugging
- * \param KeySetList Graph with Keysets and factors
+/** Stores key sets to file.
+ * \param *out output stream for debugging
+ * \param KeySetList Graph with Keysets
  * \param *path path to file
  * \return true - file written successfully, false - writing failed
  */
-bool molecule::StoreKeySetFile(ofstream *out, Graph &KeySetList, char *path)
+bool StoreKeySetFile(ofstream *out, Graph &KeySetList, char *path)
 {
   ofstream output;
@@ -190 +213 @@
   output.close();
   output.clear();
+
+  return status;
+};
+
+
+/** Stores TEFactors to file.
+ * \param *out output stream for debugging
+ * \param KeySetList Graph with factors
+ * \param *path path to file
+ * \return true - file written successfully, false - writing failed
+ */
+bool StoreTEFactorsFile(ofstream *out, Graph &KeySetList, char *path)
+{
+  ofstream output;
+  bool status =  true;
+  string line;
 
   // open TEFactors file
@@ -211 +250 @@
 };
 
+/** For a given graph, sorts KeySets into a (index, keyset) map.
+ * \param *GlobalKeySetList list of keysets with global ids (valid in "this" molecule) needed for adaptive increase
+ * \return map from index to keyset
+ */
+map<int,KeySet> * GraphToIndexedKeySet(Graph *GlobalKeySetList)
+{
+  map<int,KeySet> *IndexKeySetList = new map<int,KeySet>;
+  for(Graph::iterator runner = GlobalKeySetList->begin(); runner != GlobalKeySetList->end(); runner++) {
+    IndexKeySetList->insert( pair<int,KeySet>(runner->second.first,runner->first) );
+  }
+  return IndexKeySetList;
+};
+
+/** Inserts a (\a No, \a value) pair into the list, overwriting present one.
+ * Note if values are equal, No will decided on which is first
+ * \param *out output stream for debugging
+ * \param &AdaptiveCriteriaList list to insert into
+ * \param &IndexedKeySetList list to find key set for a given index \a No
+ * \param FragOrder current bond order of fragment
+ * \param No index of keyset
+ * \param value energy value
+ */
+void InsertIntoAdaptiveCriteriaList(ofstream *out, map<int, pair<double,int> > *AdaptiveCriteriaList, map<int,KeySet> &IndexKeySetList, int FragOrder, int No, double Value)
+{
+  map<int,KeySet>::iterator marker = IndexKeySetList.find(No);    // find keyset to Frag No.
+  if (marker != IndexKeySetList.end()) {  // if found
+    Value *= 1 + MYEPSILON*(*((*marker).second.begin()));     // in case of equal energies this makes them not equal without changing anything actually
+    // 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
+    pair <map<int, pair<double,int> >::iterator, bool> InsertedElement = AdaptiveCriteriaList->insert( make_pair(*((*marker).second.begin()), pair<double,int>( fabs(Value), FragOrder) ));
+    map<int, pair<double,int> >::iterator PresentItem = InsertedElement.first;
+    if (!InsertedElement.second) { // this root is already present
+      if ((*PresentItem).second.second < FragOrder)  // if order there is lower, update entry with higher-order term
+        //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)
+        {  // if value is smaller, update value and order
+        (*PresentItem).second.first = fabs(Value);
+        (*PresentItem).second.second = FragOrder;
+        *out << Verbose(2) << "Updated element (" <<  (*PresentItem).first << ",[" << (*PresentItem).second.first << "," << (*PresentItem).second.second << "])." << endl;
+      } else {
+        *out << Verbose(2) << "Did not update element " <<  (*PresentItem).first << " as " << FragOrder << " is less than or equal to " << (*PresentItem).second.second << "." << endl;
+      }
+    } else {
+      *out << Verbose(2) << "Inserted element (" <<  (*PresentItem).first << ",[" << (*PresentItem).second.first << "," << (*PresentItem).second.second << "])." << endl;
+    }
+  } else {
+    *out << Verbose(1) << "No Fragment under No. " << No << "found." << endl;
+  }
+};
+
+/** Scans the adaptive order file and insert (index, value) into map.
+ * \param *out output stream for debugging
+ * \param *path path to ENERGYPERFRAGMENT file (may be NULL if Order is non-negative)
+ * \param &IndexedKeySetList list to find key set for a given index \a No
+ * \return adaptive criteria list from file
+ */
+map<int, pair<double,int> > * ScanAdaptiveFileIntoMap(ofstream *out, char *path, map<int,KeySet> &IndexKeySetList)
+{
+  map<int, pair<double,int> > *AdaptiveCriteriaList = new map<int, pair<double,int> >;
+  int No = 0, FragOrder = 0;
+  double Value = 0.;
+  char *buffer = Malloc<char>(MAXSTRINGSIZE, "molecule::CheckOrderAtSite: *buffer");
+  sprintf(buffer, "%s/%s%s.dat", path, FRAGMENTPREFIX, ENERGYPERFRAGMENT);
+  ifstream InputFile(buffer, ios::in);
+
+  if (CountLinesinFile(InputFile) > 0) {
+    // each line represents a fragment root (Atom::nr) id and its energy contribution
+    InputFile.getline(buffer, MAXSTRINGSIZE); // skip comment lines
+    InputFile.getline(buffer, MAXSTRINGSIZE);
+    while(!InputFile.eof()) {
+      InputFile.getline(buffer, MAXSTRINGSIZE);
+      if (strlen(buffer) > 2) {
+        //*out << Verbose(2) << "Scanning: " << buffer << endl;
+        stringstream line(buffer);
+        line >> FragOrder;
+        line >> ws >> No;
+        line >> ws >> Value; // skip time entry
+        line >> ws >> Value;
+        No -= 1;  // indices start at 1 in file, not 0
+        //*out << Verbose(2) << " - yields (" << No << "," << Value << ", " << FragOrder << ")" << endl;
+
+        // clean the list of those entries that have been superceded by higher order terms already
+        InsertIntoAdaptiveCriteriaList(out, AdaptiveCriteriaList, IndexKeySetList, FragOrder, No, Value);
+      }
+    }
+    // close and done
+    InputFile.close();
+    InputFile.clear();
+  }
+  Free(&buffer);
+
+  return AdaptiveCriteriaList;
+};
+
+/** Maps adaptive criteria list back onto (Value, (Root Nr., Order))
+ * (i.e. sorted by value to pick the highest ones)
+ * \param *out output stream for debugging
+ * \param &AdaptiveCriteriaList list to insert into
+ * \param *mol molecule with atoms
+ * \return remapped list
+ */
+map<double, pair<int,int> >  * ReMapAdaptiveCriteriaListToValue(ofstream *out, map<int, pair<double,int> > *AdaptiveCriteriaList, molecule *mol)
+{
+  atom *Walker = mol->start;
+  map<double, pair<int,int> > *FinalRootCandidates = new map<double, pair<int,int> > ;
+  *out << Verbose(1) << "Root candidate list is: " << endl;
+  for(map<int, pair<double,int> >::iterator runner = AdaptiveCriteriaList->begin(); runner != AdaptiveCriteriaList->end(); runner++) {
+    Walker = mol->FindAtom((*runner).first);
+    if (Walker != NULL) {
+      //if ((*runner).second.second >= Walker->AdaptiveOrder) { // only insert if this is an "active" root site for the current order
+      if (!Walker->MaxOrder) {
+        *out << Verbose(2) << "(" << (*runner).first << ",[" << (*runner).second.first << "," << (*runner).second.second << "])" << endl;
+        FinalRootCandidates->insert( make_pair( (*runner).second.first, pair<int,int>((*runner).first, (*runner).second.second) ) );
+      } else {
+        *out << Verbose(2) << "Excluding (" << *Walker << ", " << (*runner).first << ",[" << (*runner).second.first << "," << (*runner).second.second << "]), as it has reached its maximum order." << endl;
+      }
+    } else {
+      cerr << "Atom No. " << (*runner).second.first << " was not found in this molecule." << endl;
+    }
+  }
+  return FinalRootCandidates;
+};
+
+/** Marks all candidate sites for update if below adaptive threshold.
+ * Picks a given number of highest values and set *AtomMask to true.
+ * \param *out output stream for debugging
+ * \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
+ * \param FinalRootCandidates list candidates to check
+ * \param Order desired order
+ * \param *mol molecule with atoms
+ * \return true - if update is necessary, false - not
+ */
+bool MarkUpdateCandidates(ofstream *out, bool *AtomMask, map<double, pair<int,int> > &FinalRootCandidates, int Order, molecule *mol)
+{
+  atom *Walker = mol->start;
+  int No = -1;
+  bool status = false;
+  for(map<double, pair<int,int> >::iterator runner = FinalRootCandidates.upper_bound(pow(10.,Order)); runner != FinalRootCandidates.end(); runner++) {
+    No = (*runner).second.first;
+    Walker = mol->FindAtom(No);
+    //if (Walker->AdaptiveOrder < MinimumRingSize[Walker->nr]) {
+      *out << Verbose(2) << "Root " << No << " is still above threshold (10^{" << Order <<"}: " << runner->first << ", setting entry " << No << " of Atom mask to true." << endl;
+      AtomMask[No] = true;
+      status = true;
+    //} else
+      //*out << 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;
+  }
+  return status;
+};
+
+/** print atom mask for debugging.
+ * \param *out output stream for debugging
+ * \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
+ * \param AtomCount number of entries in \a *AtomMask
+ */
+void PrintAtomMask(ofstream *out, bool *AtomMask, int AtomCount)
+{
+  *out << "              ";
+  for(int i=0;i<AtomCount;i++)
+    *out << (i % 10);
+  *out << endl << "Atom mask is: ";
+  for(int i=0;i<AtomCount;i++)
+    *out << (AtomMask[i] ? "t" : "f");
+  *out << endl;
+};
 
 /** Checks whether the OrderAtSite is still below \a Order at some site.
@@ -225 +427 @@
   atom *Walker = start;
   bool status = false;
-  ifstream InputFile;
 
   // initialize mask list
@@ -234 +435 @@
     if (AtomMask[AtomCount] == true)  // break after one step
       return false;
+
+    // transmorph graph keyset list into indexed KeySetList
+    if (GlobalKeySetList == NULL) {
+      cout << Verbose(1) << "ERROR: Given global key set list (graph) is NULL!" << endl;
+      return false;
+    }
+    map<int,KeySet> *IndexKeySetList = GraphToIndexedKeySet(GlobalKeySetList);
+
     // parse the EnergyPerFragment file
-    char *buffer = Malloc<char>(MAXSTRINGSIZE, "molecule::CheckOrderAtSite: *buffer");
-    sprintf(buffer, "%s/%s%s.dat", path, FRAGMENTPREFIX, ENERGYPERFRAGMENT);
-    InputFile.open(buffer, ios::in);
-    if ((InputFile != NULL) && (GlobalKeySetList != NULL)) {
-      // transmorph graph keyset list into indexed KeySetList
-      map<int,KeySet> IndexKeySetList;
-      for(Graph::iterator runner = GlobalKeySetList->begin(); runner != GlobalKeySetList->end(); runner++) {
-        IndexKeySetList.insert( pair<int,KeySet>(runner->second.first,runner->first) );
-      }
-      int lines = 0;
-      // count the number of lines, i.e. the number of fragments
-      InputFile.getline(buffer, MAXSTRINGSIZE); // skip comment lines
-      InputFile.getline(buffer, MAXSTRINGSIZE);
-      while(!InputFile.eof()) {
-        InputFile.getline(buffer, MAXSTRINGSIZE);
-        lines++;
-      }
-      //*out << Verbose(2) << "Scanned " << lines-1 << " lines." << endl;   // one endline too much
-      InputFile.clear();
-      InputFile.seekg(ios::beg);
-      map<int, pair<double,int> > AdaptiveCriteriaList;  // (Root No., (Value, Order)) !
-      int No, FragOrder;
-      double Value;
-      // each line represents a fragment root (Atom::nr) id and its energy contribution
-      InputFile.getline(buffer, MAXSTRINGSIZE); // skip comment lines
-      InputFile.getline(buffer, MAXSTRINGSIZE);
-      while(!InputFile.eof()) {
-        InputFile.getline(buffer, MAXSTRINGSIZE);
-        if (strlen(buffer) > 2) {
-          //*out << Verbose(2) << "Scanning: " << buffer << endl;
-          stringstream line(buffer);
-          line >> FragOrder;
-          line >> ws >> No;
-          line >> ws >> Value; // skip time entry
-          line >> ws >> Value;
-          No -= 1;  // indices start at 1 in file, not 0
-          //*out << Verbose(2) << " - yields (" << No << "," << Value << ", " << FragOrder << ")" << endl;
-
-          // clean the list of those entries that have been superceded by higher order terms already
-          map<int,KeySet>::iterator marker = IndexKeySetList.find(No);    // find keyset to Frag No.
-          if (marker != IndexKeySetList.end()) {  // if found
-            Value *= 1 + MYEPSILON*(*((*marker).second.begin()));     // in case of equal energies this makes em not equal without changing anything actually
-            // 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
-            pair <map<int, pair<double,int> >::iterator, bool> InsertedElement = AdaptiveCriteriaList.insert( make_pair(*((*marker).second.begin()), pair<double,int>( fabs(Value), FragOrder) ));
-            map<int, pair<double,int> >::iterator PresentItem = InsertedElement.first;
-            if (!InsertedElement.second) { // this root is already present
-              if ((*PresentItem).second.second < FragOrder)  // if order there is lower, update entry with higher-order term
-                //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)
-                {  // if value is smaller, update value and order
-                (*PresentItem).second.first = fabs(Value);
-                (*PresentItem).second.second = FragOrder;
-                *out << Verbose(2) << "Updated element (" <<  (*PresentItem).first << ",[" << (*PresentItem).second.first << "," << (*PresentItem).second.second << "])." << endl;
-              } else {
-                *out << Verbose(2) << "Did not update element " <<  (*PresentItem).first << " as " << FragOrder << " is less than or equal to " << (*PresentItem).second.second << "." << endl;
-              }
-            } else {
-              *out << Verbose(2) << "Inserted element (" <<  (*PresentItem).first << ",[" << (*PresentItem).second.first << "," << (*PresentItem).second.second << "])." << endl;
-            }
-          } else {
-            *out << Verbose(1) << "No Fragment under No. " << No << "found." << endl;
-          }
-        }
-      }
-      // then map back onto (Value, (Root Nr., Order)) (i.e. sorted by value to pick the highest ones)
-      map<double, pair<int,int> > FinalRootCandidates;
-      *out << Verbose(1) << "Root candidate list is: " << endl;
-      for(map<int, pair<double,int> >::iterator runner = AdaptiveCriteriaList.begin(); runner != AdaptiveCriteriaList.end(); runner++) {
-        Walker = FindAtom((*runner).first);
-        if (Walker != NULL) {
-          //if ((*runner).second.second >= Walker->AdaptiveOrder) { // only insert if this is an "active" root site for the current order
-          if (!Walker->MaxOrder) {
-            *out << Verbose(2) << "(" << (*runner).first << ",[" << (*runner).second.first << "," << (*runner).second.second << "])" << endl;
-            FinalRootCandidates.insert( make_pair( (*runner).second.first, pair<int,int>((*runner).first, (*runner).second.second) ) );
-          } else {
-            *out << Verbose(2) << "Excluding (" << *Walker << ", " << (*runner).first << ",[" << (*runner).second.first << "," << (*runner).second.second << "]), as it has reached its maximum order." << endl;
-          }
-        } else {
-          cerr << "Atom No. " << (*runner).second.first << " was not found in this molecule." << endl;
-        }
-      }
-      // pick the ones still below threshold and mark as to be adaptively updated
-      for(map<double, pair<int,int> >::iterator runner = FinalRootCandidates.upper_bound(pow(10.,Order)); runner != FinalRootCandidates.end(); runner++) {
-        No = (*runner).second.first;
-        Walker = FindAtom(No);
-        //if (Walker->AdaptiveOrder < MinimumRingSize[Walker->nr]) {
-          *out << Verbose(2) << "Root " << No << " is still above threshold (10^{" << Order <<"}: " << runner->first << ", setting entry " << No << " of Atom mask to true." << endl;
-          AtomMask[No] = true;
-          status = true;
-        //} else
-          //*out << 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;
-      }
-      // close and done
-      InputFile.close();
-      InputFile.clear();
-    } else {
-      cerr << "Unable to parse " << buffer << " file, incrementing all." << endl;
+    map<int, pair<double,int> > *AdaptiveCriteriaList = ScanAdaptiveFileIntoMap(out, path, *IndexKeySetList); // (Root No., (Value, Order)) !
+    if (AdaptiveCriteriaList->empty()) {
+      cerr << "Unable to parse file, incrementing all." << endl;
       while (Walker->next != end) {
         Walker = Walker->next;
@@ -342 +458 @@
       }
     }
-    Free(&buffer);
-    // pick a given number of highest values and set AtomMask
+    // then map back onto (Value, (Root Nr., Order)) (i.e. sorted by value to pick the highest ones)
+    map<double, pair<int,int> > *FinalRootCandidates = ReMapAdaptiveCriteriaListToValue(out, AdaptiveCriteriaList, this);
+
+    // pick the ones still below threshold and mark as to be adaptively updated
+    MarkUpdateCandidates(out, AtomMask, *FinalRootCandidates, Order, this);
+
+    Free(&IndexKeySetList);
+    Free(&AdaptiveCriteriaList);
+    Free(&FinalRootCandidates);
   } else { // global increase of Bond Order
     while (Walker->next != end) {
@@ -367 +490 @@
   }
 
-  // print atom mask for debugging
-  *out << "              ";
-  for(int i=0;i<AtomCount;i++)
-    *out << (i % 10);
-  *out << endl << "Atom mask is: ";
-  for(int i=0;i<AtomCount;i++)
-    *out << (AtomMask[i] ? "t" : "f");
-  *out << endl;
+  PrintAtomMask(out, AtomMask, AtomCount); // for debugging
 
   return status;
@@ -383 +499 @@
  * \param *&SortIndex Mapping array of size molecule::AtomCount
  * \return true - success, false - failure of SortIndex alloc
- * \todo do we really need this still as the IonType may appear in any order due to recent changes
  */
 bool molecule::CreateMappingLabelsToConfigSequence(ofstream *out, int *&SortIndex)
 {
-  element *runner = elemente->start;
-  int AtomNo = 0;
-  atom *Walker = NULL;
-
   if (SortIndex != NULL) {
     *out << Verbose(1) << "SortIndex is " << SortIndex << " and not NULL as expected." << endl;
     return false;
   }
-  SortIndex = Malloc<int>(AtomCount, "molecule::FragmentMolecule: *SortIndex");
+  SortIndex = Malloc<int>(AtomCount, "molecule::CreateMappingLabelsToConfigSequence: *SortIndex");
   for(int i=AtomCount;i--;)
     SortIndex[i] = -1;
-  while (runner->next != elemente->end) { // go through every element
-    runner = runner->next;
-    if (ElementsInMolecule[runner->Z]) { // if this element got atoms
-      Walker = start;
-      while (Walker->next != end) { // go through every atom of this element
-        Walker = Walker->next;
-        if (Walker->type->Z == runner->Z) // if this atom fits to element
-          SortIndex[Walker->nr] = AtomNo++;
-      }
-    }
-  }
+
+  int AtomNo = 0;
+  SetIndexedArrayForEachAtomTo( SortIndex, &atom::nr, &IncrementalAbsoluteValue, AtomNo );
+
   return true;
 };
@@ -631 +735 @@
   *out << Verbose(1) << "Writing OrderAtSite " << ORDERATSITEFILE << " ... " << endl;
   if (file != NULL) {
-    atom *Walker = start;
-    while (Walker->next != end) {
-      Walker = Walker->next;
-      file << Walker->nr << "\t" << (int)Walker->AdaptiveOrder << "\t" << (int)Walker->MaxOrder << endl;
-      *out << Verbose(2) << "Storing: " << Walker->nr << "\t" << (int)Walker->AdaptiveOrder << "\t" << (int)Walker->MaxOrder << "." << endl;
-    }
+    ActOnAllAtoms( &atom::OutputOrder, &file );
     file.close();
     *out << Verbose(1) << "done." << endl;
@@ -683 +782 @@
       }
     }
-    atom *Walker = start;
-    while (Walker->next != end) { // fill into atom classes
-      Walker = Walker->next;
-      Walker->AdaptiveOrder = OrderArray[Walker->nr];
-      Walker->MaxOrder = MaxArray[Walker->nr];
-      *out << Verbose(2) << *Walker << " gets order " << (int)Walker->AdaptiveOrder << " and is " << (!Walker->MaxOrder ? "not " : " ") << "maxed." << endl;
-    }
     file.close();
+
+    // set atom values
+    SetAtomValueToIndexedArray( OrderArray, &atom::nr, &atom::AdaptiveOrder );
+    SetAtomValueToIndexedArray( MaxArray, &atom::nr, &atom::MaxOrder );
+
     *out << Verbose(1) << "done." << endl;
     status = true;
@@ -732 +829 @@
 };
 
-/** Stores a fragment from \a KeySet into \a molecule.
- * First creates the minimal set of atoms from the KeySet, then creates the bond structure from the complete
- * molecule and adds missing hydrogen where bonds were cut.
- * \param *out output stream for debugging messages
+/** Initializes some value for putting fragment of \a *mol into \a *Leaf.
+ * \param *mol total molecule
+ * \param *Leaf fragment molecule
  * \param &Leaflet pointer to KeySet structure
- * \param IsAngstroem whether we have Ansgtroem or bohrradius
- * \return pointer to constructed molecule
- */
-molecule * molecule::StoreFragmentFromKeySet(ofstream *out, KeySet &Leaflet, bool IsAngstroem)
-{
-  atom *Runner = NULL, *FatherOfRunner = NULL, *OtherFather = NULL;
-  atom **SonList = Malloc<atom*>(AtomCount, "molecule::StoreFragmentFromStack: **SonList");
-  molecule *Leaf = new molecule(elemente);
-  bool LonelyFlag = false;
-  int size;
-
-//  *out << Verbose(1) << "Begin of StoreFragmentFromKeyset." << endl;
-
-  Leaf->BondDistance = BondDistance;
+ * \param **SonList list which atom of \a *Leaf is a son of which atom in \a *mol
+ * \return number of atoms in fragment
+ */
+int StoreFragmentFromKeySet_Init(molecule *mol, molecule *Leaf, KeySet &Leaflet, atom **SonList)
+{
+  atom *FatherOfRunner = NULL;
+
+  Leaf->BondDistance = mol->BondDistance;
   for(int i=NDIM*2;i--;)
-    Leaf->cell_size[i] = cell_size[i];
+    Leaf->cell_size[i] = mol->cell_size[i];
 
   // initialise SonList (indicates when we need to replace a bond with hydrogen instead)
-  for(int i=AtomCount;i--;)
+  for(int i=mol->AtomCount;i--;)
     SonList[i] = NULL;
 
   // first create the minimal set of atoms from the KeySet
-  size = 0;
+  int size = 0;
   for(KeySet::iterator runner = Leaflet.begin(); runner != Leaflet.end(); runner++) {
-    FatherOfRunner = FindAtom((*runner));  // find the id
+    FatherOfRunner = mol->FindAtom((*runner));  // find the id
     SonList[FatherOfRunner->nr] = Leaf->AddCopyAtom(FatherOfRunner);
     size++;
   }
-
-  // create the bonds between all: Make it an induced subgraph and add hydrogen
-//  *out << Verbose(2) << "Creating bonds from father graph (i.e. induced subgraph creation)." << endl;
-  Runner = Leaf->start;
+  return size;
+};
+
+/** Creates an induced subgraph out of a fragmental key set, adding bonds and hydrogens (if treated specially).
+ * \param *out output stream for debugging messages
+ * \param *mol total molecule
+ * \param *Leaf fragment molecule
+ * \param IsAngstroem whether we have Ansgtroem or bohrradius
+ * \param **SonList list which atom of \a *Leaf is a son of which atom in \a *mol
+ */
+void CreateInducedSubgraphOfFragment(ofstream *out, molecule *mol, molecule *Leaf, atom **SonList, bool IsAngstroem)
+{
+  bool LonelyFlag = false;
+  atom *OtherFather = NULL;
+  atom *FatherOfRunner = NULL;
+  Leaf->CountAtoms(out);
+
+  atom *Runner = Leaf->start;
   while (Runner->next != Leaf->end) {
     Runner = Runner->next;
@@ -775 +879 @@
     if (SonList[FatherOfRunner->nr] != NULL)  {  // check if this, our father, is present in list
       // create all bonds
-      for (int i=0;i<NumberOfBondsPerAtom[FatherOfRunner->nr];i++) { // go through every bond of father
-        OtherFather = ListOfBondsPerAtom[FatherOfRunner->nr][i]->GetOtherAtom(FatherOfRunner);
+      for (int i=0;i<mol->NumberOfBondsPerAtom[FatherOfRunner->nr];i++) { // go through every bond of father
+        OtherFather = mol->ListOfBondsPerAtom[FatherOfRunner->nr][i]->GetOtherAtom(FatherOfRunner);
 //        *out << Verbose(2) << "Father " << *FatherOfRunner << " of son " << *SonList[FatherOfRunner->nr] << " is bound to " << *OtherFather;
         if (SonList[OtherFather->nr] != NULL) {
@@ -783 +887 @@
 //            *out << Verbose(3) << "Adding Bond: ";
 //            *out <<
-            Leaf->AddBond(Runner, SonList[OtherFather->nr], ListOfBondsPerAtom[FatherOfRunner->nr][i]->BondDegree);
+            Leaf->AddBond(Runner, SonList[OtherFather->nr], mol->ListOfBondsPerAtom[FatherOfRunner->nr][i]->BondDegree);
 //            *out << "." << endl;
             //NumBonds[Runner->nr]++;
@@ -794 +898 @@
 #ifdef ADDHYDROGEN
           //*out << Verbose(3) << "Adding Hydrogen to " << Runner->Name << " and a bond in between." << endl;
-          if(!Leaf->AddHydrogenReplacementAtom(out, ListOfBondsPerAtom[FatherOfRunner->nr][i], Runner, FatherOfRunner, OtherFather, ListOfBondsPerAtom[FatherOfRunner->nr],NumberOfBondsPerAtom[FatherOfRunner->nr], IsAngstroem))
+          if(!Leaf->AddHydrogenReplacementAtom(out, mol->ListOfBondsPerAtom[FatherOfRunner->nr][i], Runner, FatherOfRunner, OtherFather, mol->ListOfBondsPerAtom[FatherOfRunner->nr],mol->NumberOfBondsPerAtom[FatherOfRunner->nr], IsAngstroem))
             exit(1);
 #endif
@@ -803 +907 @@
       *out << Verbose(0) << "ERROR: Son " << Runner->Name << " has father " << FatherOfRunner->Name << " but its entry in SonList is " << SonList[FatherOfRunner->nr] << "!" << endl;
     }
-    if ((LonelyFlag) && (size > 1)) {
+    if ((LonelyFlag) && (Leaf->AtomCount > 1)) {
       *out << Verbose(0) << *Runner << "has got bonds only to hydrogens!" << endl;
     }
@@ -811 +915 @@
 #endif
   }
+};
+
+/** Stores a fragment from \a KeySet into \a molecule.
+ * First creates the minimal set of atoms from the KeySet, then creates the bond structure from the complete
+ * molecule and adds missing hydrogen where bonds were cut.
+ * \param *out output stream for debugging messages
+ * \param &Leaflet pointer to KeySet structure
+ * \param IsAngstroem whether we have Ansgtroem or bohrradius
+ * \return pointer to constructed molecule
+ */
+molecule * molecule::StoreFragmentFromKeySet(ofstream *out, KeySet &Leaflet, bool IsAngstroem)
+{
+  atom **SonList = Malloc<atom*>(AtomCount, "molecule::StoreFragmentFromStack: **SonList");
+  molecule *Leaf = new molecule(elemente);
+
+//  *out << Verbose(1) << "Begin of StoreFragmentFromKeyset." << endl;
+  StoreFragmentFromKeySet_Init(this, Leaf, Leaflet, SonList);
+  // create the bonds between all: Make it an induced subgraph and add hydrogen
+//  *out << Verbose(2) << "Creating bonds from father graph (i.e. induced subgraph creation)." << endl;
+  CreateInducedSubgraphOfFragment(out, this, Leaf, SonList, IsAngstroem);
+
   Leaf->CreateListOfBondsPerAtom(out);
   //Leaflet->Leaf->ScanForPeriodicCorrection(out);
@@ -817 +942 @@
   return Leaf;
 };
-
-/** Creates \a MoleculeListClass of all unique fragments of the \a molecule containing \a Order atoms or vertices.
- * The picture to have in mind is that of a DFS "snake" of a certain length \a Order, i.e. as in the infamous
- * computer game, that winds through the connected graph representing the molecule. Color (white,
- * lightgray, darkgray, black) indicates whether a vertex has been discovered so far or not. Labels will help in
- * creating only unique fragments and not additional ones with vertices simply in different sequence.
- * The Predecessor is always the one that came before in discovering, needed on backstepping. And
- * finally, the ShortestPath is needed for removing vertices from the snake stack during the back-
- * stepping.
- * \param *out output stream for debugging
- * \param Order number of atoms in each fragment
- * \param *configuration configuration for writing config files for each fragment
- * \return List of all unique fragments with \a Order atoms
- */
-/*
-MoleculeListClass * molecule::CreateListOfUniqueFragmentsOfOrder(ofstream *out, int Order, config *configuration)
-{
-  atom **PredecessorList = Malloc<atom*>(AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: **PredecessorList");
-  int *ShortestPathList = Malloc<int>(AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *ShortestPathList");
-  int *Labels = Malloc<int>(AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *Labels");
-  enum Shading *ColorVertexList = Malloc<enum Shading>(AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *ColorList");
-  enum Shading *ColorEdgeList = Malloc<enum Shading>(BondCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *ColorBondList");
-  StackClass<atom *> *RootStack = new StackClass<atom *>(AtomCount);
-  StackClass<atom *> *TouchedStack = new StackClass<atom *>((int)pow(4,Order)+2); // number of atoms reached from one with maximal 4 bonds plus Root itself
-  StackClass<atom *> *SnakeStack = new StackClass<atom *>(Order+1); // equal to Order is not possible, as then the StackClass<atom *> cannot discern between full and empty stack!
-  MoleculeLeafClass *Leaflet = NULL, *TempLeaf = NULL;
-  MoleculeListClass *FragmentList = NULL;
-  atom *Walker = NULL, *OtherAtom = NULL, *Root = NULL, *Removal = NULL;
-  bond *Binder = NULL;
-  int RunningIndex = 0, FragmentCounter = 0;
-
-  *out << Verbose(1) << "Begin of CreateListOfUniqueFragmentsOfOrder." << endl;
-
-  // reset parent list
-  *out << Verbose(3) << "Resetting labels, parent, predecessor, color and shortest path lists." << endl;
-  for (int i=0;i<AtomCount;i++) { // reset all atom labels
-    // initialise each vertex as white with no predecessor, empty queue, color lightgray, not labelled, no sons
-    Labels[i] = -1;
-    SonList[i] = NULL;
-    PredecessorList[i] = NULL;
-    ColorVertexList[i] = white;
-    ShortestPathList[i] = -1;
-  }
-  for (int i=0;i<BondCount;i++)
-    ColorEdgeList[i] = white;
-  RootStack->ClearStack();  // clearstack and push first atom if exists
-  TouchedStack->ClearStack();
-  Walker = start->next;
-  while ((Walker != end)
-#ifdef ADDHYDROGEN
-   && (Walker->type->Z == 1)
-        }
-      }
-      *out << ", SP of " << ShortestPathList[Walker->nr]  << " and its color is " << GetColor(ColorVertexList[Walker->nr]) << "." << endl;
-
-      // then check the stack for a newly stumbled upon fragment
-      if (SnakeStack->ItemCount() == Order) { // is stack full?
-        // store the fragment if it is one and get a removal candidate
-        Removal = StoreFragmentFromStack(out, Root, Walker, Leaflet, SnakeStack, ShortestPathList, SonList, Labels, &FragmentCounter, configuration);
-        // remove the candidate if one was found
-        if (Removal != NULL) {
-          *out << Verbose(2) << "Removing item " << Removal->Name << " with SP of " << ShortestPathList[Removal->nr] << " from snake stack." << endl;
-          SnakeStack->RemoveItem(Removal);
-          ColorVertexList[Removal->nr] = lightgray; // return back to not on snake stack but explored marking
-          if (Walker == Removal) { // if the current atom is to be removed, we also have to take a step back
-            Walker = PredecessorList[Removal->nr];
-            *out << Verbose(2) << "Stepping back to " << Walker->Name << "." << endl;
-          }
-        }
-      } else
-        Removal = NULL;
-
-      // finally, look for a white neighbour as the next Walker
-      Binder = NULL;
-      if ((Removal == NULL) || (Walker != PredecessorList[Removal->nr])) {  // don't look, if a new walker has been set above
-        *out << Verbose(2) << "Snake has currently " << SnakeStack->ItemCount() << " item(s)." << endl;
-        OtherAtom = NULL; // this is actually not needed, every atom has at least one neighbour
-        if (ShortestPathList[Walker->nr] < Order) {
-          for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
-            Binder = ListOfBondsPerAtom[Walker->nr][i];
-            *out << Verbose(2) << "Current bond is " << *Binder << ": ";
-            OtherAtom = Binder->GetOtherAtom(Walker);
-            if ((Labels[OtherAtom->nr] != -1) && (Labels[OtherAtom->nr] < Labels[Root->nr])) { // we don't step up to labels bigger than us
-              *out << "Label " << Labels[OtherAtom->nr] << " is smaller than Root's " << Labels[Root->nr] << "." << endl;
-              //ColorVertexList[OtherAtom->nr] = lightgray;    // mark as explored
-            } else { // otherwise check its colour and element
-              if (
-              (OtherAtom->type->Z != 1) &&
-#endif
-                    (ColorEdgeList[Binder->nr] == white)) {  // skip hydrogen, look for unexplored vertices
-                *out << "Moving along " << GetColor(ColorEdgeList[Binder->nr]) << " bond " << Binder << " to " << ((ColorVertexList[OtherAtom->nr] == white) ? "unexplored" : "explored") << " item: " << OtherAtom->Name << "." << endl;
-                // i find it currently rather sensible to always set the predecessor in order to find one's way back
-                //if (PredecessorList[OtherAtom->nr] == NULL) {
-                PredecessorList[OtherAtom->nr] = Walker;
-                *out << Verbose(3) << "Setting Predecessor of " << OtherAtom->Name << " to " << PredecessorList[OtherAtom->nr]->Name << "." << endl;
-                //} else {
-                //  *out << Verbose(3) << "Predecessor of " << OtherAtom->Name << " is " << PredecessorList[OtherAtom->nr]->Name << "." << endl;
-                //}
-                Walker = OtherAtom;
-                break;
-              } else {
-                if (OtherAtom->type->Z == 1)
-                  *out << "Links to a hydrogen atom." << endl;
-                else
-                  *out << "Bond has not white but " << GetColor(ColorEdgeList[Binder->nr]) << " color." << endl;
-              }
-            }
-          }
-        } else {  // means we have stepped beyond the horizon: Return!
-          Walker = PredecessorList[Walker->nr];
-          OtherAtom = Walker;
-          *out << Verbose(3) << "We have gone too far, stepping back to " << Walker->Name << "." << endl;
-        }
-        if (Walker != OtherAtom) {  // if no white neighbours anymore, color it black
-          *out << Verbose(2) << "Coloring " << Walker->Name << " black." << endl;
-          ColorVertexList[Walker->nr] = black;
-          Walker = PredecessorList[Walker->nr];
-        }
-      }
-    } while ((Walker != Root) || (ColorVertexList[Root->nr] != black));
-    *out << Verbose(2) << "Inner Looping is finished." << endl;
-
-    // if we reset all AtomCount atoms, we have again technically O(N^2) ...
-    *out << Verbose(2) << "Resetting lists." << endl;
-    Walker = NULL;
-    Binder = NULL;
-    while (!TouchedStack->IsEmpty()) {
-      Walker = TouchedStack->PopLast();
-      *out << Verbose(3) << "Re-initialising entries of " << *Walker << "." << endl;
-      for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++)
-        ColorEdgeList[ListOfBondsPerAtom[Walker->nr][i]->nr] = white;
-      PredecessorList[Walker->nr] = NULL;
-      ColorVertexList[Walker->nr] = white;
-      ShortestPathList[Walker->nr] = -1;
-    }
-  }
-  *out << Verbose(1) << "Outer Looping over all vertices is done." << endl;
-
-  // copy together
-  *out << Verbose(1) << "Copying all fragments into MoleculeList structure." << endl;
-  FragmentList = new MoleculeListClass(FragmentCounter, AtomCount);
-  RunningIndex = 0;
-  while ((Leaflet != NULL) && (RunningIndex < FragmentCounter))  {
-    FragmentList->ListOfMolecules[RunningIndex++] = Leaflet->Leaf;
-    Leaflet->Leaf = NULL; // prevent molecule from being removed
-    TempLeaf = Leaflet;
-    Leaflet = Leaflet->previous;
-    delete(TempLeaf);
-  };
-
-  // free memory and exit
-  Free(&PredecessorList);
-  Free(&ShortestPathList);
-  Free(&Labels);
-  Free(&ColorVertexList);
-  delete(RootStack);
-  delete(TouchedStack);
-  delete(SnakeStack);
-
-  *out << Verbose(1) << "End of CreateListOfUniqueFragmentsOfOrder." << endl;
-  return FragmentList;
-};
-*/
 
 /** From a given set of Bond sorted by Shortest Path distance, create all possible fragments of size \a SetDimension.