Index: src/molecule.cpp
===================================================================
--- src/molecule.cpp	(revision 2fe9719f8a0aa9ecd805749ffe7ff7b10ebcc6fe)
+++ src/molecule.cpp	(revision 822f019130ed6abc1bd095b4f6652208527de7a9)
@@ -983,110 +983,4 @@
 };
 
-/** Determines whether two molecules actually contain the same atoms and coordination.
- * \param *out output stream for debugging
- * \param *OtherMolecule the molecule to compare this one to
- * \param threshold upper limit of difference when comparing the coordination.
- * \return NULL - not equal, otherwise an allocated (molecule::AtomCount) permutation map of the atom numbers (which corresponds to which)
- */
-int * molecule::IsEqualToWithinThreshold(molecule *OtherMolecule, double threshold)
-{
-  int flag;
-  double *Distances = NULL, *OtherDistances = NULL;
-  Vector CenterOfGravity, OtherCenterOfGravity;
-  size_t *PermMap = NULL, *OtherPermMap = NULL;
-  int *PermutationMap = NULL;
-  bool result = true; // status of comparison
-
-  DoLog(3) && (Log() << Verbose(3) << "Begin of IsEqualToWithinThreshold." << endl);
-  /// first count both their atoms and elements and update lists thereby ...
-  //Log() << Verbose(0) << "Counting atoms, updating list" << endl;
-
-  /// ... and compare:
-  /// -# AtomCount
-  if (result) {
-    if (getAtomCount() != OtherMolecule->getAtomCount()) {
-      DoLog(4) && (Log() << Verbose(4) << "AtomCounts don't match: " << getAtomCount() << " == " << OtherMolecule->getAtomCount() << endl);
-      result = false;
-    } else Log() << Verbose(4) << "AtomCounts match: " << getAtomCount() << " == " << OtherMolecule->getAtomCount() << endl;
-  }
-  /// -# Formula
-  if (result) {
-    if (formula != OtherMolecule->formula) {
-      DoLog(4) && (Log() << Verbose(4) << "Formulas don't match: " << formula << " == " << OtherMolecule->formula << endl);
-      result = false;
-    } else Log() << Verbose(4) << "Formulas match: " << formula << " == " << OtherMolecule->formula << endl;
-  }
-  /// then determine and compare center of gravity for each molecule ...
-  if (result) {
-    DoLog(5) && (Log() << Verbose(5) << "Calculating Centers of Gravity" << endl);
-    DeterminePeriodicCenter(CenterOfGravity);
-    OtherMolecule->DeterminePeriodicCenter(OtherCenterOfGravity);
-    DoLog(5) && (Log() << Verbose(5) << "Center of Gravity: " << CenterOfGravity << endl);
-    DoLog(5) && (Log() << Verbose(5) << "Other Center of Gravity: " << OtherCenterOfGravity << endl);
-    if (CenterOfGravity.DistanceSquared(OtherCenterOfGravity) > threshold*threshold) {
-      DoLog(4) && (Log() << Verbose(4) << "Centers of gravity don't match." << endl);
-      result = false;
-    }
-  }
-
-  /// ... then make a list with the euclidian distance to this center for each atom of both molecules
-  if (result) {
-    DoLog(5) && (Log() << Verbose(5) << "Calculating distances" << endl);
-    Distances = new double[getAtomCount()];
-    OtherDistances = new double[getAtomCount()];
-    SetIndexedArrayForEachAtomTo ( Distances, &atom::nr, &atom::DistanceSquaredToVector, (const Vector &)CenterOfGravity);
-    SetIndexedArrayForEachAtomTo ( OtherDistances, &atom::nr, &atom::DistanceSquaredToVector, (const Vector &)CenterOfGravity);
-    for(int i=0;i<getAtomCount();i++) {
-      Distances[i] = 0.;
-      OtherDistances[i] = 0.;
-    }
-
-    /// ... sort each list (using heapsort (o(N log N)) from GSL)
-    DoLog(5) && (Log() << Verbose(5) << "Sorting distances" << endl);
-    PermMap = new size_t[getAtomCount()];
-    OtherPermMap = new size_t[getAtomCount()];
-    for(int i=0;i<getAtomCount();i++) {
-      PermMap[i] = 0;
-      OtherPermMap[i] = 0;
-    }
-    gsl_heapsort_index (PermMap, Distances, getAtomCount(), sizeof(double), CompareDoubles);
-    gsl_heapsort_index (OtherPermMap, OtherDistances, getAtomCount(), sizeof(double), CompareDoubles);
-    PermutationMap = new int[getAtomCount()];
-    for(int i=0;i<getAtomCount();i++)
-      PermutationMap[i] = 0;
-    DoLog(5) && (Log() << Verbose(5) << "Combining Permutation Maps" << endl);
-    for(int i=getAtomCount();i--;)
-      PermutationMap[PermMap[i]] = (int) OtherPermMap[i];
-
-    /// ... and compare them step by step, whether the difference is individually(!) below \a threshold for all
-    DoLog(4) && (Log() << Verbose(4) << "Comparing distances" << endl);
-    flag = 0;
-    for (int i=0;i<getAtomCount();i++) {
-      DoLog(5) && (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)
-        flag = 1;
-    }
-
-    // free memory
-    delete[](PermMap);
-    delete[](OtherPermMap);
-    delete[](Distances);
-    delete[](OtherDistances);
-    if (flag) { // if not equal
-      delete[](PermutationMap);
-      result = false;
-    }
-  }
-  /// return pointer to map if all distances were below \a threshold
-  DoLog(3) && (Log() << Verbose(3) << "End of IsEqualToWithinThreshold." << endl);
-  if (result) {
-    DoLog(3) && (Log() << Verbose(3) << "Result: Equal." << endl);
-    return PermutationMap;
-  } else {
-    DoLog(3) && (Log() << Verbose(3) << "Result: Not equal." << endl);
-    return NULL;
-  }
-};
-
 /** Returns an index map for two father-son-molecules.
  * The map tells which atom in this molecule corresponds to which one in the other molecul with their fathers.
