Changeset c1a9d6 for src/Analysis

Timestamp:

Feb 22, 2012, 11:29:05 AM (13 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

b49c5d

Parents:

8fc1a6

git-author:

Frederik Heber <heber@…> (01/13/12 12:39:04)

git-committer:

Frederik Heber <heber@…> (02/22/12 11:29:05)

Message:

Rewrote PairCorrelation() to make use of new LinkedCell construct.

• removed PeriodicPairCorrelation() as this case is now contained.
• TESTFIX: Modified regression test Analysis/PairCorrelation-RangeTest because with new LinkedCell the given range can now actually be heeded.

Location:

src/Analysis

Files:

• analysis_correlation.cpp (modified) (4 diffs)
• analysis_correlation.hpp (modified) (4 diffs)
• unittests/AnalysisPairCorrelationUnitTest.cpp (modified) (3 diffs)

src/Analysis/analysis_correlation.cpp

@@ -20 +20 @@
 #include "CodePatterns/MemDebug.hpp"
 
+#include <algorithm>
 #include <iostream>
 #include <iomanip>
@@ -38 +39 @@
 #include "LinearAlgebra/Vector.hpp"
 #include "LinearAlgebra/RealSpaceMatrix.hpp"
+#include "LinkedCell/LinkedCell_View.hpp"
 #include "molecule.hpp"
 #include "Tesselation/tesselation.hpp"
@@ -275 +277 @@
 };
 
-
-/** Calculates the pair correlation between given elements.
- * Note given element order is unimportant (i.e. g(Si, O) === g(O, Si))
- * \param *molecules vector of molecules
- * \param &elements vector of elements to correlate
+/** Calculates the pair correlation between given atom sets.
+ *
+ * Note we correlate each of the \a &atomsfirst with each of the second set
+ * \a &atoms_second. However, we are aware of double counting. If an atom is
+ * in either set, the pair is counted only once.
+ *
+ * \param &atoms_first vector of atoms
+ * \param &atoms_second vector of atoms
+ * \param max_distance maximum distance for the correlation
  * \return Map of doubles with values the pair of the two atoms.
  */
-PairCorrelationMap *PairCorrelation(std::vector<molecule *> &molecules, const std::vector<const element *> &elements)
+PairCorrelationMap *PairCorrelation(
+    const World::AtomComposite &atoms_first,
+    const World::AtomComposite &atoms_second,
+    const double max_distance)
 {
   Info FunctionInfo(__func__);
@@ -289 +298 @@
   Box &domain = World::getInstance().getDomain();
 
-  if (molecules.empty()) {
-    ELOG(1, "No molecule given.");
+  if (atoms_first.empty() || atoms_second.empty()) {
+    ELOG(1, "No atoms given.");
     return outmap;
   }
-  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++)
-    (*MolWalker)->doCountAtoms();
-
-  // create all possible pairs of elements
-  set <pair<const element *,const element *> > PairsOfElements;
-  if (elements.size() >= 2) {
-    for (vector<const element *>::const_iterator type1 = elements.begin(); type1 != elements.end(); ++type1)
-      for (vector<const element *>::const_iterator type2 = elements.begin(); type2 != elements.end(); ++type2)
-        if (type1 != type2) {
-          PairsOfElements.insert( make_pair(*type1,*type2) );
-          LOG(1, "Creating element pair " << *(*type1) << " and " << *(*type2) << ".");
-        }
-  } else if (elements.size() == 1) { // one to all are valid
-    const element *elemental = *elements.begin();
-    PairsOfElements.insert( pair<const element *,const element*>(elemental,0) );
-    PairsOfElements.insert( pair<const element *,const element*>(0,elemental) );
-  } else { // all elements valid
-    PairsOfElements.insert( pair<element *, element*>((element *)NULL, (element *)NULL) );
-  }
-
+
+  //!> typedef for an unsorted container, (output) compatible with STL algorithms
+  typedef std::vector<const TesselPoint *> LinkedVector;
+
+  // create intersection (to know when to check for double-counting)
+  LinkedVector intersected_atoms(atoms_second.size(), NULL);
+  LinkedVector::iterator intersected_atoms_end =
+      std::set_intersection(
+          atoms_first.begin(),atoms_first.end(),
+          atoms_second.begin(), atoms_second.end(),
+          intersected_atoms.begin());
+  const LinkedCell::LinkedList intersected_atoms_set(intersected_atoms.begin(), intersected_atoms.end());
+
+  // create map
   outmap = new PairCorrelationMap;
-  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++){
-    LOG(2, "Current molecule is " << (*MolWalker)->getName() << ".");
-    for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
-      LOG(3, "Current atom is " << **iter << ".");
-      for (std::vector<molecule *>::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules.end(); MolOtherWalker++){
-        LOG(2, "Current other molecule is " << (*MolOtherWalker)->getName() << ".");
-        for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
-          LOG(3, "Current otheratom is " << **runner << ".");
-          if ((*iter)->getId() < (*runner)->getId()){
-            for (set <pair<const element *, const element *> >::iterator PairRunner = PairsOfElements.begin(); PairRunner != PairsOfElements.end(); ++PairRunner)
-              if ((PairRunner->first == (**iter).getType()) && (PairRunner->second == (**runner).getType())) {
-                distance = domain.periodicDistance((*iter)->getPosition(),(*runner)->getPosition());
-                //LOG(1, "Inserting " << *(*iter) << " and " << *(*runner));
-                outmap->insert (
-                    std::pair<double, std::pair <const atom *, const atom*> > (
-                        distance,
-                        std::pair<const atom *, const atom*> ((*iter), (*runner))
-                        )
-                    );
-              }
-          }
-        }
+
+  // get linked cell view
+  LinkedCell::LinkedCell_View LC = World::getInstance().getLinkedCell(max_distance);
+
+  // convert second to _sorted_ set
+  LinkedCell::LinkedList atoms_second_set(atoms_second.begin(), atoms_second.end());
+  LOG(2, "INFO: first set has " << atoms_first.size()
+      << " and second set has " << atoms_second_set.size() << " atoms.");
+
+  // fill map
+  for (World::AtomComposite::const_iterator iter = atoms_first.begin();
+      iter != atoms_first.end();
+      ++iter) {
+    const TesselPoint * const Walker = *iter;
+    LOG(3, "INFO: Current point is " << Walker->getName() << ".");
+    // obtain all possible neighbors (that is a sorted set)
+    LinkedCell::LinkedList ListOfNeighbors = LC.getPointsInsideSphere(
+        max_distance,
+        Walker->getPosition());
+    LOG(2, "INFO: There are " << ListOfNeighbors.size() << " neighbors.");
+
+    // create intersection with second set
+    // NOTE: STL algorithms do mostly not work on sorted container because reassignment
+    // of a value may also require changing its position.
+    LinkedVector intersected_set(atoms_second.size(), NULL);
+    LinkedVector::iterator intersected_end =
+        std::set_intersection(
+            ListOfNeighbors.begin(),ListOfNeighbors.end(),
+            atoms_second_set.begin(), atoms_second_set.end(),
+            intersected_set.begin());
+    // count remaining elements
+    LOG(2, "INFO: Intersection with second set has " << int(intersected_end - intersected_set.begin()) << " elements.");
+    // we have some possible candidates, go through each
+    for (LinkedVector::const_iterator neighboriter = intersected_set.begin();
+        neighboriter != intersected_end;
+        ++neighboriter) {
+      const TesselPoint * const OtherWalker = (*neighboriter);
+      LinkedCell::LinkedList::const_iterator equaliter = intersected_atoms_set.find(OtherWalker);
+      if ((equaliter !=  intersected_atoms_set.end()) && (OtherWalker <= Walker)) {
+        // present in both sets, assure that we are larger
+        continue;
       }
+      LOG(3, "INFO: Current other point is " << *OtherWalker << ".");
+      const double distance = domain.periodicDistance(OtherWalker->getPosition(),Walker->getPosition());
+      LOG(3, "INFO: Resulting distance is " << distance << ".");
+      outmap->insert (
+          std::pair<double, std::pair <const TesselPoint *, const TesselPoint*> > (
+              distance,
+              std::make_pair (Walker, OtherWalker)
+              )
+          );
     }
   }
-  return outmap;
-};
-
-/** 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 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(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const int ranges[NDIM] )
-{
-  Info FunctionInfo(__func__);
-  PairCorrelationMap *outmap = new PairCorrelationMap;
-  double distance = 0.;
-  int n[NDIM];
-  Vector checkX;
-  Vector periodicX;
-  int Othern[NDIM];
-  Vector checkOtherX;
-  Vector periodicOtherX;
-
-  if (molecules.empty()) {
-    ELOG(1, "No molecule given.");
-    return outmap;
-  }
-  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++)
-    (*MolWalker)->doCountAtoms();
-
-  // create all possible pairs of elements
-  set <pair<const element *,const element *> > PairsOfElements;
-  if (elements.size() >= 2) {
-    for (vector<const element *>::const_iterator type1 = elements.begin(); type1 != elements.end(); ++type1)
-      for (vector<const element *>::const_iterator type2 = elements.begin(); type2 != elements.end(); ++type2)
-        if (type1 != type2) {
-          PairsOfElements.insert( make_pair(*type1,*type2) );
-          LOG(1, "Creating element pair " << *(*type1) << " and " << *(*type2) << ".");
-        }
-  } else if (elements.size() == 1) { // one to all are valid
-    const element *elemental = *elements.begin();
-    PairsOfElements.insert( pair<const element *,const element*>(elemental,0) );
-    PairsOfElements.insert( pair<const element *,const element*>(0,elemental) );
-  } else { // all elements valid
-    PairsOfElements.insert( pair<element *, element*>((element *)NULL, (element *)NULL) );
-  }
-
-  outmap = new PairCorrelationMap;
-  for (std::vector<molecule *>::const_iterator MolWalker = molecules.begin(); MolWalker != molecules.end(); MolWalker++){
-    RealSpaceMatrix FullMatrix = World::getInstance().getDomain().getM();
-    RealSpaceMatrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
-    LOG(2, "Current molecule is " << *MolWalker << ".");
-    for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
-      LOG(3, "Current atom is " << **iter << ".");
-      periodicX = FullInverseMatrix * ((*iter)->getPosition()); // 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 (std::vector<molecule *>::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules.end(); MolOtherWalker++){
-                LOG(2, "Current other molecule is " << *MolOtherWalker << ".");
-                for (molecule::const_iterator runner = (*MolOtherWalker)->begin(); runner != (*MolOtherWalker)->end(); ++runner) {
-                  LOG(3, "Current otheratom is " << **runner << ".");
-                  if ((*iter)->getId() < (*runner)->getId()){
-                    for (set <pair<const element *,const element *> >::iterator PairRunner = PairsOfElements.begin(); PairRunner != PairsOfElements.end(); ++PairRunner)
-                      if ((PairRunner->first == (**iter).getType()) && (PairRunner->second == (**runner).getType())) {
-                        periodicOtherX = FullInverseMatrix * ((*runner)->getPosition()); // 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 = FullMatrix * (Vector(Othern[0], Othern[1], Othern[2]) + periodicOtherX);
-                              distance = checkX.distance(checkOtherX);
-                              //LOG(1, "Inserting " << *(*iter) << " and " << *(*runner));
-                              outmap->insert (
-                                  std::pair<double, std::pair <const atom *, const atom*> > (
-                                      distance,
-                                      std::pair<const atom *, const atom*> (
-                                          (*iter),
-                                          (*runner))
-                                      )
-                                  );
-                            }
-                      }
-                    }
-                  }
-                }
-      }
-    }
-  }
-
+  // and return
   return outmap;
 };

src/Analysis/analysis_correlation.hpp

@@ -34 +34 @@
 #include "CodePatterns/Verbose.hpp"
 #include "Helpers/helpers.hpp"
+#include "World.hpp"
 
 /****************************************** forward declarations *****************************/
@@ -46 +47 @@
 /********************************************** definitions *********************************/
 
-typedef multimap<double, pair<const atom *, const atom *> > PairCorrelationMap;
+typedef multimap<double, pair<const TesselPoint *, const TesselPoint *> > PairCorrelationMap;
 typedef multimap<double, const atom * > DipoleAngularCorrelationMap;
 typedef multimap<double, pair<const molecule *, const molecule *> > DipoleCorrelationMap;
@@ -65 +66 @@
 DipoleAngularCorrelationMap *DipoleAngularCorrelation(const Formula &DipoleFormula, const size_t timestep, const std::map<atomId_t, Vector> &ZeroVector, const enum ResetWorldTime DoTimeReset = DontResetTime);
 DipoleCorrelationMap *DipoleCorrelation(std::vector<molecule *> &molecules);
-PairCorrelationMap *PairCorrelation(std::vector<molecule *> &molecules, const std::vector<const element *> &elements);
+PairCorrelationMap *PairCorrelation(
+    const World::AtomComposite &atoms_first,
+    const World::AtomComposite &atoms_second,
+    const double max_distance);
 CorrelationToPointMap *CorrelationToPoint(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const Vector *point );
 CorrelationToSurfaceMap *CorrelationToSurface(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const Tesselation * const Surface, const LinkedCell_deprecated *LC );
-PairCorrelationMap *PeriodicPairCorrelation(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const int ranges[NDIM] );
 CorrelationToPointMap *PeriodicCorrelationToPoint(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const Vector *point, const int ranges[NDIM] );
 CorrelationToSurfaceMap *PeriodicCorrelationToSurface(std::vector<molecule *> &molecules, const std::vector<const element *> &elements, const Tesselation * const Surface, const LinkedCell_deprecated *LC, const int ranges[NDIM] );
@@ -158 +161 @@
     if (!BinPairMapInserter.second) {  // bin already present, increase
       BinPairMapInserter.first->second += 1;
+    } else { // bin not present, then remove again
+      outmap->erase(BinPairMapInserter.first);
     }
   }

src/Analysis/unittests/AnalysisPairCorrelationUnitTest.cpp

@@ -30 +30 @@
 
 #include "Atom/atom.hpp"
+#include "Descriptors/AtomTypeDescriptor.hpp"
 #include "Element/element.hpp"
 #include "Element/periodentafel.hpp"
@@ -61 +62 @@
   hydrogen = World::getInstance().getPeriode()->FindElement(1);
   CPPUNIT_ASSERT(hydrogen != NULL && "hydrogen element not found");
-  elements.push_back(hydrogen);
-  elements.push_back(hydrogen);
 
   // construct molecule (tetraeder of hydrogens)
@@ -87 +86 @@
 
   // init maps
-  World::getInstance().selectAllMolecules(AllMolecules());
-  allMolecules = World::getInstance().getSelectedMolecules();
-  CPPUNIT_ASSERT_EQUAL( (size_t) 1, allMolecules.size());
-  correlationmap = PairCorrelation( allMolecules, elements);
+  World::AtomComposite atoms_first = World::getInstance().getAllAtoms(AtomByType(hydrogen));
+  CPPUNIT_ASSERT_EQUAL( (size_t) 4, atoms_first.size());
+  correlationmap = PairCorrelation( atoms_first, atoms_first, 5.);
   binmap = NULL;