Index: molecuilder/src/molecules.cpp
===================================================================
--- molecuilder/src/molecules.cpp	(revision 5609952e2818195ff237d87e773818aacb6db3a5)
+++ molecuilder/src/molecules.cpp	(revision f3973579348a95ed314bf56d258365dbdc97f738)
@@ -62,4 +62,7 @@
   cell_size[0] = cell_size[2] = cell_size[5]= 20.;
   cell_size[1] = cell_size[3] = cell_size[4]= 0.;
+  strcpy(name,"none");
+  IndexNr  = -1;
+  ActiveFlag = false;
 };
 
@@ -598,4 +601,20 @@
 };
 
+/** Set molecule::name from the basename without suffix in the given \a *filename.
+ * \param *filename filename
+ */
+void molecule::SetNameFromFilename(const char *filename)
+{
+  int length = 0;
+  char *molname = strrchr(filename, '/')+sizeof(char);  // search for filename without dirs
+  char *endname = strchr(molname, '.');
+  if ((endname == NULL) || (endname < molname))
+    length = strlen(molname);
+  else
+    length = strlen(molname) - strlen(endname);
+  strncpy(name, molname, length);
+  name[length]='\0';
+};
+
 /** Sets the molecule::cell_size to the components of \a *dim (rectangular box)
  * \param *dim vector class
@@ -652,4 +671,5 @@
     max->Scale(0.5);
     Translate(max);
+    Center.Zero();
   }
 
@@ -691,4 +711,5 @@
     max->AddVector(min);
     Translate(min);
+    Center.Zero();
   }
   delete(min);
@@ -700,11 +721,10 @@
  * \param *center return vector for translation vector
  */
-void molecule::CenterOrigin(ofstream *out, Vector *center)
+void molecule::CenterOrigin(ofstream *out)
 {
   int Num = 0;
   atom *ptr = start->next;  // start at first in list
 
-  for(int i=NDIM;i--;) // zero center vector
-    center->x[i] = 0.;
+  Center.Zero();
 
   if (ptr != end) {   //list not empty?
@@ -712,8 +732,9 @@
       ptr = ptr->next;
       Num++;
-      center->AddVector(&ptr->x);
-    }
-    center->Scale(-1./Num); // divide through total number (and sign for direction)
-    Translate(center);
+      Center.AddVector(&ptr->x);
+    }
+    Center.Scale(-1./Num); // divide through total number (and sign for direction)
+    Translate(&Center);
+    Center.Zero();
   }
 };
@@ -780,13 +801,16 @@
  * \param *center return vector for translation vector
  */
-void molecule::CenterGravity(ofstream *out, Vector *center)
-{
-  if (center == NULL) {
-    DetermineCenter(*center);
-    Translate(center);
-    delete(center);
-  } else {
-    Translate(center);
-  }
+void molecule::CenterPeriodic(ofstream *out)
+{
+  DeterminePeriodicCenter(Center);
+};
+
+/** Centers the center of gravity of the atoms at (0,0,0).
+ * \param *out output stream for debugging
+ * \param *center return vector for translation vector
+ */
+void molecule::CenterAtVector(ofstream *out, Vector *newcenter)
+{
+  Center.CopyVector(newcenter);
 };
 
@@ -837,7 +861,7 @@
 
 /** Determines center of molecule (yet not considering atom masses).
- * \param Center reference to return vector
- */
-void molecule::DetermineCenter(Vector &Center)
+ * \param center reference to return vector
+ */
+void molecule::DeterminePeriodicCenter(Vector &center)
 {
   atom *Walker = start;
@@ -849,5 +873,5 @@
 
   do {
-    Center.Zero();
+    center.Zero();
     flag = true;
     while (Walker->next != end) {
@@ -876,5 +900,5 @@
         Testvector.AddVector(&Translationvector);
         Testvector.MatrixMultiplication(matrix);
-        Center.AddVector(&Testvector);
+        center.AddVector(&Testvector);
         cout << Verbose(1) << "vector is: ";
         Testvector.Output((ofstream *)&cout);
@@ -889,5 +913,5 @@
             Testvector.AddVector(&Translationvector);
             Testvector.MatrixMultiplication(matrix);
-            Center.AddVector(&Testvector);
+            center.AddVector(&Testvector);
             cout << Verbose(1) << "Hydrogen vector is: ";
             Testvector.Output((ofstream *)&cout);
@@ -900,5 +924,5 @@
   } while (!flag);
   Free((void **)&matrix, "molecule::DetermineCenter: *matrix");
-  Center.Scale(1./(double)AtomCount);
+  center.Scale(1./(double)AtomCount);
 };
 
@@ -913,5 +937,5 @@
   Vector *CenterOfGravity = DetermineCenterOfGravity(out);
 
-  CenterGravity(out, CenterOfGravity);
+  CenterAtVector(out, CenterOfGravity);
 
   // reset inertia tensor
@@ -1409,7 +1433,8 @@
 bool molecule::LinearInterpolationBetweenConfiguration(ofstream *out, int startstep, int endstep, const char *prefix, config &configuration) 
 {
+  molecule *mol = NULL;
   bool status = true;
   int MaxSteps = configuration.MaxOuterStep; 
-  MoleculeListClass *MoleculePerStep = new MoleculeListClass(MaxSteps+1, AtomCount);
+  MoleculeListClass *MoleculePerStep = new MoleculeListClass();
   // Get the Permutation Map by MinimiseConstrainedPotential
   atom **PermutationMap = NULL;
@@ -1443,10 +1468,11 @@
   *out << Verbose(1) << "Filling intermediate " << MaxSteps << " steps with MDSteps of " << MDSteps << "." << endl;
   for (int step = 0; step <= MaxSteps; step++) {
-    MoleculePerStep->ListOfMolecules[step] = new molecule(elemente);
+    mol = new molecule(elemente);
+    MoleculePerStep->insert(mol);
     Walker = start;
     while (Walker->next != end) {
       Walker = Walker->next;
       // add to molecule list
-      Sprinter = MoleculePerStep->ListOfMolecules[step]->AddCopyAtom(Walker);
+      Sprinter = mol->AddCopyAtom(Walker);
       for (int n=NDIM;n--;) {
         Sprinter->x.x[n] = Trajectories[Walker].R.at(startstep).x[n] + (Trajectories[PermutationMap[Walker->nr]].R.at(endstep).x[n] - Trajectories[Walker].R.at(startstep).x[n])*((double)step/(double)MaxSteps);
@@ -1467,5 +1493,5 @@
   for (int i=AtomCount; i--; )
     SortIndex[i] = i;
-  status = MoleculePerStep->OutputConfigForListOfFragments(out, "ConstrainedStep", &configuration, SortIndex, false, false);
+  status = MoleculePerStep->OutputConfigForListOfFragments(out, &configuration, SortIndex);
   
   // free and return
@@ -1828,5 +1854,5 @@
 };
 
-/** Removes atom from molecule list.
+/** Removes atom from molecule list and deletes it.
  * \param *pointer atom to be removed
  * \return true - succeeded, false - atom not found in list
@@ -1842,4 +1868,23 @@
   Trajectories.erase(pointer);
   return remove(pointer, start, end);
+};
+
+/** Removes atom from molecule list, but does not delete it.
+ * \param *pointer atom to be removed
+ * \return true - succeeded, false - atom not found in list
+ */
+bool molecule::UnlinkAtom(atom *pointer)
+{
+  if (pointer == NULL)
+    return false;
+  if (ElementsInMolecule[pointer->type->Z] != 0)  // this would indicate an error
+    ElementsInMolecule[pointer->type->Z]--; // decrease number of atom of this element
+  else
+    cerr << "ERROR: Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl;
+  if (ElementsInMolecule[pointer->type->Z] == 0)  // was last atom of this element?
+    ElementCount--;
+  Trajectories.erase(pointer);
+  unlink(pointer);
+  return true;
 };
 
@@ -2368,5 +2413,5 @@
             };
             AddBond(Walker, OtherWalker); //Add the bond between the two atoms with respective indices.
-            
+
           }
 
@@ -2484,14 +2529,13 @@
                         //*out << Verbose(0) << "Current comparison atom is " << *OtherWalker << "." << endl;
                         /// \todo periodic check is missing here!
-                        //*out << Verbose(1) << "Checking distance " << OtherWalker->x.PeriodicDistance(&(Walker->x), cell_size) << " against typical bond length of " << bonddistance*bonddistance << "." << endl;
+                        //*out << Verbose(1) << "Checking distance " << OtherWalker->x.PeriodicDistanceSquared(&(Walker->x), cell_size) << " against typical bond length of " << bonddistance*bonddistance << "." << endl;
                         MinDistance = OtherWalker->type->CovalentRadius + Walker->type->CovalentRadius;
                         MinDistance *= (IsAngstroem) ? 1. : 1./AtomicLengthToAngstroem;
                         MaxDistance = MinDistance + BONDTHRESHOLD;
                         MinDistance -= BONDTHRESHOLD;
-                        distance =   OtherWalker->x.PeriodicDistance(&(Walker->x), cell_size);
+                        distance = OtherWalker->x.PeriodicDistanceSquared(&(Walker->x), cell_size);
                         if ((OtherWalker->father->nr > Walker->father->nr) && (distance <= MaxDistance*MaxDistance) && (distance >= MinDistance*MinDistance)) { // create bond if distance is smaller
-                          //*out << Verbose(0) << "Adding Bond between " << *Walker << " and " << *OtherWalker << "." << endl;
+                          //*out << Verbose(1) << "Adding Bond between " << *Walker << " and " << *OtherWalker << " in distance " << sqrt(distance) << "." << endl;
                           AddBond(Walker->father, OtherWalker->father, 1);  // also increases molecule::BondCount
-                          BondCount++;
                         } else {
                           //*out << Verbose(1) << "Not Adding: Wrong label order or distance too great." << endl;
@@ -2564,5 +2608,5 @@
       *out << Verbose(1) << "BondCount is " << BondCount << ", no bonds between any of the " << AtomCount << " atoms." << endl;
     *out << Verbose(1) << "I detected " << BondCount << " bonds in the molecule with distance " << bonddistance << ", " << FalseBondDegree << " bonds could not be corrected." << endl;
-    
+
     // output bonds for debugging (if bond chain list was correctly installed)
     *out << Verbose(1) << endl << "From contents of bond chain list:";
@@ -3171,5 +3215,4 @@
           cerr << "KeySet file must be corrupt as there are two equal key sets therein!" << endl;
         }
-        //FragmentList->ListOfMolecules[NumberOfFragments++] = StoreFragmentFromKeySet(out, CurrentSet, IsAngstroem);
       }
     }
@@ -3731,21 +3774,21 @@
   //if (FragmentationToDo) {    // we should always store the fragments again as coordination might have changed slightly without changing bond structure
     // allocate memory for the pointer array and transmorph graphs into full molecular fragments
-    BondFragments = new MoleculeListClass(TotalGraph.size(), AtomCount);
+    BondFragments = new MoleculeListClass();
     int k=0;
     for(Graph::iterator runner = TotalGraph.begin(); runner != TotalGraph.end(); runner++) {
       KeySet test = (*runner).first;
       *out << "Fragment No." << (*runner).second.first << " with TEFactor " << (*runner).second.second << "." << endl;
-      BondFragments->ListOfMolecules[k] = StoreFragmentFromKeySet(out, test, configuration);
+      BondFragments->insert(StoreFragmentFromKeySet(out, test, configuration));
       k++;
     }
-    *out << k << "/" << BondFragments->NumberOfMolecules << " fragments generated from the keysets." << endl;
+    *out << k << "/" << BondFragments->ListOfMolecules.size() << " fragments generated from the keysets." << endl;
 
     // ===== 9. Save fragments' configuration and keyset files et al to disk ===
-    if (BondFragments->NumberOfMolecules != 0) {
+    if (BondFragments->ListOfMolecules.size() != 0) {
       // create the SortIndex from BFS labels to order in the config file
       CreateMappingLabelsToConfigSequence(out, SortIndex);
 
-      *out << Verbose(1) << "Writing " << BondFragments->NumberOfMolecules << " possible bond fragmentation configs" << endl;
-      if (BondFragments->OutputConfigForListOfFragments(out, FRAGMENTPREFIX, configuration, SortIndex, true, true))
+      *out << Verbose(1) << "Writing " << BondFragments->ListOfMolecules.size() << " possible bond fragmentation configs" << endl;
+      if (BondFragments->OutputConfigForListOfFragments(out, configuration, SortIndex))
         *out << Verbose(1) << "All configs written." << endl;
       else
@@ -5318,6 +5361,6 @@
   if (result) {
     *out << Verbose(5) << "Calculating Centers of Gravity" << endl;
-    DetermineCenter(CenterOfGravity);
-    OtherMolecule->DetermineCenter(OtherCenterOfGravity);
+    DeterminePeriodicCenter(CenterOfGravity);
+    OtherMolecule->DeterminePeriodicCenter(OtherCenterOfGravity);
     *out << Verbose(5) << "Center of Gravity: ";
     CenterOfGravity.Output(out);
@@ -5325,5 +5368,5 @@
     OtherCenterOfGravity.Output(out);
     *out << endl;
-    if (CenterOfGravity.Distance(&OtherCenterOfGravity) > threshold) {
+    if (CenterOfGravity.DistanceSquared(&OtherCenterOfGravity) > threshold*threshold) {
       *out << Verbose(4) << "Centers of gravity don't match." << endl;
       result = false;
@@ -5339,10 +5382,10 @@
     while (Walker->next != end) {
       Walker = Walker->next;
-      Distances[Walker->nr] = CenterOfGravity.Distance(&Walker->x);
+      Distances[Walker->nr] = CenterOfGravity.DistanceSquared(&Walker->x);
     }
     Walker = OtherMolecule->start;
     while (Walker->next != OtherMolecule->end) {
       Walker = Walker->next;
-      OtherDistances[Walker->nr] = OtherCenterOfGravity.Distance(&Walker->x);
+      OtherDistances[Walker->nr] = OtherCenterOfGravity.DistanceSquared(&Walker->x);
     }
 
@@ -5362,6 +5405,6 @@
     flag = 0;
     for (int i=0;i<AtomCount;i++) {
-      *out << Verbose(5) << "Distances: |" << Distances[PermMap[i]] << " - " << OtherDistances[OtherPermMap[i]] << "| = " << fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) << " ?<? " <<  threshold << endl;
-      if (fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) > threshold)
+      *out << Verbose(5) << "Distances squared: |" << Distances[PermMap[i]] << " - " << OtherDistances[OtherPermMap[i]] << "| = " << fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) << " ?<? " <<  threshold << endl;
+      if (fabs(Distances[PermMap[i]] - OtherDistances[OtherPermMap[i]]) > threshold*threshold)
         flag = 1;
     }
