Changes in / [5ec8e3:56fb09]


Ignore:
Files:
17 deleted
51 edited

Legend:

Unmodified
Added
Removed

  • src/Actions/AnalysisAction/PairCorrelationAction.cpp

    r5ec8e3 r56fb09  
    6767  dialog = UIFactory::getInstance().makeDialog();
    6868  if (type == "P")
    69     dialog->queryVector("position", &Point, false, MapOfActions::getInstance().getDescription("position"));
     69    dialog->queryVector("position", &Point, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription("position"));
    7070  if (type == "S")
    7171    dialog->queryMolecule("molecule-by-id", &Boundary, MapOfActions::getInstance().getDescription("molecule-by-id"));

  • src/Actions/AtomAction/AddAction.cpp

    r5ec8e3 r56fb09  
    4141
    4242  dialog->queryElement(NAME, &elements, MapOfActions::getInstance().getDescription(NAME));
    43   dialog->queryVector("position", &position, true, MapOfActions::getInstance().getDescription("position"));
     43  dialog->queryVector("position", &position, World::getInstance().getDomain(), true, MapOfActions::getInstance().getDescription("position"));
    4444  cout << "pre-dialog" << endl;
    4545

  • src/Actions/MoleculeAction/FillWithMoleculeAction.cpp

    r5ec8e3 r56fb09  
    6262
    6363  dialog->queryString(NAME, &filename, MapOfActions::getInstance().getDescription(NAME));
    64   dialog->queryVector("distances", &distances, false, MapOfActions::getInstance().getDescription("distances"));
    65   dialog->queryVector("lengths", &lengths, false, MapOfActions::getInstance().getDescription("lengths"));
     64  dialog->queryVector("distances", &distances, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription("distances"));
     65  dialog->queryVector("lengths", &lengths, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription("lengths"));
    6666  dialog->queryBoolean("DoRotate", &DoRotate, MapOfActions::getInstance().getDescription("DoRotate"));
    6767  dialog->queryDouble("MaxDistance", &MaxDistance, MapOfActions::getInstance().getDescription("MaxDistance"));

  • src/Actions/MoleculeAction/TranslateAction.cpp

    r5ec8e3 r56fb09  
    5555  bool periodic = false;
    5656
    57   dialog->queryVector(NAME, &x, false, MapOfActions::getInstance().getDescription(NAME));
     57  dialog->queryVector(NAME, &x, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription(NAME));
    5858  dialog->queryMolecule("molecule-by-id", &mol, MapOfActions::getInstance().getDescription("molecule-by-id"));
    5959  dialog->queryBoolean("periodic", &periodic, MapOfActions::getInstance().getDescription("periodic"));

  • src/Actions/WorldAction/AddEmptyBoundaryAction.cpp

    r5ec8e3 r56fb09  
    4141  int j=0;
    4242
    43   dialog->queryVector(NAME, &boundary, false, MapOfActions::getInstance().getDescription(NAME));
     43  dialog->queryVector(NAME, &boundary, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription(NAME));
    4444
    4545  if(dialog->display()) {
    … …  
    5959    }
    6060    // set new box size
    61     double * const cell_size = new double[6];
     61    double * const cell_size = World::getInstance().getDomain();
    6262    for (j=0;j<6;j++)
    6363      cell_size[j] = 0.;
    … …  
    6767      cell_size[j] = (Max[i]-Min[i]+2.*boundary[i]);
    6868    }
    69     World::getInstance().setDomain(cell_size);
    70     delete[] cell_size;
    7169    // translate all atoms, such that Min is aty (0,0,0)
    7270    AtomRunner = AllAtoms.begin();

  • src/Actions/WorldAction/CenterInBoxAction.cpp

    r5ec8e3 r56fb09  
    3434  Dialog *dialog = UIFactory::getInstance().makeDialog();
    3535
    36   Box& cell_size = World::getInstance().getDomain();
     36  double * cell_size = World::getInstance().getDomain();
    3737  dialog->queryBox(NAME, &cell_size, MapOfActions::getInstance().getDescription(NAME));
    3838

  • src/Actions/WorldAction/CenterOnEdgeAction.cpp

    r5ec8e3 r56fb09  
    1313#include "vector.hpp"
    1414#include "World.hpp"
    15 #include "Matrix.hpp"
    1615
    1716#include <iostream>
    … …  
    3837  Vector Min;
    3938  Vector Max;
     39  int j=0;
    4040
    4141  dialog->queryEmpty(NAME, MapOfActions::getInstance().getDescription(NAME));
    … …  
    5757    }
    5858    // set new box size
    59     Matrix domain;
     59    double * const cell_size = World::getInstance().getDomain();
     60    for (j=0;j<6;j++)
     61      cell_size[j] = 0.;
     62    j=-1;
    6063    for (int i=0;i<NDIM;i++) {
    61       double tmp = Max[i]-Min[i];
    62       tmp = fabs(tmp)>=1. ? tmp : 1.0;
    63       domain.at(i,i) = tmp;
     64      j += i+1;
     65      cell_size[j] = (Max[i]-Min[i]);
    6466    }
    65     cout << "new domain is: " << domain << endl;
    66     World::getInstance().setDomain(domain);
    6767    // translate all atoms, such that Min is aty (0,0,0)
    6868    for (vector<atom*>::iterator AtomRunner = AllAtoms.begin(); AtomRunner != AllAtoms.end(); ++AtomRunner)

  • src/Actions/WorldAction/ChangeBoxAction.cpp

    r5ec8e3 r56fb09  
    1212#include "verbose.hpp"
    1313#include "World.hpp"
    14 #include "Box.hpp"
    15 #include "Matrix.hpp"
    1614
    1715#include <iostream>
    … …  
    3634  Dialog *dialog = UIFactory::getInstance().makeDialog();
    3735
    38   Box& cell_size = World::getInstance().getDomain();
     36  double * cell_size = World::getInstance().getDomain();
    3937  dialog->queryBox(NAME, &cell_size, MapOfActions::getInstance().getDescription(NAME));
    4038
    4139  if(dialog->display()) {
    42     DoLog(0) && (Log() << Verbose(0) << "Setting box domain to " << cell_size.getM() << endl);
     40    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);
    4341    delete dialog;
    4442    return Action::success;

  • src/Actions/WorldAction/RemoveSphereOfAtomsAction.cpp

    r5ec8e3 r56fb09  
    4141
    4242  dialog->queryDouble(NAME, &radius, MapOfActions::getInstance().getDescription(NAME));
    43   dialog->queryVector("position", &point, false, MapOfActions::getInstance().getDescription("position"));
     43  dialog->queryVector("position", &point, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription("position"));
    4444
    4545  if(dialog->display()) {

  • src/Actions/WorldAction/RepeatBoxAction.cpp

    r5ec8e3 r56fb09  
    1313#include "molecule.hpp"
    1414#include "vector.hpp"
    15 #include "Matrix.hpp"
    1615#include "verbose.hpp"
    1716#include "World.hpp"
    18 #include "Box.hpp"
    1917
    2018#include <iostream>
    … …  
    4846  MoleculeListClass *molecules = World::getInstance().getMolecules();
    4947
    50   dialog->queryVector(NAME, &Repeater, false, MapOfActions::getInstance().getDescription(NAME));
     48  dialog->queryVector(NAME, &Repeater, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription(NAME));
    5149  //dialog->queryMolecule("molecule-by-id", &mol,MapOfActions::getInstance().getDescription("molecule-by-id"));
    5250  vector<molecule *> AllMolecules;
    … …  
    6159  if(dialog->display()) {
    6260    (cout << "Repeating box " << Repeater << " times for (x,y,z) axis." << endl);
    63     Matrix M = World::getInstance().getDomain().getM();
    64     Matrix newM = M;
     61    double * const cell_size = World::getInstance().getDomain();
     62    double *M = ReturnFullMatrixforSymmetric(cell_size);
    6563    Vector x,y;
    6664    int n[NDIM];
    67     Matrix repMat;
    6865    for (int axis = 0; axis < NDIM; axis++) {
    6966      Repeater[axis] = floor(Repeater[axis]);
    … …  
    7269        Repeater[axis] = 1;
    7370      }
    74       repMat.at(axis,axis) = Repeater[axis];
     71      cell_size[(abs(axis+1) == 2) ? 2 : ((abs(axis+2) == 3) ? 5 : 0)] *= Repeater[axis];
    7572    }
    76     newM *= repMat;
    77     World::getInstance().setDomain(newM);
    7873
    7974    molecule *newmol = NULL;
    … …  
    10398                DoeLog(1) && (eLog()<< Verbose(1) << "AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl);
    10499              x = y;
    105               x *= M;
     100              x.MatrixMultiplication(M);
    106101              newmol = World::getInstance().createMolecule();
    107102              molecules->insert(newmol);
    … …  
    123118      }
    124119    }
     120    delete(M);
    125121    delete dialog;
    126122    return Action::success;

  • src/Actions/WorldAction/ScaleBoxAction.cpp

    r5ec8e3 r56fb09  
    1414#include "verbose.hpp"
    1515#include "World.hpp"
    16 #include "Box.hpp"
    17 #include "Matrix.hpp"
    1816
    1917#include <iostream>
    … …  
    3937  Vector Scaler;
    4038  double x[NDIM];
     39  int j=0;
    4140
    42   dialog->queryVector(NAME, &Scaler, false, MapOfActions::getInstance().getDescription(NAME));
     41  dialog->queryVector(NAME, &Scaler, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription(NAME));
    4342
    4443  if(dialog->display()) {
    … …  
    5049      (*AtomRunner)->x.ScaleAll(x);
    5150    }
    52 
    53     Matrix M = World::getInstance().getDomain().getM();
    54     Matrix scale;
    55 
     51    j = -1;
     52    double * const cell_size = World::getInstance().getDomain();
    5653    for (int i=0;i<NDIM;i++) {
    57       scale.at(i,i) = x[i];
     54      j += i+1;
     55      cell_size[j]*=x[i];
    5856    }
    59     M *= scale;
    60     World::getInstance().setDomain(M);
    6157
    6258    delete dialog;

  • src/Line.cpp

    r5ec8e3 r56fb09  
    215215}
    216216
    217 std::vector<Vector> Line::getSphereIntersections() const{
    218   std::vector<Vector> res;
    219 
    220   // line is kept in normalized form, so we can skip a lot of calculations
    221   double discriminant = 1-origin->NormSquared();
    222   // we might have 2, 1 or 0 solutions, depending on discriminant
    223   if(discriminant>=0){
    224     if(discriminant==0){
    225       res.push_back(*origin);
    226     }
    227     else{
    228       Vector helper = sqrt(discriminant)*(*direction);
    229       res.push_back(*origin+helper);
    230       res.push_back(*origin-helper);
    231     }
    232   }
    233   return res;
    234 }
    235 
    236217Line makeLineThrough(const Vector &x1, const Vector &x2){
    237218  if(x1==x2){

  • src/Line.hpp

    r5ec8e3 r56fb09  
    3838  Plane getOrthogonalPlane(const Vector &origin) const;
    3939
    40   std::vector<Vector> getSphereIntersections() const;
    41 
    4240private:
    4341  std::auto_ptr<Vector> origin;

  • src/Makefile.am

    r5ec8e3 r56fb09  
    3838  vector.cpp
    3939                           
    40 LINALGHEADER = \
    41   gslmatrix.hpp \
     40LINALGHEADER = gslmatrix.hpp \
    4241  gslvector.hpp \
    4342  linearsystemofequations.hpp \
    … …  
    7675  Actions/MethodAction.hpp \
    7776  Actions/Process.hpp
    78 
    79 EXCEPTIONSOURCE = Exceptions/CustomException.cpp \
    80                                   Exceptions/LinearDependenceException.cpp \
    81                                   Exceptions/MathException.cpp \
    82                                   Exceptions/NotInvertibleException.cpp \
    83                                   Exceptions/SkewException.cpp \
    84                                   Exceptions/ZeroVectorException.cpp
    85                                  
    86 EXCEPTIONHEADER = Exceptions/CustomException.hpp \
    87                                   Exceptions/LinearDependenceException.hpp \
    88                                   Exceptions/MathException.hpp \
    89                                   Exceptions/NotInvertibleException.hpp \
    90                                   Exceptions/SkewException.hpp \
    91                                   Exceptions/ZeroVectorException.hpp
     77 
    9278
    9379PARSERSOURCE = \
    … …  
    115101  Patterns/Observer.hpp \
    116102  Patterns/Singleton.hpp
    117  
    118 SHAPESOURCE = \
    119   Shapes/BaseShapes.cpp \
    120   Shapes/Shape.cpp \
    121   Shapes/ShapeOps.cpp
    122 SHAPEHEADER = \
    123   Shapes/BaseShapes.hpp \
    124   Shapes/Shape.hpp \
    125   Shapes/ShapeOps.hpp
    126  
    127103
    128104QTUIMOC_HEADER = UIElements/QT4/QTDialog.hpp \
    … …  
    158134  Descriptors/MoleculeNameDescriptor.hpp \
    159135  Descriptors/MoleculePtrDescriptor.hpp
     136                                   
     137EXCEPTIONSOURCE = Exceptions/CustomException.cpp \
     138                                  Exceptions/LinearDependenceException.cpp \
     139                                  Exceptions/MathException.cpp \
     140                                  Exceptions/SkewException.cpp \
     141                                  Exceptions/ZeroVectorException.cpp
     142                                 
     143EXCEPTIONHEADER = Exceptions/CustomException.hpp \
     144                                  Exceptions/LinearDependenceException.hpp \
     145                                  Exceptions/MathException.hpp \
     146                                  Exceptions/SkewException.hpp \
     147                                  Exceptions/ZeroVectorException.hpp
    160148                                 
    161149QTUISOURCE = ${QTUIMOC_TARGETS} \
    … …  
    177165  ${ACTIONSSOURCE} \
    178166  ${ATOMSOURCE} \
    179   ${EXCEPTIONSOURCE} \
    180167  ${PATTERNSOURCE} \
    181168  ${PARSERSOURCE} \
    182   ${SHAPESOURCE} \
    183169  ${DESCRIPTORSOURCE} \
    184170  ${HELPERSOURCE} \
     171  ${EXCEPTIONSOURCE} \
    185172  bond.cpp \
    186173  bondgraph.cpp \
    187174  boundary.cpp \
    188   Box.cpp \
    189175  CommandLineParser.cpp \
    190176  config.cpp \
    … …  
    202188  log.cpp \
    203189  logger.cpp \
    204   Matrix.cpp \
    205190  moleculelist.cpp \
    206191  molecule.cpp \
    … …  
    227212  ${ACTIONSHEADER} \
    228213  ${ATOMHEADER} \
    229   ${EXCEPTIONHEADER} \
    230214  ${PARSERHEADER} \
    231215  ${PATTERNHEADER} \
    232   ${SHAPEHEADER} \
    233216  ${DESCRIPTORHEADER} \
     217  ${EXCEPTIONHEADER} \
    234218  bond.hpp \
    235219  bondgraph.hpp \
    236220  boundary.hpp \
    237   Box.hpp \
    238221  CommandLineParser.hpp \
    239222  config.hpp \
    … …  
    253236  log.hpp \
    254237  logger.hpp \
    255   Matrix.hpp \
    256238  molecule.hpp \
    257239  molecule_template.hpp \

  • src/Parser/PcpParser.cpp

    r5ec8e3 r56fb09  
    2020#include "verbose.hpp"
    2121#include "World.hpp"
    22 #include "Matrix.hpp"
    23 #include "Box.hpp"
    2422
    2523/** Constructor of PcpParser.
    … …  
    212210  // Unit cell and magnetic field
    213211  ParseForParameter(verbose,FileBuffer, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
    214   double *cell_size = new double[6];
     212  double * const cell_size = World::getInstance().getDomain();
    215213  cell_size[0] = BoxLength[0];
    216214  cell_size[1] = BoxLength[3];
    … …  
    219217  cell_size[4] = BoxLength[7];
    220218  cell_size[5] = BoxLength[8];
    221   World::getInstance().setDomain(cell_size);
    222   delete[] cell_size;
    223219  //if (1) fprintf(stderr,"\n");
    224220
    … …  
    331327void PcpParser::save(std::ostream* file)
    332328{
    333   const Matrix &domain = World::getInstance().getDomain().getM();
     329  const double * const cell_size = World::getInstance().getDomain();
    334330  class ThermoStatContainer *Thermostats = World::getInstance().getThermostats();
    335331  if (!file->fail()) {
    … …  
    416412    *file << endl;
    417413    *file << "BoxLength\t\t\t# (Length of a unit cell)" << endl;
    418     *file << domain.at(0,0) << "\t" << endl;
    419     *file << domain.at(1,0) << "\t" << domain.at(1,1) << "\t" << endl;
    420     *file << domain.at(2,0) << "\t" << domain.at(2,1) << "\t" << domain.at(2,2) << "\t" << endl;
     414    *file << cell_size[0] << "\t" << endl;
     415    *file << cell_size[1] << "\t" << cell_size[2] << "\t" << endl;
     416    *file << cell_size[3] << "\t" << cell_size[4] << "\t" << cell_size[5] << "\t" << endl;
    421417    // FIXME
    422418    *file << endl;

  • src/UIElements/CommandLineUI/CommandLineDialog.cpp

    r5ec8e3 r56fb09  
    2727#include "verbose.hpp"
    2828#include "World.hpp"
    29 #include "Box.hpp"
    3029
    3130#include "atom.hpp"
    … …  
    7473}
    7574
    76 void CommandLineDialog::queryVector(const char* title, Vector *target, bool check, string _description) {
    77   registerQuery(new VectorCommandLineQuery(title,target,check, _description));
    78 }
    79 
    80 void CommandLineDialog::queryBox(const char* title, Box* cellSize, string _description) {
     75void CommandLineDialog::queryVector(const char* title, Vector *target,const double *const cellSize, bool check, string _description) {
     76  registerQuery(new VectorCommandLineQuery(title,target,cellSize,check, _description));
     77}
     78
     79void CommandLineDialog::queryBox(const char* title, double ** const cellSize, string _description) {
    8180  registerQuery(new BoxCommandLineQuery(title,cellSize,_description));
    8281}
    … …  
    203202}
    204203
    205 CommandLineDialog::VectorCommandLineQuery::VectorCommandLineQuery(string title, Vector *_target, bool _check, string _description) :
    206     Dialog::VectorQuery(title,_target,_check, _description)
     204CommandLineDialog::VectorCommandLineQuery::VectorCommandLineQuery(string title, Vector *_target, const double *const _cellSize, bool _check, string _description) :
     205    Dialog::VectorQuery(title,_target,_cellSize,_check, _description)
    207206{}
    208207
    … …  
    225224
    226225
    227 CommandLineDialog::BoxCommandLineQuery::BoxCommandLineQuery(string title, Box* _cellSize, string _description) :
     226CommandLineDialog::BoxCommandLineQuery::BoxCommandLineQuery(string title, double ** const _cellSize, string _description) :
    228227    Dialog::BoxQuery(title,_cellSize, _description)
    229228{}

  • src/UIElements/CommandLineUI/CommandLineDialog.hpp

    r5ec8e3 r56fb09  
    3434  virtual void queryAtom(const char*,atom**, std::string = "");
    3535  virtual void queryMolecule(const char*,molecule**,std::string = "");
    36   virtual void queryVector(const char*,Vector *,bool, std::string = "");
    37   virtual void queryBox(const char*,Box *, std::string = "");
     36  virtual void queryVector(const char*,Vector *,const double * const,bool, std::string = "");
     37  virtual void queryBox(const char*,double ** const, std::string = "");
    3838  virtual void queryElement(const char*, std::vector<element *> *, std::string = "");
    3939
    … …  
    9191  class VectorCommandLineQuery : public Dialog::VectorQuery {
    9292  public:
    93     VectorCommandLineQuery(std::string title,Vector *_target,bool _check, std::string _description = "");
     93    VectorCommandLineQuery(std::string title,Vector *_target,const double *const _cellSize,bool _check, std::string _description = "");
    9494    virtual ~VectorCommandLineQuery();
    9595    virtual bool handle();
    … …  
    9898  class BoxCommandLineQuery : public Dialog::BoxQuery {
    9999  public:
    100     BoxCommandLineQuery(std::string title,Box* _cellSize, std::string _description = "");
     100    BoxCommandLineQuery(std::string title,double ** const _cellSize, std::string _description = "");
    101101    virtual ~BoxCommandLineQuery();
    102102    virtual bool handle();

  • src/UIElements/Dialog.cpp

    r5ec8e3 r56fb09  
    1414#include "molecule.hpp"
    1515#include "vector.hpp"
    16 #include "Matrix.hpp"
    17 #include "Box.hpp"
    1816
    1917using namespace std;
    … …  
    175173// Vector Queries
    176174
    177 Dialog::VectorQuery::VectorQuery(std::string title,Vector *_target,bool _check, std::string _description) :
     175Dialog::VectorQuery::VectorQuery(std::string title,Vector *_target,const double *const _cellSize,bool _check, std::string _description) :
    178176  Query(title, _description),
     177  cellSize(_cellSize),
    179178  check(_check),
    180179  target(_target)
    … …  
    194193// Box Queries
    195194
    196 Dialog::BoxQuery::BoxQuery(std::string title, Box* _cellSize, std::string _description) :
     195Dialog::BoxQuery::BoxQuery(std::string title, double ** const _cellSize, std::string _description) :
    197196  Query(title, _description),
    198197  target(_cellSize)
    … …  
    207206
    208207void Dialog::BoxQuery::setResult() {
    209   (*target)= ReturnFullMatrixforSymmetric(tmp);
     208  for (int i=0;i<6;i++) {
     209    (*target)[i] = tmp[i];
     210  }
    210211}
    211212

  • src/UIElements/Dialog.hpp

    r5ec8e3 r56fb09  
    1414
    1515class atom;
    16 class Box;
    1716class element;
    1817class molecule;
    … …  
    4342  virtual void queryAtom(const char*,atom**,std::string = "")=0;
    4443  virtual void queryMolecule(const char*,molecule**, std::string = "")=0;
    45   virtual void queryVector(const char*,Vector *,bool, std::string = "")=0;
    46   virtual void queryBox(const char*,Box*, std::string = "")=0;
     44  virtual void queryVector(const char*,Vector *,const double *const,bool, std::string = "")=0;
     45  virtual void queryBox(const char*,double ** const, std::string = "")=0;
    4746  virtual void queryElement(const char*, std::vector<element *> *, std::string = "")=0;
    4847
    … …  
    163162  class VectorQuery : public Query {
    164163  public:
    165       VectorQuery(std::string title,Vector *_target,bool _check, std::string _description = "");
     164      VectorQuery(std::string title,Vector *_target,const double *const _cellSize,bool _check, std::string _description = "");
    166165      virtual ~VectorQuery();
    167166      virtual bool handle()=0;
    … …  
    169168    protected:
    170169      Vector *tmp;
     170      const double *const cellSize;
    171171      bool check;
    172172    private:
    … …  
    176176  class BoxQuery : public Query {
    177177  public:
    178       BoxQuery(std::string title,Box *_cellSize, std::string _description = "");
     178      BoxQuery(std::string title,double ** const _cellSize, std::string _description = "");
    179179      virtual ~BoxQuery();
    180180      virtual bool handle()=0;
    181181      virtual void setResult();
    182182    protected:
    183       double* tmp;
     183      double *tmp;
    184184    private:
    185       Box* target;
     185      double **target;
    186186  };
    187187

  • src/UIElements/QT4/QTDialog.cpp

    r5ec8e3 r56fb09  
    2929#include "molecule.hpp"
    3030#include "Descriptors/MoleculeIdDescriptor.hpp"
    31 #include "Matrix.hpp"
    32 #include "Box.hpp"
    3331
    3432
    … …  
    9593}
    9694
    97 void QTDialog::queryBox(char const*, Box*, string){
     95void QTDialog::queryBox(char const*, double**, string){
    9896  // TODO
    9997  ASSERT(false, "Not implemented yet");
    … …  
    120118}
    121119
    122 void QTDialog::queryVector(const char* title, Vector *target, bool check,string) {
    123   registerQuery(new VectorQTQuery(title,target,check,inputLayout,this));
     120void QTDialog::queryVector(const char* title, Vector *target,const double *const cellSize, bool check,string) {
     121  registerQuery(new VectorQTQuery(title,target,cellSize,check,inputLayout,this));
    124122}
    125123
    … …  
    250248}
    251249
    252 QTDialog::VectorQTQuery::VectorQTQuery(std::string title, Vector *_target, bool _check,QBoxLayout *_parent,QTDialog *_dialog) :
    253     Dialog::VectorQuery(title,_target,_check),
    254     parent(_parent)
    255 {
    256   const Matrix& M = World::getInstance().getDomain().getM();
     250QTDialog::VectorQTQuery::VectorQTQuery(std::string title, Vector *_target, const double *const _cellSize, bool _check,QBoxLayout *_parent,QTDialog *_dialog) :
     251    Dialog::VectorQuery(title,_target,_cellSize,_check),
     252    parent(_parent)
     253{
     254  // 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
     255  int j = -1;
    257256  const char *coords[3] = {"x:","y:","z:"};
    258257  mainLayout= new QHBoxLayout();
    … …  
    262261  mainLayout->addLayout(subLayout);
    263262  for(int i=0; i<3; i++) {
     263    j+=i+1;
    264264    coordLayout[i] = new QHBoxLayout();
    265265    subLayout->addLayout(coordLayout[i]);
    … …  
    267267    coordLayout[i]->addWidget(coordLabel[i]);
    268268    coordInput[i] = new QDoubleSpinBox();
    269     coordInput[i]->setRange(0,M.at(i,i));
     269    coordInput[i]->setRange(0,cellSize[j]);
    270270    coordInput[i]->setDecimals(3);
    271271    coordLayout[i]->addWidget(coordInput[i]);

  • src/UIElements/QT4/QTDialog.hpp

    r5ec8e3 r56fb09  
    4242  virtual void queryAtom(const char*,atom**,std::string = "");
    4343  virtual void queryMolecule(const char*,molecule**,std::string = "");
    44   virtual void queryVector(const char*,Vector *,bool,std::string = "");
    45   virtual void queryBox(const char*,Box*, std::string = "");
     44  virtual void queryVector(const char*,Vector *,const double *const,bool,std::string = "");
     45  virtual void queryBox(const char*,double ** const, std::string = "");
    4646  virtual void queryElement(const char*,std::vector<element *> *_target,std::string = "");
    4747
    … …  
    109109    class VectorQTQuery : public Dialog::VectorQuery {
    110110    public:
    111       VectorQTQuery(std::string title,Vector *_target,bool _check,QBoxLayout *,QTDialog *);
     111      VectorQTQuery(std::string title,Vector *_target,const double *const _cellSize,bool _check,QBoxLayout *,QTDialog *);
    112112      virtual ~VectorQTQuery();
    113113      virtual bool handle();

  • src/UIElements/TextUI/TextDialog.cpp

    r5ec8e3 r56fb09  
    2525#include "molecule.hpp"
    2626#include "vector.hpp"
    27 #include "Matrix.hpp"
    28 #include "Box.hpp"
    2927
    3028using namespace std;
    … …  
    6866}
    6967
    70 void TextDialog::queryVector(const char* title, Vector *target, bool check, string description) {
    71   registerQuery(new VectorTextQuery(title,target,check,description));
    72 }
    73 
    74 void TextDialog::queryBox(const char* title,Box* cellSize, string description) {
     68void TextDialog::queryVector(const char* title, Vector *target,const double *const cellSize, bool check, string description) {
     69  registerQuery(new VectorTextQuery(title,target,cellSize,check,description));
     70}
     71
     72void TextDialog::queryBox(const char* title,double ** const cellSize, string description) {
    7573  registerQuery(new BoxTextQuery(title,cellSize,description));
    7674}
    … …  
    245243}
    246244
    247 TextDialog::VectorTextQuery::VectorTextQuery(std::string title, Vector *_target, bool _check, std::string _description) :
    248     Dialog::VectorQuery(title,_target,_check,_description)
     245TextDialog::VectorTextQuery::VectorTextQuery(std::string title, Vector *_target, const double *const _cellSize, bool _check, std::string _description) :
     246    Dialog::VectorQuery(title,_target,_cellSize,_check,_description)
    249247{}
    250248
    … …  
    255253  Log() << Verbose(0) << getTitle();
    256254
    257   const Matrix &M = World::getInstance().getDomain().getM();
    258255  char coords[3] = {'x','y','z'};
     256  int j = -1;
    259257  for (int i=0;i<3;i++) {
     258    j += i+1;
    260259    do {
    261       Log() << Verbose(0) << coords[i] << "[0.." << M.at(i,i) << "]: ";
     260      Log() << Verbose(0) << coords[i] << "[0.." << cellSize[j] << "]: ";
    262261      cin >> (*tmp)[i];
    263     } while ((((*tmp)[i] < 0) || ((*tmp)[i] >= M.at(i,i))) && (check));
     262    } while ((((*tmp)[i] < 0) || ((*tmp)[i] >= cellSize[j])) && (check));
    264263  }
    265264  return true;
    266265}
    267266
    268 TextDialog::BoxTextQuery::BoxTextQuery(std::string title, Box* _cellSize, std::string _description) :
     267TextDialog::BoxTextQuery::BoxTextQuery(std::string title, double ** const _cellSize, std::string _description) :
    269268    Dialog::BoxQuery(title,_cellSize,_description)
    270269{}

  • src/UIElements/TextUI/TextDialog.hpp

    r5ec8e3 r56fb09  
    3131  virtual void queryAtom(const char*,atom**,std::string = "");
    3232  virtual void queryMolecule(const char*,molecule**,std::string = "");
    33   virtual void queryVector(const char*,Vector *,bool, std::string = "");
    34   virtual void queryBox(const char*,Box*, std::string = "");
     33  virtual void queryVector(const char*,Vector *,const double * const,bool, std::string = "");
     34  virtual void queryBox(const char*,double ** const, std::string = "");
    3535  virtual void queryElement(const char*, std::vector<element *> *, std::string = "");
    3636
    … …  
    8888  class VectorTextQuery : public Dialog::VectorQuery {
    8989  public:
    90     VectorTextQuery(std::string title,Vector *_target,bool _check, std::string _description = NULL);
     90    VectorTextQuery(std::string title,Vector *_target,const double *const _cellSize,bool _check, std::string _description = NULL);
    9191    virtual ~VectorTextQuery();
    9292    virtual bool handle();
    … …  
    9595  class BoxTextQuery : public Dialog::BoxQuery {
    9696  public:
    97     BoxTextQuery(std::string title,Box* _cellSize, std::string _description = NULL);
     97    BoxTextQuery(std::string title,double ** const _cellSize, std::string _description = NULL);
    9898    virtual ~BoxTextQuery();
    9999    virtual bool handle();

  • src/VectorSet.hpp

    r5ec8e3 r56fb09  
    1212#include <functional>
    1313#include <algorithm>
    14 #include <limits>
    1514
    1615/**
    … …  
    2019 */
    2120
    22 #include "vector.hpp"
    23 #include <list>
     21class Vector;
    2422
    2523// this tests, whether we actually have a Vector
    … …  
    5048  void translate(const Vector &translater){
    5149    // this is needed to allow template lookup
    52     transform(this->begin(),this->end(),this->begin(),std::bind1st(std::plus<Vector>(),translater));
    53   }
    54 
    55   double minDistSquared(const Vector &point){
    56     if(!this->size())
    57       return std::numeric_limits<double>::infinity();
    58     std::list<double> helper;
    59     helper.resize(this->size());
    60     transform(this->begin(),this->end(),
    61               helper.begin(),
    62               std::bind2nd(std::mem_fun_ref(&Vector::DistanceSquared),point));
    63     return *min_element(helper.begin(),helper.end());
     50    transform(this->begin(),this->end(),this->begin(),bind1st(plus<Vector>(),translater));
    6451  }
    6552};
    6653
    67 // allows simpler definition of VectorSets
    68 #define VECTORSET(container_type) VectorSet<container_type<Vector> >
    69 
    7054#endif /* VECTORSET_HPP_ */

  • src/World.cpp

    r5ec8e3 r56fb09  
    2222#include "Actions/ManipulateAtomsProcess.hpp"
    2323#include "Helpers/Assert.hpp"
    24 #include "Box.hpp"
    25 #include "Matrix.hpp"
    2624
    2725#include "Patterns/Singleton_impl.hpp"
    … …  
    7674// system
    7775
    78 Box& World::getDomain() {
    79   return *cell_size;
    80 }
    81 
    82 void World::setDomain(const Matrix &mat){
    83   *cell_size = mat;
     76double * World::getDomain() {
     77  return cell_size;
    8478}
    8579
    8680void World::setDomain(double * matrix)
    8781{
    88   Matrix M = ReturnFullMatrixforSymmetric(matrix);
    89   cell_size->setM(M);
     82
    9083}
    9184
    … …  
    140133  molecules.erase(id);
    141134}
     135
     136double *World::cell_size = NULL;
    142137
    143138atom *World::createAtom(){
    … …  
    302297    molecules_deprecated(new MoleculeListClass(this))
    303298{
    304   cell_size = new Box;
    305   Matrix domain;
    306   domain.at(0,0) = 20;
    307   domain.at(1,1) = 20;
    308   domain.at(2,2) = 20;
    309   cell_size->setM(domain);
     299  cell_size = new double[6];
     300  cell_size[0] = 20.;
     301  cell_size[1] = 0.;
     302  cell_size[2] = 20.;
     303  cell_size[3] = 0.;
     304  cell_size[4] = 0.;
     305  cell_size[5] = 20.;
    310306  defaultName = "none";
    311307  molecules_deprecated->signOn(this);
    … …  
    315311{
    316312  molecules_deprecated->signOff(this);
    317   delete cell_size;
     313  delete[] cell_size;
    318314  delete molecules_deprecated;
    319315  delete periode;

  • src/World.hpp

    r5ec8e3 r56fb09  
    3434class AtomDescriptor_impl;
    3535template<typename T> class AtomsCalculation;
    36 class Box;
    3736class config;
    3837class ManipulateAtomsProcess;
    39 class Matrix;
    4038class molecule;
    4139class MoleculeDescriptor;
    … …  
    127125   * get the domain size as a symmetric matrix (6 components)
    128126   */
    129   Box& getDomain();
    130 
    131   /**
    132    * Set the domain size from a matrix object
    133    *
    134    * Matrix needs to be symmetric
    135    */
    136   void setDomain(const Matrix &mat);
     127  double * getDomain();
    137128
    138129  /**
    … …  
    266257  periodentafel *periode;
    267258  config *configuration;
    268   Box *cell_size;
     259  static double *cell_size;
    269260  std::string defaultName;
    270261  class ThermoStatContainer *Thermostats;

  • src/analysis_correlation.cpp

    r5ec8e3 r56fb09  
    1919#include "triangleintersectionlist.hpp"
    2020#include "vector.hpp"
    21 #include "Matrix.hpp"
    2221#include "verbose.hpp"
    2322#include "World.hpp"
    24 #include "Box.hpp"
    2523
    2624
    … …  
    3634  PairCorrelationMap *outmap = NULL;
    3735  double distance = 0.;
    38   Box &domain = World::getInstance().getDomain();
     36  double *domain = World::getInstance().getDomain();
    3937
    4038  if (molecules->ListOfMolecules.empty()) {
    … …  
    7775                for (set <pair<element *, element *> >::iterator PairRunner = PairsOfElements.begin(); PairRunner != PairsOfElements.end(); ++PairRunner)
    7876                  if ((PairRunner->first == (**iter).type) && (PairRunner->second == (**runner).type)) {
    79                     distance = domain.periodicDistance(*(*iter)->node,*(*runner)->node);
     77                    distance = (*iter)->node->PeriodicDistance(*(*runner)->node,  domain);
    8078                    //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
    8179                    outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> ((*iter), (*runner)) ) );
    … …  
    137135  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++){
    138136    if ((*MolWalker)->ActiveFlag) {
    139       Matrix FullMatrix = World::getInstance().getDomain().getM();
    140       Matrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
     137      double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
     138      double * FullInverseMatrix = InverseMatrix(FullMatrix);
    141139      DoeLog(2) && (eLog()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    142140      eLog() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
    143141      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    144142        DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
    145         periodicX = FullInverseMatrix * (*(*iter)->node); // x now in [0,1)^3
     143        periodicX = *(*iter)->node;
     144        periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
    146145        // go through every range in xyz and get distance
    147146        for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
    148147          for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
    149148            for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
    150               checkX = FullMatrix * (Vector(n[0], n[1], n[2]) + periodicX);
     149              checkX = Vector(n[0], n[1], n[2]) + periodicX;
     150              checkX.MatrixMultiplication(FullMatrix);
    151151              for (MoleculeList::const_iterator MolOtherWalker = MolWalker; MolOtherWalker != molecules->ListOfMolecules.end(); MolOtherWalker++){
    152152                if ((*MolOtherWalker)->ActiveFlag) {
    … …  
    157157                      for (set <pair<element *, element *> >::iterator PairRunner = PairsOfElements.begin(); PairRunner != PairsOfElements.end(); ++PairRunner)
    158158                        if ((PairRunner->first == (**iter).type) && (PairRunner->second == (**runner).type)) {
    159                           periodicOtherX = FullInverseMatrix * (*(*runner)->node); // x now in [0,1)^3
     159                          periodicOtherX = *(*runner)->node;
     160                          periodicOtherX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
    160161                          // go through every range in xyz and get distance
    161162                          for (Othern[0]=-ranges[0]; Othern[0] <= ranges[0]; Othern[0]++)
    162163                            for (Othern[1]=-ranges[1]; Othern[1] <= ranges[1]; Othern[1]++)
    163164                              for (Othern[2]=-ranges[2]; Othern[2] <= ranges[2]; Othern[2]++) {
    164                                 checkOtherX = FullMatrix * (Vector(Othern[0], Othern[1], Othern[2]) + periodicOtherX);
     165                                checkOtherX = Vector(Othern[0], Othern[1], Othern[2]) + periodicOtherX;
     166                                checkOtherX.MatrixMultiplication(FullMatrix);
    165167                                distance = checkX.distance(checkOtherX);
    166168                                //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
    … …  
    174176            }
    175177      }
     178      delete[](FullMatrix);
     179      delete[](FullInverseMatrix);
    176180    }
    177181  }
    … …  
    191195  CorrelationToPointMap *outmap = NULL;
    192196  double distance = 0.;
    193   Box &domain = World::getInstance().getDomain();
     197  double *cell_size = World::getInstance().getDomain();
    194198
    195199  if (molecules->ListOfMolecules.empty()) {
    … …  
    207211        for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
    208212          if ((*type == NULL) || ((*iter)->type == *type)) {
    209             distance = domain.periodicDistance(*(*iter)->node,*point);
     213            distance = (*iter)->node->PeriodicDistance(*point, cell_size);
    210214            DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
    211215            outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> ((*iter), point) ) );
    … …  
    242246  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    243247    if ((*MolWalker)->ActiveFlag) {
    244       Matrix FullMatrix = World::getInstance().getDomain().getM();
    245       Matrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
     248      double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
     249      double * FullInverseMatrix = InverseMatrix(FullMatrix);
    246250      DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    247251      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    … …  
    249253        for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
    250254          if ((*type == NULL) || ((*iter)->type == *type)) {
    251             periodicX = FullInverseMatrix * (*(*iter)->node); // x now in [0,1)^3
     255            periodicX = *(*iter)->node;
     256            periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
    252257            // go through every range in xyz and get distance
    253258            for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
    254259              for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
    255260                for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
    256                   checkX = FullMatrix * (Vector(n[0], n[1], n[2]) + periodicX);
     261                  checkX = Vector(n[0], n[1], n[2]) + periodicX;
     262                  checkX.MatrixMultiplication(FullMatrix);
    257263                  distance = checkX.distance(*point);
    258264                  DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
    … …  
    261267          }
    262268      }
     269      delete[](FullMatrix);
     270      delete[](FullInverseMatrix);
    263271    }
    264272
    … …  
    344352  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    345353    if ((*MolWalker)->ActiveFlag) {
    346       Matrix FullMatrix = World::getInstance().getDomain().getM();
    347       Matrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
     354      double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
     355      double * FullInverseMatrix = InverseMatrix(FullMatrix);
    348356      DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    349357      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    … …  
    351359        for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
    352360          if ((*type == NULL) || ((*iter)->type == *type)) {
    353             periodicX = FullInverseMatrix * (*(*iter)->node); // x now in [0,1)^3
     361            periodicX = *(*iter)->node;
     362            periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
    354363            // go through every range in xyz and get distance
    355364            ShortestDistance = -1.;
    … …  
    357366              for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
    358367                for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
    359                   checkX = FullMatrix * (Vector(n[0], n[1], n[2]) + periodicX);
     368                  checkX = Vector(n[0], n[1], n[2]) + periodicX;
     369                  checkX.MatrixMultiplication(FullMatrix);
    360370                  TriangleIntersectionList Intersections(&checkX,Surface,LC);
    361371                  distance = Intersections.GetSmallestDistance();
    … …  
    371381          }
    372382      }
     383      delete[](FullMatrix);
     384      delete[](FullInverseMatrix);
    373385    }
    374386

  • src/atom.cpp

    r5ec8e3 r56fb09  
    1717#include "World.hpp"
    1818#include "molecule.hpp"
    19 #include "Shapes/Shape.hpp"
    2019
    2120/************************************* Functions for class atom *************************************/
    … …  
    113112 * \return true - is inside, false - is not
    114113 */
    115 bool atom::IsInShape(const Shape& shape) const
    116 {
    117   return shape.isInside(*node);
     114bool atom::IsInParallelepiped(const Vector offset, const double *parallelepiped) const
     115{
     116  return (node->IsInParallelepiped(offset, parallelepiped));
    118117};
    119118

  • src/atom.hpp

    r5ec8e3 r56fb09  
    3535class World;
    3636class molecule;
    37 class Shape;
    3837
    3938/********************************************** declarations *******************************/
    … …  
    6766  double DistanceToVector(const Vector &origin) const;
    6867  double DistanceSquaredToVector(const Vector &origin) const;
    69   bool IsInShape(const Shape&) const;
     68  bool IsInParallelepiped(const Vector offset, const double *parallelepiped) const;
    7069
    7170  // getter and setter

  • src/boundary.cpp

    r5ec8e3 r56fb09  
    2222#include "World.hpp"
    2323#include "Plane.hpp"
    24 #include "Matrix.hpp"
    25 #include "Box.hpp"
    2624
    2725#include<gsl/gsl_poly.h>
    … …  
    766764  int N[NDIM];
    767765  int n[NDIM];
    768   const Matrix &M = World::getInstance().getDomain().getM();
    769   Matrix Rotations;
    770   const Matrix &MInverse = World::getInstance().getDomain().getMinv();
     766  double *M =  ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
     767  double Rotations[NDIM*NDIM];
     768  double *MInverse = InverseMatrix(M);
    771769  Vector AtomTranslations;
    772770  Vector FillerTranslations;
    … …  
    801799
    802800  // calculate filler grid in [0,1]^3
    803   FillerDistance = MInverse * Vector(distance[0], distance[1], distance[2]);
     801  FillerDistance = Vector(distance[0], distance[1], distance[2]);
     802  FillerDistance.InverseMatrixMultiplication(M);
    804803  for(int i=0;i<NDIM;i++)
    805804    N[i] = (int) ceil(1./FillerDistance[i]);
    … …  
    814813      for (n[2] = 0; n[2] < N[2]; n[2]++) {
    815814        // calculate position of current grid vector in untransformed box
    816         CurrentPosition = M * Vector((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
     815        CurrentPosition = Vector((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
     816        CurrentPosition.MatrixMultiplication(M);
    817817        // create molecule random translation vector ...
    818818        for (int i=0;i<NDIM;i++)
    … …  
    835835            }
    836836
    837             Rotations.set(0,0,  cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2])));
    838             Rotations.set(0,1,  sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2])));
    839             Rotations.set(0,2,              cos(phi[1])*sin(phi[2])                                        );
    840             Rotations.set(1,0, -sin(phi[0])*cos(phi[1])                                                    );
    841             Rotations.set(1,1,  cos(phi[0])*cos(phi[1])                                                    );
    842             Rotations.set(1,2,              sin(phi[1])                                                    );
    843             Rotations.set(2,0, -cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2])));
    844             Rotations.set(2,1, -sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2])));
    845             Rotations.set(2,2,              cos(phi[1])*cos(phi[2])                                        );
     837            Rotations[0] =   cos(phi[0])            *cos(phi[2]) + (sin(phi[0])*sin(phi[1])*sin(phi[2]));
     838            Rotations[3] =   sin(phi[0])            *cos(phi[2]) - (cos(phi[0])*sin(phi[1])*sin(phi[2]));
     839            Rotations[6] =               cos(phi[1])*sin(phi[2])                                     ;
     840            Rotations[1] = - sin(phi[0])*cos(phi[1])                                                ;
     841            Rotations[4] =   cos(phi[0])*cos(phi[1])                                                ;
     842            Rotations[7] =               sin(phi[1])                                                ;
     843            Rotations[3] = - cos(phi[0])            *sin(phi[2]) + (sin(phi[0])*sin(phi[1])*cos(phi[2]));
     844            Rotations[5] = - sin(phi[0])            *sin(phi[2]) - (cos(phi[0])*sin(phi[1])*cos(phi[2]));
     845            Rotations[8] =               cos(phi[1])*cos(phi[2])                                     ;
    846846          }
    847847
    … …  
    849849          Inserter = (*iter)->x;
    850850          if (DoRandomRotation)
    851             Inserter *= Rotations;
     851            Inserter.MatrixMultiplication(Rotations);
    852852          Inserter += AtomTranslations + FillerTranslations + CurrentPosition;
    853853
    854854          // check whether inserter is inside box
    855           Inserter *= MInverse;
     855          Inserter.MatrixMultiplication(MInverse);
    856856          FillIt = true;
    857857          for (int i=0;i<NDIM;i++)
    858858            FillIt = FillIt && (Inserter[i] >= -MYEPSILON) && ((Inserter[i]-1.) <= MYEPSILON);
    859           Inserter *= M;
     859          Inserter.MatrixMultiplication(M);
    860860
    861861          // Check whether point is in- or outside
    … …  
    892892            }
    893893      }
     894  delete[](M);
     895  delete[](MInverse);
    894896
    895897  return Filling;

  • src/config.cpp

    r5ec8e3 r56fb09  
    2626#include "ThermoStatContainer.hpp"
    2727#include "World.hpp"
    28 #include "Matrix.hpp"
    29 #include "Box.hpp"
    3028
    3129/************************************* Functions for class config ***************************/
    … …  
    681679  // Unit cell and magnetic field
    682680  ParseForParameter(verbose,FileBuffer, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
    683   double * cell_size = new double[6];
     681  double * const cell_size = World::getInstance().getDomain();
    684682  cell_size[0] = BoxLength[0];
    685683  cell_size[1] = BoxLength[3];
    … …  
    688686  cell_size[4] = BoxLength[7];
    689687  cell_size[5] = BoxLength[8];
    690   World::getInstance().setDomain(cell_size);
    691   delete cell_size;
    692688  //if (1) fprintf(stderr,"\n");
    693689
    … …  
    887883
    888884  ParseForParameter(verbose,file, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
    889   double * cell_size = new double[6];
     885  double * const cell_size = World::getInstance().getDomain();
    890886  cell_size[0] = BoxLength[0];
    891887  cell_size[1] = BoxLength[3];
    … …  
    894890  cell_size[4] = BoxLength[7];
    895891  cell_size[5] = BoxLength[8];
    896   World::getInstance().setDomain(cell_size);
    897   delete[] cell_size;
    898892  if (1) fprintf(stderr,"\n");
    899893  config::DoPerturbation = 0;
    … …  
    10331027  // bring MaxTypes up to date
    10341028  mol->CountElements();
    1035   const Matrix &domain = World::getInstance().getDomain().getM();
     1029  const double * const cell_size = World::getInstance().getDomain();
    10361030  ofstream * const output = new ofstream(filename, ios::out);
    10371031  if (output != NULL) {
    … …  
    11041098    *output << endl;
    11051099    *output << "BoxLength\t\t\t# (Length of a unit cell)" << endl;
    1106     *output << domain.at(0,0) << "\t" << endl;
    1107     *output << domain.at(1,0) << "\t" << domain.at(1,1) << "\t" << endl;
    1108     *output << domain.at(2,0) << "\t" << domain.at(2,1) << "\t" << domain.at(2,2) << "\t" << endl;
     1100    *output << cell_size[0] << "\t" << endl;
     1101    *output << cell_size[1] << "\t" << cell_size[2] << "\t" << endl;
     1102    *output << cell_size[3] << "\t" << cell_size[4] << "\t" << cell_size[5] << "\t" << endl;
    11091103    // FIXME
    11101104    *output << endl;

  • src/ellipsoid.cpp

    r5ec8e3 r56fb09  
    2121#include "tesselation.hpp"
    2222#include "vector.hpp"
    23 #include "Matrix.hpp"
    2423#include "verbose.hpp"
    2524
    … …  
    3534  Vector helper, RefPoint;
    3635  double distance = -1.;
    37   Matrix Matrix;
     36  double Matrix[NDIM*NDIM];
    3837  double InverseLength[3];
    3938  double psi,theta,phi; // euler angles in ZX'Z'' convention
    … …  
    5251  theta = EllipsoidAngle[1];
    5352  phi = EllipsoidAngle[2];
    54   Matrix.set(0,0, cos(psi)*cos(phi) - sin(psi)*cos(theta)*sin(phi));
    55   Matrix.set(1,0, -cos(psi)*sin(phi) - sin(psi)*cos(theta)*cos(phi));
    56   Matrix.set(2,0, sin(psi)*sin(theta));
    57   Matrix.set(0,1, sin(psi)*cos(phi) + cos(psi)*cos(theta)*sin(phi));
    58   Matrix.set(1,1, cos(psi)*cos(theta)*cos(phi) - sin(psi)*sin(phi));
    59   Matrix.set(2,1, -cos(psi)*sin(theta));
    60   Matrix.set(0,2, sin(theta)*sin(phi));
    61   Matrix.set(1,2, sin(theta)*cos(phi));
    62   Matrix.set(2,2, cos(theta));
    63   helper *= Matrix;
     53  Matrix[0] = cos(psi)*cos(phi) - sin(psi)*cos(theta)*sin(phi);
     54  Matrix[1] = -cos(psi)*sin(phi) - sin(psi)*cos(theta)*cos(phi);
     55  Matrix[2] = sin(psi)*sin(theta);
     56  Matrix[3] = sin(psi)*cos(phi) + cos(psi)*cos(theta)*sin(phi);
     57  Matrix[4] = cos(psi)*cos(theta)*cos(phi) - sin(psi)*sin(phi);
     58  Matrix[5] = -cos(psi)*sin(theta);
     59  Matrix[6] = sin(theta)*sin(phi);
     60  Matrix[7] = sin(theta)*cos(phi);
     61  Matrix[8] = cos(theta);
     62  helper.MatrixMultiplication(Matrix);
    6463  helper.ScaleAll(InverseLength);
    6564  //Log() << Verbose(4) << "Transformed RefPoint is at " << helper << "." << endl;
    … …  
    7473  phi = -EllipsoidAngle[2];
    7574  helper.ScaleAll(EllipsoidLength);
    76   Matrix.set(0,0, cos(psi)*cos(phi) - sin(psi)*cos(theta)*sin(phi));
    77   Matrix.set(1,0, -cos(psi)*sin(phi) - sin(psi)*cos(theta)*cos(phi));
    78   Matrix.set(2,0, sin(psi)*sin(theta));
    79   Matrix.set(0,1, sin(psi)*cos(phi) + cos(psi)*cos(theta)*sin(phi));
    80   Matrix.set(1,1, cos(psi)*cos(theta)*cos(phi) - sin(psi)*sin(phi));
    81   Matrix.set(2,1, -cos(psi)*sin(theta));
    82   Matrix.set(0,2, sin(theta)*sin(phi));
    83   Matrix.set(1,2, sin(theta)*cos(phi));
    84   Matrix.set(2,2, cos(theta));
    85   helper *= Matrix;
     75  Matrix[0] = cos(psi)*cos(phi) - sin(psi)*cos(theta)*sin(phi);
     76  Matrix[1] = -cos(psi)*sin(phi) - sin(psi)*cos(theta)*cos(phi);
     77  Matrix[2] = sin(psi)*sin(theta);
     78  Matrix[3] = sin(psi)*cos(phi) + cos(psi)*cos(theta)*sin(phi);
     79  Matrix[4] = cos(psi)*cos(theta)*cos(phi) - sin(psi)*sin(phi);
     80  Matrix[5] = -cos(psi)*sin(theta);
     81  Matrix[6] = sin(theta)*sin(phi);
     82  Matrix[7] = sin(theta)*cos(phi);
     83  Matrix[8] = cos(theta);
     84  helper.MatrixMultiplication(Matrix);
    8685  //Log() << Verbose(4) << "Intersection is at " << helper << "." << endl;
    8786

  • src/helpers.cpp

    r5ec8e3 r56fb09  
    117117};
    118118
     119/** Blows the 6-dimensional \a cell_size array up to a full NDIM by NDIM matrix.
     120 * \param *symm 6-dim array of unique symmetric matrix components
     121 * \return allocated NDIM*NDIM array with the symmetric matrix
     122 */
     123double * ReturnFullMatrixforSymmetric(const double * const symm)
     124{
     125  double *matrix = new double[NDIM * NDIM];
     126  matrix[0] = symm[0];
     127  matrix[1] = symm[1];
     128  matrix[2] = symm[3];
     129  matrix[3] = symm[1];
     130  matrix[4] = symm[2];
     131  matrix[5] = symm[4];
     132  matrix[6] = symm[3];
     133  matrix[7] = symm[4];
     134  matrix[8] = symm[5];
     135  return matrix;
     136};
     137
     138/** Calculate the inverse of a 3x3 matrix.
     139 * \param *matrix NDIM_NDIM array
     140 */
     141double * InverseMatrix( const double * const A)
     142{
     143  double *B = new double[NDIM * NDIM];
     144  double detA = RDET3(A);
     145  double detAReci;
     146
     147  for (int i=0;i<NDIM*NDIM;++i)
     148    B[i] = 0.;
     149  // calculate the inverse B
     150  if (fabs(detA) > MYEPSILON) {;  // RDET3(A) yields precisely zero if A irregular
     151    detAReci = 1./detA;
     152    B[0] =  detAReci*RDET2(A[4],A[5],A[7],A[8]);    // A_11
     153    B[1] = -detAReci*RDET2(A[1],A[2],A[7],A[8]);    // A_12
     154    B[2] =  detAReci*RDET2(A[1],A[2],A[4],A[5]);    // A_13
     155    B[3] = -detAReci*RDET2(A[3],A[5],A[6],A[8]);    // A_21
     156    B[4] =  detAReci*RDET2(A[0],A[2],A[6],A[8]);    // A_22
     157    B[5] = -detAReci*RDET2(A[0],A[2],A[3],A[5]);    // A_23
     158    B[6] =  detAReci*RDET2(A[3],A[4],A[6],A[7]);    // A_31
     159    B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]);    // A_32
     160    B[8] =  detAReci*RDET2(A[0],A[1],A[3],A[4]);    // A_33
     161  }
     162  return B;
     163};
     164
     165
     166
    119167/** Comparison function for GSL heapsort on distances in two molecules.
    120168 * \param *a

  • src/helpers.hpp

    r5ec8e3 r56fb09  
    5252bool IsValidNumber( const char *string);
    5353int CompareDoubles (const void * a, const void * b);
     54double * ReturnFullMatrixforSymmetric(const double * const cell_size);
     55double * InverseMatrix(const double * const A);
    5456void performCriticalExit();
    5557

  • src/molecule.cpp

    r5ec8e3 r56fb09  
    99#include <cstring>
    1010#include <boost/bind.hpp>
    11 #include <boost/foreach.hpp>
    1211
    1312#include "World.hpp"
    … …  
    2827#include "tesselation.hpp"
    2928#include "vector.hpp"
    30 #include "Matrix.hpp"
    3129#include "World.hpp"
    32 #include "Box.hpp"
    3330#include "Plane.hpp"
    3431#include "Exceptions/LinearDependenceException.hpp"
    … …  
    287284  Vector Orthovector1, Orthovector2;  // temporary vectors in coordination construction
    288285  Vector InBondvector;    // vector in direction of *Bond
    289   const Matrix &matrix =  World::getInstance().getDomain().getM();
     286  double *matrix = NULL;
    290287  bond *Binder = NULL;
     288  double * const cell_size = World::getInstance().getDomain();
    291289
    292290//  Log() << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl;
    … …  
    309307      } // (signs are correct, was tested!)
    310308    }
    311     Orthovector1 *= matrix;
     309    matrix = ReturnFullMatrixforSymmetric(cell_size);
     310    Orthovector1.MatrixMultiplication(matrix);
    312311    InBondvector -= Orthovector1; // subtract just the additional translation
     312    delete[](matrix);
    313313    bondlength = InBondvector.Norm();
    314314//    Log() << Verbose(4) << "Corrected InBondvector is now: ";
    … …  
    541541      break;
    542542  }
     543  delete[](matrix);
    543544
    544545//  Log() << Verbose(3) << "End of AddHydrogenReplacementAtom." << endl;
    … …  
    659660 * @param three vectors forming the matrix that defines the shape of the parallelpiped
    660661 */
    661 molecule* molecule::CopyMoleculeFromSubRegion(const Shape &region) const {
     662molecule* molecule::CopyMoleculeFromSubRegion(const Vector offset, const double *parallelepiped) const {
    662663  molecule *copy = World::getInstance().createMolecule();
    663664
    664   BOOST_FOREACH(atom *iter,atoms){
    665     if(iter->IsInShape(region)){
    666       copy->AddCopyAtom(iter);
    667     }
    668   }
     665  ActOnCopyWithEachAtomIfTrue ( &molecule::AddCopyAtom, copy, &atom::IsInParallelepiped, offset, parallelepiped );
    669666
    670667  //TODO: copy->BuildInducedSubgraph(this);
    … …  
    753750void molecule::SetBoxDimension(Vector *dim)
    754751{
    755   Matrix domain;
    756   for(int i =0; i<NDIM;++i)
    757     domain.at(i,i) = dim->at(i);
    758   World::getInstance().setDomain(domain);
     752  double * const cell_size = World::getInstance().getDomain();
     753  cell_size[0] = dim->at(0);
     754  cell_size[1] = 0.;
     755  cell_size[2] = dim->at(1);
     756  cell_size[3] = 0.;
     757  cell_size[4] = 0.;
     758  cell_size[5] = dim->at(2);
    759759};
    760760
    … …  
    849849bool molecule::CheckBounds(const Vector *x) const
    850850{
    851   const Matrix &domain = World::getInstance().getDomain().getM();
     851  double * const cell_size = World::getInstance().getDomain();
    852852  bool result = true;
     853  int j =-1;
    853854  for (int i=0;i<NDIM;i++) {
    854     result = result && ((x->at(i) >= 0) && (x->at(i) < domain.at(i,i)));
     855    j += i+1;
     856    result = result && ((x->at(i) >= 0) && (x->at(i) < cell_size[j]));
    855857  }
    856858  //return result;

  • src/molecule.hpp

    r5ec8e3 r56fb09  
    5353class periodentafel;
    5454class Vector;
    55 class Shape;
    5655
    5756/******************************** Some definitions for easier reading **********************************/
    … …  
    317316
    318317  molecule *CopyMolecule();
    319   molecule* CopyMoleculeFromSubRegion(const Shape&) const;
     318  molecule* CopyMoleculeFromSubRegion(const Vector offset, const double *parallelepiped) const;
    320319
    321320  /// Fragment molecule by two different approaches:

  • src/molecule_fragmentation.cpp

    r5ec8e3 r56fb09  
    2222#include "periodentafel.hpp"
    2323#include "World.hpp"
    24 #include "Matrix.hpp"
    25 #include "Box.hpp"
    2624
    2725/************************************* Functions for class molecule *********************************/
    … …  
    17101708  atom *Walker = NULL;
    17111709  atom *OtherWalker = NULL;
    1712   Matrix matrix = World::getInstance().getDomain().getM();
     1710  double * const cell_size = World::getInstance().getDomain();
     1711  double *matrix = ReturnFullMatrixforSymmetric(cell_size);
    17131712  enum Shading *ColorList = NULL;
    17141713  double tmp;
    … …  
    17501749          Translationvector[i] = (tmp < 0) ? +1. : -1.;
    17511750      }
    1752       Translationvector *= matrix;
     1751      Translationvector.MatrixMultiplication(matrix);
    17531752      //Log() << Verbose(3) << "Translation vector is ";
    17541753      Log() << Verbose(0) << Translationvector <<  endl;
    … …  
    17811780  delete(AtomStack);
    17821781  delete[](ColorList);
     1782  delete[](matrix);
    17831783  DoLog(2) && (Log() << Verbose(2) << "End of ScanForPeriodicCorrection." << endl);
    17841784};

  • src/molecule_geometry.cpp

    r5ec8e3 r56fb09  
    1919#include "World.hpp"
    2020#include "Plane.hpp"
    21 #include "Matrix.hpp"
    22 #include "Box.hpp"
    2321#include <boost/foreach.hpp>
    2422
    … …  
    3533  const Vector *Center = DetermineCenterOfAll();
    3634  const Vector *CenterBox = DetermineCenterOfBox();
    37   Box &domain = World::getInstance().getDomain();
     35  double * const cell_size = World::getInstance().getDomain();
     36  double *M = ReturnFullMatrixforSymmetric(cell_size);
     37  double *Minv = InverseMatrix(M);
    3838
    3939  // go through all atoms
    4040  ActOnAllVectors( &Vector::SubtractVector, *Center);
    4141  ActOnAllVectors( &Vector::SubtractVector, *CenterBox);
    42   BOOST_FOREACH(atom* iter, atoms){
    43     *iter->node = domain.WrapPeriodically(*iter->node);
    44   }
    45 
     42  ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
     43
     44  delete[](M);
     45  delete[](Minv);
    4646  delete(Center);
    4747  delete(CenterBox);
    … …  
    5656{
    5757  bool status = true;
    58   Box &domain = World::getInstance().getDomain();
     58  double * const cell_size = World::getInstance().getDomain();
     59  double *M = ReturnFullMatrixforSymmetric(cell_size);
     60  double *Minv = InverseMatrix(M);
    5961
    6062  // go through all atoms
    61   BOOST_FOREACH(atom* iter, atoms){
    62     *iter->node = domain.WrapPeriodically(*iter->node);
    63   }
    64 
     63  ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
     64
     65  delete[](M);
     66  delete[](Minv);
    6567  return status;
    6668};
    … …  
    151153{
    152154  Vector *a = new Vector(0.5,0.5,0.5);
    153   const Matrix &M = World::getInstance().getDomain().getM();
    154   (*a) *= M;
     155
     156  const double *cell_size = World::getInstance().getDomain();
     157  double *M = ReturnFullMatrixforSymmetric(cell_size);
     158  a->MatrixMultiplication(M);
     159  delete[](M);
     160
    155161  return a;
    156162};
    … …  
    238244void molecule::TranslatePeriodically(const Vector *trans)
    239245{
    240   Box &domain = World::getInstance().getDomain();
     246  double * const cell_size = World::getInstance().getDomain();
     247  double *M = ReturnFullMatrixforSymmetric(cell_size);
     248  double *Minv = InverseMatrix(M);
    241249
    242250  // go through all atoms
    243251  ActOnAllVectors( &Vector::AddVector, *trans);
    244   BOOST_FOREACH(atom* iter, atoms){
    245     *iter->node = domain.WrapPeriodically(*iter->node);
    246   }
    247 
     252  ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
     253
     254  delete[](M);
     255  delete[](Minv);
    248256};
    249257
    … …  
    256264  OBSERVE;
    257265  Plane p(*n,0);
    258   BOOST_FOREACH(atom* iter, atoms ){
     266  BOOST_FOREACH( atom* iter, atoms ){
    259267    (*iter->node) = p.mirrorVector(*iter->node);
    260268  }
    … …  
    266274void molecule::DeterminePeriodicCenter(Vector &center)
    267275{
    268   const Matrix &matrix = World::getInstance().getDomain().getM();
    269   const Matrix &inversematrix = World::getInstance().getDomain().getM();
     276  double * const cell_size = World::getInstance().getDomain();
     277  double *matrix = ReturnFullMatrixforSymmetric(cell_size);
     278  double *inversematrix = InverseMatrix(matrix);
    270279  double tmp;
    271280  bool flag;
    … …  
    279288      if ((*iter)->type->Z != 1) {
    280289#endif
    281         Testvector = inversematrix * (*iter)->x;
     290        Testvector = (*iter)->x;
     291        Testvector.MatrixMultiplication(inversematrix);
    282292        Translationvector.Zero();
    283293        for (BondList::const_iterator Runner = (*iter)->ListOfBonds.begin(); Runner != (*iter)->ListOfBonds.end(); (++Runner)) {
    … …  
    296306        }
    297307        Testvector += Translationvector;
    298         Testvector *= matrix;
     308        Testvector.MatrixMultiplication(matrix);
    299309        Center += Testvector;
    300310        Log() << Verbose(1) << "vector is: " << Testvector << endl;
    … …  
    303313        for (BondList::const_iterator Runner = (*iter)->ListOfBonds.begin(); Runner != (*iter)->ListOfBonds.end(); (++Runner)) {
    304314          if ((*Runner)->GetOtherAtom((*iter))->type->Z == 1) {
    305             Testvector = inversematrix * (*Runner)->GetOtherAtom((*iter))->x;
     315            Testvector = (*Runner)->GetOtherAtom((*iter))->x;
     316            Testvector.MatrixMultiplication(inversematrix);
    306317            Testvector += Translationvector;
    307             Testvector *= matrix;
     318            Testvector.MatrixMultiplication(matrix);
    308319            Center += Testvector;
    309320            Log() << Verbose(1) << "Hydrogen vector is: " << Testvector << endl;
    … …  
    314325    }
    315326  } while (!flag);
     327  delete[](matrix);
     328  delete[](inversematrix);
    316329
    317330  Center.Scale(1./static_cast<double>(getAtomCount()));
    … …  
    375388    DoLog(1) && (Log() << Verbose(1) << "Transforming molecule into PAS ... ");
    376389    // the eigenvectors specify the transformation matrix
    377     Matrix M = Matrix(evec->data);
    378     BOOST_FOREACH(atom* iter, atoms){
    379       (*iter->node) *= M;
    380     }
     390    ActOnAllVectors( &Vector::MatrixMultiplication, (const double *) evec->data );
    381391    DoLog(0) && (Log() << Verbose(0) << "done." << endl);
    382392

  • src/molecule_graph.cpp

    r5ec8e3 r56fb09  
    2424#include "Helpers/fast_functions.hpp"
    2525#include "Helpers/Assert.hpp"
    26 #include "Matrix.hpp"
    27 #include "Box.hpp"
    2826
    2927
    … …  
    123121  LinkedCell *LC = NULL;
    124122  bool free_BG = false;
    125   Box &domain = World::getInstance().getDomain();
     123  double * const cell_size = World::getInstance().getDomain();
    126124
    127125  if (BG == NULL) {
    … …  
    180178                          //Log() << Verbose(1) << "Checking distance " << OtherWalker->x.PeriodicDistanceSquared(&(Walker->x), cell_size) << " against typical bond length of " << bonddistance*bonddistance << "." << endl;
    181179                          (BG->*minmaxdistance)(Walker, OtherWalker, MinDistance, MaxDistance, IsAngstroem);
    182                           const double distance = domain.periodicDistanceSquared(OtherWalker->x,Walker->x);
     180                          const double distance = OtherWalker->x.PeriodicDistanceSquared(Walker->x,cell_size);
    183181                          const bool status = (distance <= MaxDistance * MaxDistance) && (distance >= MinDistance * MinDistance);
    184182//                          Log() << Verbose(1) << "MinDistance is " << MinDistance << " and MaxDistance is " << MaxDistance << "." << endl;

  • src/moleculelist.cpp

    r5ec8e3 r56fb09  
    2424#include "periodentafel.hpp"
    2525#include "Helpers/Assert.hpp"
    26 #include "Matrix.hpp"
    27 #include "Box.hpp"
    2826
    2927#include "Helpers/Assert.hpp"
    … …  
    633631  int FragmentCounter = 0;
    634632  ofstream output;
    635   Matrix cell_size = World::getInstance().getDomain().getM();
    636   Matrix cell_size_backup = cell_size;
    637 
     633  double cell_size_backup[6];
     634  double * const cell_size = World::getInstance().getDomain();
     635
     636  // backup cell_size
     637  for (int i=0;i<6;i++)
     638    cell_size_backup[i] = cell_size[i];
    638639  // store the fragments as config and as xyz
    639640  for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
    … …  
    673674    (*ListRunner)->CenterEdge(&BoxDimension);
    674675    (*ListRunner)->SetBoxDimension(&BoxDimension); // update Box of atoms by boundary
     676    int j = -1;
    675677    for (int k = 0; k < NDIM; k++) {
     678      j += k + 1;
    676679      BoxDimension[k] = 2.5 * (World::getInstance().getConfig()->GetIsAngstroem() ? 1. : 1. / AtomicLengthToAngstroem);
    677       cell_size.at(k,k) = BoxDimension[k] * 2.;
    678     }
    679     World::getInstance().setDomain(cell_size);
     680      cell_size[j] = BoxDimension[k] * 2.;
     681    }
    680682    (*ListRunner)->Translate(&BoxDimension);
    681683
    … …  
    723725
    724726  // restore cell_size
    725   World::getInstance().setDomain(cell_size_backup);
     727  for (int i=0;i<6;i++)
     728    cell_size[i] = cell_size_backup[i];
    726729
    727730  return result;
    … …  
    884887  // center at set box dimensions
    885888  mol->CenterEdge(&center);
    886   Matrix domain;
    887   for(int i =0;i<NDIM;++i)
    888     domain.at(i,i) = center[i];
    889   World::getInstance().setDomain(domain);
     889  World::getInstance().getDomain()[0] = center[0];
     890  World::getInstance().getDomain()[1] = 0;
     891  World::getInstance().getDomain()[2] = center[1];
     892  World::getInstance().getDomain()[3] = 0;
     893  World::getInstance().getDomain()[4] = 0;
     894  World::getInstance().getDomain()[5] = center[2];
    890895  insert(mol);
    891896}

  • src/unittests/LineUnittest.cpp

    r5ec8e3 r56fb09  
    1717
    1818#include <iostream>
    19 #include <cmath>
    2019
    2120using namespace std;
    … …  
    353352  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
    354353}
    355 
    356 void LineUnittest::sphereIntersectionTest(){
    357   {
    358     std::vector<Vector> res = la1->getSphereIntersections();
    359     CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
    360     CPPUNIT_ASSERT(testDirection(res[0],e1));
    361     CPPUNIT_ASSERT(testDirection(res[1],e1));
    362     CPPUNIT_ASSERT(res[0]!=res[1]);
    363   }
    364 
    365   {
    366     std::vector<Vector> res = la2->getSphereIntersections();
    367     CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
    368     CPPUNIT_ASSERT(testDirection(res[0],e2));
    369     CPPUNIT_ASSERT(testDirection(res[1],e2));
    370     CPPUNIT_ASSERT(res[0]!=res[1]);
    371   }
    372 
    373   {
    374     std::vector<Vector> res = la3->getSphereIntersections();
    375     CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
    376     CPPUNIT_ASSERT(testDirection(res[0],e3));
    377     CPPUNIT_ASSERT(testDirection(res[1],e3));
    378     CPPUNIT_ASSERT(res[0]!=res[1]);
    379   }
    380 
    381   {
    382     std::vector<Vector> res = lp1->getSphereIntersections();
    383     CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
    384     CPPUNIT_ASSERT((res[0]==e1) || (res[0]==e2));
    385     CPPUNIT_ASSERT((res[1]==e1) || (res[1]==e2));
    386     CPPUNIT_ASSERT(res[0]!=res[1]);
    387   }
    388 
    389   {
    390     std::vector<Vector> res = lp2->getSphereIntersections();
    391     CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
    392     CPPUNIT_ASSERT((res[0]==e2) || (res[0]==e3));
    393     CPPUNIT_ASSERT((res[1]==e2) || (res[1]==e3));
    394     CPPUNIT_ASSERT(res[0]!=res[1]);
    395   }
    396 
    397   {
    398     std::vector<Vector> res = lp3->getSphereIntersections();
    399     CPPUNIT_ASSERT_EQUAL(res.size(),(size_t)2);
    400     CPPUNIT_ASSERT((res[0]==e3) || (res[0]==e1));
    401     CPPUNIT_ASSERT((res[1]==e3) || (res[1]==e1));
    402     CPPUNIT_ASSERT(res[0]!=res[1]);
    403   }
    404 }

  • src/unittests/LineUnittest.hpp

    r5ec8e3 r56fb09  
    2222  CPPUNIT_TEST ( intersectionTest );
    2323  CPPUNIT_TEST ( rotationTest );
    24   CPPUNIT_TEST ( sphereIntersectionTest );
    2524  CPPUNIT_TEST_SUITE_END();
    2625
    … …  
    3433  void intersectionTest();
    3534  void rotationTest();
    36   void sphereIntersectionTest();
    3735
    3836private:

  • src/unittests/Makefile.am

    r5ec8e3 r56fb09  
    3838  periodentafelTest \
    3939  PlaneUnittest \
    40   ShapeUnittest \
    4140  SingletonTest \
    4241  StackClassUnitTest \
    … …  
    8483  periodentafelTest.cpp \
    8584  PlaneUnittest.cpp \
    86   ShapeUnittest.cpp \
    8785  SingletonTest.cpp \
    8886  stackclassunittest.cpp \
    … …  
    212210PlaneUnittest_LDADD = ${ALLLIBS}
    213211
    214 ShapeUnittest_SOURCES = UnitTestMain.cpp ShapeUnittest.cpp ShapeUnittest.hpp
    215 ShapeUnittest_LDADD = ${ALLLIBS}
    216 
    217212SingletonTest_SOURCES = UnitTestMain.cpp SingletonTest.cpp SingletonTest.hpp
    218213SingletonTest_LDADD = ${ALLLIBS} $(BOOST_LIB) ${BOOST_THREAD_LIB}

  • src/unittests/PlaneUnittest.cpp

    r5ec8e3 r56fb09  
    1111#include <cppunit/extensions/TestFactoryRegistry.h>
    1212#include <cppunit/ui/text/TestRunner.h>
    13 
    14 #include <cmath>
    1513
    1614#ifdef HAVE_TESTRUNNER

  • src/unittests/linearsystemofequationsunittest.cpp

    r5ec8e3 r56fb09  
    1212#include <cppunit/extensions/TestFactoryRegistry.h>
    1313#include <cppunit/ui/text/TestRunner.h>
    14 #include <cmath>
    1514
    1615#include "linearsystemofequationsunittest.hpp"

  • src/unittests/vectorunittest.cpp

    r5ec8e3 r56fb09  
    2121#include "Plane.hpp"
    2222#include "Exceptions/LinearDependenceException.hpp"
    23 #include "Matrix.hpp"
    2423
    2524#ifdef HAVE_TESTRUNNER
    … …  
    215214  CPPUNIT_ASSERT(testVector.ScalarProduct(three) < MYEPSILON);
    216215}
     216
     217
     218/**
     219 * UnitTest for Vector::IsInParallelepiped().
     220 */
     221void VectorTest::IsInParallelepipedTest()
     222{
     223  double parallelepiped[NDIM*NDIM];
     224  parallelepiped[0] = 1;
     225  parallelepiped[1] = 0;
     226  parallelepiped[2] = 0;
     227  parallelepiped[3] = 0;
     228  parallelepiped[4] = 1;
     229  parallelepiped[5] = 0;
     230  parallelepiped[6] = 0;
     231  parallelepiped[7] = 0;
     232  parallelepiped[8] = 1;
     233
     234  fixture = zero;
     235  CPPUNIT_ASSERT_EQUAL( false, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
     236  fixture = Vector(2.5,2.5,2.5);
     237  CPPUNIT_ASSERT_EQUAL( true, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
     238  fixture = Vector(1.,1.,1.);
     239  CPPUNIT_ASSERT_EQUAL( false, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
     240  fixture = Vector(3.5,3.5,3.5);
     241  CPPUNIT_ASSERT_EQUAL( false, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
     242  fixture = Vector(2.,2.,2.);
     243  CPPUNIT_ASSERT_EQUAL( true, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
     244  fixture = Vector(2.,3.,2.);
     245  CPPUNIT_ASSERT_EQUAL( true, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
     246  fixture = Vector(-2.,2.,-1.);
     247  CPPUNIT_ASSERT_EQUAL( false, fixture.IsInParallelepiped(Vector(2.,2.,2.), parallelepiped) );
     248}
     249

  • src/unittests/vectorunittest.hpp

    r5ec8e3 r56fb09  
    2727    CPPUNIT_TEST ( ProjectionTest );
    2828    CPPUNIT_TEST ( NormalsTest );
     29    CPPUNIT_TEST ( IsInParallelepipedTest );
    2930    CPPUNIT_TEST_SUITE_END();
    3031
    … …  
    4445    void LineIntersectionTest();
    4546    void VectorRotationTest();
     47    void IsInParallelepipedTest();
    4648
    4749private:

  • src/vector.cpp

    r5ec8e3 r56fb09  
    1212#include "Helpers/Assert.hpp"
    1313#include "Helpers/fast_functions.hpp"
    14 #include "Exceptions/MathException.hpp"
    1514
    1615#include <iostream>
    17 #include <gsl/gsl_blas.h>
    18 
    1916
    2017using namespace std;
    … …  
    3734{
    3835  content = gsl_vector_alloc(NDIM);
    39   gsl_vector_memcpy(content, src.content);
     36  gsl_vector_set(content,0,src[0]);
     37  gsl_vector_set(content,1,src[1]);
     38  gsl_vector_set(content,2,src[2]);
    4039}
    4140
    … …  
    5049};
    5150
    52 Vector::Vector(gsl_vector *_content) :
    53   content(_content)
    54 {}
    55 
    5651/**
    5752 * Assignment operator
    … …  
    6055  // check for self assignment
    6156  if(&src!=this){
    62     gsl_vector_memcpy(content, src.content);
     57    gsl_vector_set(content,0,src[0]);
     58    gsl_vector_set(content,1,src[1]);
     59    gsl_vector_set(content,2,src[2]);
    6360  }
    6461  return *this;
    … …  
    9794}
    9895
     96/** Calculates distance between this and another vector in a periodic cell.
     97 * \param *y array to second vector
     98 * \param *cell_size 6-dimensional array with (xx, xy, yy, xz, yz, zz) entries specifying the periodic cell
     99 * \return \f$| x - y |\f$
     100 */
     101double Vector::PeriodicDistance(const Vector &y, const double * const cell_size) const
     102{
     103  double res = distance(y), tmp, matrix[NDIM*NDIM];
     104    Vector Shiftedy, TranslationVector;
     105    int N[NDIM];
     106    matrix[0] = cell_size[0];
     107    matrix[1] = cell_size[1];
     108    matrix[2] = cell_size[3];
     109    matrix[3] = cell_size[1];
     110    matrix[4] = cell_size[2];
     111    matrix[5] = cell_size[4];
     112    matrix[6] = cell_size[3];
     113    matrix[7] = cell_size[4];
     114    matrix[8] = cell_size[5];
     115    // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
     116    for (N[0]=-1;N[0]<=1;N[0]++)
     117      for (N[1]=-1;N[1]<=1;N[1]++)
     118        for (N[2]=-1;N[2]<=1;N[2]++) {
     119          // create the translation vector
     120          TranslationVector.Zero();
     121          for (int i=NDIM;i--;)
     122            TranslationVector[i] = (double)N[i];
     123          TranslationVector.MatrixMultiplication(matrix);
     124          // add onto the original vector to compare with
     125          Shiftedy = y + TranslationVector;
     126          // get distance and compare with minimum so far
     127          tmp = distance(Shiftedy);
     128          if (tmp < res) res = tmp;
     129        }
     130    return (res);
     131};
     132
     133/** Calculates distance between this and another vector in a periodic cell.
     134 * \param *y array to second vector
     135 * \param *cell_size 6-dimensional array with (xx, xy, yy, xz, yz, zz) entries specifying the periodic cell
     136 * \return \f$| x - y |^2\f$
     137 */
     138double Vector::PeriodicDistanceSquared(const Vector &y, const double * const cell_size) const
     139{
     140  double res = DistanceSquared(y), tmp, matrix[NDIM*NDIM];
     141    Vector Shiftedy, TranslationVector;
     142    int N[NDIM];
     143    matrix[0] = cell_size[0];
     144    matrix[1] = cell_size[1];
     145    matrix[2] = cell_size[3];
     146    matrix[3] = cell_size[1];
     147    matrix[4] = cell_size[2];
     148    matrix[5] = cell_size[4];
     149    matrix[6] = cell_size[3];
     150    matrix[7] = cell_size[4];
     151    matrix[8] = cell_size[5];
     152    // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
     153    for (N[0]=-1;N[0]<=1;N[0]++)
     154      for (N[1]=-1;N[1]<=1;N[1]++)
     155        for (N[2]=-1;N[2]<=1;N[2]++) {
     156          // create the translation vector
     157          TranslationVector.Zero();
     158          for (int i=NDIM;i--;)
     159            TranslationVector[i] = (double)N[i];
     160          TranslationVector.MatrixMultiplication(matrix);
     161          // add onto the original vector to compare with
     162          Shiftedy = y + TranslationVector;
     163          // get distance and compare with minimum so far
     164          tmp = DistanceSquared(Shiftedy);
     165          if (tmp < res) res = tmp;
     166        }
     167    return (res);
     168};
     169
     170/** Keeps the vector in a periodic cell, defined by the symmetric \a *matrix.
     171 * \param *out ofstream for debugging messages
     172 * Tries to translate a vector into each adjacent neighbouring cell.
     173 */
     174void Vector::KeepPeriodic(const double * const matrix)
     175{
     176  //  int N[NDIM];
     177  //  bool flag = false;
     178    //vector Shifted, TranslationVector;
     179  //  Log() << Verbose(1) << "Begin of KeepPeriodic." << endl;
     180  //  Log() << Verbose(2) << "Vector is: ";
     181  //  Output(out);
     182  //  Log() << Verbose(0) << endl;
     183    InverseMatrixMultiplication(matrix);
     184    for(int i=NDIM;i--;) { // correct periodically
     185      if (at(i) < 0) {  // get every coefficient into the interval [0,1)
     186        at(i) += ceil(at(i));
     187      } else {
     188        at(i) -= floor(at(i));
     189      }
     190    }
     191    MatrixMultiplication(matrix);
     192  //  Log() << Verbose(2) << "New corrected vector is: ";
     193  //  Output(out);
     194  //  Log() << Verbose(0) << endl;
     195  //  Log() << Verbose(1) << "End of KeepPeriodic." << endl;
     196};
     197
    99198/** Calculates scalar product between this and another vector.
    100199 * \param *y array to second vector
    … …  
    104203{
    105204  double res = 0.;
    106   gsl_blas_ddot(content, y.content, &res);
     205  for (int i=NDIM;i--;)
     206    res += at(i)*y[i];
    107207  return (res);
    108208};
    … …  
    404504};
    405505
    406 void Vector::ScaleAll(const Vector &factor){
    407   gsl_vector_mul(content, factor.content);
    408 }
    409506
    410507
    411508void Vector::Scale(const double factor)
    412509{
    413   gsl_vector_scale(content,factor);
     510  for (int i=NDIM;i--;)
     511    at(i) *= factor;
     512};
     513
     514/** Given a box by its matrix \a *M and its inverse *Minv the vector is made to point within that box.
     515 * \param *M matrix of box
     516 * \param *Minv inverse matrix
     517 */
     518void Vector::WrapPeriodically(const double * const M, const double * const Minv)
     519{
     520  MatrixMultiplication(Minv);
     521  // truncate to [0,1] for each axis
     522  for (int i=0;i<NDIM;i++) {
     523    //at(i) += 0.5;  // set to center of box
     524    while (at(i) >= 1.)
     525      at(i) -= 1.;
     526    while (at(i) < 0.)
     527      at(i) += 1.;
     528  }
     529  MatrixMultiplication(M);
    414530};
    415531
    … …  
    430546  return make_pair(res,helper);
    431547}
     548
     549/** Do a matrix multiplication.
     550 * \param *matrix NDIM_NDIM array
     551 */
     552void Vector::MatrixMultiplication(const double * const M)
     553{
     554  Vector tmp;
     555  // do the matrix multiplication
     556  for(int i=NDIM;i--;)
     557    tmp[i] = M[i]*at(0)+M[i+3]*at(1)+M[i+6]*at(2);
     558
     559  (*this) = tmp;
     560};
     561
     562/** Do a matrix multiplication with the \a *A' inverse.
     563 * \param *matrix NDIM_NDIM array
     564 */
     565bool Vector::InverseMatrixMultiplication(const double * const A)
     566{
     567  double B[NDIM*NDIM];
     568  double detA = RDET3(A);
     569  double detAReci;
     570
     571  // calculate the inverse B
     572  if (fabs(detA) > MYEPSILON) {;  // RDET3(A) yields precisely zero if A irregular
     573    detAReci = 1./detA;
     574    B[0] =  detAReci*RDET2(A[4],A[5],A[7],A[8]);    // A_11
     575    B[1] = -detAReci*RDET2(A[1],A[2],A[7],A[8]);    // A_12
     576    B[2] =  detAReci*RDET2(A[1],A[2],A[4],A[5]);    // A_13
     577    B[3] = -detAReci*RDET2(A[3],A[5],A[6],A[8]);    // A_21
     578    B[4] =  detAReci*RDET2(A[0],A[2],A[6],A[8]);    // A_22
     579    B[5] = -detAReci*RDET2(A[0],A[2],A[3],A[5]);    // A_23
     580    B[6] =  detAReci*RDET2(A[3],A[4],A[6],A[7]);    // A_31
     581    B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]);    // A_32
     582    B[8] =  detAReci*RDET2(A[0],A[1],A[3],A[4]);    // A_33
     583
     584    MatrixMultiplication(B);
     585
     586    return true;
     587  } else {
     588    return false;
     589  }
     590};
     591
    432592
    433593/** Creates this vector as the b y *factors' components scaled linear combination of the given three.
    … …  
    519679void Vector::AddVector(const Vector &y)
    520680{
    521   gsl_vector_add(content, y.content);
     681  for(int i=NDIM;i--;)
     682    at(i) += y[i];
    522683}
    523684
    … …  
    527688void Vector::SubtractVector(const Vector &y)
    528689{
    529   gsl_vector_sub(content, y.content);
     690  for(int i=NDIM;i--;)
     691    at(i) -= y[i];
     692}
     693
     694/**
     695 * Checks whether this vector is within the parallelepiped defined by the given three vectors and
     696 * their offset.
     697 *
     698 * @param offest for the origin of the parallelepiped
     699 * @param three vectors forming the matrix that defines the shape of the parallelpiped
     700 */
     701bool Vector::IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const
     702{
     703  Vector a = (*this)-offset;
     704  a.InverseMatrixMultiplication(parallelepiped);
     705  bool isInside = true;
     706
     707  for (int i=NDIM;i--;)
     708    isInside = isInside && ((a[i] <= 1) && (a[i] >= 0));
     709
     710  return isInside;
    530711}
    531712

  • src/vector.hpp

    r5ec8e3 r56fb09  
    1313#include <iostream>
    1414#include <gsl/gsl_vector.h>
     15#include <gsl/gsl_multimin.h>
    1516
    1617#include <memory>
    … …  
    2324
    2425class Vector;
    25 class Matrix;
    2626
    2727typedef std::vector<Vector> pointset;
    … …  
    3131 */
    3232class Vector : public Space{
    33   friend Vector operator*(const Matrix&,const Vector&);
    3433public:
     34
    3535  Vector();
    3636  Vector(const double x1, const double x2, const double x3);
    … …  
    4242  double DistanceSquared(const Vector &y) const;
    4343  double DistanceToSpace(const Space& space) const;
     44  double PeriodicDistance(const Vector &y, const double * const cell_size) const;
     45  double PeriodicDistanceSquared(const Vector &y, const double * const cell_size) const;
    4446  double ScalarProduct(const Vector &y) const;
    4547  double Angle(const Vector &y) const;
    … …  
    5658  Vector Projection(const Vector &y) const;
    5759  void ScaleAll(const double *factor);
    58   void ScaleAll(const Vector &factor);
    5960  void Scale(const double factor);
     61  void MatrixMultiplication(const double * const M);
     62  bool InverseMatrixMultiplication(const double * const M);
     63  void KeepPeriodic(const double * const matrix);
    6064  bool GetOneNormalVector(const Vector &x1);
    6165  bool MakeNormalTo(const Vector &y1);
     66  bool IsInParallelepiped(const Vector &offset, const double * const parallelepiped) const;
     67  void WrapPeriodically(const double * const M, const double * const Minv);
    6268  std::pair<Vector,Vector> partition(const Vector&) const;
    6369  std::pair<pointset,Vector> partition(const pointset&) const;
    … …  
    98104
    99105private:
    100   Vector(gsl_vector *);
    101106  gsl_vector *content;
    102107

  • src/vector_ops.cpp

    r5ec8e3 r56fb09  
    2323#include <gsl/gsl_permutation.h>
    2424#include <gsl/gsl_vector.h>
    25 #include <gsl/gsl_multimin.h>
    2625
    2726/**

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

    r5ec8e3 r56fb09  
    4747BoxLength                       # (Length of a unit cell)
    48481       
    49 0       1       
     490       0       
    50500       0       2       
    5151
Note: See TracChangeset for help on using the changeset viewer.