Diff [ba9f5bcda4dd421dcb79d4ccff0dc34cfc072a6e:be97a841b19bbe05e88f3d3ac69c9cb5cd562f10] for / – MoleCuilder

src/Actions/ActionRegistry.hpp

rba9f5b	rbe97a8
9	9	#define ACTIONREGISTRY_HPP_
10	10
11		#include <ios~~tream~~>
	11	#include <iosfwd>
12	12	#include <string>
13	13	#include <map>

src/Actions/AnalysisAction/PairCorrelationAction.cpp

-              rba9f5b
+              rbe97a8
 #include "boundary.hpp"
 #include "linkedcell.hpp"
+#include "verbose.hpp"
 #include "log.hpp"
 #include "element.hpp"
 …
   dialog = UIFactory::getInstance().makeDialog();
   if (type == "P")
     dialog->queryVector("position", &Point, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription("position"));
+    dialog->queryVector("position", &Point, false, MapOfActions::getInstance().getDescription("position"));
   if (type == "S")
     dialog->queryMolecule("molecule-by-id", &Boundary, MapOfActions::getInstance().getDescription("molecule-by-id"));

src/Actions/AtomAction/AddAction.cpp

rba9f5b	rbe97a8
41	41
42	42	dialog->queryElement(NAME, &elements, MapOfActions::getInstance().getDescription(NAME));
43		dialog->queryVector("position", &position, ~~World::getInstance().getDomain(),~~ true, MapOfActions::getInstance().getDescription("position"));
	43	dialog->queryVector("position", &position, true, MapOfActions::getInstance().getDescription("position"));
44	44	cout << "pre-dialog" << endl;
45	45

src/Actions/AtomsCalculation_impl.hpp

-              rba9f5b
+              rbe97a8
 #include "Actions/AtomsCalculation.hpp"
 #include "Actions/Calculation_impl.hpp"
-#include <iostream>
 using namespace std;
 …
   Process::setMaxSteps(steps);
   Process::start();
+  for(World::AtomIterator iter=world->getAtomIter(descr);iter!=world->atomEnd();++iter){
+  for(World::internal_AtomIterator
+      iter=world->getAtomIter_internal(descr);
+      iter!=world->atomEnd_internal();
+      ++iter){
     Process::setCurrStep(iter.getCount());
     res->push_back(op(*iter));

src/Actions/ManipulateAtomsProcess.cpp

-              rba9f5b
+              rbe97a8
   setMaxSteps(world->numAtoms());
   start();
+  for(World::AtomIterator iter=world->getAtomIter(descr);iter!=world->atomEnd();++iter){
+  for(World::internal_AtomIterator
+      iter=world->getAtomIter_internal(descr);
+      iter!=world->atomEnd_internal();
+      ++iter){
     setCurrStep(iter.getCount());
     operation(*iter);

src/Actions/MapOfActions.cpp

-              rba9f5b
+              rbe97a8
 #include <boost/optional.hpp>
 #include <boost/program_options.hpp>
+#include <iostream>
 #include "CommandLineParser.hpp"
 …
   TypeMap["MaxDistance"] = Double;
   TypeMap["molecule-by-id"] = Molecule;
   TypeMap["molecule-by-name"] = Molecule;
+  TypeMap["molecule-by-name"] = String;
   TypeMap["nonconvex-file"] = String;
   TypeMap["order"] = Integer;

src/Actions/MoleculeAction/BondFileAction.cpp

rba9f5b	rbe97a8
24	24	#include "config.hpp"
25	25	#include "defs.hpp"
	26	#include "verbose.hpp"
26	27	#include "log.hpp"
27	28	#include "molecule.hpp"

src/Actions/MoleculeAction/FillWithMoleculeAction.cpp

-              rba9f5b
+              rbe97a8
   dialog->queryString(NAME, &filename, MapOfActions::getInstance().getDescription(NAME));
   dialog->queryVector("distances", &distances, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription("distances"));
   dialog->queryVector("lengths", &lengths, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription("lengths"));
+  dialog->queryVector("distances", &distances, false, MapOfActions::getInstance().getDescription("distances"));
+  dialog->queryVector("lengths", &lengths, false, MapOfActions::getInstance().getDescription("lengths"));
   dialog->queryBoolean("DoRotate", &DoRotate, MapOfActions::getInstance().getDescription("DoRotate"));
   dialog->queryDouble("MaxDistance", &MaxDistance, MapOfActions::getInstance().getDescription("MaxDistance"));

src/Actions/MoleculeAction/SuspendInWaterAction.cpp

rba9f5b	rbe97a8
22	22	#include "boundary.hpp"
23	23	#include "config.hpp"
	24	#include "verbose.hpp"
24	25	#include "log.hpp"
25	26	#include "config.hpp"

src/Actions/MoleculeAction/TranslateAction.cpp

rba9f5b	rbe97a8
55	55	bool periodic = false;
56	56
57		dialog->queryVector(NAME, &x, ~~World::getInstance().getDomain(),~~ false, MapOfActions::getInstance().getDescription(NAME));
	57	dialog->queryVector(NAME, &x, false, MapOfActions::getInstance().getDescription(NAME));
58	58	dialog->queryMolecule("molecule-by-id", &mol, MapOfActions::getInstance().getDescription("molecule-by-id"));
59	59	dialog->queryBoolean("periodic", &periodic, MapOfActions::getInstance().getDescription("periodic"));

src/Actions/WorldAction/AddEmptyBoundaryAction.cpp

-              rba9f5b
+              rbe97a8
   int j=0;
   dialog->queryVector(NAME, &boundary, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription(NAME));
+  dialog->queryVector(NAME, &boundary, false, MapOfActions::getInstance().getDescription(NAME));
   if(dialog->display()) {
 …
+    }
     // set new box size
     double * const cell_size = World::getInstance().getDomain();
+    double * const cell_size = new double[6];
     for (j=0;j<6;j++)
       cell_size[j] = 0.;
 …
       cell_size[j] = (Max[i]-Min[i]+2.*boundary[i]);
+    }
+    World::getInstance().setDomain(cell_size);
+    delete[] cell_size;
     // translate all atoms, such that Min is aty (0,0,0)
     AtomRunner = AllAtoms.begin();

src/Actions/WorldAction/CenterInBoxAction.cpp

-              rba9f5b
+              rbe97a8
   Dialog *dialog = UIFactory::getInstance().makeDialog();
   double * cell_size = World::getInstance().getDomain();
+  Box& cell_size = World::getInstance().getDomain();
   dialog->queryBox(NAME, &cell_size, MapOfActions::getInstance().getDescription(NAME));
   if(dialog->display()) {
-    World::getInstance().setDomain(cell_size);
     // center
     vector<molecule *> AllMolecules = World::getInstance().getAllMolecules();

src/Actions/WorldAction/CenterOnEdgeAction.cpp

-              rba9f5b
+              rbe97a8
 #include "vector.hpp"
 #include "World.hpp"
+#include "Matrix.hpp"
 #include <iostream>
 …
   Vector Min;
   Vector Max;
-  int j=0;
   dialog->queryEmpty(NAME, MapOfActions::getInstance().getDescription(NAME));
 …
+    }
     // set new box size
+    double * const cell_size = World::getInstance().getDomain();
+    for (j=0;j<6;j++)
+      cell_size[j] = 0.;
+    j=-1;
+    Matrix domain;
     for (int i=0;i<NDIM;i++) {
+      j += i+1;
+      cell_size[j] = (Max[i]-Min[i]);
+      double tmp = Max[i]-Min[i];
+      tmp = fabs(tmp)>=1. ? tmp : 1.0;
+      domain.at(i,i) = tmp;
+    }
     World::getInstance().setDomain(cell_size);
+    World::getInstance().setDomain(domain);
     // translate all atoms, such that Min is aty (0,0,0)
     for (vector<atom*>::iterator AtomRunner = AllAtoms.begin(); AtomRunner != AllAtoms.end(); ++AtomRunner)

src/Actions/WorldAction/ChangeBoxAction.cpp

-              rba9f5b
+              rbe97a8
 #include "verbose.hpp"
 #include "World.hpp"
+#include "Box.hpp"
+#include "Matrix.hpp"
 #include <iostream>
 …
   Dialog *dialog = UIFactory::getInstance().makeDialog();
   double * cell_size = World::getInstance().getDomain();
+  Box& cell_size = World::getInstance().getDomain();
   dialog->queryBox(NAME, &cell_size, MapOfActions::getInstance().getDescription(NAME));
   if(dialog->display()) {
     DoLog(0) && (Log() << Verbose(0) << "Setting box domain to " << cell_size[0] << "," << cell_size[1] << "," << cell_size[2] << "," << cell_size[3] << "," << cell_size[4] << "," << cell_size[5] << "," << endl);
     World::getInstance().setDomain(cell_size);
+    DoLog(0) && (Log() << Verbose(0) << "Setting box domain to " << cell_size.getM() << endl);
+    World::getInstance().setDomain(cell_size.getM());
     delete dialog;
     return Action::success;

src/Actions/WorldAction/RemoveSphereOfAtomsAction.cpp

rba9f5b	rbe97a8
41	41
42	42	dialog->queryDouble(NAME, &radius, MapOfActions::getInstance().getDescription(NAME));
43		dialog->queryVector("position", &point, ~~World::getInstance().getDomain(),~~ false, MapOfActions::getInstance().getDescription("position"));
	43	dialog->queryVector("position", &point, false, MapOfActions::getInstance().getDescription("position"));
44	44
45	45	if(dialog->display()) {

src/Actions/WorldAction/RepeatBoxAction.cpp

-              rba9f5b
+              rbe97a8
 #include "molecule.hpp"
 #include "vector.hpp"
+#include "Matrix.hpp"
 #include "verbose.hpp"
 #include "World.hpp"
+#include "Box.hpp"
 #include <iostream>
 …
   MoleculeListClass *molecules = World::getInstance().getMolecules();
   dialog->queryVector(NAME, &Repeater, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription(NAME));
+  dialog->queryVector(NAME, &Repeater, false, MapOfActions::getInstance().getDescription(NAME));
   //dialog->queryMolecule("molecule-by-id", &mol,MapOfActions::getInstance().getDescription("molecule-by-id"));
   vector<molecule *> AllMolecules;
 …
   if(dialog->display()) {
     (cout << "Repeating box " << Repeater << " times for (x,y,z) axis." << endl);
     double * const cell_size = World::getInstance().getDomain();
     double *M = ReturnFullMatrixforSymmetric(cell_size);
+    Matrix M = World::getInstance().getDomain().getM();
+    Matrix newM = M;
     Vector x,y;
     int n[NDIM];
+    Matrix repMat;
     for (int axis = 0; axis < NDIM; axis++) {
       Repeater[axis] = floor(Repeater[axis]);
 …
         Repeater[axis] = 1;
+      }
       cell_size[ ((axis == 1) ? 2 : (axis == 2) ? 5 : 0) ] *= Repeater[axis];
+      repMat.at(axis,axis) = Repeater[axis];
+    }
+    newM *= repMat;
+    World::getInstance().setDomain(newM);
-    World::getInstance().setDomain(cell_size);
     molecule *newmol = NULL;
     Vector ** vectors = NULL;
 …
                 DoeLog(1) && (eLog()<< Verbose(1) << "AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl);
               x = y;
               x.MatrixMultiplication(M);
+              x *= M;
               newmol = World::getInstance().createMolecule();
               molecules->insert(newmol);
 …
+      }
+    }
-    delete(M);
     delete dialog;
     return Action::success;

src/Actions/WorldAction/ScaleBoxAction.cpp

-              rba9f5b
+              rbe97a8
 #include "verbose.hpp"
 #include "World.hpp"
+#include "Box.hpp"
+#include "Matrix.hpp"
 #include <iostream>
 …
   Vector Scaler;
   double x[NDIM];
-  int j=0;
   dialog->queryVector(NAME, &Scaler, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription(NAME));
+  dialog->queryVector(NAME, &Scaler, false, MapOfActions::getInstance().getDescription(NAME));
   if(dialog->display()) {
 …
       (*AtomRunner)->x.ScaleAll(x);
+    }
+    j = -1;
+    double * const cell_size = World::getInstance().getDomain();
+    Matrix M = World::getInstance().getDomain().getM();
+    Matrix scale;
     for (int i=0;i<NDIM;i++) {
+      j += i+1;
+      cell_size[j]*=x[i];
+      scale.at(i,i) = x[i];
+    }
+    World::getInstance().setDomain(cell_size);
+    M *= scale;
+    World::getInstance().setDomain(M);
     delete dialog;

src/ConfigFileBuffer.cpp

rba9f5b	rbe97a8
9	9	#include "helpers.hpp"
10	10	#include "lists.hpp"
	11	#include "verbose.hpp"
11	12	#include "log.hpp"
12	13	#include "World.hpp"

src/ConfigFileBuffer.hpp

rba9f5b	rbe97a8
17	17	#endif
18	18
19		#include <ios~~tream~~>
	19	#include <iosfwd>
20	20
21	21	/**************************************** forward declarations ***************************/

src/Descriptors/AtomDescriptor.cpp

-              rba9f5b
+              rbe97a8
 #include "World.hpp"
 #include "atom.hpp"
+#include "Patterns/ObservedContainer_impl.hpp"
 #include <boost/bind.hpp>
 …
 using namespace std;
 typedef World::AtomSet::iterator atoms_iter_t;
+typedef World::AtomSet::internal_iterator atoms_iter_t;
 /************************ Forwarding object **************************************/
 …
 atom* AtomDescriptor_impl::find() {
   World::AtomSet atoms = getAtoms();
   atoms_iter_t res = find_if(atoms.begin(),atoms.end(),boost::bind(&AtomDescriptor_impl::predicate,this,_1));
   return (res!=atoms.end())?((*res).second):0;
+  World::AtomSet &atoms = getAtoms();
+  atoms_iter_t res = find_if(atoms.begin_internal(),atoms.end_internal(),boost::bind(&AtomDescriptor_impl::predicate,this,_1));
+  return (res!=atoms.end_internal())?((*res).second):0;
+}
 …
   vector<atom*> res;
   World::AtomSet atoms = getAtoms();
+  atoms_iter_t iter;
+  for(iter=atoms.begin();iter!=atoms.end();++iter) {
+    if(predicate(*iter)){
+      res.push_back((*iter).second);
+    }
+  for_each(atoms.begin_internal(),
+           atoms.end_internal(),
+           boost::bind(&AtomDescriptor_impl::checkAndAdd,
+                       this,&res,_1));
+  return res;
+}
+void AtomDescriptor_impl::checkAndAdd(std::vector<atom*> *v,std::pair<atomId_t,atom*> p){
+  if(predicate(p)){
+    v->push_back(p.second);
+  }
-  return res;
+}

src/Descriptors/AtomDescriptor_impl.hpp

rba9f5b	rbe97a8
50	50	*/
51	51	World::AtomSet& getAtoms();
	52
	53	void checkAndAdd(std::vector<atom>,std::pair<atomId_t,atom*>);
52	54	};
53	55

src/Descriptors/AtomIdDescriptor.cpp

-              rba9f5b
+              rbe97a8
 #include "atom.hpp"
+#include "Patterns/ObservedContainer_impl.hpp"
 using namespace std;
 …
 atom *AtomIdDescriptor_impl::find(){
   World::AtomSet atoms = getAtoms();
+  World::AtomSet &atoms = getAtoms();
   World::AtomSet::iterator res = atoms.find(id);
   return (res!=atoms.end())?((*res).second):0;

src/Descriptors/MoleculeDescriptor.cpp

-              rba9f5b
+              rbe97a8
 #include "World.hpp"
+#include "Patterns/ObservedContainer_impl.hpp"
 #include "molecule.hpp"
 …
 using namespace std;
 typedef World::MoleculeSet::iterator molecules_iter_t;
+typedef World::MoleculeSet::internal_iterator molecules_iter_t;
 /************************ Forwarding object **************************************/
 …
 molecule* MoleculeDescriptor_impl::find() {
   World::MoleculeSet molecules = getMolecules();
   molecules_iter_t res = find_if(molecules.begin(),molecules.end(),boost::bind(&MoleculeDescriptor_impl::predicate,this,_1));
   return (res!=molecules.end())?((*res).second):0;
+  World::MoleculeSet &molecules = getMolecules();
+  molecules_iter_t res = find_if(molecules.begin_internal(),molecules.end_internal(),boost::bind(&MoleculeDescriptor_impl::predicate,this,_1));
+  return (res!=molecules.end_internal())?((*res).second):0;
+}
 vector<molecule*> MoleculeDescriptor_impl::findAll() {
   vector<molecule*> res;
+  World::MoleculeSet molecules = getMolecules();
+  molecules_iter_t iter;
+  for(iter=molecules.begin();iter!=molecules.end();++iter) {
+    if(predicate(*iter)){
+      res.push_back((*iter).second);
+    }
+  World::MoleculeSet &molecules = getMolecules();
+  for_each(molecules.begin_internal(),
+           molecules.end_internal(),
+           boost::bind(&MoleculeDescriptor_impl::checkAndAdd,
+                       this,&res,_1));
+  return res;
+}
+void MoleculeDescriptor_impl::checkAndAdd(std::vector<molecule*> *v,std::pair<moleculeId_t,molecule*> p){
+  if(predicate(p)){
+    v->push_back(p.second);
+  }
-  return res;
+}

src/Descriptors/MoleculeDescriptor_impl.hpp

rba9f5b	rbe97a8
50	50	*/
51	51	World::MoleculeSet& getMolecules();
	52
	53	void checkAndAdd(std::vector<molecule>,std::pair<moleculeId_t,molecule*>);
52	54	};
53	55

src/Descriptors/MoleculeIdDescriptor.cpp

-              rba9f5b
+              rbe97a8
 #include "MoleculeIdDescriptor.hpp"
 #include "MoleculeIdDescriptor_impl.hpp"
+#include "Patterns/ObservedContainer_impl.hpp"
 #include "molecule.hpp"
 …
 molecule *MoleculeIdDescriptor_impl::find(){
   World::MoleculeSet molecules = getMolecules();
+  World::MoleculeSet &molecules = getMolecules();
   World::MoleculeSet::iterator res = molecules.find(id);
   return (res!=molecules.end())?((*res).second):0;

src/Exceptions/CustomException.cpp

rba9f5b	rbe97a8
9	9
10	10	#include "CustomException.hpp"
	11	#include <iostream>
11	12
12	13	using namespace std;

src/Exceptions/CustomException.hpp

-              rba9f5b
+              rbe97a8
 #include <exception>
+#include <iostream>
+#include <iosfwd>
+#include <string>
 /**

src/Helpers/Assert.hpp

rba9f5b	rbe97a8
11	11	#include<sstream>
12	12	#include<string>
13		#include<ios~~tream~~>
	13	#include<iosfwd>
14	14	#include<vector>
15	15	#include<map>

src/Helpers/MemDebug.cpp

-              rba9f5b
+              rbe97a8
  */
+#ifndef NDEBUG
+#ifndef NO_MEMDEBUG
+// NDEBUG implies NO_MEMDEBUG
+#ifdef NDEBUG
+# ifndef NO_MEMDEBUG
+#   define NO_MEMDEBUG
+# endif
+#endif
+// NO_MEMDEBUG and MEMDEBUG are mutually exclusive, but at least one must be set
+#ifdef NO_MEMDEBUG
+# ifdef MEMDEBUG
+#   undef MEMDEBUG
+# endif
+#else
+# ifndef MEMDEBUG
+#   define MEMDEBUG
+# endif
+#endif
+#ifdef MEMDEBUG
 #include <iostream>
 …
+}
 #endif
-#endif

src/Helpers/MemDebug.hpp

-              rba9f5b
+              rbe97a8
  * your sourcefiles.
  */
-#ifndef NDEBUG
-#ifndef NO_MEMDEBUG
+#ifndef MEMDEBUG
+#define MEMDEBUG
+// Set all flags in a way that makes sense
+// NDEBUG implies NO_MEMDEBUG
+#ifdef NDEBUG
+# ifndef NO_MEMDEBUG
+#   define NO_MEMDEBUG
+# endif
 #endif
+// NO_MEMDEBUG and MEMDEBUG are mutually exclusive, but at least one must be set
+#ifdef NO_MEMDEBUG
+# ifdef MEMDEBUG
+#   undef MEMDEBUG
+# endif
+#else
+# ifndef MEMDEBUG
+#   define MEMDEBUG
+# endif
+#endif
+#ifdef MEMDEBUG
 #include <new>
 …
 #endif
+#endif
+#endif
+#ifdef NDEBUG
+#undef MEMDEBUG
+#endif
+#ifndef MEMDEBUG
+#else
 // memory debugging was disabled
 …
 #endif
 #endif /* MEMDEBUG_HPP_ */

src/Line.cpp

-              rba9f5b
+              rbe97a8
+}
+std::vector<Vector> Line::getSphereIntersections() const{
+  std::vector<Vector> res;
+  // line is kept in normalized form, so we can skip a lot of calculations
+  double discriminant = 1-origin->NormSquared();
+  // we might have 2, 1 or 0 solutions, depending on discriminant
+  if(discriminant>=0){
+    if(discriminant==0){
+      res.push_back(*origin);
+    }
+    else{
+      Vector helper = sqrt(discriminant)*(*direction);
+      res.push_back(*origin+helper);
+      res.push_back(*origin-helper);
+    }
+  }
+  return res;
+}
 Line makeLineThrough(const Vector &x1, const Vector &x2){
   if(x1==x2){

src/Line.hpp

rba9f5b	rbe97a8
38	38	Plane getOrthogonalPlane(const Vector &origin) const;
39	39
	40	std::vector<Vector> getSphereIntersections() const;
	41
40	42	private:
41	43	std::auto_ptr<Vector> origin;

src/Makefile.am

-              rba9f5b
+              rbe97a8
   vector.cpp
+LINALGHEADER = gslmatrix.hpp \
+LINALGHEADER = \
+  gslmatrix.hpp \
   gslvector.hpp \
   linearsystemofequations.hpp \
 …
   Actions/MethodAction.hpp \
   Actions/Process.hpp
+EXCEPTIONSOURCE = Exceptions/CustomException.cpp \
+                                  Exceptions/LinearDependenceException.cpp \
+                                  Exceptions/MathException.cpp \
+                                  Exceptions/NotInvertibleException.cpp \
+                                  Exceptions/SkewException.cpp \
+                                  Exceptions/ZeroVectorException.cpp
+EXCEPTIONHEADER = Exceptions/CustomException.hpp \
+                                  Exceptions/LinearDependenceException.hpp \
+                                  Exceptions/MathException.hpp \
+                                  Exceptions/NotInvertibleException.hpp \
+                                  Exceptions/SkewException.hpp \
+                                  Exceptions/ZeroVectorException.hpp
 PARSERSOURCE = \
 …
   Patterns/Observer.hpp \
   Patterns/Singleton.hpp
+SHAPESOURCE = \
+  Shapes/BaseShapes.cpp \
+  Shapes/Shape.cpp \
+  Shapes/ShapeOps.cpp
+SHAPEHEADER = \
+  Shapes/BaseShapes.hpp \
+  Shapes/Shape.hpp \
+  Shapes/ShapeOps.hpp
 QTUIMOC_HEADER = UIElements/QT4/QTDialog.hpp \
 …
 DESCRIPTORSOURCE = Descriptors/AtomDescriptor.cpp \
   Descriptors/AtomIdDescriptor.cpp \
+  Descriptors/AtomSelectionDescriptor.cpp \
   Descriptors/AtomTypeDescriptor.cpp \
   Descriptors/MoleculeDescriptor.cpp \
   Descriptors/MoleculeIdDescriptor.cpp \
   Descriptors/MoleculeNameDescriptor.cpp \
+  Descriptors/MoleculePtrDescriptor.cpp
+  Descriptors/MoleculePtrDescriptor.cpp \
+  Descriptors/MoleculeSelectionDescriptor.cpp
 DESCRIPTORHEADER = Descriptors/AtomDescriptor.hpp \
   Descriptors/AtomIdDescriptor.hpp \
+  Descriptors/AtomSelectionDescriptor.hpp \
   Descriptors/AtomTypeDescriptor.hpp \
   Descriptors/MoleculeDescriptor.hpp \
   Descriptors/MoleculeIdDescriptor.hpp \
   Descriptors/MoleculeNameDescriptor.hpp \
+  Descriptors/MoleculePtrDescriptor.hpp
+EXCEPTIONSOURCE = Exceptions/CustomException.cpp \
+                                  Exceptions/LinearDependenceException.cpp \
+                                  Exceptions/MathException.cpp \
+                                  Exceptions/SkewException.cpp \
+                                  Exceptions/ZeroVectorException.cpp
+EXCEPTIONHEADER = Exceptions/CustomException.hpp \
+                                  Exceptions/LinearDependenceException.hpp \
+                                  Exceptions/MathException.hpp \
+                                  Exceptions/SkewException.hpp \
+                                  Exceptions/ZeroVectorException.hpp
+  Descriptors/MoleculePtrDescriptor.hpp \
+  Descriptors/MoleculeSelectionDescriptor.cpp
 QTUISOURCE = ${QTUIMOC_TARGETS} \
 …
   ${ACTIONSSOURCE} \
   ${ATOMSOURCE} \
+  ${EXCEPTIONSOURCE} \
   ${PATTERNSOURCE} \
   ${PARSERSOURCE} \
+  ${SHAPESOURCE} \
   ${DESCRIPTORSOURCE} \
   ${HELPERSOURCE} \
-  ${EXCEPTIONSOURCE} \
   bond.cpp \
   bondgraph.cpp \
   boundary.cpp \
+  Box.cpp \
   CommandLineParser.cpp \
   config.cpp \
 …
   log.cpp \
   logger.cpp \
+  Matrix.cpp \
   moleculelist.cpp \
   molecule.cpp \
 …
   ${ACTIONSHEADER} \
   ${ATOMHEADER} \
+  ${EXCEPTIONHEADER} \
   ${PARSERHEADER} \
   ${PATTERNHEADER} \
+  ${SHAPEHEADER} \
   ${DESCRIPTORHEADER} \
-  ${EXCEPTIONHEADER} \
   bond.hpp \
   bondgraph.hpp \
   boundary.hpp \
+  Box.hpp \
   CommandLineParser.hpp \
   config.hpp \
 …
   log.hpp \
   logger.hpp \
+  Matrix.hpp \
   molecule.hpp \
   molecule_template.hpp \

src/Parser/MpqcParser.hpp

rba9f5b	rbe97a8
11	11	#include "FormatParser.hpp"
12	12
13		#include <ios~~tream~~>
	13	#include <iosfwd>
14	14
15	15	/**

src/Parser/PcpParser.cpp

-              rba9f5b
+              rbe97a8
 #include <iostream>
+#include <iomanip>
 #include "atom.hpp"
 …
 #include "verbose.hpp"
 #include "World.hpp"
+#include "Matrix.hpp"
+#include "Box.hpp"
 /** Constructor of PcpParser.
 …
   // Unit cell and magnetic field
   ParseForParameter(verbose,FileBuffer, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
   double * const cell_size = World::getInstance().getDomain();
+  double *cell_size = new double[6];
   cell_size[0] = BoxLength[0];
   cell_size[1] = BoxLength[3];
 …
   cell_size[4] = BoxLength[7];
   cell_size[5] = BoxLength[8];
+  World::getInstance().setDomain(cell_size);
+  delete[] cell_size;
   //if (1) fprintf(stderr,"\n");
 …
 void PcpParser::save(std::ostream* file)
+{
   const double * const cell_size = World::getInstance().getDomain();
+  const Matrix &domain = World::getInstance().getDomain().getM();
   class ThermoStatContainer *Thermostats = World::getInstance().getThermostats();
   if (!file->fail()) {
 …
     *file << endl;
     *file << "BoxLength\t\t\t# (Length of a unit cell)" << endl;
     *file << cell_size[0] << "\t" << endl;
     *file << cell_size[1] << "\t" << cell_size[2] << "\t" << endl;
     *file << cell_size[3] << "\t" << cell_size[4] << "\t" << cell_size[5] << "\t" << endl;
+    *file << domain.at(0,0) << "\t" << endl;
+    *file << domain.at(1,0) << "\t" << domain.at(1,1) << "\t" << endl;
+    *file << domain.at(2,0) << "\t" << domain.at(2,1) << "\t" << domain.at(2,2) << "\t" << endl;
     // FIXME
     *file << endl;

src/Parser/PcpParser.hpp

rba9f5b	rbe97a8
9	9	#define PCPPARSER_HPP_
10	10
11		#include <ios~~tream~~>
	11	#include <iosfwd>
12	12	#include "Parser/FormatParser.hpp"
13	13

src/Patterns/Cacheable.hpp

rba9f5b	rbe97a8
12	12	#include <boost/function.hpp>
13	13	#include <boost/shared_ptr.hpp>
14		~~#include <iostream>~~
15	14
16	15	#include "Helpers/Assert.hpp"

src/Patterns/ObservedIterator.hpp

-              rba9f5b
+              rbe97a8
+  }
+  value_type *operator->(){
+    acquireLock();
+    return &(*iter);
+  }
   // when we turn into a const iterator we can loose our lock
   operator typename _Set::const_iterator() {

src/Patterns/Observer.hpp

-              rba9f5b
+              rbe97a8
   //! @cond
   // Structure for RAII-Style notification
 protected:
+public:
   /**
    * This structure implements the Observer-mechanism RAII-Idiom.
 …
 #define OBSERVE Observable::_Observable_protector PASTE(_scope_obs_protector_,__LINE__)(this)
 #define NOTIFY(notification) do{Observable::enque_notification_internal(this,notification);}while(0)
+#define LOCK_OBSERVABLE(observable) Observable::_Observable_protector PASTE(_scope_obs_protector_,__LINE__)(&(observable))
 #endif /* OBSERVER_HPP_ */

src/Plane.hpp

rba9f5b	rbe97a8
11	11	#include <memory>
12	12	#include <vector>
13		#include <ios~~tream~~>
	13	#include <iosfwd>
14	14	#include "Space.hpp"
15	15	#include "Exceptions/LinearDependenceException.hpp"

src/UIElements/CommandLineUI/CommandLineDialog.cpp

-              rba9f5b
+              rbe97a8
 #include "verbose.hpp"
 #include "World.hpp"
+#include "Box.hpp"
 #include "atom.hpp"
 …
+}
 void CommandLineDialog::queryVector(const char* title, Vector *target,const double *const cellSize, bool check, string _description) {
   registerQuery(new VectorCommandLineQuery(title,target,cellSize,check, _description));
+}
 void CommandLineDialog::queryBox(const char* title, double ** const cellSize, string _description) {
+void CommandLineDialog::queryVector(const char* title, Vector *target, bool check, string _description) {
+  registerQuery(new VectorCommandLineQuery(title,target,check, _description));
+}
+void CommandLineDialog::queryBox(const char* title, Box* cellSize, string _description) {
   registerQuery(new BoxCommandLineQuery(title,cellSize,_description));
+}
 …
+}
 CommandLineDialog::VectorCommandLineQuery::VectorCommandLineQuery(string title, Vector *_target, const double *const _cellSize, bool _check, string _description) :
     Dialog::VectorQuery(title,_target,_cellSize,_check, _description)
+CommandLineDialog::VectorCommandLineQuery::VectorCommandLineQuery(string title, Vector *_target, bool _check, string _description) :
+    Dialog::VectorQuery(title,_target,_check, _description)
 {}
 …
 CommandLineDialog::BoxCommandLineQuery::BoxCommandLineQuery(string title, double ** const _cellSize, string _description) :
+CommandLineDialog::BoxCommandLineQuery::BoxCommandLineQuery(string title, Box* _cellSize, string _description) :
     Dialog::BoxQuery(title,_cellSize, _description)
 {}

src/UIElements/CommandLineUI/CommandLineDialog.hpp

-              rba9f5b
+              rbe97a8
   virtual void queryAtom(const char*,atom**, std::string = "");
   virtual void queryMolecule(const char*,molecule**,std::string = "");
   virtual void queryVector(const char*,Vector *,const double * const,bool, std::string = "");
   virtual void queryBox(const char*,double ** const, std::string = "");
+  virtual void queryVector(const char*,Vector *,bool, std::string = "");
+  virtual void queryBox(const char*,Box *, std::string = "");
   virtual void queryElement(const char*, std::vector<element *> *, std::string = "");
 …
   class VectorCommandLineQuery : public Dialog::VectorQuery {
   public:
     VectorCommandLineQuery(std::string title,Vector *_target,const double *const _cellSize,bool _check, std::string _description = "");
+    VectorCommandLineQuery(std::string title,Vector *_target,bool _check, std::string _description = "");
     virtual ~VectorCommandLineQuery();
     virtual bool handle();
 …
   class BoxCommandLineQuery : public Dialog::BoxQuery {
   public:
     BoxCommandLineQuery(std::string title,double ** const _cellSize, std::string _description = "");
+    BoxCommandLineQuery(std::string title,Box* _cellSize, std::string _description = "");
     virtual ~BoxCommandLineQuery();
     virtual bool handle();

src/UIElements/Dialog.cpp

-              rba9f5b
+              rbe97a8
 #include "Dialog.hpp"
+#include "verbose.hpp"
 #include "atom.hpp"
 #include "element.hpp"
 #include "molecule.hpp"
 #include "vector.hpp"
+#include "Matrix.hpp"
+#include "Box.hpp"
 using namespace std;
 …
 // Vector Queries
 Dialog::VectorQuery::VectorQuery(std::string title,Vector *_target,const double *const _cellSize,bool _check, std::string _description) :
+Dialog::VectorQuery::VectorQuery(std::string title,Vector *_target,bool _check, std::string _description) :
   Query(title, _description),
-  cellSize(_cellSize),
   check(_check),
   target(_target)
 …
 // Box Queries
 Dialog::BoxQuery::BoxQuery(std::string title, double ** const _cellSize, std::string _description) :
+Dialog::BoxQuery::BoxQuery(std::string title, Box* _cellSize, std::string _description) :
   Query(title, _description),
   target(_cellSize)
 …
 void Dialog::BoxQuery::setResult() {
+  for (int i=0;i<6;i++) {
+    (*target)[i] = tmp[i];
+  }
+  (*target)= ReturnFullMatrixforSymmetric(tmp);
+}

src/UIElements/Dialog.hpp

-              rba9f5b
+              rbe97a8
 class atom;
+class Box;
 class element;
 class molecule;
 …
   virtual void queryAtom(const char*,atom**,std::string = "")=0;
   virtual void queryMolecule(const char*,molecule**, std::string = "")=0;
   virtual void queryVector(const char*,Vector *,const double *const,bool, std::string = "")=0;
   virtual void queryBox(const char*,double ** const, std::string = "")=0;
+  virtual void queryVector(const char*,Vector *,bool, std::string = "")=0;
+  virtual void queryBox(const char*,Box*, std::string = "")=0;
   virtual void queryElement(const char*, std::vector<element *> *, std::string = "")=0;
 …
   class VectorQuery : public Query {
   public:
       VectorQuery(std::string title,Vector *_target,const double *const _cellSize,bool _check, std::string _description = "");
+      VectorQuery(std::string title,Vector *_target,bool _check, std::string _description = "");
       virtual ~VectorQuery();
       virtual bool handle()=0;
 …
     protected:
       Vector *tmp;
-      const double *const cellSize;
       bool check;
     private:
 …
   class BoxQuery : public Query {
   public:
       BoxQuery(std::string title,double ** const _cellSize, std::string _description = "");
+      BoxQuery(std::string title,Box *_cellSize, std::string _description = "");
       virtual ~BoxQuery();
       virtual bool handle()=0;
       virtual void setResult();
     protected:
       double *tmp;
+      double* tmp;
     private:
       double **target;
+      Box* target;
   };

src/UIElements/QT4/QTDialog.cpp

-              rba9f5b
+              rbe97a8
 #include "molecule.hpp"
 #include "Descriptors/MoleculeIdDescriptor.hpp"
+#include "Matrix.hpp"
+#include "Box.hpp"
 …
+}
 void QTDialog::queryBox(char const*, double**, string){
+void QTDialog::queryBox(char const*, Box*, string){
   // TODO
   ASSERT(false, "Not implemented yet");
 …
+}
 void QTDialog::queryVector(const char* title, Vector *target,const double *const cellSize, bool check,string) {
   registerQuery(new VectorQTQuery(title,target,cellSize,check,inputLayout,this));
+void QTDialog::queryVector(const char* title, Vector *target, bool check,string) {
+  registerQuery(new VectorQTQuery(title,target,check,inputLayout,this));
+}
 …
+}
+QTDialog::VectorQTQuery::VectorQTQuery(std::string title, Vector *_target, const double *const _cellSize, bool _check,QBoxLayout *_parent,QTDialog *_dialog) :
+    Dialog::VectorQuery(title,_target,_cellSize,_check),
+    parent(_parent)
+{
+  // About the j: I don't know why it was done this way, but it was used this way in Vector::AskPosition, so I reused it
+  int j = -1;
+QTDialog::VectorQTQuery::VectorQTQuery(std::string title, Vector *_target, bool _check,QBoxLayout *_parent,QTDialog *_dialog) :
+    Dialog::VectorQuery(title,_target,_check),
+    parent(_parent)
+{
+  const Matrix& M = World::getInstance().getDomain().getM();
   const char *coords[3] = {"x:","y:","z:"};
   mainLayout= new QHBoxLayout();
 …
   mainLayout->addLayout(subLayout);
   for(int i=0; i<3; i++) {
-    j+=i+1;
     coordLayout[i] = new QHBoxLayout();
     subLayout->addLayout(coordLayout[i]);
 …
     coordLayout[i]->addWidget(coordLabel[i]);
     coordInput[i] = new QDoubleSpinBox();
     coordInput[i]->setRange(0,cellSize[j]);
+    coordInput[i]->setRange(0,M.at(i,i));
     coordInput[i]->setDecimals(3);
     coordLayout[i]->addWidget(coordInput[i]);

src/UIElements/QT4/QTDialog.hpp

-              rba9f5b
+              rbe97a8
   virtual void queryAtom(const char*,atom**,std::string = "");
   virtual void queryMolecule(const char*,molecule**,std::string = "");
   virtual void queryVector(const char*,Vector *,const double *const,bool,std::string = "");
   virtual void queryBox(const char*,double ** const, std::string = "");
+  virtual void queryVector(const char*,Vector *,bool,std::string = "");
+  virtual void queryBox(const char*,Box*, std::string = "");
   virtual void queryElement(const char*,std::vector<element *> *_target,std::string = "");
 …
     class VectorQTQuery : public Dialog::VectorQuery {
     public:
       VectorQTQuery(std::string title,Vector *_target,const double *const _cellSize,bool _check,QBoxLayout *,QTDialog *);
+      VectorQTQuery(std::string title,Vector *_target,bool _check,QBoxLayout *,QTDialog *);
       virtual ~VectorQTQuery();
       virtual bool handle();

src/UIElements/QT4/QTMainWindow.cpp

rba9f5b	rbe97a8
21	21	#include "atom.hpp"
22	22	#include "molecule.hpp"
	23	#include "verbose.hpp"
23	24	#include "Actions/Action.hpp"
24	25	#include "Actions/ActionRegistry.hpp"

src/UIElements/TextUI/TextDialog.cpp

-              rba9f5b
+              rbe97a8
 #include "molecule.hpp"
 #include "vector.hpp"
+#include "Matrix.hpp"
+#include "Box.hpp"
 using namespace std;
 …
+}
 void TextDialog::queryVector(const char* title, Vector *target,const double *const cellSize, bool check, string description) {
   registerQuery(new VectorTextQuery(title,target,cellSize,check,description));
+}
 void TextDialog::queryBox(const char* title,double ** const cellSize, string description) {
+void TextDialog::queryVector(const char* title, Vector *target, bool check, string description) {
+  registerQuery(new VectorTextQuery(title,target,check,description));
+}
+void TextDialog::queryBox(const char* title,Box* cellSize, string description) {
   registerQuery(new BoxTextQuery(title,cellSize,description));
+}
 …
+}
 TextDialog::VectorTextQuery::VectorTextQuery(std::string title, Vector *_target, const double *const _cellSize, bool _check, std::string _description) :
     Dialog::VectorQuery(title,_target,_cellSize,_check,_description)
+TextDialog::VectorTextQuery::VectorTextQuery(std::string title, Vector *_target, bool _check, std::string _description) :
+    Dialog::VectorQuery(title,_target,_check,_description)
 {}
 …
   Log() << Verbose(0) << getTitle();
+  const Matrix &M = World::getInstance().getDomain().getM();
   char coords[3] = {'x','y','z'};
-  int j = -1;
   for (int i=0;i<3;i++) {
-    j += i+1;
     do {
       Log() << Verbose(0) << coords[i] << "[0.." << cellSize[j] << "]: ";
+      Log() << Verbose(0) << coords[i] << "[0.." << M.at(i,i) << "]: ";
       cin >> (*tmp)[i];
     } while ((((*tmp)[i] < 0) || ((*tmp)[i] >= cellSize[j])) && (check));
+    } while ((((*tmp)[i] < 0) || ((*tmp)[i] >= M.at(i,i))) && (check));
+  }
   return true;
+}
 TextDialog::BoxTextQuery::BoxTextQuery(std::string title, double ** const _cellSize, std::string _description) :
+TextDialog::BoxTextQuery::BoxTextQuery(std::string title, Box* _cellSize, std::string _description) :
     Dialog::BoxQuery(title,_cellSize,_description)
 {}

src/UIElements/TextUI/TextDialog.hpp

-              rba9f5b
+              rbe97a8
   virtual void queryAtom(const char*,atom**,std::string = "");
   virtual void queryMolecule(const char*,molecule**,std::string = "");
   virtual void queryVector(const char*,Vector *,const double * const,bool, std::string = "");
   virtual void queryBox(const char*,double ** const, std::string = "");
+  virtual void queryVector(const char*,Vector *,bool, std::string = "");
+  virtual void queryBox(const char*,Box*, std::string = "");
   virtual void queryElement(const char*, std::vector<element *> *, std::string = "");
 …
   class VectorTextQuery : public Dialog::VectorQuery {
   public:
     VectorTextQuery(std::string title,Vector *_target,const double *const _cellSize,bool _check, std::string _description = NULL);
+    VectorTextQuery(std::string title,Vector *_target,bool _check, std::string _description = NULL);
     virtual ~VectorTextQuery();
     virtual bool handle();
 …
   class BoxTextQuery : public Dialog::BoxQuery {
   public:
     BoxTextQuery(std::string title,double ** const _cellSize, std::string _description = NULL);
+    BoxTextQuery(std::string title,Box* _cellSize, std::string _description = NULL);
     virtual ~BoxTextQuery();
     virtual bool handle();

src/UIElements/TextUI/TextWindow.hpp

rba9f5b	rbe97a8
11	11	#include "MainWindow.hpp"
12	12
	13	#include <string>
13	14	#include <set>
14	15

src/UIElements/Views/StreamStringView.hpp

rba9f5b	rbe97a8
10	10
11	11	#include <boost/function.hpp>
12		#include <ios~~tream~~>
	12	#include <iosfwd>
13	13
14	14	#include "Views/StringView.hpp"

src/VectorSet.hpp

-              rba9f5b
+              rbe97a8
 #include <functional>
 #include <algorithm>
+#include <limits>
 /**
 …
  */
+class Vector;
+#include "vector.hpp"
+#include <list>
 // this tests, whether we actually have a Vector
 …
   void translate(const Vector &translater){
     // this is needed to allow template lookup
+    transform(this->begin(),this->end(),this->begin(),bind1st(plus<Vector>(),translater));
+    transform(this->begin(),this->end(),this->begin(),std::bind1st(std::plus<Vector>(),translater));
+  }
+  double minDistSquared(const Vector &point){
+    if(!this->size())
+      return std::numeric_limits<double>::infinity();
+    std::list<double> helper;
+    helper.resize(this->size());
+    transform(this->begin(),this->end(),
+              helper.begin(),
+              std::bind2nd(std::mem_fun_ref(&Vector::DistanceSquared),point));
+    return *min_element(helper.begin(),helper.end());
+  }
 };
+// allows simpler definition of VectorSets
+#define VECTORSET(container_type) VectorSet<container_type<Vector> >
 #endif /* VECTORSET_HPP_ */

src/World.cpp

-              rba9f5b
+              rbe97a8
 #include "World.hpp"
+#include <functional>
 #include "atom.hpp"
 …
 #include "Actions/ManipulateAtomsProcess.hpp"
 #include "Helpers/Assert.hpp"
+#include "Box.hpp"
+#include "Matrix.hpp"
+#include "defs.hpp"
 #include "Patterns/Singleton_impl.hpp"
+#include "Patterns/ObservedContainer_impl.hpp"
 using namespace std;
 …
 // system
+double * World::getDomain() {
+  return cell_size;
+Box& World::getDomain() {
+  return *cell_size;
+}
+void World::setDomain(const Matrix &mat){
+  OBSERVE;
+  *cell_size = mat;
+}
 …
+{
   OBSERVE;
+  cell_size[0] = matrix[0];
+  cell_size[1] = matrix[1];
+  cell_size[2] = matrix[2];
+  cell_size[3] = matrix[3];
+  cell_size[4] = matrix[4];
+  cell_size[5] = matrix[5];
+  Matrix M = ReturnFullMatrixforSymmetric(matrix);
+  cell_size->setM(M);
+}
 …
 void World::setDefaultName(std::string name)
+{
+  OBSERVE;
   defaultName = name;
 };
 …
   molecule *mol = NULL;
   mol = NewMolecule();
+  ASSERT(!molecules.count(currMoleculeId),"currMoleculeId did not specify an unused ID");
+  mol->setId(currMoleculeId++);
+  moleculeId_t id = getNextMoleculeId();
+  ASSERT(!molecules.count(id),"proposed id did not specify an unused ID");
+  mol->setId(id);
   // store the molecule by ID
   molecules[mol->getId()] = mol;
 …
   DeleteMolecule(mol);
   molecules.erase(id);
+}
+double *World::cell_size = NULL;
+  releaseMoleculeId(id);
+}
 atom *World::createAtom(){
   OBSERVE;
   atomId_t id = getNextAtomId();
+  ASSERT(!atoms.count(id),"proposed id did not specify an unused ID");
   atom *res = NewAtom(id);
   res->setWorld(this);
 …
 atomId_t World::getNextAtomId(){
+  // see if we can reuse some Id
+  if(atomIdPool.empty()){
+    return currAtomId++;
+  }
+  else{
+    // we give out the first ID from the pool
+    atomId_t id = *(atomIdPool.begin());
+    atomIdPool.erase(id);
+  // try to find an Id in the pool;
+  if(!atomIdPool.empty()){
+    atomIdPool_t::iterator iter=atomIdPool.begin();
+    atomId_t id = iter->first;
+    pair<atomId_t,atomId_t> newRange = make_pair(id+1,iter->second);
+    // we wont use this iterator anymore, so we don't care about invalidating
+    atomIdPool.erase(iter);
+    if(newRange.first<newRange.second){
+      atomIdPool.insert(newRange);
+    }
     return id;
+  }
+  // Nothing in the pool... we are out of luck
+  return currAtomId++;
+}
 void World::releaseAtomId(atomId_t id){
+  atomIdPool.insert(id);
+  // defragmentation of the pool
+  set<atomId_t>::reverse_iterator iter;
+  // go through all Ids in the pool that lie immediately below the border
+  while(!atomIdPool.empty() && *(atomIdPool.rbegin())==(currAtomId-1)){
+    atomIdPool.erase(--currAtomId);
+  }
+  atomIdPool.insert(make_pair(id,id+1));
+  defragAtomIdPool();
+}
 bool World::reserveAtomId(atomId_t id){
   if(id>=currAtomId ){
     // add all ids between the new one and current border as available
     for(atomId_t pos=currAtomId; pos<id; ++pos){
       atomIdPool.insert(pos);
+    pair<atomId_t,atomId_t> newRange = make_pair(currAtomId,id);
+    if(newRange.first<newRange.second){
+      atomIdPool.insert(newRange);
+    }
     currAtomId=id+1;
+    defragAtomIdPool();
     return true;
+  }
+  else if(atomIdPool.count(id)){
+    atomIdPool.erase(id);
+  // look for a range that matches the request
+  for(atomIdPool_t::iterator iter=atomIdPool.begin();iter!=atomIdPool.end();++iter){
+    if(iter->first>id){
+      // we have coverd all available ranges... nothing to be found here
+      break;
+    }
+    // no need to check first, since it has to be <=id, since otherwise we would have broken out
+    if(iter->second > id){
+      // we found a matching range... get the id from this range
+      // split up this range at the point of id
+      pair<atomId_t,atomId_t> bottomRange = make_pair(iter->first,id);
+      pair<atomId_t,atomId_t> topRange = make_pair(id+1,iter->second);
+      // remove this range
+      atomIdPool.erase(iter);
+      if(bottomRange.first<bottomRange.second){
+        atomIdPool.insert(bottomRange);
+      }
+      if(topRange.first<topRange.second){
+        atomIdPool.insert(topRange);
+      }
+      defragAtomIdPool();
+      return true;
+    }
+  }
+  // this ID could not be reserved
+  return false;
+}
+void World::defragAtomIdPool(){
+  // check if the situation is bad enough to make defragging neccessary
+  if((numAtomDefragSkips<MAX_FRAGMENTATION_SKIPS) &&
+     (atomIdPool.size()<lastAtomPoolSize+MAX_POOL_FRAGMENTATION)){
+    ++numAtomDefragSkips;
+    return;
+  }
+  for(atomIdPool_t::iterator iter = atomIdPool.begin();iter!=atomIdPool.end();){
+    // see if this range is adjacent to the next one
+    atomIdPool_t::iterator next = iter;
+    next++;
+    if(next!=atomIdPool.end() && (next->first==iter->second)){
+      // merge the two ranges
+      pair<atomId_t,atomId_t> newRange = make_pair(iter->first,next->second);
+      atomIdPool.erase(iter);
+      atomIdPool.erase(next);
+      pair<atomIdPool_t::iterator,bool> res = atomIdPool.insert(newRange);
+      ASSERT(res.second,"Id-Pool was confused");
+      iter=res.first;
+      continue;
+    }
+    ++iter;
+  }
+  if(!atomIdPool.empty()){
+    // check if the last range is at the border
+    atomIdPool_t::iterator iter = atomIdPool.end();
+    iter--;
+    if(iter->second==currAtomId){
+      currAtomId=iter->first;
+      atomIdPool.erase(iter);
+    }
+  }
+  lastAtomPoolSize=atomIdPool.size();
+  numAtomDefragSkips=0;
+}
+// Molecules
+moleculeId_t World::getNextMoleculeId(){
+  // try to find an Id in the pool;
+  if(!moleculeIdPool.empty()){
+    moleculeIdPool_t::iterator iter=moleculeIdPool.begin();
+    moleculeId_t id = iter->first;
+    pair<moleculeId_t,moleculeId_t> newRange = make_pair(id+1,iter->second);
+    // we wont use this iterator anymore, so we don't care about invalidating
+    moleculeIdPool.erase(iter);
+    if(newRange.first<newRange.second){
+      moleculeIdPool.insert(newRange);
+    }
+    return id;
+  }
+  // Nothing in the pool... we are out of luck
+  return currMoleculeId++;
+}
+void World::releaseMoleculeId(moleculeId_t id){
+  moleculeIdPool.insert(make_pair(id,id+1));
+  defragMoleculeIdPool();
+}
+bool World::reserveMoleculeId(moleculeId_t id){
+  if(id>=currMoleculeId ){
+    pair<moleculeId_t,moleculeId_t> newRange = make_pair(currMoleculeId,id);
+    if(newRange.first<newRange.second){
+      moleculeIdPool.insert(newRange);
+    }
+    currMoleculeId=id+1;
+    defragMoleculeIdPool();
     return true;
+  }
+  else{
+    // this ID could not be reserved
+    return false;
+  }
+  // look for a range that matches the request
+  for(moleculeIdPool_t::iterator iter=moleculeIdPool.begin();iter!=moleculeIdPool.end();++iter){
+    if(iter->first>id){
+      // we have coverd all available ranges... nothing to be found here
+      break;
+    }
+    // no need to check first, since it has to be <=id, since otherwise we would have broken out
+    if(iter->second > id){
+      // we found a matching range... get the id from this range
+      // split up this range at the point of id
+      pair<moleculeId_t,moleculeId_t> bottomRange = make_pair(iter->first,id);
+      pair<moleculeId_t,moleculeId_t> topRange = make_pair(id+1,iter->second);
+      // remove this range
+      moleculeIdPool.erase(iter);
+      if(bottomRange.first<bottomRange.second){
+        moleculeIdPool.insert(bottomRange);
+      }
+      if(topRange.first<topRange.second){
+        moleculeIdPool.insert(topRange);
+      }
+      defragMoleculeIdPool();
+      return true;
+    }
+  }
+  // this ID could not be reserved
+  return false;
+}
+void World::defragMoleculeIdPool(){
+  // check if the situation is bad enough to make defragging neccessary
+  if((numMoleculeDefragSkips<MAX_FRAGMENTATION_SKIPS) &&
+     (moleculeIdPool.size()<lastMoleculePoolSize+MAX_POOL_FRAGMENTATION)){
+    ++numMoleculeDefragSkips;
+    return;
+  }
+  for(moleculeIdPool_t::iterator iter = moleculeIdPool.begin();iter!=moleculeIdPool.end();){
+    // see if this range is adjacent to the next one
+    moleculeIdPool_t::iterator next = iter;
+    next++;
+    if(next!=moleculeIdPool.end() && (next->first==iter->second)){
+      // merge the two ranges
+      pair<moleculeId_t,moleculeId_t> newRange = make_pair(iter->first,next->second);
+      moleculeIdPool.erase(iter);
+      moleculeIdPool.erase(next);
+      pair<moleculeIdPool_t::iterator,bool> res = moleculeIdPool.insert(newRange);
+      ASSERT(res.second,"Id-Pool was confused");
+      iter=res.first;
+      continue;
+    }
+    ++iter;
+  }
+  if(!moleculeIdPool.empty()){
+    // check if the last range is at the border
+    moleculeIdPool_t::iterator iter = moleculeIdPool.end();
+    iter--;
+    if(iter->second==currMoleculeId){
+      currMoleculeId=iter->first;
+      moleculeIdPool.erase(iter);
+    }
+  }
+  lastMoleculePoolSize=moleculeIdPool.size();
+  numMoleculeDefragSkips=0;
+}
+/******************************* Iterators ********************************/
+// external parts with observers
+CONSTRUCT_SELECTIVE_ITERATOR(atom*,World::AtomSet,AtomDescriptor);
+World::AtomIterator
+World::getAtomIter(AtomDescriptor descr){
+    return AtomIterator(descr,atoms);
+}
+World::AtomIterator
+World::getAtomIter(){
+    return AtomIterator(AllAtoms(),atoms);
+}
+World::AtomIterator
+World::atomEnd(){
+  return AtomIterator(AllAtoms(),atoms,atoms.end());
+}
+CONSTRUCT_SELECTIVE_ITERATOR(molecule*,World::MoleculeSet,MoleculeDescriptor);
+World::MoleculeIterator
+World::getMoleculeIter(MoleculeDescriptor descr){
+    return MoleculeIterator(descr,molecules);
+}
+World::MoleculeIterator
+World::getMoleculeIter(){
+    return MoleculeIterator(AllMolecules(),molecules);
+}
+World::MoleculeIterator
+World::moleculeEnd(){
+  return MoleculeIterator(AllMolecules(),molecules,molecules.end());
+}
+// Internal parts, without observers
+// Build the AtomIterator from template
+CONSTRUCT_SELECTIVE_ITERATOR(atom*,World::AtomSet::set_t,AtomDescriptor);
+World::internal_AtomIterator
+World::getAtomIter_internal(AtomDescriptor descr){
+  return internal_AtomIterator(descr,atoms.getContent());
+}
+World::internal_AtomIterator
+World::atomEnd_internal(){
+  return internal_AtomIterator(AllAtoms(),atoms.getContent(),atoms.end_internal());
+}
+// build the MoleculeIterator from template
+CONSTRUCT_SELECTIVE_ITERATOR(molecule*,World::MoleculeSet::set_t,MoleculeDescriptor);
+World::internal_MoleculeIterator World::getMoleculeIter_internal(MoleculeDescriptor descr){
+  return internal_MoleculeIterator(descr,molecules.getContent());
+}
+World::internal_MoleculeIterator World::moleculeEnd_internal(){
+  return internal_MoleculeIterator(AllMolecules(),molecules.getContent(),molecules.end_internal());
+}
+/************************** Selection of Atoms and molecules ******************/
+// Atoms
+void World::clearAtomSelection(){
+  selectedAtoms.clear();
+}
+void World::selectAtom(atom *atom){
+  ASSERT(atom,"Invalid pointer in selection of atom");
+  selectedAtoms[atom->getId()]=atom;
+}
+void World::selectAtom(atomId_t id){
+  ASSERT(atoms.count(id),"Atom Id selected that was not in the world");
+  selectedAtoms[id]=atoms[id];
+}
+void World::selectAllAtoms(AtomDescriptor descr){
+  internal_AtomIterator begin = getAtomIter_internal(descr);
+  internal_AtomIterator end = atomEnd_internal();
+  void (World::*func)(atom*) = &World::selectAtom; // needed for type resolution of overloaded function
+  for_each(begin,end,bind1st(mem_fun(func),this)); // func is select... see above
+}
+void World::selectAtomsOfMolecule(molecule *_mol){
+  ASSERT(_mol,"Invalid pointer to molecule in selection of Atoms of Molecule");
+  // need to make it const to get the fast iterators
+  const molecule *mol = _mol;
+  void (World::*func)(atom*) = &World::selectAtom; // needed for type resolution of overloaded function
+  for_each(mol->begin(),mol->end(),bind1st(mem_fun(func),this)); // func is select... see above
+}
+void World::selectAtomsOfMolecule(moleculeId_t id){
+  ASSERT(molecules.count(id),"No molecule with the given id upon Selection of atoms from molecule");
+  selectAtomsOfMolecule(molecules[id]);
+}
+void World::unselectAtom(atom *atom){
+  ASSERT(atom,"Invalid pointer in unselection of atom");
+  unselectAtom(atom->getId());
+}
+void World::unselectAtom(atomId_t id){
+  ASSERT(atoms.count(id),"Atom Id unselected that was not in the world");
+  selectedAtoms.erase(id);
+}
+void World::unselectAllAtoms(AtomDescriptor descr){
+  internal_AtomIterator begin = getAtomIter_internal(descr);
+  internal_AtomIterator end = atomEnd_internal();
+  void (World::*func)(atom*) = &World::unselectAtom; // needed for type resolution of overloaded function
+  for_each(begin,end,bind1st(mem_fun(func),this)); // func is unselect... see above
+}
+void World::unselectAtomsOfMolecule(molecule *_mol){
+  ASSERT(_mol,"Invalid pointer to molecule in selection of Atoms of Molecule");
+  // need to make it const to get the fast iterators
+  const molecule *mol = _mol;
+  void (World::*func)(atom*) = &World::unselectAtom; // needed for type resolution of overloaded function
+  for_each(mol->begin(),mol->end(),bind1st(mem_fun(func),this)); // func is unsselect... see above
+}
+void World::unselectAtomsOfMolecule(moleculeId_t id){
+  ASSERT(molecules.count(id),"No molecule with the given id upon Selection of atoms from molecule");
+  unselectAtomsOfMolecule(molecules[id]);
+}
 // Molecules
+/******************************* Iterators ********************************/
+// Build the AtomIterator from template
+CONSTRUCT_SELECTIVE_ITERATOR(atom*,World::AtomSet,AtomDescriptor);
+World::AtomIterator World::getAtomIter(AtomDescriptor descr){
+  return AtomIterator(descr,atoms);
+}
+World::AtomIterator World::atomEnd(){
+  return AtomIterator(AllAtoms(),atoms,atoms.end());
+}
+// build the MoleculeIterator from template
+CONSTRUCT_SELECTIVE_ITERATOR(molecule*,World::MoleculeSet,MoleculeDescriptor);
+World::MoleculeIterator World::getMoleculeIter(MoleculeDescriptor descr){
+  return MoleculeIterator(descr,molecules);
+}
+World::MoleculeIterator World::moleculeEnd(){
+  return MoleculeIterator(AllMolecules(),molecules,molecules.end());
+void World::clearMoleculeSelection(){
+  selectedMolecules.clear();
+}
+void World::selectMolecule(molecule *mol){
+  ASSERT(mol,"Invalid pointer to molecule in selection");
+  selectedMolecules[mol->getId()]=mol;
+}
+void World::selectMolecule(moleculeId_t id){
+  ASSERT(molecules.count(id),"Molecule Id selected that was not in the world");
+  selectedMolecules[id]=molecules[id];
+}
+void World::selectAllMoleculess(MoleculeDescriptor descr){
+  internal_MoleculeIterator begin = getMoleculeIter_internal(descr);
+  internal_MoleculeIterator end = moleculeEnd_internal();
+  void (World::*func)(molecule*) = &World::selectMolecule; // needed for type resolution of overloaded function
+  for_each(begin,end,bind1st(mem_fun(func),this)); // func is select... see above
+}
+void World::selectMoleculeOfAtom(atom *atom){
+  ASSERT(atom,"Invalid atom pointer in selection of MoleculeOfAtom");
+  molecule *mol=atom->getMolecule();
+  // the atom might not be part of a molecule
+  if(mol){
+    selectMolecule(mol);
+  }
+}
+void World::selectMoleculeOfAtom(atomId_t id){
+  ASSERT(atoms.count(id),"No such atom with given ID in selection of Molecules of Atom");\
+  selectMoleculeOfAtom(atoms[id]);
+}
+void World::unselectMolecule(molecule *mol){
+  ASSERT(mol,"invalid pointer in unselection of molecule");
+  unselectMolecule(mol->getId());
+}
+void World::unselectMolecule(moleculeId_t id){
+  ASSERT(molecules.count(id),"No such molecule with ID in unselection");
+  selectedMolecules.erase(id);
+}
+void World::unselectAllMoleculess(MoleculeDescriptor descr){
+  internal_MoleculeIterator begin = getMoleculeIter_internal(descr);
+  internal_MoleculeIterator end = moleculeEnd_internal();
+  void (World::*func)(molecule*) = &World::unselectMolecule; // needed for type resolution of overloaded function
+  for_each(begin,end,bind1st(mem_fun(func),this)); // func is unselect... see above
+}
+void World::unselectMoleculeOfAtom(atom *atom){
+  ASSERT(atom,"Invalid atom pointer in selection of MoleculeOfAtom");
+  molecule *mol=atom->getMolecule();
+  // the atom might not be part of a molecule
+  if(mol){
+    unselectMolecule(mol);
+  }
+}
+void World::unselectMoleculeOfAtom(atomId_t id){
+  ASSERT(atoms.count(id),"No such atom with given ID in selection of Molecules of Atom");\
+  unselectMoleculeOfAtom(atoms[id]);
+}
+/******************* Iterators over Selection *****************************/
+World::AtomSelectionIterator World::beginAtomSelection(){
+  return selectedAtoms.begin();
+}
+World::AtomSelectionIterator World::endAtomSelection(){
+  return selectedAtoms.end();
+}
+World::MoleculeSelectionIterator World::beginMoleculeSelection(){
+  return selectedMolecules.begin();
+}
+World::MoleculeSelectionIterator World::endMoleculeSelection(){
+  return selectedMolecules.end();
+}
 …
     Thermostats(new ThermoStatContainer),
     ExitFlag(0),
+    atoms(),
+    atoms(this),
+    selectedAtoms(this),
     currAtomId(0),
+    molecules(),
+    lastAtomPoolSize(0),
+    numAtomDefragSkips(0),
+    molecules(this),
+    selectedMolecules(this),
     currMoleculeId(0),
     molecules_deprecated(new MoleculeListClass(this))
+{
+  cell_size = new double[6];
+  cell_size[0] = 20.;
+  cell_size[1] = 0.;
+  cell_size[2] = 20.;
+  cell_size[3] = 0.;
+  cell_size[4] = 0.;
+  cell_size[5] = 20.;
+  cell_size = new Box;
+  Matrix domain;
+  domain.at(0,0) = 20;
+  domain.at(1,1) = 20;
+  domain.at(2,2) = 20;
+  cell_size->setM(domain);
   defaultName = "none";
   molecules_deprecated->signOn(this);
 …
+{
   molecules_deprecated->signOff(this);
   delete[] cell_size;
+  delete cell_size;
   delete molecules_deprecated;
   delete periode;

src/World.hpp

-              rba9f5b
+              rbe97a8
 #include "Patterns/Cacheable.hpp"
 #include "Patterns/Singleton.hpp"
+#include "Patterns/ObservedContainer.hpp"
 // include config.h
 …
 class AtomDescriptor_impl;
 template<typename T> class AtomsCalculation;
+class Box;
 class config;
 class ManipulateAtomsProcess;
+class Matrix;
 class molecule;
 class MoleculeDescriptor;
 …
 class ThermoStatContainer;
 /****************************************** forward declarations *****************************/
 …
 friend class MoleculeDescriptor_impl;
 friend class MoleculeDescriptor;
+// coupling with descriptors over selection
+friend class AtomSelectionDescriptor_impl;
+friend class MoleculeSelectionDescriptor_impl;
 // Actions, calculations etc associated with the World
 …
   // Types for Atom and Molecule structures
   typedef std::map<atomId_t,atom*> AtomSet;
   typedef std::map<moleculeId_t,molecule*> MoleculeSet;
+  typedef ObservedContainer<std::map<atomId_t,atom*> > AtomSet;
+  typedef ObservedContainer<std::map<moleculeId_t,molecule*> > MoleculeSet;
   /***** getter and setter *****/
 …
    * get the domain size as a symmetric matrix (6 components)
    */
+  double * getDomain();
+  Box& getDomain();
+  /**
+   * Set the domain size from a matrix object
+   *
+   * Matrix needs to be symmetric
+   */
+  void setDomain(const Matrix &mat);
   /**
 …
   ManipulateAtomsProcess* manipulateAtoms(boost::function<void(atom*)>,std::string);
+  /****
+   * Iterators to use internal data structures
+   * All these iterators are observed to track changes.
+   * There is a corresponding protected section with unobserved iterators,
+   * which can be used internally when the extra speed is needed
+   */
+  typedef SelectiveIterator<atom*,AtomSet,AtomDescriptor>       AtomIterator;
+  /**
+   * returns an iterator over all Atoms matching a given descriptor.
+   * This iterator is observed, so don't keep it around unnecessary to
+   * avoid unintended blocking.
+   */
+  AtomIterator getAtomIter(AtomDescriptor descr);
+  AtomIterator getAtomIter();
+  AtomIterator atomEnd();
+  typedef SelectiveIterator<molecule*,MoleculeSet,MoleculeDescriptor>   MoleculeIterator;
+  /**
+   * returns an iterator over all Molecules matching a given descriptor.
+   * This iterator is observed, so don't keep it around unnecessary to
+   * avoid unintended blocking.
+   */
+  MoleculeIterator getMoleculeIter(MoleculeDescriptor descr);
+  MoleculeIterator getMoleculeIter();
+  MoleculeIterator moleculeEnd();
+  /******** Selections of molecules and Atoms *************/
+  void clearAtomSelection();
+  void selectAtom(atom*);
+  void selectAtom(atomId_t);
+  void selectAllAtoms(AtomDescriptor);
+  void selectAtomsOfMolecule(molecule*);
+  void selectAtomsOfMolecule(moleculeId_t);
+  void unselectAtom(atom*);
+  void unselectAtom(atomId_t);
+  void unselectAllAtoms(AtomDescriptor);
+  void unselectAtomsOfMolecule(molecule*);
+  void unselectAtomsOfMolecule(moleculeId_t);
+  void clearMoleculeSelection();
+  void selectMolecule(molecule*);
+  void selectMolecule(moleculeId_t);
+  void selectAllMoleculess(MoleculeDescriptor);
+  void selectMoleculeOfAtom(atom*);
+  void selectMoleculeOfAtom(atomId_t);
+  void unselectMolecule(molecule*);
+  void unselectMolecule(moleculeId_t);
+  void unselectAllMoleculess(MoleculeDescriptor);
+  void unselectMoleculeOfAtom(atom*);
+  void unselectMoleculeOfAtom(atomId_t);
+  /******************** Iterators to selections *****************/
+  typedef AtomSet::iterator AtomSelectionIterator;
+  AtomSelectionIterator beginAtomSelection();
+  AtomSelectionIterator endAtomSelection();
+  typedef MoleculeSet::iterator MoleculeSelectionIterator;
+  MoleculeSelectionIterator beginMoleculeSelection();
+  MoleculeSelectionIterator endMoleculeSelection();
 protected:
+  /**** Iterators to use internal data structures */
+  /****
+   * Iterators to use internal data structures
+   * All these iterators are unobserved for speed reasons.
+   * There is a corresponding public section to these methods,
+   * which produce observed iterators.*/
   // Atoms
   typedef SelectiveIterator<atom*,AtomSet,AtomDescriptor> AtomIterator;
+  typedef SelectiveIterator<atom*,AtomSet::set_t,AtomDescriptor>        internal_AtomIterator;
   /**
 …
    * used for internal purposes, like AtomProcesses and AtomCalculations.
    */
   AtomIterator getAtomIter(AtomDescriptor descr);
+  internal_AtomIterator getAtomIter_internal(AtomDescriptor descr);
   /**
 …
    * used for internal purposes, like AtomProcesses and AtomCalculations.
    */
   AtomIterator atomEnd();
+  internal_AtomIterator atomEnd_internal();
   // Molecules
+  typedef SelectiveIterator<molecule*,MoleculeSet,MoleculeDescriptor> MoleculeIterator;
+  typedef SelectiveIterator<molecule*,MoleculeSet::set_t,MoleculeDescriptor>   internal_MoleculeIterator;
   /**
 …
    * used for internal purposes, like MoleculeProcesses and MoleculeCalculations.
    */
   MoleculeIterator getMoleculeIter(MoleculeDescriptor descr);
+  internal_MoleculeIterator getMoleculeIter_internal(MoleculeDescriptor descr);
   /**
 …
    * used for internal purposes, like MoleculeProcesses and MoleculeCalculations.
    */
   MoleculeIterator moleculeEnd();
+  internal_MoleculeIterator moleculeEnd_internal();
 …
   void releaseAtomId(atomId_t);
   bool reserveAtomId(atomId_t);
+  void defragAtomIdPool();
+  moleculeId_t getNextMoleculeId();
+  void releaseMoleculeId(moleculeId_t);
+  bool reserveMoleculeId(moleculeId_t);
+  void defragMoleculeIdPool();
   periodentafel *periode;
   config *configuration;
   static double *cell_size;
+  Box *cell_size;
   std::string defaultName;
   class ThermoStatContainer *Thermostats;
   int ExitFlag;
+public:
+private:
   AtomSet atoms;
+private:
+  std::set<atomId_t> atomIdPool; //<!stores the pool for all available AtomIds below currAtomId
+  AtomSet selectedAtoms;
+  typedef std::set<std::pair<atomId_t, atomId_t> > atomIdPool_t;
+  /**
+   * stores the pool for all available AtomIds below currAtomId
+   *
+   * The pool contains ranges of free ids in the form [bottom,top).
+   */
+  atomIdPool_t atomIdPool;
   atomId_t currAtomId; //!< stores the next available Id for atoms
+  size_t lastAtomPoolSize; //!< size of the pool after last defrag, to skip some defrags
+  unsigned int numAtomDefragSkips;
   MoleculeSet molecules;
+  MoleculeSet selectedMolecules;
+  typedef std::set<std::pair<moleculeId_t, moleculeId_t> > moleculeIdPool_t;
+  /**
+   * stores the pool for all available AtomIds below currAtomId
+   *
+   * The pool contains ranges of free ids in the form [bottom,top).
+   */
+  moleculeIdPool_t moleculeIdPool;
   moleculeId_t currMoleculeId;
+  size_t lastMoleculePoolSize; //!< size of the pool after last defrag, to skip some defrags
+  unsigned int numMoleculeDefragSkips;
 private:
   /**

src/analysis_bonds.cpp

rba9f5b	rbe97a8
13	13	#include "element.hpp"
14	14	#include "info.hpp"
	15	#include "verbose.hpp"
15	16	#include "log.hpp"
16	17	#include "molecule.hpp"

src/analysis_correlation.cpp

-              rba9f5b
+              rbe97a8
 #include <iostream>
+#include <iomanip>
 #include "analysis_correlation.hpp"
 …
 #include "triangleintersectionlist.hpp"
 #include "vector.hpp"
+#include "Matrix.hpp"
 #include "verbose.hpp"
 #include "World.hpp"
+#include "Box.hpp"
 …
   PairCorrelationMap *outmap = NULL;
   double distance = 0.;
   double *domain = World::getInstance().getDomain();
+  Box &domain = World::getInstance().getDomain();
   if (molecules->ListOfMolecules.empty()) {
 …
                 for (set <pair<element *, element *> >::iterator PairRunner = PairsOfElements.begin(); PairRunner != PairsOfElements.end(); ++PairRunner)
                   if ((PairRunner->first == (**iter).type) && (PairRunner->second == (**runner).type)) {
                     distance = (*iter)->node->PeriodicDistance(*(*runner)->node,  domain);
+                    distance = domain.periodicDistance(*(*iter)->node,*(*runner)->node);
                     //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
                     outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> ((*iter), (*runner)) ) );
 …
   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++){
     if ((*MolWalker)->ActiveFlag) {
       double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
       double * FullInverseMatrix = InverseMatrix(FullMatrix);
+      Matrix FullMatrix = World::getInstance().getDomain().getM();
+      Matrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
       DoeLog(2) && (eLog()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
       eLog() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
       for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
         DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
+        periodicX = *(*iter)->node;
+        periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
+        periodicX = FullInverseMatrix * (*(*iter)->node); // x now in [0,1)^3
         // go through every range in xyz and get distance
         for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
           for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
             for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
+              checkX = Vector(n[0], n[1], n[2]) + periodicX;
+              checkX.MatrixMultiplication(FullMatrix);
+              checkX = FullMatrix * (Vector(n[0], n[1], n[2]) + periodicX);
               for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++){
                 if ((*MolOtherWalker)->ActiveFlag) {
 …
                       for (set <pair<element *, element *> >::iterator PairRunner = PairsOfElements.begin(); PairRunner != PairsOfElements.end(); ++PairRunner)
                         if ((PairRunner->first == (**iter).type) && (PairRunner->second == (**runner).type)) {
+                          periodicOtherX = *(*runner)->node;
+                          periodicOtherX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
+                          periodicOtherX = FullInverseMatrix * (*(*runner)->node); // x now in [0,1)^3
                           // go through every range in xyz and get distance
                           for (Othern[0]=-ranges[0]; Othern[0] <= ranges[0]; Othern[0]++)
                             for (Othern[1]=-ranges[1]; Othern[1] <= ranges[1]; Othern[1]++)
                               for (Othern[2]=-ranges[2]; Othern[2] <= ranges[2]; Othern[2]++) {
+                                checkOtherX = Vector(Othern[0], Othern[1], Othern[2]) + periodicOtherX;
+                                checkOtherX.MatrixMultiplication(FullMatrix);
+                                checkOtherX = FullMatrix * (Vector(Othern[0], Othern[1], Othern[2]) + periodicOtherX);
                                 distance = checkX.distance(checkOtherX);
                                 //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
 …
+            }
+      }
-      delete[](FullMatrix);
-      delete[](FullInverseMatrix);
+    }
+  }
 …
   CorrelationToPointMap *outmap = NULL;
   double distance = 0.;
   double *cell_size = World::getInstance().getDomain();
+  Box &domain = World::getInstance().getDomain();
   if (molecules->ListOfMolecules.empty()) {
 …
         for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
           if ((*type == NULL) || ((*iter)->type == *type)) {
             distance = (*iter)->node->PeriodicDistance(*point, cell_size);
+            distance = domain.periodicDistance(*(*iter)->node,*point);
             DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
             outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> ((*iter), point) ) );
 …
   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
     if ((*MolWalker)->ActiveFlag) {
       double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
       double * FullInverseMatrix = InverseMatrix(FullMatrix);
+      Matrix FullMatrix = World::getInstance().getDomain().getM();
+      Matrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
       DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
       for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
 …
         for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
           if ((*type == NULL) || ((*iter)->type == *type)) {
+            periodicX = *(*iter)->node;
+            periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
+            periodicX = FullInverseMatrix * (*(*iter)->node); // x now in [0,1)^3
             // go through every range in xyz and get distance
             for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
               for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
                 for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
+                  checkX = Vector(n[0], n[1], n[2]) + periodicX;
+                  checkX.MatrixMultiplication(FullMatrix);
+                  checkX = FullMatrix * (Vector(n[0], n[1], n[2]) + periodicX);
                   distance = checkX.distance(*point);
                   DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
 …
+          }
+      }
-      delete[](FullMatrix);
-      delete[](FullInverseMatrix);
+    }
 …
   for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
     if ((*MolWalker)->ActiveFlag) {
       double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
       double * FullInverseMatrix = InverseMatrix(FullMatrix);
+      Matrix FullMatrix = World::getInstance().getDomain().getM();
+      Matrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
       DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
       for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
 …
         for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
           if ((*type == NULL) || ((*iter)->type == *type)) {
+            periodicX = *(*iter)->node;
+            periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
+            periodicX = FullInverseMatrix * (*(*iter)->node); // x now in [0,1)^3
             // go through every range in xyz and get distance
             ShortestDistance = -1.;
 …
               for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
                 for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
+                  checkX = Vector(n[0], n[1], n[2]) + periodicX;
+                  checkX.MatrixMultiplication(FullMatrix);
+                  checkX = FullMatrix * (Vector(n[0], n[1], n[2]) + periodicX);
                   TriangleIntersectionList Intersections(&checkX,Surface,LC);
                   distance = Intersections.GetSmallestDistance();
 …
+          }
+      }
-      delete[](FullMatrix);
-      delete[](FullInverseMatrix);
+    }

src/analysis_correlation.hpp

rba9f5b	rbe97a8
26	26
27	27	#include "atom.hpp"
	28	#include "verbose.hpp"
28	29
29	30	/**************************************** forward declarations ***************************/

src/analyzer.cpp

-              rba9f5b
+              rbe97a8
 #include "datacreator.hpp"
 #include "helpers.hpp"
-#include "memoryallocator.hpp"
 #include "parser.hpp"
 #include "periodentafel.hpp"
+#include "verbose.hpp"
 // include config.h

src/atom.cpp

-              rba9f5b
+              rbe97a8
 #include "element.hpp"
 #include "lists.hpp"
-#include "memoryallocator.hpp"
 #include "parser.hpp"
 #include "vector.hpp"
 #include "World.hpp"
 #include "molecule.hpp"
+#include "Shapes/Shape.hpp"
+#include <iomanip>
 /************************************* Functions for class atom *************************************/
 …
  * \return true - is inside, false - is not
  */
 bool atom::IsInParallelepiped(const Vector offset, const double *parallelepiped) const
+{
   return (node->IsInParallelepiped(offset, parallelepiped));
+bool atom::IsInShape(const Shape& shape) const
+{
+  return shape.isInside(*node);
 };
 …
+}
+molecule* atom::getMolecule(){
+  return mol;
+}
 void atom::removeFromMolecule(){
   if(mol){

src/atom.hpp

-              rba9f5b
+              rbe97a8
 #endif
 #include <iostream>
+#include <iosfwd>
 #include <list>
 #include <vector>
 …
 class World;
 class molecule;
+class Shape;
 /********************************************** declarations *******************************/
 …
   double DistanceToVector(const Vector &origin) const;
   double DistanceSquaredToVector(const Vector &origin) const;
   bool IsInParallelepiped(const Vector offset, const double *parallelepiped) const;
+  bool IsInShape(const Shape&) const;
   // getter and setter
 …
    void setMolecule(molecule*);
+   molecule* getMolecule();
    void removeFromMolecule();

src/atom_graphnodeinfo.cpp

-              rba9f5b
+              rbe97a8
 #include "atom_graphnodeinfo.hpp"
-#include "memoryallocator.hpp"
 /** Constructor of class GraphNodeInfo.
  */
 GraphNodeInfo::GraphNodeInfo() : GraphNr(-1), ComponentNr(NULL), LowpointNr(-1), SeparationVertex(false), IsCyclic(false), Ancestor(NULL) {};
+GraphNodeInfo::GraphNodeInfo() : GraphNr(-1), ComponentNr(0), LowpointNr(-1), SeparationVertex(false), IsCyclic(false), Ancestor(0) {};
 /** Destructor of class GraphNodeInfo.

src/atom_particleinfo.cpp

-              rba9f5b
+              rbe97a8
 #include "atom_particleinfo.hpp"
+#include "memoryallocator.hpp"
+#include <iostream>
 /** Constructor of ParticleInfo.

src/atom_particleinfo.hpp

-              rba9f5b
+              rbe97a8
 #endif
+#include <iostream>
+#include <iosfwd>
+#include <string>
 /****************************************** forward declarations *****************************/

src/atom_trajectoryparticle.hpp

rba9f5b	rbe97a8
20	20	#include <fstream>
21	21
	22	#include <gsl/gsl_inline.h>
22	23	#include <gsl/gsl_randist.h>
23	24

src/bond.cpp

rba9f5b	rbe97a8
7	7	#include "Helpers/MemDebug.hpp"
8	8
	9	#include "verbose.hpp"
9	10	#include "atom.hpp"
10	11	#include "bond.hpp"

src/bondgraph.cpp

rba9f5b	rbe97a8
15	15	#include "element.hpp"
16	16	#include "info.hpp"
	17	#include "verbose.hpp"
17	18	#include "log.hpp"
18	19	#include "molecule.hpp"

src/bondgraph.hpp

rba9f5b	rbe97a8
18	18	#endif
19	19
20		#include <ios~~tream~~>
	20	#include <iosfwd>
21	21
22	22	/********************************************* defines *********************************/

src/boundary.cpp

-              rba9f5b
+              rbe97a8
 #include "info.hpp"
 #include "linkedcell.hpp"
+#include "verbose.hpp"
 #include "log.hpp"
-#include "memoryallocator.hpp"
 #include "molecule.hpp"
 #include "tesselation.hpp"
 …
 #include "World.hpp"
 #include "Plane.hpp"
+#include "Matrix.hpp"
+#include "Box.hpp"
+#include <iostream>
+#include <iomanip>
 #include<gsl/gsl_poly.h>
 …
   int N[NDIM];
   int n[NDIM];
   double *M =  ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
   double Rotations[NDIM*NDIM];
   double *MInverse = InverseMatrix(M);
+  const Matrix &M = World::getInstance().getDomain().getM();
+  Matrix Rotations;
+  const Matrix &MInverse = World::getInstance().getDomain().getMinv();
   Vector AtomTranslations;
   Vector FillerTranslations;
 …
   // calculate filler grid in [0,1]^3
+  FillerDistance = Vector(distance[0], distance[1], distance[2]);
+  FillerDistance.InverseMatrixMultiplication(M);
+  FillerDistance = MInverse * Vector(distance[0], distance[1], distance[2]);
   for(int i=0;i<NDIM;i++)
     N[i] = (int) ceil(1./FillerDistance[i]);
 …
       for (n[2] = 0; n[2] < N[2]; n[2]++) {
         // calculate position of current grid vector in untransformed box
+        CurrentPosition = Vector((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
+        CurrentPosition.MatrixMultiplication(M);
+        CurrentPosition = M * Vector((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
         // create molecule random translation vector ...
         for (int i=0;i<NDIM;i++)
 …
+            }
             Rotations[0] =   cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2]));
             Rotations[3] =   sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2]));
             Rotations[6] =               cos(phi[1])*sin(phi[2])                                     ;
             Rotations[1] = - sin(phi[0])*cos(phi[1])                                                ;
             Rotations[4] =   cos(phi[0])*cos(phi[1])                                                ;
             Rotations[7] =               sin(phi[1])                                                ;
             Rotations[3] = - cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2]));
             Rotations[5] = - sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2]));
             Rotations[8] =               cos(phi[1])*cos(phi[2])                                     ;
+            Rotations.set(0,0,  cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2])));
+            Rotations.set(0,1,  sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2])));
+            Rotations.set(0,2,              cos(phi[1])*sin(phi[2])                                        );
+            Rotations.set(1,0, -sin(phi[0])*cos(phi[1])                                                    );
+            Rotations.set(1,1,  cos(phi[0])*cos(phi[1])                                                    );
+            Rotations.set(1,2,              sin(phi[1])                                                    );
+            Rotations.set(2,0, -cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2])));
+            Rotations.set(2,1, -sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2])));
+            Rotations.set(2,2,              cos(phi[1])*cos(phi[2])                                        );
+          }
 …
           Inserter = (*iter)->x;
           if (DoRandomRotation)
             Inserter.MatrixMultiplication(Rotations);
+            Inserter *= Rotations;
           Inserter += AtomTranslations + FillerTranslations + CurrentPosition;
           // check whether inserter is inside box
           Inserter.MatrixMultiplication(MInverse);
+          Inserter *= MInverse;
           FillIt = true;
           for (int i=0;i<NDIM;i++)
             FillIt = FillIt && (Inserter[i] >= -MYEPSILON) && ((Inserter[i]-1.) <= MYEPSILON);
           Inserter.MatrixMultiplication(M);
+          Inserter *= M;
           // Check whether point is in- or outside
 …
     delete TesselStruct[(*ListRunner)];
+  }
-  delete[](M);
-  delete[](MInverse);
   return Filling;

src/boundary.hpp

rba9f5b	rbe97a8
12	12
13	13	#include <fstream>
14		#include <ios~~tream~~>
	14	#include <iosfwd>
15	15
16	16	// STL headers

src/builder.cpp

rba9f5b	rbe97a8
60	60	#include "config.hpp"
61	61	#include "log.hpp"
62		~~#include "memoryusageobserver.hpp"~~
63	62	#include "molecule.hpp"
64	63	#include "periodentafel.hpp"

src/config.cpp

-              rba9f5b
+              rbe97a8
 #include "info.hpp"
 #include "lists.hpp"
+#include "verbose.hpp"
 #include "log.hpp"
 #include "molecule.hpp"
-#include "memoryallocator.hpp"
 #include "molecule.hpp"
 #include "periodentafel.hpp"
 #include "ThermoStatContainer.hpp"
 #include "World.hpp"
+#include "Matrix.hpp"
+#include "Box.hpp"
 /************************************* Functions for class config ***************************/
 …
   // Unit cell and magnetic field
   ParseForParameter(verbose,FileBuffer, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
   double * const cell_size = World::getInstance().getDomain();
+  double * cell_size = new double[6];
   cell_size[0] = BoxLength[0];
   cell_size[1] = BoxLength[3];
 …
   cell_size[4] = BoxLength[7];
   cell_size[5] = BoxLength[8];
+  World::getInstance().setDomain(cell_size);
+  delete cell_size;
   //if (1) fprintf(stderr,"\n");
 …
   ParseForParameter(verbose,file, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
   double * const cell_size = World::getInstance().getDomain();
+  double * cell_size = new double[6];
   cell_size[0] = BoxLength[0];
   cell_size[1] = BoxLength[3];
 …
   cell_size[4] = BoxLength[7];
   cell_size[5] = BoxLength[8];
+  World::getInstance().setDomain(cell_size);
+  delete[] cell_size;
   if (1) fprintf(stderr,"\n");
   config::DoPerturbation = 0;
 …
   // bring MaxTypes up to date
   mol->CountElements();
   const double * const cell_size = World::getInstance().getDomain();
+  const Matrix &domain = World::getInstance().getDomain().getM();
   ofstream * const output = new ofstream(filename, ios::out);
   if (output != NULL) {
 …
     *output << endl;
     *output << "BoxLength\t\t\t# (Length of a unit cell)" << endl;
     *output << cell_size[0] << "\t" << endl;
     *output << cell_size[1] << "\t" << cell_size[2] << "\t" << endl;
     *output << cell_size[3] << "\t" << cell_size[4] << "\t" << cell_size[5] << "\t" << endl;
+    *output << domain.at(0,0) << "\t" << endl;
+    *output << domain.at(1,0) << "\t" << domain.at(1,1) << "\t" << endl;
+    *output << domain.at(2,0) << "\t" << domain.at(2,1) << "\t" << domain.at(2,2) << "\t" << endl;
     // FIXME
     *output << endl;

src/datacreator.cpp

-              rba9f5b
+              rbe97a8
 #include "helpers.hpp"
 #include "parser.hpp"
+#include "verbose.hpp"
+#include <iomanip>
 //=========================== FUNCTIONS============================

src/datacreator.hpp

rba9f5b	rbe97a8
10	10	using namespace std;
11	11
12		#include <ios~~tream~~>
	12	#include <iosfwd>
13	13
14	14	/**************************************** forward declarations ***************************/

src/defs.hpp

rba9f5b	rbe97a8
84	84	#define MOLECUILDER_NAME "Molecuilder"
85	85
	86	const unsigned int MAX_POOL_FRAGMENTATION=20;
	87	const unsigned int MAX_FRAGMENTATION_SKIPS=100;
	88
86	89	#endif /DEFS_HPP_/

src/element.hpp

rba9f5b	rbe97a8
16	16	#endif
17	17
18		#include <ios~~tream~~>
	18	#include <iosfwd>
19	19	#include <string>
20	20

src/ellipsoid.cpp

-              rba9f5b
+              rbe97a8
 #include "tesselation.hpp"
 #include "vector.hpp"
+#include "Matrix.hpp"
 #include "verbose.hpp"
 …
   Vector helper, RefPoint;
   double distance = -1.;
   double Matrix[NDIM*NDIM];
+  Matrix Matrix;
   double InverseLength[3];
   double psi,theta,phi; // euler angles in ZX'Z'' convention
 …
   theta = EllipsoidAngle[1];
   phi = EllipsoidAngle[2];
   Matrix[0] = cos(psi)*cos(phi) - sin(psi)*cos(theta)*sin(phi);
   Matrix[1] = -cos(psi)*sin(phi) - sin(psi)*cos(theta)*cos(phi);
   Matrix[2] = sin(psi)*sin(theta);
   Matrix[3] = sin(psi)*cos(phi) + cos(psi)*cos(theta)*sin(phi);
   Matrix[4] = cos(psi)*cos(theta)*cos(phi) - sin(psi)*sin(phi);
   Matrix[5] = -cos(psi)*sin(theta);
   Matrix[6] = sin(theta)*sin(phi);
   Matrix[7] = sin(theta)*cos(phi);
   Matrix[8] = cos(theta);
   helper.MatrixMultiplication(Matrix);
+  Matrix.set(0,0, cos(psi)*cos(phi) - sin(psi)*cos(theta)*sin(phi));
+  Matrix.set(1,0, -cos(psi)*sin(phi) - sin(psi)*cos(theta)*cos(phi));
+  Matrix.set(2,0, sin(psi)*sin(theta));
+  Matrix.set(0,1, sin(psi)*cos(phi) + cos(psi)*cos(theta)*sin(phi));
+  Matrix.set(1,1, cos(psi)*cos(theta)*cos(phi) - sin(psi)*sin(phi));
+  Matrix.set(2,1, -cos(psi)*sin(theta));
+  Matrix.set(0,2, sin(theta)*sin(phi));
+  Matrix.set(1,2, sin(theta)*cos(phi));
+  Matrix.set(2,2, cos(theta));
+  helper *= Matrix;
   helper.ScaleAll(InverseLength);
   //Log() << Verbose(4) << "Transformed RefPoint is at " << helper << "." << endl;
 …
   phi = -EllipsoidAngle[2];
   helper.ScaleAll(EllipsoidLength);
   Matrix[0] = cos(psi)*cos(phi) - sin(psi)*cos(theta)*sin(phi);
   Matrix[1] = -cos(psi)*sin(phi) - sin(psi)*cos(theta)*cos(phi);
   Matrix[2] = sin(psi)*sin(theta);
   Matrix[3] = sin(psi)*cos(phi) + cos(psi)*cos(theta)*sin(phi);
   Matrix[4] = cos(psi)*cos(theta)*cos(phi) - sin(psi)*sin(phi);
   Matrix[5] = -cos(psi)*sin(theta);
   Matrix[6] = sin(theta)*sin(phi);
   Matrix[7] = sin(theta)*cos(phi);
   Matrix[8] = cos(theta);
   helper.MatrixMultiplication(Matrix);
+  Matrix.set(0,0, cos(psi)*cos(phi) - sin(psi)*cos(theta)*sin(phi));
+  Matrix.set(1,0, -cos(psi)*sin(phi) - sin(psi)*cos(theta)*cos(phi));
+  Matrix.set(2,0, sin(psi)*sin(theta));
+  Matrix.set(0,1, sin(psi)*cos(phi) + cos(psi)*cos(theta)*sin(phi));
+  Matrix.set(1,1, cos(psi)*cos(theta)*cos(phi) - sin(psi)*sin(phi));
+  Matrix.set(2,1, -cos(psi)*sin(theta));
+  Matrix.set(0,2, sin(theta)*sin(phi));
+  Matrix.set(1,2, sin(theta)*cos(phi));
+  Matrix.set(2,2, cos(theta));
+  helper *= Matrix;
   //Log() << Verbose(4) << "Intersection is at " << helper << "." << endl;

src/errorlogger.cpp

rba9f5b	rbe97a8
9	9
10	10	#include <fstream>
	11	#include <iostream>
11	12	#include "errorlogger.hpp"
12	13	#include "verbose.hpp"

src/errorlogger.hpp

rba9f5b	rbe97a8
9	9	#define ERRORLOGGER_HPP_
10	10
11		#include <ios~~tream~~>
	11	#include <iosfwd>
12	12
13	13	#include "Patterns/Singleton.hpp"

src/graph.cpp

rba9f5b	rbe97a8
13	13	#include "config.hpp"
14	14	#include "graph.hpp"
	15	#include "verbose.hpp"
15	16	#include "log.hpp"
16	17	#include "molecule.hpp"

src/graph.hpp

-              rba9f5b
+              rbe97a8
 class molecule;
-class Graph;
 class SubGraph;
 class Node;
 …
 /********************************************** definitions *********************************/
+#define NodeMap pair < int, class Node* >
+#define EdgeMap multimap < class Node*, class Edge* >
+typedef std::pair < int, class Node* > NodeMap;
+typedef std::multimap < class Node*, class Edge* > EdgeMap;
+#define KeyStack deque<int>
+#define KeySet set<int>
+#define NumberValuePair pair<int, double>
+#define Graph map <KeySet, NumberValuePair, KeyCompare >
+#define GraphPair pair <KeySet, NumberValuePair >
+#define KeySetTestPair pair<KeySet::iterator, bool>
+#define GraphTestPair pair<Graph::iterator, bool>
+typedef std::deque<int> KeyStack;
+typedef std::set<int> KeySet;
+typedef std::pair<int, double> NumberValuePair;
+/******************************** Some small functions and/or structures **********************************/
+// needed for definition of Graph and GraphTestPair
 struct KeyCompare
+{
   bool operator() (const KeySet SubgraphA, const KeySet SubgraphB) const;
 };
+typedef std::map <KeySet, NumberValuePair, KeyCompare > Graph;
+typedef std::pair <KeySet, NumberValuePair > GraphPair;
+typedef std::pair<KeySet::iterator, bool> KeySetTestPair;
+typedef std::pair<Graph::iterator, bool> GraphTestPair;
+/******************************** Some small functions and/or structures **********************************/
 //bool operator < (KeySet SubgraphA, KeySet SubgraphB);   //note: this declaration is important, otherwise normal < is used (producing wrong order)

src/gslvector.cpp

-              rba9f5b
+              rbe97a8
 #include <cassert>
 #include <cmath>
+#include <iostream>
 #include "gslvector.hpp"
 #include "defs.hpp"
 #include "vector.hpp"
+#include "VectorContent.hpp"
 /** Constructor of class GSLVector.
 …
  */
 void GSLVector::SetFromVector(Vector &v){
   gsl_vector_memcpy (vector, v.get());
+  gsl_vector_memcpy (vector, v.get()->content);
+}

src/gslvector.hpp

rba9f5b	rbe97a8
18	18	#endif
19	19
20		#include <ios~~tream~~>
	20	#include <iosfwd>
21	21	#include <gsl/gsl_vector.h>
22	22

src/helpers.cpp

-              rba9f5b
+              rbe97a8
 #include "helpers.hpp"
 #include "Helpers/fast_functions.hpp"
+#include "verbose.hpp"
 #include "log.hpp"
+#include "memoryusageobserver.hpp"
+#include <iostream>
 /********************************************** helpful functions *********************************/
 …
 };
-/** Blows the 6-dimensional \a cell_size array up to a full NDIM by NDIM matrix.
- * \param *symm 6-dim array of unique symmetric matrix components
- * \return allocated NDIM*NDIM array with the symmetric matrix
- */
-double * ReturnFullMatrixforSymmetric(const double * const symm)
+{
-  double *matrix = new double[NDIM * NDIM];
-  matrix[0] = symm[0];
-  matrix[1] = symm[1];
-  matrix[2] = symm[3];
-  matrix[3] = symm[1];
-  matrix[4] = symm[2];
-  matrix[5] = symm[4];
-  matrix[6] = symm[3];
-  matrix[7] = symm[4];
-  matrix[8] = symm[5];
-  return matrix;
-};
-/** Calculate the inverse of a 3x3 matrix.
- * \param *matrix NDIM_NDIM array
- */
-double * InverseMatrix( const double * const A)
+{
-  double *B = new double[NDIM * NDIM];
-  double detA = RDET3(A);
-  double detAReci;
-  for (int i=0;i<NDIM*NDIM;++i)
-    B[i] = 0.;
-  // calculate the inverse B
-  if (fabs(detA) > MYEPSILON) {;  // RDET3(A) yields precisely zero if A irregular
-    detAReci = 1./detA;
-    B[0] =  detAReci*RDET2(A[4],A[5],A[7],A[8]);    // A_11
-    B[1] = -detAReci*RDET2(A[1],A[2],A[7],A[8]);    // A_12
-    B[2] =  detAReci*RDET2(A[1],A[2],A[4],A[5]);    // A_13
-    B[3] = -detAReci*RDET2(A[3],A[5],A[6],A[8]);    // A_21
-    B[4] =  detAReci*RDET2(A[0],A[2],A[6],A[8]);    // A_22
-    B[5] = -detAReci*RDET2(A[0],A[2],A[3],A[5]);    // A_23
-    B[6] =  detAReci*RDET2(A[3],A[4],A[6],A[7]);    // A_31
-    B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]);    // A_32
-    B[8] =  detAReci*RDET2(A[0],A[1],A[3],A[4]);    // A_33
+  }
-  return B;
-};
 /** Comparison function for GSL heapsort on distances in two molecules.
  * \param *a

src/helpers.hpp

-              rba9f5b
+              rbe97a8
 #include "defs.hpp"
 #include "log.hpp"
-#include "memoryallocator.hpp"
 /********************************************** definitions *********************************/
 …
 bool IsValidNumber( const char *string);
 int CompareDoubles (const void * a, const void * b);
-double * ReturnFullMatrixforSymmetric(const double * const cell_size);
-double * InverseMatrix(const double * const A);
 void performCriticalExit();

src/joiner.cpp

-              rba9f5b
+              rbe97a8
 #include "datacreator.hpp"
 #include "helpers.hpp"
-#include "memoryallocator.hpp"
 #include "parser.hpp"
 #include "periodentafel.hpp"
+#include "verbose.hpp"
 //============================== MAIN =============================

src/linkedcell.cpp

rba9f5b	rbe97a8
10	10	#include "helpers.hpp"
11	11	#include "linkedcell.hpp"
	12	#include "verbose.hpp"
12	13	#include "log.hpp"
13	14	#include "molecule.hpp"

src/logger.cpp

rba9f5b	rbe97a8
9	9
10	10	#include <fstream>
	11	#include <iostream>
11	12	#include "logger.hpp"
12	13	#include "verbose.hpp"

src/logger.hpp

rba9f5b	rbe97a8
9	9	#define LOGGER_HPP_
10	10
11		#include <ios~~tream~~>
	11	#include <iosfwd>
12	12
13	13	#include "Patterns/Singleton.hpp"

src/molecule.cpp

-              rba9f5b
+              rbe97a8
  */
+#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"
 …
 #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"
 …
   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;
 …
       } // (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: ";
 …
       break;
+  }
-  delete[](matrix);
 //  Log() << Verbose(3) << "End of AddHydrogenReplacementAtom." << endl;
 …
  * @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);
 …
 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);
 };
 …
 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;

src/molecule.hpp

-              rba9f5b
+              rbe97a8
 #define MOLECULES_HPP_
-using namespace std;
 /*********************************************** includes ***********************************/
+// GSL headers
+#include <gsl/gsl_eigen.h>
+#include <gsl/gsl_heapsort.h>
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_multimin.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_randist.h>
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
 //// STL headers
 …
 #include "types.hpp"
 #include "graph.hpp"
-#include "stackclass.hpp"
 #include "tesselation.hpp"
 #include "Patterns/Observer.hpp"
 …
 class periodentafel;
 class Vector;
+class Shape;
+template <class> class StackClass;
 /******************************** Some definitions for easier reading **********************************/
 …
   molecule *CopyMolecule();
   molecule* CopyMoleculeFromSubRegion(const Vector offset, const double *parallelepiped) const;
+  molecule* CopyMoleculeFromSubRegion(const Shape&) const;
   /// Fragment molecule by two different approaches:

src/molecule_dynamics.cpp

-              rba9f5b
+              rbe97a8
 #include "element.hpp"
 #include "info.hpp"
+#include "verbose.hpp"
 #include "log.hpp"
-#include "memoryallocator.hpp"
 #include "molecule.hpp"
 #include "parser.hpp"
 #include "Plane.hpp"
 #include "ThermoStatContainer.hpp"
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_linalg.h>
 /************************************* Functions for class molecule *********************************/

src/molecule_fragmentation.cpp

-              rba9f5b
+              rbe97a8
 #include "helpers.hpp"
 #include "lists.hpp"
+#include "verbose.hpp"
 #include "log.hpp"
-#include "memoryallocator.hpp"
 #include "molecule.hpp"
 #include "periodentafel.hpp"
 #include "World.hpp"
+#include "Matrix.hpp"
+#include "Box.hpp"
+#include "stackclass.hpp"
 /************************************* Functions for class molecule *********************************/
 …
   atom *Walker = NULL;
   atom *OtherWalker = NULL;
+  double * const cell_size = World::getInstance().getDomain();
+  double *matrix = ReturnFullMatrixforSymmetric(cell_size);
+  Matrix matrix = World::getInstance().getDomain().getM();
   enum Shading *ColorList = NULL;
   double tmp;
 …
           Translationvector[i] = (tmp < 0) ? +1. : -1.;
+      }
       Translationvector.MatrixMultiplication(matrix);
+      Translationvector *= matrix;
       //Log() << Verbose(3) << "Translation vector is ";
       Log() << Verbose(0) << Translationvector <<  endl;
 …
   delete(AtomStack);
   delete[](ColorList);
-  delete[](matrix);
   DoLog(2) && (Log() << Verbose(2) << "End of ScanForPeriodicCorrection." << endl);
 };

src/molecule_geometry.cpp

-              rba9f5b
+              rbe97a8
  *      Author: heber
  */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
 #include "Helpers/MemDebug.hpp"
 …
 #include "helpers.hpp"
 #include "leastsquaremin.hpp"
+#include "verbose.hpp"
 #include "log.hpp"
-#include "memoryallocator.hpp"
 #include "molecule.hpp"
 #include "World.hpp"
 #include "Plane.hpp"
+#include "Matrix.hpp"
+#include "Box.hpp"
 #include <boost/foreach.hpp>
+#include <gsl/gsl_eigen.h>
+#include <gsl/gsl_multimin.h>
 …
   const Vector *Center = DetermineCenterOfAll();
   const Vector *CenterBox = DetermineCenterOfBox();
+  double * const cell_size = World::getInstance().getDomain();
+  double *M = ReturnFullMatrixforSymmetric(cell_size);
+  double *Minv = InverseMatrix(M);
+  Box &domain = World::getInstance().getDomain();
   // go through all atoms
   ActOnAllVectors( &Vector::SubtractVector, *Center);
   ActOnAllVectors( &Vector::SubtractVector, *CenterBox);
   ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
   delete[](M);
+  delete[](Minv);
+  BOOST_FOREACH(atom* iter, atoms){
+    *iter->node = domain.WrapPeriodically(*iter->node);
+  }
   delete(Center);
   delete(CenterBox);
 …
+{
   bool status = true;
+  double * const cell_size = World::getInstance().getDomain();
+  double *M = ReturnFullMatrixforSymmetric(cell_size);
+  double *Minv = InverseMatrix(M);
+  Box &domain = World::getInstance().getDomain();
   // go through all atoms
   ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
   delete[](M);
+  delete[](Minv);
+  BOOST_FOREACH(atom* iter, atoms){
+    *iter->node = domain.WrapPeriodically(*iter->node);
+  }
   return status;
 };
 …
+{
   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);
+  const Matrix &M = World::getInstance().getDomain().getM();
+  (*a) *= M;
   return a;
 };
 …
 void molecule::TranslatePeriodically(const Vector *trans)
+{
+  double * const cell_size = World::getInstance().getDomain();
+  double *M = ReturnFullMatrixforSymmetric(cell_size);
+  double *Minv = InverseMatrix(M);
+  Box &domain = World::getInstance().getDomain();
   // go through all atoms
   ActOnAllVectors( &Vector::AddVector, *trans);
   ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
   delete[](M);
+  delete[](Minv);
+  BOOST_FOREACH(atom* iter, atoms){
+    *iter->node = domain.WrapPeriodically(*iter->node);
+  }
 };
 …
   OBSERVE;
   Plane p(*n,0);
   BOOST_FOREACH( atom* iter, atoms ){
+  BOOST_FOREACH(atom* iter, atoms ){
     (*iter->node) = p.mirrorVector(*iter->node);
+  }
 …
 void molecule::DeterminePeriodicCenter(Vector &center)
+{
+  double * const cell_size = World::getInstance().getDomain();
+  double *matrix = ReturnFullMatrixforSymmetric(cell_size);
+  double *inversematrix = InverseMatrix(matrix);
+  const Matrix &matrix = World::getInstance().getDomain().getM();
+  const Matrix &inversematrix = World::getInstance().getDomain().getM();
   double tmp;
   bool flag;
 …
       if ((*iter)->type->Z != 1) {
 #endif
+        Testvector = (*iter)->x;
+        Testvector.MatrixMultiplication(inversematrix);
+        Testvector = inversematrix * (*iter)->x;
         Translationvector.Zero();
         for (BondList::const_iterator Runner = (*iter)->ListOfBonds.begin(); Runner != (*iter)->ListOfBonds.end(); (++Runner)) {
 …
+        }
         Testvector += Translationvector;
         Testvector.MatrixMultiplication(matrix);
+        Testvector *= matrix;
         Center += Testvector;
         Log() << Verbose(1) << "vector is: " << Testvector << endl;
 …
         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 = inversematrix * (*Runner)->GetOtherAtom((*iter))->x;
             Testvector += Translationvector;
             Testvector.MatrixMultiplication(matrix);
+            Testvector *= matrix;
             Center += Testvector;
             Log() << Verbose(1) << "Hydrogen vector is: " << Testvector << endl;
 …
+    }
   } while (!flag);
-  delete[](matrix);
-  delete[](inversematrix);
   Center.Scale(1./static_cast<double>(getAtomCount()));
 …
     DoLog(1) && (Log() << Verbose(1) << "Transforming molecule into PAS ... ");
     // the eigenvectors specify the transformation matrix
+    ActOnAllVectors( &Vector::MatrixMultiplication, (const double *) evec->data );
+    Matrix M = Matrix(evec->data);
+    BOOST_FOREACH(atom* iter, atoms){
+      (*iter->node) *= M;
+    }
     DoLog(0) && (Log() << Verbose(0) << "done." << endl);

src/molecule_graph.cpp

-              rba9f5b
+              rbe97a8
  *      Author: heber
  */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
 #include "Helpers/MemDebug.hpp"
 …
 #include "linkedcell.hpp"
 #include "lists.hpp"
+#include "verbose.hpp"
 #include "log.hpp"
-#include "memoryallocator.hpp"
 #include "molecule.hpp"
 #include "World.hpp"
 #include "Helpers/fast_functions.hpp"
 #include "Helpers/Assert.hpp"
+#include "Matrix.hpp"
+#include "Box.hpp"
+#include "stackclass.hpp"
 struct BFSAccounting
 …
   LinkedCell *LC = NULL;
   bool free_BG = false;
   double * const cell_size = World::getInstance().getDomain();
+  Box &domain = World::getInstance().getDomain();
   if (BG == NULL) {
 …
                           //Log() << Verbose(1) << "Checking distance " << OtherWalker->x.PeriodicDistanceSquared(&(Walker->x), cell_size) << " against typical bond length of " << bonddistance*bonddistance << "." << endl;
                           (BG->*minmaxdistance)(Walker, OtherWalker, MinDistance, MaxDistance, IsAngstroem);
                           const double distance = OtherWalker->x.PeriodicDistanceSquared(Walker->x,cell_size);
+                          const double distance = domain.periodicDistanceSquared(OtherWalker->x,Walker->x);
                           const bool status = (distance <= MaxDistance * MaxDistance) && (distance >= MinDistance * MinDistance);
 //                          Log() << Verbose(1) << "MinDistance is " << MinDistance << " and MaxDistance is " << MaxDistance << "." << endl;
 …
     // From OldGraphNr to CurrentGraphNr ranges an disconnected subgraph
     DoLog(0) && (Log() << Verbose(0) << "Disconnected subgraph ranges from " << OldGraphNr << " to " << DFS.CurrentGraphNr << "." << endl);
     LeafWalker->Leaf->Output((ofstream *)&cout);
+    LeafWalker->Leaf->Output((ofstream *)&(Log() << Verbose(0)));
     DoLog(0) && (Log() << Verbose(0) << endl);

src/molecule_pointcloud.cpp

rba9f5b	rbe97a8
11	11	#include "config.hpp"
12	12	#include "info.hpp"
13		~~#include "memoryallocator.hpp"~~
14	13	#include "molecule.hpp"
15	14

src/moleculelist.cpp

-              rba9f5b
+              rbe97a8
  */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
 #include "Helpers/MemDebug.hpp"
 #include <cstring>
+#include <gsl/gsl_inline.h>
+#include <gsl/gsl_heapsort.h>
 #include "World.hpp"
 …
 #include "linkedcell.hpp"
 #include "lists.hpp"
+#include "verbose.hpp"
 #include "log.hpp"
 #include "molecule.hpp"
-#include "memoryallocator.hpp"
 #include "periodentafel.hpp"
 #include "Helpers/Assert.hpp"
+#include "Matrix.hpp"
+#include "Box.hpp"
+#include "stackclass.hpp"
 #include "Helpers/Assert.hpp"
 …
   int FragmentCounter = 0;
   ofstream output;
+  double cell_size_backup[6];
+  double * const cell_size = World::getInstance().getDomain();
+  // backup cell_size
+  for (int i=0;i<6;i++)
+    cell_size_backup[i] = cell_size[i];
+  Matrix cell_size = World::getInstance().getDomain().getM();
+  Matrix cell_size_backup = cell_size;
   // store the fragments as config and as xyz
   for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
 …
     (*ListRunner)->CenterEdge(&BoxDimension);
     (*ListRunner)->SetBoxDimension(&BoxDimension); // update Box of atoms by boundary
-    int j = -1;
     for (int k = 0; k < NDIM; k++) {
-      j += k + 1;
       BoxDimension[k] = 2.5 * (World::getInstance().getConfig()->GetIsAngstroem() ? 1. : 1. / AtomicLengthToAngstroem);
+      cell_size[j] = BoxDimension[k] * 2.;
+    }
+      cell_size.at(k,k) = BoxDimension[k] * 2.;
+    }
+    World::getInstance().setDomain(cell_size);
     (*ListRunner)->Translate(&BoxDimension);
 …
   // restore cell_size
+  for (int i=0;i<6;i++)
+    cell_size[i] = cell_size_backup[i];
+  World::getInstance().setDomain(cell_size_backup);
   return result;
 …
   // center at set box dimensions
   mol->CenterEdge(&center);
+  World::getInstance().getDomain()[0] = center[0];
+  World::getInstance().getDomain()[1] = 0;
+  World::getInstance().getDomain()[2] = center[1];
+  World::getInstance().getDomain()[3] = 0;
+  World::getInstance().getDomain()[4] = 0;
+  World::getInstance().getDomain()[5] = center[2];
+  Matrix domain;
+  for(int i =0;i<NDIM;++i)
+    domain.at(i,i) = center[i];
+  World::getInstance().setDomain(domain);
   insert(mol);
+}

src/parser.cpp

-              rba9f5b
+              rbe97a8
 #include "helpers.hpp"
-#include "memoryallocator.hpp"
 #include "parser.hpp"
+#include "verbose.hpp"
 // include config.h

src/periodentafel.hpp

-              rba9f5b
+              rbe97a8
 #endif
 #include <iostream>
+#include <iosfwd>
 #include <map>
-#include <iterator>
 #include "unittests/periodentafelTest.hpp"

src/stackclass.hpp

rba9f5b	rbe97a8
12	12
13	13	#include "verbose.hpp"
14		#include "~~memoryallocator~~.hpp"
	14	#include "log.hpp"
15	15
16	16	/**************************************** forward declarations ***************************/

src/tesselation.cpp

rba9f5b	rbe97a8
9	9
10	10	#include <fstream>
	11	#include <iomanip>
11	12
12	13	#include "helpers.hpp"

src/tesselationhelpers.cpp

-              rba9f5b
+              rbe97a8
 #include "verbose.hpp"
 #include "Plane.hpp"
+double DetGet(gsl_matrix * const A, const int inPlace)
+{
+        Info FunctionInfo(__func__);
+  /*
+  inPlace = 1 => A is replaced with the LU decomposed copy.
+  inPlace = 0 => A is retained, and a copy is used for LU.
+  */
+  double det;
+  int signum;
+  gsl_permutation *p = gsl_permutation_alloc(A->size1);
+  gsl_matrix *tmpA=0;
+  if (inPlace)
+  tmpA = A;
+  else {
+  gsl_matrix *tmpA = gsl_matrix_alloc(A->size1, A->size2);
+  gsl_matrix_memcpy(tmpA , A);
+  }
+  gsl_linalg_LU_decomp(tmpA , p , &signum);
+  det = gsl_linalg_LU_det(tmpA , signum);
+  gsl_permutation_free(p);
+  if (! inPlace)
+  gsl_matrix_free(tmpA);
+  return det;
+};
+#include "Matrix.hpp"
 void GetSphere(Vector * const center, const Vector &a, const Vector &b, const Vector &c, const double RADIUS)
+{
         Info FunctionInfo(__func__);
   gsl_matrix *A = gsl_matrix_calloc(3,3);
+  Matrix mat;
   double m11, m12, m13, m14;
   for(int i=0;i<3;i++) {
     gsl_matrix_set(A, i, 0, a[i]);
     gsl_matrix_set(A, i, 1, b[i]);
     gsl_matrix_set(A, i, 2, c[i]);
+  }
   m11 = DetGet(A, 1);
+    mat.set(i, 0, a[i]);
+    mat.set(i, 1, b[i]);
+    mat.set(i, 2, c[i]);
+  }
+  m11 = mat.determinant();
   for(int i=0;i<3;i++) {
     gsl_matrix_set(A, i, 0, a[i]*a[i] + b[i]*b[i] + c[i]*c[i]);
     gsl_matrix_set(A, i, 1, b[i]);
     gsl_matrix_set(A, i, 2, c[i]);
+  }
   m12 = DetGet(A, 1);
+    mat.set(i, 0, a[i]*a[i] + b[i]*b[i] + c[i]*c[i]);
+    mat.set(i, 1, b[i]);
+    mat.set(i, 2, c[i]);
+  }
+  m12 = mat.determinant();
   for(int i=0;i<3;i++) {
     gsl_matrix_set(A, i, 0, a[i]*a[i] + b[i]*b[i] + c[i]*c[i]);
     gsl_matrix_set(A, i, 1, a[i]);
     gsl_matrix_set(A, i, 2, c[i]);
+  }
   m13 = DetGet(A, 1);
+    mat.set(i, 0, a[i]*a[i] + b[i]*b[i] + c[i]*c[i]);
+    mat.set(i, 1, a[i]);
+    mat.set(i, 2, c[i]);
+  }
+  m13 = mat.determinant();
   for(int i=0;i<3;i++) {
     gsl_matrix_set(A, i, 0, a[i]*a[i] + b[i]*b[i] + c[i]*c[i]);
     gsl_matrix_set(A, i, 1, a[i]);
     gsl_matrix_set(A, i, 2, b[i]);
+  }
   m14 = DetGet(A, 1);
+    mat.set(i, 0, a[i]*a[i] + b[i]*b[i] + c[i]*c[i]);
+    mat.set(i, 1, a[i]);
+    mat.set(i, 2, b[i]);
+  }
+  m14 = mat.determinant();
   if (fabs(m11) < MYEPSILON)
 …
   if (fabs(a.distance(*center) - RADIUS) > MYEPSILON)
     DoeLog(1) && (eLog()<< Verbose(1) << "The given center is further way by " << fabs(a.distance(*center) - RADIUS) << " from a than RADIUS." << endl);
-  gsl_matrix_free(A);
 };
 …
   Vector x3;
   Vector x4;
-};
-/**
- * Intersection calculation function.
+ *
- * @param x to find the result for
- * @param function parameter
- */
-double MinIntersectDistance(const gsl_vector * x, void *params)
+{
-        Info FunctionInfo(__func__);
-  double retval = 0;
-  struct Intersection *I = (struct Intersection *)params;
-  Vector intersection;
-  for (int i=0;i<NDIM;i++)
-    intersection[i] = gsl_vector_get(x, i);
-  Vector SideA = I->x1 -I->x2 ;
-  Vector HeightA = intersection - I->x1;
-  HeightA.ProjectOntoPlane(SideA);
-  Vector SideB = I->x3 - I->x4;
-  Vector HeightB = intersection - I->x3;
-  HeightB.ProjectOntoPlane(SideB);
-  retval = HeightA.ScalarProduct(HeightA) + HeightB.ScalarProduct(HeightB);
-  //Log() << Verbose(1) << "MinIntersectDistance called, result: " << retval << endl;
-  return retval;
-};
-/**
- * Calculates whether there is an intersection between two lines. The first line
- * always goes through point 1 and point 2 and the second line is given by the
- * connection between point 4 and point 5.
+ *
- * @param point 1 of line 1
- * @param point 2 of line 1
- * @param point 1 of line 2
- * @param point 2 of line 2
+ *
- * @return true if there is an intersection between the given lines, false otherwise
- */
-bool existsIntersection(const Vector &point1, const Vector &point2, const Vector &point3, const Vector &point4)
+{
-        Info FunctionInfo(__func__);
-  bool result;
-  struct Intersection par;
-    par.x1 = point1;
-    par.x2 = point2;
-    par.x3 = point3;
-    par.x4 = point4;
-    const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
-    gsl_multimin_fminimizer *s = NULL;
-    gsl_vector *ss, *x;
-    gsl_multimin_function minexFunction;
-    size_t iter = 0;
-    int status;
-    double size;
-    /* Starting point */
-    x = gsl_vector_alloc(NDIM);
-    gsl_vector_set(x, 0, point1[0]);
-    gsl_vector_set(x, 1, point1[1]);
-    gsl_vector_set(x, 2, point1[2]);
-    /* Set initial step sizes to 1 */
-    ss = gsl_vector_alloc(NDIM);
-    gsl_vector_set_all(ss, 1.0);
-    /* Initialize method and iterate */
-    minexFunction.n = NDIM;
-    minexFunction.f = &MinIntersectDistance;
-    minexFunction.params = (void *)&par;
-    s = gsl_multimin_fminimizer_alloc(T, NDIM);
-    gsl_multimin_fminimizer_set(s, &minexFunction, x, ss);
-    do {
-        iter++;
-        status = gsl_multimin_fminimizer_iterate(s);
-        if (status) {
-          break;
+        }
-        size = gsl_multimin_fminimizer_size(s);
-        status = gsl_multimin_test_size(size, 1e-2);
-        if (status == GSL_SUCCESS) {
-          DoLog(1) && (Log() << Verbose(1) << "converged to minimum" <<  endl);
+        }
-    } while (status == GSL_CONTINUE && iter < 100);
-    // check whether intersection is in between or not
-  Vector intersection;
-  double t1, t2;
-  for (int i = 0; i < NDIM; i++) {
-    intersection[i] = gsl_vector_get(s->x, i);
+  }
-  Vector SideA = par.x2 - par.x1;
-  Vector HeightA = intersection - par.x1;
-  t1 = HeightA.ScalarProduct(SideA)/SideA.ScalarProduct(SideA);
-  Vector SideB = par.x4 - par.x3;
-  Vector HeightB = intersection - par.x3;
-  t2 = HeightB.ScalarProduct(SideB)/SideB.ScalarProduct(SideB);
-  Log() << Verbose(1) << "Intersection " << intersection << " is at "
-    << t1 << " for (" << point1 << "," << point2 << ") and at "
-    << t2 << " for (" << point3 << "," << point4 << "): ";
-  if (((t1 >= 0) && (t1 <= 1)) && ((t2 >= 0) && (t2 <= 1))) {
-    DoLog(1) && (Log() << Verbose(1) << "true intersection." << endl);
-    result = true;
-  } else {
-    DoLog(1) && (Log() << Verbose(1) << "intersection out of region of interest." << endl);
-    result = false;
+  }
-  // free minimizer stuff
-    gsl_vector_free(x);
-    gsl_vector_free(ss);
-    gsl_multimin_fminimizer_free(s);
-  return result;
 };

src/tesselationhelpers.hpp

-              rba9f5b
+              rbe97a8
 #endif
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_multimin.h>
+#include <gsl/gsl_permutation.h>
+#include <gsl/gsl_vector.h>
+#include <iostream>
+#include <iosfwd>
 #include "defs.hpp"
 …
 /********************************************** declarations *******************************/
-double DetGet(gsl_matrix * const A, const int inPlace);
 void GetSphere(Vector * const Center, const Vector &a, const Vector &b, const Vector &c, const double RADIUS);
 void GetCenterOfSphere(Vector* const Center, const Vector &a, const Vector &b, const Vector &c, Vector * const NewUmkreismittelpunkt, const Vector* const Direction, const Vector* const AlternativeDirection, const double HalfplaneIndicator, const double AlternativeIndicator, const double alpha, const double beta, const double gamma, const double RADIUS, const double Umkreisradius);
 void GetCenterofCircumcircle(Vector * const Center, const Vector &a, const Vector &b, const Vector &c);
 double GetPathLengthonCircumCircle(const Vector &CircleCenter, const Vector &CirclePlaneNormal, const double CircleRadius, const Vector &NewSphereCenter, const Vector &OldSphereCenter, const Vector &NormalVector, const Vector &SearchDirection);
-double MinIntersectDistance(const gsl_vector * x, void *params);
-bool existsIntersection(const Vector &point1, const Vector &point2, const Vector &point3, const Vector &point4);
 double CalculateVolumeofGeneralTetraeder(const Vector &a, const Vector &b, const Vector &c, const Vector &d);
 double CalculateAreaofGeneralTriangle(const Vector &A, const Vector &B, const Vector &C);

src/triangleintersectionlist.cpp

rba9f5b	rbe97a8
18	18	#include "tesselation.hpp"
19	19	#include "vector.hpp"
	20	#include "verbose.hpp"
20	21
21	22	/** Constructor for class TriangleIntersectionList.

src/unittests/ActOnAllUnitTest.cpp

rba9f5b	rbe97a8
14	14	#include "../test/ActOnAlltest.hpp"
15	15	#include "ActOnAllUnitTest.hpp"
16		~~#include "memoryallocator.hpp"~~
17	16	#include "vector.hpp"
18	17

src/unittests/LineUnittest.cpp

-              rba9f5b
+              rbe97a8
 #include <iostream>
+#include <cmath>
 using namespace std;
 …
   CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
+}
+void LineUnittest::sphereIntersectionTest(){
+  {
+    std::vector<Vector> res = la1->getSphereIntersections();
+    CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
+    CPPUNIT_ASSERT(testDirection(res[0],e1));
+    CPPUNIT_ASSERT(testDirection(res[1],e1));
+    CPPUNIT_ASSERT(res[0]!=res[1]);
+  }
+  {
+    std::vector<Vector> res = la2->getSphereIntersections();
+    CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
+    CPPUNIT_ASSERT(testDirection(res[0],e2));
+    CPPUNIT_ASSERT(testDirection(res[1],e2));
+    CPPUNIT_ASSERT(res[0]!=res[1]);
+  }
+  {
+    std::vector<Vector> res = la3->getSphereIntersections();
+    CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
+    CPPUNIT_ASSERT(testDirection(res[0],e3));
+    CPPUNIT_ASSERT(testDirection(res[1],e3));
+    CPPUNIT_ASSERT(res[0]!=res[1]);
+  }
+  {
+    std::vector<Vector> res = lp1->getSphereIntersections();
+    CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
+    CPPUNIT_ASSERT((res[0]==e1) || (res[0]==e2));
+    CPPUNIT_ASSERT((res[1]==e1) || (res[1]==e2));
+    CPPUNIT_ASSERT(res[0]!=res[1]);
+  }
+  {
+    std::vector<Vector> res = lp2->getSphereIntersections();
+    CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
+    CPPUNIT_ASSERT((res[0]==e2) || (res[0]==e3));
+    CPPUNIT_ASSERT((res[1]==e2) || (res[1]==e3));
+    CPPUNIT_ASSERT(res[0]!=res[1]);
+  }
+  {
+    std::vector<Vector> res = lp3->getSphereIntersections();
+    CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
+    CPPUNIT_ASSERT((res[0]==e3) || (res[0]==e1));
+    CPPUNIT_ASSERT((res[1]==e3) || (res[1]==e1));
+    CPPUNIT_ASSERT(res[0]!=res[1]);
+  }
+}

src/unittests/LineUnittest.hpp

-              rba9f5b
+              rbe97a8
   CPPUNIT_TEST ( intersectionTest );
   CPPUNIT_TEST ( rotationTest );
+  CPPUNIT_TEST ( sphereIntersectionTest );
   CPPUNIT_TEST_SUITE_END();
 …
   void intersectionTest();
   void rotationTest();
+  void sphereIntersectionTest();
 private:

src/unittests/Makefile.am

-              rba9f5b
+              rbe97a8
   LogUnitTest \
   manipulateAtomsTest \
+  MemoryUsageObserverUnitTest \
+  MemoryAllocatorUnitTest \
+  MatrixUnittest \
   MoleculeDescriptorTest \
   ObserverTest \
 …
   periodentafelTest \
   PlaneUnittest \
+  ShapeUnittest \
   SingletonTest \
   StackClassUnitTest \
 …
   listofbondsunittest.cpp \
   logunittest.cpp \
+  MatrixUnittest.cpp \
   manipulateAtomsTest.cpp \
-  memoryallocatorunittest.cpp  \
-  memoryusageobserverunittest.cpp \
   MoleculeDescriptorTest.cpp \
   ObserverTest.cpp \
 …
   periodentafelTest.cpp \
   PlaneUnittest.cpp \
+  ShapeUnittest.cpp \
   SingletonTest.cpp \
   stackclassunittest.cpp \
 …
   logunittest.hpp \
   manipulateAtomsTest.hpp \
+  memoryallocatorunittest.hpp  \
+  memoryusageobserverunittest.hpp \
+  MatrixUnittest.hpp \
   MoleculeDescriptorTest.hpp \
   periodentafelTest.hpp \
 …
 manipulateAtomsTest_LDADD = ${ALLLIBS}
+MemoryAllocatorUnitTest_SOURCES = UnitTestMain.cpp ../memoryallocator.hpp ../memoryallocator.cpp ../memoryusageobserver.cpp ../memoryusageobserver.hpp memoryallocatorunittest.cpp memoryallocatorunittest.hpp
+MemoryAllocatorUnitTest_LDADD = ${ALLLIBS}
+MemoryUsageObserverUnitTest_SOURCES = UnitTestMain.cpp ../memoryallocator.hpp ../memoryusageobserver.cpp ../memoryusageobserver.hpp memoryusageobserverunittest.cpp memoryusageobserverunittest.hpp
+MemoryUsageObserverUnitTest_LDADD = ${ALLLIBS}
+MatrixUnittest_SOURCES = UnitTestMain.cpp MatrixUnittest.cpp MatrixUnittest.hpp
+MatrixUnittest_LDADD = ${ALLLIBS}
 MoleculeDescriptorTest_SOURCES = UnitTestMain.cpp MoleculeDescriptorTest.cpp MoleculeDescriptorTest.hpp
 …
 PlaneUnittest_LDADD = ${ALLLIBS}
+ShapeUnittest_SOURCES = UnitTestMain.cpp ShapeUnittest.cpp ShapeUnittest.hpp
+ShapeUnittest_LDADD = ${ALLLIBS}
 SingletonTest_SOURCES = UnitTestMain.cpp SingletonTest.cpp SingletonTest.hpp
 SingletonTest_LDADD = ${ALLLIBS} $(BOOST_LIB) ${BOOST_THREAD_LIB}
 …
 Tesselation_InOutsideUnitTest_LDADD = ${ALLLIBS}
 TestRunner_SOURCES = TestRunnerMain.cpp ../memoryallocator.hpp ../memoryallocator.cpp ../memoryusageobserver.cpp ../memoryusageobserver.hpp $(TESTSOURCES) $(TESTHEADERS)
+TestRunner_SOURCES = TestRunnerMain.cpp $(TESTSOURCES) $(TESTHEADERS)
 TestRunner_LDADD = ${ALLLIBS}

src/unittests/ObserverTest.cpp

-              rba9f5b
+              rbe97a8
 };
 class ObservableCollection : public Observable {
+class ObservableSet : public Observable {
 public:
   typedef std::set<SimpleObservable*> set;
 …
   typedef set::const_iterator const_iterator;
   ObservableCollection(int _num) :
+  ObservableSet(int _num) :
     Observable("ObservableCollection"),
     num(_num)
 …
+  }
   ~ObservableCollection(){
+  ~ObservableSet(){
     set::iterator iter;
     for(iter=theSet.begin(); iter!=theSet.end(); ++iter ){
       delete (*iter);
+    }
+  }
+  iterator begin(){
+    return iterator(theSet.begin(),this);
+  }
+  iterator end(){
+    return iterator(theSet.end(),this);
+  }
+  const int num;
+private:
+  set theSet;
+};
+class ObservableMap : public Observable {
+public:
+  typedef std::map<int,SimpleObservable*> set;
+  typedef ObservedIterator<set> iterator;
+  typedef set::const_iterator const_iterator;
+  ObservableMap(int _num) :
+    Observable("ObservableCollection"),
+    num(_num)
+  {
+    for(int i=0; i<num; ++i){
+      SimpleObservable *content = new SimpleObservable();
+      content->signOn(this);
+      theSet.insert(make_pair(i,content));
+    }
+  }
+  ~ObservableMap(){
+    set::iterator iter;
+    for(iter=theSet.begin(); iter!=theSet.end(); ++iter ){
+      delete iter->second;
+    }
+  }
 …
   blockObservable = new BlockObservable();
   notificationObservable = new NotificationObservable();
+  collection = new ObservableCollection(5);
+  obsset = new ObservableSet(5);
+  obsmap = new ObservableMap(5);
   observer1 = new UpdateCountObserver();
 …
   delete blockObservable;
   delete notificationObservable;
+  delete collection;
+  delete obsset;
+  delete obsmap;
   delete observer1;
 …
   simpleObservable2->signOn(observer4);
+  CPPUNIT_ASSERT_EQUAL( 0, observer1->updates );
+  CPPUNIT_ASSERT_EQUAL( 0, observer2->updates );
+  CPPUNIT_ASSERT_EQUAL( 0, observer3->updates );
+  CPPUNIT_ASSERT_EQUAL( 0, observer4->updates );
   simpleObservable1->changeMethod();
   CPPUNIT_ASSERT_EQUAL( 1, observer1->updates );
 …
 void ObserverTest::doesBlockUpdateTest() {
   callObservable->signOn(observer1);
+  CPPUNIT_ASSERT_EQUAL( 0, observer1->updates );
   callObservable->changeMethod1();
 …
   CPPUNIT_ASSERT_EQUAL( 2, observer1->updates );
   CPPUNIT_ASSERT_EQUAL( 2, observer2->updates );
+}
+void ObserverTest::outsideLockTest(){
+  callObservable->signOn(observer1);
+  CPPUNIT_ASSERT_EQUAL( 0, observer1->updates );
+  {
+    LOCK_OBSERVABLE(*callObservable);
+    CPPUNIT_ASSERT_EQUAL( 0, observer1->updates );
+  }
+  // lock is gone now, observer should have notified
+  CPPUNIT_ASSERT_EQUAL( 1, observer1->updates );
+}
 …
   int i = 0;
   // test the general iterator methods
   for(ObservableCollection::iterator iter=collection->begin(); iter!=collection->end();++iter){
     CPPUNIT_ASSERT(i< collection->num);
+  for(ObservableSet::iterator iter=obsset->begin(); iter!=obsset->end();++iter){
+    CPPUNIT_ASSERT(i< obsset->num);
     i++;
+  }
   i=0;
   for(ObservableCollection::const_iterator iter=collection->begin(); iter!=collection->end();++iter){
     CPPUNIT_ASSERT(i<collection->num);
     i++;
+  }
   collection->signOn(observer1);
+  for(ObservableSet::const_iterator iter=obsset->begin(); iter!=obsset->end();++iter){
+    CPPUNIT_ASSERT(i<obsset->num);
+    i++;
+  }
+  obsset->signOn(observer1);
+  {
     // we construct this out of the loop, so the iterator dies at the end of
     // the scope and not the end of the loop (allows more testing)
     ObservableCollection::iterator iter;
     for(iter=collection->begin(); iter!=collection->end(); ++iter){
+    ObservableSet::iterator iter;
+    for(iter=obsset->begin(); iter!=obsset->end(); ++iter){
       (*iter)->changeMethod();
+    }
 …
   // when using a const_iterator no changes should be propagated
   for(ObservableCollection::const_iterator iter = collection->begin(); iter!=collection->end();++iter);
+  for(ObservableSet::const_iterator iter = obsset->begin(); iter!=obsset->end();++iter);
   CPPUNIT_ASSERT_EQUAL( 1, observer1->updates);
+  collection->signOff(observer1);
+  // we need to test the operator-> as well
+  obsmap->signOn(observer2);
+  {
+    // we construct this out of the loop, so the iterator dies at the end of
+    // the scope and not the end of the loop (allows more testing)
+    ObservableMap::iterator iter;
+    for(iter=obsmap->begin(); iter!=obsmap->end(); ++iter){
+      iter->second->changeMethod();
+    }
+    // At this point no change should have been propagated
+    CPPUNIT_ASSERT_EQUAL( 0, observer2->updates);
+  }
+  // After the Iterator has died the propagation should take place
+  CPPUNIT_ASSERT_EQUAL( 1, observer2->updates);
+  obsset->signOff(observer1);
+  obsmap->signOff(observer2);
+}

src/unittests/ObserverTest.hpp

-              rba9f5b
+              rbe97a8
 class CallObservable;
 class SuperObservable;
+class ObservableCollection;
+class ObservableSet;
+class ObservableMap;
 class BlockObservable;
 class NotificationObservable;
 …
   CPPUNIT_TEST ( doesBlockUpdateTest );
   CPPUNIT_TEST ( doesSubObservableTest );
+  CPPUNIT_TEST ( outsideLockTest );
   CPPUNIT_TEST ( doesNotifyTest );
   CPPUNIT_TEST ( doesReportTest );
 …
   void doesBlockUpdateTest();
   void doesSubObservableTest();
+  void outsideLockTest();
   void doesNotifyTest();
   void doesReportTest();
 …
   SuperObservable *superObservable;
   NotificationObservable *notificationObservable;
+  ObservableCollection *collection;
+  ObservableSet *obsset;
+  ObservableMap *obsmap;
 };

src/unittests/PlaneUnittest.cpp

rba9f5b	rbe97a8
11	11	#include <cppunit/extensions/TestFactoryRegistry.h>
12	12	#include <cppunit/ui/text/TestRunner.h>
	13
	14	#include <cmath>
13	15
14	16	#ifdef HAVE_TESTRUNNER

src/unittests/linearsystemofequationsunittest.cpp

rba9f5b	rbe97a8
12	12	#include <cppunit/extensions/TestFactoryRegistry.h>
13	13	#include <cppunit/ui/text/TestRunner.h>
	14	#include <cmath>
14	15
15	16	#include "linearsystemofequationsunittest.hpp"

src/unittests/tesselation_insideoutsideunittest.cpp

rba9f5b	rbe97a8
17	17	#include "tesselation.hpp"
18	18	#include "tesselation_insideoutsideunittest.hpp"
	19	#include "verbose.hpp"
19	20
20	21	#ifdef HAVE_TESTRUNNER

src/unittests/vectorunittest.cpp

-              rba9f5b
+              rbe97a8
 #include "defs.hpp"
 #include "log.hpp"
-#include "memoryusageobserver.hpp"
 #include "vector.hpp"
 #include "vector_ops.hpp"
 …
 #include "Plane.hpp"
 #include "Exceptions/LinearDependenceException.hpp"
+#include "Matrix.hpp"
 #ifdef HAVE_TESTRUNNER
 …
   CPPUNIT_ASSERT(testVector.ScalarProduct(three) < MYEPSILON);
+}
-/**
- * UnitTest for Vector::IsInParallelepiped().
- */
-void VectorTest::IsInParallelepipedTest()
+{
-  double parallelepiped[NDIM*NDIM];
-  parallelepiped[0] = 1;
-  parallelepiped[1] = 0;
-  parallelepiped[2] = 0;
-  parallelepiped[3] = 0;
-  parallelepiped[4] = 1;
-  parallelepiped[5] = 0;
-  parallelepiped[6] = 0;
-  parallelepiped[7] = 0;
-  parallelepiped[8] = 1;
-  fixture = zero;
-  CPPUNIT_ASSERT_EQUAL( false, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
-  fixture = Vector(2.5,2.5,2.5);
-  CPPUNIT_ASSERT_EQUAL( true, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
-  fixture = Vector(1.,1.,1.);
-  CPPUNIT_ASSERT_EQUAL( false, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
-  fixture = Vector(3.5,3.5,3.5);
-  CPPUNIT_ASSERT_EQUAL( false, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
-  fixture = Vector(2.,2.,2.);
-  CPPUNIT_ASSERT_EQUAL( true, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
-  fixture = Vector(2.,3.,2.);
-  CPPUNIT_ASSERT_EQUAL( true, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
-  fixture = Vector(-2.,2.,-1.);
-  CPPUNIT_ASSERT_EQUAL( false, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
+}

src/unittests/vectorunittest.hpp

-              rba9f5b
+              rbe97a8
     CPPUNIT_TEST ( ProjectionTest );
     CPPUNIT_TEST ( NormalsTest );
-    CPPUNIT_TEST ( IsInParallelepipedTest );
     CPPUNIT_TEST_SUITE_END();
 …
     void LineIntersectionTest();
     void VectorRotationTest();
-    void IsInParallelepipedTest();
 private:

src/vector.cpp

-              rba9f5b
+              rbe97a8
 #include "vector.hpp"
+#include "VectorContent.hpp"
 #include "verbose.hpp"
 #include "World.hpp"
 #include "Helpers/Assert.hpp"
 #include "Helpers/fast_functions.hpp"
+#include "Exceptions/MathException.hpp"
 #include <iostream>
+#include <gsl/gsl_blas.h>
 using namespace std;
 …
 Vector::Vector()
+{
   content = gsl_vector_calloc (NDIM);
+  content = new VectorContent();
 };
 …
 Vector::Vector(const Vector& src)
+{
+  content = gsl_vector_alloc(NDIM);
+  gsl_vector_set(content,0,src[0]);
+  gsl_vector_set(content,1,src[1]);
+  gsl_vector_set(content,2,src[2]);
+  content = new VectorContent();
+  gsl_vector_memcpy(content->content, src.content->content);
+}
 …
 Vector::Vector(const double x1, const double x2, const double x3)
+{
+  content = gsl_vector_alloc(NDIM);
+  gsl_vector_set(content,0,x1);
+  gsl_vector_set(content,1,x2);
+  gsl_vector_set(content,2,x3);
+};
+  content = new VectorContent();
+  gsl_vector_set(content->content,0,x1);
+  gsl_vector_set(content->content,1,x2);
+  gsl_vector_set(content->content,2,x3);
+};
+Vector::Vector(VectorContent *_content) :
+  content(_content)
+{}
 /**
 …
   // check for self assignment
   if(&src!=this){
+    gsl_vector_set(content,0,src[0]);
+    gsl_vector_set(content,1,src[1]);
+    gsl_vector_set(content,2,src[2]);
+    gsl_vector_memcpy(content->content, src.content->content);
+  }
   return *this;
 …
  */
 Vector::~Vector() {
   gsl_vector_free(content);
+  delete content;
 };
 …
+}
-/** Calculates distance between this and another vector in a periodic cell.
- * \param *y array to second vector
- * \param *cell_size 6-dimensional array with (xx, xy, yy, xz, yz, zz) entries specifying the periodic cell
- * \return \f$| x - y |\f$
- */
-double Vector::PeriodicDistance(const Vector &y, const double * const cell_size) const
+{
-  double res = distance(y), tmp, matrix[NDIM*NDIM];
-    Vector Shiftedy, TranslationVector;
-    int N[NDIM];
-    matrix[0] = cell_size[0];
-    matrix[1] = cell_size[1];
-    matrix[2] = cell_size[3];
-    matrix[3] = cell_size[1];
-    matrix[4] = cell_size[2];
-    matrix[5] = cell_size[4];
-    matrix[6] = cell_size[3];
-    matrix[7] = cell_size[4];
-    matrix[8] = cell_size[5];
-    // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
-    for (N[0]=-1;N[0]<=1;N[0]++)
-      for (N[1]=-1;N[1]<=1;N[1]++)
-        for (N[2]=-1;N[2]<=1;N[2]++) {
-          // create the translation vector
-          TranslationVector.Zero();
-          for (int i=NDIM;i--;)
-            TranslationVector[i] = (double)N[i];
-          TranslationVector.MatrixMultiplication(matrix);
-          // add onto the original vector to compare with
-          Shiftedy = y + TranslationVector;
-          // get distance and compare with minimum so far
-          tmp = distance(Shiftedy);
-          if (tmp < res) res = tmp;
+        }
-    return (res);
-};
-/** Calculates distance between this and another vector in a periodic cell.
- * \param *y array to second vector
- * \param *cell_size 6-dimensional array with (xx, xy, yy, xz, yz, zz) entries specifying the periodic cell
- * \return \f$| x - y |^2\f$
- */
-double Vector::PeriodicDistanceSquared(const Vector &y, const double * const cell_size) const
+{
-  double res = DistanceSquared(y), tmp, matrix[NDIM*NDIM];
-    Vector Shiftedy, TranslationVector;
-    int N[NDIM];
-    matrix[0] = cell_size[0];
-    matrix[1] = cell_size[1];
-    matrix[2] = cell_size[3];
-    matrix[3] = cell_size[1];
-    matrix[4] = cell_size[2];
-    matrix[5] = cell_size[4];
-    matrix[6] = cell_size[3];
-    matrix[7] = cell_size[4];
-    matrix[8] = cell_size[5];
-    // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
-    for (N[0]=-1;N[0]<=1;N[0]++)
-      for (N[1]=-1;N[1]<=1;N[1]++)
-        for (N[2]=-1;N[2]<=1;N[2]++) {
-          // create the translation vector
-          TranslationVector.Zero();
-          for (int i=NDIM;i--;)
-            TranslationVector[i] = (double)N[i];
-          TranslationVector.MatrixMultiplication(matrix);
-          // add onto the original vector to compare with
-          Shiftedy = y + TranslationVector;
-          // get distance and compare with minimum so far
-          tmp = DistanceSquared(Shiftedy);
-          if (tmp < res) res = tmp;
+        }
-    return (res);
-};
-/** Keeps the vector in a periodic cell, defined by the symmetric \a *matrix.
- * \param *out ofstream for debugging messages
- * Tries to translate a vector into each adjacent neighbouring cell.
- */
-void Vector::KeepPeriodic(const double * const matrix)
+{
-  //  int N[NDIM];
-  //  bool flag = false;
-    //vector Shifted, TranslationVector;
-  //  Log() << Verbose(1) << "Begin of KeepPeriodic." << endl;
-  //  Log() << Verbose(2) << "Vector is: ";
-  //  Output(out);
-  //  Log() << Verbose(0) << endl;
-    InverseMatrixMultiplication(matrix);
-    for(int i=NDIM;i--;) { // correct periodically
-      if (at(i) < 0) {  // get every coefficient into the interval [0,1)
-        at(i) += ceil(at(i));
-      } else {
-        at(i) -= floor(at(i));
+      }
+    }
-    MatrixMultiplication(matrix);
-  //  Log() << Verbose(2) << "New corrected vector is: ";
-  //  Output(out);
-  //  Log() << Verbose(0) << endl;
-  //  Log() << Verbose(1) << "End of KeepPeriodic." << endl;
-};
 /** Calculates scalar product between this and another vector.
  * \param *y array to second vector
 …
+{
   double res = 0.;
+  for (int i=NDIM;i--;)
+    res += at(i)*y[i];
+  gsl_blas_ddot(content->content, y.content->content, &res);
   return (res);
 };
 …
 double& Vector::operator[](size_t i){
   ASSERT(i<=NDIM && i>=0,"Vector Index out of Range");
   return *gsl_vector_ptr (content, i);
+  return *gsl_vector_ptr (content->content, i);
+}
 const double& Vector::operator[](size_t i) const{
   ASSERT(i<=NDIM && i>=0,"Vector Index out of Range");
   return *gsl_vector_ptr (content, i);
+  return *gsl_vector_ptr (content->content, i);
+}
 …
+}
 gsl_vector* Vector::get(){
+VectorContent* Vector::get(){
   return content;
+}
 …
 };
+void Vector::ScaleAll(const Vector &factor){
+  gsl_vector_mul(content->content, factor.content->content);
+}
 void Vector::Scale(const double factor)
+{
+  for (int i=NDIM;i--;)
+    at(i) *= factor;
+};
+/** Given a box by its matrix \a *M and its inverse *Minv the vector is made to point within that box.
+ * \param *M matrix of box
+ * \param *Minv inverse matrix
+ */
+void Vector::WrapPeriodically(const double * const M, const double * const Minv)
+{
+  MatrixMultiplication(Minv);
+  // truncate to [0,1] for each axis
+  for (int i=0;i<NDIM;i++) {
+    //at(i) += 0.5;  // set to center of box
+    while (at(i) >= 1.)
+      at(i) -= 1.;
+    while (at(i) < 0.)
+      at(i) += 1.;
+  }
+  MatrixMultiplication(M);
+  gsl_vector_scale(content->content,factor);
 };
 …
   return make_pair(res,helper);
+}
-/** Do a matrix multiplication.
- * \param *matrix NDIM_NDIM array
- */
-void Vector::MatrixMultiplication(const double * const M)
+{
-  Vector tmp;
-  // do the matrix multiplication
-  for(int i=NDIM;i--;)
-    tmp[i] = M[i]*at(0)+M[i+3]*at(1)+M[i+6]*at(2);
-  (*this) = tmp;
-};
-/** Do a matrix multiplication with the \a *A' inverse.
- * \param *matrix NDIM_NDIM array
- */
-bool Vector::InverseMatrixMultiplication(const double * const A)
+{
-  double B[NDIM*NDIM];
-  double detA = RDET3(A);
-  double detAReci;
-  // calculate the inverse B
-  if (fabs(detA) > MYEPSILON) {;  // RDET3(A) yields precisely zero if A irregular
-    detAReci = 1./detA;
-    B[0] =  detAReci*RDET2(A[4],A[5],A[7],A[8]);    // A_11
-    B[1] = -detAReci*RDET2(A[1],A[2],A[7],A[8]);    // A_12
-    B[2] =  detAReci*RDET2(A[1],A[2],A[4],A[5]);    // A_13
-    B[3] = -detAReci*RDET2(A[3],A[5],A[6],A[8]);    // A_21
-    B[4] =  detAReci*RDET2(A[0],A[2],A[6],A[8]);    // A_22
-    B[5] = -detAReci*RDET2(A[0],A[2],A[3],A[5]);    // A_23
-    B[6] =  detAReci*RDET2(A[3],A[4],A[6],A[7]);    // A_31
-    B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]);    // A_32
-    B[8] =  detAReci*RDET2(A[0],A[1],A[3],A[4]);    // A_33
-    MatrixMultiplication(B);
-    return true;
-  } else {
-    return false;
+  }
-};
 /** Creates this vector as the b y *factors' components scaled linear combination of the given three.
 …
 void Vector::AddVector(const Vector &y)
+{
+  for(int i=NDIM;i--;)
+    at(i) += y[i];
+  gsl_vector_add(content->content, y.content->content);
+}
 …
 void Vector::SubtractVector(const Vector &y)
+{
+  for(int i=NDIM;i--;)
+    at(i) -= y[i];
+}
+/**
+ * Checks whether this vector is within the parallelepiped defined by the given three vectors and
+ * their offset.
+ *
+ * @param offest for the origin of the parallelepiped
+ * @param three vectors forming the matrix that defines the shape of the parallelpiped
+ */
+bool Vector::IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const
+{
+  Vector a = (*this)-offset;
+  a.InverseMatrixMultiplication(parallelepiped);
+  bool isInside = true;
+  for (int i=NDIM;i--;)
+    isInside = isInside && ((a[i] <= 1) && (a[i] >= 0));
+  return isInside;
+  gsl_vector_sub(content->content, y.content->content);
+}

src/vector.hpp

-              rba9f5b
+              rbe97a8
 #endif
+#include <iostream>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_multimin.h>
+#include <iosfwd>
 #include <memory>
 …
 class Vector;
+class Matrix;
+struct VectorContent;
 typedef std::vector<Vector> pointset;
 …
  */
 class Vector : public Space{
+  friend Vector operator*(const Matrix&,const Vector&);
+  friend class Matrix;
 public:
   Vector();
   Vector(const double x1, const double x2, const double x3);
 …
   double DistanceSquared(const Vector &y) const;
   double DistanceToSpace(const Space& space) const;
-  double PeriodicDistance(const Vector &y, const double * const cell_size) const;
-  double PeriodicDistanceSquared(const Vector &y, const double * const cell_size) const;
   double ScalarProduct(const Vector &y) const;
   double Angle(const Vector &y) const;
 …
   Vector Projection(const Vector &y) const;
   void ScaleAll(const double *factor);
+  void ScaleAll(const Vector &factor);
   void Scale(const double factor);
-  void MatrixMultiplication(const double * const M);
-  bool InverseMatrixMultiplication(const double * const M);
-  void KeepPeriodic(const double * const matrix);
   bool GetOneNormalVector(const Vector &x1);
   bool MakeNormalTo(const Vector &y1);
-  bool IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const;
-  void WrapPeriodically(const double * const M, const double * const Minv);
   std::pair<Vector,Vector> partition(const Vector&) const;
   std::pair<pointset,Vector> partition(const pointset&) const;
 …
   // Access to internal structure
   gsl_vector* get();
+  VectorContent* get();
   // Methods that are derived directly from other methods
 …
 private:
+  gsl_vector *content;
+  Vector(VectorContent *);
+  VectorContent *content;
 };

src/vector_ops.cpp

rba9f5b	rbe97a8
23	23	#include <gsl/gsl_permutation.h>
24	24	#include <gsl/gsl_vector.h>
	25	#include <gsl/gsl_multimin.h>
25	26
26	27	/**

src/verbose.cpp

rba9f5b	rbe97a8
5	5	#include "info.hpp"
6	6	#include "verbose.hpp"
	7	#include <iostream>
7	8
8	9	/** Prints the tabs according to verbosity stored in the temporary constructed class.

src/verbose.hpp

rba9f5b	rbe97a8
18	18	#endif
19	19
20		#include <ios~~tream~~>
	20	#include <iosfwd>
21	21
22	22	/*********************************** Class Verbose & Binary *****************************/

test_all.sh

-              rba9f5b
+              rbe97a8
 docheck=0;
 docheck_mem=0;
+if [ -n "$TMPDIR" ]
+then
+  tmpdir="$TMPDIR";
+else
+  tmpdir="/tmp";
+fi
+tmppattern="MolecuilderTest";
 function usage(){
 …
   echo "  -c                    Only configure and compile (implies -s)";
   echo "  -O <opt-level>        Only compile this optimization level";
+}
+while getopts âho:f:scO:â OPTION
+  echo "  -t <tmpDir>           Use tmpDir as temporary directory";
+  echo "  -p <prefix>           Prefix to use for directory names (standart MolecuilderTest)";
+}
+while getopts âho:f:scO:t:p:â OPTION
 do
   case $OPTION in
 …
       ;;
    o)
      outfile=$OPTARG;
+     outfile="$OPTARG";
      ;;
    f)
      logfile=$OPTARG;
+     logfile="$OPTARG";
      ;;
    s)
 …
      optimizations=("-O$OPTARG");
      ;;
+   t)
+     tmpdir="$OPTARG";
+     ;;
+   p)
+     tmppattern="$OPTARG";
+     ;;
    ?)
      usage;
 …
 done
+if [[ -z $outfile ]] || [[ -z $logfile ]]
+# test if all arguments were provided
+if [[ -z "$outfile" ]] || [[ -z "$logfile" ]] || [[ -z "$tmpdir" ]] || [[ -z "$tmppattern" ]]
 then
   usage;
 …
 fi
+# turn all relative paths into absolutes
 outfile=`realpath -s $outfile`;
 logfile=`realpath -s $logfile`;
+tmpdir=`realpath -s $tmpdir`;
 …
 function run(){
   echo "Testing with \"$1\":" >> $outfile;
   testdir=`mktemp -d --tmpdir MolecuilderTest.XXXXXXXXXX`;
+  testdir=`mktemp -d --tmpdir=$tmpdir $tmppattern.XXXXXXXXXX`;
   basedir=$PWD;
   cd $testdir;

tests/Tesselations/defs.in

-              rba9f5b
+              rbe97a8
         fi
         #echo "Molecuilder done with exitcode $exitcode."
+        cd ../..
         #cat stderr
         #cat stdout
         grep -E "^[0-9]* [0-9]* [0-9]*$" ../../@srcdir@/$mol/$2/${FILENAME}-$mol.dat | sort -n >reference-triangles.dat
         grep -E "^[0-9]* [0-9]* [0-9]*$" ${FILENAME}.dat | sort -n >new-triangles.dat
         diff reference-triangles.dat new-triangles.dat 2>diffstderr >diffstdout || exitcode=$?
+        grep -E "^[0-9]* [0-9]* [0-9]*$" @srcdir@/$mol/$2/${FILENAME}-$mol.dat | sort -n >$testdir/$RADIUS/reference-triangles.dat
+        grep -E "^[0-9]* [0-9]* [0-9]*$" $testdir/$RADIUS/${FILENAME}.dat | sort -n >$testdir/$RADIUS/new-triangles.dat
+        diff $testdir/$RADIUS/reference-triangles.dat $testdir/$RADIUS/new-triangles.dat 2>$testdir/$RADIUS/diffstderr >$testdir/$RADIUS/diffstdout || exitcode=$?
         #echo "Diff done with exitcode $exitcode."
         #cat diffstderr
         #cat diffstdout
-        cd ../..
         test $exitcode = $expected_exitcode || exit 1
+}

tests/Tesselations/heptan/1.5/NonConvexEnvelope-heptan.dat

-              rba9f5b
+              rbe97a8
 TITLE = "3D CONVEX SHELL"
 VARIABLES = "X" "Y" "Z" "U"
 ZONE T="heptan", N=23, E=64, DATAPACKING=POINT, ZONETYPE=FETRIANGLE
+ZONE T="none", N=23, E=64, DATAPACKING=POINT, ZONETYPE=FETRIANGLE
 -7.27e-07 -1.22006 0.930455 18.7229
 -7.27e-07 -1.22006 -0.849545 18.7229
 …
 -1.2492 0.921941 -0.849545 18.7227
 .2492 0.921941 -0.849545 18.7222
 .2492 0.921941 0.930455 27.0641
 -2.4985 -1.22006 -0.849545 19.9769
 -2.4985 -1.22006 0.930455 27.4727
 .4985 -1.22006 0.930455 19.9769
 .4985 -1.22006 -0.849545 27.4727
+.2492 0.921941 0.930455 18.7222
+-2.4985 -1.22006 -0.849545 27.4727
+-2.4985 -1.22006 0.930455 19.9769
+.4985 -1.22006 0.930455 27.4727
+.4985 -1.22006 -0.849545 19.9769
 -4.6377 -0.336759 0.0404545 21.541
 -3.7477 0.921941 0.930455 18.8853
 …
 .6377 -0.336759 0.0404545 10.6618
 .7477 0.921941 -0.849545 18.5406
 .7477 0.921941 0.930455 12.1988
+.7477 0.921941 0.930455 18.5406
 -7.27e-07 -0.590759 0.0404545 23.0174
 -1.2492 0.292641 0.0404545 23.0167
 .2492 0.292641 0.0404545 20.1632
+.2492 0.292641 0.0404545 23.0172
 -2.4985 -0.590759 0.0404545 18.9798
 .4985 -0.590759 0.0404545 18.9798
 -3.7477 0.292641 0.0404545 39.5267
 .7477 0.292641 0.0404545 23.2512
+.7477 0.292641 0.0404545 20.5497
 15 23
 …
 19 23
 15 19
+19 23
+16 19
 6 19
 6 19
 …
 4 18
 4 18
+18 22
+12 22
 18 22
 13 22
-18 22
-12 22
 13 22
 13 22
+19 23
+16 23
+12 22
+12 22
+11 22
+12 22
 13 22
 13 22
 11 22
 13 22
-12 22
-12 22
-11 22
-12 22
 16 23
 16 23
 …
 10 21
 14 21
 17 21
 10 17
+17 21
+9 17
 2 17
 2 17
 17 21
 9 21
 17 20
 7 20
+17 21
+10 21
+17 20
+8 20
 8 20
 8 20
+11 20
 11 20
+11 20
+17 20
+8 17
+17 20
+7 17

tests/regression/Domain/4/post/test.conf

rba9f5b	rbe97a8
47	47	BoxLength # (Length of a unit cell)
48	48	1
49		0 0
	49	0 1
50	50	0 0 2
51	51

tests/regression/testsuite-analysis.at

rba9f5b	rbe97a8
39	39	AT_KEYWORDS([analysis])
40	40	AT_CHECK([/bin/cp -f ${abs_top_srcdir}/${AUTOTEST_PATH}/Analysis/4/pre/test.conf .], 0)
41		AT_CHECK([../../molecuilder -i test.conf -e ${abs_top_srcdir}/src/ -v 3 -I -C S --elements 1 --output-file output.csv --bin-output-file bin_output.csv --bin-start 0 --bin-width 1. --bin-end 20 --molecule-by-id 208], 0, [stdout], [stderr])
	41	AT_CHECK([../../molecuilder -i test.conf -e ${abs_top_srcdir}/src/ -v 3 -I -C S --elements 1 --output-file output.csv --bin-output-file bin_output.csv --bin-start 0 --bin-width 1. --bin-end 20 --molecule-by-id 207], 0, [stdout], [stderr])
42	42	AT_CHECK([fgrep "Begin of CorrelationToSurface" stdout], 0, [ignore], [ignore])
43	43	#AT_CHECK([file=output.csv; diff $file ${abs_top_srcdir}/${AUTOTEST_PATH}/Analysis/4/post/$file], 0, [ignore], [ignore])

Context Navigation

Changes in / [ba9f5b:be97a8]

Legend:

Download in other formats: