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

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