Index: src/builder.cpp
===================================================================
--- src/builder.cpp	(revision a5b2c3a1b29eef3d3e6b57dfbcea27b22ec153cc)
+++ src/builder.cpp	(revision 375b4589d3335eeefe2cda072f13207108281969)
@@ -15,33 +15,33 @@
  * \section about About the Program
  *
- *	Molecuilder is a short program, written in C++, that enables the construction of a coordinate set for the
- *	atoms making up an molecule by the successive statement of binding angles and distances and referencing to
- *	already constructed atoms.
+ *  Molecuilder is a short program, written in C++, that enables the construction of a coordinate set for the
+ *  atoms making up an molecule by the successive statement of binding angles and distances and referencing to
+ *  already constructed atoms.
  *
- *	A configuration file may be written that is compatible to the format used by PCP - a parallel Car-Parrinello
- *	molecular dynamics implementation.
+ *  A configuration file may be written that is compatible to the format used by PCP - a parallel Car-Parrinello
+ *  molecular dynamics implementation.
  *
  * \section install Installation
  *
- *	Installation should without problems succeed as follows:
- *	-# ./configure (or: mkdir build;mkdir run;cd build; ../configure --bindir=../run)
- *	-# make
- *	-# make install
+ *  Installation should without problems succeed as follows:
+ *  -# ./configure (or: mkdir build;mkdir run;cd build; ../configure --bindir=../run)
+ *  -# make
+ *  -# make install
  *
- *	Further useful commands are
- *	-# make clean uninstall: deletes .o-files and removes executable from the given binary directory\n
- *	-# make doxygen-doc: Creates these html pages out of the documented source
+ *  Further useful commands are
+ *  -# make clean uninstall: deletes .o-files and removes executable from the given binary directory\n
+ *  -# make doxygen-doc: Creates these html pages out of the documented source
  *
  * \section run Running
  *
- *	The program can be executed by running: ./molecuilder
+ *  The program can be executed by running: ./molecuilder
  *
- *	Note, that it uses a database, called "elements.db", in the executable's directory. If the file is not found,
- *	it is created and any given data on elements of the periodic table will be stored therein and re-used on
- *	later re-execution.
+ *  Note, that it uses a database, called "elements.db", in the executable's directory. If the file is not found,
+ *  it is created and any given data on elements of the periodic table will be stored therein and re-used on
+ *  later re-execution.
  *
  * \section ref References
  *
- *	For the special configuration file format, see the documentation of pcp.
+ *  For the special configuration file format, see the documentation of pcp.
  *
  */
@@ -63,36 +63,36 @@
 static void AddAtoms(periodentafel *periode, molecule *mol)
 {
-	atom *first, *second, *third, *fourth;
-	Vector **atoms;
-	Vector x,y,z,n;	// coordinates for absolute point in cell volume
-	double a,b,c;
-	char choice;	// menu choice char
-	bool valid;
-
-	cout << Verbose(0) << "===========ADD ATOM============================" << endl;
-	cout << Verbose(0) << " a - state absolute coordinates of atom" << endl;
-	cout << Verbose(0) << " b - state relative coordinates of atom wrt to reference point" << endl;
-	cout << Verbose(0) << " c - state relative coordinates of atom wrt to already placed atom" << endl;
-	cout << Verbose(0) << " d - state two atoms, two angles and a distance" << endl;
-	cout << Verbose(0) << " e - least square distance position to a set of atoms" << endl;
-	cout << Verbose(0) << "all else - go back" << endl;
-	cout << Verbose(0) << "===============================================" << endl;
-	cout << Verbose(0) << "Note: Specifiy angles in degrees not multiples of Pi!" << endl;
-	cout << Verbose(0) << "INPUT: ";
-	cin >> choice;
-
-	switch (choice) {
+  atom *first, *second, *third, *fourth;
+  Vector **atoms;
+  Vector x,y,z,n;  // coordinates for absolute point in cell volume
+  double a,b,c;
+  char choice;  // menu choice char
+  bool valid;
+
+  cout << Verbose(0) << "===========ADD ATOM============================" << endl;
+  cout << Verbose(0) << " a - state absolute coordinates of atom" << endl;
+  cout << Verbose(0) << " b - state relative coordinates of atom wrt to reference point" << endl;
+  cout << Verbose(0) << " c - state relative coordinates of atom wrt to already placed atom" << endl;
+  cout << Verbose(0) << " d - state two atoms, two angles and a distance" << endl;
+  cout << Verbose(0) << " e - least square distance position to a set of atoms" << endl;
+  cout << Verbose(0) << "all else - go back" << endl;
+  cout << Verbose(0) << "===============================================" << endl;
+  cout << Verbose(0) << "Note: Specifiy angles in degrees not multiples of Pi!" << endl;
+  cout << Verbose(0) << "INPUT: ";
+  cin >> choice;
+
+  switch (choice) {
     default:
       cout << Verbose(0) << "Not a valid choice." << endl;
       break;
-			case 'a': // absolute coordinates of atom
+      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
-				break;
-
-			case 'b': // relative coordinates of atom wrt to reference point
+        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;
@@ -106,9 +106,9 @@
           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->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;
@@ -122,6 +122,6 @@
           }
         } while (!(valid = mol->CheckBounds((const Vector *)&first->x)));
-        first->type = periode->AskElement();	// give type
-        mol->AddAtom(first);	// add to molecule
+        first->type = periode->AskElement();  // give type
+        mol->AddAtom(first);  // add to molecule
         break;
 
@@ -224,9 +224,9 @@
           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->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]);
@@ -248,12 +248,12 @@
 
           first->x.Output((ofstream *)&cout);
-          first->type = periode->AskElement();	// give type
-          mol->AddAtom(first);	// add to molecule
+          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;
-	};
+        break;
+  };
 };
 
@@ -263,54 +263,54 @@
 static void CenterAtoms(molecule *mol)
 {
-	Vector x, y, helper;
-	char choice;	// menu choice char
-
-	cout << Verbose(0) << "===========CENTER ATOMS=========================" << endl;
-	cout << Verbose(0) << " a - on origin" << endl;
-	cout << Verbose(0) << " b - on center of gravity" << endl;
-	cout << Verbose(0) << " c - within box with additional boundary" << endl;
-	cout << Verbose(0) << " d - within given simulation box" << 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 'a':
-			cout << Verbose(0) << "Centering atoms in config file on origin." << endl;
-			mol->CenterOrigin((ofstream *)&cout, &x);
-			break;
-		case 'b':
-			cout << Verbose(0) << "Centering atoms in config file on center of gravity." << endl;
-			mol->CenterGravity((ofstream *)&cout, &x);
-			break;
-		case 'c':
-			cout << Verbose(0) << "Centering atoms in config file within given additional boundary." << endl;
-			for (int i=0;i<NDIM;i++) {
-				cout << Verbose(0) << "Enter axis " << i << " boundary: ";
-				cin >> y.x[i];
-			}
-			mol->CenterEdge((ofstream *)&cout, &x);	// make every coordinate positive
-			mol->Translate(&y); // translate by boundary
-			helper.CopyVector(&y);
-			helper.Scale(2.);
-			helper.AddVector(&x);
-			mol->SetBoxDimension(&helper);	// update Box of atoms by boundary
-			break;
-		case 'd':
-			cout << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
-			for (int i=0;i<NDIM;i++) {
-				cout << Verbose(0) << "Enter axis " << i << " boundary: ";
-				cin >> x.x[i];
-			}
-			// center
-			mol->CenterInBox((ofstream *)&cout, &x);
-			// update Box of atoms by boundary
-			mol->SetBoxDimension(&x);
-			break;
-	}
+  Vector x, y, helper;
+  char choice;  // menu choice char
+
+  cout << Verbose(0) << "===========CENTER ATOMS=========================" << endl;
+  cout << Verbose(0) << " a - on origin" << endl;
+  cout << Verbose(0) << " b - on center of gravity" << endl;
+  cout << Verbose(0) << " c - within box with additional boundary" << endl;
+  cout << Verbose(0) << " d - within given simulation box" << 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 'a':
+      cout << Verbose(0) << "Centering atoms in config file on origin." << endl;
+      mol->CenterOrigin((ofstream *)&cout, &x);
+      break;
+    case 'b':
+      cout << Verbose(0) << "Centering atoms in config file on center of gravity." << endl;
+      mol->CenterGravity((ofstream *)&cout, &x);
+      break;
+    case 'c':
+      cout << Verbose(0) << "Centering atoms in config file within given additional boundary." << endl;
+      for (int i=0;i<NDIM;i++) {
+        cout << Verbose(0) << "Enter axis " << i << " boundary: ";
+        cin >> y.x[i];
+      }
+      mol->CenterEdge((ofstream *)&cout, &x);  // make every coordinate positive
+      mol->Translate(&y); // translate by boundary
+      helper.CopyVector(&y);
+      helper.Scale(2.);
+      helper.AddVector(&x);
+      mol->SetBoxDimension(&helper);  // update Box of atoms by boundary
+      break;
+    case 'd':
+      cout << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
+      for (int i=0;i<NDIM;i++) {
+        cout << Verbose(0) << "Enter axis " << i << " boundary: ";
+        cin >> x.x[i];
+      }
+      // center
+      mol->CenterInBox((ofstream *)&cout, &x);
+      // update Box of atoms by boundary
+      mol->SetBoxDimension(&x);
+      break;
+  }
 };
 
@@ -321,65 +321,65 @@
 static void AlignAtoms(periodentafel *periode, molecule *mol)
 {
-	atom *first, *second, *third;
-	Vector x,n;
-	char choice;	// menu choice char
-
-	cout << Verbose(0) << "===========ALIGN ATOMS=========================" << endl;
-	cout << Verbose(0) << " a - state three atoms defining align plane" << endl;
-	cout << Verbose(0) << " b - state alignment vector" << endl;
-	cout << Verbose(0) << " c - state two atoms in alignment direction" << endl;
-	cout << Verbose(0) << " d - align automatically by least square fit" << endl;
-	cout << Verbose(0) << "all else - go back" << endl;
-	cout << Verbose(0) << "===============================================" << endl;
-	cout << Verbose(0) << "INPUT: ";
-	cin >> choice;
-
-	switch (choice) {
-		default:
-		case 'a': // three atoms defining mirror plane
-			first = mol->AskAtom("Enter first atom: ");
-			second = mol->AskAtom("Enter second atom: ");
-			third = mol->AskAtom("Enter third atom: ");
-
-			n.MakeNormalVector((const Vector *)&first->x,(const Vector *)&second->x,(const Vector *)&third->x);
-			break;
-		case 'b': // normal vector of mirror plane
-			cout << Verbose(0) << "Enter normal vector of mirror plane." << endl;
-			n.AskPosition(mol->cell_size,0);
-			n.Normalize();
-			break;
-		case 'c': // three atoms defining mirror plane
-			first = mol->AskAtom("Enter first atom: ");
-			second = mol->AskAtom("Enter second atom: ");
-
-			n.CopyVector((const Vector *)&first->x);
-			n.SubtractVector((const Vector *)&second->x);
-			n.Normalize();
-			break;
-		case 'd':
-			char shorthand[4];
-			Vector a;
-			struct lsq_params param;
-			do {
-				fprintf(stdout, "Enter the element of atoms to be chosen: ");
-				fscanf(stdin, "%3s", shorthand);
-			} while ((param.type = periode->FindElement(shorthand)) == NULL);
-			cout << Verbose(0) << "Element is " << param.type->name << endl;
-			mol->GetAlignvector(&param);
-			for (int i=NDIM;i--;) {
-				x.x[i] = gsl_vector_get(param.x,i);
-				n.x[i] = gsl_vector_get(param.x,i+NDIM);
-			}
-			gsl_vector_free(param.x);
-			cout << Verbose(0) << "Offset vector: ";
-			x.Output((ofstream *)&cout);
-			cout << Verbose(0) << endl;
-			n.Normalize();
-			break;
-	};
-	cout << Verbose(0) << "Alignment vector: ";
-	n.Output((ofstream *)&cout);
-	cout << Verbose(0) << endl;
-	mol->Align(&n);
+  atom *first, *second, *third;
+  Vector x,n;
+  char choice;  // menu choice char
+
+  cout << Verbose(0) << "===========ALIGN ATOMS=========================" << endl;
+  cout << Verbose(0) << " a - state three atoms defining align plane" << endl;
+  cout << Verbose(0) << " b - state alignment vector" << endl;
+  cout << Verbose(0) << " c - state two atoms in alignment direction" << endl;
+  cout << Verbose(0) << " d - align automatically by least square fit" << endl;
+  cout << Verbose(0) << "all else - go back" << endl;
+  cout << Verbose(0) << "===============================================" << endl;
+  cout << Verbose(0) << "INPUT: ";
+  cin >> choice;
+
+  switch (choice) {
+    default:
+    case 'a': // three atoms defining mirror plane
+      first = mol->AskAtom("Enter first atom: ");
+      second = mol->AskAtom("Enter second atom: ");
+      third = mol->AskAtom("Enter third atom: ");
+
+      n.MakeNormalVector((const Vector *)&first->x,(const Vector *)&second->x,(const Vector *)&third->x);
+      break;
+    case 'b': // normal vector of mirror plane
+      cout << Verbose(0) << "Enter normal vector of mirror plane." << endl;
+      n.AskPosition(mol->cell_size,0);
+      n.Normalize();
+      break;
+    case 'c': // three atoms defining mirror plane
+      first = mol->AskAtom("Enter first atom: ");
+      second = mol->AskAtom("Enter second atom: ");
+
+      n.CopyVector((const Vector *)&first->x);
+      n.SubtractVector((const Vector *)&second->x);
+      n.Normalize();
+      break;
+    case 'd':
+      char shorthand[4];
+      Vector a;
+      struct lsq_params param;
+      do {
+        fprintf(stdout, "Enter the element of atoms to be chosen: ");
+        fscanf(stdin, "%3s", shorthand);
+      } while ((param.type = periode->FindElement(shorthand)) == NULL);
+      cout << Verbose(0) << "Element is " << param.type->name << endl;
+      mol->GetAlignvector(&param);
+      for (int i=NDIM;i--;) {
+        x.x[i] = gsl_vector_get(param.x,i);
+        n.x[i] = gsl_vector_get(param.x,i+NDIM);
+      }
+      gsl_vector_free(param.x);
+      cout << Verbose(0) << "Offset vector: ";
+      x.Output((ofstream *)&cout);
+      cout << Verbose(0) << endl;
+      n.Normalize();
+      break;
+  };
+  cout << Verbose(0) << "Alignment vector: ";
+  n.Output((ofstream *)&cout);
+  cout << Verbose(0) << endl;
+  mol->Align(&n);
 };
 
@@ -389,44 +389,44 @@
 static void MirrorAtoms(molecule *mol)
 {
-	atom *first, *second, *third;
-	Vector n;
-	char choice;	// menu choice char
-
-	cout << Verbose(0) << "===========MIRROR ATOMS=========================" << endl;
-	cout << Verbose(0) << " a - state three atoms defining mirror plane" << endl;
-	cout << Verbose(0) << " b - state normal vector of mirror plane" << endl;
-	cout << Verbose(0) << " c - state two atoms in normal direction" << endl;
-	cout << Verbose(0) << "all else - go back" << endl;
-	cout << Verbose(0) << "===============================================" << endl;
-	cout << Verbose(0) << "INPUT: ";
-	cin >> choice;
-
-	switch (choice) {
-		default:
-		case 'a': // three atoms defining mirror plane
-			first = mol->AskAtom("Enter first atom: ");
-			second = mol->AskAtom("Enter second atom: ");
-			third = mol->AskAtom("Enter third atom: ");
-
-			n.MakeNormalVector((const Vector *)&first->x,(const Vector *)&second->x,(const Vector *)&third->x);
-			break;
-		case 'b': // normal vector of mirror plane
-			cout << Verbose(0) << "Enter normal vector of mirror plane." << endl;
-			n.AskPosition(mol->cell_size,0);
-			n.Normalize();
-			break;
-		case 'c': // three atoms defining mirror plane
-			first = mol->AskAtom("Enter first atom: ");
-			second = mol->AskAtom("Enter second atom: ");
-
-			n.CopyVector((const Vector *)&first->x);
-			n.SubtractVector((const Vector *)&second->x);
-			n.Normalize();
-			break;
-	};
-	cout << Verbose(0) << "Normal vector: ";
-	n.Output((ofstream *)&cout);
-	cout << Verbose(0) << endl;
-	mol->Mirror((const Vector *)&n);
+  atom *first, *second, *third;
+  Vector n;
+  char choice;  // menu choice char
+
+  cout << Verbose(0) << "===========MIRROR ATOMS=========================" << endl;
+  cout << Verbose(0) << " a - state three atoms defining mirror plane" << endl;
+  cout << Verbose(0) << " b - state normal vector of mirror plane" << endl;
+  cout << Verbose(0) << " c - state two atoms in normal direction" << endl;
+  cout << Verbose(0) << "all else - go back" << endl;
+  cout << Verbose(0) << "===============================================" << endl;
+  cout << Verbose(0) << "INPUT: ";
+  cin >> choice;
+
+  switch (choice) {
+    default:
+    case 'a': // three atoms defining mirror plane
+      first = mol->AskAtom("Enter first atom: ");
+      second = mol->AskAtom("Enter second atom: ");
+      third = mol->AskAtom("Enter third atom: ");
+
+      n.MakeNormalVector((const Vector *)&first->x,(const Vector *)&second->x,(const Vector *)&third->x);
+      break;
+    case 'b': // normal vector of mirror plane
+      cout << Verbose(0) << "Enter normal vector of mirror plane." << endl;
+      n.AskPosition(mol->cell_size,0);
+      n.Normalize();
+      break;
+    case 'c': // three atoms defining mirror plane
+      first = mol->AskAtom("Enter first atom: ");
+      second = mol->AskAtom("Enter second atom: ");
+
+      n.CopyVector((const Vector *)&first->x);
+      n.SubtractVector((const Vector *)&second->x);
+      n.Normalize();
+      break;
+  };
+  cout << Verbose(0) << "Normal vector: ";
+  n.Output((ofstream *)&cout);
+  cout << Verbose(0) << endl;
+  mol->Mirror((const Vector *)&n);
 };
 
@@ -436,60 +436,60 @@
 static void RemoveAtoms(molecule *mol)
 {
-	atom *first, *second, *third;
-	int axis;
-	double tmp1, tmp2;
-	char choice;	// menu choice char
-
-	cout << Verbose(0) << "===========REMOVE ATOMS=========================" << endl;
-	cout << Verbose(0) << " a - state atom for removal by number" << endl;
-	cout << Verbose(0) << " b - keep only in radius around atom" << endl;
-	cout << Verbose(0) << " c - remove this with one axis greater value" << endl;
-	cout << Verbose(0) << "all else - go back" << endl;
-	cout << Verbose(0) << "===============================================" << endl;
-	cout << Verbose(0) << "INPUT: ";
-	cin >> choice;
-
-	switch (choice) {
-		default:
-		case 'a':
-			if (mol->RemoveAtom(mol->AskAtom("Enter number of atom within molecule: ")))
-				cout << Verbose(1) << "Atom removed." << endl;
-			else
-				cout << Verbose(1) << "Atom not found." << endl;
-			break;
-		case 'b':
-		  third = mol->AskAtom("Enter number of atom as reference point: ");
-			cout << Verbose(0) << "Enter radius: ";
-			cin >> tmp1;
-			first = mol->start;
-      second = first->next;
-			while(second != mol->end) {
-				first = second;
-        second = first->next;
-				if (first->x.DistanceSquared((const Vector *)&third->x) > tmp1*tmp1) // distance to first above radius ...
-					mol->RemoveAtom(first);
-			}
-			break;
-		case 'c':
-			cout << Verbose(0) << "Which axis is it: ";
-			cin >> axis;
-			cout << Verbose(0) << "Lower boundary: ";
-			cin >> tmp1;
-			cout << Verbose(0) << "Upper boundary: ";
-			cin >> tmp2;
-			first = mol->start;
+  atom *first, *second, *third;
+  int axis;
+  double tmp1, tmp2;
+  char choice;  // menu choice char
+
+  cout << Verbose(0) << "===========REMOVE ATOMS=========================" << endl;
+  cout << Verbose(0) << " a - state atom for removal by number" << endl;
+  cout << Verbose(0) << " b - keep only in radius around atom" << endl;
+  cout << Verbose(0) << " c - remove this with one axis greater value" << endl;
+  cout << Verbose(0) << "all else - go back" << endl;
+  cout << Verbose(0) << "===============================================" << endl;
+  cout << Verbose(0) << "INPUT: ";
+  cin >> choice;
+
+  switch (choice) {
+    default:
+    case 'a':
+      if (mol->RemoveAtom(mol->AskAtom("Enter number of atom within molecule: ")))
+        cout << Verbose(1) << "Atom removed." << endl;
+      else
+        cout << Verbose(1) << "Atom not found." << endl;
+      break;
+    case 'b':
+      third = mol->AskAtom("Enter number of atom as reference point: ");
+      cout << Verbose(0) << "Enter radius: ";
+      cin >> tmp1;
+      first = mol->start;
       second = first->next;
       while(second != mol->end) {
         first = second;
         second = first->next;
-				if ((first->x.x[axis] < tmp1) || (first->x.x[axis] > tmp2)) {// out of boundary ...
-				  //cout << "Atom " << *first << " with " << first->x.x[axis] << " on axis " << axis << " is out of bounds [" << tmp1 << "," << tmp2 << "]." << endl;
-					mol->RemoveAtom(first);
-				}
-			}
-			break;
-	};
-	//mol->Output((ofstream *)&cout);
-	choice = 'r';
+        if (first->x.DistanceSquared((const Vector *)&third->x) > tmp1*tmp1) // distance to first above radius ...
+          mol->RemoveAtom(first);
+      }
+      break;
+    case 'c':
+      cout << Verbose(0) << "Which axis is it: ";
+      cin >> axis;
+      cout << Verbose(0) << "Lower boundary: ";
+      cin >> tmp1;
+      cout << Verbose(0) << "Upper boundary: ";
+      cin >> tmp2;
+      first = mol->start;
+      second = first->next;
+      while(second != mol->end) {
+        first = second;
+        second = first->next;
+        if ((first->x.x[axis] < tmp1) || (first->x.x[axis] > tmp2)) {// out of boundary ...
+          //cout << "Atom " << *first << " with " << first->x.x[axis] << " on axis " << axis << " is out of bounds [" << tmp1 << "," << tmp2 << "]." << endl;
+          mol->RemoveAtom(first);
+        }
+      }
+      break;
+  };
+  //mol->Output((ofstream *)&cout);
+  choice = 'r';
 };
 
@@ -500,107 +500,107 @@
 static void MeasureAtoms(periodentafel *periode, molecule *mol, config *configuration)
 {
-	atom *first, *second, *third;
-	Vector x,y;
-	double min[256], tmp1, tmp2, tmp3;
-	int Z;
-	char choice;	// menu choice char
-
-	cout << Verbose(0) << "===========MEASURE ATOMS=========================" << endl;
-	cout << Verbose(0) << " a - calculate bond length between one atom and all others" << endl;
-	cout << Verbose(0) << " b - calculate bond length between two atoms" << endl;
-	cout << Verbose(0) << " c - calculate bond angle" << endl;
-	cout << Verbose(0) << " d - calculate principal axis of the system" << endl;
-	cout << Verbose(0) << " e - calculate volume of the convex envelope" << endl;
-	cout << Verbose(0) << " f - calculate temperature from current velocity" << endl;
-	cout << Verbose(0) << " g - output all temperatures per step from velocities" << endl;
-	cout << Verbose(0) << "all else - go back" << endl;
-	cout << Verbose(0) << "===============================================" << endl;
-	cout << Verbose(0) << "INPUT: ";
-	cin >> choice;
-
-	switch(choice) {
-		default:
-			cout << Verbose(1) << "Not a valid choice." << endl;
-			break;
-		case 'a':
-			first = mol->AskAtom("Enter first atom: ");
-			for (int i=MAX_ELEMENTS;i--;)
-				min[i] = 0.;
-
-			second = mol->start;
-			while ((second->next != mol->end)) {
-				second = second->next; // advance
-				Z = second->type->Z;
-				tmp1 = 0.;
-				if (first != second) {
-					x.CopyVector((const Vector *)&first->x);
-					x.SubtractVector((const Vector *)&second->x);
-					tmp1 = x.Norm();
-				}
-				if ((tmp1 != 0.) && ((min[Z] == 0.) || (tmp1 < min[Z]))) min[Z] = tmp1;
-				//cout << Verbose(0) << "Bond length between Atom " << first->nr << " and " << second->nr << ": " << tmp1 << " a.u." << endl;
-			}
-			for (int i=MAX_ELEMENTS;i--;)
-				if (min[i] != 0.) cout << Verbose(0) << "Minimum Bond length between " << first->type->name << " Atom " << first->nr << " and next Ion of type " << (periode->FindElement(i))->name << ": " << min[i] << " a.u." << endl;
-			break;
-
-		case 'b':
-			first = mol->AskAtom("Enter first atom: ");
-			second = mol->AskAtom("Enter second atom: ");
-			for (int i=NDIM;i--;)
-				min[i] = 0.;
-			x.CopyVector((const Vector *)&first->x);
-			x.SubtractVector((const Vector *)&second->x);
-			tmp1 = x.Norm();
-			cout << Verbose(1) << "Distance vector is ";
-			x.Output((ofstream *)&cout);
-			cout << "." << endl << "Norm of distance is " << tmp1 << "." << endl;
-			break;
-
-		case 'c':
-			cout << Verbose(0) << "Evaluating bond angle between three - first, central, last - atoms." << endl;
-			first = mol->AskAtom("Enter first atom: ");
-			second = mol->AskAtom("Enter central atom: ");
-			third	= mol->AskAtom("Enter last atom: ");
-			tmp1 = tmp2 = tmp3 = 0.;
-			x.CopyVector((const Vector *)&first->x);
-			x.SubtractVector((const Vector *)&second->x);
-			y.CopyVector((const Vector *)&third->x);
-			y.SubtractVector((const Vector *)&second->x);
-			cout << Verbose(0) << "Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ";
-			cout << Verbose(0) << (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.) << " degrees" << endl;
-			break;
-		case 'd':
-			cout << Verbose(0) << "Evaluating prinicipal axis." << endl;
-			cout << Verbose(0) << "Shall we rotate? [0/1]: ";
-			cin >> Z;
-			if ((Z >=0) && (Z <=1))
-				mol->PrincipalAxisSystem((ofstream *)&cout, (bool)Z);
-			else
-				mol->PrincipalAxisSystem((ofstream *)&cout, false);
-			break;
-		case 'e':
-			cout << Verbose(0) << "Evaluating volume of the convex envelope.";
-			VolumeOfConvexEnvelope((ofstream *)&cout, NULL, configuration, NULL, mol);
-			break;
-		case 'f':
-			mol->OutputTemperatureFromTrajectories((ofstream *)&cout, mol->MDSteps-1, mol->MDSteps, (ofstream *)&cout);
-			break;
-		case 'g':
-			{
-				char filename[255];
-				cout << "Please enter filename: " << endl;
-				cin >> filename;
-				cout << Verbose(1) << "Storing temperatures in " << filename << "." << endl;
-				ofstream *output = new ofstream(filename, ios::trunc);
-				if (!mol->OutputTemperatureFromTrajectories((ofstream *)&cout, 0, mol->MDSteps, output))
-					cout << Verbose(2) << "File could not be written." << endl;
-				else
-					cout << Verbose(2) << "File stored." << endl;
-				output->close();
-				delete(output);
-			}
-			break;
-	}
+  atom *first, *second, *third;
+  Vector x,y;
+  double min[256], tmp1, tmp2, tmp3;
+  int Z;
+  char choice;  // menu choice char
+
+  cout << Verbose(0) << "===========MEASURE ATOMS=========================" << endl;
+  cout << Verbose(0) << " a - calculate bond length between one atom and all others" << endl;
+  cout << Verbose(0) << " b - calculate bond length between two atoms" << endl;
+  cout << Verbose(0) << " c - calculate bond angle" << endl;
+  cout << Verbose(0) << " d - calculate principal axis of the system" << endl;
+  cout << Verbose(0) << " e - calculate volume of the convex envelope" << endl;
+  cout << Verbose(0) << " f - calculate temperature from current velocity" << endl;
+  cout << Verbose(0) << " g - output all temperatures per step from velocities" << endl;
+  cout << Verbose(0) << "all else - go back" << endl;
+  cout << Verbose(0) << "===============================================" << endl;
+  cout << Verbose(0) << "INPUT: ";
+  cin >> choice;
+
+  switch(choice) {
+    default:
+      cout << Verbose(1) << "Not a valid choice." << endl;
+      break;
+    case 'a':
+      first = mol->AskAtom("Enter first atom: ");
+      for (int i=MAX_ELEMENTS;i--;)
+        min[i] = 0.;
+
+      second = mol->start;
+      while ((second->next != mol->end)) {
+        second = second->next; // advance
+        Z = second->type->Z;
+        tmp1 = 0.;
+        if (first != second) {
+          x.CopyVector((const Vector *)&first->x);
+          x.SubtractVector((const Vector *)&second->x);
+          tmp1 = x.Norm();
+        }
+        if ((tmp1 != 0.) && ((min[Z] == 0.) || (tmp1 < min[Z]))) min[Z] = tmp1;
+        //cout << Verbose(0) << "Bond length between Atom " << first->nr << " and " << second->nr << ": " << tmp1 << " a.u." << endl;
+      }
+      for (int i=MAX_ELEMENTS;i--;)
+        if (min[i] != 0.) cout << Verbose(0) << "Minimum Bond length between " << first->type->name << " Atom " << first->nr << " and next Ion of type " << (periode->FindElement(i))->name << ": " << min[i] << " a.u." << endl;
+      break;
+
+    case 'b':
+      first = mol->AskAtom("Enter first atom: ");
+      second = mol->AskAtom("Enter second atom: ");
+      for (int i=NDIM;i--;)
+        min[i] = 0.;
+      x.CopyVector((const Vector *)&first->x);
+      x.SubtractVector((const Vector *)&second->x);
+      tmp1 = x.Norm();
+      cout << Verbose(1) << "Distance vector is ";
+      x.Output((ofstream *)&cout);
+      cout << "." << endl << "Norm of distance is " << tmp1 << "." << endl;
+      break;
+
+    case 'c':
+      cout << Verbose(0) << "Evaluating bond angle between three - first, central, last - atoms." << endl;
+      first = mol->AskAtom("Enter first atom: ");
+      second = mol->AskAtom("Enter central atom: ");
+      third  = mol->AskAtom("Enter last atom: ");
+      tmp1 = tmp2 = tmp3 = 0.;
+      x.CopyVector((const Vector *)&first->x);
+      x.SubtractVector((const Vector *)&second->x);
+      y.CopyVector((const Vector *)&third->x);
+      y.SubtractVector((const Vector *)&second->x);
+      cout << Verbose(0) << "Bond angle between first atom Nr." << first->nr << ", central atom Nr." << second->nr << " and last atom Nr." << third->nr << ": ";
+      cout << Verbose(0) << (acos(x.ScalarProduct((const Vector *)&y)/(y.Norm()*x.Norm()))/M_PI*180.) << " degrees" << endl;
+      break;
+    case 'd':
+      cout << Verbose(0) << "Evaluating prinicipal axis." << endl;
+      cout << Verbose(0) << "Shall we rotate? [0/1]: ";
+      cin >> Z;
+      if ((Z >=0) && (Z <=1))
+        mol->PrincipalAxisSystem((ofstream *)&cout, (bool)Z);
+      else
+        mol->PrincipalAxisSystem((ofstream *)&cout, false);
+      break;
+    case 'e':
+      cout << Verbose(0) << "Evaluating volume of the convex envelope.";
+      VolumeOfConvexEnvelope((ofstream *)&cout, NULL, configuration, NULL, mol);
+      break;
+    case 'f':
+      mol->OutputTemperatureFromTrajectories((ofstream *)&cout, mol->MDSteps-1, mol->MDSteps, (ofstream *)&cout);
+      break;
+    case 'g':
+      {
+        char filename[255];
+        cout << "Please enter filename: " << endl;
+        cin >> filename;
+        cout << Verbose(1) << "Storing temperatures in " << filename << "." << endl;
+        ofstream *output = new ofstream(filename, ios::trunc);
+        if (!mol->OutputTemperatureFromTrajectories((ofstream *)&cout, 0, mol->MDSteps, output))
+          cout << Verbose(2) << "File could not be written." << endl;
+        else
+          cout << Verbose(2) << "File stored." << endl;
+        output->close();
+        delete(output);
+      }
+      break;
+  }
 };
 
@@ -611,17 +611,17 @@
 static void FragmentAtoms(molecule *mol, config *configuration)
 {
-	int Order1;
-	clock_t start, end;
-
-	cout << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;
-	cout << Verbose(0) << "What's the desired bond order: ";
-	cin >> Order1;
-	if (mol->first->next != mol->last) {	// there are bonds
-		start = clock();
-		mol->FragmentMolecule((ofstream *)&cout, Order1, configuration);
-		end = clock();
-		cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
-	} else
-		cout << Verbose(0) << "Connection matrix has not yet been generated!" << endl;
+  int Order1;
+  clock_t start, end;
+
+  cout << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;
+  cout << Verbose(0) << "What's the desired bond order: ";
+  cin >> Order1;
+  if (mol->first->next != mol->last) {  // there are bonds
+    start = clock();
+    mol->FragmentMolecule((ofstream *)&cout, Order1, configuration);
+    end = clock();
+    cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
+  } else
+    cout << Verbose(0) << "Connection matrix has not yet been generated!" << endl;
 };
 
@@ -1053,6 +1053,6 @@
 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
+  // the current test routine checks the functionality of the KeySet&Graph concept:
+  // We want to have a multiindex (the KeySet) describing a unique subgraph
   int i, comp, counter=0;
 
@@ -1067,63 +1067,63 @@
   atom *Walker = mol->start;
 
-	// generate some KeySets
-	cout << "Generating KeySets." << endl;
-	KeySet TestSets[mol->AtomCount+1];
-	i=1;
-	while (Walker->next != mol->end) {
-		Walker = Walker->next;
-		for (int j=0;j<i;j++) {
-			TestSets[j].insert(Walker->nr);
-		}
-		i++;
-	}
-	cout << "Testing insertion of already present item in KeySets." << endl;
-	KeySetTestPair test;
-	test = TestSets[mol->AtomCount-1].insert(Walker->nr);
-	if (test.second) {
-		cout << Verbose(1) << "Insertion worked?!" << endl;
-	} else {
-		cout << Verbose(1) << "Insertion rejected: Present object is " << (*test.first) << "." << endl;
-	}
-	TestSets[mol->AtomCount].insert(mol->end->previous->nr);
-	TestSets[mol->AtomCount].insert(mol->end->previous->previous->previous->nr);
-
-	// constructing Graph structure
-	cout << "Generating Subgraph class." << endl;
-	Graph Subgraphs;
-
-	// insert KeySets into Subgraphs
-	cout << "Inserting KeySets into Subgraph class." << endl;
-	for (int j=0;j<mol->AtomCount;j++) {
-		Subgraphs.insert(GraphPair (TestSets[j],pair<int, double>(counter++, 1.)));
-	}
-	cout << "Testing insertion of already present item in Subgraph." << endl;
-	GraphTestPair test2;
-	test2 = Subgraphs.insert(GraphPair (TestSets[mol->AtomCount],pair<int, double>(counter++, 1.)));
-	if (test2.second) {
-		cout << Verbose(1) << "Insertion worked?!" << endl;
-	} else {
-		cout << Verbose(1) << "Insertion rejected: Present object is " << (*(test2.first)).second.first << "." << endl;
-	}
-
-	// show graphs
-	cout << "Showing Subgraph's contents, checking that it's sorted." << endl;
-	Graph::iterator A = Subgraphs.begin();
-	while (A !=	Subgraphs.end()) {
-		cout << (*A).second.first << ": ";
-		KeySet::iterator key = (*A).first.begin();
-		comp = -1;
-		while (key != (*A).first.end()) {
-			if ((*key) > comp)
-				cout << (*key) << " ";
-			else
-				cout << (*key) << "! ";
-			comp = (*key);
-			key++;
-		}
-		cout << endl;
-		A++;
-	}
-	delete(mol);
+  // generate some KeySets
+  cout << "Generating KeySets." << endl;
+  KeySet TestSets[mol->AtomCount+1];
+  i=1;
+  while (Walker->next != mol->end) {
+    Walker = Walker->next;
+    for (int j=0;j<i;j++) {
+      TestSets[j].insert(Walker->nr);
+    }
+    i++;
+  }
+  cout << "Testing insertion of already present item in KeySets." << endl;
+  KeySetTestPair test;
+  test = TestSets[mol->AtomCount-1].insert(Walker->nr);
+  if (test.second) {
+    cout << Verbose(1) << "Insertion worked?!" << endl;
+  } else {
+    cout << Verbose(1) << "Insertion rejected: Present object is " << (*test.first) << "." << endl;
+  }
+  TestSets[mol->AtomCount].insert(mol->end->previous->nr);
+  TestSets[mol->AtomCount].insert(mol->end->previous->previous->previous->nr);
+
+  // constructing Graph structure
+  cout << "Generating Subgraph class." << endl;
+  Graph Subgraphs;
+
+  // insert KeySets into Subgraphs
+  cout << "Inserting KeySets into Subgraph class." << endl;
+  for (int j=0;j<mol->AtomCount;j++) {
+    Subgraphs.insert(GraphPair (TestSets[j],pair<int, double>(counter++, 1.)));
+  }
+  cout << "Testing insertion of already present item in Subgraph." << endl;
+  GraphTestPair test2;
+  test2 = Subgraphs.insert(GraphPair (TestSets[mol->AtomCount],pair<int, double>(counter++, 1.)));
+  if (test2.second) {
+    cout << Verbose(1) << "Insertion worked?!" << endl;
+  } else {
+    cout << Verbose(1) << "Insertion rejected: Present object is " << (*(test2.first)).second.first << "." << endl;
+  }
+
+  // show graphs
+  cout << "Showing Subgraph's contents, checking that it's sorted." << endl;
+  Graph::iterator A = Subgraphs.begin();
+  while (A !=  Subgraphs.end()) {
+    cout << (*A).second.first << ": ";
+    KeySet::iterator key = (*A).first.begin();
+    comp = -1;
+    while (key != (*A).first.end()) {
+      if ((*key) > comp)
+        cout << (*key) << " ";
+      else
+        cout << (*key) << "! ";
+      comp = (*key);
+      key++;
+    }
+    cout << endl;
+    A++;
+  }
+  delete(mol);
 };
 
@@ -1136,79 +1136,79 @@
 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;
-	// get correct valence orbitals
-	mol->CalculateOrbitals(*configuration);
-	configuration->InitMaxMinStopStep = configuration->MaxMinStopStep = configuration->MaxPsiDouble;
-	strcpy(filename, ConfigFileName);
-	if (ConfigFileName != NULL) { // test the file name
-		output.open(ConfigFileName, ios::trunc);
-	} else if (strlen(configuration->configname) != 0) {
-		strcpy(filename, configuration->configname);
-		output.open(configuration->configname, ios::trunc);
-		} else {
-			strcpy(filename, DEFAULTCONFIG);
-			output.open(DEFAULTCONFIG, ios::trunc);
-		}
-	output.close();
-	output.clear();
-	cout << Verbose(0) << "Saving of config file ";
-	if (configuration->Save(filename, periode, mol))
-		cout << "successful." << endl;
-	else
-		cout << "failed." << endl;
-
-	// and save to xyz file
-	if (ConfigFileName != NULL) {
-		strcpy(filename, ConfigFileName);
-		strcat(filename, ".xyz");
-		output.open(filename, ios::trunc);
-	}
-	if (output == NULL) {
-		strcpy(filename,"main_pcp_linux");
-		strcat(filename, ".xyz");
-		output.open(filename, ios::trunc);
-	}
-	cout << Verbose(0) << "Saving of XYZ file ";
-	if (mol->MDSteps <= 1) {
-		if (mol->OutputXYZ(&output))
-			cout << "successful." << endl;
-		else
-			cout << "failed." << endl;
-	} else {
-		if (mol->OutputTrajectoriesXYZ(&output))
-			cout << "successful." << endl;
-		else
-			cout << "failed." << endl;
-	}
-	output.close();
-	output.clear();
-
-	// and save as MPQC configuration
-	if (ConfigFileName != NULL)
-		strcpy(filename, ConfigFileName);
-	if (output == NULL)
-		strcpy(filename,"main_pcp_linux");
-	cout << Verbose(0) << "Saving as mpqc input ";
-	if (configuration->SaveMPQC(filename, mol))
-		cout << "done." << endl;
-	else
-		cout << "failed." << endl;
-
-	if (!strcmp(configuration->configpath, configuration->GetDefaultPath())) {
-		cerr << "WARNING: config is found under different path then stated in config file::defaultpath!" << endl;
-	}
-	delete(mol);
+  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;
+  // get correct valence orbitals
+  mol->CalculateOrbitals(*configuration);
+  configuration->InitMaxMinStopStep = configuration->MaxMinStopStep = configuration->MaxPsiDouble;
+  strcpy(filename, ConfigFileName);
+  if (ConfigFileName != NULL) { // test the file name
+    output.open(ConfigFileName, ios::trunc);
+  } else if (strlen(configuration->configname) != 0) {
+    strcpy(filename, configuration->configname);
+    output.open(configuration->configname, ios::trunc);
+    } else {
+      strcpy(filename, DEFAULTCONFIG);
+      output.open(DEFAULTCONFIG, ios::trunc);
+    }
+  output.close();
+  output.clear();
+  cout << Verbose(0) << "Saving of config file ";
+  if (configuration->Save(filename, periode, mol))
+    cout << "successful." << endl;
+  else
+    cout << "failed." << endl;
+
+  // and save to xyz file
+  if (ConfigFileName != NULL) {
+    strcpy(filename, ConfigFileName);
+    strcat(filename, ".xyz");
+    output.open(filename, ios::trunc);
+  }
+  if (output == NULL) {
+    strcpy(filename,"main_pcp_linux");
+    strcat(filename, ".xyz");
+    output.open(filename, ios::trunc);
+  }
+  cout << Verbose(0) << "Saving of XYZ file ";
+  if (mol->MDSteps <= 1) {
+    if (mol->OutputXYZ(&output))
+      cout << "successful." << endl;
+    else
+      cout << "failed." << endl;
+  } else {
+    if (mol->OutputTrajectoriesXYZ(&output))
+      cout << "successful." << endl;
+    else
+      cout << "failed." << endl;
+  }
+  output.close();
+  output.clear();
+
+  // and save as MPQC configuration
+  if (ConfigFileName != NULL)
+    strcpy(filename, ConfigFileName);
+  if (output == NULL)
+    strcpy(filename,"main_pcp_linux");
+  cout << Verbose(0) << "Saving as mpqc input ";
+  if (configuration->SaveMPQC(filename, mol))
+    cout << "done." << endl;
+  else
+    cout << "failed." << endl;
+
+  if (!strcmp(configuration->configpath, configuration->GetDefaultPath())) {
+    cerr << "WARNING: config is found under different path then stated in config file::defaultpath!" << endl;
+  }
+  delete(mol);
 };
 
@@ -1225,315 +1225,315 @@
 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
-	double *factor; // unit factor if desired
-	ifstream test;
-	ofstream output;
-	string line;
-	atom *first;
-	bool SaveFlag = false;
-	int ExitFlag = 0;
-	int j;
-	double volume = 0.;
-	enum ConfigStatus config_present = absent;
-	clock_t start,end;
-	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
-		// 1. : Parse options that just set variables or print help
-		argptr = 1;
-		do {
-			if (argv[argptr][0] == '-') {
-				cout << Verbose(0) << "Recognized command line argument: " << argv[argptr][1] << ".\n";
-				argptr++;
-				switch(argv[argptr-1][1]) {
-					case 'h':
-					case 'H':
-					case '?':
-						cout << "MoleCuilder suite" << endl << "==================" << endl << endl;
-						cout << "Usage: " << argv[0] << "[config file] [-{acefpsthH?vfrp}] [further arguments]" << endl;
-						cout << "or simply " << argv[0] << " without arguments for interactive session." << endl;
-						cout << "\t-a Z x1 x2 x3\tAdd new atom of element Z at coordinates (x1,x2,x3)." << endl;
-						cout << "\t-A <source>\tCreate adjacency list from bonds parsed from 'dbond'-style file." <<endl;
-						cout << "\t-b x1 x2 x3\tCenter atoms in domain with given edge lengths of (x1,x2,x3)." << endl;
-						cout << "\t-B <basis>\tSetting basis to store to MPQC config files." << endl;
-						cout << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl;
-						cout << "\t-D <bond distance>\tDepth-First-Search Analysis of the molecule, giving cycles and tree/back edges." << endl;
-						cout << "\t-O\tCenter atoms in origin." << endl;
-						cout << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl;
-						cout << "\t-e <file>\tSets the databases path to be parsed (default: ./)." << endl;
-						cout << "\t-E <id> <Z>\tChange atom <id>'s element to <Z>, <id> begins at 0." << endl;
-						cout << "\t-f/F <dist> <order>\tFragments the molecule in BOSSANOVA manner (with/out rings compressed) and stores config files in same dir as config (return code 0 - fragmented, 2 - no fragmentation necessary)." << endl;
-						cout << "\t-h/-H/-?\tGive this help screen." << endl;
-						cout << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl;
-						cout << "\t-n\tFast parsing (i.e. no trajectories are looked for)." << endl;
-						cout << "\t-N <radius> <file>\tGet non-convex-envelope." << endl;
-						cout << "\t-o <out>\tGet volume of the convex envelope (and store to tecplot file)." << endl;
-						cout << "\t-p <file>\tParse given xyz file and create raw config file from it." << endl;
-						cout << "\t-P <file>\tParse given forces file and append as an MD step to config file via Verlet." << endl;
-						cout << "\t-r\t\tConvert file from an old pcp syntax." << endl;
-						cout << "\t-R\t\tRemove all atoms out of sphere around a given one." << endl;
-						cout << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << endl;
-						cout << "\t-T <file> Store temperatures from the config file in <file>." << endl;
-						cout << "\t-s x1 x2 x3\tScale all atom coordinates by this vector (x1,x2,x3)." << endl;
-						cout << "\t-u rho\tsuspend in water solution and output necessary cell lengths, average density rho and repetition." << endl;
-						cout << "\t-v/-V\t\tGives version information." << endl;
-						cout << "Note: config files must not begin with '-' !" << endl;
-						delete(mol);
-						delete(periode);
-						return (1);
-						break;
-					case 'v':
-					case 'V':
-						cout << argv[0] << " " << VERSIONSTRING << endl;
-						cout << "Build your own molecule position set." << endl;
-						delete(mol);
-						delete(periode);
-						return (1);
-						break;
-					case 'e':
-						if ((argptr >= argc) || (argv[argptr][0] == '-')) {
-							cerr << "Not enough or invalid arguments for specifying element db: -e <db file>" << endl;
-						} else {
-							cout << "Using " << argv[argptr] << " as elements database." << endl;
-							PathToDatabases = argv[argptr];
-							argptr+=1;
-						}
-						break;
-					case 'n':
-						cout << "I won't parse trajectories." << endl;
-						configuration.FastParsing = true;
-						break;
-					default:	 // no match? Step on
-						argptr++;
-						break;
-				}
-			} else
-				argptr++;
-		} while (argptr < argc);
-
-		// 2. Parse the element database
-		if (periode->LoadPeriodentafel(PathToDatabases)) {
-			cout << Verbose(0) << "Element list loaded successfully." << endl;
-			//periode->Output((ofstream *)&cout);
-		} else {
-			cout << Verbose(0) << "Element list loading failed." << endl;
-			return 1;
-		}
-		// 3. Find config file name and parse if possible
-		if (argv[1][0] != '-') {
-			cout << Verbose(0) << "Config file given." << endl;
-			test.open(argv[1], ios::in);
-			if (test == NULL) {
-				//return (1);
-				output.open(argv[1], ios::out);
-				if (output == NULL) {
-					cout << Verbose(1) << "Specified config file " << argv[1] << " not found." << endl;
-					config_present = absent;
-				} else {
-					cout << "Empty configuration file." << endl;
-					ConfigFileName = argv[1];
-					config_present = empty;
-					output.close();
-				}
-			} else {
-				test.close();
-				ConfigFileName = argv[1];
-				cout << Verbose(1) << "Specified config file found, parsing ... ";
-				switch (configuration.TestSyntax(ConfigFileName, periode, mol)) {
-					case 1:
-						cout << "new syntax." << endl;
-						configuration.Load(ConfigFileName, periode, mol);
-						config_present = present;
-						break;
-					case 0:
-						cout << "old syntax." << endl;
-						configuration.LoadOld(ConfigFileName, periode, mol);
-						config_present = present;
-						break;
-					default:
-						cout << "Unknown syntax or empty, yet present file." << endl;
-						config_present = empty;
-			 }
-			}
-		} else
-			config_present = absent;
-		// 4. parse again through options, now for those depending on elements db and config presence
-		argptr = 1;
-		do {
-			cout << "Current Command line argument: " << argv[argptr] << "." << endl;
-			if (argv[argptr][0] == '-') {
-				argptr++;
-				if ((config_present == present) || (config_present == empty)) {
-					switch(argv[argptr-1][1]) {
-						case 'p':
-							ExitFlag = 1;
-							if ((argptr >= argc) || (argv[argptr][0] == '-')) {
-								ExitFlag = 255;
-								cerr << "Not enough arguments for parsing: -p <xyz file>" << endl;
-							} else {
-								SaveFlag = true;
-								cout << Verbose(1) << "Parsing xyz file for new atoms." << endl;
-								if (!mol->AddXYZFile(argv[argptr]))
-									cout << Verbose(2) << "File not found." << endl;
-								else {
-									cout << Verbose(2) << "File found and parsed." << endl;
-									config_present = present;
-								}
-							}
-							break;
-						case 'a':
-							ExitFlag = 1;
-							if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr+1]))) {
-								ExitFlag = 255;
-								cerr << "Not enough or invalid arguments for adding atom: -a <element> <x> <y> <z>" << endl;
-							} else {
-								SaveFlag = true;
-								cout << Verbose(1) << "Adding new atom with element " << argv[argptr] << " at (" << argv[argptr+1] << "," << argv[argptr+2] << "," << argv[argptr+3] << "), ";
-								first = new atom;
-								first->type = periode->FindElement(atoi(argv[argptr]));
-								if (first->type != NULL)
-									cout << Verbose(2) << "found element " << first->type->name << endl;
-								for (int i=NDIM;i--;)
-									first->x.x[i] = atof(argv[argptr+1+i]);
-								if (first->type != NULL) {
-									mol->AddAtom(first);	// add to molecule
-									if ((config_present == empty) && (mol->AtomCount != 0))
-										config_present = present;
-								} else
-									cerr << Verbose(1) << "Could not find the specified element." << endl;
-								argptr+=4;
-							}
-							break;
-						default:	 // no match? Don't step on (this is done in next switch's default)
-							break;
-					}
-				}
-				if (config_present == present) {
-					switch(argv[argptr-1][1]) {
-						case 'B':
-							if ((argptr >= argc) || (argv[argptr][0] == '-')) {
-								ExitFlag = 255;
-								cerr << "Not enough or invalid arguments given for setting MPQC basis: -B <basis name>" << endl;
-							} else {
-								configuration.basis = argv[argptr];
-								cout << Verbose(1) << "Setting MPQC basis to " << configuration.basis << "." << endl;
-								argptr+=1;
-							}
-							break;
-						case 'D':
-							ExitFlag = 1;
-							{
-								cout << Verbose(1) << "Depth-First-Search Analysis." << endl;
-								MoleculeLeafClass *Subgraphs = NULL;			// list of subgraphs from DFS analysis
-								int *MinimumRingSize = new int[mol->AtomCount];
-								atom ***ListOfLocalAtoms = NULL;
-								int FragmentCounter = 0;
-								class StackClass<bond *> *BackEdgeStack = NULL;
-								class StackClass<bond *> *LocalBackEdgeStack = NULL;
-								mol->CreateAdjacencyList((ofstream *)&cout, atof(argv[argptr]), configuration.GetIsAngstroem());
-								mol->CreateListOfBondsPerAtom((ofstream *)&cout);
-								Subgraphs = mol->DepthFirstSearchAnalysis((ofstream *)&cout, BackEdgeStack);
-								if (Subgraphs != NULL) {
-									Subgraphs->next->FillBondStructureFromReference((ofstream *)&cout, mol, (FragmentCounter = 0), ListOfLocalAtoms, false);	// we want to keep the created ListOfLocalAtoms
-									while (Subgraphs->next != NULL) {
-										Subgraphs = Subgraphs->next;
-										LocalBackEdgeStack = new StackClass<bond *> (Subgraphs->Leaf->BondCount);
-										Subgraphs->Leaf->PickLocalBackEdges((ofstream *)&cout, ListOfLocalAtoms[FragmentCounter++], BackEdgeStack, LocalBackEdgeStack);
-										Subgraphs->Leaf->CyclicStructureAnalysis((ofstream *)&cout, BackEdgeStack, MinimumRingSize);
-										delete(LocalBackEdgeStack);
-										delete(Subgraphs->previous);
-									}
-									delete(Subgraphs);
-									for (int i=0;i<FragmentCounter;i++)
-										Free((void **)&ListOfLocalAtoms[FragmentCounter], "ParseCommandLineOptions: **ListOfLocalAtoms[]");
-									Free((void **)&ListOfLocalAtoms, "ParseCommandLineOptions: ***ListOfLocalAtoms");
-								}
-								delete(BackEdgeStack);
-								delete[](MinimumRingSize);
-							}
-							//argptr+=1;
-							break;
-						case 'E':
-							ExitFlag = 1;
-							if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (argv[argptr+1][0] == '-')) {
-								ExitFlag = 255;
-								cerr << "Not enough or invalid arguments given for changing element: -E <atom nr.> <element>" << endl;
-							} else {
-								SaveFlag = true;
-								cout << Verbose(1) << "Changing atom " << argv[argptr] << " to element " << argv[argptr+1] << "." << endl;
-								first = mol->FindAtom(atoi(argv[argptr]));
-								first->type = periode->FindElement(atoi(argv[argptr+1]));
-								argptr+=2;
-							}
-							break;
-						case 'A':
-							ExitFlag = 1;
-							if ((argptr >= argc) || (argv[argptr][0] == '-')) {
-								ExitFlag =255;
-								cerr << "Missing source file for bonds in molecule: -A <bond sourcefile>" << endl;
-							} else {
-								cout << "Parsing bonds from " << argv[argptr] << "." << endl;
-								ifstream *input = new ifstream(argv[argptr]);
-								mol->CreateAdjacencyList2((ofstream *)&cout, input);
-								input->close();
-								argptr+=1;
-							}
-							break;
-						case 'N':
-							ExitFlag = 1;
-							if ((argptr+1 >= argc) || (argv[argptr+1][0] == '-')){
-								ExitFlag = 255;
-								cerr << "Not enough or invalid arguments given for non-convex envelope: -o <radius> <tecplot output file>" << endl;
-							} else {
-								class Tesselation T;
-								int N = 15;
-								int number = 100;
-								string filename(argv[argptr+1]);
-								filename.append(".csv");
-								cout << Verbose(0) << "Evaluating non-convex envelope.";
-								cout << Verbose(1) << "Using rolling ball of radius " << atof(argv[argptr]) << " and storing tecplot data in " << argv[argptr+1] << "." << endl;
-								LinkedCell LCList(mol, atof(argv[argptr]));		// \NOTE not twice the radius??
-								Find_non_convex_border((ofstream *)&cout, mol, &T, &LCList, argv[argptr+1], atof(argv[argptr]));
-								//FindDistributionOfEllipsoids((ofstream *)&cout, &T, &LCList, N, number, filename.c_str());
-								argptr+=2;
-							}
-							break;
-						case 'T':
-							ExitFlag = 1;
-							if ((argptr >= argc) || (argv[argptr][0] == '-')) {
-								ExitFlag = 255;
-								cerr << "Not enough or invalid arguments given for storing tempature: -T <temperature file>" << endl;
-							} else {
-								cout << Verbose(1) << "Storing temperatures in " << argv[argptr] << "." << endl;
-								ofstream *output = new ofstream(argv[argptr], ios::trunc);
-								if (!mol->OutputTemperatureFromTrajectories((ofstream *)&cout, 0, mol->MDSteps, output))
-									cout << Verbose(2) << "File could not be written." << endl;
-								else
-									cout << Verbose(2) << "File stored." << endl;
-								output->close();
-								delete(output);
-								argptr+=1;
-							}
-							break;
-						case 'P':
-							ExitFlag = 1;
-							if ((argptr >= argc) || (argv[argptr][0] == '-')) {
-								ExitFlag = 255;
-								cerr << "Not enough or invalid arguments given for parsing and integrating forces: -P <forces file>" << endl;
-							} else {
-								SaveFlag = true;
-								cout << Verbose(1) << "Parsing forces file and Verlet integrating." << endl;
-								if (!mol->VerletForceIntegration(argv[argptr], configuration.Deltat, configuration.GetIsAngstroem()))
-									cout << Verbose(2) << "File not found." << endl;
-								else
-									cout << Verbose(2) << "File found and parsed." << endl;
-								argptr+=1;
-							}
-							break;
+  Vector x,y,z,n;  // coordinates for absolute point in cell volume
+  double *factor; // unit factor if desired
+  ifstream test;
+  ofstream output;
+  string line;
+  atom *first;
+  bool SaveFlag = false;
+  int ExitFlag = 0;
+  int j;
+  double volume = 0.;
+  enum ConfigStatus config_present = absent;
+  clock_t start,end;
+  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
+    // 1. : Parse options that just set variables or print help
+    argptr = 1;
+    do {
+      if (argv[argptr][0] == '-') {
+        cout << Verbose(0) << "Recognized command line argument: " << argv[argptr][1] << ".\n";
+        argptr++;
+        switch(argv[argptr-1][1]) {
+          case 'h':
+          case 'H':
+          case '?':
+            cout << "MoleCuilder suite" << endl << "==================" << endl << endl;
+            cout << "Usage: " << argv[0] << "[config file] [-{acefpsthH?vfrp}] [further arguments]" << endl;
+            cout << "or simply " << argv[0] << " without arguments for interactive session." << endl;
+            cout << "\t-a Z x1 x2 x3\tAdd new atom of element Z at coordinates (x1,x2,x3)." << endl;
+            cout << "\t-A <source>\tCreate adjacency list from bonds parsed from 'dbond'-style file." <<endl;
+            cout << "\t-b x1 x2 x3\tCenter atoms in domain with given edge lengths of (x1,x2,x3)." << endl;
+            cout << "\t-B <basis>\tSetting basis to store to MPQC config files." << endl;
+            cout << "\t-c x1 x2 x3\tCenter atoms in domain with a minimum distance to boundary of (x1,x2,x3)." << endl;
+            cout << "\t-D <bond distance>\tDepth-First-Search Analysis of the molecule, giving cycles and tree/back edges." << endl;
+            cout << "\t-O\tCenter atoms in origin." << endl;
+            cout << "\t-d x1 x2 x3\tDuplicate cell along each axis by given factor." << endl;
+            cout << "\t-e <file>\tSets the databases path to be parsed (default: ./)." << endl;
+            cout << "\t-E <id> <Z>\tChange atom <id>'s element to <Z>, <id> begins at 0." << endl;
+            cout << "\t-f/F <dist> <order>\tFragments the molecule in BOSSANOVA manner (with/out rings compressed) and stores config files in same dir as config (return code 0 - fragmented, 2 - no fragmentation necessary)." << endl;
+            cout << "\t-h/-H/-?\tGive this help screen." << endl;
+            cout << "\t-m <0/1>\tCalculate (0)/ Align in(1) PAS with greatest EV along z axis." << endl;
+            cout << "\t-n\tFast parsing (i.e. no trajectories are looked for)." << endl;
+            cout << "\t-N <radius> <file>\tGet non-convex-envelope." << endl;
+            cout << "\t-o <out>\tGet volume of the convex envelope (and store to tecplot file)." << endl;
+            cout << "\t-p <file>\tParse given xyz file and create raw config file from it." << endl;
+            cout << "\t-P <file>\tParse given forces file and append as an MD step to config file via Verlet." << endl;
+            cout << "\t-r\t\tConvert file from an old pcp syntax." << endl;
+            cout << "\t-R\t\tRemove all atoms out of sphere around a given one." << endl;
+            cout << "\t-t x1 x2 x3\tTranslate all atoms by this vector (x1,x2,x3)." << endl;
+            cout << "\t-T <file> Store temperatures from the config file in <file>." << endl;
+            cout << "\t-s x1 x2 x3\tScale all atom coordinates by this vector (x1,x2,x3)." << endl;
+            cout << "\t-u rho\tsuspend in water solution and output necessary cell lengths, average density rho and repetition." << endl;
+            cout << "\t-v/-V\t\tGives version information." << endl;
+            cout << "Note: config files must not begin with '-' !" << endl;
+            delete(mol);
+            delete(periode);
+            return (1);
+            break;
+          case 'v':
+          case 'V':
+            cout << argv[0] << " " << VERSIONSTRING << endl;
+            cout << "Build your own molecule position set." << endl;
+            delete(mol);
+            delete(periode);
+            return (1);
+            break;
+          case 'e':
+            if ((argptr >= argc) || (argv[argptr][0] == '-')) {
+              cerr << "Not enough or invalid arguments for specifying element db: -e <db file>" << endl;
+            } else {
+              cout << "Using " << argv[argptr] << " as elements database." << endl;
+              PathToDatabases = argv[argptr];
+              argptr+=1;
+            }
+            break;
+          case 'n':
+            cout << "I won't parse trajectories." << endl;
+            configuration.FastParsing = true;
+            break;
+          default:   // no match? Step on
+            argptr++;
+            break;
+        }
+      } else
+        argptr++;
+    } while (argptr < argc);
+
+    // 2. Parse the element database
+    if (periode->LoadPeriodentafel(PathToDatabases)) {
+      cout << Verbose(0) << "Element list loaded successfully." << endl;
+      //periode->Output((ofstream *)&cout);
+    } else {
+      cout << Verbose(0) << "Element list loading failed." << endl;
+      return 1;
+    }
+    // 3. Find config file name and parse if possible
+    if (argv[1][0] != '-') {
+      cout << Verbose(0) << "Config file given." << endl;
+      test.open(argv[1], ios::in);
+      if (test == NULL) {
+        //return (1);
+        output.open(argv[1], ios::out);
+        if (output == NULL) {
+          cout << Verbose(1) << "Specified config file " << argv[1] << " not found." << endl;
+          config_present = absent;
+        } else {
+          cout << "Empty configuration file." << endl;
+          ConfigFileName = argv[1];
+          config_present = empty;
+          output.close();
+        }
+      } else {
+        test.close();
+        ConfigFileName = argv[1];
+        cout << Verbose(1) << "Specified config file found, parsing ... ";
+        switch (configuration.TestSyntax(ConfigFileName, periode, mol)) {
+          case 1:
+            cout << "new syntax." << endl;
+            configuration.Load(ConfigFileName, periode, mol);
+            config_present = present;
+            break;
+          case 0:
+            cout << "old syntax." << endl;
+            configuration.LoadOld(ConfigFileName, periode, mol);
+            config_present = present;
+            break;
+          default:
+            cout << "Unknown syntax or empty, yet present file." << endl;
+            config_present = empty;
+       }
+      }
+    } else
+      config_present = absent;
+    // 4. parse again through options, now for those depending on elements db and config presence
+    argptr = 1;
+    do {
+      cout << "Current Command line argument: " << argv[argptr] << "." << endl;
+      if (argv[argptr][0] == '-') {
+        argptr++;
+        if ((config_present == present) || (config_present == empty)) {
+          switch(argv[argptr-1][1]) {
+            case 'p':
+              ExitFlag = 1;
+              if ((argptr >= argc) || (argv[argptr][0] == '-')) {
+                ExitFlag = 255;
+                cerr << "Not enough arguments for parsing: -p <xyz file>" << endl;
+              } else {
+                SaveFlag = true;
+                cout << Verbose(1) << "Parsing xyz file for new atoms." << endl;
+                if (!mol->AddXYZFile(argv[argptr]))
+                  cout << Verbose(2) << "File not found." << endl;
+                else {
+                  cout << Verbose(2) << "File found and parsed." << endl;
+                  config_present = present;
+                }
+              }
+              break;
+            case 'a':
+              ExitFlag = 1;
+              if ((argptr >= argc) || (argv[argptr][0] == '-') || (!IsValidNumber(argv[argptr+1]))) {
+                ExitFlag = 255;
+                cerr << "Not enough or invalid arguments for adding atom: -a <element> <x> <y> <z>" << endl;
+              } else {
+                SaveFlag = true;
+                cout << Verbose(1) << "Adding new atom with element " << argv[argptr] << " at (" << argv[argptr+1] << "," << argv[argptr+2] << "," << argv[argptr+3] << "), ";
+                first = new atom;
+                first->type = periode->FindElement(atoi(argv[argptr]));
+                if (first->type != NULL)
+                  cout << Verbose(2) << "found element " << first->type->name << endl;
+                for (int i=NDIM;i--;)
+                  first->x.x[i] = atof(argv[argptr+1+i]);
+                if (first->type != NULL) {
+                  mol->AddAtom(first);  // add to molecule
+                  if ((config_present == empty) && (mol->AtomCount != 0))
+                    config_present = present;
+                } else
+                  cerr << Verbose(1) << "Could not find the specified element." << endl;
+                argptr+=4;
+              }
+              break;
+            default:   // no match? Don't step on (this is done in next switch's default)
+              break;
+          }
+        }
+        if (config_present == present) {
+          switch(argv[argptr-1][1]) {
+            case 'B':
+              if ((argptr >= argc) || (argv[argptr][0] == '-')) {
+                ExitFlag = 255;
+                cerr << "Not enough or invalid arguments given for setting MPQC basis: -B <basis name>" << endl;
+              } else {
+                configuration.basis = argv[argptr];
+                cout << Verbose(1) << "Setting MPQC basis to " << configuration.basis << "." << endl;
+                argptr+=1;
+              }
+              break;
+            case 'D':
+              ExitFlag = 1;
+              {
+                cout << Verbose(1) << "Depth-First-Search Analysis." << endl;
+                MoleculeLeafClass *Subgraphs = NULL;      // list of subgraphs from DFS analysis
+                int *MinimumRingSize = new int[mol->AtomCount];
+                atom ***ListOfLocalAtoms = NULL;
+                int FragmentCounter = 0;
+                class StackClass<bond *> *BackEdgeStack = NULL;
+                class StackClass<bond *> *LocalBackEdgeStack = NULL;
+                mol->CreateAdjacencyList((ofstream *)&cout, atof(argv[argptr]), configuration.GetIsAngstroem());
+                mol->CreateListOfBondsPerAtom((ofstream *)&cout);
+                Subgraphs = mol->DepthFirstSearchAnalysis((ofstream *)&cout, BackEdgeStack);
+                if (Subgraphs != NULL) {
+                  Subgraphs->next->FillBondStructureFromReference((ofstream *)&cout, mol, (FragmentCounter = 0), ListOfLocalAtoms, false);  // we want to keep the created ListOfLocalAtoms
+                  while (Subgraphs->next != NULL) {
+                    Subgraphs = Subgraphs->next;
+                    LocalBackEdgeStack = new StackClass<bond *> (Subgraphs->Leaf->BondCount);
+                    Subgraphs->Leaf->PickLocalBackEdges((ofstream *)&cout, ListOfLocalAtoms[FragmentCounter++], BackEdgeStack, LocalBackEdgeStack);
+                    Subgraphs->Leaf->CyclicStructureAnalysis((ofstream *)&cout, BackEdgeStack, MinimumRingSize);
+                    delete(LocalBackEdgeStack);
+                    delete(Subgraphs->previous);
+                  }
+                  delete(Subgraphs);
+                  for (int i=0;i<FragmentCounter;i++)
+                    Free((void **)&ListOfLocalAtoms[FragmentCounter], "ParseCommandLineOptions: **ListOfLocalAtoms[]");
+                  Free((void **)&ListOfLocalAtoms, "ParseCommandLineOptions: ***ListOfLocalAtoms");
+                }
+                delete(BackEdgeStack);
+                delete[](MinimumRingSize);
+              }
+              //argptr+=1;
+              break;
+            case 'E':
+              ExitFlag = 1;
+              if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (argv[argptr+1][0] == '-')) {
+                ExitFlag = 255;
+                cerr << "Not enough or invalid arguments given for changing element: -E <atom nr.> <element>" << endl;
+              } else {
+                SaveFlag = true;
+                cout << Verbose(1) << "Changing atom " << argv[argptr] << " to element " << argv[argptr+1] << "." << endl;
+                first = mol->FindAtom(atoi(argv[argptr]));
+                first->type = periode->FindElement(atoi(argv[argptr+1]));
+                argptr+=2;
+              }
+              break;
+            case 'A':
+              ExitFlag = 1;
+              if ((argptr >= argc) || (argv[argptr][0] == '-')) {
+                ExitFlag =255;
+                cerr << "Missing source file for bonds in molecule: -A <bond sourcefile>" << endl;
+              } else {
+                cout << "Parsing bonds from " << argv[argptr] << "." << endl;
+                ifstream *input = new ifstream(argv[argptr]);
+                mol->CreateAdjacencyList2((ofstream *)&cout, input);
+                input->close();
+                argptr+=1;
+              }
+              break;
+            case 'N':
+              ExitFlag = 1;
+              if ((argptr+1 >= argc) || (argv[argptr+1][0] == '-')){
+                ExitFlag = 255;
+                cerr << "Not enough or invalid arguments given for non-convex envelope: -o <radius> <tecplot output file>" << endl;
+              } else {
+                class Tesselation T;
+                int N = 15;
+                int number = 100;
+                string filename(argv[argptr+1]);
+                filename.append(".csv");
+                cout << Verbose(0) << "Evaluating non-convex envelope.";
+                cout << Verbose(1) << "Using rolling ball of radius " << atof(argv[argptr]) << " and storing tecplot data in " << argv[argptr+1] << "." << endl;
+                LinkedCell LCList(mol, atof(argv[argptr]));    // \NOTE not twice the radius??
+                Find_non_convex_border((ofstream *)&cout, mol, &T, &LCList, argv[argptr+1], atof(argv[argptr]));
+                //FindDistributionOfEllipsoids((ofstream *)&cout, &T, &LCList, N, number, filename.c_str());
+                argptr+=2;
+              }
+              break;
+            case 'T':
+              ExitFlag = 1;
+              if ((argptr >= argc) || (argv[argptr][0] == '-')) {
+                ExitFlag = 255;
+                cerr << "Not enough or invalid arguments given for storing tempature: -T <temperature file>" << endl;
+              } else {
+                cout << Verbose(1) << "Storing temperatures in " << argv[argptr] << "." << endl;
+                ofstream *output = new ofstream(argv[argptr], ios::trunc);
+                if (!mol->OutputTemperatureFromTrajectories((ofstream *)&cout, 0, mol->MDSteps, output))
+                  cout << Verbose(2) << "File could not be written." << endl;
+                else
+                  cout << Verbose(2) << "File stored." << endl;
+                output->close();
+                delete(output);
+                argptr+=1;
+              }
+              break;
+            case 'P':
+              ExitFlag = 1;
+              if ((argptr >= argc) || (argv[argptr][0] == '-')) {
+                ExitFlag = 255;
+                cerr << "Not enough or invalid arguments given for parsing and integrating forces: -P <forces file>" << endl;
+              } else {
+                SaveFlag = true;
+                cout << Verbose(1) << "Parsing forces file and Verlet integrating." << endl;
+                if (!mol->VerletForceIntegration(argv[argptr], configuration.Deltat, configuration.GetIsAngstroem()))
+                  cout << Verbose(2) << "File not found." << endl;
+                else
+                  cout << Verbose(2) << "File found and parsed." << endl;
+                argptr+=1;
+              }
+              break;
             case 'R':
               ExitFlag = 1;
@@ -1561,264 +1561,264 @@
               }
               break;
-						case 't':
-							ExitFlag = 1;
-							if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
-								ExitFlag = 255;
-								cerr << "Not enough or invalid arguments given for translation: -t <x> <y> <z>" << endl;
-							} else {
-								ExitFlag = 1;
-								SaveFlag = true;
-								cout << Verbose(1) << "Translating all ions to new origin." << endl;
-								for (int i=NDIM;i--;)
-									x.x[i] = atof(argv[argptr+i]);
-								mol->Translate((const Vector *)&x);
-								argptr+=3;
-							}
-							break;
-						case 's':
-							ExitFlag = 1;
-							if ((argptr >= argc) || (!IsValidNumber(argv[argptr])) ) {
-								ExitFlag = 255;
-								cerr << "Not enough or invalid arguments given for scaling: -s <factor/[factor_x]> [factor_y] [factor_z]" << endl;
-							} else {
-								SaveFlag = true;
-								j = -1;
-								cout << Verbose(1) << "Scaling all ion positions by factor." << endl;
-								factor = new double[NDIM];
-								factor[0] = atof(argv[argptr]);
-								if ((argptr < argc) && (IsValidNumber(argv[argptr])))
-									argptr++;
-								factor[1] = atof(argv[argptr]);
-								if ((argptr < argc) && (IsValidNumber(argv[argptr])))
-									argptr++;
-								factor[2] = atof(argv[argptr]);
-								mol->Scale(&factor);
-								for (int i=0;i<NDIM;i++) {
-									j += i+1;
-									x.x[i] = atof(argv[NDIM+i]);
-									mol->cell_size[j]*=factor[i];
-								}
-								delete[](factor);
-								argptr+=1;
-							}
-							break;
-						case 'b':
-							ExitFlag = 1;
-							if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
-								ExitFlag = 255;
-								cerr << "Not enough or invalid arguments given for centering in box: -b <length_x> <length_y> <length_z>" << endl;
-							} else {
-								SaveFlag = true;
-								j = -1;
-								cout << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
-								j=-1;
-								for (int i=0;i<NDIM;i++) {
-									j += i+1;
-									x.x[i] = atof(argv[argptr++]);
-									mol->cell_size[j] += x.x[i]*2.;
-								}
-								// center
-								mol->CenterInBox((ofstream *)&cout, &x);
-								// update Box of atoms by boundary
-								mol->SetBoxDimension(&x);
-							}
-							break;
-						case 'c':
-							ExitFlag = 1;
-							if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
-								ExitFlag = 255;
-								cerr << "Not enough or invalid arguments given for centering with boundary: -c <boundary_x> <boundary_y> <boundary_z>" << endl;
-							} else {
-								SaveFlag = true;
-								j = -1;
-								cout << Verbose(1) << "Centering atoms in config file within given additional boundary." << endl;
-								// make every coordinate positive
-								mol->CenterEdge((ofstream *)&cout, &x);
-								// update Box of atoms by boundary
-								mol->SetBoxDimension(&x);
-								// translate each coordinate by boundary
-								j=-1;
-								for (int i=0;i<NDIM;i++) {
-									j += i+1;
-									x.x[i] = atof(argv[argptr++]);
-									mol->cell_size[j] += x.x[i]*2.;
-								}
-								mol->Translate((const Vector *)&x);
-							}
-							break;
-						case 'O':
-							ExitFlag = 1;
-							SaveFlag = true;
-							cout << Verbose(1) << "Centering atoms in origin." << endl;
-							mol->CenterOrigin((ofstream *)&cout, &x);
-							mol->SetBoxDimension(&x);
-							break;
-						case 'r':
-							ExitFlag = 1;
-							SaveFlag = true;
-							cout << Verbose(1) << "Converting config file from supposed old to new syntax." << endl;
-							break;
-						case 'F':
-						case 'f':
-							ExitFlag = 1;
-							if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1]))) {
-								ExitFlag = 255;
-								cerr << "Not enough or invalid arguments for fragmentation: -f <max. bond distance> <bond order>" << endl;
-							} else {
-								cout << "Fragmenting molecule with bond distance " << argv[argptr] << " angstroem, order of " << argv[argptr+1] << "." << endl;
-								cout << Verbose(0) << "Creating connection matrix..." << endl;
-								start = clock();
-								mol->CreateAdjacencyList((ofstream *)&cout, atof(argv[argptr++]), configuration.GetIsAngstroem());
-								cout << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;
-								if (mol->first->next != mol->last) {
-									ExitFlag = mol->FragmentMolecule((ofstream *)&cout, atoi(argv[argptr]), &configuration);
-								}
-								end = clock();
-								cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
-								argptr+=2;
-							}
-							break;
-						case 'm':
-							ExitFlag = 1;
-							j = atoi(argv[argptr++]);
-							if ((j<0) || (j>1)) {
-								cerr << Verbose(1) << "ERROR: Argument of '-m' should be either 0 for no-rotate or 1 for rotate." << endl;
-								j = 0;
-							}
-							if (j) {
-								SaveFlag = true;
-								cout << Verbose(0) << "Converting to prinicipal axis system." << endl;
-							} else
-								cout << Verbose(0) << "Evaluating prinicipal axis." << endl;
-							mol->PrincipalAxisSystem((ofstream *)&cout, (bool)j);
-							break;
-						case 'o':
-							ExitFlag = 1;
-							if ((argptr >= argc) || (argv[argptr][0] == '-')){
-								ExitFlag = 255;
-								cerr << "Not enough or invalid arguments given for convex envelope: -o <tecplot output file>" << endl;
-							} else {
-								cout << Verbose(0) << "Evaluating volume of the convex envelope.";
-								cout << Verbose(1) << "Storing tecplot data in " << argv[argptr] << "." << endl;
-								VolumeOfConvexEnvelope((ofstream *)&cout, argv[argptr], &configuration, NULL, mol);
-								argptr+=1;
-							}
-							break;
-						case 'U':
-							ExitFlag = 1;
-							if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) ) {
-								ExitFlag = 255;
-								cerr << "Not enough or invalid arguments given for suspension with specified volume: -U <volume> <density>" << endl;
-								volume = -1; // for case 'u': don't print error again
-							} else {
-								volume = atof(argv[argptr++]);
-								cout << Verbose(0) << "Using " << volume << " angstrom^3 as the volume instead of convex envelope one's." << endl;
-							}
-						case 'u':
-							ExitFlag = 1;
-							if ((argptr >= argc) || (!IsValidNumber(argv[argptr])) ) {
-								if (volume != -1)
-									ExitFlag = 255;
-									cerr << "Not enough arguments given for suspension: -u <density>" << endl;
-							} else {
-								double density;
-								SaveFlag = true;
-								cout << Verbose(0) << "Evaluating necessary cell volume for a cluster suspended in water.";
-								density = atof(argv[argptr++]);
-								if (density < 1.0) {
-									cerr << Verbose(0) << "Density must be greater than 1.0g/cm^3 !" << endl;
-									density = 1.3;
-								}
-//								for(int i=0;i<NDIM;i++) {
-//									repetition[i] = atoi(argv[argptr++]);
-//									if (repetition[i] < 1)
-//										cerr << Verbose(0) << "ERROR: repetition value must be greater 1!" << endl;
-//									repetition[i] = 1;
-//								}
-								PrepareClustersinWater((ofstream *)&cout, &configuration, mol, volume, density);	// if volume == 0, will calculate from ConvexEnvelope
-							}
-							break;
-						case 'd':
-							ExitFlag = 1;
-							if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
-								ExitFlag = 255;
-								cerr << "Not enough or invalid arguments given for repeating cells: -d <repeat_x> <repeat_y> <repeat_z>" << endl;
-							} else {
-								SaveFlag = true;
-								for (int axis = 1; axis <= NDIM; axis++) {
-									int faktor = atoi(argv[argptr++]);
-									int count;
-									element ** Elements;
-									Vector ** vectors;
-									if (faktor < 1) {
-										cerr << Verbose(0) << "ERROR: Repetition faktor mus be greater than 1!" << endl;
-										faktor = 1;
-									}
-									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, ...)
-											}
-										}
-										// 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;
-						default:	 // no match? Step on
-							if ((argptr < argc) && (argv[argptr][0] != '-')) // if it started with a '-' we've already made a step!
-								argptr++;
-							break;
-					}
-				}
-			} else argptr++;
-		} while (argptr < argc);
-		if (SaveFlag)
-			SaveConfig(ConfigFileName, &configuration, periode, molecules);
-		if ((ExitFlag >= 1)) {
-			delete(mol);
-			delete(periode);
-			return (ExitFlag);
-		}
-	} else {	// no arguments, hence scan the elements db
-		if (periode->LoadPeriodentafel(PathToDatabases))
-			cout << Verbose(0) << "Element list loaded successfully." << endl;
-		else
-			cout << Verbose(0) << "Element list loading failed." << endl;
-		configuration.RetrieveConfigPathAndName("main_pcp_linux");
-	}
-	return(0);
+            case 't':
+              ExitFlag = 1;
+              if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
+                ExitFlag = 255;
+                cerr << "Not enough or invalid arguments given for translation: -t <x> <y> <z>" << endl;
+              } else {
+                ExitFlag = 1;
+                SaveFlag = true;
+                cout << Verbose(1) << "Translating all ions to new origin." << endl;
+                for (int i=NDIM;i--;)
+                  x.x[i] = atof(argv[argptr+i]);
+                mol->Translate((const Vector *)&x);
+                argptr+=3;
+              }
+              break;
+            case 's':
+              ExitFlag = 1;
+              if ((argptr >= argc) || (!IsValidNumber(argv[argptr])) ) {
+                ExitFlag = 255;
+                cerr << "Not enough or invalid arguments given for scaling: -s <factor/[factor_x]> [factor_y] [factor_z]" << endl;
+              } else {
+                SaveFlag = true;
+                j = -1;
+                cout << Verbose(1) << "Scaling all ion positions by factor." << endl;
+                factor = new double[NDIM];
+                factor[0] = atof(argv[argptr]);
+                if ((argptr < argc) && (IsValidNumber(argv[argptr])))
+                  argptr++;
+                factor[1] = atof(argv[argptr]);
+                if ((argptr < argc) && (IsValidNumber(argv[argptr])))
+                  argptr++;
+                factor[2] = atof(argv[argptr]);
+                mol->Scale(&factor);
+                for (int i=0;i<NDIM;i++) {
+                  j += i+1;
+                  x.x[i] = atof(argv[NDIM+i]);
+                  mol->cell_size[j]*=factor[i];
+                }
+                delete[](factor);
+                argptr+=1;
+              }
+              break;
+            case 'b':
+              ExitFlag = 1;
+              if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
+                ExitFlag = 255;
+                cerr << "Not enough or invalid arguments given for centering in box: -b <length_x> <length_y> <length_z>" << endl;
+              } else {
+                SaveFlag = true;
+                j = -1;
+                cout << Verbose(1) << "Centering atoms in config file within given simulation box." << endl;
+                j=-1;
+                for (int i=0;i<NDIM;i++) {
+                  j += i+1;
+                  x.x[i] = atof(argv[argptr++]);
+                  mol->cell_size[j] += x.x[i]*2.;
+                }
+                // center
+                mol->CenterInBox((ofstream *)&cout, &x);
+                // update Box of atoms by boundary
+                mol->SetBoxDimension(&x);
+              }
+              break;
+            case 'c':
+              ExitFlag = 1;
+              if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
+                ExitFlag = 255;
+                cerr << "Not enough or invalid arguments given for centering with boundary: -c <boundary_x> <boundary_y> <boundary_z>" << endl;
+              } else {
+                SaveFlag = true;
+                j = -1;
+                cout << Verbose(1) << "Centering atoms in config file within given additional boundary." << endl;
+                // make every coordinate positive
+                mol->CenterEdge((ofstream *)&cout, &x);
+                // update Box of atoms by boundary
+                mol->SetBoxDimension(&x);
+                // translate each coordinate by boundary
+                j=-1;
+                for (int i=0;i<NDIM;i++) {
+                  j += i+1;
+                  x.x[i] = atof(argv[argptr++]);
+                  mol->cell_size[j] += x.x[i]*2.;
+                }
+                mol->Translate((const Vector *)&x);
+              }
+              break;
+            case 'O':
+              ExitFlag = 1;
+              SaveFlag = true;
+              cout << Verbose(1) << "Centering atoms in origin." << endl;
+              mol->CenterOrigin((ofstream *)&cout, &x);
+              mol->SetBoxDimension(&x);
+              break;
+            case 'r':
+              ExitFlag = 1;
+              SaveFlag = true;
+              cout << Verbose(1) << "Converting config file from supposed old to new syntax." << endl;
+              break;
+            case 'F':
+            case 'f':
+              ExitFlag = 1;
+              if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1]))) {
+                ExitFlag = 255;
+                cerr << "Not enough or invalid arguments for fragmentation: -f <max. bond distance> <bond order>" << endl;
+              } else {
+                cout << "Fragmenting molecule with bond distance " << argv[argptr] << " angstroem, order of " << argv[argptr+1] << "." << endl;
+                cout << Verbose(0) << "Creating connection matrix..." << endl;
+                start = clock();
+                mol->CreateAdjacencyList((ofstream *)&cout, atof(argv[argptr++]), configuration.GetIsAngstroem());
+                cout << Verbose(0) << "Fragmenting molecule with current connection matrix ..." << endl;
+                if (mol->first->next != mol->last) {
+                  ExitFlag = mol->FragmentMolecule((ofstream *)&cout, atoi(argv[argptr]), &configuration);
+                }
+                end = clock();
+                cout << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl;
+                argptr+=2;
+              }
+              break;
+            case 'm':
+              ExitFlag = 1;
+              j = atoi(argv[argptr++]);
+              if ((j<0) || (j>1)) {
+                cerr << Verbose(1) << "ERROR: Argument of '-m' should be either 0 for no-rotate or 1 for rotate." << endl;
+                j = 0;
+              }
+              if (j) {
+                SaveFlag = true;
+                cout << Verbose(0) << "Converting to prinicipal axis system." << endl;
+              } else
+                cout << Verbose(0) << "Evaluating prinicipal axis." << endl;
+              mol->PrincipalAxisSystem((ofstream *)&cout, (bool)j);
+              break;
+            case 'o':
+              ExitFlag = 1;
+              if ((argptr >= argc) || (argv[argptr][0] == '-')){
+                ExitFlag = 255;
+                cerr << "Not enough or invalid arguments given for convex envelope: -o <tecplot output file>" << endl;
+              } else {
+                cout << Verbose(0) << "Evaluating volume of the convex envelope.";
+                cout << Verbose(1) << "Storing tecplot data in " << argv[argptr] << "." << endl;
+                VolumeOfConvexEnvelope((ofstream *)&cout, argv[argptr], &configuration, NULL, mol);
+                argptr+=1;
+              }
+              break;
+            case 'U':
+              ExitFlag = 1;
+              if ((argptr+1 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) ) {
+                ExitFlag = 255;
+                cerr << "Not enough or invalid arguments given for suspension with specified volume: -U <volume> <density>" << endl;
+                volume = -1; // for case 'u': don't print error again
+              } else {
+                volume = atof(argv[argptr++]);
+                cout << Verbose(0) << "Using " << volume << " angstrom^3 as the volume instead of convex envelope one's." << endl;
+              }
+            case 'u':
+              ExitFlag = 1;
+              if ((argptr >= argc) || (!IsValidNumber(argv[argptr])) ) {
+                if (volume != -1)
+                  ExitFlag = 255;
+                  cerr << "Not enough arguments given for suspension: -u <density>" << endl;
+              } else {
+                double density;
+                SaveFlag = true;
+                cout << Verbose(0) << "Evaluating necessary cell volume for a cluster suspended in water.";
+                density = atof(argv[argptr++]);
+                if (density < 1.0) {
+                  cerr << Verbose(0) << "Density must be greater than 1.0g/cm^3 !" << endl;
+                  density = 1.3;
+                }
+//                for(int i=0;i<NDIM;i++) {
+//                  repetition[i] = atoi(argv[argptr++]);
+//                  if (repetition[i] < 1)
+//                    cerr << Verbose(0) << "ERROR: repetition value must be greater 1!" << endl;
+//                  repetition[i] = 1;
+//                }
+                PrepareClustersinWater((ofstream *)&cout, &configuration, mol, volume, density);  // if volume == 0, will calculate from ConvexEnvelope
+              }
+              break;
+            case 'd':
+              ExitFlag = 1;
+              if ((argptr+2 >= argc) || (!IsValidNumber(argv[argptr])) || (!IsValidNumber(argv[argptr+1])) || (!IsValidNumber(argv[argptr+2])) ) {
+                ExitFlag = 255;
+                cerr << "Not enough or invalid arguments given for repeating cells: -d <repeat_x> <repeat_y> <repeat_z>" << endl;
+              } else {
+                SaveFlag = true;
+                for (int axis = 1; axis <= NDIM; axis++) {
+                  int faktor = atoi(argv[argptr++]);
+                  int count;
+                  element ** Elements;
+                  Vector ** vectors;
+                  if (faktor < 1) {
+                    cerr << Verbose(0) << "ERROR: Repetition faktor mus be greater than 1!" << endl;
+                    faktor = 1;
+                  }
+                  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, ...)
+                      }
+                    }
+                    // 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;
+            default:   // no match? Step on
+              if ((argptr < argc) && (argv[argptr][0] != '-')) // if it started with a '-' we've already made a step!
+                argptr++;
+              break;
+          }
+        }
+      } else argptr++;
+    } while (argptr < argc);
+    if (SaveFlag)
+      SaveConfig(ConfigFileName, &configuration, periode, molecules);
+    if ((ExitFlag >= 1)) {
+      delete(mol);
+      delete(periode);
+      return (ExitFlag);
+    }
+  } else {  // no arguments, hence scan the elements db
+    if (periode->LoadPeriodentafel(PathToDatabases))
+      cout << Verbose(0) << "Element list loaded successfully." << endl;
+    else
+      cout << Verbose(0) << "Element list loading failed." << endl;
+    configuration.RetrieveConfigPathAndName("main_pcp_linux");
+  }
+  return(0);
 };
 
@@ -1827,28 +1827,28 @@
 int main(int argc, char **argv)
 {
-	periodentafel *periode = new periodentafel; // and a period table of all elements
-	MoleculeListClass *molecules = new MoleculeListClass;  // list of all molecules
-	molecule *mol = NULL;
-	config configuration;
-	double tmp1;
-	atom *first, *second;
-	char choice;	// menu choice char
-	Vector x,y,z,n;	// coordinates for absolute point in cell volume
-	bool valid; // flag if input was valid or not
-	ifstream test;
-	ofstream output;
-	string line;
-	char *ConfigFileName = NULL;
-	char *ElementsFileName = NULL;
-	int Z;
-	int j, axis, count, faktor;
-
-	// =========================== PARSE COMMAND LINE OPTIONS ====================================
-	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 (molecules->ListOfMolecules.size() == 0) {
+  periodentafel *periode = new periodentafel; // and a period table of all elements
+  MoleculeListClass *molecules = new MoleculeListClass;  // list of all molecules
+  molecule *mol = NULL;
+  config configuration;
+  double tmp1;
+  atom *first, *second;
+  char choice;  // menu choice char
+  Vector x,y,z,n;  // coordinates for absolute point in cell volume
+  bool valid; // flag if input was valid or not
+  ifstream test;
+  ofstream output;
+  string line;
+  char *ConfigFileName = NULL;
+  char *ElementsFileName = NULL;
+  int Z;
+  int j, axis, count, faktor;
+
+  // =========================== PARSE COMMAND LINE OPTIONS ====================================
+  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 (molecules->ListOfMolecules.size() == 0) {
     mol = new molecule(periode);
     if (mol->cell_size[0] == 0.) {
@@ -1860,15 +1860,15 @@
     }
     molecules->insert(mol);
-	}
-
-	// =========================== START INTERACTIVE SESSION ====================================
-
-	// now the main construction loop
-	cout << Verbose(0) << endl << "Now comes the real construction..." << endl;
-	do {
-		cout << Verbose(0) << endl << endl;
-		cout << Verbose(0) << "============Molecule list=======================" << endl;
-		molecules->Enumerate((ofstream *)&cout);
-		cout << Verbose(0) << "============Menu===============================" << endl;
+  }
+
+  // =========================== START INTERACTIVE SESSION ====================================
+
+  // now the main construction loop
+  cout << Verbose(0) << endl << "Now comes the real construction..." << endl;
+  do {
+    cout << Verbose(0) << endl << endl;
+    cout << Verbose(0) << "============Molecule list=======================" << endl;
+    molecules->Enumerate((ofstream *)&cout);
+    cout << Verbose(0) << "============Menu===============================" << endl;
     cout << Verbose(0) << "a - set molecule (in)active" << endl;
     cout << Verbose(0) << "e - edit new molecules" << endl;
@@ -1886,6 +1886,6 @@
     cin >> choice;
 
-		switch (choice) {
-			case 'a':  // (in)activate molecule
+    switch (choice) {
+      case 'a':  // (in)activate molecule
         {
           cout << "Enter index of molecule: ";
@@ -1897,9 +1897,9 @@
             (*ListRunner)->ActiveFlag = !(*ListRunner)->ActiveFlag;
         }
-			  break;
-
-			case 'c': // edit each field of the configuration
-			 configuration.Edit();
-			 break;
+        break;
+
+      case 'c': // edit each field of the configuration
+       configuration.Edit();
+       break;
 
       case 'e': // create molecule
@@ -1919,29 +1919,29 @@
         break;
 
-			case 'q': // quit
-				break;
-
-			case 's': // save to config file
-				SaveConfig(ConfigFileName, &configuration, periode, molecules);
-				break;
-
-			case 'T':
-				testroutine(molecules);
-				break;
-
-			default:
-			  break;
-		};
-	} while (choice != 'q');
-
-	// save element data base
-	if (periode->StorePeriodentafel(ElementsFileName)) //ElementsFileName
-		cout << Verbose(0) << "Saving of elements.db successful." << endl;
-	else
-		cout << Verbose(0) << "Saving of elements.db failed." << endl;
-
-	delete(molecules);
-	delete(periode);
-	return (0);
+      case 'q': // quit
+        break;
+
+      case 's': // save to config file
+        SaveConfig(ConfigFileName, &configuration, periode, molecules);
+        break;
+
+      case 'T':
+        testroutine(molecules);
+        break;
+
+      default:
+        break;
+    };
+  } while (choice != 'q');
+
+  // save element data base
+  if (periode->StorePeriodentafel(ElementsFileName)) //ElementsFileName
+    cout << Verbose(0) << "Saving of elements.db successful." << endl;
+  else
+    cout << Verbose(0) << "Saving of elements.db failed." << endl;
+
+  delete(molecules);
+  delete(periode);
+  return (0);
 }
 
