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

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"
@@ -22 +16 @@
 #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>
 
@@ -48 +45 @@
 
   // go through all atoms
-  ActOnAllVectors( &Vector::SubtractVector, *Center);
-  ActOnAllVectors( &Vector::SubtractVector, *CenterBox);
-  atoms.transformNodes(boost::bind(&Box::WrapPeriodically,domain,_1));
+  BOOST_FOREACH(atom* iter, atoms){
+    *iter -= *Center;
+    *iter -= *CenterBox;
+    iter->setPosition(domain.WrapPeriodically(iter->getPosition()));
+  }
 
   delete(Center);
@@ -66 +65 @@
   Box &domain = World::getInstance().getDomain();
 
-  atoms.transformNodes(boost::bind(&Box::WrapPeriodically,domain,_1));
+  // go through all atoms
+  BOOST_FOREACH(atom* iter, atoms){
+    iter->setPosition(domain.WrapPeriodically(iter->getPosition()));
+  }
 
   return status;
@@ -83 +85 @@
   if (iter != end()) { //list not empty?
     for (int i=NDIM;i--;) {
-      max->at(i) = (*iter)->x[i];
-      min->at(i) = (*iter)->x[i];
+      max->at(i) = (*iter)->at(i);
+      min->at(i) = (*iter)->at(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)->x[i]) ? (*iter)->x[i] : max->at(i);
-        min->at(i) = (min->at(i) > (*iter)->x[i]) ? (*iter)->x[i] : min->at(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);
       }
     }
@@ -121 +123 @@
     for (; iter != end(); ++iter) {  // continue with second if present
       Num++;
-      Center += (*iter)->x;
+      Center += (*iter)->getPosition();
     }
     Center.Scale(-1./(double)Num); // divide through total number (and sign for direction)
@@ -143 +145 @@
     for (; iter != end(); ++iter) {  // continue with second if present
       Num++;
-      (*a) += (*iter)->x;
+      (*a) += (*iter)->getPosition();
     }
     a->Scale(1./(double)Num); // divide through total mass (and sign for direction)
@@ -165 +167 @@
  * \return pointer to center of gravity vector
  */
-Vector * molecule::DetermineCenterOfGravity() const
+Vector * molecule::DetermineCenterOfGravity()
 {
   molecule::const_iterator iter = begin();  // start at first in list
@@ -176 +178 @@
   if (iter != end()) {   //list not empty?
     for (; iter != end(); ++iter) {  // continue with second if present
-      Num += (*iter)->type->mass;
-      tmp = (*iter)->type->mass * (*iter)->x;
+      Num += (*iter)->getType()->mass;
+      tmp = (*iter)->getType()->mass * (*iter)->getPosition();
       (*a) += tmp;
     }
@@ -221 +223 @@
     for (int j=0;j<MDSteps;j++)
       (*iter)->Trajectory.R.at(j).ScaleAll(*factor);
-    (*iter)->x.ScaleAll(*factor);
+    (*iter)->ScaleAll(*factor);
   }
 };
@@ -233 +235 @@
     for (int j=0;j<MDSteps;j++)
       (*iter)->Trajectory.R.at(j) += (*trans);
-    (*iter)->x += (*trans);
+    *(*iter) += (*trans);
   }
 };
@@ -246 +248 @@
 
   // go through all atoms
-  ActOnAllVectors( &Vector::AddVector, *trans);
-  atoms.transformNodes(boost::bind(&Box::WrapPeriodically,domain,_1));
+  BOOST_FOREACH(atom* iter, atoms){
+    *iter += *trans;
+    iter->setPosition(domain.WrapPeriodically(iter->getPosition()));
+  }
 
 };
@@ -259 +263 @@
   OBSERVE;
   Plane p(*n,0);
-  atoms.transformNodes(boost::bind(&Plane::mirrorVector,p,_1));
+  BOOST_FOREACH(atom* iter, atoms ){
+    iter->setPosition(p.mirrorVector(iter->getPosition()));
+  }
 };
 
@@ -278 +284 @@
     for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
 #ifdef ADDHYDROGEN
-      if ((*iter)->type->Z != 1) {
+      if ((*iter)->getType()->Z != 1) {
 #endif
-        Testvector = inversematrix * (*iter)->x;
+        Testvector = inversematrix * (*iter)->getPosition();
         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)->x[j] - (*Runner)->GetOtherAtom(*iter)->x[j];
+              tmp = (*iter)->at(j) - (*Runner)->GetOtherAtom(*iter)->at(j);
               if ((fabs(tmp)) > BondDistance) {
                 flag = false;
@@ -303 +309 @@
         // now also change all hydrogens
         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;
+          if ((*Runner)->GetOtherAtom((*iter))->getType()->Z == 1) {
+            Testvector = inversematrix * (*Runner)->GetOtherAtom((*iter))->getPosition();
             Testvector += Translationvector;
             Testvector *= matrix;
@@ -319 +325 @@
 };
 
+/** 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.
@@ -335 +442 @@
   DoLog(1) && (Log() << Verbose(1) << "Z-X-angle: " << alpha << " ... ");
   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
-    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];
+    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));
     for (int j=0;j<MDSteps;j++) {
       tmp = (*iter)->Trajectory.R.at(j)[0];
@@ -354 +461 @@
   DoLog(1) && (Log() << Verbose(1) << "Z-Y-angle: " << alpha << " ... ");
   for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
-    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];
+    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));
     for (int j=0;j<MDSteps;j++) {
       tmp = (*iter)->Trajectory.R.at(j)[1];
@@ -394 +501 @@
   // go through all atoms
   for (molecule::const_iterator iter = par->mol->begin(); iter != par->mol->end(); ++iter) {
-    if ((*iter)->type == ((struct lsq_params *)params)->type) { // for specific type
-      c = (*iter)->x - a;
+    if ((*iter)->getType() == ((struct lsq_params *)params)->type) { // for specific type
+      c = (*iter)->getPosition() - a;
       t = c.ScalarProduct(b);           // get direction parameter
       d = t*b;       // and create vector