Changeset 99593f

Timestamp:

Nov 4, 2009, 5:34:05 PM (16 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:

7ea9e6

Parents:

c9bce3e

Message:

Extension to the periodic boundary case for analysis_correlation.cpp

• Distance() -> PeriodicDistance(): PairCorrelation(), CorrelationToPoint(), CorrelationToSurface()
• CorrelationToSurface() - goes through periodic cells given by a ranges variable and puts into map as well.
• const: ReturnFullMatrixforSymmetric() and InverseMatrix() now have const * const parameter
• InverseMatrix() shifted from class Vector into helpers.cpp
• adapted due to InverseMatrix(): molecule::CenterInBox(), molecule::BoundInBox(), molecule::TranslatePeriodically()
• new function: RDET3(), RDET2() (former #define's)
• FIX: molecule::CorrectBondDegree() would loop eternally, if bonds cannot be fully corrected.
• FIX for NULL return: Tesselation::FindClosestTrianglesToPoint(), Tesselation::GetCircleOfConnectedPoints()
• FIX: Vector::DistanceToPlane() - did not give correct sign before according to PlaneNormal orientation

other stuff:

• less output: BondedParticle::CorrectBondDegree(), molecule::CyclicStructureAnalysis()

Location:

src

Files:

• analysis_correlation.cpp (modified) (5 diffs)
• atom_bondedparticle.cpp (modified) (3 diffs)
• helpers.cpp (modified) (2 diffs)
• helpers.hpp (modified) (2 diffs)
• molecule_geometry.cpp (modified) (3 diffs)
• molecule_graph.cpp (modified) (3 diffs)
• tesselation.cpp (modified) (2 diffs)
• vector.cpp (modified) (2 diffs)
• vector.hpp (modified) (2 diffs)

src/analysis_correlation.cpp

@@ -52 +52 @@
                 if (Walker->nr < OtherWalker->nr)
                   if ((type2 == NULL) || (OtherWalker->type == type2)) {
-                    distance = Walker->node->Distance(OtherWalker->node);
+                    distance = Walker->node->PeriodicDistance(OtherWalker->node, (*MolWalker)->cell_size);
                     //*out << Verbose(1) <<"Inserting " << *Walker << " and " << *OtherWalker << endl;
                     outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> (Walker, OtherWalker) ) );
@@ -90 +90 @@
         *out << Verbose(3) << "Current atom is " << *Walker << "." << endl;
         if ((type == NULL) || (Walker->type == type)) {
-          distance = Walker->node->Distance(point);
+          distance = Walker->node->PeriodicDistance(point, (*MolWalker)->cell_size);
           *out << Verbose(4) << "Current distance is " << distance << "." << endl;
           outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> (Walker, point) ) );
@@ -111 +111 @@
 {
   CorrelationToSurfaceMap *outmap = NULL;
-  double distance = 0.;
+  double distance = 0;
   class BoundaryTriangleSet *triangle = NULL;
   Vector centroid;
+  int ranges[NDIM] = {1, 1, 1};
+  int n[NDIM];
+  Vector periodicX;
+  Vector checkX;
 
   if ((Surface == NULL) || (LC == NULL) || (molecules->ListOfMolecules.empty())) {
@@ -122 +126 @@
   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;
@@ -128 +134 @@
         *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) ) );
+          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);
+                triangle = Surface->FindClosestTriangleToPoint(out, &checkX, LC );
+                if (triangle != NULL) {
+                  distance = DistanceToTrianglePlane(out, &checkX, triangle);
+                  outmap->insert ( pair<double, pair<atom *, BoundaryTriangleSet*> >(distance, pair<atom *, BoundaryTriangleSet*> (Walker, triangle) ) );
+                }
           }
         }

src/atom_bondedparticle.cpp

@@ -126 +126 @@
   bond *CandidateBond = NULL;
 
-  *out << Verbose(3) << "Walker " << *this << ": " << (int)this->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
+  //*out << Verbose(3) << "Walker " << *this << ": " << (int)this->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
   NoBonds = CountBonds();
   if ((int)(type->NoValenceOrbitals) > NoBonds) { // we have a mismatch, check all bonding partners for mismatch
@@ -132 +132 @@
       OtherWalker = (*Runner)->GetOtherAtom(this);
       OtherNoBonds = OtherWalker->CountBonds();
-      *out << Verbose(3) << "OtherWalker " << *OtherWalker << ": " << (int)OtherWalker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
+      //*out << Verbose(3) << "OtherWalker " << *OtherWalker << ": " << (int)OtherWalker->type->NoValenceOrbitals << " > " << NoBonds << "?" << endl;
       if ((int)(OtherWalker->type->NoValenceOrbitals) > NoBonds) { // check if possible candidate
         if ((CandidateBond == NULL) || (ListOfBonds.size() > OtherWalker->ListOfBonds.size())) { // pick the one with fewer number of bonds first
           CandidateBond = (*Runner);
-          *out << Verbose(3) << "New candidate is " << *CandidateBond << "." << endl;
+          //*out << Verbose(3) << "New candidate is " << *CandidateBond << "." << endl;
         }
       }
@@ -144 +144 @@
       *out << Verbose(2) << "Increased bond degree for bond " << *CandidateBond << "." << endl;
     } else {
-      *out << Verbose(2) << "Could not find correct degree for atom " << *this << "." << endl;
+      //*out << Verbose(2) << "Could not find correct degree for atom " << *this << "." << endl;
       FalseBondDegree++;
     }

src/helpers.cpp

@@ -133 +133 @@
  * \return allocated NDIM*NDIM array with the symmetric matrix
  */
-double * ReturnFullMatrixforSymmetric(double *symm)
+double * ReturnFullMatrixforSymmetric(const double * const symm)
 {
   double *matrix = Malloc<double>(NDIM * NDIM, "molecule::ReturnFullMatrixforSymmetric: *matrix");
@@ -148 +148 @@
 };
 
+/** Calculate the inverse of a 3x3 matrix.
+ * \param *matrix NDIM_NDIM array
+ */
+double * InverseMatrix( const double * const A)
+{
+  double *B = Malloc<double>(NDIM * NDIM, "Vector::InverseMatrix: *B");
+  double detA = RDET3(A);
+  double detAReci;
+
+  for (int i=0;i<NDIM*NDIM;++i)
+    B[i] = 0.;
+  // calculate the inverse B
+  if (fabs(detA) > MYEPSILON) {;  // RDET3(A) yields precisely zero if A irregular
+    detAReci = 1./detA;
+    B[0] =  detAReci*RDET2(A[4],A[5],A[7],A[8]);    // A_11
+    B[1] = -detAReci*RDET2(A[1],A[2],A[7],A[8]);    // A_12
+    B[2] =  detAReci*RDET2(A[1],A[2],A[4],A[5]);    // A_13
+    B[3] = -detAReci*RDET2(A[3],A[5],A[6],A[8]);    // A_21
+    B[4] =  detAReci*RDET2(A[0],A[2],A[6],A[8]);    // A_22
+    B[5] = -detAReci*RDET2(A[0],A[2],A[3],A[5]);    // A_23
+    B[6] =  detAReci*RDET2(A[3],A[4],A[6],A[7]);    // A_31
+    B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]);    // A_32
+    B[8] =  detAReci*RDET2(A[0],A[1],A[3],A[4]);    // A_33
+  }
+  return B;
+};
+
+/** hard-coded determinant of a 3x3 matrix.
+ * \param a[9] matrix
+ * \return \f$det(a)\f$
+ */
+double RDET3(const double a[NDIM*NDIM])
+{
+  return ((a)[0]*(a)[4]*(a)[8] + (a)[3]*(a)[7]*(a)[2] + (a)[6]*(a)[1]*(a)[5] - (a)[2]*(a)[4]*(a)[6] - (a)[5]*(a)[7]*(a)[0] - (a)[8]*(a)[1]*(a)[3]);
+};
+
+/** hard-coded determinant of a 2x2 matrix.
+ * \param a[4] matrix
+ * \return \f$det(a)\f$
+ */
+double RDET2(const double a[4])
+{
+  return ((a[0])*(a[3])-(a[1])*(a[2]));
+};
+
+/** hard-coded determinant of a 2x2 matrix.
+ * \param a0 (0,0) entry of matrix
+ * \param a1 (0,1) entry of matrix
+ * \param a2 (1,0) entry of matrix
+ * \param a3 (1,1) entry of matrix
+ * \return \f$det(a)\f$
+ */
+double RDET2(const double a0, const double a1, const double a2, const double a3)
+{
+  return ((a0)*(a3)-(a1)*(a2));
+};
+
 /** Comparison function for GSL heapsort on distances in two molecules.
  * \param *a

src/helpers.hpp

@@ -18 +18 @@
 #include <fstream>
 
+#include "defs.hpp"
 #include "memoryallocator.hpp"
+
+/********************************************** definitions *********************************/
+
+// some algebraic matrix stuff
+double RDET3(const double a[NDIM*NDIM]);
+double RDET2(const double a[4]);
+double RDET2(const double a0, const double a1, const double a2, const double a3);
 
 /********************************************** helpful functions *********************************/
@@ -49 +57 @@
 bool IsValidNumber( const char *string);
 int CompareDoubles (const void * a, const void * b);
-double * ReturnFullMatrixforSymmetric(double *cell_size);
+double * ReturnFullMatrixforSymmetric(const double * const cell_size);
+double * InverseMatrix(const double * const A);
 void performCriticalExit();
 

src/molecule_geometry.cpp

@@ -24 +24 @@
 {
   bool status = true;
-  Vector x;
   const Vector *Center = DetermineCenterOfAll(out);
   double *M = ReturnFullMatrixforSymmetric(cell_size);
-  double *Minv = x.InverseMatrix(M);
+  double *Minv = InverseMatrix(M);
 
   // go through all atoms
@@ -46 +45 @@
 {
   bool status = true;
-  Vector x;
   double *M = ReturnFullMatrixforSymmetric(cell_size);
-  double *Minv = x.InverseMatrix(M);
+  double *Minv = InverseMatrix(M);
 
   // go through all atoms
@@ -231 +229 @@
 void molecule::TranslatePeriodically(const Vector *trans)
 {
-  Vector x;
   double *M = ReturnFullMatrixforSymmetric(cell_size);
-  double *Minv = x.InverseMatrix(M);
+  double *Minv = InverseMatrix(M);
 
   // go through all atoms

src/molecule_graph.cpp

@@ -221 +221 @@
 int molecule::CorrectBondDegree(ofstream *out) const
 {
-  int No = 0;
+  int No = 0, OldNo = -1;
 
   if (BondCount != 0) {
     *out << Verbose(1) << "Correcting Bond degree of each bond ... " << endl;
     do {
+      OldNo = No;
       No = SumPerAtom(&atom::CorrectBondDegree, out);
-    } while (No);
+    } while (OldNo != No);
     *out << " done." << endl;
   } else {
@@ -418 +419 @@
   if (!DFS.BackStepping) { // coming from (8) want to go to (3)
     // (9) remove all from stack till Walker (including), these and Root form a component
-    DFS.AtomStack->Output(out);
+    //DFS.AtomStack->Output(out);
     mol->SetNextComponentNumber(DFS.Root, DFS.ComponentNumber);
     *out << Verbose(3) << "(9) Root[" << DFS.Root->Name << "]'s Component is " << DFS.ComponentNumber << "." << endl;
@@ -881 +882 @@
   InitializeBFSAccounting(out, BFS, AtomCount);
 
-  *out << Verbose(1) << "Back edge list - ";
-  BackEdgeStack->Output(out);
+  //*out << Verbose(1) << "Back edge list - ";
+  //BackEdgeStack->Output(out);
 
   *out << Verbose(1) << "Analysing cycles ... " << endl;

src/tesselation.cpp

@@ -2538 +2538 @@
   } else {
     list<TesselPoint*> *connectedClosestPoints = GetCircleOfConnectedPoints(out, trianglePoints[0], x);
-    trianglePoints[1] = connectedClosestPoints->front();
-    trianglePoints[2] = connectedClosestPoints->back();
-    for (int i=0;i<3;i++) {
-      if (trianglePoints[i] == NULL) {
-        *out << Verbose(1) << "ERROR: IsInnerPoint encounters serious error, point " << i << " not found." << endl;
+    if (connectedClosestPoints != NULL) {
+      trianglePoints[1] = connectedClosestPoints->front();
+      trianglePoints[2] = connectedClosestPoints->back();
+      for (int i=0;i<3;i++) {
+        if (trianglePoints[i] == NULL) {
+          *out << Verbose(1) << "ERROR: IsInnerPoint encounters serious error, point " << i << " not found." << endl;
+        }
+        //*out << Verbose(2) << "List of triangle points:" << endl;
+        //*out << Verbose(3) << *trianglePoints[i] << endl;
       }
-      //*out << Verbose(2) << "List of triangle points:" << endl;
-      //*out << Verbose(3) << *trianglePoints[i] << endl;
-    }
-
-    triangles = FindTriangles(trianglePoints);
-    *out << Verbose(2) << "List of possible triangles:" << endl;
-    for(list<BoundaryTriangleSet*>::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++)
-      *out << Verbose(3) << **Runner << endl;
-
-    delete(connectedClosestPoints);
+
+      triangles = FindTriangles(trianglePoints);
+      *out << Verbose(2) << "List of possible triangles:" << endl;
+      for(list<BoundaryTriangleSet*>::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++)
+        *out << Verbose(3) << **Runner << endl;
+
+      delete(connectedClosestPoints);
+    } else {
+      triangles = NULL;
+      *out << Verbose(1) << "There is no circle of connected points!" << endl;
+    }
   }
   
-  if (triangles->empty()) {
+  if ((triangles == NULL) || (triangles->empty())) {
     *out << Verbose(0) << "ERROR: There is no nearest triangle. Please check the tesselation structure.";
     delete(triangles);
@@ -2724 +2729 @@
   Vector helper;
 
+  if (connectedPoints == NULL) {
+    *out << Verbose(2) << "Could not find any connected points!" << endl;
+    delete(connectedCircle);
+    return NULL;
+  }
   *out << Verbose(2) << "Begin of GetCircleOfConnectedPoints" << endl;
 

src/vector.cpp

@@ -275 +275 @@
   temp.AddVector(this);
   temp.SubtractVector(PlaneOffset);
-
-  return temp.Norm();
+  double sign = temp.ScalarProduct(PlaneNormal);
+  sign /= fabs(sign);
+
+  return (temp.Norm()*sign);
 };
 
@@ -713 +715 @@
   for (int i=NDIM;i--;)
     x[i] = C.x[i];
-};
-
-/** Calculate the inverse of a 3x3 matrix.
- * \param *matrix NDIM_NDIM array
- */
-double * Vector::InverseMatrix( const double * const A)
-{
-  double *B = Malloc<double>(NDIM * NDIM, "Vector::InverseMatrix: *B");
-  double detA = RDET3(A);
-  double detAReci;
-
-  for (int i=0;i<NDIM*NDIM;++i)
-    B[i] = 0.;
-  // calculate the inverse B
-  if (fabs(detA) > MYEPSILON) {;  // RDET3(A) yields precisely zero if A irregular
-    detAReci = 1./detA;
-    B[0] =  detAReci*RDET2(A[4],A[5],A[7],A[8]);    // A_11
-    B[1] = -detAReci*RDET2(A[1],A[2],A[7],A[8]);    // A_12
-    B[2] =  detAReci*RDET2(A[1],A[2],A[4],A[5]);    // A_13
-    B[3] = -detAReci*RDET2(A[3],A[5],A[6],A[8]);    // A_21
-    B[4] =  detAReci*RDET2(A[0],A[2],A[6],A[8]);    // A_22
-    B[5] = -detAReci*RDET2(A[0],A[2],A[3],A[5]);    // A_23
-    B[6] =  detAReci*RDET2(A[3],A[4],A[6],A[7]);    // A_31
-    B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]);    // A_32
-    B[8] =  detAReci*RDET2(A[0],A[1],A[3],A[4]);    // A_33
-  }
-  return B;
 };
 

src/vector.hpp

@@ -62 +62 @@
   void Scale(const double factor);
   void MatrixMultiplication(const double * const M);
-  double * InverseMatrix(const double * const A);
   void InverseMatrixMultiplication(const double * const M);
   void KeepPeriodic(ofstream *out, const double * const matrix);
@@ -91 +90 @@
 Vector& operator-(const Vector& a, const Vector& b);
 
-// some algebraic matrix stuff
-#define RDET3(a) ((a)[0]*(a)[4]*(a)[8] + (a)[3]*(a)[7]*(a)[2] + (a)[6]*(a)[1]*(a)[5] - (a)[2]*(a)[4]*(a)[6] - (a)[5]*(a)[7]*(a)[0] - (a)[8]*(a)[1]*(a)[3])  //!< hard-coded determinant of a 3x3 matrix
-#define RDET2(a0,a1,a2,a3) ((a0)*(a3)-(a1)*(a2))                      //!< hard-coded determinant of a 2x2 matrix
-
-
 
 #endif /*VECTOR_HPP_*/