Changes in src/molecule_geometry.cpp [aafd77:bdc91e]

File:

• src/molecule_geometry.cpp (modified) (16 diffs)

src/molecule_geometry.cpp

@@ -5 +5 @@
  *      Author: heber
  */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
 
 #include "Helpers/MemDebug.hpp"
@@ -18 +14 @@
 #include "helpers.hpp"
 #include "leastsquaremin.hpp"
-#include "verbose.hpp"
 #include "log.hpp"
+#include "memoryallocator.hpp"
 #include "molecule.hpp"
 #include "World.hpp"
 #include "Plane.hpp"
-#include "Matrix.hpp"
-#include "Box.hpp"
 #include <boost/foreach.hpp>
-
-#include <gsl/gsl_eigen.h>
-#include <gsl/gsl_multimin.h>
 
 
@@ -42 +33 @@
   const Vector *Center = DetermineCenterOfAll();
   const Vector *CenterBox = DetermineCenterOfBox();
-  Box &domain = World::getInstance().getDomain();
+  double * const cell_size = World::getInstance().getDomain();
+  double *M = ReturnFullMatrixforSymmetric(cell_size);
+  double *Minv = InverseMatrix(M);
 
   // go through all atoms
   ActOnAllVectors( &Vector::SubtractVector, *Center);
   ActOnAllVectors( &Vector::SubtractVector, *CenterBox);
-  BOOST_FOREACH(atom* iter, atoms){
-    *iter->node = domain.WrapPeriodically(*iter->node);
-  }
-
+  ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
+
+  delete[](M);
+  delete[](Minv);
   delete(Center);
   delete(CenterBox);
@@ -63 +56 @@
 {
   bool status = true;
-  Box &domain = World::getInstance().getDomain();
+  double * const cell_size = World::getInstance().getDomain();
+  double *M = ReturnFullMatrixforSymmetric(cell_size);
+  double *Minv = InverseMatrix(M);
 
   // go through all atoms
-  BOOST_FOREACH(atom* iter, atoms){
-    *iter->node = domain.WrapPeriodically(*iter->node);
-  }
-
+  ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
+
+  delete[](M);
+  delete[](Minv);
   return status;
 };
@@ -125 +120 @@
       Center += (*iter)->x;
     }
-    Center.Scale(-1./Num); // divide through total number (and sign for direction)
+    Center.Scale(-1./(double)Num); // divide through total number (and sign for direction)
     Translate(&Center);
     Center.Zero();
@@ -147 +142 @@
       (*a) += (*iter)->x;
     }
-    a->Scale(1./Num); // divide through total mass (and sign for direction)
+    a->Scale(1./(double)Num); // divide through total mass (and sign for direction)
   }
   return a;
@@ -158 +153 @@
 {
   Vector *a = new Vector(0.5,0.5,0.5);
-  const Matrix &M = World::getInstance().getDomain().getM();
-  (*a) *= M;
+
+  const double *cell_size = World::getInstance().getDomain();
+  double *M = ReturnFullMatrixforSymmetric(cell_size);
+  a->MatrixMultiplication(M);
+  delete[](M);
+
   return a;
 };
@@ -182 +181 @@
       (*a) += tmp;
     }
-    a->Scale(1./Num); // divide through total mass (and sign for direction)
+    a->Scale(1./Num); // divide through total mass
   }
 //  Log() << Verbose(1) << "Resulting center of gravity: ";
@@ -245 +244 @@
 void molecule::TranslatePeriodically(const Vector *trans)
 {
-  Box &domain = World::getInstance().getDomain();
+  double * const cell_size = World::getInstance().getDomain();
+  double *M = ReturnFullMatrixforSymmetric(cell_size);
+  double *Minv = InverseMatrix(M);
 
   // go through all atoms
   ActOnAllVectors( &Vector::AddVector, *trans);
-  BOOST_FOREACH(atom* iter, atoms){
-    *iter->node = domain.WrapPeriodically(*iter->node);
-  }
-
+  ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
+
+  delete[](M);
+  delete[](Minv);
 };
 
@@ -263 +264 @@
   OBSERVE;
   Plane p(*n,0);
-  BOOST_FOREACH(atom* iter, atoms ){
+  BOOST_FOREACH( atom* iter, atoms ){
     (*iter->node) = p.mirrorVector(*iter->node);
   }
@@ -273 +274 @@
 void molecule::DeterminePeriodicCenter(Vector &center)
 {
-  const Matrix &matrix = World::getInstance().getDomain().getM();
-  const Matrix &inversematrix = World::getInstance().getDomain().getM();
+  double * const cell_size = World::getInstance().getDomain();
+  double *matrix = ReturnFullMatrixforSymmetric(cell_size);
+  double *inversematrix = InverseMatrix(matrix);
   double tmp;
   bool flag;
@@ -286 +288 @@
       if ((*iter)->type->Z != 1) {
 #endif
-        Testvector = inversematrix * (*iter)->x;
+        Testvector = (*iter)->x;
+        Testvector.MatrixMultiplication(inversematrix);
         Translationvector.Zero();
         for (BondList::const_iterator Runner = (*iter)->ListOfBonds.begin(); Runner != (*iter)->ListOfBonds.end(); (++Runner)) {
@@ -303 +306 @@
         }
         Testvector += Translationvector;
-        Testvector *= matrix;
+        Testvector.MatrixMultiplication(matrix);
         Center += Testvector;
         Log() << Verbose(1) << "vector is: " << Testvector << endl;
@@ -310 +313 @@
         for (BondList::const_iterator Runner = (*iter)->ListOfBonds.begin(); Runner != (*iter)->ListOfBonds.end(); (++Runner)) {
           if ((*Runner)->GetOtherAtom((*iter))->type->Z == 1) {
-            Testvector = inversematrix * (*Runner)->GetOtherAtom((*iter))->x;
+            Testvector = (*Runner)->GetOtherAtom((*iter))->x;
+            Testvector.MatrixMultiplication(inversematrix);
             Testvector += Translationvector;
-            Testvector *= matrix;
+            Testvector.MatrixMultiplication(matrix);
             Center += Testvector;
             Log() << Verbose(1) << "Hydrogen vector is: " << Testvector << endl;
@@ -321 +325 @@
     }
   } while (!flag);
+  delete[](matrix);
+  delete[](inversematrix);
 
   Center.Scale(1./static_cast<double>(getAtomCount()));
@@ -382 +388 @@
     DoLog(1) && (Log() << Verbose(1) << "Transforming molecule into PAS ... ");
     // the eigenvectors specify the transformation matrix
-    Matrix M = Matrix(evec->data);
-    BOOST_FOREACH(atom* iter, atoms){
-      (*iter->node) *= M;
-    }
+    ActOnAllVectors( &Vector::MatrixMultiplication, (const double *) evec->data );
     DoLog(0) && (Log() << Verbose(0) << "done." << endl);