Changes in src/molecule.cpp [822f01:bf3817]

File:

• src/molecule.cpp (modified) (26 diffs)

src/molecule.cpp

@@ -5 +5 @@
  */
 
+// include config.h
 #ifdef HAVE_CONFIG_H
 #include <config.h>
@@ -24 +25 @@
 #include "element.hpp"
 #include "graph.hpp"
-#include "helpers.hpp"
+#include "Helpers/helpers.hpp"
 #include "leastsquaremin.hpp"
 #include "linkedcell.hpp"
 #include "lists.hpp"
-#include "log.hpp"
+#include "Helpers/Log.hpp"
 #include "molecule.hpp"
 
@@ -34 +35 @@
 #include "stackclass.hpp"
 #include "tesselation.hpp"
-#include "vector.hpp"
-#include "Matrix.hpp"
+#include "LinearAlgebra/Vector.hpp"
+#include "LinearAlgebra/Matrix.hpp"
 #include "World.hpp"
 #include "Box.hpp"
-#include "Plane.hpp"
+#include "LinearAlgebra/Plane.hpp"
 #include "Exceptions/LinearDependenceException.hpp"
 
@@ -157 +158 @@
   atomIds.erase( atom->getId() );
   atoms.remove( atom );
-  formula-=atom->type;
+  formula-=atom->getType();
   atom->removeFromMolecule();
   return iter;
@@ -169 +170 @@
     atomIds.erase( key->getId() );
     atoms.remove( key );
-    formula-=key->type;
+    formula-=key->getType();
     key->removeFromMolecule();
   }
@@ -191 +192 @@
   if (res.second) { // push atom if went well
     atoms.push_back(key);
-    formula+=key->type;
+    formula+=key->getType();
     return pair<iterator,bool>(molecule::iterator(--end()),res.second);
   } else {
@@ -212 +213 @@
   if (pointer != NULL) {
     pointer->sort = &pointer->nr;
-    if (pointer->type != NULL) {
-      formula += pointer->type;
-      if (pointer->type->Z != 1)
+    if (pointer->getType() != NULL) {
+      formula += pointer->getType();
+      if (pointer->getType()->Z != 1)
         NoNonHydrogen++;
       if(pointer->getName() == "Unknown"){
         stringstream sstr;
-        sstr << pointer->type->getSymbol() << pointer->nr+1;
+        sstr << pointer->getType()->getSymbol() << pointer->nr+1;
         pointer->setName(sstr.str());
       }
@@ -239 +240 @@
   if (pointer != NULL) {
     atom *walker = pointer->clone();
-    formula += walker->type;
+    formula += walker->getType();
     walker->setName(pointer->getName());
     walker->nr = last_atom++;  // increase number within molecule
     insert(walker);
-    if ((pointer->type != NULL) && (pointer->type->Z != 1))
+    if ((pointer->getType() != NULL) && (pointer->getType()->Z != 1))
       NoNonHydrogen++;
     retval=walker;
@@ -301 +302 @@
 //  Log() << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl;
   // create vector in direction of bond
-  InBondvector = TopReplacement->x - TopOrigin->x;
+  InBondvector = TopReplacement->getPosition() - TopOrigin->getPosition();
   bondlength = InBondvector.Norm();
 
@@ -313 +314 @@
     Orthovector1.Zero();
     for (int i=NDIM;i--;) {
-      l = TopReplacement->x[i] - TopOrigin->x[i];
+      l = TopReplacement->at(i) - TopOrigin->at(i);
       if (fabs(l) > BondDistance) { // is component greater than bond distance
         Orthovector1[i] = (l < 0) ? -1. : +1.;
@@ -328 +329 @@
   InBondvector.Normalize();
   // get typical bond length and store as scale factor for later
-  ASSERT(TopOrigin->type != NULL, "AddHydrogenReplacementAtom: element of TopOrigin is not given.");
-  BondRescale = TopOrigin->type->HBondDistance[TopBond->BondDegree-1];
+  ASSERT(TopOrigin->getType() != NULL, "AddHydrogenReplacementAtom: element of TopOrigin is not given.");
+  BondRescale = TopOrigin->getType()->HBondDistance[TopBond->BondDegree-1];
   if (BondRescale == -1) {
     DoeLog(1) && (eLog()<< Verbose(1) << "There is no typical hydrogen bond distance in replacing bond (" << TopOrigin->getName() << "<->" << TopReplacement->getName() << ") of degree " << TopBond->BondDegree << "!" << endl);
@@ -343 +344 @@
     case 1:
       FirstOtherAtom = World::getInstance().createAtom();    // new atom
-      FirstOtherAtom->type = elemente->FindElement(1);  // element is Hydrogen
-      FirstOtherAtom->v = TopReplacement->v; // copy velocity
+      FirstOtherAtom->setType(1);  // element is Hydrogen
+      FirstOtherAtom->AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity
       FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
-      if (TopReplacement->type->Z == 1) { // neither rescale nor replace if it's already hydrogen
+      if (TopReplacement->getType()->Z == 1) { // neither rescale nor replace if it's already hydrogen
         FirstOtherAtom->father = TopReplacement;
         BondRescale = bondlength;
@@ -353 +354 @@
       }
       InBondvector *= BondRescale;   // rescale the distance vector to Hydrogen bond length
-      FirstOtherAtom->x = TopOrigin->x; // set coordination to origin ...
-      FirstOtherAtom->x += InBondvector;  // ... and add distance vector to replacement atom
+      FirstOtherAtom->setPosition(TopOrigin->getPosition() + InBondvector); // set coordination to origin and add distance vector to replacement atom
       AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
 //      Log() << Verbose(4) << "Added " << *FirstOtherAtom << " at: ";
@@ -387 +387 @@
         // determine the plane of these two with the *origin
         try {
-          Orthovector1 =Plane(TopOrigin->x, FirstOtherAtom->x, SecondOtherAtom->x).getNormal();
+          Orthovector1 =Plane(TopOrigin->getPosition(), FirstOtherAtom->getPosition(), SecondOtherAtom->getPosition()).getNormal();
         }
         catch(LinearDependenceException &excp){
@@ -408 +408 @@
       FirstOtherAtom = World::getInstance().createAtom();
       SecondOtherAtom = World::getInstance().createAtom();
-      FirstOtherAtom->type = elemente->FindElement(1);
-      SecondOtherAtom->type = elemente->FindElement(1);
-      FirstOtherAtom->v = TopReplacement->v; // copy velocity
+      FirstOtherAtom->setType(1);
+      SecondOtherAtom->setType(1);
+      FirstOtherAtom->AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity
       FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
-      SecondOtherAtom->v = TopReplacement->v; // copy velocity
+      SecondOtherAtom->AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity
       SecondOtherAtom->FixedIon = TopReplacement->FixedIon;
       FirstOtherAtom->father = NULL;  // we are just an added hydrogen with no father
       SecondOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
-      bondangle = TopOrigin->type->HBondAngle[1];
+      bondangle = TopOrigin->getType()->HBondAngle[1];
       if (bondangle == -1) {
         DoeLog(1) && (eLog()<< Verbose(1) << "There is no typical hydrogen bond angle in replacing bond (" << TopOrigin->getName() << "<->" << TopReplacement->getName() << ") of degree " << TopBond->BondDegree << "!" << endl);
@@ -430 +430 @@
 //      Log() << Verbose(0) << endl;
 //      Log() << Verbose(3) << "Half the bond angle is " << bondangle << ", sin and cos of it: " << sin(bondangle) << ", " << cos(bondangle) << endl;
-      FirstOtherAtom->x.Zero();
-      SecondOtherAtom->x.Zero();
+      FirstOtherAtom->Zero();
+      SecondOtherAtom->Zero();
       for(int i=NDIM;i--;) { // rotate by half the bond angle in both directions (InBondvector is bondangle = 0 direction)
-        FirstOtherAtom->x[i] = InBondvector[i] * cos(bondangle) + Orthovector1[i] * (sin(bondangle));
-        SecondOtherAtom->x[i] = InBondvector[i] * cos(bondangle) + Orthovector1[i] * (-sin(bondangle));
-      }
-      FirstOtherAtom->x *= BondRescale;  // rescale by correct BondDistance
-      SecondOtherAtom->x *= BondRescale;
+        FirstOtherAtom->set(i, InBondvector[i] * cos(bondangle) + Orthovector1[i] * (sin(bondangle)));
+        SecondOtherAtom->set(i, InBondvector[i] * cos(bondangle) + Orthovector1[i] * (-sin(bondangle)));
+      }
+      FirstOtherAtom->Scale(BondRescale);  // rescale by correct BondDistance
+      SecondOtherAtom->Scale(BondRescale);
       //Log() << Verbose(3) << "ReScaleCheck: " << FirstOtherAtom->x.Norm() << " and " << SecondOtherAtom->x.Norm() << "." << endl;
-      for(int i=NDIM;i--;) { // and make relative to origin atom
-        FirstOtherAtom->x[i] += TopOrigin->x[i];
-        SecondOtherAtom->x[i] += TopOrigin->x[i];
-      }
+      *FirstOtherAtom += TopOrigin->getPosition();
+      *SecondOtherAtom += TopOrigin->getPosition();
       // ... and add to molecule
       AllWentWell = AllWentWell && AddAtom(FirstOtherAtom);
@@ -464 +462 @@
       SecondOtherAtom = World::getInstance().createAtom();
       ThirdOtherAtom = World::getInstance().createAtom();
-      FirstOtherAtom->type = elemente->FindElement(1);
-      SecondOtherAtom->type = elemente->FindElement(1);
-      ThirdOtherAtom->type = elemente->FindElement(1);
-      FirstOtherAtom->v = TopReplacement->v; // copy velocity
+      FirstOtherAtom->setType(1);
+      SecondOtherAtom->setType(1);
+      ThirdOtherAtom->setType(1);
+      FirstOtherAtom->AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity
       FirstOtherAtom->FixedIon = TopReplacement->FixedIon;
-      SecondOtherAtom->v = TopReplacement->v; // copy velocity
+      SecondOtherAtom->AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity
       SecondOtherAtom->FixedIon = TopReplacement->FixedIon;
-      ThirdOtherAtom->v = TopReplacement->v; // copy velocity
+      ThirdOtherAtom->AtomicVelocity = TopReplacement->AtomicVelocity; // copy velocity
       ThirdOtherAtom->FixedIon = TopReplacement->FixedIon;
       FirstOtherAtom->father = NULL;  //  we are just an added hydrogen with no father
@@ -494 +492 @@
 
       // create correct coordination for the three atoms
-      alpha = (TopOrigin->type->HBondAngle[2])/180.*M_PI/2.;  // retrieve triple bond angle from database
+      alpha = (TopOrigin->getType()->HBondAngle[2])/180.*M_PI/2.;  // retrieve triple bond angle from database
       l = BondRescale;        // desired bond length
       b = 2.*l*sin(alpha);    // base length of isosceles triangle
@@ -505 +503 @@
       factors[1] = f;
       factors[2] = 0.;
-      FirstOtherAtom->x.LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
+      FirstOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
       factors[1] = -0.5*f;
       factors[2] = g;
-      SecondOtherAtom->x.LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
+      SecondOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
       factors[2] = -g;
-      ThirdOtherAtom->x.LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
+      ThirdOtherAtom->LinearCombinationOfVectors(InBondvector, Orthovector1, Orthovector2, factors);
 
       // rescale each to correct BondDistance
@@ -518 +516 @@
 
       // and relative to *origin atom
-      FirstOtherAtom->x += TopOrigin->x;
-      SecondOtherAtom->x += TopOrigin->x;
-      ThirdOtherAtom->x += TopOrigin->x;
+      *FirstOtherAtom += TopOrigin->getPosition();
+      *SecondOtherAtom += TopOrigin->getPosition();
+      *ThirdOtherAtom += TopOrigin->getPosition();
 
       // ... and add to molecule
@@ -596 +594 @@
     *item >> x[1];
     *item >> x[2];
-    Walker->type = elemente->FindElement(shorthand);
-    if (Walker->type == NULL) {
+    Walker->setType(elemente->FindElement(shorthand));
+    if (Walker->getType() == NULL) {
       DoeLog(1) && (eLog()<< Verbose(1) << "Could not parse the element at line: '" << line << "', setting to H.");
-      Walker->type = elemente->FindElement(1);
+      Walker->setType(1);
     }
     if (Walker->Trajectory.R.size() <= (unsigned int)MDSteps) {
@@ -606 +604 @@
       Walker->Trajectory.F.resize(MDSteps+10);
     }
+    Walker->setPosition(Vector(x));
     for(j=NDIM;j--;) {
-      Walker->x[j] = x[j];
       Walker->Trajectory.R.at(MDSteps-1)[j] = x[j];
       Walker->Trajectory.U.at(MDSteps-1)[j] = 0;
@@ -702 +700 @@
   atom1->RegisterBond(Binder);
   atom2->RegisterBond(Binder);
-  if ((atom1->type != NULL) && (atom1->type->Z != 1) && (atom2->type != NULL) && (atom2->type->Z != 1))
+  if ((atom1->getType() != NULL) && (atom1->getType()->Z != 1) && (atom2->getType() != NULL) && (atom2->getType()->Z != 1))
     NoNonBonds++;
 
@@ -776 +774 @@
   ASSERT(pointer, "Null pointer passed to molecule::RemoveAtom().");
   OBSERVE;
-  formula-=pointer->type;
+  formula-=pointer->getType();
   RemoveBonds(pointer);
   erase(pointer);
@@ -790 +788 @@
   if (pointer == NULL)
     return false;
-  formula-=pointer->type;
+  formula-=pointer->getType();
   erase(pointer);
   return true;
@@ -972 +970 @@
   for (molecule::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
     (*iter)->nr = i;   // update number in molecule (for easier referencing in FragmentMolecule lateron)
-    if ((*iter)->type->Z != 1) // count non-hydrogen atoms whilst at it
+    if ((*iter)->getType()->Z != 1) // count non-hydrogen atoms whilst at it
       NoNonHydrogen++;
     stringstream sstr;
-    sstr << (*iter)->type->getSymbol() << (*iter)->nr+1;
+    sstr << (*iter)->getType()->getSymbol() << (*iter)->nr+1;
     (*iter)->setName(sstr.str());
     DoLog(3) && (Log() << Verbose(3) << "Naming atom nr. " << (*iter)->nr << " " << (*iter)->getName() << "." << endl);