Index: src/Makefile.am
===================================================================
--- src/Makefile.am	(revision f054074bc0e2ebd43408e39a6e2ae0cd96597630)
+++ src/Makefile.am	(revision 362b0e5abaff6113ef0ca4989b9074448eccb8b8)
@@ -1,4 +1,4 @@
-SOURCE = atom.cpp bond.cpp boundary.cpp builder.cpp config.cpp element.cpp helpers.cpp molecules.cpp moleculelist.cpp periodentafel.cpp vector.cpp verbose.cpp
-HEADER = boundary.hpp defs.hpp helpers.hpp molecules.hpp periodentafel.hpp stackclass.hpp vector.hpp
+SOURCE = atom.cpp bond.cpp boundary.cpp builder.cpp config.cpp element.cpp helpers.cpp molecules.cpp moleculelist.cpp parser.cpp periodentafel.cpp vector.cpp verbose.cpp
+HEADER = boundary.hpp defs.hpp helpers.hpp molecules.hpp parser.hpp periodentafel.hpp stackclass.hpp vector.hpp
 
 bin_PROGRAMS = molecuilder joiner analyzer
Index: src/atom.cpp
===================================================================
--- src/atom.cpp	(revision f054074bc0e2ebd43408e39a6e2ae0cd96597630)
+++ src/atom.cpp	(revision 362b0e5abaff6113ef0ca4989b9074448eccb8b8)
@@ -28,4 +28,5 @@
   LowpointNr = -1;
   AdaptiveOrder = 0;
+  MaxOrder = false;
 };
 
Index: src/builder.cpp
===================================================================
--- src/builder.cpp	(revision f054074bc0e2ebd43408e39a6e2ae0cd96597630)
+++ src/builder.cpp	(revision 362b0e5abaff6113ef0ca4989b9074448eccb8b8)
@@ -703,8 +703,15 @@
     output.open(filename, ios::trunc);
   } 
-  if (mol->OutputXYZ(&output))
-    cout << Verbose(0) << "Saving of XYZ file successful." << endl;
-  else
-    cout << Verbose(0) << "Saving of XYZ file failed." << endl;
+  if (mol->MDSteps <= 1) {
+    if (mol->OutputXYZ(&output))
+      cout << Verbose(0) << "Saving of XYZ file successful." << endl;
+    else
+      cout << Verbose(0) << "Saving of XYZ file failed." << endl;
+  } else {
+    if (mol->OutputTrajectoriesXYZ(&output))
+      cout << Verbose(0) << "Saving of XYZ file successful." << endl;
+    else
+      cout << Verbose(0) << "Saving of XYZ file failed." << endl;
+  }
   output.close();
   output.clear();
@@ -733,5 +740,6 @@
   string line;
   atom *first;
-  int ExitFlag = 0;
+  bool SaveFlag = false;
+  int ExitFlag = 1;
   int j;
   double volume = 0.;
@@ -758,8 +766,9 @@
             cout << "\t-b x1 x2 x3\tCenter atoms in domain with given edge lengths of (x1,x2,x3)." << endl;
             cout << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl;
+            cout << "\t-D <bond distance>\tDepth-First-Search Analysis of the molecule, giving cycles and tree/back edges." << endl;
             cout << "\t-O\tCenter atoms in origin." << endl;
             cout << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl;
             cout << "\t-e <file>\tSets the databases path to be parsed (default: ./)." << endl;
-            cout << "\t-f/F <dist> <order>\tFragments the molecule in BOSSANOVA manner (with/out rings compressed) and stores config files in same dir as config." << endl;
+            cout << "\t-f/F <dist> <order>\tFragments the molecule in BOSSANOVA manner (with/out rings compressed) and stores config files in same dir as config (return code 0 - fragmented, 2 - no fragmentation necessary)." << endl;
             cout << "\t-h/-H/-?\tGive this help screen." << endl;
             cout << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl;
@@ -767,5 +776,5 @@
             cout << "\t-o\tGet volume of the convex envelope (and store to tecplot file)." << endl;
             cout << "\t-p <file>\tParse given xyz file and create raw config file from it." << endl;
-            cout << "\t-P <file> <delta_t>\tParse given forces file and append as an MD step with width delta_t to config file via Verlet." << endl;
+            cout << "\t-P <file>\tParse given forces file and append as an MD step to config file via Verlet." << endl;
             cout << "\t-r\t\tConvert file from an old pcp syntax." << endl;
             cout << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << endl;
@@ -860,5 +869,5 @@
           switch(argv[argptr-1][1]) {
             case 'p':
-              ExitFlag = 1;
+              SaveFlag = true;
               cout << Verbose(1) << "Parsing xyz file for new atoms." << endl;
               if (!mol->AddXYZFile(argv[argptr]))
@@ -875,15 +884,34 @@
         if (config_present == present) {
           switch(argv[argptr-1][1]) {
+            case 'D': 
+              {
+                cout << Verbose(1) << "Depth-First-Search Analysis." << endl;
+                MoleculeLeafClass *Subgraphs = NULL;      // list of subgraphs from DFS analysis
+                int *MinimumRingSize = NULL;
+                mol->CreateAdjacencyList((ofstream *)&cout, atof(argv[argptr++]), configuration.GetIsAngstroem());
+                mol->CreateListOfBondsPerAtom((ofstream *)&cout);
+                Subgraphs = mol->DepthFirstSearchAnalysis((ofstream *)&cout, MinimumRingSize);
+                delete[](MinimumRingSize);
+                if (Subgraphs != NULL) {
+                  while (Subgraphs->next != NULL) {
+                    Subgraphs = Subgraphs->next;
+                    delete(Subgraphs->previous);
+                  }
+                  delete(Subgraphs);
+                }
+              }
+              argptr+=1;
+              break;
             case 'P':
-              ExitFlag = 1;
+              SaveFlag = true;
               cout << Verbose(1) << "Parsing forces file and Verlet integrating." << endl;
-              if (!mol->VerletForceIntegration(argv[argptr], atof(argv[argptr+1])))
+              if (!mol->VerletForceIntegration(argv[argptr], configuration.Deltat, configuration.GetIsAngstroem()))
                 cout << Verbose(2) << "File not found." << endl;
               else
                 cout << Verbose(2) << "File found and parsed." << endl;
-              argptr+=2;
+              argptr+=1;
               break;
             case 't':
-              ExitFlag = 1;
+              SaveFlag = true;
               cout << Verbose(1) << "Translating all ions to new origin." << endl;
               for (int i=NDIM;i--;)
@@ -893,5 +921,5 @@
               break;
             case 'a':
-              ExitFlag = 1;
+              SaveFlag = true;
               cout << Verbose(1) << "Adding new atom with element " << argv[argptr] << " at (" << argv[argptr+1] << "," << argv[argptr+2] << "," << argv[argptr+3] << "), ";
               first = new atom;
@@ -908,5 +936,5 @@
               break;
             case 's':
-              ExitFlag = 1;
+              SaveFlag = true;
               j = -1;
               cout << Verbose(1) << "Scaling all ion positions by factor." << endl;
@@ -929,5 +957,5 @@
               break;
             case 'b':
-              ExitFlag = 1;
+              SaveFlag = true;
               j = -1;
               cout << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
@@ -944,5 +972,5 @@
               break;
             case 'c':
-              ExitFlag = 1;
+              SaveFlag = true;
               j = -1;
               cout << Verbose(1) << "Centering atoms in config file within given additional boundary." << endl;
@@ -961,5 +989,5 @@
               break;
             case 'O':
-              ExitFlag = 1;
+              SaveFlag = true;
               cout << Verbose(1) << "Centering atoms in origin." << endl;
               mol->CenterOrigin((ofstream *)&cout, &x);
@@ -967,10 +995,9 @@
               break;
             case 'r':
-              ExitFlag = 1;
+              SaveFlag = true;
               cout << Verbose(1) << "Converting config file from supposed old to new syntax." << endl;
               break;
             case 'F':
             case 'f':
-              if (ExitFlag ==0) ExitFlag = 2;  // only set if not already by other command line switch
               cout << "Fragmenting molecule with bond distance " << argv[argptr] << " angstroem, order of " << argv[argptr+1] << "." << endl;
               if (argc >= argptr+2) {
@@ -980,5 +1007,5 @@
                 cout << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;
                 if (mol->first->next != mol->last) {
-                  mol->FragmentMolecule((ofstream *)&cout, atoi(argv[argptr]), &configuration);
+                  ExitFlag = mol->FragmentMolecule((ofstream *)&cout, atoi(argv[argptr]), &configuration);
                 }
                 end = clock();
@@ -990,5 +1017,4 @@
               break;
             case 'm':
-              ExitFlag = 1;
               j = atoi(argv[argptr++]);
               if ((j<0) || (j>1)) {
@@ -996,12 +1022,13 @@
                 j = 0;
               }
-              if (j)
+              if (j) {
+                SaveFlag = true;
                 cout << Verbose(0) << "Converting to prinicipal axis system." << endl;
-              else
+              } else
                 cout << Verbose(0) << "Evaluating prinicipal axis." << endl;
               mol->PrincipalAxisSystem((ofstream *)&cout, (bool)j);
               break;
             case 'o':
-              ExitFlag = 1;
+              SaveFlag = true;
               cout << Verbose(0) << "Evaluating volume of the convex envelope.";
               VolumeOfConvexEnvelope((ofstream *)&cout, &configuration, NULL, mol);
@@ -1013,5 +1040,5 @@
               {
                 double density;
-                ExitFlag = 1;
+                SaveFlag = true;
                 cout << Verbose(0) << "Evaluating necessary cell volume for a cluster suspended in water.";
                 density = atof(argv[argptr++]);
@@ -1030,4 +1057,5 @@
               break;
             case 'd':
+              SaveFlag = true;
               for (int axis = 1; axis <= NDIM; axis++) {
                 int faktor = atoi(argv[argptr++]);
@@ -1089,10 +1117,10 @@
       } else argptr++;
     } while (argptr < argc);
-    if (ExitFlag == 1) // 1 means save and exit
+    if (SaveFlag)
       SaveConfig(ConfigFileName, &configuration, periode, mol);
-    if (ExitFlag >= 1) { // 2 means just exit
+    if ((ExitFlag >= 1)) {
       delete(mol);
       delete(periode);
-      return (1);
+      return (ExitFlag);
     }
   } else {  // no arguments, hence scan the elements db
Index: src/config.cpp
===================================================================
--- src/config.cpp	(revision f054074bc0e2ebd43408e39a6e2ae0cd96597630)
+++ src/config.cpp	(revision 362b0e5abaff6113ef0ca4989b9074448eccb8b8)
@@ -42,5 +42,5 @@
   
   MaxOuterStep=0;
-  Deltat=0;
+  Deltat=1;
   OutVisStep=10;
   OutSrcStep=5;
@@ -461,5 +461,5 @@
   int Z, No[MAX_ELEMENTS];
   int verbose = 0;
-  int MinSteps = -1;
+  double value[3];
   
   /* Namen einlesen */
@@ -511,5 +511,6 @@
   
   ParseForParameter(verbose,file,"MaxOuterStep", 0, 1, 1, int_type, &(config::MaxOuterStep), 1, critical);
-  ParseForParameter(verbose,file,"Deltat", 0, 1, 1, double_type, &(config::Deltat), 1, optional);
+  if (!ParseForParameter(verbose,file,"Deltat", 0, 1, 1, double_type, &(config::Deltat), 1, optional))
+    config::Deltat = 1;
   ParseForParameter(verbose,file,"OutVisStep", 0, 1, 1, int_type, &(config::OutVisStep), 1, optional);
   ParseForParameter(verbose,file,"OutSrcStep", 0, 1, 1, int_type, &(config::OutSrcStep), 1, optional);
@@ -643,75 +644,110 @@
     cout << Verbose(1) << i << ". Z = " << elementhash[i]->Z << " with " << No[i] << " ions." << endl; 
   }
-  for (int i=0; i < config::MaxTypes; i++) {
-    sprintf(name,"Ion_Type%i",i+1);
-    for(int j=0;j<No[i];j++) {
-      int repetition = 1; // which repeated keyword shall be read
-      sprintf(keyword,"%s_%i",name, j+1);
-      atom *neues = new atom();
-      // then parse for each atom the coordinates as often as present
-      if (!FastParsing) {
-        while (ParseForParameter(verbose,file, keyword, 0, 1, 1, double_type, &neues->x.x[0], repetition, (repetition == 1) ? critical : optional) &&
-          ParseForParameter(verbose,file, keyword, 0, 2, 1, double_type, &neues->x.x[1], repetition, (repetition == 1) ? critical : optional) &&
-          ParseForParameter(verbose,file, keyword, 0, 3, 1, double_type, &neues->x.x[2], repetition, (repetition == 1) ? critical : optional)) {
+  int repetition = 0; // which repeated keyword shall be read
+
+  map<int, atom *> AtomList[config::MaxTypes];
+  if (!FastParsing) {
+    // parse in trajectories
+    bool status = true;
+    atom *neues = NULL;
+    while (status) {
+      for (int i=0; i < config::MaxTypes; i++) {
+        sprintf(name,"Ion_Type%i",i+1);
+        for(int j=0;j<No[i];j++) {
+          sprintf(keyword,"%s_%i",name, j+1);
+          if (repetition == 0) {
+            neues = new atom();
+            AtomList[i][j] = neues;
+            neues->type = elementhash[i]; // find element type
+            mol->AddAtom(neues);
+          } else
+            neues = AtomList[i][j];
+          status = (status && 
+                  ParseForParameter(verbose,file, keyword, 0, 1, 1, double_type, &neues->x.x[0], 1, (repetition == 0) ? critical : optional) &&
+                  ParseForParameter(verbose,file, keyword, 0, 2, 1, double_type, &neues->x.x[1], 1, (repetition == 0) ? critical : optional) &&
+                  ParseForParameter(verbose,file, keyword, 0, 3, 1, double_type, &neues->x.x[2], 1, (repetition == 0) ? critical : optional) &&
+                  ParseForParameter(verbose,file, keyword, 0, 4, 1, int_type, &neues->FixedIon, 1, (repetition == 0) ? critical : optional));
+          if (!status) break;
+  
           // check size of vectors
-          if (mol->Trajectories[neues].R.size() <= (unsigned int)(repetition-1)) {
-            cout << "Increasing size for trajectory array of " << keyword << " to " << (repetition+10) << "." << endl;
+          if (mol->Trajectories[neues].R.size() <= (unsigned int)(repetition)) {
+            //cout << "Increasing size for trajectory array of " << keyword << " to " << (repetition+10) << "." << endl;
             mol->Trajectories[neues].R.resize(repetition+10);
             mol->Trajectories[neues].U.resize(repetition+10);
             mol->Trajectories[neues].F.resize(repetition+10);
           }
-            
+        
           // put into trajectories list
           for (int d=0;d<NDIM;d++)
-            mol->Trajectories[neues].R.at(repetition-1).x[d] = neues->x.x[d];
+            mol->Trajectories[neues].R.at(repetition).x[d] = neues->x.x[d];
           
           // parse velocities if present
-          if(!ParseForParameter(verbose,file, keyword, 0, 5, 1, double_type, &neues->v.x[0], repetition,optional))
+          if(!ParseForParameter(verbose,file, keyword, 0, 5, 1, double_type, &neues->v.x[0], 1,optional))
             neues->v.x[0] = 0.;
-          if(!ParseForParameter(verbose,file, keyword, 0, 6, 1, double_type, &neues->v.x[1], repetition,optional))
+          if(!ParseForParameter(verbose,file, keyword, 0, 6, 1, double_type, &neues->v.x[1], 1,optional))
             neues->v.x[1] = 0.;
-          if(!ParseForParameter(verbose,file, keyword, 0, 7, 1, double_type, &neues->v.x[2], repetition,optional))
+          if(!ParseForParameter(verbose,file, keyword, 0, 7, 1, double_type, &neues->v.x[2], 1,optional))
             neues->v.x[2] = 0.;
           for (int d=0;d<NDIM;d++)
-            mol->Trajectories[neues].U[repetition-1].x[d] = neues->v.x[d];
-          
+            mol->Trajectories[neues].U.at(repetition).x[d] = neues->v.x[d];
+    
           // parse forces if present
-          if(!ParseForParameter(verbose,file, keyword, 0, 8, 1, double_type, &neues->v.x[0], repetition,optional))
-            neues->v.x[0] = 0.;
-          if(!ParseForParameter(verbose,file, keyword, 0, 9, 1, double_type, &neues->v.x[1], repetition,optional))
-            neues->v.x[1] = 0.;
-          if(!ParseForParameter(verbose,file, keyword, 0, 10, 1, double_type, &neues->v.x[2], repetition,optional))
-            neues->v.x[2] = 0.;
+          if(!ParseForParameter(verbose,file, keyword, 0, 8, 1, double_type, &value[0], 1,optional))
+            value[0] = 0.;
+          if(!ParseForParameter(verbose,file, keyword, 0, 9, 1, double_type, &value[1], 1,optional))
+            value[1] = 0.;
+          if(!ParseForParameter(verbose,file, keyword, 1, 10, 1, double_type, &value[2], 1,optional))
+            value[2] = 0.;
           for (int d=0;d<NDIM;d++)
-            mol->Trajectories[neues].F[repetition-1].x[d] = neues->v.x[d];
-          
-          cout << "Parsed position of step " << (repetition-1) << " "; 
-          for (int d=0;d<NDIM;d++)
-            cout << mol->Trajectories[neues].R.at(repetition-1).x[d] << "  ";          // next step
-          cout << endl;
-          
-          repetition++;
+            mol->Trajectories[neues].F.at(repetition).x[d] = value[d];
+
+//            cout << "Parsed position of step " << (repetition) << ": ("; 
+//            for (int d=0;d<NDIM;d++)
+//              cout << mol->Trajectories[neues].R.at(repetition).x[d] << " ";          // next step
+//            cout << ")\t(";
+//            for (int d=0;d<NDIM;d++)
+//              cout << mol->Trajectories[neues].U.at(repetition).x[d] << " ";          // next step
+//            cout << ")\t(";
+//            for (int d=0;d<NDIM;d++)
+//              cout << mol->Trajectories[neues].F.at(repetition).x[d] << " ";          // next step
+//            cout << ")" << endl;
         }
-        repetition--;
-        cout << "Found " << repetition << " times of the keyword " << keyword << "." << endl;
       }
-      if ((repetition < MinSteps) || (MinSteps == -1))
-        MinSteps = repetition;
-      ParseForParameter(verbose,file, keyword, 0, 1, 1, double_type, &neues->x.x[0], repetition,critical);
-      ParseForParameter(verbose,file, keyword, 0, 2, 1, double_type, &neues->x.x[1], repetition,critical);
-      ParseForParameter(verbose,file, keyword, 0, 3, 1, double_type, &neues->x.x[2], repetition,critical);
-      ParseForParameter(verbose,file, keyword, 0, 4, 1, int_type, &neues->FixedIon, repetition,critical);
-      if(!ParseForParameter(verbose,file, keyword, 0, 5, 1, double_type, &neues->v.x[0], repetition,optional))
-        neues->v.x[0] = 0.;
-      if(!ParseForParameter(verbose,file, keyword, 0, 6, 1, double_type, &neues->v.x[1], repetition,optional))
-        neues->v.x[1] = 0.;
-      if(!ParseForParameter(verbose,file, keyword, (int)FastParsing, 7, 1, double_type, &neues->v.x[2], repetition,optional))
-        neues->v.x[2] = 0.;
-      // here we don't care if forces are present (last in trajectories is always equal to current position)
-      neues->type = elementhash[i]; // find element type
-      mol->AddAtom(neues);
+      repetition++;
     }
-  }
-  mol->MDSteps = MinSteps;    // set equal to smallest number of parsed steps
+    repetition--;
+    cout << "Found " << repetition << " trajectory steps." << endl;
+    mol->MDSteps = repetition;
+  } else {
+    // find the maximum number of MD steps so that we may parse last one (Ion_Type1_1 must always be present, because is the first atom)
+    repetition = 0;
+    while ( ParseForParameter(verbose,file, "Ion_Type1_1", 0, 1, 1, double_type, &value[0], repetition, (repetition == 0) ? critical : optional) &&
+            ParseForParameter(verbose,file, "Ion_Type1_1", 0, 2, 1, double_type, &value[1], repetition, (repetition == 0) ? critical : optional) &&
+            ParseForParameter(verbose,file, "Ion_Type1_1", 0, 3, 1, double_type, &value[2], repetition, (repetition == 0) ? critical : optional))
+      repetition++;
+    cout << "I found " << repetition << " times the keyword Ion_Type1_1." << endl;
+    // parse in molecule coordinates
+    for (int i=0; i < config::MaxTypes; i++) {
+      sprintf(name,"Ion_Type%i",i+1);
+      for(int j=0;j<No[i];j++) {
+        sprintf(keyword,"%s_%i",name, j+1);
+        atom *neues = new atom();
+        // then parse for each atom the coordinates as often as present
+        ParseForParameter(verbose,file, keyword, 0, 1, 1, double_type, &neues->x.x[0], repetition,critical);
+        ParseForParameter(verbose,file, keyword, 0, 2, 1, double_type, &neues->x.x[1], repetition,critical);
+        ParseForParameter(verbose,file, keyword, 0, 3, 1, double_type, &neues->x.x[2], repetition,critical);
+        ParseForParameter(verbose,file, keyword, 0, 4, 1, int_type, &neues->FixedIon, repetition,critical);
+        if(!ParseForParameter(verbose,file, keyword, 0, 5, 1, double_type, &neues->v.x[0], repetition,optional))
+          neues->v.x[0] = 0.;
+        if(!ParseForParameter(verbose,file, keyword, 0, 6, 1, double_type, &neues->v.x[1], repetition,optional))
+          neues->v.x[1] = 0.;
+        if(!ParseForParameter(verbose,file, keyword, 0, 7, 1, double_type, &neues->v.x[2], repetition,optional))
+          neues->v.x[2] = 0.;
+        // here we don't care if forces are present (last in trajectories is always equal to current position)
+        neues->type = elementhash[i]; // find element type
+        mol->AddAtom(neues);
+      }
+    }
+  }
   file->close();
   delete(file);
@@ -1015,5 +1051,8 @@
     *output << endl;
     result = result && mol->Checkout(output);
-    result = result && mol->OutputTrajectories(output);
+    if (mol->MDSteps <=1 )
+      result = result && mol->Output(output);
+    else
+      result = result && mol->OutputTrajectories(output);
     return result;
   } else
@@ -1108,5 +1147,5 @@
     // ... and check if it is the keyword!
     //fprintf(stderr,"name %p, dummy %i/%c, dummy1 %i/%c, strlen(name) %i\n", &name, dummy, *dummy, dummy1, *dummy1, strlen(name));
-    if ((name == NULL) || ((dummy-dummy1 >= 3) && (strncmp(dummy1, name, strlen(name)) == 0)) && ((unsigned int)(dummy-dummy1) == strlen(name))) {
+    if ((name == NULL) || (((dummy-dummy1 >= 3) && (strncmp(dummy1, name, strlen(name)) == 0)) && ((unsigned int)(dummy-dummy1) == strlen(name)))) {
       found++; // found the parameter!
       //fprintf(stderr,"found %s at line %i between %i and %i\n", name, line, dummy1, dummy);    
@@ -1163,9 +1202,8 @@
               dummy1 = dummy;
               dummy = strchr(dummy1, '\t'); // seek for tab or space
-              if (dummy == NULL) {
+              if (dummy == NULL)
                 dummy = strchr(dummy1, ' ');  // if not found seek for space
-                while ((dummy != NULL) && ((*dummy == '\t') || (*dummy == ' ')))    // skip some more tabs and spaces if necessary
-                  dummy++;
-              }
+              while ((dummy != NULL) && ((*dummy == '\t') || (*dummy == ' ')))    // skip some more tabs and spaces if necessary
+                dummy++;
 /*              while ((dummy != NULL) && ((*dummy == '\t') || (*dummy == ' ')))    // skip some more tabs and spaces if necessary
                 dummy++;*/
@@ -1189,4 +1227,13 @@
               } else {
                 //fprintf(stderr,"found tab at %i\n",(char *)dummy-(char *)free_dummy);
+              }
+              if (*dummy1 == '#') {
+                // found comment, skipping rest of line
+                //if (verbose) fprintf(stderr,"Error: '#' at %i and still missing %i value(s) for parameter %s\n", line, yth-j, name);
+                if (!sequential) { // here we need it! 
+                  file->seekg(file_position, ios::beg);  // rewind to start position
+                }
+                Free((void **)&free_dummy, "config::ParseForParameter: *free_dummy");
+                return 0;
               }
               //fprintf(stderr,"value from %i to %i\n",(char *)dummy1-(char *)free_dummy,(char *)dummy-(char *)free_dummy);
Index: src/moleculelist.cpp
===================================================================
--- src/moleculelist.cpp	(revision f054074bc0e2ebd43408e39a6e2ae0cd96597630)
+++ src/moleculelist.cpp	(revision 362b0e5abaff6113ef0ca4989b9074448eccb8b8)
@@ -427,5 +427,5 @@
     *out << endl;
     // prepare output of config file
-    sprintf(FragmentName, "%s/%s%s.conf", PathBackup, FRAGMENTPREFIX, FragmentNumber);
+    sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
     outputFragment.open(FragmentName, ios::out);
     //strcpy(PathBackup, configuration->configpath);
@@ -739,4 +739,5 @@
 void MoleculeLeafClass::TranslateIndicesToGlobalIDs(ofstream *out, Graph **FragmentList, int &FragmentCounter, int &TotalNumberOfKeySets, Graph &TotalGraph)
 {
+  *out << Verbose(1) << "Begin of TranslateIndicesToGlobalIDs." << endl;
   KeySet *TempSet = new KeySet;
   if (FragmentList[FragmentCounter] != NULL) {
@@ -754,4 +755,5 @@
     next->TranslateIndicesToGlobalIDs(out, FragmentList, ++FragmentCounter, TotalNumberOfKeySets, TotalGraph);
   FragmentCounter--;
+  *out << Verbose(1) << "End of TranslateIndicesToGlobalIDs." << endl;
 };
 
Index: src/molecules.cpp
===================================================================
--- src/molecules.cpp	(revision f054074bc0e2ebd43408e39a6e2ae0cd96597630)
+++ src/molecules.cpp	(revision 362b0e5abaff6113ef0ca4989b9074448eccb8b8)
@@ -455,5 +455,7 @@
   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,7 +473,16 @@
       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,24 +976,22 @@
 
 /** 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,10 +1001,31 @@
     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,66 +1033,62 @@
           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,5 +1398,5 @@
               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,6 +1439,35 @@
 };
 
+/** 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,5 +2670,5 @@
         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,5 +2682,5 @@
         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,11 +2689,12 @@
           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,7 +2702,9 @@
                 //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,7 +2722,10 @@
         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,9 +2774,13 @@
       }
     }
-    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,5 +2797,5 @@
 };
 
-/** 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,5 +2846,5 @@
  * \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,4 +2856,5 @@
   fstream File;
   bool FragmentationToDo = true;
+  bool CheckOrder = false;
   Graph **FragmentList = NULL;
   Graph *ParsedFragmentList = NULL;
@@ -2837,5 +2903,8 @@
   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,52 +2958,56 @@
 
   // ===== 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,6 +3030,6 @@
     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,9 +3052,9 @@
 {
   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,11 +3064,18 @@
   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,5 +3083,6 @@
       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,4 +3094,5 @@
   }
   Free((void **)&OrderArray, "molecule::ParseOrderAtSiteFromFile - *OrderArray");
+  Free((void **)&MaxArray, "molecule::ParseOrderAtSiteFromFile - *MaxArray");
   
   *out << Verbose(1) << "End of ParseOrderAtSiteFromFile" << endl;
@@ -3602,5 +3684,4 @@
   int FragmentCounter;
   int CurrentIndex;
-  int *Labels;
   double TEFactor;
   int *ShortestPathList;
@@ -3670,6 +3751,4 @@
         	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,5 +3809,5 @@
         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,10 +3900,10 @@
 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,8 +3913,5 @@
 
   // 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,5 +3947,5 @@
       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,11 +3961,5 @@
           *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,5 +3969,10 @@
             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,9 +4078,9 @@
       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,10 +4270,8 @@
   // 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,5 +4322,5 @@
       // 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,46 +4333,50 @@
       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,5 +4397,4 @@
   // cleanup FragmentSearch structure
   Free((void **)&FragmentSearch.ShortestPathList, "molecule::PowerSetGenerator: *ShortestPathList");
-  Free((void **)&FragmentSearch.Labels, "molecule::PowerSetGenerator: *Labels");
   delete(FragmentSearch.FragmentSet);
   
Index: src/molecules.hpp
===================================================================
--- src/molecules.hpp	(revision f054074bc0e2ebd43408e39a6e2ae0cd96597630)
+++ src/molecules.hpp	(revision 362b0e5abaff6113ef0ca4989b9074448eccb8b8)
@@ -24,4 +24,5 @@
 
 #include "helpers.hpp"
+#include "parser.hpp"
 #include "periodentafel.hpp"
 #include "stackclass.hpp"
@@ -113,4 +114,5 @@
     bool IsCyclic;        //!< whether atom belong to as cycle or not, determined in DepthFirstSearchAnalysis()
     unsigned char AdaptiveOrder;  //!< current present bond order at site (0 means "not set")
+    bool MaxOrder;  //!< whether this atom as a root in fragmentation still creates more fragments on higher orders or not
   
   atom();
@@ -231,5 +233,5 @@
 	void PrincipalAxisSystem(ofstream *out, bool DoRotate);
 	double VolumeOfConvexEnvelope(ofstream *out, bool IsAngstroem);
-	bool VerletForceIntegration(char *file, double delta_t);
+	bool VerletForceIntegration(char *file, double delta_t, bool IsAngstroem);
 	
   bool CheckBounds(const Vector *x) const;
@@ -256,5 +258,5 @@
   
   /// Fragment molecule by two different approaches:
-  void FragmentMolecule(ofstream *out, int Order, config *configuration);
+  int FragmentMolecule(ofstream *out, int Order, config *configuration);
   bool CheckOrderAtSite(ofstream *out, bool *AtomMask, Graph *GlobalKeySetList, int Order, int *MinimumRingSize, char *path = NULL);
   bool StoreAdjacencyToFile(ofstream *out, char *path);
@@ -286,4 +288,5 @@
   void OutputListOfBonds(ofstream *out) const;
   bool OutputXYZ(ofstream *out) const;
+  bool OutputTrajectoriesXYZ(ofstream *out);
   bool Checkout(ofstream *out) const;
 
@@ -354,4 +357,5 @@
     char *configname;
     bool FastParsing;
+    double Deltat;
     
     private:
@@ -375,5 +379,4 @@
     
     int MaxOuterStep;
-    double Deltat;
     int OutVisStep;
     int OutSrcStep;
