Changeset 644ba1 for molecuilder/src/molecules.cpp

Timestamp:

Aug 18, 2008, 8:57:58 AM (17 years ago)

Author:

Frederik Heber <heber@…>

Children:

04c868

Parents:

e6971b

Message:

Adaptivity fixes, MD by VerletForceIntegration introduced, MD molecule::Trajectories, atom Max::Order, no more recursive going down the fragmentation level

MD
==
molecule::Trajectories is now a map to a struct trajectory list of all the MD steps.
struct Trajectory contains STL vectors of coordinates, velocities and forces. Both are needed for the new VerletForceIntegration.
Parsing of coordinates, velocities and forces from the config file was completely rewritten:

• in the FastParsing case, we just scan for IonType1_1 to find the last step, set the file pointer there and scan all remaining ones
• in the other case, we create the atoms in the first step and put them in a hash for lookup on later steps and read in sequentially (with file pointer moving on).
• This is a lot faster than the old variant.

VerletForceIntegration() implemented in a working manner with force smoothing (mean of actual and last one).
OutputTrajectoriesXYZ() now concatenates the single MD steps into one xyz file, so that the animation can be viewed with e.g. jmol or vmd
config:Deltat is now public (lazy me) and set to 1 instead of 0 initially, also Deltat is parsed accordingly (if not present, defaults to 1)
MatrixContainer::ParseMatrix from parser.cpp is used during VerletForceIntegration() to parse the forces file. Consequently, we have included parser.* in the Makefile.am.
Fix: MoleculeListClass::OutputConfigForListOfFragments() stores config file under config::configpath, before it backup'd the path twice to PathBackup.

Adaptivity
==========
Adaptivity (CheckOrderAtSite()) had compared not absolute, but real values which caused sign problems and faulty behaviour.
Adapatvity (CheckOrderAtSite()) had included atoms by Order (desired one) not by FragOrder (current one) into the list of candidates, which caused faulty behaviour.
CheckOrderAtSite() did single stepping wrong as the mask bit (last in AtomMask) was checked for true not false! Also bit was not set to false initially in FragmentMolecule().
Adaptivity: FragmentMolecule now returns 1 - continue, 2 - don't ... to tell whether we still have to continue with the adaptive cycle (is given as return value from molecuilder)
introduced atom::MaxOrder
StoreOrderAtSiteFile() now also stores the MaxOrder and ParseOrderAtSiteFromFile() parses it back into atom class

Removed Fragmentation Recursion
===============================
As we switched in the prelude of the adaptivity to a monotonous increase from order 1 to the desired one, we don't need to recursively go down each level from a given current bond order, as all these fragments have already been created on the lower orders. Consequently, FragmentBOSSANOVA does not contain NumLevels or alike anymore. This simplified the code a bit, but probably is not yet completely done. (but is working, checked on heptan).
SPFragmentGenerator() does not need Labels anymore, global ones are used for checks. Consequently, PowerSetGenerator() and FragmentSearch structure don't initialise/contain them anymore. We always compare against ...->GetTrueFather()->nr.

File:

• molecuilder/src/molecules.cpp (modified) (36 diffs)

molecuilder/src/molecules.cpp

@@ -455 +455 @@
   getline(xyzfile,line,'\n'); // Read comment
   cout << Verbose(1) << "Comment: " << line << endl;
-
+  
+  if (MDSteps == 0) // no atoms yet present
+    MDSteps++;
   for(i=0;i<NumberOfAtoms;i++){
     Walker = new atom;
@@ -471 +473 @@
       Walker->type = elemente->FindElement(1);
     }
-    for(j=NDIM;j--;)
+    if (Trajectories[Walker].R.size() <= (unsigned int)MDSteps) {
+      Trajectories[Walker].R.resize(MDSteps+10);
+      Trajectories[Walker].U.resize(MDSteps+10);
+      Trajectories[Walker].F.resize(MDSteps+10);
+    }
+    for(j=NDIM;j--;) {
       Walker->x.x[j] = x[j];
-     AddAtom(Walker);  // add to molecule
+      Trajectories[Walker].R.at(MDSteps-1).x[j] = x[j];
+      Trajectories[Walker].U.at(MDSteps-1).x[j] = 0;
+      Trajectories[Walker].F.at(MDSteps-1).x[j] = 0;
+    }
+    AddAtom(Walker);  // add to molecule
     delete(item);
   }
@@ -965 +976 @@
 
 /** Parses nuclear forces from file and performs Verlet integration.
+ * Note that we assume the parsed forces to be in atomic units (hence, if coordinates are in angstroem, we
+ * have to transform them).
  * This adds a new MD step to the config file.
  * \param *file filename
  * \param delta_t time step width in atomic units
+ * \param IsAngstroem whether coordinates are in angstroem (true) or bohrradius (false)
  * \return true - file found and parsed, false - file not found or imparsable
  */
-bool molecule::VerletForceIntegration(char *file, double delta_t)
-{
-  bool status = true;
+bool molecule::VerletForceIntegration(char *file, double delta_t, bool IsAngstroem)
+{
   element *runner = elemente->start;
   atom *walker = NULL;
-  int ElementNo, AtomNo;
+  int AtomNo;
   ifstream input(file);
-  double Forcesvector[3];
-  double dummy;
   string token;
-  size_t oldpos, newpos;
   stringstream item;
-  int i, Ion[2];
-  double a, IonMass;
-  unsigned int size;
+  double a, IonMass, Vector[NDIM];
+  ForceMatrix Force;
 
   CountElements();  // make sure ElementsInMolecule is up to date
@@ -992 +1001 @@
     return false;
   } else {
-    ElementNo = 0;
+    // parse file into ForceMatrix
+    if (!Force.ParseMatrix(file, 0,0,0)) {
+      cerr << "Could not parse Force Matrix file " << file << "." << endl;
+      return false;
+    }
+    if (Force.RowCounter[0] != AtomCount) {
+      cerr << "Mismatch between number of atoms in file " << Force.RowCounter[0] << " and in molecule " << AtomCount << "." << endl;
+      return false;
+    }
+    // correct Forces
+    for(int d=0;d<NDIM;d++)
+      Vector[d] = 0.;
+    for(int i=0;i<AtomCount;i++)
+      for(int d=0;d<NDIM;d++) {
+        Vector[d] += Force.Matrix[0][i][d+5];
+      }
+    for(int i=0;i<AtomCount;i++)
+      for(int d=0;d<NDIM;d++) {
+        Force.Matrix[0][i][d+5] -= Vector[d]/(double)AtomCount;
+      }
+    // and perform Verlet integration for each atom with position, velocity and force vector
+    runner = elemente->start;
     while (runner->next != elemente->end) { // go through every element
       runner = runner->next;
       IonMass = runner->mass;
+      a = delta_t*0.5/IonMass;        // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
       if (ElementsInMolecule[runner->Z]) { // if this element got atoms
-        ElementNo++;
         AtomNo = 0;
         walker = start;
@@ -1003 +1033 @@
           walker = walker->next;
           if (walker->type == runner) { // if this atom fits to element
-            AtomNo++;
-            // parse Forcevector for this ion
-            getline (input, token, '\n');
-            newpos = oldpos = i = 0;
-            while ((newpos = token.find( '\t', oldpos)) != oldpos) {
-              i++;
-              item.str();
-              switch (i) {
-                case 1:
-                case 2:
-                  item >> Ion[i-1];
-                  break;
-                case 6:
-                case 7:
-                case 8:
-                  item >> Forcesvector[i-6];
-                  break;
-                case 3:
-                  if ((Ion[0] != ElementNo) || (Ion[1] != AtomNo))
-                    cout << "ERROR: Expected " << ElementNo << ", " << AtomNo << "but parsed " << Ion[0] << ", " << Ion[1] << "." << endl;
-                  // still parse into dummy! Hence no break here ...
-                default:
-                  item >> dummy;
-              }
-              oldpos = newpos;
+            // check size of vectors
+            if (Trajectories[walker].R.size() <= (unsigned int)(MDSteps)) {
+              //cout << "Increasing size for trajectory array of " << *walker << " to " << (size+10) << "." << endl;
+              Trajectories[walker].R.resize(MDSteps+10);
+              Trajectories[walker].U.resize(MDSteps+10);
+              Trajectories[walker].F.resize(MDSteps+10);
             }
             
-            // perform Verlet integration for this atom with position, velocity and force vector
-            
-            // UpdateR
-//            for (int d=0; d<NDIM; d++) {
-//              R_old_old[d] = R_old[d];        // shift old values
-//              R_old[d] = R[d];
-//              R[d] += delta_t*(U[d]+a*FIon[d]);     // F = m * a and s = 0.5 * a * t^2
-//              FIon_old[d] = FIon[d];          // store old force
-//            }
-//            // UpdateU
-//            a = delta_t*0.5/IonMass;        // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
-//            for (int d=0; d<NDIM; d++) {
-//                U[d] += a * (FIon[d]+FIon_old[d]);
-//            }
-            
-            
-            // check size of vectors
-            size = Trajectories[walker].R.size(); 
-            if (size <= (unsigned int)(MDSteps)) {
-              cout << "Increasing size for trajectory array of " << *walker << " to " << (size+10) << "." << endl;
-              Trajectories[walker].R.resize(size+10);
-              Trajectories[walker].U.resize(size+10);
-              Trajectories[walker].F.resize(size+10);
+            // Update R (and F)
+            for (int d=0; d<NDIM; d++) {
+              Trajectories[walker].F.at(MDSteps).x[d] = -Force.Matrix[0][AtomNo][d+5]*(IsAngstroem ? AtomicLengthToAngstroem : 1.);
+              Trajectories[walker].R.at(MDSteps).x[d] = Trajectories[walker].R.at(MDSteps-1).x[d];
+              Trajectories[walker].R.at(MDSteps).x[d] += delta_t*(Trajectories[walker].U.at(MDSteps-1).x[d]);
+              Trajectories[walker].R.at(MDSteps).x[d] += delta_t*a*(Trajectories[walker].F.at(MDSteps).x[d]);     // F = m * a and s = 0.5 * F/m * t^2 = F * a * t
             }
-            // add trajectory point
-
+            // Update U
+            for (int d=0; d<NDIM; d++) {
+              Trajectories[walker].U.at(MDSteps).x[d] = Trajectories[walker].U.at(MDSteps-1).x[d]; 
+              Trajectories[walker].U.at(MDSteps).x[d] += a * (Trajectories[walker].F.at(MDSteps).x[d]+Trajectories[walker].F.at(MDSteps-1).x[d]);
+            }
+//            cout << "Integrated position&velocity of step " << (MDSteps) << ": ("; 
+//            for (int d=0;d<NDIM;d++)
+//              cout << Trajectories[walker].R.at(MDSteps).x[d] << " ";          // next step
+//            cout << ")\t(";
+//            for (int d=0;d<NDIM;d++)
+//              cout << Trajectories[walker].U.at(MDSteps).x[d] << " ";          // next step
+//            cout << ")" << endl;
+            // next atom
+            AtomNo++;
           }
         }
       }
     }
-    return true;
-  }
+  }
+//  // correct velocities (rather momenta) so that center of mass remains motionless
+//  for(int d=0;d<NDIM;d++)
+//    Vector[d] = 0.;
+//  IonMass = 0.;
+//  walker = start;
+//  while (walker->next != end) { // go through every atom
+//    walker = walker->next;
+//    IonMass += walker->type->mass;  // sum up total mass
+//    for(int d=0;d<NDIM;d++) {
+//      Vector[d] += Trajectories[walker].U.at(MDSteps).x[d]*walker->type->mass;
+//    }
+//  }
+//  walker = start;
+//  while (walker->next != end) { // go through every atom of this element
+//    walker = walker->next;
+//    for(int d=0;d<NDIM;d++) {
+//      Trajectories[walker].U.at(MDSteps).x[d] -= Vector[d]*walker->type->mass/IonMass;
+//    }
+//  }
+  MDSteps++;
+
   
   // exit
-  return status;
+  return true;
 };
 
@@ -1372 +1398 @@
               if (Trajectories[walker].F.at(step).Norm() > MYEPSILON)
                 *out << "\t" << scientific << setprecision(6) << Trajectories[walker].F.at(step).x[0] << "\t" << Trajectories[walker].F.at(step).x[1] << "\t" << Trajectories[walker].F.at(step).x[2] << "\t";
-              *out << " # Number in molecule " << walker->nr << endl;
+              *out << "\t# Number in molecule " << walker->nr << endl;
             }
           }
@@ -1413 +1439 @@
 };
 
+/** Prints molecule with all its trajectories to *out as xyz file.
+ * \param *out output stream
+ */
+bool molecule::OutputTrajectoriesXYZ(ofstream *out)
+{
+  atom *walker = NULL;
+  int No = 0;
+  time_t now;
+  
+  now = time((time_t *)NULL);   // Get the system time and put it into 'now' as 'calender time'
+  walker = start;
+  while (walker->next != end) { // go through every atom and count
+    walker = walker->next;
+    No++;
+  }
+  if (out != NULL) {
+    for (int step=0;step<MDSteps;step++) {
+      *out << No << "\n\tCreated by molecuilder, step " << step << ", on " << ctime(&now);
+      walker = start;
+      while (walker->next != end) { // go through every atom of this element
+        walker = walker->next;
+        *out << walker->type->symbol << "\t" << Trajectories[walker].R.at(step).x[0] << "\t" << Trajectories[walker].R.at(step).x[1] << "\t" << Trajectories[walker].R.at(step).x[2] << endl;
+      }
+    }
+    return true;
+  } else
+    return false;
+};
+
 /** Prints molecule to *out as xyz file.
- * \param *out output stream
+* \param *out output stream
  */
 bool molecule::OutputXYZ(ofstream *out) const
@@ -2615 +2670 @@
         lines++;
       }
-      *out << Verbose(2) << "Scanned " << lines-1 << " lines." << endl;   // one endline too much
+      //*out << Verbose(2) << "Scanned " << lines-1 << " lines." << endl;   // one endline too much
       InputFile.clear();
       InputFile.seekg(ios::beg);
@@ -2627 +2682 @@
         InputFile.getline(buffer, MAXSTRINGSIZE);
         if (strlen(buffer) > 2) {
-          //*out << Verbose(2) << "Scanning: " << buffer;
+          //*out << Verbose(2) << "Scanning: " << buffer << endl;
           stringstream line(buffer);
           line >> FragOrder;
@@ -2634 +2689 @@
           line >> ws >> Value;
           No -= 1;  // indices start at 1 in file, not 0 
-          //*out << Verbose(2) << " - yields (" << No << "," << Value << ")" << endl;
+          //*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>( Value, Order) ));
+            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
@@ -2646 +2702 @@
                 //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 = Value;
+                (*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 {
@@ -2664 +2722 @@
         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 ((*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 {
@@ -2713 +2774 @@
       }
     }
-    if ((Order == 0) && (AtomMask[AtomCount] == true))  // single stepping, just check
+    if ((Order == 0) && (AtomMask[AtomCount] == false))  // single stepping, just check
       status = true;
       
-    if (!status)
-      *out << Verbose(1) << "Order at every site is already equal or above desired order " << Order << "." << endl;
+    if (!status) {
+      if (Order == 0)
+        *out << Verbose(1) << "Single stepping done." << endl;
+      else
+        *out << Verbose(1) << "Order at every site is already equal or above desired order " << Order << "." << endl;
+    }
   }
   
@@ -2732 +2797 @@
 };
 
-/** Create a SortIndex to map from BFS labels to the sequence in which the atoms are given in the config file.
+/** Create a SortIndex to map from atomic labels to the sequence in which the atoms are given in the config file.
  * \param *out output stream for debugging
  * \param *&SortIndex Mapping array of size molecule::AtomCount
@@ -2781 +2846 @@
  * \param *configuration configuration for writing config files for each fragment
  */
-void molecule::FragmentMolecule(ofstream *out, int Order, config *configuration)
+int molecule::FragmentMolecule(ofstream *out, int Order, config *configuration)
 {
         MoleculeListClass *BondFragments = NULL;
@@ -2791 +2856 @@
   fstream File;
   bool FragmentationToDo = true;
+  bool CheckOrder = false;
   Graph **FragmentList = NULL;
   Graph *ParsedFragmentList = NULL;
@@ -2837 +2903 @@
   KeyStack *RootStack = new KeyStack[Subgraphs->next->Count()];
   AtomMask = new bool[AtomCount+1];
-  while (CheckOrderAtSite(out, AtomMask, ParsedFragmentList, Order, MinimumRingSize, configuration->configpath)) {
+  AtomMask[AtomCount] = false;
+  FragmentationToDo = false;  // if CheckOrderAtSite just ones recommends fragmentation, we will save fragments afterwards
+  while ((CheckOrder = CheckOrderAtSite(out, AtomMask, ParsedFragmentList, Order, MinimumRingSize, configuration->configpath))) {
+    FragmentationToDo = FragmentationToDo || CheckOrder;
     AtomMask[AtomCount] = true;   // last plus one entry is used as marker that we have been through this loop once already in CheckOrderAtSite()
     // ===== 6b. fill RootStack for each subgraph (second adaptivity check) ===== 
@@ -2889 +2958 @@
 
   // ===== 8b. gather keyset lists (graphs) from all subgraphs and transform into MoleculeListClass =====
-  // allocate memory for the pointer array and transmorph graphs into full molecular fragments
-  BondFragments = new MoleculeListClass(TotalGraph.size(), AtomCount);
-  int k=0;
-  for(Graph::iterator runner = TotalGraph.begin(); runner != TotalGraph.end(); runner++) {
-    KeySet test = (*runner).first;
-    *out << "Fragment No." << (*runner).second.first << " with TEFactor " << (*runner).second.second << "." << endl;
-    BondFragments->ListOfMolecules[k] = StoreFragmentFromKeySet(out, test, configuration);
-    k++;
-  }
-  *out << k << "/" << BondFragments->NumberOfMolecules << " fragments generated from the keysets." << endl;
-  
-  // ===== 9. Save fragments' configuration and keyset files et al to disk ===
-  if (BondFragments->NumberOfMolecules != 0) {
-    // create the SortIndex from BFS labels to order in the config file
-    CreateMappingLabelsToConfigSequence(out, SortIndex);
+  //if (FragmentationToDo) {    // we should always store the fragments again as coordination might have changed slightly without changing bond structure
+    // allocate memory for the pointer array and transmorph graphs into full molecular fragments
+    BondFragments = new MoleculeListClass(TotalGraph.size(), AtomCount);
+    int k=0;
+    for(Graph::iterator runner = TotalGraph.begin(); runner != TotalGraph.end(); runner++) {
+      KeySet test = (*runner).first;
+      *out << "Fragment No." << (*runner).second.first << " with TEFactor " << (*runner).second.second << "." << endl;
+      BondFragments->ListOfMolecules[k] = StoreFragmentFromKeySet(out, test, configuration);
+      k++;
+    }
+    *out << k << "/" << BondFragments->NumberOfMolecules << " fragments generated from the keysets." << endl;
     
-    *out << Verbose(1) << "Writing " << BondFragments->NumberOfMolecules << " possible bond fragmentation configs" << endl;
-    if (BondFragments->OutputConfigForListOfFragments(out, configuration, SortIndex))
-      *out << Verbose(1) << "All configs written." << endl;
-    else
-      *out << Verbose(1) << "Some config writing failed." << endl;
-
-    // store force index reference file
-    BondFragments->StoreForcesFile(out, configuration->configpath, SortIndex);
-    
-    // store keysets file 
-    StoreKeySetFile(out, TotalGraph, configuration->configpath);
-
-    // store Adjacency file 
-    StoreAdjacencyToFile(out, configuration->configpath);
-
-    // store Hydrogen saturation correction file
-    BondFragments->AddHydrogenCorrection(out, configuration->configpath);
-    
-    // store adaptive orders into file
-    StoreOrderAtSiteFile(out, configuration->configpath);
-    
-    // restore orbital and Stop values
-    CalculateOrbitals(*configuration);
-    
-    // free memory for bond part
-    *out << Verbose(1) << "Freeing bond memory" << endl;
-    delete(FragmentList); // remove bond molecule from memory
-    Free((void **)&SortIndex, "molecule::FragmentMolecule: *SortIndex"); 
-  } else
-    *out << Verbose(1) << "FragmentList is zero on return, splitting failed." << endl;
-  
+    // ===== 9. Save fragments' configuration and keyset files et al to disk ===
+    if (BondFragments->NumberOfMolecules != 0) {
+      // create the SortIndex from BFS labels to order in the config file
+      CreateMappingLabelsToConfigSequence(out, SortIndex);
+      
+      *out << Verbose(1) << "Writing " << BondFragments->NumberOfMolecules << " possible bond fragmentation configs" << endl;
+      if (BondFragments->OutputConfigForListOfFragments(out, configuration, SortIndex))
+        *out << Verbose(1) << "All configs written." << endl;
+      else
+        *out << Verbose(1) << "Some config writing failed." << endl;
+  
+      // store force index reference file
+      BondFragments->StoreForcesFile(out, configuration->configpath, SortIndex);
+      
+      // store keysets file 
+      StoreKeySetFile(out, TotalGraph, configuration->configpath);
+  
+      // store Adjacency file 
+      StoreAdjacencyToFile(out, configuration->configpath);
+  
+      // store Hydrogen saturation correction file
+      BondFragments->AddHydrogenCorrection(out, configuration->configpath);
+      
+      // store adaptive orders into file
+      StoreOrderAtSiteFile(out, configuration->configpath);
+      
+      // restore orbital and Stop values
+      CalculateOrbitals(*configuration);
+      
+      // free memory for bond part
+      *out << Verbose(1) << "Freeing bond memory" << endl;
+      delete(FragmentList); // remove bond molecule from memory
+      Free((void **)&SortIndex, "molecule::FragmentMolecule: *SortIndex"); 
+    } else
+      *out << Verbose(1) << "FragmentList is zero on return, splitting failed." << endl;
+  //} else 
+  //  *out << Verbose(1) << "No fragments to store." << endl;
   *out << Verbose(0) << "End of bond fragmentation." << endl;
+
+  return ((int)(!FragmentationToDo)+1);    // 1 - continue, 2 - stop (no fragmentation occured)
 };
 
@@ -2957 +3030 @@
     while (Walker->next != end) {
       Walker = Walker->next;
-      file << Walker->nr << "\t" << (int)Walker->AdaptiveOrder << endl;
-      *out << Verbose(2) << "Storing: " << Walker->nr << "\t" << (int)Walker->AdaptiveOrder << "." << endl;
+      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;
     }
     file.close();
@@ -2979 +3052 @@
 {
   unsigned char *OrderArray = (unsigned char *) Malloc(sizeof(unsigned char)*AtomCount, "molecule::ParseOrderAtSiteFromFile - *OrderArray");
+  bool *MaxArray = (bool *) Malloc(sizeof(bool)*AtomCount, "molecule::ParseOrderAtSiteFromFile - *MaxArray");
   bool status;
-  int AtomNr;
+  int AtomNr, value;
   stringstream line;
   ifstream file;
-  int Order;
 
   *out << Verbose(1) << "Begin of ParseOrderAtSiteFromFile" << endl;
@@ -2991 +3064 @@
   file.open(line.str().c_str());
   if (file != NULL) {
-    for (int i=AtomCount;i--;) // initialise with 0
+    for (int i=AtomCount;i--;) { // initialise with 0
       OrderArray[i] = 0;
+      MaxArray[i] = 0;
+    }
     while (!file.eof()) { // parse from file
+      AtomNr = -1;
       file >> AtomNr;
-      file >> Order;
-      OrderArray[AtomNr] = (unsigned char) Order;
-      //*out << Verbose(2) << "AtomNr " << AtomNr << " with order " << (int)OrderArray[AtomNr] << "." << endl;
+      if (AtomNr != -1) {   // test whether we really parsed something (this is necessary, otherwise last atom is set twice and to 0 on second time)
+        file >> value;
+        OrderArray[AtomNr] = value;
+        file >> value;
+        MaxArray[AtomNr] = value;
+        //*out << Verbose(2) << "AtomNr " << AtomNr << " with order " << (int)OrderArray[AtomNr] << " and max order set to " << (int)MaxArray[AtomNr] << "." << endl;
+      }
     }
     atom *Walker = start;
@@ -3003 +3083 @@
       Walker = Walker->next;
       Walker->AdaptiveOrder = OrderArray[Walker->nr];
-      *out << Verbose(2) << *Walker << " gets order " << (int)Walker->AdaptiveOrder << "." << endl; 
+      Walker->MaxOrder = MaxArray[Walker->nr];
+      *out << Verbose(2) << *Walker << " gets order " << (int)Walker->AdaptiveOrder << " and is " << (!Walker->MaxOrder ? "not " : " ") << "maxed." << endl; 
     }
     file.close();
@@ -3013 +3094 @@
   }
   Free((void **)&OrderArray, "molecule::ParseOrderAtSiteFromFile - *OrderArray");
+  Free((void **)&MaxArray, "molecule::ParseOrderAtSiteFromFile - *MaxArray");
   
   *out << Verbose(1) << "End of ParseOrderAtSiteFromFile" << endl;
@@ -3602 +3684 @@
   int FragmentCounter;
   int CurrentIndex;
-  int *Labels;
   double TEFactor;
   int *ShortestPathList;
@@ -3670 +3751 @@
                 OtherWalker = BondsSet[j]->rightatom;   // rightatom is always the one more distant, i.e. the one to add 
           //*out << Verbose(1+verbosity) << "Current Bond is " << ListOfBondsPerAtom[Walker->nr][i] << ", checking on " << *OtherWalker << "." << endl;
-          //if ((!FragmentSearch->UsedList[OtherWalker->nr][i]) && (FragmentSearch->Labels[OtherWalker->nr] > FragmentSearch->Labels[FragmentSearch->Root->nr])) { 
-            //*out << Verbose(2+verbosity) << "Not used so far, label " << FragmentSearch->Labels[OtherWalker->nr] << " is bigger than Root's " << FragmentSearch->Labels[FragmentSearch->Root->nr] << "." << endl;
           *out << Verbose(2+verbosity) << "Adding " << *OtherWalker << " with nr " << OtherWalker->nr << "." << endl;
           TestKeySetInsert = FragmentSearch->FragmentSet->insert(OtherWalker->nr); 
@@ -3730 +3809 @@
         InsertFragmentIntoGraph(out, FragmentSearch);
         //Removal = LookForRemovalCandidate(out, FragmentSearch->FragmentSet, FragmentSearch->ShortestPathList);
-        //Removal = StoreFragmentFromStack(out, FragmentSearch->Root, FragmentSearch->Leaflet, FragmentSearch->FragmentStack, FragmentSearch->ShortestPathList,FragmentSearch->Labels, &FragmentSearch->FragmentCounter, FragmentSearch->configuration);
+        //Removal = StoreFragmentFromStack(out, FragmentSearch->Root, FragmentSearch->Leaflet, FragmentSearch->FragmentStack, FragmentSearch->ShortestPathList, &FragmentSearch->FragmentCounter, FragmentSearch->configuration);
       }
       
@@ -3821 +3900 @@
 int molecule::PowerSetGenerator(ofstream *out, int Order, struct UniqueFragments &FragmentSearch, KeySet RestrictedKeySet) 
 {
-  int SP, UniqueIndex, AtomKeyNr;
+  int SP, AtomKeyNr;
   atom *Walker = NULL, *OtherWalker = NULL, *Predecessor = NULL;
   bond *Binder = NULL;
   bond *CurrentEdge = NULL;
   bond **BondsList = NULL;
-  int RootKeyNr = FragmentSearch.Root->nr;
+  int RootKeyNr = FragmentSearch.Root->GetTrueFather()->nr;
   int Counter = FragmentSearch.FragmentCounter;
   int RemainingWalkers;
@@ -3834 +3913 @@
 
   // prepare Label and SP arrays of the BFS search
-  UniqueIndex = 0;
-  if (FragmentSearch.Labels[RootKeyNr] == -1)
-    FragmentSearch.Labels[RootKeyNr] = UniqueIndex++;
-  FragmentSearch.ShortestPathList[RootKeyNr] = 0;
+  FragmentSearch.ShortestPathList[FragmentSearch.Root->nr] = 0;
 
   // prepare root level (SP = 0) and a loop bond denoting Root
@@ -3871 +3947 @@
       Predecessor = CurrentEdge->leftatom;    // ... and leftatom is predecessor
       AtomKeyNr = Walker->nr;
-      *out << Verbose(0) << "Current Walker is: " << *Walker << " with nr " << Walker->nr << " and label " << FragmentSearch.Labels[AtomKeyNr] << " and SP of " << SP << ", with " << RemainingWalkers << " remaining walkers on this level." << endl;
+      *out << Verbose(0) << "Current Walker is: " << *Walker << " with nr " << Walker->nr << " and SP of " << SP << ", with " << RemainingWalkers << " remaining walkers on this level." << endl;
       // check for new sp level
       // go through all its bonds
@@ -3885 +3961 @@
           *out << Verbose(2) << "Current partner is " << *OtherWalker << " with nr " << OtherWalker->nr << " in bond " << *Binder << "." << endl;
           // set the label if not set (and push on root stack as well)
-          if (FragmentSearch.Labels[OtherWalker->nr] == -1) {
-            FragmentSearch.Labels[OtherWalker->nr] = UniqueIndex++;
-            *out << Verbose(3) << "Set label to " << FragmentSearch.Labels[OtherWalker->nr] << "." << endl;
-          } else {
-            *out << Verbose(3) << "Label is already " << FragmentSearch.Labels[OtherWalker->nr] << "." << endl;
-          }
-          if ((OtherWalker != Predecessor) && (OtherWalker->nr > RootKeyNr)) { // only pass through those with label bigger than Root's
+          if ((OtherWalker != Predecessor) && (OtherWalker->GetTrueFather()->nr > RootKeyNr)) { // only pass through those with label bigger than Root's
             FragmentSearch.ShortestPathList[OtherWalker->nr] = SP+1;
             *out << Verbose(3) << "Set Shortest Path to " << FragmentSearch.ShortestPathList[OtherWalker->nr] << "." << endl;
@@ -3899 +3969 @@
             FragmentSearch.BondsPerSPCount[SP+1]++;
             *out << Verbose(3) << "Added its bond to SP list, having now " << FragmentSearch.BondsPerSPCount[SP+1] << " item(s)." << endl;
-          } else *out << Verbose(3) << "Not passing on, as index of " << *OtherWalker << " " << OtherWalker->nr << " is smaller than that of Root " << RootKeyNr << " or this is my predecessor " << *Predecessor << "." << endl;
+          } else {
+            if (OtherWalker != Predecessor)
+              *out << Verbose(3) << "Not passing on, as index of " << *OtherWalker << " " << OtherWalker->GetTrueFather()->nr << " is smaller than that of Root " << RootKeyNr << "." << endl;
+            else
+              *out << Verbose(3) << "This is my predecessor " << *Predecessor << "." << endl;
+          }
         } else *out << Verbose(2) << "Is not in the restricted keyset or skipping hydrogen " << *OtherWalker << "." << endl;
       }
@@ -4003 +4078 @@
       for (int i=NDIM;i--;) {
         tmp = fabs(Binder->leftatom->x.x[i] - Binder->rightatom->x.x[i]);
-        //*out << Verbose(3) << "Checking " << i << "th distance of " << *Binder->leftatom << " to " << *Binder->rightatom << ": " << tmp << "." << endl;
+        *out << Verbose(3) << "Checking " << i << "th distance of " << *Binder->leftatom << " to " << *Binder->rightatom << ": " << tmp << "." << endl;
         if (tmp > BondDistance) {
           OtherBinder = Binder->next; // note down binding partner for later re-insertion
           unlink(Binder);   // unlink bond
-          //*out << Verbose(2) << "Correcting at bond " << *Binder << "." << endl;
+          *out << Verbose(2) << "Correcting at bond " << *Binder << "." << endl;
           flag = true;
           break;
@@ -4195 +4270 @@
   // initialise the fragments structure
   FragmentSearch.ShortestPathList = (int *) Malloc(sizeof(int)*AtomCount, "molecule::PowerSetGenerator: *ShortestPathList");
-  FragmentSearch.Labels = (int *) Malloc(sizeof(int)*AtomCount, "molecule::PowerSetGenerator: *Labels");
   FragmentSearch.FragmentCounter = 0;
   FragmentSearch.FragmentSet = new KeySet;
   FragmentSearch.Root = FindAtom(RootKeyNr);
   for (int i=AtomCount;i--;) {
-    FragmentSearch.Labels[i] = -1;
     FragmentSearch.ShortestPathList[i] = -1;
   }
@@ -4249 +4322 @@
       // create top order where nothing is reduced
       *out << Verbose(0) << "==============================================================================================================" << endl;
-      *out << Verbose(0) << "Creating KeySets of Bond Order " << Order << " for " << *Walker << ", NumLevels is " << NumLevels << ", " << (RootStack.size()-RootNr) << " Roots remaining." << endl;
+      *out << Verbose(0) << "Creating KeySets of Bond Order " << Order << " for " << *Walker << ", " << (RootStack.size()-RootNr) << " Roots remaining." << endl; // , NumLevels is " << NumLevels << "
   
       // Create list of Graphs of current Bond Order (i.e. F_{ij})
@@ -4260 +4333 @@
       if (NumMoleculesOfOrder[RootNr] != 0) {
         NumMolecules = 0;
+
+        // we don't have to dive into suborders! These keysets are all already created on lower orders!
+        // this was all ancient stuff, when we still depended on the TEFactors (and for those the suborders were needed)
         
-        if ((NumLevels >> 1) > 0) {
-          // create lower order fragments
-          *out << Verbose(0) << "Creating list of unique fragments of lower Bond Order terms to be subtracted." << endl;
-          Order = Walker->AdaptiveOrder;
-          for (int source=0;source<(NumLevels >> 1);source++) { // 1-terms don't need any more splitting, that's why only half is gone through (shift again)
-            // step down to next order at (virtual) boundary of powers of 2 in array
-            while (source >= (1 << (Walker->AdaptiveOrder-Order))) // (int)pow(2,Walker->AdaptiveOrder-Order))    
-              Order--;
-            *out << Verbose(0) << "Current Order is: " << Order << "." << endl;
-            for (int SubOrder=Order-1;SubOrder>0;SubOrder--) {
-              int dest = source + (1 << (Walker->AdaptiveOrder-(SubOrder+1)));
-              *out << Verbose(0) << "--------------------------------------------------------------------------------------------------------------" << endl;
-              *out << Verbose(0) << "Current SubOrder is: " << SubOrder << " with source " << source << " to destination " << dest << "." << endl;
-        
-              // every molecule is split into a list of again (Order - 1) molecules, while counting all molecules
-              //*out << Verbose(1) << "Splitting the " << (*FragmentLowerOrdersList[RootNr][source]).size() << " molecules of the " << source << "th cell in the array." << endl;
-              //NumMolecules = 0;
-              FragmentLowerOrdersList[RootNr][dest] = new Graph;
-              for(Graph::iterator runner = (*FragmentLowerOrdersList[RootNr][source]).begin();runner != (*FragmentLowerOrdersList[RootNr][source]).end(); runner++) {
-                for (KeySet::iterator sprinter = (*runner).first.begin();sprinter != (*runner).first.end(); sprinter++) {
-                  Graph TempFragmentList;
-                  FragmentSearch.TEFactor = -(*runner).second.second;
-                  FragmentSearch.Leaflet = &TempFragmentList;      // set to insertion graph
-                  FragmentSearch.Root = FindAtom(*sprinter);
-                  NumMoleculesOfOrder[RootNr] += PowerSetGenerator(out, SubOrder, FragmentSearch, (*runner).first);
-                  // insert new keysets FragmentList into FragmentLowerOrdersList[Walker->AdaptiveOrder-1][dest]
-                  *out << Verbose(1) << "Merging resulting key sets with those present in destination " << dest << "." << endl;
-                  InsertGraphIntoGraph(out, *FragmentLowerOrdersList[RootNr][dest], TempFragmentList, &NumMolecules);
-                }
-              }
-              *out << Verbose(1) << "Number of resulting molecules for SubOrder " << SubOrder << " is: " << NumMolecules << "." << endl;
-            }
-          }
-        }
+//        if ((NumLevels >> 1) > 0) {
+//          // create lower order fragments
+//          *out << Verbose(0) << "Creating list of unique fragments of lower Bond Order terms to be subtracted." << endl;
+//          Order = Walker->AdaptiveOrder;
+//          for (int source=0;source<(NumLevels >> 1);source++) { // 1-terms don't need any more splitting, that's why only half is gone through (shift again)
+//            // step down to next order at (virtual) boundary of powers of 2 in array
+//            while (source >= (1 << (Walker->AdaptiveOrder-Order))) // (int)pow(2,Walker->AdaptiveOrder-Order))    
+//              Order--;
+//            *out << Verbose(0) << "Current Order is: " << Order << "." << endl;
+//            for (int SubOrder=Order-1;SubOrder>0;SubOrder--) {
+//              int dest = source + (1 << (Walker->AdaptiveOrder-(SubOrder+1)));
+//              *out << Verbose(0) << "--------------------------------------------------------------------------------------------------------------" << endl;
+//              *out << Verbose(0) << "Current SubOrder is: " << SubOrder << " with source " << source << " to destination " << dest << "." << endl;
+//        
+//              // every molecule is split into a list of again (Order - 1) molecules, while counting all molecules
+//              //*out << Verbose(1) << "Splitting the " << (*FragmentLowerOrdersList[RootNr][source]).size() << " molecules of the " << source << "th cell in the array." << endl;
+//              //NumMolecules = 0;
+//              FragmentLowerOrdersList[RootNr][dest] = new Graph;
+//              for(Graph::iterator runner = (*FragmentLowerOrdersList[RootNr][source]).begin();runner != (*FragmentLowerOrdersList[RootNr][source]).end(); runner++) {
+//                for (KeySet::iterator sprinter = (*runner).first.begin();sprinter != (*runner).first.end(); sprinter++) {
+//                  Graph TempFragmentList;
+//                  FragmentSearch.TEFactor = -(*runner).second.second;
+//                  FragmentSearch.Leaflet = &TempFragmentList;      // set to insertion graph
+//                  FragmentSearch.Root = FindAtom(*sprinter);
+//                  NumMoleculesOfOrder[RootNr] += PowerSetGenerator(out, SubOrder, FragmentSearch, (*runner).first);
+//                  // insert new keysets FragmentList into FragmentLowerOrdersList[Walker->AdaptiveOrder-1][dest]
+//                  *out << Verbose(1) << "Merging resulting key sets with those present in destination " << dest << "." << endl;
+//                  InsertGraphIntoGraph(out, *FragmentLowerOrdersList[RootNr][dest], TempFragmentList, &NumMolecules);
+//                }
+//              }
+//              *out << Verbose(1) << "Number of resulting molecules for SubOrder " << SubOrder << " is: " << NumMolecules << "." << endl;
+//            }
+//          }
+//        }
       } else {
-        *out << Verbose(1) << "Hence, we don't dive into SubOrders ... " << endl;
+        Walker->GetTrueFather()->MaxOrder = true;
+//        *out << Verbose(1) << "Hence, we don't dive into SubOrders ... " << endl;
       }
       // now, we have completely filled each cell of FragmentLowerOrdersList[] for the current Walker->AdaptiveOrder
       //NumMoleculesOfOrder[Walker->AdaptiveOrder-1] = NumMolecules;
       TotalNumMolecules += NumMoleculesOfOrder[RootNr];
-      *out << Verbose(1) << "Number of resulting molecules for Order " << (int)Walker->GetTrueFather()->AdaptiveOrder << " is: " << NumMoleculesOfOrder[RootNr] << "." << endl;
+//      *out << Verbose(1) << "Number of resulting molecules for Order " << (int)Walker->GetTrueFather()->AdaptiveOrder << " is: " << NumMoleculesOfOrder[RootNr] << "." << endl;
       RootStack.push_back(RootKeyNr); // put back on stack
       RootNr++;
@@ -4320 +4397 @@
   // cleanup FragmentSearch structure
   Free((void **)&FragmentSearch.ShortestPathList, "molecule::PowerSetGenerator: *ShortestPathList");
-  Free((void **)&FragmentSearch.Labels, "molecule::PowerSetGenerator: *Labels");
   delete(FragmentSearch.FragmentSet);