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

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
  *
@@ -12 +5 @@
  */
 
-// include config.h
 #ifdef HAVE_CONFIG_H
 #include <config.h>
@@ -97 +89 @@
 }
 
-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;
@@ -172 +151 @@
 molecule::const_iterator molecule::erase( const_iterator loc )
 {
-  OBSERVE;
   molecule::const_iterator iter = loc;
   iter--;
@@ -178 +156 @@
   atomIds.erase( atom->getId() );
   atoms.remove( atom );
-  formula-=atom->getType();
   atom->removeFromMolecule();
   return iter;
@@ -185 +162 @@
 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();
   }
@@ -208 +183 @@
 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 {
@@ -233 +206 @@
   if (pointer != NULL) {
     pointer->sort = &pointer->nr;
-    if (pointer->getType() != NULL) {
-      if (pointer->getType()->getAtomicNumber() != 1)
+    if (pointer->type != NULL) {
+      formula += pointer->type;
+      if (pointer->type->Z != 1)
         NoNonHydrogen++;
       if(pointer->getName() == "Unknown"){
         stringstream sstr;
-        sstr << pointer->getType()->getSymbol() << pointer->nr+1;
+        sstr << pointer->type->symbol << pointer->nr+1;
         pointer->setName(sstr.str());
       }
@@ -262 +236 @@
     walker->nr = last_atom++;  // increase number within molecule
     insert(walker);
-    if ((pointer->getType() != NULL) && (pointer->getType()->getAtomicNumber() != 1))
+    if ((pointer->type != NULL) && (pointer->type->Z != 1))
       NoNonHydrogen++;
     retval=walker;
@@ -320 +294 @@
 //  Log() << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl;
   // create vector in direction of bond
-  InBondvector = TopReplacement->getPosition() - TopOrigin->getPosition();
+  InBondvector = TopReplacement->x - TopOrigin->x;
   bondlength = InBondvector.Norm();
 
@@ -332 +306 @@
     Orthovector1.Zero();
     for (int i=NDIM;i--;) {
-      l = TopReplacement->at(i) - TopOrigin->at(i);
+      l = TopReplacement->x[i] - TopOrigin->x[i];
       if (fabs(l) > BondDistance) { // is component greater than bond distance
         Orthovector1[i] = (l < 0) ? -1. : +1.;
@@ -347 +321 @@
   InBondvector.Normalize();
   // get typical bond length and store as scale factor for later
-  ASSERT(TopOrigin->getType() != NULL, "AddHydrogenReplacementAtom: element of TopOrigin is not given.");
-  BondRescale = TopOrigin->getType()->getHBondDistance(TopBond->BondDegree-1);
+  ASSERT(TopOrigin->type != NULL, "AddHydrogenReplacementAtom: element of TopOrigin is not given.");
+  BondRescale = TopOrigin->type->HBondDistance[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);
@@ -362 +336 @@
     case 1:
       FirstOtherAtom = World::getInstance().createAtom();    // new atom
-      FirstOtherAtom->setType(1);  // element is Hydrogen
-      FirstOtherAtom->AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity
+      FirstOtherAtom->type = elemente->FindElement(1);  // element is Hydrogen
+      FirstOtherAtom->v = TopReplacement->v; // copy velocity
       FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
-      if (TopReplacement->getType()->getAtomicNumber() == 1) { // neither rescale nor replace if it's already hydrogen
+      if (TopReplacement->type->Z == 1) { // neither rescale nor replace if it's already hydrogen
         FirstOtherAtom->father = TopReplacement;
         BondRescale = bondlength;
@@ -372 +346 @@
       }
       InBondvector *= BondRescale;   // rescale the distance vector to Hydrogen bond length
-      FirstOtherAtom->setPosition(TopOrigin->getPosition() + InBondvector); // set coordination to origin and add distance vector to replacement atom
+      FirstOtherAtom->x = TopOrigin->x; // set coordination to origin ...
+      FirstOtherAtom->x += InBondvector;  // ... and add distance vector to replacement atom
       AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
 //      Log() << Verbose(4) << "Added " << *FirstOtherAtom << " at: ";
@@ -405 +380 @@
         // determine the plane of these two with the *origin
         try {
-          Orthovector1 =Plane(TopOrigin->getPosition(), FirstOtherAtom->getPosition(), SecondOtherAtom->getPosition()).getNormal();
+          Orthovector1 =Plane(TopOrigin->x, FirstOtherAtom->x, SecondOtherAtom->x).getNormal();
         }
         catch(LinearDependenceException &excp){
@@ -426 +401 @@
       FirstOtherAtom = World::getInstance().createAtom();
       SecondOtherAtom = World::getInstance().createAtom();
-      FirstOtherAtom->setType(1);
-      SecondOtherAtom->setType(1);
-      FirstOtherAtom->AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity
+      FirstOtherAtom->type = elemente->FindElement(1);
+      SecondOtherAtom->type = elemente->FindElement(1);
+      FirstOtherAtom->v = TopReplacement->v; // copy velocity
       FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
-      SecondOtherAtom->AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity
+      SecondOtherAtom->v = TopReplacement->v; // 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->getType()->getHBondAngle(1);
+      bondangle = TopOrigin->type->HBondAngle[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);
@@ -448 +423 @@
 //      Log() << Verbose(0) << endl;
 //      Log() << Verbose(3) << "Half the bond angle is " << bondangle << ", sin and cos of it: " << sin(bondangle) << ", " << cos(bondangle) << endl;
-      FirstOtherAtom->Zero();
-      SecondOtherAtom->Zero();
+      FirstOtherAtom->x.Zero();
+      SecondOtherAtom->x.Zero();
       for(int i=NDIM;i--;) { // rotate by half the bond angle in both directions (InBondvector is bondangle = 0 direction)
-        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);
+        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;
       //Log() << Verbose(3) << "ReScaleCheck: " << FirstOtherAtom->x.Norm() << " and " << SecondOtherAtom->x.Norm() << "." << endl;
-      *FirstOtherAtom += TopOrigin->getPosition();
-      *SecondOtherAtom += TopOrigin->getPosition();
+      for(int i=NDIM;i--;) { // and make relative to origin atom
+        FirstOtherAtom->x[i] += TopOrigin->x[i];
+        SecondOtherAtom->x[i] += TopOrigin->x[i];
+      }
       // ... and add to molecule
       AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
@@ -480 +457 @@
       SecondOtherAtom = World::getInstance().createAtom();
       ThirdOtherAtom = World::getInstance().createAtom();
-      FirstOtherAtom->setType(1);
-      SecondOtherAtom->setType(1);
-      ThirdOtherAtom->setType(1);
-      FirstOtherAtom->AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity
+      FirstOtherAtom->type = elemente->FindElement(1);
+      SecondOtherAtom->type = elemente->FindElement(1);
+      ThirdOtherAtom->type = elemente->FindElement(1);
+      FirstOtherAtom->v = TopReplacement->v; // copy velocity
       FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
-      SecondOtherAtom->AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity
+      SecondOtherAtom->v = TopReplacement->v; // copy velocity
       SecondOtherAtom->FixedIon = TopReplacement->FixedIon;
-      ThirdOtherAtom->AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity
+      ThirdOtherAtom->v = TopReplacement->v; // copy velocity
       ThirdOtherAtom->FixedIon = TopReplacement->FixedIon;
       FirstOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
@@ -510 +487 @@
 
       // create correct coordination for the three atoms
-      alpha = (TopOrigin->getType()->getHBondAngle(2))/180.*M_PI/2.;  // retrieve triple bond angle from database
+      alpha = (TopOrigin->type->HBondAngle[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
@@ -521 +498 @@
       factors[1] = f;
       factors[2] = 0.;
-      FirstOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
+      FirstOtherAtom->x.LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
       factors[1] = -0.5*f;
       factors[2] = g;
-      SecondOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
+      SecondOtherAtom->x.LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
       factors[2] = -g;
-      ThirdOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
+      ThirdOtherAtom->x.LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
 
       // rescale each to correct BondDistance
@@ -534 +511 @@
 
       // and relative to *origin atom
-      *FirstOtherAtom += TopOrigin->getPosition();
-      *SecondOtherAtom += TopOrigin->getPosition();
-      *ThirdOtherAtom += TopOrigin->getPosition();
+      FirstOtherAtom->x += TopOrigin->x;
+      SecondOtherAtom->x += TopOrigin->x;
+      ThirdOtherAtom->x += TopOrigin->x;
 
       // ... and add to molecule
@@ -612 +589 @@
     *item >> x[1];
     *item >> x[2];
-    Walker->setType(elemente->FindElement(shorthand));
-    if (Walker->getType() == NULL) {
+    Walker->type = elemente->FindElement(shorthand);
+    if (Walker->type == NULL) {
       DoeLog(1) && (eLog()<< Verbose(1) << "Could not parse the element at line: '" << line << "', setting to H.");
-      Walker->setType(1);
+      Walker->type = elemente->FindElement(1);
     }
     if (Walker->Trajectory.R.size() <= (unsigned int)MDSteps) {
@@ -622 +599 @@
       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;
@@ -642 +619 @@
 {
   molecule *copy = World::getInstance().createMolecule();
+  atom *LeftAtom = NULL, *RightAtom = NULL;
 
   // copy all atoms
-  for_each(atoms.begin(),atoms.end(),bind1st(mem_fun(&molecule::AddCopyAtom),copy));
+  ActOnCopyWithEachAtom ( &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) {
-        bond *Binder = (*BondRunner);
+        Binder = (*BondRunner);
 
         // get the pendant atoms of current bond in the copy molecule
-        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);
+        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);
         NewBond->Cyclic = Binder->Cyclic;
         if (Binder->Cyclic)
@@ -667 +643 @@
       }
   // correct fathers
-  for_each(atoms.begin(),atoms.end(),mem_fun(&atom::CorrectFather));
+  ActOnAllAtoms( &atom::CorrectFather );
 
   // copy values
@@ -718 +694 @@
   atom1->RegisterBond(Binder);
   atom2->RegisterBond(Binder);
-  if ((atom1->getType() != NULL) && (atom1->getType()->getAtomicNumber() != 1) && (atom2->getType() != NULL) && (atom2->getType()->getAtomicNumber() != 1))
+  if ((atom1->type != NULL) && (atom1->type->Z != 1) && (atom2->type != NULL) && (atom2->type->Z != 1))
     NoNonBonds++;
 
@@ -784 +760 @@
 };
 
-/** Removes atom from molecule list and removes all of its bonds.
+/** Removes atom from molecule list and deletes it.
  * \param *pointer atom to be removed
  * \return true - succeeded, false - atom not found in list
@@ -792 +768 @@
   ASSERT(pointer, "Null pointer passed to molecule::RemoveAtom().");
   OBSERVE;
+  formula-=pointer->type;
   RemoveBonds(pointer);
   erase(pointer);
@@ -805 +782 @@
   if (pointer == NULL)
     return false;
+  formula-=pointer->type;
   erase(pointer);
   return true;
@@ -815 +793 @@
 {
   for (molecule::iterator iter = begin(); !empty(); iter = begin())
-    erase(*iter);
+      erase(iter);
   return empty();
 };
@@ -874 +852 @@
  * \param *out output stream
  */
-bool molecule::Output(ostream * const output)
-{
+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;
+  }
   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;
-    for_each(atoms.begin(),atoms.end(),boost::bind(&atom::OutputArrayIndexed,_1,output,elementLookup,AtomNo,(const char*)0));
+    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 );
     return true;
   }
@@ -912 +895 @@
         ElementNo[i] = 0;
       }
-      for(molecule::iterator iter = begin(); iter != end(); ++iter) {
-        ElementNo[(*iter)->getType()->getAtomicNumber()] = 1;
-      }
+      SetIndexedArrayForEachAtomTo ( ElementNo, &element::Z, &AbsoluteValue, 1);
       int current=1;
       for (int i=0;i<MAX_ELEMENTS;++i) {
@@ -932 +913 @@
 {
   DoLog(2) && (Log() << Verbose(2) << endl << "From Contents of ListOfBonds, all non-hydrogen atoms:" << endl);
-  for_each(atoms.begin(),atoms.end(),mem_fun(&atom::OutputBondOfAtom));
+  ActOnAllAtoms (&atom::OutputBondOfAtom );
   DoLog(0) && (Log() << Verbose(0) << endl);
 };
@@ -955 +936 @@
     for (int step=0;step<MDSteps;step++) {
       *output << getAtomCount() << "\n\tCreated by molecuilder, step " << step << ", on " << ctime(&now);
-      for_each(atoms.begin(),atoms.end(),boost::bind(&atom::OutputTrajectoryXYZ,_1,output,step));
+      ActOnAllAtoms( &atom::OutputTrajectoryXYZ, output, step );
     }
     return true;
@@ -972 +953 @@
     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);
-    for_each(atoms.begin(),atoms.end(),bind2nd(mem_fun(&atom::OutputXYZLine),output));
+    ActOnAllAtoms( &atom::OutputXYZLine, output );
     return true;
   } else
@@ -988 +969 @@
   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)->getType()->getAtomicNumber() != 1) // count non-hydrogen atoms whilst at it
+    if ((*iter)->type->Z != 1) // count non-hydrogen atoms whilst at it
       NoNonHydrogen++;
     stringstream sstr;
-    sstr << (*iter)->getType()->getSymbol() << (*iter)->nr+1;
+    sstr << (*iter)->type->symbol << (*iter)->nr+1;
     (*iter)->setName(sstr.str());
     DoLog(3) && (Log() << Verbose(3) << "Naming atom nr. " << (*iter)->nr << " " << (*iter)->getName() << "." << endl);
@@ -997 +978 @@
   }
   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;
+  }
 };