Index: src/Actions/CommandAction/BondLengthTableAction.cpp
===================================================================
--- src/Actions/CommandAction/BondLengthTableAction.cpp	(revision fd378853756ce2d629fc774e7c82d6739ae1be78)
+++ src/Actions/CommandAction/BondLengthTableAction.cpp	(revision 4e855e867fa84e16d3c418b917eb36c8f1b45f90)
@@ -45,17 +45,26 @@
   DoLog(0) && (Log() << Verbose(0) << "Using " << params.BondGraphFileName << " as bond length table." << endl);
   config *configuration = World::getInstance().getConfig();
-  if (configuration->BG == NULL) {
-    configuration->BG = new BondGraph(configuration->GetIsAngstroem());
-    if ((!params.BondGraphFileName.empty())
-        && boost::filesystem::exists(params.BondGraphFileName)
-        && (configuration->BG->LoadBondLengthTable(params.BondGraphFileName.string()))) {
-      DoLog(0) && (Log() << Verbose(0) << "Bond length table loaded successfully." << endl);
-      return Action::success;
+  BondGraph *OldBG;
+  if (configuration->BG != NULL) {
+    OldBG = configuration->BG;
+    DoLog(0) && (Log() << Verbose(0) << "There is a bond length table already present." << endl);
+  }
+  configuration->BG = new BondGraph(configuration->GetIsAngstroem());
+  if ((!params.BondGraphFileName.empty())
+      && boost::filesystem::exists(params.BondGraphFileName)) {
+    std::ifstream input(params.BondGraphFileName.string().c_str());
+    if ((input.good()) && (configuration->BG->LoadBondLengthTable(input))) {
+      DoLog(0) && (Log() << Verbose(0) << "Bond length table parsed successfully." << endl);
+      input.close();
+      return Action::state_ptr(new CommandBondLengthTableState(*OldBG, params));
     } else {
-      DoeLog(1) && (eLog()<< Verbose(1) << "Bond length table loading failed." << endl);
+      DoeLog(1) && (eLog()<< Verbose(1) << "Bond length table parsing failed." << endl);
+      input.close();
+      configuration->BG = new BondGraph(*OldBG);
       return Action::failure;
     }
   } else {
-    DoLog(0) && (Log() << Verbose(0) << "Bond length table already present." << endl);
+    DoeLog(1) && (eLog()<< Verbose(1) << "Bond length table loading failed." << endl);
+    configuration->BG = new BondGraph(*OldBG);
     return Action::failure;
   }
@@ -63,15 +72,21 @@
 
 Action::state_ptr CommandBondLengthTableAction::performUndo(Action::state_ptr _state) {
-//  ParserLoadXyzState *state = assert_cast<ParserLoadXyzState*>(_state.get());
+  CommandBondLengthTableState *state = assert_cast<CommandBondLengthTableState*>(_state.get());
 
-  return Action::failure;
-//  string newName = state->mol->getName();
-//  state->mol->setName(state->lastName);
-//
-//  return Action::state_ptr(new ParserLoadXyzState(state->mol,newName));
+  config *configuration = World::getInstance().getConfig();
+  BondGraph *OldBG = new BondGraph(*configuration->BG);
+  configuration->BG = new BondGraph(state->OldBG);
+
+  return Action::state_ptr(new CommandBondLengthTableState(*OldBG,params));
 }
 
 Action::state_ptr CommandBondLengthTableAction::performRedo(Action::state_ptr _state){
-  return Action::failure;
+  CommandBondLengthTableState *state = assert_cast<CommandBondLengthTableState*>(_state.get());
+
+  config *configuration = World::getInstance().getConfig();
+  BondGraph *OldBG = new BondGraph(*configuration->BG);
+  configuration->BG = new BondGraph(state->OldBG);
+
+  return Action::state_ptr(new CommandBondLengthTableState(*OldBG,params));
 }
 
Index: src/Actions/CommandAction/BondLengthTableAction.def
===================================================================
--- src/Actions/CommandAction/BondLengthTableAction.def	(revision fd378853756ce2d629fc774e7c82d6739ae1be78)
+++ src/Actions/CommandAction/BondLengthTableAction.def	(revision 4e855e867fa84e16d3c418b917eb36c8f1b45f90)
@@ -8,4 +8,5 @@
 // all includes and forward declarations necessary for non-integral types below
 
+#include "bondgraph.hpp"
 
 // i.e. there is an integer with variable name Z that can be found in
@@ -18,6 +19,6 @@
 #define paramreferences (BondGraphFileName)
 
-#undef statetypes
-#undef statereferences
+#define statetypes (BondGraph)
+#define statereferences (OldBG)
 
 // some defines for all the names, you may use ACTION, STATE and PARAMS
Index: src/Actions/FragmentationAction/DepthFirstSearchAction.cpp
===================================================================
--- src/Actions/FragmentationAction/DepthFirstSearchAction.cpp	(revision fd378853756ce2d629fc774e7c82d6739ae1be78)
+++ src/Actions/FragmentationAction/DepthFirstSearchAction.cpp	(revision 4e855e867fa84e16d3c418b917eb36c8f1b45f90)
@@ -52,5 +52,9 @@
   std::deque<bond *> *BackEdgeStack = NULL;
   std::deque<bond *> *LocalBackEdgeStack = NULL;
-  mol->CreateAdjacencyList(params.distance, World::getInstance().getConfig()->GetIsAngstroem(), &BondGraph::CovalentMinMaxDistance, NULL);
+  BondGraph *BG = World::getInstance().getConfig()->BG;
+  if (BG != NULL)
+    mol->CreateAdjacencyList(params.distance, World::getInstance().getConfig()->GetIsAngstroem(), &BondGraph::BondLengthMatrixMinMaxDistance, BG);
+  else
+    mol->CreateAdjacencyList(params.distance, World::getInstance().getConfig()->GetIsAngstroem(), &BondGraph::CovalentMinMaxDistance, BG);
   Subgraphs = mol->DepthFirstSearchAnalysis(BackEdgeStack);
   if (Subgraphs != NULL) {
Index: src/bondgraph.cpp
===================================================================
--- src/bondgraph.cpp	(revision fd378853756ce2d629fc774e7c82d6739ae1be78)
+++ src/bondgraph.cpp	(revision 4e855e867fa84e16d3c418b917eb36c8f1b45f90)
@@ -57,9 +57,8 @@
  * Allocates \a MatrixContainer for BondGraph::BondLengthMatrix, using MatrixContainer::ParseMatrix(),
  * but only if parsing is successful. Otherwise variable is left as NULL.
- * \param *out output stream for debugging
- * \param filename file with bond lengths to parse
+ * \param &input input stream to parse table from
  * \return true - success in parsing file, false - failed to parse the file
  */
-bool BondGraph::LoadBondLengthTable(const string &filename)
+bool BondGraph::LoadBondLengthTable(std::istream &input)
 {
   Info FunctionInfo(__func__);
@@ -75,8 +74,9 @@
 
   // parse in matrix
-  if ((status = TempContainer->ParseMatrix(filename.c_str(), 0, 1, 0))) {
+  if ((input.good()) && (status = TempContainer->ParseMatrix(input, 0, 1, 0))) {
     DoLog(1) && (Log() << Verbose(1) << "Parsing bond length matrix successful." << endl);
   } else {
     DoeLog(1) && (eLog()<< Verbose(1) << "Parsing bond length matrix failed." << endl);
+    status = false;
   }
 
@@ -129,8 +129,12 @@
 double BondGraph::GetBondLength(int firstZ, int secondZ)
 {
-  if (BondLengthMatrix == NULL)
+  std::cout << "Request for length between " << firstZ << " and " << secondZ << ": ";
+  if (BondLengthMatrix == NULL) {
+    std::cout << "-1." << std::endl;
     return( -1. );
-  else
+  } else {
+    std::cout << BondLengthMatrix->Matrix[0][firstZ][secondZ] << std::endl;
     return (BondLengthMatrix->Matrix[0][firstZ][secondZ]);
+  }
 };
 
Index: src/bondgraph.hpp
===================================================================
--- src/bondgraph.hpp	(revision fd378853756ce2d629fc774e7c82d6739ae1be78)
+++ src/bondgraph.hpp	(revision 4e855e867fa84e16d3c418b917eb36c8f1b45f90)
@@ -36,5 +36,5 @@
   BondGraph(bool IsA);
   ~BondGraph();
-  bool LoadBondLengthTable(const string &filename);
+  bool LoadBondLengthTable(std::istream &input);
   bool ConstructBondGraph(molecule * const mol);
   double GetBondLength(int firstelement, int secondelement);
Index: src/builder_init.cpp
===================================================================
--- src/builder_init.cpp	(revision fd378853756ce2d629fc774e7c82d6739ae1be78)
+++ src/builder_init.cpp	(revision 4e855e867fa84e16d3c418b917eb36c8f1b45f90)
@@ -19,4 +19,6 @@
 
 #include "CodePatterns/MemDebug.hpp"
+
+#include <iostream>
 
 #include "bondgraph.hpp"
@@ -107,9 +109,11 @@
     World::getInstance().getConfig()->BG = new BondGraph(World::getInstance().getConfig()->GetIsAngstroem());
     if (boost::filesystem::exists(BondGraphFileName)) {
-      if (World::getInstance().getConfig()->BG->LoadBondLengthTable(BondGraphFileName)) {
+      std::ifstream input(BondGraphFileName.c_str());
+      if ((input.good()) && (World::getInstance().getConfig()->BG->LoadBondLengthTable(input))) {
         DoLog(0) && (Log() << Verbose(0) << "Bond length table loaded successfully." << endl);
       } else {
         DoeLog(1) && (eLog()<< Verbose(1) << "Bond length table loading failed." << endl);
       }
+      input.close();
     }
   }
Index: src/molecule_dynamics.cpp
===================================================================
--- src/molecule_dynamics.cpp	(revision fd378853756ce2d629fc774e7c82d6739ae1be78)
+++ src/molecule_dynamics.cpp	(revision 4e855e867fa84e16d3c418b917eb36c8f1b45f90)
@@ -574,5 +574,4 @@
 {
   Info FunctionInfo(__func__);
-  ifstream input(file);
   string token;
   stringstream item;
@@ -582,44 +581,42 @@
 
   const int AtomCount = getAtomCount();
-  // check file
-  if (input == NULL) {
+  // parse file into ForceMatrix
+  std::ifstream input(file);
+  if ((input.good()) && (!Force.ParseMatrix(input, 0,0,0))) {
+    DoeLog(0) && (eLog()<< Verbose(0) << "Could not parse Force Matrix file " << file << "." << endl);
+    performCriticalExit();
     return false;
-  } else {
-    // parse file into ForceMatrix
-    if (!Force.ParseMatrix(file, 0,0,0)) {
-      DoeLog(0) && (eLog()<< Verbose(0) << "Could not parse Force Matrix file " << file << "." << endl);
-      performCriticalExit();
-      return false;
-    }
-    if (Force.RowCounter[0] != AtomCount) {
-      DoeLog(0) && (eLog()<< Verbose(0) << "Mismatch between number of atoms in file " << Force.RowCounter[0] << " and in molecule " << getAtomCount() << "." << endl);
-      performCriticalExit();
-      return false;
-    }
-    // correct Forces
-    Velocity.Zero();
-    for(int i=0;i<AtomCount;i++)
-      for(int d=0;d<NDIM;d++) {
-        Velocity[d] += Force.Matrix[0][i][d+offset];
-      }
-    for(int i=0;i<AtomCount;i++)
-      for(int d=0;d<NDIM;d++) {
-        Force.Matrix[0][i][d+offset] -= Velocity[d]/static_cast<double>(AtomCount);
-      }
-    // solve a constrained potential if we are meant to
-    if (configuration.DoConstrainedMD) {
-      // calculate forces and potential
-      atom **PermutationMap = NULL;
-      ConstrainedPotentialEnergy = MinimiseConstrainedPotential(PermutationMap,configuration.DoConstrainedMD, 0, configuration.GetIsAngstroem());
-      EvaluateConstrainedForces(configuration.DoConstrainedMD, 0, PermutationMap, &Force);
-      delete[](PermutationMap);
-    }
-
-    // and perform Verlet integration for each atom with position, velocity and force vector
-    // check size of vectors
-    for_each(atoms.begin(),
-             atoms.end(),
-             boost::bind(&atom::VelocityVerletUpdate,_1,MDSteps+1, &configuration, &Force, (const size_t) 0));
-  }
+  }
+  input.close();
+  if (Force.RowCounter[0] != AtomCount) {
+    DoeLog(0) && (eLog()<< Verbose(0) << "Mismatch between number of atoms in file " << Force.RowCounter[0] << " and in molecule " << getAtomCount() << "." << endl);
+    performCriticalExit();
+    return false;
+  }
+  // correct Forces
+  Velocity.Zero();
+  for(int i=0;i<AtomCount;i++)
+    for(int d=0;d<NDIM;d++) {
+      Velocity[d] += Force.Matrix[0][i][d+offset];
+    }
+  for(int i=0;i<AtomCount;i++)
+    for(int d=0;d<NDIM;d++) {
+      Force.Matrix[0][i][d+offset] -= Velocity[d]/static_cast<double>(AtomCount);
+    }
+  // solve a constrained potential if we are meant to
+  if (configuration.DoConstrainedMD) {
+    // calculate forces and potential
+    atom **PermutationMap = NULL;
+    ConstrainedPotentialEnergy = MinimiseConstrainedPotential(PermutationMap,configuration.DoConstrainedMD, 0, configuration.GetIsAngstroem());
+    EvaluateConstrainedForces(configuration.DoConstrainedMD, 0, PermutationMap, &Force);
+    delete[](PermutationMap);
+  }
+
+  // and perform Verlet integration for each atom with position, velocity and force vector
+  // check size of vectors
+  for_each(atoms.begin(),
+           atoms.end(),
+           boost::bind(&atom::VelocityVerletUpdate,_1,MDSteps+1, &configuration, &Force, (const size_t) 0));
+
   // correct velocities (rather momenta) so that center of mass remains motionless
   Velocity = atoms.totalMomentumAtStep(MDSteps+1);
Index: src/molecule_graph.cpp
===================================================================
--- src/molecule_graph.cpp	(revision fd378853756ce2d629fc774e7c82d6739ae1be78)
+++ src/molecule_graph.cpp	(revision 4e855e867fa84e16d3c418b917eb36c8f1b45f90)
@@ -174,10 +174,10 @@
         for (LC->n[2] = 0; LC->n[2] < LC->N[2]; LC->n[2]++) {
           const LinkedCell::LinkedNodes *List = LC->GetCurrentCell();
-//          Log() << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << " containing " << List->size() << " points." << endl;
+          Log() << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << " containing " << List->size() << " points." << endl;
           if (List != NULL) {
             for (LinkedCell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
               Walker = dynamic_cast<atom*>(*Runner);
               ASSERT(Walker,"Tesselpoint that was not an atom retrieved from LinkedNode");
-              //Log() << Verbose(0) << "Current Atom is " << *Walker << "." << endl;
+              Log() << Verbose(0) << "Current Atom is " << *Walker << "." << endl;
               // 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]++)
@@ -185,24 +185,24 @@
                   for (n[2] = -1; n[2] <= 1; n[2]++) {
                     const LinkedCell::LinkedNodes *OtherList = LC->GetRelativeToCurrentCell(n);
-//                    Log() << Verbose(2) << "Current relative cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << " containing " << List->size() << " points." << endl;
                     if (OtherList != NULL) {
+                      Log() << Verbose(2) << "Current relative cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << " containing " << List->size() << " points." << endl;
                       for (LinkedCell::LinkedNodes::const_iterator OtherRunner = OtherList->begin(); OtherRunner != OtherList->end(); OtherRunner++) {
                         if ((*OtherRunner)->nr > Walker->nr) {
                           OtherWalker = dynamic_cast<atom*>(*OtherRunner);
                           ASSERT(OtherWalker,"TesselPoint that was not an atom retrieved from LinkedNode");
-                          //Log() << Verbose(1) << "Checking distance " << OtherWalker->x.PeriodicDistanceSquared(&(Walker->x), cell_size) << " against typical bond length of " << bonddistance*bonddistance << "." << endl;
                           (BG->*minmaxdistance)(Walker, OtherWalker, MinDistance, MaxDistance, IsAngstroem);
                           const double distance = domain.periodicDistanceSquared(OtherWalker->getPosition(),Walker->getPosition());
+                          Log() << Verbose(1) << "Checking distance " << distance << " against typical bond length of " << bonddistance*bonddistance << "." << endl;
                           const bool status = (distance <= MaxDistance * MaxDistance) && (distance >= MinDistance * MinDistance);
-//                          Log() << Verbose(1) << "MinDistance is " << MinDistance << " and MaxDistance is " << MaxDistance << "." << endl;
+                          Log() << Verbose(1) << "MinDistance is " << MinDistance << " and MaxDistance is " << MaxDistance << "." << endl;
                           if (OtherWalker->father->nr > Walker->father->nr) {
                             if (status) { // create bond if distance is smaller
-//                              Log() << Verbose(1) << "Adding Bond between " << *Walker << " and " << *OtherWalker << " in distance " << sqrt(distance) << "." << endl;
+                              Log() << Verbose(1) << "Adding Bond between " << *Walker << " and " << *OtherWalker << " in distance " << sqrt(distance) << "." << endl;
                               AddBond(Walker->father, OtherWalker->father, 1); // also increases molecule::BondCount
                             } else {
-//                              Log() << Verbose(1) << "Not Adding: distance too great." << endl;
+                              Log() << Verbose(1) << "Not Adding: distance too great." << endl;
                             }
                           } else {
-//                            Log() << Verbose(1) << "Not Adding: Wrong order of labels." << endl;
+                            Log() << Verbose(1) << "Not Adding: Wrong order of labels." << endl;
                           }
                         }
Index: src/parser.cpp
===================================================================
--- src/parser.cpp	(revision fd378853756ce2d629fc774e7c82d6739ae1be78)
+++ src/parser.cpp	(revision 4e855e867fa84e16d3c418b917eb36c8f1b45f90)
@@ -22,4 +22,5 @@
 
 #include <cstring>
+#include <iomanip>
 
 #include "Helpers/helpers.hpp"
@@ -158,5 +159,5 @@
  *    -# allocate matrix
  *    -# loop over found column and row counts and parse in each entry
- * \param *name directory with files
+ * \param &input input stream
  * \param skiplines number of inital lines to skip
  * \param skiplines number of inital columns to skip
@@ -164,15 +165,13 @@
  * \return parsing successful
  */
-bool MatrixContainer::ParseMatrix(const char *name, int skiplines, int skipcolumns, int MatrixNr)
-{
-  ifstream input;
+bool MatrixContainer::ParseMatrix(std::istream &input, int skiplines, int skipcolumns, int MatrixNr)
+{
   stringstream line;
   string token;
   char filename[1023];
 
-  input.open(name, ios::in);
   //Log() << Verbose(1) << "Opening " << name << " ... "  << input << endl;
-  if (input == NULL) {
-    DoeLog(1) && (eLog()<< Verbose(1) << endl << "MatrixContainer::ParseMatrix: Unable to open " << name << ", is the directory correct?" << endl);
+  if (input.bad()) {
+    DoeLog(1) && (eLog()<< Verbose(1) << endl << "MatrixContainer::ParseMatrix: Unable to parse istream." << endl);
     //performCriticalExit();
     return false;
@@ -190,5 +189,5 @@
   for(int k=skipcolumns;k--;)
     line >> Header[MatrixNr];
-  //Log() << Verbose(0) << line.str() << endl;
+  Log() << Verbose(0) << line.str() << endl;
   ColumnCounter[MatrixNr]=0;
   while ( getline(line,token, '\t') ) {
@@ -196,8 +195,8 @@
       ColumnCounter[MatrixNr]++;
   }
-  //Log() << Verbose(0) << line.str() << endl;
-  //Log() << Verbose(1) << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << "." << endl;
+  Log() << Verbose(0) << line.str() << endl;
+  Log() << Verbose(1) << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << "." << endl;
   if (ColumnCounter[MatrixNr] == 0) {
-    DoeLog(0) && (eLog()<< Verbose(0) << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << " from file " << name << ", this is probably an error!" << endl);
+    DoeLog(0) && (eLog()<< Verbose(0) << "ColumnCounter[" << MatrixNr << "]: " << ColumnCounter[MatrixNr] << " from ostream." << endl);
     performCriticalExit();
   }
@@ -207,5 +206,5 @@
   while (!input.eof()) {
     input.getline(filename, 1023);
-    //Log() << Verbose(0) << "Comparing: " << strncmp(filename,"MeanForce",9) << endl;
+    Log() << Verbose(0) << "Comparing: " << strncmp(filename,"MeanForce",9) << endl;
     RowCounter[MatrixNr]++; // then line was not last MeanForce
     if (strncmp(filename,"MeanForce", 9) == 0) {
@@ -213,7 +212,7 @@
     }
   }
-  //Log() << Verbose(1) << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << "." << endl;
+  Log() << Verbose(1) << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from input stream." << endl;
   if (RowCounter[MatrixNr] == 0) {
-    DoeLog(0) && (eLog()<< Verbose(0) << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from file " << name << ", this is probably an error!" << endl);
+    DoeLog(0) && (eLog()<< Verbose(0) << "RowCounter[" << MatrixNr << "]: " << RowCounter[MatrixNr] << " from input stream." << endl);
     performCriticalExit();
   }
@@ -233,4 +232,5 @@
       input.getline(Header[MatrixNr], 1023);    // skip header
     line.str(Header[MatrixNr]);
+    Log() << Verbose(0) << "Header: " << line.str() << endl;
     for(int k=skipcolumns;k--;)  // skip columns in header too
       line >> filename;
@@ -245,10 +245,10 @@
       input.getline(filename, 1023);
       stringstream lines(filename);
-      //Log() << Verbose(2) << "Matrix at level " << j << ":";// << filename << endl;
+      Log() << Verbose(2) << "Matrix at level " << j << ":";// << filename << endl;
       for(int k=skipcolumns;k--;)
         lines >> filename;
       for(int k=0;(k<ColumnCounter[MatrixNr]) && (!lines.eof());k++) {
         lines >> Matrix[MatrixNr][j][k];
-        //Log() << Verbose(1) << " " << setprecision(2) << Matrix[MatrixNr][j][k] << endl;
+        Log() << Verbose(1) << " " << setprecision(2) << Matrix[MatrixNr][j][k] << endl;
       }
       if (Matrix[MatrixNr][ RowCounter[MatrixNr] ] != NULL)
@@ -261,5 +261,4 @@
     DoeLog(1) && (eLog()<< Verbose(1) << "Matrix nr. " << MatrixNr << " has column and row count of (" << ColumnCounter[MatrixNr] << "," << RowCounter[MatrixNr] << "), could not allocate nor parse!" << endl); 
   }
-  input.close();
   return true;
 };
@@ -327,6 +326,10 @@
     file.str(" ");
     file << name << FRAGMENTPREFIX << FragmentNumber << prefix << suffix;
-    if (!ParseMatrix(file.str().c_str(), skiplines, skipcolumns, i))
+    std::ifstream input(file.str().c_str());
+    if (!ParseMatrix(input, skiplines, skipcolumns, i)) {
+      input.close();
       return false;
+    }
+    input.close();
     delete[](FragmentNumber);
   }
@@ -628,5 +631,7 @@
     strncat(filename, prefix, 1023-strlen(filename));
     strncat(filename, suffix.c_str(), 1023-strlen(filename));
-    ParseMatrix(filename, skiplines, skipcolumns, MatrixCounter);
+    std::ifstream input(filename);
+    ParseMatrix(input, skiplines, skipcolumns, MatrixCounter);
+    input.close();
   }
   return status;
@@ -764,5 +769,7 @@
     strncat(filename, prefix, 1023-strlen(filename));
     strncat(filename, suffix.c_str(), 1023-strlen(filename));
-    ParseMatrix(filename, skiplines, skipcolumns, MatrixCounter);
+    std::ifstream input(filename);
+    ParseMatrix(input, skiplines, skipcolumns, MatrixCounter);
+    input.close();
   }
  
@@ -992,5 +999,7 @@
     strncat(filename, prefix, 1023-strlen(filename));
     strncat(filename, suffix.c_str(), 1023-strlen(filename));
-    ParseMatrix(filename, skiplines, skipcolumns, MatrixCounter);
+    std::ifstream input(filename);
+    ParseMatrix(input, skiplines, skipcolumns, MatrixCounter);
+    input.close();
   }
 
Index: src/parser.hpp
===================================================================
--- src/parser.hpp	(revision fd378853756ce2d629fc774e7c82d6739ae1be78)
+++ src/parser.hpp	(revision 4e855e867fa84e16d3c418b917eb36c8f1b45f90)
@@ -86,5 +86,5 @@
   
   bool InitialiseIndices(class MatrixContainer *Matrix = NULL);
-  bool ParseMatrix(const char *name, int skiplines, int skipcolumns, int MatrixNr);
+  bool ParseMatrix(std::istream &input, int skiplines, int skipcolumns, int MatrixNr);
   virtual bool ParseFragmentMatrix(const char *name, const char *prefix, string suffix, int skiplines, int skipcolumns);
   bool AllocateMatrix(char **GivenHeader, int MCounter, int *RCounter, int *CCounter);
Index: src/unittests/AnalysisBondsUnitTest.cpp
===================================================================
--- src/unittests/AnalysisBondsUnitTest.cpp	(revision fd378853756ce2d629fc774e7c82d6739ae1be78)
+++ src/unittests/AnalysisBondsUnitTest.cpp	(revision 4e855e867fa84e16d3c418b917eb36c8f1b45f90)
@@ -92,8 +92,6 @@
   CPPUNIT_ASSERT_EQUAL( TestMolecule->getAtomCount(), 5 );
 
-  // create a small file with table
-  filename = new string("test.dat");
-  CPPUNIT_ASSERT(filename != NULL && "could not create string");
-  ofstream test(filename->c_str());
+  // create stream with table
+  std::stringstream test;
   test << ".\tH\tHe\tLi\tBe\tB\tC\n";
   test << "H\t1.\t1.\t1.\t1.\t1.\t1.2\n";
@@ -103,9 +101,8 @@
   test << "B\t1.\t1.\t1.\t1.\t1.\t1.\n";
   test << "C\t1.2\t1.\t1.\t1.\t1.\t1.5\n";
-  test.close();
   BG = new BondGraph(true);
   CPPUNIT_ASSERT(BG != NULL && "could not create BondGraph");
 
-  CPPUNIT_ASSERT_EQUAL( true , BG->LoadBondLengthTable(*filename) );
+  CPPUNIT_ASSERT_EQUAL( true , BG->LoadBondLengthTable(test) );
   CPPUNIT_ASSERT_EQUAL( 1., BG->GetBondLength(0,0) );
   CPPUNIT_ASSERT_EQUAL( 1.2, BG->GetBondLength(0,5) );
@@ -119,6 +116,4 @@
 {
   // remove the file
-  remove(filename->c_str());
-  delete(filename);
   delete(BG);
 
Index: src/unittests/BondGraphUnitTest.cpp
===================================================================
--- src/unittests/BondGraphUnitTest.cpp	(revision fd378853756ce2d629fc774e7c82d6739ae1be78)
+++ src/unittests/BondGraphUnitTest.cpp	(revision 4e855e867fa84e16d3c418b917eb36c8f1b45f90)
@@ -92,10 +92,5 @@
   CPPUNIT_ASSERT_EQUAL( TestMolecule->getAtomCount(), 4 );
 
-  // create a small file with table
-  dummyname = new string("dummy.dat");
-  CPPUNIT_ASSERT(dummyname != NULL && "could not create string");
-  filename = new string("test.dat");
-  CPPUNIT_ASSERT(filename != NULL && "could not create string");
-  ofstream test(filename->c_str());
+  // create stream with table
   test << ".\tH\tHe\tLi\tBe\tB\tC\n";
   test << "H\t1.\t1.\t1.\t1.\t1.\t1.2\n";
@@ -105,5 +100,7 @@
   test << "B\t1.\t1.\t1.\t1.\t1.\t1.\n";
   test << "C\t1.2\t1.\t1.\t1.\t1.\t1.5\n";
-  test.close();
+  // created bad stream (i.e. non-present file)
+  dummy.setstate(ios::eofbit);
+  CPPUNIT_ASSERT(dummy.eof());
   BG = new BondGraph(true);
   CPPUNIT_ASSERT(BG != NULL && "could not create BondGraph");
@@ -114,7 +111,4 @@
 {
   // remove the file
-  remove(filename->c_str());
-  delete(filename);
-  delete(dummyname);
   delete(BG);
 
@@ -139,5 +133,5 @@
 void BondGraphTest::LoadTableTest()
 {
-  CPPUNIT_ASSERT_EQUAL( true , BG->LoadBondLengthTable(*filename) );
+  CPPUNIT_ASSERT_EQUAL( true , BG->LoadBondLengthTable(test) );
   CPPUNIT_ASSERT_EQUAL( 1., BG->GetBondLength(0,0) );
   CPPUNIT_ASSERT_EQUAL( 1.2, BG->GetBondLength(0,5) );
@@ -152,5 +146,5 @@
   molecule::iterator Runner = TestMolecule->begin();
   Runner++;
-  CPPUNIT_ASSERT_EQUAL( true , BG->LoadBondLengthTable(*filename) );
+  CPPUNIT_ASSERT_EQUAL( true , BG->LoadBondLengthTable(test) );
   CPPUNIT_ASSERT_EQUAL( true , BG->ConstructBondGraph(TestMolecule) );
   CPPUNIT_ASSERT_EQUAL( true , (*Walker)->IsBondedTo((*Runner)) );
@@ -165,5 +159,5 @@
   //atom *Runner = TestMolecule->end->previous;
   //CPPUNIT_ASSERT( TestMolecule->end != Walker );
-  CPPUNIT_ASSERT_EQUAL( false , BG->LoadBondLengthTable(*dummyname) );
+  CPPUNIT_ASSERT_EQUAL( false , BG->LoadBondLengthTable(dummy) );
   CPPUNIT_ASSERT_EQUAL( true , BG->ConstructBondGraph(TestMolecule) );
 
Index: src/unittests/BondGraphUnitTest.hpp
===================================================================
--- src/unittests/BondGraphUnitTest.hpp	(revision fd378853756ce2d629fc774e7c82d6739ae1be78)
+++ src/unittests/BondGraphUnitTest.hpp	(revision 4e855e867fa84e16d3c418b917eb36c8f1b45f90)
@@ -49,6 +49,6 @@
 
       BondGraph *BG;
-      string *filename;
-      string *dummyname;
+      std::stringstream dummy;
+      std::stringstream test;
 };
 
Index: tests/regression/Graph/1/pre/test2.conf
===================================================================
--- tests/regression/Graph/1/pre/test2.conf	(revision 4e855e867fa84e16d3c418b917eb36c8f1b45f90)
+++ tests/regression/Graph/1/pre/test2.conf	(revision 4e855e867fa84e16d3c418b917eb36c8f1b45f90)
@@ -0,0 +1,73 @@
+# ParallelCarParinello - main configuration file - created with molecuilder
+
+mainname	pcp	# programm name (for runtime files)
+defaultpath	not specified	# where to put files during runtime
+pseudopotpath	not specified	# where to find pseudopotentials
+
+ProcPEGamma	8	# for parallel computing: share constants
+ProcPEPsi	1	# for parallel computing: share wave functions
+DoOutVis	0	# Output data for OpenDX
+DoOutMes	1	# Output data for measurements
+DoOutOrbitals	0	# Output all Orbitals
+DoOutCurr	0	# Ouput current density for OpenDx
+DoOutNICS	0	# Output Nucleus independent current shieldings
+DoPerturbation	0	# Do perturbation calculate and determine susceptibility and shielding
+DoFullCurrent	0	# Do full perturbation
+DoConstrainedMD	0	# Do perform a constrained (>0, relating to current MD step) instead of unconstrained (0) MD
+Thermostat	Berendsen	2.5	# Which Thermostat and its parameters to use in MD case.
+CommonWannier	0	# Put virtual centers at indivual orbits, all common, merged by variance, to grid point, to cell center
+SawtoothStart	0.01	# Absolute value for smooth transition at cell border 
+VectorPlane	0	# Cut plane axis (x, y or z: 0,1,2) for two-dim current vector plot
+VectorCut	0	# Cut plane axis value
+AddGramSch	1	# Additional GramSchmidtOrtogonalization to be safe
+Seed		1	# initial value for random seed for Psi coefficients
+
+MaxOuterStep	0	# number of MolecularDynamics/Structure optimization steps
+Deltat	0.01	# time per MD step
+OutVisStep	10	# Output visual data every ...th step
+OutSrcStep	5	# Output "restart" data every ..th step
+TargetTemp	0.000950045	# Target temperature
+MaxPsiStep	0	# number of Minimisation steps per state (0 - default)
+EpsWannier	1e-07	# tolerance value for spread minimisation of orbitals
+
+# Values specifying when to stop
+MaxMinStep	100	# Maximum number of steps
+RelEpsTotalE	1e-07	# relative change in total energy
+RelEpsKineticE	1e-05	# relative change in kinetic energy
+MaxMinStopStep	0	# check every ..th steps
+MaxMinGapStopStep	0	# check every ..th steps
+
+# Values specifying when to stop for INIT, otherwise same as above
+MaxInitMinStep	100	# Maximum number of steps
+InitRelEpsTotalE	1e-05	# relative change in total energy
+InitRelEpsKineticE	0.0001	# relative change in kinetic energy
+InitMaxMinStopStep	0	# check every ..th steps
+InitMaxMinGapStopStep	0	# check every ..th steps
+
+BoxLength			# (Length of a unit cell)
+20	
+0	20	
+0	0	20	
+
+ECut		128	# energy cutoff for discretization in Hartrees
+MaxLevel	5	# number of different levels in the code, >=2
+Level0Factor	2	# factor by which node number increases from S to 0 level
+RiemannTensor	0	# (Use metric)
+PsiType		0	# 0 - doubly occupied, 1 - SpinUp,SpinDown
+MaxPsiDouble	0	# here: specifying both maximum number of SpinUp- and -Down-states
+PsiMaxNoUp	0	# here: specifying maximum number of SpinUp-states
+PsiMaxNoDown	0	# here: specifying maximum number of SpinDown-states
+AddPsis		0	# Additional unoccupied Psis for bandgap determination
+
+RCut		20	# R-cut for the ewald summation
+StructOpt	0	# Do structure optimization beforehand
+IsAngstroem	1	# 0 - Bohr, 1 - Angstroem
+RelativeCoord	0	# whether ion coordinates are relative (1) or absolute (0)
+MaxTypes	1	# maximum number of different ion types
+
+# Ion type data (PP = PseudoPotential, Z = atomic number)
+#Ion_TypeNr.	Amount	Z	RGauss	L_Max(PP)L_Loc(PP)IonMass	# chemical name, symbol
+Ion_Type1	2	1	1.0	3	3	1.00800000000	Hydrogen	H
+#Ion_TypeNr._Nr.R[0]    R[1]    R[2]    MoveType (0 MoveIon, 1 FixedIon)
+Ion_Type1_1	5.000000000	5.000000000	5.000000000	0 # molecule nr 3
+Ion_Type1_2	6.500000000	5.000000000	5.000000000	0 # molecule nr 4
Index: tests/regression/testsuite-graph.at
===================================================================
--- tests/regression/testsuite-graph.at	(revision fd378853756ce2d629fc774e7c82d6739ae1be78)
+++ tests/regression/testsuite-graph.at	(revision 4e855e867fa84e16d3c418b917eb36c8f1b45f90)
@@ -6,4 +6,23 @@
 AT_CHECK([../../molecuilder -i test.conf  -v 3 -D 2.], 0, [stdout], [stderr])
 AT_CHECK([fgrep -c "No rings were detected in the molecular structure." stdout], 0, [1
+], [ignore])
+AT_CLEANUP
+
+AT_SETUP([Graph - DFS analysis with bond length table])
+AT_KEYWORDS([graph])
+AT_CHECK([/bin/cp -f ${abs_top_srcdir}/${AUTOTEST_PATH}/Graph/1/pre/test2.conf .], 0)
+AT_DATA([table1.dat],[.	H
+H	1.
+])
+AT_DATA([table2.dat],[.	H
+H	1.5
+])
+AT_CHECK([../../molecuilder -i test2.conf -v 4 -g table1.dat -D 2.], 0, [stdout], [stderr])
+AT_CHECK([fgrep -c "Ion_Type1_1" stdout], 0, [2
+], [ignore])
+AT_CHECK([../../molecuilder -i test2.conf -v 4 -g table2.dat -D 3.], 0, [stdout], [stderr])
+AT_CHECK([fgrep -c "Ion_Type1_1" stdout], 0, [1
+], [ignore])
+AT_CHECK([fgrep -c "Ion_Type1_2" stdout], 0, [1
 ], [ignore])
 AT_CLEANUP
