Changes in src/molecule_geometry.cpp [68f03d:a67d19]

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::getInstance().getDomain();
+  double * const cell_size = World::get()->cell_size;
   double *M = ReturnFullMatrixforSymmetric(cell_size);
   double *Minv = InverseMatrix(M);
 
   // go through all atoms
-  ActOnAllVectors( &Vector::SubtractVector, *Center);
+  ActOnAllVectors( &Vector::SubtractVector, Center);
   ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
 
@@ -48 +48 @@
 {
   bool status = true;
-  double * const cell_size = World::getInstance().getDomain();
+  double * const cell_size = World::get()->cell_size;
   double *M = ReturnFullMatrixforSymmetric(cell_size);
   double *Minv = InverseMatrix(M);
@@ -72 +72 @@
   if (ptr != end) {   //list not empty?
     for (int i=NDIM;i--;) {
-      max->at(i) = ptr->x[i];
-      min->at(i) = ptr->x[i];
+      max->x[i] = ptr->x.x[i];
+      min->x[i] = ptr->x.x[i];
     }
     while (ptr->next != end) {  // continue with second if present
@@ -79 +79 @@
       //ptr->Output(1,1,out);
       for (int i=NDIM;i--;) {
-        max->at(i) = (max->at(i) < ptr->x[i]) ? ptr->x[i] : max->at(i);
-        min->at(i) = (min->at(i) > ptr->x[i]) ? ptr->x[i] : min->at(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];
       }
     }
@@ -89 +89 @@
 //    Log() << Verbose(0) << endl;
     min->Scale(-1.);
-    (*max) += (*min);
+    max->AddVector(min);
     Translate(min);
     Center.Zero();
@@ -112 +112 @@
       ptr = ptr->next;
       Num++;
-      Center += ptr->x;
+      Center.AddVector(&ptr->x);
     }
     Center.Scale(-1./Num); // divide through total number (and sign for direction)
@@ -136 +136 @@
       ptr = ptr->next;
       Num += 1.;
-      tmp = ptr->x;
-      (*a) += tmp;
+      tmp.CopyVector(&ptr->x);
+      a->AddVector(&tmp);
     }
     a->Scale(1./Num); // divide through total mass (and sign for direction)
@@ -161 +161 @@
       ptr = ptr->next;
       Num += ptr->type->mass;
-      tmp = ptr->type->mass * ptr->x;
-      (*a) += tmp;
+      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)
@@ -188 +189 @@
 void molecule::CenterAtVector(Vector *newcenter)
 {
-  Center = *newcenter;
+  Center.CopyVector(newcenter);
 };
 
@@ -194 +195 @@
 /** 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)
@@ -207 +203 @@
     ptr = ptr->next;
     for (int j=0;j<MDSteps;j++)
-      ptr->Trajectory.R.at(j).ScaleAll(*factor);
-    ptr->x.ScaleAll(*factor);
+      ptr->Trajectory.R.at(j).Scale(factor);
+    ptr->x.Scale(factor);
   }
 };
@@ -222 +218 @@
     ptr = ptr->next;
     for (int j=0;j<MDSteps;j++)
-      ptr->Trajectory.R.at(j) += (*trans);
-    ptr->x += (*trans);
+      ptr->Trajectory.R.at(j).Translate(trans);
+    ptr->x.Translate(trans);
   }
 };
@@ -233 +229 @@
 void molecule::TranslatePeriodically(const Vector *trans)
 {
-  double * const cell_size = World::getInstance().getDomain();
+  double * const cell_size = World::get()->cell_size;
   double *M = ReturnFullMatrixforSymmetric(cell_size);
   double *Minv = InverseMatrix(M);
 
   // go through all atoms
-  ActOnAllVectors( &Vector::SubtractVector, *trans);
+  ActOnAllVectors( &Vector::SubtractVector, trans);
   ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
 
@@ -251 +247 @@
 void molecule::Mirror(const Vector *n)
 {
-  ActOnAllVectors( &Vector::Mirror, *n );
+  ActOnAllVectors( &Vector::Mirror, n );
 };
 
@@ -260 +256 @@
 {
   atom *Walker = start;
-  double * const cell_size = World::getInstance().getDomain();
+  double * const cell_size = World::get()->cell_size;
   double *matrix = ReturnFullMatrixforSymmetric(cell_size);
   double *inversematrix = InverseMatrix(cell_size);
@@ -275 +271 @@
       if (Walker->type->Z != 1) {
 #endif
-        Testvector = Walker->x;
+        Testvector.CopyVector(&Walker->x);
         Testvector.MatrixMultiplication(inversematrix);
         Translationvector.Zero();
@@ -281 +277 @@
          if (Walker->nr < (*Runner)->GetOtherAtom(Walker)->nr) // otherwise we shift one to, the other fro and gain nothing
             for (int j=0;j<NDIM;j++) {
-              tmp = Walker->x[j] - (*Runner)->GetOtherAtom(Walker)->x[j];
+              tmp = Walker->x.x[j] - (*Runner)->GetOtherAtom(Walker)->x.x[j];
               if ((fabs(tmp)) > BondDistance) {
                 flag = false;
-                DoLog(0) && (Log() << Verbose(0) << "Hit: atom " << Walker->getName() << " in bond " << *(*Runner) << " has to be shifted due to " << tmp << "." << endl);
+                DoLog(0) && (Log() << Verbose(0) << "Hit: atom " << Walker->Name << " in bond " << *(*Runner) << " has to be shifted due to " << tmp << "." << endl);
                 if (tmp > 0)
-                  Translationvector[j] -= 1.;
+                  Translationvector.x[j] -= 1.;
                 else
-                  Translationvector[j] += 1.;
+                  Translationvector.x[j] += 1.;
               }
             }
         }
-        Testvector += Translationvector;
+        Testvector.AddVector(&Translationvector);
         Testvector.MatrixMultiplication(matrix);
-        Center += Testvector;
-        Log() << Verbose(1) << "vector is: " << Testvector << endl;
+        Center.AddVector(&Testvector);
+        DoLog(1) && (Log() << Verbose(1) << "vector is: ");
+        Testvector.Output();
+        DoLog(0) && (Log() << Verbose(0) << 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 = (*Runner)->GetOtherAtom(Walker)->x;
+            Testvector.CopyVector(&(*Runner)->GetOtherAtom(Walker)->x);
             Testvector.MatrixMultiplication(inversematrix);
-            Testvector += Translationvector;
+            Testvector.AddVector(&Translationvector);
             Testvector.MatrixMultiplication(matrix);
-            Center += Testvector;
-            Log() << Verbose(1) << "Hydrogen vector is: " << Testvector << endl;
+            Center.AddVector(&Testvector);
+            DoLog(1) && (Log() << Verbose(1) << "Hydrogen vector is: ");
+            Testvector.Output();
+            DoLog(0) && (Log() << Verbose(0) << endl);
           }
         }
@@ -338 +338 @@
   while (ptr->next != end) {
     ptr = ptr->next;
-    Vector x = ptr->x;
+    Vector x;
+    x.CopyVector(&ptr->x);
     //x.SubtractVector(CenterOfGravity);
-    InertiaTensor[0] += ptr->type->mass*(x[1]*x[1] + x[2]*x[2]);
-    InertiaTensor[1] += ptr->type->mass*(-x[0]*x[1]);
-    InertiaTensor[2] += ptr->type->mass*(-x[0]*x[2]);
-    InertiaTensor[3] += ptr->type->mass*(-x[1]*x[0]);
-    InertiaTensor[4] += ptr->type->mass*(x[0]*x[0] + x[2]*x[2]);
-    InertiaTensor[5] += ptr->type->mass*(-x[1]*x[2]);
-    InertiaTensor[6] += ptr->type->mass*(-x[2]*x[0]);
-    InertiaTensor[7] += ptr->type->mass*(-x[2]*x[1]);
-    InertiaTensor[8] += ptr->type->mass*(x[0]*x[0] + x[1]*x[1]);
+    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]);
   }
   // print InertiaTensor for debugging
@@ -389 +390 @@
     while (ptr->next != end) {
       ptr = ptr->next;
-      Vector x = ptr->x;
+      Vector x;
+      x.CopyVector(&ptr->x);
       //x.SubtractVector(CenterOfGravity);
-      InertiaTensor[0] += ptr->type->mass*(x[1]*x[1] + x[2]*x[2]);
-      InertiaTensor[1] += ptr->type->mass*(-x[0]*x[1]);
-      InertiaTensor[2] += ptr->type->mass*(-x[0]*x[2]);
-      InertiaTensor[3] += ptr->type->mass*(-x[1]*x[0]);
-      InertiaTensor[4] += ptr->type->mass*(x[0]*x[0] + x[2]*x[2]);
-      InertiaTensor[5] += ptr->type->mass*(-x[1]*x[2]);
-      InertiaTensor[6] += ptr->type->mass*(-x[2]*x[0]);
-      InertiaTensor[7] += ptr->type->mass*(-x[2]*x[1]);
-      InertiaTensor[8] += ptr->type->mass*(x[0]*x[0] + x[1]*x[1]);
+      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]);
     }
     // print InertiaTensor for debugging
@@ -426 +428 @@
   double alpha, tmp;
   Vector z_axis;
-  z_axis[0] = 0.;
-  z_axis[1] = 0.;
-  z_axis[2] = 1.;
+  z_axis.x[0] = 0.;
+  z_axis.x[1] = 0.;
+  z_axis.x[2] = 1.;
 
   // rotate on z-x plane
   DoLog(0) && (Log() << Verbose(0) << "Begin of Aligning all atoms." << endl);
-  alpha = atan(-n->at(0)/n->at(2));
+  alpha = atan(-n->x[0]/n->x[2]);
   DoLog(1) && (Log() << Verbose(1) << "Z-X-angle: " << alpha << " ... ");
   while (ptr->next != end) {
     ptr = ptr->next;
-    tmp = ptr->x[0];
-    ptr->x[0] =  cos(alpha) * tmp + sin(alpha) * ptr->x[2];
-    ptr->x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x[2];
+    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 (int j=0;j<MDSteps;j++) {
-      tmp = ptr->Trajectory.R.at(j)[0];
-      ptr->Trajectory.R.at(j)[0] =  cos(alpha) * tmp + sin(alpha) * ptr->Trajectory.R.at(j)[2];
-      ptr->Trajectory.R.at(j)[2] = -sin(alpha) * tmp + cos(alpha) * ptr->Trajectory.R.at(j)[2];
+      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];
     }
   }
   // 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);
-  DoLog(1) && (Log() << Verbose(1) << "alignment vector after first rotation: " << n << endl);
+  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);
 
   // rotate on z-y plane
   ptr = start;
-  alpha = atan(-n->at(1)/n->at(2));
+  alpha = atan(-n->x[1]/n->x[2]);
   DoLog(1) && (Log() << Verbose(1) << "Z-Y-angle: " << alpha << " ... ");
   while (ptr->next != end) {
     ptr = ptr->next;
-    tmp = ptr->x[1];
-    ptr->x[1] =  cos(alpha) * tmp + sin(alpha) * ptr->x[2];
-    ptr->x[2] = -sin(alpha) * tmp + cos(alpha) * ptr->x[2];
+    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 (int j=0;j<MDSteps;j++) {
-      tmp = ptr->Trajectory.R.at(j)[1];
-      ptr->Trajectory.R.at(j)[1] =  cos(alpha) * tmp + sin(alpha) * ptr->Trajectory.R.at(j)[2];
-      ptr->Trajectory.R.at(j)[2] = -sin(alpha) * tmp + cos(alpha) * ptr->Trajectory.R.at(j)[2];
+      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];
     }
   }
   // 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);
-
-
-  DoLog(1) && (Log() << Verbose(1) << "alignment vector after second rotation: " << n << endl);
+  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);
   DoLog(0) && (Log() << Verbose(0) << "End of Aligning all atoms." << endl);
 };
@@ -490 +495 @@
 
   // initialize vectors
-  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);
+  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);
   // go through all atoms
   while (ptr != par->mol->end) {
     ptr = ptr->next;
     if (ptr->type == ((struct lsq_params *)params)->type) { // for specific type
-      c = ptr->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
+      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
     }
   }