Changes in src/molecule_geometry.cpp [d74077:4bb63c]

File:

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

src/molecule_geometry.cpp

@@ -10 +10 @@
 #endif
 
+#include "Helpers/helpers.hpp"
+#include "Helpers/Log.hpp"
 #include "Helpers/MemDebug.hpp"
+#include "Helpers/Verbose.hpp"
+#include "LinearAlgebra/Line.hpp"
+#include "LinearAlgebra/Matrix.hpp"
+#include "LinearAlgebra/Plane.hpp"
 
 #include "atom.hpp"
@@ -16 +22 @@
 #include "config.hpp"
 #include "element.hpp"
-#include "helpers.hpp"
 #include "leastsquaremin.hpp"
-#include "verbose.hpp"
-#include "log.hpp"
 #include "molecule.hpp"
 #include "World.hpp"
-#include "Plane.hpp"
-#include "Matrix.hpp"
+
 #include "Box.hpp"
+
 #include <boost/foreach.hpp>
 
@@ -45 +48 @@
 
   // go through all atoms
-  BOOST_FOREACH(atom* iter, atoms){
-    *iter -= *Center;
-    *iter -= *CenterBox;
-    iter->setPosition(domain.WrapPeriodically(iter->getPosition()));
-  }
+  ActOnAllVectors( &Vector::SubtractVector, *Center);
+  ActOnAllVectors( &Vector::SubtractVector, *CenterBox);
+  atoms.transformNodes(boost::bind(&Box::WrapPeriodically,domain,_1));
 
   delete(Center);
@@ -65 +66 @@
   Box &domain = World::getInstance().getDomain();
 
-  // go through all atoms
-  BOOST_FOREACH(atom* iter, atoms){
-    iter->setPosition(domain.WrapPeriodically(iter->getPosition()));
-  }
+  atoms.transformNodes(boost::bind(&Box::WrapPeriodically,domain,_1));
 
   return status;
@@ -85 +83 @@
   if (iter != end()) { //list not empty?
     for (int i=NDIM;i--;) {
-      max->at(i) = (*iter)->at(i);
-      min->at(i) = (*iter)->at(i);
+      max->at(i) = (*iter)->x[i];
+      min->at(i) = (*iter)->x[i];
     }
     for (; iter != end(); ++iter) {// continue with second if present
       //(*iter)->Output(1,1,out);
       for (int i=NDIM;i--;) {
-        max->at(i) = (max->at(i) < (*iter)->at(i)) ? (*iter)->at(i) : max->at(i);
-        min->at(i) = (min->at(i) > (*iter)->at(i)) ? (*iter)->at(i) : min->at(i);
+        max->at(i) = (max->at(i) < (*iter)->x[i]) ? (*iter)->x[i] : max->at(i);
+        min->at(i) = (min->at(i) > (*iter)->x[i]) ? (*iter)->x[i] : min->at(i);
       }
     }
@@ -123 +121 @@
     for (; iter != end(); ++iter) {  // continue with second if present
       Num++;
-      Center += (*iter)->getPosition();
+      Center += (*iter)->x;
     }
     Center.Scale(-1./(double)Num); // divide through total number (and sign for direction)
@@ -145 +143 @@
     for (; iter != end(); ++iter) {  // continue with second if present
       Num++;
-      (*a) += (*iter)->getPosition();
+      (*a) += (*iter)->x;
     }
     a->Scale(1./(double)Num); // divide through total mass (and sign for direction)
@@ -167 +165 @@
  * \return pointer to center of gravity vector
  */
-Vector * molecule::DetermineCenterOfGravity()
+Vector * molecule::DetermineCenterOfGravity() const
 {
   molecule::const_iterator iter = begin();  // start at first in list
@@ -178 +176 @@
   if (iter != end()) {   //list not empty?
     for (; iter != end(); ++iter) {  // continue with second if present
-      Num += (*iter)->getType()->mass;
-      tmp = (*iter)->getType()->mass * (*iter)->getPosition();
+      Num += (*iter)->type->mass;
+      tmp = (*iter)->type->mass * (*iter)->x;
       (*a) += tmp;
     }
@@ -223 +221 @@
     for (int j=0;j<MDSteps;j++)
       (*iter)->Trajectory.R.at(j).ScaleAll(*factor);
-    (*iter)->ScaleAll(*factor);
+    (*iter)->x.ScaleAll(*factor);
   }
 };
@@ -235 +233 @@
     for (int j=0;j<MDSteps;j++)
       (*iter)->Trajectory.R.at(j) += (*trans);
-    *(*iter) += (*trans);
+    (*iter)->x += (*trans);
   }
 };
@@ -248 +246 @@
 
   // go through all atoms
-  BOOST_FOREACH(atom* iter, atoms){
-    *iter += *trans;
-    iter->setPosition(domain.WrapPeriodically(iter->getPosition()));
-  }
+  ActOnAllVectors( &Vector::AddVector, *trans);
+  atoms.transformNodes(boost::bind(&Box::WrapPeriodically,domain,_1));
 
 };
@@ -263 +259 @@
   OBSERVE;
   Plane p(*n,0);
-  BOOST_FOREACH(atom* iter, atoms ){
-    iter->setPosition(p.mirrorVector(iter->getPosition()));
-  }
+  atoms.transformNodes(boost::bind(&Plane::mirrorVector,p,_1));
 };
 
@@ -284 +278 @@
     for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
 #ifdef ADDHYDROGEN
-      if ((*iter)->getType()->Z != 1) {
+      if ((*iter)->type->Z != 1) {
 #endif
-        Testvector = inversematrix * (*iter)->getPosition();
+        Testvector = inversematrix * (*iter)->x;
         Translationvector.Zero();
         for (BondList::const_iterator Runner = (*iter)->ListOfBonds.begin(); Runner != (*iter)->ListOfBonds.end(); (++Runner)) {
          if ((*iter)->nr < (*Runner)->GetOtherAtom((*iter))->nr) // otherwise we shift one to, the other fro and gain nothing
             for (int j=0;j<NDIM;j++) {
-              tmp = (*iter)->at(j) - (*Runner)->GetOtherAtom(*iter)->at(j);
+              tmp = (*iter)->x[j] - (*Runner)->GetOtherAtom(*iter)->x[j];
               if ((fabs(tmp)) > BondDistance) {
                 flag = false;
@@ -309 +303 @@
         // now also change all hydrogens
         for (BondList::const_iterator Runner = (*iter)->ListOfBonds.begin(); Runner != (*iter)->ListOfBonds.end(); (++Runner)) {
-          if ((*Runner)->GetOtherAtom((*iter))->getType()->Z == 1) {
-            Testvector = inversematrix * (*Runner)->GetOtherAtom((*iter))->getPosition();
+          if ((*Runner)->GetOtherAtom((*iter))->type->Z == 1) {
+            Testvector = inversematrix * (*Runner)->GetOtherAtom((*iter))->x;
             Testvector += Translationvector;
             Testvector *= matrix;
@@ -325 +319 @@
 };
 
-/** Transforms/Rotates the given molecule into its principal axis system.
- * \param *out output stream for debugging
- * \param DoRotate whether to rotate (true) or only to determine the PAS.
- * TODO treatment of trajetories missing
- */
-void molecule::PrincipalAxisSystem(bool DoRotate)
-{
-  double InertiaTensor[NDIM*NDIM];
-  Vector *CenterOfGravity = DetermineCenterOfGravity();
-
-  CenterPeriodic();
-
-  // reset inertia tensor
-  for(int i=0;i<NDIM*NDIM;i++)
-    InertiaTensor[i] = 0.;
-
-  // sum up inertia tensor
-  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
-    Vector x = (*iter)->getPosition();
-    //x.SubtractVector(CenterOfGravity);
-    InertiaTensor[0] += (*iter)->getType()->mass*(x[1]*x[1] + x[2]*x[2]);
-    InertiaTensor[1] += (*iter)->getType()->mass*(-x[0]*x[1]);
-    InertiaTensor[2] += (*iter)->getType()->mass*(-x[0]*x[2]);
-    InertiaTensor[3] += (*iter)->getType()->mass*(-x[1]*x[0]);
-    InertiaTensor[4] += (*iter)->getType()->mass*(x[0]*x[0] + x[2]*x[2]);
-    InertiaTensor[5] += (*iter)->getType()->mass*(-x[1]*x[2]);
-    InertiaTensor[6] += (*iter)->getType()->mass*(-x[2]*x[0]);
-    InertiaTensor[7] += (*iter)->getType()->mass*(-x[2]*x[1]);
-    InertiaTensor[8] += (*iter)->getType()->mass*(x[0]*x[0] + x[1]*x[1]);
-  }
-  // print InertiaTensor for debugging
-  DoLog(0) && (Log() << Verbose(0) << "The inertia tensor is:" << endl);
-  for(int i=0;i<NDIM;i++) {
-    for(int j=0;j<NDIM;j++)
-      DoLog(0) && (Log() << Verbose(0) << InertiaTensor[i*NDIM+j] << " ");
-    DoLog(0) && (Log() << Verbose(0) << endl);
-  }
-  DoLog(0) && (Log() << Verbose(0) << endl);
-
-  // diagonalize to determine principal axis system
-  gsl_eigen_symmv_workspace *T = gsl_eigen_symmv_alloc(NDIM);
-  gsl_matrix_view m = gsl_matrix_view_array(InertiaTensor, NDIM, NDIM);
-  gsl_vector *eval = gsl_vector_alloc(NDIM);
-  gsl_matrix *evec = gsl_matrix_alloc(NDIM, NDIM);
-  gsl_eigen_symmv(&m.matrix, eval, evec, T);
-  gsl_eigen_symmv_free(T);
-  gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_ABS_DESC);
-
-  for(int i=0;i<NDIM;i++) {
-    DoLog(1) && (Log() << Verbose(1) << "eigenvalue = " << gsl_vector_get(eval, i));
-    DoLog(0) && (Log() << Verbose(0) << ", eigenvector = (" << evec->data[i * evec->tda + 0] << "," << evec->data[i * evec->tda + 1] << "," << evec->data[i * evec->tda + 2] << ")" << endl);
-  }
-
-  // check whether we rotate or not
-  if (DoRotate) {
-    DoLog(1) && (Log() << Verbose(1) << "Transforming molecule into PAS ... ");
-    // the eigenvectors specify the transformation matrix
-    Matrix M = Matrix(evec->data);
-    Vector tempVector;
-    BOOST_FOREACH(atom* iter, atoms){
-      tempVector = iter->getPosition();
-      tempVector *= M;
-      iter->setPosition(tempVector);
-    }
-    DoLog(0) && (Log() << Verbose(0) << "done." << endl);
-
-    // summing anew for debugging (resulting matrix has to be diagonal!)
-    // reset inertia tensor
-    for(int i=0;i<NDIM*NDIM;i++)
-      InertiaTensor[i] = 0.;
-
-    // sum up inertia tensor
-    for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
-      Vector x = (*iter)->getPosition();
-      InertiaTensor[0] += (*iter)->getType()->mass*(x[1]*x[1] + x[2]*x[2]);
-      InertiaTensor[1] += (*iter)->getType()->mass*(-x[0]*x[1]);
-      InertiaTensor[2] += (*iter)->getType()->mass*(-x[0]*x[2]);
-      InertiaTensor[3] += (*iter)->getType()->mass*(-x[1]*x[0]);
-      InertiaTensor[4] += (*iter)->getType()->mass*(x[0]*x[0] + x[2]*x[2]);
-      InertiaTensor[5] += (*iter)->getType()->mass*(-x[1]*x[2]);
-      InertiaTensor[6] += (*iter)->getType()->mass*(-x[2]*x[0]);
-      InertiaTensor[7] += (*iter)->getType()->mass*(-x[2]*x[1]);
-      InertiaTensor[8] += (*iter)->getType()->mass*(x[0]*x[0] + x[1]*x[1]);
-    }
-    // print InertiaTensor for debugging
-    DoLog(0) && (Log() << Verbose(0) << "The inertia tensor is:" << endl);
-    for(int i=0;i<NDIM;i++) {
-      for(int j=0;j<NDIM;j++)
-        DoLog(0) && (Log() << Verbose(0) << InertiaTensor[i*NDIM+j] << " ");
-      DoLog(0) && (Log() << Verbose(0) << endl);
-    }
-    DoLog(0) && (Log() << Verbose(0) << endl);
-  }
-
-  // free everything
-  delete(CenterOfGravity);
-  gsl_vector_free(eval);
-  gsl_matrix_free(evec);
-};
-
-
 /** Align all atoms in such a manner that given vector \a *n is along z axis.
  * \param n[] alignment vector.
@@ -442 +335 @@
   DoLog(1) && (Log() << Verbose(1) << "Z-X-angle: " << alpha << " ... ");
   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
-    tmp = (*iter)->at(0);
-    (*iter)->set(0,  cos(alpha) * tmp + sin(alpha) * (*iter)->at(2));
-    (*iter)->set(2, -sin(alpha) * tmp + cos(alpha) * (*iter)->at(2));
+    tmp = (*iter)->x[0];
+    (*iter)->x[0] =  cos(alpha) * tmp + sin(alpha) * (*iter)->x[2];
+    (*iter)->x[2] = -sin(alpha) * tmp + cos(alpha) * (*iter)->x[2];
     for (int j=0;j<MDSteps;j++) {
       tmp = (*iter)->Trajectory.R.at(j)[0];
@@ -461 +354 @@
   DoLog(1) && (Log() << Verbose(1) << "Z-Y-angle: " << alpha << " ... ");
   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
-    tmp = (*iter)->at(1);
-    (*iter)->set(1,  cos(alpha) * tmp + sin(alpha) * (*iter)->at(2));
-    (*iter)->set(2, -sin(alpha) * tmp + cos(alpha) * (*iter)->at(2));
+    tmp = (*iter)->x[1];
+    (*iter)->x[1] =  cos(alpha) * tmp + sin(alpha) * (*iter)->x[2];
+    (*iter)->x[2] = -sin(alpha) * tmp + cos(alpha) * (*iter)->x[2];
     for (int j=0;j<MDSteps;j++) {
       tmp = (*iter)->Trajectory.R.at(j)[1];
@@ -501 +394 @@
   // go through all atoms
   for (molecule::const_iterator iter = par->mol->begin(); iter != par->mol->end(); ++iter) {
-    if ((*iter)->getType() == ((struct lsq_params *)params)->type) { // for specific type
-      c = (*iter)->getPosition() - a;
+    if ((*iter)->type == ((struct lsq_params *)params)->type) { // for specific type
+      c = (*iter)->x - a;
       t = c.ScalarProduct(b);           // get direction parameter
       d = t*b;       // and create vector