Changes in src/molecule.cpp [906822:aafd77]

File:

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

src/molecule.cpp

@@ -5 +5 @@
  */
 
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
 #include "Helpers/MemDebug.hpp"
 
 #include <cstring>
 #include <boost/bind.hpp>
+#include <boost/foreach.hpp>
+
+#include <gsl/gsl_inline.h>
+#include <gsl/gsl_heapsort.h>
 
 #include "World.hpp"
@@ -22 +30 @@
 #include "log.hpp"
 #include "molecule.hpp"
-#include "memoryallocator.hpp"
+
 #include "periodentafel.hpp"
 #include "stackclass.hpp"
 #include "tesselation.hpp"
 #include "vector.hpp"
+#include "Matrix.hpp"
 #include "World.hpp"
+#include "Box.hpp"
 #include "Plane.hpp"
 #include "Exceptions/LinearDependenceException.hpp"
@@ -79 +89 @@
 void molecule::setName(const std::string _name){
   OBSERVE;
+  cout << "Set name of molecule " << getId() << " to " << _name << endl;
   strncpy(name,_name.c_str(),MAXSTRINGSIZE);
 }
@@ -283 +294 @@
   Vector Orthovector1, Orthovector2;  // temporary vectors in coordination construction
   Vector InBondvector;    // vector in direction of *Bond
-  double *matrix = NULL;
+  const Matrix &matrix =  World::getInstance().getDomain().getM();
   bond *Binder = NULL;
-  double * const cell_size = World::getInstance().getDomain();
 
 //  Log() << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl;
@@ -306 +316 @@
       } // (signs are correct, was tested!)
     }
-    matrix = ReturnFullMatrixforSymmetric(cell_size);
-    Orthovector1.MatrixMultiplication(matrix);
+    Orthovector1 *= matrix;
     InBondvector -= Orthovector1; // subtract just the additional translation
-    delete[](matrix);
     bondlength = InBondvector.Norm();
 //    Log() << Verbose(4) << "Corrected InBondvector is now: ";
@@ -540 +548 @@
       break;
   }
-  delete[](matrix);
 
 //  Log() << Verbose(3) << "End of AddHydrogenReplacementAtom." << endl;
@@ -659 +666 @@
  * @param three vectors forming the matrix that defines the shape of the parallelpiped
  */
-molecule* molecule::CopyMoleculeFromSubRegion(const Vector offset, const double *parallelepiped) const {
+molecule* molecule::CopyMoleculeFromSubRegion(const Shape &region) const {
   molecule *copy = World::getInstance().createMolecule();
 
-  ActOnCopyWithEachAtomIfTrue ( &molecule::AddCopyAtom, copy, &atom::IsInParallelepiped, offset, parallelepiped );
+  BOOST_FOREACH(atom *iter,atoms){
+    if(iter->IsInShape(region)){
+      copy->AddCopyAtom(iter);
+    }
+  }
 
   //TODO: copy->BuildInducedSubgraph(this);
@@ -739 +750 @@
   else
     length = strlen(molname) - strlen(endname);
+  cout << "Set name of molecule " << getId() << " to " << molname << endl;
   strncpy(name, molname, length);
   name[length]='\0';
@@ -748 +760 @@
 void molecule::SetBoxDimension(Vector *dim)
 {
-  double * const cell_size = World::getInstance().getDomain();
-  cell_size[0] = dim->at(0);
-  cell_size[1] = 0.;
-  cell_size[2] = dim->at(1);
-  cell_size[3] = 0.;
-  cell_size[4] = 0.;
-  cell_size[5] = dim->at(2);
+  Matrix domain;
+  for(int i =0; i<NDIM;++i)
+    domain.at(i,i) = dim->at(i);
+  World::getInstance().setDomain(domain);
 };
 
@@ -847 +856 @@
 bool molecule::CheckBounds(const Vector *x) const
 {
-  double * const cell_size = World::getInstance().getDomain();
+  const Matrix &domain = World::getInstance().getDomain().getM();
   bool result = true;
-  int j =-1;
   for (int i=0;i<NDIM;i++) {
-    j += i+1;
-    result = result && ((x->at(i) >= 0) && (x->at(i) < cell_size[j]));
+    result = result && ((x->at(i) >= 0) && (x->at(i) < domain.at(i,i)));
   }
   //return result;
@@ -880 +887 @@
         ElementNo[i] = current++;
     }
-    ActOnAllAtoms( &atom::OutputArrayIndexed, output, (const int *)ElementNo, (int *)AtomNo, (const char *) NULL );
+    ActOnAllAtoms( &atom::OutputArrayIndexed, (ostream * const) output, (const int *)ElementNo, (int *)AtomNo, (const char *) NULL );
     return true;
   }
@@ -1003 +1010 @@
   for(int i=MAX_ELEMENTS;i--;)
     ElementCount += (ElementsInMolecule[i] != 0 ? 1 : 0);
-};
-
-
-/** Counts necessary number of valence electrons and returns number and SpinType.
- * \param configuration containing everything
- */
-void molecule::CalculateOrbitals(class config &configuration)
-{
-  configuration.MaxPsiDouble = configuration.PsiMaxNoDown = configuration.PsiMaxNoUp = configuration.PsiType = 0;
-  for(int i=MAX_ELEMENTS;i--;) {
-    if (ElementsInMolecule[i] != 0) {
-      //Log() << Verbose(0) << "CalculateOrbitals: " << elemente->FindElement(i)->name << " has a valence of " << (int)elemente->FindElement(i)->Valence << " and there are " << ElementsInMolecule[i] << " of it." << endl;
-      configuration.MaxPsiDouble += ElementsInMolecule[i]*((int)elemente->FindElement(i)->Valence);
-    }
-  }
-  configuration.PsiMaxNoDown = configuration.MaxPsiDouble/2 + (configuration.MaxPsiDouble % 2);
-  configuration.PsiMaxNoUp = configuration.MaxPsiDouble/2;
-  configuration.MaxPsiDouble /= 2;
-  configuration.PsiType = (configuration.PsiMaxNoDown == configuration.PsiMaxNoUp) ? 0 : 1;
-  if ((configuration.PsiType == 1) && (configuration.ProcPEPsi < 2) && ((configuration.PsiMaxNoDown != 1) || (configuration.PsiMaxNoUp != 0))) {
-    configuration.ProcPEGamma /= 2;
-    configuration.ProcPEPsi *= 2;
-  } else {
-    configuration.ProcPEGamma *= configuration.ProcPEPsi;
-    configuration.ProcPEPsi = 1;
-  }
-  cout << configuration.PsiMaxNoDown << ">" << configuration.PsiMaxNoUp << endl;
-  if (configuration.PsiMaxNoDown > configuration.PsiMaxNoUp) {
-    configuration.InitMaxMinStopStep = configuration.MaxMinStopStep = configuration.PsiMaxNoDown;
-    cout << configuration.PsiMaxNoDown << " " << configuration.InitMaxMinStopStep << endl;
-  } else {
-    configuration.InitMaxMinStopStep = configuration.MaxMinStopStep = configuration.PsiMaxNoUp;
-    cout << configuration.PsiMaxNoUp << " " << configuration.InitMaxMinStopStep << endl;
-  }
 };