Changes in src/molecule_geometry.cpp [a67d19:7fd416]

File:

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

src/molecule_geometry.cpp

@@ -27 +27 @@
   bool status = true;
   const Vector *Center = DetermineCenterOfAll();
-  double * const cell_size = World::get()->cell_size;
+  const Vector *CenterBox = DetermineCenterOfBox();
+  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, *Center);
+  ActOnAllVectors( &Vector::SubtractVector, *CenterBox);
   ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
 
-  Free(&M);
-  Free(&Minv);
+  delete[](M);
+  delete[](Minv);
   delete(Center);
   return status;
@@ -48 +50 @@
 {
   bool status = true;
-  double * const cell_size = World::get()->cell_size;
+  double * const cell_size = World::getInstance().getDomain();
   double *M = ReturnFullMatrixforSymmetric(cell_size);
   double *Minv = InverseMatrix(M);
@@ -55 +57 @@
   ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
 
-  Free(&M);
-  Free(&Minv);
+  delete[](M);
+  delete[](Minv);
   return status;
 };
@@ -69 +71 @@
 
 //  Log() << Verbose(3) << "Begin of CenterEdge." << endl;
-  atom *ptr = start->next;  // start at first in list
-  if (ptr != end) {   //list not empty?
+  molecule::const_iterator iter = begin();  // start at first in list
+  if (iter != end()) { //list not empty?
     for (int i=NDIM;i--;) {
-      max->x[i] = ptr->x.x[i];
-      min->x[i] = ptr->x.x[i];
-    }
-    while (ptr->next != end) {  // continue with second if present
-      ptr = ptr->next;
-      //ptr->Output(1,1,out);
+      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->x[i] = (max->x[i] < ptr->x.x[i]) ? ptr->x.x[i] : max->x[i];
-        min->x[i] = (min->x[i] > ptr->x.x[i]) ? ptr->x.x[i] : min->x[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);
       }
     }
@@ -89 +90 @@
 //    Log() << Verbose(0) << endl;
     min->Scale(-1.);
-    max->AddVector(min);
+    (*max) += (*min);
     Translate(min);
     Center.Zero();
@@ -104 +105 @@
 {
   int Num = 0;
-  atom *ptr = start;  // start at first in list
+  molecule::const_iterator iter = begin();  // start at first in list
 
   Center.Zero();
 
-  if (ptr->next != end) {   //list not empty?
-    while (ptr->next != end) {  // continue with second if present
-      ptr = ptr->next;
+  if (iter != end()) {   //list not empty?
+    for (; iter != end(); ++iter) {  // continue with second if present
       Num++;
-      Center.AddVector(&ptr->x);
+      Center += (*iter)->x;
     }
     Center.Scale(-1./Num); // divide through total number (and sign for direction)
@@ -125 +125 @@
 Vector * molecule::DetermineCenterOfAll() const
 {
-  atom *ptr = start->next;  // start at first in list
+  molecule::const_iterator iter = begin();  // start at first in list
+  Vector *a = new Vector();
+  double Num = 0;
+
+  a->Zero();
+
+  if (iter != end()) {   //list not empty?
+    for (; iter != end(); ++iter) {  // continue with second if present
+      Num++;
+      (*a) += (*iter)->x;
+    }
+    a->Scale(1./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 double *cell_size = World::getInstance().getDomain();
+  double *M = ReturnFullMatrixforSymmetric(cell_size);
+  a->MatrixMultiplication(M);
+  delete[](M);
+
+  return a;
+};
+
+/** Returns vector pointing to center of gravity.
+ * \param *out output stream for debugging
+ * \return pointer to center of gravity vector
+ */
+Vector * molecule::DetermineCenterOfGravity()
+{
+  molecule::const_iterator iter = begin();  // start at first in list
   Vector *a = new Vector();
   Vector tmp;
@@ -132 +169 @@
   a->Zero();
 
-  if (ptr != end) {   //list not empty?
-    while (ptr->next != end) {  // continue with second if present
-      ptr = ptr->next;
-      Num += 1.;
-      tmp.CopyVector(&ptr->x);
-      a->AddVector(&tmp);
+  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;
+      (*a) += tmp;
     }
     a->Scale(1./Num); // divide through total mass (and sign for direction)
-  }
-  return a;
-};
-
-/** Returns vector pointing to center of gravity.
- * \param *out output stream for debugging
- * \return pointer to center of gravity vector
- */
-Vector * molecule::DetermineCenterOfGravity()
-{
-  atom *ptr = start->next;  // start at first in list
-  Vector *a = new Vector();
-  Vector tmp;
-  double Num = 0;
-
-  a->Zero();
-
-  if (ptr != end) {   //list not empty?
-    while (ptr->next != end) {  // continue with second if present
-      ptr = ptr->next;
-      Num += ptr->type->mass;
-      tmp.CopyVector(&ptr->x);
-      tmp.Scale(ptr->type->mass);  // scale by mass
-      a->AddVector(&tmp);
-    }
-    a->Scale(-1./Num); // divide through total mass (and sign for direction)
   }
 //  Log() << Verbose(1) << "Resulting center of gravity: ";
@@ -189 +199 @@
 void molecule::CenterAtVector(Vector *newcenter)
 {
-  Center.CopyVector(newcenter);
+  Center = *newcenter;
 };
 
@@ -195 +205 @@
 /** Scales all atoms by \a *factor.
  * \param *factor pointer to scaling factor
+ *
+ * TODO: Is this realy what is meant, i.e.
+ * x=(x[0]*factor[0],x[1]*factor[1],x[2]*factor[2]) (current impl)
+ * or rather
+ * x=(**factor) * x (as suggested by comment)
  */
 void molecule::Scale(const double ** const factor)
 {
-  atom *ptr = start;
-
-  while (ptr->next != end) {
-    ptr = ptr->next;
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
     for (int j=0;j<MDSteps;j++)
-      ptr->Trajectory.R.at(j).Scale(factor);
-    ptr->x.Scale(factor);
+      (*iter)->Trajectory.R.at(j).ScaleAll(*factor);
+    (*iter)->x.ScaleAll(*factor);
   }
 };
@@ -213 +225 @@
 void molecule::Translate(const Vector *trans)
 {
-  atom *ptr = start;
-
-  while (ptr->next != end) {
-    ptr = ptr->next;
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
     for (int j=0;j<MDSteps;j++)
-      ptr->Trajectory.R.at(j).Translate(trans);
-    ptr->x.Translate(trans);
+      (*iter)->Trajectory.R.at(j) += (*trans);
+    (*iter)->x += (*trans);
   }
 };
@@ -229 +238 @@
 void molecule::TranslatePeriodically(const Vector *trans)
 {
-  double * const cell_size = World::get()->cell_size;
+  double * const cell_size = World::getInstance().getDomain();
   double *M = ReturnFullMatrixforSymmetric(cell_size);
   double *Minv = InverseMatrix(M);
 
   // go through all atoms
-  ActOnAllVectors( &Vector::SubtractVector, trans);
+  ActOnAllVectors( &Vector::AddVector, *trans);
   ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
 
-  Free(&M);
-  Free(&Minv);
+  delete[](M);
+  delete[](Minv);
 };
 
@@ -247 +256 @@
 void molecule::Mirror(const Vector *n)
 {
-  ActOnAllVectors( &Vector::Mirror, n );
+  ActOnAllVectors( &Vector::Mirror, *n );
 };
 
@@ -255 +264 @@
 void molecule::DeterminePeriodicCenter(Vector &center)
 {
-  atom *Walker = start;
-  double * const cell_size = World::get()->cell_size;
+  double * const cell_size = World::getInstance().getDomain();
   double *matrix = ReturnFullMatrixforSymmetric(cell_size);
-  double *inversematrix = InverseMatrix(cell_size);
+  double *inversematrix = InverseMatrix(matrix);
   double tmp;
   bool flag;
@@ -266 +274 @@
     Center.Zero();
     flag = true;
-    while (Walker->next != end) {
-      Walker = Walker->next;
+    for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
 #ifdef ADDHYDROGEN
-      if (Walker->type->Z != 1) {
+      if ((*iter)->type->Z != 1) {
 #endif
-        Testvector.CopyVector(&Walker->x);
+        Testvector = (*iter)->x;
         Testvector.MatrixMultiplication(inversematrix);
         Translationvector.Zero();
-        for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
-         if (Walker->nr < (*Runner)->GetOtherAtom(Walker)->nr) // otherwise we shift one to, the other fro and gain nothing
+        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 = Walker->x.x[j] - (*Runner)->GetOtherAtom(Walker)->x.x[j];
+              tmp = (*iter)->x[j] - (*Runner)->GetOtherAtom(*iter)->x[j];
               if ((fabs(tmp)) > BondDistance) {
                 flag = false;
-                DoLog(0) && (Log() << Verbose(0) << "Hit: atom " << Walker->Name << " in bond " << *(*Runner) << " has to be shifted due to " << tmp << "." << endl);
+                DoLog(0) && (Log() << Verbose(0) << "Hit: atom " << (*iter)->getName() << " in bond " << *(*Runner) << " has to be shifted due to " << tmp << "." << endl);
                 if (tmp > 0)
-                  Translationvector.x[j] -= 1.;
+                  Translationvector[j] -= 1.;
                 else
-                  Translationvector.x[j] += 1.;
+                  Translationvector[j] += 1.;
               }
             }
         }
-        Testvector.AddVector(&Translationvector);
+        Testvector += Translationvector;
         Testvector.MatrixMultiplication(matrix);
-        Center.AddVector(&Testvector);
-        DoLog(1) && (Log() << Verbose(1) << "vector is: ");
-        Testvector.Output();
-        DoLog(0) && (Log() << Verbose(0) << endl);
+        Center += Testvector;
+        Log() << Verbose(1) << "vector is: " << Testvector << endl;
 #ifdef ADDHYDROGEN
         // now also change all hydrogens
-        for (BondList::const_iterator Runner = Walker->ListOfBonds.begin(); Runner != Walker->ListOfBonds.end(); (++Runner)) {
-          if ((*Runner)->GetOtherAtom(Walker)->type->Z == 1) {
-            Testvector.CopyVector(&(*Runner)->GetOtherAtom(Walker)->x);
+        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.AddVector(&Translationvector);
+            Testvector += Translationvector;
             Testvector.MatrixMultiplication(matrix);
-            Center.AddVector(&Testvector);
-            DoLog(1) && (Log() << Verbose(1) << "Hydrogen vector is: ");
-            Testvector.Output();
-            DoLog(0) && (Log() << Verbose(0) << endl);
+            Center += Testvector;
+            Log() << Verbose(1) << "Hydrogen vector is: " << Testvector << endl;
           }
         }
@@ -312 +315 @@
     }
   } while (!flag);
-  Free(&matrix);
-  Free(&inversematrix);
-
-  Center.Scale(1./(double)AtomCount);
+  delete[](matrix);
+  delete[](inversematrix);
+
+  Center.Scale(1./static_cast<double>(getAtomCount()));
 };
 
@@ -325 +328 @@
 void molecule::PrincipalAxisSystem(bool DoRotate)
 {
-  atom *ptr = start;  // start at first in list
   double InertiaTensor[NDIM*NDIM];
   Vector *CenterOfGravity = DetermineCenterOfGravity();
@@ -336 +338 @@
 
   // sum up inertia tensor
-  while (ptr->next != end) {
-    ptr = ptr->next;
-    Vector x;
-    x.CopyVector(&ptr->x);
+  for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+    Vector x = (*iter)->x;
     //x.SubtractVector(CenterOfGravity);
-    InertiaTensor[0] += ptr->type->mass*(x.x[1]*x.x[1] + x.x[2]*x.x[2]);
-    InertiaTensor[1] += ptr->type->mass*(-x.x[0]*x.x[1]);
-    InertiaTensor[2] += ptr->type->mass*(-x.x[0]*x.x[2]);
-    InertiaTensor[3] += ptr->type->mass*(-x.x[1]*x.x[0]);
-    InertiaTensor[4] += ptr->type->mass*(x.x[0]*x.x[0] + x.x[2]*x.x[2]);
-    InertiaTensor[5] += ptr->type->mass*(-x.x[1]*x.x[2]);
-    InertiaTensor[6] += ptr->type->mass*(-x.x[2]*x.x[0]);
-    InertiaTensor[7] += ptr->type->mass*(-x.x[2]*x.x[1]);
-    InertiaTensor[8] += ptr->type->mass*(x.x[0]*x.x[0] + x.x[1]*x.x[1]);
+    InertiaTensor[0] += (*iter)->type->mass*(x[1]*x[1] + x[2]*x[2]);
+    InertiaTensor[1] += (*iter)->type->mass*(-x[0]*x[1]);
+    InertiaTensor[2] += (*iter)->type->mass*(-x[0]*x[2]);
+    InertiaTensor[3] += (*iter)->type->mass*(-x[1]*x[0]);
+    InertiaTensor[4] += (*iter)->type->mass*(x[0]*x[0] + x[2]*x[2]);
+    InertiaTensor[5] += (*iter)->type->mass*(-x[1]*x[2]);
+    InertiaTensor[6] += (*iter)->type->mass*(-x[2]*x[0]);
+    InertiaTensor[7] += (*iter)->type->mass*(-x[2]*x[1]);
+    InertiaTensor[8] += (*iter)->type->mass*(x[0]*x[0] + x[1]*x[1]);
   }
   // print InertiaTensor for debugging
@@ -387 +387 @@
 
     // sum up inertia tensor
-    ptr = start;
-    while (ptr->next != end) {
-      ptr = ptr->next;
-      Vector x;
-      x.CopyVector(&ptr->x);
-      //x.SubtractVector(CenterOfGravity);
-      InertiaTensor[0] += ptr->type->mass*(x.x[1]*x.x[1] + x.x[2]*x.x[2]);
-      InertiaTensor[1] += ptr->type->mass*(-x.x[0]*x.x[1]);
-      InertiaTensor[2] += ptr->type->mass*(-x.x[0]*x.x[2]);
-      InertiaTensor[3] += ptr->type->mass*(-x.x[1]*x.x[0]);
-      InertiaTensor[4] += ptr->type->mass*(x.x[0]*x.x[0] + x.x[2]*x.x[2]);
-      InertiaTensor[5] += ptr->type->mass*(-x.x[1]*x.x[2]);
-      InertiaTensor[6] += ptr->type->mass*(-x.x[2]*x.x[0]);
-      InertiaTensor[7] += ptr->type->mass*(-x.x[2]*x.x[1]);
-      InertiaTensor[8] += ptr->type->mass*(x.x[0]*x.x[0] + x.x[1]*x.x[1]);
+    for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
+      Vector x = (*iter)->x;
+      InertiaTensor[0] += (*iter)->type->mass*(x[1]*x[1] + x[2]*x[2]);
+      InertiaTensor[1] += (*iter)->type->mass*(-x[0]*x[1]);
+      InertiaTensor[2] += (*iter)->type->mass*(-x[0]*x[2]);
+      InertiaTensor[3] += (*iter)->type->mass*(-x[1]*x[0]);
+      InertiaTensor[4] += (*iter)->type->mass*(x[0]*x[0] + x[2]*x[2]);
+      InertiaTensor[5] += (*iter)->type->mass*(-x[1]*x[2]);
+      InertiaTensor[6] += (*iter)->type->mass*(-x[2]*x[0]);
+      InertiaTensor[7] += (*iter)->type->mass*(-x[2]*x[1]);
+      InertiaTensor[8] += (*iter)->type->mass*(x[0]*x[0] + x[1]*x[1]);
     }
     // print InertiaTensor for debugging
@@ -425 +421 @@
 void molecule::Align(Vector *n)
 {
-  atom *ptr = start;
   double alpha, tmp;
   Vector z_axis;
-  z_axis.x[0] = 0.;
-  z_axis.x[1] = 0.;
-  z_axis.x[2] = 1.;
+  z_axis[0] = 0.;
+  z_axis[1] = 0.;
+  z_axis[2] = 1.;
 
   // rotate on z-x plane
   DoLog(0) && (Log() << Verbose(0) << "Begin of Aligning all atoms." << endl);
-  alpha = atan(-n->x[0]/n->x[2]);
+  alpha = atan(-n->at(0)/n->at(2));
   DoLog(1) && (Log() << Verbose(1) << "Z-X-angle: " << alpha << " ... ");
-  while (ptr->next != end) {
-    ptr = ptr->next;
-    tmp = ptr->x.x[0];
-    ptr->x.x[0] =  cos(alpha) * tmp + sin(alpha) * ptr->x.x[2];
-    ptr->x.x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x.x[2];
+  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];
     for (int j=0;j<MDSteps;j++) {
-      tmp = ptr->Trajectory.R.at(j).x[0];
-      ptr->Trajectory.R.at(j).x[0] =  cos(alpha) * tmp + sin(alpha) * ptr->Trajectory.R.at(j).x[2];
-      ptr->Trajectory.R.at(j).x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->Trajectory.R.at(j).x[2];
+      tmp = (*iter)->Trajectory.R.at(j)[0];
+      (*iter)->Trajectory.R.at(j)[0] =  cos(alpha) * tmp + sin(alpha) * (*iter)->Trajectory.R.at(j)[2];
+      (*iter)->Trajectory.R.at(j)[2] = -sin(alpha) * tmp + cos(alpha) * (*iter)->Trajectory.R.at(j)[2];
     }
   }
   // rotate n vector
-  tmp = n->x[0];
-  n->x[0] =  cos(alpha) * tmp +  sin(alpha) * n->x[2];
-  n->x[2] = -sin(alpha) * tmp +  cos(alpha) * n->x[2];
-  DoLog(1) && (Log() << Verbose(1) << "alignment vector after first rotation: ");
-  n->Output();
-  DoLog(0) && (Log() << Verbose(0) << endl);
+  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);
+  DoLog(1) && (Log() << Verbose(1) << "alignment vector after first rotation: " << n << endl);
 
   // rotate on z-y plane
-  ptr = start;
-  alpha = atan(-n->x[1]/n->x[2]);
+  alpha = atan(-n->at(1)/n->at(2));
   DoLog(1) && (Log() << Verbose(1) << "Z-Y-angle: " << alpha << " ... ");
-  while (ptr->next != end) {
-    ptr = ptr->next;
-    tmp = ptr->x.x[1];
-    ptr->x.x[1] =  cos(alpha) * tmp + sin(alpha) * ptr->x.x[2];
-    ptr->x.x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x.x[2];
+  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];
     for (int j=0;j<MDSteps;j++) {
-      tmp = ptr->Trajectory.R.at(j).x[1];
-      ptr->Trajectory.R.at(j).x[1] =  cos(alpha) * tmp + sin(alpha) * ptr->Trajectory.R.at(j).x[2];
-      ptr->Trajectory.R.at(j).x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->Trajectory.R.at(j).x[2];
+      tmp = (*iter)->Trajectory.R.at(j)[1];
+      (*iter)->Trajectory.R.at(j)[1] =  cos(alpha) * tmp + sin(alpha) * (*iter)->Trajectory.R.at(j)[2];
+      (*iter)->Trajectory.R.at(j)[2] = -sin(alpha) * tmp + cos(alpha) * (*iter)->Trajectory.R.at(j)[2];
     }
   }
   // rotate n vector (for consistency check)
-  tmp = n->x[1];
-  n->x[1] =  cos(alpha) * tmp +  sin(alpha) * n->x[2];
-  n->x[2] = -sin(alpha) * tmp +  cos(alpha) * n->x[2];
-
-  DoLog(1) && (Log() << Verbose(1) << "alignment vector after second rotation: ");
-  n->Output();
-  DoLog(1) && (Log() << Verbose(1) << endl);
+  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);
+
+
+  DoLog(1) && (Log() << Verbose(1) << "alignment vector after second rotation: " << n << endl);
   DoLog(0) && (Log() << Verbose(0) << "End of Aligning all atoms." << endl);
 };
@@ -492 +481 @@
   Vector a,b,c,d;
   struct lsq_params *par = (struct lsq_params *)params;
-  atom *ptr = par->mol->start;
 
   // initialize vectors
-  a.x[0] = gsl_vector_get(x,0);
-  a.x[1] = gsl_vector_get(x,1);
-  a.x[2] = gsl_vector_get(x,2);
-  b.x[0] = gsl_vector_get(x,3);
-  b.x[1] = gsl_vector_get(x,4);
-  b.x[2] = gsl_vector_get(x,5);
+  a[0] = gsl_vector_get(x,0);
+  a[1] = gsl_vector_get(x,1);
+  a[2] = gsl_vector_get(x,2);
+  b[0] = gsl_vector_get(x,3);
+  b[1] = gsl_vector_get(x,4);
+  b[2] = gsl_vector_get(x,5);
   // go through all atoms
-  while (ptr != par->mol->end) {
-    ptr = ptr->next;
-    if (ptr->type == ((struct lsq_params *)params)->type) { // for specific type
-      c.CopyVector(&ptr->x);  // copy vector to temporary one
-      c.SubtractVector(&a);   // subtract offset vector
-      t = c.ScalarProduct(&b);           // get direction parameter
-      d.CopyVector(&b);       // and create vector
-      d.Scale(&t);
-      c.SubtractVector(&d);   // ... yielding distance vector
-      res += d.ScalarProduct((const Vector *)&d);        // add squared distance
+  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;
+      t = c.ScalarProduct(b);           // get direction parameter
+      d = t*b;       // and create vector
+      c -= d;   // ... yielding distance vector
+      res += d.ScalarProduct(d);        // add squared distance
     }
   }