Changes in src/vector.cpp [ce5ac3:042f82]

File:

• src/vector.cpp (modified) (8 diffs)

src/vector.cpp

@@ -49 +49 @@
  * \param *y array to second vector
  * \param *cell_size 6-dimensional array with (xx, xy, yy, xz, yz, zz) entries specifying the periodic cell
- * \return \f$| x - y |^2\f$
+ * \return \f$| x - y |\f$
  */
 double Vector::PeriodicDistance(const Vector *y, const double *cell_size) const
@@ -84 +84 @@
 };
 
+/** Calculates distance between this and another vector in a periodic cell.
+ * \param *y array to second vector
+ * \param *cell_size 6-dimensional array with (xx, xy, yy, xz, yz, zz) entries specifying the periodic cell
+ * \return \f$| x - y |^2\f$
+ */
+double Vector::PeriodicDistanceSquared(const Vector *y, const double *cell_size) const
+{
+  double res = DistanceSquared(y), tmp, matrix[NDIM*NDIM];
+  Vector Shiftedy, TranslationVector;
+  int N[NDIM];
+  matrix[0] = cell_size[0];
+  matrix[1] = cell_size[1];
+  matrix[2] = cell_size[3];
+  matrix[3] = cell_size[1];
+  matrix[4] = cell_size[2];
+  matrix[5] = cell_size[4];
+  matrix[6] = cell_size[3];
+  matrix[7] = cell_size[4];
+  matrix[8] = cell_size[5];
+  // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
+  for (N[0]=-1;N[0]<=1;N[0]++)
+    for (N[1]=-1;N[1]<=1;N[1]++)
+      for (N[2]=-1;N[2]<=1;N[2]++) {
+        // create the translation vector
+        TranslationVector.Zero();
+        for (int i=NDIM;i--;)
+          TranslationVector.x[i] = (double)N[i];
+        TranslationVector.MatrixMultiplication(matrix);
+        // add onto the original vector to compare with
+        Shiftedy.CopyVector(y);
+        Shiftedy.AddVector(&TranslationVector);
+        // get distance and compare with minimum so far
+        tmp = DistanceSquared(&Shiftedy);
+        if (tmp < res) res = tmp;
+      }
+  return (res);
+};
+
 /** Keeps the vector in a periodic cell, defined by the symmetric \a *matrix.
  * \param *out ofstream for debugging messages
@@ -221 +259 @@
 double Vector::Angle(const Vector *y) const
 {
-  return acos(this->ScalarProduct(y)/Norm()/y->Norm());
+  double angle = this->ScalarProduct(y)/Norm()/y->Norm();
+  // -1-MYEPSILON occured due to numerical imprecision, catch ...
+  //cout << Verbose(2) << "INFO: acos(-1) = " << acos(-1) << ", acos(-1+MYEPSILON) = " << acos(-1+MYEPSILON) << ", acos(-1-MYEPSILON) = " << acos(-1-MYEPSILON) << "." << endl;
+  if (angle < -1)
+    angle = -1;
+  if (angle > 1)
+    angle = 1;
+  return acos(angle);
 };
 
@@ -313 +358 @@
 };
 
-/** Prints a 3dim vector to a stream.
- * \param ost output stream
- * \param v Vector to be printed
- * \return output stream
- */
 ostream& operator<<(ostream& ost,Vector& m)
 {
@@ -375 +415 @@
 };
 
-/** Do a matrix multiplication with \a *matrix' inverse.
+/** Calculate the inverse of a 3x3 matrix.
  * \param *matrix NDIM_NDIM array
  */
-void Vector::InverseMatrixMultiplication(double *A)
-{
-  Vector C;
-  double B[NDIM*NDIM];
+double * Vector::InverseMatrix(double *A)
+{
+  double *B = (double *) Malloc(sizeof(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
@@ -397 +438 @@
     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;
+};
+
+/** Do a matrix multiplication with \a *matrix' inverse.
+ * \param *matrix NDIM_NDIM array
+ */
+void Vector::InverseMatrixMultiplication(double *A)
+{
+  Vector C;
+  double B[NDIM*NDIM];
+  double detA = RDET3(A);
+  double detAReci;
+
+  // 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
 
     // do the matrix multiplication
@@ -457 +524 @@
   x2.CopyVector(y3);
   x2.SubtractVector(y2);
-  if ((x1.Norm()==0) || (x2.Norm()==0)) {
+  if ((fabs(x1.Norm()) < MYEPSILON) || (fabs(x2.Norm()) < MYEPSILON) || (fabs(x1.Angle(&x2)) < MYEPSILON)) {
     cout << Verbose(4) << "Given vectors are linear dependent." << endl;
     return false;
@@ -491 +558 @@
   x2.CopyVector(y2);
   Zero();
-  if ((x1.Norm()==0) || (x2.Norm()==0)) {
+  if ((fabs(x1.Norm()) < MYEPSILON) || (fabs(x2.Norm()) < MYEPSILON) || (fabs(x1.Angle(&x2)) < MYEPSILON)) {
     cout << Verbose(4) << "Given vectors are linear dependent." << endl;
     return false;