Changeset 6fb785 for molecuilder/src/vector.cpp

Timestamp:

Jul 10, 2009, 12:48:05 PM (16 years ago)

Author:

Frederik Heber <heber@…>

Children:

22f587

Parents:

a8c8ba

Message:

molecule::CenterInBox puts atoms now periodically into the given box, new function molecule::TranslatePeriodically, BUGFIX: molecule::ReturnFullMatrixforSymmetrical()

• molecule::CenterInBox() has no more a vector as a parameter, but instead enforces the periodicity of the simulation box, i.e. all atoms out of bounds are put back in with wrap-around at boundaries. Call of function was changed in everywhere.
• in ParseCommandLineParameters() a SetBoxDimension was missing in certain Center...() commands.
• new function molecule::TranslatePeriodically translates all atoms of a molecule while adhering to the periodicity of the domain
• new function vector::InverseMatrix() returns the hard-encoded inverse of 3x3 real matrix
• BUGFIX: molecule::ReturnFullMatrixforSymmetrical()'s assignment from 6-doubles to 9-doubles was all wrong (symmetric to full 3x3 matrix)

File:

• molecuilder/src/vector.cpp (modified) (1 diff)

molecuilder/src/vector.cpp

@@ -413 +413 @@
         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(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
+    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;
 };