Index: src/molecule.cpp
===================================================================
--- src/molecule.cpp	(revision 8d27720451e7340d82da7537d832ffb204273b10)
+++ src/molecule.cpp	(revision ea71760f1257d93e6396b87b269d71544d16a690)
@@ -32,10 +32,9 @@
  */
 molecule::molecule(const periodentafel * const teil) : elemente(teil),
-  first(new bond(0, 0, 1, -1)), last(new bond(0, 0, 1, -1)), MDSteps(0), AtomCount(0),
+  first(new bond(0, 0, 1, -1)), last(new bond(0, 0, 1, -1)), MDSteps(0),
   BondCount(0), ElementCount(0), NoNonHydrogen(0), NoNonBonds(0), NoCyclicBonds(0), BondDistance(0.), 
   ActiveFlag(false), IndexNr(-1),
   formula(this,boost::bind(&molecule::calcFormula,this)),
-  last_atom(0),
-  InternalPointer(begin())
+  AtomCount(this,boost::bind(&molecule::doCountAtoms,this)), last_atom(0),  InternalPointer(begin())
 {
   // init bond chain list
@@ -72,4 +71,8 @@
 const std::string molecule::getName(){
   return std::string(name);
+}
+
+int molecule::getAtomCount() const{
+  return *AtomCount;
 }
 
@@ -180,5 +183,4 @@
   if (pointer != NULL) {
     pointer->sort = &pointer->nr;
-    AtomCount++;
     if (pointer->type != NULL) {
       if (ElementsInMolecule[pointer->type->Z] == 0)
@@ -214,5 +216,4 @@
     if ((pointer->type != NULL) && (pointer->type->Z != 1))
       NoNonHydrogen++;
-    AtomCount++;
     retval=walker;
   }
@@ -611,5 +612,4 @@
 
   // copy values
-  copy->CountAtoms();
   copy->CountElements();
   if (first->next != last) {  // if adjaceny list is present
@@ -728,7 +728,7 @@
 bool molecule::RemoveAtom(atom *pointer)
 {
+  OBSERVE;
   if (ElementsInMolecule[pointer->type->Z] != 0)  { // this would indicate an error
     ElementsInMolecule[pointer->type->Z]--;  // decrease number of atom of this element
-    AtomCount--;
   } else
     eLog() << Verbose(1) << "Atom " << pointer->Name << " is of element " << pointer->type->Z << " but the entry in the table of the molecule is 0!" << endl;
@@ -909,5 +909,5 @@
     now = time((time_t *)NULL);   // Get the system time and put it into 'now' as 'calender time'
     for (int step=0;step<MDSteps;step++) {
-      *output << AtomCount << "\n\tCreated by molecuilder, step " << step << ", on " << ctime(&now);
+      *output << getAtomCount() << "\n\tCreated by molecuilder, step " << step << ", on " << ctime(&now);
       ActOnAllAtoms( &atom::OutputTrajectoryXYZ, output, step );
     }
@@ -926,5 +926,5 @@
   if (output != NULL) {
     now = time((time_t *)NULL);   // Get the system time and put it into 'now' as 'calender time'
-    *output << AtomCount << "\n\tCreated by molecuilder on " << ctime(&now);
+    *output << getAtomCount() << "\n\tCreated by molecuilder on " << ctime(&now);
     ActOnAllAtoms( &atom::OutputXYZLine, output );
     return true;
@@ -936,29 +936,20 @@
  * \param *out output stream for debugging
  */
-void molecule::CountAtoms()
-{
-  int i = size();
-  if ((AtomCount == 0) || (i != AtomCount)) {
-    cout << "!!!!!!!!! Counting needed" << endl;
-    Log() << Verbose(3) << "Mismatch in AtomCount " << AtomCount << " and recounted number " << i << ", renaming all." << endl;
-    AtomCount = i;
-
-    // count NonHydrogen atoms and give each atom a unique name
-    if (AtomCount != 0) {
-      i=0;
-      NoNonHydrogen = 0;
-      for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
-        (*iter)->nr = i;   // update number in molecule (for easier referencing in FragmentMolecule lateron)
-        if ((*iter)->type->Z != 1) // count non-hydrogen atoms whilst at it
-          NoNonHydrogen++;
-        Free(&(*iter)->Name);
-        (*iter)->Name = Malloc<char>(6, "molecule::CountAtoms: *walker->Name");
-        sprintf((*iter)->Name, "%2s%02d", (*iter)->type->symbol, (*iter)->nr+1);
-        Log() << Verbose(3) << "Naming atom nr. " << (*iter)->nr << " " << (*iter)->Name << "." << endl;
-        i++;
-      }
-    } else
-      Log() << Verbose(3) << "AtomCount is still " << AtomCount << ", thus counting nothing." << endl;
-  }
+int molecule::doCountAtoms()
+{
+  int res = size();
+  int i = 0;
+  NoNonHydrogen = 0;
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    (*iter)->nr = i;   // update number in molecule (for easier referencing in FragmentMolecule lateron)
+    if ((*iter)->type->Z != 1) // count non-hydrogen atoms whilst at it
+      NoNonHydrogen++;
+    Free(&(*iter)->Name);
+    (*iter)->Name = Malloc<char>(6, "molecule::CountAtoms: *walker->Name");
+    sprintf((*iter)->Name, "%2s%02d", (*iter)->type->symbol, (*iter)->nr+1);
+    Log() << Verbose(3) << "Naming atom nr. " << (*iter)->nr << " " << (*iter)->Name << "." << endl;
+    i++;
+  }
+  return res;
 };
 
@@ -1022,6 +1013,4 @@
   /// first count both their atoms and elements and update lists thereby ...
   //Log() << Verbose(0) << "Counting atoms, updating list" << endl;
-  CountAtoms();
-  OtherMolecule->CountAtoms();
   CountElements();
   OtherMolecule->CountElements();
@@ -1030,8 +1019,8 @@
   /// -# AtomCount
   if (result) {
-    if (AtomCount != OtherMolecule->AtomCount) {
-      Log() << Verbose(4) << "AtomCounts don't match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl;
+    if (getAtomCount() != OtherMolecule->getAtomCount()) {
+      Log() << Verbose(4) << "AtomCounts don't match: " << getAtomCount() << " == " << OtherMolecule->getAtomCount() << endl;
       result = false;
-    } else Log() << Verbose(4) << "AtomCounts match: " << AtomCount << " == " << OtherMolecule->AtomCount << endl;
+    } else Log() << Verbose(4) << "AtomCounts match: " << getAtomCount() << " == " << OtherMolecule->getAtomCount() << endl;
   }
   /// -# ElementCount
@@ -1073,6 +1062,6 @@
   if (result) {
     Log() << Verbose(5) << "Calculating distances" << endl;
-    Distances = Calloc<double>(AtomCount, "molecule::IsEqualToWithinThreshold: Distances");
-    OtherDistances = Calloc<double>(AtomCount, "molecule::IsEqualToWithinThreshold: OtherDistances");
+    Distances = Calloc<double>(getAtomCount(), "molecule::IsEqualToWithinThreshold: Distances");
+    OtherDistances = Calloc<double>(getAtomCount(), "molecule::IsEqualToWithinThreshold: OtherDistances");
     SetIndexedArrayForEachAtomTo ( Distances, &atom::nr, &atom::DistanceSquaredToVector, (const Vector &)CenterOfGravity);
     SetIndexedArrayForEachAtomTo ( OtherDistances, &atom::nr, &atom::DistanceSquaredToVector, (const Vector &)CenterOfGravity);
@@ -1080,11 +1069,11 @@
     /// ... sort each list (using heapsort (o(N log N)) from GSL)
     Log() << Verbose(5) << "Sorting distances" << endl;
-    PermMap = Calloc<size_t>(AtomCount, "molecule::IsEqualToWithinThreshold: *PermMap");
-    OtherPermMap = Calloc<size_t>(AtomCount, "molecule::IsEqualToWithinThreshold: *OtherPermMap");
-    gsl_heapsort_index (PermMap, Distances, AtomCount, sizeof(double), CompareDoubles);
-    gsl_heapsort_index (OtherPermMap, OtherDistances, AtomCount, sizeof(double), CompareDoubles);
-    PermutationMap = Calloc<int>(AtomCount, "molecule::IsEqualToWithinThreshold: *PermutationMap");
+    PermMap = Calloc<size_t>(getAtomCount(), "molecule::IsEqualToWithinThreshold: *PermMap");
+    OtherPermMap = Calloc<size_t>(getAtomCount(), "molecule::IsEqualToWithinThreshold: *OtherPermMap");
+    gsl_heapsort_index (PermMap, Distances, getAtomCount(), sizeof(double), CompareDoubles);
+    gsl_heapsort_index (OtherPermMap, OtherDistances, getAtomCount(), sizeof(double), CompareDoubles);
+    PermutationMap = Calloc<int>(getAtomCount(), "molecule::IsEqualToWithinThreshold: *PermutationMap");
     Log() << Verbose(5) << "Combining Permutation Maps" << endl;
-    for(int i=AtomCount;i--;)
+    for(int i=getAtomCount();i--;)
       PermutationMap[PermMap[i]] = (int) OtherPermMap[i];
 
@@ -1092,5 +1081,5 @@
     Log() << Verbose(4) << "Comparing distances" << endl;
     flag = 0;
-    for (int i=0;i<AtomCount;i++) {
+    for (int i=0;i<getAtomCount();i++) {
       Log() << 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)
@@ -1129,9 +1118,9 @@
 {
   Log() << Verbose(3) << "Begin of GetFatherAtomicMap." << endl;
-  int *AtomicMap = Malloc<int>(AtomCount, "molecule::GetAtomicMap: *AtomicMap");
-  for (int i=AtomCount;i--;)
+  int *AtomicMap = Malloc<int>(getAtomCount(), "molecule::GetAtomicMap: *AtomicMap");
+  for (int i=getAtomCount();i--;)
     AtomicMap[i] = -1;
   if (OtherMolecule == this) {  // same molecule
-    for (int i=AtomCount;i--;) // no need as -1 means already that there is trivial correspondence
+    for (int i=getAtomCount();i--;) // no need as -1 means already that there is trivial correspondence
       AtomicMap[i] = i;
     Log() << Verbose(4) << "Map is trivial." << endl;
