Changes in src/molecule_geometry.cpp [9291d04:833b15]

File:

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

src/molecule_geometry.cpp

@@ -58 +58 @@
 /************************************* Functions for class molecule *********************************/
 
+/** Returns vector pointing to center of the domain.
+ * \return pointer to center of the domain
+ */
+#ifdef HAVE_INLINE
+inline
+#else
+static
+#endif
+const Vector DetermineCenterOfBox()
+{
+  Vector a(0.5,0.5,0.5);
+  const RealSpaceMatrix &M = World::getInstance().getDomain().getM();
+  a *= M;
+  return a;
+}
 
 /** Centers the molecule in the box whose lengths are defined by vector \a *BoxLengths.
@@ -65 +80 @@
 {
   bool status = true;
-  const Vector *Center = DetermineCenterOfAll();
-  const Vector *CenterBox = DetermineCenterOfBox();
+  const Vector Center = DetermineCenterOfAll();
+  const Vector CenterBox = DetermineCenterOfBox();
   Box &domain = World::getInstance().getDomain();
 
   // go through all atoms
-  for (iterator iter = begin(); iter != end(); ++iter) {
-    if (DoLog(4) && (*Center != *CenterBox))
-      LOG(4, "INFO: atom before is at " << **iter);
-    **iter -= *Center;
-    **iter += *CenterBox;
-    if (DoLog(4) && (*Center != *CenterBox))
-      LOG(4, "INFO: atom after is at " << **iter);
-  }
+  Translate(CenterBox - Center);
   getAtomSet().transformNodes(boost::bind(&Box::enforceBoundaryConditions,domain,_1));
 
-  delete(Center);
-  delete(CenterBox);
   return status;
-};
+}
 
 
@@ -98 +104 @@
 
   return status;
-};
+}
 
 /** Centers the edge of the atoms at (0,0,0).
- * \param *out output stream for debugging
- * \param *max coordinates of other edge, specifying box dimensions.
- */
-void molecule::CenterEdge(Vector *max)
-{
-//  Info info(__func__);
-  Vector *min = new Vector;
-
-  const_iterator iter = begin();  // start at first in list
-  if (iter != end()) { //list not empty?
-    for (int i=NDIM;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)->at(i)) ? (*iter)->at(i) : max->at(i);
-        min->at(i) = (min->at(i) > (*iter)->at(i)) ? (*iter)->at(i) : min->at(i);
-      }
-    }
-    LOG(1, "INFO: Maximum is " << *max << ", Minimum is " << *min << ".");
-    min->Scale(-1.);
-    (*max) += (*min);
-    Translate(min);
-  }
-  delete(min);
-};
+ */
+void molecule::CenterEdge()
+{
+  const_iterator iter = begin();
+  if (iter != end()) {   //list not empty?
+    Vector min = (*begin())->getPosition();
+    for (;iter != end(); ++iter) { // continue with second if present
+      const Vector &currentPos = (*iter)->getPosition();
+      for (size_t i=0;i<NDIM;++i)
+        if (min[i] > currentPos[i])
+          min[i] = currentPos[i];
+    }
+    Translate(-1.*min);
+  }
+}
 
 /** Centers the center of the atoms at (0,0,0).
@@ -147 +140 @@
     }
     Center.Scale(-1./(double)Num); // divide through total number (and sign for direction)
-    Translate(&Center);
-  }
-};
+    Translate(Center);
+  }
+}
 
 /** Returns vector pointing to center of all atoms.
  * \return pointer to center of all vector
  */
-Vector * molecule::DetermineCenterOfAll() const
+const Vector molecule::DetermineCenterOfAll() const
 {
   const_iterator iter = begin();  // start at first in list
-  Vector *a = new Vector();
+  Vector a;
   double Num = 0;
 
-  a->Zero();
+  a.Zero();
 
   if (iter != end()) {   //list not empty?
     for (; iter != end(); ++iter) {  // continue with second if present
       Num++;
-      (*a) += (*iter)->getPosition();
-    }
-    a->Scale(1./(double)Num); // divide through total mass (and sign for direction)
+      a += (*iter)->getPosition();
+    }
+    a.Scale(1./(double)Num); // divide through total mass (and sign for direction)
   }
   return a;
-};
-
-/** Returns vector pointing to center of the domain.
- * \return pointer to center of the domain
- */
-Vector * molecule::DetermineCenterOfBox() const
-{
-  Vector *a = new Vector(0.5,0.5,0.5);
-  const RealSpaceMatrix &M = World::getInstance().getDomain().getM();
-  (*a) *= M;
-  return a;
-};
+}
+
 
 /** Returns vector pointing to center of gravity.
@@ -187 +170 @@
  * \return pointer to center of gravity vector
  */
-Vector * molecule::DetermineCenterOfGravity() const
+const Vector molecule::DetermineCenterOfGravity() const
 {
   const_iterator iter = begin();  // start at first in list
-  Vector *a = new Vector();
+  Vector a;
   Vector tmp;
   double Num = 0;
 
-  a->Zero();
+  a.Zero();
 
   if (iter != end()) {   //list not empty?
@@ -200 +183 @@
       Num += (*iter)->getType()->getMass();
       tmp = (*iter)->getType()->getMass() * (*iter)->getPosition();
-      (*a) += tmp;
-    }
-    a->Scale(1./Num); // divide through total mass
-  }
-  LOG(1, "INFO: Resulting center of gravity: " << *a << ".");
+      a += tmp;
+    }
+    a.Scale(1./Num); // divide through total mass
+  }
+  LOG(1, "INFO: Resulting center of gravity: " << a << ".");
   return a;
-};
+}
 
 /** Centers the center of gravity of the atoms at (0,0,0).
@@ -216 +199 @@
   Vector NewCenter;
   DeterminePeriodicCenter(NewCenter);
-  // go through all atoms
-  for (iterator iter = begin(); iter != end(); ++iter) {
-    **iter -= NewCenter;
-  }
-};
+  Translate(-1.*NewCenter);
+}
 
 
@@ -227 +207 @@
  * \param *center return vector for translation vector
  */
-void molecule::CenterAtVector(Vector *newcenter)
-{
-  // go through all atoms
-  for (iterator iter = begin(); iter != end(); ++iter) {
-    **iter -= *newcenter;
-  }
-};
+void molecule::CenterAtVector(const Vector &newcenter)
+{
+  Translate(-1.*newcenter);
+}
 
 /** Calculate the inertia tensor of a the molecule.
@@ -242 +219 @@
 {
   RealSpaceMatrix InertiaTensor;
-  Vector *CenterOfGravity = DetermineCenterOfGravity();
+  const Vector CenterOfGravity = DetermineCenterOfGravity();
 
   // reset inertia tensor
@@ -250 +227 @@
   for (const_iterator iter = begin(); iter != end(); ++iter) {
     Vector x = (*iter)->getPosition();
-    x -= *CenterOfGravity;
+    x -= CenterOfGravity;
     const double mass = (*iter)->getType()->getMass();
     InertiaTensor.at(0,0) += mass*(x[1]*x[1] + x[2]*x[2]);
@@ -265 +242 @@
   LOG(1, "INFO: The inertia tensor of molecule " << getName() <<  " is:" << InertiaTensor);
 
-  delete CenterOfGravity;
   return InertiaTensor;
 }
@@ -276 +252 @@
 void molecule::RotateToPrincipalAxisSystem(const Vector &Axis)
 {
-  Vector *CenterOfGravity = DetermineCenterOfGravity();
+  const Vector CenterOfGravity = DetermineCenterOfGravity();
   RealSpaceMatrix InertiaTensor = getInertiaTensor();
 
@@ -288 +264 @@
 
   // obtain first column, eigenvector to biggest eigenvalue
-  Vector BiggestEigenvector(InertiaTensor.column(Eigenvalues.SmallestComponent()));
+  const Vector BiggestEigenvector(InertiaTensor.column(Eigenvalues.SmallestComponent()));
   Vector DesiredAxis(Axis.getNormalized());
 
@@ -302 +278 @@
   // and rotate
   for (iterator iter = begin(); iter != end(); ++iter) {
-    *(*iter) -= *CenterOfGravity;
+    *(*iter) -= CenterOfGravity;
     (*iter)->setPosition(RotationAxis.rotateVector((*iter)->getPosition(), alpha));
-    *(*iter) += *CenterOfGravity;
+    *(*iter) += CenterOfGravity;
   }
   LOG(0, "STATUS: done.");
-
-  delete CenterOfGravity;
 }
 
@@ -319 +293 @@
  * x=(**factor) * x (as suggested by comment)
  */
-void molecule::Scale(const double ** const factor)
+void molecule::Scale(const double *factor)
 {
   for (iterator iter = begin(); iter != end(); ++iter) {
     for (size_t j=0;j<(*iter)->getTrajectorySize();j++) {
       Vector temp = (*iter)->getPositionAtStep(j);
-      temp.ScaleAll(*factor);
+      temp.ScaleAll(factor);
       (*iter)->setPositionAtStep(j,temp);
     }
@@ -333 +307 @@
  * \param trans[] translation vector.
  */
-void molecule::Translate(const Vector *trans)
-{
-  for (iterator iter = begin(); iter != end(); ++iter) {
-    for (size_t j=0;j<(*iter)->getTrajectorySize();j++) {
-      (*iter)->setPositionAtStep(j, (*iter)->getPositionAtStep(j) + (*trans));
-    }
-  }
+void molecule::Translate(const Vector &trans)
+{
+  getAtomSet().translate(trans);
 };
 
@@ -346 +316 @@
  * TODO treatment of trajectories missing
  */
-void molecule::TranslatePeriodically(const Vector *trans)
-{
+void molecule::TranslatePeriodically(const Vector &trans)
+{
+  Translate(trans);
   Box &domain = World::getInstance().getDomain();
-
-  // go through all atoms
-  for (iterator iter = begin(); iter != end(); ++iter) {
-    **iter += *trans;
-  }
   getAtomSet().transformNodes(boost::bind(&Box::enforceBoundaryConditions,domain,_1));
-
 };
 
@@ -362 +327 @@
  * \param n[] normal vector of mirror plane.
  */
-void molecule::Mirror(const Vector *n)
-{
-  OBSERVE;
-  Plane p(*n,0);
+void molecule::Mirror(const Vector &n)
+{
+  Plane p(n,0);
   getAtomSet().transformNodes(boost::bind(&Plane::mirrorVector,p,_1));
 };
@@ -381 +345 @@
   Vector Testvector, Translationvector;
   Vector Center;
-  BondGraph *BG = World::getInstance().getBondGraph();
+  const BondGraph * const BG = World::getInstance().getBondGraph();
 
   do {
@@ -432 +396 @@
 
   Center.Scale(1./static_cast<double>(getAtomCount()));
-  CenterAtVector(&Center);
+  CenterAtVector(Center);
 };
 
@@ -438 +402 @@
  * \param n[] alignment vector.
  */
-void molecule::Align(Vector *n)
+void molecule::Align(const Vector &n)
 {
   double alpha, tmp;
   Vector z_axis;
+  Vector alignment(n);
   z_axis[0] = 0.;
   z_axis[1] = 0.;
@@ -448 +413 @@
   // rotate on z-x plane
   LOG(0, "Begin of Aligning all atoms.");
-  alpha = atan(-n->at(0)/n->at(2));
+  alpha = atan(-alignment.at(0)/alignment.at(2));
   LOG(1, "INFO: Z-X-angle: " << alpha << " ... ");
   for (iterator iter = begin(); iter != end(); ++iter) {
@@ -462 +427 @@
   }
   // rotate n vector
-  tmp = n->at(0);
-  n->at(0) =  cos(alpha) * tmp +  sin(alpha) * n->at(2);
-  n->at(2) = -sin(alpha) * tmp +  cos(alpha) * n->at(2);
-  LOG(1, "alignment vector after first rotation: " << n);
+  tmp = alignment.at(0);
+  alignment.at(0) =  cos(alpha) * tmp +  sin(alpha) * alignment.at(2);
+  alignment.at(2) = -sin(alpha) * tmp +  cos(alpha) * alignment.at(2);
+  LOG(1, "alignment vector after first rotation: " << alignment);
 
   // rotate on z-y plane
-  alpha = atan(-n->at(1)/n->at(2));
+  alpha = atan(-alignment.at(1)/alignment.at(2));
   LOG(1, "INFO: Z-Y-angle: " << alpha << " ... ");
   for (iterator iter = begin(); iter != end(); ++iter) {
@@ -482 +447 @@
   }
   // rotate n vector (for consistency check)
-  tmp = n->at(1);
-  n->at(1) =  cos(alpha) * tmp +  sin(alpha) * n->at(2);
-  n->at(2) = -sin(alpha) * tmp +  cos(alpha) * n->at(2);
-
-
-  LOG(1, "alignment vector after second rotation: " << n);
+  tmp = alignment.at(1);
+  alignment.at(1) =  cos(alpha) * tmp +  sin(alpha) * alignment.at(2);
+  alignment.at(2) = -sin(alpha) * tmp +  cos(alpha) * alignment.at(2);
+
+  LOG(1, "alignment vector after second rotation: " << alignment);
   LOG(0, "End of Aligning all atoms.");
 };