Changes in src/molecule.cpp [952f38:83f176]

File:

• src/molecule.cpp (modified) (37 diffs)

src/molecule.cpp

@@ +1 @@
+/*
+ * Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2010 University of Bonn. All rights reserved.
+ * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
+ */
+
 /** \file molecules.cpp
  *
@@ -5 +12 @@
  */
 
+// include config.h
 #ifdef HAVE_CONFIG_H
 #include <config.h>
@@ -89 +97 @@
 }
 
+bool molecule::changeId(moleculeId_t newId){
+  // first we move ourselves in the world
+  // the world lets us know if that succeeded
+  if(World::getInstance().changeMoleculeId(id,newId,this)){
+    id = newId;
+    return true;
+  }
+  else{
+    return false;
+  }
+}
+
+
 moleculeId_t molecule::getId(){
   return id;
@@ -151 +172 @@
 molecule::const_iterator molecule::erase( const_iterator loc )
 {
+  OBSERVE;
   molecule::const_iterator iter = loc;
   iter--;
@@ -156 +178 @@
   atomIds.erase( atom->getId() );
   atoms.remove( atom );
+  formula-=atom->getType();
   atom->removeFromMolecule();
   return iter;
@@ -162 +185 @@
 molecule::const_iterator molecule::erase( atom * key )
 {
+  OBSERVE;
   molecule::const_iterator iter = find(key);
   if (iter != end()){
     atomIds.erase( key->getId() );
     atoms.remove( key );
+    formula-=key->getType();
     key->removeFromMolecule();
   }
@@ -183 +208 @@
 pair<molecule::iterator,bool> molecule::insert ( atom * const key )
 {
+  OBSERVE;
   pair<atomIdSet::iterator,bool> res = atomIds.insert(key->getId());
   if (res.second) { // push atom if went well
     atoms.push_back(key);
+    formula+=key->getType();
     return pair<iterator,bool>(molecule::iterator(--end()),res.second);
   } else {
@@ -206 +233 @@
   if (pointer != NULL) {
     pointer->sort = &pointer->nr;
-    if (pointer->type != NULL) {
-      formula += pointer->type;
-      if (pointer->type->Z != 1)
+    if (pointer->getType() != NULL) {
+      if (pointer->getType()->getAtomicNumber() != 1)
         NoNonHydrogen++;
       if(pointer->getName() == "Unknown"){
         stringstream sstr;
-        sstr << pointer->type->symbol << pointer->nr+1;
+        sstr << pointer->getType()->getSymbol() << pointer->nr+1;
         pointer->setName(sstr.str());
       }
@@ -236 +262 @@
     walker->nr = last_atom++;  // increase number within molecule
     insert(walker);
-    if ((pointer->type != NULL) && (pointer->type->Z != 1))
+    if ((pointer->getType() != NULL) && (pointer->getType()->getAtomicNumber() != 1))
       NoNonHydrogen++;
     retval=walker;
@@ -294 +320 @@
 //  Log() << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl;
   // create vector in direction of bond
-  InBondvector = TopReplacement->x - TopOrigin->x;
+  InBondvector = TopReplacement->getPosition() - TopOrigin->getPosition();
   bondlength = InBondvector.Norm();
 
@@ -306 +332 @@
     Orthovector1.Zero();
     for (int i=NDIM;i--;) {
-      l = TopReplacement->x[i] - TopOrigin->x[i];
+      l = TopReplacement->at(i) - TopOrigin->at(i);
       if (fabs(l) > BondDistance) { // is component greater than bond distance
         Orthovector1[i] = (l < 0) ? -1. : +1.;
@@ -321 +347 @@
   InBondvector.Normalize();
   // get typical bond length and store as scale factor for later
-  ASSERT(TopOrigin->type != NULL, "AddHydrogenReplacementAtom: element of TopOrigin is not given.");
-  BondRescale = TopOrigin->type->HBondDistance[TopBond->BondDegree-1];
+  ASSERT(TopOrigin->getType() != NULL, "AddHydrogenReplacementAtom: element of TopOrigin is not given.");
+  BondRescale = TopOrigin->getType()->getHBondDistance(TopBond->BondDegree-1);
   if (BondRescale == -1) {
     DoeLog(1) && (eLog()<< Verbose(1) << "There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->getName() << "<->" << TopReplacement->getName() << ") of degree " << TopBond->BondDegree << "!" << endl);
@@ -336 +362 @@
     case 1:
       FirstOtherAtom = World::getInstance().createAtom();    // new atom
-      FirstOtherAtom->type = elemente->FindElement(1);  // element is Hydrogen
-      FirstOtherAtom->v = TopReplacement->v; // copy velocity
+      FirstOtherAtom->setType(1);  // element is Hydrogen
+      FirstOtherAtom->AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity
       FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
-      if (TopReplacement->type->Z == 1) { // neither rescale nor replace if it's already hydrogen
+      if (TopReplacement->getType()->getAtomicNumber() == 1) { // neither rescale nor replace if it's already hydrogen
         FirstOtherAtom->father = TopReplacement;
         BondRescale = bondlength;
@@ -346 +372 @@
       }
       InBondvector *= BondRescale;   // rescale the distance vector to Hydrogen bond length
-      FirstOtherAtom->x = TopOrigin->x; // set coordination to origin ...
-      FirstOtherAtom->x += InBondvector;  // ... and add distance vector to replacement atom
+      FirstOtherAtom->setPosition(TopOrigin->getPosition() + InBondvector); // set coordination to origin and add distance vector to replacement atom
       AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
 //      Log() << Verbose(4) << "Added " << *FirstOtherAtom << " at: ";
@@ -380 +405 @@
         // determine the plane of these two with the *origin
         try {
-          Orthovector1 =Plane(TopOrigin->x, FirstOtherAtom->x, SecondOtherAtom->x).getNormal();
+          Orthovector1 =Plane(TopOrigin->getPosition(), FirstOtherAtom->getPosition(), SecondOtherAtom->getPosition()).getNormal();
         }
         catch(LinearDependenceException &excp){
@@ -401 +426 @@
       FirstOtherAtom = World::getInstance().createAtom();
       SecondOtherAtom = World::getInstance().createAtom();
-      FirstOtherAtom->type = elemente->FindElement(1);
-      SecondOtherAtom->type = elemente->FindElement(1);
-      FirstOtherAtom->v = TopReplacement->v; // copy velocity
+      FirstOtherAtom->setType(1);
+      SecondOtherAtom->setType(1);
+      FirstOtherAtom->AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity
       FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
-      SecondOtherAtom->v = TopReplacement->v; // copy velocity
+      SecondOtherAtom->AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity
       SecondOtherAtom->FixedIon = TopReplacement->FixedIon;
       FirstOtherAtom->father = NULL;  // we are just an added hydrogen with no father
       SecondOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
-      bondangle = TopOrigin->type->HBondAngle[1];
+      bondangle = TopOrigin->getType()->getHBondAngle(1);
       if (bondangle == -1) {
         DoeLog(1) && (eLog()<< Verbose(1) << "There is no typical hydrogen bond angle in replacing bond (" << TopOrigin->getName() << "<->" << TopReplacement->getName() << ") of degree " << TopBond->BondDegree << "!" << endl);
@@ -423 +448 @@
 //      Log() << Verbose(0) << endl;
 //      Log() << Verbose(3) << "Half the bond angle is " << bondangle << ", sin and cos of it: " << sin(bondangle) << ", " << cos(bondangle) << endl;
-      FirstOtherAtom->x.Zero();
-      SecondOtherAtom->x.Zero();
+      FirstOtherAtom->Zero();
+      SecondOtherAtom->Zero();
       for(int i=NDIM;i--;) { // rotate by half the bond angle in both directions (InBondvector is bondangle = 0 direction)
-        FirstOtherAtom->x[i] = InBondvector[i] * cos(bondangle) + Orthovector1[i] * (sin(bondangle));
-        SecondOtherAtom->x[i] = InBondvector[i] * cos(bondangle) + Orthovector1[i] * (-sin(bondangle));
-      }
-      FirstOtherAtom->x *= BondRescale;  // rescale by correct BondDistance
-      SecondOtherAtom->x *= BondRescale;
+        FirstOtherAtom->set(i, InBondvector[i] * cos(bondangle) + Orthovector1[i] * (sin(bondangle)));
+        SecondOtherAtom->set(i, InBondvector[i] * cos(bondangle) + Orthovector1[i] * (-sin(bondangle)));
+      }
+      FirstOtherAtom->Scale(BondRescale);  // rescale by correct BondDistance
+      SecondOtherAtom->Scale(BondRescale);
       //Log() << Verbose(3) << "ReScaleCheck: " << FirstOtherAtom->x.Norm() << " and " << SecondOtherAtom->x.Norm() << "." << endl;
-      for(int i=NDIM;i--;) { // and make relative to origin atom
-        FirstOtherAtom->x[i] += TopOrigin->x[i];
-        SecondOtherAtom->x[i] += TopOrigin->x[i];
-      }
+      *FirstOtherAtom += TopOrigin->getPosition();
+      *SecondOtherAtom += TopOrigin->getPosition();
       // ... and add to molecule
       AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
@@ -457 +480 @@
       SecondOtherAtom = World::getInstance().createAtom();
       ThirdOtherAtom = World::getInstance().createAtom();
-      FirstOtherAtom->type = elemente->FindElement(1);
-      SecondOtherAtom->type = elemente->FindElement(1);
-      ThirdOtherAtom->type = elemente->FindElement(1);
-      FirstOtherAtom->v = TopReplacement->v; // copy velocity
+      FirstOtherAtom->setType(1);
+      SecondOtherAtom->setType(1);
+      ThirdOtherAtom->setType(1);
+      FirstOtherAtom->AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity
       FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
-      SecondOtherAtom->v = TopReplacement->v; // copy velocity
+      SecondOtherAtom->AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity
       SecondOtherAtom->FixedIon = TopReplacement->FixedIon;
-      ThirdOtherAtom->v = TopReplacement->v; // copy velocity
+      ThirdOtherAtom->AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity
       ThirdOtherAtom->FixedIon = TopReplacement->FixedIon;
       FirstOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
@@ -487 +510 @@
 
       // create correct coordination for the three atoms
-      alpha = (TopOrigin->type->HBondAngle[2])/180.*M_PI/2.;  // retrieve triple bond angle from database
+      alpha = (TopOrigin->getType()->getHBondAngle(2))/180.*M_PI/2.;  // retrieve triple bond angle from database
       l = BondRescale;        // desired bond length
       b = 2.*l*sin(alpha);    // base length of isosceles triangle
@@ -498 +521 @@
       factors[1] = f;
       factors[2] = 0.;
-      FirstOtherAtom->x.LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
+      FirstOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
       factors[1] = -0.5*f;
       factors[2] = g;
-      SecondOtherAtom->x.LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
+      SecondOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
       factors[2] = -g;
-      ThirdOtherAtom->x.LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
+      ThirdOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
 
       // rescale each to correct BondDistance
@@ -511 +534 @@
 
       // and relative to *origin atom
-      FirstOtherAtom->x += TopOrigin->x;
-      SecondOtherAtom->x += TopOrigin->x;
-      ThirdOtherAtom->x += TopOrigin->x;
+      *FirstOtherAtom += TopOrigin->getPosition();
+      *SecondOtherAtom += TopOrigin->getPosition();
+      *ThirdOtherAtom += TopOrigin->getPosition();
 
       // ... and add to molecule
@@ -589 +612 @@
     *item >> x[1];
     *item >> x[2];
-    Walker->type = elemente->FindElement(shorthand);
-    if (Walker->type == NULL) {
+    Walker->setType(elemente->FindElement(shorthand));
+    if (Walker->getType() == NULL) {
       DoeLog(1) && (eLog()<< Verbose(1) << "Could not parse the element at line: '" << line << "', setting to H.");
-      Walker->type = elemente->FindElement(1);
+      Walker->setType(1);
     }
     if (Walker->Trajectory.R.size() <= (unsigned int)MDSteps) {
@@ -599 +622 @@
       Walker->Trajectory.F.resize(MDSteps+10);
     }
+    Walker->setPosition(Vector(x));
     for(j=NDIM;j--;) {
-      Walker->x[j] = x[j];
       Walker->Trajectory.R.at(MDSteps-1)[j] = x[j];
       Walker->Trajectory.U.at(MDSteps-1)[j] = 0;
@@ -619 +642 @@
 {
   molecule *copy = World::getInstance().createMolecule();
-  atom *LeftAtom = NULL, *RightAtom = NULL;
 
   // copy all atoms
-  ActOnCopyWithEachAtom ( &molecule::AddCopyAtom, copy );
+  for_each(atoms.begin(),atoms.end(),bind1st(mem_fun(&molecule::AddCopyAtom),copy));
 
   // copy all bonds
-  bond *Binder = NULL;
-  bond *NewBond = NULL;
   for(molecule::iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner)
     for(BondList::iterator BondRunner = (*AtomRunner)->ListOfBonds.begin(); !(*AtomRunner)->ListOfBonds.empty(); BondRunner = (*AtomRunner)->ListOfBonds.begin())
       if ((*BondRunner)->leftatom == *AtomRunner) {
-        Binder = (*BondRunner);
+        bond *Binder = (*BondRunner);
 
         // get the pendant atoms of current bond in the copy molecule
-        copy->ActOnAllAtoms( &atom::EqualsFather, (const atom *)Binder->leftatom, (const atom **)&LeftAtom );
-        copy->ActOnAllAtoms( &atom::EqualsFather, (const atom *)Binder->rightatom, (const atom **)&RightAtom );
-
-        NewBond = copy->AddBond(LeftAtom, RightAtom, Binder->BondDegree);
+        atomSet::iterator leftiter=find_if(atoms.begin(),atoms.end(),bind2nd(mem_fun(&atom::isFather),Binder->leftatom));
+        atomSet::iterator rightiter=find_if(atoms.begin(),atoms.end(),bind2nd(mem_fun(&atom::isFather),Binder->rightatom));
+        ASSERT(leftiter!=atoms.end(),"No original left atom for bondcopy found");
+        ASSERT(leftiter!=atoms.end(),"No original right atom for bondcopy found");
+        atom *LeftAtom = *leftiter;
+        atom *RightAtom = *rightiter;
+
+        bond *NewBond = copy->AddBond(LeftAtom, RightAtom, Binder->BondDegree);
         NewBond->Cyclic = Binder->Cyclic;
         if (Binder->Cyclic)
@@ -643 +667 @@
       }
   // correct fathers
-  ActOnAllAtoms( &atom::CorrectFather );
+  for_each(atoms.begin(),atoms.end(),mem_fun(&atom::CorrectFather));
 
   // copy values
@@ -694 +718 @@
   atom1->RegisterBond(Binder);
   atom2->RegisterBond(Binder);
-  if ((atom1->type != NULL) && (atom1->type->Z != 1) && (atom2->type != NULL) && (atom2->type->Z != 1))
+  if ((atom1->getType() != NULL) && (atom1->getType()->getAtomicNumber() != 1) && (atom2->getType() != NULL) && (atom2->getType()->getAtomicNumber() != 1))
     NoNonBonds++;
 
@@ -760 +784 @@
 };
 
-/** Removes atom from molecule list and deletes it.
+/** Removes atom from molecule list and removes all of its bonds.
  * \param *pointer atom to be removed
  * \return true - succeeded, false - atom not found in list
@@ -768 +792 @@
   ASSERT(pointer, "Null pointer passed to molecule::RemoveAtom().");
   OBSERVE;
-  formula-=pointer->type;
   RemoveBonds(pointer);
   erase(pointer);
@@ -782 +805 @@
   if (pointer == NULL)
     return false;
-  formula-=pointer->type;
   erase(pointer);
   return true;
@@ -793 +815 @@
 {
   for (molecule::iterator iter = begin(); !empty(); iter = begin())
-      erase(iter);
+    erase(*iter);
   return empty();
 };
@@ -852 +874 @@
  * \param *out output stream
  */
-bool molecule::Output(ofstream * const output)
-{
-  int ElementNo[MAX_ELEMENTS], AtomNo[MAX_ELEMENTS];
-
-  for (int i=0;i<MAX_ELEMENTS;++i) {
-    AtomNo[i] = 0;
-    ElementNo[i] = 0;
-  }
+bool molecule::Output(ostream * const output)
+{
   if (output == NULL) {
     return false;
   } else {
+    int AtomNo[MAX_ELEMENTS];
+    memset(AtomNo,0,(MAX_ELEMENTS-1)*sizeof(*AtomNo));
+    enumeration<const element*> elementLookup = formula.enumerateElements();
+    for(map<const element*,unsigned int>::iterator iter=elementLookup.there.begin();
+        iter!=elementLookup.there.end();++iter){
+      cout << "Enumerated element " << *iter->first << " with number " << iter->second << endl;
+    }
     *output << "#Ion_TypeNr._Nr.R[0]    R[1]    R[2]    MoveType (0 MoveIon, 1 FixedIon)" << endl;
-    SetIndexedArrayForEachAtomTo ( ElementNo, &element::Z, &AbsoluteValue, 1);
-    int current=1;
-    for (int i=0;i<MAX_ELEMENTS;++i) {
-      if (ElementNo[i] == 1)
-        ElementNo[i] = current++;
-    }
-    ActOnAllAtoms( &atom::OutputArrayIndexed, (ostream * const) output, (const int *)ElementNo, (int *)AtomNo, (const char *) NULL );
+    for_each(atoms.begin(),atoms.end(),boost::bind(&atom::OutputArrayIndexed,_1,output,elementLookup,AtomNo,(const char*)0));
     return true;
   }
@@ -895 +912 @@
         ElementNo[i] = 0;
       }
-      SetIndexedArrayForEachAtomTo ( ElementNo, &element::Z, &AbsoluteValue, 1);
+      for(molecule::iterator iter = begin(); iter != end(); ++iter) {
+        ElementNo[(*iter)->getType()->getAtomicNumber()] = 1;
+      }
       int current=1;
       for (int i=0;i<MAX_ELEMENTS;++i) {
@@ -913 +932 @@
 {
   DoLog(2) && (Log() << Verbose(2) << endl << "From Contents of ListOfBonds, all non-hydrogen atoms:" << endl);
-  ActOnAllAtoms (&atom::OutputBondOfAtom );
+  for_each(atoms.begin(),atoms.end(),mem_fun(&atom::OutputBondOfAtom));
   DoLog(0) && (Log() << Verbose(0) << endl);
 };
@@ -936 +955 @@
     for (int step=0;step<MDSteps;step++) {
       *output << getAtomCount() << "\n\tCreated by molecuilder, step " << step << ", on " << ctime(&now);
-      ActOnAllAtoms( &atom::OutputTrajectoryXYZ, output, step );
+      for_each(atoms.begin(),atoms.end(),boost::bind(&atom::OutputTrajectoryXYZ,_1,output,step));
     }
     return true;
@@ -953 +972 @@
     now = time((time_t *)NULL);   // Get the system time and put it into 'now' as 'calender time'
     *output << getAtomCount() << "\n\tCreated by molecuilder on " << ctime(&now);
-    ActOnAllAtoms( &atom::OutputXYZLine, output );
+    for_each(atoms.begin(),atoms.end(),bind2nd(mem_fun(&atom::OutputXYZLine),output));
     return true;
   } else
@@ -969 +988 @@
   for (molecule::const_iterator iter = atoms.begin(); iter != atoms.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
+    if ((*iter)->getType()->getAtomicNumber() != 1) // count non-hydrogen atoms whilst at it
       NoNonHydrogen++;
     stringstream sstr;
-    sstr << (*iter)->type->symbol << (*iter)->nr+1;
+    sstr << (*iter)->getType()->getSymbol() << (*iter)->nr+1;
     (*iter)->setName(sstr.str());
     DoLog(3) && (Log() << Verbose(3) << "Naming atom nr. " << (*iter)->nr << " " << (*iter)->getName() << "." << endl);
@@ -978 +997 @@
   }
   return res;
-};
-
-/** 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;
-  }
 };