Changes in molecuilder/src/molecules.cpp [d01f4d:848729]

File:

• molecuilder/src/molecules.cpp (modified) (211 diffs, 1 prop)

molecuilder/src/molecules.cpp

@@ -1 +1 @@
 /** \file molecules.cpp
- * 
+ *
  * Functions for the class molecule.
- * 
+ *
  */
 
@@ -25 +25 @@
       sum += (gsl_vector_get(x,j) - (vectors[i])->x[j])*(gsl_vector_get(x,j) - (vectors[i])->x[j]);
   }
-  
+
   return sum;
 };
@@ -34 +34 @@
  * Initialises molecule list with correctly referenced start and end, and sets molecule::last_atom to zero.
  */
-molecule::molecule(periodentafel *teil) 
-{ 
+molecule::molecule(periodentafel *teil)
+{
   // init atom chain list
-  start = new atom; 
+  start = new atom;
   end = new atom;
-  start->father = NULL; 
+  start->father = NULL;
   end->father = NULL;
   link(start,end);
@@ -46 +46 @@
   last = new bond(start, end, 1, -1);
   link(first,last);
-  // other stuff 
+  // other stuff
   MDSteps = 0;
-  last_atom = 0; 
+  last_atom = 0;
   elemente = teil;
   AtomCount = 0;
@@ -67 +67 @@
  * Initialises molecule list with correctly referenced start and end, and sets molecule::last_atom to zero.
  */
-molecule::~molecule() 
+molecule::~molecule()
 {
   if (ListOfBondsPerAtom != NULL)
@@ -78 +78 @@
   delete(last);
   delete(end);
-  delete(start); 
+  delete(start);
 };
 
 /** Adds given atom \a *pointer from molecule list.
- * Increases molecule::last_atom and gives last number to added atom and names it according to its element::abbrev and molecule::AtomCount 
+ * Increases molecule::last_atom and gives last number to added atom and names it according to its element::abbrev and molecule::AtomCount
  * \param *pointer allocated and set atom
  * \return true - succeeded, false - atom not found in list
  */
 bool molecule::AddAtom(atom *pointer)
-{ 
+{
   if (pointer != NULL) {
-    pointer->sort = &pointer->nr; 
+    pointer->sort = &pointer->nr;
     pointer->nr = last_atom++;  // increase number within molecule
     AtomCount++;
@@ -106 +106 @@
     return add(pointer, end);
   } else
-    return false; 
+    return false;
 };
 
@@ -115 +115 @@
  */
 atom * molecule::AddCopyAtom(atom *pointer)
-{ 
+{
   if (pointer != NULL) {
         atom *walker = new atom();
@@ -122 +122 @@
     walker->v.CopyVector(&pointer->v); // copy velocity
     walker->FixedIon = pointer->FixedIon;
-    walker->sort = &walker->nr; 
+    walker->sort = &walker->nr;
     walker->nr = last_atom++;  // increase number within molecule
     walker->father = pointer; //->GetTrueFather();
@@ -133 +133 @@
     return walker;
   } else
-    return NULL; 
+    return NULL;
 };
 
@@ -156 +156 @@
  *    The lengths for these are \f$f\f$ and \f$g\f$ (from triangle H1-H2-(center of H1-H2-H3)) with knowledge that
  *    the median lines in an isosceles triangle meet in the center point with a ratio 2:1.
- *    \f[ f = \frac{b}{\sqrt{3}} \qquad g = \frac{b}{2} 
+ *    \f[ f = \frac{b}{\sqrt{3}} \qquad g = \frac{b}{2}
  *    \f]
- *    as the coordination of all three atoms in the coordinate system of these three vectors: 
+ *    as the coordination of all three atoms in the coordinate system of these three vectors:
  *    \f$\pmatrix{d & f & 0}\f$, \f$\pmatrix{d & -0.5 \cdot f & g}\f$ and \f$\pmatrix{d & -0.5 \cdot f & -g}\f$.
- * 
+ *
  * \param *out output stream for debugging
- * \param *Bond pointer to bond between \a *origin and \a *replacement 
- * \param *TopOrigin son of \a *origin of upper level molecule (the atom added to this molecule as a copy of \a *origin) 
+ * \param *Bond pointer to bond between \a *origin and \a *replacement
+ * \param *TopOrigin son of \a *origin of upper level molecule (the atom added to this molecule as a copy of \a *origin)
  * \param *origin pointer to atom which acts as the origin for scaling the added hydrogen to correct bond length
  * \param *replacement pointer to the atom which shall be copied as a hydrogen atom in this molecule
@@ -191 +191 @@
   InBondvector.SubtractVector(&TopOrigin->x);
   bondlength = InBondvector.Norm();
-   
+
    // is greater than typical bond distance? Then we have to correct periodically
    // the problem is not the H being out of the box, but InBondvector have the wrong direction
-   // due to TopReplacement or Origin being on the wrong side! 
-  if (bondlength > BondDistance) { 
+   // due to TopReplacement or Origin being on the wrong side!
+  if (bondlength > BondDistance) {
 //    *out << Verbose(4) << "InBondvector is: ";
 //    InBondvector.Output(out);
@@ -215 +215 @@
 //    *out << endl;
   } // periodic correction finished
-  
+
   InBondvector.Normalize();
   // get typical bond length and store as scale factor for later
   BondRescale = TopOrigin->type->HBondDistance[TopBond->BondDegree-1];
   if (BondRescale == -1) {
-    cerr << Verbose(3) << "WARNING: There is no typical bond distance for bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl;
-    BondRescale = bondlength; 
+    cerr << Verbose(3) << "ERROR: There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl;
+    return false;
+    BondRescale = bondlength;
   } else {
     if (!IsAngstroem)
@@ -272 +273 @@
       if (FirstOtherAtom != NULL) { // then we just have this double bond and the plane does not matter at all
 //        *out << Verbose(3) << "Regarding the double bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") to be constructed: Taking " << FirstOtherAtom->Name << " and " << SecondOtherAtom->Name << " along with " << TopOrigin->Name << " to determine orthogonal plane." << endl;
-        
+
         // determine the plane of these two with the *origin
         AllWentWell = AllWentWell && Orthovector1.MakeNormalVector(&TopOrigin->x, &FirstOtherAtom->x, &SecondOtherAtom->x);
@@ -285 +286 @@
       Orthovector1.Normalize();
       //*out << Verbose(3) << "ReScaleCheck: " << Orthovector1.Norm() << " and " << InBondvector.Norm() << "." << endl;
-      
+
       // create the two Hydrogens ...
       FirstOtherAtom = new atom();
@@ -299 +300 @@
       bondangle = TopOrigin->type->HBondAngle[1];
       if (bondangle == -1) {
-        *out << Verbose(3) << "WARNING: There is no typical bond angle for bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl;
+        *out << Verbose(3) << "ERROR: There is no typical hydrogen bond angle in replacing bond (" << TopOrigin->Name << "<->" << TopReplacement->Name << ") of degree " << TopBond->BondDegree << "!" << endl;
+        return false;
         bondangle = 0;
       }
@@ -316 +318 @@
         SecondOtherAtom->x.x[i] = InBondvector.x[i] * cos(bondangle) + Orthovector1.x[i] * (-sin(bondangle));
       }
-      FirstOtherAtom->x.Scale(&BondRescale);  // rescale by correct BondDistance 
+      FirstOtherAtom->x.Scale(&BondRescale);  // rescale by correct BondDistance
       SecondOtherAtom->x.Scale(&BondRescale);
       //*out << Verbose(3) << "ReScaleCheck: " << FirstOtherAtom->x.Norm() << " and " << SecondOtherAtom->x.Norm() << "." << endl;
       for(int i=NDIM;i--;) { // and make relative to origin atom
-        FirstOtherAtom->x.x[i] += TopOrigin->x.x[i]; 
+        FirstOtherAtom->x.x[i] += TopOrigin->x.x[i];
         SecondOtherAtom->x.x[i] += TopOrigin->x.x[i];
       }
@@ -363 +365 @@
 //      *out << endl;
       AllWentWell = AllWentWell && Orthovector2.MakeNormalVector(&InBondvector, &Orthovector1);
-//      *out << Verbose(3) << "Orthovector2: "; 
+//      *out << Verbose(3) << "Orthovector2: ";
 //      Orthovector2.Output(out);
 //      *out << endl;
-      
+
       // create correct coordination for the three atoms
       alpha = (TopOrigin->type->HBondAngle[2])/180.*M_PI/2.;  // retrieve triple bond angle from database
@@ -378 +380 @@
       factors[0] = d;
       factors[1] = f;
-      factors[2] = 0.; 
+      factors[2] = 0.;
       FirstOtherAtom->x.LinearCombinationOfVectors(&InBondvector, &Orthovector1, &Orthovector2, factors);
       factors[1] = -0.5*f;
-      factors[2] = g; 
+      factors[2] = g;
       SecondOtherAtom->x.LinearCombinationOfVectors(&InBondvector, &Orthovector1, &Orthovector2, factors);
-      factors[2] = -g; 
+      factors[2] = -g;
       ThirdOtherAtom->x.LinearCombinationOfVectors(&InBondvector, &Orthovector1, &Orthovector2, factors);
 
@@ -435 +437 @@
  */
 bool molecule::AddXYZFile(string filename)
-{ 
+{
   istringstream *input = NULL;
   int NumberOfAtoms = 0; // atom number in xyz read
@@ -444 +446 @@
   string line;    // currently parsed line
   double x[3];    // atom coordinates
-  
+
   xyzfile.open(filename.c_str());
   if (!xyzfile)
@@ -452 +454 @@
   input = new istringstream(line);
   *input >> NumberOfAtoms;
-  cout << Verbose(0) << "Parsing " << NumberOfAtoms << " atoms in file." << endl; 
+  cout << Verbose(0) << "Parsing " << NumberOfAtoms << " atoms in file." << endl;
   getline(xyzfile,line,'\n'); // Read comment
   cout << Verbose(1) << "Comment: " << line << endl;
-  
+
   if (MDSteps == 0) // no atoms yet present
     MDSteps++;
@@ -489 +491 @@
   xyzfile.close();
   delete(input);
-  return true; 
+  return true;
 };
 
@@ -501 +503 @@
   atom *LeftAtom = NULL, *RightAtom = NULL;
   atom *Walker = NULL;
-  
+
   // copy all atoms
   Walker = start;
@@ -508 +510 @@
     CurrentAtom = copy->AddCopyAtom(Walker);
   }
-  
+
   // copy all bonds
   bond *Binder = first;
@@ -532 +534 @@
       copy->NoCyclicBonds++;
     NewBond->Type = Binder->Type;
-  }  
+  }
   // correct fathers
   Walker = copy->start;
@@ -549 +551 @@
     copy->CreateListOfBondsPerAtom((ofstream *)&cout);
   }
-  
+
   return copy;
 };
@@ -574 +576 @@
 
 /** Remove bond from bond chain list.
- * \todo Function not implemented yet 
+ * \todo Function not implemented yet
  * \param *pointer bond pointer
  * \return true - bound found and removed, false - bond not found/removed
@@ -586 +588 @@
 
 /** Remove every bond from bond chain list that atom \a *BondPartner is a constituent of.
- * \todo Function not implemented yet 
+ * \todo Function not implemented yet
  * \param *BondPartner atom to be removed
  * \return true - bounds found and removed, false - bonds not found/removed
@@ -619 +621 @@
   Vector *min = new Vector;
   Vector *max = new Vector;
-  
+
   // gather min and max for each axis
   ptr = start->next;  // start at first in list
@@ -665 +667 @@
 {
   Vector *min = new Vector;
-  
+
 //  *out << Verbose(3) << "Begin of CenterEdge." << endl;
   atom *ptr = start->next;  // start at first in list
@@ -681 +683 @@
       }
     }
-//    *out << Verbose(4) << "Maximum is "; 
+//    *out << Verbose(4) << "Maximum is ";
 //    max->Output(out);
 //    *out << ", Minimum is ";
@@ -689 +691 @@
     max->AddVector(min);
     Translate(min);
-  } 
+  }
   delete(min);
 //  *out << Verbose(3) << "End of CenterEdge." << endl;
-}; 
+};
 
 /** Centers the center of the atoms at (0,0,0).
@@ -702 +704 @@
   int Num = 0;
   atom *ptr = start->next;  // start at first in list
- 
+
   for(int i=NDIM;i--;) // zero center vector
     center->x[i] = 0.;
-    
+
   if (ptr != end) {   //list not empty?
     while (ptr->next != end) {  // continue with second if present
       ptr = ptr->next;
       Num++;
-      center->AddVector(&ptr->x);       
+      center->AddVector(&ptr->x);
     }
     center->Scale(-1./Num); // divide through total number (and sign for direction)
     Translate(center);
   }
-}; 
+};
 
 /** Returns vector pointing to center of gravity.
@@ -727 +729 @@
   Vector tmp;
   double Num = 0;
-  
+
   a->Zero();
 
@@ -735 +737 @@
       Num += 1.;
       tmp.CopyVector(&ptr->x);
-      a->AddVector(&tmp);       
+      a->AddVector(&tmp);
     }
     a->Scale(-1./Num); // divide through total mass (and sign for direction)
@@ -755 +757 @@
         Vector tmp;
   double Num = 0;
-        
+
         a->Zero();
 
@@ -764 +766 @@
       tmp.CopyVector(&ptr->x);
       tmp.Scale(ptr->type->mass);  // scale by mass
-      a->AddVector(&tmp);       
+      a->AddVector(&tmp);
     }
     a->Scale(-1./Num); // divide through total mass (and sign for direction)
@@ -787 +789 @@
     Translate(center);
   }
-}; 
+};
 
 /** Scales all atoms by \a *factor.
@@ -801 +803 @@
       Trajectories[ptr].R.at(j).Scale(factor);
     ptr->x.Scale(factor);
-  }      
-};
-
-/** Translate all atoms by given vector. 
+  }
+};
+
+/** Translate all atoms by given vector.
  * \param trans[] translation vector.
  */
@@ -816 +818 @@
       Trajectories[ptr].R.at(j).Translate(trans);
     ptr->x.Translate(trans);
-  }      
-};
-
-/** Mirrors all atoms against a given plane. 
+  }
+};
+
+/** Mirrors all atoms against a given plane.
  * \param n[] normal vector of mirror plane.
  */
@@ -831 +833 @@
       Trajectories[ptr].R.at(j).Mirror(n);
     ptr->x.Mirror(n);
-  }      
+  }
 };
 
@@ -845 +847 @@
   bool flag;
   Vector Testvector, Translationvector;
-  
+
   do {
     Center.Zero();
@@ -861 +863 @@
           if (Walker->nr < Binder->GetOtherAtom(Walker)->nr) // otherwise we shift one to, the other fro and gain nothing
             for (int j=0;j<NDIM;j++) {
-              tmp = Walker->x.x[j] - Binder->GetOtherAtom(Walker)->x.x[j]; 
+              tmp = Walker->x.x[j] - Binder->GetOtherAtom(Walker)->x.x[j];
               if ((fabs(tmp)) > BondDistance) {
                 flag = false;
@@ -877 +879 @@
         cout << Verbose(1) << "vector is: ";
         Testvector.Output((ofstream *)&cout);
-        cout << endl;     
+        cout << endl;
 #ifdef ADDHYDROGEN
         // now also change all hydrogens
@@ -890 +892 @@
             cout << Verbose(1) << "Hydrogen vector is: ";
             Testvector.Output((ofstream *)&cout);
-            cout << endl;     
+            cout << endl;
           }
         }
@@ -912 +914 @@
 
         CenterGravity(out, CenterOfGravity);
-        
-        // reset inertia tensor 
+
+        // reset inertia tensor
         for(int i=0;i<NDIM*NDIM;i++)
                 InertiaTensor[i] = 0.;
-        
+
         // sum up inertia tensor
         while (ptr->next != end) {
@@ -941 +943 @@
         }
         *out << endl;
-        
+
         // diagonalize to determine principal axis system
         gsl_eigen_symmv_workspace *T = gsl_eigen_symmv_alloc(NDIM);
@@ -950 +952 @@
         gsl_eigen_symmv_free(T);
         gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_ABS_DESC);
-        
+
         for(int i=0;i<NDIM;i++) {
                 *out << Verbose(1) << "eigenvalue = " << gsl_vector_get(eval, i);
                 *out << ", eigenvector = (" << evec->data[i * evec->tda + 0] << "," << evec->data[i * evec->tda + 1] << "," << evec->data[i * evec->tda + 2] << ")" << endl;
         }
-        
+
         // check whether we rotate or not
         if (DoRotate) {
-          *out << Verbose(1) << "Transforming molecule into PAS ... "; 
+          *out << Verbose(1) << "Transforming molecule into PAS ... ";
           // the eigenvectors specify the transformation matrix
           ptr = start;
@@ -970 +972 @@
 
           // summing anew for debugging (resulting matrix has to be diagonal!)
-          // reset inertia tensor 
+          // reset inertia tensor
     for(int i=0;i<NDIM*NDIM;i++)
       InertiaTensor[i] = 0.;
-    
+
     // sum up inertia tensor
     ptr = start;
@@ -1000 +1002 @@
     *out << endl;
         }
-        
+
         // free everything
         delete(CenterOfGravity);
@@ -1007 +1009 @@
 };
 
-/** Evaluates the potential energy used for constrained molecular dynamics.
- * \f$V_i^{con} = c^{bond} \cdot | r_{P(i)} - R_i | + sum_{i \neq j} C^{min} \cdot \frac{1}{C_{ij}} + C^{inj} \Bigl (1 - \theta \bigl (\prod_{i \neq j} (P(i) - P(j)) \bigr ) \Bigr )\f$
- *     where the first term points to the target in minimum distance, the second is a penalty for trajectories lying too close to each other (\f$C_{ij}$ is minimum distance between 
- *     trajectories i and j) and the third term is a penalty for two atoms trying to each the same target point.
- * Note that for the second term we have to solve the following linear system:
- * \f$-c_1 \cdot n_1 + c_2 \cdot n_2 + C \cdot n_3 = - p_2 + p_1\f$, where \f$c_1\f$, \f$c_2\f$ and \f$C\f$ are constants, 
- * offset vector \f$p_1\f$ in direction \f$n_1\f$, offset vector \f$p_2\f$ in direction \f$n_2\f$,
- * \f$n_3\f$ is the normal vector to both directions. \f$C\f$ would be the minimum distance between the two lines.
- * \sa molecule::MinimiseConstrainedPotential(), molecule::VerletForceIntegration()
- * \param *out output stream for debugging
- * \param *PermutationMap gives target ptr for each atom, array of size molecule::AtomCount (this is "x" in \f$ V^{con}(x) \f$ )
- * \param startstep start configuration (MDStep in molecule::trajectories)
- * \param endstep end configuration (MDStep in molecule::trajectories)
- * \param *constants constant in front of each term
- * \param IsAngstroem whether coordinates are in angstroem (true) or bohrradius (false)
- * \return potential energy
- * \note This routine is scaling quadratically which is not optimal.
- * \todo There's a bit double counting going on for the first time, bu nothing to worry really about.
- */
-double molecule::ConstrainedPotential(ofstream *out, atom **PermutationMap, int startstep, int endstep, double *constants, bool IsAngstroem)
-{
-  double result = 0., tmp, Norm1, Norm2;
-  atom *Walker = NULL, *Runner = NULL, *Sprinter = NULL;
-  Vector trajectory1, trajectory2, normal, TestVector;
-  gsl_matrix *A = gsl_matrix_alloc(NDIM,NDIM);
-  gsl_vector *x = gsl_vector_alloc(NDIM);
-
-  // go through every atom
-  Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    // first term: distance to target
-    Runner = PermutationMap[Walker->nr];   // find target point
-    tmp = (Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Runner].R.at(endstep)));
-    tmp *= IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
-    result += constants[0] * tmp;
-    //*out << Verbose(4) << "Adding " << tmp*constants[0] << "." << endl;
-    
-    // second term: sum of distances to other trajectories
-    Runner = start;
-    while (Runner->next != end) {
-      Runner = Runner->next;
-      if (Runner == Walker) // hence, we only go up to the Walker, not beyond (similar to i=0; i<j; i++)
-        break;
-      // determine normalized trajectories direction vector (n1, n2)
-      Sprinter = PermutationMap[Walker->nr];   // find first target point
-      trajectory1.CopyVector(&Trajectories[Sprinter].R.at(endstep));
-      trajectory1.SubtractVector(&Trajectories[Walker].R.at(startstep));
-      trajectory1.Normalize();
-      Norm1 = trajectory1.Norm();
-      Sprinter = PermutationMap[Runner->nr];   // find second target point
-      trajectory2.CopyVector(&Trajectories[Sprinter].R.at(endstep));
-      trajectory2.SubtractVector(&Trajectories[Runner].R.at(startstep));
-      trajectory2.Normalize();
-      Norm2 = trajectory1.Norm();
-      // check whether either is zero()
-      if ((Norm1 < MYEPSILON) && (Norm2 < MYEPSILON)) {
-        tmp = Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Runner].R.at(startstep));
-      } else if (Norm1 < MYEPSILON) {
-        Sprinter = PermutationMap[Walker->nr];   // find first target point
-        trajectory1.CopyVector(&Trajectories[Sprinter].R.at(endstep));  // copy first offset
-        trajectory1.SubtractVector(&Trajectories[Runner].R.at(startstep));  // subtract second offset
-        trajectory2.Scale( trajectory1.ScalarProduct(&trajectory2) ); // trajectory2 is scaled to unity, hence we don't need to divide by anything 
-        trajectory1.SubtractVector(&trajectory2);   // project the part in norm direction away
-        tmp = trajectory1.Norm();  // remaining norm is distance
-      } else if (Norm2 < MYEPSILON) {
-        Sprinter = PermutationMap[Runner->nr];   // find second target point
-        trajectory2.CopyVector(&Trajectories[Sprinter].R.at(endstep));  // copy second offset
-        trajectory2.SubtractVector(&Trajectories[Walker].R.at(startstep));  // subtract first offset
-        trajectory1.Scale( trajectory2.ScalarProduct(&trajectory1) ); // trajectory1 is scaled to unity, hence we don't need to divide by anything 
-        trajectory2.SubtractVector(&trajectory1);   // project the part in norm direction away
-        tmp = trajectory2.Norm();  // remaining norm is distance
-      } else if ((fabs(trajectory1.ScalarProduct(&trajectory2)/Norm1/Norm2) - 1.) < MYEPSILON) { // check whether they're linear dependent
-//        *out << Verbose(3) << "Both trajectories of " << *Walker << " and " << *Runner << " are linear dependent: ";
-//        *out << trajectory1;
-//        *out << " and ";
-//        *out << trajectory2;
-        tmp = Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Runner].R.at(startstep));
-//        *out << " with distance " << tmp << "." << endl;
-      } else { // determine distance by finding minimum distance
-//        *out << Verbose(3) << "Both trajectories of " << *Walker << " and " << *Runner << " are linear independent ";
-//        *out << endl;
-//        *out << "First Trajectory: ";
-//        *out << trajectory1 << endl;
-//        *out << "Second Trajectory: ";
-//        *out << trajectory2 << endl;
-        // determine normal vector for both
-        normal.MakeNormalVector(&trajectory1, &trajectory2);
-        // print all vectors for debugging
-//        *out << "Normal vector in between: ";
-//        *out << normal << endl;
-        // setup matrix
-        for (int i=NDIM;i--;) {
-          gsl_matrix_set(A, 0, i, trajectory1.x[i]);
-          gsl_matrix_set(A, 1, i, trajectory2.x[i]);
-          gsl_matrix_set(A, 2, i, normal.x[i]);
-          gsl_vector_set(x,i, (Trajectories[Walker].R.at(startstep).x[i] - Trajectories[Runner].R.at(startstep).x[i]));
-        }
-        // solve the linear system by Householder transformations
-        gsl_linalg_HH_svx(A, x);
-        // distance from last component
-        tmp = gsl_vector_get(x,2);
-//        *out << " with distance " << tmp << "." << endl;
-        // test whether we really have the intersection (by checking on c_1 and c_2)
-        TestVector.CopyVector(&Trajectories[Runner].R.at(startstep));
-        trajectory2.Scale(gsl_vector_get(x,1));
-        TestVector.AddVector(&trajectory2);
-        normal.Scale(gsl_vector_get(x,2));
-        TestVector.AddVector(&normal);
-        TestVector.SubtractVector(&Trajectories[Walker].R.at(startstep));
-        trajectory1.Scale(gsl_vector_get(x,0));
-        TestVector.SubtractVector(&trajectory1);
-        if (TestVector.Norm() < MYEPSILON) {
-//          *out << Verbose(2) << "Test: ok.\tDistance of " << tmp << " is correct." << endl;  
-        } else {
-//          *out << Verbose(2) << "Test: failed.\tIntersection is off by ";
-//          *out << TestVector;
-//          *out << "." << endl;  
-        }
-      }
-      // add up
-      tmp *= IsAngstroem ? 1. : 1./AtomicLengthToAngstroem;
-      if (fabs(tmp) > MYEPSILON) {
-        result += constants[1] * 1./tmp;
-        //*out << Verbose(4) << "Adding " << 1./tmp*constants[1] << "." << endl;
-      }
-    }
-      
-    // third term: penalty for equal targets
-    Runner = start;
-    while (Runner->next != end) {
-      Runner = Runner->next;
-      if ((PermutationMap[Walker->nr] == PermutationMap[Runner->nr]) && (Walker->nr < Runner->nr)) {
-        Sprinter = PermutationMap[Walker->nr];
-//        *out << *Walker << " and " << *Runner << " are heading to the same target at ";
-//        *out << Trajectories[Sprinter].R.at(endstep);
-//        *out << ", penalting." << endl;
-        result += constants[2];
-        //*out << Verbose(4) << "Adding " << constants[2] << "." << endl;
-      }
-    }
-  }
-  
-  return result;
-};
-
-void PrintPermutationMap(ofstream *out, atom **PermutationMap, int Nr) 
-{
-  stringstream zeile1, zeile2;
-  int *DoubleList = (int *) Malloc(Nr*sizeof(int), "PrintPermutationMap: *DoubleList");
-  int doubles = 0;
-  for (int i=0;i<Nr;i++)
-    DoubleList[i] = 0;
-  zeile1 << "PermutationMap: ";
-  zeile2 << "                ";
-  for (int i=0;i<Nr;i++) {
-    DoubleList[PermutationMap[i]->nr]++;
-    zeile1 << i << " ";
-    zeile2 << PermutationMap[i]->nr << " ";
-  }
-  for (int i=0;i<Nr;i++)
-    if (DoubleList[i] > 1)
-    doubles++;
- // *out << "Found " << doubles << " Doubles." << endl;
-  Free((void **)&DoubleList, "PrintPermutationMap: *DoubleList");
-//  *out << zeile1.str() << endl << zeile2.str() << endl;
-};
-
-/** Minimises the extra potential for constrained molecular dynamics and gives forces and the constrained potential energy.
- * We do the following:
- *  -# Generate a distance list from all source to all target points
- *  -# Sort this per source point
- *  -# Take for each source point the target point with minimum distance, use this as initial permutation
- *  -# check whether molecule::ConstrainedPotential() is greater than injective penalty
- *     -# If so, we go through each source point, stepping down in the sorted target point distance list and re-checking potential.
- *  -# Next, we only apply transformations that keep the injectivity of the permutations list.
- *  -# Hence, for one source point we step down the ladder and seek the corresponding owner of this new target 
- *     point and try to change it for one with lesser distance, or for the next one with greater distance, but only
- *     if this decreases the conditional potential.
- *  -# finished. 
- *  -# Then, we calculate the forces by taking the spatial derivative, where we scale the potential to such a degree,
- *     that the total force is always pointing in direction of the constraint force (ensuring that we move in the
- *     right direction).
- *  -# Finally, we calculate the potential energy and return.
- * \param *out output stream for debugging
- * \param **PermutationMap on return: mapping between the atom label of the initial and the final configuration
- * \param startstep current MD step giving initial position between which and \a endstep we perform the constrained MD (as further steps are always concatenated)
- * \param endstep step giving final position in constrained MD 
- * \param IsAngstroem whether coordinates are in angstroem (true) or bohrradius (false)
- * \sa molecule::VerletForceIntegration()
- * \return potential energy (and allocated **PermutationMap (array of molecule::AtomCount ^2)
- * \todo The constrained potential's constants are set to fixed values right now, but they should scale based on checks of the system in order
- *       to ensure they're properties (e.g. constants[2] always greater than the energy of the system).
- * \bug this all is not O(N log N) but O(N^2)
- */
-double molecule::MinimiseConstrainedPotential(ofstream *out, atom **&PermutationMap, int startstep, int endstep, bool IsAngstroem)
-{
-  double Potential, OldPotential, OlderPotential;
-  PermutationMap = (atom **) Malloc(AtomCount*sizeof(atom *), "molecule::MinimiseConstrainedPotential: **PermutationMap");
-  DistanceMap **DistanceList = (DistanceMap **) Malloc(AtomCount*sizeof(DistanceMap *), "molecule::MinimiseConstrainedPotential: **DistanceList");
-  DistanceMap::iterator *DistanceIterators = (DistanceMap::iterator *) Malloc(AtomCount*sizeof(DistanceMap::iterator), "molecule::MinimiseConstrainedPotential: *DistanceIterators");
-  int *DoubleList = (int *) Malloc(AtomCount*sizeof(int), "molecule::MinimiseConstrainedPotential: *DoubleList");
-  DistanceMap::iterator *StepList = (DistanceMap::iterator *) Malloc(AtomCount*sizeof(DistanceMap::iterator), "molecule::MinimiseConstrainedPotential: *StepList");
-  double constants[3];
-  int round;
-  atom *Walker = NULL, *Runner = NULL, *Sprinter = NULL;
-  DistanceMap::iterator Rider, Strider;
-  
-  /// Minimise the potential
-  // set Lagrange multiplier constants
-  constants[0] = 10.;
-  constants[1] = 1.;
-  constants[2] = 1e+7;    // just a huge penalty
-  // generate the distance list
-  *out << Verbose(1) << "Creating the distance list ... " << endl;
-  for (int i=AtomCount; i--;) {
-    DoubleList[i] = 0;  // stores for how many atoms in startstep this atom is a target in endstep
-    DistanceList[i] = new DistanceMap;    // is the distance sorted target list per atom
-    DistanceList[i]->clear();
-  }
-  *out << Verbose(1) << "Filling the distance list ... " << endl;
-  Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    Runner = start;
-    while(Runner->next != end) {
-      Runner = Runner->next;
-      DistanceList[Walker->nr]->insert( DistancePair(Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Runner].R.at(endstep)), Runner) );
-    }
-  }
-  // create the initial PermutationMap (source -> target)
-  Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    StepList[Walker->nr] = DistanceList[Walker->nr]->begin();    // stores the step to the next iterator that could be a possible next target
-    PermutationMap[Walker->nr] = DistanceList[Walker->nr]->begin()->second;   // always pick target with the smallest distance
-    DoubleList[DistanceList[Walker->nr]->begin()->second->nr]++;            // increase this target's source count (>1? not injective)
-    DistanceIterators[Walker->nr] = DistanceList[Walker->nr]->begin();    // and remember which one we picked
-    *out << *Walker << " starts with distance " << DistanceList[Walker->nr]->begin()->first << "." << endl;
-  }
-  *out << Verbose(1) << "done." << endl;
-  // make the PermutationMap injective by checking whether we have a non-zero constants[2] term in it
-  *out << Verbose(1) << "Making the PermutationMap injective ... " << endl;
-  Walker = start;
-  DistanceMap::iterator NewBase;
-  OldPotential = fabs(ConstrainedPotential(out, PermutationMap, startstep, endstep, constants, IsAngstroem));
-  while ((OldPotential) > constants[2]) {
-    PrintPermutationMap(out, PermutationMap, AtomCount);
-    Walker = Walker->next;
-    if (Walker == end) // round-robin at the end
-      Walker = start->next;
-    if (DoubleList[DistanceIterators[Walker->nr]->second->nr] <= 1)  // no need to make those injective that aren't
-      continue;
-    // now, try finding a new one
-    NewBase = DistanceIterators[Walker->nr];  // store old base
-    do {
-      NewBase++;  // take next further distance in distance to targets list that's a target of no one
-    } while ((DoubleList[NewBase->second->nr] != 0) && (NewBase != DistanceList[Walker->nr]->end()));
-    if (NewBase != DistanceList[Walker->nr]->end()) {
-      PermutationMap[Walker->nr] = NewBase->second;
-      Potential = fabs(ConstrainedPotential(out, PermutationMap, startstep, endstep, constants, IsAngstroem));
-      if (Potential > OldPotential) { // undo
-        PermutationMap[Walker->nr] = DistanceIterators[Walker->nr]->second;
-      } else {  // do
-        DoubleList[DistanceIterators[Walker->nr]->second->nr]--;  // decrease the old entry in the doubles list
-        DoubleList[NewBase->second->nr]++;    // increase the old entry in the doubles list
-        DistanceIterators[Walker->nr] = NewBase;
-        OldPotential = Potential;
-        *out << Verbose(3) << "Found a new permutation, new potential is " << OldPotential << "." << endl;
-      }
-    }
-  } 
-  for (int i=AtomCount; i--;) // now each single entry in the DoubleList should be <=1
-    if (DoubleList[i] > 1) { 
-      cerr << "Failed to create an injective PermutationMap!" << endl;
-      exit(1);
-    }
-  *out << Verbose(1) << "done." << endl;
-  Free((void **)&DoubleList, "molecule::MinimiseConstrainedPotential: *DoubleList");
-  // argument minimise the constrained potential in this injective PermutationMap
-  *out << Verbose(1) << "Argument minimising the PermutationMap, at current potential " << OldPotential << " ... " << endl;
-  OldPotential = 1e+10;
-  round = 0;
-  do {
-    *out << "Starting round " << ++round << " ... " << endl;
-    OlderPotential = OldPotential;
-    do {
-      Walker = start;
-      while (Walker->next != end) { // pick one
-        Walker = Walker->next;
-        PrintPermutationMap(out, PermutationMap, AtomCount);
-        Sprinter = DistanceIterators[Walker->nr]->second;   // store initial partner
-        Strider = DistanceIterators[Walker->nr];  //remember old iterator
-        DistanceIterators[Walker->nr] = StepList[Walker->nr];  
-        if (DistanceIterators[Walker->nr] == DistanceList[Walker->nr]->end()) {// stop, before we run through the list and still on
-          DistanceIterators[Walker->nr] == DistanceList[Walker->nr]->begin();
-          break;
-        }
-        //*out << Verbose(2) << "Current Walker: " << *Walker << " with old/next candidate " << *Sprinter << "/" << *DistanceIterators[Walker->nr]->second << "." << endl;
-        // find source of the new target
-        Runner = start->next;
-        while(Runner != end) { // find the source whose toes we might be stepping on (Walker's new target should be in use by another already)
-          if (PermutationMap[Runner->nr] == DistanceIterators[Walker->nr]->second) {
-            //*out << Verbose(2) << "Found the corresponding owner " << *Runner << " to " << *PermutationMap[Runner->nr] << "." << endl;
-            break;
-          }
-          Runner = Runner->next;
-        }
-        if (Runner != end) { // we found the other source
-          // then look in its distance list for Sprinter
-          Rider = DistanceList[Runner->nr]->begin();
-          for (; Rider != DistanceList[Runner->nr]->end(); Rider++)
-            if (Rider->second == Sprinter)
-              break;
-          if (Rider != DistanceList[Runner->nr]->end()) { // if we have found one
-            //*out << Verbose(2) << "Current Other: " << *Runner << " with old/next candidate " << *PermutationMap[Runner->nr] << "/" << *Rider->second << "." << endl; 
-            // exchange both 
-            PermutationMap[Walker->nr] = DistanceIterators[Walker->nr]->second; // put next farther distance into PermutationMap 
-            PermutationMap[Runner->nr] = Sprinter;  // and hand the old target to its respective owner
-            PrintPermutationMap(out, PermutationMap, AtomCount);
-            // calculate the new potential
-            //*out << Verbose(2) << "Checking new potential ..." << endl;
-            Potential = ConstrainedPotential(out, PermutationMap, startstep, endstep, constants, IsAngstroem);
-            if (Potential > OldPotential) { // we made everything worse! Undo ...
-              //*out << Verbose(3) << "Nay, made the potential worse: " << Potential << " vs. " << OldPotential << "!" << endl;
-              //*out << Verbose(3) << "Setting " << *Runner << "'s source to " << *DistanceIterators[Runner->nr]->second << "." << endl; 
-              // Undo for Runner (note, we haven't moved the iteration yet, we may use this)
-              PermutationMap[Runner->nr] = DistanceIterators[Runner->nr]->second;
-              // Undo for Walker
-              DistanceIterators[Walker->nr] = Strider;  // take next farther distance target
-              //*out << Verbose(3) << "Setting " << *Walker << "'s source to " << *DistanceIterators[Walker->nr]->second << "." << endl; 
-              PermutationMap[Walker->nr] = DistanceIterators[Walker->nr]->second;
-            } else {
-              DistanceIterators[Runner->nr] = Rider;  // if successful also move the pointer in the iterator list
-              *out << Verbose(3) << "Found a better permutation, new potential is " << Potential << " vs." << OldPotential << "." << endl;
-              OldPotential = Potential;
-            }
-            if (Potential > constants[2]) {
-              cerr << "ERROR: The two-step permutation procedure did not maintain injectivity!" << endl;
-              exit(255);
-            }
-            //*out << endl;
-          } else {
-            cerr << "ERROR: " << *Runner << " was not the owner of " << *Sprinter << "!" << endl;
-            exit(255);
-          }
-        } else {
-          PermutationMap[Walker->nr] = DistanceIterators[Walker->nr]->second; // new target has no source!
-        }
-        StepList[Walker->nr]++; // take next farther distance target
-      }
-    } while (Walker->next != end);
-  } while ((OlderPotential - OldPotential) > 1e-3);
-  *out << Verbose(1) << "done." << endl;
-
-  
-  /// free memory and return with evaluated potential 
-  for (int i=AtomCount; i--;)
-    DistanceList[i]->clear();
-  Free((void **)&DistanceList, "molecule::MinimiseConstrainedPotential: **DistanceList");
-  Free((void **)&DistanceIterators, "molecule::MinimiseConstrainedPotential: *DistanceIterators");
-  return ConstrainedPotential(out, PermutationMap, startstep, endstep, constants, IsAngstroem);
-};
-
-/** Evaluates the (distance-related part) of the constrained potential for the constrained forces.
- * \param *out output stream for debugging
- * \param startstep current MD step giving initial position between which and \a endstep we perform the constrained MD (as further steps are always concatenated)
- * \param endstep step giving final position in constrained MD 
- * \param **PermutationMap mapping between the atom label of the initial and the final configuration
- * \param *Force ForceMatrix containing force vectors from the external energy functional minimisation.
- * \todo the constant for the constrained potential distance part is hard-coded independently of the hard-coded value in MinimiseConstrainedPotential()
- */
-void molecule::EvaluateConstrainedForces(ofstream *out, int startstep, int endstep, atom **PermutationMap, ForceMatrix *Force)
-{
-  double constant = 10.;
-  atom *Walker = NULL, *Sprinter = NULL; 
-  
-  /// evaluate forces (only the distance to target dependent part) with the final PermutationMap
-  *out << Verbose(1) << "Calculating forces and adding onto ForceMatrix ... " << endl;
-  Walker = start;
-  while (Walker->next != NULL) {
-    Walker = Walker->next;
-    Sprinter = PermutationMap[Walker->nr];
-    // set forces
-    for (int i=NDIM;i++;)
-      Force->Matrix[0][Walker->nr][5+i] += 2.*constant*sqrt(Trajectories[Walker].R.at(startstep).Distance(&Trajectories[Sprinter].R.at(endstep)));
-  }  
-  *out << Verbose(1) << "done." << endl;
-};
-
-/** Performs a linear interpolation between two desired atomic configurations with a given number of steps.
- * Note, step number is config::MaxOuterStep
- * \param *out output stream for debugging
- * \param startstep stating initial configuration in molecule::Trajectories
- * \param endstep stating final configuration in molecule::Trajectories
- * \param &config configuration structure
- * \return true - success in writing step files, false - error writing files or only one step in molecule::Trajectories
- */
-bool molecule::LinearInterpolationBetweenConfiguration(ofstream *out, int startstep, int endstep, const char *prefix, config &configuration) 
-{
-  bool status = true;
-  int MaxSteps = configuration.MaxOuterStep; 
-  MoleculeListClass *MoleculePerStep = new MoleculeListClass(MaxSteps+1, AtomCount);
-  // Get the Permutation Map by MinimiseConstrainedPotential
-  atom **PermutationMap = NULL;
-  atom *Walker = NULL, *Sprinter = NULL;
-  MinimiseConstrainedPotential(out, PermutationMap, startstep, endstep, configuration.GetIsAngstroem());
-
-  // check whether we have sufficient space in Trajectories for each atom
-  Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    if (Trajectories[Walker].R.size() <= (unsigned int)(MaxSteps)) {
-      //cout << "Increasing size for trajectory array of " << keyword << " to " << (MaxSteps+1) << "." << endl;
-      Trajectories[Walker].R.resize(MaxSteps+1);
-      Trajectories[Walker].U.resize(MaxSteps+1);
-      Trajectories[Walker].F.resize(MaxSteps+1);
-    }
-  }
-  // push endstep to last one
-  Walker = start;
-  while (Walker->next != end) { // remove the endstep (is now the very last one)
-    Walker = Walker->next;
-    for (int n=NDIM;n--;) {
-      Trajectories[Walker].R.at(MaxSteps).x[n] = Trajectories[Walker].R.at(endstep).x[n];
-      Trajectories[Walker].U.at(MaxSteps).x[n] = Trajectories[Walker].U.at(endstep).x[n];
-      Trajectories[Walker].F.at(MaxSteps).x[n] = Trajectories[Walker].F.at(endstep).x[n];
-    }
-  }
-  endstep = MaxSteps;
-  
-  // go through all steps and add the molecular configuration to the list and to the Trajectories of \a this molecule
-  *out << Verbose(1) << "Filling intermediate " << MaxSteps << " steps with MDSteps of " << MDSteps << "." << endl;
-  for (int step = 0; step <= MaxSteps; step++) {
-    MoleculePerStep->ListOfMolecules[step] = new molecule(elemente);
-    Walker = start;
-    while (Walker->next != end) {
-      Walker = Walker->next;
-      // add to molecule list
-      Sprinter = MoleculePerStep->ListOfMolecules[step]->AddCopyAtom(Walker);
-      for (int n=NDIM;n--;) {
-        Sprinter->x.x[n] = Trajectories[Walker].R.at(startstep).x[n] + (Trajectories[PermutationMap[Walker->nr]].R.at(endstep).x[n] - Trajectories[Walker].R.at(startstep).x[n])*((double)step/(double)MaxSteps);
-        // add to Trajectories
-        //*out << Verbose(3) << step << ">=" << MDSteps-1 << endl;
-        if (step < MaxSteps) {
-          Trajectories[Walker].R.at(step).x[n] = Trajectories[Walker].R.at(startstep).x[n] + (Trajectories[PermutationMap[Walker->nr]].R.at(endstep).x[n] - Trajectories[Walker].R.at(startstep).x[n])*((double)step/(double)MaxSteps);
-          Trajectories[Walker].U.at(step).x[n] = 0.;
-          Trajectories[Walker].F.at(step).x[n] = 0.;
-        }
-      }
-    }
-  }
-  MDSteps = MaxSteps+1;   // otherwise new Trajectories' points aren't stored on save&exit
-  
-  // store the list to single step files
-  int *SortIndex = (int *) Malloc(AtomCount*sizeof(int), "molecule::LinearInterpolationBetweenConfiguration: *SortIndex");
-  for (int i=AtomCount; i--; )
-    SortIndex[i] = i;
-  status = MoleculePerStep->OutputConfigForListOfFragments(out, "ConstrainedStep", &configuration, SortIndex, false, false);
-  
-  // free and return
-  Free((void **)&PermutationMap, "molecule::MinimiseConstrainedPotential: *PermutationMap");
-  delete(MoleculePerStep);
-  return status;
-};
-
 /** 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 *out output stream for debugging
  * \param *file filename
- * \param config structure with config::Deltat, config::IsAngstroem, config::DoConstrained
  * \param delta_t time step width in atomic units
  * \param IsAngstroem whether coordinates are in angstroem (true) or bohrradius (false)
- * \param DoConstrained whether we perform a constrained (>0, target step in molecule::trajectories) or unconstrained (0) molecular dynamics, \sa molecule::MinimiseConstrainedPotential() 
  * \return true - file found and parsed, false - file not found or imparsable
- * \todo This is not yet checked if it is correctly working with DoConstrained set to true.
- */
-bool molecule::VerletForceIntegration(ofstream *out, char *file, config &configuration)
-{
+ */
+bool molecule::VerletForceIntegration(char *file, double delta_t, bool IsAngstroem)
+{
+  element *runner = elemente->start;
   atom *walker = NULL;
   int AtomNo;
@@ -1493 +1026 @@
   string token;
   stringstream item;
-  double a, IonMass, Vector[NDIM], ConstrainedPotentialEnergy, ActualTemp;
+  double a, IonMass;
   ForceMatrix Force;
+  Vector tmpvector;
 
   CountElements();  // make sure ElementsInMolecule is up to date
-  
+
   // check file
   if (input == NULL) {
@@ -1512 +1046 @@
     }
     // 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;
-      }
-    // solve a constrained potential if we are meant to
-    if (configuration.DoConstrainedMD) {
-      // calculate forces and potential
-      atom **PermutationMap = NULL;
-      ConstrainedPotentialEnergy = MinimiseConstrainedPotential(out, PermutationMap,configuration.DoConstrainedMD, 0, configuration.GetIsAngstroem());
-      EvaluateConstrainedForces(out, configuration.DoConstrainedMD, 0, PermutationMap, &Force);
-      Free((void **)&PermutationMap, "molecule::MinimiseConstrainedPotential: *PermutationMap");
-    }
-    
+//    for(int d=0;d<NDIM;d++)
+//      tmpvector.x[d] = 0.;
+//    for(int i=0;i<AtomCount;i++)
+//      for(int d=0;d<NDIM;d++) {
+//        tmpvector.x[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] -= tmpvector.x[d]/(double)AtomCount;
+//      }
     // and perform Verlet integration for each atom with position, velocity and force vector
-    walker = start;
-    while (walker->next != end) { // go through every atom of this element
-      walker = walker->next;
-      //a = configuration.Deltat*0.5/walker->type->mass;        // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
-      // check size of vectors
-      if (Trajectories[walker].R.size() <= (unsigned int)(MDSteps)) {
-        //out << "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);
-      }
-      
-      // Update R (and F)
-      for (int d=0; d<NDIM; d++) {
-        Trajectories[walker].F.at(MDSteps).x[d] = -Force.Matrix[0][walker->nr][d+5]*(configuration.GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
-        Trajectories[walker].R.at(MDSteps).x[d] = Trajectories[walker].R.at(MDSteps-1).x[d];
-        Trajectories[walker].R.at(MDSteps).x[d] += configuration.Deltat*(Trajectories[walker].U.at(MDSteps-1).x[d]);     // s(t) = s(0) + v * deltat + 1/2 a * deltat^2 
-        Trajectories[walker].R.at(MDSteps).x[d] += 0.5*configuration.Deltat*configuration.Deltat*(Trajectories[walker].F.at(MDSteps).x[d]/walker->type->mass);     // F = m * a and s = 0.5 * F/m * t^2 = F * a * t
-      }
-      // 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] += configuration.Deltat * (Trajectories[walker].F.at(MDSteps).x[d]+Trajectories[walker].F.at(MDSteps-1).x[d]/walker->type->mass); // v = F/m * t
-      }
-//      out << "Integrated position&velocity of step " << (MDSteps) << ": ("; 
-//      for (int d=0;d<NDIM;d++)
-//        out << Trajectories[walker].R.at(MDSteps).x[d] << " ";          // next step
-//      out << ")\t(";
-//      for (int d=0;d<NDIM;d++)
-//        cout << Trajectories[walker].U.at(MDSteps).x[d] << " ";          // next step
-//      out << ")" << endl;
-            // next atom
-    }
-  }
-  // 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;
-    }
-  }
-  for(int d=0;d<NDIM;d++)
-    Vector[d] /= IonMass;
-  ActualTemp = 0.;
-  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];
-      ActualTemp += 0.5 * walker->type->mass * Trajectories[walker].U.at(MDSteps).x[d] * Trajectories[walker].U.at(MDSteps).x[d];
-    }
-  }
-  Thermostats(configuration, ActualTemp, Berendsen);
-  MDSteps++;
-
-  
-  // exit
-  return true;
-};
-
-
-/** Implementation of various thermostats.
- * All these thermostats apply an additional force which has the following forms:
- * -# Woodcock
- *  \f$p_i \rightarrow \sqrt{\frac{T_0}{T}} \cdot p_i\f$
- * -# Gaussian
- *  \f$ \frac{ \sum_i \frac{p_i}{m_i} \frac{\partial V}{\partial q_i}} {\sum_i \frac{p^2_i}{m_i}} \cdot p_i\f$
- * -# Langevin
- *  \f$p_{i,n} \rightarrow \sqrt{1-\alpha^2} p_{i,0} + \alpha p_r\f$  
- * -# Berendsen
- *  \f$p_i \rightarrow \left [ 1+ \frac{\delta t}{\tau_T} \left ( \frac{T_0}{T} \right ) \right ]^{\frac{1}{2}} \cdot p_i\f$
- * -# Nose-Hoover
- *  \f$\zeta p_i \f$ with \f$\frac{\partial \zeta}{\partial t} = \frac{1}{M_s} \left ( \sum^N_{i=1} \frac{p_i^2}{m_i} - g k_B T \right )\f$
- * These Thermostats either simply rescale the velocities, thus this function should be called after ion velocities have been updated, and/or
- * have a constraint force acting additionally on the ions. In the latter case, the ion speeds have to be modified
- * belatedly and the constraint force set.
- * \param *P Problem at hand
- * \param i which of the thermostats to take: 0 - none, 1 - Woodcock, 2 - Gaussian, 3 - Langevin, 4 - Berendsen, 5 - Nose-Hoover
- * \sa InitThermostat()
- */
-void molecule::Thermostats(config &configuration, double ActualTemp, int Thermostat)
-{
-  double ekin = 0.;
-  double E = 0., G = 0.;
-  double delta_alpha = 0.;
-  double ScaleTempFactor;
-  double sigma;
-  double IonMass;
-  int d;
-  gsl_rng * r;
-  const gsl_rng_type * T;
-  double *U = NULL, *F = NULL, FConstraint[NDIM];
-  atom *walker = NULL;
-  
-  // calculate scale configuration
-  ScaleTempFactor = configuration.TargetTemp/ActualTemp;
-    
-  // differentating between the various thermostats
-  switch(Thermostat) {
-     case None:
-      cout << Verbose(2) <<  "Applying no thermostat..." << endl;
-      break;
-     case Woodcock:
-      if ((configuration.ScaleTempStep > 0) && ((MDSteps-1) % configuration.ScaleTempStep == 0)) {
-        cout << Verbose(2) <<  "Applying Woodcock thermostat..." << endl;
+    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
+        AtomNo = 0;
         walker = start;
         while (walker->next != end) { // go through every atom of this element
           walker = walker->next;
-          IonMass = walker->type->mass;
-          U = Trajectories[walker].U.at(MDSteps).x;
-          if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
-            for (d=0; d<NDIM; d++) {
-              U[d] *= sqrt(ScaleTempFactor);
-              ekin += 0.5*IonMass * U[d]*U[d];
+          if (walker->type == runner) { // if this atom fits to element
+            // 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);
             }
-        }
-      }
-      break;
-     case Gaussian:
-      cout << Verbose(2) <<  "Applying Gaussian thermostat..." << endl;
-      walker = start;
-      while (walker->next != end) { // go through every atom of this element
-        walker = walker->next;
-        IonMass = walker->type->mass;
-        U = Trajectories[walker].U.at(MDSteps).x;
-        F = Trajectories[walker].F.at(MDSteps).x;
-        if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
-          for (d=0; d<NDIM; d++) {
-            G += U[d] * F[d]; 
-            E += U[d]*U[d]*IonMass;
-          }
-      }
-      cout << Verbose(1) << "Gaussian Least Constraint constant is " << G/E << "." << endl;
-      walker = start;
-      while (walker->next != end) { // go through every atom of this element
-        walker = walker->next;
-        IonMass = walker->type->mass;
-        U = Trajectories[walker].U.at(MDSteps).x;
-        F = Trajectories[walker].F.at(MDSteps).x;
-        if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
-          for (d=0; d<NDIM; d++) {
-            FConstraint[d] = (G/E) * (U[d]*IonMass);
-            U[d] += configuration.Deltat/IonMass * (FConstraint[d]);
-            ekin += IonMass * U[d]*U[d];
-          }
-      }
-      break;
-     case Langevin:
-      cout << Verbose(2) <<  "Applying Langevin thermostat..." << endl;
-      // init random number generator
-      gsl_rng_env_setup();
-      T = gsl_rng_default;
-      r = gsl_rng_alloc (T);
-      // Go through each ion
-      walker = start;
-      while (walker->next != end) { // go through every atom of this element
-        walker = walker->next;
-        IonMass = walker->type->mass;
-        sigma  = sqrt(configuration.TargetTemp/IonMass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)
-        U = Trajectories[walker].U.at(MDSteps).x;
-        F = Trajectories[walker].F.at(MDSteps).x;
-        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
-          // throw a dice to determine whether it gets hit by a heat bath particle
-          if (((((rand()/(double)RAND_MAX))*configuration.TempFrequency) < 1.)) {  
-            cout << Verbose(3) << "Particle " << *walker << " was hit (sigma " << sigma << "): " << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << " -> ";
-            // pick three random numbers from a Boltzmann distribution around the desired temperature T for each momenta axis
-            for (d=0; d<NDIM; d++) {
-              U[d] = gsl_ran_gaussian (r, sigma);
+            // 1. calculate x(t+\delta t)
+            for (int d=0; d<NDIM; d++) {
+              Trajectories[walker].F.at(MDSteps).x[d] = -Force.Matrix[0][AtomNo][d+5];
+              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] += 0.5*delta_t*delta_t*(Trajectories[walker].F.at(MDSteps-1).x[d])/IonMass;     // F = m * a and s = 0.5 * F/m * t^2 = F * a * t
             }
-            cout << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << endl;
-          }
-          for (d=0; d<NDIM; d++)
-            ekin += 0.5*IonMass * U[d]*U[d];
-        }
-      }
-      break;
-     case Berendsen:
-      cout << Verbose(2) <<  "Applying Berendsen-VanGunsteren thermostat..." << endl;
-      walker = start;
-      while (walker->next != end) { // go through every atom of this element
-        walker = walker->next;
-        IonMass = walker->type->mass;
-        U = Trajectories[walker].U.at(MDSteps).x;
-        F = Trajectories[walker].F.at(MDSteps).x;
-        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
-          for (d=0; d<NDIM; d++) {
-            U[d] *= sqrt(1+(configuration.Deltat/configuration.TempFrequency)*(ScaleTempFactor-1));
-            ekin += 0.5*IonMass * U[d]*U[d];
+            // 2. Calculate v(t+\delta t)
+            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] += 0.5*delta_t*(Trajectories[walker].F.at(MDSteps-1).x[d]+Trajectories[walker].F.at(MDSteps).x[d])/IonMass;
+            }
+//            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++;
           }
         }
       }
-      break;
-     case NoseHoover:
-      cout << Verbose(2) <<  "Applying Nose-Hoover thermostat..." << endl;
-      // dynamically evolve alpha (the additional degree of freedom)
-      delta_alpha = 0.;
-      walker = start;
-      while (walker->next != end) { // go through every atom of this element
-        walker = walker->next;
-        IonMass = walker->type->mass;
-        U = Trajectories[walker].U.at(MDSteps).x;
-        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
-          for (d=0; d<NDIM; d++) {
-            delta_alpha += U[d]*U[d]*IonMass;
-          }
-        }
-      }
-      delta_alpha = (delta_alpha - (3.*AtomCount+1.) * configuration.TargetTemp)/(configuration.HooverMass*Units2Electronmass);
-      configuration.alpha += delta_alpha*configuration.Deltat;
-      cout << Verbose(3) << "alpha = " << delta_alpha << " * " << configuration.Deltat << " = " << configuration.alpha << "." << endl;
-      // apply updated alpha as additional force
-      walker = start;
-      while (walker->next != end) { // go through every atom of this element
-        walker = walker->next;
-        IonMass = walker->type->mass;
-        U = Trajectories[walker].U.at(MDSteps).x;
-        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
-          for (d=0; d<NDIM; d++) {
-              FConstraint[d] = - configuration.alpha * (U[d] * IonMass);
-              U[d] += configuration.Deltat/IonMass * (FConstraint[d]);
-              ekin += (0.5*IonMass) * U[d]*U[d];
-            }
-        }
-      }
-      break;
-  }    
-  cout << Verbose(1) << "Kinetic energy is " << ekin << "." << endl;
-};
-
-/** Align all atoms in such a manner that given vector \a *n is along z axis. 
+    }
+  }
+//  // correct velocities (rather momenta) so that center of mass remains motionless
+//  tmpvector.zero()
+//  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++) {
+//      tmpvector.x[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] -= tmpvector.x[d]*walker->type->mass/IonMass;
+//    }
+//  }
+  MDSteps++;
+
+
+  // exit
+  return true;
+};
+
+/** Align all atoms in such a manner that given vector \a *n is along z axis.
  * \param n[] alignment vector.
  */
@@ -1785 +1145 @@
     ptr = ptr->next;
     tmp = ptr->x.x[0];
-    ptr->x.x[0] =  cos(alpha) * tmp + sin(alpha) * ptr->x.x[2]; 
+    ptr->x.x[0] =  cos(alpha) * tmp + sin(alpha) * ptr->x.x[2];
     ptr->x.x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x.x[2];
     for (int j=0;j<MDSteps;j++) {
       tmp = Trajectories[ptr].R.at(j).x[0];
-      Trajectories[ptr].R.at(j).x[0] =  cos(alpha) * tmp + sin(alpha) * Trajectories[ptr].R.at(j).x[2]; 
+      Trajectories[ptr].R.at(j).x[0] =  cos(alpha) * tmp + sin(alpha) * Trajectories[ptr].R.at(j).x[2];
       Trajectories[ptr].R.at(j).x[2] = -sin(alpha) * tmp + cos(alpha) * Trajectories[ptr].R.at(j).x[2];
     }
-  }      
+  }
   // rotate n vector
   tmp = n->x[0];
@@ -1800 +1160 @@
   n->Output((ofstream *)&cout);
   cout << endl;
-  
+
   // rotate on z-y plane
   ptr = start;
@@ -1808 +1168 @@
     ptr = ptr->next;
     tmp = ptr->x.x[1];
-    ptr->x.x[1] =  cos(alpha) * tmp + sin(alpha) * ptr->x.x[2]; 
+    ptr->x.x[1] =  cos(alpha) * tmp + sin(alpha) * ptr->x.x[2];
     ptr->x.x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x.x[2];
     for (int j=0;j<MDSteps;j++) {
       tmp = Trajectories[ptr].R.at(j).x[1];
-      Trajectories[ptr].R.at(j).x[1] =  cos(alpha) * tmp + sin(alpha) * Trajectories[ptr].R.at(j).x[2]; 
+      Trajectories[ptr].R.at(j).x[1] =  cos(alpha) * tmp + sin(alpha) * Trajectories[ptr].R.at(j).x[2];
       Trajectories[ptr].R.at(j).x[2] = -sin(alpha) * tmp + cos(alpha) * Trajectories[ptr].R.at(j).x[2];
     }
-  }      
+  }
   // rotate n vector (for consistency check)
   tmp = n->x[1];
   n->x[1] =  cos(alpha) * tmp +  sin(alpha) * n->x[2];
   n->x[2] = -sin(alpha) * tmp +  cos(alpha) * n->x[2];
-  
+
   cout << Verbose(1) << "alignment vector after second rotation: ";
   n->Output((ofstream *)&cout);
@@ -1831 +1191 @@
  * \return true - succeeded, false - atom not found in list
  */
-bool molecule::RemoveAtom(atom *pointer)  
-{ 
+bool molecule::RemoveAtom(atom *pointer)
+{
   if (ElementsInMolecule[pointer->type->Z] != 0)  // this would indicate an error
     ElementsInMolecule[pointer->type->Z]--;  // decrease number of atom of this element
   else
-    cerr << "ERROR: Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl; 
+    cerr << "ERROR: Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl;
   if (ElementsInMolecule[pointer->type->Z] == 0)  // was last atom of this element?
     ElementCount--;
@@ -1846 +1206 @@
  * \return true - succeeded, false - atom not found in list
  */
-bool molecule::CleanupMolecule() 
-{ 
-  return (cleanup(start,end) && cleanup(first,last)); 
+bool molecule::CleanupMolecule()
+{
+  return (cleanup(start,end) && cleanup(first,last));
 };
 
@@ -1862 +1222 @@
   } else {
     cout << Verbose(0) << "Atom not found in list." << endl;
-    return NULL;  
+    return NULL;
   }
 };
@@ -1911 +1271 @@
   struct lsq_params *par = (struct lsq_params *)params;
   atom *ptr = par->mol->start;
-  
+
   // initialize vectors
   a.x[0] = gsl_vector_get(x,0);
@@ -1941 +1301 @@
 {
     int np = 6;
-    
+
    const gsl_multimin_fminimizer_type *T =
      gsl_multimin_fminimizer_nmsimplex;
@@ -1947 +1307 @@
    gsl_vector *ss;
    gsl_multimin_function minex_func;
- 
+
    size_t iter = 0, i;
    int status;
    double size;
- 
+
    /* Initial vertex size vector */
    ss = gsl_vector_alloc (np);
- 
+
    /* Set all step sizes to 1 */
    gsl_vector_set_all (ss, 1.0);
- 
+
    /* Starting point */
    par->x = gsl_vector_alloc (np);
    par->mol = this;
- 
+
    gsl_vector_set (par->x, 0, 0.0);  // offset
    gsl_vector_set (par->x, 1, 0.0);
@@ -1968 +1328 @@
    gsl_vector_set (par->x, 4, 0.0);
    gsl_vector_set (par->x, 5, 1.0);
- 
+
    /* Initialize method and iterate */
    minex_func.f = &LeastSquareDistance;
    minex_func.n = np;
    minex_func.params = (void *)par;
- 
+
    s = gsl_multimin_fminimizer_alloc (T, np);
    gsl_multimin_fminimizer_set (s, &minex_func, par->x, ss);
- 
+
    do
      {
        iter++;
        status = gsl_multimin_fminimizer_iterate(s);
- 
+
        if (status)
          break;
- 
+
        size = gsl_multimin_fminimizer_size (s);
        status = gsl_multimin_test_size (size, 1e-2);
- 
+
        if (status == GSL_SUCCESS)
          {
            printf ("converged to minimum at\n");
          }
- 
+
        printf ("%5d ", (int)iter);
        for (i = 0; i < (size_t)np; i++)
@@ -2001 +1361 @@
      }
    while (status == GSL_CONTINUE && iter < 100);
- 
+
   for (i=0;i<(size_t)np;i++)
     gsl_vector_set(par->x, i, gsl_vector_get(s->x, i));
@@ -2018 +1378 @@
   int ElementNo, AtomNo;
   CountElements();
-  
+
   if (out == NULL) {
     return false;
@@ -2053 +1413 @@
   int ElementNo, AtomNo;
   CountElements();
-  
+
   if (out == NULL) {
     return false;
@@ -2100 +1460 @@
   atom *Walker = start;
   while (Walker->next != end) {
-    Walker = Walker->next; 
+    Walker = Walker->next;
 #ifdef ADDHYDROGEN
     if (Walker->type->Z != 1) {   // regard only non-hydrogen
@@ -2131 +1491 @@
   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;
@@ -2158 +1518 @@
 {
   atom *walker = NULL;
-  int No = 0;
+  int AtomNo = 0, ElementNo;
   time_t now;
-  
+  element *runner = NULL;
+
   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++;
+    AtomNo++;
   }
   if (out != NULL) {
-    *out << No << "\n\tCreated by molecuilder on " << ctime(&now);
-    walker = start;
-    while (walker->next != end) { // go through every atom of this element
-      walker = walker->next;
-      walker->OutputXYZLine(out);
+    *out << AtomNo << "\n\tCreated by molecuilder on " << ctime(&now);
+    ElementNo = 0;
+    runner = elemente->start;
+    while (runner->next != elemente->end) { // go through every element
+                runner = runner->next;
+      if (ElementsInMolecule[runner->Z]) { // if this element got atoms
+        ElementNo++;
+        walker = start;
+        while (walker->next != end) { // go through every atom of this element
+          walker = walker->next;
+          if (walker->type == runner) { // if this atom fits to element
+            walker->OutputXYZLine(out);
+          }
+        }
+      }
     }
     return true;
@@ -2203 +1574 @@
               Walker->nr = i;   // update number in molecule (for easier referencing in FragmentMolecule lateron)
               if (Walker->type->Z != 1) // count non-hydrogen atoms whilst at it
-                NoNonHydrogen++; 
+                NoNonHydrogen++;
               Free((void **)&Walker->Name, "molecule::CountAtoms: *walker->Name");
               Walker->Name = (char *) Malloc(sizeof(char)*6, "molecule::CountAtoms: *walker->Name");
@@ -2211 +1582 @@
             }
     } else
-        *out << Verbose(3) << "AtomCount is still " << AtomCount << ", thus counting nothing." << endl; 
+        *out << Verbose(3) << "AtomCount is still " << AtomCount << ", thus counting nothing." << endl;
   }
 };
@@ -2223 +1594 @@
         ElementsInMolecule[i] = 0;
         ElementCount = 0;
-        
+
   atom *walker = start;
   while (walker->next != end) {
@@ -2259 +1630 @@
     Binder = Binder->next;
     if (Binder->Cyclic)
-      No++;    
+      No++;
   }
   delete(BackEdgeStack);
@@ -2317 +1688 @@
 
 /** Creates an adjacency list of the molecule.
+ * We obtain an outside file with the indices of atoms which are bondmembers.
+ */
+void molecule::CreateAdjacencyList2(ofstream *out, ifstream *input)
+{
+
+        // 1 We will parse bonds out of the dbond file created by tremolo.
+                        int atom1, atom2, temp;
+                        atom *Walker, *OtherWalker;
+
+                if (!input)
+                {
+                        cout << Verbose(1) << "Opening silica failed \n";
+                };
+
+                        *input >> ws >> atom1;
+                        *input >> ws >> atom2;
+                cout << Verbose(1) << "Scanning file\n";
+                while (!input->eof()) // Check whether we read everything already
+                {
+                                *input >> ws >> atom1;
+                                *input >> ws >> atom2;
+                        if(atom2<atom1) //Sort indices of atoms in order
+                        {
+                                temp=atom1;
+                                atom1=atom2;
+                                atom2=temp;
+                        };
+
+                        Walker=start;
+                        while(Walker-> nr != atom1) // Find atom corresponding to first index
+                        {
+                                Walker = Walker->next;
+                        };
+                        OtherWalker = Walker->next;
+                        while(OtherWalker->nr != atom2) // Find atom corresponding to second index
+                        {
+                                OtherWalker= OtherWalker->next;
+                        };
+                        AddBond(Walker, OtherWalker); //Add the bond between the two atoms with respective indices.
+                        
+                }
+
+                CreateListOfBondsPerAtom(out);
+
+};
+
+
+/** Creates an adjacency list of the molecule.
  * Generally, we use the CSD approach to bond recognition, that is the the distance
  * between two atoms A and B must be within [Rcov(A)+Rcov(B)-t,Rcov(A)+Rcov(B)+t] with
- * a threshold t = 0.4 Angstroem. 
+ * a threshold t = 0.4 Angstroem.
  * To make it O(N log N) the function uses the linked-cell technique as follows:
  * The procedure is step-wise:
@@ -2336 +1755 @@
 void molecule::CreateAdjacencyList(ofstream *out, double bonddistance, bool IsAngstroem)
 {
+
   atom *Walker = NULL, *OtherWalker = NULL, *Candidate = NULL;
   int No, NoBonds, CandidateBondNo;
@@ -2344 +1764 @@
   Vector x;
   int FalseBondDegree = 0;
- 
+
   BondDistance = bonddistance; // * ((IsAngstroem) ? 1. : 1./AtomicLengthToAngstroem);
   *out << Verbose(0) << "Begin of CreateAdjacencyList." << endl;
@@ -2351 +1771 @@
     cleanup(first,last);
   }
-        
+
   // count atoms in molecule = dimension of matrix (also give each unique name and continuous numbering)
   CountAtoms(out);
@@ -2370 +1790 @@
     for (int i=NumberCells;i--;)
       CellList[i] = NULL;
-  
+
     // 2b. put all atoms into its corresponding list
     Walker = start;
@@ -2391 +1811 @@
       if (CellList[index] == NULL)  // allocate molecule if not done
         CellList[index] = new molecule(elemente);
-      OtherWalker = CellList[index]->AddCopyAtom(Walker); // add a copy of walker to this atom, father will be walker for later reference 
-      //*out << Verbose(1) << "Copy Atom is " << *OtherWalker << "." << endl; 
+      OtherWalker = CellList[index]->AddCopyAtom(Walker); // add a copy of walker to this atom, father will be walker for later reference
+      //*out << Verbose(1) << "Copy Atom is " << *OtherWalker << "." << endl;
     }
     //for (int i=0;i<NumberCells;i++)
       //*out << Verbose(1) << "Cell number " << i << ": " << CellList[i] << "." << endl;
-      
+
+
     // 3a. go through every cell
     for (N[0]=divisor[0];N[0]--;)
@@ -2409 +1830 @@
               Walker = Walker->next;
               //*out << Verbose(0) << "Current Atom is " << *Walker << "." << endl;
-              // 3c. check for possible bond between each atom in this and every one in the 27 cells 
+              // 3c. check for possible bond between each atom in this and every one in the 27 cells
               for (n[0]=-1;n[0]<=1;n[0]++)
                 for (n[1]=-1;n[1]<=1;n[1]++)
@@ -2441 +1862 @@
           }
         }
+
+
+
     // 4. free the cell again
     for (int i=NumberCells;i--;)
@@ -2447 +1871 @@
       }
     Free((void **)&CellList, "molecule::CreateAdjacencyList - ** CellList");
-    
+
     // create the adjacency list per atom
     CreateListOfBondsPerAtom(out);
-                
+
     // correct Bond degree of each bond by checking both bond partners for a mismatch between valence and current sum of bond degrees,
     // iteratively increase the one first where the other bond partner has the fewest number of bonds (i.e. in general bonds oxygene
@@ -2499 +1923 @@
                         *out << Verbose(1) << "BondCount is " << BondCount << ", no bonds between any of the " << AtomCount << " atoms." << endl;
           *out << Verbose(1) << "I detected " << BondCount << " bonds in the molecule with distance " << bonddistance << ", " << FalseBondDegree << " bonds could not be corrected." << endl;
-                
+
           // output bonds for debugging (if bond chain list was correctly installed)
           *out << Verbose(1) << endl << "From contents of bond chain list:";
@@ -2509 +1933 @@
     *out << endl;
   } else
-        *out << Verbose(1) << "AtomCount is " << AtomCount << ", thus no bonds, no connections!." << endl; 
+        *out << Verbose(1) << "AtomCount is " << AtomCount << ", thus no bonds, no connections!." << endl;
   *out << Verbose(0) << "End of CreateAdjacencyList." << endl;
   Free((void **)&matrix, "molecule::CreateAdjacencyList: *matrix");
-};
+
+};
+
+
 
 /** Performs a Depth-First search on this molecule.
@@ -2533 +1960 @@
   bond *Binder = NULL;
   bool BackStepping = false;
-  
+
   *out << Verbose(0) << "Begin of DepthFirstSearchAnalysis" << endl;
-  
+
   ResetAllBondsToUnused();
   ResetAllAtomNumbers();
@@ -2548 +1975 @@
     LeafWalker->Leaf = new molecule(elemente);
     LeafWalker->Leaf->AddCopyAtom(Root);
-    
+
     OldGraphNr = CurrentGraphNr;
     Walker = Root;
@@ -2559 +1986 @@
           AtomStack->Push(Walker);
           CurrentGraphNr++;
-        }      
+        }
         do { // (3) if Walker has no unused egdes, go to (5)
           BackStepping = false; // reset backstepping flag for (8)
@@ -2593 +2020 @@
           Binder = NULL;
       } while (1);  // (2)
-      
+
       // if we came from backstepping, yet there were no more unused bonds, we end up here with no Ancestor, because Walker is Root! Then we are finished!
       if ((Walker == Root) && (Binder == NULL))
         break;
-        
-      // (5) if Ancestor of Walker is ... 
+
+      // (5) if Ancestor of Walker is ...
       *out << Verbose(1) << "(5) Number of Walker[" << Walker->Name << "]'s Ancestor[" << Walker->Ancestor->Name << "] is " << Walker->Ancestor->GraphNr << "." << endl;
       if (Walker->Ancestor->GraphNr != Root->GraphNr) {
@@ -2641 +2068 @@
         } while (OtherAtom != Walker);
         ComponentNumber++;
-    
+
         // (11) Root is separation vertex,  set Walker to Root and go to (4)
         Walker = Root;
@@ -2654 +2081 @@
 
     // From OldGraphNr to CurrentGraphNr ranges an disconnected subgraph
-    *out << Verbose(0) << "Disconnected subgraph ranges from " << OldGraphNr << " to " << CurrentGraphNr << "." << endl;    
+    *out << Verbose(0) << "Disconnected subgraph ranges from " << OldGraphNr << " to " << CurrentGraphNr << "." << endl;
     LeafWalker->Leaf->Output(out);
     *out << endl;
@@ -2662 +2089 @@
       //*out << Verbose(1) << "Current next subgraph root candidate is " << Root->Name << "." << endl;
       if (Root->GraphNr != -1) // if already discovered, step on
-        Root = Root->next; 
+        Root = Root->next;
     }
   }
@@ -2684 +2111 @@
     *out << " with Lowpoint " << Walker->LowpointNr << " and Graph Nr. " << Walker->GraphNr << "." << endl;
   }
-  
+
   *out << Verbose(1) << "Final graph info for each bond is:" << endl;
   Binder = first;
@@ -2695 +2122 @@
     *out << ((Binder->rightatom->SeparationVertex) ? "SP," : "") << "L" << Binder->rightatom->LowpointNr << " G" << Binder->rightatom->GraphNr << " Comp.";
     OutputComponentNumber(out, Binder->rightatom);
-    *out << ">." << endl; 
+    *out << ">." << endl;
     if (Binder->Cyclic) // cyclic ??
       *out << Verbose(3) << "Lowpoint at each side are equal: CYCLIC!" << endl;
@@ -2710 +2137 @@
  * the other our initial Walker - and do a Breadth First Search for the Root. We mark down each Predecessor and as soon as
  * we have found the Root via BFS, we may climb back the closed cycle via the Predecessors. Thereby we mark atoms and bonds
- * as cyclic and print out the cycles. 
+ * as cyclic and print out the cycles.
  * \param *out output stream for debugging
  * \param *BackEdgeStack stack with all back edges found during DFS scan. Beware: This stack contains the bonds from the total molecule, not from the subgraph!
@@ -2721 +2148 @@
   int *ShortestPathList = (int *) Malloc(sizeof(int)*AtomCount, "molecule::CyclicStructureAnalysis: *ShortestPathList");
   enum Shading *ColorList = (enum Shading *) Malloc(sizeof(enum Shading)*AtomCount, "molecule::CyclicStructureAnalysis: *ColorList");
-  class StackClass<atom *> *BFSStack = new StackClass<atom *> (AtomCount);   // will hold the current ring 
+  class StackClass<atom *> *BFSStack = new StackClass<atom *> (AtomCount);   // will hold the current ring
   class StackClass<atom *> *TouchedStack = new StackClass<atom *> (AtomCount);   // contains all "touched" atoms (that need to be reset after BFS loop)
   atom *Walker = NULL, *OtherAtom = NULL, *Root = NULL;
@@ -2733 +2160 @@
     ColorList[i] = white;
   }
-  
+
   *out << Verbose(1) << "Back edge list - ";
   BackEdgeStack->Output(out);
-  
+
   *out << Verbose(1) << "Analysing cycles ... " << endl;
   NumCycles = 0;
@@ -2742 +2169 @@
     BackEdge = BackEdgeStack->PopFirst();
     // this is the target
-    Root = BackEdge->leftatom;  
+    Root = BackEdge->leftatom;
     // this is the source point
-    Walker = BackEdge->rightatom; 
+    Walker = BackEdge->rightatom;
     ShortestPathList[Walker->nr] = 0;
     BFSStack->ClearStack();  // start with empty BFS stack
@@ -2758 +2185 @@
         if (Binder != BackEdge) { // only walk along DFS spanning tree (otherwise we always find SP of one being backedge Binder)
           OtherAtom = Binder->GetOtherAtom(Walker);
-#ifdef ADDHYDROGEN          
+#ifdef ADDHYDROGEN
           if (OtherAtom->type->Z != 1) {
 #endif
@@ -2767 +2194 @@
               PredecessorList[OtherAtom->nr] = Walker;  // Walker is the predecessor
               ShortestPathList[OtherAtom->nr] = ShortestPathList[Walker->nr]+1;
-              *out << Verbose(2) << "Coloring OtherAtom " << OtherAtom->Name << " lightgray, its predecessor is " << Walker->Name << " and its Shortest Path is " << ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl; 
+              *out << Verbose(2) << "Coloring OtherAtom " << OtherAtom->Name << " lightgray, its predecessor is " << Walker->Name << " and its Shortest Path is " << ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl;
               //if (ShortestPathList[OtherAtom->nr] < MinimumRingSize[Walker->GetTrueFather()->nr]) { // Check for maximum distance
                 *out << Verbose(3) << "Putting OtherAtom into queue." << endl;
@@ -2777 +2204 @@
             if (OtherAtom == Root)
               break;
-#ifdef ADDHYDROGEN          
+#ifdef ADDHYDROGEN
           } else {
             *out << Verbose(2) << "Skipping hydrogen atom " << *OtherAtom << "." << endl;
@@ -2815 +2242 @@
       }
     } while ((!BFSStack->IsEmpty()) && (OtherAtom != Root) && (OtherAtom != NULL)); // || (ShortestPathList[OtherAtom->nr] < MinimumRingSize[Walker->GetTrueFather()->nr])));
-    
+
     if (OtherAtom == Root) {
       // now climb back the predecessor list and thus find the cycle members
@@ -2843 +2270 @@
       *out << Verbose(1) << "No ring containing " << *Root << " with length equal to or smaller than " << MinimumRingSize[Walker->GetTrueFather()->nr] << " found." << endl;
     }
-    
+
     // now clean the lists
     while (!TouchedStack->IsEmpty()){
@@ -2853 +2280 @@
   }
   if (MinRingSize != -1) {
-    // go over all atoms 
+    // go over all atoms
     Root = start;
     while(Root->next != end) {
       Root = Root->next;
-      
+
       if (MinimumRingSize[Root->GetTrueFather()->nr] == AtomCount) { // check whether MinimumRingSize is set, if not BFS to next where it is
         Walker = Root;
@@ -2894 +2321 @@
           }
           ColorList[Walker->nr] = black;
-          //*out << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl; 
+          //*out << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl;
         }
-    
+
         // now clean the lists
         while (!TouchedStack->IsEmpty()){
@@ -2945 +2372 @@
 void molecule::OutputComponentNumber(ofstream *out, atom *vertex)
 {
-  for(int i=0;i<NumberOfBondsPerAtom[vertex->nr];i++) 
+  for(int i=0;i<NumberOfBondsPerAtom[vertex->nr];i++)
     *out << vertex->ComponentNr[i] << "  ";
 };
@@ -3023 +2450 @@
 {
   int c = 0;
-  int FragmentCount; 
+  int FragmentCount;
   // get maximum bond degree
   atom *Walker = start;
@@ -3033 +2460 @@
   *out << Verbose(1) << "Upper limit for this subgraph is " << FragmentCount << " for " << NoNonHydrogen << " non-H atoms with maximum bond degree of " << c << "." << endl;
   return FragmentCount;
-};  
+};
 
 /** Scans a single line for number and puts them into \a KeySet.
  * \param *out output stream for debugging
  * \param *buffer buffer to scan
- * \param &CurrentSet filled KeySet on return 
+ * \param &CurrentSet filled KeySet on return
  * \return true - at least one valid atom id parsed, false - CurrentSet is empty
  */
@@ -3046 +2473 @@
   int AtomNr;
   int status = 0;
-  
+
   line.str(buffer);
   while (!line.eof()) {
@@ -3082 +2509 @@
   double TEFactor;
   char *filename = (char *) Malloc(sizeof(char)*MAXSTRINGSIZE, "molecule::ParseKeySetFile - filename");
-  
+
   if (FragmentList == NULL) { // check list pointer
     FragmentList = new Graph;
   }
-  
+
   // 1st pass: open file and read
   *out << Verbose(1) << "Parsing the KeySet file ... " << endl;
@@ -3115 +2542 @@
     status = false;
   }
-  
+
   // 2nd pass: open TEFactors file and read
   *out << Verbose(1) << "Parsing the TEFactors file ... " << endl;
@@ -3127 +2554 @@
         InputFile >> TEFactor;
         (*runner).second.second = TEFactor;
-        *out << Verbose(2) << "Setting " << ++NumberOfFragments << " fragment's TEFactor to " << (*runner).second.second << "." << endl; 
+        *out << Verbose(2) << "Setting " << ++NumberOfFragments << " fragment's TEFactor to " << (*runner).second.second << "." << endl;
       } else {
         status = false;
@@ -3168 +2595 @@
   if(output != NULL) {
     for(Graph::iterator runner = KeySetList.begin(); runner != KeySetList.end(); runner++) {
-      for (KeySet::iterator sprinter = (*runner).first.begin();sprinter != (*runner).first.end(); sprinter++) { 
+      for (KeySet::iterator sprinter = (*runner).first.begin();sprinter != (*runner).first.end(); sprinter++) {
         if (sprinter != (*runner).first.begin())
           output << "\t";
@@ -3234 +2661 @@
     status = false;
   }
-  
+
   return status;
 };
@@ -3243 +2670 @@
  * \param **ListOfAtoms allocated (molecule::AtomCount) and filled lookup table for ids (Atom::nr) to *Atom
  * \return true - structure is equal, false - not equivalence
- */ 
+ */
 bool molecule::CheckAdjacencyFileAgainstMolecule(ofstream *out, char *path, atom **ListOfAtoms)
 {
@@ -3250 +2677 @@
   bool status = true;
   char *buffer = (char *) Malloc(sizeof(char)*MAXSTRINGSIZE, "molecule::CheckAdjacencyFileAgainstMolecule: *buffer");
-  
+
   filename << path << "/" << FRAGMENTPREFIX << ADJACENCYFILE;
   File.open(filename.str().c_str(), ios::out);
@@ -3309 +2736 @@
   *out << endl;
   Free((void **)&buffer, "molecule::CheckAdjacencyFileAgainstMolecule: *buffer");
-  
+
   return status;
 };
@@ -3331 +2758 @@
   for(int i=AtomCount;i--;)
     AtomMask[i] = false;
-  
+
   if (Order < 0) { // adaptive increase of BondOrder per site
     if (AtomMask[AtomCount] == true)  // break after one step
@@ -3371 +2798 @@
           line >> ws >> Value; // skip time entry
           line >> ws >> Value;
-          No -= 1;  // indices start at 1 in file, not 0 
+          No -= 1;  // indices start at 1 in file, not 0
           //*out << Verbose(2) << " - yields (" << No << "," << Value << ", " << FragOrder << ")" << endl;
 
@@ -3380 +2807 @@
             // as the smallest number in each set has always been the root (we use global id to keep the doubles away), seek smallest and insert into AtomMask
             pair <map<int, pair<double,int> >::iterator, bool> InsertedElement = AdaptiveCriteriaList.insert( make_pair(*((*marker).second.begin()), pair<double,int>( fabs(Value), FragOrder) ));
-            map<int, pair<double,int> >::iterator PresentItem = InsertedElement.first; 
+            map<int, pair<double,int> >::iterator PresentItem = InsertedElement.first;
             if (!InsertedElement.second) { // this root is already present
-              if ((*PresentItem).second.second < FragOrder)  // if order there is lower, update entry with higher-order term 
-                //if ((*PresentItem).second.first < (*runner).first)    // as higher-order terms are not always better, we skip this part (which would always include this site into adaptive increase) 
+              if ((*PresentItem).second.second < FragOrder)  // if order there is lower, update entry with higher-order term
+                //if ((*PresentItem).second.first < (*runner).first)    // as higher-order terms are not always better, we skip this part (which would always include this site into adaptive increase)
                 {  // if value is smaller, update value and order
                 (*PresentItem).second.first = fabs(Value);
@@ -3421 +2848 @@
         Walker = FindAtom(No);
         //if (Walker->AdaptiveOrder < MinimumRingSize[Walker->nr]) {
-          *out << Verbose(2) << "Root " << No << " is still above threshold (10^{" << Order <<"}: " << runner->first << ", setting entry " << No << " of Atom mask to true." << endl; 
+          *out << Verbose(2) << "Root " << No << " is still above threshold (10^{" << Order <<"}: " << runner->first << ", setting entry " << No << " of Atom mask to true." << endl;
           AtomMask[No] = true;
           status = true;
         //} else
-          //*out << Verbose(2) << "Root " << No << " is still above threshold (10^{" << Order <<"}: " << runner->first << ", however MinimumRingSize of " << MinimumRingSize[Walker->nr] << " does not allow further adaptive increase." << endl; 
+          //*out << Verbose(2) << "Root " << No << " is still above threshold (10^{" << Order <<"}: " << runner->first << ", however MinimumRingSize of " << MinimumRingSize[Walker->nr] << " does not allow further adaptive increase." << endl;
       }
       // close and done
@@ -3459 +2886 @@
     if ((Order == 0) && (AtomMask[AtomCount] == false))  // single stepping, just check
       status = true;
-      
+
     if (!status) {
       if (Order == 0)
@@ -3467 +2894 @@
     }
   }
-  
+
   // print atom mask for debugging
   *out << "              ";
@@ -3476 +2903 @@
     *out << (AtomMask[i] ? "t" : "f");
   *out << endl;
- 
+
   return status;
 };
@@ -3490 +2917 @@
   int AtomNo = 0;
   atom *Walker = NULL;
-  
+
   if (SortIndex != NULL) {
     *out << Verbose(1) << "SortIndex is " << SortIndex << " and not NULL as expected." << endl;
@@ -3548 +2975 @@
   atom **ListOfAtoms = NULL;
   atom ***ListOfLocalAtoms = NULL;
-  bool *AtomMask = NULL; 
-  
+  bool *AtomMask = NULL;
+
   *out << endl;
 #ifdef ADDHYDROGEN
@@ -3558 +2985 @@
 
   // ++++++++++++++++++++++++++++ INITIAL STUFF: Bond structure analysis, file parsing, ... ++++++++++++++++++++++++++++++++++++++++++
-  
+
   // ===== 1. Check whether bond structure is same as stored in files ====
-  
+
   // fill the adjacency list
   CreateListOfBondsPerAtom(out);
@@ -3566 +2993 @@
   // create lookup table for Atom::nr
   FragmentationToDo = FragmentationToDo && CreateFatherLookupTable(out, start, end, ListOfAtoms, AtomCount);
-  
+
   // === compare it with adjacency file ===
-  FragmentationToDo = FragmentationToDo && CheckAdjacencyFileAgainstMolecule(out, configuration->configpath, ListOfAtoms); 
+  FragmentationToDo = FragmentationToDo && CheckAdjacencyFileAgainstMolecule(out, configuration->configpath, ListOfAtoms);
   Free((void **)&ListOfAtoms, "molecule::FragmentMolecule - **ListOfAtoms");
 
@@ -3582 +3009 @@
   while (MolecularWalker->next != NULL) {
     MolecularWalker = MolecularWalker->next;
-    *out << Verbose(0) << "Analysing the cycles of subgraph " << MolecularWalker->Leaf << " with nr. " << FragmentCounter << "." << endl;
     LocalBackEdgeStack = new StackClass<bond *> (MolecularWalker->Leaf->BondCount);
 //    // check the list of local atoms for debugging
@@ -3591 +3017 @@
 //      else
 //        *out << "\t" << ListOfLocalAtoms[FragmentCounter][i]->Name;
+    *out << Verbose(0) << "Gathering local back edges for subgraph " << MolecularWalker->Leaf << " with nr. " << FragmentCounter << "." << endl;
     MolecularWalker->Leaf->PickLocalBackEdges(out, ListOfLocalAtoms[FragmentCounter++], BackEdgeStack, LocalBackEdgeStack);
+    *out << Verbose(0) << "Analysing the cycles of subgraph " << MolecularWalker->Leaf << " with nr. " << FragmentCounter << "." << endl;
     MolecularWalker->Leaf->CyclicStructureAnalysis(out, LocalBackEdgeStack, MinimumRingSize);
+    *out << Verbose(0) << "Done with Analysing the cycles of subgraph " << MolecularWalker->Leaf << " with nr. " << FragmentCounter << "." << endl;
     delete(LocalBackEdgeStack);
   }
-  
+
   // ===== 3. if structure still valid, parse key set file and others =====
   FragmentationToDo = FragmentationToDo && ParseKeySetFile(out, configuration->configpath, ParsedFragmentList);
@@ -3601 +3030 @@
   // ===== 4. check globally whether there's something to do actually (first adaptivity check)
   FragmentationToDo = FragmentationToDo && ParseOrderAtSiteFromFile(out, configuration->configpath);
-  
-  // =================================== Begin of FRAGMENTATION =============================== 
-  // ===== 6a. assign each keyset to its respective subgraph ===== 
+
+  // =================================== Begin of FRAGMENTATION ===============================
+  // ===== 6a. assign each keyset to its respective subgraph =====
   Subgraphs->next->AssignKeySetsToFragment(out, this, ParsedFragmentList, ListOfLocalAtoms, FragmentList, (FragmentCounter = 0), true);
 
@@ -3614 +3043 @@
     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) ===== 
+    // ===== 6b. fill RootStack for each subgraph (second adaptivity check) =====
     Subgraphs->next->FillRootStackForSubgraphs(out, RootStack, AtomMask, (FragmentCounter = 0));
 
@@ -3623 +3052 @@
       MolecularWalker = MolecularWalker->next;
       *out << Verbose(1) << "Fragmenting subgraph " << MolecularWalker << "." << endl;
-      // output ListOfBondsPerAtom for debugging
-      MolecularWalker->Leaf->OutputListOfBonds(out);
+      //MolecularWalker->Leaf->OutputListOfBonds(out);  // output ListOfBondsPerAtom for debugging
       if (MolecularWalker->Leaf->first->next != MolecularWalker->Leaf->last) {
-      
         // call BOSSANOVA method
         *out << Verbose(0) << endl << " ========== BOND ENERGY of subgraph " << FragmentCounter << " ========================= " << endl;
@@ -3640 +3067 @@
   delete(ParsedFragmentList);
   delete[](MinimumRingSize);
-  
+
 
   // ==================================== End of FRAGMENTATION ============================================
@@ -3646 +3073 @@
   // ===== 8a. translate list into global numbers (i.e. ones that are valid in "this" molecule, not in MolecularWalker->Leaf)
   Subgraphs->next->TranslateIndicesToGlobalIDs(out, FragmentList, (FragmentCounter = 0), TotalNumberOfKeySets, TotalGraph);
-  
+
   // free subgraph memory again
   FragmentCounter = 0;
@@ -3671 +3098 @@
     }
     *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);
-      
+
       *out << Verbose(1) << "Writing " << BondFragments->NumberOfMolecules << " possible bond fragmentation configs" << endl;
-      if (BondFragments->OutputConfigForListOfFragments(out, FRAGMENTPREFIX, configuration, SortIndex, true, true))
+      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 
+
+      // store keysets file
       StoreKeySetFile(out, TotalGraph, configuration->configpath);
-  
-      // store Adjacency file 
+
+      // 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"); 
+      Free((void **)&SortIndex, "molecule::FragmentMolecule: *SortIndex");
     } else
       *out << Verbose(1) << "FragmentList is zero on return, splitting failed." << endl;
-  //} else 
+  //} else
   //  *out << Verbose(1) << "No fragments to store." << endl;
   *out << Verbose(0) << "End of bond fragmentation." << endl;
@@ -3733 +3160 @@
   atom *Walker = NULL, *OtherAtom = NULL;
   ReferenceStack->Push(Binder);
-  
+
   do {  // go through all bonds and push local ones
     Walker = ListOfLocalAtoms[Binder->leftatom->nr];  // get one atom in the reference molecule
-    if (Walker == NULL) // if this Walker exists in the subgraph ... 
-      continue;
-    for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {    // go through the local list of bonds
-      OtherAtom = ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
-      if (OtherAtom == ListOfLocalAtoms[Binder->rightatom->nr]) { // found the bond
-        LocalStack->Push(ListOfBondsPerAtom[Walker->nr][i]);
-        break;
-      }
-    }
+    if (Walker != NULL) // if this Walker exists in the subgraph ...
+        for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {    // go through the local list of bonds
+        OtherAtom = ListOfBondsPerAtom[Walker->nr][i]->GetOtherAtom(Walker);
+              if (OtherAtom == ListOfLocalAtoms[Binder->rightatom->nr]) { // found the bond
+              LocalStack->Push(ListOfBondsPerAtom[Walker->nr][i]);
+                                        *out << Verbose(3) << "Found local edge " << *(ListOfBondsPerAtom[Walker->nr][i]) << "." << endl;
+                break;
+            }
+        }
     Binder = ReferenceStack->PopFirst();  // loop the stack for next item
+                *out << Verbose(3) << "Current candidate edge " << Binder << "." << endl;
     ReferenceStack->Push(Binder);
   } while (FirstBond != Binder);
-  
+
   return status;
 };
@@ -3824 +3252 @@
       Walker->AdaptiveOrder = OrderArray[Walker->nr];
       Walker->MaxOrder = MaxArray[Walker->nr];
-      *out << Verbose(2) << *Walker << " gets order " << (int)Walker->AdaptiveOrder << " and is " << (!Walker->MaxOrder ? "not " : " ") << "maxed." << endl; 
+      *out << Verbose(2) << *Walker << " gets order " << (int)Walker->AdaptiveOrder << " and is " << (!Walker->MaxOrder ? "not " : " ") << "maxed." << endl;
     }
     file.close();
@@ -3835 +3263 @@
   Free((void **)&OrderArray, "molecule::ParseOrderAtSiteFromFile - *OrderArray");
   Free((void **)&MaxArray, "molecule::ParseOrderAtSiteFromFile - *MaxArray");
-  
+
   *out << Verbose(1) << "End of ParseOrderAtSiteFromFile" << endl;
   return status;
@@ -3893 +3321 @@
   Walker = start;
   while (Walker->next != end) {
-    Walker = Walker->next; 
+    Walker = Walker->next;
     *out << Verbose(4) << "Atom " << Walker->Name << "/" << Walker->nr << " with " << NumberOfBondsPerAtom[Walker->nr] << " bonds: ";
     TotalDegree = 0;
@@ -3901 +3329 @@
     }
     *out << " -- TotalDegree: " << TotalDegree << endl;
-  }      
+  }
   *out << Verbose(1) << "End of Creating ListOfBondsPerAtom." << endl << endl;
 };
@@ -3907 +3335 @@
 /** Adds atoms up to \a BondCount distance from \a *Root and notes them down in \a **AddedAtomList.
  * Gray vertices are always enqueued in an StackClass<atom *> FIFO queue, the rest is usual BFS with adding vertices found was
- * white and putting into queue. 
+ * white and putting into queue.
  * \param *out output stream for debugging
  * \param *Mol Molecule class to add atoms to
@@ -3916 +3344 @@
  * \param BondOrder maximum distance for vertices to add
  * \param IsAngstroem lengths are in angstroem or bohrradii
- */  
+ */
 void molecule::BreadthFirstSearchAdd(ofstream *out, molecule *Mol, atom **&AddedAtomList, bond **&AddedBondList, atom *Root, bond *Bond, int BondOrder, bool IsAngstroem)
 {
@@ -3942 +3370 @@
   }
   ShortestPathList[Root->nr] = 0;
-  
+
   // and go on ... Queue always contains all lightgray vertices
   while (!AtomStack->IsEmpty()) {
@@ -3950 +3378 @@
     // followed by n+1 till top of stack.
     Walker = AtomStack->PopFirst(); // pop oldest added
-    *out << Verbose(1) << "Current Walker is: " << Walker->Name << ", and has " << NumberOfBondsPerAtom[Walker->nr] << " bonds." << endl; 
+    *out << Verbose(1) << "Current Walker is: " << Walker->Name << ", and has " << NumberOfBondsPerAtom[Walker->nr] << " bonds." << endl;
     for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
       Binder = ListOfBondsPerAtom[Walker->nr][i];
@@ -3957 +3385 @@
         *out << Verbose(2) << "Current OtherAtom is: " << OtherAtom->Name << " for bond " << *Binder << "." << endl;
         if (ColorList[OtherAtom->nr] == white) {
-          if (Binder != Bond) // let other atom white if it's via Root bond. In case it's cyclic it has to be reached again (yet Root is from OtherAtom already black, thus no problem) 
+          if (Binder != Bond) // let other atom white if it's via Root bond. In case it's cyclic it has to be reached again (yet Root is from OtherAtom already black, thus no problem)
             ColorList[OtherAtom->nr] = lightgray;
           PredecessorList[OtherAtom->nr] = Walker;  // Walker is the predecessor
           ShortestPathList[OtherAtom->nr] = ShortestPathList[Walker->nr]+1;
-          *out << Verbose(2) << "Coloring OtherAtom " << OtherAtom->Name << " " << ((ColorList[OtherAtom->nr] == white) ? "white" : "lightgray") << ", its predecessor is " << Walker->Name << " and its Shortest Path is " << ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl; 
+          *out << Verbose(2) << "Coloring OtherAtom " << OtherAtom->Name << " " << ((ColorList[OtherAtom->nr] == white) ? "white" : "lightgray") << ", its predecessor is " << Walker->Name << " and its Shortest Path is " << ShortestPathList[OtherAtom->nr] << " egde(s) long." << endl;
           if ((((ShortestPathList[OtherAtom->nr] < BondOrder) && (Binder != Bond))) ) { // Check for maximum distance
             *out << Verbose(3);
@@ -3999 +3427 @@
               } else {
 #ifdef ADDHYDROGEN
-                Mol->AddHydrogenReplacementAtom(out, Binder, AddedAtomList[Walker->nr], Walker, OtherAtom, ListOfBondsPerAtom[Walker->nr], NumberOfBondsPerAtom[Walker->nr], IsAngstroem);
+                if (!Mol->AddHydrogenReplacementAtom(out, Binder, AddedAtomList[Walker->nr], Walker, OtherAtom, ListOfBondsPerAtom[Walker->nr], NumberOfBondsPerAtom[Walker->nr], IsAngstroem))
+                  exit(1);
 #endif
               }
@@ -4008 +3437 @@
           // This has to be a cyclic bond, check whether it's present ...
           if (AddedBondList[Binder->nr] == NULL) {
-            if ((Binder != Bond) && (Binder->Cyclic) && (((ShortestPathList[Walker->nr]+1) < BondOrder))) { 
+            if ((Binder != Bond) && (Binder->Cyclic) && (((ShortestPathList[Walker->nr]+1) < BondOrder))) {
               AddedBondList[Binder->nr] = Mol->AddBond(AddedAtomList[Walker->nr], AddedAtomList[OtherAtom->nr], Binder->BondDegree);
               AddedBondList[Binder->nr]->Cyclic = Binder->Cyclic;
@@ -4014 +3443 @@
             } else { // if it's root bond it has to broken (otherwise we would not create the fragments)
 #ifdef ADDHYDROGEN
-              Mol->AddHydrogenReplacementAtom(out, Binder, AddedAtomList[Walker->nr], Walker, OtherAtom, ListOfBondsPerAtom[Walker->nr], NumberOfBondsPerAtom[Walker->nr], IsAngstroem);
+              if(!Mol->AddHydrogenReplacementAtom(out, Binder, AddedAtomList[Walker->nr], Walker, OtherAtom, ListOfBondsPerAtom[Walker->nr], NumberOfBondsPerAtom[Walker->nr], IsAngstroem))
+                exit(1);
 #endif
             }
@@ -4022 +3452 @@
     }
     ColorList[Walker->nr] = black;
-    *out << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl; 
+    *out << Verbose(1) << "Coloring Walker " << Walker->Name << " black." << endl;
   }
   Free((void **)&PredecessorList, "molecule::BreadthFirstSearchAdd: **PredecessorList");
@@ -4046 +3476 @@
 
   *out << Verbose(2) << "Begin of BuildInducedSubgraph." << endl;
-  
+
   // reset parent list
   *out << Verbose(3) << "Resetting ParentList." << endl;
   for (int i=Father->AtomCount;i--;)
     ParentList[i] = NULL;
-  
+
   // fill parent list with sons
   *out << Verbose(3) << "Filling Parent List." << endl;
@@ -4092 +3522 @@
  * \param *&Leaf KeySet to look through
  * \param *&ShortestPathList list of the shortest path to decide which atom to suggest as removal candidate in the end
- * \param index of the atom suggested for removal 
+ * \param index of the atom suggested for removal
  */
 int molecule::LookForRemovalCandidate(ofstream *&out, KeySet *&Leaf, int *&ShortestPathList)
@@ -4098 +3528 @@
   atom *Runner = NULL;
   int SP, Removal;
-  
+
   *out << Verbose(2) << "Looking for removal candidate." << endl;
   SP = -1; //0;  // not -1, so that Root is never removed
@@ -4116 +3546 @@
 /** Stores a fragment from \a KeySet into \a molecule.
  * First creates the minimal set of atoms from the KeySet, then creates the bond structure from the complete
- * molecule and adds missing hydrogen where bonds were cut. 
+ * molecule and adds missing hydrogen where bonds were cut.
  * \param *out output stream for debugging messages
- * \param &Leaflet pointer to KeySet structure 
+ * \param &Leaflet pointer to KeySet structure
  * \param IsAngstroem whether we have Ansgtroem or bohrradius
  * \return pointer to constructed molecule
@@ -4129 +3559 @@
   bool LonelyFlag = false;
   int size;
-  
+
 //  *out << Verbose(1) << "Begin of StoreFragmentFromKeyset." << endl;
-  
+
   Leaf->BondDistance = BondDistance;
   for(int i=NDIM*2;i--;)
-    Leaf->cell_size[i] = cell_size[i]; 
+    Leaf->cell_size[i] = cell_size[i];
 
   // initialise SonList (indicates when we need to replace a bond with hydrogen instead)
@@ -4147 +3577 @@
     size++;
   }
-  
+
   // create the bonds between all: Make it an induced subgraph and add hydrogen
 //  *out << Verbose(2) << "Creating bonds from father graph (i.e. induced subgraph creation)." << endl;
@@ -4157 +3587 @@
     if (SonList[FatherOfRunner->nr] != NULL)  {  // check if this, our father, is present in list
       // create all bonds
-      for (int i=0;i<NumberOfBondsPerAtom[FatherOfRunner->nr];i++) { // go through every bond of father 
+      for (int i=0;i<NumberOfBondsPerAtom[FatherOfRunner->nr];i++) { // go through every bond of father
         OtherFather = ListOfBondsPerAtom[FatherOfRunner->nr][i]->GetOtherAtom(FatherOfRunner);
 //        *out << Verbose(2) << "Father " << *FatherOfRunner << " of son " << *SonList[FatherOfRunner->nr] << " is bound to " << *OtherFather;
         if (SonList[OtherFather->nr] != NULL) {
-//          *out << ", whose son is " << *SonList[OtherFather->nr] << "." << endl; 
+//          *out << ", whose son is " << *SonList[OtherFather->nr] << "." << endl;
           if (OtherFather->nr > FatherOfRunner->nr) { // add bond (nr check is for adding only one of both variants: ab, ba)
 //            *out << Verbose(3) << "Adding Bond: ";
-//            *out << 
+//            *out <<
             Leaf->AddBond(Runner, SonList[OtherFather->nr], ListOfBondsPerAtom[FatherOfRunner->nr][i]->BondDegree);
 //            *out << "." << endl;
             //NumBonds[Runner->nr]++;
-          } else { 
+          } else {
 //            *out << Verbose(3) << "Not adding bond, labels in wrong order." << endl;
           }
           LonelyFlag = false;
         } else {
-//          *out << ", who has no son in this fragment molecule." << endl; 
+//          *out << ", who has no son in this fragment molecule." << endl;
 #ifdef ADDHYDROGEN
           //*out << Verbose(3) << "Adding Hydrogen to " << Runner->Name << " and a bond in between." << endl;
-          Leaf->AddHydrogenReplacementAtom(out, ListOfBondsPerAtom[FatherOfRunner->nr][i], Runner, FatherOfRunner, OtherFather, ListOfBondsPerAtom[FatherOfRunner->nr],NumberOfBondsPerAtom[FatherOfRunner->nr], IsAngstroem);
+          if(!Leaf->AddHydrogenReplacementAtom(out, ListOfBondsPerAtom[FatherOfRunner->nr][i], Runner, FatherOfRunner, OtherFather, ListOfBondsPerAtom[FatherOfRunner->nr],NumberOfBondsPerAtom[FatherOfRunner->nr], IsAngstroem))
+            exit(1);
 #endif
           //NumBonds[Runner->nr] += ListOfBondsPerAtom[FatherOfRunner->nr][i]->BondDegree;
@@ -4190 +3621 @@
     while ((Runner->next != Leaf->end) && (Runner->next->type->Z == 1)) // skip added hydrogen
       Runner = Runner->next;
-#endif       
+#endif
   }
   Leaf->CreateListOfBondsPerAtom(out);
@@ -4223 +3654 @@
   StackClass<atom *> *TouchedStack = new StackClass<atom *>((int)pow(4,Order)+2); // number of atoms reached from one with maximal 4 bonds plus Root itself
   StackClass<atom *> *SnakeStack = new StackClass<atom *>(Order+1); // equal to Order is not possible, as then the StackClass<atom *> cannot discern between full and empty stack!
-  MoleculeLeafClass *Leaflet = NULL, *TempLeaf = NULL; 
+  MoleculeLeafClass *Leaflet = NULL, *TempLeaf = NULL;
   MoleculeListClass *FragmentList = NULL;
   atom *Walker = NULL, *OtherAtom = NULL, *Root = NULL, *Removal = NULL;
@@ -4273 +3704 @@
                 // clear snake stack
           SnakeStack->ClearStack();
-    //SnakeStack->TestImplementation(out, start->next);    
+    //SnakeStack->TestImplementation(out, start->next);
 
     ///// INNER LOOP ////////////
@@ -4294 +3725 @@
         }
         if (ShortestPathList[Walker->nr] == -1) {
-          ShortestPathList[Walker->nr] = ShortestPathList[PredecessorList[Walker->nr]->nr]+1; 
+          ShortestPathList[Walker->nr] = ShortestPathList[PredecessorList[Walker->nr]->nr]+1;
           TouchedStack->Push(Walker); // mark every atom for lists cleanup later, whose shortest path has been changed
         }
@@ -4332 +3763 @@
                                 OtherAtom = Binder->GetOtherAtom(Walker);
             if ((Labels[OtherAtom->nr] != -1) && (Labels[OtherAtom->nr] < Labels[Root->nr])) { // we don't step up to labels bigger than us
-              *out << "Label " << Labels[OtherAtom->nr] << " is smaller than Root's " << Labels[Root->nr] << "." << endl; 
+              *out << "Label " << Labels[OtherAtom->nr] << " is smaller than Root's " << Labels[Root->nr] << "." << endl;
               //ColorVertexList[OtherAtom->nr] = lightgray;    // mark as explored
             } else { // otherwise check its colour and element
                                 if (
 #ifdef ADDHYDROGEN
-              (OtherAtom->type->Z != 1) && 
+              (OtherAtom->type->Z != 1) &&
 #endif
                     (ColorEdgeList[Binder->nr] == white)) {  // skip hydrogen, look for unexplored vertices
@@ -4347 +3778 @@
                 //} else {
                 //  *out << Verbose(3) << "Predecessor of " << OtherAtom->Name << " is " << PredecessorList[OtherAtom->nr]->Name << "." << endl;
-                //} 
+                //}
                 Walker = OtherAtom;
                 break;
               } else {
-                if (OtherAtom->type->Z == 1) 
+                if (OtherAtom->type->Z == 1)
                   *out << "Links to a hydrogen atom." << endl;
-                else                 
+                else
                   *out << "Bond has not white but " << GetColor(ColorEdgeList[Binder->nr]) << " color." << endl;
               }
@@ -4363 +3794 @@
           *out << Verbose(3) << "We have gone too far, stepping back to " << Walker->Name << "." << endl;
         }
-                        if (Walker != OtherAtom) {      // if no white neighbours anymore, color it black 
+                        if (Walker != OtherAtom) {      // if no white neighbours anymore, color it black
           *out << Verbose(2) << "Coloring " << Walker->Name << " black." << endl;
                                 ColorVertexList[Walker->nr] = black;
@@ -4370 +3801 @@
       }
         } while ((Walker != Root) || (ColorVertexList[Root->nr] != black));
-    *out << Verbose(2) << "Inner Looping is finished." << endl;    
+    *out << Verbose(2) << "Inner Looping is finished." << endl;
 
     // if we reset all AtomCount atoms, we have again technically O(N^2) ...
@@ -4386 +3817 @@
     }
   }
-  *out << Verbose(1) << "Outer Looping over all vertices is done." << endl; 
+  *out << Verbose(1) << "Outer Looping over all vertices is done." << endl;
 
   // copy together
-  *out << Verbose(1) << "Copying all fragments into MoleculeList structure." << endl; 
+  *out << Verbose(1) << "Copying all fragments into MoleculeList structure." << endl;
   FragmentList = new MoleculeListClass(FragmentCounter, AtomCount);
   RunningIndex = 0;
@@ -4460 +3891 @@
 
   NumCombinations = 1 << SetDimension;
-  
+
   // Hier muessen von 1 bis NumberOfBondsPerAtom[Walker->nr] alle Kombinationen
   // von Endstuecken (aus den Bonds) hinzugefￜￜgt werden und fￜￜr verbleibende ANOVAOrder
   // rekursiv GraphCrawler in der nￜￜchsten Ebene aufgerufen werden
-  
+
   *out << Verbose(1+verbosity) << "Begin of SPFragmentGenerator." << endl;
   *out << Verbose(1+verbosity) << "We are " << RootDistance << " away from Root, which is " << *FragmentSearch->Root << ", SubOrder is " << SubOrder << ", SetDimension is " << SetDimension << " and this means " <<  NumCombinations-1 << " combination(s)." << endl;
 
-  // initialised touched list (stores added atoms on this level) 
+  // initialised touched list (stores added atoms on this level)
   *out << Verbose(1+verbosity) << "Clearing touched list." << endl;
   for (TouchedIndex=SubOrder+1;TouchedIndex--;)  // empty touched list
     TouchedList[TouchedIndex] = -1;
   TouchedIndex = 0;
-  
+
   // create every possible combination of the endpieces
   *out << Verbose(1+verbosity) << "Going through all combinations of the power set." << endl;
@@ -4481 +3912 @@
     for (int j=SetDimension;j--;)
       bits += (i & (1 << j)) >> j;
-      
+
     *out << Verbose(1+verbosity) << "Current set is " << Binary(i | (1 << SetDimension)) << ", number of bits is " << bits << "." << endl;
     if (bits <= SubOrder) { // if not greater than additional atoms allowed on stack, continue
@@ -4489 +3920 @@
         bit = ((i & (1 << j)) != 0);  // mask the bit for the j-th bond
         if (bit) {  // if bit is set, we add this bond partner
-                OtherWalker = BondsSet[j]->rightatom;   // rightatom is always the one more distant, i.e. the one to add 
+                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;
           *out << Verbose(2+verbosity) << "Adding " << *OtherWalker << " with nr " << OtherWalker->nr << "." << endl;
-          TestKeySetInsert = FragmentSearch->FragmentSet->insert(OtherWalker->nr); 
+          TestKeySetInsert = FragmentSearch->FragmentSet->insert(OtherWalker->nr);
           if (TestKeySetInsert.second) {
             TouchedList[TouchedIndex++] = OtherWalker->nr;  // note as added
@@ -4505 +3936 @@
         }
       }
-      
-      SpaceLeft = SubOrder - Added ;// SubOrder - bits; // due to item's maybe being already present, this does not work anymore 
+
+      SpaceLeft = SubOrder - Added ;// SubOrder - bits; // due to item's maybe being already present, this does not work anymore
       if (SpaceLeft > 0) {
         *out << Verbose(1+verbosity) << "There's still some space left on stack: " << SpaceLeft << "." << endl;
@@ -4534 +3965 @@
           }
           *out << Verbose(2+verbosity) << "Calling subset generator " << SP << " away from root " << *FragmentSearch->Root << " with sub set dimension " << SubSetDimension << "." << endl;
-          SPFragmentGenerator(out, FragmentSearch, SP, BondsList, SubSetDimension, SubOrder-bits);  
+          SPFragmentGenerator(out, FragmentSearch, SP, BondsList, SubSetDimension, SubOrder-bits);
           Free((void **)&BondsList, "molecule::SPFragmentGenerator: **BondsList");
         }
@@ -4543 +3974 @@
         *out << Verbose(2) << "Found a new fragment[" << FragmentSearch->FragmentCounter << "], local nr.s are: ";
         for(KeySet::iterator runner = FragmentSearch->FragmentSet->begin(); runner != FragmentSearch->FragmentSet->end(); runner++)
-          *out << (*runner) << " "; 
+          *out << (*runner) << " ";
         *out << endl;
         //if (!CheckForConnectedSubgraph(out, FragmentSearch->FragmentSet))
@@ -4551 +3982 @@
         //Removal = StoreFragmentFromStack(out, FragmentSearch->Root, FragmentSearch->Leaflet, FragmentSearch->FragmentStack, FragmentSearch->ShortestPathList, &FragmentSearch->FragmentCounter, FragmentSearch->configuration);
       }
-      
+
       // --3-- remove all added items in this level from snake stack
       *out << Verbose(1+verbosity) << "Removing all items that were added on this SP level " << RootDistance << "." << endl;
@@ -4562 +3993 @@
       *out << Verbose(2) << "Remaining local nr.s on snake stack are: ";
       for(KeySet::iterator runner = FragmentSearch->FragmentSet->begin(); runner != FragmentSearch->FragmentSet->end(); runner++)
-        *out << (*runner) << " "; 
+        *out << (*runner) << " ";
       *out << endl;
       TouchedIndex = 0; // set Index to 0 for list of atoms added on this level
@@ -4569 +4000 @@
     }
   }
-  Free((void **)&TouchedList, "molecule::SPFragmentGenerator: *TouchedList");  
+  Free((void **)&TouchedList, "molecule::SPFragmentGenerator: *TouchedList");
   *out << Verbose(1+verbosity) << "End of SPFragmentGenerator, " << RootDistance << " away from Root " << *FragmentSearch->Root << " and SubOrder is " << SubOrder << "." << endl;
 };
@@ -4578 +4009 @@
  * \return true - connected, false - disconnected
  * \note this is O(n^2) for it's just a bug checker not meant for permanent use!
- */ 
+ */
 bool molecule::CheckForConnectedSubgraph(ofstream *out, KeySet *Fragment)
 {
@@ -4584 +4015 @@
   bool BondStatus = false;
   int size;
-  
+
   *out << Verbose(1) << "Begin of CheckForConnectedSubgraph" << endl;
   *out << Verbose(2) << "Disconnected atom: ";
-  
+
   // count number of atoms in graph
   size = 0;
@@ -4633 +4064 @@
  * \param *out output stream for debugging
  * \param Order bond order (limits BFS exploration and "number of digits" in power set generation
- * \param FragmentSearch UniqueFragments structure containing TEFactor, root atom and so on 
+ * \param FragmentSearch UniqueFragments structure containing TEFactor, root atom and so on
  * \param RestrictedKeySet Restricted vertex set to use in context of molecule
  * \return number of inserted fragments
  * \note ShortestPathList in FragmentSearch structure is probably due to NumberOfAtomsSPLevel and SP not needed anymore
  */
-int molecule::PowerSetGenerator(ofstream *out, int Order, struct UniqueFragments &FragmentSearch, KeySet RestrictedKeySet) 
+int molecule::PowerSetGenerator(ofstream *out, int Order, struct UniqueFragments &FragmentSearch, KeySet RestrictedKeySet)
 {
   int SP, AtomKeyNr;
@@ -4659 +4090 @@
     FragmentSearch.BondsPerSPCount[i] = 0;
   FragmentSearch.BondsPerSPCount[0] = 1;
-  Binder = new bond(FragmentSearch.Root, FragmentSearch.Root); 
+  Binder = new bond(FragmentSearch.Root, FragmentSearch.Root);
   add(Binder, FragmentSearch.BondsPerSPList[1]);
-  
+
   // do a BFS search to fill the SP lists and label the found vertices
   // Actually, we should construct a spanning tree vom the root atom and select all edges therefrom and put them into
   // according shortest path lists. However, we don't. Rather we fill these lists right away, as they do form a spanning
   // tree already sorted into various SP levels. That's why we just do loops over the depth (CurrentSP) and breadth
-  // (EdgeinSPLevel) of this tree ... 
+  // (EdgeinSPLevel) of this tree ...
   // In another picture, the bonds always contain a direction by rightatom being the one more distant from root and hence
   // naturally leftatom forming its predecessor, preventing the BFS"seeker" from continuing in the wrong direction.
@@ -4719 +4150 @@
     }
   }
-  
+
   // outputting all list for debugging
   *out << Verbose(0) << "Printing all found lists." << endl;
@@ -4728 +4159 @@
       Binder = Binder->next;
       *out << Verbose(2) << *Binder << endl;
-    } 
-  }
-  
+    }
+  }
+
   // creating fragments with the found edge sets  (may be done in reverse order, faster)
   SP = -1;  // the Root <-> Root edge must be subtracted!
@@ -4737 +4168 @@
     while (Binder->next != FragmentSearch.BondsPerSPList[2*i+1]) {
       Binder = Binder->next;
-      SP ++; 
+      SP ++;
     }
   }
@@ -4744 +4175 @@
     // start with root (push on fragment stack)
     *out << Verbose(0) << "Starting fragment generation with " << *FragmentSearch.Root << ", local nr is " << FragmentSearch.Root->nr << "." << endl;
-    FragmentSearch.FragmentSet->clear();  
+    FragmentSearch.FragmentSet->clear();
     *out << Verbose(0) << "Preparing subset for this root and calling generator." << endl;
     // prepare the subset and call the generator
     BondsList = (bond **) Malloc(sizeof(bond *)*FragmentSearch.BondsPerSPCount[0], "molecule::PowerSetGenerator: **BondsList");
     BondsList[0] = FragmentSearch.BondsPerSPList[0]->next;  // on SP level 0 there's only the root bond
-    
+
     SPFragmentGenerator(out, &FragmentSearch, 0, BondsList, FragmentSearch.BondsPerSPCount[0], Order);
-    
+
     Free((void **)&BondsList, "molecule::PowerSetGenerator: **BondsList");
   } else {
@@ -4760 +4191 @@
   // remove root from stack
   *out << Verbose(0) << "Removing root again from stack." << endl;
-  FragmentSearch.FragmentSet->erase(FragmentSearch.Root->nr);    
+  FragmentSearch.FragmentSet->erase(FragmentSearch.Root->nr);
 
   // free'ing the bonds lists
@@ -4779 +4210 @@
   }
 
-  // return list  
+  // return list
   *out << Verbose(0) << "End of PowerSetGenerator." << endl;
   return (FragmentSearch.FragmentCounter - Counter);
@@ -4810 +4241 @@
     // remove bonds that are beyond bonddistance
     for(int i=NDIM;i--;)
-      Translationvector.x[i] = 0.; 
+      Translationvector.x[i] = 0.;
     // scan all bonds
     Binder = first;
@@ -4857 +4288 @@
         }
       }
-      // re-add bond    
+      // re-add bond
       link(Binder, OtherBinder);
     } else {
@@ -4911 +4342 @@
         IteratorB++;
       } // end of while loop
-    }// end of check in case of equal sizes 
+    }// end of check in case of equal sizes
   }
   return false; // if we reach this point, they are equal
@@ -4955 +4386 @@
  * \param graph1 first (dest) graph
  * \param graph2 second (source) graph
- * \param *counter keyset counter that gets increased 
+ * \param *counter keyset counter that gets increased
  */
 inline void InsertGraphIntoGraph(ofstream *out, Graph &graph1, Graph &graph2, int *counter)
@@ -5000 +4431 @@
   int RootKeyNr, RootNr;
   struct UniqueFragments FragmentSearch;
-  
+
   *out << Verbose(0) << "Begin of FragmentBOSSANOVA." << endl;
 
@@ -5023 +4454 @@
     Walker = Walker->next;
     CompleteMolecule.insert(Walker->GetTrueFather()->nr);
-  } 
+  }
 
   // this can easily be seen: if Order is 5, then the number of levels for each lower order is the total sum of the number of levels above, as
@@ -5037 +4468 @@
     //if ((MinimumRingSize[Walker->GetTrueFather()->nr] != -1) && (Walker->GetTrueFather()->AdaptiveOrder+1 > MinimumRingSize[Walker->GetTrueFather()->nr])) {
     //  *out << Verbose(0) << "Bond order " << Walker->GetTrueFather()->AdaptiveOrder << " of Root " << *Walker << " greater than or equal to Minimum Ring size of " << MinimumRingSize << " found is not allowed." << endl;
-    //} else 
+    //} else
     {
       // increase adaptive order by one
       Walker->GetTrueFather()->AdaptiveOrder++;
       Order = Walker->AdaptiveOrder = Walker->GetTrueFather()->AdaptiveOrder;
-  
+
       // initialise Order-dependent entries of UniqueFragments structure
       FragmentSearch.BondsPerSPList = (bond **) Malloc(sizeof(bond *)*Order*2, "molecule::PowerSetGenerator: ***BondsPerSPList");
@@ -5049 +4480 @@
         FragmentSearch.BondsPerSPList[2*i] = new bond();    // start node
         FragmentSearch.BondsPerSPList[2*i+1] = new bond();  // end node
-        FragmentSearch.BondsPerSPList[2*i]->next = FragmentSearch.BondsPerSPList[2*i+1];     // intertwine these two 
+        FragmentSearch.BondsPerSPList[2*i]->next = FragmentSearch.BondsPerSPList[2*i+1];     // intertwine these two
         FragmentSearch.BondsPerSPList[2*i+1]->previous = FragmentSearch.BondsPerSPList[2*i];
         FragmentSearch.BondsPerSPCount[i] = 0;
-      }  
-  
+      }
+
       // allocate memory for all lower level orders in this 1D-array of ptrs
       NumLevels = 1 << (Order-1); // (int)pow(2,Order);
@@ -5059 +4490 @@
       for (int i=0;i<NumLevels;i++)
         FragmentLowerOrdersList[RootNr][i] = NULL;
-      
+
       // create top order where nothing is reduced
       *out << Verbose(0) << "==============================================================================================================" << 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})
       FragmentLowerOrdersList[RootNr][0] =  new Graph;
@@ -5076 +4507 @@
         // 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
@@ -5083 +4514 @@
 //          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))    
+//            while (source >= (1 << (Walker->AdaptiveOrder-Order))) // (int)pow(2,Walker->AdaptiveOrder-Order))
 //              Order--;
 //            *out << Verbose(0) << "Current Order is: " << Order << "." << endl;
@@ -5090 +4521 @@
 //              *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;
@@ -5121 +4552 @@
       RootStack.push_back(RootKeyNr); // put back on stack
       RootNr++;
-  
+
       // free Order-dependent entries of UniqueFragments structure for next loop cycle
       Free((void **)&FragmentSearch.BondsPerSPCount, "molecule::PowerSetGenerator: *BondsPerSPCount");
@@ -5127 +4558 @@
         delete(FragmentSearch.BondsPerSPList[2*i]);
         delete(FragmentSearch.BondsPerSPList[2*i+1]);
-      }  
+      }
       Free((void **)&FragmentSearch.BondsPerSPList, "molecule::PowerSetGenerator: ***BondsPerSPList");
     }
@@ -5138 +4569 @@
   Free((void **)&FragmentSearch.ShortestPathList, "molecule::PowerSetGenerator: *ShortestPathList");
   delete(FragmentSearch.FragmentSet);
-  
-  // now, FragmentLowerOrdersList is complete, it looks - for BondOrder 5 - as this (number is the ANOVA Order of the terms therein) 
+
+  // now, FragmentLowerOrdersList is complete, it looks - for BondOrder 5 - as this (number is the ANOVA Order of the terms therein)
   // 5433222211111111
   // 43221111
@@ -5159 +4590 @@
     RootKeyNr = RootStack.front();
     RootStack.pop_front();
-    Walker = FindAtom(RootKeyNr); 
+    Walker = FindAtom(RootKeyNr);
     NumLevels = 1 << (Walker->AdaptiveOrder - 1);
     for(int i=0;i<NumLevels;i++) {
@@ -5172 +4603 @@
   Free((void **)&FragmentLowerOrdersList, "molecule::FragmentBOSSANOVA: ***FragmentLowerOrdersList");
   Free((void **)&NumMoleculesOfOrder, "molecule::FragmentBOSSANOVA: *NumMoleculesOfOrder");
-  
+
   *out << Verbose(0) << "End of FragmentBOSSANOVA." << endl;
 };
@@ -5206 +4637 @@
   atom *Walker = NULL;
   bool result = true; // status of comparison
-  
-  *out << Verbose(3) << "Begin of IsEqualToWithinThreshold." << endl; 
+
+  *out << Verbose(3) << "Begin of IsEqualToWithinThreshold." << endl;
   /// first count both their atoms and elements and update lists thereby ...
   //*out << Verbose(0) << "Counting atoms, updating list" << endl;
@@ -5254 +4685 @@
     if (CenterOfGravity.Distance(&OtherCenterOfGravity) > threshold) {
       *out << Verbose(4) << "Centers of gravity don't match." << endl;
-      result = false;  
-    }
-  }
-  
+      result = false;
+    }
+  }
+
   /// ... then make a list with the euclidian distance to this center for each atom of both molecules
   if (result) {
@@ -5273 +4704 @@
       OtherDistances[Walker->nr] = OtherCenterOfGravity.Distance(&Walker->x);
     }
-  
+
     /// ... sort each list (using heapsort (o(N log N)) from GSL)
     *out << Verbose(5) << "Sorting distances" << endl;
@@ -5284 +4715 @@
     for(int i=AtomCount;i--;)
       PermutationMap[PermMap[i]] = (int) OtherPermMap[i];
-    
+
     /// ... and compare them step by step, whether the difference is individiually(!) below \a threshold for all
     *out << Verbose(4) << "Comparing distances" << endl;
@@ -5295 +4726 @@
     Free((void **)&PermMap, "molecule::IsEqualToWithinThreshold: *PermMap");
     Free((void **)&OtherPermMap, "molecule::IsEqualToWithinThreshold: *OtherPermMap");
-      
+
     /// free memory
     Free((void **)&Distances, "molecule::IsEqualToWithinThreshold: Distances");
@@ -5303 +4734 @@
       result = false;
     }
-  } 
+  }
   /// return pointer to map if all distances were below \a threshold
   *out << Verbose(3) << "End of IsEqualToWithinThreshold." << endl;
@@ -5312 +4743 @@
     *out << Verbose(3) << "Result: Not equal." << endl;
     return NULL;
-  } 
+  }
 };
 
@@ -5367 +4798 @@
  * \param *output output stream of temperature file
  * \return file written (true), failure on writing file (false)
- */ 
+ */
 bool molecule::OutputTemperatureFromTrajectories(ofstream *out, int startstep, int endstep, ofstream *output)
 {
@@ -5375 +4806 @@
   if (output == NULL)
     return false;
-  else 
+  else
     *output << "# Step Temperature [K] Temperature [a.u.]" << endl;
   for (int step=startstep;step < endstep; step++) { // loop over all time steps