Index: src/atom.cpp
===================================================================
--- src/atom.cpp	(revision d8b94a86944c4679c88bc250bff510f242d63309)
+++ src/atom.cpp	(revision 1907a789609146bbdebfeff09de44cb4ba21b808)
@@ -1,15 +1,15 @@
 /** \file atom.cpp
- * 
+ *
  * Function implementations for the class atom.
- * 
+ *
  */
 
 #include "molecules.hpp"
- 
+
 /************************************* Functions for class atom *************************************/
 
 /** Constructor of class atom.
  */
-atom::atom() 
+atom::atom()
 {
 	Name = NULL;
@@ -33,5 +33,5 @@
 /** Destructor of class atom.
  */
-atom::~atom() 
+atom::~atom()
 {
 	Free((void **)&Name, "atom::~atom: *Name");
@@ -58,6 +58,7 @@
  * \param AtomNo cardinal number among these atoms of the same element
  * \param *out stream to output to
+ * \param *comment commentary after '#' sign
  */
-bool atom::Output(int ElementNo, int AtomNo, ofstream *out) const
+bool atom::Output(int ElementNo, int AtomNo, ofstream *out, const char *comment) const
 {
 	if (out != NULL) {
@@ -67,5 +68,8 @@
 		if (v.Norm() > MYEPSILON)
 			*out << "\t" << scientific << setprecision(6) << v.x[0] << "\t" << v.x[1] << "\t" << v.x[2] << "\t";
-		*out << " # Number in molecule " << nr << endl;
+		if (comment != NULL)
+		  *out << " # " << comment << endl;
+		else
+		  *out << " # molecule nr " << nr << endl;
 		return true;
 	} else
@@ -85,5 +89,5 @@
 };
 
-ostream & operator << (ostream &ost, atom &a) 
+ostream & operator << (ostream &ost, atom &a)
 {
 	ost << "[" << a.Name << "|" << &a << "]";
@@ -94,5 +98,5 @@
  * \param ptr atom to compare index against
  * \return true - this one's is smaller, false - not
- */ 
+ */
 bool atom::Compare(atom &ptr)
 {
@@ -103,5 +107,5 @@
 };
 
-bool operator < (atom &a, atom &b) 
+bool operator < (atom &a, atom &b)
 {
 	return a.Compare(b);
Index: src/builder.cpp
===================================================================
--- src/builder.cpp	(revision d8b94a86944c4679c88bc250bff510f242d63309)
+++ src/builder.cpp	(revision 1907a789609146bbdebfeff09de44cb4ba21b808)
@@ -55,9 +55,9 @@
 #include "molecules.hpp"
 
-/********************************************** Submenu routine **************************************/
+/********************************************* Subsubmenu routine ************************************/
 
 /** Submenu for adding atoms to the molecule.
  * \param *periode periodentafel
- * \param *mol the molecule to add to
+ * \param *molecule molecules with atoms
  */
 static void AddAtoms(periodentafel *periode, molecule *mol)
@@ -83,172 +83,175 @@
 
 	switch (choice) {
+    default:
+      cout << Verbose(0) << "Not a valid choice." << endl;
+      break;
 			case 'a': // absolute coordinates of atom
-				cout << Verbose(0) << "Enter absolute coordinates." << endl;
-				first = new atom;
-				first->x.AskPosition(mol->cell_size, false);
-				first->type = periode->AskElement();	// give type
-				mol->AddAtom(first);	// add to molecule
+        cout << Verbose(0) << "Enter absolute coordinates." << endl;
+        first = new atom;
+        first->x.AskPosition(mol->cell_size, false);
+        first->type = periode->AskElement();	// give type
+        mol->AddAtom(first);	// add to molecule
 				break;
 
 			case 'b': // relative coordinates of atom wrt to reference point
-				first = new atom;
-				valid = true;
-				do {
-					if (!valid) cout << Verbose(0) << "Resulting position out of cell." << endl;
-					cout << Verbose(0) << "Enter reference coordinates." << endl;
-					x.AskPosition(mol->cell_size, true);
-					cout << Verbose(0) << "Enter relative coordinates." << endl;
-					first->x.AskPosition(mol->cell_size, false);
-					first->x.AddVector((const Vector *)&x);
-					cout << Verbose(0) << "\n";
-				} while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
-				first->type = periode->AskElement();	// give type
-				mol->AddAtom(first);	// add to molecule
+        first = new atom;
+        valid = true;
+        do {
+          if (!valid) cout << Verbose(0) << "Resulting position out of cell." << endl;
+          cout << Verbose(0) << "Enter reference coordinates." << endl;
+          x.AskPosition(mol->cell_size, true);
+          cout << Verbose(0) << "Enter relative coordinates." << endl;
+          first->x.AskPosition(mol->cell_size, false);
+          first->x.AddVector((const Vector *)&x);
+          cout << Verbose(0) << "\n";
+        } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
+        first->type = periode->AskElement();	// give type
+        mol->AddAtom(first);	// add to molecule
 				break;
 
 			case 'c': // relative coordinates of atom wrt to already placed atom
-				first = new atom;
-				valid = true;
-				do {
-					if (!valid) cout << Verbose(0) << "Resulting position out of cell." << endl;
-					second = mol->AskAtom("Enter atom number: ");
-					cout << Verbose(0) << "Enter relative coordinates." << endl;
-					first->x.AskPosition(mol->cell_size, false);
-					for (int i=NDIM;i--;) {
-						first->x.x[i] += second->x.x[i];
-					}
-				} while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
-				first->type = periode->AskElement();	// give type
-				mol->AddAtom(first);	// add to molecule
+        first = new atom;
+        valid = true;
+        do {
+          if (!valid) cout << Verbose(0) << "Resulting position out of cell." << endl;
+          second = mol->AskAtom("Enter atom number: ");
+          cout << Verbose(0) << "Enter relative coordinates." << endl;
+          first->x.AskPosition(mol->cell_size, false);
+          for (int i=NDIM;i--;) {
+            first->x.x[i] += second->x.x[i];
+          }
+        } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
+        first->type = periode->AskElement();	// give type
+        mol->AddAtom(first);	// add to molecule
+        break;
+
+    case 'd': // two atoms, two angles and a distance
+        first = new atom;
+        valid = true;
+        do {
+          if (!valid) {
+            cout << Verbose(0) << "Resulting coordinates out of cell - ";
+            first->x.Output((ofstream *)&cout);
+            cout << Verbose(0) << endl;
+          }
+          cout << Verbose(0) << "First, we need two atoms, the first atom is the central, while the second is the outer one." << endl;
+          second = mol->AskAtom("Enter central atom: ");
+          third = mol->AskAtom("Enter second atom (specifying the axis for first angle): ");
+          fourth = mol->AskAtom("Enter third atom (specifying a plane for second angle): ");
+          a = ask_value("Enter distance between central (first) and new atom: ");
+          b = ask_value("Enter angle between new, first and second atom (degrees): ");
+          b *= M_PI/180.;
+          bound(&b, 0., 2.*M_PI);
+          c = ask_value("Enter second angle between new and normal vector of plane defined by first, second and third atom (degrees): ");
+          c *= M_PI/180.;
+          bound(&c, -M_PI, M_PI);
+          cout << Verbose(0) << "radius: " << a << "\t phi: " << b*180./M_PI << "\t theta: " << c*180./M_PI << endl;
+/*
+          second->Output(1,1,(ofstream *)&cout);
+          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);
+          x.SubtractVector(&third->x);
+          x.Copyvector(&fourth->x);
+          x.SubtractVector(&third->x);
+
+          if (!z.SolveSystem(&x,&y,&n, b, c, a)) {
+            cout << Verbose(0) << "Failure solving self-dependent linear system!" << endl;
+            continue;
+          }
+          cout << Verbose(0) << "resulting relative coordinates: ";
+          z.Output((ofstream *)&cout);
+          cout << Verbose(0) << endl;
+          */
+          // calc axis vector
+          x.CopyVector(&second->x);
+          x.SubtractVector(&third->x);
+          x.Normalize();
+          cout << "x: ",
+          x.Output((ofstream *)&cout);
+          cout << endl;
+          z.MakeNormalVector(&second->x,&third->x,&fourth->x);
+          cout << "z: ",
+          z.Output((ofstream *)&cout);
+          cout << endl;
+          y.MakeNormalVector(&x,&z);
+          cout << "y: ",
+          y.Output((ofstream *)&cout);
+          cout << endl;
+
+          // rotate vector around first angle
+          first->x.CopyVector(&x);
+          first->x.RotateVector(&z,b - M_PI);
+          cout << "Rotated vector: ",
+          first->x.Output((ofstream *)&cout);
+          cout << endl;
+          // remove the projection onto the rotation plane of the second angle
+          n.CopyVector(&y);
+          n.Scale(first->x.Projection(&y));
+          cout << "N1: ",
+          n.Output((ofstream *)&cout);
+          cout << endl;
+          first->x.SubtractVector(&n);
+          cout << "Subtracted vector: ",
+          first->x.Output((ofstream *)&cout);
+          cout << endl;
+          n.CopyVector(&z);
+          n.Scale(first->x.Projection(&z));
+          cout << "N2: ",
+          n.Output((ofstream *)&cout);
+          cout << endl;
+          first->x.SubtractVector(&n);
+          cout << "2nd subtracted vector: ",
+          first->x.Output((ofstream *)&cout);
+          cout << endl;
+
+          // rotate another vector around second angle
+          n.CopyVector(&y);
+          n.RotateVector(&x,c - M_PI);
+          cout << "2nd Rotated vector: ",
+          n.Output((ofstream *)&cout);
+          cout << endl;
+
+          // add the two linear independent vectors
+          first->x.AddVector(&n);
+          first->x.Normalize();
+          first->x.Scale(a);
+          first->x.AddVector(&second->x);
+
+          cout << Verbose(0) << "resulting coordinates: ";
+          first->x.Output((ofstream *)&cout);
+          cout << Verbose(0) << endl;
+        } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
+        first->type = periode->AskElement();	// give type
+        mol->AddAtom(first);	// add to molecule
 				break;
 
-			case 'd': // two atoms, two angles and a distance
-				first = new atom;
-				valid = true;
-				do {
-					if (!valid) {
-						cout << Verbose(0) << "Resulting coordinates out of cell - ";
-						first->x.Output((ofstream *)&cout);
-						cout << Verbose(0) << endl;
-					}
-					cout << Verbose(0) << "First, we need two atoms, the first atom is the central, while the second is the outer one." << endl;
-					second = mol->AskAtom("Enter central atom: ");
-					third = mol->AskAtom("Enter second atom (specifying the axis for first angle): ");
-					fourth = mol->AskAtom("Enter third atom (specifying a plane for second angle): ");
-					a = ask_value("Enter distance between central (first) and new atom: ");
-					b = ask_value("Enter angle between new, first and second atom (degrees): ");
-					b *= M_PI/180.;
-					bound(&b, 0., 2.*M_PI);
-					c = ask_value("Enter second angle between new and normal vector of plane defined by first, second and third atom (degrees): ");
-					c *= M_PI/180.;
-					bound(&c, -M_PI, M_PI);
-					cout << Verbose(0) << "radius: " << a << "\t phi: " << b*180./M_PI << "\t theta: " << c*180./M_PI << endl;
-/*
-					second->Output(1,1,(ofstream *)&cout);
-					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);
-					x.SubtractVector(&third->x);
-					x.Copyvector(&fourth->x);
-					x.SubtractVector(&third->x);
-
-					if (!z.SolveSystem(&x,&y,&n, b, c, a)) {
-						cout << Verbose(0) << "Failure solving self-dependent linear system!" << endl;
-						continue;
-					}
-					cout << Verbose(0) << "resulting relative coordinates: ";
-					z.Output((ofstream *)&cout);
-					cout << Verbose(0) << endl;
-					*/
-					// calc axis vector
-					x.CopyVector(&second->x);
-					x.SubtractVector(&third->x);
-					x.Normalize();
-					cout << "x: ",
-					x.Output((ofstream *)&cout);
-					cout << endl;
-					z.MakeNormalVector(&second->x,&third->x,&fourth->x);
-					cout << "z: ",
-					z.Output((ofstream *)&cout);
-					cout << endl;
-					y.MakeNormalVector(&x,&z);
-					cout << "y: ",
-					y.Output((ofstream *)&cout);
-					cout << endl;
-
-					// rotate vector around first angle
-					first->x.CopyVector(&x);
-					first->x.RotateVector(&z,b - M_PI);
-					cout << "Rotated vector: ",
-					first->x.Output((ofstream *)&cout);
-					cout << endl;
-					// remove the projection onto the rotation plane of the second angle
-					n.CopyVector(&y);
-					n.Scale(first->x.Projection(&y));
-					cout << "N1: ",
-					n.Output((ofstream *)&cout);
-					cout << endl;
-					first->x.SubtractVector(&n);
-					cout << "Subtracted vector: ",
-					first->x.Output((ofstream *)&cout);
-					cout << endl;
-					n.CopyVector(&z);
-					n.Scale(first->x.Projection(&z));
-					cout << "N2: ",
-					n.Output((ofstream *)&cout);
-					cout << endl;
-					first->x.SubtractVector(&n);
-					cout << "2nd subtracted vector: ",
-					first->x.Output((ofstream *)&cout);
-					cout << endl;
-
-					// rotate another vector around second angle
-					n.CopyVector(&y);
-					n.RotateVector(&x,c - M_PI);
-					cout << "2nd Rotated vector: ",
-					n.Output((ofstream *)&cout);
-					cout << endl;
-
-					// add the two linear independent vectors
-					first->x.AddVector(&n);
-					first->x.Normalize();
-					first->x.Scale(a);
-					first->x.AddVector(&second->x);
-
-					cout << Verbose(0) << "resulting coordinates: ";
-					first->x.Output((ofstream *)&cout);
-					cout << Verbose(0) << endl;
-				} while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
-				first->type = periode->AskElement();	// give type
-				mol->AddAtom(first);	// add to molecule
-				break;
-
 			case 'e': // least square distance position to a set of atoms
-				first = new atom;
-				atoms = new (Vector*[128]);
-				valid = true;
-				for(int i=128;i--;)
-					atoms[i] = NULL;
-				int i=0, j=0;
-				cout << Verbose(0) << "Now we need at least three molecules.\n";
-				do {
-					cout << Verbose(0) << "Enter " << i+1 << "th atom: ";
-					cin >> j;
-					if (j != -1) {
-						second = mol->FindAtom(j);
-						atoms[i++] = &(second->x);
-					}
-				} while ((j != -1) && (i<128));
-				if (i >= 2) {
-					first->x.LSQdistance(atoms, i);
-
-					first->x.Output((ofstream *)&cout);
-					first->type = periode->AskElement();	// give type
-					mol->AddAtom(first);	// add to molecule
-				} else {
-					delete first;
-					cout << Verbose(0) << "Please enter at least two vectors!\n";
-				}
+        first = new atom;
+        atoms = new (Vector*[128]);
+        valid = true;
+        for(int i=128;i--;)
+          atoms[i] = NULL;
+        int i=0, j=0;
+        cout << Verbose(0) << "Now we need at least three molecules.\n";
+        do {
+          cout << Verbose(0) << "Enter " << i+1 << "th atom: ";
+          cin >> j;
+          if (j != -1) {
+            second = mol->FindAtom(j);
+            atoms[i++] = &(second->x);
+          }
+        } while ((j != -1) && (i<128));
+        if (i >= 2) {
+          first->x.LSQdistance(atoms, i);
+
+          first->x.Output((ofstream *)&cout);
+          first->type = periode->AskElement();	// give type
+          mol->AddAtom(first);	// add to molecule
+        } else {
+          delete first;
+          cout << Verbose(0) << "Please enter at least two vectors!\n";
+        }
 				break;
 	};
@@ -256,5 +259,5 @@
 
 /** Submenu for centering the atoms in the molecule.
- * \param *mol the molecule with all the atoms
+ * \param *mol molecule with all the atoms
  */
 static void CenterAtoms(molecule *mol)
@@ -314,5 +317,5 @@
 /** Submenu for aligning the atoms in the molecule.
  * \param *periode periodentafel
- * \param *mol the molecule with all the atoms
+ * \param *mol molecule with all the atoms
  */
 static void AlignAtoms(periodentafel *periode, molecule *mol)
@@ -382,5 +385,5 @@
 
 /** Submenu for mirroring the atoms in the molecule.
- * \param *mol the molecule with all the atoms
+ * \param *mol molecule with all the atoms
  */
 static void MirrorAtoms(molecule *mol)
@@ -429,5 +432,5 @@
 
 /** Submenu for removing the atoms from the molecule.
- * \param *mol the molecule with all the atoms
+ * \param *mol molecule with all the atoms
  */
 static void RemoveAtoms(molecule *mol)
@@ -487,5 +490,5 @@
 /** Submenu for measuring out the atoms in the molecule.
  * \param *periode periodentafel
- * \param *mol the molecule with all the atoms
+ * \param *mol molecule with all the atoms
  */
 static void MeasureAtoms(periodentafel *periode, molecule *mol, config *configuration)
@@ -597,5 +600,5 @@
 
 /** Submenu for measuring out the atoms in the molecule.
- * \param *mol the molecule with all the atoms
+ * \param *mol molecule with all the atoms
  * \param *configuration configuration structure for the to be written config files of all fragments
  */
@@ -617,14 +620,444 @@
 };
 
+/********************************************** Submenu routine **************************************/
+
+/** Submenu for manipulating atoms.
+ * \param *periode periodentafel
+ * \param *molecules list of molecules whose atoms are to be manipulated
+ */
+static void ManipulateAtoms(periodentafel *periode, MoleculeListClass *molecules, config *configuration)
+{
+  atom *first, *second, *third, *fourth;
+  Vector **atoms;
+  molecule *mol = NULL;
+  Vector x,y,z,n; // coordinates for absolute point in cell volume
+  double *factor; // unit factor if desired
+  double a,b,c;
+  double bond, min_bond;
+  char choice;  // menu choice char
+  bool valid;
+
+  cout << Verbose(0) << "=========MANIPULATE ATOMS======================" << endl;
+  cout << Verbose(0) << "a - add an atom" << endl;
+  cout << Verbose(0) << "r - remove an atom" << endl;
+  cout << Verbose(0) << "b - scale a bond between atoms" << endl;
+  cout << Verbose(0) << "u - change an atoms element" << endl;
+  cout << Verbose(0) << "l - measure lengths, angles, ... for an atom" << endl;
+  cout << Verbose(0) << "all else - go back" << endl;
+  cout << Verbose(0) << "===============================================" << endl;
+  if (molecules->NumberOfActiveMolecules() > 0)
+    cout << Verbose(0) << "WARNING: There is more than one molecule active! Atoms will be added to each." << endl;
+  cout << Verbose(0) << "INPUT: ";
+  cin >> choice;
+
+  switch (choice) {
+    default:
+      cout << Verbose(0) << "Not a valid choice." << endl;
+      break;
+
+    case 'a': // add atom
+      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
+        mol = *ListRunner;
+        cout << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
+        AddAtoms(periode, mol);
+      }
+      break;
+
+    case 'b': // scale a bond
+      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
+        mol = *ListRunner;
+        cout << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
+        cout << Verbose(0) << "Scaling bond length between two atoms." << endl;
+        first = mol->AskAtom("Enter first (fixed) atom: ");
+        second = mol->AskAtom("Enter second (shifting) atom: ");
+        min_bond = 0.;
+        for (int i=NDIM;i--;)
+          min_bond += (first->x.x[i]-second->x.x[i])*(first->x.x[i] - second->x.x[i]);
+        min_bond = sqrt(min_bond);
+        cout << Verbose(0) << "Current Bond length between " << first->type->name << " Atom " << first->nr << " and " << second->type->name << " Atom " << second->nr << ": " << min_bond << " a.u." << endl;
+        cout << Verbose(0) << "Enter new bond length [a.u.]: ";
+        cin >> bond;
+        for (int i=NDIM;i--;) {
+          second->x.x[i] -= (second->x.x[i]-first->x.x[i])/min_bond*(min_bond-bond);
+        }
+        //cout << Verbose(0) << "New coordinates of Atom " << second->nr << " are: ";
+        //second->Output(second->type->No, 1, (ofstream *)&cout);
+      }
+      break;
+
+    case 'c': // unit scaling of the metric
+      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
+        mol = *ListRunner;
+        cout << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
+       cout << Verbose(0) << "Angstroem -> Bohrradius: 1.8897261\t\tBohrradius -> Angstroem: 0.52917721" << endl;
+       cout << Verbose(0) << "Enter three factors: ";
+       factor = new double[NDIM];
+       cin >> factor[0];
+       cin >> factor[1];
+       cin >> factor[2];
+       valid = true;
+       mol->Scale(&factor);
+       delete[](factor);
+      }
+     break;
+
+    case 'l': // measure distances or angles
+      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
+        mol = *ListRunner;
+        cout << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
+        MeasureAtoms(periode, mol, configuration);
+      }
+      break;
+
+    case 'r': // remove atom
+      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
+        mol = *ListRunner;
+        cout << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
+        RemoveAtoms(mol);
+      }
+      break;
+
+    case 'u': // change an atom's element
+      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
+        int Z;
+        mol = *ListRunner;
+        cout << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
+        first = NULL;
+        do {
+          cout << Verbose(0) << "Change the element of which atom: ";
+          cin >> Z;
+        } while ((first = mol->FindAtom(Z)) == NULL);
+        cout << Verbose(0) << "New element by atomic number Z: ";
+        cin >> Z;
+        first->type = periode->FindElement(Z);
+        cout << Verbose(0) << "Atom " << first->nr << "'s element is " << first->type->name << "." << endl;
+      }
+      break;
+  }
+};
+
+/** Submenu for manipulating molecules.
+ * \param *periode periodentafel
+ * \param *molecules list of molecule to manipulate
+ */
+static void ManipulateMolecules(periodentafel *periode, MoleculeListClass *molecules, config *configuration)
+{
+  atom *first, *second, *third, *fourth;
+  Vector **atoms;
+  Vector x,y,z,n; // coordinates for absolute point in cell volume
+  double a,b,c;
+  int j, axis, count, faktor;
+  char choice;  // menu choice char
+  bool valid;
+  molecule *mol = NULL;
+  element **Elements;
+  Vector **vectors;
+  MoleculeLeafClass *Subgraphs = NULL;
+
+  cout << Verbose(0) << "=========MANIPULATE GLOBALLY===================" << endl;
+  cout << Verbose(0) << "c - scale by unit transformation" << endl;
+  cout << Verbose(0) << "d - duplicate molecule/periodic cell" << endl;
+  cout << Verbose(0) << "f - fragment molecule many-body bond order style" << endl;
+  cout << Verbose(0) << "g - center atoms in box" << endl;
+  cout << Verbose(0) << "i - realign molecule" << endl;
+  cout << Verbose(0) << "m - mirror all molecules" << endl;
+  cout << Verbose(0) << "o - create connection matrix" << endl;
+  cout << Verbose(0) << "t - translate molecule by vector" << endl;
+  cout << Verbose(0) << "all else - go back" << endl;
+  cout << Verbose(0) << "===============================================" << endl;
+  if (molecules->NumberOfActiveMolecules() > 0)
+    cout << Verbose(0) << "WARNING: There is more than one molecule active! Atoms will be added to each." << endl;
+  cout << Verbose(0) << "INPUT: ";
+  cin >> choice;
+
+  switch (choice) {
+    default:
+      cout << Verbose(0) << "Not a valid choice." << endl;
+      break;
+
+    case 'd': // duplicate the periodic cell along a given axis, given times
+      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
+        mol = *ListRunner;
+        cout << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
+        cout << Verbose(0) << "State the axis [(+-)123]: ";
+        cin >> axis;
+        cout << Verbose(0) << "State the factor: ";
+        cin >> faktor;
+
+        mol->CountAtoms((ofstream *)&cout); // recount atoms
+        if (mol->AtomCount != 0) {  // if there is more than none
+          count = mol->AtomCount;  // is changed becausing of adding, thus has to be stored away beforehand
+          Elements = new element *[count];
+          vectors = new Vector *[count];
+          j = 0;
+          first = mol->start;
+          while (first->next != mol->end) { // make a list of all atoms with coordinates and element
+            first = first->next;
+            Elements[j] = first->type;
+            vectors[j] = &first->x;
+            j++;
+          }
+          if (count != j)
+            cout << Verbose(0) << "ERROR: AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl;
+          x.Zero();
+          y.Zero();
+          y.x[abs(axis)-1] = mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude
+          for (int i=1;i<faktor;i++) {  // then add this list with respective translation factor times
+            x.AddVector(&y); // per factor one cell width further
+            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.AddVector(&x);     // translate the coordinates
+              first->type = Elements[k];  // insert original element
+              mol->AddAtom(first);        // and add to the molecule (which increments ElementsInMolecule, AtomCount, ...)
+            }
+          }
+          if (mol->first->next != mol->last) // if connect matrix is present already, redo it
+            mol->CreateAdjacencyList((ofstream *)&cout, mol->BondDistance, configuration->GetIsAngstroem());
+          // free memory
+          delete[](Elements);
+          delete[](vectors);
+          // correct cell size
+          if (axis < 0) { // if sign was negative, we have to translate everything
+            x.Zero();
+            x.AddVector(&y);
+            x.Scale(-(faktor-1));
+            mol->Translate(&x);
+          }
+          mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
+        }
+      }
+      break;
+
+    case 'f':
+      FragmentAtoms(mol, configuration);
+      break;
+
+    case 'g': // center the atoms
+      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
+        mol = *ListRunner;
+        cout << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
+        CenterAtoms(mol);
+      }
+      break;
+
+    case 'i': // align all atoms
+      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
+        mol = *ListRunner;
+        cout << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
+        AlignAtoms(periode, mol);
+      }
+      break;
+
+    case 'm': // mirror atoms along a given axis
+      for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
+        mol = *ListRunner;
+        cout << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
+        MirrorAtoms(mol);
+      }
+      break;
+
+    case 'o': // create the connection matrix
+      {
+        double bonddistance;
+        clock_t start,end;
+        cout << Verbose(0) << "What's the maximum bond distance: ";
+        cin >> bonddistance;
+        start = clock();
+        mol->CreateAdjacencyList((ofstream *)&cout, bonddistance, configuration->GetIsAngstroem());
+        mol->CreateListOfBondsPerAtom((ofstream *)&cout);
+        end = clock();
+        cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
+      }
+      break;
+
+    case 't': // translate all atoms
+     for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) {
+        mol = *ListRunner;
+        cout << Verbose(0) << "Current molecule is: " << mol->IndexNr << "\t" << mol->name << endl;
+        cout << Verbose(0) << "Enter translation vector." << endl;
+        x.AskPosition(mol->cell_size,0);
+        mol->Translate((const Vector *)&x);
+     }
+     break;
+  }
+  // Free all
+  if (Subgraphs != NULL) {  // free disconnected subgraph list of DFS analysis was performed
+    while (Subgraphs->next != NULL) {
+      Subgraphs = Subgraphs->next;
+      delete(Subgraphs->previous);
+    }
+    delete(Subgraphs);
+  }
+};
+
+
+/** Submenu for creating new molecules.
+ * \param *periode periodentafel
+ * \param *molecules list of molecules to add to
+ */
+static void EditMolecules(periodentafel *periode, MoleculeListClass *molecules)
+{
+  char choice;  // menu choice char
+  bool valid;
+  Vector Center;
+  int nr, count;
+  molecule *mol = NULL;
+  char *molname = NULL;
+  int length;
+  char filename[MAXSTRINGSIZE];
+
+  cout << Verbose(0) << "==========Edit MOLECULES=====================" << endl;
+  cout << Verbose(0) << "c - create new molecule" << endl;
+  cout << Verbose(0) << "l - load molecule from xyz file" << endl;
+  cout << Verbose(0) << "n - change molecule's name" << endl;
+  cout << Verbose(0) << "N - give molecules filename" << endl;
+  cout << Verbose(0) << "p - parse xyz file into molecule" << endl;
+  cout << Verbose(0) << "r - remove a molecule" << endl;
+  cout << Verbose(0) << "all else - go back" << endl;
+  cout << Verbose(0) << "===============================================" << endl;
+  cout << Verbose(0) << "INPUT: ";
+  cin >> choice;
+
+  switch (choice) {
+    default:
+      cout << Verbose(0) << "Not a valid choice." << endl;
+      break;
+    case 'c':
+      mol = new molecule(periode);
+      molecules->insert(mol);
+      break;
+
+    case 'l': // laod from XYZ file
+      cout << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl;
+      mol = new molecule(periode);
+      do {
+        cout << Verbose(0) << "Enter file name: ";
+        cin >> filename;
+      } while (!mol->AddXYZFile(filename));
+      mol->SetNameFromFilename(filename);
+      // center at set box dimensions
+      mol->CenterEdge((ofstream *)&cout, &Center);
+      mol->cell_size[0] = Center.x[0];
+      mol->cell_size[1] = 0;
+      mol->cell_size[2] = Center.x[1];
+      mol->cell_size[3] = 0;
+      mol->cell_size[4] = 0;
+      mol->cell_size[5] = Center.x[2];
+      molecules->insert(mol);
+      break;
+
+    case 'n':
+      do {
+        cout << Verbose(0) << "Enter index of molecule: ";
+        cin >> nr;
+        mol = molecules->ReturnIndex(nr);
+      } while (mol != NULL);
+      cout << Verbose(0) << "Enter name: ";
+      cin >> filename;
+      strcpy(mol->name, filename);
+      break;
+
+    case 'N':
+      do {
+        cout << Verbose(0) << "Enter index of molecule: ";
+        cin >> nr;
+        mol = molecules->ReturnIndex(nr);
+      } while (mol != NULL);
+      cout << Verbose(0) << "Enter name: ";
+      cin >> filename;
+      mol->SetNameFromFilename(filename);
+      break;
+
+    case 'p': // parse XYZ file
+      mol = NULL;
+      do {
+        cout << Verbose(0) << "Enter index of molecule: ";
+        cin >> nr;
+        mol = molecules->ReturnIndex(nr);
+      } while (mol == NULL);
+      cout << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl;
+      do {
+        cout << Verbose(0) << "Enter file name: ";
+        cin >> filename;
+      } while (!mol->AddXYZFile(filename));
+      mol->SetNameFromFilename(filename);
+      break;
+
+    case 'r':
+      cout << Verbose(0) << "Enter index of molecule: ";
+      cin >> nr;
+      count = 1;
+      MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin();
+      for(; ((ListRunner != molecules->ListOfMolecules.end()) && (count < nr)); ListRunner++);
+      mol = *ListRunner;
+      if (count == nr) {
+        molecules->ListOfMolecules.erase(ListRunner);
+        delete(mol);
+      }
+      break;
+  }
+};
+
+
+/** Submenu for merging molecules.
+ * \param *periode periodentafel
+ * \param *molecules list of molecules to add to
+ */
+static void MergeMolecules(periodentafel *periode, MoleculeListClass *molecules)
+{
+  char choice;  // menu choice char
+  bool valid;
+
+  cout << Verbose(0) << "===========MERGE MOLECULES=====================" << endl;
+  cout << Verbose(0) << "e - embedding merge of two molecules" << endl;
+  cout << Verbose(0) << "m - multi-merge of all molecules" << endl;
+  cout << Verbose(0) << "s - scatter merge of two molecules" << endl;
+  cout << Verbose(0) << "t - simple merge of two molecules" << endl;
+  cout << Verbose(0) << "all else - go back" << endl;
+  cout << Verbose(0) << "===============================================" << endl;
+  cout << Verbose(0) << "INPUT: ";
+  cin >> choice;
+
+  switch (choice) {
+    default:
+      cout << Verbose(0) << "Not a valid choice." << endl;
+      break;
+
+    case 'e':
+      break;
+
+    case 'm':
+      break;
+
+    case 's':
+      break;
+
+    case 't':
+      break;
+  }
+};
+
+
 /********************************************** Test routine **************************************/
 
 /** Is called always as option 'T' in the menu.
+ * \param *molecules list of molecules
  */
-static void testroutine(molecule *mol)
+static void testroutine(MoleculeListClass *molecules)
 {
 	// the current test routine checks the functionality of the KeySet&Graph concept:
 	// We want to have a multiindex (the KeySet) describing a unique subgraph
-	atom *Walker = mol->start;
-	int i, comp, counter=0;
+  int i, comp, counter=0;
+
+  // create a clone
+  molecule *mol = NULL;
+  if (molecules->ListOfMolecules.size() != 0) // clone
+    mol = (molecules->ListOfMolecules.front())->CopyMolecule();
+  else {
+    cerr << "I don't have anything to test on ... ";
+    return;
+  }
+  atom *Walker = mol->start;
 
 	// generate some KeySets
@@ -686,4 +1119,5 @@
 		A++;
 	}
+	delete(mol);
 };
 
@@ -692,10 +1126,19 @@
  * \param *configuration pointer to configuration structure with all the values
  * \param *periode pointer to periodentafel structure with all the elements
- * \param *mol pointer to molecule structure with all the atoms and coordinates
+ * \param *molecules list of molecules structure with all the atoms and coordinates
  */
-static void SaveConfig(char *ConfigFileName, config *configuration, periodentafel *periode, molecule *mol)
+static void SaveConfig(char *ConfigFileName, config *configuration, periodentafel *periode, MoleculeListClass *molecules)
 {
 	char filename[MAXSTRINGSIZE];
 	ofstream output;
+	molecule *mol = new molecule(periode);
+
+	// merge all molecules in MoleculeListClass into this molecule
+	int N = molecules->ListOfMolecules.size();
+	int *src = new int(N);
+	N=0;
+	for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++)
+	  src[N++] = (*ListRunner)->IndexNr;
+	molecules->SimpleMultiAdd(mol, src, N);
 
 	cout << Verbose(0) << "Storing configuration ... " << endl;
@@ -761,4 +1204,5 @@
 		cerr << "WARNING: config is found under different path then stated in config file::defaultpath!" << endl;
 	}
+	delete(mol);
 };
 
@@ -766,5 +1210,5 @@
  * \param argc argument count
  * \param **argv arguments array
- * \param *mol molecule structure
+ * \param *molecules list of molecules structure
  * \param *periode elements structure
  * \param configuration config file structure
@@ -773,5 +1217,5 @@
  * \return exit code (0 - successful, all else - something's wrong)
  */
-static int ParseCommandLineOptions(int argc, char **argv, molecule *&mol, periodentafel *&periode, config& configuration, char *&ConfigFileName, char *&PathToDatabases)
+static int ParseCommandLineOptions(int argc, char **argv, MoleculeListClass *&molecules, periodentafel *&periode, config& configuration, char *&ConfigFileName, char *&PathToDatabases)
 {
 	Vector x,y,z,n;	// coordinates for absolute point in cell volume
@@ -789,4 +1233,8 @@
 	int argptr;
 	PathToDatabases = LocalPath;
+
+	// simply create a new molecule, wherein the config file is loaded and the manipulation takes place
+	molecule *mol = new molecule(periode);
+	molecules->insert(mol);
 
 	if (argc > 1) { // config file specified as option
@@ -1327,5 +1775,5 @@
 		} while (argptr < argc);
 		if (SaveFlag)
-			SaveConfig(ConfigFileName, &configuration, periode, mol);
+			SaveConfig(ConfigFileName, &configuration, periode, molecules);
 		if ((ExitFlag >= 1)) {
 			delete(mol);
@@ -1348,40 +1796,36 @@
 {
 	periodentafel *periode = new periodentafel; // and a period table of all elements
-	molecule *mol = new molecule(periode);		// first we need an empty molecule
+	MoleculeListClass *molecules = new MoleculeListClass;  // list of all molecules
+	molecule *mol = NULL;
 	config configuration;
 	double tmp1;
-	double bond, min_bond;
 	atom *first, *second;
 	char choice;	// menu choice char
 	Vector x,y,z,n;	// coordinates for absolute point in cell volume
-	double *factor;	// unit factor if desired
 	bool valid; // flag if input was valid or not
 	ifstream test;
 	ofstream output;
 	string line;
-	char filename[MAXSTRINGSIZE];
 	char *ConfigFileName = NULL;
 	char *ElementsFileName = NULL;
 	int Z;
 	int j, axis, count, faktor;
-	clock_t start,end;
-//	int *MinimumRingSize = NULL;
-	MoleculeLeafClass *Subgraphs = NULL;
-//	class StackClass<bond *> *BackEdgeStack = NULL;
-	element **Elements;
-	Vector **vectors;
 
 	// =========================== PARSE COMMAND LINE OPTIONS ====================================
-	j = ParseCommandLineOptions(argc, argv, mol, periode, configuration, ConfigFileName, ElementsFileName);
+	j = ParseCommandLineOptions(argc, argv, molecules, periode, configuration, ConfigFileName, ElementsFileName);
 	if (j == 1) return 0; // just for -v and -h options
 	if (j) return j;	// something went wrong
 
 	// General stuff
-	if (mol->cell_size[0] == 0.) {
-		cout << Verbose(0) << "enter lower triadiagonal form of basis matrix" << endl << endl;
-		for (int i=0;i<6;i++) {
-			cout << Verbose(1) << "Cell size" << i << ": ";
-			cin >> mol->cell_size[i];
-		}
+	if (molecules->ListOfMolecules.size() == 0) {
+    mol = new molecule(periode);
+    if (mol->cell_size[0] == 0.) {
+      cout << Verbose(0) << "enter lower tridiagonal form of basis matrix" << endl << endl;
+      for (int i=0;i<6;i++) {
+        cout << Verbose(1) << "Cell size" << i << ": ";
+        cin >> mol->cell_size[i];
+      }
+    }
+    molecules->insert(mol);
 	}
 
@@ -1392,206 +1836,68 @@
 	do {
 		cout << Verbose(0) << endl << endl;
-		cout << Verbose(0) << "============Element list=======================" << endl;
-		mol->Checkout((ofstream *)&cout);
-		cout << Verbose(0) << "============Atom list==========================" << endl;
-		mol->Output((ofstream *)&cout);
+		cout << Verbose(0) << "============Molecule list=======================" << endl;
+		molecules->Enumerate((ofstream *)&cout);
 		cout << Verbose(0) << "============Menu===============================" << endl;
-		cout << Verbose(0) << "a - add an atom" << endl;
-		cout << Verbose(0) << "r - remove an atom" << endl;
-		cout << Verbose(0) << "b - scale a bond between atoms" << endl;
-		cout << Verbose(0) << "u - change an atoms element" << endl;
-		cout << Verbose(0) << "l - measure lengths, angles, ... for an atom" << endl;
-		cout << Verbose(0) << "-----------------------------------------------" << endl;
-		cout << Verbose(0) << "p - Parse xyz file" << endl;
-		cout << Verbose(0) << "e - edit the current configuration" << endl;
-		cout << Verbose(0) << "o - create connection matrix" << endl;
-		cout << Verbose(0) << "f - fragment molecule many-body bond order style" << endl;
-		cout << Verbose(0) << "-----------------------------------------------" << endl;
-		cout << Verbose(0) << "d - duplicate molecule/periodic cell" << endl;
-		cout << Verbose(0) << "i - realign molecule" << endl;
-		cout << Verbose(0) << "m - mirror all molecules" << endl;
-		cout << Verbose(0) << "t - translate molecule by vector" << endl;
-		cout << Verbose(0) << "c - scale by unit transformation" << endl;
-		cout << Verbose(0) << "g - center atoms in box" << endl;
-		cout << Verbose(0) << "-----------------------------------------------" << endl;
-		cout << Verbose(0) << "s - save current setup to config file" << endl;
-		cout << Verbose(0) << "T - call the current test routine" << endl;
-		cout << Verbose(0) << "q - quit" << endl;
-		cout << Verbose(0) << "===============================================" << endl;
-		cout << Verbose(0) << "Input: ";
-		cin >> choice;
+    cout << Verbose(0) << "a - set molecule (in)active" << endl;
+    cout << Verbose(0) << "e - edit new molecules" << endl;
+    cout << Verbose(0) << "g - globally manipulate atoms in molecule" << endl;
+    cout << Verbose(0) << "M - Merge molecules" << endl;
+    cout << Verbose(0) << "m - manipulate atoms" << endl;
+    cout << Verbose(0) << "-----------------------------------------------" << endl;
+    cout << Verbose(0) << "c - edit the current configuration" << endl;
+    cout << Verbose(0) << "-----------------------------------------------" << endl;
+    cout << Verbose(0) << "s - save current setup to config file" << endl;
+    cout << Verbose(0) << "T - call the current test routine" << endl;
+    cout << Verbose(0) << "q - quit" << endl;
+    cout << Verbose(0) << "===============================================" << endl;
+    cout << Verbose(0) << "Input: ";
+    cin >> choice;
 
 		switch (choice) {
-			default:
-			case 'a': // add atom
-				AddAtoms(periode, mol);
-				choice = 'a';
-				break;
-
-			case 'b': // scale a bond
-				cout << Verbose(0) << "Scaling bond length between two atoms." << endl;
-				first = mol->AskAtom("Enter first (fixed) atom: ");
-				second = mol->AskAtom("Enter second (shifting) atom: ");
-				min_bond = 0.;
-				for (int i=NDIM;i--;)
-					min_bond += (first->x.x[i]-second->x.x[i])*(first->x.x[i] - second->x.x[i]);
-				min_bond = sqrt(min_bond);
-				cout << Verbose(0) << "Current Bond length between " << first->type->name << " Atom " << first->nr << " and " << second->type->name << " Atom " << second->nr << ": " << min_bond << " a.u." << endl;
-				cout << Verbose(0) << "Enter new bond length [a.u.]: ";
-				cin >> bond;
-				for (int i=NDIM;i--;) {
-					second->x.x[i] -= (second->x.x[i]-first->x.x[i])/min_bond*(min_bond-bond);
-				}
-				//cout << Verbose(0) << "New coordinates of Atom " << second->nr << " are: ";
-				//second->Output(second->type->No, 1, (ofstream *)&cout);
-				break;
-
-			case 'c': // unit scaling of the metric
-			 cout << Verbose(0) << "Angstroem -> Bohrradius: 1.8897261\t\tBohrradius -> Angstroem: 0.52917721" << endl;
-			 cout << Verbose(0) << "Enter three factors: ";
-			 factor = new double[NDIM];
-			 cin >> factor[0];
-			 cin >> factor[1];
-			 cin >> factor[2];
-			 valid = true;
-			 mol->Scale(&factor);
-			 delete[](factor);
+			case 'a':  // (in)activate molecule
+        {
+          cout << "Enter index of molecule: ";
+          cin >> j;
+          count = 1;
+          MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin();
+          for(; ((ListRunner != molecules->ListOfMolecules.end()) && (count < j)); ListRunner++);
+          if (count == j)
+            (*ListRunner)->ActiveFlag = !(*ListRunner)->ActiveFlag;
+        }
+			  break;
+
+			case 'c': // edit each field of the configuration
+			 configuration.Edit();
 			 break;
 
-			case 'd': // duplicate the periodic cell along a given axis, given times
-				cout << Verbose(0) << "State the axis [(+-)123]: ";
-				cin >> axis;
-				cout << Verbose(0) << "State the factor: ";
-				cin >> faktor;
-
-				mol->CountAtoms((ofstream *)&cout);	// recount atoms
-				if (mol->AtomCount != 0) {	// if there is more than none
-					count = mol->AtomCount;	 // is changed becausing of adding, thus has to be stored away beforehand
-					Elements = new element *[count];
-					vectors = new Vector *[count];
-					j = 0;
-					first = mol->start;
-					while (first->next != mol->end) {	// make a list of all atoms with coordinates and element
-						first = first->next;
-						Elements[j] = first->type;
-						vectors[j] = &first->x;
-						j++;
-					}
-					if (count != j)
-						cout << Verbose(0) << "ERROR: AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl;
-					x.Zero();
-					y.Zero();
-					y.x[abs(axis)-1] = mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] * abs(axis)/axis; // last term is for sign, first is for magnitude
-					for (int i=1;i<faktor;i++) {	// then add this list with respective translation factor times
-						x.AddVector(&y); // per factor one cell width further
-						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.AddVector(&x);			// translate the coordinates
-							first->type = Elements[k];	// insert original element
-							mol->AddAtom(first);				// and add to the molecule (which increments ElementsInMolecule, AtomCount, ...)
-						}
-					}
-					if (mol->first->next != mol->last) // if connect matrix is present already, redo it
-						mol->CreateAdjacencyList((ofstream *)&cout, mol->BondDistance, configuration.GetIsAngstroem());
-					// free memory
-					delete[](Elements);
-					delete[](vectors);
-					// correct cell size
-					if (axis < 0) { // if sign was negative, we have to translate everything
-						x.Zero();
-						x.AddVector(&y);
-						x.Scale(-(faktor-1));
-						mol->Translate(&x);
-					}
-					mol->cell_size[(abs(axis) == 2) ? 2 : ((abs(axis) == 3) ? 5 : 0)] *= faktor;
-				}
-				break;
-
-			case 'e': // edit each field of the configuration
-			 configuration.Edit(mol);
-			 break;
-
-			case 'f':
-				FragmentAtoms(mol, &configuration);
-				break;
-
-			case 'g': // center the atoms
-				CenterAtoms(mol);
-				break;
-
-			case 'i': // align all atoms
-				AlignAtoms(periode, mol);
-				break;
-
-			case 'l': // measure distances or angles
-				MeasureAtoms(periode, mol, &configuration);
-				break;
-
-			case 'm': // mirror atoms along a given axis
-				MirrorAtoms(mol);
-				break;
-
-			case 'o': // create the connection matrix
-				{
-					cout << Verbose(0) << "What's the maximum bond distance: ";
-					cin >> tmp1;
-					start = clock();
-					mol->CreateAdjacencyList((ofstream *)&cout, tmp1, configuration.GetIsAngstroem());
-					mol->CreateListOfBondsPerAtom((ofstream *)&cout);
-//					Subgraphs = mol->DepthFirstSearchAnalysis((ofstream *)&cout, BackEdgeStack);
-//					while (Subgraphs->next != NULL) {
-//						Subgraphs = Subgraphs->next;
-//						Subgraphs->Leaf->CyclicStructureAnalysis((ofstream *)&cout, BackEdgeStack, MinimumRingSize);
-//						delete(Subgraphs->previous);
-//					}
-//					delete(Subgraphs);		// we don't need the list here, so free everything
-//					delete[](MinimumRingSize);
-//					Subgraphs = NULL;
-					end = clock();
-					cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
-				}
-				break;
-
-			case 'p': // parse and XYZ file
-				cout << Verbose(0) << "Format should be XYZ with: ShorthandOfElement\tX\tY\tZ" << endl;
-				do {
-					cout << Verbose(0) << "Enter file name: ";
-					cin >> filename;
-				} while (!mol->AddXYZFile(filename));
-				break;
+      case 'g': // manipulate molecules
+        ManipulateMolecules(periode, molecules, &configuration);
+        break;
+
+      case 'M':  // merge molecules
+        MergeMolecules(periode, molecules);
+        break;
+
+      case 'm': // manipulate atoms
+        ManipulateAtoms(periode, molecules, &configuration);
+        break;
+
+      case 'e': // create molecule
+        EditMolecules(periode, molecules);
+        break;
 
 			case 'q': // quit
 				break;
 
-			case 'r': // remove atom
-				RemoveAtoms(mol);
+			case 's': // save to config file
+				SaveConfig(ConfigFileName, &configuration, periode, molecules);
 				break;
 
-			case 's': // save to config file
-				SaveConfig(ConfigFileName, &configuration, periode, mol);
+			case 'T':
+				testroutine(molecules);
 				break;
 
-			case 't': // translate all atoms
-			 cout << Verbose(0) << "Enter translation vector." << endl;
-			 x.AskPosition(mol->cell_size,0);
-			 mol->Translate((const Vector *)&x);
-			 break;
-
-			case 'T':
-				testroutine(mol);
-				break;
-
-			case 'u': // change an atom's element
-				first = NULL;
-				do {
-					cout << Verbose(0) << "Change the element of which atom: ";
-					cin >> Z;
-				} while ((first = mol->FindAtom(Z)) == NULL);
-				cout << Verbose(0) << "New element by atomic number Z: ";
-				cin >> Z;
-				first->type = periode->FindElement(Z);
-				cout << Verbose(0) << "Atom " << first->nr << "'s element is " << first->type->name << "." << endl;
-				break;
+			default:
+			  break;
 		};
 	} while (choice != 'q');
@@ -1603,13 +1909,5 @@
 		cout << Verbose(0) << "Saving of elements.db failed." << endl;
 
-	// Free all
-	if (Subgraphs != NULL) {	// free disconnected subgraph list of DFS analysis was performed
-		while (Subgraphs->next != NULL) {
-			Subgraphs = Subgraphs->next;
-			delete(Subgraphs->previous);
-		}
-		delete(Subgraphs);
-	}
-	delete(mol);
+	delete(molecules);
 	delete(periode);
 	return (0);
Index: src/config.cpp
===================================================================
--- src/config.cpp	(revision d8b94a86944c4679c88bc250bff510f242d63309)
+++ src/config.cpp	(revision 1907a789609146bbdebfeff09de44cb4ba21b808)
@@ -1,6 +1,6 @@
 /** \file config.cpp
- * 
+ *
  * Function implementations for the class config.
- * 
+ *
  */
 
@@ -24,5 +24,5 @@
 	configname[0]='\0';
 	basis = "3-21G";
-	
+
 	FastParsing = false;
 	ProcPEGamma=8;
@@ -42,5 +42,5 @@
 	UseAddGramSch=1;
 	Seed=1;
-	
+
 	MaxOuterStep=0;
 	Deltat=1;
@@ -51,5 +51,5 @@
 	MaxPsiStep=0;
 	EpsWannier=1e-7;
-	
+
 	MaxMinStep=100;
 	RelEpsTotalEnergy=1e-7;
@@ -62,7 +62,7 @@
 	InitMaxMinStopStep=1;
 	InitMaxMinGapStopStep=0;
-	
+
 	//BoxLength[NDIM*NDIM];
-	
+
 	ECut=128.;
 	MaxLevel=5;
@@ -77,5 +77,5 @@
 	PsiMaxNoDown=0;
 	AddPsis=0;
-	
+
 	RCut=20.;
 	StructOpt=0;
@@ -101,8 +101,8 @@
  * for each entry of the config file structure.
  */
-void config::Edit(molecule *mol)
+void config::Edit()
 {
 	char choice;
-	
+
 	do {
 		cout << Verbose(0) << "===========EDIT CONFIGURATION============================" << endl;
@@ -137,5 +137,5 @@
 		cout << Verbose(0) << " g - Relative change in kinetic energy to stop min. iteration during initial level" << endl;
 		cout << Verbose(0) << " h - Check stop conditions every ..th step during min. iteration during initial level" << endl;
-		cout << Verbose(0) << " j - six lower diagonal entries of matrix, defining the unit cell" << endl;
+//		cout << Verbose(0) << " j - six lower diagonal entries of matrix, defining the unit cell" << endl;
 		cout << Verbose(0) << " k - Energy cutoff of plane wave basis in Hartree" << endl;
 		cout << Verbose(0) << " l - Maximum number of levels in multi-level-ansatz" << endl;
@@ -157,5 +157,5 @@
 		cout << Verbose(0) << "INPUT: ";
 		cin >> choice;
-		
+
 		switch (choice) {
 				case 'A': // mainname
@@ -282,13 +282,13 @@
 					cin >> config::InitMaxMinStopStep;
 					break;
-				
-				case 'j': // BoxLength
-					cout << Verbose(0) << "enter lower triadiagonalo form of basis matrix" << endl << endl;
-					for (int i=0;i<6;i++) {
-						cout << Verbose(0) << "Cell size" << i << ": ";
-						cin >> mol->cell_size[i];
-					}
-					break;
-				
+
+//				case 'j': // BoxLength
+//					cout << Verbose(0) << "enter lower triadiagonalo form of basis matrix" << endl << endl;
+//					for (int i=0;i<6;i++) {
+//						cout << Verbose(0) << "Cell size" << i << ": ";
+//						cin >> mol->cell_size[i];
+//					}
+//					break;
+
 				case 'k': // ECut
 					cout << Verbose(0) << "Old: " << config::ECut << "\t new: ";
@@ -364,5 +364,5 @@
  * \param *periode pointer to a periodentafel class with all elements
  * \param *mol pointer to molecule containing all atoms of the molecule
- * \return 0 - old syntax, 1 - new syntax, -1 - unknown syntax 
+ * \return 0 - old syntax, 1 - new syntax, -1 - unknown syntax
  */
 int config::TestSyntax(char *filename, periodentafel *periode, molecule *mol)
@@ -370,5 +370,5 @@
 	int test;
 	ifstream file(filename);
-	
+
 	// search file for keyword: ProcPEGamma (new syntax)
 	if (ParseForParameter(1,&file,"ProcPEGamma", 0, 1, 1, int_type, &test, 1, optional)) {
@@ -441,5 +441,5 @@
  * \param *file input file stream being the opened config file
  * \param *periode pointer to a periodentafel class with all elements
- * \param *mol pointer to molecule containing all atoms of the molecule 
+ * \param *mol pointer to molecule containing all atoms of the molecule
  */
 void config::Load(char *filename, periodentafel *periode, molecule *mol)
@@ -452,5 +452,5 @@
 	RetrieveConfigPathAndName(filename);
 	// ParseParameters
-	
+
 	/* Oeffne Hauptparameterdatei */
 	int di;
@@ -464,5 +464,5 @@
 	int verbose = 0;
 	double value[3];
-	
+
 	/* Namen einlesen */
 
@@ -495,5 +495,5 @@
 		config::DoOutCurrent = 0;
 	if (config::DoOutCurrent < 0) config::DoOutCurrent = 0;
-	if (config::DoOutCurrent > 1) config::DoOutCurrent = 1; 
+	if (config::DoOutCurrent > 1) config::DoOutCurrent = 1;
 	ParseForParameter(verbose,file,"AddGramSch", 0, 1, 1, int_type, &(config::UseAddGramSch), 1, critical);
 	if (config::UseAddGramSch < 0) config::UseAddGramSch = 0;
@@ -517,5 +517,5 @@
 		if (config::SawtoothStart > 1.) config::SawtoothStart = 1.;
 	}
-	
+
 	ParseForParameter(verbose,file,"MaxOuterStep", 0, 1, 1, int_type, &(config::MaxOuterStep), 1, critical);
 	if (!ParseForParameter(verbose,file,"Deltat", 0, 1, 1, double_type, &(config::Deltat), 1, optional))
@@ -527,10 +527,10 @@
 	if (!ParseForParameter(verbose,file,"EpsWannier", 0, 1, 1, double_type, &(config::EpsWannier), 1, optional))
 		config::EpsWannier = 1e-8;
-	
+
 	// stop conditions
 	//if (config::MaxOuterStep <= 0) config::MaxOuterStep = 1;
 	ParseForParameter(verbose,file,"MaxPsiStep", 0, 1, 1, int_type, &(config::MaxPsiStep), 1, critical);
 	if (config::MaxPsiStep <= 0) config::MaxPsiStep = 3;
-	
+
 	ParseForParameter(verbose,file,"MaxMinStep", 0, 1, 1, int_type, &(config::MaxMinStep), 1, critical);
 	ParseForParameter(verbose,file,"RelEpsTotalE", 0, 1, 1, double_type, &(config::RelEpsTotalEnergy), 1, critical);
@@ -541,5 +541,5 @@
 	if (config::MaxMinStopStep < 1) config::MaxMinStopStep = 1;
 	if (config::MaxMinGapStopStep < 1) config::MaxMinGapStopStep = 1;
-	
+
 	ParseForParameter(verbose,file,"MaxInitMinStep", 0, 1, 1, int_type, &(config::MaxInitMinStep), 1, critical);
 	ParseForParameter(verbose,file,"InitRelEpsTotalE", 0, 1, 1, double_type, &(config::InitRelEpsTotalEnergy), 1, critical);
@@ -550,5 +550,5 @@
 	if (config::InitMaxMinStopStep < 1) config::InitMaxMinStopStep = 1;
 	if (config::InitMaxMinGapStopStep < 1) config::InitMaxMinGapStopStep = 1;
-	
+
 	// Unit cell and magnetic field
 	ParseForParameter(verbose,file, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
@@ -566,7 +566,7 @@
 		config::DoFullCurrent = 0;
 	if (config::DoFullCurrent < 0) config::DoFullCurrent = 0;
-	if (config::DoFullCurrent > 2) config::DoFullCurrent = 2; 
+	if (config::DoFullCurrent > 2) config::DoFullCurrent = 2;
 	if (config::DoOutNICS < 0) config::DoOutNICS = 0;
-	if (config::DoOutNICS > 2) config::DoOutNICS = 2; 
+	if (config::DoOutNICS > 2) config::DoOutNICS = 2;
 	if (config::DoPerturbation == 0) {
 		config::DoFullCurrent = 0;
@@ -585,5 +585,5 @@
 	} else {
 		fprintf(stderr, "0 <= RiemanTensor < 2: 0 UseNotRT, 1 UseRT");
-		exit(1); 
+		exit(1);
 	}
 	switch (config::RiemannTensor) {
@@ -602,5 +602,5 @@
 			if (config::RiemannLevel < 2) {
 				config::RiemannLevel = 2;
-			} 
+			}
 			if (config::RiemannLevel > config::MaxLevel-1) {
 				config::RiemannLevel = config::MaxLevel-1;
@@ -609,5 +609,5 @@
 			if (config::LevRFactor < 2) {
 				config::LevRFactor = 2;
-			} 
+			}
 			config::Lev0Factor = 2;
 			config::RTActualUse = 2;
@@ -619,5 +619,5 @@
 	} else {
 		fprintf(stderr, "0 <= PsiType < 2: 0 UseSpinDouble, 1 UseSpinUpDown");
-		exit(1); 
+		exit(1);
 	}
 	switch (config::PsiType) {
@@ -633,7 +633,7 @@
 		break;
 	}
-	
+
 	// IonsInitRead
-	
+
 	ParseForParameter(verbose,file,"RCut", 0, 1, 1, double_type, &(config::RCut), 1, critical);
 	ParseForParameter(verbose,file,"IsAngstroem", 0, 1, 1, int_type, &(config::IsAngstroem), 1, critical);
@@ -653,8 +653,8 @@
 			ParseForParameter(verbose,file, name, 0, 2, 1, int_type, &Z, 1, critical);
 			elementhash[i] = periode->FindElement(Z);
-			cout << Verbose(1) << i << ". Z = " << elementhash[i]->Z << " with " << No[i] << " ions." << endl; 
+			cout << Verbose(1) << i << ". Z = " << elementhash[i]->Z << " with " << No[i] << " ions." << endl;
 		}
 		int repetition = 0; // which repeated keyword shall be read
-	
+
 		map<int, atom *> AtomList[config::MaxTypes];
 		if (!FastParsing) {
@@ -675,5 +675,5 @@
 						} else
 							neues = AtomList[i][j];
-						status = (status && 
+						status = (status &&
 										ParseForParameter(verbose,file, keyword, 0, 1, 1, double_type, &neues->x.x[0], 1, (repetition == 0) ? critical : optional) &&
 										ParseForParameter(verbose,file, keyword, 0, 2, 1, double_type, &neues->x.x[1], 1, (repetition == 0) ? critical : optional) &&
@@ -681,5 +681,5 @@
 										ParseForParameter(verbose,file, keyword, 0, 4, 1, int_type, &neues->FixedIon, 1, (repetition == 0) ? critical : optional));
 						if (!status) break;
-		
+
 						// check size of vectors
 						if (mol->Trajectories[neues].R.size() <= (unsigned int)(repetition)) {
@@ -689,9 +689,9 @@
 							mol->Trajectories[neues].F.resize(repetition+10);
 						}
-					
+
 						// put into trajectories list
 						for (int d=0;d<NDIM;d++)
 							mol->Trajectories[neues].R.at(repetition).x[d] = neues->x.x[d];
-						
+
 						// parse velocities if present
 						if(!ParseForParameter(verbose,file, keyword, 0, 5, 1, double_type, &neues->v.x[0], 1,optional))
@@ -703,5 +703,5 @@
 						for (int d=0;d<NDIM;d++)
 							mol->Trajectories[neues].U.at(repetition).x[d] = neues->v.x[d];
-			
+
 						// parse forces if present
 						if(!ParseForParameter(verbose,file, keyword, 0, 8, 1, double_type, &value[0], 1,optional))
@@ -713,6 +713,6 @@
 						for (int d=0;d<NDIM;d++)
 							mol->Trajectories[neues].F.at(repetition).x[d] = value[d];
-	
-	//						cout << "Parsed position of step " << (repetition) << ": ("; 
+
+	//						cout << "Parsed position of step " << (repetition) << ": (";
 	//						for (int d=0;d<NDIM;d++)
 	//							cout << mol->Trajectories[neues].R.at(repetition).x[d] << " ";					// next step
@@ -770,5 +770,5 @@
  * \param *file input file stream being the opened config file with old pcp syntax
  * \param *periode pointer to a periodentafel class with all elements
- * \param *mol pointer to molecule containing all atoms of the molecule 
+ * \param *mol pointer to molecule containing all atoms of the molecule
  */
 void config::LoadOld(char *filename, periodentafel *periode, molecule *mol)
@@ -781,5 +781,5 @@
 	RetrieveConfigPathAndName(filename);
 	// ParseParameters
-	
+
 	/* Oeffne Hauptparameterdatei */
 	int l, i, di;
@@ -791,5 +791,5 @@
 	int Z, No, AtomNo, found;
 	int verbose = 0;
-	
+
 	/* Namen einlesen */
 
@@ -823,10 +823,10 @@
 	ParseForParameter(verbose,file,"ScaleTempStep", 0, 1, 1, int_type, &(config::ScaleTempStep), 1, optional);
 	config::EpsWannier = 1e-8;
-	
+
 	// stop conditions
 	//if (config::MaxOuterStep <= 0) config::MaxOuterStep = 1;
 	ParseForParameter(verbose,file,"MaxPsiStep", 0, 1, 1, int_type, &(config::MaxPsiStep), 1, critical);
 	if (config::MaxPsiStep <= 0) config::MaxPsiStep = 3;
-	
+
 	ParseForParameter(verbose,file,"MaxMinStep", 0, 1, 1, int_type, &(config::MaxMinStep), 1, critical);
 	ParseForParameter(verbose,file,"MaxMinStep", 0, 2, 1, double_type, &(config::RelEpsTotalEnergy), 1, critical);
@@ -836,5 +836,5 @@
 	if (config::MaxMinStopStep < 1) config::MaxMinStopStep = 1;
 	config::MaxMinGapStopStep = 1;
-	
+
 	ParseForParameter(verbose,file,"MaxInitMinStep", 0, 1, 1, int_type, &(config::MaxInitMinStep), 1, critical);
 	ParseForParameter(verbose,file,"MaxInitMinStep", 0, 2, 1, double_type, &(config::InitRelEpsTotalEnergy), 1, critical);
@@ -867,5 +867,5 @@
 	} else {
 		fprintf(stderr, "0 <= RiemanTensor < 2: 0 UseNotRT, 1 UseRT");
-		exit(1); 
+		exit(1);
 	}
 	switch (config::RiemannTensor) {
@@ -884,5 +884,5 @@
 			if (config::RiemannLevel < 2) {
 				config::RiemannLevel = 2;
-			} 
+			}
 			if (config::RiemannLevel > config::MaxLevel-1) {
 				config::RiemannLevel = config::MaxLevel-1;
@@ -891,5 +891,5 @@
 			if (config::LevRFactor < 2) {
 				config::LevRFactor = 2;
-			} 
+			}
 			config::Lev0Factor = 2;
 			config::RTActualUse = 2;
@@ -901,5 +901,5 @@
 	} else {
 		fprintf(stderr, "0 <= PsiType < 2: 0 UseSpinDouble, 1 UseSpinUpDown");
-		exit(1); 
+		exit(1);
 	}
 	switch (config::PsiType) {
@@ -915,7 +915,7 @@
 		break;
 	}
-	
+
 	// IonsInitRead
-	
+
 	ParseForParameter(verbose,file,"RCut", 0, 1, 1, double_type, &(config::RCut), 1, critical);
 	ParseForParameter(verbose,file,"IsAngstroem", 0, 1, 1, int_type, &(config::IsAngstroem), 1, critical);
@@ -924,7 +924,7 @@
 
 	// Routine from builder.cpp
-	
-	
-	for (i=MAX_ELEMENTS;i--;) 
+
+
+	for (i=MAX_ELEMENTS;i--;)
 		elementhash[i] = NULL;
 	cout << Verbose(0) << "Parsing Ions ..." << endl;
@@ -937,5 +937,5 @@
 		}
 		if (found > 0) {
-			if (zeile.find("Ions_Data") == 0) 
+			if (zeile.find("Ions_Data") == 0)
 				getline(*file,zeile,'\n'); // read next line and parse this one
 			istringstream input(zeile);
@@ -966,5 +966,5 @@
 			No++;
 		}
-	}	 
+	}
 	file->close();
 	delete(file);
@@ -974,5 +974,5 @@
  * \param *filename name of file
  * \param *periode pointer to a periodentafel class with all elements
- * \param *mol pointer to molecule containing all atoms of the molecule 
+ * \param *mol pointer to molecule containing all atoms of the molecule
  */
 bool config::Save(const char *filename, periodentafel *periode, molecule *mol) const
@@ -1081,5 +1081,5 @@
  * Note that this format cannot be parsed again.
  * \param *filename name of file (without ".in" suffix!)
- * \param *mol pointer to molecule containing all atoms of the molecule 
+ * \param *mol pointer to molecule containing all atoms of the molecule
  */
 bool config::SaveMPQC(const char *filename, molecule *mol) const
@@ -1092,5 +1092,5 @@
 	ofstream *output = NULL;
 	stringstream *fname = NULL;
-	
+
 	// first without hessian
 	fname = new stringstream;
@@ -1191,5 +1191,5 @@
 	delete(output);
 	delete(fname);
-	
+
 	return true;
 };
@@ -1203,5 +1203,5 @@
  * \param file file to be parsed
  * \param name Name of value in file (at least 3 chars!)
- * \param sequential 1 - do not reset file pointer to begin of file, 0 - set to beginning 
+ * \param sequential 1 - do not reset file pointer to begin of file, 0 - set to beginning
  *				(if file is sequentially parsed this can be way faster! However, beware of multiple values per line, as whole line is read -
  *				 best approach: 0 0 0 1 (not resetted only on last value of line) - and of yth, which is now
@@ -1223,5 +1223,5 @@
 	dummy1 = free_dummy = (char *) Malloc(256 * sizeof(char), "config::ParseForParameter: *free_dummy");
 
-	//fprintf(stderr,"Parsing for %s\n",name);	
+	//fprintf(stderr,"Parsing for %s\n",name);
 	if (repetition == 0)
 		//Error(SomeError, "ParseForParameter(): argument repetition must not be 0!");
@@ -1244,5 +1244,5 @@
 					file->clear();
 					file->seekg(file_position, ios::beg);	// rewind to start position
-					Free((void **)&free_dummy, "config::ParseForParameter: *free_dummy");					
+					Free((void **)&free_dummy, "config::ParseForParameter: *free_dummy");
 					return 0;
 				}
@@ -1250,5 +1250,5 @@
 			line++;
 		} while (dummy != NULL && dummy1 != NULL && ((dummy1[0] == '#') || (dummy1[0] == '\0'))); // skip commentary and empty lines
-		
+
 		// C++ getline removes newline at end, thus re-add
 		if ((dummy1 != NULL) && (strchr(dummy1,'\n') == NULL)) {
@@ -1276,5 +1276,5 @@
 				//fprintf(stderr,"Error: Cannot find tabs or spaces on line %i in search for %s\n", line, name);
 				//Free((void **)&free_dummy);
-				//Error(FileOpenParams, NULL);			
+				//Error(FileOpenParams, NULL);
 			} else {
 				//fprintf(stderr,"found tab at %i\n",(char *)dummy-(char *)dummy1);
@@ -1285,6 +1285,6 @@
 		if ((name == NULL) || (((dummy-dummy1 >= 3) && (strncmp(dummy1, name, strlen(name)) == 0)) && ((unsigned int)(dummy-dummy1) == strlen(name)))) {
 			found++; // found the parameter!
-			//fprintf(stderr,"found %s at line %i between %i and %i\n", name, line, dummy1, dummy);		
-			
+			//fprintf(stderr,"found %s at line %i between %i and %i\n", name, line, dummy1, dummy);
+
 			if (found == repetition) {
 				for (i=0;i<xth;i++) { // i = rows
@@ -1325,5 +1325,5 @@
 						dummy1[j+1] = '\0';
 					}
-	
+
 					int start = (type >= grid) ? 0 : yth-1 ;
 					for (j=start;j<yth;j++) { // j = columns
@@ -1348,5 +1348,5 @@
 										//return 0;
 										exit(255);
-										//Error(FileOpenParams, NULL);			
+										//Error(FileOpenParams, NULL);
 									} else {
 										//if (!sequential)
@@ -1363,5 +1363,5 @@
 								// found comment, skipping rest of line
 								//if (verbose) fprintf(stderr,"Error: '#' at %i and still missing %i value(s) for parameter %s\n", line, yth-j, name);
-								if (!sequential) { // here we need it! 
+								if (!sequential) { // here we need it!
 									file->seekg(file_position, ios::beg);	// rewind to start position
 								}
@@ -1378,5 +1378,5 @@
 										//value += sizeof(int);
 									break;
-								case(row_double): 
+								case(row_double):
 								case(grid):
 								case(lower_trigrid):
@@ -1419,5 +1419,5 @@
 			}
 		}
-	} 
+	}
 	if ((type >= row_int) && (verbose)) fprintf(stderr,"\n");
 	Free((void **)&free_dummy, "config::ParseForParameter: *free_dummy");
@@ -1427,5 +1427,5 @@
 	}
 	//fprintf(stderr, "End of Parsing\n\n");
-	
+
 	return (found); // true if found, false if not
 }
Index: src/helpers.cpp
===================================================================
--- src/helpers.cpp	(revision d8b94a86944c4679c88bc250bff510f242d63309)
+++ src/helpers.cpp	(revision 1907a789609146bbdebfeff09de44cb4ba21b808)
@@ -170,5 +170,5 @@
 	}
 	// allocate string
-	returnstring = (char *) Malloc(sizeof(char)*(order+2), "MoleculeListClass::CreateFragmentNumberForOutput: *returnstring");
+	returnstring = (char *) Malloc(sizeof(char)*(order+2), "FixedDigitNumber: *returnstring");
 	// terminate	and fill string array from end backward
 	returnstring[order] = '\0';
Index: src/moleculelist.cpp
===================================================================
--- src/moleculelist.cpp	(revision d8b94a86944c4679c88bc250bff510f242d63309)
+++ src/moleculelist.cpp	(revision 1907a789609146bbdebfeff09de44cb4ba21b808)
@@ -1,6 +1,6 @@
 /** \file MoleculeListClass.cpp
- * 
+ *
  * Function implementations for the class MoleculeListClass.
- * 
+ *
  */
 
@@ -13,19 +13,8 @@
 MoleculeListClass::MoleculeListClass()
 {
-};
-
-/** constructor for MoleculeListClass.
- * \param NumMolecules number of molecules to allocate for
- * \param NumAtoms number of atoms to allocate for
- */
-MoleculeListClass::MoleculeListClass(int NumMolecules, int NumAtoms)
-{
-	ListOfMolecules = (molecule **) Malloc(sizeof(molecule *)*NumMolecules, "MoleculeListClass:MoleculeListClass: **ListOfMolecules");
-	for (int i=NumMolecules;i--;)
-		ListOfMolecules[i] = NULL;
-	NumberOfMolecules = NumMolecules;
-	NumberOfTopAtoms = NumAtoms;
-};
-
+  // empty lists
+  ListOfMolecules.clear();
+  MaxIndex = 1;
+};
 
 /** Destructor for MoleculeListClass.
@@ -33,14 +22,21 @@
 MoleculeListClass::~MoleculeListClass()
 {
-	cout << Verbose(3) << this << ": Freeing ListOfMolcules." << endl;
-	for (int i=NumberOfMolecules;i--;) {
-		if (ListOfMolecules[i] != NULL) { // if NULL don't free
-			cout << Verbose(4) << "ListOfMolecules: Freeing " << ListOfMolecules[i] << "." << endl;
-			delete(ListOfMolecules[i]);
-			ListOfMolecules[i] = NULL;
-		}
-	}
-	cout << Verbose(4) << "Freeing ListOfMolecules." << endl;
-	Free((void **)&ListOfMolecules, "MoleculeListClass:MoleculeListClass: **ListOfMolecules");
+  cout << Verbose(3) << this << ": Freeing ListOfMolcules." << endl;
+  for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
+    cout << Verbose(4) << "ListOfMolecules: Freeing " << *ListRunner << "." << endl;
+    delete (*ListRunner);
+  }
+  cout << Verbose(4) << "Freeing ListOfMolecules." << endl;
+  ListOfMolecules.clear(); // empty list
+};
+
+/** Insert a new molecule into the list and set its number.
+ * \param *mol molecule to add to list.
+ * \return true - add successful
+ */
+bool MoleculeListClass::insert(molecule *mol)
+{
+  mol->IndexNr = MaxIndex++;
+  ListOfMolecules.push_back(mol);
 };
 
@@ -52,72 +48,259 @@
 int MolCompare(const void *a, const void *b)
 {
-	int *aList = NULL, *bList = NULL;
-	int Count, Counter, aCounter, bCounter;
-	int flag;
-	atom *aWalker = NULL;
-	atom *bWalker = NULL;
-	
-	// sort each atom list and put the numbers into a list, then go through
-	//cout << "Comparing fragment no. " << *(molecule **)a << " to " << *(molecule **)b << "." << endl;
-	if ( (**(molecule **)a).AtomCount < (**(molecule **)b).AtomCount ) {
-		return -1;
-	} else { if ((**(molecule **)a).AtomCount > (**(molecule **)b).AtomCount)
-		return +1;
-		else {
-			Count = (**(molecule **)a).AtomCount;
-			aList = new int[Count];
-			bList = new int[Count];
-	
-			// fill the lists
-			aWalker = (**(molecule **)a).start;
-			bWalker = (**(molecule **)b).start;
-			Counter = 0;
-			aCounter = 0;
-			bCounter = 0;
-			while ((aWalker->next != (**(molecule **)a).end) && (bWalker->next != (**(molecule **)b).end)) {
-				aWalker = aWalker->next;
-				bWalker = bWalker->next;
-				if (aWalker->GetTrueFather() == NULL)
-					aList[Counter] = Count + (aCounter++);
-				else
-					aList[Counter] = aWalker->GetTrueFather()->nr;
-				if (bWalker->GetTrueFather() == NULL)
-					bList[Counter] = Count + (bCounter++);
-				else
-					bList[Counter] = bWalker->GetTrueFather()->nr;
-				Counter++;
-			}
-			// check if AtomCount was for real
-			flag = 0;
-			if ((aWalker->next == (**(molecule **)a).end) && (bWalker->next != (**(molecule **)b).end)) {
-				flag = -1;
-			} else {
-				if ((aWalker->next != (**(molecule **)a).end) && (bWalker->next == (**(molecule **)b).end))
-					flag = 1;
-			}
-			if (flag == 0) {
-				// sort the lists
-				gsl_heapsort(aList,Count, sizeof(int), CompareDoubles); 
-				gsl_heapsort(bList,Count, sizeof(int), CompareDoubles);
-				// compare the lists 
-				
-				flag = 0;
-				for(int i=0;i<Count;i++) {
-					if (aList[i] < bList[i]) {
-						flag = -1;
-					} else {
-						if (aList[i] > bList[i])
-							flag = 1;
-					}
-					if (flag != 0)
-						break;
-				}
-			}
-			delete[](aList);
-			delete[](bList);
-			return flag;
-		}
-	}
-	return	-1;
+  int *aList = NULL, *bList = NULL;
+  int Count, Counter, aCounter, bCounter;
+  int flag;
+  atom *aWalker = NULL;
+  atom *bWalker = NULL;
+
+  // sort each atom list and put the numbers into a list, then go through
+  //cout << "Comparing fragment no. " << *(molecule **)a << " to " << *(molecule **)b << "." << endl;
+  if ((**(molecule **) a).AtomCount < (**(molecule **) b).AtomCount) {
+    return -1;
+  } else {
+    if ((**(molecule **) a).AtomCount > (**(molecule **) b).AtomCount)
+      return +1;
+    else {
+      Count = (**(molecule **) a).AtomCount;
+      aList = new int[Count];
+      bList = new int[Count];
+
+      // fill the lists
+      aWalker = (**(molecule **) a).start;
+      bWalker = (**(molecule **) b).start;
+      Counter = 0;
+      aCounter = 0;
+      bCounter = 0;
+      while ((aWalker->next != (**(molecule **) a).end) && (bWalker->next != (**(molecule **) b).end)) {
+        aWalker = aWalker->next;
+        bWalker = bWalker->next;
+        if (aWalker->GetTrueFather() == NULL)
+          aList[Counter] = Count + (aCounter++);
+        else
+          aList[Counter] = aWalker->GetTrueFather()->nr;
+        if (bWalker->GetTrueFather() == NULL)
+          bList[Counter] = Count + (bCounter++);
+        else
+          bList[Counter] = bWalker->GetTrueFather()->nr;
+        Counter++;
+      }
+      // check if AtomCount was for real
+      flag = 0;
+      if ((aWalker->next == (**(molecule **) a).end) && (bWalker->next != (**(molecule **) b).end)) {
+        flag = -1;
+      } else {
+        if ((aWalker->next != (**(molecule **) a).end) && (bWalker->next == (**(molecule **) b).end))
+          flag = 1;
+      }
+      if (flag == 0) {
+        // sort the lists
+        gsl_heapsort(aList, Count, sizeof(int), CompareDoubles);
+        gsl_heapsort(bList, Count, sizeof(int), CompareDoubles);
+        // compare the lists
+
+        flag = 0;
+        for (int i = 0; i < Count; i++) {
+          if (aList[i] < bList[i]) {
+            flag = -1;
+          } else {
+            if (aList[i] > bList[i])
+              flag = 1;
+          }
+          if (flag != 0)
+            break;
+        }
+      }
+      delete[] (aList);
+      delete[] (bList);
+      return flag;
+    }
+  }
+  return -1;
+};
+
+/** Output of a list of all molecules.
+ * \param *out output stream
+ */
+void MoleculeListClass::Enumerate(ofstream *out)
+{
+  int i=1;
+  element* Elemental = NULL;
+  atom *Walker = NULL;
+  int Counts[MAX_ELEMENTS];
+
+  // header
+  *out << "Index\tName\tNo.Atoms\tformula" << endl;
+  cout << Verbose(0) << "-----------------------------------------------" << endl;
+  if (ListOfMolecules.size() == 0)
+    *out << "\tNone" << endl;
+  else {
+    for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
+      // reset element counts
+      for (int j = 0; j<MAX_ELEMENTS;j++)
+        Counts[j] = 0;
+      // count atoms per element
+      Walker = (*ListRunner)->start;
+      while (Walker->next != (*ListRunner)->end) {
+        Walker = Walker->next;
+        Counts[Walker->type->Z]++;
+      }
+      // output Index, Name, number of atoms, chemical formula
+      *out << ((*ListRunner)->ActiveFlag ? "*" : " ") << (*ListRunner)->IndexNr << "\t" << (*ListRunner)->name << "\t\t" << (*ListRunner)->AtomCount << "\t";
+      Elemental = (*ListRunner)->elemente->end;
+      while(Elemental != (*ListRunner)->elemente->start) {
+        Elemental = Elemental->previous;
+        if (Counts[Elemental->Z] != 0)
+          *out << Elemental->symbol << Counts[Elemental->Z];
+      }
+      *out << endl;
+    }
+  }
+};
+
+/** Returns the molecule with the given index \a index.
+ * \param index index of the desired molecule
+ * \return pointer to molecule structure, NULL if not found
+ */
+molecule * MoleculeListClass::ReturnIndex(int index)
+{
+  int count = 1;
+  MoleculeList::iterator ListRunner = ListOfMolecules.begin();
+  for(; ((ListRunner != ListOfMolecules.end()) && (count < index)); ListRunner++);
+  if (count == index)
+    return (*ListRunner);
+  else
+    return NULL;
+};
+
+/** Simple merge of two molecules into one.
+ * \param *mol destination molecule
+ * \param *srcmol source molecule
+ * \return true - merge successful, false - merge failed (probably due to non-existant indices
+ */
+bool MoleculeListClass::SimpleMerge(molecule *mol, molecule *srcmol)
+{
+  if (srcmol == NULL)
+    return false;
+
+  // put all molecules of src into mol
+  atom *Walker = srcmol->start;
+  atom *NextAtom = Walker->next;
+  while (NextAtom != srcmol->end) {
+    Walker = NextAtom;
+    NextAtom = Walker->next;
+    srcmol->UnlinkAtom(Walker);
+    mol->AddAtom(Walker);
+  }
+
+  // remove src
+  ListOfMolecules.remove(srcmol);
+  delete(srcmol);
+  return true;
+};
+
+/** Simple add of one molecules into another.
+ * \param *mol destination molecule
+ * \param *srcmol source molecule
+ * \return true - merge successful, false - merge failed (probably due to non-existant indices
+ */
+bool MoleculeListClass::SimpleAdd(molecule *mol, molecule *srcmol)
+{
+  if (srcmol == NULL)
+    return false;
+
+  // put all molecules of src into mol
+  atom *Walker = srcmol->start;
+  atom *NextAtom = Walker->next;
+  while (NextAtom != srcmol->end) {
+    Walker = NextAtom;
+    NextAtom = Walker->next;
+    Walker = mol->AddCopyAtom(Walker);
+    Walker->father = Walker;
+  }
+
+  return true;
+};
+
+/** Simple merge of a given set of molecules into one.
+ * \param *mol destination molecule
+ * \param *src index of set of source molecule
+ * \param N number of source molecules
+ * \return true - merge successful, false - some merges failed (probably due to non-existant indices)
+ */
+bool MoleculeListClass::SimpleMultiMerge(molecule *mol, int *src, int N)
+{
+  bool status = true;
+  // check presence of all source molecules
+  for (int i=0;i<N;i++) {
+    molecule *srcmol = ReturnIndex(src[i]);
+    status = status && SimpleMerge(mol, srcmol);
+  }
+  return status;
+};
+
+/** Simple add of a given set of molecules into one.
+ * \param *mol destination molecule
+ * \param *src index of set of source molecule
+ * \param N number of source molecules
+ * \return true - merge successful, false - some merges failed (probably due to non-existant indices)
+ */
+bool MoleculeListClass::SimpleMultiAdd(molecule *mol, int *src, int N)
+{
+  bool status = true;
+  // check presence of all source molecules
+  for (int i=0;i<N;i++) {
+    molecule *srcmol = ReturnIndex(src[i]);
+    status = status && SimpleAdd(mol, srcmol);
+  }
+  return status;
+};
+
+/** Scatter merge of a given set of molecules into one.
+ * Scatter merge distributes the molecules in such a manner that they don't overlap.
+ * \param *mol destination molecule
+ * \param *src index of set of source molecule
+ * \param N number of source molecules
+ * \return true - merge successful, false - merge failed (probably due to non-existant indices
+ * \TODO find scatter center for each src molecule
+ */
+bool MoleculeListClass::ScatterMerge(molecule *mol, int *src, int N)
+{
+  // check presence of all source molecules
+  for (int i=0;i<N;i++) {
+    // get pointer to src molecule
+    molecule *srcmol = ReturnIndex(src[i]);
+    if (srcmol == NULL)
+      return false;
+  }
+  // adapt each Center
+  for (int i=0;i<N;i++) {
+    // get pointer to src molecule
+    molecule *srcmol = ReturnIndex(src[i]);
+    //srcmol->Center.Zero();
+    srcmol->Translate(&srcmol->Center);
+  }
+  // perform a simple multi merge
+  SimpleMultiMerge(mol, src, N);
+  return true;
+};
+
+/** Embedding merge of a given set of molecules into one.
+ * Embedding merge inserts one molecule into the other.
+ * \param *mol destination molecule
+ * \param *srcmol source molecule
+ * \return true - merge successful, false - merge failed (probably due to non-existant indices
+ * \TODO find embedding center
+ */
+bool MoleculeListClass::EmbedMerge(molecule *mol, molecule *srcmol)
+{
+  if (srcmol == NULL)
+    return false;
+
+  // calculate center for merge
+  //srcmol->Center.Zero();
+
+  // perform simple merge
+  SimpleMerge(mol, srcmol);
+  return true;
 };
 
@@ -127,8 +310,8 @@
 void MoleculeListClass::Output(ofstream *out)
 {
-	*out<< Verbose(1) << "MoleculeList: ";
-	for (int i=0;i<NumberOfMolecules;i++)
-		*out << ListOfMolecules[i] << "\t";
-	*out << endl;
+  *out << Verbose(1) << "MoleculeList: ";
+  for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++)
+    *out << *ListRunner << "\t";
+  *out << endl;
 };
 
@@ -142,181 +325,183 @@
 bool MoleculeListClass::AddHydrogenCorrection(ofstream *out, char *path)
 {
-	atom *Walker = NULL;
-	atom *Runner = NULL;
-	double ***FitConstant = NULL, **correction = NULL;
-	int a,b;
-	ofstream output;
-	ifstream input;
-	string line;
-	stringstream zeile;
-	double distance;
-	char ParsedLine[1023];
-	double tmp;
-	char *FragmentNumber = NULL;
-
-	cout << Verbose(1) << "Saving hydrogen saturation correction ... ";
-	// 0. parse in fit constant files that should have the same dimension as the final energy files
-	// 0a. find dimension of matrices with constants
-	line = path;
-	line.append("/");
-	line += FRAGMENTPREFIX;
-	line += "1";
-	line += FITCONSTANTSUFFIX;
-	input.open(line.c_str());
-	if (input == NULL) {
-		cerr << endl << "Unable to open " << line << ", is the directory correct?" << endl;
-		return false;
-	}
-	a=0;
-	b=-1; // we overcount by one 
-	while (!input.eof()) {
-		input.getline(ParsedLine, 1023);
-		zeile.str(ParsedLine);
-		int i=0;
-		while (!zeile.eof()) {
-			zeile >> distance;
-			i++;
-		}
-		if (i > a) 
-			a = i; 
-		b++;
-	}
-	cout << "I recognized " << a << " columns and " << b << " rows, ";
-	input.close();
-	
-	// 0b. allocate memory for constants
-	FitConstant = (double ***) Malloc(sizeof(double **)*3, "MoleculeListClass::AddHydrogenCorrection: ***FitConstant");
-	for (int k=0;k<3;k++) {
-		FitConstant[k] = (double **) Malloc(sizeof(double *)*a, "MoleculeListClass::AddHydrogenCorrection: **FitConstant[]");
-		for (int i=a;i--;) {
-			FitConstant[k][i] = (double *) Malloc(sizeof(double)*b, "MoleculeListClass::AddHydrogenCorrection: *FitConstant[][]");
-		}
-	}
-	// 0c. parse in constants
-	for (int i=0;i<3;i++) {
-		line = path;
-		line.append("/");
-		line += FRAGMENTPREFIX;
-		sprintf(ParsedLine, "%d", i+1);
-		line += ParsedLine;
-		line += FITCONSTANTSUFFIX;
-		input.open(line.c_str());
-		if (input == NULL) {
-			cerr << endl << "Unable to open " << line << ", is the directory correct?" << endl;
-			return false;
-		}
-		int k = 0,l;
-		while ((!input.eof()) && (k < b)) {
-			input.getline(ParsedLine, 1023);
-			//cout << "Current Line: " << ParsedLine << endl;
-			zeile.str(ParsedLine);
-			zeile.clear();
-			l = 0;
-			while ((!zeile.eof()) && (l < a)) {
-				zeile >> FitConstant[i][l][k];
-				//cout << FitConstant[i][l][k] << "\t";
-				l++;
-			}
-			//cout << endl;
-			k++;
-		}
-		input.close();
-	}
-	for(int k=0;k<3;k++) {
-		cout << "Constants " << k << ":" << endl;
-		for (int j=0;j<b;j++) {
-			for (int i=0;i<a;i++) {
-				cout << FitConstant[k][i][j] << "\t";
-			}
-			cout << endl;
-		}
-		cout << endl;
-	}
-	
-	// 0d. allocate final correction matrix
-	correction = (double **) Malloc(sizeof(double *)*a, "MoleculeListClass::AddHydrogenCorrection: **correction");
-	for (int i=a;i--;)
-		correction[i] = (double *) Malloc(sizeof(double)*b, "MoleculeListClass::AddHydrogenCorrection: *correction[]");
-				
-	// 1a. go through every molecule in the list
-	for(int i=NumberOfMolecules;i--;) {
-		// 1b. zero final correction matrix
-		for (int k=a;k--;)
-			for (int j=b;j--;)
-				correction[k][j] = 0.;
-		// 2. take every hydrogen that is a saturated one
-		Walker = ListOfMolecules[i]->start;
-		while (Walker->next != ListOfMolecules[i]->end) {
-			Walker = Walker->next;
-			//cout << Verbose(1) << "Walker: " << *Walker << " with first bond " << *ListOfMolecules[i]->ListOfBondsPerAtom[Walker->nr][0] << "." << endl;
-			if ((Walker->type->Z == 1) && ((Walker->father == NULL) || (Walker->father->type->Z != 1))) { // if it's a hydrogen
-				Runner = ListOfMolecules[i]->start;
-				while (Runner->next != ListOfMolecules[i]->end) {
-					Runner = Runner->next;
-					//cout << Verbose(2) << "Runner: " << *Runner << " with first bond " << *ListOfMolecules[i]->ListOfBondsPerAtom[Runner->nr][0] << "." << endl;
-					// 3. take every other hydrogen that is the not the first and not bound to same bonding partner
-					if ((Runner->type->Z == 1) && (Runner->nr > Walker->nr) && (ListOfMolecules[i]->ListOfBondsPerAtom[Runner->nr][0]->GetOtherAtom(Runner) != ListOfMolecules[i]->ListOfBondsPerAtom[Walker->nr][0]->GetOtherAtom(Walker))) {	// (hydrogens have only one bonding partner!)
-						// 4. evaluate the morse potential for each matrix component and add up
-						distance = Runner->x.Distance(&Walker->x);
-						//cout << "Fragment " << i << ": " << *Runner << "<= " << distance << "=>" << *Walker << ":" << endl;
-						for(int k=0;k<a;k++) {
-							for (int j=0;j<b;j++) {
-								switch(k) {
-								case 1:
-								case 7:
-								case 11:
-									tmp = pow(FitConstant[0][k][j] * ( 1. - exp(-FitConstant[1][k][j] * (distance - FitConstant[2][k][j]) ) ),2);
-									break;
-								default:
-									tmp = FitConstant[0][k][j] * pow( distance,	FitConstant[1][k][j]) + FitConstant[2][k][j];
-								};
-								correction[k][j] -= tmp;		// ground state is actually lower (disturbed by additional interaction)
-								//cout << tmp << "\t";
-							}
-							//cout << endl;
-						}
-						//cout << endl;
-					}
-				}
-			}
-		}
-		// 5. write final matrix to file
-		line = path;
-		line.append("/");
-		line += FRAGMENTPREFIX;
-		FragmentNumber = FixedDigitNumber(NumberOfMolecules, i);
-		line += FragmentNumber;
-		delete(FragmentNumber);
-		line += HCORRECTIONSUFFIX;
-		output.open(line.c_str());
-		output << "Time\t\tTotal\t\tKinetic\t\tNonLocal\tCorrelation\tExchange\tPseudo\t\tHartree\t\t-Gauss\t\tEwald\t\tIonKin\t\tETotal" << endl;
-		for (int j=0;j<b;j++) {
-			for(int i=0;i<a;i++)
-				output << correction[i][j] << "\t";
-			output << endl;
-		}
-		output.close();
-	}
-	line = path;
-	line.append("/");
-	line += HCORRECTIONSUFFIX;
-	output.open(line.c_str());
-	output << "Time\t\tTotal\t\tKinetic\t\tNonLocal\tCorrelation\tExchange\tPseudo\t\tHartree\t\t-Gauss\t\tEwald\t\tIonKin\t\tETotal" << endl;
-	for (int j=0;j<b;j++) {
-		for(int i=0;i<a;i++)
-			output << 0 << "\t";
-		output << endl;
-	}
-	output.close();
-	// 6. free memory of parsed matrices
-	FitConstant = (double ***) Malloc(sizeof(double **)*a, "MoleculeListClass::AddHydrogenCorrection: ***FitConstant");
-	for (int k=0;k<3;k++) {
-		FitConstant[k] = (double **) Malloc(sizeof(double *)*a, "MoleculeListClass::AddHydrogenCorrection: **FitConstant[]");
-		for (int i=a;i--;) {
-			FitConstant[k][i] = (double *) Malloc(sizeof(double)*b, "MoleculeListClass::AddHydrogenCorrection: *FitConstant[][]");
-		}
-	}
-	cout << "done." << endl;
-	return true;
+  atom *Walker = NULL;
+  atom *Runner = NULL;
+  double ***FitConstant = NULL, **correction = NULL;
+  int a, b;
+  ofstream output;
+  ifstream input;
+  string line;
+  stringstream zeile;
+  double distance;
+  char ParsedLine[1023];
+  double tmp;
+  char *FragmentNumber = NULL;
+
+  cout << Verbose(1) << "Saving hydrogen saturation correction ... ";
+  // 0. parse in fit constant files that should have the same dimension as the final energy files
+  // 0a. find dimension of matrices with constants
+  line = path;
+  line.append("/");
+  line += FRAGMENTPREFIX;
+  line += "1";
+  line += FITCONSTANTSUFFIX;
+  input.open(line.c_str());
+  if (input == NULL) {
+    cerr << endl << "Unable to open " << line << ", is the directory correct?"
+        << endl;
+    return false;
+  }
+  a = 0;
+  b = -1; // we overcount by one
+  while (!input.eof()) {
+    input.getline(ParsedLine, 1023);
+    zeile.str(ParsedLine);
+    int i = 0;
+    while (!zeile.eof()) {
+      zeile >> distance;
+      i++;
+    }
+    if (i > a)
+      a = i;
+    b++;
+  }
+  cout << "I recognized " << a << " columns and " << b << " rows, ";
+  input.close();
+
+  // 0b. allocate memory for constants
+  FitConstant = (double ***) Malloc(sizeof(double **) * 3, "MoleculeListClass::AddHydrogenCorrection: ***FitConstant");
+  for (int k = 0; k < 3; k++) {
+    FitConstant[k] = (double **) Malloc(sizeof(double *) * a, "MoleculeListClass::AddHydrogenCorrection: **FitConstant[]");
+    for (int i = a; i--;) {
+      FitConstant[k][i] = (double *) Malloc(sizeof(double) * b, "MoleculeListClass::AddHydrogenCorrection: *FitConstant[][]");
+    }
+  }
+  // 0c. parse in constants
+  for (int i = 0; i < 3; i++) {
+    line = path;
+    line.append("/");
+    line += FRAGMENTPREFIX;
+    sprintf(ParsedLine, "%d", i + 1);
+    line += ParsedLine;
+    line += FITCONSTANTSUFFIX;
+    input.open(line.c_str());
+    if (input == NULL) {
+      cerr << endl << "Unable to open " << line << ", is the directory correct?" << endl;
+      return false;
+    }
+    int k = 0, l;
+    while ((!input.eof()) && (k < b)) {
+      input.getline(ParsedLine, 1023);
+      //cout << "Current Line: " << ParsedLine << endl;
+      zeile.str(ParsedLine);
+      zeile.clear();
+      l = 0;
+      while ((!zeile.eof()) && (l < a)) {
+        zeile >> FitConstant[i][l][k];
+        //cout << FitConstant[i][l][k] << "\t";
+        l++;
+      }
+      //cout << endl;
+      k++;
+    }
+    input.close();
+  }
+  for (int k = 0; k < 3; k++) {
+    cout << "Constants " << k << ":" << endl;
+    for (int j = 0; j < b; j++) {
+      for (int i = 0; i < a; i++) {
+        cout << FitConstant[k][i][j] << "\t";
+      }
+      cout << endl;
+    }
+    cout << endl;
+  }
+
+  // 0d. allocate final correction matrix
+  correction = (double **) Malloc(sizeof(double *) * a, "MoleculeListClass::AddHydrogenCorrection: **correction");
+  for (int i = a; i--;)
+    correction[i] = (double *) Malloc(sizeof(double) * b, "MoleculeListClass::AddHydrogenCorrection: *correction[]");
+
+  // 1a. go through every molecule in the list
+  for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
+    // 1b. zero final correction matrix
+    for (int k = a; k--;)
+      for (int j = b; j--;)
+        correction[k][j] = 0.;
+    // 2. take every hydrogen that is a saturated one
+    Walker = (*ListRunner)->start;
+    while (Walker->next != (*ListRunner)->end) {
+      Walker = Walker->next;
+      //cout << Verbose(1) << "Walker: " << *Walker << " with first bond " << *(*Runner)->ListOfBondsPerAtom[Walker->nr][0] << "." << endl;
+      if ((Walker->type->Z == 1) && ((Walker->father == NULL)
+          || (Walker->father->type->Z != 1))) { // if it's a hydrogen
+        Runner = (*ListRunner)->start;
+        while (Runner->next != (*ListRunner)->end) {
+          Runner = Runner->next;
+          //cout << Verbose(2) << "Runner: " << *Runner << " with first bond " << *(*Runner)->ListOfBondsPerAtom[Runner->nr][0] << "." << endl;
+          // 3. take every other hydrogen that is the not the first and not bound to same bonding partner
+          if ((Runner->type->Z == 1) && (Runner->nr > Walker->nr) && ((*ListRunner)->ListOfBondsPerAtom[Runner->nr][0]->GetOtherAtom(Runner) != (*ListRunner)->ListOfBondsPerAtom[Walker->nr][0]->GetOtherAtom(Walker))) { // (hydrogens have only one bonding partner!)
+            // 4. evaluate the morse potential for each matrix component and add up
+            distance = Runner->x.Distance(&Walker->x);
+            //cout << "Fragment " << (*ListRunner)->name << ": " << *Runner << "<= " << distance << "=>" << *Walker << ":" << endl;
+            for (int k = 0; k < a; k++) {
+              for (int j = 0; j < b; j++) {
+                switch (k) {
+                  case 1:
+                  case 7:
+                  case 11:
+                    tmp = pow(FitConstant[0][k][j] * (1. - exp(-FitConstant[1][k][j] * (distance - FitConstant[2][k][j]))), 2);
+                    break;
+                  default:
+                    tmp = FitConstant[0][k][j] * pow(distance, FitConstant[1][k][j]) + FitConstant[2][k][j];
+                };
+                correction[k][j] -= tmp; // ground state is actually lower (disturbed by additional interaction)
+                //cout << tmp << "\t";
+              }
+              //cout << endl;
+            }
+            //cout << endl;
+          }
+        }
+      }
+    }
+    // 5. write final matrix to file
+    line = path;
+    line.append("/");
+    line += FRAGMENTPREFIX;
+    FragmentNumber = FixedDigitNumber(ListOfMolecules.size(), (*ListRunner)->IndexNr);
+    line += FragmentNumber;
+    delete (FragmentNumber);
+    line += HCORRECTIONSUFFIX;
+    output.open(line.c_str());
+    output << "Time\t\tTotal\t\tKinetic\t\tNonLocal\tCorrelation\tExchange\tPseudo\t\tHartree\t\t-Gauss\t\tEwald\t\tIonKin\t\tETotal" << endl;
+    for (int j = 0; j < b; j++) {
+      for (int i = 0; i < a; i++)
+        output << correction[i][j] << "\t";
+      output << endl;
+    }
+    output.close();
+  }
+  line = path;
+  line.append("/");
+  line += HCORRECTIONSUFFIX;
+  output.open(line.c_str());
+  output << "Time\t\tTotal\t\tKinetic\t\tNonLocal\tCorrelation\tExchange\tPseudo\t\tHartree\t\t-Gauss\t\tEwald\t\tIonKin\t\tETotal" << endl;
+  for (int j = 0; j < b; j++) {
+    for (int i = 0; i < a; i++)
+      output << 0 << "\t";
+    output << endl;
+  }
+  output.close();
+  // 6. free memory of parsed matrices
+  FitConstant = (double ***) Malloc(sizeof(double **) * a, "MoleculeListClass::AddHydrogenCorrection: ***FitConstant");
+  for (int k = 0; k < 3; k++) {
+    FitConstant[k] = (double **) Malloc(sizeof(double *) * a, "MoleculeListClass::AddHydrogenCorrection: **FitConstant[]");
+    for (int i = a; i--;) {
+      FitConstant[k][i] = (double *) Malloc(sizeof(double) * b, "MoleculeListClass::AddHydrogenCorrection: *FitConstant[][]");
+    }
+  }
+  cout << "done." << endl;
+  return true;
 };
 
@@ -327,49 +512,50 @@
  * \return true - file written successfully, false - writing failed
  */
-bool MoleculeListClass::StoreForcesFile(ofstream *out, char *path, int *SortIndex)
-{
-	bool status = true;
-	ofstream ForcesFile;
-	stringstream line;
-	atom *Walker = NULL;
-	element *runner = NULL;
-
-	// open file for the force factors
-	*out << Verbose(1) << "Saving	force factors ... ";
-	line << path << "/" << FRAGMENTPREFIX << FORCESFILE;
-	ForcesFile.open(line.str().c_str(), ios::out);
-	if (ForcesFile != NULL) {
-		//cout << Verbose(1) << "Final AtomicForcesList: ";
-		//output << prefix << "Forces" << endl;
-		for(int j=0;j<NumberOfMolecules;j++) {
-			//if (TEList[j] != 0) {
-			runner = ListOfMolecules[j]->elemente->start;
-			while (runner->next != ListOfMolecules[j]->elemente->end) { // go through every element
-				runner = runner->next;
-				if (ListOfMolecules[j]->ElementsInMolecule[runner->Z]) { // if this element got atoms
-					Walker = ListOfMolecules[j]->start;
-					while (Walker->next != ListOfMolecules[j]->end) { // go through every atom of this element
-						Walker = Walker->next;
-						if (Walker->type->Z == runner->Z) {
-							if ((Walker->GetTrueFather() != NULL) && (Walker->GetTrueFather() != Walker)) {// if there is a rea
-								//cout << "Walker is " << *Walker << " with true father " << *( Walker->GetTrueFather()) << ", it
-								ForcesFile <<	SortIndex[Walker->GetTrueFather()->nr] << "\t";
-								} else	// otherwise a -1 to indicate an added saturation hydrogen
-									ForcesFile << "-1\t";
-							}
-						}
-					}
-			}
-			ForcesFile << endl;
-		}
-		ForcesFile.close();
-		*out << Verbose(1) << "done." << endl;
-	} else {
-		status = false;
-		*out << Verbose(1) << "failed to open file " << line.str() << "." << endl;
-	}
-	ForcesFile.close();
-	
-	return status;
+bool MoleculeListClass::StoreForcesFile(ofstream *out, char *path,
+    int *SortIndex)
+{
+  bool status = true;
+  ofstream ForcesFile;
+  stringstream line;
+  atom *Walker = NULL;
+  element *runner = NULL;
+
+  // open file for the force factors
+  *out << Verbose(1) << "Saving	force factors ... ";
+  line << path << "/" << FRAGMENTPREFIX << FORCESFILE;
+  ForcesFile.open(line.str().c_str(), ios::out);
+  if (ForcesFile != NULL) {
+    //cout << Verbose(1) << "Final AtomicForcesList: ";
+    //output << prefix << "Forces" << endl;
+    for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
+      runner = (*ListRunner)->elemente->start;
+      while (runner->next != (*ListRunner)->elemente->end) { // go through every element
+        runner = runner->next;
+        if ((*ListRunner)->ElementsInMolecule[runner->Z]) { // if this element got atoms
+          Walker = (*ListRunner)->start;
+          while (Walker->next != (*ListRunner)->end) { // go through every atom of this element
+            Walker = Walker->next;
+            if (Walker->type->Z == runner->Z) {
+              if ((Walker->GetTrueFather() != NULL) && (Walker->GetTrueFather() != Walker)) {// if there is a rea
+                //cout << "Walker is " << *Walker << " with true father " << *( Walker->GetTrueFather()) << ", it
+                ForcesFile << SortIndex[Walker->GetTrueFather()->nr] << "\t";
+              } else
+                // otherwise a -1 to indicate an added saturation hydrogen
+                ForcesFile << "-1\t";
+            }
+          }
+        }
+      }
+      ForcesFile << endl;
+    }
+    ForcesFile.close();
+    *out << Verbose(1) << "done." << endl;
+  } else {
+    status = false;
+    *out << Verbose(1) << "failed to open file " << line.str() << "." << endl;
+  }
+  ForcesFile.close();
+
+  return status;
 };
 
@@ -380,106 +566,118 @@
  * \return true - success (each file was written), false - something went wrong.
  */
-bool MoleculeListClass::OutputConfigForListOfFragments(ofstream *out, config *configuration, int *SortIndex)
-{
-	ofstream outputFragment;
-	char FragmentName[MAXSTRINGSIZE];
-	char PathBackup[MAXSTRINGSIZE];
-	bool result = true;
-	bool intermediateResult = true;
-	atom *Walker = NULL;
-	Vector BoxDimension;
-	char *FragmentNumber = NULL;
-	char *path = NULL;
-	int FragmentCounter = 0;
-	ofstream output;
-	
-	// store the fragments as config and as xyz
-	for(int i=0;i<NumberOfMolecules;i++) {
-		// save default path as it is changed for each fragment
-		path = configuration->GetDefaultPath();
-		if (path != NULL)
-			strcpy(PathBackup, path);
-		else
-			cerr << "OutputConfigForListOfFragments: NULL default path obtained from config!" << endl;
-
-		// correct periodic
-		ListOfMolecules[i]->ScanForPeriodicCorrection(out);
-
-		// output xyz file
-		FragmentNumber = FixedDigitNumber(NumberOfMolecules, FragmentCounter++);
-		sprintf(FragmentName, "%s/%s%s.conf.xyz", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
-		outputFragment.open(FragmentName, ios::out);
-		*out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter-1 << " as XYZ ...";
-		if ((intermediateResult = ListOfMolecules[i]->OutputXYZ(&outputFragment)))
-			*out << " done." << endl;
-		else
-			*out << " failed." << endl;
-		result = result && intermediateResult;
-		outputFragment.close();
-		outputFragment.clear();
-
-		// list atoms in fragment for debugging
-		*out << Verbose(2) << "Contained atoms: ";
-		Walker = ListOfMolecules[i]->start;
-		while (Walker->next != ListOfMolecules[i]->end) {
-			Walker = Walker->next;
-			*out << Walker->Name << " ";
-		}
-		*out << endl;
-		
-		// center on edge
-		ListOfMolecules[i]->CenterEdge(out, &BoxDimension);
-		ListOfMolecules[i]->SetBoxDimension(&BoxDimension);	// update Box of atoms by boundary
-		int j = -1;
-		for (int k=0;k<NDIM;k++) {
-			j += k+1;
-			BoxDimension.x[k] = 2.5* (configuration->GetIsAngstroem() ? 1. : 1./AtomicLengthToAngstroem);
-			ListOfMolecules[i]->cell_size[j] += BoxDimension.x[k]*2.;
-		}
-		ListOfMolecules[i]->Translate(&BoxDimension);
-
-		// also calculate necessary orbitals
-		ListOfMolecules[i]->CountElements();	// this is a bugfix, atoms should should actually be added correctly to this fragment 
-		ListOfMolecules[i]->CalculateOrbitals(*configuration);
-		
-		// change path in config
-		//strcpy(PathBackup, configuration->configpath);
-		sprintf(FragmentName, "%s/%s%s/", PathBackup, FRAGMENTPREFIX, FragmentNumber);
-		configuration->SetDefaultPath(FragmentName);
-		
-		// and save as config
-		sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
-		*out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter-1 << " as config ...";
-		if ((intermediateResult = configuration->Save(FragmentName, ListOfMolecules[i]->elemente, ListOfMolecules[i])))
-			*out << " done." << endl;
-		else
-			*out << " failed." << endl;
-		result = result && intermediateResult;
-
-		// restore old config
-		configuration->SetDefaultPath(PathBackup);
-
-
-		// and save as mpqc input file
-		sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
-		*out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter-1 << " as mpqc input ...";
-		if ((intermediateResult = configuration->SaveMPQC(FragmentName, ListOfMolecules[i])))
-			*out << " done." << endl;
-		else
-			*out << " failed." << endl;
-			 
-		result = result && intermediateResult;
-		//outputFragment.close();
-		//outputFragment.clear();
-		delete(FragmentNumber);
-		//Free((void **)&FragmentNumber, "MoleculeListClass::OutputConfigForListOfFragments: *FragmentNumber");
-	}
-	cout << " done." << endl;
-	
-	// printing final number
-	*out << "Final number of fragments: " << FragmentCounter << "." << endl;
-			
-	return result;
-};
+bool MoleculeListClass::OutputConfigForListOfFragments(ofstream *out,
+    config *configuration, int *SortIndex)
+{
+  ofstream outputFragment;
+  char FragmentName[MAXSTRINGSIZE];
+  char PathBackup[MAXSTRINGSIZE];
+  bool result = true;
+  bool intermediateResult = true;
+  atom *Walker = NULL;
+  Vector BoxDimension;
+  char *FragmentNumber = NULL;
+  char *path = NULL;
+  int FragmentCounter = 0;
+  ofstream output;
+
+  // store the fragments as config and as xyz
+  for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
+    // save default path as it is changed for each fragment
+    path = configuration->GetDefaultPath();
+    if (path != NULL)
+      strcpy(PathBackup, path);
+    else
+      cerr << "OutputConfigForListOfFragments: NULL default path obtained from config!" << endl;
+
+    // correct periodic
+    (*ListRunner)->ScanForPeriodicCorrection(out);
+
+    // output xyz file
+    FragmentNumber = FixedDigitNumber(ListOfMolecules.size(), FragmentCounter++);
+    sprintf(FragmentName, "%s/%s%s.conf.xyz", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
+    outputFragment.open(FragmentName, ios::out);
+    *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter - 1 << " as XYZ ...";
+    if ((intermediateResult = (*ListRunner)->OutputXYZ(&outputFragment)))
+      *out << " done." << endl;
+    else
+      *out << " failed." << endl;
+    result = result && intermediateResult;
+    outputFragment.close();
+    outputFragment.clear();
+
+    // list atoms in fragment for debugging
+    *out << Verbose(2) << "Contained atoms: ";
+    Walker = (*ListRunner)->start;
+    while (Walker->next != (*ListRunner)->end) {
+      Walker = Walker->next;
+      *out << Walker->Name << " ";
+    }
+    *out << endl;
+
+    // center on edge
+    (*ListRunner)->CenterEdge(out, &BoxDimension);
+    (*ListRunner)->SetBoxDimension(&BoxDimension); // update Box of atoms by boundary
+    int j = -1;
+    for (int k = 0; k < NDIM; k++) {
+      j += k + 1;
+      BoxDimension.x[k] = 2.5 * (configuration->GetIsAngstroem() ? 1. : 1. / AtomicLengthToAngstroem);
+      (*ListRunner)->cell_size[j] += BoxDimension.x[k] * 2.;
+    }
+    (*ListRunner)->Translate(&BoxDimension);
+
+    // also calculate necessary orbitals
+    (*ListRunner)->CountElements(); // this is a bugfix, atoms should shoulds actually be added correctly to this fragment
+    (*ListRunner)->CalculateOrbitals(*configuration);
+
+    // change path in config
+    //strcpy(PathBackup, configuration->configpath);
+    sprintf(FragmentName, "%s/%s%s/", PathBackup, FRAGMENTPREFIX, FragmentNumber);
+    configuration->SetDefaultPath(FragmentName);
+
+    // and save as config
+    sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
+    *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter - 1 << " as config ...";
+    if ((intermediateResult = configuration->Save(FragmentName, (*ListRunner)->elemente, (*ListRunner))))
+      *out << " done." << endl;
+    else
+      *out << " failed." << endl;
+    result = result && intermediateResult;
+
+    // restore old config
+    configuration->SetDefaultPath(PathBackup);
+
+    // and save as mpqc input file
+    sprintf(FragmentName, "%s/%s%s.conf", configuration->configpath, FRAGMENTPREFIX, FragmentNumber);
+    *out << Verbose(2) << "Saving bond fragment No. " << FragmentNumber << "/" << FragmentCounter - 1 << " as mpqc input ...";
+    if ((intermediateResult = configuration->SaveMPQC(FragmentName, (*ListRunner))))
+      *out << " done." << endl;
+    else
+      *out << " failed." << endl;
+
+    result = result && intermediateResult;
+    //outputFragment.close();
+    //outputFragment.clear();
+    delete (FragmentNumber);
+    //Free((void **)&FragmentNumber, "MoleculeListClass::OutputConfigForListOfFragments: *FragmentNumber");
+  }
+  cout << " done." << endl;
+
+  // printing final number
+  *out << "Final number of fragments: " << FragmentCounter << "." << endl;
+
+  return result;
+};
+
+/** Counts the number of molecules with the molecule::ActiveFlag set.
+ * \return number of molecules with ActiveFlag set to true.
+ */
+int MoleculeListClass::NumberOfActiveMolecules()
+{
+  int count = 0;
+  for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++)
+    count += ((*ListRunner)->ActiveFlag ? 1 : 0);
+  return count;
+};
+
 
 /******************************************* Class MoleculeLeafClass ************************************************/
@@ -492,18 +690,18 @@
 MoleculeLeafClass::MoleculeLeafClass(MoleculeLeafClass *PreviousLeaf = NULL)
 {
-//	if (Up != NULL)
-//		if (Up->DownLeaf == NULL) // are we the first down leaf for the upper leaf?
-//			Up->DownLeaf = this;
-//	UpLeaf = Up;
-//	DownLeaf = NULL;
-	Leaf = NULL;
-	previous = PreviousLeaf;
-	if (previous != NULL) {
-		MoleculeLeafClass *Walker = previous->next;
-		previous->next = this;
-		next = Walker;
-	} else {
-		next = NULL;
-	}
+  //	if (Up != NULL)
+  //		if (Up->DownLeaf == NULL) // are we the first down leaf for the upper leaf?
+  //			Up->DownLeaf = this;
+  //	UpLeaf = Up;
+  //	DownLeaf = NULL;
+  Leaf = NULL;
+  previous = PreviousLeaf;
+  if (previous != NULL) {
+    MoleculeLeafClass *Walker = previous->next;
+    previous->next = this;
+    next = Walker;
+  } else {
+    next = NULL;
+  }
 };
 
@@ -512,33 +710,33 @@
 MoleculeLeafClass::~MoleculeLeafClass()
 {
-//	if (DownLeaf != NULL) {// drop leaves further down
-//		MoleculeLeafClass *Walker = DownLeaf;
-//		MoleculeLeafClass *Next;
-//		do {
-//			Next = Walker->NextLeaf;
-//			delete(Walker);
-//			Walker = Next;
-//		} while (Walker != NULL);
-//		// Last Walker sets DownLeaf automatically to NULL
-//	}
-	// remove the leaf itself
-	if (Leaf != NULL) {
-		delete(Leaf);
-		Leaf = NULL;
-	}
-	// remove this Leaf from level list
-	if (previous != NULL)	 
-		previous->next = next;
-//	} else { // we are first in list (connects to UpLeaf->DownLeaf)
-//		if ((NextLeaf != NULL) && (NextLeaf->UpLeaf == NULL))
-//			NextLeaf->UpLeaf = UpLeaf;	// either null as we are top level or the upleaf of the first node
-//		if (UpLeaf != NULL)
-//			UpLeaf->DownLeaf = NextLeaf;	// either null as we are only leaf or NextLeaf if we are just the first
-//	}
-//	UpLeaf = NULL;
-	if (next != NULL) // are we last in list
-		next->previous = previous;
-	next = NULL;
-	previous = NULL;
+  //	if (DownLeaf != NULL) {// drop leaves further down
+  //		MoleculeLeafClass *Walker = DownLeaf;
+  //		MoleculeLeafClass *Next;
+  //		do {
+  //			Next = Walker->NextLeaf;
+  //			delete(Walker);
+  //			Walker = Next;
+  //		} while (Walker != NULL);
+  //		// Last Walker sets DownLeaf automatically to NULL
+  //	}
+  // remove the leaf itself
+  if (Leaf != NULL) {
+    delete (Leaf);
+    Leaf = NULL;
+  }
+  // remove this Leaf from level list
+  if (previous != NULL)
+    previous->next = next;
+  //	} else { // we are first in list (connects to UpLeaf->DownLeaf)
+  //		if ((NextLeaf != NULL) && (NextLeaf->UpLeaf == NULL))
+  //			NextLeaf->UpLeaf = UpLeaf;	// either null as we are top level or the upleaf of the first node
+  //		if (UpLeaf != NULL)
+  //			UpLeaf->DownLeaf = NextLeaf;	// either null as we are only leaf or NextLeaf if we are just the first
+  //	}
+  //	UpLeaf = NULL;
+  if (next != NULL) // are we last in list
+    next->previous = previous;
+  next = NULL;
+  previous = NULL;
 };
 
@@ -550,5 +748,5 @@
 bool MoleculeLeafClass::AddLeaf(molecule *ptr, MoleculeLeafClass *Previous)
 {
-	return false;
+  return false;
 };
 
@@ -564,50 +762,51 @@
 bool MoleculeLeafClass::FillBondStructureFromReference(ofstream *out, molecule *reference, int &FragmentCounter, atom ***&ListOfLocalAtoms, bool FreeList)
 {
-	atom *Walker = NULL, *OtherWalker = NULL;
-	bond *Binder = NULL;
-	bool status = true;
-	int AtomNo;
-
-	*out << Verbose(1) << "Begin of FillBondStructureFromReference." << endl;
-	// fill ListOfLocalAtoms if NULL was given
-	if (!FillListOfLocalAtoms(out, ListOfLocalAtoms, FragmentCounter, reference->AtomCount, FreeList)) {
-		*out << Verbose(1) << "Filling of ListOfLocalAtoms failed." << endl;
-		return false;
-	}
-	
-	if (status) {
-		*out << Verbose(1) << "Creating adjacency list for subgraph " << this << "." << endl;
-		Walker = Leaf->start;
-		while (Walker->next != Leaf->end) {
-			Walker = Walker->next;
-			AtomNo = Walker->GetTrueFather()->nr;	// global id of the current walker
-			for(int i=0;i<reference->NumberOfBondsPerAtom[AtomNo];i++) { // go through father's bonds and copy them all
-				Binder = reference->ListOfBondsPerAtom[AtomNo][i];
-				OtherWalker = ListOfLocalAtoms[FragmentCounter][Binder->GetOtherAtom(Walker->GetTrueFather())->nr];		// local copy of current bond partner of walker
-				if (OtherWalker != NULL) {
-					if (OtherWalker->nr > Walker->nr)
-					Leaf->AddBond(Walker, OtherWalker, Binder->BondDegree);
-				} else {
-					*out << Verbose(1) << "OtherWalker = ListOfLocalAtoms[" << FragmentCounter << "][" << Binder->GetOtherAtom(Walker->GetTrueFather())->nr << "] is NULL!" << endl;
-					status = false;
-				}
-			}
-		}
-		Leaf->CreateListOfBondsPerAtom(out);
-		FragmentCounter++;
-		if (next != NULL)
-			status = next->FillBondStructureFromReference(out, reference, FragmentCounter, ListOfLocalAtoms);
-		FragmentCounter--;
-	}
-	
-	if ((FreeList) && (ListOfLocalAtoms != NULL)) {
-		// free the index lookup list
-		Free((void **)&ListOfLocalAtoms[FragmentCounter], "MoleculeLeafClass::FillBondStructureFromReference - **ListOfLocalAtoms[]");
-		if (FragmentCounter == 0) // first fragments frees the initial pointer to list
-			Free((void **)&ListOfLocalAtoms, "MoleculeLeafClass::FillBondStructureFromReference - ***ListOfLocalAtoms");
-	}
-	FragmentCounter--;
-	*out << Verbose(1) << "End of FillBondStructureFromReference." << endl;
-	return status;
+  atom *Walker = NULL, *OtherWalker = NULL;
+  bond *Binder = NULL;
+  bool status = true;
+  int AtomNo;
+
+  *out << Verbose(1) << "Begin of FillBondStructureFromReference." << endl;
+  // fill ListOfLocalAtoms if NULL was given
+  if (!FillListOfLocalAtoms(out, ListOfLocalAtoms, FragmentCounter, reference->AtomCount, FreeList)) {
+    *out << Verbose(1) << "Filling of ListOfLocalAtoms failed." << endl;
+    return false;
+  }
+
+  if (status) {
+    *out << Verbose(1) << "Creating adjacency list for subgraph " << this
+        << "." << endl;
+    Walker = Leaf->start;
+    while (Walker->next != Leaf->end) {
+      Walker = Walker->next;
+      AtomNo = Walker->GetTrueFather()->nr; // global id of the current walker
+      for (int i = 0; i < reference->NumberOfBondsPerAtom[AtomNo]; i++) { // go through father's bonds and copy them all
+        Binder = reference->ListOfBondsPerAtom[AtomNo][i];
+        OtherWalker = ListOfLocalAtoms[FragmentCounter][Binder->GetOtherAtom(Walker->GetTrueFather())->nr]; // local copy of current bond partner of walker
+        if (OtherWalker != NULL) {
+          if (OtherWalker->nr > Walker->nr)
+            Leaf->AddBond(Walker, OtherWalker, Binder->BondDegree);
+        } else {
+          *out << Verbose(1) << "OtherWalker = ListOfLocalAtoms[" << FragmentCounter << "][" << Binder->GetOtherAtom(Walker->GetTrueFather())->nr << "] is NULL!" << endl;
+          status = false;
+        }
+      }
+    }
+    Leaf->CreateListOfBondsPerAtom(out);
+    FragmentCounter++;
+    if (next != NULL)
+      status = next->FillBondStructureFromReference(out, reference, FragmentCounter, ListOfLocalAtoms);
+    FragmentCounter--;
+  }
+
+  if ((FreeList) && (ListOfLocalAtoms != NULL)) {
+    // free the index lookup list
+    Free((void **) &ListOfLocalAtoms[FragmentCounter], "MoleculeLeafClass::FillBondStructureFromReference - **ListOfLocalAtoms[]");
+    if (FragmentCounter == 0) // first fragments frees the initial pointer to list
+      Free((void **) &ListOfLocalAtoms, "MoleculeLeafClass::FillBondStructureFromReference - ***ListOfLocalAtoms");
+  }
+  FragmentCounter--;
+  *out << Verbose(1) << "End of FillBondStructureFromReference." << endl;
+  return status;
 };
 
@@ -620,34 +819,35 @@
  * \return true - stack is non-empty, fragmentation necessary, false - stack is empty, no more sites to update
  */
-bool MoleculeLeafClass::FillRootStackForSubgraphs(ofstream *out, KeyStack *&RootStack, bool *AtomMask, int &FragmentCounter)
-{
-	atom *Walker = NULL, *Father = NULL;
-
-	if (RootStack != NULL) {
-		// find first root candidates 
-		if (&(RootStack[FragmentCounter]) != NULL) {
-			RootStack[FragmentCounter].clear();	
-			Walker = Leaf->start;
-			while (Walker->next != Leaf->end) { // go through all (non-hydrogen) atoms
-				Walker = Walker->next;
-				Father = Walker->GetTrueFather();
-				if (AtomMask[Father->nr]) // apply mask
-		#ifdef ADDHYDROGEN
-					if (Walker->type->Z != 1) // skip hydrogen
-		#endif
-						RootStack[FragmentCounter].push_front(Walker->nr);
-			}
-			if (next != NULL)
-				next->FillRootStackForSubgraphs(out, RootStack, AtomMask, ++FragmentCounter);
-		}	else {
-			*out << Verbose(1) << "Rootstack[" << FragmentCounter	<< "] is NULL." << endl;
-			return false;
-		}
-		FragmentCounter--;
-		return true;
-	} else {
-		*out << Verbose(1) << "Rootstack is NULL." << endl;
-		return false;
-	}
+bool MoleculeLeafClass::FillRootStackForSubgraphs(ofstream *out,
+    KeyStack *&RootStack, bool *AtomMask, int &FragmentCounter)
+{
+  atom *Walker = NULL, *Father = NULL;
+
+  if (RootStack != NULL) {
+    // find first root candidates
+    if (&(RootStack[FragmentCounter]) != NULL) {
+      RootStack[FragmentCounter].clear();
+      Walker = Leaf->start;
+      while (Walker->next != Leaf->end) { // go through all (non-hydrogen) atoms
+        Walker = Walker->next;
+        Father = Walker->GetTrueFather();
+        if (AtomMask[Father->nr]) // apply mask
+#ifdef ADDHYDROGEN
+          if (Walker->type->Z != 1) // skip hydrogen
+#endif
+          RootStack[FragmentCounter].push_front(Walker->nr);
+      }
+      if (next != NULL)
+        next->FillRootStackForSubgraphs(out, RootStack, AtomMask, ++FragmentCounter);
+    } else {
+      *out << Verbose(1) << "Rootstack[" << FragmentCounter << "] is NULL." << endl;
+      return false;
+    }
+    FragmentCounter--;
+    return true;
+  } else {
+    *out << Verbose(1) << "Rootstack is NULL." << endl;
+    return false;
+  }
 };
 
@@ -662,24 +862,24 @@
 bool MoleculeLeafClass::FillListOfLocalAtoms(ofstream *out, atom ***&ListOfLocalAtoms, const int FragmentCounter, const int GlobalAtomCount, bool &FreeList)
 {
-	bool status = true;
-	
-	int Counter = Count();
-	if (ListOfLocalAtoms == NULL) { // allocated initial pointer
-		// allocate and set each field to NULL
-		ListOfLocalAtoms = (atom ***) Malloc(sizeof(atom **)*Counter, "MoleculeLeafClass::FillBondStructureFromReference - ***ListOfLocalAtoms");
-		if (ListOfLocalAtoms != NULL) {
-			for (int i=Counter;i--;)
-				ListOfLocalAtoms[i] = NULL;
-			FreeList = FreeList && true;
-		} else
-			status = false;
-	}
-	
-	if ((ListOfLocalAtoms != NULL) && (ListOfLocalAtoms[FragmentCounter] == NULL)) { // allocate and fill list of this fragment/subgraph
-		status = status && CreateFatherLookupTable(out, Leaf->start, Leaf->end, ListOfLocalAtoms[FragmentCounter], GlobalAtomCount);
-		FreeList = FreeList && true;
-	}
-	
-	return status;
+  bool status = true;
+
+  int Counter = Count();
+  if (ListOfLocalAtoms == NULL) { // allocated initial pointer
+    // allocate and set each field to NULL
+    ListOfLocalAtoms = (atom ***) Malloc(sizeof(atom **) * Counter, "MoleculeLeafClass::FillBondStructureFromReference - ***ListOfLocalAtoms");
+    if (ListOfLocalAtoms != NULL) {
+      for (int i = Counter; i--;)
+        ListOfLocalAtoms[i] = NULL;
+      FreeList = FreeList && true;
+    } else
+      status = false;
+  }
+
+  if ((ListOfLocalAtoms != NULL) && (ListOfLocalAtoms[FragmentCounter] == NULL)) { // allocate and fill list of this fragment/subgraph
+    status = status && CreateFatherLookupTable(out, Leaf->start, Leaf->end, ListOfLocalAtoms[FragmentCounter], GlobalAtomCount);
+    FreeList = FreeList && true;
+  }
+
+  return status;
 };
 
@@ -694,61 +894,63 @@
  * \retuen true - success, false - failure
  */
-bool MoleculeLeafClass::AssignKeySetsToFragment(ofstream *out, molecule *reference, Graph *KeySetList, atom ***&ListOfLocalAtoms, Graph **&FragmentList, int &FragmentCounter, bool FreeList)
-{
-	bool status = true;
-	int KeySetCounter = 0;
-	
-	*out << Verbose(1) << "Begin of AssignKeySetsToFragment." << endl;
-	// fill ListOfLocalAtoms if NULL was given
-	if (!FillListOfLocalAtoms(out, ListOfLocalAtoms, FragmentCounter, reference->AtomCount, FreeList)) {
-		*out << Verbose(1) << "Filling of ListOfLocalAtoms failed." << endl;
-		return false;
-	}
-
-	// allocate fragment list
-	if (FragmentList == NULL) {
-		KeySetCounter = Count();
-		FragmentList = (Graph **) Malloc(sizeof(Graph *)*KeySetCounter, "MoleculeLeafClass::AssignKeySetsToFragment - **FragmentList");
-		for(int i=KeySetCounter;i--;)
-			FragmentList[i] = NULL;
-		KeySetCounter = 0;
-	}
-	
-	if ((KeySetList != NULL) && (KeySetList->size() != 0)) { // if there are some scanned keysets at all
-		// assign scanned keysets
-		if (FragmentList[FragmentCounter] == NULL)
-			FragmentList[FragmentCounter] = new Graph;
-		KeySet *TempSet = new KeySet;
-		for(Graph::iterator runner = KeySetList->begin();runner != KeySetList->end(); runner++) { // key sets contain global numbers!
-			if ( ListOfLocalAtoms[FragmentCounter][reference->FindAtom(*((*runner).first.begin()))->nr] != NULL) {// as we may assume that that bond structure is unchanged, we only test the first key in each set
-				// translate keyset to local numbers
-				for(KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++)
-					TempSet->insert(ListOfLocalAtoms[FragmentCounter][reference->FindAtom(*sprinter)->nr]->nr);
-				// insert into FragmentList
-				FragmentList[FragmentCounter]->insert(GraphPair (*TempSet, pair<int,double>(KeySetCounter++, (*runner).second.second)));
-			}
-			TempSet->clear();
-		}
-		delete(TempSet);
-		if (KeySetCounter == 0) {// if there are no keysets, delete the list
-			*out << Verbose(1) << "KeySetCounter is zero, deleting FragmentList." << endl;
-			delete(FragmentList[FragmentCounter]);
-		} else
-			*out << Verbose(1) << KeySetCounter << " keysets were assigned to subgraph " << FragmentCounter << "." << endl;
-		FragmentCounter++;
-		if (next != NULL)
-			next->AssignKeySetsToFragment(out, reference, KeySetList, ListOfLocalAtoms, FragmentList, FragmentCounter, FreeList);
-		FragmentCounter--;
-	} else
-		*out << Verbose(1) << "KeySetList is NULL or empty." << endl;
-	
-	if ((FreeList) && (ListOfLocalAtoms != NULL)) {
-		// free the index lookup list
-		Free((void **)&ListOfLocalAtoms[FragmentCounter], "MoleculeLeafClass::AssignKeySetsToFragment - **ListOfLocalAtoms[]");
-		if (FragmentCounter == 0) // first fragments frees the initial pointer to list
-			Free((void **)&ListOfLocalAtoms, "MoleculeLeafClass::AssignKeySetsToFragment - ***ListOfLocalAtoms");
-	}
-	*out << Verbose(1) << "End of AssignKeySetsToFragment." << endl;
-	return status;
+bool MoleculeLeafClass::AssignKeySetsToFragment(ofstream *out,
+    molecule *reference, Graph *KeySetList, atom ***&ListOfLocalAtoms,
+    Graph **&FragmentList, int &FragmentCounter, bool FreeList)
+{
+  bool status = true;
+  int KeySetCounter = 0;
+
+  *out << Verbose(1) << "Begin of AssignKeySetsToFragment." << endl;
+  // fill ListOfLocalAtoms if NULL was given
+  if (!FillListOfLocalAtoms(out, ListOfLocalAtoms, FragmentCounter, reference->AtomCount, FreeList)) {
+    *out << Verbose(1) << "Filling of ListOfLocalAtoms failed." << endl;
+    return false;
+  }
+
+  // allocate fragment list
+  if (FragmentList == NULL) {
+    KeySetCounter = Count();
+    FragmentList = (Graph **) Malloc(sizeof(Graph *) * KeySetCounter, "MoleculeLeafClass::AssignKeySetsToFragment - **FragmentList");
+    for (int i = KeySetCounter; i--;)
+      FragmentList[i] = NULL;
+    KeySetCounter = 0;
+  }
+
+  if ((KeySetList != NULL) && (KeySetList->size() != 0)) { // if there are some scanned keysets at all
+    // assign scanned keysets
+    if (FragmentList[FragmentCounter] == NULL)
+      FragmentList[FragmentCounter] = new Graph;
+    KeySet *TempSet = new KeySet;
+    for (Graph::iterator runner = KeySetList->begin(); runner != KeySetList->end(); runner++) { // key sets contain global numbers!
+      if (ListOfLocalAtoms[FragmentCounter][reference->FindAtom(*((*runner).first.begin()))->nr] != NULL) {// as we may assume that that bond structure is unchanged, we only test the first key in each set
+        // translate keyset to local numbers
+        for (KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++)
+          TempSet->insert(ListOfLocalAtoms[FragmentCounter][reference->FindAtom(*sprinter)->nr]->nr);
+        // insert into FragmentList
+        FragmentList[FragmentCounter]->insert(GraphPair(*TempSet, pair<int, double> (KeySetCounter++, (*runner).second.second)));
+      }
+      TempSet->clear();
+    }
+    delete (TempSet);
+    if (KeySetCounter == 0) {// if there are no keysets, delete the list
+      *out << Verbose(1) << "KeySetCounter is zero, deleting FragmentList." << endl;
+      delete (FragmentList[FragmentCounter]);
+    } else
+      *out << Verbose(1) << KeySetCounter << " keysets were assigned to subgraph " << FragmentCounter << "." << endl;
+    FragmentCounter++;
+    if (next != NULL)
+      next->AssignKeySetsToFragment(out, reference, KeySetList, ListOfLocalAtoms, FragmentList, FragmentCounter, FreeList);
+    FragmentCounter--;
+  } else
+    *out << Verbose(1) << "KeySetList is NULL or empty." << endl;
+
+  if ((FreeList) && (ListOfLocalAtoms != NULL)) {
+    // free the index lookup list
+    Free((void **) &ListOfLocalAtoms[FragmentCounter], "MoleculeLeafClass::AssignKeySetsToFragment - **ListOfLocalAtoms[]");
+    if (FragmentCounter == 0) // first fragments frees the initial pointer to list
+      Free((void **) &ListOfLocalAtoms, "MoleculeLeafClass::AssignKeySetsToFragment - ***ListOfLocalAtoms");
+  }
+  *out << Verbose(1) << "End of AssignKeySetsToFragment." << endl;
+  return status;
 };
 
@@ -759,24 +961,26 @@
  * \param &TotalNumberOfKeySets global key set counter
  * \param &TotalGraph Graph to be filled with global numbers
- */ 
-void MoleculeLeafClass::TranslateIndicesToGlobalIDs(ofstream *out, Graph **FragmentList, int &FragmentCounter, int &TotalNumberOfKeySets, Graph &TotalGraph)
-{
-	*out << Verbose(1) << "Begin of TranslateIndicesToGlobalIDs." << endl;
-	KeySet *TempSet = new KeySet;
-	if (FragmentList[FragmentCounter] != NULL) {
-		for(Graph::iterator runner = FragmentList[FragmentCounter]->begin(); runner != FragmentList[FragmentCounter]->end(); runner++) {
-			for(KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++)
-				TempSet->insert((Leaf->FindAtom(*sprinter))->GetTrueFather()->nr);
-			TotalGraph.insert(GraphPair(*TempSet, pair<int,double>(TotalNumberOfKeySets++, (*runner).second.second)));
-			TempSet->clear();
-		}
-		delete(TempSet);
-	} else {
-		*out << Verbose(1) << "FragmentList is NULL." << endl;
-	}
-	if (next != NULL)
-		next->TranslateIndicesToGlobalIDs(out, FragmentList, ++FragmentCounter, TotalNumberOfKeySets, TotalGraph);
-	FragmentCounter--;
-	*out << Verbose(1) << "End of TranslateIndicesToGlobalIDs." << endl;
+ */
+void MoleculeLeafClass::TranslateIndicesToGlobalIDs(ofstream *out,
+    Graph **FragmentList, int &FragmentCounter, int &TotalNumberOfKeySets,
+    Graph &TotalGraph)
+{
+  *out << Verbose(1) << "Begin of TranslateIndicesToGlobalIDs." << endl;
+  KeySet *TempSet = new KeySet;
+  if (FragmentList[FragmentCounter] != NULL) {
+    for (Graph::iterator runner = FragmentList[FragmentCounter]->begin(); runner != FragmentList[FragmentCounter]->end(); runner++) {
+      for (KeySet::iterator sprinter = (*runner).first.begin(); sprinter != (*runner).first.end(); sprinter++)
+        TempSet->insert((Leaf->FindAtom(*sprinter))->GetTrueFather()->nr);
+      TotalGraph.insert(GraphPair(*TempSet, pair<int, double> (TotalNumberOfKeySets++, (*runner).second.second)));
+      TempSet->clear();
+    }
+    delete (TempSet);
+  } else {
+    *out << Verbose(1) << "FragmentList is NULL." << endl;
+  }
+  if (next != NULL)
+    next->TranslateIndicesToGlobalIDs(out, FragmentList, ++FragmentCounter, TotalNumberOfKeySets, TotalGraph);
+  FragmentCounter--;
+  *out << Verbose(1) << "End of TranslateIndicesToGlobalIDs." << endl;
 };
 
@@ -786,7 +990,8 @@
 int MoleculeLeafClass::Count() const
 {
-	if (next != NULL)
-		return next->Count()+1;
-	else
-		return 1; 
-};
+  if (next != NULL)
+    return next->Count() + 1;
+  else
+    return 1;
+};
+
Index: src/molecules.cpp
===================================================================
--- src/molecules.cpp	(revision d8b94a86944c4679c88bc250bff510f242d63309)
+++ src/molecules.cpp	(revision 1907a789609146bbdebfeff09de44cb4ba21b808)
@@ -62,4 +62,5 @@
 	cell_size[0] = cell_size[2] = cell_size[5]= 20.;
 	cell_size[1] = cell_size[3] = cell_size[4]= 0.;
+	strcpy(name,"none");
 };
 
@@ -596,4 +597,19 @@
 	cerr << Verbose(1) << "molecule::RemoveBond: Function not implemented yet." << endl;
 	return false;
+};
+
+/** Set molecule::name from the basename without suffix in the given \a *filename.
+ * \param *filename filename
+ */
+void molecule::SetNameFromFilename(char *filename)
+{
+  int length = 0;
+  char *molname = strrchr(filename, '/')+sizeof(char);  // search for filename without dirs
+  char *endname = strrchr(filename, '.');
+  if ((endname == NULL) || (endname < molname))
+    length = strlen(molname);
+  else
+    length = strlen(molname) - strlen(endname);
+  strncpy(name, molname, length);
 };
 
@@ -1187,5 +1203,5 @@
 };
 
-/** Removes atom from molecule list.
+/** Removes atom from molecule list and deletes it.
  * \param *pointer atom to be removed
  * \return true - succeeded, false - atom not found in list
@@ -1201,4 +1217,23 @@
 	Trajectories.erase(pointer);
 	return remove(pointer, start, end);
+};
+
+/** Removes atom from molecule list, but does not delete it.
+ * \param *pointer atom to be removed
+ * \return true - succeeded, false - atom not found in list
+ */
+bool molecule::UnlinkAtom(atom *pointer)
+{
+  if (pointer == NULL)
+    return false;
+  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;
+  if (ElementsInMolecule[pointer->type->Z] == 0)  // was last atom of this element?
+    ElementCount--;
+  Trajectories.erase(pointer);
+  unlink(pointer);
+  return true;
 };
 
@@ -1727,5 +1762,5 @@
 						};
 						AddBond(Walker, OtherWalker); //Add the bond between the two atoms with respective indices.
-		
+
 					}
 
@@ -2529,5 +2564,4 @@
 					cerr << "KeySet file must be corrupt as there are two equal key sets therein!" << endl;
 				}
-				//FragmentList->ListOfMolecules[NumberOfFragments++] = StoreFragmentFromKeySet(out, CurrentSet, IsAngstroem);
 			}
 		}
@@ -3088,20 +3122,20 @@
 	//if (FragmentationToDo) {		// we should always store the fragments again as coordination might have changed slightly without changing bond structure
 		// allocate memory for the pointer array and transmorph graphs into full molecular fragments
-		BondFragments = new MoleculeListClass(TotalGraph.size(), AtomCount);
+		BondFragments = new MoleculeListClass();
 		int k=0;
 		for(Graph::iterator runner = TotalGraph.begin(); runner != TotalGraph.end(); runner++) {
 			KeySet test = (*runner).first;
 			*out << "Fragment No." << (*runner).second.first << " with TEFactor " << (*runner).second.second << "." << endl;
-			BondFragments->ListOfMolecules[k] = StoreFragmentFromKeySet(out, test, configuration);
+			BondFragments->insert(StoreFragmentFromKeySet(out, test, configuration));
 			k++;
 		}
-		*out << k << "/" << BondFragments->NumberOfMolecules << " fragments generated from the keysets." << endl;
+		*out << k << "/" << BondFragments->ListOfMolecules.size() << " fragments generated from the keysets." << endl;
 
 		// ===== 9. Save fragments' configuration and keyset files et al to disk ===
-		if (BondFragments->NumberOfMolecules != 0) {
+		if (BondFragments->ListOfMolecules.size() != 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;
+			*out << Verbose(1) << "Writing " << BondFragments->ListOfMolecules.size() << " possible bond fragmentation configs" << endl;
 			if (BondFragments->OutputConfigForListOfFragments(out, configuration, SortIndex))
 				*out << Verbose(1) << "All configs written." << endl;
@@ -3629,219 +3663,4 @@
 };
 
-/** Creates \a MoleculeListClass of all unique fragments of the \a molecule containing \a Order atoms or vertices.
- * The picture to have in mind is that of a DFS "snake" of a certain length \a Order, i.e. as in the infamous
- * computer game, that winds through the connected graph representing the molecule. Color (white,
- * lightgray, darkgray, black) indicates whether a vertex has been discovered so far or not. Labels will help in
- * creating only unique fragments and not additional ones with vertices simply in different sequence.
- * The Predecessor is always the one that came before in discovering, needed on backstepping. And
- * finally, the ShortestPath is needed for removing vertices from the snake stack during the back-
- * stepping.
- * \param *out output stream for debugging
- * \param Order number of atoms in each fragment
- * \param *configuration configuration for writing config files for each fragment
- * \return List of all unique fragments with \a Order atoms
- */
-/*
-MoleculeListClass * molecule::CreateListOfUniqueFragmentsOfOrder(ofstream *out, int Order, config *configuration)
-{
-	atom **PredecessorList = (atom **) Malloc(sizeof(atom *)*AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: **PredecessorList");
-	int *ShortestPathList = (int *) Malloc(sizeof(int)*AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *ShortestPathList");
-	int *Labels = (int *) Malloc(sizeof(int)*AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *Labels");
-	enum Shading *ColorVertexList = (enum Shading *) Malloc(sizeof(enum Shading)*AtomCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *ColorList");
-	enum Shading *ColorEdgeList = (enum Shading *) Malloc(sizeof(enum Shading)*BondCount, "molecule::CreateListOfUniqueFragmentsOfOrder: *ColorBondList");
-	StackClass<atom *> *RootStack = new StackClass<atom *>(AtomCount);
-	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;
-	MoleculeListClass *FragmentList = NULL;
-	atom *Walker = NULL, *OtherAtom = NULL, *Root = NULL, *Removal = NULL;
-	bond *Binder = NULL;
-	int RunningIndex = 0, FragmentCounter = 0;
-
-	*out << Verbose(1) << "Begin of CreateListOfUniqueFragmentsOfOrder." << endl;
-
-	// reset parent list
-	*out << Verbose(3) << "Resetting labels, parent, predecessor, color and shortest path lists." << endl;
-	for (int i=0;i<AtomCount;i++) { // reset all atom labels
-		// initialise each vertex as white with no predecessor, empty queue, color lightgray, not labelled, no sons
-		Labels[i] = -1;
-		SonList[i] = NULL;
-		PredecessorList[i] = NULL;
-		ColorVertexList[i] = white;
-		ShortestPathList[i] = -1;
-	}
-	for (int i=0;i<BondCount;i++)
-		ColorEdgeList[i] = white;
-	RootStack->ClearStack();	// clearstack and push first atom if exists
-	TouchedStack->ClearStack();
-	Walker = start->next;
-	while ((Walker != end)
-#ifdef ADDHYDROGEN
-	 && (Walker->type->Z == 1)
-#endif
-												) { // search for first non-hydrogen atom
-		*out << Verbose(4) << "Current Root candidate is " << Walker->Name << "." << endl;
-		Walker = Walker->next;
-	}
-	if (Walker != end)
-		RootStack->Push(Walker);
-	else
-		*out << Verbose(0) << "ERROR: Could not find an appropriate Root atom!" << endl;
-	*out << Verbose(3) << "Root " << Walker->Name << " is on AtomStack, beginning loop through all vertices ..." << endl;
-
-	///// OUTER LOOP ////////////
-	while (!RootStack->IsEmpty()) {
-		// get new root vertex from atom stack
-		Root = RootStack->PopFirst();
-		ShortestPathList[Root->nr] = 0;
-		if (Labels[Root->nr] == -1)
-			Labels[Root->nr] = RunningIndex++; // prevent it from getting again on AtomStack
-		PredecessorList[Root->nr] = Root;
-		TouchedStack->Push(Root);
-		*out << Verbose(0) << "Root for this loop is: " << Root->Name << ".\n";
-
-		// clear snake stack
-		SnakeStack->ClearStack();
-		//SnakeStack->TestImplementation(out, start->next);
-
-		///// INNER LOOP ////////////
-		// Problems:
-		// - what about cyclic bonds?
-		Walker = Root;
-		do {
-			*out << Verbose(1) << "Current Walker is: " << Walker->Name;
-			// initial setting of the new Walker: label, color, shortest path and put on stacks
-			if (Labels[Walker->nr] == -1)	{	// give atom a unique, monotonely increasing number
-				Labels[Walker->nr] = RunningIndex++;
-				RootStack->Push(Walker);
-			}
-			*out << ", has label " << Labels[Walker->nr];
-			if ((ColorVertexList[Walker->nr] == white) || ((Binder != NULL) && (ColorEdgeList[Binder->nr] == white))) {	// color it if newly discovered and push on stacks (and if within reach!)
-				if ((Binder != NULL) && (ColorEdgeList[Binder->nr] == white)) {
-					// Binder ought to be set still from last neighbour search
-					*out << ", coloring bond " << *Binder << " black";
-					ColorEdgeList[Binder->nr] = black; // mark this bond as used
-				}
-				if (ShortestPathList[Walker->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
-				}
-				if ((ShortestPathList[Walker->nr] < Order) && (ColorVertexList[Walker->nr] != darkgray)) {	// if not already on snake stack
-					SnakeStack->Push(Walker);
-					ColorVertexList[Walker->nr] = darkgray; // mark as dark gray of on snake stack
-				}
-			}
-			*out << ", SP of " << ShortestPathList[Walker->nr]	<< " and its color is " << GetColor(ColorVertexList[Walker->nr]) << "." << endl;
-
-			// then check the stack for a newly stumbled upon fragment
-			if (SnakeStack->ItemCount() == Order) { // is stack full?
-				// store the fragment if it is one and get a removal candidate
-				Removal = StoreFragmentFromStack(out, Root, Walker, Leaflet, SnakeStack, ShortestPathList, SonList, Labels, &FragmentCounter, configuration);
-				// remove the candidate if one was found
-				if (Removal != NULL) {
-					*out << Verbose(2) << "Removing item " << Removal->Name << " with SP of " << ShortestPathList[Removal->nr] << " from snake stack." << endl;
-					SnakeStack->RemoveItem(Removal);
-					ColorVertexList[Removal->nr] = lightgray; // return back to not on snake stack but explored marking
-					if (Walker == Removal) { // if the current atom is to be removed, we also have to take a step back
-						Walker = PredecessorList[Removal->nr];
-						*out << Verbose(2) << "Stepping back to " << Walker->Name << "." << endl;
-					}
-				}
-			} else
-				Removal = NULL;
-
-			// finally, look for a white neighbour as the next Walker
-			Binder = NULL;
-			if ((Removal == NULL) || (Walker != PredecessorList[Removal->nr])) {	// don't look, if a new walker has been set above
-				*out << Verbose(2) << "Snake has currently " << SnakeStack->ItemCount() << " item(s)." << endl;
-				OtherAtom = NULL; // this is actually not needed, every atom has at least one neighbour
-				if (ShortestPathList[Walker->nr] < Order) {
-					for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++) {
-						Binder = ListOfBondsPerAtom[Walker->nr][i];
-						*out << Verbose(2) << "Current bond is " << *Binder << ": ";
-						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;
-							//ColorVertexList[OtherAtom->nr] = lightgray;		// mark as explored
-						} else { // otherwise check its colour and element
-							if (
-#ifdef ADDHYDROGEN
-							(OtherAtom->type->Z != 1) &&
-#endif
-										(ColorEdgeList[Binder->nr] == white)) {	// skip hydrogen, look for unexplored vertices
-								*out << "Moving along " << GetColor(ColorEdgeList[Binder->nr]) << " bond " << Binder << " to " << ((ColorVertexList[OtherAtom->nr] == white) ? "unexplored" : "explored") << " item: " << OtherAtom->Name << "." << endl;
-								// i find it currently rather sensible to always set the predecessor in order to find one's way back
-								//if (PredecessorList[OtherAtom->nr] == NULL) {
-								PredecessorList[OtherAtom->nr] = Walker;
-								*out << Verbose(3) << "Setting Predecessor of " << OtherAtom->Name << " to " << PredecessorList[OtherAtom->nr]->Name << "." << endl;
-								//} else {
-								//	*out << Verbose(3) << "Predecessor of " << OtherAtom->Name << " is " << PredecessorList[OtherAtom->nr]->Name << "." << endl;
-								//}
-								Walker = OtherAtom;
-								break;
-							} else {
-								if (OtherAtom->type->Z == 1)
-									*out << "Links to a hydrogen atom." << endl;
-								else
-									*out << "Bond has not white but " << GetColor(ColorEdgeList[Binder->nr]) << " color." << endl;
-							}
-						}
-					}
-				} else {	// means we have stepped beyond the horizon: Return!
-					Walker = PredecessorList[Walker->nr];
-					OtherAtom = Walker;
-					*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
-					*out << Verbose(2) << "Coloring " << Walker->Name << " black." << endl;
-					ColorVertexList[Walker->nr] = black;
-					Walker = PredecessorList[Walker->nr];
-				}
-			}
-		} while ((Walker != Root) || (ColorVertexList[Root->nr] != black));
-		*out << Verbose(2) << "Inner Looping is finished." << endl;
-
-		// if we reset all AtomCount atoms, we have again technically O(N^2) ...
-		*out << Verbose(2) << "Resetting lists." << endl;
-		Walker = NULL;
-		Binder = NULL;
-		while (!TouchedStack->IsEmpty()) {
-			Walker = TouchedStack->PopLast();
-			*out << Verbose(3) << "Re-initialising entries of " << *Walker << "." << endl;
-			for(int i=0;i<NumberOfBondsPerAtom[Walker->nr];i++)
-				ColorEdgeList[ListOfBondsPerAtom[Walker->nr][i]->nr] = white;
-			PredecessorList[Walker->nr] = NULL;
-			ColorVertexList[Walker->nr] = white;
-			ShortestPathList[Walker->nr] = -1;
-		}
-	}
-	*out << Verbose(1) << "Outer Looping over all vertices is done." << endl;
-
-	// copy together
-	*out << Verbose(1) << "Copying all fragments into MoleculeList structure." << endl;
-	FragmentList = new MoleculeListClass(FragmentCounter, AtomCount);
-	RunningIndex = 0;
-	while ((Leaflet != NULL) && (RunningIndex < FragmentCounter))	{
-		FragmentList->ListOfMolecules[RunningIndex++] = Leaflet->Leaf;
-		Leaflet->Leaf = NULL; // prevent molecule from being removed
-		TempLeaf = Leaflet;
-		Leaflet = Leaflet->previous;
-		delete(TempLeaf);
-	};
-
-	// free memory and exit
-	Free((void **)&PredecessorList, "molecule::CreateListOfUniqueFragmentsOfOrder: **PredecessorList");
-	Free((void **)&ShortestPathList, "molecule::CreateListOfUniqueFragmentsOfOrder: *ShortestPathList");
-	Free((void **)&Labels, "molecule::CreateListOfUniqueFragmentsOfOrder: *Labels");
-	Free((void **)&ColorVertexList, "molecule::CreateListOfUniqueFragmentsOfOrder: *ColorList");
-	delete(RootStack);
-	delete(TouchedStack);
-	delete(SnakeStack);
-
-	*out << Verbose(1) << "End of CreateListOfUniqueFragmentsOfOrder." << endl;
-	return FragmentList;
-};
-*/
-
 /** Structure containing all values in power set combination generation.
  */
Index: src/molecules.hpp
===================================================================
--- src/molecules.hpp	(revision d8b94a86944c4679c88bc250bff510f242d63309)
+++ src/molecules.hpp	(revision 1907a789609146bbdebfeff09de44cb4ba21b808)
@@ -68,4 +68,7 @@
 #define DistanceMultiMap multimap <double, pair < PointMap::iterator, PointMap::iterator> >
 #define DistanceMultiMapPair pair <double, pair < PointMap::iterator, PointMap::iterator> >
+
+#define MoleculeList list <molecule *>
+#define MoleculeListTest pair <MoleculeList::iterator, bool>
 
 /******************************** Some small functions and/or structures **********************************/
@@ -142,5 +145,5 @@
 	~atom();
 
-	bool Output(int ElementNo, int AtomNo, ofstream *out) const;
+	bool Output(int ElementNo, int AtomNo, ofstream *out, const char *comment = NULL) const;
 	bool OutputXYZLine(ofstream *out) const;
 	atom *GetTrueFather();
@@ -216,4 +219,8 @@
 		int NoCyclicBonds;	//!< number of cyclic bonds in molecule, by DepthFirstSearchAnalysis()
 		double BondDistance;	//!< typical bond distance used in CreateAdjacencyList() and furtheron
+		bool ActiveFlag;    //!< in a MoleculeListClass used to discern active from inactive molecules
+		Vector Center;      //!< Center of molecule in a global box
+		char name[MAXSTRINGSIZE];         //!< arbitrary name
+		int IndexNr;        //!< index of molecule in a MoleculeListClass
 
 	molecule(periodentafel *teil);
@@ -223,4 +230,5 @@
 	bool AddAtom(atom *pointer);
 	bool RemoveAtom(atom *pointer);
+	bool UnlinkAtom(atom *pointer);
 	bool CleanupMolecule();
 
@@ -252,4 +260,5 @@
 	Vector * DetermineCenterOfGravity(ofstream *out);
 	Vector * DetermineCenterOfAll(ofstream *out);
+	void SetNameFromFilename(char *filename);
 	void SetBoxDimension(Vector *dim);
 	double * ReturnFullMatrixforSymmetric(double *cell_size);
@@ -326,17 +335,26 @@
 class MoleculeListClass {
 	public:
-		molecule **ListOfMolecules;	 //!< pointer list of fragment molecules to check for equality
-		int NumberOfMolecules;				//!< Number of entries in \a **FragmentList and of to be returned one.
-		int NumberOfTopAtoms;				 //!< Number of atoms in the molecule from which all fragments originate
+	  MoleculeList ListOfMolecules; //!< List of the contained molecules
+	  int MaxIndex;
 
 	MoleculeListClass();
-	MoleculeListClass(int Num, int NumAtoms);
 	~MoleculeListClass();
 
-	/// Output configs.
 	bool AddHydrogenCorrection(ofstream *out, char *path);
 	bool StoreForcesFile(ofstream *out, char *path, int *SortIndex);
+	bool insert(molecule *mol);
+	molecule * ReturnIndex(int index);
 	bool OutputConfigForListOfFragments(ofstream *out, config *configuration, int *SortIndex);
+	int NumberOfActiveMolecules();
+	void Enumerate(ofstream *out);
 	void Output(ofstream *out);
+
+	// merging of molecules
+  bool SimpleMerge(molecule *mol, molecule *srcmol);
+  bool SimpleAdd(molecule *mol, molecule *srcmol);
+  bool SimpleMultiMerge(molecule *mol, int *src, int N);
+  bool SimpleMultiAdd(molecule *mol, int *src, int N);
+  bool ScatterMerge(molecule *mol, int *src, int N);
+  bool EmbedMerge(molecule *mol, molecule *srcmol);
 
 	private:
@@ -454,5 +472,5 @@
 	bool Save(const char *filename, periodentafel *periode, molecule *mol) const;
 	bool SaveMPQC(const char *filename, molecule *mol) const;
-	void Edit(molecule *mol);
+	void Edit();
 	bool GetIsAngstroem() const;
 	char *GetDefaultPath() const;
