Changeset 7ea9e6

Timestamp:

Nov 4, 2009, 6:15:10 PM (15 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:

7326b2

Parents:

99593f

git-author:

Frederik Heber <heber@…> (11/04/09 18:11:09)

git-committer:

Frederik Heber <heber@…> (11/04/09 18:15:10)

Message:

Periodic variants of AnalysisPair...() implemented.

• BUGFIX: Vector::DistanceToPlane() - we did not check whether sign is 0.
• NOTE: Due to changes in Vector::DistanceToPlane() with sign, Unit test AnalysisCorrelationToSurfaceUnitTest had to be changed:
  • we now have bin -0.5 filled with 1 (instead of 0.) and -0.288 instead of 0.288
  • find() replaced by lower_bound
• new functions: PeriodicPairCorrelation(), PeriodicCorrelationToPoint(), PeriodicCorrelationToSurface()
  • each has a ranges[NDIM] argument with specifies the neighbours to scan: [ -ranges[i], ranges[i] ]
  • the atom::node is periodically translated to each periodic cell and the distance calculated.
• NOTE: make check was broken before due to implementation of Periodic variants as normal ones (yielding more points that expected in unit test).

Signed-off-by: Frederik Heber <heber@…>

Location:

src

Files:

• analysis_correlation.cpp (modified) (3 diffs)
• analysis_correlation.hpp (modified) (1 diff)
• builder.cpp (modified) (1 diff)
• unittests/AnalysisCorrelationToSurfaceUnitTest.cpp (modified) (2 diffs)
• vector.cpp (modified) (1 diff)

src/analysis_correlation.cpp

@@ -65 +65 @@
 };
 
+/** 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 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( ofstream * const out, MoleculeListClass * const &molecules, const element * const type1, const element * const type2, const int ranges[NDIM] )
+{
+  PairCorrelationMap *outmap = NULL;
+  double distance = 0.;
+  int n[NDIM];
+  Vector checkX;
+  Vector periodicX;
+  int Othern[NDIM];
+  Vector checkOtherX;
+  Vector periodicOtherX;
+
+  if (molecules->ListOfMolecules.empty()) {
+    cerr << Verbose(1) <<"No molecule given." << endl;
+    return outmap;
+  }
+  outmap = new PairCorrelationMap;
+  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
+    if ((*MolWalker)->ActiveFlag) {
+      const double * const FullMatrix = ReturnFullMatrixforSymmetric((*MolWalker)->cell_size);
+      const double * const FullInverseMatrix = InverseMatrix(FullMatrix);
+      cerr << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
+      atom *Walker = (*MolWalker)->start;
+      while (Walker->next != (*MolWalker)->end) {
+        Walker = Walker->next;
+        *out << Verbose(3) << "Current atom is " << *Walker << "." << endl;
+        if ((type1 == NULL) || (Walker->type == type1)) {
+          periodicX.CopyVector(Walker->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.Init(n[0], n[1], n[2]);
+                checkX.AddVector(&periodicX);
+                checkX.MatrixMultiplication(FullMatrix);
+                for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++)
+                  if ((*MolOtherWalker)->ActiveFlag) {
+                    *out << Verbose(2) << "Current other molecule is " << *MolOtherWalker << "." << endl;
+                    atom *OtherWalker = (*MolOtherWalker)->start;
+                    while (OtherWalker->next != (*MolOtherWalker)->end) { // only go up to Walker
+                      OtherWalker = OtherWalker->next;
+                      *out << Verbose(3) << "Current otheratom is " << *OtherWalker << "." << endl;
+                      if (Walker->nr < OtherWalker->nr)
+                        if ((type2 == NULL) || (OtherWalker->type == type2)) {
+                          periodicOtherX.CopyVector(OtherWalker->node);
+                          periodicOtherX.MatrixMultiplication(FullInverseMatrix);  // 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.Init(Othern[0], Othern[1], Othern[2]);
+                                checkOtherX.AddVector(&periodicOtherX);
+                                checkOtherX.MatrixMultiplication(FullMatrix);
+                                distance = checkX.Distance(&checkOtherX);
+                                //*out << Verbose(1) <<"Inserting " << *Walker << " and " << *OtherWalker << endl;
+                                outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> (Walker, OtherWalker) ) );
+                              }
+                        }
+                  }
+              }
+            }
+        }
+      }
+    }
+
+  return outmap;
+};
+
 /** Calculates the distance (pair) correlation between a given element and a point.
  * \param *out output stream for debugging
@@ -100 +177 @@
 };
 
+/** 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 *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(  ofstream * const out, MoleculeListClass * const &molecules, const element * const type, const Vector *point, const int ranges[NDIM] )
+{
+  CorrelationToPointMap *outmap = NULL;
+  double distance = 0.;
+  int n[NDIM];
+  Vector periodicX;
+  Vector checkX;
+
+  if (molecules->ListOfMolecules.empty()) {
+    *out << Verbose(1) <<"No molecule given." << endl;
+    return outmap;
+  }
+  outmap = new CorrelationToPointMap;
+  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
+    if ((*MolWalker)->ActiveFlag) {
+      const double * const FullMatrix = ReturnFullMatrixforSymmetric((*MolWalker)->cell_size);
+      const double * const FullInverseMatrix = InverseMatrix(FullMatrix);
+      *out << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
+      atom *Walker = (*MolWalker)->start;
+      while (Walker->next != (*MolWalker)->end) {
+        Walker = Walker->next;
+        *out << Verbose(3) << "Current atom is " << *Walker << "." << endl;
+        if ((type == NULL) || (Walker->type == type)) {
+          periodicX.CopyVector(Walker->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.Init(n[0], n[1], n[2]);
+                checkX.AddVector(&periodicX);
+                checkX.MatrixMultiplication(FullMatrix);
+                distance = checkX.Distance(point);
+                *out << Verbose(4) << "Current distance is " << distance << "." << endl;
+                outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (Walker, point) ) );
+              }
+        }
+      }
+    }
+
+  return outmap;
+};
+
 /** Calculates the distance (pair) correlation between a given element and a surface.
  * \param *out output stream for debugging
@@ -114 +242 @@
   class BoundaryTriangleSet *triangle = NULL;
   Vector centroid;
-  int ranges[NDIM] = {1, 1, 1};
+
+  if ((Surface == NULL) || (LC == NULL) || (molecules->ListOfMolecules.empty())) {
+    *out << Verbose(1) <<"No Tesselation, no LinkedCell or no molecule given." << endl;
+    return outmap;
+  }
+  outmap = new CorrelationToSurfaceMap;
+  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
+    if ((*MolWalker)->ActiveFlag) {
+      *out << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
+      atom *Walker = (*MolWalker)->start;
+      while (Walker->next != (*MolWalker)->end) {
+        Walker = Walker->next;
+        *out << Verbose(3) << "Current atom is " << *Walker << "." << endl;
+        if ((type == NULL) || (Walker->type == type)) {
+          triangle = Surface->FindClosestTriangleToPoint(out, Walker->node, LC );
+          if (triangle != NULL) {
+            distance = DistanceToTrianglePlane(out, Walker->node, triangle);
+            outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> (Walker, triangle) ) );
+          }
+        }
+      }
+    }
+
+  return outmap;
+};
+
+/** Calculates the distance (pair) correlation between a given element, all its periodic images and and a surface.
+ * Note that we also put all periodic images found in the cells given by [ -ranges[i], ranges[i] ] and i=0,...,NDIM-1.
+ * I.e. We multiply the atom::node with the inverse of the domain matrix, i.e. transform it to \f$[0,0^3\f$, then add per
+ * 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 *Surface pointer to Tesselation class surface
+ * \param *LC LinkedCell structure to quickly find neighbouring atoms
+ * \param ranges[NDIM] interval boundaries for the periodic images to scan also
+ * \return Map of doubles with values as pairs of atom and the BoundaryTriangleSet that's closest
+ */
+CorrelationToSurfaceMap *PeriodicCorrelationToSurface( ofstream * const out, MoleculeListClass * const &molecules, const element * const type, const Tesselation * const Surface, const LinkedCell *LC, const int ranges[NDIM] )
+{
+  CorrelationToSurfaceMap *outmap = NULL;
+  double distance = 0;
+  class BoundaryTriangleSet *triangle = NULL;
+  Vector centroid;
   int n[NDIM];
   Vector periodicX;

src/analysis_correlation.hpp

@@ -48 +48 @@
 CorrelationToPointMap *CorrelationToPoint( ofstream * const out, MoleculeListClass * const &molecules, const element * const type, const Vector *point );
 CorrelationToSurfaceMap *CorrelationToSurface( ofstream * const out, MoleculeListClass * const &molecules, const element * const type, const Tesselation * const Surface, const LinkedCell *LC );
+PairCorrelationMap *PeriodicPairCorrelation( ofstream * const out, MoleculeListClass * const &molecules, const element * const type1, const element * const  type2, const int ranges[NDIM] );
+CorrelationToPointMap *PeriodicCorrelationToPoint( ofstream * const out, MoleculeListClass * const &molecules, const element * const type, const Vector *point, const int ranges[NDIM] );
+CorrelationToSurfaceMap *PeriodicCorrelationToSurface( ofstream * const out, MoleculeListClass * const &molecules, const element * const type, const Tesselation * const Surface, const LinkedCell *LC, const int ranges[NDIM] );
 double GetBin ( const double value, const double BinWidth, const double BinStart );
 void OutputCorrelation( ofstream * const file, const BinPairMap * const map );

src/builder.cpp

@@ -1632 +1632 @@
                 element *elemental = periode->FindElement((const int) atoi(argv[argptr]));
                 FindNonConvexBorder((ofstream *)&cout, Boundary, TesselStruct, LCList, radius, NULL);
-                CorrelationToSurfaceMap *surfacemap = CorrelationToSurface( (ofstream *)&cout, molecules, elemental, TesselStruct, LCList );
+                int ranges[NDIM] = {1,1,1};
+                CorrelationToSurfaceMap *surfacemap = PeriodicCorrelationToSurface( (ofstream *)&cout, molecules, elemental, TesselStruct, LCList, ranges );
                 BinPairMap *binmap = BinData( (ofstream *)&cout, surfacemap, 0.5, 0., 0. );
                 OutputCorrelation ( &binoutput, binmap );

src/unittests/AnalysisCorrelationToSurfaceUnitTest.cpp

@@ -186 +186 @@
   CPPUNIT_ASSERT_EQUAL( 3, tester->second );
   // inner point
-  tester = binmap->lower_bound(0.);
+  tester = binmap->lower_bound(-0.5);
   CPPUNIT_ASSERT( tester != binmap->end() );
   CPPUNIT_ASSERT_EQUAL( 1, tester->second );
@@ -200 +200 @@
   OutputCorrelation ( (ofstream *)&cout, binmap );
   // three outside points
-  tester = binmap->find(3.);
+  tester = binmap->lower_bound(3.);
   CPPUNIT_ASSERT( tester != binmap->end() );
   CPPUNIT_ASSERT_EQUAL( 3, tester->second );
   // inner point
-  tester = binmap->find(0.);
+  tester = binmap->lower_bound(-0.5);
   CPPUNIT_ASSERT( tester != binmap->end() );
   CPPUNIT_ASSERT_EQUAL( 1, tester->second );

src/vector.cpp

@@ -276 +276 @@
   temp.SubtractVector(PlaneOffset);
   double sign = temp.ScalarProduct(PlaneNormal);
-  sign /= fabs(sign);
+  if (fabs(sign) > MYEPSILON)
+    sign /= fabs(sign);
+  else
+    sign = 0.;
 
   return (temp.Norm()*sign);