Changes in src/analysis_correlation.cpp [58bbd3:112b09]

File:

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

src/analysis_correlation.cpp

@@ -5 +5 @@
  *      Author: heber
  */
+
+#include "Helpers/MemDebug.hpp"
 
 #include <iostream>
@@ -23 +25 @@
 /** Calculates the pair correlation between given elements.
  * Note given element order is unimportant (i.e. g(Si, O) === g(O, Si))
- * \param *molecules list of molecules structure
- * \param &elements vector of elements to correlate
+ * \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)
  * \return Map of doubles with values the pair of the two atoms.
  */
-PairCorrelationMap *PairCorrelation(MoleculeListClass * const &molecules, const std::vector<element *> &elements)
+PairCorrelationMap *PairCorrelation(MoleculeListClass * const &molecules, const element * const type1, const element * const type2 )
 {
   Info FunctionInfo(__func__);
   PairCorrelationMap *outmap = NULL;
   double distance = 0.;
-  double *domain = World::getInstance().getDomain();
 
   if (molecules->ListOfMolecules.empty()) {
@@ -40 +43 @@
   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++){
@@ -65 +50 @@
       for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
         DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
-        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 = (*iter)->node->PeriodicDistance(*(*runner)->node,  domain);
+        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());
                     //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
                     outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> ((*iter), (*runner)) ) );
                   }
+                }
               }
             }
@@ -89 +75 @@
 /** Calculates the pair correlation between given elements.
  * Note given element order is unimportant (i.e. g(Si, O) === g(O, Si))
- * \param *molecules list of molecules structure
- * \param &elements vector of elements to correlate
+ * \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 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 std::vector<element *> &elements, const int ranges[NDIM] )
+PairCorrelationMap *PeriodicPairCorrelation(MoleculeListClass * const &molecules, const element * const type1, const element * const type2, const int ranges[NDIM] )
 {
   Info FunctionInfo(__func__);
@@ -112 +100 @@
   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++){
+  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);
       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);
-        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)->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)) {
+        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
@@ -168 +137 @@
                               }
                         }
-                    }
                   }
-                }
               }
             }
+        }
       }
       delete[](FullMatrix);
       delete[](FullInverseMatrix);
     }
-  }
 
   return outmap;
@@ -183 +150 @@
 
 /** Calculates the distance (pair) correlation between a given element and a point.
- * \param *molecules list of molecules structure
- * \param &elements vector of elements to correlate with point
+ * \param *out output stream for debugging
+ * \param *molecules list of molecules structure
+ * \param *type element or NULL (if any element)
  * \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 std::vector<element *> &elements, const Vector *point )
+CorrelationToPointMap *CorrelationToPoint(MoleculeListClass * const &molecules, const element * const type, const Vector *point )
 {
   Info FunctionInfo(__func__);
   CorrelationToPointMap *outmap = NULL;
   double distance = 0.;
-  double *cell_size = World::getInstance().getDomain();
 
   if (molecules->ListOfMolecules.empty()) {
@@ -207 +174 @@
       for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
         DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
-        for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
-          if ((*type == NULL) || ((*iter)->type == *type)) {
-            distance = (*iter)->node->PeriodicDistance(*point, cell_size);
-            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) ) );
-          }
+        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) ) );
+        }
       }
     }
@@ -220 +186 @@
 
 /** Calculates the distance (pair) correlation between a given element, all its periodic images and a point.
- * \param *molecules list of molecules structure
- * \param &elements vector of elements to correlate to point
+ * \param *out output stream for debugging
+ * \param *molecules list of molecules structure
+ * \param *type element or NULL (if any element)
  * \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 std::vector<element *> &elements, const Vector *point, const int ranges[NDIM] )
+CorrelationToPointMap *PeriodicCorrelationToPoint(MoleculeListClass * const &molecules, const element * const type, const Vector *point, const int ranges[NDIM] )
 {
   Info FunctionInfo(__func__);
@@ -249 +216 @@
       for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
         DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
-        for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
-          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) ) );
-                }
-          }
+        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);
@@ -273 +239 @@
 
 /** Calculates the distance (pair) correlation between a given element and a surface.
- * \param *molecules list of molecules structure
- * \param &elements vector of elements to correlate to surface
+ * \param *out output stream for debugging
+ * \param *molecules list of molecules structure
+ * \param *type element or NULL (if any element)
  * \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 std::vector<element *> &elements, const Tesselation * const Surface, const LinkedCell *LC )
+CorrelationToSurfaceMap *CorrelationToSurface(MoleculeListClass * const &molecules, const element * const type, const Tesselation * const Surface, const LinkedCell *LC )
 {
   Info FunctionInfo(__func__);
@@ -301 +268 @@
       for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
         DoLog(1) && (Log() << Verbose(1) << "\tCurrent atom is " << *(*iter) << "." << endl);
-        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) ) );
-          }
+        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 {
@@ -321 +287 @@
  * 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 *molecules list of molecules structure
- * \param &elements vector of elements to correlate to surface
+ * \param *out output stream for debugging
+ * \param *molecules list of molecules structure
+ * \param *type element or NULL (if any element)
  * \param *Surface pointer to Tesselation class surface
  * \param *LC LinkedCell structure to quickly find neighbouring atoms
@@ -328 +295 @@
  * \return Map of doubles with values as pairs of atom and the BoundaryTriangleSet that's closest
  */
-CorrelationToSurfaceMap *PeriodicCorrelationToSurface(MoleculeListClass * const &molecules, const std::vector<element *> &elements, const Tesselation * const Surface, const LinkedCell *LC, const int ranges[NDIM] )
+CorrelationToSurfaceMap *PeriodicCorrelationToSurface(MoleculeListClass * const &molecules, const element * const type, const Tesselation * const Surface, const LinkedCell *LC, const int ranges[NDIM] )
 {
   Info FunctionInfo(__func__);
@@ -355 +322 @@
       for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
         DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
-        for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
-          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;
-                  }
+        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;
                 }
-            // 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;
-          }
+              }
+          // 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);