Index: src/builder.cpp
===================================================================
--- src/builder.cpp	(revision b4b7c37aea37c05e26de1669061a747720f5d760)
+++ src/builder.cpp	(revision d7e30ca04fa0adf2762a650a9f855a8ded1885f5)
@@ -147,8 +147,8 @@
           third->Output(1,2,(ofstream *)&cout);
           fourth->Output(1,3,(ofstream *)&cout);
-          n.MakeNormalVector((const vector *)&second->x, (const vector *)&third->x, (const vector *)&fourth->x);
-          x.CopyVector(&second->x);
+          n.MakeNormalvector((const vector *)&second->x, (const vector *)&third->x, (const vector *)&fourth->x);
+          x.Copyvector(&second->x);
           x.SubtractVector(&third->x);
-          x.CopyVector(&fourth->x);
+          x.Copyvector(&fourth->x);
           x.SubtractVector(&third->x);
           
@@ -359,5 +359,5 @@
       } while ((param.type = periode->FindElement(shorthand)) == NULL);
       cout << Verbose(0) << "Element is " << param.type->name << endl;
-      mol->GetAlignVector(&param);
+      mol->GetAlignvector(&param);
       for (int i=NDIM;i--;) {
         x.x[i] = gsl_vector_get(param.x,i);
@@ -761,5 +761,5 @@
             cout << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl;
             cout << "\t-e <file>\tSets the databases path to be parsed (default: ./)." << endl;
-            cout << "\t-f <dist> <order>\tFragments the molecule in BOSSANOVA manner and stores config files in same dir as config." << endl;
+            cout << "\t-f/F <dist> <order>\tFragments the molecule in BOSSANOVA manner (with/out rings compressed) and stores config files in same dir as config." << endl;
             cout << "\t-h/-H/-?\tGive this help screen." << endl;
             cout << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl;
@@ -767,4 +767,5 @@
             cout << "\t-o\tGet volume of the convex envelope (and store to tecplot file)." << endl;
             cout << "\t-p <file>\tParse given xyz file and create raw config file from it." << endl;
+            cout << "\t-P <file> <delta_t>\tParse given forces file and append as an MD step with width delta_t to config file via Verlet." << endl;
             cout << "\t-r\t\tConvert file from an old pcp syntax." << endl;
             cout << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << endl;
@@ -863,7 +864,8 @@
               if (!mol->AddXYZFile(argv[argptr]))
                 cout << Verbose(2) << "File not found." << endl;
-              else
+              else {
                 cout << Verbose(2) << "File found and parsed." << endl;
                 config_present = present;
+              }
               break;
             default:   // no match? Don't step on (this is done in next switch's default)
@@ -873,4 +875,13 @@
         if (config_present == present) {
           switch(argv[argptr-1][1]) {
+            case 'P':
+              ExitFlag = 1;
+              cout << Verbose(1) << "Parsing forces file and Verlet integrating." << endl;
+              if (!mol->VerletForceIntegration(argv[argptr], atof(argv[argptr+1])))
+                cout << Verbose(2) << "File not found." << endl;
+              else
+                cout << Verbose(2) << "File found and parsed." << endl;
+              argptr+=2;
+              break;
             case 't':
               ExitFlag = 1;
@@ -959,4 +970,5 @@
               cout << Verbose(1) << "Converting config file from supposed old to new syntax." << endl;
               break;
+            case 'F':
             case 'f':
               if (ExitFlag ==0) ExitFlag = 2;  // only set if not already by other command line switch
@@ -1022,5 +1034,5 @@
                 int count;
                 element ** Elements;
-                vector ** Vectors;
+                vector ** vectors;
                 if (faktor < 1) {
                   cerr << Verbose(0) << "ERROR: Repetition faktor mus be greater than 1!" << endl;
@@ -1031,5 +1043,5 @@
                   count = mol->AtomCount;   // is changed becausing of adding, thus has to be stored away beforehand
                   Elements = new element *[count];
-                  Vectors = new vector *[count];
+                  vectors = new vector *[count];
                   j = 0;
                   first = mol->start;
@@ -1037,5 +1049,5 @@
                     first = first->next;
                     Elements[j] = first->type;
-                    Vectors[j] = &first->x;
+                    vectors[j] = &first->x;
                     j++;
                   }
@@ -1049,5 +1061,5 @@
                     for (int k=count;k--;) { // go through every atom of the original cell
                       first = new atom(); // create a new body
-                      first->x.CopyVector(Vectors[k]);  // use coordinate of original atom
+                      first->x.CopyVector(vectors[k]);  // use coordinate of original atom
                       first->x.AddVector(&x);      // translate the coordinates
                       first->type = Elements[k];  // insert original element
@@ -1057,5 +1069,5 @@
                   // free memory
                   delete[](Elements);
-                  delete[](Vectors);
+                  delete[](vectors);
                   // correct cell size
                   if (axis < 0) { // if sign was negative, we have to translate everything
@@ -1120,5 +1132,5 @@
   clock_t start,end;
   element **Elements;
-  vector **Vectors;
+  vector **vectors;
 
   // =========================== PARSE COMMAND LINE OPTIONS ====================================
@@ -1219,5 +1231,5 @@
           count = mol->AtomCount;   // is changed becausing of adding, thus has to be stored away beforehand
           Elements = new element *[count];
-          Vectors = new vector *[count];
+          vectors = new vector *[count];
           j = 0;
           first = mol->start;
@@ -1225,5 +1237,5 @@
             first = first->next;
             Elements[j] = first->type;
-            Vectors[j] = &first->x;
+            vectors[j] = &first->x;
             j++;
           }
@@ -1237,5 +1249,5 @@
             for (int k=count;k--;) { // go through every atom of the original cell
               first = new atom(); // create a new body
-              first->x.CopyVector(Vectors[k]);  // use coordinate of original atom
+              first->x.CopyVector(vectors[k]);  // use coordinate of original atom
               first->x.AddVector(&x);      // translate the coordinates
               first->type = Elements[k];  // insert original element
@@ -1247,5 +1259,5 @@
           // free memory
           delete[](Elements);
-          delete[](Vectors);
+          delete[](vectors);
           // correct cell size
           if (axis < 0) { // if sign was negative, we have to translate everything
Index: src/molecules.cpp
===================================================================
--- src/molecules.cpp	(revision b4b7c37aea37c05e26de1669061a747720f5d760)
+++ src/molecules.cpp	(revision d7e30ca04fa0adf2762a650a9f855a8ded1885f5)
@@ -140,5 +140,5 @@
  * -# Single Bond: Simply add new atom with bond distance rescaled to typical hydrogen one
  * -# Double Bond: Here, we need the **BondList of the \a *origin atom, by scanning for the other bonds instead of
- *    *Bond, we use the through these connected atoms to determine the plane they lie in, vector::MakeNormalVector().
+ *    *Bond, we use the through these connected atoms to determine the plane they lie in, vector::MakeNormalvector().
  *    The orthonormal vector to this plane along with the vector in *Bond direction determines the plane the two
  *    replacing hydrogens shall lie in. Now, all remains to do is take the usual hydrogen double bond angle for the
@@ -151,5 +151,5 @@
  *    \f[ h = l \cdot \cos{\left (\frac{\alpha}{2} \right )} \qquad b = 2l \cdot \sin{\left (\frac{\alpha}{2} \right)} \quad \rightarrow \quad d = l \cdot \sqrt{\cos^2{\left (\frac{\alpha}{2} \right)}-\frac{1}{3}\cdot\sin^2{\left (\frac{\alpha}{2}\right )}}
  *    \f]
- *    vector::GetNormalVector() creates one orthonormal vector from this *Bond vector and vector::MakeNormalVector creates
+ *    vector::GetNormalvector() creates one orthonormal vector from this *Bond vector and vector::MakeNormalvector creates
  *    the third one from the former two vectors. The latter ones form the plane of the base triangle mentioned above.
  *    The lengths for these are \f$f\f$ and \f$g\f$ (from triangle H1-H2-(center of H1-H2-H3)) with knowledge that
@@ -180,6 +180,6 @@
   atom *FirstOtherAtom = NULL, *SecondOtherAtom = NULL, *ThirdOtherAtom = NULL; // pointer to hydrogen atoms to be added
   double b,l,d,f,g, alpha, factors[NDIM];    // hold temporary values in triple bond case for coordination determination
-  vector OrthoVector1, OrthoVector2;  // temporary vectors in coordination construction
-  vector InBondVector;    // vector in direction of *Bond
+  vector Orthovector1, Orthovector2;  // temporary vectors in coordination construction
+  vector InBondvector;    // vector in direction of *Bond
   bond *Binder = NULL;
   double *matrix;
@@ -187,33 +187,33 @@
 //	*out << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl;
   // create vector in direction of bond
-  InBondVector.CopyVector(&TopReplacement->x);
-  InBondVector.SubtractVector(&TopOrigin->x);
-  bondlength = InBondVector.Norm();
+  InBondvector.CopyVector(&TopReplacement->x);
+  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
+   // 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) { 
-//    *out << Verbose(4) << "InBondVector is: ";
-//    InBondVector.Output(out);
+//    *out << Verbose(4) << "InBondvector is: ";
+//    InBondvector.Output(out);
 //    *out << endl;
-    OrthoVector1.Zero();
+    Orthovector1.Zero();
     for (int i=NDIM;i--;) {
       l = TopReplacement->x.x[i] - TopOrigin->x.x[i];
       if (fabs(l) > BondDistance) { // is component greater than bond distance
-        OrthoVector1.x[i] = (l < 0) ? -1. : +1.;
+        Orthovector1.x[i] = (l < 0) ? -1. : +1.;
       } // (signs are correct, was tested!)
     }
     matrix = ReturnFullMatrixforSymmetric(cell_size);
-    OrthoVector1.MatrixMultiplication(matrix);
-    InBondVector.SubtractVector(&OrthoVector1); // subtract just the additional translation
+    Orthovector1.MatrixMultiplication(matrix);
+    InBondvector.SubtractVector(&Orthovector1); // subtract just the additional translation
     Free((void **)&matrix, "molecule::AddHydrogenReplacementAtom: *matrix");
-    bondlength = InBondVector.Norm();
-//    *out << Verbose(4) << "Corrected InBondVector is now: ";
-//    InBondVector.Output(out);
+    bondlength = InBondvector.Norm();
+//    *out << Verbose(4) << "Corrected InBondvector is now: ";
+//    InBondvector.Output(out);
 //    *out << endl;
   } // periodic correction finished
   
-  InBondVector.Normalize();
+  InBondvector.Normalize();
   // get typical bond length and store as scale factor for later
   BondRescale = TopOrigin->type->HBondDistance[TopBond->BondDegree-1];
@@ -239,7 +239,7 @@
         FirstOtherAtom->father = NULL;  // if we replace hydrogen, we mark it as our father, otherwise we are just an added hydrogen with no father
       }
-      InBondVector.Scale(&BondRescale);   // rescale the distance vector to Hydrogen bond length
+      InBondvector.Scale(&BondRescale);   // rescale the distance vector to Hydrogen bond length
       FirstOtherAtom->x.CopyVector(&TopOrigin->x); // set coordination to origin ...
-      FirstOtherAtom->x.AddVector(&InBondVector);  // ... and add distance vector to replacement atom
+      FirstOtherAtom->x.AddVector(&InBondvector);  // ... and add distance vector to replacement atom
       AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
 //      *out << Verbose(4) << "Added " << *FirstOtherAtom << " at: ";
@@ -273,15 +273,15 @@
         
         // determine the plane of these two with the *origin
-        AllWentWell = AllWentWell && OrthoVector1.MakeNormalVector(&TopOrigin->x, &FirstOtherAtom->x, &SecondOtherAtom->x);
+        AllWentWell = AllWentWell && Orthovector1.MakeNormalVector(&TopOrigin->x, &FirstOtherAtom->x, &SecondOtherAtom->x);
       } else {
-        OrthoVector1.GetOneNormalVector(&InBondVector);
+        Orthovector1.GetOneNormalVector(&InBondvector);
       }
       //*out << Verbose(3)<< "Orthovector1: ";
-      //OrthoVector1.Output(out);
+      //Orthovector1.Output(out);
       //*out << endl;
       // orthogonal vector and bond vector between origin and replacement form the new plane
-      OrthoVector1.MakeNormalVector(&InBondVector);
-      OrthoVector1.Normalize();
-      //*out << Verbose(3) << "ReScaleCheck: " << OrthoVector1.Norm() << " and " << InBondVector.Norm() << "." << endl;
+      Orthovector1.MakeNormalVector(&InBondvector);
+      Orthovector1.Normalize();
+      //*out << Verbose(3) << "ReScaleCheck: " << Orthovector1.Norm() << " and " << InBondvector.Norm() << "." << endl;
       
       // create the two Hydrogens ...
@@ -302,16 +302,16 @@
       }
       bondangle *= M_PI/180./2.;
-//      *out << Verbose(3) << "ReScaleCheck: InBondVector ";
-//      InBondVector.Output(out);
+//      *out << Verbose(3) << "ReScaleCheck: InBondvector ";
+//      InBondvector.Output(out);
 //      *out << endl;
 //      *out << Verbose(3) << "ReScaleCheck: Orthovector ";
-//      OrthoVector1.Output(out);
+//      Orthovector1.Output(out);
 //      *out << endl;
 //      *out << Verbose(3) << "Half the bond angle is " << bondangle << ", sin and cos of it: " << sin(bondangle) << ", " << cos(bondangle) << endl;
       FirstOtherAtom->x.Zero();
       SecondOtherAtom->x.Zero();
-      for(int i=NDIM;i--;) { // rotate by half the bond angle in both directions (InBondVector is bondangle = 0 direction)
-        FirstOtherAtom->x.x[i] = InBondVector.x[i] * cos(bondangle) + OrthoVector1.x[i] * (sin(bondangle));
-        SecondOtherAtom->x.x[i] = InBondVector.x[i] * cos(bondangle) + OrthoVector1.x[i] * (-sin(bondangle));
+      for(int i=NDIM;i--;) { // rotate by half the bond angle in both directions (InBondvector is bondangle = 0 direction)
+        FirstOtherAtom->x.x[i] = InBondvector.x[i] * cos(bondangle) + Orthovector1.x[i] * (sin(bondangle));
+        SecondOtherAtom->x.x[i] = InBondvector.x[i] * cos(bondangle) + Orthovector1.x[i] * (-sin(bondangle));
       }
       FirstOtherAtom->x.Scale(&BondRescale);  // rescale by correct BondDistance 
@@ -356,12 +356,12 @@
       ThirdOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
 
-      // we need to vectors orthonormal the InBondVector
-      AllWentWell = AllWentWell && OrthoVector1.GetOneNormalVector(&InBondVector);
+      // we need to vectors orthonormal the InBondvector
+      AllWentWell = AllWentWell && Orthovector1.GetOneNormalVector(&InBondvector);
 //      *out << Verbose(3) << "Orthovector1: ";
-//      OrthoVector1.Output(out);
+//      Orthovector1.Output(out);
 //      *out << endl;
-      AllWentWell = AllWentWell && OrthoVector2.MakeNormalVector(&InBondVector, &OrthoVector1);
+      AllWentWell = AllWentWell && Orthovector2.MakeNormalVector(&InBondvector, &Orthovector1);
 //      *out << Verbose(3) << "Orthovector2: "; 
-//      OrthoVector2.Output(out);
+//      Orthovector2.Output(out);
 //      *out << endl;
       
@@ -371,6 +371,6 @@
       b = 2.*l*sin(alpha);    // base length of isosceles triangle
       d = l*sqrt(cos(alpha)*cos(alpha) - sin(alpha)*sin(alpha)/3.);   // length for InBondvector
-      f = b/sqrt(3.);   // length for OrthVector1
-      g = b/2.;         // length for OrthVector2
+      f = b/sqrt(3.);   // length for Orthvector1
+      g = b/2.;         // length for Orthvector2
 //      *out << Verbose(3) << "Bond length and half-angle: " << l << ", " << alpha << "\t (b,d,f,g) = " << b << ", " << d << ", " << f << ", " << g << ", " << endl;
 //      *out << Verbose(3) << "The three Bond lengths: " << sqrt(d*d+f*f) << ", " << sqrt(d*d+(-0.5*f)*(-0.5*f)+g*g) << ", "  << sqrt(d*d+(-0.5*f)*(-0.5*f)+g*g) << endl;
@@ -378,10 +378,10 @@
       factors[1] = f;
       factors[2] = 0.; 
-      FirstOtherAtom->x.LinearCombinationOfVectors(&InBondVector, &OrthoVector1, &OrthoVector2, factors);
+      FirstOtherAtom->x.LinearCombinationOfVectors(&InBondvector, &Orthovector1, &Orthovector2, factors);
       factors[1] = -0.5*f;
       factors[2] = g; 
-      SecondOtherAtom->x.LinearCombinationOfVectors(&InBondVector, &OrthoVector1, &OrthoVector2, factors);
+      SecondOtherAtom->x.LinearCombinationOfVectors(&InBondvector, &Orthovector1, &Orthovector2, factors);
       factors[2] = -g; 
-      ThirdOtherAtom->x.LinearCombinationOfVectors(&InBondVector, &OrthoVector1, &OrthoVector2, factors);
+      ThirdOtherAtom->x.LinearCombinationOfVectors(&InBondvector, &Orthovector1, &Orthovector2, factors);
 
       // rescale each to correct BondDistance
@@ -802,5 +802,5 @@
   double tmp;
   bool flag;
-  vector TestVector, TranslationVector;
+  vector Testvector, Translationvector;
   
   do {
@@ -812,7 +812,7 @@
       if (Walker->type->Z != 1) {
 #endif
-        TestVector.CopyVector(&Walker->x);
-        TestVector.InverseMatrixMultiplication(matrix);
-        TranslationVector.Zero();
+        Testvector.CopyVector(&Walker->x);
+        Testvector.InverseMatrixMultiplication(matrix);
+        Translationvector.Zero();
         for (int i=0;i<NumberOfBondsPerAtom[Walker->nr]; i++) {
           Binder = ListOfBondsPerAtom[Walker->nr][i];
@@ -824,15 +824,15 @@
                 cout << Verbose(0) << "Hit: atom " << Walker->Name << " in bond " << *Binder << " has to be shifted due to " << tmp << "." << endl;
                 if (tmp > 0)
-                  TranslationVector.x[j] -= 1.;
+                  Translationvector.x[j] -= 1.;
                 else
-                  TranslationVector.x[j] += 1.;
+                  Translationvector.x[j] += 1.;
               }
             }
         }
-        TestVector.AddVector(&TranslationVector);
-        TestVector.MatrixMultiplication(matrix);
-        Center.AddVector(&TestVector);
-        cout << Verbose(1) << "Vector is: ";
-        TestVector.Output((ofstream *)&cout);
+        Testvector.AddVector(&Translationvector);
+        Testvector.MatrixMultiplication(matrix);
+        Center.AddVector(&Testvector);
+        cout << Verbose(1) << "vector is: ";
+        Testvector.Output((ofstream *)&cout);
         cout << endl;     
 #ifdef ADDHYDROGEN
@@ -841,11 +841,11 @@
           Binder = ListOfBondsPerAtom[Walker->nr][i];
           if (Binder->GetOtherAtom(Walker)->type->Z == 1) {
-            TestVector.CopyVector(&Binder->GetOtherAtom(Walker)->x);
-            TestVector.InverseMatrixMultiplication(matrix);
-            TestVector.AddVector(&TranslationVector);
-            TestVector.MatrixMultiplication(matrix);
-            Center.AddVector(&TestVector);
-            cout << Verbose(1) << "Hydrogen Vector is: ";
-            TestVector.Output((ofstream *)&cout);
+            Testvector.CopyVector(&Binder->GetOtherAtom(Walker)->x);
+            Testvector.InverseMatrixMultiplication(matrix);
+            Testvector.AddVector(&Translationvector);
+            Testvector.MatrixMultiplication(matrix);
+            Center.AddVector(&Testvector);
+            cout << Verbose(1) << "Hydrogen vector is: ";
+            Testvector.Output((ofstream *)&cout);
             cout << endl;     
           }
@@ -963,4 +963,98 @@
 };
 
+/** Parses nuclear forces from file and performs Verlet integration.
+ * This adds a new MD step to the config file.
+ * \param *file filename
+ * \param delta_t time step width in atomic units
+ * \return true - file found and parsed, false - file not found or imparsable
+ */
+bool molecule::VerletForceIntegration(char *file, double delta_t)
+{
+  bool status = true;
+  element *runner = elemente->start;
+  atom *walker = NULL;
+  int ElementNo, AtomNo;
+  ifstream input(file);
+  double Forcesvector[3];
+  double dummy;
+  string token;
+  size_t oldpos, newpos;
+  stringstream item;
+  int i, Ion[2];
+  double a, IonMass;
+
+  CountElements();  // make sure ElementsInMolecule is up to date
+  
+  // check file
+  if (input == NULL) {
+    return false;
+  } else {
+    ElementNo = 0;
+    while (runner->next != elemente->end) { // go through every element
+      runner = runner->next;
+      IonMass = runner->mass;
+      if (ElementsInMolecule[runner->Z]) { // if this element got atoms
+        ElementNo++;
+        AtomNo = 0;
+        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
+            AtomNo++;
+            // parse Forcevector for this ion
+            getline (input, token, '\n');
+            newpos = oldpos = i = 0;
+            while ((newpos = token.find( '\t', oldpos)) != oldpos) {
+              i++;
+              item.str();
+              switch (i) {
+                case 1:
+                case 2:
+                  item >> Ion[i-1];
+                  break;
+                case 6:
+                case 7:
+                case 8:
+                  item >> Forcesvector[i-6];
+                  break;
+                case 3:
+                  if ((Ion[0] != ElementNo) || (Ion[1] != AtomNo))
+                    cout << "ERROR: Expected " << ElementNo << ", " << AtomNo << "but parsed " << Ion[0] << ", " << Ion[1] << "." << endl;
+                  // still parse into dummy! Hence no break here ...
+                default:
+                  item >> dummy;
+              }
+              oldpos = newpos;
+            }
+            
+            // perform Verlet integration for this atom with position, velocity and force vector
+            
+            // UpdateR
+//            for (int d=0; d<NDIM; d++) {
+//              R_old_old[d] = R_old[d];        // shift old values
+//              R_old[d] = R[d];
+//              R[d] += delta_t*(U[d]+a*FIon[d]);     // F = m * a and s = 0.5 * a * t^2
+//              FIon_old[d] = FIon[d];          // store old force
+//            }
+//            // UpdateU
+//            a = delta_t*0.5/IonMass;        // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
+//            for (int d=0; d<NDIM; d++) {
+//                U[d] += a * (FIon[d]+FIon_old[d]);
+//            }
+            
+            
+            // add trajectory point
+
+          }
+        }
+      }
+    }
+    return true;
+  }
+  
+  // exit
+  return status;
+};
+
 /** Align all atoms in such a manner that given vector \a *n is along z axis. 
  * \param n[] alignment vector.
@@ -1124,5 +1218,5 @@
  * \bug this is not yet working properly it seems
  */
-void molecule::GetAlignVector(struct lsq_params * par) const
+void molecule::GetAlignvector(struct lsq_params * par) const
 {
     int np = 6;
@@ -3573,6 +3667,6 @@
           *out << (*runner) << " "; 
         *out << endl;
-        if (!CheckForConnectedSubgraph(out, FragmentSearch->FragmentSet))
-          *out << Verbose(0) << "ERROR: The found fragment is not a connected subgraph!" << endl;
+        //if (!CheckForConnectedSubgraph(out, FragmentSearch->FragmentSet))
+          //*out << Verbose(0) << "ERROR: The found fragment is not a connected subgraph!" << endl;
         InsertFragmentIntoGraph(out, FragmentSearch);
         //Removal = LookForRemovalCandidate(out, FragmentSearch->FragmentSet, FragmentSearch->ShortestPathList);
@@ -3831,5 +3925,5 @@
   enum Shading *ColorList = NULL;
   double tmp;
-  vector TranslationVector;
+  vector Translationvector;
   //class StackClass<atom *> *CompStack = NULL;
   class StackClass<atom *> *AtomStack = new StackClass<atom *>(AtomCount);
@@ -3842,5 +3936,5 @@
     // remove bonds that are beyond bonddistance
     for(int i=NDIM;i--;)
-      TranslationVector.x[i] = 0.; 
+      Translationvector.x[i] = 0.; 
     // scan all bonds
     Binder = first;
@@ -3865,9 +3959,9 @@
         tmp = Binder->leftatom->x.x[i] - Binder->rightatom->x.x[i];
         if (fabs(tmp) > BondDistance)
-          TranslationVector.x[i] = (tmp < 0) ? +1. : -1.;
-      }
-      TranslationVector.MatrixMultiplication(matrix);
+          Translationvector.x[i] = (tmp < 0) ? +1. : -1.;
+      }
+      Translationvector.MatrixMultiplication(matrix);
       *out << Verbose(3) << "Translation vector is ";
-      TranslationVector.Output(out);
+      Translationvector.Output(out);
       *out << endl;
       // apply to all atoms of first component via BFS
@@ -3879,5 +3973,5 @@
         //*out << Verbose (3) << "Current Walker is: " << *Walker << "." << endl;
         ColorList[Walker->nr] = black;    // mark as explored
-        Walker->x.AddVector(&TranslationVector); // translate
+        Walker->x.AddVector(&Translationvector); // translate
         for (int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {  // go through all binding partners
           if (ListOfBondsPerAtom[Walker->nr][i] != Binder) {
Index: src/molecules.hpp
===================================================================
--- src/molecules.hpp	(revision b4b7c37aea37c05e26de1669061a747720f5d760)
+++ src/molecules.hpp	(revision d7e30ca04fa0adf2762a650a9f855a8ded1885f5)
@@ -20,4 +20,5 @@
 #include <set>
 #include <deque>
+#include <list>
 
 #include "helpers.hpp"
@@ -46,4 +47,11 @@
 {
   bool operator() (const KeySet SubgraphA, const KeySet SubgraphB) const;
+};
+
+struct Trajectory
+{
+  double R[3];  //!< position vector
+  double U[3];    //!< velocity vector
+  double F[3];    //!< last force vector
 };
 
@@ -219,7 +227,8 @@
 	void PrincipalAxisSystem(ofstream *out, bool DoRotate);
 	double VolumeOfConvexEnvelope(ofstream *out, bool IsAngstroem);
+	bool VerletForceIntegration(char *file, double delta_t);
 	
   bool CheckBounds(const vector *x) const;
-  void GetAlignVector(struct lsq_params * par) const;
+  void GetAlignvector(struct lsq_params * par) const;
 
   /// Initialising routines in fragmentation  
