Changes in / [56fb09:5ec8e3]


Ignore:
Files:
17 added
51 edited

Legend:

Unmodified
Added
Removed

  • src/Actions/AnalysisAction/PairCorrelationAction.cpp

    r56fb09 r5ec8e3  
    6767  dialog = UIFactory::getInstance().makeDialog();
    6868  if (type == "P")
    69     dialog->queryVector("position", &Point, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription("position"));
     69    dialog->queryVector("position", &Point, 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

    r56fb09 r5ec8e3  
    4141
    4242  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"));
    4444  cout << "pre-dialog" << endl;
    4545

  • src/Actions/MoleculeAction/FillWithMoleculeAction.cpp

    r56fb09 r5ec8e3  
    6262
    6363  dialog->queryString(NAME, &filename, MapOfActions::getInstance().getDescription(NAME));
    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"));
     64  dialog->queryVector("distances", &distances, false, MapOfActions::getInstance().getDescription("distances"));
     65  dialog->queryVector("lengths", &lengths, 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

    r56fb09 r5ec8e3  
    5555  bool periodic = false;
    5656
    57   dialog->queryVector(NAME, &x, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription(NAME));
     57  dialog->queryVector(NAME, &x, 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

    r56fb09 r5ec8e3  
    4141  int j=0;
    4242
    43   dialog->queryVector(NAME, &boundary, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription(NAME));
     43  dialog->queryVector(NAME, &boundary, false, MapOfActions::getInstance().getDescription(NAME));
    4444
    4545  if(dialog->display()) {
    … …  
    5959    }
    6060    // set new box size
    61     double * const cell_size = World::getInstance().getDomain();
     61    double * const cell_size = new double[6];
    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;
    6971    // translate all atoms, such that Min is aty (0,0,0)
    7072    AtomRunner = AllAtoms.begin();

  • src/Actions/WorldAction/CenterInBoxAction.cpp

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

  • src/Actions/WorldAction/CenterOnEdgeAction.cpp

    r56fb09 r5ec8e3  
    1313#include "vector.hpp"
    1414#include "World.hpp"
     15#include "Matrix.hpp"
    1516
    1617#include <iostream>
    … …  
    3738  Vector Min;
    3839  Vector Max;
    39   int j=0;
    4040
    4141  dialog->queryEmpty(NAME, MapOfActions::getInstance().getDescription(NAME));
    … …  
    5757    }
    5858    // set new box size
    59     double * const cell_size = World::getInstance().getDomain();
    60     for (j=0;j<6;j++)
    61       cell_size[j] = 0.;
    62     j=-1;
     59    Matrix domain;
    6360    for (int i=0;i<NDIM;i++) {
    64       j += i+1;
    65       cell_size[j] = (Max[i]-Min[i]);
     61      double tmp = Max[i]-Min[i];
     62      tmp = fabs(tmp)>=1. ? tmp : 1.0;
     63      domain.at(i,i) = tmp;
    6664    }
     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

    r56fb09 r5ec8e3  
    1212#include "verbose.hpp"
    1313#include "World.hpp"
     14#include "Box.hpp"
     15#include "Matrix.hpp"
    1416
    1517#include <iostream>
    … …  
    3436  Dialog *dialog = UIFactory::getInstance().makeDialog();
    3537
    36   double * cell_size = World::getInstance().getDomain();
     38  Box& cell_size = World::getInstance().getDomain();
    3739  dialog->queryBox(NAME, &cell_size, MapOfActions::getInstance().getDescription(NAME));
    3840
    3941  if(dialog->display()) {
    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);
     42    DoLog(0) && (Log() << Verbose(0) << "Setting box domain to " << cell_size.getM() << endl);
    4143    delete dialog;
    4244    return Action::success;

  • src/Actions/WorldAction/RemoveSphereOfAtomsAction.cpp

    r56fb09 r5ec8e3  
    4141
    4242  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"));
    4444
    4545  if(dialog->display()) {

  • src/Actions/WorldAction/RepeatBoxAction.cpp

    r56fb09 r5ec8e3  
    1313#include "molecule.hpp"
    1414#include "vector.hpp"
     15#include "Matrix.hpp"
    1516#include "verbose.hpp"
    1617#include "World.hpp"
     18#include "Box.hpp"
    1719
    1820#include <iostream>
    … …  
    4648  MoleculeListClass *molecules = World::getInstance().getMolecules();
    4749
    48   dialog->queryVector(NAME, &Repeater, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription(NAME));
     50  dialog->queryVector(NAME, &Repeater, false, MapOfActions::getInstance().getDescription(NAME));
    4951  //dialog->queryMolecule("molecule-by-id", &mol,MapOfActions::getInstance().getDescription("molecule-by-id"));
    5052  vector<molecule *> AllMolecules;
    … …  
    5961  if(dialog->display()) {
    6062    (cout << "Repeating box " << Repeater << " times for (x,y,z) axis." << endl);
    61     double * const cell_size = World::getInstance().getDomain();
    62     double *M = ReturnFullMatrixforSymmetric(cell_size);
     63    Matrix M = World::getInstance().getDomain().getM();
     64    Matrix newM = M;
    6365    Vector x,y;
    6466    int n[NDIM];
     67    Matrix repMat;
    6568    for (int axis = 0; axis < NDIM; axis++) {
    6669      Repeater[axis] = floor(Repeater[axis]);
    … …  
    6972        Repeater[axis] = 1;
    7073      }
    71       cell_size[(abs(axis+1) == 2) ? 2 : ((abs(axis+2) == 3) ? 5 : 0)] *= Repeater[axis];
     74      repMat.at(axis,axis) = Repeater[axis];
    7275    }
     76    newM *= repMat;
     77    World::getInstance().setDomain(newM);
    7378
    7479    molecule *newmol = NULL;
    … …  
    98103                DoeLog(1) && (eLog()<< Verbose(1) << "AtomCount " << count << " is not equal to number of atoms in molecule " << j << "!" << endl);
    99104              x = y;
    100               x.MatrixMultiplication(M);
     105              x *= M;
    101106              newmol = World::getInstance().createMolecule();
    102107              molecules->insert(newmol);
    … …  
    118123      }
    119124    }
    120     delete(M);
    121125    delete dialog;
    122126    return Action::success;

  • src/Actions/WorldAction/ScaleBoxAction.cpp

    r56fb09 r5ec8e3  
    1414#include "verbose.hpp"
    1515#include "World.hpp"
     16#include "Box.hpp"
     17#include "Matrix.hpp"
    1618
    1719#include <iostream>
    … …  
    3739  Vector Scaler;
    3840  double x[NDIM];
    39   int j=0;
    4041
    41   dialog->queryVector(NAME, &Scaler, World::getInstance().getDomain(), false, MapOfActions::getInstance().getDescription(NAME));
     42  dialog->queryVector(NAME, &Scaler, false, MapOfActions::getInstance().getDescription(NAME));
    4243
    4344  if(dialog->display()) {
    … …  
    4950      (*AtomRunner)->x.ScaleAll(x);
    5051    }
    51     j = -1;
    52     double * const cell_size = World::getInstance().getDomain();
     52
     53    Matrix M = World::getInstance().getDomain().getM();
     54    Matrix scale;
     55
    5356    for (int i=0;i<NDIM;i++) {
    54       j += i+1;
    55       cell_size[j]*=x[i];
     57      scale.at(i,i) = x[i];
    5658    }
     59    M *= scale;
     60    World::getInstance().setDomain(M);
    5761
    5862    delete dialog;

  • src/Line.cpp

    r56fb09 r5ec8e3  
    215215}
    216216
     217std::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
    217236Line makeLineThrough(const Vector &x1, const Vector &x2){
    218237  if(x1==x2){

  • src/Line.hpp

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

  • src/Makefile.am

    r56fb09 r5ec8e3  
    3838  vector.cpp
    3939                           
    40 LINALGHEADER = gslmatrix.hpp \
     40LINALGHEADER = \
     41  gslmatrix.hpp \
    4142  gslvector.hpp \
    4243  linearsystemofequations.hpp \
    … …  
    7576  Actions/MethodAction.hpp \
    7677  Actions/Process.hpp
    77  
     78
     79EXCEPTIONSOURCE = Exceptions/CustomException.cpp \
     80                                  Exceptions/LinearDependenceException.cpp \
     81                                  Exceptions/MathException.cpp \
     82                                  Exceptions/NotInvertibleException.cpp \
     83                                  Exceptions/SkewException.cpp \
     84                                  Exceptions/ZeroVectorException.cpp
     85                                 
     86EXCEPTIONHEADER = Exceptions/CustomException.hpp \
     87                                  Exceptions/LinearDependenceException.hpp \
     88                                  Exceptions/MathException.hpp \
     89                                  Exceptions/NotInvertibleException.hpp \
     90                                  Exceptions/SkewException.hpp \
     91                                  Exceptions/ZeroVectorException.hpp
    7892
    7993PARSERSOURCE = \
    … …  
    101115  Patterns/Observer.hpp \
    102116  Patterns/Singleton.hpp
     117 
     118SHAPESOURCE = \
     119  Shapes/BaseShapes.cpp \
     120  Shapes/Shape.cpp \
     121  Shapes/ShapeOps.cpp
     122SHAPEHEADER = \
     123  Shapes/BaseShapes.hpp \
     124  Shapes/Shape.hpp \
     125  Shapes/ShapeOps.hpp
     126 
    103127
    104128QTUIMOC_HEADER = UIElements/QT4/QTDialog.hpp \
    … …  
    134158  Descriptors/MoleculeNameDescriptor.hpp \
    135159  Descriptors/MoleculePtrDescriptor.hpp
    136                                    
    137 EXCEPTIONSOURCE = Exceptions/CustomException.cpp \
    138                                   Exceptions/LinearDependenceException.cpp \
    139                                   Exceptions/MathException.cpp \
    140                                   Exceptions/SkewException.cpp \
    141                                   Exceptions/ZeroVectorException.cpp
    142                                  
    143 EXCEPTIONHEADER = Exceptions/CustomException.hpp \
    144                                   Exceptions/LinearDependenceException.hpp \
    145                                   Exceptions/MathException.hpp \
    146                                   Exceptions/SkewException.hpp \
    147                                   Exceptions/ZeroVectorException.hpp
    148160                                 
    149161QTUISOURCE = ${QTUIMOC_TARGETS} \
    … …  
    165177  ${ACTIONSSOURCE} \
    166178  ${ATOMSOURCE} \
     179  ${EXCEPTIONSOURCE} \
    167180  ${PATTERNSOURCE} \
    168181  ${PARSERSOURCE} \
     182  ${SHAPESOURCE} \
    169183  ${DESCRIPTORSOURCE} \
    170184  ${HELPERSOURCE} \
    171   ${EXCEPTIONSOURCE} \
    172185  bond.cpp \
    173186  bondgraph.cpp \
    174187  boundary.cpp \
     188  Box.cpp \
    175189  CommandLineParser.cpp \
    176190  config.cpp \
    … …  
    188202  log.cpp \
    189203  logger.cpp \
     204  Matrix.cpp \
    190205  moleculelist.cpp \
    191206  molecule.cpp \
    … …  
    212227  ${ACTIONSHEADER} \
    213228  ${ATOMHEADER} \
     229  ${EXCEPTIONHEADER} \
    214230  ${PARSERHEADER} \
    215231  ${PATTERNHEADER} \
     232  ${SHAPEHEADER} \
    216233  ${DESCRIPTORHEADER} \
    217   ${EXCEPTIONHEADER} \
    218234  bond.hpp \
    219235  bondgraph.hpp \
    220236  boundary.hpp \
     237  Box.hpp \
    221238  CommandLineParser.hpp \
    222239  config.hpp \
    … …  
    236253  log.hpp \
    237254  logger.hpp \
     255  Matrix.hpp \
    238256  molecule.hpp \
    239257  molecule_template.hpp \

  • src/Parser/PcpParser.cpp

    r56fb09 r5ec8e3  
    2020#include "verbose.hpp"
    2121#include "World.hpp"
     22#include "Matrix.hpp"
     23#include "Box.hpp"
    2224
    2325/** Constructor of PcpParser.
    … …  
    210212  // Unit cell and magnetic field
    211213  ParseForParameter(verbose,FileBuffer, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
    212   double * const cell_size = World::getInstance().getDomain();
     214  double *cell_size = new double[6];
    213215  cell_size[0] = BoxLength[0];
    214216  cell_size[1] = BoxLength[3];
    … …  
    217219  cell_size[4] = BoxLength[7];
    218220  cell_size[5] = BoxLength[8];
     221  World::getInstance().setDomain(cell_size);
     222  delete[] cell_size;
    219223  //if (1) fprintf(stderr,"\n");
    220224
    … …  
    327331void PcpParser::save(std::ostream* file)
    328332{
    329   const double * const cell_size = World::getInstance().getDomain();
     333  const Matrix &domain = World::getInstance().getDomain().getM();
    330334  class ThermoStatContainer *Thermostats = World::getInstance().getThermostats();
    331335  if (!file->fail()) {
    … …  
    412416    *file << endl;
    413417    *file << "BoxLength\t\t\t# (Length of a unit cell)" << 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;
     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;
    417421    // FIXME
    418422    *file << endl;

  • src/UIElements/CommandLineUI/CommandLineDialog.cpp

    r56fb09 r5ec8e3  
    2727#include "verbose.hpp"
    2828#include "World.hpp"
     29#include "Box.hpp"
    2930
    3031#include "atom.hpp"
    … …  
    7374}
    7475
    75 void 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 
    79 void CommandLineDialog::queryBox(const char* title, double ** const cellSize, string _description) {
     76void CommandLineDialog::queryVector(const char* title, Vector *target, bool check, string _description) {
     77  registerQuery(new VectorCommandLineQuery(title,target,check, _description));
     78}
     79
     80void CommandLineDialog::queryBox(const char* title, Box* cellSize, string _description) {
    8081  registerQuery(new BoxCommandLineQuery(title,cellSize,_description));
    8182}
    … …  
    202203}
    203204
    204 CommandLineDialog::VectorCommandLineQuery::VectorCommandLineQuery(string title, Vector *_target, const double *const _cellSize, bool _check, string _description) :
    205     Dialog::VectorQuery(title,_target,_cellSize,_check, _description)
     205CommandLineDialog::VectorCommandLineQuery::VectorCommandLineQuery(string title, Vector *_target, bool _check, string _description) :
     206    Dialog::VectorQuery(title,_target,_check, _description)
    206207{}
    207208
    … …  
    224225
    225226
    226 CommandLineDialog::BoxCommandLineQuery::BoxCommandLineQuery(string title, double ** const _cellSize, string _description) :
     227CommandLineDialog::BoxCommandLineQuery::BoxCommandLineQuery(string title, Box* _cellSize, string _description) :
    227228    Dialog::BoxQuery(title,_cellSize, _description)
    228229{}

  • src/UIElements/CommandLineUI/CommandLineDialog.hpp

    r56fb09 r5ec8e3  
    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 *,const double * const,bool, std::string = "");
    37   virtual void queryBox(const char*,double ** const, std::string = "");
     36  virtual void queryVector(const char*,Vector *,bool, std::string = "");
     37  virtual void queryBox(const char*,Box *, 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,const double *const _cellSize,bool _check, std::string _description = "");
     93    VectorCommandLineQuery(std::string title,Vector *_target,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,double ** const _cellSize, std::string _description = "");
     100    BoxCommandLineQuery(std::string title,Box* _cellSize, std::string _description = "");
    101101    virtual ~BoxCommandLineQuery();
    102102    virtual bool handle();

  • src/UIElements/Dialog.cpp

    r56fb09 r5ec8e3  
    1414#include "molecule.hpp"
    1515#include "vector.hpp"
     16#include "Matrix.hpp"
     17#include "Box.hpp"
    1618
    1719using namespace std;
    … …  
    173175// Vector Queries
    174176
    175 Dialog::VectorQuery::VectorQuery(std::string title,Vector *_target,const double *const _cellSize,bool _check, std::string _description) :
     177Dialog::VectorQuery::VectorQuery(std::string title,Vector *_target,bool _check, std::string _description) :
    176178  Query(title, _description),
    177   cellSize(_cellSize),
    178179  check(_check),
    179180  target(_target)
    … …  
    193194// Box Queries
    194195
    195 Dialog::BoxQuery::BoxQuery(std::string title, double ** const _cellSize, std::string _description) :
     196Dialog::BoxQuery::BoxQuery(std::string title, Box* _cellSize, std::string _description) :
    196197  Query(title, _description),
    197198  target(_cellSize)
    … …  
    206207
    207208void Dialog::BoxQuery::setResult() {
    208   for (int i=0;i<6;i++) {
    209     (*target)[i] = tmp[i];
    210   }
     209  (*target)= ReturnFullMatrixforSymmetric(tmp);
    211210}
    212211

  • src/UIElements/Dialog.hpp

    r56fb09 r5ec8e3  
    1414
    1515class atom;
     16class Box;
    1617class element;
    1718class molecule;
    … …  
    4243  virtual void queryAtom(const char*,atom**,std::string = "")=0;
    4344  virtual void queryMolecule(const char*,molecule**, 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;
     45  virtual void queryVector(const char*,Vector *,bool, std::string = "")=0;
     46  virtual void queryBox(const char*,Box*, std::string = "")=0;
    4647  virtual void queryElement(const char*, std::vector<element *> *, std::string = "")=0;
    4748
    … …  
    162163  class VectorQuery : public Query {
    163164  public:
    164       VectorQuery(std::string title,Vector *_target,const double *const _cellSize,bool _check, std::string _description = "");
     165      VectorQuery(std::string title,Vector *_target,bool _check, std::string _description = "");
    165166      virtual ~VectorQuery();
    166167      virtual bool handle()=0;
    … …  
    168169    protected:
    169170      Vector *tmp;
    170       const double *const cellSize;
    171171      bool check;
    172172    private:
    … …  
    176176  class BoxQuery : public Query {
    177177  public:
    178       BoxQuery(std::string title,double ** const _cellSize, std::string _description = "");
     178      BoxQuery(std::string title,Box *_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       double **target;
     185      Box* target;
    186186  };
    187187

  • src/UIElements/QT4/QTDialog.cpp

    r56fb09 r5ec8e3  
    2929#include "molecule.hpp"
    3030#include "Descriptors/MoleculeIdDescriptor.hpp"
     31#include "Matrix.hpp"
     32#include "Box.hpp"
    3133
    3234
    … …  
    9395}
    9496
    95 void QTDialog::queryBox(char const*, double**, string){
     97void QTDialog::queryBox(char const*, Box*, string){
    9698  // TODO
    9799  ASSERT(false, "Not implemented yet");
    … …  
    118120}
    119121
    120 void QTDialog::queryVector(const char* title, Vector *target,const double *const cellSize, bool check,string) {
    121   registerQuery(new VectorQTQuery(title,target,cellSize,check,inputLayout,this));
     122void QTDialog::queryVector(const char* title, Vector *target, bool check,string) {
     123  registerQuery(new VectorQTQuery(title,target,check,inputLayout,this));
    122124}
    123125
    … …  
    248250}
    249251
    250 QTDialog::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;
     252QTDialog::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();
    256257  const char *coords[3] = {"x:","y:","z:"};
    257258  mainLayout= new QHBoxLayout();
    … …  
    261262  mainLayout->addLayout(subLayout);
    262263  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,cellSize[j]);
     269    coordInput[i]->setRange(0,M.at(i,i));
    270270    coordInput[i]->setDecimals(3);
    271271    coordLayout[i]->addWidget(coordInput[i]);

  • src/UIElements/QT4/QTDialog.hpp

    r56fb09 r5ec8e3  
    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 *,const double *const,bool,std::string = "");
    45   virtual void queryBox(const char*,double ** const, std::string = "");
     44  virtual void queryVector(const char*,Vector *,bool,std::string = "");
     45  virtual void queryBox(const char*,Box*, 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,const double *const _cellSize,bool _check,QBoxLayout *,QTDialog *);
     111      VectorQTQuery(std::string title,Vector *_target,bool _check,QBoxLayout *,QTDialog *);
    112112      virtual ~VectorQTQuery();
    113113      virtual bool handle();

  • src/UIElements/TextUI/TextDialog.cpp

    r56fb09 r5ec8e3  
    2525#include "molecule.hpp"
    2626#include "vector.hpp"
     27#include "Matrix.hpp"
     28#include "Box.hpp"
    2729
    2830using namespace std;
    … …  
    6668}
    6769
    68 void 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 
    72 void TextDialog::queryBox(const char* title,double ** const cellSize, string description) {
     70void TextDialog::queryVector(const char* title, Vector *target, bool check, string description) {
     71  registerQuery(new VectorTextQuery(title,target,check,description));
     72}
     73
     74void TextDialog::queryBox(const char* title,Box* cellSize, string description) {
    7375  registerQuery(new BoxTextQuery(title,cellSize,description));
    7476}
    … …  
    243245}
    244246
    245 TextDialog::VectorTextQuery::VectorTextQuery(std::string title, Vector *_target, const double *const _cellSize, bool _check, std::string _description) :
    246     Dialog::VectorQuery(title,_target,_cellSize,_check,_description)
     247TextDialog::VectorTextQuery::VectorTextQuery(std::string title, Vector *_target, bool _check, std::string _description) :
     248    Dialog::VectorQuery(title,_target,_check,_description)
    247249{}
    248250
    … …  
    253255  Log() << Verbose(0) << getTitle();
    254256
     257  const Matrix &M = World::getInstance().getDomain().getM();
    255258  char coords[3] = {'x','y','z'};
    256   int j = -1;
    257259  for (int i=0;i<3;i++) {
    258     j += i+1;
    259260    do {
    260       Log() << Verbose(0) << coords[i] << "[0.." << cellSize[j] << "]: ";
     261      Log() << Verbose(0) << coords[i] << "[0.." << M.at(i,i) << "]: ";
    261262      cin >> (*tmp)[i];
    262     } while ((((*tmp)[i] < 0) || ((*tmp)[i] >= cellSize[j])) && (check));
     263    } while ((((*tmp)[i] < 0) || ((*tmp)[i] >= M.at(i,i))) && (check));
    263264  }
    264265  return true;
    265266}
    266267
    267 TextDialog::BoxTextQuery::BoxTextQuery(std::string title, double ** const _cellSize, std::string _description) :
     268TextDialog::BoxTextQuery::BoxTextQuery(std::string title, Box* _cellSize, std::string _description) :
    268269    Dialog::BoxQuery(title,_cellSize,_description)
    269270{}

  • src/UIElements/TextUI/TextDialog.hpp

    r56fb09 r5ec8e3  
    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 *,const double * const,bool, std::string = "");
    34   virtual void queryBox(const char*,double ** const, std::string = "");
     33  virtual void queryVector(const char*,Vector *,bool, std::string = "");
     34  virtual void queryBox(const char*,Box*, 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,const double *const _cellSize,bool _check, std::string _description = NULL);
     90    VectorTextQuery(std::string title,Vector *_target,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,double ** const _cellSize, std::string _description = NULL);
     97    BoxTextQuery(std::string title,Box* _cellSize, std::string _description = NULL);
    9898    virtual ~BoxTextQuery();
    9999    virtual bool handle();

  • src/VectorSet.hpp

    r56fb09 r5ec8e3  
    1212#include <functional>
    1313#include <algorithm>
     14#include <limits>
    1415
    1516/**
    … …  
    1920 */
    2021
    21 class Vector;
     22#include "vector.hpp"
     23#include <list>
    2224
    2325// this tests, whether we actually have a Vector
    … …  
    4850  void translate(const Vector &translater){
    4951    // this is needed to allow template lookup
    50     transform(this->begin(),this->end(),this->begin(),bind1st(plus<Vector>(),translater));
     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());
    5164  }
    5265};
    5366
     67// allows simpler definition of VectorSets
     68#define VECTORSET(container_type) VectorSet<container_type<Vector> >
     69
    5470#endif /* VECTORSET_HPP_ */

  • src/World.cpp

    r56fb09 r5ec8e3  
    2222#include "Actions/ManipulateAtomsProcess.hpp"
    2323#include "Helpers/Assert.hpp"
     24#include "Box.hpp"
     25#include "Matrix.hpp"
    2426
    2527#include "Patterns/Singleton_impl.hpp"
    … …  
    7476// system
    7577
    76 double * World::getDomain() {
    77   return cell_size;
     78Box& World::getDomain() {
     79  return *cell_size;
     80}
     81
     82void World::setDomain(const Matrix &mat){
     83  *cell_size = mat;
    7884}
    7985
    8086void World::setDomain(double * matrix)
    8187{
    82 
     88  Matrix M = ReturnFullMatrixforSymmetric(matrix);
     89  cell_size->setM(M);
    8390}
    8491
    … …  
    133140  molecules.erase(id);
    134141}
    135 
    136 double *World::cell_size = NULL;
    137142
    138143atom *World::createAtom(){
    … …  
    297302    molecules_deprecated(new MoleculeListClass(this))
    298303{
    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.;
     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);
    306310  defaultName = "none";
    307311  molecules_deprecated->signOn(this);
    … …  
    311315{
    312316  molecules_deprecated->signOff(this);
    313   delete[] cell_size;
     317  delete cell_size;
    314318  delete molecules_deprecated;
    315319  delete periode;

  • src/World.hpp

    r56fb09 r5ec8e3  
    3434class AtomDescriptor_impl;
    3535template<typename T> class AtomsCalculation;
     36class Box;
    3637class config;
    3738class ManipulateAtomsProcess;
     39class Matrix;
    3840class molecule;
    3941class MoleculeDescriptor;
    … …  
    125127   * get the domain size as a symmetric matrix (6 components)
    126128   */
    127   double * getDomain();
     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);
    128137
    129138  /**
    … …  
    257266  periodentafel *periode;
    258267  config *configuration;
    259   static double *cell_size;
     268  Box *cell_size;
    260269  std::string defaultName;
    261270  class ThermoStatContainer *Thermostats;

  • src/analysis_correlation.cpp

    r56fb09 r5ec8e3  
    1919#include "triangleintersectionlist.hpp"
    2020#include "vector.hpp"
     21#include "Matrix.hpp"
    2122#include "verbose.hpp"
    2223#include "World.hpp"
     24#include "Box.hpp"
    2325
    2426
    … …  
    3436  PairCorrelationMap *outmap = NULL;
    3537  double distance = 0.;
    36   double *domain = World::getInstance().getDomain();
     38  Box &domain = World::getInstance().getDomain();
    3739
    3840  if (molecules->ListOfMolecules.empty()) {
    … …  
    7577                for (set <pair<element *, element *> >::iterator PairRunner = PairsOfElements.begin(); PairRunner != PairsOfElements.end(); ++PairRunner)
    7678                  if ((PairRunner->first == (**iter).type) && (PairRunner->second == (**runner).type)) {
    77                     distance = (*iter)->node->PeriodicDistance(*(*runner)->node,  domain);
     79                    distance = domain.periodicDistance(*(*iter)->node,*(*runner)->node);
    7880                    //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
    7981                    outmap->insert ( pair<double, pair <atom *, atom*> > (distance, pair<atom *, atom*> ((*iter), (*runner)) ) );
    … …  
    135137  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++){
    136138    if ((*MolWalker)->ActiveFlag) {
    137       double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
    138       double * FullInverseMatrix = InverseMatrix(FullMatrix);
     139      Matrix FullMatrix = World::getInstance().getDomain().getM();
     140      Matrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
    139141      DoeLog(2) && (eLog()<< Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    140142      eLog() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl;
    141143      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    142144        DoLog(3) && (Log() << Verbose(3) << "Current atom is " << **iter << "." << endl);
    143         periodicX = *(*iter)->node;
    144         periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
     145        periodicX = FullInverseMatrix * (*(*iter)->node); // x now in [0,1)^3
    145146        // go through every range in xyz and get distance
    146147        for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
    147148          for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
    148149            for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
    149               checkX = Vector(n[0], n[1], n[2]) + periodicX;
    150               checkX.MatrixMultiplication(FullMatrix);
     150              checkX = FullMatrix * (Vector(n[0], n[1], n[2]) + periodicX);
    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 = *(*runner)->node;
    160                           periodicOtherX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
     159                          periodicOtherX = FullInverseMatrix * (*(*runner)->node); // x now in [0,1)^3
    161160                          // go through every range in xyz and get distance
    162161                          for (Othern[0]=-ranges[0]; Othern[0] <= ranges[0]; Othern[0]++)
    163162                            for (Othern[1]=-ranges[1]; Othern[1] <= ranges[1]; Othern[1]++)
    164163                              for (Othern[2]=-ranges[2]; Othern[2] <= ranges[2]; Othern[2]++) {
    165                                 checkOtherX = Vector(Othern[0], Othern[1], Othern[2]) + periodicOtherX;
    166                                 checkOtherX.MatrixMultiplication(FullMatrix);
     164                                checkOtherX = FullMatrix * (Vector(Othern[0], Othern[1], Othern[2]) + periodicOtherX);
    167165                                distance = checkX.distance(checkOtherX);
    168166                                //Log() << Verbose(1) <<"Inserting " << *(*iter) << " and " << *(*runner) << endl;
    … …  
    176174            }
    177175      }
    178       delete[](FullMatrix);
    179       delete[](FullInverseMatrix);
    180176    }
    181177  }
    … …  
    195191  CorrelationToPointMap *outmap = NULL;
    196192  double distance = 0.;
    197   double *cell_size = World::getInstance().getDomain();
     193  Box &domain = World::getInstance().getDomain();
    198194
    199195  if (molecules->ListOfMolecules.empty()) {
    … …  
    211207        for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
    212208          if ((*type == NULL) || ((*iter)->type == *type)) {
    213             distance = (*iter)->node->PeriodicDistance(*point, cell_size);
     209            distance = domain.periodicDistance(*(*iter)->node,*point);
    214210            DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
    215211            outmap->insert ( pair<double, pair<atom *, const Vector*> >(distance, pair<atom *, const Vector*> ((*iter), point) ) );
    … …  
    246242  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    247243    if ((*MolWalker)->ActiveFlag) {
    248       double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
    249       double * FullInverseMatrix = InverseMatrix(FullMatrix);
     244      Matrix FullMatrix = World::getInstance().getDomain().getM();
     245      Matrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
    250246      DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    251247      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    … …  
    253249        for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
    254250          if ((*type == NULL) || ((*iter)->type == *type)) {
    255             periodicX = *(*iter)->node;
    256             periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
     251            periodicX = FullInverseMatrix * (*(*iter)->node); // x now in [0,1)^3
    257252            // go through every range in xyz and get distance
    258253            for (n[0]=-ranges[0]; n[0] <= ranges[0]; n[0]++)
    259254              for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
    260255                for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
    261                   checkX = Vector(n[0], n[1], n[2]) + periodicX;
    262                   checkX.MatrixMultiplication(FullMatrix);
     256                  checkX = FullMatrix * (Vector(n[0], n[1], n[2]) + periodicX);
    263257                  distance = checkX.distance(*point);
    264258                  DoLog(4) && (Log() << Verbose(4) << "Current distance is " << distance << "." << endl);
    … …  
    267261          }
    268262      }
    269       delete[](FullMatrix);
    270       delete[](FullInverseMatrix);
    271263    }
    272264
    … …  
    352344  for (MoleculeList::const_iterator MolWalker = molecules->ListOfMolecules.begin(); MolWalker != molecules->ListOfMolecules.end(); MolWalker++)
    353345    if ((*MolWalker)->ActiveFlag) {
    354       double * FullMatrix = ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
    355       double * FullInverseMatrix = InverseMatrix(FullMatrix);
     346      Matrix FullMatrix = World::getInstance().getDomain().getM();
     347      Matrix FullInverseMatrix = World::getInstance().getDomain().getMinv();
    356348      DoLog(2) && (Log() << Verbose(2) << "Current molecule is " << *MolWalker << "." << endl);
    357349      for (molecule::const_iterator iter = (*MolWalker)->begin(); iter != (*MolWalker)->end(); ++iter) {
    … …  
    359351        for (vector<element *>::const_iterator type = elements.begin(); type != elements.end(); ++type)
    360352          if ((*type == NULL) || ((*iter)->type == *type)) {
    361             periodicX = *(*iter)->node;
    362             periodicX.MatrixMultiplication(FullInverseMatrix);  // x now in [0,1)^3
     353            periodicX = FullInverseMatrix * (*(*iter)->node); // x now in [0,1)^3
    363354            // go through every range in xyz and get distance
    364355            ShortestDistance = -1.;
    … …  
    366357              for (n[1]=-ranges[1]; n[1] <= ranges[1]; n[1]++)
    367358                for (n[2]=-ranges[2]; n[2] <= ranges[2]; n[2]++) {
    368                   checkX = Vector(n[0], n[1], n[2]) + periodicX;
    369                   checkX.MatrixMultiplication(FullMatrix);
     359                  checkX = FullMatrix * (Vector(n[0], n[1], n[2]) + periodicX);
    370360                  TriangleIntersectionList Intersections(&checkX,Surface,LC);
    371361                  distance = Intersections.GetSmallestDistance();
    … …  
    381371          }
    382372      }
    383       delete[](FullMatrix);
    384       delete[](FullInverseMatrix);
    385373    }
    386374

  • src/atom.cpp

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

  • src/atom.hpp

    r56fb09 r5ec8e3  
    3535class World;
    3636class molecule;
     37class Shape;
    3738
    3839/********************************************** declarations *******************************/
    … …  
    6667  double DistanceToVector(const Vector &origin) const;
    6768  double DistanceSquaredToVector(const Vector &origin) const;
    68   bool IsInParallelepiped(const Vector offset, const double *parallelepiped) const;
     69  bool IsInShape(const Shape&) const;
    6970
    7071  // getter and setter

  • src/boundary.cpp

    r56fb09 r5ec8e3  
    2222#include "World.hpp"
    2323#include "Plane.hpp"
     24#include "Matrix.hpp"
     25#include "Box.hpp"
    2426
    2527#include<gsl/gsl_poly.h>
    … …  
    764766  int N[NDIM];
    765767  int n[NDIM];
    766   double *M =  ReturnFullMatrixforSymmetric(World::getInstance().getDomain());
    767   double Rotations[NDIM*NDIM];
    768   double *MInverse = InverseMatrix(M);
     768  const Matrix &M = World::getInstance().getDomain().getM();
     769  Matrix Rotations;
     770  const Matrix &MInverse = World::getInstance().getDomain().getMinv();
    769771  Vector AtomTranslations;
    770772  Vector FillerTranslations;
    … …  
    799801
    800802  // calculate filler grid in [0,1]^3
    801   FillerDistance = Vector(distance[0], distance[1], distance[2]);
    802   FillerDistance.InverseMatrixMultiplication(M);
     803  FillerDistance = MInverse * Vector(distance[0], distance[1], distance[2]);
    803804  for(int i=0;i<NDIM;i++)
    804805    N[i] = (int) ceil(1./FillerDistance[i]);
    … …  
    813814      for (n[2] = 0; n[2] < N[2]; n[2]++) {
    814815        // calculate position of current grid vector in untransformed box
    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);
     816        CurrentPosition = M * Vector((double)n[0]/(double)N[0], (double)n[1]/(double)N[1], (double)n[2]/(double)N[2]);
    817817        // create molecule random translation vector ...
    818818        for (int i=0;i<NDIM;i++)
    … …  
    835835            }
    836836
    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])                                     ;
     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])                                        );
    846846          }
    847847
    … …  
    849849          Inserter = (*iter)->x;
    850850          if (DoRandomRotation)
    851             Inserter.MatrixMultiplication(Rotations);
     851            Inserter *= Rotations;
    852852          Inserter += AtomTranslations + FillerTranslations + CurrentPosition;
    853853
    854854          // check whether inserter is inside box
    855           Inserter.MatrixMultiplication(MInverse);
     855          Inserter *= MInverse;
    856856          FillIt = true;
    857857          for (int i=0;i<NDIM;i++)
    858858            FillIt = FillIt && (Inserter[i] >= -MYEPSILON) && ((Inserter[i]-1.) <= MYEPSILON);
    859           Inserter.MatrixMultiplication(M);
     859          Inserter *= M;
    860860
    861861          // Check whether point is in- or outside
    … …  
    892892            }
    893893      }
    894   delete[](M);
    895   delete[](MInverse);
    896894
    897895  return Filling;

  • src/config.cpp

    r56fb09 r5ec8e3  
    2626#include "ThermoStatContainer.hpp"
    2727#include "World.hpp"
     28#include "Matrix.hpp"
     29#include "Box.hpp"
    2830
    2931/************************************* Functions for class config ***************************/
    … …  
    679681  // Unit cell and magnetic field
    680682  ParseForParameter(verbose,FileBuffer, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
    681   double * const cell_size = World::getInstance().getDomain();
     683  double * cell_size = new double[6];
    682684  cell_size[0] = BoxLength[0];
    683685  cell_size[1] = BoxLength[3];
    … …  
    686688  cell_size[4] = BoxLength[7];
    687689  cell_size[5] = BoxLength[8];
     690  World::getInstance().setDomain(cell_size);
     691  delete cell_size;
    688692  //if (1) fprintf(stderr,"\n");
    689693
    … …  
    883887
    884888  ParseForParameter(verbose,file, "BoxLength", 0, 3, 3, lower_trigrid, BoxLength, 1, critical); /* Lattice->RealBasis */
    885   double * const cell_size = World::getInstance().getDomain();
     889  double * cell_size = new double[6];
    886890  cell_size[0] = BoxLength[0];
    887891  cell_size[1] = BoxLength[3];
    … …  
    890894  cell_size[4] = BoxLength[7];
    891895  cell_size[5] = BoxLength[8];
     896  World::getInstance().setDomain(cell_size);
     897  delete[] cell_size;
    892898  if (1) fprintf(stderr,"\n");
    893899  config::DoPerturbation = 0;
    … …  
    10271033  // bring MaxTypes up to date
    10281034  mol->CountElements();
    1029   const double * const cell_size = World::getInstance().getDomain();
     1035  const Matrix &domain = World::getInstance().getDomain().getM();
    10301036  ofstream * const output = new ofstream(filename, ios::out);
    10311037  if (output != NULL) {
    … …  
    10981104    *output << endl;
    10991105    *output << "BoxLength\t\t\t# (Length of a unit cell)" << 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;
     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;
    11031109    // FIXME
    11041110    *output << endl;

  • src/ellipsoid.cpp

    r56fb09 r5ec8e3  
    2121#include "tesselation.hpp"
    2222#include "vector.hpp"
     23#include "Matrix.hpp"
    2324#include "verbose.hpp"
    2425
    … …  
    3435  Vector helper, RefPoint;
    3536  double distance = -1.;
    36   double Matrix[NDIM*NDIM];
     37  Matrix Matrix;
    3738  double InverseLength[3];
    3839  double psi,theta,phi; // euler angles in ZX'Z'' convention
    … …  
    5152  theta = EllipsoidAngle[1];
    5253  phi = EllipsoidAngle[2];
    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);
     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;
    6364  helper.ScaleAll(InverseLength);
    6465  //Log() << Verbose(4) << "Transformed RefPoint is at " << helper << "." << endl;
    … …  
    7374  phi = -EllipsoidAngle[2];
    7475  helper.ScaleAll(EllipsoidLength);
    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);
     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;
    8586  //Log() << Verbose(4) << "Intersection is at " << helper << "." << endl;
    8687

  • src/helpers.cpp

    r56fb09 r5ec8e3  
    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  */
    123 double * 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  */
    141 double * 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 
    167119/** Comparison function for GSL heapsort on distances in two molecules.
    168120 * \param *a

  • src/helpers.hpp

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

  • src/molecule.cpp

    r56fb09 r5ec8e3  
    99#include <cstring>
    1010#include <boost/bind.hpp>
     11#include <boost/foreach.hpp>
    1112
    1213#include "World.hpp"
    … …  
    2728#include "tesselation.hpp"
    2829#include "vector.hpp"
     30#include "Matrix.hpp"
    2931#include "World.hpp"
     32#include "Box.hpp"
    3033#include "Plane.hpp"
    3134#include "Exceptions/LinearDependenceException.hpp"
    … …  
    284287  Vector Orthovector1, Orthovector2;  // temporary vectors in coordination construction
    285288  Vector InBondvector;    // vector in direction of *Bond
    286   double *matrix = NULL;
     289  const Matrix &matrix =  World::getInstance().getDomain().getM();
    287290  bond *Binder = NULL;
    288   double * const cell_size = World::getInstance().getDomain();
    289291
    290292//  Log() << Verbose(3) << "Begin of AddHydrogenReplacementAtom." << endl;
    … …  
    307309      } // (signs are correct, was tested!)
    308310    }
    309     matrix = ReturnFullMatrixforSymmetric(cell_size);
    310     Orthovector1.MatrixMultiplication(matrix);
     311    Orthovector1 *= matrix;
    311312    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);
    544543
    545544//  Log() << Verbose(3) << "End of AddHydrogenReplacementAtom." << endl;
    … …  
    660659 * @param three vectors forming the matrix that defines the shape of the parallelpiped
    661660 */
    662 molecule* molecule::CopyMoleculeFromSubRegion(const Vector offset, const double *parallelepiped) const {
     661molecule* molecule::CopyMoleculeFromSubRegion(const Shape &region) const {
    663662  molecule *copy = World::getInstance().createMolecule();
    664663
    665   ActOnCopyWithEachAtomIfTrue ( &molecule::AddCopyAtom, copy, &atom::IsInParallelepiped, offset, parallelepiped );
     664  BOOST_FOREACH(atom *iter,atoms){
     665    if(iter->IsInShape(region)){
     666      copy->AddCopyAtom(iter);
     667    }
     668  }
    666669
    667670  //TODO: copy->BuildInducedSubgraph(this);
    … …  
    750753void molecule::SetBoxDimension(Vector *dim)
    751754{
    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);
     755  Matrix domain;
     756  for(int i =0; i<NDIM;++i)
     757    domain.at(i,i) = dim->at(i);
     758  World::getInstance().setDomain(domain);
    759759};
    760760
    … …  
    849849bool molecule::CheckBounds(const Vector *x) const
    850850{
    851   double * const cell_size = World::getInstance().getDomain();
     851  const Matrix &domain = World::getInstance().getDomain().getM();
    852852  bool result = true;
    853   int j =-1;
    854853  for (int i=0;i<NDIM;i++) {
    855     j += i+1;
    856     result = result && ((x->at(i) >= 0) && (x->at(i) < cell_size[j]));
     854    result = result && ((x->at(i) >= 0) && (x->at(i) < domain.at(i,i)));
    857855  }
    858856  //return result;

  • src/molecule.hpp

    r56fb09 r5ec8e3  
    5353class periodentafel;
    5454class Vector;
     55class Shape;
    5556
    5657/******************************** Some definitions for easier reading **********************************/
    … …  
    316317
    317318  molecule *CopyMolecule();
    318   molecule* CopyMoleculeFromSubRegion(const Vector offset, const double *parallelepiped) const;
     319  molecule* CopyMoleculeFromSubRegion(const Shape&) const;
    319320
    320321  /// Fragment molecule by two different approaches:

  • src/molecule_fragmentation.cpp

    r56fb09 r5ec8e3  
    2222#include "periodentafel.hpp"
    2323#include "World.hpp"
     24#include "Matrix.hpp"
     25#include "Box.hpp"
    2426
    2527/************************************* Functions for class molecule *********************************/
    … …  
    17081710  atom *Walker = NULL;
    17091711  atom *OtherWalker = NULL;
    1710   double * const cell_size = World::getInstance().getDomain();
    1711   double *matrix = ReturnFullMatrixforSymmetric(cell_size);
     1712  Matrix matrix = World::getInstance().getDomain().getM();
    17121713  enum Shading *ColorList = NULL;
    17131714  double tmp;
    … …  
    17491750          Translationvector[i] = (tmp < 0) ? +1. : -1.;
    17501751      }
    1751       Translationvector.MatrixMultiplication(matrix);
     1752      Translationvector *= matrix;
    17521753      //Log() << Verbose(3) << "Translation vector is ";
    17531754      Log() << Verbose(0) << Translationvector <<  endl;
    … …  
    17801781  delete(AtomStack);
    17811782  delete[](ColorList);
    1782   delete[](matrix);
    17831783  DoLog(2) && (Log() << Verbose(2) << "End of ScanForPeriodicCorrection." << endl);
    17841784};

  • src/molecule_geometry.cpp

    r56fb09 r5ec8e3  
    1919#include "World.hpp"
    2020#include "Plane.hpp"
     21#include "Matrix.hpp"
     22#include "Box.hpp"
    2123#include <boost/foreach.hpp>
    2224
    … …  
    3335  const Vector *Center = DetermineCenterOfAll();
    3436  const Vector *CenterBox = DetermineCenterOfBox();
    35   double * const cell_size = World::getInstance().getDomain();
    36   double *M = ReturnFullMatrixforSymmetric(cell_size);
    37   double *Minv = InverseMatrix(M);
     37  Box &domain = World::getInstance().getDomain();
    3838
    3939  // go through all atoms
    4040  ActOnAllVectors( &Vector::SubtractVector, *Center);
    4141  ActOnAllVectors( &Vector::SubtractVector, *CenterBox);
    42   ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
    43 
    44   delete[](M);
    45   delete[](Minv);
     42  BOOST_FOREACH(atom* iter, atoms){
     43    *iter->node = domain.WrapPeriodically(*iter->node);
     44  }
     45
    4646  delete(Center);
    4747  delete(CenterBox);
    … …  
    5656{
    5757  bool status = true;
    58   double * const cell_size = World::getInstance().getDomain();
    59   double *M = ReturnFullMatrixforSymmetric(cell_size);
    60   double *Minv = InverseMatrix(M);
     58  Box &domain = World::getInstance().getDomain();
    6159
    6260  // go through all atoms
    63   ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
    64 
    65   delete[](M);
    66   delete[](Minv);
     61  BOOST_FOREACH(atom* iter, atoms){
     62    *iter->node = domain.WrapPeriodically(*iter->node);
     63  }
     64
    6765  return status;
    6866};
    … …  
    153151{
    154152  Vector *a = new Vector(0.5,0.5,0.5);
    155 
    156   const double *cell_size = World::getInstance().getDomain();
    157   double *M = ReturnFullMatrixforSymmetric(cell_size);
    158   a->MatrixMultiplication(M);
    159   delete[](M);
    160 
     153  const Matrix &M = World::getInstance().getDomain().getM();
     154  (*a) *= M;
    161155  return a;
    162156};
    … …  
    244238void molecule::TranslatePeriodically(const Vector *trans)
    245239{
    246   double * const cell_size = World::getInstance().getDomain();
    247   double *M = ReturnFullMatrixforSymmetric(cell_size);
    248   double *Minv = InverseMatrix(M);
     240  Box &domain = World::getInstance().getDomain();
    249241
    250242  // go through all atoms
    251243  ActOnAllVectors( &Vector::AddVector, *trans);
    252   ActOnAllVectors( &Vector::WrapPeriodically, (const double *)M, (const double *)Minv);
    253 
    254   delete[](M);
    255   delete[](Minv);
     244  BOOST_FOREACH(atom* iter, atoms){
     245    *iter->node = domain.WrapPeriodically(*iter->node);
     246  }
     247
    256248};
    257249
    … …  
    264256  OBSERVE;
    265257  Plane p(*n,0);
    266   BOOST_FOREACH( atom* iter, atoms ){
     258  BOOST_FOREACH(atom* iter, atoms ){
    267259    (*iter->node) = p.mirrorVector(*iter->node);
    268260  }
    … …  
    274266void molecule::DeterminePeriodicCenter(Vector &center)
    275267{
    276   double * const cell_size = World::getInstance().getDomain();
    277   double *matrix = ReturnFullMatrixforSymmetric(cell_size);
    278   double *inversematrix = InverseMatrix(matrix);
     268  const Matrix &matrix = World::getInstance().getDomain().getM();
     269  const Matrix &inversematrix = World::getInstance().getDomain().getM();
    279270  double tmp;
    280271  bool flag;
    … …  
    288279      if ((*iter)->type->Z != 1) {
    289280#endif
    290         Testvector = (*iter)->x;
    291         Testvector.MatrixMultiplication(inversematrix);
     281        Testvector = inversematrix * (*iter)->x;
    292282        Translationvector.Zero();
    293283        for (BondList::const_iterator Runner = (*iter)->ListOfBonds.begin(); Runner != (*iter)->ListOfBonds.end(); (++Runner)) {
    … …  
    306296        }
    307297        Testvector += Translationvector;
    308         Testvector.MatrixMultiplication(matrix);
     298        Testvector *= matrix;
    309299        Center += Testvector;
    310300        Log() << Verbose(1) << "vector is: " << Testvector << endl;
    … …  
    313303        for (BondList::const_iterator Runner = (*iter)->ListOfBonds.begin(); Runner != (*iter)->ListOfBonds.end(); (++Runner)) {
    314304          if ((*Runner)->GetOtherAtom((*iter))->type->Z == 1) {
    315             Testvector = (*Runner)->GetOtherAtom((*iter))->x;
    316             Testvector.MatrixMultiplication(inversematrix);
     305            Testvector = inversematrix * (*Runner)->GetOtherAtom((*iter))->x;
    317306            Testvector += Translationvector;
    318             Testvector.MatrixMultiplication(matrix);
     307            Testvector *= matrix;
    319308            Center += Testvector;
    320309            Log() << Verbose(1) << "Hydrogen vector is: " << Testvector << endl;
    … …  
    325314    }
    326315  } while (!flag);
    327   delete[](matrix);
    328   delete[](inversematrix);
    329316
    330317  Center.Scale(1./static_cast<double>(getAtomCount()));
    … …  
    388375    DoLog(1) && (Log() << Verbose(1) << "Transforming molecule into PAS ... ");
    389376    // the eigenvectors specify the transformation matrix
    390     ActOnAllVectors( &Vector::MatrixMultiplication, (const double *) evec->data );
     377    Matrix M = Matrix(evec->data);
     378    BOOST_FOREACH(atom* iter, atoms){
     379      (*iter->node) *= M;
     380    }
    391381    DoLog(0) && (Log() << Verbose(0) << "done." << endl);
    392382

  • src/molecule_graph.cpp

    r56fb09 r5ec8e3  
    2424#include "Helpers/fast_functions.hpp"
    2525#include "Helpers/Assert.hpp"
     26#include "Matrix.hpp"
     27#include "Box.hpp"
    2628
    2729
    … …  
    121123  LinkedCell *LC = NULL;
    122124  bool free_BG = false;
    123   double * const cell_size = World::getInstance().getDomain();
     125  Box &domain = World::getInstance().getDomain();
    124126
    125127  if (BG == NULL) {
    … …  
    178180                          //Log() << Verbose(1) << "Checking distance " << OtherWalker->x.PeriodicDistanceSquared(&(Walker->x), cell_size) << " against typical bond length of " << bonddistance*bonddistance << "." << endl;
    179181                          (BG->*minmaxdistance)(Walker, OtherWalker, MinDistance, MaxDistance, IsAngstroem);
    180                           const double distance = OtherWalker->x.PeriodicDistanceSquared(Walker->x,cell_size);
     182                          const double distance = domain.periodicDistanceSquared(OtherWalker->x,Walker->x);
    181183                          const bool status = (distance <= MaxDistance * MaxDistance) && (distance >= MinDistance * MinDistance);
    182184//                          Log() << Verbose(1) << "MinDistance is " << MinDistance << " and MaxDistance is " << MaxDistance << "." << endl;

  • src/moleculelist.cpp

    r56fb09 r5ec8e3  
    2424#include "periodentafel.hpp"
    2525#include "Helpers/Assert.hpp"
     26#include "Matrix.hpp"
     27#include "Box.hpp"
    2628
    2729#include "Helpers/Assert.hpp"
    … …  
    631633  int FragmentCounter = 0;
    632634  ofstream output;
    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];
     635  Matrix cell_size = World::getInstance().getDomain().getM();
     636  Matrix cell_size_backup = cell_size;
     637
    639638  // store the fragments as config and as xyz
    640639  for (MoleculeList::iterator ListRunner = ListOfMolecules.begin(); ListRunner != ListOfMolecules.end(); ListRunner++) {
    … …  
    674673    (*ListRunner)->CenterEdge(&BoxDimension);
    675674    (*ListRunner)->SetBoxDimension(&BoxDimension); // update Box of atoms by boundary
    676     int j = -1;
    677675    for (int k = 0; k < NDIM; k++) {
    678       j += k + 1;
    679676      BoxDimension[k] = 2.5 * (World::getInstance().getConfig()->GetIsAngstroem() ? 1. : 1. / AtomicLengthToAngstroem);
    680       cell_size[j] = BoxDimension[k] * 2.;
    681     }
     677      cell_size.at(k,k) = BoxDimension[k] * 2.;
     678    }
     679    World::getInstance().setDomain(cell_size);
    682680    (*ListRunner)->Translate(&BoxDimension);
    683681
    … …  
    725723
    726724  // restore cell_size
    727   for (int i=0;i<6;i++)
    728     cell_size[i] = cell_size_backup[i];
     725  World::getInstance().setDomain(cell_size_backup);
    729726
    730727  return result;
    … …  
    887884  // center at set box dimensions
    888885  mol->CenterEdge(&center);
    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];
     886  Matrix domain;
     887  for(int i =0;i<NDIM;++i)
     888    domain.at(i,i) = center[i];
     889  World::getInstance().setDomain(domain);
    895890  insert(mol);
    896891}

  • src/unittests/LineUnittest.cpp

    r56fb09 r5ec8e3  
    1717
    1818#include <iostream>
     19#include <cmath>
    1920
    2021using namespace std;
    … …  
    352353  CPPUNIT_ASSERT_EQUAL(fixture,zeroVec);
    353354}
     355
     356void 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

    r56fb09 r5ec8e3  
    2222  CPPUNIT_TEST ( intersectionTest );
    2323  CPPUNIT_TEST ( rotationTest );
     24  CPPUNIT_TEST ( sphereIntersectionTest );
    2425  CPPUNIT_TEST_SUITE_END();
    2526
    … …  
    3334  void intersectionTest();
    3435  void rotationTest();
     36  void sphereIntersectionTest();
    3537
    3638private:

  • src/unittests/Makefile.am

    r56fb09 r5ec8e3  
    3838  periodentafelTest \
    3939  PlaneUnittest \
     40  ShapeUnittest \
    4041  SingletonTest \
    4142  StackClassUnitTest \
    … …  
    8384  periodentafelTest.cpp \
    8485  PlaneUnittest.cpp \
     86  ShapeUnittest.cpp \
    8587  SingletonTest.cpp \
    8688  stackclassunittest.cpp \
    … …  
    210212PlaneUnittest_LDADD = ${ALLLIBS}
    211213
     214ShapeUnittest_SOURCES = UnitTestMain.cpp ShapeUnittest.cpp ShapeUnittest.hpp
     215ShapeUnittest_LDADD = ${ALLLIBS}
     216
    212217SingletonTest_SOURCES = UnitTestMain.cpp SingletonTest.cpp SingletonTest.hpp
    213218SingletonTest_LDADD = ${ALLLIBS} $(BOOST_LIB) ${BOOST_THREAD_LIB}

  • src/unittests/PlaneUnittest.cpp

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

  • src/unittests/linearsystemofequationsunittest.cpp

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

  • src/unittests/vectorunittest.cpp

    r56fb09 r5ec8e3  
    2121#include "Plane.hpp"
    2222#include "Exceptions/LinearDependenceException.hpp"
     23#include "Matrix.hpp"
    2324
    2425#ifdef HAVE_TESTRUNNER
    … …  
    214215  CPPUNIT_ASSERT(testVector.ScalarProduct(three) < MYEPSILON);
    215216}
    216 
    217 
    218 /**
    219  * UnitTest for Vector::IsInParallelepiped().
    220  */
    221 void 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

    r56fb09 r5ec8e3  
    2727    CPPUNIT_TEST ( ProjectionTest );
    2828    CPPUNIT_TEST ( NormalsTest );
    29     CPPUNIT_TEST ( IsInParallelepipedTest );
    3029    CPPUNIT_TEST_SUITE_END();
    3130
    … …  
    4544    void LineIntersectionTest();
    4645    void VectorRotationTest();
    47     void IsInParallelepipedTest();
    4846
    4947private:

  • src/vector.cpp

    r56fb09 r5ec8e3  
    1212#include "Helpers/Assert.hpp"
    1313#include "Helpers/fast_functions.hpp"
     14#include "Exceptions/MathException.hpp"
    1415
    1516#include <iostream>
     17#include <gsl/gsl_blas.h>
     18
    1619
    1720using namespace std;
    … …  
    3437{
    3538  content = gsl_vector_alloc(NDIM);
    36   gsl_vector_set(content,0,src[0]);
    37   gsl_vector_set(content,1,src[1]);
    38   gsl_vector_set(content,2,src[2]);
     39  gsl_vector_memcpy(content, src.content);
    3940}
    4041
    … …  
    4950};
    5051
     52Vector::Vector(gsl_vector *_content) :
     53  content(_content)
     54{}
     55
    5156/**
    5257 * Assignment operator
    … …  
    5560  // check for self assignment
    5661  if(&src!=this){
    57     gsl_vector_set(content,0,src[0]);
    58     gsl_vector_set(content,1,src[1]);
    59     gsl_vector_set(content,2,src[2]);
     62    gsl_vector_memcpy(content, src.content);
    6063  }
    6164  return *this;
    … …  
    9497}
    9598
    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  */
    101 double 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  */
    138 double 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  */
    174 void 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 
    19899/** Calculates scalar product between this and another vector.
    199100 * \param *y array to second vector
    … …  
    203104{
    204105  double res = 0.;
    205   for (int i=NDIM;i--;)
    206     res += at(i)*y[i];
     106  gsl_blas_ddot(content, y.content, &res);
    207107  return (res);
    208108};
    … …  
    504404};
    505405
     406void Vector::ScaleAll(const Vector &factor){
     407  gsl_vector_mul(content, factor.content);
     408}
    506409
    507410
    508411void Vector::Scale(const double factor)
    509412{
    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  */
    518 void 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);
     413  gsl_vector_scale(content,factor);
    530414};
    531415
    … …  
    546430  return make_pair(res,helper);
    547431}
    548 
    549 /** Do a matrix multiplication.
    550  * \param *matrix NDIM_NDIM array
    551  */
    552 void 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  */
    565 bool 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 
    592432
    593433/** Creates this vector as the b y *factors' components scaled linear combination of the given three.
    … …  
    679519void Vector::AddVector(const Vector &y)
    680520{
    681   for(int i=NDIM;i--;)
    682     at(i) += y[i];
     521  gsl_vector_add(content, y.content);
    683522}
    684523
    … …  
    688527void Vector::SubtractVector(const Vector &y)
    689528{
    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  */
    701 bool 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;
     529  gsl_vector_sub(content, y.content);
    711530}
    712531

  • src/vector.hpp

    r56fb09 r5ec8e3  
    1313#include <iostream>
    1414#include <gsl/gsl_vector.h>
    15 #include <gsl/gsl_multimin.h>
    1615
    1716#include <memory>
    … …  
    2423
    2524class Vector;
     25class Matrix;
    2626
    2727typedef std::vector<Vector> pointset;
    … …  
    3131 */
    3232class Vector : public Space{
     33  friend Vector operator*(const Matrix&,const Vector&);
    3334public:
    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;
    4644  double ScalarProduct(const Vector &y) const;
    4745  double Angle(const Vector &y) const;
    … …  
    5856  Vector Projection(const Vector &y) const;
    5957  void ScaleAll(const double *factor);
     58  void ScaleAll(const Vector &factor);
    6059  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);
    6460  bool GetOneNormalVector(const Vector &x1);
    6561  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);
    6862  std::pair<Vector,Vector> partition(const Vector&) const;
    6963  std::pair<pointset,Vector> partition(const pointset&) const;
    … …  
    10498
    10599private:
     100  Vector(gsl_vector *);
    106101  gsl_vector *content;
    107102

  • src/vector_ops.cpp

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

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

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