Changes in src/analysis_correlation.cpp [112b09:36166d]

File:

• src/analysis_correlation.cpp (modified) (16 diffs)

src/analysis_correlation.cpp

@@ -9 +9 @@
 
 #include <iostream>
+#include <iomanip>
 
 #include "analysis_correlation.hpp"
@@ -19 +20 @@
 #include "triangleintersectionlist.hpp"
 #include "vector.hpp"
+#include "Matrix.hpp"
 #include "verbose.hpp"
 #include "World.hpp"
+#include "Box.hpp"
 
 
 /** Calculates the pair correlation between given elements.
  * Note given element order is unimportant (i.e. g(Si, O) === g(O, Si))
- * \param *out output stream for debugging
- * \param *molecules list of molecules structure
- * \param *type1 first element or NULL (if any element)
- * \param *type2 second element or NULL (if any element)
+ * \param *molecules list of molecules structure
+ * \param &elements vector of elements to correlate
  * \return Map of doubles with values the pair of the two atoms.
  */
-PairCorrelationMap *PairCorrelation(MoleculeListClass * const &molecules, const element * const type1, const element * const type2 )
+PairCorrelationMap *PairCorrelation(MoleculeListClass * const &molecules, const std::vector<element *> &elements)
 {
   Info FunctionInfo(__func__);
   PairCorrelationMap *outmap = NULL;
   double distance = 0.;
+  Box &domain = World::getInstance().getDomain();
 
   if (molecules->ListOfMolecules.empty()) {
@@ -43 +45 @@
   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
     (*MolWalker)->doCountAtoms();
+
+  // create all possible pairs of elements
+  set <pair<element *, element *> > PairsOfElements;
+  if (elements.size() >= 2) {
+    for (vector<element *>::const_iterator type1 = elements.begin(); type1 != elements.end(); ++type1)
+      for (vector<element *>::const_iterator type2 = elements.begin(); type2 != elements.end(); ++type2)
+        if (type1 != type2) {
+          PairsOfElements.insert( pair<element *, element*>(*type1,*type2) );
+          DoLog(1) && (Log() << Verbose(1) << "Creating element pair " << (*type1)->symbol << " and " << (*type2)->symbol << "." << endl);
+        }
+  } else if (elements.size() == 1) { // one to all are valid
+    element *elemental = *elements.begin();
+    PairsOfElements.insert( pair<element *, element*>(elemental,(element *)NULL) );
+    PairsOfElements.insert( pair<element *, element*>((element *)NULL,elemental) );
+  } else { // all elements valid
+    PairsOfElements.insert( pair<element *, element*>((element *)NULL, (element *)NULL) );
+  }
+
   outmap = new PairCorrelationMap;
   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++){
@@ -50 +70 @@
       for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
         DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
-        if ((type1 == NULL) || ((*iter)->type == type1)) {
-          for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++){
-            if ((*MolOtherWalker)->ActiveFlag) {
-              DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
-              for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
-                DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << **runner << "." << endl);
-                if ((*iter)->getId() < (*runner)->getId()){
-                  if ((type2 == NULL) || ((*runner)->type == type2)) {
-                    distance = (*iter)->node->PeriodicDistance(*(*runner)->node,  World::getInstance().getDomain());
+        for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++){
+          if ((*MolOtherWalker)->ActiveFlag) {
+            DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
+            for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
+              DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << **runner << "." << endl);
+              if ((*iter)->getId() < (*runner)->getId()){
+                for (set <pair<element *, element *> >::iterator PairRunner = PairsOfElements.begin(); PairRunner != PairsOfElements.end(); ++PairRunner)
+                  if ((PairRunner->first == (**iter).type) && (PairRunner->second == (**runner).type)) {
+                    distance = domain.periodicDistance(*(*iter)->node,*(*runner)->node);
                     //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
                     outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> ((*iter), (*runner)) ) );
                   }
-                }
               }
             }
@@ -75 +94 @@
 /** Calculates the pair correlation between given elements.
  * Note given element order is unimportant (i.e. g(Si, O) === g(O, Si))
- * \param *out output stream for debugging
- * \param *molecules list of molecules structure
- * \param *type1 first element or NULL (if any element)
- * \param *type2 second element or NULL (if any element)
+ * \param *molecules list of molecules structure
+ * \param &elements vector of elements to correlate
  * \param ranges[NDIM] interval boundaries for the periodic images to scan also
  * \return Map of doubles with values the pair of the two atoms.
  */
-PairCorrelationMap *PeriodicPairCorrelation(MoleculeListClass * const &molecules, const element * const type1, const element * const type2, const int ranges[NDIM] )
+PairCorrelationMap *PeriodicPairCorrelation(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const int ranges[NDIM] )
 {
   Info FunctionInfo(__func__);
@@ -100 +117 @@
   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
     (*MolWalker)->doCountAtoms();
+
+  // create all possible pairs of elements
+  set <pair<element *, element *> > PairsOfElements;
+  if (elements.size() >= 2) {
+    for (vector<element *>::const_iterator type1 = elements.begin(); type1 != elements.end(); ++type1)
+      for (vector<element *>::const_iterator type2 = elements.begin(); type2 != elements.end(); ++type2)
+        if (type1 != type2) {
+          PairsOfElements.insert( pair<element *, element*>(*type1,*type2) );
+          DoLog(1) && (Log() << Verbose(1) << "Creating element pair " << (*type1)->symbol << " and " << (*type2)->symbol << "." << endl);
+        }
+  } else if (elements.size() == 1) { // one to all are valid
+    element *elemental = *elements.begin();
+    PairsOfElements.insert( pair<element *, element*>(elemental,(element *)NULL) );
+    PairsOfElements.insert( pair<element *, element*>((element *)NULL,elemental) );
+  } else { // all elements valid
+    PairsOfElements.insert( pair<element *, element*>((element *)NULL, (element *)NULL) );
+  }
+
   outmap = new PairCorrelationMap;
-  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
-    if ((*MolWalker)->ActiveFlag) {
-      double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
-      double * FullInverseMatrix = InverseMatrix(FullMatrix);
+  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++){
+    if ((*MolWalker)->ActiveFlag) {
+      Matrix FullMatrix = World::getInstance().getDomain().getM();
+      Matrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
       DoeLog(2) && (eLog()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
+      eLog() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
       for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
         DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
-        if ((type1 == NULL) || ((*iter)->type == type1)) {
-          periodicX = *(*iter)->node;
-          periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
-          // go through every range in xyz and get distance
-          for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
-            for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
-              for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
-                checkX = Vector(n[0], n[1], n[2]) + periodicX;
-                checkX.MatrixMultiplication(FullMatrix);
-                for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++)
-                  if ((*MolOtherWalker)->ActiveFlag) {
-                    DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
-                    for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
-                      DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << **runner << "." << endl);
-                      if ((*iter)->nr < (*runner)->nr)
-                        if ((type2 == NULL) || ((*runner)->type == type2)) {
-                          periodicOtherX = *(*runner)->node;
-                          periodicOtherX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
+        periodicX = FullInverseMatrix * (*(*iter)->node); // x now in [0,1)^3
+        // go through every range in xyz and get distance
+        for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
+          for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
+            for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
+              checkX = FullMatrix * (Vector(n[0], n[1], n[2]) + periodicX);
+              for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++){
+                if ((*MolOtherWalker)->ActiveFlag) {
+                  DoLog(2) && (Log() << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl);
+                  for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
+                    DoLog(3) && (Log() << Verbose(3) << "Current otheratom is " << **runner << "." << endl);
+                    if ((*iter)->getId() < (*runner)->getId()){
+                      for (set <pair<element *, element *> >::iterator PairRunner = PairsOfElements.begin(); PairRunner != PairsOfElements.end(); ++PairRunner)
+                        if ((PairRunner->first == (**iter).type) && (PairRunner->second == (**runner).type)) {
+                          periodicOtherX = FullInverseMatrix * (*(*runner)->node); // x now in [0,1)^3
                           // go through every range in xyz and get distance
                           for (Othern[0]=-ranges[0]; Othern[0] <= ranges[0]; Othern[0]++)
                             for (Othern[1]=-ranges[1]; Othern[1] <= ranges[1]; Othern[1]++)
                               for (Othern[2]=-ranges[2]; Othern[2] <= ranges[2]; Othern[2]++) {
-                                checkOtherX = Vector(Othern[0], Othern[1], Othern[2]) + periodicOtherX;
-                                checkOtherX.MatrixMultiplication(FullMatrix);
+                                checkOtherX = FullMatrix * (Vector(Othern[0], Othern[1], Othern[2]) + periodicOtherX);
                                 distance = checkX.distance(checkOtherX);
                                 //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
@@ -137 +169 @@
                               }
                         }
+                    }
                   }
+                }
               }
             }
-        }
-      }
-      delete[](FullMatrix);
-      delete[](FullInverseMatrix);
-    }
+      }
+    }
+  }
 
   return outmap;
@@ -150 +182 @@
 
 /** Calculates the distance (pair) correlation between a given element and a point.
- * \param *out output stream for debugging
- * \param *molecules list of molecules structure
- * \param *type element or NULL (if any element)
+ * \param *molecules list of molecules structure
+ * \param &elements vector of elements to correlate with point
  * \param *point vector to the correlation point
  * \return Map of dobules with values as pairs of atom and the vector
  */
-CorrelationToPointMap *CorrelationToPoint(MoleculeListClass * const &molecules, const element * const type, const Vector *point )
+CorrelationToPointMap *CorrelationToPoint(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const Vector *point )
 {
   Info FunctionInfo(__func__);
   CorrelationToPointMap *outmap = NULL;
   double distance = 0.;
+  Box &domain = World::getInstance().getDomain();
 
   if (molecules->ListOfMolecules.empty()) {
@@ -174 +206 @@
       for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
         DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
-        if ((type == NULL) || ((*iter)->type == type)) {
-          distance = (*iter)->node->PeriodicDistance(*point, World::getInstance().getDomain());
-          DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
-          outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> ((*iter), point) ) );
-        }
+        for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
+          if ((*type == NULL) || ((*iter)->type == *type)) {
+            distance = domain.periodicDistance(*(*iter)->node,*point);
+            DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
+            outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> ((*iter), point) ) );
+          }
       }
     }
@@ -186 +219 @@
 
 /** Calculates the distance (pair) correlation between a given element, all its periodic images and a point.
- * \param *out output stream for debugging
- * \param *molecules list of molecules structure
- * \param *type element or NULL (if any element)
+ * \param *molecules list of molecules structure
+ * \param &elements vector of elements to correlate to point
  * \param *point vector to the correlation point
  * \param ranges[NDIM] interval boundaries for the periodic images to scan also
  * \return Map of dobules with values as pairs of atom and the vector
  */
-CorrelationToPointMap *PeriodicCorrelationToPoint(MoleculeListClass * const &molecules, const element * const type, const Vector *point, const int ranges[NDIM] )
+CorrelationToPointMap *PeriodicCorrelationToPoint(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const Vector *point, const int ranges[NDIM] )
 {
   Info FunctionInfo(__func__);
@@ -211 +243 @@
   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
     if ((*MolWalker)->ActiveFlag) {
-      double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
-      double * FullInverseMatrix = InverseMatrix(FullMatrix);
+      Matrix FullMatrix = World::getInstance().getDomain().getM();
+      Matrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
       DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
       for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
         DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
-        if ((type == NULL) || ((*iter)->type == type)) {
-          periodicX = *(*iter)->node;
-          periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
-          // go through every range in xyz and get distance
-          for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
-            for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
-              for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
-                checkX = Vector(n[0], n[1], n[2]) + periodicX;
-                checkX.MatrixMultiplication(FullMatrix);
-                distance = checkX.distance(*point);
-                DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
-                outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (*iter, point) ) );
-              }
-        }
-      }
-      delete[](FullMatrix);
-      delete[](FullInverseMatrix);
+        for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
+          if ((*type == NULL) || ((*iter)->type == *type)) {
+            periodicX = FullInverseMatrix * (*(*iter)->node); // x now in [0,1)^3
+            // go through every range in xyz and get distance
+            for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
+              for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
+                for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
+                  checkX = FullMatrix * (Vector(n[0], n[1], n[2]) + periodicX);
+                  distance = checkX.distance(*point);
+                  DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
+                  outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (*iter, point) ) );
+                }
+          }
+      }
     }
 
@@ -239 +268 @@
 
 /** Calculates the distance (pair) correlation between a given element and a surface.
- * \param *out output stream for debugging
- * \param *molecules list of molecules structure
- * \param *type element or NULL (if any element)
+ * \param *molecules list of molecules structure
+ * \param &elements vector of elements to correlate to surface
  * \param *Surface pointer to Tesselation class surface
  * \param *LC LinkedCell structure to quickly find neighbouring atoms
  * \return Map of doubles with values as pairs of atom and the BoundaryTriangleSet that's closest
  */
-CorrelationToSurfaceMap *CorrelationToSurface(MoleculeListClass * const &molecules, const element * const type, const Tesselation * const Surface, const LinkedCell *LC )
+CorrelationToSurfaceMap *CorrelationToSurface(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const Tesselation * const Surface, const LinkedCell *LC )
 {
   Info FunctionInfo(__func__);
@@ -268 +296 @@
       for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
         DoLog(1) && (Log() << Verbose(1) << "\tCurrent atom is " << *(*iter) << "." << endl);
-        if ((type == NULL) || ((*iter)->type == type)) {
-          TriangleIntersectionList Intersections((*iter)->node,Surface,LC);
-          distance = Intersections.GetSmallestDistance();
-          triangle = Intersections.GetClosestTriangle();
-          outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> ((*iter), triangle) ) );
-        }
+        for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
+          if ((*type == NULL) || ((*iter)->type == *type)) {
+            TriangleIntersectionList Intersections((*iter)->node,Surface,LC);
+            distance = Intersections.GetSmallestDistance();
+            triangle = Intersections.GetClosestTriangle();
+            outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> ((*iter), triangle) ) );
+          }
       }
     } else {
@@ -287 +316 @@
  * axis an integer from [ -ranges[i], ranges[i] ] onto it and multiply with the domain matrix to bring it back into
  * the real space. Then, we Tesselation::FindClosestTriangleToPoint() and DistanceToTrianglePlane().
- * \param *out output stream for debugging
- * \param *molecules list of molecules structure
- * \param *type element or NULL (if any element)
+ * \param *molecules list of molecules structure
+ * \param &elements vector of elements to correlate to surface
  * \param *Surface pointer to Tesselation class surface
  * \param *LC LinkedCell structure to quickly find neighbouring atoms
@@ -295 +323 @@
  * \return Map of doubles with values as pairs of atom and the BoundaryTriangleSet that's closest
  */
-CorrelationToSurfaceMap *PeriodicCorrelationToSurface(MoleculeListClass * const &molecules, const element * const type, const Tesselation * const Surface, const LinkedCell *LC, const int ranges[NDIM] )
+CorrelationToSurfaceMap *PeriodicCorrelationToSurface(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const Tesselation * const Surface, const LinkedCell *LC, const int ranges[NDIM] )
 {
   Info FunctionInfo(__func__);
@@ -317 +345 @@
   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
     if ((*MolWalker)->ActiveFlag) {
-      double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
-      double * FullInverseMatrix = InverseMatrix(FullMatrix);
+      Matrix FullMatrix = World::getInstance().getDomain().getM();
+      Matrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
       DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
       for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
         DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
-        if ((type == NULL) || ((*iter)->type == type)) {
-          periodicX = *(*iter)->node;
-          periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
-          // go through every range in xyz and get distance
-          ShortestDistance = -1.;
-          for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
-            for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
-              for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
-                checkX = Vector(n[0], n[1], n[2]) + periodicX;
-                checkX.MatrixMultiplication(FullMatrix);
-                TriangleIntersectionList Intersections(&checkX,Surface,LC);
-                distance = Intersections.GetSmallestDistance();
-                triangle = Intersections.GetClosestTriangle();
-                if ((ShortestDistance == -1.) || (distance < ShortestDistance)) {
-                  ShortestDistance = distance;
-                  ShortestTriangle = triangle;
+        for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
+          if ((*type == NULL) || ((*iter)->type == *type)) {
+            periodicX = FullInverseMatrix * (*(*iter)->node); // x now in [0,1)^3
+            // go through every range in xyz and get distance
+            ShortestDistance = -1.;
+            for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
+              for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
+                for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
+                  checkX = FullMatrix * (Vector(n[0], n[1], n[2]) + periodicX);
+                  TriangleIntersectionList Intersections(&checkX,Surface,LC);
+                  distance = Intersections.GetSmallestDistance();
+                  triangle = Intersections.GetClosestTriangle();
+                  if ((ShortestDistance == -1.) || (distance < ShortestDistance)) {
+                    ShortestDistance = distance;
+                    ShortestTriangle = triangle;
+                  }
                 }
-              }
-          // insert
-          outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(ShortestDistance, pair<atom *, BoundaryTriangleSet*> (*iter, ShortestTriangle) ) );
-          //Log() << Verbose(1) << "INFO: Inserting " << Walker << " with distance " << ShortestDistance << " to " << *ShortestTriangle << "." << endl;
-        }
-      }
-      delete[](FullMatrix);
-      delete[](FullInverseMatrix);
+            // insert
+            outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(ShortestDistance, pair<atom *, BoundaryTriangleSet*> (*iter, ShortestTriangle) ) );
+            //Log() << Verbose(1) << "INFO: Inserting " << Walker << " with distance " << ShortestDistance << " to " << *ShortestTriangle << "." << endl;
+          }
+      }
     }