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

File:

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

src/molecule_geometry.cpp

@@ +1 @@
+/*
+ * Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2010 University of Bonn. All rights reserved.
+ * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
+ */
+
 /*
  * molecule_geometry.cpp
@@ -6 +13 @@
  */
 
+// include config.h
 #ifdef HAVE_CONFIG_H
 #include <config.h>
 #endif
 
+#include "Helpers/MemDebug.hpp"
+
 #include "Helpers/helpers.hpp"
 #include "Helpers/Log.hpp"
-#include "Helpers/MemDebug.hpp"
 #include "Helpers/Verbose.hpp"
 #include "LinearAlgebra/Line.hpp"
@@ -48 +57 @@
 
   // go through all atoms
-  ActOnAllVectors( &Vector::SubtractVector, *Center);
-  ActOnAllVectors( &Vector::SubtractVector, *CenterBox);
+  BOOST_FOREACH(atom* iter, atoms){
+    *iter -= *Center;
+    *iter -= *CenterBox;
+  }
   atoms.transformNodes(boost::bind(&Box::WrapPeriodically,domain,_1));
 
@@ -66 +77 @@
   Box &domain = World::getInstance().getDomain();
 
+  // go through all atoms
   atoms.transformNodes(boost::bind(&Box::WrapPeriodically,domain,_1));
 
@@ -83 +95 @@
   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);
       }
     }
@@ -101 +113 @@
     (*max) += (*min);
     Translate(min);
-    Center.Zero();
   }
   delete(min);
@@ -115 +126 @@
   int Num = 0;
   molecule::const_iterator iter = begin();  // start at first in list
+  Vector Center;
 
   Center.Zero();
-
   if (iter != end()) {   //list not empty?
     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)
     Translate(&Center);
-    Center.Zero();
   }
 };
@@ -143 +153 @@
     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)
@@ -176 +186 @@
   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()->getMass();
+      tmp = (*iter)->getType()->getMass() * (*iter)->getPosition();
       (*a) += tmp;
     }
@@ -194 +204 @@
 void molecule::CenterPeriodic()
 {
-  DeterminePeriodicCenter(Center);
+  Vector NewCenter;
+  DeterminePeriodicCenter(NewCenter);
+  // go through all atoms
+  BOOST_FOREACH(atom* iter, atoms){
+    *iter -= NewCenter;
+  }
 };
 
@@ -204 +219 @@
 void molecule::CenterAtVector(Vector *newcenter)
 {
-  Center = *newcenter;
+  // go through all atoms
+  BOOST_FOREACH(atom* iter, atoms){
+    *iter -= *newcenter;
+  }
 };
 
@@ -221 +239 @@
     for (int j=0;j<MDSteps;j++)
       (*iter)->Trajectory.R.at(j).ScaleAll(*factor);
-    (*iter)->x.ScaleAll(*factor);
+    (*iter)->ScaleAll(*factor);
   }
 };
@@ -233 +251 @@
     for (int j=0;j<MDSteps;j++)
       (*iter)->Trajectory.R.at(j) += (*trans);
-    (*iter)->x += (*trans);
+    *(*iter) += (*trans);
   }
 };
@@ -246 +264 @@
 
   // go through all atoms
-  ActOnAllVectors( &Vector::AddVector, *trans);
+  BOOST_FOREACH(atom* iter, atoms){
+    *iter += *trans;
+  }
   atoms.transformNodes(boost::bind(&Box::WrapPeriodically,domain,_1));
 
@@ -272 +292 @@
   bool flag;
   Vector Testvector, Translationvector;
+  Vector Center;
 
   do {
@@ -278 +299 @@
     for (molecule::const_iterator iter = begin(); iter != end(); ++iter) {
 #ifdef ADDHYDROGEN
-      if ((*iter)->type->Z != 1) {
+      if ((*iter)->getType()->getAtomicNumber() != 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 +324 @@
         // 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()->getAtomicNumber() == 1) {
+            Testvector = inversematrix * (*Runner)->GetOtherAtom((*iter))->getPosition();
             Testvector += Translationvector;
             Testvector *= matrix;
@@ -317 +338 @@
 
   Center.Scale(1./static_cast<double>(getAtomCount()));
+  CenterAtVector(&Center);
 };
 
@@ -335 +357 @@
   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 +376 @@
   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 +416 @@
   // 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